Update halo and interior restoring (CICE-Consortium#1110)

This is a significant refactor of the restoring implementation.  Separated the exterior halo boundary update and the interior restoring capability.  Added namelist set_boundary_flds to set exterior halo values.  Added namelist restore_data, restore_flds, restore_mask, restore_width, and restore_timescale which work when restore_ice is on to computed interior restoring terms.  Calls to initialize and execute the boundary and interior restoring were updated in the code.

In summary, several model defined fields can be restored (controlled by set_boundary_flds and restore_flds) in ice_restoring.F90.  The restoring dataset is specified with restore_data.  Exterior halo values are set with set_boundary_flds.  Interior restoring is turned on with restore_ice and the restoring strength is specified with restore_mask, restore_width, and restore_timescale.  The implementation provides some basic options and a framework for adding new fields, new masks, and new datasets by application users as needed.

The current restore_data options are
   - defined, internally defined state at all time
   - initial, uses the initial model state at all time
   - restartfiles, reads in extended restart files at each timestep, filenames are somewhat hardcoded into the implementation

The fields (set_boundary_flds, restore_flds) currently supported are
   - aicen, restore ice concentration in all categories
   - state, restore aicen, vicen, vsnon, and trcrn
   - trcrn, restore all tracers in all categories
   - vicen, restore ice thickness in all categories
   - vsnon,  restore snow thickness in all categories
   - velocity, restore uvel and vvel

The interior restoring strength is set as follows.  The interior mask (restore_mask) specifies a mask that is applied to interior points based on restore_width (the number of gridcells around the exterior to restore) and restore_timescale (the restoring time in days) as follows
   - all, sets all interior points with full restoring value (restore_timescale)
   - none, no interior restoring is carried out
   - constant, restores restore_width gridpoints on the exterior of the domain with full restoring value (restore_timescale)
   - linear, restores restore_width gridpoints on the exterior of the domain with a linearly decreasing restoring value moving from the exterior (set to restore_timescale) to the interior.

The boundary and interior restoring uses the same internal arrays to hold the restoring data.  The boundary updates always just set the values on the exterior halo to the restoring values.  The interior restoring updates the interior fields with a forcing strength proportional to the model timestep (dt) over the restoring time (restore_timescale).  It's important not to do interior restoring more than once per model timestep for each field.  The halo updates are done via subroutine ice_restoring_halo while the interior restoring is done via calls to subroutine ice_restoring_interior.

A new timer was added, timer_restore, to track restoring time.  Documentation was updated, and several tests were added or updated.  There is a new script, configuration/scripts/tests/full2nest.sh, that provides a basis for extracting a smaller (nest) regional grid from a (full) run and generating the extended restart files to force the boundary and interior for the nest testcase.  The full2nest.sh script describes how to setup that test.

All results are bit-for-bit except the boxrestore test case which was modified to produce new results.

There are still several issues to investigate
  - Can the exterior halo values for grid and technical fields (like the advection metrics) be computed directly and exactly for a regional grid rather than by a "zero_gradient" type implementation?  Up to now, the only time these values were needed on the halo was in cyclic cases and those halo values can be set by the cyclic bc exactly.
  - Can/should we be able to mix "zero_gradient" and "linear_extrapolation" for regional bcs in a model run depending on fields (maybe velocity should be extrapolated and concentration should not?)
  - We should implement restoring for aice, vice, vsno, tracers (not defined by category but then distributed across categories) and maybe subsets of tracer fields (instead of all or none).
  - Will the regional grids work better with 2 or 3 ghost cells?  We need to validate the model with ghost cells greater than 1 and test whether extending the halo and specifying the halo values will improve the regional solutions.
  - Are the restoring calls in the best place in the code?  For instance, can we implement interior restoring in dynamics subcycling (ndtd or ndte) while preserving the restoring forcing strength somehow?

Several other minor modifications were made
  - bound_state was separated from ice_state into it's own module due to some circular dependencies created by adding ice_restoring_halo calls to bound_state to update the exterior halo when required.  Also added an optional argument that turns off any restoring in bound_state, needed during initialization.
  - deprecated restore_bgc feature
  - maxval_spval in runtime_diagnostics was updated to better reflect the implementation
  - the restoring calculation for restore_ocn and restore_bgc was maintained but updated to reduce flops
  - added some memory allocation checks
  - updated some abort_ice calls
  - cleaned up indentation in ice_forcing_bgc.F90
  - updated restoring calls in all drivers including unit test drivers as needed
  - removed a couple of unused variables
  - updated the Boreas port

Author: apcraig

cicecore/cicedyn/analysis/ice_diagnostics.F90

@@ -207,11 +207,11 @@ subroutine runtime_diags (dt)
          work1, work2      ! temporary
 
       real (kind=dbl_kind), parameter :: &
-         maxval_spval = -0.9_dbl_kind*HUGE(0.0_dbl_kind)  ! spval to detect
+         maxval_spval = -0.9_dbl_kind*1.0e36_dbl_kind  ! spval to detect
                  ! undefined values returned from global_maxval.  if global_maxval
                  ! is applied to a region that does not exist (for instance
                  ! southern hemisphere in box cases), global_maxval
-                 ! returns -HUGE which we want to avoid writing.  The
+                 ! returns -1.0e36 which we want to avoid writing.  The
                  ! return value is checked against maxval_spval before writing.
 
       character(len=*), parameter :: subname = '(runtime_diags)'

cicecore/cicedyn/dynamics/ice_dyn_eap.F90

@@ -121,6 +121,7 @@ subroutine eap (dt)
       use ice_grid, only: tmask, umask, dxT, dyT, dxU, dyU, &
           tarear, uarear, grid_average_X2Y, &
           grid_atm_dynu, grid_atm_dynv, grid_ocn_dynu, grid_ocn_dynv
+      use ice_restoring, only: ice_restoring_halo
       use ice_state, only: aice, aiU, vice, vsno, uvel, vvel, divu, shear, vort, &
           aice_init, aice0, aicen, vicen, strength
       use ice_timers, only: timer_dynamics, timer_bound, &
@@ -520,6 +521,8 @@ subroutine eap (dt)
          call ice_timer_stop(timer_bound)
          call unstack_fields(fld2, uvel, vvel)
 
+         call ice_restoring_halo(setfld='velocity')
+
       enddo                     ! subcycling
 
       if (maskhalo_dyn) call ice_HaloDestroy(halo_info_mask)

cicecore/cicedyn/dynamics/ice_dyn_evp.F90

@@ -285,6 +285,7 @@ subroutine evp (dt)
           dxE, dxN, dxT, dxU, dyE, dyN, dyT, dyU, &
           tarear, uarear, earear, narear, grid_average_X2Y, uarea, &
           grid_ice, grid_atm_dynu, grid_atm_dynv, grid_ocn_dynu, grid_ocn_dynv
+      use ice_restoring, only: ice_restoring_halo
       use ice_state, only: aice, aiU, vice, vsno, uvel, vvel, uvelN, vvelN, &
           uvelE, vvelE, divu, shear, vort, &
           aice_init, aice0, aicen, vicen, strength
@@ -907,6 +908,7 @@ subroutine evp (dt)
                call dyn_haloUpdate (halo_info,          halo_info_mask,    &
                                     field_loc_NEcorner, field_type_vector, &
                                     uvel, vvel)
+               call ice_restoring_halo(setfld='velocity')
 
             enddo  ! sub cycling
          endif ! evp algorithm
@@ -1094,6 +1096,7 @@ subroutine evp (dt)
             call dyn_haloUpdate (halo_info,          halo_info_mask,    &
                                  field_loc_NEcorner, field_type_vector, &
                                  uvel, vvel)
+            call ice_restoring_halo(setfld='velocity')
 
          enddo                     ! subcycling
 
@@ -1285,6 +1288,7 @@ subroutine evp (dt)
             call dyn_haloUpdate (halo_info,          halo_info_mask,    &
                                  field_loc_NEcorner, field_type_vector, &
                                  uvel, vvel)
+            call ice_restoring_halo(setfld='velocity')
 
          enddo                     ! subcycling
 

cicecore/cicedyn/dynamics/ice_dyn_vp.F90

@@ -682,6 +682,7 @@ subroutine anderson_solver (icellT  , icellU , &
       use ice_grid, only: dxT, dyT, uarear
       use ice_dyn_shared, only: DminTarea, dxhy, dyhx, cxp, cyp, cxm, cym
       use ice_state, only: uvel, vvel, strength
+      use ice_restoring, only: ice_restoring_halo
       use ice_timers, only: ice_timer_start, ice_timer_stop, timer_bound
 
       integer (kind=int_kind), intent(in) :: &
@@ -1092,6 +1093,8 @@ subroutine anderson_solver (icellT  , icellU , &
          call ice_timer_stop(timer_bound)
          call unstack_fields(fld2, uvel, vvel)
 
+         call ice_restoring_halo(setfld='velocity')
+
          ! Compute "progress" residual norm
          !$OMP PARALLEL DO PRIVATE(iblk)
          do iblk = 1, nblocks

cicecore/cicedyn/dynamics/ice_transport_driver.F90

@@ -258,7 +258,8 @@ subroutine transport_remap (dt)
       use ice_domain_size, only: ncat, max_blocks
       use ice_blocks, only: nx_block, ny_block, block, get_block, nghost
       use ice_state, only: aice0, aicen, vicen, vsnon, trcrn, &
-          uvel, vvel, bound_state, uvelE, vvelN
+          uvel, vvel, uvelE, vvelN
+      use ice_bound_state, only: bound_state
       use ice_grid, only: tarea, grid_ice
       use ice_calendar, only: istep1
       use ice_timers, only: ice_timer_start, ice_timer_stop, &
@@ -706,13 +707,15 @@ end subroutine transport_remap
 
       subroutine transport_upwind (dt)
 
-      use ice_boundary, only: ice_HaloUpdate
+      use ice_boundary, only: ice_HaloUpdate, ice_HaloExtrapolate
       use ice_blocks, only: nx_block, ny_block, block, get_block, nx_block, ny_block
-      use ice_domain, only: blocks_ice, halo_info, nblocks
+      use ice_domain, only: blocks_ice, halo_info, nblocks, distrb_info, &
+          ew_boundary_type, ns_boundary_type
       use ice_domain_size, only: ncat, max_blocks
       use ice_state, only: aice0, aicen, vicen, vsnon, trcrn, &
-          uvel, vvel, trcr_depend, bound_state, trcr_base, &
+          uvel, vvel, trcr_depend, trcr_base, &
           n_trcr_strata, nt_strata, uvelE, vvelN
+      use ice_bound_state, only: bound_state
       use ice_flux, only: Tf
       use ice_grid, only: HTE, HTN, tarea, tmask, grid_ice
       use ice_timers, only: ice_timer_start, ice_timer_stop, &
@@ -793,6 +796,10 @@ subroutine transport_upwind (dt)
          !$OMP END PARALLEL DO
 
          call ice_timer_start(timer_bound)
+         call ice_HaloExtrapolate(uee, distrb_info, &
+                                  ew_boundary_type, ns_boundary_type)
+         call ice_HaloExtrapolate(vnn, distrb_info, &
+                                  ew_boundary_type, ns_boundary_type)
          call ice_HaloUpdate (uee,             halo_info,     &
                               field_loc_Eface, field_type_vector)
          call ice_HaloUpdate (vnn,             halo_info,     &

cicecore/cicedyn/general/ice_bound_state.F90

@@ -0,0 +1,152 @@
+!=======================================================================
+!
+! Primary state variables in various configurations
+! Note: other state variables are at the end of this...
+! The primary state variable names are:
+!-------------------------------------------------------------------
+! for each category   aggregated over     units
+!                       categories
+!-------------------------------------------------------------------
+! aicen(i,j,n)         aice(i,j)           ---
+! vicen(i,j,n)         vice(i,j)           m
+! vsnon(i,j,n)         vsno(i,j)           m
+! trcrn(i,j,it,n)      trcr(i,j,it)
+!
+! Area is dimensionless because aice is the fractional area
+! (normalized so that the sum over all categories, including open
+! water, is 1.0).  That is why vice/vsno have units of m instead of m^3.
+!
+! Variable names follow these rules:
+!
+! (1) For 3D variables (indices i,j,n), write 'ice' or 'sno' or
+!     'sfc' and put an 'n' at the end.
+! (2) For 2D variables (indices i,j) aggregated over all categories,
+!     write 'ice' or 'sno' or 'sfc' without the 'n'.
+! (3) For 2D variables (indices i,j) associated with an individual
+!     category, write 'i' or 's' instead of 'ice' or 'sno' and put an 'n'
+!     at the end: e.g. hin, hsn.  These are not declared here
+!     but in individual modules (e.g., icepack_therm_vertical).
+!
+! authors C. M. Bitz, UW
+!         Elizabeth C. Hunke and William H. Lipscomb, LANL
+!
+! 2004: Block structure added by William Lipscomb
+! 2006: Converted to free form source (F90) by Elizabeth Hunke
+
+      module ice_bound_state
+
+      use ice_kinds_mod
+      use ice_constants, only: field_loc_center, field_type_scalar, c0
+      use ice_exit, only: abort_ice
+      use ice_fileunits, only: nu_diag
+!      use icepack_intfc, only: icepack_query_tracer_sizes
+!      use icepack_intfc, only: icepack_warnings_flush, icepack_warnings_aborted
+
+      implicit none
+      private
+      public :: bound_state
+
+!=======================================================================
+
+      contains
+
+!=======================================================================
+!
+! Get ghost cell values for ice state variables in each thickness category.
+! NOTE: This subroutine cannot be called from inside a block loop!
+!
+! author: William H. Lipscomb, LANL
+
+      subroutine bound_state (aicen,        &
+                              vicen, vsnon, &
+                              ntrcr, trcrn, &
+                              restore)
+
+      use ice_boundary, only: ice_halo, ice_HaloMask, ice_HaloUpdate, &
+          ice_HaloDestroy
+      use ice_domain_size, only: max_blocks, ncat
+      use ice_blocks, only: nx_block, ny_block
+      use ice_domain, only: halo_info, maskhalo_bound, nblocks
+      use ice_restoring, only: ice_restoring_halo
+
+      integer (kind=int_kind), intent(in) :: &
+         ntrcr     ! number of tracers in use
+
+      real (kind=dbl_kind), dimension(nx_block,ny_block,ncat,max_blocks), intent(inout) :: &
+         aicen , & ! fractional ice area
+         vicen , & ! volume per unit area of ice          (m)
+         vsnon     ! volume per unit area of snow         (m)
+
+      real (kind=dbl_kind), intent(inout), dimension(:,:,:,:,:) :: &  ! (nx_block,ny_block,ntrcr,ncat,max_blocks)
+         trcrn     ! ice tracers
+
+      logical (kind=log_kind), intent(in), optional :: &
+         restore
+
+      ! local variables
+
+      integer (kind=int_kind) :: i, j, n, iblk
+
+      integer (kind=int_kind), &
+         dimension(nx_block,ny_block,max_blocks) :: halomask
+
+      type (ice_halo) :: halo_info_aicemask  ! halo mask
+
+      logical (kind=log_kind) :: lrestore  ! local restart flag
+
+      character(len=*), parameter :: subname = '(bound_state)'
+
+      if (present(restore)) then
+         lrestore = restore
+      else
+         lrestore = .true.
+      endif
+
+      call ice_HaloUpdate (aicen,            halo_info, &
+                           field_loc_center, field_type_scalar)
+
+      if (maskhalo_bound) then
+         halomask(:,:,:) = 0
+
+         !$OMP PARALLEL DO PRIVATE(iblk,n,i,j)
+         do iblk = 1, nblocks
+         do n = 1, ncat
+         do j = 1, ny_block
+         do i = 1, nx_block
+            if (aicen(i,j,n,iblk) > c0) halomask(i,j,iblk) = 1
+         enddo
+         enddo
+         enddo
+         enddo
+         !$OMP END PARALLEL DO
+
+         call ice_HaloMask(halo_info_aicemask, halo_info, halomask)
+
+         call ice_HaloUpdate (trcrn(:,:,:,:,:), halo_info_aicemask, &
+                              field_loc_center, field_type_scalar)
+         call ice_HaloUpdate (vicen,            halo_info_aicemask, &
+                              field_loc_center, field_type_scalar)
+         call ice_HaloUpdate (vsnon,            halo_info_aicemask, &
+                              field_loc_center, field_type_scalar)
+         call ice_HaloDestroy(halo_info_aicemask)
+
+      else
+         call ice_HaloUpdate (trcrn(:,:,:,:,:), halo_info, &
+                              field_loc_center, field_type_scalar)
+         call ice_HaloUpdate (vicen,            halo_info, &
+                              field_loc_center, field_type_scalar)
+         call ice_HaloUpdate (vsnon,            halo_info, &
+                              field_loc_center, field_type_scalar)
+      endif
+
+      if (lrestore) then
+         call ice_restoring_halo(setfld='state')
+      endif
+
+      end subroutine bound_state
+
+!=======================================================================
+
+      end module ice_bound_state
+
+!=======================================================================

cicecore/cicedyn/general/ice_forcing.F90

@@ -19,7 +19,6 @@
       module ice_forcing
 
       use ice_kinds_mod
-      use ice_boundary, only: ice_HaloUpdate
       use ice_blocks, only: nx_block, ny_block
       use ice_domain, only: halo_info
       use ice_domain_size, only: ncat, max_blocks, nx_global, ny_global, nfreq
@@ -34,7 +33,7 @@ module ice_forcing
                                 ice_open_nc, ice_read_nc, ice_close_nc
       use ice_timers, only: ice_timer_start, ice_timer_stop, timer_readwrite, &
                             timer_bound, timer_forcing
-      use ice_arrays_column, only: oceanmixed_ice, restore_bgc
+      use ice_arrays_column, only: oceanmixed_ice
       use ice_constants, only: c0, c1, c2, c3, c4, c5, c8, c10, c12, c15, c20, &
                                c180, c360, c365, c1000, c3600
       use ice_constants, only: p001, p01, p1, p2, p25, p5, p6
@@ -160,11 +159,8 @@ module ice_forcing
       logical (kind=log_kind), public :: &
          restore_ocn                 ! restore sst if true
 
-      integer (kind=int_kind), public :: &
-         trestore                    ! restoring time scale (days)
-
       real (kind=dbl_kind), public :: &
-         trest                       ! restoring time scale (sec)
+         trestore                    ! restoring time scale (days)
 
       logical (kind=log_kind), public :: &
          debug_forcing               ! prints forcing debugging output if true
@@ -188,6 +184,9 @@ module ice_forcing
 
       ! PRIVATE:
 
+      real (dbl_kind) :: &
+         frestore        ! restoring value based on trestore (0.-1.)
+
       real (dbl_kind), parameter :: &
          mixed_layer_depth_default = c20  ! default mixed layer depth in m
 
@@ -411,11 +410,11 @@ subroutine init_forcing_ocn(dt)
 
       nbits = 64              ! double precision data
 
-      if (restore_ocn .or. restore_bgc) then
-         if (trestore == 0) then
-            trest = dt        ! use data instantaneously
+      if (restore_ocn) then
+         if (trestore == c0) then
+            frestore = c1      ! use data instantaneously
          else
-            trest = real(trestore,kind=dbl_kind) * secday ! seconds
+            frestore = max(dt/(trestore*secday),c1)
          endif
       endif
 
@@ -731,17 +730,6 @@ subroutine get_forcing_atmo
       enddo                     ! iblk
       !$OMP END PARALLEL DO
 
-      call ice_timer_start(timer_bound)
-      call ice_HaloUpdate (swvdr,             halo_info, &
-                           field_loc_center,  field_type_scalar, fillvalue=c0)
-      call ice_HaloUpdate (swvdf,             halo_info, &
-                           field_loc_center,  field_type_scalar, fillvalue=c0)
-      call ice_HaloUpdate (swidr,             halo_info, &
-                           field_loc_center,  field_type_scalar, fillvalue=c0)
-      call ice_HaloUpdate (swidf,             halo_info, &
-                           field_loc_center,  field_type_scalar, fillvalue=c0)
-      call ice_timer_stop(timer_bound)
-
       call ice_timer_stop(timer_forcing)
 
       end subroutine get_forcing_atmo
@@ -3629,7 +3617,7 @@ subroutine ocn_data_clim (dt)
             do j = 1, ny_block
             do i = 1, nx_block
                sst(i,j,iblk) = sst(i,j,iblk)  &
-                         + (sstdat(i,j,iblk)-sst(i,j,iblk))*dt/trest
+                         + (sstdat(i,j,iblk)-sst(i,j,iblk))*frestore
             enddo
             enddo
          enddo
@@ -4085,7 +4073,7 @@ subroutine ocn_data_ncar(dt)
       if (restore_ocn) then
         do j = 1, ny_block
          do i = 1, nx_block
-           sst(i,j,:) = sst(i,j,:) + (work1(i,j,:)-sst(i,j,:))*dt/trest
+           sst(i,j,:) = sst(i,j,:) + (work1(i,j,:)-sst(i,j,:))*frestore
          enddo
         enddo
 !     else sst is only updated in ice_ocean.F
@@ -4288,7 +4276,7 @@ subroutine ocn_data_hadgem(dt)
             do j = 1, ny_block
             do i = 1, nx_block
                sst(i,j,iblk) = sst(i,j,iblk)  &
-                         + (sstdat(i,j,iblk)-sst(i,j,iblk))*dt/trest
+                         + (sstdat(i,j,iblk)-sst(i,j,iblk))*frestore
             enddo
             enddo
          enddo