50 |
#include "tamc.h" |
#include "tamc.h" |
51 |
#include "tamc_keys.h" |
#include "tamc_keys.h" |
52 |
#endif /* ALLOW_AUTODIFF_TAMC */ |
#endif /* ALLOW_AUTODIFF_TAMC */ |
53 |
|
#include "SURFACE.h" |
54 |
|
|
55 |
C !INPUT/OUTPUT PARAMETERS: |
C !INPUT/OUTPUT PARAMETERS: |
56 |
C == Routine arguments == |
C == Routine arguments == |
67 |
INTEGER i,j, Kp1 |
INTEGER i,j, Kp1 |
68 |
_RL zero, one, half |
_RL zero, one, half |
69 |
_RL alphaRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
_RL alphaRho(1-OLx:sNx+OLx,1-OLy:sNy+OLy) |
70 |
_RL dRloc,dRlocKp1 |
_RL dRloc,dRlocKp1,locAlpha |
71 |
_RL ddPI, ddPIm, ddPIp, ratioRp, ratioRm |
_RL ddPI, ddPIm, ddPIp, ratioRp, ratioRm |
72 |
CEOP |
CEOP |
73 |
|
|
171 |
ENDDO |
ENDDO |
172 |
ENDDO |
ENDDO |
173 |
|
|
174 |
|
ELSEIF ( buoyancyRelation .eq. 'OCEANICP' ) THEN |
175 |
|
C This is the hydrostatic pressure calculation for the Ocean |
176 |
|
C which uses the FIND_RHO() routine to calculate density |
177 |
|
C before integrating g*rho over the current layer/interface |
178 |
|
|
179 |
|
dRloc=drC(k) |
180 |
|
IF (k.EQ.1) dRloc=drF(1) |
181 |
|
IF (k.EQ.Nr) THEN |
182 |
|
dRlocKp1=0. |
183 |
|
ELSE |
184 |
|
dRlocKp1=drC(k+1) |
185 |
|
ENDIF |
186 |
|
|
187 |
|
IF (k.EQ.1) THEN |
188 |
|
DO j=jMin,jMax |
189 |
|
DO i=iMin,iMax |
190 |
|
phiHyd(i,j,k)=0. |
191 |
|
phiHyd(i,j,k)=pload(i,j,bi,bj) |
192 |
|
c & -Ro_surf(i,j,bi,bj)*recip_rhoNil |
193 |
|
c & -(Ro_surf(i,j,bi,bj)-.5*drF( kSurfC(i,j,bi,bj) ))/1000. |
194 |
|
c & -(Ro_surf(i,j,bi,bj)-.5*drF( kSurfC(i,j,bi,bj) ))*recip_rhoNil |
195 |
|
ENDDO |
196 |
|
ENDDO |
197 |
|
ENDIF |
198 |
|
|
199 |
|
C Calculate density |
200 |
|
#ifdef ALLOW_AUTODIFF_TAMC |
201 |
|
kkey = (ikey-1)*Nr + k |
202 |
|
CADJ STORE theta(:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte |
203 |
|
CADJ STORE salt (:,:,k,bi,bj) = comlev1_bibj_k, key=kkey, byte=isbyte |
204 |
|
#endif /* ALLOW_AUTODIFF_TAMC */ |
205 |
|
CALL FIND_RHO( bi, bj, iMin, iMax, jMin, jMax, k, k, eosType, |
206 |
|
& theta, salt, |
207 |
|
& alphaRho, myThid) |
208 |
|
|
209 |
|
C Hydrostatic pressure at cell centers |
210 |
|
DO j=jMin,jMax |
211 |
|
DO i=iMin,iMax |
212 |
|
locAlpha=alphaRho(i,j)+rhoNil |
213 |
|
IF (locAlpha.NE.0.) locAlpha=maskC(i,j,k,bi,bj)/locAlpha |
214 |
|
|
215 |
|
C---------- This discretization is the "finite volume" form |
216 |
|
C which has not been used to date since it does not |
217 |
|
C conserve KE+PE exactly even though it is more natural |
218 |
|
C |
219 |
|
c IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)+ |
220 |
|
c & drF(K)*locAlpha |
221 |
|
c phiHyd(i,j,k)=phiHyd(i,j,k)+ |
222 |
|
c & 0.5*drF(K)*locAlpha |
223 |
|
C----------------------------------------------------------------------- |
224 |
|
|
225 |
|
C---------- This discretization is the "energy conserving" form |
226 |
|
C which has been used since at least Adcroft et al., MWR 1997 |
227 |
|
C |
228 |
|
phiHyd(i,j,k)=phiHyd(i,j,k)+ |
229 |
|
& 0.5*dRloc*locAlpha |
230 |
|
IF (k.LT.Nr) phiHyd(i,j,k+1)=phiHyd(i,j,k)+ |
231 |
|
& 0.5*dRlocKp1*locAlpha |
232 |
|
C----------------------------------------------------------------------- |
233 |
|
ENDDO |
234 |
|
ENDDO |
235 |
|
|
236 |
ELSEIF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN |
ELSEIF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN |
237 |
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |
C---+----1----+----2----+----3----+----4----+----5----+----6----+----7-|--+----| |