A new sea level pond scheme.

columnphysics/icepack_flux.F90

@@ -12,7 +12,7 @@ module icepack_flux
       use icepack_parameters, only: c1, emissivity, snwgrain
       use icepack_warnings, only: warnstr, icepack_warnings_add
       use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
-      use icepack_tracers, only: tr_iso
+      use icepack_tracers, only: tr_iso, tr_pond
 
       implicit none
       private
@@ -64,7 +64,12 @@ subroutine merge_fluxes (aicen,                &
                                Uref,     Urefn,      &
                                Qref_iso, Qrefn_iso,  &
                                fiso_ocn, fiso_ocnn,  &
-                               fiso_evap, fiso_evapn)
+                               fiso_evap, fiso_evapn,&
+                               flpnd,  flpndn,       &
+                               expnd,  expndn,       &
+                               frpnd,  frpndn,       &
+                               rfpnd,  rfpndn,       &
+                               ilpnd,  ilpndn)
 
       ! single category fluxes
       real (kind=dbl_kind), intent(in) :: &
@@ -98,6 +103,11 @@ subroutine merge_fluxes (aicen,                &
           dsnown  , & ! change in snow depth            (m)
           congeln , & ! congelation ice growth          (m)
           snoicen , & ! snow-ice growth                 (m)
+          flpndn  , & ! pond flushing rate due to ice permeability (m/step)
+          expndn  , & ! exponential pond drainage rate (m/step)
+          frpndn  , & ! pond drainage rate due to freeboard constraint (m/step)
+          rfpndn  , & ! runoff rate due to rfrac (m/step)
+          ilpndn  , & ! pond loss/gain due to ice lid (m/step)
           fswthrun_vdr, & ! vis dir sw radiation through ice bot    (W/m**2)
           fswthrun_vdf, & ! vis dif sw radiation through ice bot    (W/m**2)
           fswthrun_idr, & ! nir dir sw radiation through ice bot    (W/m**2)
@@ -131,10 +141,15 @@ subroutine merge_fluxes (aicen,                &
           meltsliq, & ! mass of snow melt               (kg/m^2)
           congel  , & ! congelation ice growth          (m)
           snoice  , & ! snow-ice growth                 (m)
-          fswthru_vdr, & ! vis dir sw radiation through ice bot    (W/m**2)
-          fswthru_vdf, & ! vis dif sw radiation through ice bot    (W/m**2)
-          fswthru_idr, & ! nir dir sw radiation through ice bot    (W/m**2)
-          fswthru_idf, & ! nir dif sw radiation through ice bot    (W/m**2)
+          flpnd   , & ! pond flushing rate due to ice permeability (m/step)
+          expnd   , & ! exponential pond drainage rate (m/step)
+          frpnd   , & ! pond drainage rate due to freeboard constraint (m/step)
+          rfpnd   , & ! runoff rate due to rfrac (m/step)
+          ilpnd   , & ! pond loss/gain (+/-) to ice lid freezing/melting (m/step)
+          fswthru_vdr , & ! vis dir sw radiation through ice bot    (W/m**2)
+          fswthru_vdf , & ! vis dif sw radiation through ice bot    (W/m**2)
+          fswthru_idr , & ! nir dir sw radiation through ice bot    (W/m**2)
+          fswthru_idf , & ! nir dif sw radiation through ice bot    (W/m**2)
           dsnow,    & ! change in snow depth            (m)
           Uref        ! air speed reference level       (m/s)
 
@@ -247,6 +262,19 @@ subroutine merge_fluxes (aicen,                &
          congel    = congel    + congeln   * aicen
       if (present(snoicen) .and. present(snoice)) &
          snoice    = snoice    + snoicen   * aicen
+      ! Meltwater fluxes
+      if (tr_pond) then
+         if (present(flpndn) .and. present(flpnd)) &
+            flpnd     = flpnd     + flpndn    * aicen
+         if (present(expndn) .and. present(expnd)) &
+            expnd     = expnd     + expndn    * aicen
+         if (present(frpndn) .and. present(frpnd)) &
+            frpnd     = frpnd     + frpndn    * aicen
+         if (present(rfpndn) .and. present(rfpnd)) &
+            rfpnd     = rfpnd     + rfpndn    * aicen
+         if (present(ilpndn) .and. present(ilpnd)) &
+            ilpnd     = ilpnd     + ilpndn    * aicen
+      endif
 
       end subroutine merge_fluxes
 

columnphysics/icepack_intfc.F90

@@ -125,6 +125,8 @@ module icepack_intfc
       use icepack_mushy_physics , only: icepack_mushy_liquid_fraction
       use icepack_mushy_physics , only: icepack_mushy_temperature_mush
 
+      use icepack_meltpond_sealvl, only: icepack_init_sealvlpnd
+
       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

@@ -32,6 +32,7 @@ module icepack_itd
       use icepack_tracers,    only: ncat, nilyr, nslyr, nblyr, ntrcr, nbtrcr, n_aero
       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, bio_index
+      use icepack_tracers,    only: tr_pond, tr_pond_lvl, nt_alvl
       use icepack_tracers,    only: n_iso, tr_iso, nt_smice, nt_rsnw, nt_rhos, nt_sice
       use icepack_tracers,    only: icepack_compute_tracers
       use icepack_parameters, only: skl_bgc, z_tracers, hi_min
@@ -295,28 +296,36 @@ end subroutine rebin
 
       subroutine reduce_area (hin_max,            &
                               aicen,     vicen,   &
-                              aicen_init,vicen_init)
+                              aicen_init,vicen_init, &
+                              mipnd, trcrn)
 
       real (kind=dbl_kind), intent(in) :: &
          hin_max       ! lowest category boundary
 
       real (kind=dbl_kind), intent(inout) :: &
          aicen     , & ! concentration of ice
-         vicen         ! volume per unit area of ice          (m)
+         vicen     , & ! volume per unit area of ice          (m)
+         mipnd         ! pond 'drainage' due to ice melting (m / step)
 
       real (kind=dbl_kind), intent(in) :: &
          aicen_init, & ! old ice area for category 1 (m)
          vicen_init    ! old ice volume for category 1 (m)
 
+      real (kind=dbl_kind), dimension (:,:), intent(in) :: &
+         trcrn     ! ice tracers
+
       ! local variables
 
       real (kind=dbl_kind) :: &
          hi0       , & ! initial hi
          hi1       , & ! current hi
-         dhi           ! hi1 - hi0
+         dhi       , & ! hi1 - hi0
+         da            ! total change in area complete within this subroutine
 
       character(len=*),parameter :: subname='(reduce_area)'
 
+            da = aicen ! store value of aicen at start of the subroutine
+
             hi0 = c0
             if (aicen_init > c0) &
                 hi0 = vicen_init / aicen_init
@@ -339,6 +348,15 @@ subroutine reduce_area (hin_max,            &
                endif
             endif
 
+            da = da - aicen ! -1*change in fractional area over the subroutine
+            if (tr_pond) then
+               if (tr_pond_lvl) then
+                  mipnd = da*trcrn(nt_apnd,1)*trcrn(nt_hpnd,1)*trcrn(nt_alvl,1)
+               else
+                  mipnd = da*trcrn(nt_apnd,1)*trcrn(nt_hpnd,1)
+               endif
+            endif
+
       end subroutine reduce_area
 
 !=======================================================================

columnphysics/icepack_mechred.F90

@@ -42,7 +42,8 @@ module icepack_mechred
       use icepack_parameters, only: kstrength, krdg_partic, krdg_redist, mu_rdg
       use icepack_parameters, only: conserv_check, z_tracers
       use icepack_tracers, only: ncat, nilyr, nslyr, nblyr, n_aero
-      use icepack_tracers, only: tr_pond_topo, tr_aero, tr_iso, tr_brine, ntrcr, nbtrcr
+      use icepack_tracers, only: tr_aero, tr_iso, tr_brine, ntrcr, nbtrcr
+      use icepack_tracers, only: tr_pond_lvl, tr_pond_topo, tr_pond_sealvl
       use icepack_tracers, only: nt_qice, nt_qsno, nt_fbri, nt_sice
       use icepack_tracers, only: nt_alvl, nt_vlvl, nt_aero, nt_isosno, nt_isoice
       use icepack_tracers, only: nt_apnd, nt_hpnd
@@ -108,7 +109,7 @@ subroutine ridge_ice (dt,          ndtd,       &
                             dardg1ndt,   dardg2ndt,  &
                             dvirdgndt,   Tf,         &
                             araftn,      vraftn,     &
-                            closing )
+                            closing,     rdpnd)
 
       integer (kind=int_kind), intent(in) :: &
          ndtd       ! number of dynamics subcycles
@@ -165,7 +166,8 @@ subroutine ridge_ice (dt,          ndtd,       &
          closing    , & ! rate of closing due to divergence/shear (1/s)
          fpond      , & ! fresh water flux to ponds (kg/m^2/s)
          fresh      , & ! fresh water flux to ocean (kg/m^2/s)
-         fhocn          ! net heat flux to ocean (W/m^2)
+         fhocn      , & ! net heat flux to ocean (W/m^2)
+         rdpnd          ! pond drainage due to ridging (m w.e. avg. over cell)
 
       real (kind=dbl_kind), dimension(:), intent(inout), optional :: &
          dardg1ndt  , & ! rate of fractional area loss by ridging ice (1/s)
@@ -202,7 +204,8 @@ subroutine ridge_ice (dt,          ndtd,       &
          aksum      , & ! ratio of area removed to area ridged
          msnow_mlt  , & ! mass of snow added to ocean (kg m-2)
          esnow_mlt  , & ! energy needed to melt snow in ocean (J m-2)
-         mpond      , & ! mass of pond added to ocean (kg m-2)
+         mpond      , & ! thickness of pond water (avg. over grid cell) lost
+                        ! due to ridging (m)
          closing_net, & ! net rate at which area is removed    (1/s)
                         ! (ridging ice area - area of new ridges) / dt
          divu_adv   , & ! divu as implied by transport scheme  (1/s)
@@ -605,6 +608,7 @@ subroutine ridge_ice (dt,          ndtd,       &
       if (present(fpond)) then
          fpond = fpond - mpond ! units change later
       endif
+      if (present(rdpnd)) rdpnd = mpond
 
       !-----------------------------------------------------------------
       ! Check for fractional ice area > 1.
@@ -1143,7 +1147,8 @@ subroutine ridge_shift (dt,          hin_max,         &
       real (kind=dbl_kind), intent(inout) :: &
          msnow_mlt  , & ! mass of snow added to ocean (kg m-2)
          esnow_mlt  , & ! energy needed to melt snow in ocean (J m-2)
-         mpond          ! mass of pond added to ocean (kg m-2)
+         mpond          ! thickness of pond water, averaged over entire grid
+                        ! cell area, lost to ridging (m)
 
       real (kind=dbl_kind), dimension(:), intent(inout) :: &
          maero          ! aerosol mass added to ocean (kg m-2)
@@ -1389,11 +1394,26 @@ subroutine ridge_shift (dt,          hin_max,         &
                enddo
             endif
 
-            if (tr_pond_topo) then
+            if (z_tracers .and. nbtrcr > 0) then
+               dzssl = p5/real(nslyr,kind=dbl_kind)
+               dzint = c1-dzssl
+               do it = 1, nbtrcr
+                  mbio(it) = mbio(it) + vsrdgn*(c1-fsnowrdg) &
+                           * (trcrn(bio_index(it) + nblyr + 1,n) * dzssl &
+                            + trcrn(bio_index(it) + nblyr + 2,n) * dzint)
+               enddo
+            endif
+
+            if (tr_pond_topo .or. tr_pond_sealvl) then
                mpond = mpond + ardg1n * trcrn(nt_apnd,n) &
                                       * trcrn(nt_hpnd,n)
             endif
 
+            if (tr_pond_lvl) then
+               mpond = mpond + ardg1n * trcrn(nt_apnd,n) &
+                                      * trcrn(nt_hpnd,n) * trcrn(nt_alvl,n)
+            endif
+
       !-----------------------------------------------------------------
       ! Compute quantities used to apportion ice among categories
       ! in the nr loop below
@@ -1741,7 +1761,7 @@ subroutine icepack_step_ridge(dt,           ndtd,          &
                                     aice,         fsalt,         &
                                     first_ice,    fzsal,         &
                                     flux_bio,     closing,       &
-                                    Tf,                          &
+                                    Tf,           rdpnd,         &
                                     docleanup,    dorebin)
 
       real (kind=dbl_kind), intent(in) :: &
@@ -1817,6 +1837,9 @@ subroutine icepack_step_ridge(dt,           ndtd,          &
       logical (kind=log_kind), dimension(:), intent(inout) :: &
          first_ice    ! true until ice forms
 
+      real (kind=dbl_kind), intent(inout), optional :: &
+         rdpnd        ! pond drainage due to ridging (m w.e. avg. over cell)
+
      logical (kind=log_kind), intent(in), optional ::   &
          docleanup, & ! if false, do not call cleanup_itd (default true)
          dorebin      ! if false, do not call rebin in cleanup_itd (default true)
@@ -1893,7 +1916,7 @@ subroutine icepack_step_ridge(dt,           ndtd,          &
                       dardg1ndt,    dardg2ndt,      &
                       dvirdgndt,    Tf,             &
                       araftn,       vraftn,         &
-                      closing )
+                      closing,      rdpnd )
       if (icepack_warnings_aborted(subname)) return
 
       !-----------------------------------------------------------------

columnphysics/icepack_meltpond_lvl.F90

@@ -42,7 +42,9 @@ subroutine compute_ponds_lvl(dt,                   &
                                    qicen,  sicen,        &
                                    Tsfcn,  alvl,         &
                                    apnd,   hpnd,  ipnd,  &
-                                   meltsliqn)
+                                   meltsliqn, frpndn,    &
+                                   rfpndn, ilpndn,       &
+                                   flpndn)
 
       real (kind=dbl_kind), intent(in) :: &
          dt          ! time step (s)
@@ -62,7 +64,11 @@ subroutine compute_ponds_lvl(dt,                   &
          meltsliqn   ! liquid contribution to meltponds in dt (kg/m^2)
 
       real (kind=dbl_kind), intent(inout) :: &
-         apnd, hpnd, ipnd
+         apnd, hpnd, ipnd, &
+         frpndn, &   ! pond drainage rate due to freeboard constraint (m/step)
+         rfpndn, &   ! runoff rate due to rfrac (m/step)
+         ilpndn, &   ! pond loss/gain due to ice lid (m/step)
+         flpndn      ! pond flushing rate due to ice permeability (m/s)
 
       real (kind=dbl_kind), dimension (:), intent(in) :: &
          qicen, &  ! ice layer enthalpy (J m-3)
@@ -77,7 +83,9 @@ subroutine compute_ponds_lvl(dt,                   &
       ! local temporary variables
 
       real (kind=dbl_kind) :: &
-         volpn     ! pond volume per unit area (m)
+         volpn, &     ! pond volume per unit area (m)
+         hpond_tmp, & ! local variable for hpond before flushing
+         dvn_temp     ! local variable for change in volume due to rfrac
 
       real (kind=dbl_kind), dimension (nilyr) :: &
          Tmlt      ! melting temperature (C)
@@ -149,6 +157,12 @@ subroutine compute_ponds_lvl(dt,                   &
                    +                 melts*rhos &
                    +                 frain*  dt)*aicen
             endif
+            ! Track lost meltwater dvn is volume of meltwater (m3/m2) captured
+            ! over entire grid cell area. Multiply by (1-rfrac)/rfrac to get
+            ! loss over entire area. And divide by aicen to get loss per unit
+            ! category area (for consistency with melttn, frpndn, etc)
+            rfpndn = dvn * (c1-rfrac) / (rfrac * aicen)
+            dvn_temp = dvn
 
             ! shrink pond volume under freezing conditions
             if (trim(frzpnd) == 'cesm') then
@@ -188,6 +202,9 @@ subroutine compute_ponds_lvl(dt,                   &
             endif
 
             volpn = volpn + dvn
+            ! Track lost/gained meltwater per unit category area from pond 
+            ! lid freezing/melting. Note sign flip relative to dvn convention
+            ilpndn = (dvn_temp - dvn) / aicen
 
             !-----------------------------------------------------------
             ! update pond area and depth
@@ -206,21 +223,31 @@ subroutine compute_ponds_lvl(dt,                   &
 
             elseif (alvl_tmp*aicen > c10*puny) then ! new ponds
                apondn = min (sqrt(volpn/(pndaspect*aicen)), alvl_tmp)
-               hpondn = pndaspect * apondn
+               hpondn = pndaspect * apondn ! Possible loss of meltwater if apondn == alvl_tmp
 
             else           ! melt water runs off deformed ice
                apondn = c0
-               hpondn = c0
+               hpondn = c0 ! Loss of meltwater for very deformed ice
             endif
             apondn = max(apondn, c0)
 
             ! limit pond depth to maintain nonnegative freeboard
+            hpond_tmp = hpondn
             hpondn = min(hpondn, ((rhow-rhoi)*hi - rhos*hs)/rhofresh)
+            ! The way apondn is used is very confusing but at this point
+            ! apondn is the fraction of the entire category (level + deformed)
+            ! with ponds on it. Thus, multiplying the change in hpondn (i.e.,
+            ! the meltwater lost from the ponded area) by apondn here yields
+            ! the meltwater height lost averaged over the category area
+            ! analogous to how melttn is defined.
+            frpndn = (hpond_tmp - hpondn) * apondn
 
             ! fraction of grid cell covered by ponds
             apondn = apondn * aicen
 
             volpn = hpondn*apondn
+            ! note, this implies that if ponds fully drain or freeze their
+            ! depressions cease to exist and the lid ice also ceases to exist
             if (volpn <= c0) then
                volpn = c0
                apondn = c0
@@ -250,6 +277,7 @@ subroutine compute_ponds_lvl(dt,                   &
                     + 0.5*dvn/(pndaspect*apondn), alvl_tmp*aicen))
                hpondn = c0
                if (apondn > puny) hpondn = volpn/apondn
+               flpndn = -dvn/aicen
             endif
 
          endif