! source file: /Users/oschlies/UVIC/master/source/ice/therm.F
      subroutine therm  (is, ie, js, je)
!=======================================================================
!     thermodynamic ice model

!     Note: if run with embm this routine must be called after "fluxes"
!           and before "solve"

!     calculates ice and open water growth rates based on the surface
!     energy budget. see Parkinson and Washington, JGR, Vol.84, C1,
!     311-337, 1979 and Hibler, JPO, Vol.9, 815-846, 1979

!     heat and fresh water fluxes between the ocean and atmosphere
!     are adjusted depending on ice growth or melt. ice thickness is
!     changed by the amount of growth or melt and sublimation

!     based on code by: M. Eby
!=======================================================================

      implicit none

      include "param.h"
      include "csbc.h"
      include "cembm.h"
      include "atm.h"
      include "ice.h"
      include "coord.h"
      include "grdvar.h"
      include "veg.h"
      include "mtlm.h"

      integer i, ie, iem1, imax, index, iter, is, isp1, j, je, jem1
      integer jmax, js, jsp1, maxit

      real ai, aice3, al, amin, ao, as, as_crops, ca,  delta, df, dh
      real dha, dhflxi, dhflxs, dhi, hice3, dho, dhs, hsno3, dhss
      real dhstot, dhtot, dswr, dt, dultnt, dulwr, dusens, emax, f, fa
      real fas, fb, fbot, fcond, fd, fds, fe, ff, ffs, fh, fl, fls, fm
      real fn, fptf, fpts, ftopi, ftopo, ho, hsextra, qair,qice, sla
      real sub, tair, tcdh, ti, tiold, tol, tol2, ultnt, ulwr, usens
      real wspd, zintfc

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

      fa = dts/(rhoice*flice)
      fb = 0.94*rhoatm*cpatm
      fd = rhoatm/rhoice
      fe = rhoatm*slice
      ff = rhoice*flice
      fh = 21.8746*265.5
      fas = dts/(rhosno*flice)
      fds = rhoatm/rhosno
      ffs = rhosno*flice
      ho = 1.0
      amin = 0.15
      maxit = 10
      tol = 0.01
      emax = 0.0
      imax = 0
      jmax = 0
      fptf = 0.0
      index = 1
      sla = zw(1)*secday/dampice/2.389e-8
      tol2 = tol*2.0

      do j=jsp1,jem1
        do i=isp1,iem1

          hsno3 = hsno(i,j,index)
          hice3 = hice(i,j,index)
          aice3 = 0

!-----------------------------------------------------------------------
!         set the incoming shortwave over snow and ice
!-----------------------------------------------------------------------
!         if snow is less than 25 cm linearly reduce to ice albedo
          as = min(hsno(i,j,2)*0.04/(aice(i,j,2) + epsln), c1)
          ca = ice_calb*(c1 - as) + sno_calb*as
          dswr = solins(i,j)*sbc(i,j,iaca)*pass*ca
          wspd = sbc(i,j,iws)

          aice(i,j,index) = min(c1, aice(i,j,index))
          ai = aice(i,j,2)

          if (tmsk(i,j) .lt. 0.5 .and. land_map(i,j) .eq. 0) then

!-----------------------------------------------------------------------
!           land points
!-----------------------------------------------------------------------
            tair = at(i,j,2,isat) - elev(i,j)*rlapse
     &           - hicel(i,j,2)*rlapse
            aice3 = 0.0
!           snow masking distance for different vegetation types
            as = hsno(i,j,2)/(100.0*veg_smd(iveg(i,j)))
            as_crops = hsno(i,j,2)/(100.0*veg_smd(icrops))
            as_crops = min(c1, max(c0, as_crops))
            as = (1.-crops(i,j,2))*as + crops(i,j,2)*as_crops
!           limit snow coverage between 0 and 1
            aice3 = min(c1, max(c0, as))
            aice3 = max(aice3, aicel(i,j,2))

!-----------------------------------------------------------------------
!           start loop for grid points with snow or ice
!-----------------------------------------------------------------------
            if (ai .gt. c0) then

!-----------------------------------------------------------------------
!             find snow temperature by balancing the surface heat budget
!               dwsr = ultnt + usens + ulwr
!             using Newton's method:
!               t(i+1) = t(i) - f(t(i))/df(t(i))
!             where:
!               f(t(i)) = dwsr - ultnt - usens - ulwr
!               -df(t(i)) = dultnt - dusens - dulwr
!-----------------------------------------------------------------------
              al = 1.0 - ai
              ti = tice(i,j)
              tiold = tice(i,j)
              fm = esatm*(tair + 273.15)**4
              fls = fe*dalt_i*wspd
              dusens = fb*dalt_i*wspd
              qair = rh(i,j)*cssh*exp(17.67*tair/(tair + 243.5))
              iter = 0
              delta = tol2
              do while (abs(delta) .gt. tol .and. iter .le. maxit)
                iter = iter + 1
                dt = ti - tair
                qice = cssh*exp(21.8746*ti/(ti + 265.5))
                if (qice .gt. qair) then
                  ultnt = fls*(qice - qair)
                  dultnt = fls*qice*21.8746*265.5/(ti + 265.5)**2
                else
                  ultnt = 0.0
                  dultnt = 0.0
                endif
                usens = dusens*dt
                ulwr = esice*(ti + 273.15)**4 - fm
                dulwr = 4.0*esice*(ti + 273.15)**3
                f = dswr - ultnt - usens - ulwr
                df = dultnt + dusens + dulwr
                delta = f/df
                ti = ti + delta
              enddo
              if (iter .gt. maxit) then
!               if not converged, set to last converged temperature
                if (abs(delta) .gt. emax .and. ti .lt. fptf) then
                  emax = abs(delta)
                  imax = i
                  jmax = j
                  ti = tiold
                endif
              endif

!-----------------------------------------------------------------------
!             set maximum tice to freezing and calculate fluxes
!-----------------------------------------------------------------------
              ti = min(ti, fptf)
              dt = ti - tair
              qice = cssh*exp(21.8746*ti/(ti + 265.5))
              sub = max(c0, fds*dalt_i*wspd*(qice - qair))
              usens = dusens*dt
              ulwr = esice*(ti + 273.15)**4 - fm

!             ensure that snow sublimated does not exceed hsno
              dha = -dts*sub*ai
              if (-dha .gt. hsno3) then
                dha = -hsno3
                sub = -dha/(ai*dts)
              endif
                     ultnt = rhosno*slice*sub
              sub = sub*rhosno

              ftopi = dswr - ulwr - usens - ultnt

!-----------------------------------------------------------------------
!             add ice/snow covered area fluxes to land fluxes
!-----------------------------------------------------------------------
              tice(i,j) = ti
              evap(i,j) = evap(i,j)*al + sub*ai
              dnswr(i,j) = dnswr(i,j)*al + dswr*ai
              upltnt(i,j) = upltnt(i,j)*al + ultnt*ai
              upsens(i,j) = upsens(i,j)*al + usens*ai
              uplwr(i,j) = uplwr(i,j)*al + ulwr*ai

!-----------------------------------------------------------------------
!             calculate total change in snow volume on land
!             allocate snowmelt to flux for input to bucket
!-----------------------------------------------------------------------
              dhs = 0.0
              dhtot = 0.0
              if (tice(i,j) .ge. fptf .and. ftopi .gt. 0.0)
     &          dhs = -ai*fas*ftopi
              dhs = min(0.0, max(-(hsno3 + dha), dhs))
              dhtot = dhs + dha
              hsno3 = hsno3 + dhtot

              flux(i,j,ishum) = flux(i,j,ishum) + dhs*rhosno/dts

!             ensure fluxes are balanced since land can't absorb error
              upsens(i,j) = dnswr(i,j) - upltnt(i,j) - uplwr(i,j)
     &                   + dhs*ffs/dts

!-----------------------------------------------------------------------
!     if no snow on land
!-----------------------------------------------------------------------
            else
              tice(i,j) = 0.0
            endif

          elseif (tmsk(i,j) .ge. 0.5) then

!-----------------------------------------------------------------------
!           ocean points
!-----------------------------------------------------------------------
            ao = 1.0 - ai

!-----------------------------------------------------------------------
!           calculate fluxes to and from the ocean (without ice)
!-----------------------------------------------------------------------
            ftopo = dnswr(i,j) - uplwr(i,j) - upsens(i,j) - upltnt(i,j)
            fbot = sla*(frzpt(i,j) - sbc(i,j,isst))
!           calculate growth of ice in open water areas
            dho = fa*(fbot - ftopo)

            if (ai .ne. 0.0) then

!-----------------------------------------------------------------------
!             find ice temperature by balancing the surface heat budget:
!               tcdh*(ti - fpts) = dswr - ultnt - usens - ulwr
!             using Newton's method:
!               t(i+1) = t(i) - f(t(i))/df(t(i))
!             where:
!               f(t(i)) = dswr - ultnt - usens - ulwr - tcdh*(ti - fpts)
!               -df(t(i)) = dultnt + dusens + dulwr + tcdh
!-----------------------------------------------------------------------
              tcdh = condice/(hice(i,j,2) + 6.5*hsno(i,j,2))
              ti = tice(i,j)
              tiold = tice(i,j)
              fpts = frzpt(i,j)
              fm = esatm*(at(i,j,2,isat) + 273.15)**4
              fn = 4.0*esice

              fl = fe*dalt_i*wspd
              dusens = fb*dalt_i*wspd
              qair = at(i,j,2,ishum)
              iter = 0
              delta = tol2
              do while (abs(delta) .gt. tol .and. iter .le. maxit)
                iter = iter + 1
                dt = ti - at(i,j,2,isat)
                qice = cssh*exp(21.8746*ti/(ti + 265.5))
                if (qice .gt. qair) then
                  ultnt = fl*(qice - qair)
                  dultnt = fl*qice*fh/(ti + 265.5)**2
                else
                  ultnt = 0.0
                  dultnt = 0.0
                endif
                usens = dusens*dt
                ulwr = esice*(ti + 273.15)**4 - fm
                dulwr = fn*(ti + 273.15)**3
                f = dswr - ultnt - usens - ulwr - tcdh*(ti - fpts)
                df = dultnt + dusens + dulwr + tcdh
                delta = f/df
                ti = ti + delta
              enddo
              if (iter .gt. maxit) then
!               if not converged, set to last converged temperature
                if (abs(delta) .gt. emax .and. ti .lt. fptf) then
                  emax = abs(delta)
                  imax = i
                  jmax = j
                  ti = tiold
                endif
              endif

!-----------------------------------------------------------------------
!             set maximum tice to freezing and calculate fluxes
!-----------------------------------------------------------------------
              ti = min(ti, fptf)
              dt = ti - at(i,j,2,isat)
              qice = cssh*exp(21.8746*ti/(ti + 265.5))
              sub = max(c0, dalt_i*wspd*(qice - qair))
              ultnt = fe*sub
              fcond = tcdh*(ti - fpts)
              if (hsno(i,j,index) .gt. 0.0) then
                sub = fds*sub
                dha = -dts*sub
                sub = sub*ai*rhosno
              else
                sub = fd*sub
                dha = -dts*sub
                sub = sub*ai*rhoice
              endif

              usens = dusens*dt
              ulwr = esice*(ti + 273.15)**4 - fm

!-----------------------------------------------------------------------
!             add ice covered area fluxes to ocean area fluxes
!-----------------------------------------------------------------------
              tice(i,j) = ti
              dnswr(i,j) = dnswr(i,j)*ao + dswr*ai
              upltnt(i,j) = upltnt(i,j)*ao + ultnt*ai
              upsens(i,j) = upsens(i,j)*ao + usens*ai
              uplwr(i,j) = uplwr(i,j)*ao + ulwr*ai
              ftopi = dswr - ulwr - usens - ultnt

!-----------------------------------------------------------------------
!             calculate change in ice volume due to sublimation of
!             ice (dha). adjust evaporation to the atmosphere to
!             account for sublimation from ice. subtract this
!             adjustment from the ocean freshwater flux
!-----------------------------------------------------------------------
              if (addflux) flux(i,j,ishum) = flux(i,j,ishum) + dts*sub
              evap(i,j) = evap(i,j)*ao + sub
            else
              tice(i,j) = sbc(i,j,isst)
              ftopi = 0.0
              fcond = 0.0
              dha = 0.0
            endif

!-----------------------------------------------------------------------
!           calculate total change in ice volume (dh)
!-----------------------------------------------------------------------
            dha = dha*ai
            dhflxs = 0.0
            dhs = 0.0
            if (hsno(i,j,index) .le. 0.0) then
!             total growth of ice from the ocean
              dhi = ai*fa*(fbot - ftopi) + ao*dho
!             total growth (loss + sublimation limited to total amount)
              dh = max(-hice(i,j,index), dhi + dha)
!             adjust ocean fluxes for ice growth or melt + sublimation
              dhflxi = dh - dha
            else
!             total growth of ice from the ocean
              dhi = ai*fa*(fbot - fcond) + ao*dho
!             loss of snow due to melt
              if (tice(i,j) .ge. fptf) dhs = ai*fas*(fcond - ftopi)
!             total loss of snow including sublimation
              dhs = dhs + dha
!             check if snow loss greater than total
              if (-dhs .gt. hsno(i,j,index)) then
!               take extra melt from ice
                dhi = dhi + rhosno/rhoice*(dhs + hsno(i,j,index))
!               remove all snow
                dhs = -hsno(i,j,index)
              endif
!             adjust ocean fluxes for snow melt and sublimation
              dhflxs = dhs - dha
!             ice loss limited to total
              dh = max(-hice(i,j,index), dhi)
!             adjust ocean fluxes for ice growth or melt
              dhflxi = dh
            endif

!-----------------------------------------------------------------------
!           calculate new area and thickness from thermodynamics
!-----------------------------------------------------------------------
!           use minimum area (amin) of open water
            ai = max(amin, aice(i,j,index))
            aice3 = aice(i,j,index) + ((1.0 - ai)*max(c0, dho)/ho
     &         + 0.5*min(c0, dhi)*ai/(hice(i,j,index) + epsln))
!           update ice and snow thickness
            hice3 = hice(i,j,index) + dh
            hsno3 = hsno(i,j,index) + dhs

!           lower ice area where thickness is < 1 cm
            aice3 = min(aice3, hice3)
!           max ice area where thickness is > 10 m
            aice3 = max(aice3, hice3*0.001)
            aice3 = max(c0, min(c1, aice3))
            if (aice3 .eq. 0.0) then
              dhflxs = dhflxs - hsno3
              hsno3 = 0.0
            endif

!           check if the weight of the snow pushes the ice/snow
!           interface below the waterline (if so, change snow to ice)
            zintfc = hice3 - (rhosno*hsno3 + rhoice*hice3)/rhoocn
            if (zintfc .lt. 0.0) then
              dhss = rhoice/rhosno*zintfc
              if (-dhss .gt. hsno3) then
                write(*,*) '==> Warning: dhss is too large: ',dhss
                dhss = -hsno3
              endif
              hice3 = hice3 - rhosno/rhoice*dhss
              hsno3 = hsno3 + dhss
            endif
            hsno3 = max(hsno3, c0)

!-----------------------------------------------------------------------
!           adjust fluxes to the ocean due to ice melt or growth
!-----------------------------------------------------------------------
            if (addflux) then
              flux(i,j,isat) = flux(i,j,isat) + ff*dhflxi + ffs*dhflxs
              flux(i,j,ishum) = flux(i,j,ishum) - rhoice*dhflxi
     &                        - rhosno*dhflxs
            endif

          endif

!-----------------------------------------------------------------------
!         shuffle time levels
!-----------------------------------------------------------------------
          hice(i,j,1) = hice(i,j,2)
          hice(i,j,2) = hice3
          aice(i,j,1) = aice(i,j,2)
          aice(i,j,2) = aice3
          hsno(i,j,1) = hsno(i,j,2)
          hsno(i,j,2) = hsno3

        enddo
      enddo

      if (emax .gt. 0.0) write (stdout,*)
     &  '==> Warning: ice temperature not converging: emax, i, j:'
     &, emax, imax, jmax

!-----------------------------------------------------------------------
!     set boundary conditions
!-----------------------------------------------------------------------
      call embmbc (hice(1,1,2))
      call embmbc (aice(1,1,2))
      call embmbc (hsno(1,1,2))

      return
      end