bug fixes for snow grain radius, brine conservation check (#415)

* bug fixes for snow grain radius, brine conservation check

* bug fix for snow drainage; allow drainage for snow aging when meltpond water is not used

columnphysics/icepack_snow.F90

@@ -16,7 +16,7 @@ module icepack_snow
       use icepack_parameters, only: isnw_T, isnw_Tgrd, isnw_rhos
       use icepack_parameters, only: snowage_rhos, snowage_Tgrd, snowage_T
       use icepack_parameters, only: snowage_tau, snowage_kappa, snowage_drdt0
-      use icepack_parameters, only: snw_aging_table
+      use icepack_parameters, only: snw_aging_table, use_smliq_pnd
 
       use icepack_therm_shared, only: icepack_ice_temperature
       use icepack_therm_shared, only: adjust_enthalpy
@@ -31,7 +31,8 @@ module icepack_snow
 
       real (kind=dbl_kind), parameter, public :: &
          S_r  = 0.033_dbl_kind, & ! irreducible saturation (Anderson 1976)
-         S_wet= 4.22e-5_dbl_kind  ! (um^3/s) wet metamorphism parameters
+         S_wet= 4.22e5_dbl_kind  ! wet metamorphism parameter (um^3/s)
+                                  ! = 1.e18 * 4.22e-13 (Oleson 2010)
 
       real (kind=dbl_kind) :: &
          min_rhos, &   ! snowtable axis data, assumes linear data
@@ -1079,7 +1080,7 @@ subroutine snow_wet_metamorph (dt, dr_wet, rsnw, smice, smliq)
       dr_wet = c0
       fliq = c1
       if (smice + smliq > c0 .and. rsnw > c0) then
-         fliq = min(smliq/(smice + smliq),p1)*c100
+         fliq = min(smliq/(smice + smliq),p1)
          dr_wet = S_wet * fliq**3*dt/(c4*pi*rsnw**2)
       endif
 
@@ -1160,7 +1161,8 @@ subroutine drain_snow (nslyr, vsnon, aicen, &
 
       real (kind=dbl_kind) :: &
          hslyr,  & ! snow layer thickness (m)
-         hsn       ! snow thickness (m)
+         hsn,    & ! snow thickness (m)
+         sliq      ! snow liquid content (kg/m^2)
 
       real (kind=dbl_kind), dimension(nslyr) :: &
          dlin    , & ! liquid mass into the layer from above (kg/m^2)
@@ -1171,27 +1173,29 @@ subroutine drain_snow (nslyr, vsnon, aicen, &
       character (len=*), parameter :: subname='(drain_snow)'
 
       hsn = c0
+      sliq = c0
       if (aicen > c0) hsn = vsnon/aicen
       if (hsn > puny) then
          dlin (:) = c0
          dlout(:) = c0
          hslyr    = hsn / real(nslyr,kind=dbl_kind)
-         meltsliq = c0
          do k = 1, nslyr
             massliq(k) = massliq(k) + dlin(k)   ! add liquid in from layer above
             phi_ice(k) = min(c1, massice(k) / (rhoi    *hslyr))
             phi_liq(k) =         massliq(k) / (rhofresh*hslyr)
-            dlout(k)   = max(c0, (phi_liq(k) - S_r*(c1-phi_ice(k))) / rhofresh*hslyr)
+            dlout(k)   = max(c0, (phi_liq(k) - S_r*(c1-phi_ice(k))) * rhofresh * hslyr)
             massliq(k) = massliq(k) - dlout(k)
             if (k < nslyr) then
                dlin(k+1) = dlout(k)
             else
-               meltsliq = dlout(nslyr) ! this (re)initializes meltsliq
+               sliq = dlout(nslyr) ! this (re)initializes meltsliq
             endif
          enddo
       else
-         meltsliq = meltsliq  ! computed in thickness_changes
+         sliq = meltsliq  ! computed in thickness_changes
       endif
+      meltsliq = meltsliq
+      if (use_smliq_pnd) meltsliq = sliq
 
       end subroutine drain_snow
 

columnphysics/icepack_therm_itd.F90

@@ -236,7 +236,7 @@ subroutine linear_itd (ncat,        hin_max,     &
 
          if (tr_brine) then
             vbrin(n) = vbrin(n) + trcrn(nt_fbri,n) &
-                     * vicen(n)/real(nilyr,kind=dbl_kind)
+                     * vicen(n)
          endif
 
          do k = 1, nilyr
@@ -653,7 +653,7 @@ subroutine linear_itd (ncat,        hin_max,     &
 
          if (tr_brine) then
             vbrin(n) = vbrin(n) + trcrn(nt_fbri,n) &
-                     * vicen(n)/real(nilyr,kind=dbl_kind)
+                     * vicen(n)
          endif
 
          do k = 1, nilyr

columnphysics/icepack_therm_vertical.F90

@@ -2800,7 +2800,7 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr,    &
 
          endif   ! aicen_init
 
-         if (snwgrain .and. use_smliq_pnd) then
+         if (snwgrain) then
             call drain_snow (nslyr = nslyr, &
                              vsnon = vsnon(n), &
                              aicen = aicen(n), &