/[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.13 by jmc, Fri Sep 10 17:23:17 2010 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
# Line 46  C     myThid :: Instance number for this Line 47  C     myThid :: Instance number for this
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 57  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
# Line 71  CEOP Line 71  CEOP
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,jMax          DO j=jMin,jMax
77           DO i=iMin,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
# 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,jMax          DO j=jMin,jMax
96           DO i=iMin,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
# Line 115  C-     Consistent with Phi'= Integr[ the Line 115  C-     Consistent with Phi'= Integr[ the
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 160  C--   Zonal & Meridional gradient of pot Line 160  C--   Zonal & Meridional gradient of pot
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_prime/rho0 * Grad_r(g.z)  C--    z* coordinate slope term: rho_prime/rho0 * Grad_r(g.z)
167          factorZ = gravity*recip_rhoConst*recip_rhoFacC(k)*0.5 _d 0          factorZ = gravity*recip_rhoConst*recip_rhoFacC(k)*0.5 _d 0
168          DO j=jMin,jMax          DO j=jMin,jMax
# Line 186  C--    z* coordinate slope term: rho_pri Line 187  C--    z* coordinate slope term: rho_pri
187       &             *recip_dyC(i,j,bi,bj)*recip_deepFacC(k)       &             *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_prime * 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
# Line 205  C--    p* coordinate slope term: alpha_p Line 206  C--    p* coordinate slope term: alpha_p
206       &             *rC(k)*recip_dyC(i,j,bi,bj)*recip_deepFacC(k)       &             *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_prime * 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+1,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)*recip_deepFacC(k)       &             *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           ENDDO
235          ENDDO          ENDDO
236          DO j=jMin+1,jMax          DO j=jMin+1,jMax
237           DO i=iMin,iMax           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)*recip_deepFacC(k)       &             *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.13  
changed lines
  Added in v.1.16
  ViewVC Help
Powered by ViewVC 1.1.22