model/src/calc_phi_hyd.F

9bf3441755 Jean*0001 #include "PACKAGES_CONFIG.h"
1dbaea09ee Chri*0002 #include "CPP_OPTIONS.h"
517dbdc414 Jean*0003 #ifdef ALLOW_AUTODIFF
                0004 # include "AUTODIFF_OPTIONS.h"
                0005 #endif
924557e60a Chri*0006
9366854e02 Chri*0007 CBOP
                0008 C     !ROUTINE: CALC_PHI_HYD
                0009 C     !INTERFACE:
fb481a83c2 Alis*0010       SUBROUTINE CALC_PHI_HYD(
f1f873e18c Jean*0011      I                         bi, bj, iMin, iMax, jMin, jMax, k,
                0012      U                         phiHydF,
                0013      O                         phiHydC, dPhiHydX, dPhiHydY,
842349fcb7 Jean*0014      I                         myTime, myIter, myThid )
9366854e02 Chri*0015 C     !DESCRIPTION: \bv
                0016 C     *==========================================================*
11734d097d Chri*0017 C     | SUBROUTINE CALC_PHI_HYD                                  |
0f5ba23f97 Jean*0018 C     | o Integrate the hydrostatic relation to find the Hydros. |
9366854e02 Chri*0019 C     *==========================================================*
f1f873e18c Jean*0020 C     |    Potential (ocean: Pressure/rho ; atmos = geopotential)
                0021 C     | On entry:
                0022 C     |   phiHydF(i,j) is the hydrostatic Potential anomaly
0f5ba23f97 Jean*0023 C     |                at middle between tracer points k-1,k
f1f873e18c Jean*0024 C     | On exit:
                0025 C     |   phiHydC(i,j) is the hydrostatic Potential anomaly
                0026 C     |                at cell centers (tracer points), level k
                0027 C     |   phiHydF(i,j) is the hydrostatic Potential anomaly
0f5ba23f97 Jean*0028 C     |                at middle between tracer points k,k+1
f1f873e18c Jean*0029 C     |   dPhiHydX,Y   hydrostatic Potential gradient (X&Y dir)
                0030 C     |                at cell centers (tracer points), level k
                0031 C     | integr_GeoPot allows to select one integration method
                0032 C     |    1= Finite volume form ; else= Finite difference form
9366854e02 Chri*0033 C     *==========================================================*
                0034 C     \ev
                0035 C     !USES:
924557e60a Chri*0036       IMPLICIT NONE
                0037 C     == Global variables ==
                0038 #include "SIZE.h"
                0039 #include "GRID.h"
81bc00c2f0 Chri*0040 #include "EEPARAMS.h"
924557e60a Chri*0041 #include "PARAMS.h"
7c50f07931 Mart*0042 #ifdef ALLOW_AUTODIFF_TAMC
612f1eab29 Patr*0043 #include "tamc.h"
7c50f07931 Mart*0044 #endif /* ALLOW_AUTODIFF_TAMC */
9e5f2093d3 Alis*0045 #include "SURFACE.h"
60c223928f Mart*0046 #include "DYNVARS.h"
612f1eab29 Patr*0047
9366854e02 Chri*0048 C     !INPUT/OUTPUT PARAMETERS:
924557e60a Chri*0049 C     == Routine arguments ==
0f5ba23f97 Jean*0050 C     bi, bj, k  :: tile and level indices
f1f873e18c Jean*0051 C     iMin,iMax,jMin,jMax :: computational domain
                0052 C     phiHydF    :: hydrostatic potential anomaly at middle between
                0053 C                   2 centers (entry: Interf_k ; output: Interf_k+1)
                0054 C     phiHydC    :: hydrostatic potential anomaly at cell center
                0055 C     dPhiHydX,Y :: gradient (X & Y dir.) of hydrostatic potential anom.
                0056 C     myTime     :: current time
                0057 C     myIter     :: current iteration number
                0058 C     myThid     :: thread number for this instance of the routine.
                0059       INTEGER bi,bj,iMin,iMax,jMin,jMax,k
                0060       _RL phiHydF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0061       _RL phiHydC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
2b8326b7dd Jean*0062       _RL dPhiHydX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0063       _RL dPhiHydY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
7e1bd93d4f Jean*0064       _RL myTime
                0065       INTEGER myIter, myThid
0f5ba23f97 Jean*0066
1dbaea09ee Chri*0067 #ifdef INCLUDE_PHIHYD_CALCULATION_CODE
                0068
9366854e02 Chri*0069 C     !LOCAL VARIABLES:
fb481a83c2 Alis*0070 C     == Local variables ==
3d9080e0a5 Jean*0071 C     phiHydU    :: hydrostatic potential anomaly at interface k+1 (Upper k)
                0072 C     pKappaF    :: (p/Po)^kappa at interface k
                0073 C     pKappaU    :: (p/Po)^kappa at interface k+1 (Upper k)
                0074 C     pKappaC    :: (p/Po)^kappa at cell center k
f1f873e18c Jean*0075       INTEGER i,j
fb481a83c2 Alis*0076       _RL alphaRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
3d9080e0a5 Jean*0077 #ifndef DISABLE_SIGMA_CODE
                0078       _RL phiHydU (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0079       _RL pKappaF (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0080       _RL pKappaU (1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0081       _RL pKappaC, locDepth, fullDepth
                0082 #endif /* DISABLE_SIGMA_CODE */
8fabf86828 Jean*0083       _RL thetaRef, locBuoy
cfbe180681 Jean*0084       _RL dRlocM,dRlocP, ddRloc
f1f873e18c Jean*0085       _RL ddPIm, ddPIp, rec_dRm, rec_dRp
                0086       _RL surfPhiFac
                0087       LOGICAL useDiagPhiRlow, addSurfPhiAnom
3d9080e0a5 Jean*0088       LOGICAL useFVgradPhi
d9efc637ba Mart*0089 #ifdef ALLOW_AUTODIFF_TAMC
                0090 C     tkey :: tape key (tile dependent)
                0091 C     kkey :: tape key (level and tile dependent)
                0092       INTEGER tkey, kkey
                0093 #endif
9366854e02 Chri*0094 CEOP
1adfdf6281 Jean*0095       useDiagPhiRlow = .TRUE.
2b8326b7dd Jean*0096       addSurfPhiAnom = select_rStar.EQ.0 .AND. nonlinFreeSurf.GE.4
3d9080e0a5 Jean*0097       useFVgradPhi   = selectSigmaCoord.NE.0
                0098
f1f873e18c Jean*0099       surfPhiFac = 0.
                0100       IF (addSurfPhiAnom) surfPhiFac = 1.
45adaae56c Jean*0101
                0102 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
0f5ba23f97 Jean*0103 C  Atmosphere:
463053c692 Jean*0104 C integr_GeoPot => select one option for the integration of the Geopotential:
f1f873e18c Jean*0105 C   = 0 : Energy Conserving Form, accurate with Topo full cell;
                0106 C   = 1 : Finite Volume Form, with Part-Cell, linear in P by Half level;
0f5ba23f97 Jean*0107 C   =2,3: Finite Difference Form, with Part-Cell,
f1f873e18c Jean*0108 C         linear in P between 2 Tracer levels.
0f5ba23f97 Jean*0109 C       can handle both cases: Tracer lev at the middle of InterFace_W
f1f873e18c Jean*0110 C                          and InterFace_W at the middle of Tracer lev;
45adaae56c Jean*0111 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
fb481a83c2 Alis*0112
612f1eab29 Patr*0113 #ifdef ALLOW_AUTODIFF_TAMC
edb6656069 Mart*0114       tkey = bi + (bj-1)*nSx + (ikey_dynamics-1)*nSx*nSy
7448700841 Mart*0115       kkey = k + (tkey-1)*Nr
612f1eab29 Patr*0116 #endif /* ALLOW_AUTODIFF_TAMC */
                0117
0f5ba23f97 Jean*0118 C--   Initialize phiHydF to zero :
f1f873e18c Jean*0119 C     note: atmospheric_loading or Phi_topo anomaly are incorporated
                0120 C           later in S/R calc_grad_phi_hyd
                0121       IF (k.EQ.1) THEN
2b8326b7dd Jean*0122         DO j=1-OLy,sNy+OLy
                0123          DO i=1-OLx,sNx+OLx
f1f873e18c Jean*0124            phiHydF(i,j) = 0.
                0125          ENDDO
                0126         ENDDO
                0127       ENDIF
7e1bd93d4f Jean*0128
                0129 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
f1f873e18c Jean*0130       IF ( buoyancyRelation .EQ. 'OCEANIC' ) THEN
fb481a83c2 Alis*0131 C       This is the hydrostatic pressure calculation for the Ocean
                0132 C       which uses the FIND_RHO() routine to calculate density
                0133 C       before integrating g*rho over the current layer/interface
842349fcb7 Jean*0134 #ifdef ALLOW_AUTODIFF_TAMC
7e1bd93d4f Jean*0135 CADJ GENERAL
842349fcb7 Jean*0136 #endif /* ALLOW_AUTODIFF_TAMC */
fb481a83c2 Alis*0137
842349fcb7 Jean*0138         IF ( implicitIntGravWave .OR. myIter.LT.0 ) THEN
f1f873e18c Jean*0139 C---    Calculate density
612f1eab29 Patr*0140 #ifdef ALLOW_AUTODIFF_TAMC
7448700841 Mart*0141 CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, kind = isbyte
                0142 CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, kind = isbyte
612f1eab29 Patr*0143 #endif /* ALLOW_AUTODIFF_TAMC */
842349fcb7 Jean*0144           CALL FIND_RHO_2D(
                0145      I              iMin, iMax, jMin, jMax, k,
b98355c8c5 jm-c 0146      I              theta(1-OLx,1-OLy,k,bi,bj),
                0147      I              salt(1-OLx,1-OLy,k,bi,bj),
842349fcb7 Jean*0148      O              alphaRho,
                0149      I              k, bi, bj, myThid )
                0150         ELSE
                0151           DO j=jMin,jMax
                0152            DO i=iMin,iMax
                0153              alphaRho(i,j) = rhoInSitu(i,j,k,bi,bj)
                0154            ENDDO
                0155           ENDDO
                0156         ENDIF
0f5ba23f97 Jean*0157
a6cbc7a360 Mart*0158 #ifdef ALLOW_SHELFICE
0f5ba23f97 Jean*0159 C     mask rho, so that there is no contribution of phiHyd from
a6cbc7a360 Mart*0160 C     overlying shelfice (whose density we do not know)
807d4a3bf5 Jean*0161         IF ( useShelfIce .AND. useDOWN_SLOPE ) THEN
                0162 C- note: does not work for down_slope pkg which needs rho below the bottom.
                0163 C    setting rho=0 above the ice-shelf base is enough (and works in both cases)
                0164 C    but might be slower (--> keep original masking if not using down_slope pkg)
                0165          DO j=jMin,jMax
                0166           DO i=iMin,iMax
                0167            IF ( k.LT.kSurfC(i,j,bi,bj) ) alphaRho(i,j) = 0. _d 0
                0168           ENDDO
                0169          ENDDO
                0170         ELSEIF ( useShelfIce ) THEN
a6cbc7a360 Mart*0171          DO j=jMin,jMax
                0172           DO i=iMin,iMax
                0173            alphaRho(i,j) = alphaRho(i,j)*maskC(i,j,k,bi,bj)
                0174           ENDDO
                0175          ENDDO
                0176         ENDIF
                0177 #endif /* ALLOW_SHELFICE */
fb481a83c2 Alis*0178
af53f2701c Jean*0179 #ifdef ALLOW_MOM_COMMON
0d373cd03b Jean*0180 C--     Quasi-hydrostatic terms are added in as if they modify the buoyancy
ee19c4a296 Jean*0181         IF ( quasiHydrostatic ) THEN
                0182          CALL MOM_QUASIHYDROSTATIC( bi, bj, k, uVel, vVel,
                0183      U                              alphaRho,
                0184      I                              myTime, myIter, myThid )
775ad00e06 Alis*0185         ENDIF
af53f2701c Jean*0186 #endif /* ALLOW_MOM_COMMON */
7448700841 Mart*0187 #ifdef ALLOW_AUTODIFF_TAMC
                0188 CADJ STORE alphaRho = comlev1_bibj_k, key = kkey, kind = isbyte
                0189 #endif
775ad00e06 Alis*0190
f1f873e18c Jean*0191 #ifdef NONLIN_FRSURF
2b8326b7dd Jean*0192         IF ( addSurfPhiAnom .AND.
                0193      &       uniformFreeSurfLev .AND. k.EQ.1 ) THEN
f1f873e18c Jean*0194           DO j=jMin,jMax
                0195             DO i=iMin,iMax
                0196               phiHydF(i,j) = surfPhiFac*etaH(i,j,bi,bj)
                0197      &                      *gravity*alphaRho(i,j)*recip_rhoConst
                0198             ENDDO
                0199           ENDDO
                0200         ENDIF
                0201 #endif /* NONLIN_FRSURF */
1adfdf6281 Jean*0202
f1f873e18c Jean*0203 C----  Hydrostatic pressure at cell centers
7e1bd93d4f Jean*0204
                0205        IF (integr_GeoPot.EQ.1) THEN
                0206 C  --  Finite Volume Form
                0207
                0208 C---------- This discretization is the "finite volume" form
                0209 C           which has not been used to date since it does not
                0210 C           conserve KE+PE exactly even though it is more natural
2b8326b7dd Jean*0211
                0212         IF ( uniformFreeSurfLev ) THEN
                0213          DO j=jMin,jMax
                0214           DO i=iMin,iMax
                0215             phiHydC(i,j) = phiHydF(i,j)
bd053599dc Jean*0216      &              + halfRL*drF(k)*gravFacC(k)*gravity
                0217      &                             *alphaRho(i,j)*recip_rhoConst
2b8326b7dd Jean*0218             phiHydF(i,j) = phiHydF(i,j)
bd053599dc Jean*0219      &                     + drF(k)*gravFacC(k)*gravity
                0220      &                             *alphaRho(i,j)*recip_rhoConst
7e1bd93d4f Jean*0221           ENDDO
                0222          ENDDO
2b8326b7dd Jean*0223         ELSE
                0224          DO j=jMin,jMax
                0225           DO i=iMin,iMax
                0226            IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
                0227             ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
                0228 #ifdef NONLIN_FRSURF
                0229             ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
                0230 #endif
bd053599dc Jean*0231             phiHydC(i,j) = ddRloc*gravFacC(k)*gravity
                0232      &                           *alphaRho(i,j)*recip_rhoConst
2b8326b7dd Jean*0233            ELSE
                0234             phiHydC(i,j) = phiHydF(i,j)
bd053599dc Jean*0235      &              + halfRL*drF(k)*gravFacC(k)*gravity
                0236      &                             *alphaRho(i,j)*recip_rhoConst
2b8326b7dd Jean*0237            ENDIF
bd053599dc Jean*0238            phiHydF(i,j) = phiHydC(i,j)
                0239      &              + halfRL*drF(k)*gravFacC(k)*gravity
                0240      &                             *alphaRho(i,j)*recip_rhoConst
2b8326b7dd Jean*0241           ENDDO
                0242          ENDDO
                0243         ENDIF
7e1bd93d4f Jean*0244
                0245        ELSE
                0246 C  --  Finite Difference Form
                0247
2b8326b7dd Jean*0248 C---------- This discretization is the "energy conserving" form
                0249 C           which has been used since at least Adcroft et al., MWR 1997
f1f873e18c Jean*0250
bd053599dc Jean*0251         dRlocM = halfRL*drC(k)*gravFacF(k)
                0252         IF (k.EQ.1) dRlocM = (rF(k)-rC(k))*gravFacF(k)
2b8326b7dd Jean*0253         IF (k.EQ.Nr) THEN
bd053599dc Jean*0254           dRlocP = (rC(k)-rF(k+1))*gravFacF(k+1)
2b8326b7dd Jean*0255         ELSE
bd053599dc Jean*0256           dRlocP = halfRL*drC(k+1)*gravFacF(k+1)
2b8326b7dd Jean*0257         ENDIF
                0258         IF ( uniformFreeSurfLev ) THEN
7e1bd93d4f Jean*0259          DO j=jMin,jMax
                0260           DO i=iMin,iMax
2b8326b7dd Jean*0261             phiHydC(i,j) = phiHydF(i,j)
bd053599dc Jean*0262      &             + dRlocM*gravity*alphaRho(i,j)*recip_rhoConst
2b8326b7dd Jean*0263             phiHydF(i,j) = phiHydC(i,j)
bd053599dc Jean*0264      &             + dRlocP*gravity*alphaRho(i,j)*recip_rhoConst
fb481a83c2 Alis*0265           ENDDO
7e1bd93d4f Jean*0266          ENDDO
2b8326b7dd Jean*0267         ELSE
0d373cd03b Jean*0268          rec_dRm = oneRL/(rF(k)-rC(k))
                0269          rec_dRp = oneRL/(rC(k)-rF(k+1))
2b8326b7dd Jean*0270          DO j=jMin,jMax
                0271           DO i=iMin,iMax
                0272            IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
                0273             ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
                0274 #ifdef NONLIN_FRSURF
                0275             ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
                0276 #endif
0d373cd03b Jean*0277             phiHydC(i,j) =( MAX(zeroRL,ddRloc)*rec_dRm*dRlocM
                0278      &                     +MIN(zeroRL,ddRloc)*rec_dRp*dRlocP
                0279      &                    )*gravity*alphaRho(i,j)*recip_rhoConst
2b8326b7dd Jean*0280            ELSE
                0281             phiHydC(i,j) = phiHydF(i,j)
bd053599dc Jean*0282      &             + dRlocM*gravity*alphaRho(i,j)*recip_rhoConst
2b8326b7dd Jean*0283            ENDIF
bd053599dc Jean*0284            phiHydF(i,j) = phiHydC(i,j)
                0285      &             + dRlocP*gravity*alphaRho(i,j)*recip_rhoConst
2b8326b7dd Jean*0286           ENDDO
                0287          ENDDO
                0288         ENDIF
7e1bd93d4f Jean*0289
                0290 C  --  end if integr_GeoPot = ...
                0291        ENDIF
0f5ba23f97 Jean*0292
7e1bd93d4f Jean*0293 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
f1f873e18c Jean*0294       ELSEIF ( buoyancyRelation .EQ. 'OCEANICP' ) THEN
9e5f2093d3 Alis*0295 C       This is the hydrostatic pressure calculation for the Ocean
ff02675122 Jean*0296 C       which uses the FIND_RHO() routine to calculate density before
                0297 C       integrating (1/rho)_prime*dp over the current layer/interface
e305438401 Mart*0298 #ifdef      ALLOW_AUTODIFF_TAMC
                0299 CADJ GENERAL
                0300 #endif      /* ALLOW_AUTODIFF_TAMC */
9e5f2093d3 Alis*0301
842349fcb7 Jean*0302         IF ( implicitIntGravWave .OR. myIter.LT.0 ) THEN
1adfdf6281 Jean*0303 C--     Calculate density
9e5f2093d3 Alis*0304 #ifdef ALLOW_AUTODIFF_TAMC
7448700841 Mart*0305 CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, kind = isbyte
                0306 CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, kind = isbyte
9e5f2093d3 Alis*0307 #endif /* ALLOW_AUTODIFF_TAMC */
842349fcb7 Jean*0308           CALL FIND_RHO_2D(
                0309      I              iMin, iMax, jMin, jMax, k,
b98355c8c5 jm-c 0310      I              theta(1-OLx,1-OLy,k,bi,bj),
                0311      I              salt(1-OLx,1-OLy,k,bi,bj),
842349fcb7 Jean*0312      O              alphaRho,
                0313      I              k, bi, bj, myThid )
                0314         ELSE
                0315           DO j=jMin,jMax
                0316            DO i=iMin,iMax
                0317              alphaRho(i,j) = rhoInSitu(i,j,k,bi,bj)
                0318            ENDDO
                0319           ENDDO
                0320         ENDIF
c5c831ccb0 Jean*0321
7448700841 Mart*0322 #ifdef ALLOW_AUTODIFF_TAMC
                0323 CADJ STORE alphaRho = comlev1_bibj_k, key = kkey, kind = isbyte
                0324 #endif
8fabf86828 Jean*0325 C--     Calculate specific volume anomaly "alpha_prime":
                0326 C       = 1/(rho_prime + rho_Cst) - 1/rho_Cst
                0327 C       = - (rho_prime/rho_Cst) / ( 1 + rho_prime/rho_Cst ) / rho_Cst
1adfdf6281 Jean*0328         DO j=jMin,jMax
                0329           DO i=iMin,iMax
8fabf86828 Jean*0330             locBuoy = alphaRho(i,j)*recip_rhoConst
                0331             alphaRho(i,j) = -maskC(i,j,k,bi,bj)*recip_rhoConst
                0332      &                    * locBuoy/( oneRL + locBuoy )
1adfdf6281 Jean*0333           ENDDO
                0334         ENDDO
                0335
0d373cd03b Jean*0336 #ifdef ALLOW_MOM_COMMON
                0337 C--     Quasi-hydrostatic terms are added as if they modify the specific-volume
ee19c4a296 Jean*0338         IF ( quasiHydrostatic ) THEN
                0339          CALL MOM_QUASIHYDROSTATIC( bi, bj, k, uVel, vVel,
                0340      U                              alphaRho,
                0341      I                              myTime, myIter, myThid )
0d373cd03b Jean*0342         ENDIF
                0343 #endif /* ALLOW_MOM_COMMON */
7448700841 Mart*0344 #ifdef ALLOW_AUTODIFF_TAMC
                0345 CADJ STORE alphaRho = comlev1_bibj_k, key = kkey, kind = isbyte
                0346 #endif
0d373cd03b Jean*0347
7e1bd93d4f Jean*0348 C----  Hydrostatic pressure at cell centers
                0349
                0350        IF (integr_GeoPot.EQ.1) THEN
                0351 C  --  Finite Volume Form
                0352
                0353          DO j=jMin,jMax
9e5f2093d3 Alis*0354           DO i=iMin,iMax
fb481a83c2 Alis*0355
7e1bd93d4f Jean*0356 C---------- This discretization is the "finite volume" form
                0357 C           which has not been used to date since it does not
                0358 C           conserve KE+PE exactly even though it is more natural
2b8326b7dd Jean*0359
807d4a3bf5 Jean*0360            IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
f1f873e18c Jean*0361              ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
                0362 #ifdef NONLIN_FRSURF
                0363              ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
                0364 #endif
                0365              phiHydC(i,j) = ddRloc*alphaRho(i,j)
0f5ba23f97 Jean*0366 c--to reproduce results of c48d_post: uncomment those 4+1 lines
f1f873e18c Jean*0367 c            phiHydC(i,j)=phiHydF(i,j)
0d373cd03b Jean*0368 c    &          +(hFacC(i,j,k,bi,bj)-halfRL)*drF(k)*alphaRho(i,j)
f1f873e18c Jean*0369 c            phiHydF(i,j)=phiHydF(i,j)
                0370 c    &          + hFacC(i,j,k,bi,bj)*drF(k)*alphaRho(i,j)
                0371            ELSE
0d373cd03b Jean*0372              phiHydC(i,j) = phiHydF(i,j) + halfRL*drF(k)*alphaRho(i,j)
                0373 c            phiHydF(i,j) = phiHydF(i,j) +        drF(k)*alphaRho(i,j)
f1f873e18c Jean*0374            ENDIF
                0375 c-- and comment this last one:
0d373cd03b Jean*0376              phiHydF(i,j) = phiHydC(i,j) + halfRL*drF(k)*alphaRho(i,j)
f1f873e18c Jean*0377 c-----
7e1bd93d4f Jean*0378           ENDDO
                0379          ENDDO
                0380
                0381        ELSE
f1f873e18c Jean*0382 C  --  Finite Difference Form, with Part-Cell Bathy
                0383
0d373cd03b Jean*0384          dRlocM = halfRL*drC(k)
f1f873e18c Jean*0385          IF (k.EQ.1) dRlocM=rF(k)-rC(k)
                0386          IF (k.EQ.Nr) THEN
                0387            dRlocP=rC(k)-rF(k+1)
                0388          ELSE
0d373cd03b Jean*0389            dRlocP=halfRL*drC(k+1)
f1f873e18c Jean*0390          ENDIF
0d373cd03b Jean*0391          rec_dRm = oneRL/(rF(k)-rC(k))
                0392          rec_dRp = oneRL/(rC(k)-rF(k+1))
7e1bd93d4f Jean*0393
                0394          DO j=jMin,jMax
                0395           DO i=iMin,iMax
                0396
                0397 C---------- This discretization is the "energy conserving" form
                0398
807d4a3bf5 Jean*0399            IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
f1f873e18c Jean*0400              ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
                0401 #ifdef NONLIN_FRSURF
                0402              ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
                0403 #endif
0d373cd03b Jean*0404              phiHydC(i,j) =( MAX(zeroRL,ddRloc)*rec_dRm*dRlocM
                0405      &                      +MIN(zeroRL,ddRloc)*rec_dRp*dRlocP
f1f873e18c Jean*0406      &                     )*alphaRho(i,j)
                0407            ELSE
                0408              phiHydC(i,j) = phiHydF(i,j) + dRlocM*alphaRho(i,j)
                0409            ENDIF
                0410              phiHydF(i,j) = phiHydC(i,j) + dRlocP*alphaRho(i,j)
9e5f2093d3 Alis*0411           ENDDO
7e1bd93d4f Jean*0412          ENDDO
                0413
                0414 C  --  end if integr_GeoPot = ...
                0415        ENDIF
fb481a83c2 Alis*0416
f1f873e18c Jean*0417       ELSEIF ( buoyancyRelation .EQ. 'ATMOSPHERIC' ) THEN
45adaae56c Jean*0418 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
fb481a83c2 Alis*0419 C       This is the hydrostatic geopotential calculation for the Atmosphere
                0420 C       The ideal gas law is used implicitly here rather than calculating
                0421 C       the specific volume, analogous to the oceanic case.
                0422
5edeab65c1 Jean*0423         IF ( implicitIntGravWave .OR. myIter.LT.0 ) THEN
6cdaaf6acc Jean*0424 C--     virtual potential temperature anomaly (including water vapour effect)
5edeab65c1 Jean*0425           IF ( select_rStar.GE.1 .OR. selectSigmaCoord.GE.1 ) THEN
cfbe180681 Jean*0426 C-      isothermal (theta=const) reference state
5edeab65c1 Jean*0427             thetaRef = thetaConst
                0428           ELSE
cfbe180681 Jean*0429 C-      horizontally uniform (tRef) reference state
5edeab65c1 Jean*0430             thetaRef = tRef(k)
                0431           ENDIF
                0432           DO j=jMin,jMax
                0433            DO i=iMin,iMax
b98355c8c5 jm-c 0434             alphaRho(i,j) = ( theta(i,j,k,bi,bj)
                0435      &                        *( salt(i,j,k,bi,bj)*atm_Rq + oneRL )
5edeab65c1 Jean*0436      &                      - thetaRef )*maskC(i,j,k,bi,bj)
                0437            ENDDO
                0438           ENDDO
                0439         ELSE
                0440           DO j=jMin,jMax
                0441            DO i=iMin,iMax
                0442              alphaRho(i,j) = rhoInSitu(i,j,k,bi,bj)
                0443            ENDDO
                0444           ENDDO
cfbe180681 Jean*0445         ENDIF
6cdaaf6acc Jean*0446
0d373cd03b Jean*0447 #ifdef ALLOW_MOM_COMMON
                0448 C--     Quasi-hydrostatic terms are added in as if they modify the Pot.Temp
ee19c4a296 Jean*0449         IF ( quasiHydrostatic ) THEN
                0450          CALL MOM_QUASIHYDROSTATIC( bi, bj, k, uVel, vVel,
                0451      U                              alphaRho,
                0452      I                              myTime, myIter, myThid )
0d373cd03b Jean*0453         ENDIF
                0454 #endif /* ALLOW_MOM_COMMON */
7448700841 Mart*0455 #ifdef ALLOW_AUTODIFF_TAMC
                0456 CADJ STORE alphaRho = comlev1_bibj_k, key = kkey, kind = isbyte
                0457 #endif
0d373cd03b Jean*0458
f1f873e18c Jean*0459 C---    Integrate d Phi / d pi
fb481a83c2 Alis*0460
3d9080e0a5 Jean*0461 #ifndef DISABLE_SIGMA_CODE
                0462 C  --  Finite Volume Form, integrated to r-level (cell center & upper interface)
                0463        IF ( useFVgradPhi ) THEN
                0464
                0465          fullDepth = rF(1)-rF(Nr+1)
                0466          DO j=jMin,jMax
                0467           DO i=iMin,iMax
                0468            locDepth = Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj)
                0469 #ifdef NONLIN_FRSURF
                0470            locDepth = locDepth + surfPhiFac*etaH(i,j,bi,bj)
                0471 #endif
                0472            pKappaF(i,j) = (
                0473      &           ( R_low(i,j,bi,bj) + aHybSigmF( k )*fullDepth
                0474      &                              + bHybSigmF( k )*locDepth
                0475      &           )/atm_Po )**atm_kappa
                0476            pKappaC      = (
                0477      &           ( R_low(i,j,bi,bj) + aHybSigmC( k )*fullDepth
                0478      &                              + bHybSigmC( k )*locDepth
                0479      &           )/atm_Po )**atm_kappa
                0480            pKappaU(i,j) = (
                0481      &           ( R_low(i,j,bi,bj)+ aHybSigmF(k+1)*fullDepth
                0482      &                             + bHybSigmF(k+1)*locDepth
                0483      &           )/atm_Po )**atm_kappa
                0484            phiHydC(i,j) = phiHydF(i,j)
                0485      &        + atm_Cp*( pKappaF(i,j) - pKappaC      )*alphaRho(i,j)
                0486            phiHydU(i,j) = phiHydF(i,j)
                0487      &        + atm_Cp*( pKappaF(i,j) - pKappaU(i,j) )*alphaRho(i,j)
                0488           ENDDO
                0489          ENDDO
                0490 C end: Finite Volume Form, integrated to r-level.
                0491
                0492        ELSEIF (integr_GeoPot.EQ.0) THEN
                0493 #else /* DISABLE_SIGMA_CODE */
f1f873e18c Jean*0494        IF (integr_GeoPot.EQ.0) THEN
3d9080e0a5 Jean*0495 #endif /* DISABLE_SIGMA_CODE */
f1f873e18c Jean*0496 C  --  Energy Conserving Form, accurate with Full cell topo  --
45adaae56c Jean*0497 C------------ The integration for the first level phi(k=1) is the same
                0498 C             for both the "finite volume" and energy conserving methods.
0f5ba23f97 Jean*0499 C    *NOTE* o Working with geopotential Anomaly, the geopotential boundary
e9b27c9813 Alis*0500 C             condition is simply Phi-prime(Ro_surf)=0.
45adaae56c Jean*0501 C           o convention ddPI > 0 (same as drF & drC)
                0502 C-----------------------------------------------------------------------
f1f873e18c Jean*0503          IF (k.EQ.1) THEN
                0504            ddPIm=atm_Cp*( ((rF( k )/atm_Po)**atm_kappa)
                0505      &                   -((rC( k )/atm_Po)**atm_kappa) )
                0506          ELSE
                0507            ddPIm=atm_Cp*( ((rC(k-1)/atm_Po)**atm_kappa)
0d373cd03b Jean*0508      &                   -((rC( k )/atm_Po)**atm_kappa) )*halfRL
f1f873e18c Jean*0509          ENDIF
                0510          IF (k.EQ.Nr) THEN
                0511            ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
                0512      &                   -((rF(k+1)/atm_Po)**atm_kappa) )
                0513          ELSE
                0514            ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
0d373cd03b Jean*0515      &                   -((rC(k+1)/atm_Po)**atm_kappa) )*halfRL
f1f873e18c Jean*0516          ENDIF
45adaae56c Jean*0517 C-------- This discretization is the energy conserving form
f1f873e18c Jean*0518          DO j=jMin,jMax
                0519           DO i=iMin,iMax
6cdaaf6acc Jean*0520              phiHydC(i,j) = phiHydF(i,j) +ddPIm*alphaRho(i,j)
                0521              phiHydF(i,j) = phiHydC(i,j) +ddPIp*alphaRho(i,j)
45adaae56c Jean*0522           ENDDO
f1f873e18c Jean*0523          ENDDO
45adaae56c Jean*0524 C end: Energy Conserving Form, No hFac  --
fb481a83c2 Alis*0525 C-----------------------------------------------------------------------
                0526
f1f873e18c Jean*0527        ELSEIF (integr_GeoPot.EQ.1) THEN
                0528 C  --  Finite Volume Form, with Part-Cell Topo, linear in P by Half level
45adaae56c Jean*0529 C---------
                0530 C  Finite Volume formulation consistent with Partial Cell, linear in p by piece
                0531 C   Note: a true Finite Volume form should be linear between 2 Interf_W :
                0532 C     phi_C = (phi_W_k+ phi_W_k+1)/2 ; but not accurate in Stratosphere (low p)
                0533 C   also: if Interface_W at the middle between tracer levels, this form
0f5ba23f97 Jean*0534 C     is close to the Energy Cons. form in the Interior, except for the
45adaae56c Jean*0535 C     non-linearity in PI(p)
                0536 C---------
f1f873e18c Jean*0537            ddPIm=atm_Cp*( ((rF( k )/atm_Po)**atm_kappa)
                0538      &                   -((rC( k )/atm_Po)**atm_kappa) )
                0539            ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
                0540      &                   -((rF(k+1)/atm_Po)**atm_kappa) )
                0541          DO j=jMin,jMax
                0542           DO i=iMin,iMax
807d4a3bf5 Jean*0543            IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
f1f873e18c Jean*0544              ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
                0545 #ifdef NONLIN_FRSURF
                0546              ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
                0547 #endif
                0548              phiHydC(i,j) = ddRloc*recip_drF(k)*2. _d 0
6cdaaf6acc Jean*0549      &          *ddPIm*alphaRho(i,j)
f1f873e18c Jean*0550            ELSE
6cdaaf6acc Jean*0551              phiHydC(i,j) = phiHydF(i,j) +ddPIm*alphaRho(i,j)
f1f873e18c Jean*0552            ENDIF
6cdaaf6acc Jean*0553              phiHydF(i,j) = phiHydC(i,j) +ddPIp*alphaRho(i,j)
45adaae56c Jean*0554           ENDDO
f1f873e18c Jean*0555          ENDDO
                0556 C end: Finite Volume Form, with Part-Cell Topo, linear in P by Half level
fb481a83c2 Alis*0557 C-----------------------------------------------------------------------
                0558
f1f873e18c Jean*0559        ELSEIF ( integr_GeoPot.EQ.2
                0560      &     .OR. integr_GeoPot.EQ.3 ) THEN
0f5ba23f97 Jean*0561 C  --  Finite Difference Form, with Part-Cell Topo,
f1f873e18c Jean*0562 C       works with Interface_W  at the middle between 2.Tracer_Level
                0563 C        and  with Tracer_Level at the middle between 2.Interface_W.
45adaae56c Jean*0564 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0565 C  Finite Difference formulation consistent with Partial Cell,
                0566 C   Valid & accurate if Interface_W at middle between tracer levels
0f5ba23f97 Jean*0567 C   linear in p between 2 Tracer levels ; conserve energy in the Interior
45adaae56c Jean*0568 C---------
f1f873e18c Jean*0569          IF (k.EQ.1) THEN
                0570            ddPIm=atm_Cp*( ((rF( k )/atm_Po)**atm_kappa)
                0571      &                   -((rC( k )/atm_Po)**atm_kappa) )
                0572          ELSE
                0573            ddPIm=atm_Cp*( ((rC(k-1)/atm_Po)**atm_kappa)
0d373cd03b Jean*0574      &                   -((rC( k )/atm_Po)**atm_kappa) )*halfRL
f1f873e18c Jean*0575          ENDIF
                0576          IF (k.EQ.Nr) THEN
                0577            ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
                0578      &                   -((rF(k+1)/atm_Po)**atm_kappa) )
                0579          ELSE
                0580            ddPIp=atm_Cp*( ((rC( k )/atm_Po)**atm_kappa)
0d373cd03b Jean*0581      &                   -((rC(k+1)/atm_Po)**atm_kappa) )*halfRL
f1f873e18c Jean*0582          ENDIF
0d373cd03b Jean*0583          rec_dRm = oneRL/(rF(k)-rC(k))
                0584          rec_dRp = oneRL/(rC(k)-rF(k+1))
f1f873e18c Jean*0585          DO j=jMin,jMax
                0586           DO i=iMin,iMax
807d4a3bf5 Jean*0587            IF (k.EQ.kSurfC(i,j,bi,bj)) THEN
f1f873e18c Jean*0588              ddRloc = Ro_surf(i,j,bi,bj)-rC(k)
                0589 #ifdef NONLIN_FRSURF
                0590              ddRloc = ddRloc + surfPhiFac*etaH(i,j,bi,bj)
                0591 #endif
0d373cd03b Jean*0592              phiHydC(i,j) =( MAX(zeroRL,ddRloc)*rec_dRm*ddPIm
                0593      &                      +MIN(zeroRL,ddRloc)*rec_dRp*ddPIp
                0594      &                     )*alphaRho(i,j)
f1f873e18c Jean*0595            ELSE
6cdaaf6acc Jean*0596              phiHydC(i,j) = phiHydF(i,j) +ddPIm*alphaRho(i,j)
f1f873e18c Jean*0597            ENDIF
6cdaaf6acc Jean*0598              phiHydF(i,j) = phiHydC(i,j) +ddPIp*alphaRho(i,j)
fb481a83c2 Alis*0599           ENDDO
f1f873e18c Jean*0600          ENDDO
                0601 C end: Finite Difference Form, with Part-Cell Topo
45adaae56c Jean*0602 C-----------------------------------------------------------------------
fb481a83c2 Alis*0603
f1f873e18c Jean*0604        ELSE
                0605          STOP 'CALC_PHI_HYD: Bad integr_GeoPot option !'
                0606        ENDIF
fb481a83c2 Alis*0607
45adaae56c Jean*0608 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
fb481a83c2 Alis*0609       ELSE
463053c692 Jean*0610         STOP 'CALC_PHI_HYD: Bad value of buoyancyRelation !'
fb481a83c2 Alis*0611       ENDIF
924557e60a Chri*0612
7448700841 Mart*0613 #ifdef ALLOW_AUTODIFF_TAMC
                0614 # ifdef NONLIN_FRSURF
                0615 CADJ STORE alphaRho = comlev1_bibj_k, key = kkey, kind = isbyte
                0616 CADJ STORE phiHydC  = comlev1_bibj_k, key = kkey, kind = isbyte
                0617 CADJ STORE phiHydF  = comlev1_bibj_k, key = kkey, kind = isbyte
                0618 # endif
                0619 #endif
                0620
3d9080e0a5 Jean*0621       IF ( .NOT. useFVgradPhi ) THEN
                0622 C--   r-coordinate and r*-coordinate cases:
                0623
                0624        IF ( momPressureForcing ) THEN
                0625         CALL CALC_GRAD_PHI_HYD(
                0626      I                         k, bi, bj, iMin,iMax, jMin,jMax,
b98355c8c5 jm-c 0627      I                         phiHydC, alphaRho,
3d9080e0a5 Jean*0628      O                         dPhiHydX, dPhiHydY,
                0629      I                         myTime, myIter, myThid)
                0630        ENDIF
                0631
                0632 #ifndef DISABLE_SIGMA_CODE
                0633       ELSE
                0634 C--   else (SigmaCoords part)
                0635
                0636        IF ( fluidIsWater ) THEN
                0637         STOP 'CALC_PHI_HYD: missing code for SigmaCoord'
                0638        ENDIF
                0639        IF ( momPressureForcing ) THEN
                0640         CALL CALC_GRAD_PHI_FV(
                0641      I                         k, bi, bj, iMin,iMax, jMin,jMax,
                0642      I                         phiHydF, phiHydU, pKappaF, pKappaU,
                0643      O                         dPhiHydX, dPhiHydY,
                0644      I                         myTime, myIter, myThid)
                0645        ENDIF
                0646        DO j=jMin,jMax
                0647          DO i=iMin,iMax
                0648            phiHydF(i,j) = phiHydU(i,j)
                0649          ENDDO
                0650        ENDDO
                0651
                0652 #endif /* DISABLE_SIGMA_CODE */
                0653 C--   end if-not/else useFVgradPhi
                0654       ENDIF
                0655
f1f873e18c Jean*0656 C---   Diagnose Phi at boundary r=R_low :
                0657 C       = Ocean bottom pressure (Ocean, Z-coord.)
                0658 C       = Sea-surface height    (Ocean, P-coord.)
                0659 C       = Top atmosphere height (Atmos, P-coord.)
                0660       IF (useDiagPhiRlow) THEN
                0661         CALL DIAGS_PHI_RLOW(
                0662      I                      k, bi, bj, iMin,iMax, jMin,jMax,
b98355c8c5 jm-c 0663      I                      phiHydF, phiHydC, alphaRho,
0f5ba23f97 Jean*0664      I                      myTime, myIter, myThid)
f1f873e18c Jean*0665       ENDIF
                0666
                0667 C---   Diagnose Full Hydrostatic Potential at cell center level
                0668         CALL DIAGS_PHI_HYD(
                0669      I                      k, bi, bj, iMin,iMax, jMin,jMax,
                0670      I                      phiHydC,
                0671      I                      myTime, myIter, myThid)
                0672
45adaae56c Jean*0673 #endif /* INCLUDE_PHIHYD_CALCULATION_CODE */
1dbaea09ee Chri*0674
09f06f8597 Jean*0675       RETURN
                0676       END