Merge pull request #13 from E3SM-Project/Tliquidus_limit

add namelist parameter to set maximum liquidus temperature

columnphysics/icepack_parameters.F90

@@ -129,6 +129,7 @@ module icepack_parameters
          phi_init  = 0.75_dbl_kind    ,&! initial liquid fraction of frazil
          min_salin = p1               ,&! threshold for brine pocket treatment
          salt_loss = 0.4_dbl_kind     ,&! fraction of salt retained in zsalinity
+         Tliquidus_max = c0           ,&! maximum liquidus temperature of mush (C)
          dSin0_frazil = c3            ,&! bulk salinity reduction of newly formed frazil
          dts_b     = 50._dbl_kind     ,&! zsalinity timestep
          ustar_min = 0.005_dbl_kind   ,&! minimum friction velocity for ocean heat flux (m/s)
@@ -455,6 +456,7 @@ subroutine icepack_init_parameters(   &
 #endif
          zref_in, hs_min_in, snowpatch_in, rhosi_in, sk_l_in, &
          saltmax_in, phi_init_in, min_salin_in, salt_loss_in, &
+         Tliquidus_max_in, &
          min_bgc_in, dSin0_frazil_in, hi_ssl_in, hs_ssl_in, &
          awtvdr_in, awtidr_in, awtvdf_in, awtidf_in, &
          qqqice_in, TTTice_in, qqqocn_in, TTTocn_in, &
@@ -549,6 +551,7 @@ subroutine icepack_init_parameters(   &
          phi_init_in,   & ! initial liquid fraction of frazil
          min_salin_in,  & ! threshold for brine pocket treatment
          salt_loss_in,  & ! fraction of salt retained in zsalinity
+         Tliquidus_max_in, & ! maximum liquidus temperature of mush (C)
          dSin0_frazil_in  ! bulk salinity reduction of newly formed frazil
 
       integer (kind=int_kind), intent(in), optional :: &
@@ -910,6 +913,7 @@ subroutine icepack_init_parameters(   &
       if (present(phi_init_in)          ) phi_init         = phi_init_in
       if (present(min_salin_in)         ) min_salin        = min_salin_in
       if (present(salt_loss_in)         ) salt_loss        = salt_loss_in
+      if (present(Tliquidus_max_in)     ) Tliquidus_max    = Tliquidus_max_in
       if (present(min_bgc_in)           ) min_bgc          = min_bgc_in
       if (present(dSin0_frazil_in)      ) dSin0_frazil     = dSin0_frazil_in
       if (present(hi_ssl_in)            ) hi_ssl           = hi_ssl_in
@@ -1173,6 +1177,7 @@ subroutine icepack_query_parameters(   &
 #endif
          zref_out, hs_min_out, snowpatch_out, rhosi_out, sk_l_out, &
          saltmax_out, phi_init_out, min_salin_out, salt_loss_out, &
+         Tliquidus_max_out, &
          min_bgc_out, dSin0_frazil_out, hi_ssl_out, hs_ssl_out, &
          awtvdr_out, awtidr_out, awtvdf_out, awtidf_out, &
          qqqice_out, TTTice_out, qqqocn_out, TTTocn_out, update_ocn_f_out, &
@@ -1276,6 +1281,7 @@ subroutine icepack_query_parameters(   &
          phi_init_out,   & ! initial liquid fraction of frazil
          min_salin_out,  & ! threshold for brine pocket treatment
          salt_loss_out,  & ! fraction of salt retained in zsalinity
+         Tliquidus_max_out, & ! maximum liquidus temperature of mush (C)
          dSin0_frazil_out  ! bulk salinity reduction of newly formed frazil
 
       integer (kind=int_kind), intent(out), optional :: &
@@ -1669,6 +1675,7 @@ subroutine icepack_query_parameters(   &
       if (present(phi_init_out)          ) phi_init_out     = phi_init
       if (present(min_salin_out)         ) min_salin_out    = min_salin
       if (present(salt_loss_out)         ) salt_loss_out    = salt_loss
+      if (present(Tliquidus_max_out)     ) Tliquidus_max_out= Tliquidus_max
       if (present(min_bgc_out)           ) min_bgc_out      = min_bgc
       if (present(dSin0_frazil_out)      ) dSin0_frazil_out = dSin0_frazil
       if (present(hi_ssl_out)            ) hi_ssl_out       = hi_ssl
@@ -1869,6 +1876,7 @@ subroutine icepack_write_parameters(iounit)
         write(iounit,*) "  phi_init  = ",phi_init
         write(iounit,*) "  min_salin = ",min_salin
         write(iounit,*) "  salt_loss = ",salt_loss
+        write(iounit,*) "  Tliquidus_max = ",Tliquidus_max
         write(iounit,*) "  min_bgc   = ",min_bgc
         write(iounit,*) "  dSin0_frazil = ",dSin0_frazil
         write(iounit,*) "  hi_ssl = ",hi_ssl

columnphysics/icepack_therm_itd.F90

@@ -24,6 +24,7 @@ module icepack_therm_itd
       use icepack_parameters, only: p001, p1, p333, p5, p666, puny, bignum
       use icepack_parameters, only: rhos, rhoi, Lfresh, ice_ref_salinity
       use icepack_parameters, only: phi_init, dsin0_frazil, hs_ssl, salt_loss
+      use icepack_parameters, only: Tliquidus_max
       use icepack_parameters, only: rhosi, conserv_check, rhosmin
 #ifdef UNDEPRECATE_0LAYER
       use icepack_parameters, only: kitd, ktherm, heat_capacity
@@ -1038,7 +1039,7 @@ subroutine lateral_melt (dt,         ncat,       &
          else
             Si0 = sss**2 / (c4*dSin0_frazil)
          endif
-         Ti = min(liquidus_temperature_mush(Si0/phi_init), -p1)
+         Ti = min(liquidus_temperature_mush(Si0/phi_init), Tliquidus_max)
          qi0 = enthalpy_mush(Ti, Si0)
 
          do n = 1, ncat
@@ -1569,7 +1570,7 @@ subroutine add_new_ice (ncat,      nilyr,      &
          do k = 1, nilyr
             Sprofile(k) = Si0new
          enddo
-         Ti = min(liquidus_temperature_mush(Si0new/phi_init), -p1)
+         Ti = min(liquidus_temperature_mush(Si0new/phi_init), Tliquidus_max)
          qi0new = enthalpy_mush(Ti, Si0new)
       else
          do k = 1, nilyr

configuration/driver/icedrv_init.F90

@@ -88,7 +88,7 @@ subroutine input_data
       character (len=char_len), dimension(4) :: nx_names_default
 
       real (kind=dbl_kind) :: ustar_min, albicev, albicei, albsnowv, albsnowi, &
-         ahmax, R_ice, R_pnd, R_snw, dT_mlt, rsnw_mlt, ksno, hi_min, &
+         ahmax, R_ice, R_pnd, R_snw, dT_mlt, rsnw_mlt, ksno, hi_min, Tliquidus_max, &
          mu_rdg, hs0, dpscale, rfracmin, rfracmax, pndaspect, hs1, hp1, &
          a_rapid_mode, Rac_rapid_mode, aspect_rapid_mode, dSdt_slow_mode, &
          phi_c_slow_mode, phi_i_mushy, kalg, emissivity, floediam, hfrazilmin, &
@@ -150,7 +150,7 @@ subroutine input_data
         kitd,           ktherm,          ksno,     conduct,             &
         a_rapid_mode,   Rac_rapid_mode,  aspect_rapid_mode,             &
         dSdt_slow_mode, phi_c_slow_mode, phi_i_mushy,                   &
-        floediam,       hfrazilmin
+        floediam,       hfrazilmin,      Tliquidus_max
 
       namelist /dynamics_nml/ &
         kstrength,      krdg_partic,    krdg_redist,    mu_rdg,         &
@@ -227,7 +227,7 @@ subroutine input_data
            Rac_rapid_mode_out=Rac_rapid_mode, &
            aspect_rapid_mode_out=aspect_rapid_mode, &
            dSdt_slow_mode_out=dSdt_slow_mode, &
-           phi_c_slow_mode_out=phi_c_slow_mode, &
+           phi_c_slow_mode_out=phi_c_slow_mode, Tliquidus_max_out=Tliquidus_max, &
            phi_i_mushy_out=phi_i_mushy, conserv_check_out=conserv_check, &
            tfrz_option_out=tfrz_option, kalg_out=kalg, &
            fbot_xfer_type_out=fbot_xfer_type, puny_out=puny, &
@@ -759,6 +759,7 @@ subroutine input_data
          write(nu_diag,1005) ' dSdt_slow_mode            = ', dSdt_slow_mode
          write(nu_diag,1005) ' phi_c_slow_mode           = ', phi_c_slow_mode
          write(nu_diag,1005) ' phi_i_mushy               = ', phi_i_mushy
+         write(nu_diag,1005) ' Tliquidus_max             = ', Tliquidus_max
          endif
 
          write(nu_diag,1030) ' atmbndy                   = ', trim(atmbndy)
@@ -976,7 +977,7 @@ subroutine input_data
            Rac_rapid_mode_in=Rac_rapid_mode, &
            aspect_rapid_mode_in=aspect_rapid_mode, &
            dSdt_slow_mode_in=dSdt_slow_mode, &
-           phi_c_slow_mode_in=phi_c_slow_mode, &
+           phi_c_slow_mode_in=phi_c_slow_mode, Tliquidus_max_in=Tliquidus_max, &
            phi_i_mushy_in=phi_i_mushy, conserv_check_in=conserv_check, &
            tfrz_option_in=tfrz_option, kalg_in=kalg, &
            fbot_xfer_type_in=fbot_xfer_type, &

configuration/scripts/icepack_in

@@ -45,6 +45,7 @@
     dSdt_slow_mode    = -5.0e-8
     phi_c_slow_mode   =    0.05
     phi_i_mushy       =    0.85
+    Tliquidus_max     = -0.1d0
     floediam          = 300.0d0
     hfrazilmin        = 0.05d0
 /

doc/source/icepack_index.rst

@@ -462,6 +462,7 @@ either Celsius or Kelvin units).
    "time_forc", "time of last forcing update", "s"
    "Timelt", "melting temperature of ice top surface", "0. C"
    "TLAT", "latitude of cell center", "radians"
+   "Tliquidus_max", "maximum liquidus temperature of mush", "0. C"
    "TLON", "longitude of cell center", "radians"
    "tmask", "land/boundary mask, thickness (T-cell)", ""
    "tmass", "total mass of ice and snow", "kg/m\ :math:`^2`"

doc/source/user_guide/ug_case_settings.rst

@@ -218,6 +218,7 @@ thermo_nml
    "``phi_c_slow_mode``", ":math:`0<\phi_c < 1`", "critical liquid fraction", "0.05"
    "``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"
    "", "", "", ""
 
 ..