Tmax in icepack_itd code

[DmitryDukhovskoy]

Hello,

While looking into the icepack_itd.F90 code, I'm having trouble to get right units in the line whereTmax is computed from snow enthalpy:

        if (hsn > hs_min) then
           ! zqsn < 0
           zqsn = trcrn(nt_qsno+k-1,n)
           Tmax = -zqsn*puny*rnslyr / (rhos*cp_ice*vsnon(n))
        else
           zqsn = -rhos * Lfresh
           Tmax = puny
        endif

For Tmax = -zqsn * puny * rnslyr / (rhos* cp_ice* vsnon(n))
We have: zqsn = J/m3, puny, rnslyr - no units
rhos = kg/m3
cp_ice = J/kg/degree
vsnon = m3/m2 = m

This gives Tmax = [ J/m3 / J/(degree *m2) ] = degree/m

I wonder if there is a bug in this formula and the extra term (rnslyr/vsnon(n)) was copied mistakenly from the formula that converts CICE4 snow enthalpy (J/m2) to CICE6 that is:
trcrn(i,j,nt_qsno+k-1,n) = &
min(esnon(i,j,slyr1(n)+k-1)rnslyr/vsnon(i,j,n),-rhosLfresh)

If rnslyr/vsnon(n) is removed from the equation, the units work fine:

Tmax * rhos * cp_ice = -zqsn ==> J/m3 = J/m3

What am I missing?

Thank you

[apcraig]

You could also try asking this question here, https://bb.cgd.ucar.edu/cesm/forums/cice-consortium.146/

[phil-blain]

I think the units are OK, as the comment above indicates that puny is used as a small volume, and so it has units of meters (m^3/m^2):

Icepack/columnphysics/icepack_therm_vertical.F90

Lines 752 to 763 in f5e093f

	!-----------------------------------------------------------------
	! Tmax based on the idea that dT ~ dq / (rhos*cp_ice)
	! dq ~ q dv / v
	! dv ~ puny = eps11
	! where 'd' denotes an error due to roundoff.
	!-----------------------------------------------------------------
	if (hslyr > hs_min/rnslyr) then
	! zqsn < 0
	Tmax = -zqsn(k)*puny*rnslyr / &
	(rhos*cp_ice*vsnon)
	else

so this gets rid of your remaining meter and so Tmax has units of degrees.

I think that the rnslyr factor is also OK, it was added between CICE4 and CICE5, in CICE-Consortium/CICE-svn-trunk@147e57d#diff-0dea52c09d4accde92a7a7f9a321db395f5d2e0b506cd24fb9949b1a35acfb8eL969 and the commit messages indicates:

A bug fix in
ice_therm_vertical only affects simulations using more than one snow layer
(the default is one layer).
source/ice_therm_vertical.F90

• BUG FIX: correct layer thickness for multiple snow layers

[phil-blain]

I see you are in fact referring to icepack_itd.F90:

Icepack/columnphysics/icepack_itd.F90

Lines 1513 to 1521 in f5e093f

	! snow enthalpy and max temperature
	if (hsn > hs_min) then
	! zqsn < 0
	zqsn = trcrn(nt_qsno+k-1,n)
	Tmax = -zqsn*puny*rnslyr / (rhos*cp_ice*vsnon(n))
	else
	zqsn = -rhos * Lfresh
	Tmax = puny
	endif

That code seems to do the same thing as in icepack_therm_vertical.F90 which I point out above.

[DmitryDukhovskoy]

Thank you Phil, I thought about puny as a very small volume to make the units work but was not sure if this was correct. Your explanation makes sense.

[DmitryDukhovskoy]

I'm still struggling with zqsn units in icepack_itd.F90. In the code, zqsn is defined as J/m2 (snow layer enthalpy) - units used in CICE4.
1515 zqsn , & ! snow layer enthalpy (J m-2)

For Tmax the formula seems to be ok, assuming puny = [m3] (a tiny volume of snow) and using rnslyr/vsnon to convert J/m2 --> J/m3 (internal energy of a snow layer).

Next, the code computes snow T using zqsn:
1552 ! snow temperature
1553 zTsn = (Lfresh + zqsn/rhos)/cp_ice

Here, units of zqsn should be J/m3 not J/m2 (Lfresh = [J/kg], rhos = [kg/m3], cp_ice = [J/(kg * degree)]). I thought that zqsn should be J/m3 - units for internal snow energy (that's what is saved in restart) and I could not find in the code, where initial energy is converted to J/m2 (enthalpy) before being passed to icepack_itd.
zqsn is taken from trcrn array:
1545 zqsn = trcrn(nt_qsno+k-1,n) (should be J/m3)
or computed as:
1548 zqsn = -rhos * Lfresh (J/m3)

In zap_snow subroutine, trcrn is converted to snow enthalpy (J/m2):
1447 xtmp = trcrn(nt_qsno+k-1) / dt &
1448 * vsnon/real(nslyr,kind=dbl_kind) ! < 0

What am I missing? Thank you!

[eclare108213]

Hi Dmitry,
This is confusing, I agree. Looking through the code, I believe that the qsno tracer (or zqsn) should be J/m3 everywhere. qsno is fundamentally Lfresh*rhos ~ J/kg * kg/m3 ~ J/m3. It is listed as J/m2 in the flood_ice routine in icepack_therm_mushy.F90, which is clearly wrong (they are passed in as J/m3). The only other place in Icepack that they are listed as J/m2 is in icepack_itd.F90, and I think the units comment is wrong here also. As @phil-blain notes above, the comments in icepack_therm_vertical.F90 regarding Tmax indicate that puny does have units for this approximation.
Does the code make sense if you assume the J/m2 units in the comments are wrong?

[DmitryDukhovskoy]

Hi Elizabeth, thank you for clarification and confirming the units of the qsno (zqsn) tracer. The code does make sense for J/m3 units.