Update to Icepack1.3.3 release, #8f96707, Dec 9, 2022

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, hi_min
-      use icepack_parameters, only: kcatbound, kitd
+      use icepack_parameters, only: kcatbound, kitd, saltflux_option
       use icepack_therm_shared, only: Tmin
       use icepack_warnings,   only: warnstr, icepack_warnings_add
       use icepack_warnings,   only: icepack_warnings_setabort, icepack_warnings_aborted
@@ -1101,7 +1101,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
 
@@ -1199,7 +1199,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)
@@ -1316,8 +1324,17 @@ subroutine zap_small_areas (dt,        ntrcr,        &
                  * (aice-c1)/aice / dt
             dfresh = dfresh + xtmp
 
-            xtmp = rhoi*vicen(n)*ice_ref_salinity*p001 &
-                 * (aice-c1)/aice / 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 &
+                    * (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

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
 !-----------------------------------------------------------------------
@@ -455,6 +460,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, &
@@ -561,12 +567,16 @@ subroutine icepack_init_parameters(   &
          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
-                                         ! 'constant' = Tocnfrz
-                                         ! '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
@@ -931,6 +941,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
@@ -1160,6 +1171,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, &
@@ -1282,6 +1294,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
 !-----------------------------------------------------------------------
@@ -1677,6 +1695,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
@@ -1883,6 +1902,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

@@ -28,11 +28,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, hi_min
+      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
@@ -985,7 +986,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)
 
@@ -1094,7 +1095,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
 
@@ -1597,18 +1606,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,8 @@ 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
       use icepack_tracers, only: n_aero, n_iso
@@ -237,7 +239,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)'
 
@@ -290,6 +292,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.
@@ -431,7 +441,15 @@ 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
+      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
@@ -2565,7 +2583,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

columnphysics/version.txt

@@ -1 +1 @@
-ICEPACK 1.3.2
+ICEPACK 1.3.3

configuration/driver/icedrv_InitMod.F90

@@ -44,7 +44,7 @@ subroutine icedrv_initialize
       use icedrv_flux, only: init_coupler_flux, init_history_therm, &
           init_flux_atm_ocn
       use icedrv_forcing, only: init_forcing, get_forcing, get_wave_spec
-      use icedrv_forcing_bgc, only: get_forcing_bgc, faero_default, fiso_default, init_forcing_bgc 
+      use icedrv_forcing_bgc, only: get_forcing_bgc, faero_default, fiso_default, init_forcing_bgc
       use icedrv_restart_shared, only: restart
       use icedrv_init, only: input_data, init_state, init_grid2, init_fsd
       use icedrv_init_column, only: init_thermo_vertical, init_shortwave, init_zbgc
@@ -138,7 +138,7 @@ subroutine icedrv_initialize
       if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
           file=__FILE__,line= __LINE__)
 
-      call init_forcing      ! initialize forcing (standalone)     
+      call init_forcing      ! initialize forcing (standalone)
       if (skl_bgc .or. z_tracers) call init_forcing_bgc !cn
       if (tr_fsd .and. wave_spec) call get_wave_spec ! wave spectrum in ice
       call get_forcing(istep1)       ! get forcing from data arrays

configuration/driver/icedrv_MAIN.F90

@@ -1,24 +1,24 @@
 !=======================================================================
 ! Copyright (c) 2022, Triad National Security, LLC
 ! All rights reserved.
-!                
+!
 ! Copyright 2022. Triad National Security, LLC. This software was
-! produced under U.S. Government contract DE-AC52-06NA25396 for Los 
+! produced under U.S. Government contract DE-AC52-06NA25396 for Los
 ! Alamos National Laboratory (LANL), which is operated by Triad
-! National Security, LLC for the U.S. Department of Energy. The U.S.  
-! Government has rights to use, reproduce, and distribute this software.  
-! NEITHER THE GOVERNMENT NOR TRIAD NATIONAL SECURITY, LLC MAKES ANY  
+! National Security, LLC for the U.S. Department of Energy. The U.S.
+! Government has rights to use, reproduce, and distribute this software.
+! NEITHER THE GOVERNMENT NOR TRIAD NATIONAL SECURITY, LLC MAKES ANY
 ! WARRANTY, EXPRESS OR IMPLIED, OR ASSUMES ANY LIABILITY FOR THE USE OF
-! THIS SOFTWARE. If software is modified to produce derivative works, 
-! such modified software should be clearly marked, so as not to confuse 
+! THIS SOFTWARE. If software is modified to produce derivative works,
+! such modified software should be clearly marked, so as not to confuse
 ! it with the version available from LANL.
 !
 ! The full license and distribution policy are available from
 ! https://github.com/CICE-Consortium
-! 
+!
 !=======================================================================
 
-! Main driver routine for Icepack, the column package for CICE.  
+! Main driver routine for Icepack, the column package for CICE.
 ! Initializes and steps through the model.
 !
 ! author Elizabeth C. Hunke, LANL