Implement Runge-Kutta time integration

components/mpas-albany-landice/src/Registry.xml

@@ -292,6 +292,10 @@
                         description="If true, then distribute unablatedVolumeDynCell among dynamic neighbors when converting ablation velocity to ablation thickness. This should only be used as a clean-up measure, while limiting the timestep based on ablation velocity should be used as the primary method of getting accurate ablation thickness from ablation velocity.  If you choose to set config_adaptive_timestep_calvingCFL_fraction much larger than 1.0 (which is not recommended), setting this option to true usually results in more accurate calving behavior. "
                         possible_values=".true. or .false."
                 />
+                <nml_option name="config_update_velocity_before_calving" type="logical" default_value=".false." units="unitless"
+                            description="If true, add an additional velocity solve between advection and calving. If false, use velocity field from beginning of time step to calculate calving rate. For certain calving laws, like damage threshold calving, it is not necessary to update the velocity before calving, while for von Mises stress and eigencalving, it is more accurate to have an updated velocity state before solving for calvingThickness."
+                            possible_values=".true. or .false."
+                />
 	</nml_record>
 
 	<nml_record name="thermal_solver" in_defaults="true">
@@ -485,8 +489,16 @@
 		            possible_values="Any time interval of the format 'YYYY-MM-DD_HH:MM:SS', but limited by CFL condition."
 		/>
 		<nml_option name="config_time_integration" type="character" default_value="forward_euler" units="unitless"
-		            description="Time integration method (currently, only forward Euler is supported)."
-		            possible_values="'forward_euler'"
+		            description="Time integration method."
+		            possible_values="'forward_euler' or 'runge_kutta'"
+		/>
+                <nml_option name="config_rk_order" type="integer" default_value="2" units="unitless"
+                            description="Order of Runge-Kutta time integration to use. A value of 1 would be equivalent to forward euler, but will cause an error to avoid unnecessary redundancy. Values of 2 and 3 indicate strong-stability preserving RK2 and RK3. There is currently no support for classical RK2 or RK4 methods."
+                            possible_values="2 or 3"
+                />
+		<nml_option name="config_rk3_stages" type="integer" default_value="3" units="stages"
+                        description="Number of stages for strong stability preserving RK3 time integration. If set to 3, this involves 3 velocity solves and a maximum CFL fraction of 1. If set to 4, this involves 4 velocity solves, but the maximum CFL fraction is 2."
+                        possible_values="3 or 4"
 		/>
 		<nml_option name="config_adaptive_timestep" type="logical" default_value=".false." units="unitless"
 			description="Determines if the time step should be adjusted based on the CFL condition or should be steady in time. If true, the config_dt_* options are ignored."

components/mpas-albany-landice/src/landice.cmake

@@ -57,7 +57,7 @@ list(APPEND RAW_SOURCES
   core_landice/mode_forward/mpas_li_core.F
   core_landice/mode_forward/mpas_li_core_interface.F
   core_landice/mode_forward/mpas_li_time_integration.F
-  core_landice/mode_forward/mpas_li_time_integration_fe.F
+  core_landice/mode_forward/mpas_li_time_integration_fe_rk.F
   core_landice/mode_forward/mpas_li_diagnostic_vars.F
   core_landice/mode_forward/mpas_li_advection.F
   core_landice/mode_forward/mpas_li_calving.F

components/mpas-albany-landice/src/mode_forward/Makefile

@@ -4,7 +4,7 @@
 OBJS =  mpas_li_core.o \
     mpas_li_core_interface.o \
     mpas_li_time_integration.o \
-    mpas_li_time_integration_fe.o \
+    mpas_li_time_integration_fe_rk.o \
     mpas_li_diagnostic_vars.o \
     mpas_li_advection.o \
     mpas_li_advection_fct_shared.o \
@@ -35,9 +35,9 @@ mpas_li_core.o: mpas_li_time_integration.o \
                      mpas_li_statistics.o \
                      mpas_li_calving.o
 
-mpas_li_time_integration.o: mpas_li_time_integration_fe.o
+mpas_li_time_integration.o: mpas_li_time_integration_fe_rk.o
 
-mpas_li_time_integration_fe.o: mpas_li_advection.o \
+mpas_li_time_integration_fe_rk.o: mpas_li_advection.o \
 	                       mpas_li_calving.o \
                                mpas_li_thermal.o \
 			       mpas_li_iceshelf_melt.o \

components/mpas-albany-landice/src/mode_forward/mpas_li_advection.F

@@ -46,7 +46,8 @@ module li_advection
    !--------------------------------------------------------------------
    public :: li_advection_thickness_tracers, &
              li_layer_normal_velocity, &
-             li_update_geometry
+             li_update_geometry, &
+             li_grounded_to_floating
 
    !--------------------------------------------------------------------
    !
@@ -629,10 +630,6 @@ subroutine li_advection_thickness_tracers(&
          call li_calculate_mask(meshPool, velocityPool, geometryPool, err_tmp)
          err = ior(err, err_tmp)
 
-         ! Calculate volume converted from grounded to floating
-         ! This needs to be determined after SMB/BMB are applied because those can change floating/grounded state
-         call grounded_to_floating(cellMaskTemporaryField % array, cellMask, thickness, groundedToFloatingThickness, nCells)
-
          ! Calculate flux across grounding line
          ! Do this after new thickness & mask have been calculated, including SMB/BMB
          fluxAcrossGroundingLine(:) = 0.0_RKIND
@@ -2107,7 +2104,7 @@ subroutine vertical_remap(thickness, cellMask, meshPool, layerThickness, tracers
 
     end subroutine vertical_remap
 
-    subroutine grounded_to_floating(cellMaskOrig, cellMaskNew, thicknessNew, groundedToFloatingThickness, nCells)
+    subroutine li_grounded_to_floating(cellMaskOrig, cellMaskNew, thicknessNew, groundedToFloatingThickness, nCells)
       !-----------------------------------------------------------------
       ! input variables
       !-----------------------------------------------------------------
@@ -2151,7 +2148,7 @@ subroutine grounded_to_floating(cellMaskOrig, cellMaskNew, thicknessNew, grounde
         endif
       enddo
 
-    end subroutine grounded_to_floating
+    end subroutine li_grounded_to_floating
 
 !***********************************************************************
 

components/mpas-albany-landice/src/mode_forward/mpas_li_calving.F

@@ -74,7 +74,7 @@ module li_calving
 !> (3) Calve ice based on an ice thickness threshold
 !-----------------------------------------------------------------------
 
-   subroutine li_calve_ice(domain, err)
+   subroutine li_calve_ice(domain, err, solveVeloAfterCalving)
 
       use li_advection
 
@@ -92,7 +92,7 @@ subroutine li_calve_ice(domain, err)
       ! output variables
       !-----------------------------------------------------------------
       integer, intent(out) :: err !< Output: error flag
-
+      logical, intent(out), optional :: solveVeloAfterCalving
       !-----------------------------------------------------------------
       ! local variables
       !-----------------------------------------------------------------
@@ -111,13 +111,14 @@ subroutine li_calve_ice(domain, err)
       logical, pointer :: config_apply_calving_mask
       real(kind=RKIND), pointer :: config_calving_timescale
 
-      integer, pointer :: nCells
+      integer, pointer :: nCells, nCellsSolve
       integer, pointer :: config_number_of_blocks
 
       real (kind=RKIND), pointer :: deltat  !< time step (s)
 
       integer, dimension(:), pointer :: &
-           indexToCellID       ! list of global cell IDs
+           indexToCellID, &       ! list of global cell IDs
+           cellMask
 
       real (kind=RKIND) ::  &
            calvingFraction ! fraction of ice that calves in each column; depends on calving_timescale
@@ -136,14 +137,22 @@ subroutine li_calve_ice(domain, err)
 
       real (kind=RKIND), dimension(:), pointer :: calvingVelocity
 
+      real (kind=RKIND), dimension(:), pointer :: areaCell
+
       type (field1dReal), pointer :: originalThicknessField
 
       real (kind=RKIND), dimension(:), pointer :: originalThickness
 
+      type (field1dInteger), pointer :: cellMaskTemporaryField
+
       integer :: iCell
 
+      real (kind=RKIND), parameter :: smallNumber = 1.0e-6_RKIND
+
       integer :: err_tmp
 
+      real (kind=RKIND), dimension(1) :: calvingSumLocal, calvingSumGlobal
+
       err = 0
       err_tmp = 0
 
@@ -155,6 +164,15 @@ subroutine li_calve_ice(domain, err)
       call mpas_pool_get_config(liConfigs, 'config_data_calving', config_data_calving)
       call mpas_pool_get_config(liConfigs, 'config_number_of_blocks', config_number_of_blocks)
 
+      if (present(solveVeloAfterCalving)) then
+         call mpas_pool_get_subpool(domain % blocklist % structs, 'geometry', geometryPool)
+         call mpas_pool_get_field(geometryPool, 'cellMaskTemporary', cellMaskTemporaryField)
+         call mpas_allocate_scratch_field(cellMaskTemporaryField, .true.)
+         call mpas_pool_get_array(geometryPool, 'cellMask', cellMask)
+         ! Store cellMask prior to calving
+         cellMaskTemporaryField % array(:) = cellMask(:)
+      endif
+
       ! Zero calvingThickness here instead of or in addition to in individual subroutines.
       ! This is necessary when using damage threshold calving with other calving
       ! routines. Some individual routines still set calvingThickness to zero, but
@@ -206,7 +224,7 @@ subroutine li_calve_ice(domain, err)
       ! However, the eigencalving method requires multiple applications of the calvingThickness
       ! to the thickness.  So the simplest method to apply data calving is to store the old
       ! thickness and then set it back when we are done.
-      if (config_data_calving) then
+      if ( config_data_calving .or. present(solveVeloAfterCalving) ) then
          call mpas_pool_get_subpool(domain % blocklist % structs, 'scratch', scratchPool)
          call mpas_pool_get_field(scratchPool, 'workCell2',  originalThicknessField)
          call mpas_allocate_scratch_field(originalThicknessField, single_block_in = .false.)
@@ -337,8 +355,22 @@ subroutine li_calve_ice(domain, err)
 
          block => block % next
       end do
+      if (present(solveVeloAfterCalving)) then
+         call mpas_pool_get_array(meshPool, 'areaCell', areaCell)
+         call mpas_pool_get_dimension(meshPool, 'nCellsSolve', nCellsSolve)
+         call mpas_pool_get_array(geometryPool, 'calvingThickness', calvingThickness)
+         calvingSumLocal(1) = sum(calvingThickness(1:nCellsSolve) * areaCell(1:nCellsSolve) * &
+                               real(li_mask_is_dynamic_ice_int(cellMaskTemporaryField % array(1:nCellsSolve)), RKIND))
+         call mpas_dmpar_sum_real_array(domain % dminfo, 1, calvingSumLocal(1), calvingSumGlobal(1))
+         if (calvingSumGlobal(1) > smallNumber) then
+            solveVeloAfterCalving = .true.
+         else
+            solveVeloAfterCalving = .false.
+         endif
+         call mpas_deallocate_scratch_field(cellMaskTemporaryField, .true.)
+      endif
 
-      if (config_data_calving) then
+      if ( config_data_calving .or. present(solveVeloAfterCalving) ) then
          call mpas_deallocate_scratch_field(originalThicknessField, single_block_in=.false.)
       endif
 

components/mpas-albany-landice/src/mode_forward/mpas_li_time_integration.F

@@ -26,7 +26,7 @@ module li_time_integration
    use mpas_timekeeping
    use mpas_log
 
-   use li_time_integration_fe
+   use li_time_integration_fe_rk
    use li_setup
    use li_constants
 
@@ -166,10 +166,9 @@ subroutine li_timestep(domain, err)
       !call mpas_log_write('Using ' // trim(config_time_integration) // ' time integration.')
       select case (config_time_integration)
       case ('forward_euler')
-         call li_time_integrator_forwardeuler(domain, err_tmp)
-      case ('rk4')
-         call mpas_log_write(trim(config_time_integration) // ' is not currently supported.', MPAS_LOG_ERR)
-         err_tmp = 1
+         call li_time_integrator_forwardeuler_rungekutta(domain, err_tmp)
+      case ('runge_kutta')
+         call li_time_integrator_forwardeuler_rungekutta(domain, err_tmp)
       case default
          call mpas_log_write(trim(config_time_integration) // ' is not a valid land ice time integration option.', MPAS_LOG_ERR)
          err_tmp = 1