53 |
_RL atm_cp, atm_kappa, atm_po |
_RL atm_cp, atm_kappa, atm_po |
54 |
|
|
55 |
IF ( buoyancyRelation .eq. 'OCEANIC' ) THEN |
IF ( buoyancyRelation .eq. 'OCEANIC' ) THEN |
56 |
|
C This is the hydrostatic pressure calculation for the Ocean |
57 |
|
C which uses the FIND_RHO() routine to calculate density |
58 |
|
C before integrating g*rho over the current layer/interface |
59 |
|
|
60 |
dRloc=drC(k) |
dRloc=drC(k) |
61 |
IF (k.EQ.1) dRloc=drF(1) |
IF (k.EQ.1) dRloc=drF(1) |
85 |
C Hydrostatic pressure at cell centers |
C Hydrostatic pressure at cell centers |
86 |
DO j=jMin,jMax |
DO j=jMin,jMax |
87 |
DO i=iMin,iMax |
DO i=iMin,iMax |
88 |
#ifdef ALLOW_AUTODIFF_TAMC |
#ifdef ALLOW_AUTODIFF_TAMC |
89 |
Is this directive correct or even necessary in this new code? |
Is this directive correct or even necessary in this new code? |
90 |
CADJ GENERAL |
CADJ GENERAL |
91 |
#endif |
#endif /* ALLOW_AUTODIFF_TAMC */ |
92 |
C This discretization is the "finite volume" form |
|
93 |
|
C---------- This discretization is the "finite volume" form |
94 |
C which has not been used to date since it does not |
C which has not been used to date since it does not |
95 |
C conserve KE+PE exactly even though it is more natural |
C conserve KE+PE exactly even though it is more natural |
96 |
C |
C |
97 |
C phiHyd(i,j,k)=phiHydInterface(i,j)+ |
c phiHyd(i,j,k)=phiHydInterface(i,j)+ |
98 |
C & 0.5*drF(K)*gravity*alphaRho(i,j) |
c & 0.5*drF(K)*gravity*alphaRho(i,j) |
99 |
C |
c phiHydInterface(i,j)=phiHydInterface(i,j)+ |
100 |
C This discretization is the "energy conserving" form |
c & drF(K)*gravity*alphaRho(i,j) |
101 |
|
C----------------------------------------------------------------------- |
102 |
|
|
103 |
|
C---------- This discretization is the "energy conserving" form |
104 |
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 |
105 |
C |
C |
106 |
phiHyd(i,j,k)=phiHyd(i,j,k)+ |
phiHyd(i,j,k)=phiHyd(i,j,k)+ |
107 |
& 0.5*dRloc*gravity*alphaRho(i,j) |
& 0.5*dRloc*gravity*alphaRho(i,j) |
108 |
phiHyd(i,j,k+1)=phiHyd(i,j,k)+ |
phiHyd(i,j,k+1)=phiHyd(i,j,k)+ |
109 |
& 0.5*dRlocKp1*gravity*alphaRho(i,j) |
& 0.5*dRlocKp1*gravity*alphaRho(i,j) |
110 |
|
C----------------------------------------------------------------------- |
111 |
ENDDO |
ENDDO |
112 |
ENDDO |
ENDDO |
113 |
|
|
114 |
C Hydrostatic pressure at interface below |
|
115 |
DO j=jMin,jMax |
|
|
DO i=iMin,iMax |
|
|
phiHydInterface(i,j)=phiHydInterface(i,j)+ |
|
|
& drF(K)*gravity*alphaRho(i,j) |
|
|
ENDDO |
|
|
ENDDO |
|
116 |
|
|
117 |
ELSEIF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN |
ELSEIF ( buoyancyRelation .eq. 'ATMOSPHERIC' ) THEN |
118 |
|
C This is the hydrostatic geopotential calculation for the Atmosphere |
119 |
|
C The ideal gas law is used implicitly here rather than calculating |
120 |
|
C the specific volume, analogous to the oceanic case. |
121 |
|
|
122 |
|
C Integrate d Phi / d pi |
123 |
|
|
124 |
|
C *NOTE* These constants should be in the data file and PARAMS.h |
125 |
atm_cp=1004. _d 0 |
atm_cp=1004. _d 0 |
126 |
atm_kappa=2. _d 0/7. _d 0 |
atm_kappa=2. _d 0/7. _d 0 |
127 |
atm_po=1. _d 5 |
atm_po=1. _d 5 |
128 |
IF (K.EQ.1) THEN |
IF (K.EQ.1) THEN |
|
C Integrate d Phi / d R |
|
129 |
ddRp1=atm_cp*( ((rC(K)/atm_po)**atm_kappa) |
ddRp1=atm_cp*( ((rC(K)/atm_po)**atm_kappa) |
130 |
& -((rF(K)/atm_po)**atm_kappa) ) |
& -((rF(K)/atm_po)**atm_kappa) ) |
131 |
DO j=jMin,jMax |
DO j=jMin,jMax |
132 |
DO i=iMin,iMax |
DO i=iMin,iMax |
133 |
ddRp=ddRp1 |
ddRp=ddRp1 |
134 |
IF (hFacC(I,J, K ,bi,bj).EQ.0.) ddRp=0. |
IF (hFacC(I,J, K ,bi,bj).EQ.0.) ddRp=0. |
135 |
phiHyd(i,j,K)=0. |
C------------ The integration for the first level phi(k=1) is the |
136 |
|
C same for both the "finite volume" and energy conserving |
137 |
|
C methods. |
138 |
|
C *NOTE* The geopotential boundary condition should go |
139 |
|
C here but has not been implemented yet |
140 |
|
phiHyd(i,j,K)=0. |
141 |
& -ddRp*(theta(I,J,K,bi,bj)-tRef(K)) |
& -ddRp*(theta(I,J,K,bi,bj)-tRef(K)) |
142 |
|
C----------------------------------------------------------------------- |
143 |
ENDDO |
ENDDO |
144 |
ENDDO |
ENDDO |
145 |
ELSE |
ELSE |
146 |
C Integrate d Phi / d R |
|
147 |
|
C-------- This discretization is the "finite volume" form which |
148 |
|
C integrates the hydrostatic equation of each half/sub-layer. |
149 |
|
C This seems most natural and could easily allow for lopped cells |
150 |
|
C by replacing rF(K) with the height of the surface (not implemented). |
151 |
|
C in the lower layers (e.g. at k=1). |
152 |
|
C |
153 |
|
c ddRm1=atm_cp*( ((rF( K )/atm_po)**atm_kappa) |
154 |
|
c & -((rC(K-1)/atm_po)**atm_kappa) ) |
155 |
|
c ddRp1=atm_cp*( ((rC( K )/atm_po)**atm_kappa) |
156 |
|
c & -((rF( K )/atm_po)**atm_kappa) ) |
157 |
|
C----------------------------------------------------------------------- |
158 |
|
|
159 |
|
|
160 |
|
C-------- This discretization is the energy conserving form |
161 |
ddRp1=atm_cp*( ((rC( K )/atm_po)**atm_kappa) |
ddRp1=atm_cp*( ((rC( K )/atm_po)**atm_kappa) |
162 |
& -((rF( K )/atm_po)**atm_kappa) ) |
& -((rC(K-1)/atm_po)**atm_kappa) )*0.5 |
163 |
ddRm1=atm_cp*( ((rF( K )/atm_po)**atm_kappa) |
ddRm1=ddRp1 |
164 |
& -((rC(K-1)/atm_po)**atm_kappa) ) |
C----------------------------------------------------------------------- |
165 |
|
|
166 |
DO j=jMin,jMax |
DO j=jMin,jMax |
167 |
DO i=iMin,iMax |
DO i=iMin,iMax |
168 |
ddRp=ddRp1 |
ddRp=ddRp1 |
169 |
ddRm=ddRm1 |
ddRm=ddRm1 |
170 |
IF (hFacC(I,J, K ,bi,bj).EQ.0.) ddRp=0. |
IF (hFacC(I,J, K ,bi,bj).EQ.0.) ddRp=0. |
171 |
IF (hFacC(I,J,K-1,bi,bj).EQ.0.) ddRm=0. |
IF (hFacC(I,J,K-1,bi,bj).EQ.0.) ddRm=0. |
172 |
phiHyd(i,j,K)=phiHyd(i,j,K-1) |
phiHyd(i,j,K)=phiHyd(i,j,K-1) |
173 |
& -( ddRm*(theta(I,J,K-1,bi,bj)-tRef(K-1)) |
& -( ddRm*(theta(I,J,K-1,bi,bj)-tRef(K-1)) |
174 |
& +ddRp*(theta(I,J, K ,bi,bj)-tRef( K )) ) |
& +ddRp*(theta(I,J, K ,bi,bj)-tRef( K )) ) |
175 |
ENDDO |
C Old code bug looked like this |
176 |
|
Cold phiHyd(i,j,K)=phiHyd(i,j,K-1)-(ddRm1* |
177 |
|
Cold & (theta(I,J,K-1,bi,bj)+theta(I,J,K,bi,bj))-tRef(K)) |
178 |
|
ENDDO |
179 |
ENDDO |
ENDDO |
180 |
ENDIF |
ENDIF |
181 |
|
|
182 |
|
|
183 |
ELSE |
ELSE |
184 |
STOP 'CALC_PHI_HYD: We should never reach this point!' |
STOP 'CALC_PHI_HYD: We should never reach this point!' |
185 |
ENDIF |
ENDIF |