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

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

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


Revision 1.10 - (hide annotations) (download)
Mon Feb 5 03:22:44 2007 UTC (17 years, 3 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint62c, checkpoint59, checkpoint62a, checkpoint58y_post, checkpoint60, checkpoint61, checkpoint62, checkpoint58w_post, checkpoint59q, checkpoint59p, checkpoint59r, checkpoint59e, checkpoint59d, checkpoint59g, checkpoint59f, checkpoint59a, checkpoint59c, checkpoint59b, checkpoint59m, checkpoint59l, checkpoint59o, checkpoint59n, checkpoint59i, checkpoint59h, checkpoint59k, checkpoint62b, checkpoint58v_post, checkpoint61f, checkpoint58x_post, checkpoint61n, checkpoint59j, checkpoint61q, checkpoint61e, checkpoint61g, checkpoint61d, checkpoint61b, checkpoint61c, checkpoint61a, checkpoint61o, checkpoint61l, checkpoint61m, checkpoint61j, checkpoint61k, checkpoint61h, checkpoint61i, checkpoint61v, checkpoint61w, checkpoint61t, checkpoint61u, checkpoint61r, checkpoint61s, checkpoint61p, checkpoint61z, checkpoint61x, checkpoint61y
Changes since 1.9: +33 -17 lines
calculate pressure gradient only where necessary (avoid Pb with uninitialized array)

1 jmc 1.10 C $Header: /u/gcmpack/MITgcm/model/src/calc_grad_phi_hyd.F,v 1.9 2006/12/05 05:25:08 jmc Exp $
2 jmc 1.1 C $Name: $
3    
4     #include "CPP_OPTIONS.h"
5    
6     CBOP
7     C !ROUTINE: CALC_GRAD_PHI_HYD
8     C !INTERFACE:
9 jmc 1.9 SUBROUTINE CALC_GRAD_PHI_HYD(
10 jmc 1.1 I k, bi, bj, iMin,iMax, jMin,jMax,
11 jmc 1.3 I phiHydC, alphRho, tFld, sFld,
12 jmc 1.1 O dPhiHydX, dPhiHydY,
13     I myTime, myIter, myThid)
14     C !DESCRIPTION: \bv
15     C *==========================================================*
16 jmc 1.9 C | S/R CALC_GRAD_PHI_HYD
17     C | o Calculate the gradient of Hydrostatic potential anomaly
18 jmc 1.1 C *==========================================================*
19     C \ev
20    
21     C !USES:
22     IMPLICIT NONE
23     C == Global variables ==
24     #include "SIZE.h"
25     #include "EEPARAMS.h"
26     #include "PARAMS.h"
27     #include "GRID.h"
28     #include "SURFACE.h"
29 jmc 1.3 #include "DYNVARS.h"
30 jmc 1.1
31     C !INPUT/OUTPUT PARAMETERS:
32     C == Routine Arguments ==
33 jmc 1.9 C bi,bj :: tile index
34 jmc 1.1 C iMin,iMax,jMin,jMax :: Loop counters
35 jmc 1.9 C phiHydC :: Hydrostatic Potential anomaly
36 jmc 1.1 C (atmos: =Geopotential ; ocean-z: =Pressure/rho)
37     C alphRho :: Density (z-coord) or specific volume (p-coord)
38     C tFld :: Potential temp.
39 jmc 1.9 C sFld :: Salinity
40 jmc 1.1 C dPhiHydX,Y :: Gradient (X & Y directions) of Hyd. Potential
41     C myTime :: Current time
42     C myIter :: Current iteration number
43     C myThid :: Instance number for this call of the routine.
44     INTEGER k, bi,bj, iMin,iMax, jMin,jMax
45 jmc 1.3 _RL phiHydC(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
46 jmc 1.1 _RL alphRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy)
47     _RL tFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
48     _RL sFld(1-OLx:sNx+OLx,1-OLy:sNy+OLy,Nr,nSx,nSy)
49     _RL dPhiHydX(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
50     _RL dPhiHydY(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
51     _RL myTime
52     INTEGER myIter, myThid
53    
54     #ifdef INCLUDE_PHIHYD_CALCULATION_CODE
55    
56     C !LOCAL VARIABLES:
57     C == Local variables ==
58     C i,j :: Loop counters
59     INTEGER i,j
60     _RL varLoc(1-Olx:sNx+Olx,1-Oly:sNy+Oly)
61 jmc 1.7 #ifdef NONLIN_FRSURF
62 jmc 1.3 _RL factorZ, factorP, conv_theta2T
63 jmc 1.5 _RL factPI
64     CHARACTER*(MAX_LEN_MBUF) msgBuf
65 jmc 1.7 #endif
66 jmc 1.1 CEOP
67    
68 jmc 1.3 #ifdef NONLIN_FRSURF
69 jmc 1.5 IF (select_rStar.GE.2 .AND. nonlinFreeSurf.GE.4 ) THEN
70 heimbach 1.8 # ifndef DISABLE_RSTAR_CODE
71 jmc 1.5 C- Integral of b.dr = rStarFac * Integral of b.dr* :
72 jmc 1.9 C and will add later (select_rStar=2) the contribution of
73 jmc 1.5 C the slope of the r* coordinate.
74     IF ( buoyancyRelation .EQ. 'ATMOSPHERIC' ) THEN
75     C- Consistent with Phi'= Integr[ theta'.dPi ] :
76 jmc 1.10 DO j=jMin,jMax
77     DO i=iMin,iMax
78 jmc 1.5 varLoc(i,j) = phiHydC(i,j)*rStarFacC(i,j,bi,bj)**atm_kappa
79     & + phi0surf(i,j,bi,bj)
80     ENDDO
81     ENDDO
82     ELSE
83 jmc 1.10 DO j=jMin,jMax
84     DO i=iMin,iMax
85 jmc 1.5 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 ( buoyancyRelation .EQ. 'ATMOSPHERIC' ) THEN
94     C- Consistent with Phi'= Integr[ theta'.dPi ] :
95 jmc 1.10 DO j=jMin,jMax
96     DO i=iMin,iMax
97 jmc 1.5 IF (Ro_surf(i,j,bi,bj).EQ.rC(k)) THEN
98 jmc 1.9 factPI=atm_Cp*( ((etaH(i,j,bi,bj)+rC(k))/atm_Po)**atm_kappa
99     & -( rC(k) /atm_Po)**atm_kappa
100 jmc 1.5 & )
101     varLoc(i,j) = factPI*alphRho(i,j)
102     ELSEIF (Ro_surf(i,j,bi,bj).NE.0. _d 0) THEN
103     factPI = (rC(k)/Ro_surf(i,j,bi,bj))**atm_kappa
104     varLoc(i,j) = phiHydC(i,j)
105     & *(rStarFacC(i,j,bi,bj)**atm_kappa - factPI)
106     & /(1. _d 0 -factPI)
107     & + phi0surf(i,j,bi,bj)
108     ENDIF
109     ENDDO
110     ENDDO
111     ELSE
112 jmc 1.10 DO j=jMin,jMax
113     DO i=iMin,iMax
114 jmc 1.5 IF (Ro_surf(i,j,bi,bj).EQ.rC(k)) THEN
115 jmc 1.6 WRITE(msgBuf,'(3A)') 'CALC_GRAD_PHI_HYD: ',
116     & 'Problem when Ro_surf=rC',
117     & ' with select_rStar,integr_GeoPot=1,4'
118 jmc 1.5 CALL PRINT_ERROR( msgBuf , myThid)
119     STOP 'CALC_GRAD_PHI_HYD: Pb in r* options implementation'
120     ELSE
121     varLoc(i,j) = phiHydC(i,j)
122     & *(etaH(i,j,bi,bj)+Ro_surf(i,j,bi,bj)-rC(k))
123     & / (Ro_surf(i,j,bi,bj)-rC(k))
124     & + phi0surf(i,j,bi,bj)
125     ENDIF
126     ENDDO
127 jmc 1.3 ENDDO
128 jmc 1.5 ENDIF
129 heimbach 1.8 # endif /* DISABLE_RSTAR_CODE */
130 jmc 1.3 ELSE
131     #else /* NONLIN_FRSURF */
132     IF (.TRUE.) THEN
133     #endif /* NONLIN_FRSURF */
134 jmc 1.10 DO j=jMin,jMax
135     DO i=iMin,iMax
136 jmc 1.3 varLoc(i,j) = phiHydC(i,j)+phi0surf(i,j,bi,bj)
137     ENDDO
138     ENDDO
139     ENDIF
140 jmc 1.1
141 jmc 1.3 C-- Zonal & Meridional gradient of potential anomaly
142 jmc 1.1 DO j=jMin,jMax
143 jmc 1.10 DO i=iMin+1,iMax
144 jmc 1.9 dPhiHydX(i,j) = _recip_dxC(i,j,bi,bj)*recip_deepFacC(k)
145     & *( varLoc(i,j)-varLoc(i-1,j) )*recip_rhoFacC(k)
146 jmc 1.10 ENDDO
147     ENDDO
148     DO j=jMin+1,jMax
149     DO i=iMin,iMax
150 jmc 1.9 dPhiHydY(i,j) = _recip_dyC(i,j,bi,bj)*recip_deepFacC(k)
151     & *( varLoc(i,j)-varLoc(i,j-1) )*recip_rhoFacC(k)
152 jmc 1.1 ENDDO
153     ENDDO
154 jmc 1.3
155     #ifdef NONLIN_FRSURF
156     IF (select_rStar.GE.2 .AND. nonlinFreeSurf.GE.1 ) THEN
157     IF ( buoyancyRelation .EQ. 'OCEANIC' ) THEN
158     C-- z* coordinate slope term: rho'/rho0 * Grad_r(g.z)
159 jmc 1.9 factorZ = gravity*recip_rhoConst*recip_rhoFacC(k)*0.5 _d 0
160 jmc 1.10 DO j=jMin,jMax
161     DO i=iMin,iMax
162 jmc 1.3 varLoc(i,j) = etaH(i,j,bi,bj)
163     & *(1. _d 0 + rC(k)*recip_Rcol(i,j,bi,bj))
164     ENDDO
165     ENDDO
166     DO j=jMin,jMax
167 jmc 1.10 DO i=iMin+1,iMax
168 jmc 1.3 dPhiHydX(i,j) = dPhiHydX(i,j)
169     & +factorZ*(alphRho(i-1,j)+alphRho(i,j))
170     & *(varLoc(i,j)-varLoc(i-1,j))
171 jmc 1.9 & *recip_dxC(i,j,bi,bj)*recip_deepFacC(k)
172 jmc 1.10 ENDDO
173     ENDDO
174     DO j=jMin+1,jMax
175     DO i=iMin,iMax
176 jmc 1.3 dPhiHydY(i,j) = dPhiHydY(i,j)
177     & +factorZ*(alphRho(i,j-1)+alphRho(i,j))
178     & *(varLoc(i,j)-varLoc(i,j-1))
179 jmc 1.9 & *recip_dyC(i,j,bi,bj)*recip_deepFacC(k)
180 jmc 1.3 ENDDO
181     ENDDO
182     ELSEIF (buoyancyRelation .EQ. 'OCEANICP' ) THEN
183     C-- p* coordinate slope term: alpha' * Grad_r( p )
184     factorP = 0.5 _d 0
185     DO j=jMin,jMax
186 jmc 1.10 DO i=iMin+1,iMax
187 jmc 1.3 dPhiHydX(i,j) = dPhiHydX(i,j)
188     & +factorP*(alphRho(i-1,j)+alphRho(i,j))
189     & *(rStarFacC(i,j,bi,bj)-rStarFacC(i-1,j,bi,bj))
190 jmc 1.9 & *rC(k)*recip_dxC(i,j,bi,bj)*recip_deepFacC(k)
191 jmc 1.10 ENDDO
192     ENDDO
193     DO j=jMin+1,jMax
194     DO i=iMin,iMax
195 jmc 1.3 dPhiHydY(i,j) = dPhiHydY(i,j)
196     & +factorP*(alphRho(i,j-1)+alphRho(i,j))
197     & *(rStarFacC(i,j,bi,bj)-rStarFacC(i,j-1,bi,bj))
198 jmc 1.9 & *rC(k)*recip_dyC(i,j,bi,bj)*recip_deepFacC(k)
199 jmc 1.3 ENDDO
200     ENDDO
201     ELSEIF ( buoyancyRelation .EQ. 'ATMOSPHERIC' ) THEN
202     C-- p* coordinate slope term: alpha' * Grad_r( p )
203     conv_theta2T = (rC(k)/atm_Po)**atm_kappa
204     factorP = (atm_Rd/rC(k))*conv_theta2T*0.5 _d 0
205     DO j=jMin,jMax
206 jmc 1.10 DO i=iMin+1,iMax
207 jmc 1.3 dPhiHydX(i,j) = dPhiHydX(i,j)
208 jmc 1.5 & +factorP*(alphRho(i-1,j)+alphRho(i,j))
209 jmc 1.3 & *(rStarFacC(i,j,bi,bj)-rStarFacC(i-1,j,bi,bj))
210 jmc 1.9 & *rC(k)*recip_dxC(i,j,bi,bj)*recip_deepFacC(k)
211 jmc 1.10 ENDDO
212     ENDDO
213     DO j=jMin+1,jMax
214     DO i=iMin,iMax
215 jmc 1.3 dPhiHydY(i,j) = dPhiHydY(i,j)
216 jmc 1.5 & +factorP*(alphRho(i,j-1)+alphRho(i,j))
217 jmc 1.3 & *(rStarFacC(i,j,bi,bj)-rStarFacC(i,j-1,bi,bj))
218 jmc 1.9 & *rC(k)*recip_dyC(i,j,bi,bj)*recip_deepFacC(k)
219 jmc 1.3 ENDDO
220     ENDDO
221     ENDIF
222     ENDIF
223     #endif /* NONLIN_FRSURF */
224 jmc 1.1
225     #endif /* INCLUDE_PHIHYD_CALCULATION_CODE */
226    
227     RETURN
228     END

  ViewVC Help
Powered by ViewVC 1.1.22