Add saltflux option.

columnphysics/icepack_itd.F90

@@ -31,10 +31,10 @@ module icepack_itd
       use icepack_parameters, only: rhosi, sk_l, hs_ssl, min_salin, rsnw_fall
       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, nt_bgc_S, bio_index
-      use icepack_tracers,    only: n_iso, tr_iso, tr_snow, nt_smice, nt_rsnw, nt_rhos
+      use icepack_tracers,    only: n_iso, tr_iso, tr_snow, nt_smice, nt_rsnw, nt_rhos, nt_sice
       use icepack_tracers,    only: icepack_compute_tracers
       use icepack_parameters, only: solve_zsal, skl_bgc, z_tracers
-      use icepack_parameters, only: kcatbound, kitd
+      use icepack_parameters, only: kcatbound, kitd, saltflux_option
       use icepack_therm_shared, only: Tmin, hi_min
       use icepack_warnings,   only: warnstr, icepack_warnings_add
       use icepack_warnings,   only: icepack_warnings_setabort, icepack_warnings_aborted
@@ -1089,7 +1089,7 @@ subroutine zap_small_areas (dt,        ntrcr,        &
          n, k, it, & !counting indices
          blevels
 
-      real (kind=dbl_kind) :: xtmp      ! temporary variables
+      real (kind=dbl_kind) :: xtmp, sicen      ! temporary variables
       real (kind=dbl_kind) , dimension (1):: trcr_skl
       real (kind=dbl_kind) , dimension (nblyr+1):: bvol
 
@@ -1187,7 +1187,15 @@ subroutine zap_small_areas (dt,        ntrcr,        &
             xtmp = (rhoi*vicen(n)) / dt
             dfresh = dfresh + xtmp
 
-            xtmp = rhoi*vicen(n)*ice_ref_salinity*p001 / dt
+            if (saltflux_option == 'prognostic') then
+               sicen = c0
+               do k=1,nilyr
+                  sicen = sicen + trcrn(nt_sice+k-1,n) / real(nilyr,kind=dbl_kind)
+               enddo
+               xtmp = rhoi*vicen(n)*sicen*p001 / dt
+            else
+               xtmp = rhoi*vicen(n)*ice_ref_salinity*p001 / dt
+            endif
             dfsalt = dfsalt + xtmp
 
             aice0 = aice0 + aicen(n)
@@ -1301,13 +1309,22 @@ subroutine zap_small_areas (dt,        ntrcr,        &
       !-----------------------------------------------------------------
 
             xtmp = (rhoi*vicen(n) + rhos*vsnon(n)) &
-                 * (aice-c1)/aice / dt
-            dfresh = dfresh + xtmp
-
-            xtmp = rhoi*vicen(n)*ice_ref_salinity*p001 &
-                 * (aice-c1)/aice / dt
-            dfsalt = dfsalt + xtmp
-
+                 * (aice-c1)/aice / dt 
+            dfresh = dfresh + xtmp 
+
+            if (saltflux_option == 'prognostic') then
+               sicen = c0
+               do k=1,nilyr
+                  sicen = sicen + trcrn(nt_sice+k-1,n) / real(nilyr,kind=dbl_kind)
+               enddo
+               xtmp = rhoi*vicen(n)*sicen*p001 &
+                    * (aice-c1)/aice / dt
+            else
+               xtmp = rhoi*vicen(n)*ice_ref_salinity*p001 &
+                    * (aice-c1)/aice / dt
+            endif
+            dfsalt = dfsalt + xtmp 
+
             if (solve_zsal) then
             do k = 1,nblyr
                xtmp = rhosi*trcrn(nt_fbri,n)*vicen(n)*p001&

columnphysics/icepack_parameters.F90

@@ -165,6 +165,11 @@ module icepack_parameters
                                   ! 'linear_salt' = -depressT * sss
                                   ! 'mushy' conforms with ktherm=2
 
+      character(len=char_len), public :: &
+         saltflux_option  = 'constant'! Salt flux computation
+                                      ! 'constant' reference value of ice_ref_salinity
+                                      ! 'prognostic' prognostic salt flux
+
 !-----------------------------------------------------------------------
 ! Parameters for radiation
 !-----------------------------------------------------------------------
@@ -454,6 +459,7 @@ subroutine icepack_init_parameters(   &
          atmbndy_in, calc_strair_in, formdrag_in, highfreq_in, natmiter_in, &
          atmiter_conv_in, calc_dragio_in, &
          tfrz_option_in, kitd_in, kcatbound_in, hs0_in, frzpnd_in, &
+         saltflux_option_in, &
          floeshape_in, wave_spec_in, wave_spec_type_in, nfreq_in, &
          dpscale_in, rfracmin_in, rfracmax_in, pndaspect_in, hs1_in, hp1_in, &
          bgc_flux_type_in, z_tracers_in, scale_bgc_in, solve_zbgc_in, &
@@ -557,12 +563,17 @@ subroutine icepack_init_parameters(   &
          dSdt_slow_mode_in    , & ! slow mode drainage strength (m s-1 K-1)
          phi_c_slow_mode_in   , & ! liquid fraction porosity cutoff for slow mode
          phi_i_mushy_in           ! liquid fraction of congelation ice
+
+      character(len=*), intent(in), optional :: &
+         tfrz_option_in              ! form of ocean freezing temperature
+                                     ! 'minus1p8' = -1.8 C
+                                     ! 'linear_salt' = -depressT * sss
+                                     ! 'mushy' conforms with ktherm=2
 
-        character(len=*), intent(in), optional :: &
-             tfrz_option_in              ! form of ocean freezing temperature
-                                         ! 'minus1p8' = -1.8 C
-                                         ! 'linear_salt' = -depressT * sss
-                                         ! 'mushy' conforms with ktherm=2
+      character(len=*), intent(in), optional :: &
+         saltflux_option_in         ! Salt flux computation
+                                    ! 'constant' reference value of ice_ref_salinity
+                                    ! 'prognostic' prognostic salt flux
 
 !-----------------------------------------------------------------------
 ! Parameters for radiation
@@ -917,6 +928,7 @@ subroutine icepack_init_parameters(   &
       if (present(natmiter_in)          ) natmiter         = natmiter_in
       if (present(atmiter_conv_in)      ) atmiter_conv     = atmiter_conv_in
       if (present(tfrz_option_in)       ) tfrz_option      = tfrz_option_in
+      if (present(saltflux_option_in)   ) saltflux_option  = saltflux_option_in
       if (present(kitd_in)              ) kitd             = kitd_in
       if (present(kcatbound_in)         ) kcatbound        = kcatbound_in
       if (present(floeshape_in)         ) floeshape        = floeshape_in
@@ -1145,6 +1157,7 @@ subroutine icepack_query_parameters(   &
          atmbndy_out, calc_strair_out, formdrag_out, highfreq_out, natmiter_out, &
          atmiter_conv_out, calc_dragio_out, &
          tfrz_option_out, kitd_out, kcatbound_out, hs0_out, frzpnd_out, &
+         saltflux_option_out, &
          floeshape_out, wave_spec_out, wave_spec_type_out, nfreq_out, &
          dpscale_out, rfracmin_out, rfracmax_out, pndaspect_out, hs1_out, hp1_out, &
          bgc_flux_type_out, z_tracers_out, scale_bgc_out, solve_zbgc_out, &
@@ -1264,6 +1277,12 @@ subroutine icepack_query_parameters(   &
                                       ! 'linear_salt' = -depressT * sss
                                       ! 'mushy' conforms with ktherm=2
 
+      character(len=*), intent(out), optional :: &
+         saltflux_option_out         ! Salt flux computation
+                                     ! 'constant' reference value of ice_ref_salinity
+                                     ! 'prognostic' prognostic salt flux
+
+
 !-----------------------------------------------------------------------
 ! Parameters for radiation
 !-----------------------------------------------------------------------
@@ -1653,6 +1672,7 @@ subroutine icepack_query_parameters(   &
       if (present(natmiter_out)          ) natmiter_out     = natmiter
       if (present(atmiter_conv_out)      ) atmiter_conv_out = atmiter_conv
       if (present(tfrz_option_out)       ) tfrz_option_out  = tfrz_option
+      if (present(saltflux_option_out)   ) saltflux_option_out = saltflux_option
       if (present(kitd_out)              ) kitd_out         = kitd
       if (present(kcatbound_out)         ) kcatbound_out    = kcatbound
       if (present(floeshape_out)         ) floeshape_out    = floeshape
@@ -1856,6 +1876,7 @@ subroutine icepack_write_parameters(iounit)
         write(iounit,*) "  natmiter      = ", natmiter
         write(iounit,*) "  atmiter_conv  = ", atmiter_conv
         write(iounit,*) "  tfrz_option   = ", tfrz_option
+        write(iounit,*) "  saltflux_option = ", saltflux_option
         write(iounit,*) "  kitd          = ", kitd
         write(iounit,*) "  kcatbound     = ", kcatbound
         write(iounit,*) "  floeshape     = ", floeshape

columnphysics/icepack_therm_itd.F90

@@ -27,11 +27,12 @@ module icepack_therm_itd
       use icepack_parameters, only: rhosi, conserv_check, rhosmin
       use icepack_parameters, only: kitd, ktherm
       use icepack_parameters, only: z_tracers, solve_zsal, hfrazilmin
+      use icepack_parameters, only: saltflux_option
 
       use icepack_tracers, only: ntrcr, nbtrcr
       use icepack_tracers, only: nt_qice, nt_qsno, nt_fbri, nt_sice
       use icepack_tracers, only: nt_apnd, nt_hpnd, nt_aero, nt_isosno, nt_isoice
-      use icepack_tracers, only: nt_Tsfc, nt_iage, nt_FY, nt_fsd, nt_rhos
+      use icepack_tracers, only: nt_Tsfc, nt_iage, nt_FY, nt_fsd, nt_rhos, nt_sice
       use icepack_tracers, only: nt_alvl, nt_vlvl
       use icepack_tracers, only: tr_pond_lvl, tr_pond_topo, tr_snow
       use icepack_tracers, only: tr_iage, tr_FY, tr_lvl, tr_aero, tr_iso, tr_brine, tr_fsd
@@ -984,7 +985,7 @@ subroutine lateral_melt (dt,         ncat,       &
       real (kind=dbl_kind), intent(in) :: &
          sss
       real (kind=dbl_kind) :: &
-         Ti, Si0, qi0, &
+         Ti, Si0, qi0, sicen, &
          elapsed_t,    & ! FSD subcycling
          subdt           ! FSD timestep (s)
 
@@ -1093,7 +1094,15 @@ subroutine lateral_melt (dt,         ncat,       &
             ! dfresh > 0, dfsalt > 0, dfpond > 0
 
             dfresh = (rhoi*vicen(n) + rhos*vsnon(n))      * rsiden(n) / dt
-            dfsalt =  rhoi*vicen(n)*ice_ref_salinity*p001 * rsiden(n) / dt
+            if (saltflux_option == 'prognostic') then
+               sicen = c0
+               do k=1,nilyr
+                  sicen = sicen + trcrn(nt_sice+k-1,n) / real(nilyr,kind=dbl_kind)
+               enddo
+               dfsalt = rhoi*vicen(n)*sicen*p001 * rsiden(n) / dt
+            else
+               dfsalt = rhoi*vicen(n)*ice_ref_salinity*p001 * rsiden(n) / dt
+            endif
             fresh  = fresh + dfresh
             fsalt  = fsalt + dfsalt
 
@@ -1596,18 +1605,23 @@ subroutine add_new_ice (ncat,      nilyr,      &
       !-----------------------------------------------------------------
 
       if (update_ocn_f) then
-         if (ktherm <= 1) then
-            dfresh = -rhoi*vi0new/dt
+         dfresh = -rhoi*vi0new/dt 
+         if (saltflux_option == 'prognostic') then
+            dfsalt = Si0new*p001*dfresh
+         else
             dfsalt = ice_ref_salinity*p001*dfresh
-            fresh  = fresh + dfresh
-            fsalt  = fsalt + dfsalt
-         ! elseif (ktherm == 2) the fluxes are added elsewhere
          endif
+         fresh  = fresh + dfresh
+         fsalt  = fsalt + dfsalt
       else ! update_ocn_f = false
          if (ktherm == 2) then ! return mushy-layer frazil to ocean (POP)
             vi0tmp = fnew*dt / (rhoi*Lfresh)
             dfresh = -rhoi*(vi0new - vi0tmp)/dt
-            dfsalt = ice_ref_salinity*p001*dfresh
+            if (saltflux_option == 'prognostic') then
+               dfsalt = Si0new*p001*dfresh
+            else
+               dfsalt = ice_ref_salinity*p001*dfresh
+            endif
             fresh  = fresh + dfresh
             fsalt  = fsalt + dfsalt
             frazil_diag = frazil - vi0tmp

columnphysics/icepack_therm_vertical.F90

@@ -27,6 +27,7 @@ module icepack_therm_vertical
       use icepack_parameters, only: ustar_min, fbot_xfer_type, formdrag, calc_strair
       use icepack_parameters, only: rfracmin, rfracmax, dpscale, frzpnd, snwgrain, snwlvlfac
       use icepack_parameters, only: phi_i_mushy, floeshape, floediam, use_smliq_pnd, snwredist
+      use icepack_parameters, only: saltflux_option
 
       use icepack_tracers, only: tr_iage, tr_FY, tr_aero, tr_pond, tr_fsd, tr_iso
       use icepack_tracers, only: tr_pond_lvl, tr_pond_topo
@@ -239,7 +240,7 @@ subroutine thermo_vertical (nilyr,       nslyr,     &
          einter          ! intermediate energy
 
       real (kind=dbl_kind) :: &
-         fadvocn ! advective heat flux to ocean
+         fadvocn, saltvol, dfsalt ! advective heat flux to ocean
 
       character(len=*),parameter :: subname='(thermo_vertical)'
 
@@ -292,6 +293,14 @@ subroutine thermo_vertical (nilyr,       nslyr,     &
       worki = hin
       works = hsn
 
+      ! Save initial salt volume for prognostic flux
+      if (saltflux_option == 'prognostic') then
+         saltvol = c0
+         do k=1,nilyr
+            saltvol = saltvol + rhoi*zSin(k)*hin*p001 / real(nilyr,kind=dbl_kind)
+         enddo
+      endif
+
       !-----------------------------------------------------------------
       ! Compute new surface temperature and internal ice and snow
       !  temperatures.
@@ -435,8 +444,16 @@ subroutine thermo_vertical (nilyr,       nslyr,     &
       dhs = hsn - works - hsn_new
 
       freshn = freshn + evapn - (rhoi*dhi + rhos*dhs) / dt
-      fsaltn = fsaltn - rhoi*dhi*ice_ref_salinity*p001/dt
-      fhocnn = fhocnn + fadvocn ! for ktherm=2
+      if (saltflux_option == 'prognostic') then
+         dfsalt = c0
+         do k=1,nilyr
+            dfsalt = dfsalt + rhoi*zSin(k)*hin*p001 / real(nilyr,kind=dbl_kind)
+         enddo
+         fsaltn = fsaltn - (dfsalt - saltvol) / dt
+      else
+         fsaltn = fsaltn - rhoi*dhi*ice_ref_salinity*p001/dt
+      endif
+      fhocnn = fhocnn + fadvocn ! for ktherm=2 
 
       if (hin == c0) then
          if (tr_pond_topo) fpond = fpond - aicen * apond * hpond
@@ -2569,7 +2586,6 @@ subroutine icepack_step_therm1(dt, ncat, nilyr, nslyr,    &
          enddo
       endif
 
-
       !-----------------------------------------------------------------
       ! Update the neutral drag coefficients to account for form drag
       ! Oceanic and atmospheric drag coefficients

configuration/driver/icedrv_init.F90

@@ -97,14 +97,17 @@ subroutine input_data
          natmiter, kitd, kcatbound
 
       character (len=char_len) :: shortwave, albedo_type, conduct, fbot_xfer_type, &
-         tfrz_option, frzpnd, atmbndy, wave_spec_type, snwredist, snw_aging_table
+         tfrz_option, saltflux_option, frzpnd, atmbndy, wave_spec_type, snwredist, snw_aging_table
 
       logical (kind=log_kind) :: sw_redist, use_smliq_pnd, snwgrain
       real (kind=dbl_kind)    :: sw_frac, sw_dtemp
 
       ! Flux convergence tolerance
       real (kind=dbl_kind) :: atmiter_conv
 
+      ! Ice reference salinity for fluxes
+      real (kind=dbl_kind) :: ice_ref_salinity
+
       logical (kind=log_kind) :: calc_Tsfc, formdrag, highfreq, calc_strair, calc_dragio
       logical (kind=log_kind) :: conserv_check
 
@@ -168,7 +171,8 @@ subroutine input_data
         fbot_xfer_type,  oceanmixed_ice,  emissivity,      &
         formdrag,        highfreq,        natmiter,        &
         atmiter_conv,    calc_dragio,                      &
-        tfrz_option,     default_season,  wave_spec_type,  &
+        tfrz_option,     saltflux_option, ice_ref_salinity, &
+        default_season,  wave_spec_type,  &
         precip_units,    fyear_init,      ycycle,          &
         atm_data_type,   ocn_data_type,   bgc_data_type,   &
         atm_data_file,   ocn_data_file,   bgc_data_file,   &
@@ -216,7 +220,8 @@ subroutine input_data
            dSdt_slow_mode_out=dSdt_slow_mode, &
            phi_c_slow_mode_out=phi_c_slow_mode, &
            phi_i_mushy_out=phi_i_mushy, conserv_check_out=conserv_check, &
-           tfrz_option_out=tfrz_option, kalg_out=kalg, &
+           tfrz_option_out=tfrz_option, saltflux_option_out=saltflux_option, &
+           ice_ref_salinity_out=ice_ref_salinity, kalg_out=kalg, &
            fbot_xfer_type_out=fbot_xfer_type, puny_out=puny, &
            wave_spec_type_out=wave_spec_type, &
            sw_redist_out=sw_redist, sw_frac_out=sw_frac, sw_dtemp_out=sw_dtemp, &
@@ -573,6 +578,13 @@ subroutine input_data
          'WARNING: For consistency, set tfrz_option = mushy'
       endif
 
+      if (ktherm == 1 .and. trim(saltflux_option) /= 'constant') then
+         write (nu_diag,*) &
+         'WARNING: ktherm = 1 and saltflux_option = ',trim(saltflux_option)
+         write (nu_diag,*) &
+         'WARNING: For consistency, set saltflux_option = constant'
+      endif
+
       if (ktherm == 0) then
          write (nu_diag,*) 'WARNING: ktherm = 0 zero-layer thermodynamics'
          write (nu_diag,*) 'WARNING: has been deprecated'
@@ -761,6 +773,11 @@ subroutine input_data
          write(nu_diag,1010) ' oceanmixed_ice            = ', oceanmixed_ice
          write(nu_diag,*)    ' tfrz_option               = ', &
                                trim(tfrz_option)
+         write(nu_diag,*)    ' saltflux_option           = ', &
+                               trim(saltflux_option)
+         if (trim(saltflux_option) == 'constant') then
+            write(nu_diag,1005)    ' ice_ref_salinity          = ', ice_ref_salinity
+         endif
          write(nu_diag,1010) ' restore_ocn               = ', restore_ocn
          if (restore_ocn) &
          write(nu_diag,1005) ' trestore                  = ', trestore
@@ -942,7 +959,8 @@ subroutine input_data
            dSdt_slow_mode_in=dSdt_slow_mode, &
            phi_c_slow_mode_in=phi_c_slow_mode, &
            phi_i_mushy_in=phi_i_mushy, conserv_check_in=conserv_check, &
-           tfrz_option_in=tfrz_option, kalg_in=kalg, &
+           tfrz_option_in=tfrz_option, saltflux_option_in=saltflux_option, &
+           ice_ref_salinity_in=ice_ref_salinity, kalg_in=kalg, &
            fbot_xfer_type_in=fbot_xfer_type, &
            wave_spec_type_in=wave_spec_type, wave_spec_in=wave_spec, &
            sw_redist_in=sw_redist, sw_frac_in=sw_frac, sw_dtemp_in=sw_dtemp, &

configuration/scripts/icepack_in

@@ -109,6 +109,8 @@
     update_ocn_f    = .false.
     l_mpond_fresh   = .false.
     tfrz_option     = 'mushy'
+    ice_ref_salinity = 4.0
+    saltflux_option = 'constant'
     oceanmixed_ice  = .true.
     wave_spec_type  = 'none'
     restore_ocn     = .false.

doc/source/icepack_index.rst

@@ -403,6 +403,7 @@ either Celsius or Kelvin units).  Deprecated parameters are listed at the end.
    "rsnw_sig", "standard deviation of snow grain radius", "250. :math:`\times` 10\ :math:`^{-6}` m"
    "**S**", "", ""
    "salinz", "ice salinity profile", "ppt"
+   "saltflux_option", "constant or prognostic salinity fluxes","constant"
    "saltmax", "max salinity, at ice base (:cite:`Bitz99`)", "3.2 ppt"
    "scale_factor", "scaling factor for shortwave radiation components", ""
    "sec", "seconds elasped into idate", ""

doc/source/user_guide/ug_case_settings.rst

@@ -358,7 +358,9 @@ forcing_nml
    "", "``mm_per_sec``", "(same as MKS units)", ""
    "", "``m_per_sec``", "", ""
    "``restore_ocn``", "logical", "restore sst to data", "``.false.``"
-   "``tfrz_option``", "``linear_salt``", "linear functino of salinity (ktherm=1)", "``mushy``"
+   "``saltflux_option``", "``constant``","salt flux is referenced to a constant salinity","``constant``"
+   "","``prognostic``","use actual sea ice bulk salinity in flux"
+   "``tfrz_option``", "``linear_salt``", "linear function of salinity (ktherm=1)", "``mushy``"
    "", "``minus1p8``", "constant ocean freezing temperature (:math:`-1.8^{\circ} C`)", ""
    "", "``mushy``", "matches mushy-layer thermo (ktherm=2)", ""
    "``trestore``", "integer", "sst restoring time scale (days)", "90"