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):

https://github.com/CICE-Consortium/Icepack/blob/f5e093f5148554674079d5c7fc0702a41b81f744/columnphysics/icepack_therm_vertical.F90#L752-L763

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 https://github.com/CICE-Consortium/CICE-svn-trunk/commit/147e57de0c22410352b40877d33dffa80f51583c#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:

https://github.com/CICE-Consortium/Icepack/blob/f5e093f5148554674079d5c7fc0702a41b81f744/columnphysics/icepack_itd.F90#L1513-L1521

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.