Bug fix for mushy thermodynamics in brackish seas:

To get the initial ice salinity when new ice forms, a fixed amount of 3 ppt 
was subtracted from the sea surface salinity (SSS), causing negative ice 
salinities when SSS < 3 ppt.  This fix keeps the same 3 ppt reduction in 
salinity for all SSS > 6 ppt, but it makes the initial sea ice salinity go 
to 0 continuously (and differentiably) as SSS approaches 0. We also added 
an arbitrary -0.1 C upper limit to the new ice temperature.

Author: eclare108213

source/ice_therm_itd.F90

@@ -1460,11 +1460,15 @@ subroutine add_new_ice (nx_block,  ny_block,   &
          do ij = 1, icells
             i = indxi(ij)
             j = indxj(ij)
-            Si0new(ij) = sss(i,j) - dSin0_frazil
+            if (sss(i,j) > c2 * dSin0_frazil) then
+               Si0new(ij) = sss(i,j) - dSin0_frazil
+            else
+               Si0new(ij) = sss(i,j)**2 / (c4*dSin0_frazil)
+            endif
             do k = 1, nilyr
                Sprofile(ij,k) = Si0new(ij)
             enddo
-            Ti = liquidus_temperature_mush(Si0new(ij) / phi_init)
+            Ti = min(liquidus_temperature_mush(Si0new(ij)/phi_init), -p1)
             qi0new(ij) = enthalpy_mush(Ti, Si0new(ij))
          enddo ! ij