/[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.7 by jmc, Sat Nov 5 01:00:57 2005 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                       phiHydC, 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 30  C     == Global variables == Line 31  C     == Global variables ==
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     phiHydC    :: 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
 c     _RL phiHyd(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr)  
46        _RL phiHydC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)        _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 58  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  #ifdef NONLIN_FRSURF
63        _RL factorZ, factorP, conv_theta2T        _RL factorZ, factorP, factPI
       _RL factPI  
64        CHARACTER*(MAX_LEN_MBUF) msgBuf        CHARACTER*(MAX_LEN_MBUF) msgBuf
65  #endif  #endif
66  CEOP  CEOP
67    
68  #ifdef NONLIN_FRSURF  #ifdef NONLIN_FRSURF
69        IF (select_rStar.GE.2 .AND. nonlinFreeSurf.GE.4 ) THEN        IF (select_rStar.GE.2 .AND. nonlinFreeSurf.GE.4 ) THEN
70    # ifndef DISABLE_RSTAR_CODE
71  C-    Integral of b.dr = rStarFac * Integral of b.dr* :  C-    Integral of b.dr = rStarFac * Integral of b.dr* :
72  C      and will add later (select_rStar=2) the contribution of  C      and will add later (select_rStar=2) the contribution of
73  C      the slope of the r* coordinate.  C      the slope of the r* coordinate.
74         IF ( buoyancyRelation .EQ. 'ATMOSPHERIC' ) THEN         IF ( fluidIsAir ) THEN
75  C-     Consistent with Phi'= Integr[ theta'.dPi ] :  C-     Consistent with Phi'= Integr[ theta'.dPI ] :
76          DO j=jMin-1,jMax          DO j=jMin,jMax
77           DO i=iMin-1,iMax           DO i=iMin,iMax
78            varLoc(i,j) = phiHydC(i,j)*rStarFacC(i,j,bi,bj)**atm_kappa            varLoc(i,j) = phiHydC(i,j)*pStarFacK(i,j,bi,bj)
79       &                + phi0surf(i,j,bi,bj)       &                + phi0surf(i,j,bi,bj)
80           ENDDO           ENDDO
81          ENDDO          ENDDO
82         ELSE         ELSE
83          DO j=jMin-1,jMax          DO j=jMin,jMax
84           DO i=iMin-1,iMax           DO i=iMin,iMax
85            varLoc(i,j) = phiHydC(i,j)*rStarFacC(i,j,bi,bj)            varLoc(i,j) = phiHydC(i,j)*rStarFacC(i,j,bi,bj)
86       &                + phi0surf(i,j,bi,bj)       &                + phi0surf(i,j,bi,bj)
87           ENDDO           ENDDO
# Line 90  C-     Consistent with Phi'= Integr[ the Line 90  C-     Consistent with Phi'= Integr[ the
90        ELSEIF (select_rStar.GE.1 .AND. nonlinFreeSurf.GE.4 ) THEN        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*)  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)  C      no contribution of the slope of the r* coordinate (select_rStar=1)
93         IF ( buoyancyRelation .EQ. 'ATMOSPHERIC' ) THEN         IF ( fluidIsAir ) THEN
94  C-     Consistent with Phi'= Integr[ theta'.dPi ] :  C-     Consistent with Phi'= Integr[ theta'.dPI ] :
95          DO j=jMin-1,jMax          DO j=jMin,jMax
96           DO i=iMin-1,iMax           DO i=iMin,iMax
97            IF (Ro_surf(i,j,bi,bj).EQ.rC(k)) THEN            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               factPI=atm_Cp*( ((etaH(i,j,bi,bj)+rC(k))/atm_Po)**atm_kappa
99       &                    -(                 rC(k) /atm_Po)**atm_kappa       &                    -(                 rC(k) /atm_Po)**atm_kappa
100       &                  )       &                  )
101             varLoc(i,j) = factPI*alphRho(i,j)             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            ELSEIF (Ro_surf(i,j,bi,bj).NE.0. _d 0) THEN
104             factPI = (rC(k)/Ro_surf(i,j,bi,bj))**atm_kappa             factPI = (rC(k)/Ro_surf(i,j,bi,bj))**atm_kappa
105             varLoc(i,j) = phiHydC(i,j)             varLoc(i,j) = phiHydC(i,j)
# Line 109  C-     Consistent with Phi'= Integr[ the Line 110  C-     Consistent with Phi'= Integr[ the
110           ENDDO           ENDDO
111          ENDDO          ENDDO
112         ELSE         ELSE
113          DO j=jMin-1,jMax          DO j=jMin,jMax
114           DO i=iMin-1,iMax           DO i=iMin,iMax
115            IF (Ro_surf(i,j,bi,bj).EQ.rC(k)) THEN            IF (Ro_surf(i,j,bi,bj).EQ.rC(k)) THEN
116             WRITE(msgBuf,'(3A)') 'CALC_GRAD_PHI_HYD: ',             WRITE(msgBuf,'(3A)') 'CALC_GRAD_PHI_HYD: ',
117       &      'Problem when Ro_surf=rC',       &      'Problem when Ro_surf=rC',
118       &      ' with select_rStar,integr_GeoPot=1,4'       &      ' with select_rStar,nonlinFreeSurf=1,4'
119             CALL PRINT_ERROR( msgBuf , myThid)             CALL PRINT_ERROR( msgBuf , myThid)
120             STOP 'CALC_GRAD_PHI_HYD: Pb in r* options implementation'             STOP 'CALC_GRAD_PHI_HYD: Pb in r* options implementation'
121            ELSE            ELSE
# Line 126  C-     Consistent with Phi'= Integr[ the Line 127  C-     Consistent with Phi'= Integr[ the
127           ENDDO           ENDDO
128          ENDDO          ENDDO
129         ENDIF         ENDIF
130    # endif /* DISABLE_RSTAR_CODE */
131        ELSE        ELSE
132  #else /* NONLIN_FRSURF */  #else /* NONLIN_FRSURF */
133        IF (.TRUE.) THEN        IF (.TRUE.) THEN
134  #endif /* NONLIN_FRSURF */  #endif /* NONLIN_FRSURF */
135         DO j=jMin-1,jMax         DO j=jMin,jMax
136          DO i=iMin-1,iMax          DO i=iMin,iMax
137           varLoc(i,j) = phiHydC(i,j)+phi0surf(i,j,bi,bj)           varLoc(i,j) = phiHydC(i,j)+phi0surf(i,j,bi,bj)
138          ENDDO          ENDDO
139         ENDDO         ENDDO
140        ENDIF        ENDIF
141    
142  C--   Zonal & Meridional gradient of potential anomaly  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)
         dPhiHydY(i,j) = _recip_dyC(i,j,bi,bj)  
      &                *( varLoc(i,j)-varLoc(i,j-1) )  
159         ENDDO         ENDDO
160        ENDDO        ENDDO
161    
162  #ifdef NONLIN_FRSURF  #ifdef NONLIN_FRSURF
163    # ifndef DISABLE_RSTAR_CODE
164        IF (select_rStar.GE.2 .AND. nonlinFreeSurf.GE.1 ) THEN        IF (select_rStar.GE.2 .AND. nonlinFreeSurf.GE.1 ) THEN
165         IF ( buoyancyRelation .EQ. 'OCEANIC' ) THEN         IF ( fluidIsWater .AND. usingZCoords ) THEN
166  C--    z* coordinate slope term: rho'/rho0 * Grad_r(g.z)  C--    z* coordinate slope term: rho_prime/rho0 * Grad_r(g.z)
167          factorZ = gravity*recip_rhoConst*0.5 _d 0          factorZ = gravity*recip_rhoConst*recip_rhoFacC(k)*0.5 _d 0
168          DO j=jMin-1,jMax          DO j=jMin,jMax
169           DO i=iMin-1,iMax           DO i=iMin,iMax
170            varLoc(i,j) = etaH(i,j,bi,bj)            varLoc(i,j) = etaH(i,j,bi,bj)
171       &                *(1. _d 0 + rC(k)*recip_Rcol(i,j,bi,bj))       &                *(1. _d 0 + rC(k)*recip_Rcol(i,j,bi,bj))
172           ENDDO           ENDDO
173          ENDDO          ENDDO
174          DO j=jMin,jMax          DO j=jMin,jMax
175           DO i=iMin,iMax           DO i=iMin+1,iMax
176            dPhiHydX(i,j) = dPhiHydX(i,j)            dPhiHydX(i,j) = dPhiHydX(i,j)
177       &     +factorZ*(alphRho(i-1,j)+alphRho(i,j))       &     +factorZ*(alphRho(i-1,j)+alphRho(i,j))
178       &             *(varLoc(i,j)-varLoc(i-1,j))       &             *(varLoc(i,j)-varLoc(i-1,j))
179       &             *recip_dxC(i,j,bi,bj)       &             *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)            dPhiHydY(i,j) = dPhiHydY(i,j)
185       &     +factorZ*(alphRho(i,j-1)+alphRho(i,j))       &     +factorZ*(alphRho(i,j-1)+alphRho(i,j))
186       &             *(varLoc(i,j)-varLoc(i,j-1))       &             *(varLoc(i,j)-varLoc(i,j-1))
187       &             *recip_dyC(i,j,bi,bj)         &             *recip_dyC(i,j,bi,bj)*recip_deepFacC(k)
188           ENDDO           ENDDO
189          ENDDO          ENDDO
190         ELSEIF (buoyancyRelation .EQ. 'OCEANICP' ) THEN         ELSEIF ( fluidIsWater ) THEN
191  C--    p* coordinate slope term: alpha' * Grad_r( p )  C--    p* coordinate slope term: alpha_prime * Grad_r( p )
192          factorP = 0.5 _d 0          factorP = 0.5 _d 0
193          DO j=jMin,jMax          DO j=jMin,jMax
194           DO i=iMin,iMax           DO i=iMin+1,iMax
195            dPhiHydX(i,j) = dPhiHydX(i,j)            dPhiHydX(i,j) = dPhiHydX(i,j)
196       &     +factorP*(alphRho(i-1,j)+alphRho(i,j))       &     +factorP*(alphRho(i-1,j)+alphRho(i,j))
197       &             *(rStarFacC(i,j,bi,bj)-rStarFacC(i-1,j,bi,bj))       &             *(rStarFacC(i,j,bi,bj)-rStarFacC(i-1,j,bi,bj))
198       &             *rC(k)*recip_dxC(i,j,bi,bj)         &             *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)            dPhiHydY(i,j) = dPhiHydY(i,j)
204       &     +factorP*(alphRho(i,j-1)+alphRho(i,j))       &     +factorP*(alphRho(i,j-1)+alphRho(i,j))
205       &             *(rStarFacC(i,j,bi,bj)-rStarFacC(i,j-1,bi,bj))       &             *(rStarFacC(i,j,bi,bj)-rStarFacC(i,j-1,bi,bj))
206       &             *rC(k)*recip_dyC(i,j,bi,bj)         &             *rC(k)*recip_dyC(i,j,bi,bj)*recip_deepFacC(k)
207           ENDDO           ENDDO
208          ENDDO          ENDDO
209         ELSEIF ( buoyancyRelation .EQ. 'ATMOSPHERIC' ) THEN         ELSEIF ( fluidIsAir ) THEN
210  C--    p* coordinate slope term: alpha' * Grad_r( p )  #ifdef OLD_PSTAR_SLOPE_TERM
211          conv_theta2T = (rC(k)/atm_Po)**atm_kappa  C--    p* coordinate slope term: alpha_prime * Grad_r( p ):
212          factorP = (atm_Rd/rC(k))*conv_theta2T*0.5 _d 0  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          DO j=jMin,jMax
224           DO i=iMin,iMax           DO i=iMin+1,iMax
225            dPhiHydX(i,j) = dPhiHydX(i,j)            dPhiHydX(i,j) = dPhiHydX(i,j)
226       &     +factorP*(alphRho(i-1,j)+alphRho(i,j))       &     +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))       &             *(rStarFacC(i,j,bi,bj)-rStarFacC(i-1,j,bi,bj))
229       &             *rC(k)*recip_dxC(i,j,bi,bj)         &             *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)            dPhiHydY(i,j) = dPhiHydY(i,j)
239       &     +factorP*(alphRho(i,j-1)+alphRho(i,j))       &     +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))       &             *(rStarFacC(i,j,bi,bj)-rStarFacC(i,j-1,bi,bj))
242       &             *rC(k)*recip_dyC(i,j,bi,bj)         &             *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           ENDDO
248          ENDDO          ENDDO
249         ENDIF         ENDIF
250        ENDIF        ENDIF
251    # endif /* DISABLE_RSTAR_CODE */
252  #endif /* NONLIN_FRSURF */  #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
260          ENDDO
261    
262  #endif /* INCLUDE_PHIHYD_CALCULATION_CODE */  #endif /* INCLUDE_PHIHYD_CALCULATION_CODE */
263    
264        RETURN        RETURN
Legend:
Removed from v.1.7  
changed lines
  Added in v.1.16
  ViewVC Help
Powered by ViewVC 1.1.22