LCOV - code coverage report
Current view: top level - columnphysics - icepack_meltpond_lvl.F90 (source / functions) Hit Total Coverage
Test: 200804-015008:4c42a82e3d:3:base,travis,quick Lines: 104 107 97.20 %
Date: 2020-08-03 20:10:57 Functions: 2 2 100.00 %

          Line data    Source code
       1             : !=======================================================================
       2             : 
       3             : ! Level-ice meltpond parameterization
       4             : !
       5             : ! This meltpond parameterization was developed for use with the delta-
       6             : ! Eddington radiation scheme, and only affects the radiation budget in
       7             : ! the model.  That is, although the pond volume is tracked, that liquid
       8             : ! water is not used elsewhere in the model for mass budgets or other
       9             : ! physical processes.
      10             : !
      11             : ! authors Elizabeth Hunke (LANL)
      12             : !         David Hebert (NRL Stennis)
      13             : !         Olivier Lecomte (Univ. Louvain)
      14             : 
      15             :       module icepack_meltpond_lvl
      16             : 
      17             :       use icepack_kinds
      18             :       use icepack_parameters, only: c0, c1, c2, c10, p01, p5, puny
      19             :       use icepack_parameters, only: viscosity_dyn, rhoi, rhos, rhow, Timelt, Tffresh, Lfresh
      20             :       use icepack_parameters, only: gravit, depressT, rhofresh, kice
      21             :       use icepack_warnings, only: warnstr, icepack_warnings_add
      22             :       use icepack_warnings, only: icepack_warnings_setabort, icepack_warnings_aborted
      23             : 
      24             :       implicit none
      25             : 
      26             :       private
      27             :       public :: compute_ponds_lvl
      28             : 
      29             : !=======================================================================
      30             : 
      31             :       contains
      32             : 
      33             : !=======================================================================
      34             : 
      35     4305360 :       subroutine compute_ponds_lvl(dt,     nilyr,        &
      36             :                                    ktherm,               &
      37             :                                    hi_min, dpscale,      &
      38           0 :                                    frzpnd, pndaspect,    &
      39             :                                    rfrac,  meltt, melts, &
      40             :                                    frain,  Tair,  fsurfn,&
      41             :                                    dhs,    ffrac,        &
      42             :                                    aicen,  vicen, vsnon, &
      43     4305360 :                                    qicen,  sicen,        &
      44             :                                    Tsfcn,  alvl,         &
      45             :                                    apnd,   hpnd,  ipnd)
      46             : 
      47             :       integer (kind=int_kind), intent(in) :: &
      48             :          nilyr, &    ! number of ice layers
      49             :          ktherm      ! type of thermodynamics (0 0-layer, 1 BL99, 2 mushy)
      50             : 
      51             :       real (kind=dbl_kind), intent(in) :: &
      52             :          dt,       & ! time step (s)  
      53             :          hi_min,   & ! minimum ice thickness allowed for thermo (m)
      54             :          dpscale,  & ! alter e-folding time scale for flushing 
      55             :          pndaspect   ! ratio of pond depth to pond fraction
      56             : 
      57             :       character (len=char_len), intent(in) :: &
      58             :          frzpnd      ! pond refreezing parameterization
      59             : 
      60             :       real (kind=dbl_kind), &
      61             :          intent(in) :: &
      62             :          Tsfcn, &    ! surface temperature (C)
      63             :          alvl,  &    ! fraction of level ice
      64             :          rfrac, &    ! water fraction retained for melt ponds
      65             :          meltt, &    ! top melt rate (m/s)
      66             :          melts, &    ! snow melt rate (m/s)
      67             :          frain, &    ! rainfall rate (kg/m2/s)
      68             :          Tair,  &    ! air temperature (K)
      69             :          fsurfn,&    ! atm-ice surface heat flux  (W/m2)
      70             :          aicen, &    ! ice area fraction
      71             :          vicen, &    ! ice volume (m)
      72             :          vsnon       ! snow volume (m)
      73             : 
      74             :       real (kind=dbl_kind), &
      75             :          intent(inout) :: &
      76             :          apnd, hpnd, ipnd
      77             : 
      78             :       real (kind=dbl_kind), dimension (:), intent(in) :: &
      79             :          qicen, &  ! ice layer enthalpy (J m-3)
      80             :          sicen     ! salinity (ppt)   
      81             : 
      82             :       real (kind=dbl_kind), &
      83             :          intent(in) :: &
      84             :          dhs       ! depth difference for snow on sea ice and pond ice
      85             : 
      86             :       real (kind=dbl_kind), &
      87             :          intent(out) :: &
      88             :          ffrac     ! fraction of fsurfn over pond used to melt ipond
      89             : 
      90             :       ! local temporary variables
      91             : 
      92             :       real (kind=dbl_kind) :: &
      93     1670640 :          volpn     ! pond volume per unit area (m)
      94             : 
      95             :       real (kind=dbl_kind), dimension (nilyr) :: &
      96    19341120 :          Tmlt      ! melting temperature (C)
      97             : 
      98             :       real (kind=dbl_kind) :: &
      99     1670640 :          hi                     , & ! ice thickness (m)
     100     1670640 :          hs                     , & ! snow depth (m)
     101     1670640 :          dTs                    , & ! surface temperature diff for freeze-up (C)
     102     1670640 :          Tp                     , & ! pond freezing temperature (C)
     103     1670640 :          Ts                     , & ! surface air temperature (C)
     104     1670640 :          apondn                 , & ! local pond area 
     105     1670640 :          hpondn                 , & ! local pond depth (m)
     106     1670640 :          dvn                    , & ! change in pond volume (m)
     107     3341280 :          hlid, alid             , & ! refrozen lid thickness, area
     108     1670640 :          dhlid                  , & ! change in refrozen lid thickness
     109     1670640 :          bdt                    , & ! 2 kice dT dt / (rhoi Lfresh)
     110     1670640 :          alvl_tmp               , & ! level ice fraction of ice area
     111     5011920 :          draft, deltah, pressure_head, perm, drain ! for permeability
     112             : 
     113             :       real (kind=dbl_kind), parameter :: &
     114             :          Td       = c2          , & ! temperature difference for freeze-up (C)
     115             :          rexp     = p01             ! pond contraction scaling
     116             : 
     117             :       character(len=*),parameter :: subname='(compute_ponds_lvl)'
     118             : 
     119             :       !-----------------------------------------------------------------
     120             :       ! Initialize 
     121             :       !-----------------------------------------------------------------
     122             : 
     123     4305360 :       volpn = hpnd * aicen * alvl * apnd
     124     4305360 :       ffrac = c0
     125             : 
     126             :       !-----------------------------------------------------------------
     127             :       ! Identify grid cells where ponds can be
     128             :       !-----------------------------------------------------------------
     129             : 
     130     4305360 :       if (aicen*alvl > puny**2) then
     131             :          
     132     2758086 :          hi = vicen/aicen
     133     2758086 :          hs = vsnon/aicen
     134     2758086 :          alvl_tmp = alvl
     135             : 
     136     2758086 :          if (hi < hi_min) then
     137             : 
     138             :             !--------------------------------------------------------------
     139             :             ! Remove ponds on thin ice
     140             :             !--------------------------------------------------------------
     141        8145 :             apondn = c0
     142        8145 :             hpondn = c0
     143        8145 :             volpn  = c0
     144        8145 :             hlid = c0
     145             : 
     146             :          else
     147             : 
     148             :             !-----------------------------------------------------------
     149             :             ! initialize pond area as fraction of ice
     150             :             !-----------------------------------------------------------
     151     2749941 :             apondn = apnd*alvl_tmp
     152             : 
     153             :             !-----------------------------------------------------------
     154             :             ! update pond volume
     155             :             !-----------------------------------------------------------
     156             :             ! add melt water
     157             :             dvn = rfrac/rhofresh*(meltt*rhoi &
     158             :                 +                 melts*rhos &
     159     2749941 :                 +                 frain*  dt)*aicen
     160             : 
     161             :             ! shrink pond volume under freezing conditions
     162     2749941 :             if (trim(frzpnd) == 'cesm') then
     163      261523 :                Tp = Timelt - Td
     164      261523 :                dTs = max(Tp - Tsfcn,c0)
     165      261523 :                dvn = dvn - volpn * (c1 - exp(rexp*dTs/Tp))
     166             : 
     167             :             else 
     168             :                ! trim(frzpnd) == 'hlid' Stefan approximation
     169             :                ! assumes pond is fresh (freezing temperature = 0 C)
     170             :                ! and ice grows from existing pond ice
     171     2488418 :                hlid = ipnd
     172     2488418 :                if (dvn == c0) then ! freeze pond
     173     1463019 :                   Ts = Tair - Tffresh
     174     1463019 :                   if (Ts < c0) then
     175             :                      ! if (Ts < -c2) then ! as in meltpond_cesm
     176     1463019 :                      bdt = -c2*Ts*kice*dt/(rhoi*Lfresh)
     177     1463019 :                      dhlid = p5*sqrt(bdt)                  ! open water freezing
     178     1463019 :                      if (hlid > dhlid) dhlid = p5*bdt/hlid ! existing ice
     179     1463019 :                      dhlid = min(dhlid, hpnd*rhofresh/rhoi)
     180     1463019 :                      hlid = hlid + dhlid
     181             :                   else
     182           0 :                      dhlid = c0 ! to account for surface inversions
     183             :                   endif
     184             :                else ! convert refrozen pond ice back to water
     185     1025399 :                   dhlid = max(fsurfn*dt / (rhoi*Lfresh), c0) ! > 0
     186     1025399 :                   dhlid = -min(dhlid, hlid) ! < 0
     187     1025399 :                   hlid = max(hlid + dhlid, c0)
     188     1025399 :                   if (hs - dhs < puny) then ! pond ice is snow-free
     189       52858 :                      ffrac = c1 ! fraction of fsurfn over pond used to melt ipond
     190       52858 :                      if (fsurfn > puny) &
     191       44163 :                           ffrac = min(-dhlid*rhoi*Lfresh/(dt*fsurfn), c1)
     192             :                   endif
     193             :                endif
     194     2488418 :                alid = apondn * aicen
     195     2488418 :                dvn = dvn - dhlid*alid*rhoi/rhofresh
     196             :             endif
     197             : 
     198     2749941 :             volpn = volpn + dvn
     199             : 
     200             :             !-----------------------------------------------------------
     201             :             ! update pond area and depth
     202             :             !-----------------------------------------------------------
     203     2749941 :             if (volpn <= c0) then
     204     1457026 :                volpn = c0
     205     1457026 :                apondn = c0
     206             :             endif
     207             : 
     208     2749941 :             if (apondn*aicen > puny) then ! existing ponds
     209             :                apondn = max(c0, min(alvl_tmp, &
     210      774853 :                     apondn + 0.5*dvn/(pndaspect*apondn*aicen)))
     211      774853 :                hpondn = c0
     212      774853 :                if (apondn > puny) &
     213      772955 :                     hpondn = volpn/(apondn*aicen)
     214             : 
     215     1975088 :             elseif (alvl_tmp*aicen > c10*puny) then ! new ponds
     216     1961388 :                apondn = min (sqrt(volpn/(pndaspect*aicen)), alvl_tmp)
     217     1961388 :                hpondn = pndaspect * apondn
     218             : 
     219             :             else           ! melt water runs off deformed ice      
     220       13700 :                apondn = c0
     221       13700 :                hpondn = c0
     222             :             endif
     223     2749941 :             apondn = max(apondn, c0)
     224             : 
     225             :             ! limit pond depth to maintain nonnegative freeboard
     226     2749941 :             hpondn = min(hpondn, ((rhow-rhoi)*hi - rhos*hs)/rhofresh)
     227             : 
     228             :             ! fraction of grid cell covered by ponds
     229     2749941 :             apondn = apondn * aicen
     230             : 
     231     2749941 :             volpn = hpondn*apondn
     232     2749941 :             if (volpn <= c0) then
     233     1468440 :                volpn = c0
     234     1468440 :                apondn = c0
     235     1468440 :                hpondn = c0
     236     1468440 :                hlid = c0
     237             :             endif
     238             : 
     239             :             !-----------------------------------------------------------
     240             :             ! drainage due to permeability (flushing)
     241             :             ! setting dpscale = 0 turns this off
     242             :             ! NOTE this uses the initial salinity and melting T profiles
     243             :             !-----------------------------------------------------------
     244             : 
     245     2749941 :             if (ktherm /= 2 .and. hpondn > c0 .and. dpscale > puny) then
     246      361363 :                draft = (rhos*hs + rhoi*hi)/rhow + hpondn
     247      361363 :                deltah = hpondn + hi - draft
     248      361363 :                pressure_head = gravit * rhow * max(deltah, c0)
     249     2890904 :                Tmlt(:) = -sicen(:) * depressT
     250           0 :                call brine_permeability(nilyr, qicen, &
     251      361363 :                     sicen, Tmlt, perm)
     252      361363 :                if (icepack_warnings_aborted(subname)) return
     253      361363 :                drain = perm*pressure_head*dt / (viscosity_dyn*hi) * dpscale
     254      361363 :                deltah = min(drain, hpondn)
     255      361363 :                dvn = -deltah*apondn
     256      361363 :                volpn = volpn + dvn
     257             :                apondn = max(c0, min(apondn &
     258      361363 :                     + 0.5*dvn/(pndaspect*apondn), alvl_tmp*aicen))
     259      361363 :                hpondn = c0
     260      361363 :                if (apondn > puny) hpondn = volpn/apondn
     261             :             endif
     262             : 
     263             :          endif
     264             : 
     265             :          !-----------------------------------------------------------
     266             :          ! Reload tracer array
     267             :          !-----------------------------------------------------------
     268             : 
     269     2758086 :          hpnd = hpondn
     270     2758086 :          apnd = apondn / (aicen*alvl_tmp)
     271     2758086 :          if (trim(frzpnd) == 'hlid') ipnd = hlid
     272             : 
     273             :       endif
     274             : 
     275             :       end subroutine compute_ponds_lvl
     276             : 
     277             : !=======================================================================
     278             : 
     279             : ! determine the liquid fraction of brine in the ice and the permeability
     280             : 
     281      361363 :       subroutine brine_permeability(nilyr, qicen, salin, Tmlt, perm)
     282             : 
     283             :       use icepack_therm_shared, only: calculate_Tin_from_qin
     284             : 
     285             :       integer (kind=int_kind), intent(in) :: &
     286             :          nilyr     ! number of ice layers
     287             : 
     288             :       real (kind=dbl_kind), dimension(:), intent(in) :: &
     289             :          qicen, &  ! enthalpy for each ice layer (J m-3)
     290             :          salin, &  ! salinity (ppt)   
     291             :          Tmlt      ! melting temperature (C)
     292             :     
     293             :       real (kind=dbl_kind), intent(out) :: &
     294             :          perm      ! permeability (m^2)
     295             : 
     296             :       ! local variables
     297             : 
     298             :       real (kind=dbl_kind) ::   &
     299      129283 :          Sbr       ! brine salinity
     300             : 
     301             :       real (kind=dbl_kind), dimension(nilyr) ::   &
     302     1498424 :          Tin, &    ! ice temperature (C)
     303     1627707 :          phi       ! liquid fraction
     304             : 
     305             :       integer (kind=int_kind) :: k
     306             :     
     307             :       character(len=*),parameter :: subname='(brine_permeability)'
     308             : 
     309             :       !-----------------------------------------------------------------
     310             :       ! Compute ice temperatures from enthalpies using quadratic formula
     311             :       !-----------------------------------------------------------------
     312             : 
     313     2890904 :       do k = 1,nilyr
     314     2890904 :          Tin(k) = calculate_Tin_from_qin(qicen(k),Tmlt(k))
     315             :       enddo
     316             : 
     317             :       !-----------------------------------------------------------------
     318             :       ! brine salinity and liquid fraction
     319             :       !-----------------------------------------------------------------
     320             : 
     321     2890904 :       do k = 1,nilyr
     322     2529541 :          Sbr = c1/(1.e-3_dbl_kind - depressT/Tin(k)) ! Notz thesis eq 3.6
     323     2529541 :          phi(k) = salin(k)/Sbr ! liquid fraction
     324     2890904 :          if (phi(k) < 0.05) phi(k) = c0 ! impermeable
     325             :       enddo
     326             : 
     327             :       !-----------------------------------------------------------------
     328             :       ! permeability
     329             :       !-----------------------------------------------------------------
     330             : 
     331     3252267 :       perm = 3.0e-8_dbl_kind * (minval(phi))**3
     332             :     
     333      361363 :       end subroutine brine_permeability
     334             :   
     335             : !=======================================================================
     336             : 
     337             :       end module icepack_meltpond_lvl
     338             : 
     339             : !=======================================================================

Generated by: LCOV version 1.14-6-g40580cd