Merge pull request #33 from apcraig/upd230930

Update Icepack Consortium main #390fb5518326217c2501dbf5, Sept 30, 2023

columnphysics/icepack_itd.F90

@@ -1495,7 +1495,7 @@ subroutine zap_snow_temperature(dt,         ncat,     &
       real (kind=dbl_kind) :: &
          rnslyr   , & ! real(nslyr)
          hsn      , & ! snow thickness (m)
-         zqsn     , & ! snow layer enthalpy (J m-2)
+         zqsn     , & ! snow layer enthalpy (J m-3)
          zTsn     , & ! snow layer temperature (C)
          Tmax         ! maximum allowed snow temperature
 

columnphysics/icepack_parameters.F90

@@ -108,6 +108,7 @@ module icepack_parameters
                                         ! freshwater value needed for enthalpy
          depressT  = 0.054_dbl_kind   ,&! Tf:brine salinity ratio (C/ppt)
          viscosity_dyn = 1.79e-3_dbl_kind, & ! dynamic viscosity of brine (kg/m/s)
+         tscale_pnd_drain = c10       ,&! mushy macroscopic drainage timescale (days)
          Tocnfrz   = -1.8_dbl_kind    ,&! freezing temp of seawater (C),
                                         ! used as Tsfcn for open water
          Tffresh   = 273.15_dbl_kind  ,&! freezing temp of fresh ice (K)
@@ -331,7 +332,7 @@ module icepack_parameters
          frzpnd    = 'cesm'           ! pond refreezing parameterization
 
       real (kind=dbl_kind), public :: &
-         dpscale   = c1, &            ! alter e-folding time scale for flushing
+         dpscale   = 0.001_dbl_kind,& ! alter e-folding time scale for flushing (ktherm=1)
          rfracmin  = 0.15_dbl_kind, & ! minimum retained fraction of meltwater
          rfracmax  = 0.85_dbl_kind, & ! maximum retained fraction of meltwater
          pndaspect = 0.8_dbl_kind, &  ! ratio of pond depth to area fraction
@@ -440,7 +441,7 @@ subroutine icepack_init_parameters(   &
          rhos_in, rhoi_in, rhow_in, cp_air_in, emissivity_in, &
          cp_ice_in, cp_ocn_in, hfrazilmin_in, floediam_in, &
          depressT_in, dragio_in, thickness_ocn_layer1_in, iceruf_ocn_in, &
-         albocn_in, gravit_in, viscosity_dyn_in, &
+         albocn_in, gravit_in, viscosity_dyn_in, tscale_pnd_drain_in, &
          Tocnfrz_in, rhofresh_in, zvir_in, vonkar_in, cp_wv_in, &
          stefan_boltzmann_in, ice_ref_salinity_in, &
          Tffresh_in, Lsub_in, Lvap_in, Timelt_in, Tsmelt_in, &
@@ -523,6 +524,7 @@ subroutine icepack_init_parameters(   &
          cp_ocn_in,     & ! specific heat of ocn    (J/kg/K)
          depressT_in,   & ! Tf:brine salinity ratio (C/ppt)
          viscosity_dyn_in, & ! dynamic viscosity of brine (kg/m/s)
+         tscale_pnd_drain_in,&! mushy macroscopic drainage timescale (days)
          Tocnfrz_in,    & ! freezing temp of seawater (C)
          Tffresh_in,    & ! freezing temp of fresh ice (K)
          Lsub_in,       & ! latent heat, sublimation freshwater (J/kg)
@@ -866,6 +868,7 @@ subroutine icepack_init_parameters(   &
       if (present(albocn_in)            ) albocn           = albocn_in
       if (present(gravit_in)            ) gravit           = gravit_in
       if (present(viscosity_dyn_in)     ) viscosity_dyn    = viscosity_dyn_in
+      if (present(tscale_pnd_drain_in)  ) tscale_pnd_drain = tscale_pnd_drain_in
       if (present(Tocnfrz_in)           ) Tocnfrz          = Tocnfrz_in
       if (present(rhofresh_in)          ) rhofresh         = rhofresh_in
       if (present(zvir_in)              ) zvir             = zvir_in
@@ -1165,7 +1168,8 @@ subroutine icepack_query_parameters(   &
          p333_out, p666_out, spval_const_out, pih_out, piq_out, pi2_out, &
          rhos_out, rhoi_out, rhow_out, cp_air_out, emissivity_out, &
          cp_ice_out, cp_ocn_out, hfrazilmin_out, floediam_out, &
-         depressT_out, dragio_out, thickness_ocn_layer1_out, iceruf_ocn_out, albocn_out, gravit_out, viscosity_dyn_out, &
+         depressT_out, dragio_out, thickness_ocn_layer1_out, iceruf_ocn_out, &
+         albocn_out, gravit_out, viscosity_dyn_out, tscale_pnd_drain_out, &
          Tocnfrz_out, rhofresh_out, zvir_out, vonkar_out, cp_wv_out, &
          stefan_boltzmann_out, ice_ref_salinity_out, &
          Tffresh_out, Lsub_out, Lvap_out, Timelt_out, Tsmelt_out, &
@@ -1256,6 +1260,7 @@ subroutine icepack_query_parameters(   &
          cp_ocn_out,     & ! specific heat of ocn    (J/kg/K)
          depressT_out,   & ! Tf:brine salinity ratio (C/ppt)
          viscosity_dyn_out, & ! dynamic viscosity of brine (kg/m/s)
+         tscale_pnd_drain_out, & ! mushy macroscopic drainage timescale (days)
          Tocnfrz_out,    & ! freezing temp of seawater (C)
          Tffresh_out,    & ! freezing temp of fresh ice (K)
          Lsub_out,       & ! latent heat, sublimation freshwater (J/kg)
@@ -1633,6 +1638,7 @@ subroutine icepack_query_parameters(   &
       if (present(albocn_out)            ) albocn_out       = albocn
       if (present(gravit_out)            ) gravit_out       = gravit
       if (present(viscosity_dyn_out)     ) viscosity_dyn_out= viscosity_dyn
+      if (present(tscale_pnd_drain_out)  ) tscale_pnd_drain_out = tscale_pnd_drain
       if (present(Tocnfrz_out)           ) Tocnfrz_out      = Tocnfrz
       if (present(rhofresh_out)          ) rhofresh_out     = rhofresh
       if (present(zvir_out)              ) zvir_out         = zvir
@@ -1831,6 +1837,7 @@ subroutine icepack_write_parameters(iounit)
         write(iounit,*) "  albocn     = ",albocn
         write(iounit,*) "  gravit     = ",gravit
         write(iounit,*) "  viscosity_dyn = ",viscosity_dyn
+        write(iounit,*) "  tscale_pnd_drain = ",tscale_pnd_drain
         write(iounit,*) "  Tocnfrz    = ",Tocnfrz
         write(iounit,*) "  rhofresh   = ",rhofresh
         write(iounit,*) "  zvir       = ",zvir

columnphysics/icepack_therm_mushy.F90

@@ -7,7 +7,7 @@ module icepack_therm_mushy
   use icepack_parameters, only: p01, p05, p1, p2, p5, pi, bignum, puny
   use icepack_parameters, only: viscosity_dyn, rhow, rhoi, rhos, cp_ocn, cp_ice, Lfresh, gravit
   use icepack_parameters, only: hs_min, snwgrain
-  use icepack_parameters, only: a_rapid_mode, Rac_rapid_mode
+  use icepack_parameters, only: a_rapid_mode, Rac_rapid_mode, tscale_pnd_drain
   use icepack_parameters, only: aspect_rapid_mode, dSdt_slow_mode, phi_c_slow_mode
   use icepack_parameters, only: sw_redist, sw_frac, sw_dtemp
   use icepack_mushy_physics, only: icepack_mushy_density_brine, enthalpy_brine, icepack_enthalpy_snow
@@ -3167,9 +3167,11 @@ subroutine flush_pond(w, hpond, apond, dt)
          hpond     ! melt pond thickness (m)
 
     real(kind=dbl_kind), parameter :: &
-         lambda_pond = c1 / (10.0_dbl_kind * 24.0_dbl_kind * 3600.0_dbl_kind), &
          hpond0 = 0.01_dbl_kind
 
+    real(kind=dbl_kind) :: &
+         lambda_pond ! 1 / macroscopic drainage time scale (s)
+
     character(len=*),parameter :: subname='(flush_pond)'
 
     if (tr_pond) then
@@ -3181,6 +3183,7 @@ subroutine flush_pond(w, hpond, apond, dt)
           hpond = max(hpond, c0)
 
           ! exponential decay of pond
+          lambda_pond = c1 / (tscale_pnd_drain * 24.0_dbl_kind * 3600.0_dbl_kind)
           hpond = hpond - lambda_pond * dt * (hpond + hpond0)
 
           hpond = max(hpond, c0)
@@ -3214,11 +3217,11 @@ subroutine flood_ice(hsn,    hin,      &
          hsn               , & ! snow thickness (m)
          hin               , & ! ice thickness (m)
          sss               , & ! sea surface salinity (ppt)
-         qocn                  ! ocean brine enthalpy (J m-2)
+         qocn                  ! ocean brine enthalpy (J m-3)
 
     real(kind=dbl_kind), dimension(:), intent(inout) :: &
-         zqsn              , & ! snow layer enthalpy (J m-2)
-         zqin              , & ! ice layer enthalpy (J m-2)
+         zqsn              , & ! snow layer enthalpy (J m-3)
+         zqin              , & ! ice layer enthalpy (J m-3)
          zSin              , & ! ice layer bulk salinity (ppt)
          phi               , & ! ice liquid fraction
          smice             , & ! ice mass tracer in snow (kg/m^3)
@@ -3246,8 +3249,8 @@ subroutine flood_ice(hsn,    hin,      &
          phi_snowice       , & ! liquid fraction of new snow ice
          rho_snowice       , & ! density of snowice (kg m-3)
          zSin_snowice      , & ! bulk salinity of new snowice (ppt)
-         zqin_snowice      , & ! ice enthalpy of new snowice (J m-2)
-         zqsn_snowice      , & ! snow enthalpy of snow thats becoming snowice (J m-2)
+         zqin_snowice      , & ! ice enthalpy of new snowice (J m-3)
+         zqsn_snowice      , & ! snow enthalpy of snow that is becoming snowice (J m-3)
          freeboard_density , & ! negative of ice surface freeboard times the ocean density (kg m-2)
          ice_mass          , & ! mass of the ice (kg m-2)
 !        snow_mass         , & ! mass of the ice (kg m-2)
@@ -3384,10 +3387,10 @@ subroutine enthalpy_snow_snowice(nslyr, dh, hsn, zqsn, zqsn_snowice)
          hsn          ! initial snow thickness
 
     real(kind=dbl_kind), dimension(:), intent(in) :: &
-         zqsn         ! snow layer enthalpy (J m-2)
+         zqsn         ! snow layer enthalpy (J m-3)
 
     real(kind=dbl_kind), intent(out) :: &
-         zqsn_snowice ! enthalpy of snow becoming snowice (J m-2)
+         zqsn_snowice ! enthalpy of snow becoming snowice (J m-3)
 
     real(kind=dbl_kind) :: &
          rnlyr        ! real value of number of snow layers turning to snowice

doc/source/icepack_index.rst

@@ -121,7 +121,7 @@ either Celsius or Kelvin units).  Deprecated parameters are listed at the end.
    "divu_adv", "divergence associated with advection", "1/s"
    "dms", "dimethyl sulfide concentration", "mmol/m\ :math:`^3`"
    "dmsp", "dimethyl sulfoniopropionate concentration", "mmol/m\ :math:`^3`"
-   "dpscale", ":math:`\bullet` time scale for flushing in permeable ice", ":math:`1\times 10^{-3}`"
+   "dpscale", ":math:`\bullet` scaling factor for flushing in permeable ice (ktherm=1)", ":math:`1\times 10^{-3}`"
    "dragio", "drag coefficient for water on ice", "0.00536"
    "dSdt_slow_mode", ":math:`\bullet` drainage strength parameter", ""
    "dsnow", "change in snow thickness", "m"
@@ -463,6 +463,7 @@ either Celsius or Kelvin units).  Deprecated parameters are listed at the end.
    "time", "total elapsed time", "s"
    "time_forc", "time of last forcing update", "s"
    "Timelt", "melting temperature of ice top surface", "0. C"
+   "tscale_pnd_drain", "mushy pond macroscopic drainage timescale", "10 days"
    "TLAT", "latitude of cell center", "radians"
    "Tliquidus_max", "maximum liquidus temperature of mush", "0. C"
    "TLON", "longitude of cell center", "radians"

doc/source/user_guide/ug_case_settings.rst

@@ -219,6 +219,7 @@ thermo_nml
    "``phi_i_mushy``", ":math:`0<\phi_i < 1`", "solid fraction at lower boundary", "0.85"
    "``Rac_rapid_mode``", "real", "critical Rayleigh number", "10.0"
    "``Tliquidus_max``", "real", "maximum liquidus temperature of mush (C)", "0.0"
+   "``tscale_pnd_drain``", "real", "mushy pond macroscopic drainage timescale in days", "10."
    "", "", "", ""
 
 
@@ -279,7 +280,7 @@ ponds_nml
    :widths: 15, 15, 30, 15 
 
    "", "", "", ""
-   "``dpscale``", "real", "time scale for flushing in permeable ice", "1.0"
+   "``dpscale``", "real", "scaling factor for flushing in permeable ice (ktherm=1)", "1.e-3"
    "``frzpnd``", "``cesm``", "CESM pond refreezing forumulation", "``cesm``"
    "", "``hlid``", "Stefan refreezing with pond ice thickness", ""
    "``hp1``", "real", "critical ice lid thickness for topo ponds in m", "0.01"