Merge branch 'sbrus89/ocn/rk4-time-fix' (PR #6695)

Evaluate time-dependent tendencies at correct time for each RK4 stage

This PR corrects a longstanding issue with RK4 and tendency terms which
have an explicit time dependence. Big thanks to @gcapodag and
@jeremy-lilly for providing the fix in their LTS development branch.

Here, I have migrated these changes from their MPAS-Dev branch to E3SM
and added the modifications to correct the convergence issues for the
manufactured solution test case

[BFB] - mpaso standalone

components/mpas-framework/src/framework/mpas_forcing.F

@@ -35,7 +35,8 @@ module mpas_forcing
        mpas_forcing_init_field_data, &
        mpas_forcing_get_forcing, &
        mpas_forcing_get_forcing_time, &
-       mpas_forcing_write_restart_times
+       mpas_forcing_write_restart_times, &
+       mpas_advance_forcing_clock
 
 contains
 
@@ -1193,7 +1194,7 @@ subroutine mpas_forcing_get_forcing(&!{{{
        if (trim(forcingGroup % forcingGroupName) == trim(forcingGroupName)) then
 
           ! advance the forcing time
-          call advance_forcing_clock(forcingGroup, dt)
+          call mpas_advance_forcing_clock(forcingGroup, dt)
 
           ! cycle the forcing clock
           if (forcingGroup % forcingCycleUse) then
@@ -1230,7 +1231,7 @@ end subroutine mpas_forcing_get_forcing!}}}
 
 !|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
 !
-!  advance_forcing_clock
+!  mpas_advance_forcing_clock
 !
 !> \brief set the forcing clock
 !> \author Adrian K. Turner, LANL
@@ -1241,7 +1242,7 @@ end subroutine mpas_forcing_get_forcing!}}}
 !
 !-----------------------------------------------------------------------
 
-  subroutine advance_forcing_clock(&!{{{
+  subroutine mpas_advance_forcing_clock(&!{{{
        forcingGroup, &
        dt)
 
@@ -1256,9 +1257,10 @@ subroutine advance_forcing_clock(&!{{{
 
     ! increment clock with timestep
     call mpas_set_timeInterval(timeStep, dt=dt)
+    !call mpas_log_write('time step from mpas_set_timeInterval: $i $i $i', intArgs=(/int(timeStep%ti%basetime%S), int(timeStep%ti%basetime%Sn), int(timeStep%ti%basetime%Sd) /))
     call mpas_advance_clock(forcingGroup % forcingClock, timeStep)
 
-  end subroutine advance_forcing_clock!}}}
+  end subroutine mpas_advance_forcing_clock!}}}
 
 !|||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||||
 !
@@ -1597,8 +1599,10 @@ subroutine get_interpolants_linear(interpolants, forcingStream, currentTime)!{{{
     call mpas_get_timeInterval(diff1, forcingStream % forcingTimes(1), dt=diffr1)
     call mpas_get_timeInterval(diff2, currentTime, dt=diffr2)
 
+    !call mpas_log_write('diffr2 $r, diffr, $r', realArgs=(/ diffr2, diffr /))
     interpolants(1) = diffr2 / diffr
     interpolants(2) = 1.0_RKIND - interpolants(1) !diffr1 / diffr
+
 
   end subroutine get_interpolants_linear!}}}
 
@@ -1612,7 +1616,7 @@ end subroutine get_interpolants_linear!}}}
 !> \details
 !>  Given the current time and forcing times calculate the correct
 !>  interpolants with piecewise constant interpolation
-!
+
 !-----------------------------------------------------------------------
 
   subroutine get_interpolants_constant(interpolants, forcingStream, currentTime)!{{{

components/mpas-framework/src/framework/mpas_timekeeping.F

@@ -1427,6 +1427,8 @@ subroutine mpas_set_timeInterval(interval, YY, MM, DD, H, M, S, S_n, S_d, S_i8,
             timeString_ = timeString
          end if
 
+         !call mpas_log_write('timeString_ '//trim(timeString_))
+
          numerator = 0
          denominator = 1
 
@@ -1519,6 +1521,12 @@ subroutine mpas_set_timeInterval(interval, YY, MM, DD, H, M, S, S_n, S_d, S_i8,
             return
          end if
 
+         if (sec < 0) then
+            numerator = -numerator
+         end if
+
+         !call mpas_log_write('sec $i, num $i, denom $i', intArgs=(/int(sec),numerator,denominator/))
+
          call ESMF_TimeIntervalSet(interval % ti, YY=year, MM=month, D=day, H=hour, M=min, S_i8=sec, Sn=numerator, Sd=denominator, rc=ierr)
 
       else

components/mpas-ocean/src/Registry.xml

@@ -4175,6 +4175,10 @@
 			 description="Multiplication factors to smooth ssh at coastlines for SAL caculation"
 			 packages="tidalPotentialForcingPKG"
 		/>
+               	<var name="forcingTimeIncrement" type="real" dimensions="" units="s"
+			description="Time increment to compute forcing terms at intermediate time intervals"
+			packages="forwardMode"
+		/>
 	</var_struct>
 	<var_struct name="timeVaryingForcing" time_levs="1">
 		<var name="windSpeedU" type="real" dimensions="nCells Time" units="m s^-1"

components/mpas-ocean/src/mode_forward/mpas_ocn_forward_mode.F

@@ -26,6 +26,7 @@ module ocn_forward_mode
    use mpas_stream_manager
    use mpas_timekeeping
    use mpas_dmpar
+   use mpas_forcing
    use mpas_timer
    use mpas_log
    use mpas_decomp
@@ -94,6 +95,7 @@ module ocn_forward_mode
 
    use ocn_forcing
    use ocn_time_varying_forcing
+   use ocn_framework_forcing
 
    use ocn_constants
    use ocn_config
@@ -664,7 +666,12 @@ function ocn_forward_mode_run(domain) result(ierr)!{{{
 
       ! initialize time-varying forcing
       call ocn_time_varying_forcing_init(domain)
-      call ocn_time_varying_forcing_get(domain % streamManager, domain, domain % clock)
+
+      ! if not using RK4, calculate time varying forcing terms once per
+      ! time-step as opposed at each RK substage as implemented in RK4 
+      if (timeIntegratorChoice /= timeIntRK4) then
+        call ocn_time_varying_forcing_get(domain % streamManager, domain, domain % clock)
+      endif
 
       ! During integration, time level 1 stores the model state at the beginning of the
       !   time step, and time level 2 stores the state advanced dt in time by timestep(...)
@@ -834,7 +841,17 @@ function ocn_forward_mode_run(domain) result(ierr)!{{{
          endif
 
          ! read in next time level data required for time-varying forcing
-         call ocn_time_varying_forcing_get(domain % streamManager, domain, domain % clock)
+         if (timeIntegratorChoice /= timeIntRK4) then
+           ! if not using RK4, calculate time varying forcing terms once per
+           ! time-step as opposed at each RK substage as implemented in RK4
+           call ocn_time_varying_forcing_get(domain % streamManager, domain, domain % clock)
+         else
+           if (config_use_time_varying_atmospheric_forcing .or. &
+               config_use_time_varying_land_ice_forcing) then
+             ! increment forcing clock to next time-step
+             call mpas_advance_forcing_clock(forcingGroupHead, dt)
+           endif
+         endif
 
          ! Validate that the state is OK to run with for the next timestep.
          call ocn_validate_state(domain, timeLevel=1)

components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration.F

@@ -48,6 +48,19 @@ module ocn_time_integration
    !
    !--------------------------------------------------------------------
 
+    ! Enum for selecting different time integrators
+    integer, public :: timeIntegratorChoice
+
+    integer, public, parameter :: &
+       timeIntUnknown          = 0, &! unknown or undefined
+       timeIntSplitExplicit    = 1, &! split-explicit
+       timeIntUnsplitExplicit  = 2, &! unsplit-explicit
+       timeIntSemiImplicit     = 3, &! Semi-implicit
+       timeIntRK4              = 4, &! 4th-order Runge-Kutta
+       timeIntLTS              = 5, &! local time-stepping
+       timeIntFBLTS            = 6, &! forward-backward lts
+       timeIntSplitExplicitAB2 = 7   ! split-explicit AB2 baroclinic
+
    !--------------------------------------------------------------------
    !
    ! Public member functions
@@ -63,18 +76,6 @@ module ocn_time_integration
    !
    !--------------------------------------------------------------------
 
-    ! Enum for selecting different time integrators
-    integer :: timeIntegratorChoice
-
-    integer, parameter :: &
-       timeIntUnknown          = 0, &! unknown or undefined
-       timeIntSplitExplicit    = 1, &! split-explicit
-       timeIntUnsplitExplicit  = 2, &! unsplit-explicit
-       timeIntSemiImplicit     = 3, &! Semi-implicit
-       timeIntRK4              = 4, &! 4th-order Runge-Kutta
-       timeIntLTS              = 5, &! local time-stepping
-       timeIntFBLTS            = 6, &! forward-backward lts
-       timeIntSplitExplicitAB2 = 7   ! split-explicit AB2 baroclinic
 
 !***********************************************************************
 

components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_rk4.F

@@ -43,6 +43,7 @@ module ocn_time_integration_rk4
    use ocn_effective_density_in_land_ice
    use ocn_surface_land_ice_fluxes
    use ocn_transport_tests
+   use ocn_time_varying_forcing
 
    use ocn_subgrid
 
@@ -118,7 +119,7 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
 
       type (mpas_pool_type), pointer :: nextProvisPool, prevProvisPool
 
-      real (kind=RKIND), dimension(4) :: rk_weights, rk_substep_weights
+      real (kind=RKIND), dimension(4) :: rk_weights, rk_substep_weights, forcingTimeIncrementRK4
 
       real (kind=RKIND) :: coef
       real (kind=RKIND), dimension(:,:), pointer :: &
@@ -186,6 +187,7 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
       ! Tidal boundary condition
       logical, pointer :: config_use_tidal_forcing
       character (len=StrKIND), pointer :: config_tidal_forcing_type
+      real (kind=RKIND), pointer :: forcingTimeIncrement
       real (kind=RKIND), dimension(:), pointer :: tidalInputMask, tidalBCValue
       real (kind=RKIND), dimension(:,:), pointer :: restingThickness
       real (kind=RKIND) :: totalDepth
@@ -233,6 +235,7 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
          call mpas_pool_get_subpool(statePool, 'tracers', tracersPool)
          call mpas_pool_get_subpool(block % structs, 'mesh', meshPool)
          call mpas_pool_get_subpool(block % structs, 'diagnostics', diagnosticsPool)
+         call mpas_pool_get_subpool(block % structs, 'forcing', forcingPool)
 
          call mpas_pool_get_subpool(block % structs, 'provis_state', provisStatePool)
 
@@ -260,6 +263,9 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
          call mpas_pool_get_array(meshPool, 'subgridSshCellTableRange', &
                                 subgridSshCellTableRange)
 
+         call mpas_pool_get_array(forcingPool, 'forcingTimeIncrement', forcingTimeIncrement)
+
+         forcingTimeIncrement = 0.0_RKIND
 
          ! Lower k-loop limit of 1 rather than minLevel* needed in *New = *Cur
          ! assignments below are needed to maintain bit-for-bit results
@@ -343,10 +349,10 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
       ! The coefficients of k_j are b_j = (1/6, 1/3, 1/3, 1/6) and are
       ! initialized here as delta t * b_j:
 
-      rk_weights(1) = dt/6.
-      rk_weights(2) = dt/3.
-      rk_weights(3) = dt/3.
-      rk_weights(4) = dt/6.
+      rk_weights(1) = dt/6.0_RKIND
+      rk_weights(2) = dt/3.0_RKIND
+      rk_weights(3) = dt/3.0_RKIND
+      rk_weights(4) = dt/6.0_RKIND
 
       ! The a_j coefficients of h in the computation of k_j are typically written (0, 1/2, 1/2, 1).
       ! However, in the algorithm below we pre-compute the state for the tendency one iteration early.
@@ -355,11 +361,19 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
       ! That is why the coefficients of h are one index early in the following, i.e.
       ! a = (1/2, 1/2, 1)
 
-      rk_substep_weights(1) = dt/2.
-      rk_substep_weights(2) = dt/2.
+      rk_substep_weights(1) = dt/2.0_RKIND
+      rk_substep_weights(2) = dt/2.0_RKIND
       rk_substep_weights(3) = dt
       rk_substep_weights(4) = dt ! a_4 only used for ALE step, otherwise it is skipped.
 
+      ! these are time increments to evaluate the tidal forcing at the
+      ! intermediate time-steps as required by RK4
+
+      forcingTimeIncrementRK4(1) = 0.0_RKIND
+      forcingTimeIncrementRK4(2) = dt/2.0_RKIND
+      forcingTimeIncrementRK4(3) = dt/2.0_RKIND
+      forcingTimeIncrementRK4(4) = dt
+
       call mpas_timer_start("RK4-main loop")
 
       !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
@@ -481,8 +495,17 @@ subroutine ocn_time_integrator_rk4(domain, dt)!{{{
 
         block => domain % blocklist
         do while (associated(block))
-           call ocn_time_integrator_rk4_compute_vel_tends(domain, block, dt, rk_substep_weights(rk_step), domain % dminfo, err )
+           call mpas_pool_get_subpool(block % structs, 'forcing', forcingPool)
+           call mpas_pool_get_array(forcingPool, 'forcingTimeIncrement', forcingTimeIncrement)
+           forcingTimeIncrement = forcingTimeIncrementRK4(rk_step)
+           block => block % next
+        end do
+
+        call ocn_time_varying_forcing_get(domain % streamManager, domain, domain % clock)
 
+        block => domain % blocklist
+        do while (associated(block))
+           call ocn_time_integrator_rk4_compute_vel_tends(domain, block, dt, rk_substep_weights(rk_step), domain % dminfo, err )
            call ocn_time_integrator_rk4_compute_thick_tends( block, dt, rk_substep_weights(rk_step), err )
            block => block % next
         end do

components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_si.F

@@ -380,6 +380,8 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{
 
       real (kind=RKIND), dimension(:,:,:), pointer :: activeTracersNew
 
+      real (kind=RKIND), pointer :: forcingTimeIncrement
+
       ! Remap variables
       real (kind=RKIND), dimension(:,:), pointer :: &
          layerThicknessLagNew
@@ -507,6 +509,11 @@ subroutine ocn_time_integrator_si(domain, dt)!{{{
 
       call mpas_pool_get_array(tracersPool, 'activeTracers', activeTracersNew, 2)
 
+      call mpas_pool_get_array(forcingPool, 'forcingTimeIncrement', &
+                                             forcingTimeIncrement)
+
+      forcingTimeIncrement = 0.0_RKIND
+
       allocate(bottomDepthEdge(nEdgesAll+1))
 
       if (config_transport_tests_flow_id > 0) then

components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F

@@ -280,6 +280,8 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
 
       real (kind=RKIND), dimension(:,:,:), pointer :: activeTracersNew
 
+      real (kind=RKIND), pointer :: forcingTimeIncrement
+
       ! Remap variables
       real (kind=RKIND), dimension(:,:), pointer :: &
          layerThicknessLagNew
@@ -398,6 +400,11 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
 
       allocate(baroclinicThickness(nEdgesAll+1))
 
+      call mpas_pool_get_array(forcingPool, 'forcingTimeIncrement', &
+                                             forcingTimeIncrement)
+
+      forcingTimeIncrement = 0.0_RKIND
+
       if (config_transport_tests_flow_id > 0) then
         ! This is a transport test. Write advection velocity from prescribed
         ! flow field.

components/mpas-ocean/src/shared/mpas_ocn_manufactured_solution.F

@@ -76,7 +76,13 @@ module ocn_manufactured_solution
 !
 !-----------------------------------------------------------------------
 
-   subroutine ocn_manufactured_solution_tend_thick(tend, err)!{{{
+   subroutine ocn_manufactured_solution_tend_thick(forcingPool, tend, err)!{{{
+
+      !-----------------------------------------------------------------
+      ! input variables
+      !-----------------------------------------------------------------
+
+      type (mpas_pool_type), intent(in) :: forcingPool
 
       !-----------------------------------------------------------------
       ! input/output variables
@@ -97,14 +103,17 @@ subroutine ocn_manufactured_solution_tend_thick(tend, err)!{{{
 
       integer :: iCell, kmin, kmax, k
       real (kind=RKIND) :: phase, time
+      real (kind=RKIND), pointer :: forcingTimeIncrement
 
       ! End preamble
       !-------------
       ! Begin code
 
       if (.not. config_use_manufactured_solution) return
 
-      time = daysSinceStartOfSim*86400.0_RKIND
+      call mpas_pool_get_array(forcingPool, 'forcingTimeIncrement', forcingTimeIncrement) 
+
+      time = daysSinceStartOfSim*86400.0_RKIND + forcingTimeIncrement 
 
       do iCell = 1,nCellsOwned
 
@@ -137,7 +146,13 @@ end subroutine ocn_manufactured_solution_tend_thick!}}}
 !>  This routine computes the velocity tendency for the manufactured solution
 !
 !-----------------------------------------------------------------------
-   subroutine ocn_manufactured_solution_tend_vel(tend, err)!{{{
+   subroutine ocn_manufactured_solution_tend_vel(forcingPool, tend, err)!{{{
+
+      !-----------------------------------------------------------------
+      ! input variables
+      !-----------------------------------------------------------------
+
+      type (mpas_pool_type), intent(in) :: forcingPool
 
       !-----------------------------------------------------------------
       ! input/output variables
@@ -158,14 +173,17 @@ subroutine ocn_manufactured_solution_tend_vel(tend, err)!{{{
 
       integer :: iEdge, kmin, kmax, k
       real (kind=RKIND) :: phase, u, v, time
+      real (kind=RKIND), pointer :: forcingTimeIncrement
 
       ! End preamble
       !-----------------------------------------------------------------
       ! Begin code
 
       if (.not. config_use_manufactured_solution) return
 
-      time = daysSinceStartOfSim*86400.0_RKIND
+      call mpas_pool_get_array(forcingPool, 'forcingTimeIncrement', forcingTimeIncrement) 
+
+      time = daysSinceStartOfSim*86400.0_RKIND + forcingTimeIncrement 
 
       do iEdge = 1, nEdgesOwned