!======================================================================= ! Level-ice meltpond parameterization ! ! This meltpond parameterization was developed for use with the delta- ! Eddington radiation scheme, and only affects the radiation budget in ! the model. That is, although the pond volume is tracked, that liquid ! water is not used elsewhere in the model for mass budgets or other ! physical processes. ! ! authors Elizabeth Hunke (LANL) ! David Hebert (NRL Stennis) ! Olivier Lecomte (Univ. Louvain) 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_warnings, only: warnstr, icepack_warnings_add use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted implicit none private public :: compute_ponds_lvl !======================================================================= contains !======================================================================= subroutine compute_ponds_lvl(dt, nilyr, & ktherm, & hi_min, dpscale, & frzpnd, pndaspect, & rfrac, meltt, melts, & frain, Tair, fsurfn,& dhs, ffrac, & aicen, vicen, vsnon, & qicen, sicen, & Tsfcn, alvl, & apnd, hpnd, ipnd) integer (kind=int_kind), intent(in) :: & nilyr, & ! number of ice layers ktherm ! type of thermodynamics (0 0-layer, 1 BL99, 2 mushy) 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 character (len=char_len), intent(in) :: & frzpnd ! pond refreezing parameterization real (kind=dbl_kind), & intent(in) :: & Tsfcn, & ! surface temperature (C) alvl, & ! fraction of level ice rfrac, & ! water fraction retained for melt ponds meltt, & ! top melt rate (m/s) melts, & ! snow melt rate (m/s) frain, & ! rainfall rate (kg/m2/s) Tair, & ! air temperature (K) fsurfn,& ! atm-ice surface heat flux (W/m2) aicen, & ! ice area fraction vicen, & ! ice volume (m) vsnon ! snow volume (m) real (kind=dbl_kind), & intent(inout) :: & apnd, hpnd, ipnd real (kind=dbl_kind), dimension (:), intent(in) :: & qicen, & ! ice layer enthalpy (J m-3) sicen ! salinity (ppt) real (kind=dbl_kind), & intent(in) :: & dhs ! depth difference for snow on sea ice and pond ice real (kind=dbl_kind), & intent(out) :: & ffrac ! fraction of fsurfn over pond used to melt ipond ! local temporary variables real (kind=dbl_kind) :: & volpn ! pond volume per unit area (m) real (kind=dbl_kind), dimension (nilyr) :: & Tmlt ! melting temperature (C) real (kind=dbl_kind) :: & hi , & ! ice thickness (m) hs , & ! snow depth (m) dTs , & ! surface temperature diff for freeze-up (C) Tp , & ! pond freezing temperature (C) Ts , & ! surface air temperature (C) apondn , & ! local pond area hpondn , & ! local pond depth (m) dvn , & ! change in pond volume (m) hlid, alid , & ! refrozen lid thickness, area dhlid , & ! change in refrozen lid thickness bdt , & ! 2 kice dT dt / (rhoi Lfresh) alvl_tmp , & ! level ice fraction of ice area draft, deltah, pressure_head, perm, drain ! for permeability real (kind=dbl_kind), parameter :: & Td = c2 , & ! temperature difference for freeze-up (C) rexp = p01 ! pond contraction scaling character(len=*),parameter :: subname='(compute_ponds_lvl)' !----------------------------------------------------------------- ! Initialize !----------------------------------------------------------------- volpn = hpnd * aicen * alvl * apnd ffrac = c0 !----------------------------------------------------------------- ! Identify grid cells where ponds can be !----------------------------------------------------------------- if (aicen*alvl > puny**2) then hi = vicen/aicen hs = vsnon/aicen alvl_tmp = alvl if (hi < hi_min) then !-------------------------------------------------------------- ! Remove ponds on thin ice !-------------------------------------------------------------- apondn = c0 hpondn = c0 volpn = c0 hlid = c0 else !----------------------------------------------------------- ! initialize pond area as fraction of ice !----------------------------------------------------------- apondn = apnd*alvl_tmp !----------------------------------------------------------- ! update pond volume !----------------------------------------------------------- ! add melt water dvn = rfrac/rhofresh*(meltt*rhoi & + melts*rhos & + frain* dt)*aicen ! shrink pond volume under freezing conditions if (trim(frzpnd) == 'cesm') then Tp = Timelt - Td dTs = max(Tp - Tsfcn,c0) dvn = dvn - volpn * (c1 - exp(rexp*dTs/Tp)) else ! trim(frzpnd) == 'hlid' Stefan approximation ! assumes pond is fresh (freezing temperature = 0 C) ! and ice grows from existing pond ice hlid = ipnd if (dvn == c0) then ! freeze pond Ts = Tair - Tffresh if (Ts < c0) then ! if (Ts < -c2) then ! as in meltpond_cesm bdt = -c2*Ts*kice*dt/(rhoi*Lfresh) dhlid = p5*sqrt(bdt) ! open water freezing if (hlid > dhlid) dhlid = p5*bdt/hlid ! existing ice dhlid = min(dhlid, hpnd*rhofresh/rhoi) hlid = hlid + dhlid else dhlid = c0 ! to account for surface inversions endif else ! convert refrozen pond ice back to water dhlid = max(fsurfn*dt / (rhoi*Lfresh), c0) ! > 0 dhlid = -min(dhlid, hlid) ! < 0 hlid = max(hlid + dhlid, c0) if (hs - dhs < puny) then ! pond ice is snow-free ffrac = c1 ! fraction of fsurfn over pond used to melt ipond if (fsurfn > puny) & ffrac = min(-dhlid*rhoi*Lfresh/(dt*fsurfn), c1) endif endif alid = apondn * aicen dvn = dvn - dhlid*alid*rhoi/rhofresh endif volpn = volpn + dvn !----------------------------------------------------------- ! update pond area and depth !----------------------------------------------------------- if (volpn <= c0) then volpn = c0 apondn = c0 endif if (apondn*aicen > puny) then ! existing ponds apondn = max(c0, min(alvl_tmp, & apondn + 0.5*dvn/(pndaspect*apondn*aicen))) hpondn = c0 if (apondn > puny) & hpondn = volpn/(apondn*aicen) elseif (alvl_tmp*aicen > c10*puny) then ! new ponds apondn = min (sqrt(volpn/(pndaspect*aicen)), alvl_tmp) hpondn = pndaspect * apondn else ! melt water runs off deformed ice apondn = c0 hpondn = c0 endif apondn = max(apondn, c0) ! limit pond depth to maintain nonnegative freeboard hpondn = min(hpondn, ((rhow-rhoi)*hi - rhos*hs)/rhofresh) ! fraction of grid cell covered by ponds apondn = apondn * aicen volpn = hpondn*apondn if (volpn <= c0) then volpn = c0 apondn = c0 hpondn = c0 hlid = c0 endif !----------------------------------------------------------- ! drainage due to permeability (flushing) ! setting dpscale = 0 turns this off ! NOTE this uses the initial salinity and melting T profiles !----------------------------------------------------------- if (ktherm /= 2 .and. hpondn > c0 .and. dpscale > puny) then draft = (rhos*hs + rhoi*hi)/rhow + hpondn deltah = hpondn + hi - draft pressure_head = gravit * rhow * max(deltah, c0) Tmlt(:) = -sicen(:) * depressT call brine_permeability(nilyr, qicen, & sicen, Tmlt, perm) if (icepack_warnings_aborted(subname)) return drain = perm*pressure_head*dt / (viscosity_dyn*hi) * dpscale deltah = min(drain, hpondn) dvn = -deltah*apondn volpn = volpn + dvn apondn = max(c0, min(apondn & + 0.5*dvn/(pndaspect*apondn), alvl_tmp*aicen)) hpondn = c0 if (apondn > puny) hpondn = volpn/apondn endif endif !----------------------------------------------------------- ! Reload tracer array !----------------------------------------------------------- hpnd = hpondn apnd = apondn / (aicen*alvl_tmp) if (trim(frzpnd) == 'hlid') ipnd = hlid endif end subroutine compute_ponds_lvl !======================================================================= ! determine the liquid fraction of brine in the ice and the permeability subroutine brine_permeability(nilyr, qicen, salin, Tmlt, perm) use icepack_therm_shared, only: calculate_Tin_from_qin integer (kind=int_kind), intent(in) :: & nilyr ! number of ice layers real (kind=dbl_kind), dimension(:), intent(in) :: & qicen, & ! enthalpy for each ice layer (J m-3) salin, & ! salinity (ppt) Tmlt ! melting temperature (C) real (kind=dbl_kind), intent(out) :: & perm ! permeability (m^2) ! local variables real (kind=dbl_kind) :: & Sbr ! brine salinity real (kind=dbl_kind), dimension(nilyr) :: & Tin, & ! ice temperature (C) phi ! liquid fraction integer (kind=int_kind) :: k character(len=*),parameter :: subname='(brine_permeability)' !----------------------------------------------------------------- ! Compute ice temperatures from enthalpies using quadratic formula !----------------------------------------------------------------- do k = 1,nilyr Tin(k) = calculate_Tin_from_qin(qicen(k),Tmlt(k)) enddo !----------------------------------------------------------------- ! brine salinity and liquid fraction !----------------------------------------------------------------- do k = 1,nilyr Sbr = c1/(1.e-3_dbl_kind - depressT/Tin(k)) ! Notz thesis eq 3.6 phi(k) = salin(k)/Sbr ! liquid fraction if (phi(k) < 0.05) phi(k) = c0 ! impermeable enddo !----------------------------------------------------------------- ! permeability !----------------------------------------------------------------- perm = 3.0e-8_dbl_kind * (minval(phi))**3 end subroutine brine_permeability !======================================================================= end module icepack_meltpond_lvl !=======================================================================