From a350b641ebf6ed7fa3b30ac1cc94a43779d7c8d2 Mon Sep 17 00:00:00 2001
From: Mark Petersen <mpetersen@lanl.gov>
Date: Mon, 18 Apr 2022 17:07:46 -0600
Subject: [PATCH 1/5] Add four split explicit changes. See PRs:

MPAS#728, MPAS ocean discussion #16
---
 components/mpas-ocean/src/Registry.xml        |  28 +++++
 .../mpas_ocn_time_integration_split.F         | 118 ++++++++++++++++--
 .../shared/mpas_ocn_diagnostics_variables.F   |  31 +++++
 3 files changed, 167 insertions(+), 10 deletions(-)

diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml
index 8eec749bb6ed..d6df61e02de3 100644
--- a/components/mpas-ocean/src/Registry.xml
+++ b/components/mpas-ocean/src/Registry.xml
@@ -1263,6 +1263,22 @@
 					description="number of iterations of the velocity corrector step in stage 2"
 					possible_values="any positive integer, but typically 1, 2, or 3"
 		/>
+		<nml_option name="config_btr_new1_resplit" type="logical" default_value=".false."
+					description="experimental flag"
+					possible_values=".true. or .false."
+		/>
+		<nml_option name="config_btr_new2_ubtr_foH" type="logical" default_value=".false."
+					description="experimental flag"
+					possible_values=".true. or .false."
+		/>
+		<nml_option name="config_btr_new3_sub_btrAvgbclV" type="logical" default_value=".false."
+					description="experimental flag"
+					possible_values=".true. or .false."
+		/>
+		<nml_option name="config_btr_new4_VelCor" type="logical" default_value=".false."
+					description="experimental flag"
+					possible_values=".true. or .false."
+		/>
 		<nml_option name="config_vel_correction" type="logical" default_value=".true."
 					description="If true, the velocity correction term is included in the horizontal advection of thickness and tracers"
 					possible_values=".true. or .false."
@@ -2722,6 +2738,18 @@
 			 description="horizontal velocity, tangential to an edge"
 			 packages="forwardMode;analysisMode"
 		/>
+		<var name="BtrAvgOfBaroclinicVelocity" type="real" dimensions="nEdges Time" units="m s^{-1}"
+			 description="barotropic average of baroclinic velocity"
+			 packages="splitTimeIntegrator"
+		/>
+		<var name="BtrAvgOfGMBolusVelocity" type="real" dimensions="nEdges Time" units="m s^{-1}"
+			 description="barotropic average of GM  Bolus velocity"
+			 packages="splitTimeIntegrator"
+		/>
+		<var name="normalVelocityCorrection" type="real" dimensions="nEdges Time" units="m s^{-1}"
+			 description="read it"
+			 packages="splitTimeIntegrator"
+		/>
 		<var name="layerThicknessEdgeMean" type="real" dimensions="nVertLevels nEdges Time" units="m"
 			 description="layer thickness averaged from cell center to edges"
 			 packages="forwardMode;analysisMode"
diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F
index 7f1ac1f921e8..c70174b94dfa 100644
--- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F
+++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F
@@ -179,7 +179,6 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
          flux,             &! temp for computing flux for barotropic
          sshEdge,          &! sea surface height at edge
          CoriolisTerm,     &! temp for computing coriolis term (fuperp)
-         normalVelocityCorrection, &! velocity correction
          temp,             &! temp for holding vars at new time
          temp_h,           &! temporary for phi
          temp_mask,        &! temporary mask (edgeMask)
@@ -407,6 +406,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
       ! Initialize * variables that are used to compute baroclinic
       ! tendencies below.
 
+if (.not.config_btr_new1_resplit) then
       ! The baroclinic velocity needs be recomputed at the beginning
       ! of a timestep because the implicit vertical mixing is
       ! conducted on the total u.  We keep normalBarotropicVelocity
@@ -445,11 +445,40 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
       end do
 #ifndef MPAS_OPENACC
       !$omp end do
+      !$omp end parallel
 #endif
 
+else ! config_btr_resplit
+
+      !$omp parallel
+      !$omp do schedule(runtime) private(k, thicknessSum)
+      do iEdge = 1,nEdgesAll
+
+         ! Compute barotropic component ubar at the beginning of every time step
+         thicknessSum  = layerThickEdgeMean(1,iEdge)
+         normalBarotropicVelocityCur(iEdge) = normalVelocityCur(1,iEdge) * layerThickEdgeMean(1,iEdge)
+         do k = 2, maxLevelEdgeTop(iEdge)
+            thicknessSum  =  thicknessSum + layerThickEdgeMean(k,iEdge)
+            normalBarotropicVelocityCur(iEdge) = normalBarotropicVelocityCur(iEdge) + normalVelocityCur(k,iEdge) * layerThickEdgeMean(k,iEdge)
+         enddo
+         normalBarotropicVelocityCur(iEdge) = normalBarotropicVelocityCur(iEdge) / thicknessSum
+
+         ! Compute baroclinic component u' = u - ubar at the beginning of every time step
+         do k = 1, nVertLevels !mrp, could tighten these bounds later: maxLevelEdgeTop % array(iEdge)
+            normalBaroclinicVelocityCur(k,iEdge) = normalVelocityCur(k,iEdge) - normalBarotropicVelocityCur(iEdge)
+            normalVelocityNew(k,iEdge) = normalVelocityCur(k,iEdge)
+            normalBaroclinicVelocityNew(k,iEdge) = normalBaroclinicVelocityCur(k,iEdge)
+         end do
+      end do
+      !$omp end do
+      !$omp end parallel
+
+end if ! config_btr_resplit
+
 #ifdef MPAS_OPENACC
       !$acc loop gang
 #else
+      !$omp parallel
       !$omp do schedule(runtime) private(k)
 #endif
       do iCell = 1, nCellsAll
@@ -1556,7 +1585,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
             ! efficiency: This next block of code is really a check for
             ! debugging, and can be removed later.
             call mpas_timer_start('btr se ssh verif')
-
+if (.not.config_btr_new4_VelCor) then
             ! Correction velocity
             !   normalVelocityCorrection = (Flux - Sum(h u*))/H
             ! or, for the full latex version:
@@ -1624,7 +1653,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
                                  layerThickEdgeFlux(k,iEdge)
                enddo
 
-               normalVelocityCorrection = useVelocityCorrection* &
+               normalVelocityCorrection(iEdge) = useVelocityCorrection* &
                              ((barotropicThicknessFlux(iEdge) - &
                                normalThicknessFluxSum)/thicknessSum)
 
@@ -1639,15 +1668,81 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
                   ! and tracers in tendency calls in stage 3.
                   normalTransportVelocity(k,iEdge) = edgeMask(k,iEdge)*( &
                              normalTransportVelocity(k,iEdge) + &
-                             normalVelocityCorrection)
+                             normalVelocityCorrection(iEdge))
+               enddo
+
+               k = minLevelEdgeBot(iEdge)
+               thicknessSum = 1.0e-15_RKIND
+               BtrAvgOfBaroclinicVelocity(iEdge) = 0.0_RKIND
+               BtrAvgOfGMBolusVelocity(iEdge) = 0.0_RKIND
+               ! note, BtrAvgOfGMBolusVelocity could be changed from an array
+               ! to a scalar if it is not needed to save for later inspection
+               do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge)
+                  thicknessSum  =  thicknessSum + layerThickEdgeFlux(k,iEdge)
+                  BtrAvgOfBaroclinicVelocity(iEdge) = BtrAvgOfBaroclinicVelocity(iEdge) + normalBaroclinicVelocityNew(k,iEdge) * layerThickEdgeFlux(k,iEdge)
+                  BtrAvgOfGMBolusVelocity(iEdge) = BtrAvgOfGMBolusVelocity(iEdge) + normalGMBolusVelocity(k,iEdge) * layerThickEdgeFlux(k,iEdge)
                enddo
+               BtrAvgOfBaroclinicVelocity(iEdge) = BtrAvgOfBaroclinicVelocity(iEdge) / thicknessSum
+               BtrAvgOfGMBolusVelocity(iEdge) = BtrAvgOfGMBolusVelocity(iEdge) / thicknessSum
+  if (config_btr_new2_ubtr_foH) then
+               normalBarotropicVelocityNew(iEdge) = barotropicThicknessFlux(iEdge) / thicknessSum
+  end if
 
             end do ! iEdge
             !$omp end do
             !$omp end parallel
 
             deallocate(uTemp)
+else ! if (config_btr_new4_VelCor) then
+            nEdges = nEdgesHalo(config_num_halos+1 ) !mrp check min w nEdgesPtr in PR
 
+            !$omp parallel
+            !$omp do schedule(runtime) private(k, thicknessSum)
+            do iEdge = 1, nEdges
+
+               ! thicknessSum is initialized outside the loop because on land boundaries
+               ! maxLevelEdgeTop=0, but I want to initialize thicknessSum with a
+               ! nonzero value to avoid a NaN.
+               k = minLevelEdgeBot(iEdge)
+               thicknessSum = 1.0e-15_RKIND
+               BtrAvgOfBaroclinicVelocity(iEdge) = 0.0_RKIND
+               BtrAvgOfGMBolusVelocity(iEdge) = 0.0_RKIND
+               ! note, BtrAvgOfGMBolusVelocity could be changed from an array
+               ! to a scalar if it is not needed to save for later inspection
+               do k = minLevelEdgeBot(iEdge), maxLevelEdgeTop(iEdge)
+                  thicknessSum  =  thicknessSum + layerThickEdgeFlux(k,iEdge)
+                  BtrAvgOfBaroclinicVelocity(iEdge) = BtrAvgOfBaroclinicVelocity(iEdge) + normalBaroclinicVelocityNew(k,iEdge) * layerThickEdgeFlux(k,iEdge)
+                  BtrAvgOfGMBolusVelocity(iEdge) = BtrAvgOfGMBolusVelocity(iEdge) + normalGMBolusVelocity(k,iEdge) * layerThickEdgeFlux(k,iEdge)
+               enddo
+               BtrAvgOfBaroclinicVelocity(iEdge) = BtrAvgOfBaroclinicVelocity(iEdge) / thicknessSum
+               BtrAvgOfGMBolusVelocity(iEdge) = BtrAvgOfGMBolusVelocity(iEdge) / thicknessSum
+
+               ! Compute barotropic velocity from thickness flux.
+               ! mrp changed back. With this commented out, barotropic velocity
+               ! is computed from average of subcycles above.
+               ! Before, we had
+               ! Set momentum velocity (normalVelocity) to:
+               ! u = F/H + u' - btrAvg( u' )
+               ! with this commented out, we have
+               ! u = subcycleAvg(ubar) + u' - btrAvg( u' )
+  if (config_btr_new2_ubtr_foH) then
+               normalBarotropicVelocityNew(iEdge) = barotropicThicknessFlux(iEdge) / thicknessSum
+  end if
+
+               do k = 1, nVertLevels
+
+                  ! Compute transport velocity, used to transport thickness and tracers.
+                  ! uTr = F/H + u' - btrAvg( u' ) + uBolus - btrAvg( uBolus )
+                  normalTransportVelocity(k,iEdge) = edgeMask(k,iEdge) &
+                       *(  barotropicThicknessFlux(iEdge) / thicknessSum &
+                         + normalBaroclinicVelocityNew(k,iEdge) - BtrAvgOfBaroclinicVelocity(iEdge) &
+                         + normalGMBolusVelocity(k,iEdge) - BtrAvgOfGMBolusVelocity(iEdge) )
+               enddo
+
+            end do ! iEdge
+            !$omp end do
+            !$omp end parallel
+end if 
             call mpas_timer_stop('btr se ssh verif')
 
          endif ! split_explicit
@@ -1909,9 +2004,17 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
                ! u = normalBarotropicVelocity + normalBaroclinicVelocity
                ! here normalBaroclinicVelocity is at time n+1/2
                ! This is u used in next iteration or step
+if (.not.config_btr_new3_sub_btrAvgbclV) then
                normalVelocityNew(k,iEdge) = edgeMask(k,iEdge)* &
                          ( normalBarotropicVelocityNew(iEdge) + &
                            normalBaroclinicVelocityNew(k,iEdge) )
+else !if (.not.config_btr_new3_sub_btrAvgbclV) then
+               
+               normalVelocityNew(k,iEdge) = edgeMask(k,iEdge)* &
+                         ( normalBarotropicVelocityNew(iEdge) + &
+                           normalBaroclinicVelocityNew(k,iEdge) - &
+                           BtrAvgOfBaroclinicVelocity(iEdge) )
+end if
 
             enddo
             end do ! iEdge
@@ -2794,12 +2897,7 @@ subroutine ocn_time_integration_split_init(domain)!{{{
          nBtrSubcycles = nBtrSubcycles + 1
       end if
 
-      !*** Set mask for using velocity correction
-      if (config_vel_correction) then
-         useVelocityCorrection = 1.0_RKIND
-      else
-         useVelocityCorrection = 0.0_RKIND
-      endif
+         useVelocityCorrection = 1.0_RKIND ! mrp 
 
       !*** Compute some halo needs to reduce the number
       !*** of halo updates during subcycling
diff --git a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F
index 3850155fb4a3..092df1fa8dfe 100644
--- a/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F
+++ b/components/mpas-ocean/src/shared/mpas_ocn_diagnostics_variables.F
@@ -191,6 +191,8 @@ module ocn_diagnostics_variables
    real (kind=RKIND), dimension(:), pointer :: betaEdge
    real (kind=RKIND), dimension(:), pointer :: RediHorizontalTaper, gmHorizontalTaper, RossbyRadius
 
+   real (kind=RKIND), dimension(:), pointer :: BtrAvgOfBaroclinicVelocity, BtrAvgOfGMBolusVelocity, normalVelocityCorrection
+
    ! Semi-implicit Array Pointers
    real (kind=RKIND), dimension(:), pointer :: SIvec_r0,SIvec_r00,SIvec_r1 ,SIvec_rh0,SIvec_rh1,SIvec_ph0,SIvec_ph1
    real (kind=RKIND), dimension(:), pointer :: SIvec_v0,SIvec_s0 ,SIvec_s1 ,SIvec_sh0,SIvec_sh1
@@ -611,6 +613,13 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e
       call mpas_pool_get_array(diagnosticsPool, 'surfacePressure', &
                 surfacePressure)
 
+      if (trim(config_time_integrator) == 'split_explicit' .or. &
+          trim(config_time_integrator) == 'split_implicit') then
+         call mpas_pool_get_array(diagnosticsPool, 'BtrAvgOfBaroclinicVelocity', BtrAvgOfBaroclinicVelocity)
+         call mpas_pool_get_array(diagnosticsPool, 'BtrAvgOfGMBolusVelocity', BtrAvgOfGMBolusVelocity)
+      end if
+         call mpas_pool_get_array(diagnosticsPool, 'normalVelocityCorrection', normalVelocityCorrection)
+
       ! Semi-implicit Array Pointer retrievals
       if (trim(config_time_integrator) == 'split_implicit') then
          call mpas_pool_get_array(diagnosticsPool, 'SIvec_r0', SIvec_r0)
@@ -801,6 +810,14 @@ subroutine ocn_diagnostics_variables_init(domain, jenkinsOn, hollandJenkinsOn, e
          !$acc                   salinitySurfaceRestoringTendency  &
          !$acc                   )
       end if
+      if (trim(config_time_integrator) == 'split_explicit' .or. &
+          trim(config_time_integrator) == 'split_implicit') then
+         !$acc enter data create(                       &
+         !$acc                   BtrAvgOfBaroclinicVelocity,              &
+         !$acc                   normalVelocityCorrection,              &
+         !$acc                   BtrAvgOfGMBolusVelocity              &
+         !$acc                   )
+      end if
       if ( trim(config_time_integrator) == 'split_implicit' ) then
          !$acc enter data create(                       &
          !$acc                   SIvec_r0,              &
@@ -1049,6 +1066,14 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{
          !$acc                   salinitySurfaceRestoringTendency  &
          !$acc                   )
       end if
+      if (trim(config_time_integrator) == 'split_explicit' .or. &
+          trim(config_time_integrator) == 'split_implicit') then
+         !$acc exit data delete(                        &
+         !$acc                   BtrAvgOfBaroclinicVelocity,              &
+         !$acc                   normalVelocityCorrection,              &
+         !$acc                   BtrAvgOfGMBolusVelocity              &
+         !$acc                   )
+      end if
       if ( trim(config_time_integrator) == 'split_implicit' ) then
          !$acc exit data delete(                        &
          !$acc                   SIvec_r0,              &
@@ -1249,6 +1274,12 @@ subroutine ocn_diagnostics_variables_destroy(err) !{{{
       if ( config_use_self_attraction_loading ) then
          nullify(ssh_sal)
       end if
+      if (trim(config_time_integrator) == 'split_explicit' .or. &
+          trim(config_time_integrator) == 'split_implicit') then
+         nullify(BtrAvgOfBaroclinicVelocity,              &
+                 normalVelocityCorrection,              &
+                 BtrAvgOfGMBolusVelocity)
+      end if
       if ( trim(config_time_integrator) == 'split_implicit' ) then
          nullify(SIvec_r0,              &
                  SIvec_r00,             &

From 35de4211bb3ab69dac8fd3dc36d9c15cabd2eb2c Mon Sep 17 00:00:00 2001
From: Mark Petersen <mpetersen@lanl.gov>
Date: Mon, 19 Sep 2022 10:02:27 -0600
Subject: [PATCH 2/5] turn on gm package (temporary)

---
 components/mpas-ocean/src/driver/mpas_ocn_core_interface.F | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/components/mpas-ocean/src/driver/mpas_ocn_core_interface.F b/components/mpas-ocean/src/driver/mpas_ocn_core_interface.F
index 9d9fd8b19f89..59340c911f99 100644
--- a/components/mpas-ocean/src/driver/mpas_ocn_core_interface.F
+++ b/components/mpas-ocean/src/driver/mpas_ocn_core_interface.F
@@ -356,6 +356,8 @@ function ocn_setup_packages(configPool, packagePool, iocontext) result(ierr)!{{{
       if (config_use_GM.or.config_use_Redi) then
          gmActive = .true.
       end if
+         ! mrp TEMPORARY, always turn on gm package
+         gmActive = .true.
 
       !
       ! test for use of gm

From 5fa93e9cba3c8c74508311d76989f6cc86f85aa6 Mon Sep 17 00:00:00 2001
From: Mark Petersen <mpetersen@lanl.gov>
Date: Mon, 19 Sep 2022 10:13:33 -0600
Subject: [PATCH 3/5] add config_btr_new5_NoALE

---
 components/mpas-ocean/src/Registry.xml        |  4 ++++
 .../mpas-ocean/src/shared/mpas_ocn_tendency.F | 24 +++++++++++++++++++
 2 files changed, 28 insertions(+)

diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml
index d6df61e02de3..17e645d5ec38 100644
--- a/components/mpas-ocean/src/Registry.xml
+++ b/components/mpas-ocean/src/Registry.xml
@@ -1279,6 +1279,10 @@
 					description="experimental flag"
 					possible_values=".true. or .false."
 		/>
+		<nml_option name="config_btr_new5_NoALE" type="logical" default_value=".false."
+					description="experimental flag"
+					possible_values=".true. or .false."
+		/>
 		<nml_option name="config_vel_correction" type="logical" default_value=".true."
 					description="If true, the velocity correction term is included in the horizontal advection of thickness and tracers"
 					possible_values=".true. or .false."
diff --git a/components/mpas-ocean/src/shared/mpas_ocn_tendency.F b/components/mpas-ocean/src/shared/mpas_ocn_tendency.F
index e081df8ae097..189bd8abd116 100644
--- a/components/mpas-ocean/src/shared/mpas_ocn_tendency.F
+++ b/components/mpas-ocean/src/shared/mpas_ocn_tendency.F
@@ -207,8 +207,15 @@ subroutine ocn_tend_thick(tendPool, forcingPool)!{{{
                                layerThickEdgeFlux, tendThick, err)
 
       ! Compute vertical advection term -d/dz(hw)
+if (config_btr_new5_NoALE) then
+      ! force z-level with Eulerian transport at top
+      call ocn_thick_vadv_tend(vertTransportVelocityTop, &
+                               tendThick, err)
+else
+      ! original
       call ocn_thick_vadv_tend(vertAleTransportTop, &
                                tendThick, err)
+endif
 
       ! Compute surface flux tendency
       call ocn_thick_surface_flux_tend(fractionAbsorbed, &
@@ -411,8 +418,15 @@ subroutine ocn_tend_vel(domain, tendPool, statePool, forcingPool, &
                                       kineticEnergyCell, tendVel, err)
 
       ! Add vertical advection term -w du/dz
+if (config_btr_new5_NoALE) then
+      ! force z-level with Eulerian transport at top
+      call ocn_vel_vadv_tend(normalVelocity, layerThickEdgeFlux, &
+                             vertVelocityTop, tendVel, err)
+else
+      ! original
       call ocn_vel_vadv_tend(normalVelocity, layerThickEdgeFlux, &
                              vertAleTransportTop, tendVel, err)
+endif
 
       ! Add pressure gradient
       call ocn_vel_pressure_grad_tend(ssh, pressure, surfacePressure, &
@@ -1095,10 +1109,20 @@ subroutine ocn_tend_tracer(tendPool, statePool, forcingPool, &
          !     -div( layerThickness \phi u)
          ! Tracer budget is computed and stored within the tracer adv
          ! routine
+if (config_btr_new5_NoALE) then
+      ! force z-level with Eulerian transport at top
+      ! I think this is just for diagnostics, but being safe here.
+         call ocn_tracer_advection_tend(tracerGroup, &
+                            normalThicknessFlux, vertTransportVelocityTop, &
+                            layerThickness, dt, tracerGroupTend, &
+                            computeBudgets)
+else
+      ! original
          call ocn_tracer_advection_tend(tracerGroup, &
                             normalThicknessFlux, vertAleTransportTop, &
                             layerThickness, dt, tracerGroupTend, &
                             computeBudgets)
+endif
 
          !$acc update host (activeTracerHorizontalAdvectionEdgeFlux, &
          !$acc              activeTracerHorizontalAdvectionTendency, &

From 830652e4a4fccbbbc2841cb5b6879e5ca84c0c06 Mon Sep 17 00:00:00 2001
From: Mark Petersen <mpetersen@lanl.gov>
Date: Tue, 20 Sep 2022 08:07:48 -0600
Subject: [PATCH 4/5] add config_btr_new6_EdgeThickness

---
 components/mpas-ocean/src/Registry.xml        |  4 ++
 .../mpas_ocn_time_integration_split.F         | 64 +++++++++++++------
 2 files changed, 48 insertions(+), 20 deletions(-)

diff --git a/components/mpas-ocean/src/Registry.xml b/components/mpas-ocean/src/Registry.xml
index 17e645d5ec38..101a7cfa4c1e 100644
--- a/components/mpas-ocean/src/Registry.xml
+++ b/components/mpas-ocean/src/Registry.xml
@@ -1283,6 +1283,10 @@
 					description="experimental flag"
 					possible_values=".true. or .false."
 		/>
+		<nml_option name="config_btr_new6_EdgeThickness" type="logical" default_value=".false."
+					description="experimental flag"
+					possible_values=".true. or .false."
+		/>
 		<nml_option name="config_vel_correction" type="logical" default_value=".true."
 					description="If true, the velocity correction term is included in the horizontal advection of thickness and tracers"
 					possible_values=".true. or .false."
diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F
index c70174b94dfa..92eba83b8912 100644
--- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F
+++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F
@@ -1112,12 +1112,13 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
 
                !$omp parallel
                !$omp do schedule(runtime) &
-               !$omp private(i, iEdge, cell1, cell2, sshEdge, &
+               !$omp private(i, k, iEdge, cell1, cell2, sshEdge, &
                !$omp         thicknessSum, flux)
                do iCell = 1, nCells
                   sshTend(iCell) = 0.0_RKIND
                   do i = 1, nEdgesOnCell(iCell)
                      iEdge = edgesOnCell(i, iCell)
+                     if (maxLevelEdgeTop(iEdge).eq.0) cycle
 
                      cell1 = cellsOnEdge(1, iEdge)
                      cell2 = cellsOnEdge(2, iEdge)
@@ -1129,16 +1130,22 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
                      !thicknessSum = sshEdge + &
                      ! refBottomDepthTopOfCell(maxLevelEdgeTop(iEdge)+1)
 
+if (.not.config_btr_new6_EdgeThickness) then
                      ! method 1, matches method 0 without pbcs, works with pbcs.
                      thicknessSum = sshEdge + min(bottomDepth(cell1), &
                                                   bottomDepth(cell2))
-
+else
+! mrp note, this works for z-level right now.
+! for z-star, could use -zTop(k+1) for refBottomDepth(k)
+! but that would not work for general coordinate
                      ! method 2: may be better than method 1.
                      ! Take average  of full thickness at two
                      ! neighboring cells.
-                     !thicknessSum = sshEdge + 0.5* &
-                     !      (bottomDepth(cell1) + bottomDepth(cell2))
-
+                     k = maxLevelEdgeTop(iEdge)
+                     thicknessSum = sshEdge + 0.5* &
+                        (  min(bottomDepth(cell1), refBottomDepth(k)) &
+                         + min(bottomDepth(cell2), refBottomDepth(k)) )
+endif
 
                      flux = ((1.0-config_btr_gam1_velWt1)* &
                             normalBarotropicVelocitySubcycleCur(iEdge) &
@@ -1178,16 +1185,26 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
                   !$omp do schedule(runtime) &
                   !$omp private(cell1, cell2, sshEdge,thicknessSum,flux)
                   do iEdge = 1, nEdges
+                     if (maxLevelEdgeTop(iEdge).eq.0) cycle
                      cell1 = cellsOnEdge(1,iEdge)
                      cell2 = cellsOnEdge(2,iEdge)
 
                      sshEdge = 0.5_RKIND*(sshSubcycleCur(cell1) + &
                                           sshSubcycleCur(cell2))
 
-                     ! method 1, matches method 0 without pbcs,
-                     !    works with pbcs.
+if (.not.config_btr_new6_EdgeThickness) then
+                     ! method 1, matches method 0 without pbcs, works with pbcs.
                      thicknessSum = sshEdge + min(bottomDepth(cell1), &
                                                   bottomDepth(cell2))
+else
+                     ! method 2: may be better than method 1.
+                     ! Take average  of full thickness at two
+                     ! neighboring cells.
+                     k = maxLevelEdgeTop(iEdge)
+                     thicknessSum = sshEdge + 0.5* &
+                        (  min(bottomDepth(cell1), refBottomDepth(k)) &
+                         + min(bottomDepth(cell2), refBottomDepth(k)) )
+endif
 
                      flux = ((1.0-config_btr_gam1_velWt1)* &
                             normalBarotropicVelocitySubcycleCur(iEdge) &
@@ -1399,6 +1416,7 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
                      sshTend(iCell) = 0.0_RKIND
                      do i = 1, nEdgesOnCell(iCell)
                         iEdge = edgesOnCell(i, iCell)
+                        if (maxLevelEdgeTop(iEdge).eq.0) cycle
 
                         cell1 = cellsOnEdge(1,iEdge)
                         cell2 = cellsOnEdge(2,iEdge)
@@ -1420,17 +1438,21 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
                         !thicknessSum = sshEdge + &
                         !refBottomDepthTopOfCell(maxLevelEdgeTop(iEdge)+1)
 
+if (.not.config_btr_new6_EdgeThickness) then
                         ! method 1, matches method 0 without pbcs,
                         !    but works with pbcs.
                         thicknessSum = sshEdge + &
                                        min(bottomDepth(cell1), &
                                            bottomDepth(cell2))
-
+else
                         ! method 2: may be better than method 1.
                         ! take average of full thickness at two
                         !    neighboring cells
-                        !thicknessSum = sshEdge + 0.5* &
-                        !   (bottomDepth(cell1) + bottomDepth (cell2))
+                     k = maxLevelEdgeTop(iEdge)
+                     thicknessSum = sshEdge + 0.5* &
+                        (  min(bottomDepth(cell1), refBottomDepth(k)) &
+                         + min(bottomDepth(cell2), refBottomDepth(k)) )
+endif
 
                         flux = ((1.0-config_btr_gam3_velWt2)* &
                            normalBarotropicVelocitySubcycleCur(iEdge) &
@@ -1474,17 +1496,19 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
                      !thicknessSum = sshEdge + &
                      !refBottomDepthTopOfCell(maxLevelEdgeTop(iEdge)+1)
 
-                     ! method 1, matches method 0 without pbcs,
-                     !           but works with pbcs.
-                     thicknessSum = sshEdge + min(bottomDepth(cell1),&
+if (.not.config_btr_new6_EdgeThickness) then
+                     ! method 1, matches method 0 without pbcs, works with pbcs.
+                     thicknessSum = sshEdge + min(bottomDepth(cell1), &
                                                   bottomDepth(cell2))
-
-                     ! method 2, better, I think.
-                     !           take average of full thickness at two
-                     !           neighboring cells
-                     !thicknessSum = sshEdge + 0.5* &
-                     !               (bottomDepth(cell1) + &
-                     !                bottomDepth(cell2) )
+else
+                     ! method 2: may be better than method 1.
+                     ! Take average  of full thickness at two
+                     ! neighboring cells.
+                     k = maxLevelEdgeTop(iEdge)
+                     thicknessSum = sshEdge + 0.5* &
+                        (  min(bottomDepth(cell1), refBottomDepth(k)) &
+                         + min(bottomDepth(cell2), refBottomDepth(k)) )
+endif
 
                      flux = ((1.0-config_btr_gam3_velWt2) * &
                          normalBarotropicVelocitySubcycleCur(iEdge) &

From 1a151a591dfe4968b115ba07460f577f5f7eb8a3 Mon Sep 17 00:00:00 2001
From: Mark Petersen <mpetersen@lanl.gov>
Date: Tue, 20 Sep 2022 09:53:58 -0600
Subject: [PATCH 5/5] change edge thickness from refBD to zTop

---
 .../mpas_ocn_time_integration_split.F         | 28 +++++++++++++------
 1 file changed, 20 insertions(+), 8 deletions(-)

diff --git a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F
index 92eba83b8912..d7d4b02dce8c 100644
--- a/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F
+++ b/components/mpas-ocean/src/mode_forward/mpas_ocn_time_integration_split.F
@@ -1143,8 +1143,11 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
                      ! neighboring cells.
                      k = maxLevelEdgeTop(iEdge)
                      thicknessSum = sshEdge + 0.5* &
-                        (  min(bottomDepth(cell1), refBottomDepth(k)) &
-                         + min(bottomDepth(cell2), refBottomDepth(k)) )
+                        (  min(bottomDepth(cell1), -zTop(k+1,cell1)) &
+                         + min(bottomDepth(cell2), -zTop(k+1,cell2)) )
+                        ! first attempt, for z-level only:
+                        !(  min(bottomDepth(cell1), refBottomDepth(k)) &
+                        ! + min(bottomDepth(cell2), refBottomDepth(k)) )
 endif
 
                      flux = ((1.0-config_btr_gam1_velWt1)* &
@@ -1202,8 +1205,11 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
                      ! neighboring cells.
                      k = maxLevelEdgeTop(iEdge)
                      thicknessSum = sshEdge + 0.5* &
-                        (  min(bottomDepth(cell1), refBottomDepth(k)) &
-                         + min(bottomDepth(cell2), refBottomDepth(k)) )
+                        (  min(bottomDepth(cell1), -zTop(k+1,cell1)) &
+                         + min(bottomDepth(cell2), -zTop(k+1,cell2)) )
+                        ! first attempt, for z-level only:
+                        !(  min(bottomDepth(cell1), refBottomDepth(k)) &
+                        ! + min(bottomDepth(cell2), refBottomDepth(k)) )
 endif
 
                      flux = ((1.0-config_btr_gam1_velWt1)* &
@@ -1450,8 +1456,11 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
                         !    neighboring cells
                      k = maxLevelEdgeTop(iEdge)
                      thicknessSum = sshEdge + 0.5* &
-                        (  min(bottomDepth(cell1), refBottomDepth(k)) &
-                         + min(bottomDepth(cell2), refBottomDepth(k)) )
+                        (  min(bottomDepth(cell1), -zTop(k+1,cell1)) &
+                         + min(bottomDepth(cell2), -zTop(k+1,cell2)) )
+                        ! first attempt, for z-level only:
+                        !(  min(bottomDepth(cell1), refBottomDepth(k)) &
+                        ! + min(bottomDepth(cell2), refBottomDepth(k)) )
 endif
 
                         flux = ((1.0-config_btr_gam3_velWt2)* &
@@ -1506,8 +1515,11 @@ subroutine ocn_time_integrator_split(domain, dt)!{{{
                      ! neighboring cells.
                      k = maxLevelEdgeTop(iEdge)
                      thicknessSum = sshEdge + 0.5* &
-                        (  min(bottomDepth(cell1), refBottomDepth(k)) &
-                         + min(bottomDepth(cell2), refBottomDepth(k)) )
+                        (  min(bottomDepth(cell1), -zTop(k+1,cell1)) &
+                         + min(bottomDepth(cell2), -zTop(k+1,cell2)) )
+                        ! first attempt, for z-level only:
+                        !(  min(bottomDepth(cell1), refBottomDepth(k)) &
+                        ! + min(bottomDepth(cell2), refBottomDepth(k)) )
 endif
 
                      flux = ((1.0-config_btr_gam3_velWt2) * &