| 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-|--+----| |
Parent Directory
| 
Revision Log
| 
Revision Graph
| 
Patch