More precise specific-vol. anomaly for ocean in P (MITgcm#827)

* More precise specific-vol. anomaly for ocean in P

This is a minor change in Geopotential-Anomaly calculation for Ocean in P-coord:
Use a more robust/precise (regarding floating point arithmetics) expression
for specific-volume anomaly.

* move/remove minus sign

* document minor change for Ocean in P-coord

doc/tag-index

@@ -1,6 +1,9 @@
     Notes on tags used in MITgcmUV
     ==============================
 
+o model/src:
+  - more precise specific-volume anomaly expression for Ocean in P-coord;
+    update reference output of all 3 Ocean in P test experiments.
 o pkg/generic_advdiff:
   - adjust 2 SOM (Prather advection) store directives to reduce memory overhead
   - include OS7MP (advScheme 7) in AD code ; switch AD test experiment

model/src/calc_phi_hyd.F

@@ -80,7 +80,7 @@ SUBROUTINE CALC_PHI_HYD(
       _RL pKappaU (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
       _RL pKappaC, locDepth, fullDepth
 #endif /* DISABLE_SIGMA_CODE */
-      _RL thetaRef, locAlpha
+      _RL thetaRef, locBuoy
       _RL dRlocM,dRlocP, ddRloc
       _RL ddPIm, ddPIp, rec_dRm, rec_dRp
       _RL surfPhiFac
@@ -322,12 +322,14 @@ SUBROUTINE CALC_PHI_HYD(
 #ifdef ALLOW_AUTODIFF_TAMC
 CADJ STORE alphaRho = comlev1_bibj_k, key = kkey, kind = isbyte
 #endif
-C--     Calculate specific volume anomaly : alpha_prime = 1/rho - alpha_Cst
+C--     Calculate specific volume anomaly "alpha_prime":
+C       = 1/(rho_prime + rho_Cst) - 1/rho_Cst
+C       = - (rho_prime/rho_Cst) / ( 1 + rho_prime/rho_Cst ) / rho_Cst
         DO j=jMin,jMax
           DO i=iMin,iMax
-            locAlpha=alphaRho(i,j)+rhoConst
-            alphaRho(i,j)=maskC(i,j,k,bi,bj)*
-     &              (oneRL/locAlpha - recip_rhoConst)
+            locBuoy = alphaRho(i,j)*recip_rhoConst
+            alphaRho(i,j) = -maskC(i,j,k,bi,bj)*recip_rhoConst
+     &                    * locBuoy/( oneRL + locBuoy )
           ENDDO
         ENDDO