1 |
C $Header: /u/gcmpack/MITgcm/pkg/bbl/bbl_calc_rhs.F,v 1.1 2011/08/06 03:13:22 dimitri Exp $ |
2 |
C $Name: $ |
3 |
|
4 |
#include "BBL_OPTIONS.h" |
5 |
|
6 |
CBOP |
7 |
C !ROUTINE: BBL_CALC_RHS |
8 |
|
9 |
C !INTERFACE: |
10 |
SUBROUTINE BBL_CALC_RHS( |
11 |
I myTime, myIter, myThid ) |
12 |
|
13 |
C !DESCRIPTION: |
14 |
C Calculate tendency of tracers due to bottom boundary layer. |
15 |
C Scheme is currently coded for dTtracerLev(k) .EQ. deltaT. |
16 |
|
17 |
C !USES: |
18 |
IMPLICIT NONE |
19 |
#include "SIZE.h" |
20 |
#include "EEPARAMS.h" |
21 |
#include "PARAMS.h" |
22 |
#include "GRID.h" |
23 |
#include "DYNVARS.h" |
24 |
#include "BBL.h" |
25 |
|
26 |
C !INPUT PARAMETERS: |
27 |
C myTime :: Current time in simulation |
28 |
C myIter :: Current time-step number |
29 |
C myThid :: my Thread Id number |
30 |
_RL myTime |
31 |
INTEGER myIter, myThid |
32 |
|
33 |
C !OUTPUT PARAMETERS: |
34 |
|
35 |
C !LOCAL VARIABLES: |
36 |
C bi,bj :: Tile indices |
37 |
C i,j :: Loop indices |
38 |
C kBot :: k index of bottommost wet grid box |
39 |
C kLowC1 :: k index of bottommost (i,j) cell |
40 |
C kLowC2 :: k index of bottommost (i+1,j) or (i,j+1) cell |
41 |
C kl :: k index at which to compare 2 cells |
42 |
C resThk :: thickness of bottommost wet grid box minus bbl_eta |
43 |
C resTheta :: temperature of this residual volume |
44 |
C resSalt :: salinity of this residual volume |
45 |
C deltaRho :: density change |
46 |
C deltaDpt :: depth change |
47 |
C bbl_tend :: temporary variable for tendency terms |
48 |
C sloc :: salinity of bottommost wet grid box |
49 |
C tloc :: temperature of bottommost wet grid box |
50 |
C rholoc :: in situ density of bottommost wet grid box |
51 |
C rhoBBL :: in situ density of bottom boundary layer |
52 |
C fZon :: Zonal flux |
53 |
C fMer :: Meridional flux |
54 |
C bbl_rho1 :: local (i,j) density |
55 |
C bbl_rho2 :: local (i+1, j) or (i,j+1) density |
56 |
INTEGER bi, bj |
57 |
INTEGER i, j, kBot, kLowC1, kLowC2, kl |
58 |
_RL resThk, resTheta, resSalt |
59 |
_RL deltaRho, deltaDpt, bbl_tend |
60 |
_RL bbl_rho1, bbl_rho2 |
61 |
_RL sloc ( 1-Olx:sNx+Olx, 1-Oly:sNy+Oly ) |
62 |
_RL tloc ( 1-Olx:sNx+Olx, 1-Oly:sNy+Oly ) |
63 |
_RL rholoc ( 1-Olx:sNx+Olx, 1-Oly:sNy+Oly ) |
64 |
_RL rhoBBL ( 1-Olx:sNx+Olx, 1-Oly:sNy+Oly ) |
65 |
_RL fZon ( 1-Olx:sNx+Olx, 1-Oly:sNy+Oly ) |
66 |
_RL fMer ( 1-Olx:sNx+Olx, 1-Oly:sNy+Oly ) |
67 |
CHARACTER*(MAX_LEN_MBUF) msgBuf |
68 |
CEOP |
69 |
|
70 |
C-- Loops on tile indices bi,bj |
71 |
DO bj=myByLo(myThid),myByHi(myThid) |
72 |
DO bi=myBxLo(myThid),myBxHi(myThid) |
73 |
|
74 |
C Initialize tendency terms and make local copy of |
75 |
C bottomost temperature, salinity, and in-situ desnity |
76 |
C and of in-situ BBL density |
77 |
DO j=1-OLy,sNy+OLy |
78 |
DO i=1-OLx,sNx+OLx |
79 |
bbl_TendTheta(i,j,bi,bj) = 0. _d 0 |
80 |
bbl_TendSalt (i,j,bi,bj) = 0. _d 0 |
81 |
kBot = max(1,kLowC(i,j,bi,bj)) |
82 |
tLoc(i,j) = theta(i,j,kBot,bi,bj) |
83 |
sLoc(i,j) = salt (i,j,kBot,bi,bj) |
84 |
rholoc(i,j) = rhoInSitu(i,j,kBot,bi,bj) |
85 |
IF ( kBot .EQ. Nr ) THEN |
86 |
rhoBBL(i,j) = bbl_rho_nr(i,j,bi,bj) |
87 |
ELSE |
88 |
rhoBBL(i,j) = rhoInSitu(i,j,kBot+1, bi,bj) |
89 |
ENDIF |
90 |
ENDDO |
91 |
ENDDO |
92 |
|
93 |
C==== Compute and apply vertical exchange between BBL and |
94 |
C residual volume of botommost wet grid box. |
95 |
C This operation does not change total tracer quantity |
96 |
C in botommost wet grid box. |
97 |
|
98 |
DO j=1-Oly,sNy+Oly |
99 |
DO i=1-Olx,sNx+Olx |
100 |
kBot = kLowC(i,j,bi,bj) |
101 |
IF ( kBot .GT. 0 ) THEN |
102 |
resThk = hFacC(i,j,kBot,bi,bj)*drF(kBot) - bbl_eta(i,j,bi,bj) |
103 |
|
104 |
C If bbl occupies the complete bottom model grid box or |
105 |
C if model density is higher than BBL then mix instantly. |
106 |
IF ( (resThk.LE.0) .OR. (rhoLoc(i,j).GE.rhoBBL(i,j)) ) THEN |
107 |
bbl_theta(i,j,bi,bj) = tLoc(i,j) |
108 |
bbl_salt (i,j,bi,bj) = sLoc(i,j) |
109 |
|
110 |
C If model density is lower than BBL, slowly diffuse upward. |
111 |
ELSE |
112 |
resTheta = ( tLoc(i,j) * (resThk+bbl_eta(i,j,bi,bj)) - |
113 |
& (bbl_theta(i,j,bi,bj)*bbl_eta(i,j,bi,bj)) ) / resThk |
114 |
resSalt = ( sLoc(i,j) * (resThk+bbl_eta(i,j,bi,bj)) - |
115 |
& (bbl_salt (i,j,bi,bj)*bbl_eta(i,j,bi,bj)) ) / resThk |
116 |
bbl_theta(i,j,bi,bj) = bbl_theta(i,j,bi,bj) + |
117 |
& deltaT * (resTheta-bbl_theta(i,j,bi,bj)) / bbl_RelaxR |
118 |
bbl_salt (i,j,bi,bj) = bbl_salt (i,j,bi,bj) + |
119 |
& deltaT * (resSalt -bbl_salt (i,j,bi,bj)) / bbl_RelaxR |
120 |
ENDIF |
121 |
ENDIF |
122 |
ENDDO |
123 |
ENDDO |
124 |
|
125 |
C==== Compute zonal bbl exchange. |
126 |
DO j=1-Oly,sNy+Oly |
127 |
DO i=1-Olx,sNx+Olx-1 |
128 |
kLowC1 = kLowC(i,j,bi,bj) |
129 |
kLowC2 = kLowC(i+1,j,bi,bj) |
130 |
IF ((kLowC1.GT.0).AND.(kLowC2.GT.0)) THEN |
131 |
C Compare the bbl densities at the higher pressure |
132 |
C (highest possible density of given t,s) |
133 |
C bbl in situ density is stored in kLowC + 1 index |
134 |
kl = MAX(kLowC1, kLowC2) + 1 |
135 |
IF ( kl .GT. Nr ) THEN |
136 |
bbl_rho1 = bbl_rho_nr(i,j,bi,bj) |
137 |
bbl_rho2 = bbl_rho_nr(i+1,j,bi,bj) |
138 |
else |
139 |
bbl_rho1 = rhoInSitu(i,j,kl,bi,bj) |
140 |
bbl_rho2 = rhoInSitu(i+1,j,kl,bi,bj) |
141 |
endif |
142 |
deltaRho = bbl_rho2 - bbl_rho1 |
143 |
deltaDpt = R_low(i ,j,bi,bj) + bbl_eta(i ,j,bi,bj) - |
144 |
& R_low(i+1,j,bi,bj) - bbl_eta(i+1,j,bi,bj) |
145 |
|
146 |
C If heavy BBL water is higher than light BBL water, |
147 |
C exchange properties laterally. |
148 |
IF ( (deltaRho*deltaDpt) .LE. 0. ) THEN |
149 |
bbl_TendTheta(i,j,bi,bj) = bbl_TendTheta(i,j,bi,bj) + |
150 |
& ( bbl_theta(i+1,j,bi,bj) - bbl_theta(i,j,bi,bj) ) / |
151 |
& bbl_RelaxH |
152 |
bbl_TendTheta(i+1,j,bi,bj) = bbl_TendTheta(i+1,j,bi,bj) - |
153 |
& ( bbl_theta(i+1,j,bi,bj) - bbl_theta(i,j,bi,bj) ) * |
154 |
& ( rA(i ,j,bi,bj) * bbl_eta(i ,j,bi,bj) ) / |
155 |
& ( rA(i+1,j,bi,bj) * bbl_eta(i+1,j,bi,bj) * bbl_RelaxH ) |
156 |
bbl_TendSalt(i,j,bi,bj) = bbl_TendSalt(i,j,bi,bj) + |
157 |
& ( bbl_salt(i+1,j,bi,bj) - bbl_salt(i,j,bi,bj) ) / |
158 |
& bbl_RelaxH |
159 |
bbl_TendSalt(i+1,j,bi,bj) = bbl_TendSalt(i+1,j,bi,bj) - |
160 |
& ( bbl_salt(i+1,j,bi,bj) - bbl_salt(i,j,bi,bj) ) * |
161 |
& ( rA(i ,j,bi,bj) * bbl_eta(i ,j,bi,bj) ) / |
162 |
& ( rA(i+1,j,bi,bj) * bbl_eta(i+1,j,bi,bj) * bbl_RelaxH ) |
163 |
ENDIF |
164 |
ENDIF |
165 |
ENDDO |
166 |
ENDDO |
167 |
|
168 |
C==== Compute meridional bbl exchange. |
169 |
DO j=1-Oly,sNy+Oly-1 |
170 |
DO i=1-Olx,sNx+Olx |
171 |
kLowC1 = kLowC(i,j,bi,bj) |
172 |
kLowC2 = kLowC(i,j+1, bi,bj) |
173 |
IF ((kLowC1.GT.0).AND.(kLowC2.GT.0)) THEN |
174 |
C compare the bbl densities at the higher pressure |
175 |
C (highest possible density of given t,s) |
176 |
C bbl in situ density is stored in kLowC + 1 index |
177 |
kl = MAX(kLowC1, kLowC2) + 1 |
178 |
IF ( kl .GT. Nr ) THEN |
179 |
bbl_rho1 = bbl_rho_nr(i,j,bi,bj) |
180 |
bbl_rho2 = bbl_rho_nr(i,j+1,bi,bj) |
181 |
else |
182 |
bbl_rho1 = rhoInSitu(i,j,kl,bi,bj) |
183 |
bbl_rho2 = rhoInSitu(i,j+1,kl,bi,bj) |
184 |
endif |
185 |
deltaRho = bbl_rho2 - bbl_rho1 |
186 |
deltaDpt = R_low(i,j ,bi,bj) + bbl_eta(i,j ,bi,bj) - |
187 |
& R_low(i,j+1,bi,bj) - bbl_eta(i,j+1,bi,bj) |
188 |
|
189 |
C If heavy BBL water is higher than light BBL water, |
190 |
C exchange properties laterally. |
191 |
IF ( (deltaRho*deltaDpt) .LE. 0. ) THEN |
192 |
bbl_TendTheta(i,j,bi,bj) = bbl_TendTheta(i,j,bi,bj) + |
193 |
& ( bbl_theta(i,j+1,bi,bj) - bbl_theta(i,j,bi,bj) ) / |
194 |
& bbl_RelaxH |
195 |
bbl_TendTheta(i,j+1,bi,bj) = bbl_TendTheta(i,j+1,bi,bj) - |
196 |
& ( bbl_theta(i,j+1,bi,bj) - bbl_theta(i,j,bi,bj) ) * |
197 |
& ( rA(i ,j,bi,bj) * bbl_eta(i ,j,bi,bj) ) / |
198 |
& ( rA(i,j+1,bi,bj) * bbl_eta(i,j+1,bi,bj) ) / |
199 |
& bbl_RelaxH |
200 |
bbl_TendSalt(i,j,bi,bj) = bbl_TendSalt(i,j,bi,bj) + |
201 |
& ( bbl_salt(i,j+1,bi,bj) - bbl_salt(i,j,bi,bj) ) / |
202 |
& bbl_RelaxH |
203 |
bbl_TendSalt(i,j+1,bi,bj) = bbl_TendSalt(i,j+1,bi,bj) - |
204 |
& ( bbl_salt(i,j+1,bi,bj)-bbl_salt(i,j,bi,bj)) * |
205 |
& ( rA(i ,j,bi,bj) * bbl_eta(i ,j,bi,bj) ) / |
206 |
& ( rA(i,j+1,bi,bj) * bbl_eta(i,j+1,bi,bj) * bbl_RelaxH ) |
207 |
ENDIF |
208 |
ENDIF |
209 |
ENDDO |
210 |
ENDDO |
211 |
|
212 |
C==== Apply lateral BBL exchange then scale tendency term |
213 |
C for botommost wet grid box. |
214 |
DO j=1-Oly,sNy+Oly-1 |
215 |
DO i=1-Olx,sNx+Olx-1 |
216 |
kBot = kLowC(i,j,bi,bj) |
217 |
IF ( kBot .GT. 0 ) THEN |
218 |
bbl_theta(i,j,bi,bj) = bbl_theta(i,j,bi,bj) + |
219 |
& deltaT * bbl_TendTheta(i,j,bi,bj) |
220 |
bbl_salt (i,j,bi,bj) = bbl_salt (i,j,bi,bj) + |
221 |
& deltaT * bbl_TendSalt (i,j,bi,bj) |
222 |
bbl_TendTheta(i,j,bi,bj) = bbl_TendTheta(i,j,bi,bj) * |
223 |
& bbl_eta(i,j,bi,bj) / (hFacC(i,j,kBot,bi,bj)*drF(kBot)) |
224 |
bbl_TendSalt (i,j,bi,bj) = bbl_TendSalt (i,j,bi,bj) * |
225 |
& bbl_eta(i,j,bi,bj) / (hFacC(i,j,kBot,bi,bj)*drF(kBot)) |
226 |
ENDIF |
227 |
ENDDO |
228 |
ENDDO |
229 |
|
230 |
#ifdef ALLOW_DEBUG |
231 |
IF ( debugLevel .GE. debLevB ) THEN |
232 |
C Check salinity conservation |
233 |
bbl_tend=0 |
234 |
DO j=1,sNy |
235 |
DO i=1,sNx |
236 |
kBot = kLowC(i,j,bi,bj) |
237 |
IF ( kBot .GT. 0 ) THEN |
238 |
bbl_tend = bbl_tend + bbl_TendSalt(i,j,bi,bj) * |
239 |
& hFacC(i,j,kBot,bi,bj) * drF(kBot) *rA(i,j,bi,bj) |
240 |
ENDIF |
241 |
ENDDO |
242 |
ENDDO |
243 |
_GLOBAL_SUM_RL( bbl_tend, myThid ) |
244 |
WRITE(msgBuf,'(A,E10.2)') 'total salt tendency = ', bbl_tend |
245 |
CALL PRINT_MESSAGE( msgBuf, standardMessageUnit, |
246 |
& SQUEEZE_RIGHT, myThid ) |
247 |
ENDIF |
248 |
#endif /* ALLOW_DEBUG */ |
249 |
|
250 |
C-- end bi,bj loops. |
251 |
ENDDO |
252 |
ENDDO |
253 |
|
254 |
CALL EXCH_XY_RL( bbl_theta, myThid ) |
255 |
CALL EXCH_XY_RL( bbl_salt , myThid ) |
256 |
|
257 |
RETURN |
258 |
END |