76 |
_RL ddPI, ddPIm, ddPIp, ratioRp, ratioRm |
_RL ddPI, ddPIm, ddPIp, ratioRp, ratioRm |
77 |
INTEGER iMnLoc,jMnLoc |
INTEGER iMnLoc,jMnLoc |
78 |
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 ) |
79 |
|
LOGICAL useDiagPhiRlow |
80 |
CEOP |
CEOP |
81 |
|
useDiagPhiRlow = .TRUE. |
82 |
|
|
83 |
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
84 |
C Atmosphere: |
C Atmosphere: |
151 |
CALL QUASIHYDROSTATICTERMS(bi,bj,k,alphaRho,myThid) |
CALL QUASIHYDROSTATICTERMS(bi,bj,k,alphaRho,myThid) |
152 |
ENDIF |
ENDIF |
153 |
|
|
154 |
|
C--- Diagnose Hydrostatic pressure at the bottom: |
155 |
|
IF (useDiagPhiRlow) THEN |
156 |
|
CALL DIAGS_PHI_RLOW( |
157 |
|
I k, bi, bj, iMin,iMax, jMin,jMax, |
158 |
|
I phiHyd, alphaRho, tFld, sFld, |
159 |
|
I myTime, myIter, myThid) |
160 |
|
ENDIF |
161 |
|
|
162 |
C--- Hydrostatic pressure at cell centers |
C--- Hydrostatic pressure at cell centers |
163 |
|
|
164 |
IF (integr_GeoPot.EQ.1) THEN |
IF (integr_GeoPot.EQ.1) THEN |
171 |
C which has not been used to date since it does not |
C which has not been used to date since it does not |
172 |
C conserve KE+PE exactly even though it is more natural |
C conserve KE+PE exactly even though it is more natural |
173 |
C |
C |
|
IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN |
|
|
phiHydLow(i,j,bi,bj) = phiHyd(i,j,k) |
|
|
& + hFacC(i,j,k,bi,bj) |
|
|
& *drF(K)*gravity*alphaRho(i,j)*recip_rhoConst |
|
|
& + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj) |
|
|
ENDIF |
|
174 |
IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k) |
IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k) |
175 |
& + drF(K)*gravity*alphaRho(i,j)*recip_rhoConst |
& + drF(K)*gravity*alphaRho(i,j)*recip_rhoConst |
176 |
phiHyd(i,j,k)=phiHyd(i,j,k)+ |
phiHyd(i,j,k)=phiHyd(i,j,k)+ |
188 |
C---------- This discretization is the "energy conserving" form |
C---------- This discretization is the "energy conserving" form |
189 |
C which has been used since at least Adcroft et al., MWR 1997 |
C which has been used since at least Adcroft et al., MWR 1997 |
190 |
C |
C |
|
|
|
191 |
phiHyd(i,j,k)=phiHyd(i,j,k) |
phiHyd(i,j,k)=phiHyd(i,j,k) |
192 |
& +half*dRloc*gravity*alphaRho(i,j)*recip_rhoConst |
& +half*dRloc*gravity*alphaRho(i,j)*recip_rhoConst |
193 |
IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k) |
IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k) |
194 |
& +half*dRlocKp1*gravity*alphaRho(i,j)*recip_rhoConst |
& +half*dRlocKp1*gravity*alphaRho(i,j)*recip_rhoConst |
195 |
|
|
|
C---------- Compute bottom pressure deviation from gravity*rho0*H |
|
|
C This has to be done starting from phiHyd at the current |
|
|
C tracer point and .5 of the cell's thickness has to be |
|
|
C substracted from hFacC |
|
|
IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN |
|
|
phiHydLow(i,j,bi,bj) = phiHyd(i,j,k) |
|
|
& + (hFacC(i,j,k,bi,bj)-half)*drF(K) |
|
|
& *gravity*alphaRho(i,j)*recip_rhoConst |
|
|
& + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj) |
|
|
ENDIF |
|
|
|
|
196 |
ENDDO |
ENDDO |
197 |
ENDDO |
ENDDO |
198 |
|
|
225 |
ENDDO |
ENDDO |
226 |
ENDIF |
ENDIF |
227 |
|
|
228 |
C Calculate density |
C-- Calculate density |
229 |
#ifdef ALLOW_AUTODIFF_TAMC |
#ifdef ALLOW_AUTODIFF_TAMC |
230 |
kkey = (ikey-1)*Nr + k |
kkey = (ikey-1)*Nr + k |
231 |
CADJ STORE tFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte |
CADJ STORE tFld (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte |
238 |
CADJ STORE alphaRho (:,:) = comlev1_bibj_k, key=kkey, byte=isbyte |
CADJ STORE alphaRho (:,:) = comlev1_bibj_k, key=kkey, byte=isbyte |
239 |
#endif /* ALLOW_AUTODIFF_TAMC */ |
#endif /* ALLOW_AUTODIFF_TAMC */ |
240 |
|
|
241 |
|
C-- Calculate specific volume anomaly : alpha' = 1/rho - alpha_Cst |
242 |
|
DO j=jMin,jMax |
243 |
|
DO i=iMin,iMax |
244 |
|
locAlpha=alphaRho(i,j)+rhoConst |
245 |
|
alphaRho(i,j)=maskC(i,j,k,bi,bj)* |
246 |
|
& (one/locAlpha - recip_rhoConst) |
247 |
|
ENDDO |
248 |
|
ENDDO |
249 |
|
|
250 |
|
C--- Diagnose Sea-surface height (Hydrostatic geopotential at r=Rlow): |
251 |
|
IF (useDiagPhiRlow) THEN |
252 |
|
CALL DIAGS_PHI_RLOW( |
253 |
|
I k, bi, bj, iMin,iMax, jMin,jMax, |
254 |
|
I phiHyd, alphaRho, tFld, sFld, |
255 |
|
I myTime, myIter, myThid) |
256 |
|
ENDIF |
257 |
|
|
258 |
C---- Hydrostatic pressure at cell centers |
C---- Hydrostatic pressure at cell centers |
259 |
|
|
262 |
|
|
263 |
DO j=jMin,jMax |
DO j=jMin,jMax |
264 |
DO i=iMin,iMax |
DO i=iMin,iMax |
|
locAlpha=alphaRho(i,j)+rhoConst |
|
|
locAlpha=maskC(i,j,k,bi,bj)* |
|
|
& (one/locAlpha - recip_rhoConst) |
|
|
c IF (locAlpha.NE.0.) locAlpha=maskC(i,j,k,bi,bj)/locAlpha |
|
265 |
|
|
266 |
C---------- This discretization is the "finite volume" form |
C---------- This discretization is the "finite volume" form |
267 |
C which has not been used to date since it does not |
C which has not been used to date since it does not |
268 |
C conserve KE+PE exactly even though it is more natural |
C conserve KE+PE exactly even though it is more natural |
269 |
C |
C |
|
IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN |
|
|
phiHydLow(i,j,bi,bj) = phiHyd(i,j,k) |
|
|
& + hFacC(i,j,k,bi,bj)*drF(K)*locAlpha |
|
|
& + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj) |
|
|
ENDIF |
|
270 |
IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k) |
IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k) |
271 |
& + hFacC(i,j,k,bi,bj)*drF(K)*locAlpha |
& + hFacC(i,j,k,bi,bj)*drF(K)*alphaRho(i,j) |
272 |
phiHyd(i,j,k)=phiHyd(i,j,k) |
phiHyd(i,j,k)=phiHyd(i,j,k) |
273 |
& +(hFacC(i,j,k,bi,bj)-half)*drF(K)*locAlpha |
& +(hFacC(i,j,k,bi,bj)-half)*drF(K)*alphaRho(i,j) |
274 |
|
|
275 |
ENDDO |
ENDDO |
276 |
ENDDO |
ENDDO |
280 |
|
|
281 |
DO j=jMin,jMax |
DO j=jMin,jMax |
282 |
DO i=iMin,iMax |
DO i=iMin,iMax |
|
locAlpha=alphaRho(i,j)+rhoConst |
|
|
locAlpha=maskC(i,j,k,bi,bj)* |
|
|
& (one/locAlpha - recip_rhoConst) |
|
|
c IF (locAlpha.NE.0.) locAlpha=maskC(i,j,k,bi,bj)/locAlpha |
|
283 |
|
|
284 |
C---------- This discretization is the "energy conserving" form |
C---------- This discretization is the "energy conserving" form |
285 |
|
|
286 |
phiHyd(i,j,k)=phiHyd(i,j,k) |
phiHyd(i,j,k)=phiHyd(i,j,k) |
287 |
& + half*dRloc*locAlpha |
& + half*dRloc*alphaRho(i,j) |
288 |
IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k) |
IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k) |
289 |
& + half*dRlocKp1*locAlpha |
& + half*dRlocKp1*alphaRho(i,j) |
|
|
|
|
|
|
|
C---------- Compute gravity*(sea surface elevation) first |
|
|
C This has to be done starting from phiHyd at the current |
|
|
C tracer point and .5 of the cell's thickness has to be |
|
|
C substracted from hFacC |
|
|
IF ( K .EQ. kLowC(i,j,bi,bj) ) THEN |
|
|
phiHydLow(i,j,bi,bj) = phiHyd(i,j,k) |
|
|
& + (hFacC(i,j,k,bi,bj)-half)*drF(k)*locAlpha |
|
|
& + Bo_surf(i,j,bi,bj)*etaN(i,j,bi,bj) |
|
|
ENDIF |
|
290 |
|
|
291 |
ENDDO |
ENDDO |
292 |
ENDDO |
ENDDO |