Correction for thin ice/snow treatment of enthalpy and other tracers

To be tested as a possible fix to #5124 - mpas-si crashes on small timesteps bug fix

For snow/ice thicknesses < puny but > 0, enthalpy and other adjusted tracers are spread evenly among layers.
Non-BFB

components/mpas-seaice/src/column/ice_therm_vertical.F90

@@ -1710,8 +1710,8 @@ subroutine thickness_changes (nilyr,     nslyr,    &
       !-----------------------------------------------------------------
 
       if (ktherm == 2) then
-         do k = 1, nslyr
-            if (hsn <= puny) then
+         if (hsn <= puny) then
+            do k = 1, nslyr
                fhocnn = fhocnn &
                       + zqsn(k)*hsn/(real(nslyr,kind=dbl_kind)*dt)
                zqsn(k) = -rhos*Lfresh
@@ -1720,9 +1720,12 @@ subroutine thickness_changes (nilyr,     nslyr,    &
                  smice(k) = c0
                  smliq(k) = c0
                endif
-               hslyr = c0
-            endif
-         enddo
+               if (tr_rsnw) rsnw(k) = rsnw_fall
+            enddo
+            melts = melts + hsn
+            hsn = c0
+            hslyr = c0
+         endif
       endif
 
       !-----------------------------------------------------------------
@@ -1903,7 +1906,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
@@ -1913,36 +1917,53 @@ subroutine adjust_enthalpy (nlyr,               &
       !-----------------------------------------------------------------
 
       rhlyr = c0
-      if (hn > puny) rhlyr = c1 / hlyr
+      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
+         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
+         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
+            do k = 1, nlyr
+               qn(k) = c0
+            enddo
+         endif
+
+      endif
 
       end subroutine adjust_enthalpy