Advanced snow physics

columnphysics/icepack_algae.F90

@@ -332,6 +332,7 @@ subroutine zbio   (dt,           nblyr,       &
       if (icepack_warnings_aborted(subname)) return
 
       call merge_bgc_fluxes   (dt,           nblyr,      &
+                               nslyr,                    &
                                bio_index,    n_algae,    &
                                nbtrcr,       aicen,      &    
                                vicen,        vsnon,      &

columnphysics/icepack_flux.F90

@@ -58,8 +58,9 @@ subroutine merge_fluxes (aicen,                &
                                fswthru_idr, fswthru_idf,&
                                melttn, meltsn, meltbn, congeln, snoicen, &
                                meltt,  melts,        &
-                               meltb,                &
+                               meltb,  dsnow, dsnown,&
                                congel,  snoice,      &
+                               meltsliq, meltsliqn,  &
                                Uref,     Urefn,      &
                                Qref_iso, Qrefn_iso,  &
                                fiso_ocn, fiso_ocnn,  &
@@ -95,6 +96,8 @@ subroutine merge_fluxes (aicen,                &
           melttn  , & ! top ice melt                    (m)
           meltbn  , & ! bottom ice melt                 (m)
           meltsn  , & ! snow melt                       (m)
+          meltsliqn,& ! mass of snow melt               (kg/m^2)
+          dsnown  , & ! change in snow depth            (m)
           congeln , & ! congelation ice growth          (m)
           snoicen     ! snow-ice growth                 (m)
 
@@ -125,6 +128,8 @@ subroutine merge_fluxes (aicen,                &
           meltt   , & ! top ice melt                    (m)
           meltb   , & ! bottom ice melt                 (m)
           melts   , & ! snow melt                       (m)
+          meltsliq, & ! mass of snow melt               (kg/m^2)
+          dsnow   , & ! change in snow depth            (m)
           congel  , & ! congelation ice growth          (m)
           snoice      ! snow-ice growth                 (m)
 
@@ -212,6 +217,8 @@ subroutine merge_fluxes (aicen,                &
       meltt     = meltt     + melttn    * aicen
       meltb     = meltb     + meltbn    * aicen
       melts     = melts     + meltsn    * aicen
+      meltsliq  = meltsliq  + meltsliqn * aicen
+      dsnow     = dsnow     + dsnown    * aicen
       congel    = congel    + congeln   * aicen
       snoice    = snoice    + snoicen   * aicen
 

columnphysics/icepack_intfc.F90

@@ -82,6 +82,8 @@ module icepack_intfc
       use icepack_wavefracspec, only: icepack_init_wave
       use icepack_wavefracspec, only: icepack_step_wavefracture
 
+      use icepack_snow, only: icepack_step_snow
+
       use icepack_shortwave, only: icepack_prep_radiation
       use icepack_shortwave, only: icepack_step_radiation
 
@@ -114,6 +116,8 @@ module icepack_intfc
       use icepack_mushy_physics , only: icepack_mushy_liquid_fraction
       use icepack_mushy_physics , only: icepack_mushy_temperature_mush
 
+      use icepack_snow, only: icepack_init_snow
+
       use icepack_warnings, only: icepack_warnings_clear
       use icepack_warnings, only: icepack_warnings_print
       use icepack_warnings, only: icepack_warnings_flush

columnphysics/icepack_itd.F90

@@ -28,11 +28,10 @@ module icepack_itd
       use icepack_kinds
       use icepack_parameters, only: c0, c1, c2, c3, c15, c25, c100, p1, p01, p001, p5, puny
       use icepack_parameters, only: Lfresh, rhos, ice_ref_salinity, hs_min, cp_ice, Tocnfrz, rhoi
-      use icepack_parameters, only: rhosi, sk_l, hs_ssl, min_salin
+      use icepack_parameters, only: rhosi, sk_l, hs_ssl, min_salin, rsnw_fall
       use icepack_tracers,    only: nt_Tsfc, nt_qice, nt_qsno, nt_aero, nt_isosno, nt_isoice
       use icepack_tracers,    only: nt_apnd, nt_hpnd, nt_fbri, tr_brine, nt_bgc_S, bio_index
-      use icepack_tracers,    only: n_iso
-      use icepack_tracers,    only: tr_iso
+      use icepack_tracers,    only: n_iso, tr_iso, tr_snow, nt_smice, nt_rsnw, nt_rhos
       use icepack_tracers,    only: icepack_compute_tracers
       use icepack_parameters, only: solve_zsal, skl_bgc, z_tracers
       use icepack_parameters, only: kcatbound, kitd
@@ -1230,11 +1229,15 @@ subroutine zap_small_areas (dt,        ntrcr,        &
 
             if (ntrcr >= 2) then
                do it = 2, ntrcr
-                  if (tr_brine .and. it == nt_fbri) then
-                     trcrn(it,n) = c1
-                  else
-                     trcrn(it,n) = c0
-                  endif
+                  trcrn(it,n) = c0
+               enddo
+            endif
+            if (tr_brine) trcrn(nt_fbri,n) = c1
+            if (tr_snow) then
+               do k = 1, nslyr
+                  trcrn(nt_rhos +k-1,n) = rhos
+                  trcrn(nt_smice+k-1,n) = rhos
+                  trcrn(nt_rsnw +k-1,n) = rsnw_fall
                enddo
             endif
             first_ice(n) = .true.
@@ -1424,14 +1427,14 @@ subroutine zap_snow(dt,         nslyr,    &
       endif ! tr_iso
 
       if (z_tracers) then
-            dvssl  = min(p5*vsnon, hs_ssl*aicen)   !snow surface layer
-            dvint  = vsnon- dvssl                  !snow interior
-
-            do it = 1, nbtrcr
-               xtmp = (trcrn(bio_index(it)+nblyr+1)*dvssl + &
-                       trcrn(bio_index(it)+nblyr+2)*dvint)/dt
-               dflux_bio(it) = dflux_bio(it) + xtmp
-            enddo                 ! it
+         dvssl = min(p5*vsnon/real(nslyr,kind=dbl_kind), hs_ssl*aicen) ! snow surface layer
+         dvint = vsnon - dvssl                                         ! snow interior
+
+         do it = 1, nbtrcr
+            xtmp = (trcrn(bio_index(it)+nblyr+1)*dvssl + &
+                    trcrn(bio_index(it)+nblyr+2)*dvint)/dt
+            dflux_bio(it) = dflux_bio(it) + xtmp
+         enddo                 ! it
 
       endif ! z_tracers
 

columnphysics/icepack_meltpond_cesm.F90

@@ -16,7 +16,7 @@ module icepack_meltpond_cesm
 
       use icepack_kinds
       use icepack_parameters, only: c0, c1, c2, p01, puny
-      use icepack_parameters, only: rhofresh, rhoi, rhos, Timelt
+      use icepack_parameters, only: rhofresh, rhoi, rhos, Timelt, pndaspect, use_smliq_pnd
       use icepack_warnings, only: warnstr, icepack_warnings_add
       use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
 
@@ -32,19 +32,19 @@ module icepack_meltpond_cesm
 !=======================================================================
 
       subroutine compute_ponds_cesm(dt,    hi_min,       &
-                                    pndaspect,           &
                                     rfrac, meltt,        &
                                     melts, frain,        &
-                                    aicen, vicen, &
-                                    Tsfcn, apnd,  hpnd)
+                                    aicen, vicen,        &
+                                    Tsfcn, apnd,  hpnd,  &
+                                    meltsliqn)
 
       real (kind=dbl_kind), intent(in) :: &
          dt,       & ! time step (s)
-         hi_min,   & ! minimum ice thickness allowed for thermo (m)
-         pndaspect   ! ratio of pond depth to pond fraction
+         hi_min      ! minimum ice thickness allowed for thermo (m)
 
       real (kind=dbl_kind), intent(in) :: &
-         rfrac, &    ! water fraction retained for melt ponds
+         meltsliqn, & ! liquid input from snow liquid tracer
+         rfrac, &     ! water fraction retained for melt ponds
          meltt, &
          melts, &
          frain, &
@@ -104,11 +104,18 @@ subroutine compute_ponds_cesm(dt,    hi_min,       &
             !-----------------------------------------------------------
             ! Update pond volume
             !-----------------------------------------------------------
-            volpn = volpn &
-                  + rfrac/rhofresh*(meltt*rhoi &
-                  +                 melts*rhos &
-                  +                 frain*  dt)&
-                  * aicen
+            if (use_smliq_pnd) then
+               volpn = volpn &
+                     + rfrac/rhofresh*(meltt*rhoi &
+                     +                 meltsliqn) &
+                     * aicen
+            else
+               volpn = volpn &
+                     + rfrac/rhofresh*(meltt*rhoi &
+                     +                 melts*rhos &
+                     +                 frain*  dt)&
+                     * aicen
+            endif
 
             !-----------------------------------------------------------
             ! Shrink pond volume under freezing conditions

columnphysics/icepack_meltpond_lvl.F90

@@ -17,7 +17,7 @@ module icepack_meltpond_lvl
       use icepack_kinds
       use icepack_parameters, only: c0, c1, c2, c10, p01, p5, puny
       use icepack_parameters, only: viscosity_dyn, rhoi, rhos, rhow, Timelt, Tffresh, Lfresh
-      use icepack_parameters, only: gravit, depressT, rhofresh, kice
+      use icepack_parameters, only: gravit, depressT, rhofresh, kice, pndaspect, use_smliq_pnd
       use icepack_warnings, only: warnstr, icepack_warnings_add
       use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
 
@@ -35,14 +35,15 @@ module icepack_meltpond_lvl
       subroutine compute_ponds_lvl(dt,     nilyr,        &
                                    ktherm,               &
                                    hi_min, dpscale,      &
-                                   frzpnd, pndaspect,    &
+                                   frzpnd,               &
                                    rfrac,  meltt, melts, &
                                    frain,  Tair,  fsurfn,&
                                    dhs,    ffrac,        &
                                    aicen,  vicen, vsnon, &
                                    qicen,  sicen,        &
                                    Tsfcn,  alvl,         &
-                                   apnd,   hpnd,  ipnd)
+                                   apnd,   hpnd,  ipnd,  &
+                                   meltsliqn)
 
       integer (kind=int_kind), intent(in) :: &
          nilyr, &    ! number of ice layers
@@ -51,8 +52,7 @@ subroutine compute_ponds_lvl(dt,     nilyr,        &
       real (kind=dbl_kind), intent(in) :: &
          dt,       & ! time step (s)  
          hi_min,   & ! minimum ice thickness allowed for thermo (m)
-         dpscale,  & ! alter e-folding time scale for flushing 
-         pndaspect   ! ratio of pond depth to pond fraction
+         dpscale     ! alter e-folding time scale for flushing
 
       character (len=char_len), intent(in) :: &
          frzpnd      ! pond refreezing parameterization
@@ -69,7 +69,8 @@ subroutine compute_ponds_lvl(dt,     nilyr,        &
          fsurfn,&    ! atm-ice surface heat flux  (W/m2)
          aicen, &    ! ice area fraction
          vicen, &    ! ice volume (m)
-         vsnon       ! snow volume (m)
+         vsnon, &    ! snow volume (m)
+         meltsliqn   ! liquid contribution to meltponds in dt (kg/m^2)
 
       real (kind=dbl_kind), &
          intent(inout) :: &
@@ -154,9 +155,14 @@ subroutine compute_ponds_lvl(dt,     nilyr,        &
             ! update pond volume
             !-----------------------------------------------------------
             ! add melt water
-            dvn = rfrac/rhofresh*(meltt*rhoi &
-                +                 melts*rhos &
-                +                 frain*  dt)*aicen
+            if (use_smliq_pnd) then
+               dvn = rfrac/rhofresh*(meltt*rhoi &
+                   +                 meltsliqn)*aicen
+            else
+               dvn = rfrac/rhofresh*(meltt*rhoi &
+                   +                 melts*rhos &
+                   +                 frain*  dt)*aicen
+            endif
 
             ! shrink pond volume under freezing conditions
             if (trim(frzpnd) == 'cesm') then