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

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

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

revision 1.2 by jmc, Sun Feb 9 02:00:50 2003 UTC revision 1.16 by jmc, Tue Apr 29 21:10:48 2014 UTC
# Line 2  C $Header$ Line 2  C $Header$
2  C $Name$  C $Name$
3    
4  #include "CPP_OPTIONS.h"  #include "CPP_OPTIONS.h"
5    #define OLD_PSTAR_SLOPE_TERM
6    
7  CBOP  CBOP
8  C     !ROUTINE: CALC_GRAD_PHI_HYD  C     !ROUTINE: CALC_GRAD_PHI_HYD
9  C     !INTERFACE:  C     !INTERFACE:
10        SUBROUTINE CALC_GRAD_PHI_HYD(        SUBROUTINE CALC_GRAD_PHI_HYD(
11       I                       k, bi, bj, iMin,iMax, jMin,jMax,       I                       k, bi, bj, iMin,iMax, jMin,jMax,
12       I                       phiHyd, alphRho, tFld, sFld,       I                       phiHydC, alphRho, tFld, sFld,
13       O                       dPhiHydX, dPhiHydY,       O                       dPhiHydX, dPhiHydY,
14       I                       myTime, myIter, myThid)       I                       myTime, myIter, myThid)
15  C     !DESCRIPTION: \bv  C     !DESCRIPTION: \bv
16  C     *==========================================================*  C     *==========================================================*
17  C     | S/R CALC_GRAD_PHI_HYD                                      C     | S/R CALC_GRAD_PHI_HYD
18  C     | o Calculate the gradient of Hydrostatic potential anomaly  C     | o Calculate the gradient of Hydrostatic potential anomaly
19  C     *==========================================================*  C     *==========================================================*
20  C     \ev  C     \ev
21    
# Line 26  C     == Global variables == Line 27  C     == Global variables ==
27  #include "PARAMS.h"  #include "PARAMS.h"
28  #include "GRID.h"  #include "GRID.h"
29  #include "SURFACE.h"  #include "SURFACE.h"
30    #include "DYNVARS.h"
31    
32  C     !INPUT/OUTPUT PARAMETERS:  C     !INPUT/OUTPUT PARAMETERS:
33  C     == Routine Arguments ==  C     == Routine Arguments ==
34  C     bi,bj      :: tile index              C     bi,bj      :: tile index
35  C     iMin,iMax,jMin,jMax :: Loop counters  C     iMin,iMax,jMin,jMax :: Loop counters
36  C     phiHyd     :: Hydrostatic Potential anomaly  C     phiHydC    :: Hydrostatic Potential anomaly
37  C                  (atmos: =Geopotential ; ocean-z: =Pressure/rho)  C                  (atmos: =Geopotential ; ocean-z: =Pressure/rho)
38  C     alphRho    :: Density (z-coord) or specific volume (p-coord)  C     alphRho    :: Density (z-coord) or specific volume (p-coord)
39  C     tFld       :: Potential temp.  C     tFld       :: Potential temp.
40  C     sFld       :: Salinity  C     sFld       :: Salinity
41  C     dPhiHydX,Y :: Gradient (X & Y directions) of Hyd. Potential  C     dPhiHydX,Y :: Gradient (X & Y directions) of Hyd. Potential
42  C     myTime :: Current time  C     myTime :: Current time
43  C     myIter :: Current iteration number  C     myIter :: Current iteration number
44  C     myThid :: Instance number for this call of the routine.  C     myThid :: Instance number for this call of the routine.
45        INTEGER k, bi,bj, iMin,iMax, jMin,jMax        INTEGER k, bi,bj, iMin,iMax, jMin,jMax
46        _RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)        _RL phiHydC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
47        _RL alphRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _RL alphRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
48        _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)
49        _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)
50        _RL dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL dPhiHydX(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
51        _RL dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL dPhiHydY(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
52        _RL myTime        _RL myTime
53        INTEGER myIter, myThid        INTEGER myIter, myThid
54    
# Line 56  C     !LOCAL VARIABLES: Line 58  C     !LOCAL VARIABLES:
58  C     == Local variables ==  C     == Local variables ==
59  C     i,j :: Loop counters  C     i,j :: Loop counters
60        INTEGER i,j        INTEGER i,j
61        _RL varLoc(1-Olx:sNx+Olx,1-Oly:sNy+Oly)        _RL varLoc(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
62    #ifdef NONLIN_FRSURF
63          _RL factorZ, factorP, factPI
64          CHARACTER*(MAX_LEN_MBUF) msgBuf
65    #endif
66  CEOP  CEOP
67    
68        DO j=jMin-1,jMax  #ifdef NONLIN_FRSURF
69         DO i=iMin-1,iMax        IF (select_rStar.GE.2 .AND. nonlinFreeSurf.GE.4 ) THEN
70          varLoc(i,j) = phiHyd(i,j,k)+phi0surf(i,j,bi,bj)  # ifndef DISABLE_RSTAR_CODE
71  c       varLoc(i,j) = phiHyd(i,j,k)  C-    Integral of b.dr = rStarFac * Integral of b.dr* :
72    C      and will add later (select_rStar=2) the contribution of
73    C      the slope of the r* coordinate.
74           IF ( fluidIsAir ) THEN
75    C-     Consistent with Phi'= Integr[ theta'.dPI ] :
76            DO j=jMin,jMax
77             DO i=iMin,iMax
78              varLoc(i,j) = phiHydC(i,j)*pStarFacK(i,j,bi,bj)
79         &                + phi0surf(i,j,bi,bj)
80             ENDDO
81            ENDDO
82           ELSE
83            DO j=jMin,jMax
84             DO i=iMin,iMax
85              varLoc(i,j) = phiHydC(i,j)*rStarFacC(i,j,bi,bj)
86         &                + phi0surf(i,j,bi,bj)
87             ENDDO
88            ENDDO
89           ENDIF
90          ELSEIF (select_rStar.GE.1 .AND. nonlinFreeSurf.GE.4 ) THEN
91    C-    Integral of b.dr but scaled to correspond to a fixed r-level (=r*)
92    C      no contribution of the slope of the r* coordinate (select_rStar=1)
93           IF ( fluidIsAir ) THEN
94    C-     Consistent with Phi'= Integr[ theta'.dPI ] :
95            DO j=jMin,jMax
96             DO i=iMin,iMax
97              IF (Ro_surf(i,j,bi,bj).EQ.rC(k)) THEN
98               factPI=atm_Cp*( ((etaH(i,j,bi,bj)+rC(k))/atm_Po)**atm_kappa
99         &                    -(                 rC(k) /atm_Po)**atm_kappa
100         &                  )
101               varLoc(i,j) = factPI*alphRho(i,j)
102         &                 + phi0surf(i,j,bi,bj)
103              ELSEIF (Ro_surf(i,j,bi,bj).NE.0. _d 0) THEN
104               factPI = (rC(k)/Ro_surf(i,j,bi,bj))**atm_kappa
105               varLoc(i,j) = phiHydC(i,j)
106         &                  *(rStarFacC(i,j,bi,bj)**atm_kappa - factPI)
107         &                  /(1. _d 0 -factPI)
108         &                 + phi0surf(i,j,bi,bj)
109              ENDIF
110             ENDDO
111            ENDDO
112           ELSE
113            DO j=jMin,jMax
114             DO i=iMin,iMax
115              IF (Ro_surf(i,j,bi,bj).EQ.rC(k)) THEN
116               WRITE(msgBuf,'(3A)') 'CALC_GRAD_PHI_HYD: ',
117         &      'Problem when Ro_surf=rC',
118         &      ' with select_rStar,nonlinFreeSurf=1,4'
119               CALL PRINT_ERROR( msgBuf , myThid)
120               STOP 'CALC_GRAD_PHI_HYD: Pb in r* options implementation'
121              ELSE
122               varLoc(i,j) = phiHydC(i,j)
123         &                  *(etaH(i,j,bi,bj)+Ro_surf(i,j,bi,bj)-rC(k))
124         &                  /                (Ro_surf(i,j,bi,bj)-rC(k))
125         &                 + phi0surf(i,j,bi,bj)
126              ENDIF
127             ENDDO
128            ENDDO
129           ENDIF
130    # endif /* DISABLE_RSTAR_CODE */
131          ELSE
132    #else /* NONLIN_FRSURF */
133          IF (.TRUE.) THEN
134    #endif /* NONLIN_FRSURF */
135           DO j=jMin,jMax
136            DO i=iMin,iMax
137             varLoc(i,j) = phiHydC(i,j)+phi0surf(i,j,bi,bj)
138            ENDDO
139         ENDDO         ENDDO
140        ENDDO        ENDIF
141    
142  C-    Zonal gradient  C--   Zonal & Meridional gradient of potential anomaly
143          DO j=1-OLy,sNy+OLy
144           DO i=1-OLx,sNx+OLx
145            dPhiHydX(i,j)  = 0. _d 0
146            dPhiHydY(i,j)  = 0. _d 0
147           ENDDO
148          ENDDO
149        DO j=jMin,jMax        DO j=jMin,jMax
150           DO i=iMin+1,iMax
151            dPhiHydX(i,j) = _recip_dxC(i,j,bi,bj)*recip_deepFacC(k)
152         &                *( varLoc(i,j)-varLoc(i-1,j) )*recip_rhoFacC(k)
153           ENDDO
154          ENDDO
155          DO j=jMin+1,jMax
156         DO i=iMin,iMax         DO i=iMin,iMax
157          dPhiHydX(i,j) = _recip_dxC(i,j,bi,bj)          dPhiHydY(i,j) = _recip_dyC(i,j,bi,bj)*recip_deepFacC(k)
158       &                *( varLoc(i,j)-varLoc(i-1,j) )       &                *( varLoc(i,j)-varLoc(i,j-1) )*recip_rhoFacC(k)
 c    &              *( (phiHyd(i,j,k)-phiHyd(i-1,j,k))  
 c    &                +(phi0surf(i,j,bi,bj)-phi0surf(i-1,j,bi,bj))  
 c    &              *( (phiHyd(i,j,k)+phi0surf(i,j,bi,bj))  
 c    &                -(phiHyd(i-1,j,k)+phi0surf(i-1,j,bi,bj))  
 c    &               )  
159         ENDDO         ENDDO
160        ENDDO        ENDDO
161    
162  C-    Meridional gradient  #ifdef NONLIN_FRSURF
163        DO j=jMin,jMax  # ifndef DISABLE_RSTAR_CODE
164         DO i=iMin,iMax        IF (select_rStar.GE.2 .AND. nonlinFreeSurf.GE.1 ) THEN
165          dPhiHydY(i,j) = _recip_dyC(i,j,bi,bj)         IF ( fluidIsWater .AND. usingZCoords ) THEN
166       &                *( varLoc(i,j)-varLoc(i,j-1) )  C--    z* coordinate slope term: rho_prime/rho0 * Grad_r(g.z)
167  c    &              *( (phiHyd(i,j,k)-phiHyd(i,j-1,k))          factorZ = gravity*recip_rhoConst*recip_rhoFacC(k)*0.5 _d 0
168  c    &                +(phi0surf(i,j,bi,bj)-phi0surf(i,j-1,bi,bj))          DO j=jMin,jMax
169  c    &              *( (phiHyd(i,j,k)+phi0surf(i,j,bi,bj))           DO i=iMin,iMax
170  c    &                -(phiHyd(i,j-1,k)+phi0surf(i,j-1,bi,bj))            varLoc(i,j) = etaH(i,j,bi,bj)
171  c    &               )       &                *(1. _d 0 + rC(k)*recip_Rcol(i,j,bi,bj))
172             ENDDO
173            ENDDO
174            DO j=jMin,jMax
175             DO i=iMin+1,iMax
176              dPhiHydX(i,j) = dPhiHydX(i,j)
177         &     +factorZ*(alphRho(i-1,j)+alphRho(i,j))
178         &             *(varLoc(i,j)-varLoc(i-1,j))
179         &             *recip_dxC(i,j,bi,bj)*recip_deepFacC(k)
180             ENDDO
181            ENDDO
182            DO j=jMin+1,jMax
183             DO i=iMin,iMax
184              dPhiHydY(i,j) = dPhiHydY(i,j)
185         &     +factorZ*(alphRho(i,j-1)+alphRho(i,j))
186         &             *(varLoc(i,j)-varLoc(i,j-1))
187         &             *recip_dyC(i,j,bi,bj)*recip_deepFacC(k)
188             ENDDO
189            ENDDO
190           ELSEIF ( fluidIsWater ) THEN
191    C--    p* coordinate slope term: alpha_prime * Grad_r( p )
192            factorP = 0.5 _d 0
193            DO j=jMin,jMax
194             DO i=iMin+1,iMax
195              dPhiHydX(i,j) = dPhiHydX(i,j)
196         &     +factorP*(alphRho(i-1,j)+alphRho(i,j))
197         &             *(rStarFacC(i,j,bi,bj)-rStarFacC(i-1,j,bi,bj))
198         &             *rC(k)*recip_dxC(i,j,bi,bj)*recip_deepFacC(k)
199             ENDDO
200            ENDDO
201            DO j=jMin+1,jMax
202             DO i=iMin,iMax
203              dPhiHydY(i,j) = dPhiHydY(i,j)
204         &     +factorP*(alphRho(i,j-1)+alphRho(i,j))
205         &             *(rStarFacC(i,j,bi,bj)-rStarFacC(i,j-1,bi,bj))
206         &             *rC(k)*recip_dyC(i,j,bi,bj)*recip_deepFacC(k)
207             ENDDO
208            ENDDO
209           ELSEIF ( fluidIsAir ) THEN
210    #ifdef OLD_PSTAR_SLOPE_TERM
211    C--    p* coordinate slope term: alpha_prime * Grad_r( p ):
212    C      PI_star * (Theta_eq^prime)_bar_i * kappa * delta^i( rStarFacC )
213    C- Note: factor: ( p_s / p_s^o )^(kappa - 1) = rStarFacC^(kappa -1)
214    C        is missing here.
215            factorP = (rC(k)/atm_Po)**atm_kappa
216            factorP = (atm_Rd/rC(k))*factorP*0.5 _d 0
217    #else
218    C--    p* coordinate slope term: theta_prime * Grad_r( PI ):
219    C      PI_star * (Theta_eq^prime)_bar_i * delta^i( rStarFacC^kappa )
220    C      This is also consitent with geopotential factor: rStarFacC^kappa
221            factorP = halfRL*atm_Cp*(rC(k)/atm_Po)**atm_kappa
222    #endif
223            DO j=jMin,jMax
224             DO i=iMin+1,iMax
225              dPhiHydX(i,j) = dPhiHydX(i,j)
226         &     +factorP*(alphRho(i-1,j)+alphRho(i,j))
227    #ifdef OLD_PSTAR_SLOPE_TERM
228         &             *(rStarFacC(i,j,bi,bj)-rStarFacC(i-1,j,bi,bj))
229         &             *rC(k)*recip_dxC(i,j,bi,bj)*recip_deepFacC(k)
230    #else
231         &             *(pStarFacK(i,j,bi,bj)-pStarFacK(i-1,j,bi,bj))
232         &             *recip_dxC(i,j,bi,bj)*recip_deepFacC(k)
233    #endif
234             ENDDO
235            ENDDO
236            DO j=jMin+1,jMax
237             DO i=iMin,iMax
238              dPhiHydY(i,j) = dPhiHydY(i,j)
239         &     +factorP*(alphRho(i,j-1)+alphRho(i,j))
240    #ifdef OLD_PSTAR_SLOPE_TERM
241         &             *(rStarFacC(i,j,bi,bj)-rStarFacC(i,j-1,bi,bj))
242         &             *rC(k)*recip_dyC(i,j,bi,bj)*recip_deepFacC(k)
243    #else
244         &             *(pStarFacK(i,j,bi,bj)-pStarFacK(i,j-1,bi,bj))
245         &             *recip_dyC(i,j,bi,bj)*recip_deepFacC(k)
246    #endif
247             ENDDO
248            ENDDO
249           ENDIF
250          ENDIF
251    # endif /* DISABLE_RSTAR_CODE */
252    #endif /* NONLIN_FRSURF */
253    
254    C--   Apply mask:
255          DO j=1-OLy,sNy+OLy
256           DO i=1-OLx,sNx+OLx
257             dPhiHydX(i,j) = dPhiHydX(i,j)*_maskW(i,j,k,bi,bj)
258             dPhiHydY(i,j) = dPhiHydY(i,j)*_maskS(i,j,k,bi,bj)
259         ENDDO         ENDDO
260        ENDDO        ENDDO
261    
Legend:
Removed from v.1.2  
changed lines
  Added in v.1.16
  ViewVC Help
Powered by ViewVC 1.1.22