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