Merge branch 'cbegeman/ocn/add-explicit-layerThickEdge-vars' (PR #4832)

Add explicit layerThicknessEdge variables

This PR replaces the existing layerThicknessEdge variable with two
variables:
* layerThicknessEdgeMean refers to the mean of the neighboring
  cell-centered layerThickness
* layerThicknessEdgeFlux refers to the layerThicknessEdge used to flux
  momentum and tracers.
When config_thickness_flux_type = 'centered' layerThicknessEdgeFlux =
layerThicknessEdgeMean, but when config_thickness_flux_type = 'upwind'
layerThicknessEdgeFlux is upwinded.

Prior to this PR, when config_thickness_flux_type = 'upwind'
layerThicknessEdge was upwinded everywhere. This PR uses
layerThicknessEdgeFlux for horizontal advection terms and
layerThicknessEdgeMean elsewhere, such as for mixing parameterizations.
This change is in preparation for adding additional options for
thickness fluxes.

This PR also removes layerThicknessEdge from global stats, since it did
not offer much information not given by the global stats for
layerThickness.

[BFB]

components/mpas-ocean/src/Registry.xml

@@ -1014,7 +1014,7 @@
 					possible_values=".true. or .false."
 		/>
 		<nml_option name="config_thickness_flux_type" type="character" default_value="centered" units="unitless"
-					description="If 'upwind', use upwind to evaluate the edge-value for layerThickness, i.e., layerThicknessEdge.  The standard MPAS-O approach is 'centered'."
+					description="If 'upwind', use upwind to evaluate the edge-value for layerThickness, i.e., layerThickEdgeFlux.  The standard MPAS-O approach is 'centered'."
 					possible_values="'upwind', 'centered'"
 		/>
 		<nml_option name="config_drying_safety_height" type="real" default_value="1.0e-10" units="m"
@@ -2563,10 +2563,14 @@
 			 description="horizontal velocity, tangential to an edge"
 			 packages="forwardMode;analysisMode"
 		/>
-		<var name="layerThicknessEdge" type="real" dimensions="nVertLevels nEdges Time" units="m"
+		<var name="layerThicknessEdgeMean" type="real" dimensions="nVertLevels nEdges Time" units="m"
 			 description="layer thickness averaged from cell center to edges"
 			 packages="forwardMode;analysisMode"
 		/>
+		<var name="layerThicknessEdgeFlux" type="real" dimensions="nVertLevels nEdges Time" units="m"
+			 description="layer thickness used for fluxes through edges. May be centered, upwinded, or a combination of the two."
+			 packages="forwardMode;analysisMode"
+		/>
 		<var name="layerThicknessVertex" type="real" dimensions="nVertLevels nVertices Time" units="m"
 			 description="layer thickness averaged from cell center to vertices"
 			 packages="forwardMode;analysisMode"

components/mpas-ocean/src/analysis_members/Registry_global_stats.xml

@@ -71,7 +71,7 @@
 				 description="Minimum global value of tangentialVelocity on ocean edges."
 			/>
 			<var name="layerThicknessEdgeMin" array_group="mins" units="m"
-				 description="Minimum global value of layerThicknessEdge on ocean edges."
+				 description="Minimum global value of layerThicknessEdgeMean on ocean edges."
 			/>
 			<var name="relativeVorticityMin" array_group="mins" units="s^{-1}"
 				 description="Minimum global value of relativeVorticity on ocean vertices."
@@ -170,7 +170,7 @@
 				 description="Maximum global value of tangentialVelocity on ocean edges."
 			/>
 			<var name="layerThicknessEdgeMax" array_group="maxes" units="m"
-				 description="Maximum global value of layerThicknessEdge on ocean edges."
+				 description="Maximum global value of layerThicknessEdgeMean on ocean edges."
 			/>
 			<var name="relativeVorticityMax" array_group="maxes" units="s^{-1}"
 				 description="Maximum global value of relativeVorticity on ocean vertices."
@@ -269,7 +269,7 @@
 				 description="Accumulated global value of tangentialVelocity on ocean edges."
 			/>
 			<var name="layerThicknessEdgeSum" array_group="sums" units="m^3"
-				 description="Accumulated global value of layerThicknessEdge on ocean edges."
+				 description="Accumulated global value of layerThicknessEdgeMean on ocean edges."
 			/>
 			<var name="relativeVorticitySum" array_group="sums" units="m^2 s^{-1}"
 				 description="Accumulated global value of relativeVorticity on ocean vertices."
@@ -368,7 +368,7 @@
 				 description="Global root mean square value of tangentialVelocity on ocean edges."
 			/>
 			<var name="layerThicknessEdgeRms" array_group="rms" units="m"
-				 description="Global root mean square value of layerThicknessEdge on ocean edges."
+				 description="Global root mean square value of layerThicknessEdgeMean on ocean edges."
 			/>
 			<var name="relativeVorticityRms" array_group="rms" units="s^{-1}"
 				 description="Global root mean square value of relativeVorticity on ocean vertices."
@@ -467,7 +467,7 @@
 				 description="Average value of tangentialVelocity on ocean edges."
 			/>
 			<var name="layerThicknessEdgeAvg" array_group="avg" units="m"
-				 description="Average value of layerThicknessEdge on ocean edges."
+				 description="Average value of layerThicknessEdgeMean on ocean edges."
 			/>
 			<var name="relativeVorticityAvg" array_group="avg" units="s^{-1}"
 				 description="Average value of relativeVorticity on ocean vertices."
@@ -566,7 +566,7 @@
 				 description="Minimum vertical sum of tangentialVelocity on ocean edges."
 			/>
 			<var name="layerThicknessEdgeMinVertSum" array_group="vertSumMin" units="m"
-				 description="Minimum vertical sum of layerThicknessEdge on ocean edges."
+				 description="Minimum vertical sum of layerThicknessEdgeMean on ocean edges."
 			/>
 			<var name="relativeVorticityMinVertSum" array_group="vertSumMin" units="s^{-1}"
 				 description="Minimum vertical sum of relativeVorticity on ocean vertices."
@@ -665,7 +665,7 @@
 				 description="Maximum vertical sum of tangentialVelocity on ocean edges."
 			/>
 			<var name="layerThicknessEdgeMaxVertSum" array_group="vertSumMax" units="m"
-				 description="Maximum vertical sum of layerThicknessEdge on ocean edges."
+				 description="Maximum vertical sum of layerThicknessEdgeMean on ocean edges."
 			/>
 			<var name="relativeVorticityMaxVertSum" array_group="vertSumMax" units="s^{-1}"
 				 description="Maximum vertical sum of relativeVorticity on ocean vertices."

components/mpas-ocean/src/analysis_members/mpas_ocn_global_stats.F

@@ -125,7 +125,7 @@ subroutine ocn_init_global_stats(domain, err)!{{{
             write (fileID,'(i5,a)') i,'. layerThickness'; i=i+1
             write (fileID,'(i5,a)') i,'. normalVelocity'; i=i+1
             write (fileID,'(i5,a)') i,'. tangentialVelocity'; i=i+1
-            write (fileID,'(i5,a)') i,'. layerThicknessEdge'; i=i+1
+            write (fileID,'(i5,a)') i,'. layerThicknessEdgeMean'; i=i+1
             write (fileID,'(i5,a)') i,'. relativeVorticity'; i=i+1
             write (fileID,'(i5,a)') i,'. enstrophy = relativeVorticity**2'; i=i+1
             write (fileID,'(i5,a)') i,'. kineticEnergyCell'; i=i+1
@@ -409,7 +409,7 @@ subroutine ocn_compute_global_stats(domain, timeLevel, err)!{{{
          variableIndex = variableIndex + 1
          call ocn_compute_field_volume_weighted_local_stats_max_level(dminfo, nVertLevels, nEdgesSolve, &
                               minLevelEdgeBot(1:nEdgesSolve), maxLevelEdgeTop(1:nEdgesSolve), areaEdge(1:nEdgesSolve), &
-                              layerThickEdge(:,1:nEdgesSolve), &
+                              layerThickEdgeMean(:,1:nEdgesSolve), &
                               normalVelocity(:,1:nEdgesSolve), sums_tmp(variableIndex), &
                               sumSquares_tmp(variableIndex), mins_tmp(variableIndex), &
                               maxes_tmp(variableIndex), &
@@ -426,7 +426,7 @@ subroutine ocn_compute_global_stats(domain, timeLevel, err)!{{{
          variableIndex = variableIndex + 1
          call ocn_compute_field_volume_weighted_local_stats_max_level(dminfo, nVertLevels, nEdgesSolve, &
                               minLevelEdgeBot(1:nEdgesSolve), maxLevelEdgeTop(1:nEdgesSolve), areaEdge(1:nEdgesSolve), &
-                              layerThickEdge(:,1:nEdgesSolve), &
+                              layerThickEdgeMean(:,1:nEdgesSolve), &
                               tangentialVelocity(:,1:nEdgesSolve), &
                               sums_tmp(variableIndex), sumSquares_tmp(variableIndex), &
                               mins_tmp(variableIndex), maxes_tmp(variableIndex), &
@@ -439,11 +439,11 @@ subroutine ocn_compute_global_stats(domain, timeLevel, err)!{{{
          verticalSumMins(variableIndex) = min(verticalSumMins(variableIndex), verticalSumMins_tmp(variableIndex))
          verticalSumMaxes(variableIndex) = max(verticalSumMaxes(variableIndex), verticalSumMaxes_tmp(variableIndex))
 
-         ! layerThicknessEdge
+         ! layerThicknessEdgeMean
          variableIndex = variableIndex + 1
          call ocn_compute_field_area_weighted_local_stats_max_level(dminfo, nVertLevels, nEdgesSolve, &
                             minLevelEdgeBot(1:nEdgesSolve), maxLevelEdgeTop(1:nEdgesSolve), areaEdge(1:nEdgesSolve), &
-                            layerThickEdge(:,1:nEdgesSolve), sums_tmp(variableIndex), &
+                            layerThickEdgeMean(:,1:nEdgesSolve), sums_tmp(variableIndex), &
                             sumSquares_tmp(variableIndex), mins_tmp(variableIndex), &
                             maxes_tmp(variableIndex), verticalSumMins_tmp(variableIndex), &
                             verticalSumMaxes_tmp(variableIndex))
@@ -514,7 +514,7 @@ subroutine ocn_compute_global_stats(domain, timeLevel, err)!{{{
          variableIndex = variableIndex + 1
          call ocn_compute_field_volume_weighted_local_stats_max_level(dminfo, nVertLevels, nEdgesSolve, &
                               minLevelEdgeBot(1:nEdgesSolve), maxLevelEdgeTop(1:nEdgesSolve), areaEdge(1:nEdgesSolve), &
-                              layerThickEdge(:,1:nEdgesSolve), &
+                              layerThickEdgeMean(:,1:nEdgesSolve), &
                               normalizedAbsoluteVorticity(:,1:nEdgesSolve), &
                               sums_tmp(variableIndex), sumSquares_tmp(variableIndex), &
                               mins_tmp(variableIndex), maxes_tmp(variableIndex), &
@@ -1128,7 +1128,7 @@ subroutine ocn_compute_global_stats(domain, timeLevel, err)!{{{
       averages(variableIndex) = sums(variableIndex)/volumeEdgeGlobal
       rms(variableIndex) = sqrt(sumSquares(variableIndex)/volumeEdgeGlobal)
 
-      ! layerThicknessEdge
+      ! layerThicknessEdgeMean
       variableIndex = variableIndex + 1
       averages(variableIndex) = sums(variableIndex)/(areaEdgeGlobal*nVertLevels)
       rms(variableIndex) = sqrt(sumSquares(variableIndex)/(areaEdgeGlobal*nVertLevels))

components/mpas-ocean/src/analysis_members/mpas_ocn_meridional_heat_transport.F

@@ -460,7 +460,7 @@ subroutine ocn_compute_meridional_heat_transport(domain, timeLevel, err)!{{{
 
                   do k = 1, kMax
 
-                     ! Compute divergence of huT, i.e. layerThicknessEdge * normalTransportVelocity * temperature, at an edge
+                     ! Compute divergence of huT, i.e. layerThicknessEdgeFlux * normalTransportVelocity * temperature, at an edge
                      ! for meridional heat transport.  Here we use a centered difference to compute the temperature at
                      ! the edge, which is an approximation to the actual edge temperature used in the horizontal
                      ! advection scheme (for example, FCT).  We expect that the error in this approximation is small.
@@ -470,7 +470,7 @@ subroutine ocn_compute_meridional_heat_transport(domain, timeLevel, err)!{{{
                      div_huT = 0.0_RKIND
                      do i = 1, nEdgesOnCell(iCell)
                         iEdge = edgesOnCell(i, iCell)
-                        div_huT = div_huT - layerThickEdge(k, iEdge) * normalTransportVelocity(k, iEdge) &
+                        div_huT = div_huT - layerThickEdgeFlux(k, iEdge) * normalTransportVelocity(k, iEdge) &
                              * 0.5_RKIND * (activeTracers(indexTemperature,k,cellsOnEdge(1,iEdge)) &
                              + activeTracers(indexTemperature,k,cellsOnEdge(2,iEdge))) &
                              * edgeSignOnCell(i, iCell) * dvEdge(iEdge)

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

@@ -887,12 +887,12 @@ subroutine ocn_time_integrator_rk4_compute_vel_tends(block, dt, &
       ! advection of u uses u, while advection of layerThickness and tracers use normalTransportVelocity.
       if (associated(highFreqThicknessProvis)) then
          call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, &
-            layerThicknessCur,layerThickEdge, normalVelocityProvis, &
+            layerThicknessCur,layerThickEdgeFlux, normalVelocityProvis, &
             sshCur, rkSubstepWeight, &
             vertAleTransportTop, err, highFreqThicknessProvis)
       else
          call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, &
-            layerThicknessCur,layerThickEdge, normalVelocityProvis, &
+            layerThicknessCur,layerThickEdgeFlux, normalVelocityProvis, &
             sshCur, rkSubstepWeight, &
             vertAleTransportTop, err)
       endif
@@ -955,12 +955,12 @@ subroutine ocn_time_integrator_rk4_compute_thick_tends(block, dt, rkSubstepWeigh
       ! advection of u uses u, while advection of layerThickness and tracers use normalTransportVelocity.
       if (associated(highFreqThicknessProvis)) then
          call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, &
-            layerThicknessCur, layerThickEdge, normalTransportVelocity, &
+            layerThicknessCur, layerThickEdgeFlux, normalTransportVelocity, &
             sshCur, rkSubstepWeight, &
             vertAleTransportTop, err, highFreqThicknessProvis)
       else
          call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, &
-            layerThicknessCur, layerThickEdge, normalTransportVelocity, &
+            layerThicknessCur, layerThickEdgeFlux, normalTransportVelocity, &
             sshCur, rkSubstepWeight, &
             vertAleTransportTop, err)
       endif
@@ -1025,12 +1025,12 @@ subroutine ocn_time_integrator_rk4_compute_tracer_tends(block, dt, rkSubstepWeig
       ! advection of u uses u, while advection of layerThickness and tracers use normalTransportVelocity.
       if (associated(highFreqThicknessProvis)) then
          call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, &
-            layerThicknessCur, layerThickEdge, normalTransportVelocity, &
+            layerThicknessCur, layerThickEdgeFlux, normalTransportVelocity, &
             sshCur, rkSubstepWeight, &
             vertAleTransportTop, err, highFreqThicknessProvis)
       else
          call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, &
-            layerThicknessCur, layerThickEdge, normalTransportVelocity, &
+            layerThicknessCur, layerThickEdgeFlux, normalTransportVelocity, &
             sshCur, rkSubstepWeight, &
             vertAleTransportTop, err)
       endif
@@ -1086,12 +1086,12 @@ subroutine ocn_time_integrator_rk4_compute_tends(block, dt, rkWeight, err)!{{{
       ! advection of u uses u, while advection of layerThickness and tracers use normalTransportVelocity.
       if (associated(highFreqThicknessProvis)) then
          call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, &
-            layerThicknessCur,layerThickEdge, normalVelocityProvis, &
+            layerThicknessCur,layerThickEdgeFlux, normalVelocityProvis, &
             sshCur, rkWeight, &
             vertAleTransportTop, err, highFreqThicknessProvis)
       else
          call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, &
-            layerThicknessCur,layerThickEdge, normalVelocityProvis, &
+            layerThicknessCur,layerThickEdgeFlux, normalVelocityProvis, &
             sshCur, rkWeight, &
             vertAleTransportTop, err)
       endif
@@ -1100,12 +1100,12 @@ subroutine ocn_time_integrator_rk4_compute_tends(block, dt, rkWeight, err)!{{{
 
       if (associated(highFreqThicknessProvis)) then
          call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, &
-            layerThicknessCur, layerThickEdge, normalTransportVelocity, &
+            layerThicknessCur, layerThickEdgeFlux, normalTransportVelocity, &
             sshCur, rkWeight, &
             vertAleTransportTop, err, highFreqThicknessProvis)
       else
          call ocn_vert_transport_velocity_top(meshPool, verticalMeshPool, &
-            layerThicknessCur, layerThickEdge, normalTransportVelocity, &
+            layerThicknessCur, layerThickEdgeFlux, normalTransportVelocity, &
             sshCur, rkWeight, &
             vertAleTransportTop, err)
       endif