Corrects thin ice/snow treatment of enthalpy and other tracers (CICE-…

…Consortium#454)

* Corrects thin ice/snow treatment of enthalpy and other tracers

Adopted from E3SM-Project/E3SM#5630

"This fix redistributes enthalpy and other tracers evenly in snow and ice when their
respective thicknesses are < 1e-15 . Previously, these tracers were zero'd non-conservatively.
Also corrects a bug in the zeroing of snow thicknesses, and removes snow in thickness_changes
if the ice vanishes."

* Clean up modifications

columnphysics/icepack_therm_shared.F90

@@ -497,7 +497,8 @@ subroutine adjust_enthalpy (nlyr,               &
          hovlp           ! overlap between old and new layers (m)
 
       real (kind=dbl_kind) :: &
-         rhlyr           ! 1./hlyr
+         rhlyr,        & ! 1./hlyr
+         qtot            ! total h*q in the column
 
       real (kind=dbl_kind), dimension (nlyr) :: &
          hq              ! h * q for a layer
@@ -509,36 +510,55 @@ subroutine adjust_enthalpy (nlyr,               &
       !-----------------------------------------------------------------
 
       rhlyr = c0
-      if (hn > puny) rhlyr = c1 / hlyr
-
-      !-----------------------------------------------------------------
-      ! Compute h*q for new layers (k2) given overlap with old layers (k1)
-      !-----------------------------------------------------------------
-
-      do k2 = 1, nlyr
-         hq(k2) = c0
-      enddo                     ! k
-      k1 = 1
-      k2 = 1
-      do while (k1 <= nlyr .and. k2 <= nlyr)
-         hovlp = min (z1(k1+1), z2(k2+1)) &
-               - max (z1(k1),   z2(k2))
-         hovlp = max (hovlp, c0)
-         hq(k2) = hq(k2) + hovlp*qn(k1)
-         if (z1(k1+1) > z2(k2+1)) then
-            k2 = k2 + 1
+      if (hn > puny) then
+         rhlyr = c1 / hlyr
+
+         !-----------------------------------------------------------------
+         ! Compute h*q for new layers (k2) given overlap with old layers (k1)
+         !-----------------------------------------------------------------
+
+         do k2 = 1, nlyr
+            hq(k2) = c0
+         enddo                     ! k
+         k1 = 1
+         k2 = 1
+         do while (k1 <= nlyr .and. k2 <= nlyr)
+            hovlp = min (z1(k1+1), z2(k2+1)) &
+                  - max (z1(k1),   z2(k2))
+            hovlp = max (hovlp, c0)
+            hq(k2) = hq(k2) + hovlp*qn(k1)
+            if (z1(k1+1) > z2(k2+1)) then
+               k2 = k2 + 1
+            else
+               k1 = k1 + 1
+            endif
+         enddo                  ! while
+
+         !-----------------------------------------------------------------
+         ! Compute new enthalpies.
+         !-----------------------------------------------------------------
+
+         do k = 1, nlyr
+            qn(k) = hq(k) * rhlyr
+         enddo                     ! k
+
+      else
+
+         qtot = c0
+         do k = 1, nlyr
+            qtot = qtot + qn(k) * (z1(k+1)-z1(k))
+         enddo
+         if (hn > c0) then
+            do k = 1, nlyr
+               qn(k) = qtot/hn
+            enddo
          else
-            k1 = k1 + 1
+            do k = 1, nlyr
+               qn(k) = c0
+            enddo
          endif
-      enddo                  ! while
 
-      !-----------------------------------------------------------------
-      ! Compute new enthalpies.
-      !-----------------------------------------------------------------
-
-      do k = 1, nlyr
-         qn(k) = hq(k) * rhlyr
-      enddo                     ! k
+      endif
 
       end subroutine adjust_enthalpy
 

columnphysics/icepack_therm_vertical.F90

@@ -1706,7 +1706,7 @@ subroutine thickness_changes (nilyr,     nslyr,    &
       !-----------------------------------------------------------------
 
       if (ktherm == 2) then
-         if (hsn <= puny) then
+         if (hsn <= puny .or. hin <= c0) then
             do k = 1, nslyr
                fhocnn = fhocnn &
                       + zqsn(k)*hsn/(real(nslyr,kind=dbl_kind)*dt)
@@ -1715,8 +1715,11 @@ subroutine thickness_changes (nilyr,     nslyr,    &
                   meltsliq = meltsliq + massice(k)  ! add to meltponds
                   smice(k) = rhos
                   smliq(k) = c0
+                  rsnw(k) = rsnw_fall
                endif
             enddo
+            melts = melts + hsn
+            hsn = c0
             hslyr = c0
          endif
       endif