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

@@ -1910,7 +1910,7 @@ subroutine icepack_step_ridge (dt,           ndtd,         &
       end subroutine icepack_step_ridge
 
 !=======================================================================
-
+      
       end module icepack_mechred
 
 !=======================================================================

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
 

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

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             = ', &

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)'
 

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/'

configuration/scripts/options/set_nml.fluxopenw

@@ -0,0 +1 @@
+lateral_flux_type = 'open_water'

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
 

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", ""

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", ""

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
 ---------------