Back to home page

MITgcm

 
 

    


File indexing completed on 2019-03-20 05:09:58 UTC

view on githubraw file Latest commit b98355c8 on 2019-03-19 19:03:21 UTC
1adfdf6281 Jean*0001 #include "CPP_OPTIONS.h"
                0002 
                0003 CBOP
                0004 C     !ROUTINE: DIAGS_PHI_RLOW
                0005 C     !INTERFACE:
e8c0c8afb8 Jean*0006       SUBROUTINE DIAGS_PHI_RLOW(
1adfdf6281 Jean*0007      I                       k, bi, bj, iMin,iMax, jMin,jMax,
b98355c8c5 jm-c 0008      I                       phiHydF, phiHydC, alphRho,
1adfdf6281 Jean*0009      I                       myTime, myIter, myThid)
                0010 C     !DESCRIPTION: \bv
                0011 C     *==========================================================*
e8c0c8afb8 Jean*0012 C     | S/R DIAGS_PHI_RLOW
1adfdf6281 Jean*0013 C     | o Diagnose Phi-Hydrostatic at r-lower boundary
e8c0c8afb8 Jean*0014 C     |   = bottom pressure (ocean in z-coord) ;
1adfdf6281 Jean*0015 C     |   = sea surface elevation (ocean in p-coord) ;
                0016 C     |   = height at the top of atmosphere (in p-coord) ;
                0017 C     *==========================================================*
                0018 C     \ev
                0019 
                0020 C     !USES:
                0021       IMPLICIT NONE
                0022 C     == Global variables ==
                0023 #include "SIZE.h"
                0024 #include "EEPARAMS.h"
                0025 #include "PARAMS.h"
                0026 #include "GRID.h"
                0027 #include "SURFACE.h"
                0028 #include "DYNVARS.h"
                0029 
                0030 C     !INPUT/OUTPUT PARAMETERS:
                0031 C     == Routine Arguments ==
e8c0c8afb8 Jean*0032 C     bi,bj      :: tile index
1adfdf6281 Jean*0033 C     iMin,iMax,jMin,jMax :: Loop counters
e8c0c8afb8 Jean*0034 C     phiHydF    :: hydrostatic potential anomaly at middle between
f1f873e18c Jean*0035 C                   2 centers k & k+1 (interface k+1)
                0036 C     phiHydC    :: hydrostatic potential anomaly at cell center
1adfdf6281 Jean*0037 C                  (atmos: =Geopotential ; ocean-z: =Pressure/rho)
                0038 C     alphRho    :: Density (z-coord) or specific volume (p-coord)
e8c0c8afb8 Jean*0039 C     myTime     :: Current time
                0040 C     myIter     :: Current iteration number
                0041 C     myThid     :: my Thread Id number
1adfdf6281 Jean*0042       INTEGER k, bi,bj, iMin,iMax, jMin,jMax
f1f873e18c Jean*0043       _RL phiHydF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0044       _RL phiHydC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
1adfdf6281 Jean*0045       _RL alphRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0046       _RL myTime
                0047       INTEGER myIter, myThid
                0048 
                0049 #ifdef INCLUDE_PHIHYD_CALCULATION_CODE
                0050 
                0051 C     !LOCAL VARIABLES:
                0052 C     == Local variables ==
                0053 C     i,j :: Loop counters
                0054       INTEGER i,j
f1f873e18c Jean*0055       _RL ddRloc, ratioRm, ratioRp
c75189e180 Jean*0056 #ifdef NONLIN_FRSURF
                0057       _RL facP, dPhiRef
                0058 #endif /* NONLIN_FRSURF */
1adfdf6281 Jean*0059 CEOP
                0060 
fb2b773b61 Jean*0061       IF ( usingZCoords ) THEN
1adfdf6281 Jean*0062 
f1f873e18c Jean*0063 C----- Compute bottom pressure deviation from gravity*rho0*H
ff02675122 Jean*0064 C      Start from phiHyd at the (bottom) tracer point and add Del_h*g*rho_prime
f1f873e18c Jean*0065 C      with Del_h = distance from the bottom up to tracer point
1adfdf6281 Jean*0066 
fb2b773b61 Jean*0067 C--    Initialise to zero (otherwise phi0surf accumulate over land)
                0068        IF ( k.EQ.1 ) THEN
a41ace1808 Jean*0069          DO j=1-OLy,sNy+OLy
                0070           DO i=1-OLx,sNx+OLx
fb2b773b61 Jean*0071             phiHydLow(i,j,bi,bj) = 0. _d 0
                0072           ENDDO
                0073          ENDDO
                0074        ENDIF
                0075 
1adfdf6281 Jean*0076        IF (integr_GeoPot.EQ.1) THEN
                0077 C  --  Finite Volume Form
                0078 
                0079          DO j=jMin,jMax
                0080           DO i=iMin,iMax
                0081            IF ( k .EQ. kLowC(i,j,bi,bj) ) THEN
f1f873e18c Jean*0082              ddRloc = rC(k)-R_low(i,j,bi,bj)
                0083              phiHydLow(i,j,bi,bj) = phiHydC(i,j)
bd053599dc Jean*0084      &        + ddRloc*gravFacC(k)*gravity*alphRho(i,j)*recip_rhoConst
1adfdf6281 Jean*0085            ENDIF
                0086           ENDDO
                0087          ENDDO
                0088 
                0089        ELSE
                0090 C  --  Finite Difference Form
                0091 
a41ace1808 Jean*0092          ratioRm = oneRL
                0093          ratioRp = oneRL
                0094          IF (k.GT.1 ) ratioRm = halfRL*drC(k)/(rF(k)-rC(k))
                0095          IF (k.LT.Nr) ratioRp = halfRL*drC(k+1)/(rC(k)-rF(k+1))
bd053599dc Jean*0096          ratioRm = ratioRm*gravFacF(k)
                0097          ratioRp = ratioRp*gravFacF(k+1)
1adfdf6281 Jean*0098 
                0099          DO j=jMin,jMax
                0100           DO i=iMin,iMax
f1f873e18c Jean*0101            IF ( k .EQ. kLowC(i,j,bi,bj) ) THEN
                0102              ddRloc = rC(k)-R_low(i,j,bi,bj)
                0103              phiHydLow(i,j,bi,bj) = phiHydC(i,j)
a41ace1808 Jean*0104      &                  +( MIN(zeroRL,ddRloc)*ratioRm
                0105      &                    +MAX(zeroRL,ddRloc)*ratioRp
f1f873e18c Jean*0106      &                   )*gravity*alphRho(i,j)*recip_rhoConst
1adfdf6281 Jean*0107            ENDIF
                0108           ENDDO
                0109          ENDDO
                0110 
                0111 C  --  end if integr_GeoPot = ...
                0112        ENDIF
e8c0c8afb8 Jean*0113 
fb2b773b61 Jean*0114 C  -- end usingZCoords
f1f873e18c Jean*0115       ENDIF
1adfdf6281 Jean*0116 
fb2b773b61 Jean*0117       IF ( k.EQ.Nr ) THEN
f1f873e18c Jean*0118 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
                0119 C  --  last level (bottom): rescale (r*) and add surface contribution
                0120 
fb2b773b61 Jean*0121        IF ( usingPCoords ) THEN
f1f873e18c Jean*0122 C  -- P coordinate : Phi(R_low) is simply at the top :
                0123         DO j=jMin,jMax
                0124          DO i=iMin,iMax
                0125            phiHydLow(i,j,bi,bj) = phiHydF(i,j)
1adfdf6281 Jean*0126          ENDDO
f1f873e18c Jean*0127         ENDDO
                0128        ENDIF
1adfdf6281 Jean*0129 
e8c0c8afb8 Jean*0130 #ifdef NONLIN_FRSURF
                0131 c      IF ( select_rStar.GE.2 .AND. nonlinFreeSurf.GE.4 ) THEN
4efe92625a Jean*0132        IF ( select_rStar.GE.1 .AND. nonlinFreeSurf.GE.4 ) THEN
                0133 C-    Integral of b.dr = rStarFac * Integral of b.dr* :
c75189e180 Jean*0134         IF ( fluidIsAir ) THEN
4efe92625a Jean*0135 C-     Consistent with Phi'= Integr[ theta'.dPi ] :
                0136          DO j=jMin,jMax
                0137           DO i=iMin,iMax
5332d6319a Jean*0138            facP = pStarFacK(i,j,bi,bj)
c75189e180 Jean*0139            dPhiRef = phiRef(2*k+1) - gravity*topoZ(i,j,bi,bj)
                0140      &                             - phi0surf(i,j,bi,bj)
                0141            phiHydLow(i,j,bi,bj) =
                0142      &              phiHydLow(i,j,bi,bj)*facP
a41ace1808 Jean*0143      &            + MAX( dPhiRef, zeroRL )*( facP - 1. _d 0 )
c75189e180 Jean*0144      &            + phi0surf(i,j,bi,bj)
                0145           ENDDO
                0146          ENDDO
                0147         ELSEIF ( usingPCoords ) THEN
                0148 C- Note: dPhiRef*(rStarFacC -1) = 1/rho*PSo*((eta+PSo)/PSo -1) = Bo_surf*etaN
                0149 C     so that this expression is the same as before (same as in "ELSE" block below)
                0150          DO j=jMin,jMax
                0151           DO i=iMin,iMax
                0152            dPhiRef = ( Ro_surf(i,j,bi,bj) - R_low(i,j,bi,bj) )
                0153      &              *recip_rhoConst
                0154            phiHydLow(i,j,bi,bj) =
                0155      &              phiHydLow(i,j,bi,bj)*rStarFacC(i,j,bi,bj)
                0156      &            + dPhiRef*( rStarFacC(i,j,bi,bj) - 1. _d 0 )
                0157      &            + phi0surf(i,j,bi,bj)
4efe92625a Jean*0158           ENDDO
e8c0c8afb8 Jean*0159          ENDDO
4efe92625a Jean*0160         ELSE
c75189e180 Jean*0161 C- Note: dPhiRef*(rStarFacC -1) = g*H*((eta+H)/H -1) = Bo_surf*etaN
                0162 C     so that this expression is the same as before (same as in "ELSE" block below)
4efe92625a Jean*0163          DO j=jMin,jMax
                0164           DO i=iMin,iMax
c75189e180 Jean*0165            dPhiRef = ( Ro_surf(i,j,bi,bj)-R_low(i,j,bi,bj) )
                0166      &              *gravity
                0167            phiHydLow(i,j,bi,bj) =
                0168      &              phiHydLow(i,j,bi,bj)*rStarFacC(i,j,bi,bj)
                0169      &            + dPhiRef*( rStarFacC(i,j,bi,bj) - 1. _d 0 )
                0170      &            + phi0surf(i,j,bi,bj)
4efe92625a Jean*0171           ENDDO
                0172          ENDDO
                0173         ENDIF
c75189e180 Jean*0174        ELSE
                0175 #else /* NONLIN_FRSURF */
                0176        IF ( .TRUE. ) THEN
e8c0c8afb8 Jean*0177 #endif /* NONLIN_FRSURF */
                0178 
c75189e180 Jean*0179          DO j=jMin,jMax
                0180           DO i=iMin,iMax
f1f873e18c Jean*0181            phiHydLow(i,j,bi,bj) = phiHydLow(i,j,bi,bj)
                0182      &            + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
e8c0c8afb8 Jean*0183      &            + phi0surf(i,j,bi,bj)
c75189e180 Jean*0184           ENDDO
                0185          ENDDO
                0186        ENDIF
1adfdf6281 Jean*0187 
                0188 C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
f1f873e18c Jean*0189 C  -- end if k=Nr.
1adfdf6281 Jean*0190       ENDIF
                0191 
                0192 #endif /* INCLUDE_PHIHYD_CALCULATION_CODE */
                0193 
                0194       RETURN
                0195       END