From e6727c13f1caef5fb8c79db46416cffb6c4668bf Mon Sep 17 00:00:00 2001
From: davidclemenssewall
 <51925832+davidclemenssewall@users.noreply.github.com>
Date: Fri, 12 May 2023 12:59:44 -0600
Subject: [PATCH 1/2] Fixes ice area inconsistency under convergence in single
 column mode (#433)

* started work on rescaling, stopped and commented code to fix divu history issue

* working advection changes, commit before deleting unneeded comments

* Fix creation of open water on ice closing

* restricted advection to ocn_data_type == SHEBA to match step_dyn_ridge

* Fixed bug in ocn_data_type if block and floating point error with small aicen

* changes requested by ECH: changed to 'uniform_ice', namelist name, code cleanup

* couple more documentation changes

* inlined uniform ice advection step

* removed icepack_step_advection_scm

* added tests advecting open water

* undoing changes to icepack_mechred

* removed unused category index variable

* updated advection step such advection divergence in consistent with opening and closing

* adjusted such that when ice_advection_type is none and the forcing is net opening, ice is exported from the cell

* minor code cleanup (per eclare108213's comments)

* changed ice_advection_type to lateral_flux_type per conversation with eclare108213

* trivial change to documentation to check RtD build
---
 columnphysics/icepack_mechred.F90             |   2 +-
 configuration/driver/icedrv_RunMod.F90        |  15 ++-
 configuration/driver/icedrv_forcing.F90       |   1 +
 configuration/driver/icedrv_init.F90          |   7 ++
 configuration/driver/icedrv_step.F90          | 112 +++++++++++++++++-
 configuration/scripts/icepack_in              |   1 +
 .../scripts/options/set_nml.fluxopenw         |   1 +
 configuration/scripts/tests/base_suite.ts     |   3 +
 doc/source/icepack_index.rst                  |   1 +
 doc/source/user_guide/ug_case_settings.rst    |   2 +
 doc/source/user_guide/ug_running.rst          |  11 ++
 11 files changed, 146 insertions(+), 10 deletions(-)
 create mode 100644 configuration/scripts/options/set_nml.fluxopenw

diff --git a/columnphysics/icepack_mechred.F90 b/columnphysics/icepack_mechred.F90
index ea049761f..314dcc7e1 100644
--- a/columnphysics/icepack_mechred.F90
+++ b/columnphysics/icepack_mechred.F90
@@ -1910,7 +1910,7 @@ subroutine icepack_step_ridge (dt,           ndtd,         &
       end subroutine icepack_step_ridge
 
 !=======================================================================
-
+      
       end module icepack_mechred
 
 !=======================================================================
diff --git a/configuration/driver/icedrv_RunMod.F90 b/configuration/driver/icedrv_RunMod.F90
index ff963673f..a51a9e72f 100644
--- a/configuration/driver/icedrv_RunMod.F90
+++ b/configuration/driver/icedrv_RunMod.F90
@@ -106,7 +106,7 @@ subroutine ice_step
       use icedrv_restart_bgc, only: write_restart_bgc
       use icedrv_step, only: prep_radiation, step_therm1, step_therm2, &
           update_state, step_dyn_ridge, step_snow, step_radiation, &
-          biogeochemistry, step_dyn_wave
+          biogeochemistry, step_dyn_wave, step_lateral_flux_scm
 
       integer (kind=int_kind) :: &
          k               ! dynamics supercycling index
@@ -177,13 +177,16 @@ subroutine ice_step
       if (tr_fsd .and. wave_spec) call step_dyn_wave(dt)
 
       do k = 1, ndtd
+         
+         ! horizontal advection of ice or open water into the single column
+         call step_lateral_flux_scm(dt_dyn)
 
-        ! ridging
-        call step_dyn_ridge (dt_dyn, ndtd)
+         ! ridging
+         call step_dyn_ridge (dt_dyn, ndtd)
 
-        ! clean up, update tendency diagnostics
-        offset = c0
-        call update_state (dt_dyn, daidtd, dvidtd, dagedtd, offset)
+         ! clean up, update tendency diagnostics
+         offset = c0
+         call update_state (dt_dyn, daidtd, dvidtd, dagedtd, offset)
 
       enddo
 
diff --git a/configuration/driver/icedrv_forcing.F90 b/configuration/driver/icedrv_forcing.F90
index d417d5eb5..4c7196aac 100644
--- a/configuration/driver/icedrv_forcing.F90
+++ b/configuration/driver/icedrv_forcing.F90
@@ -77,6 +77,7 @@ module icedrv_forcing
          atm_data_type,   & ! 'default', 'clim', 'CFS'
          ocn_data_type,   & ! 'default', 'SHEBA'
          bgc_data_type,   & ! 'default', 'ISPOL', 'NICE'
+         lateral_flux_type,   & ! 'uniform_ice', 'open_water'
          atm_data_file,   & ! atmospheric forcing data file
          ocn_data_file,   & ! ocean forcing data file
          ice_data_file,   & ! ice forcing data file
diff --git a/configuration/driver/icedrv_init.F90 b/configuration/driver/icedrv_init.F90
index 8e4958447..5a1dbac77 100644
--- a/configuration/driver/icedrv_init.F90
+++ b/configuration/driver/icedrv_init.F90
@@ -73,6 +73,7 @@ subroutine input_data
       use icedrv_forcing, only: atm_data_format, ocn_data_format, bgc_data_format
       use icedrv_forcing, only: data_dir
       use icedrv_forcing, only: oceanmixed_ice, restore_ocn, trestore
+      use icedrv_forcing, only: lateral_flux_type
 
       ! local variables
 
@@ -176,6 +177,7 @@ subroutine input_data
         default_season,  wave_spec_type,  &
         precip_units,    fyear_init,      ycycle,          &
         atm_data_type,   ocn_data_type,   bgc_data_type,   &
+        lateral_flux_type,                                &
         atm_data_file,   ocn_data_file,   bgc_data_file,   &
         ice_data_file,                                     &
         atm_data_format, ocn_data_format, bgc_data_format, &
@@ -276,6 +278,8 @@ subroutine input_data
       ocn_data_type   = 'default' ! source of ocean forcing data
       ocn_data_file   = ' '       ! ocean forcing data file
       ice_data_file   = ' '       ! ice forcing data file (opening, closing)
+      lateral_flux_type   = 'uniform_ice' ! if 'uniform_ice' assume closing 
+                                           ! fluxes in uniform ice
       bgc_data_format = 'bin'     ! file format ('bin'=binary or 'nc'=netcdf)
       bgc_data_type   = 'default' ! source of BGC forcing data
       bgc_data_file   = ' '       ! biogeochemistry forcing data file
@@ -764,6 +768,9 @@ subroutine input_data
          write(nu_diag,*)    ' bgc_data_type             = ', &
                                trim(bgc_data_type)
 
+         write(nu_diag,*)    '  lateral_flux_type       = ', &
+                               trim(lateral_flux_type)
+
          write(nu_diag,*)    ' atm_data_file             = ', &
                                trim(atm_data_file)
          write(nu_diag,*)    ' ocn_data_file             = ', &
diff --git a/configuration/driver/icedrv_step.F90 b/configuration/driver/icedrv_step.F90
index 0b5e74dda..add380006 100644
--- a/configuration/driver/icedrv_step.F90
+++ b/configuration/driver/icedrv_step.F90
@@ -6,10 +6,11 @@
 
       module icedrv_step
 
-      use icedrv_constants, only: c0, nu_diag, c4
+      use icedrv_constants, only: c0, nu_diag, c4, c1
       use icedrv_kinds
 !      use icedrv_calendar, only: istep1
       use icedrv_forcing, only: ocn_data_type
+      use icedrv_forcing, only: lateral_flux_type
       use icedrv_system, only: icedrv_system_abort
       use icepack_intfc, only: icepack_warnings_flush
       use icepack_intfc, only: icepack_warnings_aborted
@@ -23,7 +24,8 @@ module icedrv_step
 
       public :: step_therm1, step_therm2, step_dyn_ridge, &
                 prep_radiation, step_radiation, ocean_mixed_layer, &
-                update_state, biogeochemistry, step_dyn_wave, step_snow
+                update_state, biogeochemistry, step_dyn_wave, step_snow, &
+                step_lateral_flux_scm
 
 !=======================================================================
 
@@ -698,6 +700,110 @@ subroutine step_dyn_wave (dt)
 
       end subroutine step_dyn_wave
 
+
+!=======================================================================
+!
+! Run one time step of horizontal advection (if there is closing)
+!
+! authors: David Clemens-Sewall, NCAR
+
+      subroutine step_lateral_flux_scm (dt)
+
+         use icedrv_domain_size, only: ncat, nx
+         use icedrv_flux, only: closing, opening
+         use icedrv_init, only: tmask
+         use icedrv_state, only: vsnon, aicen, vicen, aice0
+      
+         real (kind=dbl_kind), intent(in) :: &
+            dt      ! time step
+   
+         ! local variables
+   
+         integer (kind=int_kind) :: &
+            i,            & ! horizontal indices
+            n               ! ice thickness category index         !
+         
+         real (kind=dbl_kind) :: &
+            expansion_ratio  ! how much the ice area will change
+                  
+         character(len=*), parameter :: subname='(step_lateral_flux_scm)'
+   
+         !-----------------------------------------------------------------
+         ! query icepack values
+         !-----------------------------------------------------------------
+            call icepack_warnings_flush(nu_diag)
+            if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
+                file=__FILE__,line= __LINE__)
+   
+         !-----------------------------------------------------------------
+         ! Ice advection
+         !-----------------------------------------------------------------
+            ! Currently we only do ridging for the SHEBA ocean data type (in step_dyn_ridge)
+            if (trim(ocn_data_type) == "SHEBA") then
+               ! Currently only uniform_ice (and none) advection is implemented
+               if (trim(lateral_flux_type) == "uniform_ice") then
+      
+                  do i = 1, nx
+         
+                  if (tmask(i)) then                     
+                     ! We assume that this single column grid cell is surrounded by
+                     ! identical ice. If so, ice closing implies the flux of
+                     ! this surrounding ice into the single column grid cell.
+                     ! Equivalently, one can think of this step as expanding the
+                     ! domain of the grid cell before the ridging step will
+                     ! contract the domain.
+                     expansion_ratio = c1 + (closing(i) - opening(i)) * dt
+                     aice0(i) = aice0(i) * expansion_ratio
+                     do n = 1, ncat
+                        ! Scale up state variables
+                        aicen(i,n) = aicen(i,n) * expansion_ratio
+                        vicen(i,n) = vicen(i,n) * expansion_ratio
+                        vsnon(i,n) = vsnon(i,n) * expansion_ratio
+                     enddo ! n
+                  endif ! tmask
+         
+                  enddo ! i
+               elseif (trim(lateral_flux_type) == "open_water") then
+                  do i = 1, nx
+         
+                     if (tmask(i)) then                     
+                        ! We assume that this single column grid cell is surrounded by
+                        ! open water. If so, net ice closing implies the flux of
+                        ! this surrounding open water into the single column grid cell.
+                        ! To accomplish this without modifying the icepack 
+                        ! columnphysics code, we do nothing at this step. Within the
+                        ! ridge_ice subroutine, icepack will ridge ice by the amount
+                        ! given in the forcing. This will drop the cell area (asum)
+                        ! below 1 and then in the second ridging iteration loop
+                        ! 'opning' will be set such that it adds enough open water
+                        ! to return the cell area to 1.
+                        ! If the forcing is net opening, we still need to flux
+                        ! ice out of the grid cell as above.
+                        expansion_ratio = c1 + (closing(i) - opening(i)) * dt
+                        if (expansion_ratio < 1) then ! net opening
+                           aice0(i) = aice0(i) * expansion_ratio
+                           do n = 1, ncat
+                              ! Remove ice from cell
+                              aicen(i,n) = aicen(i,n) * expansion_ratio
+                              vicen(i,n) = vicen(i,n) * expansion_ratio
+                              vsnon(i,n) = vsnon(i,n) * expansion_ratio
+                           enddo ! n
+                        endif ! expansion ratio < 1
+                     endif ! tmask
+            
+                     enddo ! i
+               else
+                  call icedrv_system_abort(string=subname//' ERROR: unknown lateral_flux_type: '&
+                  //trim(lateral_flux_type),file=__FILE__,line=__LINE__)
+               endif
+            endif
+
+            call icepack_warnings_flush(nu_diag)
+            if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
+                file=__FILE__, line=__LINE__)
+   
+         end subroutine step_lateral_flux_scm
+   
 !=======================================================================
 !
 ! Run one time step of ridging.
@@ -732,7 +838,7 @@ subroutine step_dyn_ridge (dt, ndtd)
       integer (kind=int_kind) :: &
          i,            & ! horizontal indices
          ntrcr,        & !
-         nbtrcr          !
+         nbtrcr
 
       character(len=*), parameter :: subname='(step_dyn_ridge)'
 
diff --git a/configuration/scripts/icepack_in b/configuration/scripts/icepack_in
index bc9e2942c..d989f2a83 100644
--- a/configuration/scripts/icepack_in
+++ b/configuration/scripts/icepack_in
@@ -121,6 +121,7 @@
     atm_data_type   = 'clim'
     ocn_data_type   = 'SHEBA'
     bgc_data_type   = 'default'
+    lateral_flux_type   = 'uniform_ice'
     fyear_init      = 2015
     ycycle          = 1
     data_dir        = '/Users/ftuser/Desktop/CICE-Consortium/ICEPACK_DATA/'
diff --git a/configuration/scripts/options/set_nml.fluxopenw b/configuration/scripts/options/set_nml.fluxopenw
new file mode 100644
index 000000000..fc2c2e088
--- /dev/null
+++ b/configuration/scripts/options/set_nml.fluxopenw
@@ -0,0 +1 @@
+lateral_flux_type = 'open_water'
diff --git a/configuration/scripts/tests/base_suite.ts b/configuration/scripts/tests/base_suite.ts
index 771381355..ce98270fd 100644
--- a/configuration/scripts/tests/base_suite.ts
+++ b/configuration/scripts/tests/base_suite.ts
@@ -19,6 +19,8 @@ smoke          col     1x1        debug,run1year,snw30percent,snwgrain
 smoke          col     1x1        debug,run1year,snwitdrdg
 smoke          col     1x1        debug,run1year,calcdragio
 smoke          col     1x1        debug,run1year,modal
+smoke          col     1x1        debug,run1year,fluxopenw
+smoke          col     1x1        debug,run1year,dyn,fluxopenw
 restart        col     1x1        debug
 restart        col     1x1        diag1
 restart        col     1x1        pondlvl
@@ -36,4 +38,5 @@ restart        col     1x1        dyn
 restart        col     1x1        fsd12,short
 restart        col     1x1        snwitdrdg,snwgrain
 restart        col     1x1        modal
+restart        col     1x1        fluxopenw
 
diff --git a/doc/source/icepack_index.rst b/doc/source/icepack_index.rst
index 08a3b0d94..ab7d582c8 100755
--- a/doc/source/icepack_index.rst
+++ b/doc/source/icepack_index.rst
@@ -224,6 +224,7 @@ either Celsius or Kelvin units).  Deprecated parameters are listed at the end.
    "Hstar", "determines mean thickness of ridged ice", "25. m"
    "**I**", "", ""
    "i0vis","fraction of penetrating visible solar radiation", "0.70"
+   "lateral_flux_type", ":math:`\bullet` laterally flux ice or open water into grid cell when closing", ""
    "ice_ic", ":math:`\bullet` choice of initial conditions", ""
    "ice_stdout", "unit number for standard output", ""
    "ice_stderr", "unit number for standard error output", ""
diff --git a/doc/source/user_guide/ug_case_settings.rst b/doc/source/user_guide/ug_case_settings.rst
index 6d4f877ca..5fcc188d9 100644
--- a/doc/source/user_guide/ug_case_settings.rst
+++ b/doc/source/user_guide/ug_case_settings.rst
@@ -345,6 +345,8 @@ forcing_nml
    "``formdrag``", "logical", "calculate form drag", "``.false.``"
    "``fyear_init``", "integer", "first year of atmospheric forcing data", "1998"
    "``highfreq``", "logical", "high-frequency atmo coupling", "``.false.``"
+   "``lateral_flux_type``", "``uniform_ice``", "flux ice with identical properties into the cell when closing (Icepack only)", ""
+   "", "``none``", "advect open water into the cell when closing (Icepack only)", ""
    "``ice_data_file``", "string", "file containing ice opening, closing data", "' '"
    "``l_mpond_fresh``", "``.false.``", "release pond water immediately to ocean", "``.false.``"
    "", "``true``", "retain (topo) pond water until ponds drain", ""
diff --git a/doc/source/user_guide/ug_running.rst b/doc/source/user_guide/ug_running.rst
index d30045d66..782ef56b8 100755
--- a/doc/source/user_guide/ug_running.rst
+++ b/doc/source/user_guide/ug_running.rst
@@ -747,6 +747,17 @@ cases.  Current filenames can be found in the options scripts in
    about 4 :math:`W/m^2` for biases in the forcings or the model. In the namelist, set ``atm_data_type = clim`` 
    to use climatological atmospheric forcing.
 
+Horizontal ice advection
+------------------------
+
+When Icepack is run in standalone mode with a dynamical forcing (e.g., ``ocn_data_type = SHEBA``),
+closing implies the lateral flux of ice or open water area into the grid cell. The default assumption
+(in the namelist, ``lateral_flux_type = 'uniform_ice'``) is that the active grid cell is surrounded by grid cells 
+with identical ice properties to the active grid cell, i.e. the ice is uniform across all of those cells, and 
+when the dynamical forcing is net convergence, this uniform ice is fluxed into the grid cell. Alternatively, one may
+assume that the active grid cell is surrounded by open water (in the namelist ``lateral_flux_type = 'open_water'``), in which case closing (i.e., ice convergence) will produce open water in 
+the grid cell. In either case, when the forcing is net divergence, ice area and volume are removed from the grid 
+cell to accommodate the formation of open water implied by the net divergence.
 
 Run Directories
 ---------------

From 6703bc533c96802235e2f20de5fffc0bc6cc4c97 Mon Sep 17 00:00:00 2001
From: Tony Craig <apcraig@users.noreply.github.com>
Date: Mon, 22 May 2023 04:55:36 -0700
Subject: [PATCH 2/2] Fix RTD, minor updates to documentation (#437)

* - Add .readthedocs.yaml to fix readthedocs
- Minor update to icepack_mushy_physics inline documentation
- Update Icepack interfaces documentation, add icepack_mushy public interfaces
- Remove trailing blank spaces in source code

* Add cpp CCSMCOUPLED to trigger CESMCOUPLED cpp in icepack_kinds.
icepack_kinds is used by all files in columnphysics, so should be comprehensive.
---
 .readthedocs.yaml                        | 29 +++++++++
 columnphysics/icepack_kinds.F90          |  3 +
 columnphysics/icepack_mechred.F90        |  2 +-
 columnphysics/icepack_mushy_physics.F90  | 68 ++++++++++----------
 columnphysics/icepack_therm_itd.F90      | 10 +--
 configuration/driver/icedrv_RunMod.F90   |  2 +-
 configuration/driver/icedrv_init.F90     |  2 +-
 configuration/driver/icedrv_restart.F90  | 62 +++++++++----------
 configuration/driver/icedrv_step.F90     | 34 +++++-----
 doc/source/user_guide/interfaces.include | 79 ++++++++++++++++++++++--
 10 files changed, 193 insertions(+), 98 deletions(-)
 create mode 100644 .readthedocs.yaml

diff --git a/.readthedocs.yaml b/.readthedocs.yaml
new file mode 100644
index 000000000..f83760cce
--- /dev/null
+++ b/.readthedocs.yaml
@@ -0,0 +1,29 @@
+# .readthedocs.yaml
+# Read the Docs configuration file
+# See https://docs.readthedocs.io/en/stable/config-file/v2.html for details
+
+# Required
+version: 2
+
+# Set the version of Python and other tools you might need
+build:
+  os: ubuntu-22.04
+  tools:
+    python: "3.7"
+    # You can also specify other tool versions:
+    # nodejs: "19"
+    # rust: "1.64"
+    # golang: "1.19"
+
+# Build documentation in the docs/ directory with Sphinx
+sphinx:
+   configuration: doc/source/conf.py
+
+# If using Sphinx, optionally build your docs in additional formats such as PDF
+# formats:
+#    - pdf
+
+# Optionally declare the Python requirements required to build your docs
+python:
+   install:
+   - requirements: doc/requirements.txt
diff --git a/columnphysics/icepack_kinds.F90 b/columnphysics/icepack_kinds.F90
index 664643fa3..832f45466 100644
--- a/columnphysics/icepack_kinds.F90
+++ b/columnphysics/icepack_kinds.F90
@@ -1,3 +1,6 @@
+#ifdef CCSMCOUPLED
+#define CESMCOUPLED
+#endif
 !=======================================================================
 
 ! Defines variable precision for all common data types
diff --git a/columnphysics/icepack_mechred.F90 b/columnphysics/icepack_mechred.F90
index 314dcc7e1..ea049761f 100644
--- a/columnphysics/icepack_mechred.F90
+++ b/columnphysics/icepack_mechred.F90
@@ -1910,7 +1910,7 @@ subroutine icepack_step_ridge (dt,           ndtd,         &
       end subroutine icepack_step_ridge
 
 !=======================================================================
-      
+
       end module icepack_mechred
 
 !=======================================================================
diff --git a/columnphysics/icepack_mushy_physics.F90 b/columnphysics/icepack_mushy_physics.F90
index d100ebd2e..466c5d727 100644
--- a/columnphysics/icepack_mushy_physics.F90
+++ b/columnphysics/icepack_mushy_physics.F90
@@ -68,11 +68,10 @@ module icepack_mushy_physics
 !=======================================================================
 ! Physical Quantities
 !=======================================================================
+! Detemine the conductivity of the mush from enthalpy and salinity
 
   subroutine conductivity_mush_array(nilyr, zqin, zSin, km)
 
-    ! detemine the conductivity of the mush from enthalpy and salinity
-
     integer (kind=int_kind), intent(in) :: &
          nilyr   ! number of ice layers
 
@@ -101,17 +100,19 @@ subroutine conductivity_mush_array(nilyr, zqin, zSin, km)
   end subroutine conductivity_mush_array
 
 !=======================================================================
+!autodocument_start icepack_mushy_density_brine
+! Compute density of brine from brine salinity
 
   function icepack_mushy_density_brine(Sbr) result(rho)
 
-    ! density of brine from brine salinity
-
     real(kind=dbl_kind), intent(in) :: &
          Sbr ! brine salinity (ppt)
 
     real(kind=dbl_kind) :: &
          rho ! brine density (kg m-3)
 
+!autodocument_end
+
     real(kind=dbl_kind), parameter :: &
          a = 1000.3_dbl_kind    , & ! zeroth empirical coefficient
          b = 0.78237_dbl_kind   , & ! linear empirical coefficient
@@ -126,11 +127,10 @@ end function icepack_mushy_density_brine
 !=======================================================================
 ! Snow
 !=======================================================================
+! Heat conductivity of the snow
 
   subroutine conductivity_snow_array(ks)
 
-    ! heat conductivity of the snow
-
     real(kind=dbl_kind), dimension(:), intent(out) :: &
          ks ! snow layer conductivity (W m-1 K-1)
 
@@ -141,11 +141,10 @@ subroutine conductivity_snow_array(ks)
   end subroutine conductivity_snow_array
 
 !=======================================================================
+! Enthalpy of snow from snow temperature
 
   function enthalpy_snow(zTsn) result(zqsn)
 
-    ! enthalpy of snow from snow temperature
-
     real(kind=dbl_kind), intent(in) :: &
          zTsn ! snow layer temperature (C)
 
@@ -159,11 +158,10 @@ function enthalpy_snow(zTsn) result(zqsn)
   end function enthalpy_snow
 
 !=======================================================================
+! Temperature of snow from the snow enthalpy
 
   function temperature_snow(zqsn) result(zTsn)
 
-    ! temperature of snow from the snow enthalpy
-
     real(kind=dbl_kind), intent(in) :: &
          zqsn ! snow layer enthalpy (J m-3)
 
@@ -182,12 +180,11 @@ end function temperature_snow
 !=======================================================================
 ! Mushy Layer Formulation - Assur (1958) liquidus
 !=======================================================================
+! Liquidus relation: equilibrium brine salinity as function of temperature
+! based on empirical data from Assur (1958)
 
   function liquidus_brine_salinity_mush(zTin) result(Sbr)
 
-    ! liquidus relation: equilibrium brine salinity as function of temperature
-    ! based on empirical data from Assur (1958)
-
     real(kind=dbl_kind), intent(in) :: &
          zTin         ! ice layer temperature (C)
 
@@ -223,12 +220,11 @@ function liquidus_brine_salinity_mush(zTin) result(Sbr)
   end function liquidus_brine_salinity_mush
 
 !=======================================================================
+! Liquidus relation: equilibrium temperature as function of brine salinity
+! based on empirical data from Assur (1958)
 
   function liquidus_temperature_mush(Sbr) result(zTin)
 
-    ! liquidus relation: equilibrium temperature as function of brine salinity
-    ! based on empirical data from Assur (1958)
-
     real(kind=dbl_kind), intent(in) :: &
          Sbr    ! ice brine salinity (ppt)
 
@@ -264,11 +260,10 @@ function liquidus_temperature_mush(Sbr) result(zTin)
   end function liquidus_temperature_mush
 
 !=======================================================================
+! Enthalpy of mush from mush temperature and bulk salinity
 
   function enthalpy_mush(zTin, zSin) result(zqin)
 
-    ! enthalpy of mush from mush temperature and bulk salinity
-
     real(kind=dbl_kind), intent(in) :: &
          zTin, & ! ice layer temperature (C)
          zSin    ! ice layer bulk salinity (ppt)
@@ -289,11 +284,10 @@ function enthalpy_mush(zTin, zSin) result(zqin)
   end function enthalpy_mush
 
 !=======================================================================
+! Enthalpy of mush from mush temperature and liquid fraction
 
   function enthalpy_mush_liquid_fraction(zTin, phi) result(zqin)
 
-    ! enthalpy of mush from mush temperature and bulk salinity
-
     real(kind=dbl_kind), intent(in) :: &
          zTin, & ! ice layer temperature (C)
          phi     ! liquid fraction
@@ -309,12 +303,11 @@ function enthalpy_mush_liquid_fraction(zTin, phi) result(zqin)
   end function enthalpy_mush_liquid_fraction
 
 !=======================================================================
+! Enthalpy of melting of mush from bulk salinity
+! Energy needed to fully melt mush (T < 0)
 
   function enthalpy_of_melting(zSin) result(qm)
 
-    ! enthalpy of melting of mush
-    ! energy needed to fully melt mush (T < 0)
-
     real(kind=dbl_kind), intent(in) :: &
          zSin ! ice layer bulk salinity (ppt)
 
@@ -328,11 +321,10 @@ function enthalpy_of_melting(zSin) result(qm)
   end function enthalpy_of_melting
 
 !=======================================================================
+! Enthalpy of brine (fully liquid) from temperature
 
   function enthalpy_brine(zTin) result(qbr)
 
-    ! enthalpy of brine (fully liquid)
-
     real(kind=dbl_kind), intent(in) :: &
          zTin ! ice layer temperature (C)
 
@@ -346,11 +338,11 @@ function enthalpy_brine(zTin) result(qbr)
   end function enthalpy_brine
 
 !=======================================================================
+!autodocument_start icepack_mushy_temperature_mush
+! Temperature of mush from mush enthalpy and bulk salinity
 
   function icepack_mushy_temperature_mush(zqin, zSin) result(zTin)
 
-    ! temperature of mush from mush enthalpy
-
     real(kind=dbl_kind), intent(in) :: &
          zqin   , & ! ice enthalpy (J m-3)
          zSin       ! ice layer bulk salinity (ppt)
@@ -358,6 +350,8 @@ function icepack_mushy_temperature_mush(zqin, zSin) result(zTin)
     real(kind=dbl_kind) :: &
          zTin       ! ice layer temperature (C)
 
+!autodocument_end
+
     real(kind=dbl_kind) :: &
          qb     , & ! liquidus break enthalpy
          q0     , & ! fully melted enthalpy
@@ -454,11 +448,10 @@ function icepack_mushy_temperature_mush(zqin, zSin) result(zTin)
   end function icepack_mushy_temperature_mush
 
 !=======================================================================
+! Temperature of mush from mush enthalpy and liquid fraction
 
   function temperature_mush_liquid_fraction(zqin, phi) result(zTin)
 
-    ! temperature of mush from mush enthalpy
-
     real(kind=dbl_kind), intent(in) :: &
          zqin   , & ! ice enthalpy (J m-3)
          phi        ! liquid fraction
@@ -474,11 +467,10 @@ function temperature_mush_liquid_fraction(zqin, phi) result(zTin)
   end function temperature_mush_liquid_fraction
 
 !=======================================================================
+! Mush heat conductivity from mush temperature and bulk salinity
 
   function heat_conductivity(zTin, zSin) result(km)
 
-    ! msuh heat conductivity from mush temperature and bulk salinity
-
     real(kind=dbl_kind), intent(in) :: &
          zTin              , & ! ice layer temperature (C)
          zSin                  ! ice layer bulk salinity (ppt)
@@ -497,18 +489,22 @@ function heat_conductivity(zTin, zSin) result(km)
 
   end function heat_conductivity
 
-  !=======================================================================
+!=======================================================================
+!autodocument_start icepack_mushy_liquid_fraction
+! Liquid fraction of mush from mush temperature and bulk salinity
 
   function icepack_mushy_liquid_fraction(zTin, zSin) result(phi)
 
-    ! liquid fraction of mush from mush temperature and bulk salinity
-
     real(kind=dbl_kind), intent(in) :: &
          zTin, & ! ice layer temperature (C)
          zSin    ! ice layer bulk salinity (ppt)
 
     real(kind=dbl_kind) :: &
-         phi , & ! liquid fraction
+         phi     ! liquid fraction
+
+!autodocument_end
+
+    real(kind=dbl_kind) :: &
          Sbr     ! brine salinity (ppt)
 
     character(len=*),parameter :: subname='(icepack_mushy_liquid_fraction)'
@@ -521,5 +517,3 @@ end function icepack_mushy_liquid_fraction
 !=======================================================================
 
 end module icepack_mushy_physics
-
-
diff --git a/columnphysics/icepack_therm_itd.F90 b/columnphysics/icepack_therm_itd.F90
index 3aac1ac70..c16e38fe7 100644
--- a/columnphysics/icepack_therm_itd.F90
+++ b/columnphysics/icepack_therm_itd.F90
@@ -915,7 +915,7 @@ subroutine lateral_melt (dt,         ncat,       &
       real (kind=dbl_kind), intent(in) :: &
          rside   , & ! fraction of ice that melts laterally
          wlat        ! lateral melt rate (m/s)
-         
+
       real (kind=dbl_kind), intent(inout) :: &
          fside       ! lateral heat flux (W/m^2)
 
@@ -982,7 +982,7 @@ subroutine lateral_melt (dt,         ncat,       &
          f_flx         !
 
       real (kind=dbl_kind) :: &
-         sicen,        & 
+         sicen,        &
          etot,         & ! column energy per itd cat, for FSD code
          elapsed_t,    & ! FSD subcycling
          subdt           ! FSD timestep (s)
@@ -1064,8 +1064,8 @@ subroutine lateral_melt (dt,         ncat,       &
             do k = 1, nilyr
                etot = etot + trcrn(nt_qice+k-1,n) * vicen(n)/real(nilyr,kind=dbl_kind)
             enddo                  ! nilyr
-               
-            ! lateral heat flux, fside < 0        
+
+            ! lateral heat flux, fside < 0
             fside = fside + rsiden(n)*etot/dt
 
          enddo ! ncat
@@ -1245,7 +1245,7 @@ subroutine lateral_melt (dt,         ncat,       &
       if (tr_fsd) then
 
          trcrn(nt_fsd:nt_fsd+nfsd-1,:) =  afsdn
-         
+
          call icepack_cleanup_fsd (ncat, nfsd, trcrn(nt_fsd:nt_fsd+nfsd-1,:) )
          if (icepack_warnings_aborted(subname)) return
 
diff --git a/configuration/driver/icedrv_RunMod.F90 b/configuration/driver/icedrv_RunMod.F90
index a51a9e72f..db97c79c6 100644
--- a/configuration/driver/icedrv_RunMod.F90
+++ b/configuration/driver/icedrv_RunMod.F90
@@ -177,7 +177,7 @@ subroutine ice_step
       if (tr_fsd .and. wave_spec) call step_dyn_wave(dt)
 
       do k = 1, ndtd
-         
+
          ! horizontal advection of ice or open water into the single column
          call step_lateral_flux_scm(dt_dyn)
 
diff --git a/configuration/driver/icedrv_init.F90 b/configuration/driver/icedrv_init.F90
index 5a1dbac77..62b852565 100644
--- a/configuration/driver/icedrv_init.F90
+++ b/configuration/driver/icedrv_init.F90
@@ -278,7 +278,7 @@ subroutine input_data
       ocn_data_type   = 'default' ! source of ocean forcing data
       ocn_data_file   = ' '       ! ocean forcing data file
       ice_data_file   = ' '       ! ice forcing data file (opening, closing)
-      lateral_flux_type   = 'uniform_ice' ! if 'uniform_ice' assume closing 
+      lateral_flux_type   = 'uniform_ice' ! if 'uniform_ice' assume closing
                                            ! fluxes in uniform ice
       bgc_data_format = 'bin'     ! file format ('bin'=binary or 'nc'=netcdf)
       bgc_data_type   = 'default' ! source of BGC forcing data
diff --git a/configuration/driver/icedrv_restart.F90 b/configuration/driver/icedrv_restart.F90
index 745ec9210..7202da455 100644
--- a/configuration/driver/icedrv_restart.F90
+++ b/configuration/driver/icedrv_restart.F90
@@ -30,9 +30,9 @@ module icedrv_restart
           write_restart_fsd,       read_restart_fsd, &
           write_restart_iso,       read_restart_iso, &
           write_restart_aero,      read_restart_aero
-      
+
       character (len=3), private :: nchar !
-      
+
       logical (kind=log_kind), private :: diag ! netCDF diagnostics flag
 
 #ifdef USE_NETCDF
@@ -59,7 +59,7 @@ module icedrv_restart
 !=======================================================================
 !---! these subroutines write/read restart data files ..
 !     Which can either be in unformatted Fortran format
-!     (restart_format = 'bin') or NetCDF (restart_format = 'nc')      
+!     (restart_format = 'bin') or NetCDF (restart_format = 'nc')
 !=======================================================================
 
 ! Dumps all values needed for a restart
@@ -96,7 +96,7 @@ subroutine dumpfile
 
       character(len=char_len_long) :: filename
       character(len=*), parameter :: subname='(dumpfile)'
-      
+
 #ifdef USE_NETCDF
       ! local variables if we're writing in NetCDF
       integer (kind=int_kind) :: &
@@ -108,7 +108,7 @@ subroutine dumpfile
 
       ! get year
       iyear = nyr + year_init - 1
-      
+
       ! query which tracers are active and their indices
       call icepack_query_tracer_indices(nt_Tsfc_out=nt_Tsfc, nt_sice_out=nt_sice, &
          nt_qice_out=nt_qice, nt_qsno_out=nt_qsno)
@@ -228,14 +228,14 @@ subroutine dumpfile
    ! called in icedrv_RunMod.F90 to prevent circular dependencies
    !      if (solve_zsal .or. skl_bgc .or. z_tracers) &
    !                        call write_restart_bgc         ! biogeochemistry
-      
+
       ! Close netcdf file
       if (restart_format == 'nc') then
 #ifdef USE_NETCDF
          status = nf90_close(ncid)
          if (status /= nf90_noerr) call icedrv_system_abort(string=subname, &
             file=__FILE__,line= __LINE__)
-#endif         
+#endif
       endif
 
       end subroutine dumpfile
@@ -284,7 +284,7 @@ subroutine restartfile (ice_ic)
       ! local variable for reading from a netcdf file
       integer (kind=int_kind) :: &
          status
-      
+
       ! Query tracers
       call icepack_query_tracer_indices(nt_Tsfc_out=nt_Tsfc, nt_sice_out=nt_sice, &
           nt_qice_out=nt_qice, nt_qsno_out=nt_qsno)
@@ -455,7 +455,7 @@ subroutine read_restart_field(nu,work,ndim,name)
 
       real (kind=dbl_kind), dimension(nx,ndim), intent(inout) :: &
          work              ! input array (real, 8-byte)
-      
+
       character (len=*), intent(in), optional :: &
          name
 
@@ -514,10 +514,10 @@ subroutine write_restart_field(nu,work,ndim,vname,dims)
 
       real (kind=dbl_kind), dimension(nx,ndim), intent(in) :: &
          work              ! input array (real, 8-byte)
-      
+
       character (len=*), intent(in), optional :: &
          vname             ! variable name for netcdf output
-      
+
       integer (kind=int_kind), intent(in), optional :: &
          dims(:)           ! netcdf dimension IDs
 
@@ -540,7 +540,7 @@ subroutine write_restart_field(nu,work,ndim,vname,dims)
          status, &          ! Status variable from netCDF routine
          varid              ! Variable ID
 #endif
-      
+
       if (restart_format == 'bin') then
          do n = 1, ndim
          work2(:) = work(:,n)
@@ -724,10 +724,10 @@ end subroutine read_restart_snow
 ! author Elizabeth C. Hunke, LANL
 
       subroutine write_restart_age(dims)
-      
+
       use icedrv_state, only: trcrn
       use icedrv_domain_size, only: ncat
-      
+
       integer (kind=int_kind), intent(in), optional :: &
          dims(:)           ! netcdf dimension IDs
 
@@ -761,7 +761,7 @@ subroutine read_restart_age()
       call icepack_warnings_flush(nu_diag)
       if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
           file=__FILE__,line= __LINE__)
-      
+
       call read_restart_field(nu_restart,trcrn(:,nt_iage,:),ncat,'iage')
 
       end subroutine read_restart_age
@@ -775,7 +775,7 @@ subroutine write_restart_fsd(dims)
 
       use icedrv_state, only: trcrn
       use icedrv_domain_size, only: ncat, nfsd
-      
+
       integer (kind=int_kind), intent(in), optional :: &
          dims(:)           ! netcdf dimension IDs
 
@@ -831,10 +831,10 @@ subroutine write_restart_FY(dims)
       use icedrv_flux, only: frz_onset
       use icedrv_state, only: trcrn
       use icedrv_domain_size, only: ncat
-      
+
       integer (kind=int_kind), intent(in), optional :: &
          dims(:)           ! netcdf dimension IDs
-      
+
       integer (kind=int_kind) :: nt_FY
       character(len=*), parameter :: subname='(write_restart_FY)'
 
@@ -940,7 +940,7 @@ subroutine write_restart_pond_lvl(dims)
       use icedrv_flux, only: fsnow
       use icedrv_state, only: trcrn
       use icedrv_domain_size, only: ncat
-      
+
       integer (kind=int_kind), intent(in), optional :: &
          dims(:)           ! netcdf dimension IDs
 
@@ -1007,7 +1007,7 @@ subroutine write_restart_aero(dims)
       use icedrv_domain_size, only: n_aero
       use icedrv_state, only: trcrn
       use icedrv_domain_size, only: ncat
-      
+
       integer (kind=int_kind), intent(in), optional :: &
          dims(:)           ! netcdf dimension IDs
 
@@ -1092,7 +1092,7 @@ subroutine write_restart_iso(dims)
       use icedrv_domain_size, only: n_iso
       use icedrv_state, only: trcrn
       use icedrv_domain_size, only: ncat
-      
+
       integer (kind=int_kind), intent(in), optional :: &
          dims(:)           ! netcdf dimension IDs
 
@@ -1331,7 +1331,7 @@ subroutine write_restart_field_nc1D(ncid,nrec,work,atype,vname,ndim3,diag)
 
       status = nf90_inq_varid(ncid,trim(vname),varid)
       if (status /= 0) then
-         write(nu_diag,*) 'Writing out ',trim(vname) 
+         write(nu_diag,*) 'Writing out ',trim(vname)
          write(nu_diag,*) 'Erros Status ',status
          call icedrv_system_abort(string='Write out restart', &
             file=__FILE__,line= __LINE__)
@@ -1346,7 +1346,7 @@ subroutine ice_write_nc2D(fid,  nrec,  varid, work,  diag)
 
 ! Write a 2D field (nx, ncat) in NetCDF restart
 ! author Chris Riedel, NCAR
-   
+
       use icedrv_domain_size, only: ncat,nx
 
       integer (kind=int_kind), intent(in) :: &
@@ -1368,7 +1368,7 @@ subroutine ice_write_nc2D(fid,  nrec,  varid, work,  diag)
          id,              & ! dimension index
          dimlen             ! size of dimension
       integer (kind=int_kind) :: start2(2),count2(2)
-      
+
       real (kind=dbl_kind) :: &
          amin, amax, asum   ! min, max values and sum of input array
 
@@ -1382,7 +1382,7 @@ subroutine ice_write_nc2D(fid,  nrec,  varid, work,  diag)
       count2(1) = nx
       start2(2) = 1
       count2(2) = ncat
-   
+
       status = nf90_put_var( fid, varid, work, &
                start=start2, &
                count=count2)
@@ -1462,7 +1462,7 @@ end subroutine ice_write_nc1D
 !=======================================================================
       subroutine read_restart_field_net2D(nu,nrec,work,vname,ndim3, &
                      diag)
-      
+
 ! author Chris Riedel, NCAR
 
       use icedrv_domain_size, only: ncat,nx
@@ -1493,7 +1493,7 @@ end subroutine read_restart_field_net2D
 !=======================================================================
       subroutine read_restart_field_net1D(nu,nrec,work,vname,ndim3, &
                   diag)
-            
+
 ! author Chris Riedel, NCAR
 
       use icedrv_domain_size, only: ncat,nx
@@ -1519,11 +1519,11 @@ subroutine read_restart_field_net1D(nu,nrec,work,vname,ndim3, &
       status        ! status variable from netCDF routine
 
       call ice_read_nc1D(ncid, 1, vname, work, diag)
-      
+
       end subroutine read_restart_field_net1D
 !=======================================================================
       subroutine ice_read_nc2D(fid,  nrec,  varname, work,  diag)
-      
+
 ! author Chris Riedel, NCAR
 
       use icedrv_domain_size, only: ncat,nx
@@ -1582,11 +1582,11 @@ subroutine ice_read_nc2D(fid,  nrec,  varname, work,  diag)
          asum = sum   (work)
          write(nu_diag,*) ' min, max, sum =', amin, amax, asum
       endif
-   
+
       end subroutine ice_read_nc2D
 !=======================================================================
       subroutine ice_read_nc1D(fid,  nrec,  varname, work,  diag)
-      
+
 ! author Chris Riedel, NCAR
 
       use icedrv_domain_size, only: ncat,nx
diff --git a/configuration/driver/icedrv_step.F90 b/configuration/driver/icedrv_step.F90
index add380006..f21bc1677 100644
--- a/configuration/driver/icedrv_step.F90
+++ b/configuration/driver/icedrv_step.F90
@@ -713,28 +713,28 @@ subroutine step_lateral_flux_scm (dt)
          use icedrv_flux, only: closing, opening
          use icedrv_init, only: tmask
          use icedrv_state, only: vsnon, aicen, vicen, aice0
-      
+
          real (kind=dbl_kind), intent(in) :: &
             dt      ! time step
-   
+
          ! local variables
-   
+
          integer (kind=int_kind) :: &
             i,            & ! horizontal indices
             n               ! ice thickness category index         !
-         
+
          real (kind=dbl_kind) :: &
             expansion_ratio  ! how much the ice area will change
-                  
+
          character(len=*), parameter :: subname='(step_lateral_flux_scm)'
-   
+
          !-----------------------------------------------------------------
          ! query icepack values
          !-----------------------------------------------------------------
             call icepack_warnings_flush(nu_diag)
             if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
                 file=__FILE__,line= __LINE__)
-   
+
          !-----------------------------------------------------------------
          ! Ice advection
          !-----------------------------------------------------------------
@@ -742,10 +742,10 @@ subroutine step_lateral_flux_scm (dt)
             if (trim(ocn_data_type) == "SHEBA") then
                ! Currently only uniform_ice (and none) advection is implemented
                if (trim(lateral_flux_type) == "uniform_ice") then
-      
+
                   do i = 1, nx
-         
-                  if (tmask(i)) then                     
+
+                  if (tmask(i)) then
                      ! We assume that this single column grid cell is surrounded by
                      ! identical ice. If so, ice closing implies the flux of
                      ! this surrounding ice into the single column grid cell.
@@ -761,16 +761,16 @@ subroutine step_lateral_flux_scm (dt)
                         vsnon(i,n) = vsnon(i,n) * expansion_ratio
                      enddo ! n
                   endif ! tmask
-         
+
                   enddo ! i
                elseif (trim(lateral_flux_type) == "open_water") then
                   do i = 1, nx
-         
-                     if (tmask(i)) then                     
+
+                     if (tmask(i)) then
                         ! We assume that this single column grid cell is surrounded by
                         ! open water. If so, net ice closing implies the flux of
                         ! this surrounding open water into the single column grid cell.
-                        ! To accomplish this without modifying the icepack 
+                        ! To accomplish this without modifying the icepack
                         ! columnphysics code, we do nothing at this step. Within the
                         ! ridge_ice subroutine, icepack will ridge ice by the amount
                         ! given in the forcing. This will drop the cell area (asum)
@@ -790,7 +790,7 @@ subroutine step_lateral_flux_scm (dt)
                            enddo ! n
                         endif ! expansion ratio < 1
                      endif ! tmask
-            
+
                      enddo ! i
                else
                   call icedrv_system_abort(string=subname//' ERROR: unknown lateral_flux_type: '&
@@ -801,9 +801,9 @@ subroutine step_lateral_flux_scm (dt)
             call icepack_warnings_flush(nu_diag)
             if (icepack_warnings_aborted()) call icedrv_system_abort(string=subname, &
                 file=__FILE__, line=__LINE__)
-   
+
          end subroutine step_lateral_flux_scm
-   
+
 !=======================================================================
 !
 ! Run one time step of ridging.
diff --git a/doc/source/user_guide/interfaces.include b/doc/source/user_guide/interfaces.include
index 95db22ba6..c45625551 100644
--- a/doc/source/user_guide/interfaces.include
+++ b/doc/source/user_guide/interfaces.include
@@ -601,6 +601,65 @@ icepack_step_ridge
   
  
 
+icepack_mushy_physics.F90
+~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+
+.. _icepack_mushy_density_brine:
+
+icepack_mushy_density_brine
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+.. code-block:: fortran
+
+  ! Compute density of brine from brine salinity
+  
+    function icepack_mushy_density_brine(Sbr) result(rho)
+  
+      real(kind=dbl_kind), intent(in) :: &
+           Sbr ! brine salinity (ppt)
+  
+      real(kind=dbl_kind) :: &
+           rho ! brine density (kg m-3)
+  
+ 
+
+.. _icepack_mushy_temperature_mush:
+
+icepack_mushy_temperature_mush
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+.. code-block:: fortran
+
+  ! Temperature of mush from mush enthalpy and bulk salinity
+  
+    function icepack_mushy_temperature_mush(zqin, zSin) result(zTin)
+  
+      real(kind=dbl_kind), intent(in) :: &
+           zqin   , & ! ice enthalpy (J m-3)
+           zSin       ! ice layer bulk salinity (ppt)
+  
+      real(kind=dbl_kind) :: &
+           zTin       ! ice layer temperature (C)
+  
+ 
+
+.. _icepack_mushy_liquid_fraction:
+
+icepack_mushy_liquid_fraction
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+.. code-block:: fortran
+
+  ! Liquid fraction of mush from mush temperature and bulk salinity
+  
+    function icepack_mushy_liquid_fraction(zTin, zSin) result(phi)
+  
+      real(kind=dbl_kind), intent(in) :: &
+           zTin, & ! ice layer temperature (C)
+           zSin    ! ice layer bulk salinity (ppt)
+  
+      real(kind=dbl_kind) :: &
+           phi     ! liquid fraction
+  
+ 
+
 icepack_ocean.F90
 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 
@@ -1908,7 +1967,7 @@ icepack_step_therm2
                                        Tf,           sss,           &
                                        salinz,                      &
                                        rside,        meltl,         &
-                                       fside,                       &
+                                       fside,        wlat,          &
                                        frzmlt,       frazil,        &
                                        frain,        fpond,         &
                                        fresh,        fsalt,         &
@@ -1948,10 +2007,12 @@ icepack_step_therm2
            Tf       , & ! freezing temperature (C)
            sss      , & ! sea surface salinity (ppt)
            rside    , & ! fraction of ice that melts laterally
-           fside    , & ! lateral heat flux (W/m^2)
            frzmlt   , & ! freezing/melting potential (W/m^2)
            wave_sig_ht ! significant height of waves in ice (m)
   
+        real (kind=dbl_kind), intent(in), optional :: &
+           wlat         ! lateral melt rate (m/s)
+  
         real (kind=dbl_kind), dimension(:), intent(in)  :: &
            wave_spectrum  ! ocean surface wave spectrum E(f) (m^2 s)
   
@@ -1990,6 +2051,7 @@ icepack_step_therm2
         real (kind=dbl_kind), intent(inout) :: &
            aice     , & ! sea ice concentration
            aice0    , & ! concentration of open water
+           fside    , & ! lateral heat flux (W/m^2)
            frain    , & ! rainfall rate (kg/m^2 s)
            fpond    , & ! fresh water flux to ponds (kg/m^2/s)
            fresh    , & ! fresh water flux to ocean (kg/m^2/s)
@@ -2226,7 +2288,7 @@ icepack_step_therm1
                                       fbot        ,               &
                                       Tbot        , Tsnice      , &
                                       frzmlt      , rside       , &
-                                      fside       ,               &
+                                      fside       , wlat        , &
                                       fsnow       , frain       , &
                                       fpond       , fsloss      , &
                                       fsurf       , fsurfn      , &
@@ -2367,21 +2429,28 @@ icepack_step_therm1
            mlt_onset   , & ! day of year that sfc melting begins
            frz_onset       ! day of year that freezing begins (congel or frazil)
   
+        real (kind=dbl_kind), intent(out), optional :: &
+           wlat            ! lateral melt rate                    (m/s)
+  
         real (kind=dbl_kind), intent(inout), optional :: &
            fswthru_vdr , & ! vis dir shortwave penetrating to ocean (W/m^2)
            fswthru_vdf , & ! vis dif shortwave penetrating to ocean (W/m^2)
            fswthru_idr , & ! nir dir shortwave penetrating to ocean (W/m^2)
            fswthru_idf , & ! nir dif shortwave penetrating to ocean (W/m^2)
            dsnow       , & ! change in snow depth     (m/step-->cm/day)
-           meltsliq    , & ! mass of snow melt                 (kg/m^2)
            fsloss          ! rate of snow loss to leads      (kg/m^2/s)
   
+        real (kind=dbl_kind), intent(out), optional :: &
+           meltsliq        ! mass of snow melt                 (kg/m^2)
+  
         real (kind=dbl_kind), dimension(:), intent(inout), optional :: &
            Qa_iso      , & ! isotope specific humidity          (kg/kg)
            Qref_iso    , & ! isotope 2m atm ref spec humidity   (kg/kg)
            fiso_atm    , & ! isotope deposition rate         (kg/m^2 s)
            fiso_ocn    , & ! isotope flux to ocean           (kg/m^2/s)
-           fiso_evap   , & ! isotope evaporation             (kg/m^2/s)
+           fiso_evap       ! isotope evaporation             (kg/m^2/s)
+  
+        real (kind=dbl_kind), dimension(:), intent(inout), optional :: &
            meltsliqn       ! mass of snow melt                 (kg/m^2)
   
         real (kind=dbl_kind), dimension(:,:), intent(inout), optional :: &