/[MITgcm]/MITgcm/model/src/diags_phi_rlow.F
ViewVC logotype

Diff of /MITgcm/model/src/diags_phi_rlow.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph | View Patch Patch

revision 1.1 by jmc, Sun Feb 9 02:58:39 2003 UTC revision 1.2 by jmc, Tue Feb 18 15:30:47 2003 UTC
# Line 8  C     !ROUTINE: DIAGS_PHI_RLOW Line 8  C     !ROUTINE: DIAGS_PHI_RLOW
8  C     !INTERFACE:  C     !INTERFACE:
9        SUBROUTINE DIAGS_PHI_RLOW(        SUBROUTINE DIAGS_PHI_RLOW(
10       I                       k, bi, bj, iMin,iMax, jMin,jMax,       I                       k, bi, bj, iMin,iMax, jMin,jMax,
11       I                       phiHyd, alphRho, tFld, sFld,       I                       phiHydF, phiHydC, alphRho, tFld, sFld,
12       I                       myTime, myIter, myThid)       I                       myTime, myIter, myThid)
13  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
14  C     *==========================================================*  C     *==========================================================*
# Line 34  C     !INPUT/OUTPUT PARAMETERS: Line 34  C     !INPUT/OUTPUT PARAMETERS:
34  C     == Routine Arguments ==  C     == Routine Arguments ==
35  C     bi,bj      :: tile index              C     bi,bj      :: tile index            
36  C     iMin,iMax,jMin,jMax :: Loop counters  C     iMin,iMax,jMin,jMax :: Loop counters
37  C     phiHyd     :: Hydrostatic Potential anomaly  C     phiHydF    :: hydrostatic potential anomaly at middle between
38    C                   2 centers k & k+1 (interface k+1)
39    C     phiHydC    :: hydrostatic potential anomaly at cell center
40  C                  (atmos: =Geopotential ; ocean-z: =Pressure/rho)  C                  (atmos: =Geopotential ; ocean-z: =Pressure/rho)
41  C     alphRho    :: Density (z-coord) or specific volume (p-coord)  C     alphRho    :: Density (z-coord) or specific volume (p-coord)
42  C     tFld       :: Potential temp.  C     tFld       :: Potential temp.
# Line 43  C     myTime :: Current time Line 45  C     myTime :: Current time
45  C     myIter :: Current iteration number  C     myIter :: Current iteration number
46  C     myThid :: Instance number for this call of the routine.  C     myThid :: Instance number for this call of the routine.
47        INTEGER k, bi,bj, iMin,iMax, jMin,jMax        INTEGER k, bi,bj, iMin,iMax, jMin,jMax
48        _RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL phiHydF(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
49          _RL phiHydC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
50        _RL alphRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL alphRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
51        _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
52        _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)        _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
# Line 57  C     == Local variables == Line 60  C     == Local variables ==
60  C     i,j :: Loop counters  C     i,j :: Loop counters
61        INTEGER i,j        INTEGER i,j
62        _RL zero, one, half        _RL zero, one, half
63        _RL dRloc        _RL ddRloc, ratioRm, ratioRp
64        PARAMETER ( zero= 0. _d 0 , one= 1. _d 0 , half= .5 _d 0 )        PARAMETER ( zero= 0. _d 0 , one= 1. _d 0 , half= .5 _d 0 )
65  CEOP  CEOP
66    
67          dRloc=drC(k)        IF ( buoyancyRelation .EQ. 'OCEANIC' ) THEN
         IF (k.EQ.1) dRloc=drF(1)  
68    
69        IF ( buoyancyRelation .eq. 'OCEANIC' ) THEN  C----- Compute bottom pressure deviation from gravity*rho0*H
70    C      Start from phiHyd at the (bottom) tracer point and add Del_h*g*rho'
71    C      with Del_h = distance from the bottom up to tracer point
72    
73         IF (integr_GeoPot.EQ.1) THEN         IF (integr_GeoPot.EQ.1) THEN
74  C  --  Finite Volume Form  C  --  Finite Volume Form
# Line 72  C  --  Finite Volume Form Line 76  C  --  Finite Volume Form
76           DO j=jMin,jMax           DO j=jMin,jMax
77            DO i=iMin,iMax            DO i=iMin,iMax
78             IF ( k .EQ. kLowC(i,j,bi,bj) ) THEN             IF ( k .EQ. kLowC(i,j,bi,bj) ) THEN
79              phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)               ddRloc = rC(k)-R_low(i,j,bi,bj)
80       &          + hFacC(i,j,k,bi,bj)               phiHydLow(i,j,bi,bj) = phiHydC(i,j)
81       &             *drF(K)*gravity*alphRho(i,j)*recip_rhoConst       &            + ddRloc*gravity*alphRho(i,j)*recip_rhoConst
      &          + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)  
      &          + phi0surf(i,j,bi,bj)  
82             ENDIF             ENDIF
83            ENDDO            ENDDO
84           ENDDO           ENDDO
# Line 84  C  --  Finite Volume Form Line 86  C  --  Finite Volume Form
86         ELSE         ELSE
87  C  --  Finite Difference Form  C  --  Finite Difference Form
88    
89  C---------- Compute bottom pressure deviation from gravity*rho0*H           ratioRm = one
90  C           This has to be done starting from phiHyd at the current           ratioRp = one
91  C           tracer point and .5 of the cell's thickness has to be           IF (k.GT.1 ) ratioRm = half*drC(k)/(rF(k)-rC(k))
92  C           substracted from hFacC           IF (k.LT.Nr) ratioRp = half*drC(k+1)/(rC(k)-rF(k+1))
93    
94           DO j=jMin,jMax           DO j=jMin,jMax
95            DO i=iMin,iMax            DO i=iMin,iMax
96             IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN             IF ( k .EQ. kLowC(i,j,bi,bj) ) THEN
97              phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)               ddRloc = rC(k)-R_low(i,j,bi,bj)
98       &          + (half*dRloc+(hFacC(i,j,k,bi,bj)-half)*drF(k))               phiHydLow(i,j,bi,bj) = phiHydC(i,j)
99       &                    *gravity*alphRho(i,j)*recip_rhoConst       &                  +( MIN(zero,ddRloc)*ratioRm
100       &          + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)       &                    +MAX(zero,ddRloc)*ratioRp
101       &          + phi0surf(i,j,bi,bj)       &                   )*gravity*alphRho(i,j)*recip_rhoConst
102             ENDIF             ENDIF
103            ENDDO            ENDDO
104           ENDDO           ENDDO
# Line 104  C           substracted from hFacC Line 106  C           substracted from hFacC
106  C  --  end if integr_GeoPot = ...  C  --  end if integr_GeoPot = ...
107         ENDIF         ENDIF
108                    
109  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C  -- end buoyancyR = Oceanic (z)
110        ELSEIF ( buoyancyRelation .eq. 'OCEANICP' ) THEN        ENDIF
111    
112         IF (integr_GeoPot.EQ.1) THEN        IF (k.EQ.Nr) THEN
113  C  --  Finite Volume Form  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
114    C  --  last level (bottom): rescale (r*) and add surface contribution
115    
116           DO j=jMin,jMax         IF ( buoyancyRelation .EQ. 'OCEANICP' .OR.
117            DO i=iMin,iMax       &      buoyancyRelation .EQ. 'ATMOSPHERIC' ) THEN
118             IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN  C  -- P coordinate : Phi(R_low) is simply at the top :
119              phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)          DO j=jMin,jMax
120       &          + hFacC(i,j,k,bi,bj)*drF(K)*alphRho(i,j)           DO i=iMin,iMax
121       &          + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)             phiHydLow(i,j,bi,bj) = phiHydF(i,j)
      &          + phi0surf(i,j,bi,bj)  
            ENDIF  
           ENDDO  
122           ENDDO           ENDDO
123            ENDDO
124           ENDIF
125    
126         ELSE          DO j=jMin,jMax
127  C  --  Finite Difference Form           DO i=iMin,iMax
128               phiHydLow(i,j,bi,bj) = phiHydLow(i,j,bi,bj)
129  C---------- Compute gravity*(sea surface elevation) first       &            + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)
130  C           This has to be done starting from phiHyd at the current       &            + phi0surf(i,j,bi,bj)
 C           tracer point and .5 of the cell's thickness has to be  
 C           substracted from hFacC  
   
          DO j=jMin,jMax  
           DO i=iMin,iMax  
            IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN  
             phiHydLow(i,j,bi,bj) = phiHyd(i,j,k)  
      &          + ( half*dRloc+(hFacC(i,j,k,bi,bj)-half)*drF(k)  
      &            )*alphRho(i,j)  
      &          + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj)  
      &          + phi0surf(i,j,bi,bj)  
            ENDIF  
           ENDDO  
131           ENDDO           ENDDO
132            ENDDO
133    
 C  --  end if integr_GeoPot = ...  
        ENDIF  
   
 c     ELSEIF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN  
134  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|  C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----|
135    C  -- end if k=Nr.
136        ENDIF        ENDIF
137    
138  #endif /* INCLUDE_PHIHYD_CALCULATION_CODE */  #endif /* INCLUDE_PHIHYD_CALCULATION_CODE */
Legend:
Removed from v.1.1  
changed lines
  Added in v.1.2
  ViewVC Help
Powered by ViewVC 1.1.22