From 7e0821dcfe824ae99a3bf0cb957497a8e10c4da7 Mon Sep 17 00:00:00 2001
From: apcraig <anthony.p.craig@gmail.com>
Date: Mon, 28 Aug 2023 17:24:46 -0600
Subject: [PATCH] Corrects thin ice/snow treatment of enthalpy and other
 tracers

Adopted from https://github.com/E3SM-Project/E3SM/pull/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."
---
 columnphysics/icepack_therm_shared.F90   | 76 +++++++++++++++---------
 columnphysics/icepack_therm_vertical.F90 | 14 ++++-
 2 files changed, 61 insertions(+), 29 deletions(-)

diff --git a/columnphysics/icepack_therm_shared.F90 b/columnphysics/icepack_therm_shared.F90
index ad00db2ac..53113c08d 100644
--- a/columnphysics/icepack_therm_shared.F90
+++ b/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
 
diff --git a/columnphysics/icepack_therm_vertical.F90 b/columnphysics/icepack_therm_vertical.F90
index 09328c5d3..98acbdcb9 100644
--- a/columnphysics/icepack_therm_vertical.F90
+++ b/columnphysics/icepack_therm_vertical.F90
@@ -1706,17 +1706,29 @@ 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)
                zqsn(k) = -rhos*Lfresh
+!tcx, tcraig, in columnphysics, this is 
+! is it correct that now everything is "if snwgrain"?
+!               if (tr_snow) then
+!                 meltsliq = meltsliq + smicetot(k)  ! add to meltponds
+!                 smice(k) = rhos
+!                 smliq(k) = c0
+!               endif
+!               if (tr_rsnw) rsnw(k) = rsnw_fall
                if (snwgrain) then
                   meltsliq = meltsliq + massice(k)  ! add to meltponds
                   smice(k) = rhos
                   smliq(k) = c0
+                  rsnw(k) = rsnw_fall
                endif
+!tcx
             enddo
+            melts = melts + hsn
+            hsn = c0
             hslyr = c0
          endif
       endif