subroutine dh (sal1d, Fx, Fmix, Fneti, ultnt, condb, Tw, Ts, Ti
     &,              Tbot, Tf, qi, hi, hs, dhib, subi, subs, dhit, dhs
     &,              dhif, dhsf, ni)

#if defined O_ice && defined O_ice_cpts && defined O_embm

! Energy-conserving sea ice model
! Routines to grow/melt ice and adjust temperature profile

! See Bitz, C.M., and W.H. Lipscomb, 1999:
! A new energy-conserving sea ice model for climate study,
! J. Geophys. Res., in press.

! The author grants permission to the public to copy and use this
! software without charge, provided that this Notice and any statement
! of authorship are reproduced on all copies.

!     Computes the thickness changes at the top and bottom
!     and adjusts layer specific enthalpy
!     This algorithm does not allow h<0

!     Fx= actual flux of heat from the mixed layer under sea ice
!     (equal to Fmix unless all the ice melts away)
!     parameterized according to the
!     Tw as described in papers by McPhee/Steel
!     compensates for rare case of melting entire slab through

      implicit none

      include "cpts.h"
      include "thermo.h"
      include "cembm.h"

! input variables
      real sal1d(NMAX) ! ice salinity (ppt)
      real Ts          ! top srf temp (C)
      real Ti(0:NMAX)  ! snow/ice internal temp (C)
      real Tbot        ! ice bottom in (C)
      real Tw          ! water temperature below ice (C)
      real Tf          ! freezing temperature of sea water
      real hi          ! initial ice thickness (m)
      real hs          ! initial snow thickness (m)
      real Fmix        ! Flux from mixed layer from relaxation (W/m**2)
      real Fneti       ! net flux at upper surface (W/m**2)
      real condb       ! conductive flux at bottom (w/m**2)
      real ultnt       ! latent heat flux, for sublimation (w/m**2)
      integer ni       ! number of layers

! output variables
!     thickness changes from grow/melt (default) or sublimate/flooding
      real dhib         ! ice bot, dhib<0 if melt (m)
      real dhit         ! ice top, dhit<=0 (m)
      real dhs          ! snow top, dhit<=0 (m)
      real subi         ! ice top, subi<0 if sublimating (m)
      real subs         ! snow, subs<0 if sublimating (m)
      real dhif         ! ice top from flooding, dhif>0 (m)
      real dhsf         ! snow from flooding, dhsf>0 (m)

      real qi(NMAX)     ! energy of melt of ice per unit vol. (J/m**3)
      real Fx           ! actual flx of heat used from ocn (w/m**2)

! local variables
      real delti(NMAX)  ! evolving ice layer thickness
      real delts        ! evolving snow thickness
      real sumr         ! dummy for adjusting delti
      real rmvi         ! another dummy for adjusting delti
      real dtop         ! dhit+subi for adjust

      real qigrow       ! energy of melt of ice that grows
      real Ebot         ! heat available to grow/melt at bottom
      real Etop         ! heat available to melt at top

      real qiflood      ! energy of melt of flooded ice (W/m**2)
      real zintfc       ! height of snow/ice interf. wrt ocn (m)
      real qs           ! energy of melting of snow per unit volume
      real Enet         ! sum of the energy of melting relative to melt
      real Evnet        ! sum of the energy of melting relative to vapour
      real energ        ! function
      integer layer
      logical alarm2,alarm

      alarm2=.false.
      alarm=.false.

      dhib  = 0.
      dhit  = 0.
      dhs   = 0.
      dhif  = 0.
      dhsf  = 0.
      subi  = 0.
      subs  = 0.

! thickness of snow and each ice layer
      delti(1) = hi/ni
      do layer = 2,ni
        delti(layer) = delti(1)
      enddo
      delts = hs

! energy of melt per unit vol for snow and ice for each layer and sum
      qs = -RFLSNO
      enet = 0.
      do layer = 1,ni
        qi(layer) = energ(Ti(layer),sal1d(layer))
        enet = enet + qi(layer)
      enddo
      enet = enet*delti(1) + qs*hs

      if (ultnt.gt.0.0) then
! etop is positive is sublimating and negative if condensing
         Etop  = ultnt*dtau
         Evnet = -RVLSNO*hs-RVLICE*hi+Enet+Etop
         if ( Evnet .ge. 0.0 ) then
            subi=-hi
            subs=-hs
            Fx=condb-Evnet/dtau  ! maybe wrong, irrelevant though
# if defined O_ice_cpts_warnings
            print*,'sublimated though all snow and ice',Fx,hi,hs
            stop 'sublimate too much is not allowed, model will stop'
!            return
# endif
         endif
         call srfsub(    qi,   qs, delti, delts,    ni,
     $                 subi, subs,  etop,  enet,  alarm2 )

! adjust the layer thickness to reflect subl/cond changes
         delts = hs + subs
         rmvi = subi
         do layer = 1,ni
           sumr = max( -delti(layer), rmvi )
           rmvi = rmvi-sumr
           delti(layer) = delti(layer) + sumr
         enddo
      endif

      if ( Fneti .gt. 0.0 ) then
!     it is possible to melt at srf when Ts is below melting
!     but it never melts more than a negligible amount
         Etop=Fneti*dtau
         Enet=Enet+Etop
         if (Enet.ge.0.0) then
            dhit=-(hi+subi)
            dhs=-(hs+subs)
            Fx=condb-Enet/dtau
            return
         endif
         call srfmelt (qi, qs, delti, delts, ni,
     $                 dhit, dhs, etop, alarm2)
!       adjust the layer thickness to reflect melt changes
        delts = delts + dhs
        rmvi = dhit
        do layer = 1,ni
          sumr = max( -delti(layer), rmvi )
          rmvi = rmvi - sumr
          delti(layer) = delti(layer) + sumr
        enddo
      endif

      dtop=dhit+subi

      Fx=Fmix
      Ebot=dtau*(Fmix-condb)

      if (Ebot.le.0.0) then
         qigrow=energ(Tbot,SALNEW)
         dhib=Ebot/qigrow
      else
!     melt at bottom, melting temperature = Tbot
         qigrow=0.              ! on purpose
         if ((Enet+Ebot).ge.0.0) then
            dhib=-(hi+dtop)
            dhs=-(hs+subs)
            Fx=condb-Enet/dtau
            return
         endif
         call botmelt (qi, qs, delti, delts, ni,
     $                 dhib, dhs, ebot, alarm2)
      endif

! print out helpful variables if something went wrong
      if (alarm2) then
         print*,'Ti ',(Ti(layer),layer=1,ni)
         print*,'Tw ',Tw,condb,Fx,Fneti
         print*,'Enet, Ebot',Enet,Ebot,Fneti*dtau,dtau*(Fmix-condb)
         Enet=0.
         do layer=1,ni
           Enet=Enet+qi(layer)
         enddo
         Enet=Enet*hi/ni+qs*hs
         print*,'orig Enet',Enet
         print*,'delti',delts, (delti(layer),layer=1,ni)
         print*,dhib,dhit,subs,subi,dhs
         alarm=.true.
      endif

! flooding
      zintfc = hi - (rhosno*hs +  rhoice*hi)/rhoocn
      qiflood = 0.
      if ((zintfc .lt. 0.) .and. (hs+dhs.gt.0.)) then
         dhsf   = -rhoice/rhosno*zintfc
         dhsf   = min(dhsf,hs+dhs)
         dhif   = rhosno/rhoice*dhsf
         if (dhif .gt. 0) qiflood = qs*dhsf/dhif
      endif

      call adjust(hi,dhib,dtop,dhif,dhsf,qiflood,qigrow,ni,qi)

      return
      end

      subroutine srfsub( qi   ,qs   ,delti ,delts ,ni    ,
     $                   subi ,subs ,etop  ,enet  ,alarm )

!     compute the sea ice and snow thickness changes from
!     sublimation/condensation
! Etop>0     qi<0  qs<=0 and   Enet<0

      implicit none

      include "cpts.h"
      include "thermo.h"
      include "cembm.h"

! input variables
      real qi (1:nmax), qs    ! energy of melt per vol  (J/m**3)
      real delti(NMAX), delts ! thickness of ice/snow layer (m)
      integer ni              ! number of layers

! output variables
      real subi, subs ! subl/cond. amount for ice/snow (m)

! input/output variables
      real etop ! energy avail to sub/cond ice/snow   (J/m**2)
      real enet ! energy needed to melt all ice/snow  (J/m**2)

! local variables
      real u     ! energy of melt + vapour for ice/snow (J/m**3)
      real subil ! subl from this layer
      integer layer
      logical alarm

      if ( delts .gt. 0. ) then
!       convert etop into snow thickness
!       positive if cond/negative if subl
        u = qs - rvlsno
        subs = etop/u
        if ( (delts + subs) .ge. 0. ) then
!         either all condensation becomes snow
!         or sublimate some snow but no ice
          etop = 0.0
          enet = enet + subs * qs
          return
        else
!         sublimate all of the snow and some ice too
          subs = -delts
          etop = etop + delts * u
          enet = enet + subs  * qs
        endif
      endif

      do layer = 1,ni
!       convert etop into ice thickness
!       positive if cond/negative if subl
        u = qi(layer) - rvlice
        subil = etop/u
        if ( (delti(layer)+subil) .ge. 0. ) then
!         either all condensation becomes ice
!         or sublimate some ice from this layer
          subi = subi + subil
          etop = 0.0
          enet = enet + subil * qi(layer)
          return
        else
!         sublimate all of the layer
          subi = subi - delti(layer)
          etop = etop + delti(layer) * u
          enet = enet - delti(layer) * qi(layer)
        endif
      enddo

! model should never get here, probably a only minor error if it does
      print*, 'ERROR in srfsub',Etop
      alarm=.true.

      return
      end

      subroutine srfmelt(   qi,   qs, delti, delts, ni,
     $                    dhit,  dhs,  etop, alarm  )

! melt ice/snow from the top srf
! Etop>0     qi<0  and  qs<=0

      implicit none

      include "cpts.h"
      include "thermo.h"
      include "cembm.h"

! input variables
      real qi (1:nmax), qs    ! energy of melt per vol (J/m**3)
      real delti(NMAX), delts ! thickness of ice/snow layer (m)
      integer ni              ! number of layers

! output variables
      real dhit, dhs ! ice/snow thickness change(m)

! input/output variables
      real etop ! energy avail to melt ice and snow (J/m**2)

! local variables
      real dhitl ! melt from this layer (m)
      integer layer
      logical alarm

      if ( delts .gt. 0. ) then
!       convert etop into snow thickness
        dhs  =  etop/qs
        if ( (delts + dhs) .ge. 0. ) then
!         melt only some of the snow
          etop = 0.0
          return
        else
!         melt all of the snow and some ice too
          dhs  = -delts
          etop = etop+qs*delts
        endif
      endif

      do layer = 1,ni
!       convert etop into ice thickness
        dhitl = etop/qi(layer)
        if ( (delti(layer)+dhitl) .ge. 0. ) then
!         melt some ice from this layer
          dhit = dhit + dhitl
          etop = 0.0
          return
        else
!         melt all of the ice in this layer
          dhit = dhit - delti(layer)
          etop = etop + delti(layer) * qi(layer)
        endif
      enddo

! model should never get here, probably a only minor error if it does
      print*, 'ERROR in srfmelt',Etop
      alarm=.true.
      return
      end

      subroutine botmelt(   qi,   qs, delti, delts, ni,
     $                    dhib,  dhs,  ebot, alarm  )

! melt from bottom
! Ebot>0     qi<0  and  qs<=0

      implicit none

      include "cpts.h"
      include "thermo.h"
      include "cembm.h"

! input variables
      real qi (1:nmax), qs    ! energy of melt per vol (J/m**3)
      real delti(NMAX), delts ! thickness of ice/snow layer (m)
      integer ni              ! number of layers

! output variables
      real dhib ! ice thickness change (m)

! input/output variables
      real ebot ! energy avail to melt ice and snow (J/m**2)
      real dhs  ! snow thickness change (m)

! local variables
      real dhibl ! melt from this ice layer (m)
      real dhsl  ! melt from this snow layer (m)
      integer layer
      logical alarm

      do layer = ni,1,-1
         dhibl = ebot/qi(layer)
         if ( (delti(layer)+dhibl) .ge. 0. ) then
            dhib = dhib + dhibl
            ebot = 0.0
            return
         else
            dhib = dhib - delti(layer)
            ebot = ebot + delti(layer) * qi(layer)
         endif
      enddo

!     finally melt snow if necessary
      dhsl = ebot/qs
      if ( (delts + dhsl) .ge. 0. ) then
         dhs = dhs + dhsl
         ebot = 0.0
         return
      endif

! model should never get here, probably a only minor error if it does
      print*,'melted completely through all ice and snow layer',delts
      print*,'some diagnostics',(delti(layer),layer=1,ni),(qi(layer)
     $      ,layer=1,ni)
      print*,'ERROR in botmelt',Ebot
      alarm=.true.

      return
      end

      subroutine adjust(hi,dhib,dhit,dhif,dhsf,qiflood,qigrow,ni,qi_tw)

! Adjusts temperature profile after melting/growing.

! E_tw = E at start
! E is the energy density after updating Ti from
! the heat equation without regard dhit and dhib

! h is the thickness from previous time step
! h_tw is the NEW thickness h(i,j)
! dhib is negative if there is melt at the bottom
! dhit is negative if there is melt at the top

! generally _tw is a suffix to label the new layer spacing variables

      implicit none

      include "cpts.h"
      include "thermo.h"
      include "cembm.h"

! input variables
      real hi   ! initial ice thickness (m)
      real dhib ! ice bot, dhib<0 if melt (m)
      real dhit ! ice top, dhit<=0 (m)
      real dhif ! ice top from flooding, dhif>0 (m)
      real dhsf ! snow from flooding, dhsf>0 (m)

      real qiflood ! qi for flooded ice (J/m**3)
      real qigrow  ! qi for ice growing on bot (J/m**3)
      integer ni

! input and output variables
      real qi_tw(nmax) ! qi for adjusted layers (J/m**3)

! local variables
!     the following pairs are before & after adjusting
      real     h, h_tw      ! ice thickness (m)
      real     z(nmax+1)    ! vertical layer position (m)
      real     z_tw(nmax+1) ! vertical layer position (m)
      real     D, D_tw      ! layer thickness (m)
      integer  k, k_tw      ! index of layer

      real qi(nmax)         ! qi of initial layers (J/m**3)
      real fract(nmax,nmax) ! fract of layer k that makes up layer k_tw
      integer layer

! first check if there is any ice left after top/bottom melt/growth
      h = hi
      h_tw = hi+dhib+dhit

      if ( h_tw .le. 0. ) then
         do layer = 1,ni
            qi_tw(layer) = 0.
         enddo
         dhsf = 0.
         dhif = 0.
         return
      endif

! there is ice left so allow snow-ice conversion from flooding
      h_tw = h_tw+dhif

      if ( (abs(dhib) .gt. tiny) .or. (abs(dhit) .gt. tiny)
     $     .or. (dhif.gt.0)) then
         do layer = 1,ni
            qi(layer) = qi_tw(layer)
         enddo

! layer thickness
         D = h/ni
         D_tw = h_tw/ni

! z is positive down and zero is relative to the top
! of the ice from the old time step
         z(1) = -dhit          ! necessary
         z_tw(1) = -dhit-dhif
         do layer = 2,ni
           z(layer) = D*(layer-1)
           z_tw(layer) = z_tw(1)+D_tw*(layer-1)
         enddo
         z(ni+1) = h + min(dhib,0.)
         z_tw(ni+1) = z_tw(1)+h_tw

         do k_tw = 1,ni
           qi_tw(k_tw) = 0.
           do k = 1,ni
             fract(k,k_tw) = ( min( z_tw(k_tw+1), z(k+1)) -
     $                         max( z_tw(k_tw)  , z(k)  )   )
             if (fract(k,k_tw).gt.0.)
     $            qi_tw(k_tw) = qi_tw(k_tw)+fract(k,k_tw)*qi(k)
           enddo
         enddo

         if (dhif.gt.D_tw) then
           do k = 1,ni
             qi_tw(k) = qi_tw(k) + qiflood*max(0.,
     $                   min(D_tw,dhif-D_tw*(k-1)))
           enddo
         else
           qi_tw(1) = qi_tw(1) + dhif*qiflood
         endif

         if (dhib.gt.D_tw) then
           z(ni+2) = h + dhib
           do k = 1,ni
             qi_tw(k) = qi_tw(k) +  qigrow*max(0.,
     $                   (min(z(ni+2),Z_tw(k+1))-max(h,Z_tw(k))))
           enddo
         else
           qi_tw(ni) = qi_tw(ni) + max(dhib,0.)*qigrow
         endif

         do k_tw = 1,ni
           qi_tw(k_tw) = qi_tw(k_tw)/D_tw
         enddo

      endif

      return
      end

      subroutine init_ice_cpts (is, ie, js, je)

! if the cpts ice option is used without the embm, figure out where to
! call this routine

      implicit none

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
      include "cembm.h"
      include "csbc.h"
      include "cpts.h"
      include "ice.h"

      integer i, ie, is, j, je, js, nc, layers

! use double/single Euler steps
      dtau=dts
      idx=1
      if (mod(nats,2).eq.0) idx=2

      if (nats.eq.1) then   ! mixing step, mode 2 is kept
        do nc=1,ncat
          do j=js,je
            do i=is,ie
              A(i,j,1,nc)=A(i,j,2,nc)
              heff(i,j,1,nc)=heff(i,j,2,nc)
              hseff(i,j,1,nc)=hseff(i,j,2,nc)
              Ts(i,j,1,nc)=Ts(i,j,2,nc)
            enddo
          enddo
        enddo
        do layers=1,ntilay
          do j=js,je
            do i=is,ie
              E(i,j,1,layers)=E(i,j,2,layers)
            enddo
          enddo
        enddo
      endif

      return
      end

      subroutine adv_ridge_cpts (is, ie, js, je)

      implicit none

# if defined O_ice_evp
      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
      include "cembm.h"
      include "cpts.h"
      include "ice.h"

      integer i, ie, is, j, je, js, layer, nc

      real A0(is:ie,js:je), ah, A_Ts(is:ie,js:je,ncat)

      do nc=1,ncat
         call embmbc (heff(1,1,idx,nc))
         call embmbc (A(1,1,idx,nc))
         call embmbc (Ts(1,1,idx,nc))
         call embmbc (hseff(1,1,idx,nc))
      enddo
      do layer=1,ntilay
         call embmbc (E(1,1,idx,layer))
      enddo

!     compute open water area so it can be advected
      do j=js,je
        do i=is,ie
          A0(i,j)=1.0
          do nc=1,ncat
            A0(i,j)=A0(i,j)-A(i,j,idx,nc)
          enddo
          A0(i,j)=max(0.,A0(i,j))
        enddo
      enddo

!     proper surface temperature advection is done by advecting
!     the area-temperature product
      do nc=1,ncat
        do j=js,je
          do i=is,ie
            A_Ts(i,j,nc) = min(A(i,j,idx,nc)*Ts(i,j,idx,nc),0.)
          enddo
        enddo
      enddo

! At this time mechanical redistribution for this model assumes
! you are using upstream advection.
      call advupb(uice, vice ,A0, is, ie, js, je)
      do nc=1,ncat
        call advupb(uice, vice, A(1,1,idx,nc), is, ie, js, je)
      enddo
      do nc=1,ncat
        call advupb(uice, vice, A_Ts(1,1,nc), is, ie, js, je)
      enddo
      do nc=1,ncat
        call advupb(uice, vice, heff(1,1,idx,nc), is, ie, js, je)
      enddo
      do nc=1,ncat
        call advupb(uice, vice, hseff(1,1,idx,nc), is, ie, js, je)
      enddo
      do layer=1,ntilay
        call advupb(uice, vice, E(1,1,idx,layer), is, ie, js, je)
      enddo

      do nc=1,ncat
        do j=js,je
          do i=is,ie
            Ts(i,j,idx,nc)=frzpt(i,j)
            if (A(i,j,idx,nc).gt.tiny) then
              Ts(i,j,idx,nc)=A_Ts(i,j,nc)/A(i,j,idx,nc)
              Ts(i,j,idx,nc)=min(Ts(i,j,idx,nc),0.)
            endif
          enddo
        enddo
      enddo

      call mechred(A0, is, ie, js, je)

# endif

      return
      end

! Energy-conserving sea ice model
! Routines to do basic energy/temperature conversion

! See Bitz, C.M., and W.H. Lipscomb, 1999:
! A new energy-conserving sea ice model for climate study,
! J. Geophys. Res., in press.

! The author grants permission to the public to copy and use this
! software without charge, provided that this Notice and any statement
! of authorship are reproduced on all copies.

      function energ(Tmp,sal)
!     compute the energy of melting for non-new ice
!     relative to melting (negative quantity)
!     assuming Tmp is in Celsius
      implicit none
      include "cpts.h"
      include "thermo.h"

      real Tmp,energ,sal

      energ=-rflice
     $        -rcpice*(-alpha*sal-Tmp)-gamma*sal/Tmp

      return
      end

      subroutine getTmp(Tmp,E,sal1d,ni)

!     compute the midpoint temperature from the
!     energy of melting for each layer

      implicit none

      include "cpts.h"
      include "thermo.h"

      real Tmp(nmax),E(nmax),sal1d(nmax),quad
      integer ni,i

      real(kind=8) q, B, C

      do i=1,ni
        q=E(i)+rflice-rcpice*alpha*sal1d(i)
        B=-q/rcpice
        C=-gamma*sal1d(i)/rcpice
        Tmp(i)= quad(B,C)
      enddo

      return
      end

      function quad(B,C)
      implicit none
      include "cpts.h"
      include "thermo.h"

      real quad,root,T1,T2
      real(kind=8) B, C, B_2

      B_2=B/2.0
      root=sqrt(B_2*B_2-C)
!      T1=-B_2+root
      T2=-B_2-root
      quad=T2
!      print*,'quad T1=',T1,T2
!      if (T1.le.(TMELT-0.1)) print*,'quad T1=',T1,T2
      return
      end

      subroutine grownew(fsh_new,Io_new,flo,focean,Fbot,
     $            hnew,tnew,tair,qair,ug,dtp,fice,frzpt,evp,
     $            dswr,ultnt,usens,ulwr)

!     Uses a Newton_Raphson (non-linear)
!     method to compute new ice surface temp
!     and thickness.
!     This algorithm tends to fail when the correct solution is
!     to grow a very thin layer of ice. In this case I just set
!     the new ice surface temperature to 0.5*(Tair+frzpt) and
!     grow accordingly. It still conserves energy.

      implicit none

      include "cpts.h"
      include "thermo.h"
      include "cembm.h"
      include "calendar.h"

      real  coltnt,cosens
      real  fsh_new,flo,Io_new
      real  Fbot,focean,hnew,tnew,frzpt
      real  dswr,ultnt,usens,ulwr
      real  tair,qair,evp,qnew,rday
      real  ug,dt
      real  dtp
      real  b,fice,e1,e2,e3,e4,e5,terror,herror
      real  u
      real  J11,J12,J21,J22,F1,F2,detJ
      real  rand,energ,kthin
      real  ult_new,use_new,ulw_new,evp_new
      integer iter

      kthin=kappai+beta*salnew/(frzpt-5.)
      Tnew=0.5*(frzpt+Tair) ! initial guess
      Tnew=min(frzpt,Tnew)

# if defined O_ice_cpts_simple_growth
       fice=focean
       u=0.5*(frzpt+Tnew)
       u=energ(u,salnew)
       hnew=(focean+Fbot)*dtp/u
# else
      dt=Tnew-Tair
      cosens=0.94*rcpatm*dalt_i*ug

      ult_new=ultnt
      e1=fsh_new-ult_new+flo+cosens*tair

      e2=-cosens*Tnew-ESICE*(Tnew+TFFRESH)**4
      fice=e1+e2

      if ((fice+focean+2*Fbot).lt.0.0) then
         rday=daylen/dtau
         hnew=centi*0.0133/rday*(frzpt-Tnew)

         iter=1

 100     continue

            e2=-cosens*Tnew-esice*Tnew**4
            e3=-cosens-4.*ESICE*(Tnew+TFFRESH)**3

            if (Tnew.gt.frzpt) then
               u=frzpt
            else
               u=0.5*(frzpt+Tnew)
            endif
            e4=-0.5*(rcpice-gamma*salnew/(u**2))
            u=-energ(u,salnew)

            J11=hnew*e3-kthin
            J12=e1+e2
            J21=0.5*e3+hnew/dtp*e4
            J22=u/dtp
            F1=(J12-Io_new)*hnew+kthin*(frzpt-Tnew)
            F2=0.5*(focean+J12)+Fbot+hnew*J22
            detJ=J11*J22-J12*J21

            terror=(J22*F1-J12*F2)/detJ
            herror=(-J21*F1+J11*F2)/detJ
            Tnew= Tnew-terror
            hnew=hnew-herror
            iter=iter+1
            if ((iter.gt.10).or.(hnew.lt.0).or.(Tnew.gt.0.)) then
!     when cannot converge, assume a temperature and continue
               Tnew=0.5*(frzpt+Tair)
               Tnew=min(frzpt,Tnew)
               go to 200
            endif
            if (abs(terror).gt.0.005) go to 100

 200        continue

            u=0.5*(frzpt+Tnew)
            u=energ(u,salnew)

            ulw_new=esice*(Tnew+TFFRESH)**4-Flo
            use_new=cosens*(Tnew-tair)

            fice=fsh_new-ulw_new-use_new-ult_new
            hnew=(0.5*(fice+focean)+Fbot)*dtp/u

            dswr = 0.5*(dswr+fsh_new)
            ulwr = 0.5*(ulwr +ulw_new)
            usens = 0.5*(usens+use_new)

      else
!     it is not clear what to do with the negative
!     ocean energy balance if in the presence of
!     open water the energy balance is negative and
!     would allow ice to grow, yet once it grows a little
!     ice the energy balance swings to positive

!     I am going to set h=0 in this case and balance the fluxes
!     by messing with the sensible heat

         usens = dswr-ultnt-ulwr+Fbot
         fice=-focean-2*Fbot
         Tnew=frzpt
         hnew=0
      endif

# endif
      return
      end

      subroutine mechred (A0, is, ie, js, je)

!  it is perfectly fine to call this routine when NCAT=1
!  in fact, it will provide a source of open water from shear
!  deformation. See Stern et al, 1995 for information about
!  how and why this is done

      implicit none

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
      include "evp.h"
      include "atm.h"
      include "cpts.h"
      include "ice.h"
      include "grdvar.h"

!   input variables
      real A0(is:ie,js:je)

!    local variables

      integer i, ie, iem1, is, isp1, j, je, jem1, js, jsp1, n, nc
      integer layer, l

      real En(NMAX,ncat)
      real Esn(ncat)
      real heffn(ncat)
      real hseffn(ncat)
      real An(ncat)
      real Tsn(ncat)
      real closing       ! lead closing rate
      real M(ncat,ncat)  ! area fraction of ice from i into j
      real HN(ncat,ncat) ! H*volume fraction of ice from i into j
      real rpress

      isp1 = is + 1
      iem1 = ie - 1
      jsp1 = js + 1
      jem1 = je - 1

# if defined O_ice_cpts_roth_press
      do j=js,je
        do i=isp1,ie
          strength(i,j,idx)=0.
        enddo
      enddo
# endif

      do n=1,nseg
        do j=jsi(n),jei(n)+1
          do i=isp1,iem1
            if (tmsk(i,j) .ge. 0.5) then

!    ice velocity divergence and lead closing rate, following FH95
              closing=0.25*(del(i,j)-abs(eI(i,j)))-min(eI(i,j),0.)
              closing=max(closing,-eI(i,j))

              l=1
              do nc=1,ncat
                do layer=1,nilay(nc)
                   En(layer,nc)=E(i,j,idx,l)
                   l=l+1
                enddo
                Tsn(nc)=Ts(i,j,idx,nc)
                An(nc)=A(i,j,idx,nc)
                heffn(nc)=heff(i,j,idx,nc)
                hseffn(nc)=hseff(i,j,idx,nc)
              enddo

              call ridge(heffn,hseffn,A0(i,j),An,En,Esn,Tsn,
     $            closing,eI(i,j),M,HN,frzpt(i,j))

# if defined O_ice_cpts_roth_press
              strength(i,j,idx)=rpress(heffn,A0(i,j),An,HN,M)
# endif
              l=1
              do nc=1,ncat
                do layer=1,nilay(nc)
                   E(i,j,idx,l)=En(layer,nc)
                   l=l+1
                 enddo
                 Ts(i,j,idx,nc)=Tsn(nc)
                 A(i,j,idx,nc)=An(nc)
                 heff(i,j,idx,nc)=heffn(nc)
                 hseff(i,j,idx,nc)=hseffn(nc)
              enddo

            endif
          enddo
        enddo
      enddo

      do nc=1,ncat
         call embmbc (heff(1,1,idx,nc))
         call embmbc (A(1,1,idx,nc))
         call embmbc (Ts(1,1,idx,nc))
         call embmbc (hseff(1,1,idx,nc))
      enddo
      do layer=1,ntilay
         call embmbc (E(1,1,idx,layer))
      enddo
      return
      end

      subroutine ridge(heffn,hseffn,A0,An,En,Esn,Tsn,closed,divu,
     $     M,HN,frzpt)

!     note the time available for ridging between the cats is
!     dtau, however have to subcycle if a category runs
!     out of ice then dtn is less than dtau

!     M(nc,k) is the fraction of ice that participates in ridging

      implicit none

      include "cpts.h"
      include "thermo.h"

      real heffn(ncat)
      real hseffn(ncat)
      real An(ncat)
      real Tsn(ncat),dtf(ncat),Atsn(ncat)
      real En(nmax,ncat)
      real Esn(ncat)
      real closed        ! the closing as defined by Stern et al
      real divu, A0
      real frzpt
      real strg          ! ice strength

      real Asum
      real W(0:ncat)     ! ridging for each category
      real Wa(0:ncat)    ! the area that participates
      real M(ncat,ncat)  ! area fraction of ice from i into j
      real N(ncat,ncat)  ! volume fraction of ice from i into j
      real HN(ncat,ncat) ! dummy
      real Hmean(ncat)   ! Hi+sqrt(cK*Hi)

      real dtn,dtnn,dtp
      integer k,layer,index,isub,nc,layk,lay,i,j

      real Qi(nmax,ncat) ! energy of melting ice per unit volume
      real Qs(ncat)      ! energy of melting of snow per unit area
      real delt(ncat),delts(ncat),deltEs(ncat)
      real deltE(nmax,ncat)

      real Psi, closdtn

      real dsum,ssum,esum,fsum
      real Hi(ncat),Hs(ncat), eta

      dtn=dtau

      dtp=0

      A0=A0
      Asum=A0
      do nc=1,ncat
         Asum = Asum + An(nc)
      enddo

! unfortunately the second order accurate advection
! schemes do not have Asum=1-divu*dt after advection
! owing to diffusion.
! Adjust opening rate here to compensate.
      divu=(1-Asum)/dtn

      closed=-min(divu,0.) ! from Bill Hibler suggestions
      closed=max(closed,-divu)

      A0=A0+dtn*(divu+closed)
      Asum=Asum+dtn*(divu+closed)

      if (Asum.le.1.0) then
        A0=A0+1.-Asum
# if defined O_ice_cpts_roth_press
        call ridge_matrices(heffn,An,Hi,Hmean,M,N,HN)
# endif
        return
      endif

 50   continue

      call ridge_matrices(heffn,An,Hi,Hmean,M,N,HN)
      call ridging_mode(A0,An,W,Wa,M)

# if defined O_ice_cpts3 || defined O_ice_cpts5 || defined O_ice_cpts10
      do nc=1,ncat
         Hs(nc)=0.0
         Qs(nc)=0.0
         do layer=1,nilay(nc)
            Qi(layer,nc)=0.0
         enddo

         if ((An(nc).gt.TINY) .and. (heffn(nc).gt.0.))  then
            Hs(nc)=hseffn(nc)/An(nc)
            Qs(nc)=Esn(nc)/An(nc)
            do layer=1,nilay(nc)
               Qi(layer,nc)=En(layer,nc)*nilay(nc)/heffn(nc)
            enddo
         endif
      enddo

      dsum=0.
      ssum=0.
      esum=0.
      do nc=1,ncat-1
         delt(nc)  =-Hi(nc)*Wa(nc)
         delts(nc) =-Hs(nc)*Wa(nc)
         deltEs(nc)=-Qs(nc)*Wa(nc)
         do k=1,nc
            delt(nc)  =delt(nc)+Wa(k)*N(k,nc)
            delts(nc) =delts(nc)+Wa(k)*M(k,nc)*Hs(k)*Hmean(k)/hi(k)
            deltEs(nc)=deltEs(nc)+Wa(k)*M(k,nc)*Qs(k)*Hmean(k)/hi(k)
         enddo
         dsum=dsum+delt(nc)
         ssum=ssum+delts(nc)
         esum=esum+deltEs(nc)
      enddo

      delt(ncat)  =-dsum
      delts(ncat) =-ssum
      deltEs(ncat)=-esum

      fsum=0.
      do nc=1,ncat-1
         do layer=1,nilay(nc)
            deltE(layer,nc)=-Wa(nc)*Hi(nc)*Qi(layer,nc) /nilay(nc)
            do k=1,nc
               layk=(layer+ncrel(k,nc)-1)/ncrel(k,nc)
               deltE(layer,nc)=deltE(layer,nc)+Qi(layk,k)*Wa(k)*N(k,nc)
     $                        /nilay(nc)
            enddo
            fsum=fsum+deltE(layer,nc)
         enddo
      enddo
      do layer=1,nilay(ncat)
         deltE(layer,ncat)=0.
         do k=1,ncat-1
            layk=(layer+ncrel(k,ncat)-1)/ncrel(k,ncat)
            deltE(layer,ncat)=deltE(layer,ncat)+Qi(layk,k)*Wa(k)
     $                       *N(k,ncat)/nilay(ncat)
         enddo
         fsum=fsum+deltE(layer,ncat)
      enddo

!     geometrical factor for treating ice/snow surface temp when ridging
      do nc = 1,ncat
         dtf(nc) = - tsn(nc)*Wa(nc)
         do k=1,nc
           dtf(nc) = dtf(nc)+ tsn(k)*Wa(k)*M(k,nc)
         enddo
      enddo

# endif

!     make sure that only ridge as much ice as is available, otherwise
!     reduce timestep
      Psi=closed*Wa(0)
      if ((Psi*dtn.gt.A0) .and. (Psi.gt.0.)) then
         dtnn=A0/(Psi)
         dtn=min(dtn,dtnn)
      endif
      do nc=1,ncat
         Psi=closed*Wa(nc)
         if ((Psi*dtn.gt.An(nc)) .and. (Psi.gt.0.)) then
            dtnn=An(nc)/(Psi)
            dtn=min(dtn,dtnn)
         endif
      enddo
      dtp=dtp+dtn

      closdtn=closed*dtn

      A0=A0+closdtn*W(0)
      Asum=A0
      do nc=1,ncat
         Atsn(nc) = An(nc)*tsn(nc) + closdtn*dtf(nc)
         An(nc)=An(nc)+closdtn*W(nc)
         Asum=Asum+An(nc)
# if defined O_ice_cpts3 || defined O_ice_cpts5 || defined O_ice_cpts10
         heffn(nc)=heffn(nc)+closdtn*delt(nc)
         hseffn(nc)=hseffn(nc)+closdtn*delts(nc)
         Esn(nc)=Esn(nc)+closdtn*deltEs(nc)
         do layer=1,nilay(nc)
            En(layer,nc)=En(layer,nc)+closdtn*deltE(layer,nc)
         enddo
# endif
      enddo
      do nc = 1,ncat
        if (An(nc).gt.tiny) then
          tsn(nc) = Atsn(nc)/An(nc)
        else
          tsn(nc) = frzpt
        endif
      enddo

      dtn=dtau-dtp

      if (dtn.gt.TINY) go to 50

      return
      end

      subroutine ridging_mode(A0,An,W,Wa,M)

!     compute W and Wa
!     M(nc,k) is the fraction of ice that participates in ridging
!     from category nc that goes into category k

      implicit none

      include "cpts.h"

      real An(ncat),A0

      real W(0:ncat)       ! ridging for each category
      real Wa(0:ncat)      ! the area that participates
      real M(ncat,ncat)    ! area fraction of ice from i into j

      real V(1:ncat+1)     ! the area that participates
      real Asum,eta,Wsum
      integer k,layer,index,isub,nc,layk,lay,i,j

      Asum=0.
      Wa(0)=min(gstar-Asum,A0)
      Asum=Asum+Wa(0)
!      print*,0,Wa(0),Asum,(1.0 + 0.5*Wa(0)/gstar - Asum/gstar)
      Wa(0)=Wa(0)*(1.0 + 0.5*Wa(0)/gstar - Asum/gstar)
      Wa(0)=max(Wa(0),0.)
      do nc=1,ncat
         if (An(nc).gt.0.) then
            Wa(nc)=min(gstar-Asum,An(nc))
            Asum=Asum+Wa(nc)
            Wa(nc)=Wa(nc)*(1.0 + 0.5*Wa(nc)/gstar - Asum/gstar)
            Wa(nc)=max(Wa(nc),0.)
         else
            Wa(nc)=0.
         endif
      enddo

      Wsum=-Wa(0)
      do nc=1,ncat
         W(nc)=-Wa(nc)
         do k=1,nc
            W(nc)=W(nc)+Wa(k)*M(k,nc)
         enddo
         Wsum=Wsum+W(nc)
      enddo

      eta=-Wsum
      do nc=1,ncat
         W(nc)=W(nc)/eta
         Wa(nc)=Wa(nc)/eta
      enddo
      Wa(0)=Wa(0)/eta
      W(0)=-Wa(0)

      return
      end

      subroutine ridge_matrices(heffn,An,Hi,Hmean,M,N,HN)

!     compute M,N,HN for ridge and press for true thicknesses

      implicit none

      include "cpts.h"

      real heffn(ncat)
      real An(ncat)

      real Hi(ncat)
      real M(ncat,ncat)    ! area fraction of ice from i into j
!     N is used in ridge and HN is used in rpress
      real N(ncat,ncat)    ! volume fraction of ice from i into j
      real HN(ncat,ncat)   ! H*volume fraction of ice from i into j

      real rcKH            ! sqrt(cK*Hi)
      real Hmean(ncat)     ! Hi+sqrt(cK*Hi)
      real bet             ! 2*(K-Hi)
      integer n1,n2        ! see above
      integer k,nc
      logical row_flag(ncat)

      do nc=1,ncat
         if ((An(nc).gt.TINY) .and. (heffn(nc).gt.0.))  then
            Hi(nc)=heffn(nc)/An(nc)
            row_flag(nc)=.true.
            do k=1,ncat
              M(nc,k)=0.
              N(nc,k)=0.
              HN(nc,k)=0.
            enddo
         else
            Hi(nc)=hstar(nc-1)
            row_flag(nc)=.false.
            do k=1,ncat
              M(nc,k)=M_def(nc,k)
              N(nc,k)=N_def(nc,k)
              HN(nc,k)=HN_def(nc,k)
            enddo
         endif
      enddo

      call comp_matrices(row_flag,Hi,Hmean,M,N,HN)

      return
      end

      subroutine comp_matrices(row_flag,Hi,Hmean,M,N,HN)

!     compute M,N,HN for ridge and rpress for true thicknesses
!     N is used in ridge and HN is used in rpress

      implicit none

      include "cpts.h"

      logical row_flag(ncat)
      real Hi(ncat)
      real M(ncat,ncat)    ! area fraction of ice from i into j
      real N(ncat,ncat)    ! volume fraction of ice from i into j
      real HN(ncat,ncat)   ! H*volume fraction of ice from i into j

      real rcKH            ! sqrt(cK*Hi)
      real Hmean(ncat)     ! Hi+sqrt(cK*Hi)
      real bet             ! 2*(K-Hi)
      integer n1,n2        ! see above
      integer k,layer,index,isub,nc,layk,lay,i,j

      do nc=1,ncat
         rcKH=sqrt(cK*Hi(nc))
         Hmean(nc)=Hi(nc)+rcKH
         bet=2*(cK-Hi(nc))

         if (row_flag(nc)) then
! ice from category nc is ridging into a distribution from n1 to n2
            n1=nc  ! first guess
            n2=n1  ! first guess
            do k=nc+1,ncat
               if (2*Hi(nc).ge.hstar(k-1)) n1=k
            enddo
            if (n1.eq.ncat) then
!              case when all ice ridges into thickest cat
               M(nc,n1)=Hi(nc)/Hmean(nc)
               N(nc,n1)=Hi(nc)
               HN(nc,n1)=Hi(nc)*Hmean(nc)
            elseif (2*rcKH.le.hstar(n1)) then
!              case when all ice ridges into just one category
               M(nc,n1)=Hi(nc)/Hmean(nc)
               N(nc,n1)=Hi(nc)
               HN(nc,n1)=Hi(nc)*Hmean(nc)
            else
               n2=ncat
               do k=ncat-1,n1+1,-1
                  if (2*rcKH.lt.hstar(k)) n2=k
               enddo
               M(nc,n1)=(hstar(n1)-2*Hi(nc))/bet
               N(nc,n1)=M(nc,n1)*(Hi(nc)+hstar(n1)/2.)
               HN(nc,n1)=N(nc,n1)*(Hi(nc)+hstar(n1)/2.)

               M(nc,n2)=(2.*rcKH-hstar(n2-1))/bet
               N(nc,n2)=M(nc,n2)*(hstar(n2-1)/2.+rcKH)
               HN(nc,n2)=N(nc,n2)*(hstar(n2-1)/2.+rcKH)

               do k=n1+1,n2-1
                  M(nc,k)=(hstar(k)-hstar(k-1))/bet
                  N(nc,k)=M(nc,k)*(hstar(k)+hstar(k-1))/2.
                  HN(nc,k)=N(nc,k)*(hstar(k)+hstar(k-1))/2.
               enddo
            endif
         endif
      enddo

      return
      end

      function rpress(heffn,A0,An,HN,M)
!     Compute the ice strength based on Rothrock, 1975 and also see FH95
!     does not make sense to do this unless the ice that participates in
!     ridging is well resolved -- must have about 5 categories, 10 would
!     be better

!     M(nc,k) is the fraction of ice that participates in ridging
!     from category nc that goes into category k
!     n1 is the smallest k with nonzero fraction
!     n2 is the largest  k with nonzero fraction
!     the ice that participates has thickness Hi(nc)
!     and it ridges up to linear distribution between
!     2*Hi(nc) and 2*sqrt(cK*Hi(nc))

      implicit none
      include "cpts.h"

      real rpress
      real heffn(ncat)
      real An(ncat)

      real W(0:ncat)         ! ridging for each category
      real Wa(0:ncat)        ! the area that participates
      real M(ncat,ncat)    ! area fraction of ice from i into j
      real HN(ncat,ncat)   ! H*volume fraction of ice from i into j

      real A0

      integer n1,n2        ! see above
      integer k,nc,i,j

      real Hi(ncat),Hs(ncat)

      real peterm, frterm

!     ideally cpe should be a function of density (a very sensitive one)
!     see FH95, I set it equal to a number so you are more likely to get
!     a reasonable ice strength.
!     It has units (kg/m**3)*(m/s**2) = (kg/m**2/s)

      rpress=0.

# if defined O_ice_cpts_roth_press
      if (A0.ge.gstar) return

      do nc=1,ncat
         if ((An(nc).gt.TINY) .and. (heffn(nc).gt.0.))  then
            Hi(nc)=heffn(nc)/An(nc)
         else
            Hi(nc)=hstar(nc-1)
         endif
      enddo

!     strictly speaking, should recompute the M and HN matrices here,
!     but they would be only very slightly different from what they were
!     before ridging and it does not affect the mass conservation
      call ridging_mode(A0,An,W,Wa,M)

      peterm=0.
      do nc=1,ncat
         peterm=peterm-Wa(nc)*Hi(nc)*Hi(nc)
         do k=1,nc
            peterm=peterm+Wa(k)*HN(k,nc)
         enddo
      enddo
      peterm=max(0.,cpe*peterm)

      rpress=peterm*17.

# endif
      return
      end

      subroutine movedown(frzpt,heff1,hseff1,A1,hsnow1,hice1,h1min,Es1,
     $     ts1,nn1,heff2,hseff2,A2,hsnow2,hice2,h2min,Es2,ts2,nn2,E1,E2)

!     hice and hsnow are outputs (not needed for input)

      implicit none

      include "cpts.h"
      real heff1,hseff1,A1,hsnow1,hice1,h1min,Es1,E1(NMAX),ts1
      real heff2,hseff2,A2,hsnow2,hice2,h2min,Es2,E2(NMAX),ts2
      real frzpt
      integer nn1,nn2,M,layer,index,isub

      real x(NMAX)
      real Atemp

      Atemp=A1+A2
      if (Atemp.le.0.0) return

      heff1=heff1+heff2
      hseff1=hseff1+hseff2
      ts1=ts1*A1+ts2*A2

      A1=Atemp

      hsnow1=hseff1/A1
      hice1=heff1/A1
      ts1=ts1/A1

      Es1=Es1+Es2

      M=nn2/nn1
!     combine the N2 layers of energy in category 2
!     into N1 portions
      do 1000 layer=1,nn1
         x(layer)=0.
         do 1010 isub=1,M
            index=(layer-1)*M+isub
            x(layer)=x(layer)+E2(index)
 1010    continue
 1000 continue

!     sum the energies
      do layer=1,nn1
         E1(layer)=E1(layer)+x(layer)
      enddo

      call zerocat(heff2,hseff2,A2,hsnow2,hice2,h2min,
     $     Es2,E2,ts2,frzpt,nn2)

      return
      end

      subroutine moveup(frzpt,heff1,hseff1,A1,hsnow1,hice1,h1min,Es1,
     $   ts1,nn1,heff2,hseff2,A2,hsnow2,hice2,h2min,Es2,ts2,nn2,E1,E2)

!     hice and hsnow are outputs (not needed for input)

      implicit none

      include "cpts.h"
      real heff1,hseff1,A1,hsnow1,hice1,h1min,Es1,E1(NMAX),ts1
      real heff2,hseff2,A2,hsnow2,hice2,h2min,Es2,E2(NMAX),ts2
      real frzpt
      integer nn1,nn2,M,layer,index,isub

      real x(NMAX)
      real Atemp

      Atemp=A1+A2
      if (Atemp.le.0.0) return

      heff2=heff1+heff2
      hseff2=hseff1+hseff2
      ts2=ts1*A1+ts2*A2

      A2=Atemp

      hsnow2=hseff2/A2
      hice2=heff2/A2
      ts2=ts2/A2

      Es2=Es1+Es2

      M=nn2/nn1
!     move energy of nn1 layers from 1 into the nn2 of 2
      do index=1,nn1
         do isub =1,M
            layer=M*(index-1)+isub
            x(layer)=E1(index)/M
         enddo
      enddo

!     sum the energies
      do layer=1,nn2
         E2(layer)=E2(layer)+x(layer)
      enddo

      call zerocat(heff1,hseff1,A1,hsnow1,hice1,h1min,Es1,E1,ts1,frzpt,
     $             nn1)

      return
      end

      subroutine zerocat(heff,hseff,A,hsnow,hice,hmin,Es,E,ts,frzpt,nng)

      implicit none

      include "cpts.h"
      real heff,hseff,A,hsnow,hice,hmin,Es,E(NMAX)
      real ts, frzpt
      integer nng,layer

      heff=0.0
      hseff=0.0
      A=0.0
      hsnow=0.0
      hice=hmin
      Es=0.0
      do layer=1,nng
         E(layer)=0.0
      enddo
      ts=frzpt
      return
      end

      subroutine thermo (is, ie, js, je)

!     heat budget

!     compute fluxes, temp and thickness  over open water and ice

!     flux is positive toward the atm/ice or atm/ocean surface even if
!     the flux is below the surface

      implicit none

      include "size.h"
      include "param.h"
      include "pconst.h"
      include "stdunits.h"
      include "csbc.h"
      include "cembm.h"
      include "atm.h"
      include "cpts.h"
      include "ice.h"
      include "thermo.h"
      include "coord.h"
      include "grdvar.h"

      real  tocn           ! surface ocean temperature in Celsius
      real  tair           ! surface air temperature in celcius
      real  qair           ! surface air humidity
      real  ug             ! geostrophic wind speed

      real  evp(0:ncat)    ! evapourative water flux for each category
      real  wff            ! ice-ocean water flux for each category

      real  dswr(0:ncat)   ! downward shortwave for each category
      real  ultnt(0:ncat)  ! upward latent heat for each category
      real  usens(0:ncat)  ! upward sensible heat for each category
      real  ulwr(0:ncat)   ! upward longwave for each category

      real  Io             ! solar transmitted through top surface
      real  Ib             ! solar transmitted through bottom surface
      real  absorb         ! sum of Ib*A for each category

      real  Focean         ! net flux into ocean from atmosphere
      real  Flwatm         ! downward longwave from atmosphere (not net)
      real  Fbot           ! flux from ocean due to thermal relaxation
      real  Flat           ! flux from ocean due to lateral melt

      real  Fice           ! net flux into ice from atmosphere
      real  Fneti          ! Fice + conductive flux
      real  Fx             ! adjusted Fbot for when ice melts away
      real  condb          ! conductive flux in ice at bottom surface

      real  A0             ! open water fraction
      real  Fsh_new        ! solar incident on new ice (over open water)
      real  Io_new         ! solar transmitted through top of new ice
      real  hnew           ! thickness of new ice
      real  Tnew           ! surface temperature of new ice

      real  hi(ncat)       ! ice thickness for each category
      real  hs(ncat)       ! snow thickness for each category
      real  En(1:NMAX,ncat)! energy of melt of each layer and category
      real  qi(1:NMAX,ncat)! like Ei but per unit volume
      real  Ti(0:NMAX,ncat)! Temperature of each layer and category (C)
      real  Tbot(ncat)     ! Temperature of ice bottom and category (C)

      real  delhib         ! ice thickness change at bottom
      real  delhit         ! ice thickness change (no sublimation)
      real  delhs          ! snow thickness change (no sublimation)
      real  subi           ! ice thickness change from sublimation
      real  subs           ! snow thickness change from sublimation
      real  delf           ! ice thickness change from flooding
      real  delfs          ! snow thickness change from flooding

      integer i,j,l,layer,nc,nr  ! counters

      real sblmt            ! heat to ice-snow from sublimation (J/m**2)
      real htatm            ! heat to ice-snow from atm (J/m**2)
      real wtatm            ! water to ice-snow from atm (m/s)
      real htocn            ! heat to ice-snow from ocn (J/m**2)
      real wtocn            ! waterto ice-snow from ocn (m/s)

      real einit            ! initial energy in ice/snow (J/m**2)
      real winit            ! initial water equiv. in ice/snow (m)
      real efinl            ! final energy in ice/snow (J/m**2)
      real wfinl            ! final water equiv. in ice/snow (m)
      real ediff            ! energy difference in ice/snow (J/m**2)
      real wdiff            ! water difference  in ice/snow (m)

      logical alarm

! vars for small stuff check and computing Ti, hi, hs and albedo
      real Aleft, tio, Enet

! lateral melt vars
      real Esum   ! sum of energy in ice and snow (J/m**2)
      real mrate  ! lateral melt rate (m/s)
      real arate  ! rate ice concentration is melting laterally (s)
      real DAice  ! change in concentration from lateral melt
      real Vfrac  ! ratio of volume after to before melting
      real Dhlatm, Dslatm, sla, fracsnow

! redistribute the ice categories vars
      real h0prime,A0prime,h1prime,A1prime,E0,Ad0
      real energ,Tmp,epsilon

      real ca, as, aicet, coaloff, C2K

      integer ie, iem1, is, isp1, je, jem1, js, jsp1

!     dummy energy of melting in snow for calls to movedown/up
      real esdum1, esdum2
!-----------------------------------------------------------------------

      sla = zw(1)*secday/dampice/2.389e-8
      isp1 = is + 1
      iem1 = ie - 1
      jsp1 = js + 1
      jem1 = je - 1
      C2K = 273.15

      do j=jsp1,jem1
        do i=isp1,iem1
          if (tmsk(i,j) .ge. 0.5) then

            Io_new=0.0
!     adjust albedo for new ice growth
# if defined O_sulphate_data || defined O_sulphate_data_transient
            fsh_new = solins(i,j)*sbc(i,j,iaca)*pass
     &                *max(c0, sbc(i,j,isca) - sulph(i,j,2))
# else
            fsh_new = solins(i,j)*sbc(i,j,iaca)*pass*sbc(i,j,isca)
# endif
            tocn = sbc(i,j,isst)
            tair = at(i,j,2,isat)
# if defined O_sealev || defined O_sealev_data
     &           - elev_sealev(i,j)*rlapse
# endif
            qair = at(i,j,2,ishum)
            Fbot = -sla*(frzpt(i,j) - tocn)
            ug = sbc(i,j,iws)
            Flwatm = esatm*(tair + C2K)**4
            Focean = dnswr(i,j) - uplwr(i,j) - upsens(i,j) - upltnt(i,j)
            dswr(0) = dnswr(i,j)
            ultnt(0) = upltnt(i,j)
            usens(0) = upsens(i,j)
            ulwr(0) = uplwr(i,j)
            evp(0) = evap(i,j)

!-----------------------------------------------------------------------
!     make simple energy of melt matrix for speed
!-----------------------------------------------------------------------
            do nc=1,ncat
              do layer=1,nilay(nc)
                En(layer,nc)=E(i,j,idx,layer1(nc)+layer-1)
              enddo
            enddo

!-----------------------------------------------------------------------
!     initial energy and water diagnostics
!     could be eliminated to save time
!-----------------------------------------------------------------------
            einit = 0.
            winit = 0.
            do nc=1,ncat
              do layer=1,nilay(nc)
                einit=einit+En(layer,nc)
              enddo
              einit = einit-hseff(i,j,idx,nc)*rflsno
              winit = winit+(heff(i,j,idx,nc)*rhoice
     $                        +hseff(i,j,idx,nc)*rhosno)
            enddo

# if defined O_ice_cpts3 || defined O_ice_cpts5 || defined O_ice_cpts10
!-----------------------------------------------------------------------
!     get rid of small area by moving to other category (if possible)
!     does not fix negative ice area/volumes (that comes later)

!     push small bit of ice to thicker category
!     hi and hs are outputs of moveup (not inputs)
!-----------------------------------------------------------------------
            do nc=1,(ncat-1)
              if ((A(i,j,idx,nc) .le. ASMALL(nc))
     $             .and.(A(i,j,idx,nc) .gt. 0.0) ) then

                A(i,j,idx,nc)=heff(i,j,idx,nc)/hstar(nc)
                call moveup(frzpt(i,j),heff(i,j,idx,nc),
     $               hseff(i,j,idx,nc),
     $               A(i,j,idx,nc),hs(nc),hi(nc),hstar(nc-1),
     $               esdum1,Ts(i,j,idx,nc),nilay(nc),
     $               heff(i,j,idx,nc+1),hseff(i,j,idx,nc+1),
     $               A(i,j,idx,nc+1),hs(nc+1),hi(nc+1),hstar(nc),
     $               esdum2,Ts(i,j,idx,nc+1),nilay(nc+1),
     $               En(1,nc),En(1,nc+1))
              endif
            enddo

!     push small bit of ice to thinner category
            do nc=ncat,2,-1
              if ((A(i,j,idx,nc) .le. ASMALL(nc))
     $             .and.(A(i,j,idx,nc) .gt. 0.0) ) then
!     reshape the thicker ice to thickness of thinner category
!     while making certain not to make area > 1
                Aleft=1.0
                do  nr=1,(nc-1)
                  Aleft=Aleft-A(i,j,idx,nr)
                enddo
                A(i,j,idx,nc)=min(heff(i,j,idx,nc)/hstar(nc-1),Aleft)
                call movedown(frzpt(i,j),heff(i,j,idx,nc-1),
     $               hseff(i,j,idx,nc-1),
     $               A(i,j,idx,nc-1),hs(nc-1),hi(nc-1),hstar(nc-2),
     $               esdum2,Ts(i,j,idx,nc-1),nilay(nc-1),
     $               heff(i,j,idx,nc),hseff(i,j,idx,nc),A(i,j,idx,nc),
     $               hs(nc),hi(nc),hstar(nc-1),
     $               esdum1,Ts(i,j,idx,nc),nilay(nc),
     $               En(1,nc-1),En(1,nc))
              endif
            enddo
# endif

!-----------------------------------------------------------------------
!     get rid of small open water area
!-----------------------------------------------------------------------

            A0=1.0
            do nc=1,ncat
              A0=A0-max(A(i,j,idx,nc),0.0)
            enddo
            if (A0.lt.A0SMALL) then
!     put tiny open water areas out of misery by reshaping
!     the ice starting with the thin category
              if (A0.ge.0.) then
                do nc=1,ncat
                  if (A(i,j,idx,nc).gt.0.) then
                    A(i,j,idx,nc)=A(i,j,idx,nc)+A0
                    A0=0.0
                  endif
                enddo
              else
!     if the open water fraction is negative then reduce the area
!     of the category with the thin ice first
                nr=1
                do nc=2,ncat
                  if ((A(i,j,idx,nr)+A0).lt.ASMALL(nr)) nr=nc
                enddo
                A(i,j,idx,nr)=A(i,j,idx,nr)+A0
              endif
            endif

!-----------------------------------------------------------------------
!     fix negative ice/snow area/volumes and compute hi, hs, ti
!     if ice volume is too small
!     energy and water in ice-snow will be transferred into open water
!     negative ice or snow volume will increase the temperature
!     of the open water
!-----------------------------------------------------------------------

            htatm = 0.
            wtatm = 0.
            htocn = 0.
            wtocn = 0.

            Aicet=0.0

            do nc=1,ncat
              if ((A(i,j,idx,nc) .le. ASMALL(nc))
     $             .or.(heff(i,j,idx,nc) .le. 0.0)) then

                Enet=0.0
                do layer=1,nilay(nc)
                  Enet=Enet-En(layer,nc)
                enddo

                htocn=htocn+(Enet+hseff(i,j,idx,nc)*rflsno)
                wtocn=wtocn+(heff(i,j,idx,nc)*rhoice
     $                      +hseff(i,j,idx,nc)*rhosno)

                groice(i,j,nc)=-heff(i,j,idx,nc)
                A(i,j,idx,nc)=0.0
                heff(i,j,idx,nc)=0.0
                hseff(i,j,idx,nc)=0.0
                Ts(i,j,idx,nc)=frzpt(i,j)
                do layer=0,nilay(nc)
                  Ti(layer,nc)=frzpt(i,j)
                enddo
                hi(nc)=HSTAR(nc-1)
                hs(nc)=0.0

              else

                groice(i,j,nc)=0.
                hi(nc)=heff(i,j,idx,nc)/A(i,j,idx,nc)
                hi(nc)=max(HSTAR(nc-1),hi(nc))
                A(i,j,idx,nc)=heff(i,j,idx,nc)/hi(nc)

                do layer=1,nilay(nc)
                  qi(layer,nc)=En(layer,nc)*nilay(nc)/heff(i,j,idx,nc)
                enddo
                call getTmp(Ti(1,nc),qi(1,nc),saltz(1,nc),nilay(nc))
                Ti(0,nc)=Ti(1,nc) ! need a starting value for tstm

                if (hseff(i,j,idx,nc) .le. HSSTAR ) then
                  htocn=htocn+hseff(i,j,idx,nc)*rflsno
                  wtocn=wtocn+hseff(i,j,idx,nc)*rhosno
                  hseff(i,j,idx,nc)=0.0
                  hs(nc)=0.0
                else
                  hs(nc)=hseff(i,j,idx,nc)/A(i,j,idx,nc)
                endif
              endif
              Aicet=Aicet+A(i,j,idx,nc)
            enddo

!-----------------------------------------------------------------------
!     finally see if aicet is greater than 1, it never should
!     be after all the trouble above
!-----------------------------------------------------------------------
            if (Aicet.gt.(1.0+TINY)) then
# if defined O_ice_cpts_warnings
              print*,'warning too much ice area by',Aicet
              print*,nats,i,j,(A(i,j,idx,nc),nc=1,ncat)
# endif
              A(i,j,idx,ncat)=min(1.0,A(i,j,idx,ncat))
              Aleft=1.0-A(i,j,idx,ncat)
              do nc=(ncat-1),1,-1
                A(i,j,idx,nc)=min(A(i,j,idx,nc),Aleft)
                Aleft=Aleft-A(i,j,idx,nc)
              enddo
            endif

            A0=1.0
            do nc=1,ncat
              A0=A0-A(i,j,idx,nc)
            enddo

!-----------------------------------------------------------------------
!     heat conduction in ice and change thickness from growth/melt
!     this is only section that needs Ti
!-----------------------------------------------------------------------

            sblmt = 0.
            do nc=1,ncat
              if (A(i,j,idx,nc).ge.TINY) then
!               if snow is less than 25 cm linearly reduce to ice albedo
                as = min(hseff(i,j,idx,nc)*0.04, c1)
                ca = ice_calb*(c1 - as) + sno_calb*as
# if defined O_sulphate_data || defined O_sulphate_data_transient
                dswr = solins(i,j)*sbc(i,j,iaca)*pass
     &                 *max(c0, ca - sulph(i,j,2))
# else
                dswr = solins(i,j)*sbc(i,j,iaca)*pass*ca
# endif
                delhit=0.
                delhib=0.
                delhs=0.
                subi=0.
                subs=0.
                delf=0.
                delfs=0.

                fracsnow = hs(nc)/(hs(nc) + 0.1*centi) ! ~ snow param
                Io=dswr(nc)*0.3*(1-fracsnow) ! solar penetrating surf

                call tstm(tmelz(1,nc),saltz(1,nc),
     $               ts(i,j,idx,nc),ti(0,nc),tocn,tbot(nc),frzpt(i,j),
     $               hi(nc),hs(nc),tair,Qair,Flwatm,
     $               Io,Ib,ug,A(i,j,idx,nc),
     $               dswr(nc),ultnt(nc),usens(nc),ulwr(nc),
     $               nilay(nc),Fneti,condb,i,j)

                call dh(saltz(1,nc),Fx,Fbot,Fneti,ultnt(nc),condb,
     $               tocn,ts(i,j,idx,nc),ti(0,nc),tbot(nc),frzpt(i,j),
     $               qi(1,nc),hi(nc),hs(nc),delhib,
     $               subi,subs,
     $               delhit,delhs,delf,delfs,nilay(nc))

!     note that Ti is no longer current because of layer adjustment

                hi(nc)=hi(nc)+delhit+delhib+subi+delf
                hs(nc)=hs(nc)+delhs+subs-delfs

                groice(i,j,nc)=groice(i,j,nc)+(delhit+delhib)
     $                        *A(i,j,idx,nc)

!     atm gains lwe from sublimation
                evp(nc) = -(rhosno*subs+rhoice*subi)
!     ocn gains lwe from melt/growth
# if defined O_convect_brine
                if (delhib .le. 0.) then
                  wff = -(rhosno*delhs+rhoice*(delhit+delhib))
                else
                  wff = -(rhosno*delhs + rhoice*delhit)
                  if (addflxa) then
                    cbf(i,j,nc) = cbf(i,j,nc)
     $                          - A(i,j,idx,nc)*rhoice*delhib
                    cba(i,j,nc) = cba(i,j,nc) + A(i,j,idx,nc)*dtau
                   endif
                endif
# else
                wff = -(rhosno*delhs+rhoice*(delhit+delhib))
# endif
!     heat flux from atmosphere to ice
                Fice = dswr(nc)-ultnt(nc)-usens(nc)-ulwr(nc)

                htatm = htatm + Fice * dtau * A(i,j,idx,nc)
                sblmt = sblmt + (ultnt(nc) * dtau +
     $                  (rvlsno*subs+rvlice*subi) ) * A(i,j,idx,nc)
                wtatm = wtatm + evp(nc) * A(i,j,idx,nc)
                htocn = htocn + (Fx-Ib)* dtau * A(i,j,idx,nc)
                wtocn = wtocn + wff * A(i,j,idx,nc)

              else
!     give a reasonable value for ts, even if there is no ice
                ts(i,j,idx,nc)=frzpt(i,j)
                dswr(nc) = 0.
                ultnt(nc) = 0.
                usens(nc) = 0.
                ulwr(nc) = 0.
                evp(nc) = 0.
              endif

            enddo               ! loop over bins
            dnswr(i,j)  = 0.
            upltnt(i,j) = 0.
            upsens(i,j) = 0.
            uplwr(i,j)  = 0.
            evap(i,j)   = 0.
            do nc=1,ncat
              dnswr(i,j)  = dnswr(i,j)  + A(i,j,idx,nc)*dswr(nc)
              upltnt(i,j) = upltnt(i,j) + A(i,j,idx,nc)*ultnt(nc)
              upsens(i,j) = upsens(i,j) + A(i,j,idx,nc)*usens(nc)
              uplwr(i,j)  = uplwr(i,j)  + A(i,j,idx,nc)*ulwr(nc)
              evap(i,j)   = evap(i,j)   + A(i,j,idx,nc)*evp(nc)/dtau
            enddo

!-----------------------------------------------------------------------
!     recalc heff hseff and energy of melting for each layer
!-----------------------------------------------------------------------

            do nc=1,ncat
              heff(I,J,idx,nc)=hi(nc)*A(I,J,idx,nc)
              hseff(I,J,idx,nc)=hs(nc)*A(I,J,idx,nc)
              do layer=1,nilay(nc)
                En(layer,nc)=qi(layer,nc)*heff(i,j,idx,nc)/nilay(nc)
              enddo
            enddo

!-----------------------------------------------------------------------
!     ice growth in lead and upper layer water temp
!-----------------------------------------------------------------------
            hnew=0.0
            Tnew=frzpt(i,j)

            if (A0.le.TINY) then                !   do nothing
            elseif ((Focean+Fbot).lt.0.0) then  !   grow new ice
              call grownew(fsh_new,Io_new,Flwatm,Focean,Fbot,
     $             hnew,Tnew,tair,Qair,ug,
     $             dtau,fice,frzpt(i,j),evp(0),
     $             dswr(0),ultnt(0),usens(0),ulwr(0))

              Focean=0.5*(Focean+Fice) ! avg. flx w/ and w/o new ice

              htatm = htatm + Focean*A0*dtau
              htocn = htocn + A0*Fbot*dtau
# if defined O_convect_brine
              if (hnew .le. 0.) then
                wtocn = wtocn - A0*hnew*rhoice
              else
                if (addflxa) then
                  cbf(i,j,0) = cbf(i,j,0) - A0*hnew*rhoice
                  cba(i,j,0) = cba(i,j,0) + A0*dtau
                endif
              endif
# else
              wtocn = wtocn - A0*hnew*rhoice
# endif
              groice(i,j,1)=groice(i,j,1)+hnew*A0

            elseif ((A0.lt.1.0)) then           ! lateral ice melt
!     melt ice lateral
              mrate=m1*(max(tocn-frzpt(i,j),0.0))**m2
              arate=pi_eta*mrate/wclead
              arate=min(arate,1./dtau)
!     concentration that melts laterally per ice conc
              DAice=dtau*arate
              Vfrac=1.0-DAice
              Flat=0.
              Dhlatm=0.
              Dslatm=0.

              do nc=1,ncat
                if (A(i,j,idx,nc).ge.TINY) then

!     diagnostic for heat flux, negative when lateral melt
                  Esum=0.
                  do layer=1,nilay(nc)
                    Esum=Esum+En(layer,nc)
                  enddo
                  Esum=Esum-rflsno*hseff(i,j,idx,nc)
                  Flat=Flat+Esum*arate

!     diagnostic for water flux, positive when lateral melt
                  Dhlatm=Dhlatm+heff(i,j,idx,nc)
                  Dslatm=Dslatm+hseff(i,j,idx,nc)

!     update the ice state
                  heff(i,j,idx,nc)=heff(i,j,idx,nc)*Vfrac
                  hseff(i,j,idx,nc)=hseff(i,j,idx,nc)*Vfrac
                  A(i,j,idx,nc)=A(i,j,idx,nc)*Vfrac
                  do layer=1,nilay(nc)
                    En(layer,nc)=En(layer,nc)*Vfrac
                  enddo
                  groice(i,j,nc)=groice(i,j,nc)-heff(i,j,idx,nc)*DAice
                endif
              enddo
              Dhlatm=Dhlatm*DAice
              Dslatm=Dslatm*DAice

              htocn = htocn - Flat*dtau
              wtocn = wtocn + (Dhlatm*rhoice+Dslatm*rhosno)
            endif

!-----------------------------------------------------------------------
!     budgets:   ice-ocn fluxes of heat/water for coupled model
!     subtract htatm and add wtatm to compensate
!     for adding them back later in embm
!-----------------------------------------------------------------------
! do not move this block
            if (addflxa) then
              flux(i,j,isat) = flux(i,j,isat) - (htocn + htatm)
              flux(i,j,ishum) = flux(i,j,ishum) + (wtatm + wtocn)
            endif

            dnswr(i,j)  = dnswr(i,j)  + dswr(0)*A0
            upltnt(i,j) = upltnt(i,j) + ultnt(0)*A0
            upsens(i,j) = upsens(i,j) + usens(0)*A0
            uplwr(i,j)  = uplwr(i,j)  + ulwr(0)*A0
            evap(i,j)   = evap(i,j)   + evp(0)*A0

!-----------------------------------------------------------------------
!     add lead ice to first category
!-----------------------------------------------------------------------
!     compute the energy in the new ice,
!     divide it into nilay(1) equal portions
            A0=1                       ! must recompute after meltlat
            do nc=1,ncat
              A0=A0-A(i,j,idx,nc)
            enddo
            H0prime=max(hstar(0),hnew) ! new ice h, made >=HOSTAR
            A0prime=A0*hnew/H0prime    ! area covered by H0prime ice

            A1prime=A0prime+A(i,j,idx,1)
            h1prime=0.0
            if (A1prime.gt.0.0) then
              h1prime=(hnew*A0+heff(i,j,idx,1))/A1prime
              Ts(i,j,idx,1)=(Ts(i,j,idx,1)*A(i,j,idx,1)+Tnew*A0prime)
     $                     /A1prime
              A(i,j,idx,1)=A1prime
              hi(1)=h1prime
              heff(i,j,idx,1)=hi(1)*A(i,j,idx,1)
              Ad0=A0*hnew/nilay(1)
              Tmp=0.5*(Tnew+frzpt(i,j))
              E0=Ad0*energ(Tmp,salnew)

              epsilon=Ad0*(energ(Tnew,salnew)-energ(frzpt(i,j),salnew))
              epsilon=0.0
              E0=E0+epsilon/2.

!     sum the energies
              do layer=1,nilay(1)
                En(layer,1)=E0+En(layer,1)-epsilon*(layer-0.5)/nilay(1)
              enddo
            endif               ! finished adding new ice to category 1

!     make sure thickness of category 1 is at least hstar(0)
            heff(i,j,idx,1)=A(i,j,idx,1)*hi(1)
            if (hi(1).gt.hstar(0)) then
              A(i,j,idx,1)=heff(i,j,idx,1)/hi(1)
            else
              A(i,j,idx,1)=heff(i,j,idx,1)/hstar(0)
              hi(1)=hstar(0)
            endif

# if defined O_ice_cpts3 || defined O_ice_cpts5 || defined O_ice_cpts10
!-----------------------------------------------------------------------
!     redistribute the ice categories
!-----------------------------------------------------------------------
!     see if category 1 has outgrown its boundary
!     if so, add it to category 2, and so forth
            do nc=1,(ncat-1)    ! only does loop if ncat>1
              if (hi(nc).ge.hstar(nc)) then
                call moveup(frzpt(i,j),heff(i,j,idx,nc),
     $               hseff(i,j,idx,nc),
     $               A(i,j,idx,nc),hs(nc),hi(nc),hstar(nc-1),
     $               esdum1,Ts(i,j,idx,nc),nilay(nc),
     $               heff(i,j,idx,nc+1),hseff(i,j,idx,nc+1),
     $               A(i,j,idx,nc+1),hs(nc+1),hi(nc+1),hstar(nc),
     $               esdum2,Ts(i,j,idx,nc+1),nilay(nc+1),En(1,nc),
     $               En(1,nc+1))
              endif
            enddo

!     see if the thickest category has become too thin (and so forth)
            do nc=ncat,2,-1
              if (hi(nc).lt.hstar(nc-1)) then
                call movedown(frzpt(i,j),heff(i,j,idx,nc-1),
     $               hseff(i,j,idx,nc-1),
     $               A(i,j,idx,nc-1),hs(nc-1),hi(nc-1),hstar(nc-2),
     $               esdum2,Ts(i,j,idx,nc-1),nilay(nc-1),
     $               heff(i,j,idx,nc),hseff(i,j,idx,nc),A(i,j,idx,nc),
     $               hs(nc),hi(nc),hstar(nc-1),
     $               esdum1,Ts(i,j,idx,nc),nilay(nc),En(1,nc-1),
     $               En(1,nc))
              endif
            enddo

!     cannot hurt to do this again
            heff(i,j,idx,1)=A(i,j,idx,1)*hi(1)
            hi(1)=max(hstar(0),hi(1))
            A(i,j,idx,1)=heff(i,j,idx,1)/hi(1)
# endif

!-----------------------------------------------------------------------
!     reconstruct 3d energy of melting matrix
!-----------------------------------------------------------------------
            l=1
            do nc=1,ncat
              do layer=1,nilay(nc)
                E(i,j,idx,l)=En(layer,nc)
                l=l+1
              enddo
            enddo

!-----------------------------------------------------------------------
!     final energy and water diagnostic
!     could be eliminated to save time
!-----------------------------------------------------------------------
            efinl = 0.
            wfinl = 0.
            do nc = 1,ncat
              do layer = 1,nilay(nc)
                efinl = efinl+En(layer,nc)
              enddo
              efinl = efinl-hseff(i,j,idx,nc)*rflsno
              wfinl = wfinl+(heff(i,j,idx,nc)*rhoice
     $                     +hseff(i,j,idx,nc)*rhosno)
            enddo

            ediff = (efinl-einit)-(htatm+htocn+sblmt)
            wdiff = (wfinl-winit)+(wtatm+wtocn)
# if defined O_ice_cpts_warnings
            if (abs(ediff).gt.maxedif)
     $           print*,'big energy difference in ice model at',i,j,
     $                  ediff*1.e-9,
     $           (htatm+htocn+sblmt)*1.e-9, (efinl-einit)*1.e-9
#  if !defined O_convect_brine
            if (abs(wdiff).gt.maxwdif)
     $           print*,'big water difference in ice model at',i,j,
     $                  wdiff,(wtatm+wtocn), (wfinl-winit)
#  endif
# endif

          endif
        enddo
      enddo

      do nc=1,ncat
        call embmbc (heff(1,1,idx,nc))
        call embmbc (A(1,1,idx,nc))
        call embmbc (Ts(1,1,idx,nc))
        call embmbc (hseff(1,1,idx,nc))
      enddo
      do layer=1,ntilay
        call embmbc (E(1,1,idx,layer))
      enddo

      return
      end

      subroutine tstm(tmz1d, sal1d, ts,ti,tw,tbot,tfrz,
     $     h,hs,Tair,qair,Flo,
     $     Io,Ib,ug,area,
     $     dswr,ultnt,usens,ulwr,
     $     ni,F,condb,i,j)

      implicit none

!     This subroutine calculates the evolution of the ice interior and
!     surface temperature from the heat equation and surface energy
!     balance

!     the albedo is fixed for this calculation
!     if you think you want it to be computed implicitly, then Iabs
!     will have to depend on ts too and that could get nasty
!     I think it would be alot of trouble because the albedo affects
!     boundary conditions too

!     solves the heat equation which is non-linear for saline ice
!     by linearizing and then iterating a set of equations.
!     Scheme is backward Euler giving a tridiagonal
!     set of equations implicit in ti
!     (tried crank-nicholson but it behaves poorly when
!     thin ice has a weird initial temperature profile. This can happen
!     when I redistribute between cats so I always use backward

!     Must have a sufficient amount of snow to solve heat equation in
!     snow hsmin is the minimum depth of snow in order to solve for
!     ti(0) if snow thickness < hsmin then do not change ti(0)

!     The number of equations that must be solved by the tridiagonal
!     solver depends on whether the surface is melting and whether
!     there is snow. Four cases are possible (see variable case below):
!     1 = freezing w/ snow, 2 = freezing w/ no snow,
!     3 = melting w/ snow, and 4 = melting w/ no snow

      include "cpts.h"
      include "thermo.h"
      include "cembm.h"
      include "tmngr.h"

! input variables
      integer ni       ! number of vertical layers in the ice
      integer i,j      ! grid points
      real tmz1d(NMAX) ! melting temp of each layer (C)
      real sal1d(NMAX) ! salinity of each layer (ppt)
      real tfrz        ! freezing temperature of ocean (C)
      real h,hs        ! ice and snow thickness (m)
      real area        ! area of the ice/snow
      real Tw          ! water temp below ice (C)
      real ug          ! geostrophic wind speed (m)
      real qice,qair   ! vapour at ice/snow srfc and in air above
      real tair        ! air temperature (C)
      real Flo         ! long wave down at srfc
      real dswr        ! above srfc net dnwd shortwave, + down (W/m**2)
      real Io          ! solar penetrating top srfc, + down (W/m**2)

! output variables
      real Tbot        ! bottom ice temp (C)
      real condb       ! conductive flux at bottom srfc, + up (W/m**2)
      real Ib          ! solar passing through bottom of ice (W/m**2)
      real F           ! net flux at top srfc including conductive flux
                       ! in ice/snow, + down (W/m**2)
      real usens       ! upward sensible flux, + up (W/m**2)
      real ultnt       ! upward latent flux, + up (W/m**2)
      real ulwr        ! upward longwave flux, + up (W/m**2)

! input/output variables
      real ts          ! srfc temp of snow or ice (C)
      real ti(0:NMAX)  ! snow(0) and ice(1:NMAX) interior temp (C)

! local variables
      real ts_old      ! initial temperature of the ice/snow srfc (C)
      real ti_old(0:NMAX) ! initial temperature in the ice and snow (C)
      real tinext(-1:NMAX) ! incremented ts and ti (C)
      real dti(-1:NMAX) ! incremental changes to ts and ti (C)
      real Ts_kelv     ! srfc temperature in kelvin (K)
      real z           ! vertical coordinate (m)

      real coltnt,cosens ! bulk coef for latent and sensible heat flux

      real condt       ! conductive flux at top srfc, + down (W/m**2)
      real Fo          ! net flux at top srfc excluding conductive flux
                       ! in ice/snow (F-condt), + down (W/m**2)
      real dFo_dt      ! derivative of Fo wrt temperature (W/m**2/K)
      real iru         ! net upward longwave flux, + up (W/m**2)
      real dultnt      ! derivative of ultnt wrt temperature (W/m**2/K)
      real ultnti      ! latent heat from initial temperature (W/m**2)
      real dultnti     ! derivative of ultnti wrt temperature (W/m**2/K)
      logical zult     ! flag that is true if latent heat is zero
      real Iabs(NMAX+1) ! solar absorbed in each layer (W/m**2)
      real absorb      ! sum of Iabs (W/m**2)

      real melts       ! the srfc melting temp (either TSMELT or TMELT)
      real alph, bet   ! parameters for maintaining 2nd order accurate
                       ! diff at boundary

      real tinterf     ! snow and ice interf temp (C)

! a,b,c are vectors that describe the diagonal and off-diagonal elements
! of the matrix [A], such that [A] ti = r
      real a(-1:NMAX),b(-1:NMAX),c(-1:NMAX),d(-1:1),r(-1:NMAX)
      real ki(0:NMAX+1) ! layer conductivity divided by layer thickness
      real zeta(0:NMAX) ! the terms in heat equation independent of ti
      real eta(0:NMAX) ! time step / ice layer thickness / fresh ice cp
      real cpi(1:NMAX) ! time step / ice layer thickness / fresh ice cp
      real dh          ! ice layer thickness
      real dt_dh,dt_hs ! time step / ice or snow thickness
      integer  N       ! number of equations solved by tridiag solver
      integer layer    ! counter for ice layers
      integer ipc      ! counter for iterations of dti
      real errit       ! the absolute value of the maximum dti

      integer case     ! indicates the case of snow vs. no snow and
                       ! melting vs. freezing

      real cond2i,dheat
      logical convrg   ! flag that is true if temperature converges
      logical complt   ! flag that is true if a layer melts internally

      real hsmin
      parameter ( hsmin = 0.01*centi )

      real EXPMAX
      parameter(EXPMAX=10.)

! -------------------------------------------------------------------------
!     setup helpful parameters
! -------------------------------------------------------------------------
      dh    = h/(ni + 1.e-20)
      dt_dh = dtau/(dh + 1.e-20)
      dt_hs = 0.0
      if (hs.gt.hsmin) dt_hs=dtau/(hs + 1.e-20)

      ts_old = ts
      do layer = 0,ni
         ti_old(layer) = ti(layer)
      enddo
      tbot = tfrz     ! min(tw,TMELT)

      call conductiv(sal1d,ki,ti,tbot,hs,dh,ni)

!     the solar radiation absorbed internally
      z = 0.0
      Iabs(1) = 1.0
      Iabs(1) = 1.0
      do layer = 1,ni
        z = z+dh
        if (z.lt.0.10*centi) then ! no absorption in the first 10 cm
          Iabs(layer+1) = 1.0
        elseif (z.lt.0.20*centi) then ! ramp up the absorption
          Iabs(layer+1) = 1.0-(1.0-exp(-kappa*0.20*centi))*
     $         (z/(0.10*centi)-1.0)**2
        else
          Iabs(layer+1) = exp(-kappa*z)
        endif
        Iabs(layer) = Io*(Iabs(layer)-Iabs(layer+1))
      enddo

      coltnt = rslatm*dalt_i*ug
      cosens = 0.94*rcpatm*dalt_i*ug
      qice = QS1*exp(min(AI*ts/(ts-BI+1.e-20),EXPMAX))

      if (qice.lt.qair) then
        zult     = .true.
        qice     = qair
        dultnti  = 0.
        ultnti   = 0.
      else
        zult    = .false.
        dultnti = -coltnt*qice*ai*bi/(ts-bi+1.e-20)**2
        ultnti  = coltnt*(qice-qair)
      endif
      ultnt  = ultnti
      dultnt = dultnti

      usens   = cosens*(ts-tair)
      ts_kelv = ts+TFFRESH
      iru     = ESICE*(Ts_kelv)**4

      dFo_dt  = -cosens-dultnt-4*ESICE*Ts_kelv**3
      Fo      = dswr-Io+Flo-ultnt-usens-iru

      if (hs.gt.hsmin) then
         alph  = 2.0*(2.*hs + dh)/(hs+dh+1.e-20)
         bet   = -2.0*hs*hs/(2.*hs+dh+1.e-20)/(hs+dh+1.e-20)
         melts = tsmelt
      else
         alph  =  3.
         bet   =  -1./3.
         melts = tmelt
      endif

      ts_old = ts
      ts     = min(ts,melts-0.01)   ! absolutely necessary
! -------------------------------------------------------------------------
!     beginning of iterative procedure
! -------------------------------------------------------------------------

      ipc = 1
 1000 continue

! -------------------------------------------------------------------------
!     setup terms that depend on the iterated temperature
! -------------------------------------------------------------------------
      call get_Cpice(sal1d,tmz1d,Cpi,Ti(1),Ti_old(1),ni)
      eta(0) = 0.0
      if (hs.gt.hsmin) eta(0) = dtau/(hs+1.e-20)
      do layer = 1,ni
         eta(layer) = dt_dh/(cpi(layer)+1.e-20)
      enddo

      zeta(0) = rcpsno*Ti_old(0)
      do layer = 1,ni
        zeta(layer) = Ti_old(layer)+eta(layer)*Iabs(layer)
      enddo

      if (ts.lt.melts-tiny) then
!        solve heat equation for ice/snow using flux BC

         if (hs.gt.hsmin) then
! -------------------------------------------------------------------------
            case = 1      ! case of freezing with snow layer
! -------------------------------------------------------------------------

            call getabc(a,b,c,r,ti,tbot,zeta,ki,eta,ni,1)

            a(0) = -eta(0)*ki(0)*(alph+bet)
            c(0) = eta(0)*(bet*ki(0)-ki(1))
            b(0) = rcpsno+eta(0)*(alph*ki(0)+ki(1))
            r(0) = zeta(0)

            a(-1) = 0.
            c(-1) = ki(0)*alph
            d(-1) = ki(0)*bet
            b(-1) = dFo_dt-c(-1)-d(-1)
            r(-1) = -Fo+dFo_dt*ts_old

!           row operation to get rid of d(-1)
            b(-1) = c(0)*b(-1)-d(-1)*a(0)
            c(-1) = c(0)*c(-1)-d(-1)*b(0)
            r(-1) = c(0)*r(-1)-d(-1)*r(0)

            N = ni+2
            call tridag(a(-1),b(-1),c(-1),r(-1),tinext(-1),N)
            dti(-1) = tinext(-1)-ts
            ts = tinext(-1)
            do layer = 0,ni
               dti(layer) = tinext(layer)-ti(layer)
               ti(layer)  = tinext(layer)
            enddo

!            write(*,2000) 'I. iterations',ipc*1.0,ts,(ti(layer),layer = 0,4)
!     $           ,(dti(layer),layer = -1,4)

         else
! -------------------------------------------------------------------------
            case = 2      ! case of freezing with no snow layer
! -------------------------------------------------------------------------

            call getabc(a,b,c,r,ti,tbot,zeta,ki,eta,ni,2)

            a(1) = -eta(1)*ki(1)*(alph+bet)
            c(1) = -eta(1)*(ki(2)-bet*ki(1))
            b(1) = 1+eta(1)*(ki(2)+alph*ki(1))
            r(1) = zeta(1)

            a(0) = 0.
            c(0) = ki(1)*alph
            d(0) = ki(1)*bet
            b(0) = dFo_dt-c(0)-d(0)
            r(0) = -Fo+dFo_dt*ts_old

!           row operation to get rid of d(0)
            b(0) = c(1)*b(0)-d(0)*a(1)
            c(0) = c(1)*c(0)-d(0)*b(1)
            r(0) = c(1)*r(0)-d(0)*r(1)

            N = ni+1
            call tridag(a(0),b(0),c(0),r(0),tinext(0),N)

            dti(-1) = 0.
            dti(0)  = tinext(0)-ts
            ts = tinext(0)
            do layer = 1,ni
               dti(layer) = tinext(layer)-ti(layer)
               ti(layer)  = tinext(layer)
            enddo
            ti(0) = ts

         endif
         ts = min(ts,melts)
      else
         ts = melts
         if (hs.gt.hsmin) then
! -------------------------------------------------------------------------
            case = 3      ! case of melting with snow layer
! -------------------------------------------------------------------------

            a(0) = -eta(0)*ki(0)*(alph+bet)
            c(0) = eta(0)*(bet*ki(0)-ki(1))
            b(0) = rcpsno+eta(0)*(alph*ki(0)+ki(1))
            r(0) = zeta(0)-a(0)*ts
            a(0) = 0.0

            call getabc(a,b,c,r,ti,tbot,zeta,ki,eta,ni,1)

            N = ni+1
            call tridag(a(0),b(0),c(0),r(0),tinext(0),N)

            dti(-1) = 0.
            do layer = 0,ni
               dti(layer) = tinext(layer)-ti(layer)
               ti(layer)  = tinext(layer)
            enddo

!            write(*,2000)'III. T =  ', ipc*1.0,(ti(layer),layer = 0,4)
!     $           ,(dti(layer),layer = 1,4)

         else
! -------------------------------------------------------------------------
            case = 4      ! case of melting with no snow layer
! -------------------------------------------------------------------------

            call getabc(a,b,c,r,ti,tbot,zeta,ki,eta,ni,2)

            a(1) = -eta(1)*ki(1)*(alph+bet)
            c(1) = -eta(1)*(ki(2)-bet*ki(1))
            b(1) = 1+eta(1)*(ki(2)+alph*ki(1))
            r(1) = zeta(1)-a(1)*Ts
            a(1) = 0.0

            N = ni
            call tridag(a(1),b(1),c(1),r(1),tinext(1),N)

            dti(-1) = 0.
            dti(0) = 0.
            do layer = 1,ni
               dti(layer) = tinext(layer)-ti(layer)
               ti(layer)  = tinext(layer)
            enddo

            ti(0) = ts
         endif
      endif

! -------------------------------------------------------------------------
!     end iterative procedure, see if need to reiterate
! -------------------------------------------------------------------------
! -------------------------------------------------------------------------
!     make certain that latent heat is always positive
!     be sure to iterate one more time if latent heat becomes
!     less than zero during the iteration (and was positive before)
! -------------------------------------------------------------------------
      ultnt  =  ultnti + dultnti*(ts-ts_old)
      if (ultnt.lt.0. .and. .not.zult) then
        zult    = .true.
        ultnt   = 0.
        ultnti  = 0.
        dultnti = 0.
        dFo_dt  = -cosens-4*ESICE*Ts_kelv**3
        Fo      = dswr-Io+Flo-usens-iru
        ipc     = ipc+1
        go to 1000
      endif

      errit = 0.0
      do layer = -1,ni
         errit = max(errit,abs(dti(layer)))
      enddo
      if ((errit.lt.errmax).or.(ipc.gt.20)) go to 6000 ! done iterating

      ipc = ipc+1
      go to 1000    ! to beginning of iterative routine

! -------------------------------------------------------------------------
!     continue from here when done iterative
!     begin with error checking
! -------------------------------------------------------------------------
 6000 continue

      complt = .false.
      do layer = 1,ni
         if (ti(layer).gt.tmz1d(layer)) complt = .true.
      enddo

      convrg = .true.
      if (errit.gt.errmax) then
        if (errit.gt.0.005) then
          print*,'WARNING No CONVERGENCE in icemodel',errit,i,j
          convrg = .false.
        else
          if (errit.gt.5.e-4)
     $      print*,'WARNING POOR CONVERGENCE in icemodel',errit,i,j
        endif
      endif

      if (complt .or. .not.convrg) then
        if (complt .and. convrg)
     $       print*,'WARNING CONVERGES to ti>TMELT in icemodel',i,j
        if (complt .and. .not.convrg)
     $       print*,'& final iteration has ti>TMELT'
        print*,'diagnostics useful for debugging'
        print*,area,h,hs,dswr,Flo
        print*,'tstmnew computed T(z):'
        print*,ts,(ti(layer),layer = 0,ni)
        print*,'tstmnew started with T(z):'
        print*,ts_old,(ti_old(layer),layer = 0,ni),tbot

        ts = min(ts_old,melts)
        tinterf=(h*kappas*Ts+hs*kappai*tbot)/(h*kappas+hs*kappai+1.e-20)
        Ti(0) = 0.5*(Ts+tinterf)
        do layer = 1,ni
          Ti(layer) = tinterf+(layer-0.5)*(tbot-tinterf)/(ni+1.e-20)
        enddo
        print*,'setting the temperature profile to be linear T(z):'
        print*,ts,(ti(layer),layer = 0,ni)
      elseif (ti(0).gt.(TSMELT+TINY)) then
        print*,'WARNING SNOW has ti>TSMELT in icemodel',i,j
        print*,'  likely no problem if hs or area is tiny'
        print*,'  area=',area,' hs=',hs
        print*,'setting snow temp to be TSMELT'
        ti(0) = tsmelt
       endif

! -------------------------------------------------------------------------
!     finish up by updating the surface fluxes, etc.
! -------------------------------------------------------------------------

      if (hs.gt.hsmin) then
         condt = ki(0)*(alph*(ti(0)-ts)+bet*(ti(1)-ts))
      else
         condt = ki(1)*(alph*(ti(1)-ts)+bet*(ti(2)-ts))
      endif
      condb = ki(ni+1)*( 3.*(tbot-ti(ni)) - (tbot-ti(ni-1))/3.0  )

      iru   =  iru + (4*ESICE*(ts_old+tffresh)**3)*(Ts-Ts_old)
      ulwr  =  -Flo + iru
      usens =  usens + cosens*(Ts-Ts_old)

      Fo    =  Fo+dFo_dt*(Ts-Ts_old)
      F     =  Fo + condt

!     in the rare event that F<0 set it equal to 0, adjust sensible heat
      if (F .lt. 0.) then
        Fo = -condt
        F  = 0.
        usens = condt+dswr-Io-ulwr-ultnt
      endif

      absorb = 0.
      do layer = 1,ni
         absorb = absorb+Iabs(layer)
      enddo
      Ib = Io-absorb  ! solar passing through the bottom of the ice

! 2000 format(A15, 50(f8.3))

      return
      end

      subroutine getabc(a,b,c,r,ti,tbot,zeta,k,eta,ni,lfirst)

!     compute elements of tridiagonal matrix

      implicit none
      include "cpts.h"
      include "thermo.h"

!  parameters for maintaining 2nd order accurate diff at bottom boundary
      real       alph, bet
      parameter( alph = 3.0 ,bet = -1.0/3.0 )

! input variables
      integer ni              ! number of layers
     $       ,lfirst          ! start with this layer
      real    ti   (0:nmax)   ! temperature of ice-snow layers
     $       ,tbot            ! temperature of ice bottom surface
     $       ,zeta (0:nmax)
     $       ,k    (0:nmax+1) ! ice-snow conductivity
     $       ,eta  (0:nmax)
! output variables
      real    a    (-1:nmax)  ! sub-diagonal elements
     $       ,b    (-1:nmax)  ! diagonal elements
     $       ,c    (-1:nmax)  ! super-diagonal elements
     $       ,r    (-1:nmax)  ! constants (indep. of ti)

! local variable
      integer layer

! if there is snow lfirst = 1 otherwise it is 2
      do layer = lfirst,(ni-1)
        a(layer) = -eta(layer)*k(layer)
        c(layer) = -eta(layer)*k(layer+1)
        b(layer) = 1-c(layer)-a(layer)
        r(layer) = zeta(layer)
      enddo
      a(ni) = -eta(ni)*(k(ni)-bet*k(ni+1))
      c(ni) = 0.0
      b(ni) = 1+eta(ni)*(k(ni)+alph*k(ni+1))
      r(ni) = zeta(ni)+eta(ni)*(alph+bet)*k(ni+1)*tbot

      return
      end

      subroutine tridag(a,b,c,r,u,n)

      implicit none

      include "cpts.h"
      integer nnmax
      parameter (nnmax = nmax+2) ! don't change this

! input variables
      integer n               ! number of layers
      real    a      (nnmax)  ! sub-diagonal elements
     $       ,b      (nnmax)  ! diagonal elements
     $       ,c      (nnmax)  ! super-diagonal elements
     $       ,r      (nnmax)  ! constants (indep. of ti)

! output variables
      real    u      (nnmax)  ! solution

! local variables
      integer layer
      real    bet             ! work variable
     $       ,gam    (nnmax)  ! work array

      bet = b(1)
      u(1) = r(1)/bet
      do layer = 2,n
        gam(layer) = c(layer-1)/bet
        bet = b(layer)-a(layer)*gam(layer)
        u(layer) = (r(layer)-a(layer)*u(layer-1))/bet
      enddo
      do layer = n-1,1,-1
        u(layer) = u(layer)-gam(layer+1)*u(layer+1)
      enddo

      return
      end

      subroutine conductiv(sal1d,ki,ti,tbot,hs,dh,ni)

      implicit none

      include "cpts.h"
      include "thermo.h"

      real sal1d    (NMAX) ! salinity of each layer
      real ki(0:nmax+1),ti(0:nmax),tbot,hs,dh
      real tmax
      parameter (tmax = -0.1)
      integer ni,layer

      real hsmin
      parameter ( hsmin = 0.01*centi )

      do layer = 2,ni
         ki(layer) = kappai+BETA*(sal1d(layer)+sal1d(layer+1))/2.
     $        /min(tmax,(ti(layer-1)+ti(layer))/2.)
         ki(layer) = max(ki(layer),KIMIN)
         ki(layer) = ki(layer)/dh
      enddo
      ki(ni+1) = kappai+BETA*sal1d(ni+1)/tbot
      ki(ni+1) = max(ki(ni+1),KIMIN)
      ki(ni+1) = ki(ni+1)/dh
      ki(1) = kappai+BETA*sal1d(1)/min(tmax,ti(1))
      ki(1) = max(ki(1),KIMIN)
      ki(1) = ki(1)/dh
!cc      ki(1) = kappai/dh
      if (hs.gt.hsmin) then
         ki(0) = kappas/hs
         ki(1) = 2.0*ki(1)*ki(0)/(ki(1)+ki(0))
      else
         ki(0) = 0.0
      endif

      return
      end

      subroutine get_Cpice(sal1d,tmz1d,cpi,Tnew,Told,ni)

!     compute the heat capacity of each layer in the ICE
!     Tnew is assumed to be the midpoint temperature of each layer

      implicit none

      include "cpts.h"
      include "thermo.h"

      real sal1d    (NMAX) ! salinity of each layer
      real tmz1d    (NMAX) ! melting temp of each layer
      real cpi(nmax),Tnew(nmax),Told(nmax)
      integer ni,layer

      do layer = 1,ni
        cpi(layer) = rcpice+gamma*sal1d(layer)/Told(layer)/
     $       min(Tnew(layer),tmz1d(layer))
      enddo

#endif

      return
      end