Back to home page

MITgcm
 
 

    

File indexing completed on 2022-10-25 05:09:03 UTC

view on githubraw file Latest commit dbb54b7d on 2022-10-24 17:07:16 UTC
7e1bd93d4f Jean*0001 #include "CPP_OPTIONS.h"
64b9454396 Jean*0002 #undef OLD_PSTAR_SLOPE_TERM
7e1bd93d4f Jean*0003 
                0004 CBOP
                0005 C     !ROUTINE: CALC_GRAD_PHI_HYD
                0006 C     !INTERFACE:
4606c28752 Jean*0007       SUBROUTINE CALC_GRAD_PHI_HYD(
7e1bd93d4f Jean*0008      I                       k, bi, bj, iMin,iMax, jMin,jMax,
b98355c8c5 jm-c 0009      I                       phiHydC, alphRho,
7e1bd93d4f Jean*0010      O                       dPhiHydX, dPhiHydY,
                0011      I                       myTime, myIter, myThid)
                0012 C     !DESCRIPTION: \bv
                0013 C     *==========================================================*
4606c28752 Jean*0014 C     | S/R CALC_GRAD_PHI_HYD
                0015 C     | o Calculate the gradient of Hydrostatic potential anomaly
7e1bd93d4f Jean*0016 C     *==========================================================*
                0017 C     \ev
                0018 
                0019 C     !USES:
                0020       IMPLICIT NONE
                0021 C     == Global variables ==
                0022 #include "SIZE.h"
                0023 #include "EEPARAMS.h"
                0024 #include "PARAMS.h"
                0025 #include "GRID.h"
                0026 #include "SURFACE.h"
f1f873e18c Jean*0027 #include "DYNVARS.h"
7e1bd93d4f Jean*0028 
                0029 C     !INPUT/OUTPUT PARAMETERS:
                0030 C     == Routine Arguments ==
4606c28752 Jean*0031 C     bi,bj      :: tile index
7e1bd93d4f Jean*0032 C     iMin,iMax,jMin,jMax :: Loop counters
4606c28752 Jean*0033 C     phiHydC    :: Hydrostatic Potential anomaly
7e1bd93d4f Jean*0034 C                  (atmos: =Geopotential ; ocean-z: =Pressure/rho)
                0035 C     alphRho    :: Density (z-coord) or specific volume (p-coord)
                0036 C     dPhiHydX,Y :: Gradient (X & Y directions) of Hyd. Potential
                0037 C     myTime :: Current time
                0038 C     myIter :: Current iteration number
                0039 C     myThid :: Instance number for this call of the routine.
                0040       INTEGER k, bi,bj, iMin,iMax, jMin,jMax
f1f873e18c Jean*0041       _RL phiHydC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
7e1bd93d4f Jean*0042       _RL alphRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
08b3a88a6b Jean*0043       _RL dPhiHydX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
                0044       _RL dPhiHydY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
7e1bd93d4f Jean*0045       _RL myTime
                0046       INTEGER myIter, myThid
                0047 
                0048 #ifdef INCLUDE_PHIHYD_CALCULATION_CODE
                0049 
                0050 C     !LOCAL VARIABLES:
                0051 C     == Local variables ==
                0052 C     i,j :: Loop counters
                0053       INTEGER i,j
08b3a88a6b Jean*0054       _RL varLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
8996cf5a3c Jean*0055 #ifdef NONLIN_FRSURF
dbb54b7d7a Jean*0056       LOGICAL generalForm
                0057       _RL factorP, factPI
6cdaaf6acc Jean*0058       CHARACTER*(MAX_LEN_MBUF) msgBuf
8996cf5a3c Jean*0059 #endif
7e1bd93d4f Jean*0060 CEOP
                0061 
f1f873e18c Jean*0062 #ifdef NONLIN_FRSURF
6cdaaf6acc Jean*0063       IF (select_rStar.GE.2 .AND. nonlinFreeSurf.GE.4 ) THEN
cdc9f269ae Patr*0064 # ifndef DISABLE_RSTAR_CODE
6cdaaf6acc Jean*0065 C-    Integral of b.dr = rStarFac * Integral of b.dr* :
4606c28752 Jean*0066 C      and will add later (select_rStar=2) the contribution of
6cdaaf6acc Jean*0067 C      the slope of the r* coordinate.
45e2ae5569 Jean*0068        IF ( fluidIsAir ) THEN
                0069 C-     Consistent with Phi'= Integr[ theta'.dPI ] :
1657c55ee0 Jean*0070         DO j=jMin,jMax
                0071          DO i=iMin,iMax
45e2ae5569 Jean*0072           varLoc(i,j) = phiHydC(i,j)*pStarFacK(i,j,bi,bj)
6cdaaf6acc Jean*0073      &                + phi0surf(i,j,bi,bj)
                0074          ENDDO
f1f873e18c Jean*0075         ENDDO
6cdaaf6acc Jean*0076        ELSE
1657c55ee0 Jean*0077         DO j=jMin,jMax
                0078          DO i=iMin,iMax
6cdaaf6acc Jean*0079           varLoc(i,j) = phiHydC(i,j)*rStarFacC(i,j,bi,bj)
                0080      &                + phi0surf(i,j,bi,bj)
                0081          ENDDO
                0082         ENDDO
                0083        ENDIF
                0084       ELSEIF (select_rStar.GE.1 .AND. nonlinFreeSurf.GE.4 ) THEN
                0085 C-    Integral of b.dr but scaled to correspond to a fixed r-level (=r*)
                0086 C      no contribution of the slope of the r* coordinate (select_rStar=1)
45e2ae5569 Jean*0087        IF ( fluidIsAir ) THEN
                0088 C-     Consistent with Phi'= Integr[ theta'.dPI ] :
1657c55ee0 Jean*0089         DO j=jMin,jMax
                0090          DO i=iMin,iMax
6cdaaf6acc Jean*0091           IF (Ro_surf(i,j,bi,bj).EQ.rC(k)) THEN
4606c28752 Jean*0092            factPI=atm_Cp*( ((etaH(i,j,bi,bj)+rC(k))/atm_Po)**atm_kappa
                0093      &                    -(                 rC(k) /atm_Po)**atm_kappa
6cdaaf6acc Jean*0094      &                  )
                0095            varLoc(i,j) = factPI*alphRho(i,j)
0bf1d01915 Jean*0096      &                 + phi0surf(i,j,bi,bj)
6cdaaf6acc Jean*0097           ELSEIF (Ro_surf(i,j,bi,bj).NE.0. _d 0) THEN
                0098            factPI = (rC(k)/Ro_surf(i,j,bi,bj))**atm_kappa
                0099            varLoc(i,j) = phiHydC(i,j)
64b9454396 Jean*0100      &                  *(pStarFacK(i,j,bi,bj) - factPI)
6cdaaf6acc Jean*0101      &                  /(1. _d 0 -factPI)
                0102      &                 + phi0surf(i,j,bi,bj)
                0103           ENDIF
                0104          ENDDO
                0105         ENDDO
                0106        ELSE
1657c55ee0 Jean*0107         DO j=jMin,jMax
                0108          DO i=iMin,iMax
6cdaaf6acc Jean*0109           IF (Ro_surf(i,j,bi,bj).EQ.rC(k)) THEN
a7400ee2ba Jean*0110            WRITE(msgBuf,'(3A)') 'CALC_GRAD_PHI_HYD: ',
                0111      &      'Problem when Ro_surf=rC',
806349206d Jean*0112      &      ' with select_rStar,nonlinFreeSurf=1,4'
6cdaaf6acc Jean*0113            CALL PRINT_ERROR( msgBuf , myThid)
                0114            STOP 'CALC_GRAD_PHI_HYD: Pb in r* options implementation'
                0115           ELSE
                0116            varLoc(i,j) = phiHydC(i,j)
                0117      &                  *(etaH(i,j,bi,bj)+Ro_surf(i,j,bi,bj)-rC(k))
                0118      &                  /                (Ro_surf(i,j,bi,bj)-rC(k))
                0119      &                 + phi0surf(i,j,bi,bj)
                0120           ENDIF
                0121          ENDDO
                0122         ENDDO
                0123        ENDIF
cdc9f269ae Patr*0124 # endif /* DISABLE_RSTAR_CODE */
f1f873e18c Jean*0125       ELSE
                0126 #else /* NONLIN_FRSURF */
                0127       IF (.TRUE.) THEN
                0128 #endif /* NONLIN_FRSURF */
1657c55ee0 Jean*0129        DO j=jMin,jMax
                0130         DO i=iMin,iMax
f1f873e18c Jean*0131          varLoc(i,j) = phiHydC(i,j)+phi0surf(i,j,bi,bj)
                0132         ENDDO
                0133        ENDDO
                0134       ENDIF
7e1bd93d4f Jean*0135 
f1f873e18c Jean*0136 C--   Zonal & Meridional gradient of potential anomaly
d7fe92a587 Jean*0137       DO j=1-OLy,sNy+OLy
                0138        DO i=1-OLx,sNx+OLx
                0139         dPhiHydX(i,j)  = 0. _d 0
                0140         dPhiHydY(i,j)  = 0. _d 0
                0141        ENDDO
                0142       ENDDO
7e1bd93d4f Jean*0143       DO j=jMin,jMax
1657c55ee0 Jean*0144        DO i=iMin+1,iMax
4606c28752 Jean*0145         dPhiHydX(i,j) = _recip_dxC(i,j,bi,bj)*recip_deepFacC(k)
                0146      &                *( varLoc(i,j)-varLoc(i-1,j) )*recip_rhoFacC(k)
1657c55ee0 Jean*0147        ENDDO
                0148       ENDDO
                0149       DO j=jMin+1,jMax
                0150        DO i=iMin,iMax
4606c28752 Jean*0151         dPhiHydY(i,j) = _recip_dyC(i,j,bi,bj)*recip_deepFacC(k)
                0152      &                *( varLoc(i,j)-varLoc(i,j-1) )*recip_rhoFacC(k)
7e1bd93d4f Jean*0153        ENDDO
                0154       ENDDO
                0155 
f1f873e18c Jean*0156 #ifdef NONLIN_FRSURF
45e2ae5569 Jean*0157 # ifndef DISABLE_RSTAR_CODE
f1f873e18c Jean*0158       IF (select_rStar.GE.2 .AND. nonlinFreeSurf.GE.1 ) THEN
dbb54b7d7a Jean*0159 C-     need to use general form anytime r* @ Top is not uniformly zero:
                0160        generalForm = useShelfIce .OR.
                0161      &  ( usingPCoords .AND. ( rF(Nr+1).NE.zeroRS ) ) .OR.
                0162      &  ( usingZCoords .AND. ( topoFile.NE.' ' .OR. rF(1).NE.zeroRS ) )
                0163 c      generalForm = .TRUE.
                0164        IF ( fluidIsWater .AND. ( usingZCoords .OR. generalForm ) ) THEN
                0165         IF ( usingZCoords ) THEN
ff02675122 Jean*0166 C--    z* coordinate slope term: rho_prime/rho0 * Grad_r(g.z)
dbb54b7d7a Jean*0167          factorP = gravity*recip_rhoConst*recip_rhoFacC(k)*0.5 _d 0
                0168         ELSE
                0169 C--    p* coordinate slope term: alpha_prime * Grad_r( p )
                0170          factorP = 0.5 _d 0
                0171         ENDIF
                0172         IF ( generalForm ) THEN
                0173 C-     general case, works for both P & Z coordinates:
                0174          DO j=jMin,jMax
                0175           DO i=iMin,iMax
                0176            varLoc(i,j) = etaH(i,j,bi,bj)*recip_Rcol(i,j,bi,bj)
                0177      &                 *( rC(k) - R_low(i,j,bi,bj) )
                0178           ENDDO
f1f873e18c Jean*0179          ENDDO
dbb54b7d7a Jean*0180         ELSE
                0181 C-     Z-coordinate with flat top at z = 0:
                0182          DO j=jMin,jMax
                0183           DO i=iMin,iMax
                0184            varLoc(i,j) = etaH(i,j,bi,bj)
                0185      &                 *(1. _d 0 + rC(k)*recip_Rcol(i,j,bi,bj))
                0186           ENDDO
                0187          ENDDO
                0188         ENDIF
f1f873e18c Jean*0189         DO j=jMin,jMax
1657c55ee0 Jean*0190          DO i=iMin+1,iMax
f1f873e18c Jean*0191           dPhiHydX(i,j) = dPhiHydX(i,j)
dbb54b7d7a Jean*0192      &     +factorP*(alphRho(i-1,j)+alphRho(i,j))
f1f873e18c Jean*0193      &             *(varLoc(i,j)-varLoc(i-1,j))
4606c28752 Jean*0194      &             *recip_dxC(i,j,bi,bj)*recip_deepFacC(k)
1657c55ee0 Jean*0195          ENDDO
                0196         ENDDO
                0197         DO j=jMin+1,jMax
                0198          DO i=iMin,iMax
f1f873e18c Jean*0199           dPhiHydY(i,j) = dPhiHydY(i,j)
dbb54b7d7a Jean*0200      &     +factorP*(alphRho(i,j-1)+alphRho(i,j))
f1f873e18c Jean*0201      &             *(varLoc(i,j)-varLoc(i,j-1))
4606c28752 Jean*0202      &             *recip_dyC(i,j,bi,bj)*recip_deepFacC(k)
f1f873e18c Jean*0203          ENDDO
                0204         ENDDO
45e2ae5569 Jean*0205        ELSEIF ( fluidIsWater ) THEN
ff02675122 Jean*0206 C--    p* coordinate slope term: alpha_prime * Grad_r( p )
dbb54b7d7a Jean*0207 C-     requires top to be at p = 0:
f1f873e18c Jean*0208         factorP = 0.5 _d 0
                0209         DO j=jMin,jMax
1657c55ee0 Jean*0210          DO i=iMin+1,iMax
f1f873e18c Jean*0211           dPhiHydX(i,j) = dPhiHydX(i,j)
                0212      &     +factorP*(alphRho(i-1,j)+alphRho(i,j))
                0213      &             *(rStarFacC(i,j,bi,bj)-rStarFacC(i-1,j,bi,bj))
4606c28752 Jean*0214      &             *rC(k)*recip_dxC(i,j,bi,bj)*recip_deepFacC(k)
1657c55ee0 Jean*0215          ENDDO
                0216         ENDDO
                0217         DO j=jMin+1,jMax
                0218          DO i=iMin,iMax
f1f873e18c Jean*0219           dPhiHydY(i,j) = dPhiHydY(i,j)
                0220      &     +factorP*(alphRho(i,j-1)+alphRho(i,j))
                0221      &             *(rStarFacC(i,j,bi,bj)-rStarFacC(i,j-1,bi,bj))
4606c28752 Jean*0222      &             *rC(k)*recip_dyC(i,j,bi,bj)*recip_deepFacC(k)
f1f873e18c Jean*0223          ENDDO
                0224         ENDDO
45e2ae5569 Jean*0225        ELSEIF ( fluidIsAir ) THEN
                0226 #ifdef OLD_PSTAR_SLOPE_TERM
                0227 C--    p* coordinate slope term: alpha_prime * Grad_r( p ):
                0228 C      PI_star * (Theta_eq^prime)_bar_i * kappa * delta^i( rStarFacC )
                0229 C- Note: factor: ( p_s / p_s^o )^(kappa - 1) = rStarFacC^(kappa -1)
                0230 C        is missing here.
                0231         factorP = (rC(k)/atm_Po)**atm_kappa
                0232         factorP = (atm_Rd/rC(k))*factorP*0.5 _d 0
                0233 #else
                0234 C--    p* coordinate slope term: theta_prime * Grad_r( PI ):
                0235 C      PI_star * (Theta_eq^prime)_bar_i * delta^i( rStarFacC^kappa )
                0236 C      This is also consitent with geopotential factor: rStarFacC^kappa
                0237         factorP = halfRL*atm_Cp*(rC(k)/atm_Po)**atm_kappa
                0238 #endif
f1f873e18c Jean*0239         DO j=jMin,jMax
1657c55ee0 Jean*0240          DO i=iMin+1,iMax
f1f873e18c Jean*0241           dPhiHydX(i,j) = dPhiHydX(i,j)
6cdaaf6acc Jean*0242      &     +factorP*(alphRho(i-1,j)+alphRho(i,j))
45e2ae5569 Jean*0243 #ifdef OLD_PSTAR_SLOPE_TERM
f1f873e18c Jean*0244      &             *(rStarFacC(i,j,bi,bj)-rStarFacC(i-1,j,bi,bj))
4606c28752 Jean*0245      &             *rC(k)*recip_dxC(i,j,bi,bj)*recip_deepFacC(k)
45e2ae5569 Jean*0246 #else
                0247      &             *(pStarFacK(i,j,bi,bj)-pStarFacK(i-1,j,bi,bj))
                0248      &             *recip_dxC(i,j,bi,bj)*recip_deepFacC(k)
                0249 #endif
1657c55ee0 Jean*0250          ENDDO
                0251         ENDDO
                0252         DO j=jMin+1,jMax
                0253          DO i=iMin,iMax
f1f873e18c Jean*0254           dPhiHydY(i,j) = dPhiHydY(i,j)
6cdaaf6acc Jean*0255      &     +factorP*(alphRho(i,j-1)+alphRho(i,j))
45e2ae5569 Jean*0256 #ifdef OLD_PSTAR_SLOPE_TERM
f1f873e18c Jean*0257      &             *(rStarFacC(i,j,bi,bj)-rStarFacC(i,j-1,bi,bj))
4606c28752 Jean*0258      &             *rC(k)*recip_dyC(i,j,bi,bj)*recip_deepFacC(k)
45e2ae5569 Jean*0259 #else
                0260      &             *(pStarFacK(i,j,bi,bj)-pStarFacK(i,j-1,bi,bj))
                0261      &             *recip_dyC(i,j,bi,bj)*recip_deepFacC(k)
                0262 #endif
f1f873e18c Jean*0263          ENDDO
                0264         ENDDO
                0265        ENDIF
                0266       ENDIF
45e2ae5569 Jean*0267 # endif /* DISABLE_RSTAR_CODE */
f1f873e18c Jean*0268 #endif /* NONLIN_FRSURF */
                0269 
08b3a88a6b Jean*0270 C--   Apply mask:
                0271       DO j=1-OLy,sNy+OLy
                0272        DO i=1-OLx,sNx+OLx
                0273          dPhiHydX(i,j) = dPhiHydX(i,j)*_maskW(i,j,k,bi,bj)
                0274          dPhiHydY(i,j) = dPhiHydY(i,j)*_maskS(i,j,k,bi,bj)
                0275        ENDDO
                0276       ENDDO
                0277 
7e1bd93d4f Jean*0278 #endif /* INCLUDE_PHIHYD_CALCULATION_CODE */
                0279 
                0280       RETURN
                0281       END
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated from https://github.com/MITgcm/MITgcm by the 2.2.1-MITgcm-0.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!