Change displacedDensity to be relative displacement of one level.

Add potentialDensity for displacement to top level.

src/core_atmosphere/Registry.xml

@@ -1,5 +1,5 @@
 <?xml version="1.0"?>
-<registry model="mpas" core="atmosphere" version="1.0">
+<registry model="mpas" core="atmosphere" version="1.1">
 
 <!-- **************************************************************************************** -->
 <!-- ************************************** Dimensions ************************************** -->

src/core_init_atmosphere/Registry.xml

@@ -1,5 +1,5 @@
 <?xml version="1.0"?>
-<registry model="mpas" core="init_atmosphere" version="1.0">
+<registry model="mpas" core="init_atmosphere" version="1.1">
 
 <!-- **************************************************************************************** -->
 <!-- ************************************** Dimensions ************************************** -->

src/core_ocean/Registry.xml

@@ -1,5 +1,5 @@
 <?xml version="1.0"?>
-<registry model="mpas" core="ocean" version="1.0">
+<registry model="mpas" core="ocean" version="1.1">
 	<dims>
 		<dim name="nCells" units="unitless"
 		     description="The number of polygons in the primary grid."
@@ -683,7 +683,10 @@
 		     description="vertical velocity defined at center (horizonally) and top (vertically) of cell"
 		/>
 		<var name="displacedDensity" type="real" dimensions="nVertLevels nCells Time" streams="o" units="kg m^{-3}"
-		     description="potential density displaced to the mid-depth of top layer"
+		     description="Density displaced adiabatically to the mid-depth one layer deeper.  That is, layer k has been displaced to the depth of layer k+1."
+		/>
+		<var name="potentialDensity" type="real" dimensions="nVertLevels nCells Time" streams="o" units="kg m^{-3}"
+		     description="potential density: density displaced adiabatically to the mid-depth of top layer"
 		/>
 		<var name="BruntVaisalaFreqTop" type="real" dimensions="nVertLevels nCells Time" streams="o" units="s^{-2}"
 		     description="Brunt Vaisala frequency defined at the center (horizontally) and top (vertically) of cell"
@@ -1015,7 +1018,7 @@
 		     description="vertical viscosity defined at the edge (horizontally) and top (vertically)"
 		/>
 		<var name="vertDiffTopOfCell" type="real" dimensions="nVertLevelsP1 nCells Time" units="m^2 s^{-1}"
-		     description="vertical diffusion defined at the edge (horizontally) and top (vertically)"
+		     description="vertical diffusion defined at the cell center (horizontally) and top (vertically)"
 		/>
         <var name="windStressZonal" type="real" dimensions="nCells Time" streams="o" units="N m^{-2}"
              description="reconstructed surface wind stress in the eastward direction."

src/core_ocean/mpas_ocn_diagnostics.F

@@ -110,7 +110,7 @@ subroutine ocn_diagnostic_solve(dt, s, grid, scratch)!{{{
         relativeVorticity, relativeVorticityCell, &
         normalizedPlanetaryVorticityEdge, normalizedPlanetaryVorticityVertex, &
         normalizedRelativeVorticityEdge, normalizedRelativeVorticityVertex, normalizedRelativeVorticityCell, &
-        density, displacedDensity, temperature, salinity, kineticEnergyVertex, kineticEnergyVertexOnCells, uBolusGM, uTransport, &
+        density, displacedDensity, potentialDensity, temperature, salinity, kineticEnergyVertex, kineticEnergyVertexOnCells, uBolusGM, uTransport, &
         vertVelocityTop, BruntVaisalaFreqTop, &
         vorticityGradientNormalComponent, vorticityGradientTangentialComponent
       real (kind=RKIND), dimension(:,:,:), pointer :: tracers, deriv_two
@@ -132,6 +132,7 @@ subroutine ocn_diagnostic_solve(dt, s, grid, scratch)!{{{
       normalizedRelativeVorticityCell    => s % normalizedRelativeVorticityCell % array
       density         => s % density % array
       displacedDensity=> s % displacedDensity % array
+      potentialDensity=> s % potentialDensity % array
       montgomeryPotential     => s % montgomeryPotential % array
       pressure    => s % pressure % array
       zMid        => s % zMid % array
@@ -413,17 +414,21 @@ subroutine ocn_diagnostic_solve(dt, s, grid, scratch)!{{{
          call mpas_timer_start("equation of state", .false., diagEOSTimer)
 
          ! compute in-place density
-         call ocn_equation_of_state_density(s, grid, 0, 'relative', err)
+         call ocn_equation_of_state_density(s, grid, 0, 'relative', density, err)
 
-         ! compute displacedDensity, the potential density referenced to the top layer
-         call ocn_equation_of_state_density(s, grid, 1, 'relative', err)
+         ! compute potentialDensity, the density displaced adiabatically to the mid-depth of top layer.
+         call ocn_equation_of_state_density(s, grid, 1, 'absolute', potentialDensity, err)
+
+         ! compute displacedDensity, density displaced adiabatically to the mid-depth one layer deeper.  
+         ! That is, layer k has been displaced to the depth of layer k+1.
+         call ocn_equation_of_state_density(s, grid, 1, 'relative', displacedDensity, err)
 
          call mpas_timer_stop("equation of state", diagEOSTimer)
       endif
 
       !
       ! Pressure
-      ! This section must be after computing density
+      ! This section must be placed in the code after computing the density.
       !
       if (config_pressure_gradient_type.eq.'MontgomeryPotential') then
 
@@ -491,7 +496,7 @@ subroutine ocn_diagnostic_solve(dt, s, grid, scratch)!{{{
       do iCell=1,nCells
          BruntVaisalaFreqTop(1,iCell) = 0.0
          do k=2,maxLevelCell(iCell)
-            BruntVaisalaFreqTop(k,iCell) = coef * (displacedDensity(k-1,iCell) - displacedDensity(k,iCell)) & 
+            BruntVaisalaFreqTop(k,iCell) = coef * (displacedDensity(k-1,iCell) - density(k,iCell)) & 
               / (zMid(k-1,iCell) - zMid(k,iCell))
           end do
       end do

src/core_ocean/mpas_ocn_equation_of_state.F

@@ -71,7 +71,7 @@ module ocn_equation_of_state
 !
 !-----------------------------------------------------------------------
 
-   subroutine ocn_equation_of_state_density(s, grid, k_displaced, displacement_type, err)!{{{
+   subroutine ocn_equation_of_state_density(s, grid, k_displaced, displacement_type, density, err)!{{{
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !  This module contains routines necessary for computing the density
    !  from model temperature and salinity using an equation of state.
@@ -80,9 +80,9 @@ subroutine ocn_equation_of_state_density(s, grid, k_displaced, displacement_type
    !        s - state: tracers
    !        k_displaced 
    !
-   !  If k_displaced==0, state % density is returned with no displacement 
+   !  If k_displaced==0, density is returned with no displacement 
    !
-   !  If k_displaced~=0,the state % displacedDensity is returned, and is for
+   !  If k_displaced~=0, density is returned, and is for
    !  a parcel adiabatically displaced from its original level to level 
    !  k_displaced.  When using the linear EOS, state % displacedDensity is 
    !  still filled, but depth (i.e. pressure) does not modify the output.
@@ -93,12 +93,12 @@ subroutine ocn_equation_of_state_density(s, grid, k_displaced, displacement_type
 
       type (state_type), intent(inout) :: s
       type (mesh_type), intent(in) :: grid
+      real (kind=RKIND), dimension(:,:), intent(out) :: density
       integer, intent(out) :: err
       integer :: k_displaced
       character(len=*), intent(in) :: displacement_type
 
       integer, dimension(:), pointer :: maxLevelCell
-      real (kind=RKIND), dimension(:,:), pointer :: density
       real (kind=RKIND), dimension(:,:,:), pointer :: tracers
       integer :: nCells, iCell, k, indexT, indexS
       type (dm_info) :: dminfo
@@ -111,13 +111,6 @@ subroutine ocn_equation_of_state_density(s, grid, k_displaced, displacement_type
       indexT = s % index_temperature
       indexS = s % index_salinity
 
-      !  Choose to fill the array density or displacedDensity
-      if (k_displaced == 0) then
-         density => s % density % array
-      else
-         density => s % displacedDensity % array
-      endif
-
       if (linearEos) then
 
          call ocn_equation_of_state_linear_density(grid, indexT, indexS, tracers, density, err)

src/core_ocean/mpas_ocn_equation_of_state_jm.F

@@ -102,9 +102,8 @@ subroutine ocn_equation_of_state_jm_density(grid, k_displaced, displacement_type
 
       real (kind=RKIND), dimension(:), pointer :: &
         refBottomDepth, pRefEOS
-      real (kind=RKIND), dimension(:,:), intent(inout) :: &
-        density
       real (kind=RKIND), dimension(:,:,:), intent(in) :: tracers
+      real (kind=RKIND), dimension(:,:), intent(out) :: density
 
       integer, dimension(:), pointer :: maxLevelCell
 

src/core_ocean/mpas_ocn_vmix_coefs_rich.F

@@ -146,8 +146,14 @@ subroutine ocn_vmix_coefs_rich_build(grid, s, d, err)!{{{
       tracers => s % tracers % array
 
       call mpas_timer_start("eos rich", .false., richEOSTimer)
-      call ocn_equation_of_state_density(s, grid, 0,'relative', err)
-      call ocn_equation_of_state_density(s, grid, 1,'relative', err)
+
+         ! compute in-place density
+         call ocn_equation_of_state_density(s, grid, 0, 'relative', density, err)
+
+         ! compute displacedDensity, density displaced adiabatically to the mid-depth one layer deeper.  
+         ! That is, layer k has been displaced to the depth of layer k+1.
+         call ocn_equation_of_state_density(s, grid, 1, 'relative', displacedDensity, err)
+
       call mpas_timer_stop("eos rich", richEOSTimer)
 
       call ocn_vmix_get_rich_numbers(grid, indexT, indexS, normalVelocity, layerThickness, layerThicknessEdge, & 
@@ -430,11 +436,11 @@ subroutine ocn_vmix_get_rich_numbers(grid, indexT, indexS, normalVelocity, layer
          ddensityTopOfCell(nVertLevels+1,nCells+1), ddensityTopOfEdge(nVertLevels+1,nEdges), &
          du2TopOfCell(nVertLevels+1,nCells+1), du2TopOfEdge(nVertLevels+1,nEdges))
 
-      ! ddensityTopOfCell(k) = $\density^*_{k-1}-\density^*_k$
+      ! ddensityTopOfCell(k) = $\rho^*_{k-1}-\rho_k$, where $\rho^*$ has been adiabatically displaced to level k.
       ddensityTopOfCell = 0.0
       do iCell=1,nCells
          do k=2,maxLevelCell(iCell)
-            ddensityTopOfCell(k,iCell) = displacedDensity(k-1,iCell) - displacedDensity(k,iCell)
+            ddensityTopOfCell(k,iCell) = displacedDensity(k-1,iCell) - density(k,iCell)
           end do
       end do
 

src/core_sw/Registry.xml

@@ -1,5 +1,5 @@
 <?xml version="1.0"?>
-<registry model="mpas" core="sw" version="1.0">
+<registry model="mpas" core="sw" version="1.1">
 	<dims>
 		<dim name="nCells"/>
 		<dim name="nEdges"/>