1 |
C $Id: g_calc.F,v 1.3 1998/06/23 13:36:41 adcroft Exp $ |
2 |
#include "CPP_OPTIONS.h" |
3 |
#include "CPP_MACROS.h" |
4 |
|
5 |
#define _DB0 0. |
6 |
#define _DB1 1. |
7 |
#define _wdvdz 1.D0 |
8 |
|
9 |
C=================================================================== |
10 |
C Procedure name: G_CALC | |
11 |
C Function: Caculates the tendency terms for the momentum| |
12 |
C equations. | |
13 |
C Comments: | |
14 |
C=================================================================== |
15 |
CStartofinterface |
16 |
SUBROUTINE G_CALC( |
17 |
I U, V, W, PH, PS, K, |
18 |
U RHS2D, divH, |
19 |
O GU, GV ) |
20 |
IMPLICIT NONE |
21 |
C============== Global data ========================================== |
22 |
#include "SIZE.h" |
23 |
#include "AJAINF.h" |
24 |
#include "OPERATORS.h" |
25 |
#include "GRID.h" |
26 |
#include "PARAMS.h" |
27 |
#include "OLDG.h" |
28 |
#include "FORCING.h" |
29 |
#include "MASKS.h" |
30 |
C============= Routine arguments ===================================== |
31 |
C /--------------------------------------------------------------\ |
32 |
C |U, V, W - X,Y,Z Velocity ( m/s, m/s, Pa/s ). | |
33 |
C |PH - Hydrostatic pressure perturbation (m). | |
34 |
C |PS - Surface pressure (m). | |
35 |
C |GU - d/dt U (m/s/s). | |
36 |
C |GU - d/dt V (m/s/s). | |
37 |
C |RHS2D - Contribution to righ-hand side of | |
38 |
C | pressure eqaution. | |
39 |
C |divH - Contribution to barotropic divergence. | |
40 |
C \--------------------------------------------------------------/ |
41 |
REAL U (nx,ny,nz) |
42 |
REAL V (nx,ny,nz) |
43 |
REAL W (nx,ny,nz) |
44 |
REAL PH(nx,ny,nz) |
45 |
REAL PS(nx,ny ) |
46 |
REAL GU(nx,ny,nz) |
47 |
REAL GV(nx,ny,nz) |
48 |
REAL RHS2D(nx,ny) |
49 |
REAL divH (nx,ny) |
50 |
INTEGER K |
51 |
CEndofinterface |
52 |
C============= Local variables ====================================== |
53 |
C /--------------------------------------------------------------\ |
54 |
C |Slice Notation | |
55 |
C |phiK - XY slice of variable "phi" at level K. | |
56 |
C |phiKP1 - XY slice of variable "phi" at level K+1. | |
57 |
C |phiKM1 - XY slice of variable "phi" at level K-1. | |
58 |
C |phiNM1 - Variable "phi" at time level N-1. | |
59 |
C |phiU - Varaible "phi" interpolated to U grid. | |
60 |
C |phiV - Varaible "phi" interpolated to V grid. | |
61 |
C \--------------------------------------------------------------/ |
62 |
C /--------------------------------------------------------------\ |
63 |
C |Geometry terms | |
64 |
C | o xa,ya,za - X, Y, Z face areas ( m.Pa, m.Pa, m**2 ). | |
65 |
C | o rV - 1/Volume ( m**2.Pa ) | |
66 |
C | o ?Grad?Op - Generic gradient operator A[XYZ][uvwp].d/d[xyz] | |
67 |
C | uGradxOp: AXu/DXu, uGradyOp: AYu/DYu etc... | |
68 |
C \--------------------------------------------------------------/ |
69 |
REAL xaK (nx+1,0:ny) |
70 |
REAL yaK (0:nx,ny+1) |
71 |
REAL uGradxOp (nx ,ny ) |
72 |
REAL uGradyOp (nx ,ny ) |
73 |
REAL uGradzOpK (nx ,ny ) |
74 |
REAL uGradzOpKP1(nx ,ny ) |
75 |
REAL vGradxOp (nx ,ny ) |
76 |
REAL vGradyOp (nx ,ny ) |
77 |
REAL vGradzOpK (nx ,ny ) |
78 |
REAL vGradzOpKP1(nx ,ny ) |
79 |
REAL zaK (0:nx,0:ny) |
80 |
REAL zaKP1 (0:nx,0:ny) |
81 |
REAL rVk (nx ,ny ) |
82 |
C /--------------------------------------------------------------\ |
83 |
C |State variables | |
84 |
C | o u, v, w - X,Y,Z velocity slices. | |
85 |
C | o ph - Hydrostatic pressure. | |
86 |
C \--------------------------------------------------------------/ |
87 |
REAL uK (nx+1,0:ny) |
88 |
REAL uMskK (nx+1,0:ny) |
89 |
REAL uKNM1 (nx+1,0:ny) |
90 |
REAL wK (0:nx,0:ny) |
91 |
REAL wkP1 (0:nx,0:ny) |
92 |
REAL phK (nx ,ny ) |
93 |
REAL vK (0:nx,ny+1) |
94 |
REAL vMskK (nx+1,0:ny) |
95 |
REAL vKNM1 (0:nx,ny+1) |
96 |
REAL vKP1 (nx ,ny ) |
97 |
REAL vKM1 (nx ,ny ) |
98 |
C /--------------------------------------------------------------\ |
99 |
C |Work space arrays | |
100 |
C \--------------------------------------------------------------/ |
101 |
REAL tmp (0:nx+1,0:ny+1) |
102 |
REAL tmp2 (0:nx+1,0:ny+1) |
103 |
REAL stmp |
104 |
C /--------------------------------------------------------------\ |
105 |
C |Arrays for accumulating fluxes | |
106 |
C | o fX, fY, fZ - Flux in X, Y and Z. | |
107 |
C \--------------------------------------------------------------/ |
108 |
REAL fX (0:nx+1,ny) |
109 |
REAL fY (nx,0:ny+1) |
110 |
REAL fZK (nx,ny) |
111 |
REAL fZKP1 (nx,ny) |
112 |
C Loop counters: |
113 |
INTEGER I, J |
114 |
C Flags: |
115 |
LOGICAL TOP_LAYER |
116 |
LOGICAL BOTTOM_LAYER |
117 |
C === Set up controlling flags ======================================== |
118 |
TOP_LAYER = K .EQ. 1 |
119 |
BOTTOM_LAYER = K .EQ. nz |
120 |
C === Set up the slices =============================================== |
121 |
C /---------------------------------------------------------------\ |
122 |
C | xaK <- xa(K) | |
123 |
C | yaK <- ya(K) | |
124 |
C | za[K,KP1] <- za[(K),(K+1)] | |
125 |
C \---------------------------------------------------------------/ |
126 |
xaK (1:nx,1:ny) = xa (_I3(K,:,:)) |
127 |
xaK (nx+1,1:ny) = xaK (1,1:ny ) |
128 |
xaK (:,0 ) = xaK (:,ny ) |
129 |
yaK (1:nx,1:ny) = ya (_I3(K,:,:)) |
130 |
yaK (0,1:ny ) = yaK (nx,1:ny ) |
131 |
yaK (:,ny+1 ) = yaK (:,1 ) |
132 |
zaK (1:nx,1:ny) = za (_I3(K,:,:)) |
133 |
zaK (0,1:ny ) = zaK (nx,1:ny ) |
134 |
zaK (:,0 ) = zaK (:,ny ) |
135 |
IF ( .NOT. BOTTOM_LAYER ) THEN |
136 |
zaKP1(1:nx,1:ny) = za (_I3(K+1,:,:)) |
137 |
zaKP1(0,1:ny ) = zaKP1(nx,1:ny ) |
138 |
zaKP1(:,0 ) = zaKP1(:,ny ) |
139 |
ENDIF |
140 |
C /---------------------------------------------------------------\ |
141 |
C | uGradXOp <- 1/DXu*Bx{xa} | |
142 |
C | uGradYOp <- 1/DYu*Bx{ya} | |
143 |
C | uGradZOp[K,KP1] <- 1/DZu*Bz{ZA[K,K+1]} | |
144 |
C \---------------------------------------------------------------/ |
145 |
DO J = 1, ny |
146 |
DO I = 1, nx |
147 |
uGradxOp(I,J) = 0.5*( xaK(I,J)+xaK(I+1,J) )*uDdxOp(_I3(K,I ,J)) |
148 |
ENDDO |
149 |
ENDDO |
150 |
DO J = 1, ny |
151 |
DO I = 2, nx |
152 |
uGradyOp(I,J) = 0.5*( yaK(I,J)+yaK(I-1,J) )*uDdyOp(_I3(K,I,J)) |
153 |
uGradzOpK(I,J) = |
154 |
& 0.5*(ZA(_I3(K,I,J))+ZA(_I3(K,I-1,J)))*uDdzOp(_I3(K,I,J)) |
155 |
ENDDO |
156 |
uGradyOp(1,J) = 0.5*( yaK(1,J)+yaK(nx,J) )*uDdyOp(_I3(K,1,J)) |
157 |
uGradzOpK(1,J) = |
158 |
& 0.5*(ZA(_I3(K,1,J))+ZA(_I3(K,nx,J)))*uDdzOp(_I3(K,1,J)) |
159 |
ENDDO |
160 |
IF ( .NOT. BOTTOM_LAYER ) THEN |
161 |
DO J = 1, ny |
162 |
DO I = 2, nx |
163 |
uGradzOpKP1(I,J) = |
164 |
& 0.5*(ZA(_I3(K+1,I,J))+ZA(_I3(K+1,I-1,J)))*uDdzOp(_I3(K+1,I,J)) |
165 |
ENDDO |
166 |
uGradzOpKP1(1,J) = |
167 |
& 0.5*( ZA(_I3(K+1,1,J))+ZA(_I3(K+1,nx,J)) )*uDdzOp(_I3(K+1,1,J)) |
168 |
ENDDO |
169 |
ENDIF |
170 |
C /---------------------------------------------------------------\ |
171 |
C | vGradXOp <- 1/DXu*By{xa} | |
172 |
C | vGradYOp <- 1/DYu*By{ya} | |
173 |
C | vGradZOp[K,KP1] <- 1/DZu*Bz{ZA[(K),(K+1)]} | |
174 |
C \---------------------------------------------------------------/ |
175 |
DO J = 2, ny |
176 |
DO I = 1, nx |
177 |
vGradxOp(I,J) = 0.5*( xaK(I,J)+xaK(I,J-1) )*vDdxOp(_I3(K,I,J)) |
178 |
vGradzOpK(I,J)= |
179 |
& 0.5*(ZA(_I3(K,I,J))+ZA(_I3(K,I,J-1)))*vDdzOp(_I3(K,I,J)) |
180 |
ENDDO |
181 |
ENDDO |
182 |
DO I = 1, nx |
183 |
vGradxOp (I,1) = 0.5*( xaK(I,1)+xaK(I,ny) )*vDdxOp(_I3(K,I,1 )) |
184 |
vGradzOpK(I,1) = 0.5*(ZA(_I3(K,I,1))+ZA(_I3(K,I,ny)))*vDdzOp(_I3(K,I,1)) |
185 |
ENDDO |
186 |
IF ( .NOT. BOTTOM_LAYER ) THEN |
187 |
DO J = 2, ny |
188 |
DO I = 1, nx |
189 |
vGradzOpKP1(I,J) = |
190 |
& 0.5*(ZA(_I3(K+1,I,J))+ZA(_I3(K+1,I,J-1)))*vDdzOp(_I3(K+1,I,J)) |
191 |
ENDDO |
192 |
ENDDO |
193 |
DO I = 1, nx |
194 |
vGradzOpKP1(I,1) = |
195 |
& 0.5*(ZA(_I3(K+1,I,1))+ZA(_I3(K+1,I,ny)))*vDdzOp(_I3(K+1,I,1)) |
196 |
ENDDO |
197 |
ENDIF |
198 |
DO J = 1, ny-1 |
199 |
DO I = 1, nx |
200 |
vGradyOp(I,J) = 0.5*( yaK(I,J)+yaK(I,J+1) )*vDdyOp(_I3(K,I,J)) |
201 |
ENDDO |
202 |
ENDDO |
203 |
DO I = 1, nx |
204 |
vGradyOp(I,ny) = 0.5*( yaK(I,ny)+yaK(I,1) )*vDdyOp(_I3(K,I,ny)) |
205 |
ENDDO |
206 |
C /---------------------------------------------------------------\ |
207 |
C | w[K,KP1] <- W[(K),(K+1)] | |
208 |
C \---------------------------------------------------------------/ |
209 |
wKP1 = 0. |
210 |
wK(1:nx,1:ny) = W(_I3(K,:,:)) |
211 |
wK(0,1:ny) = wK(nx,1:ny) |
212 |
wK(:,0) = wK(:,ny) |
213 |
IF ( .NOT. BOTTOM_LAYER ) wKP1(1:nx,1:ny) = W(_I3(K+1,:,:)) |
214 |
wKP1(0,1:ny) = wKP1(nx,1:ny) |
215 |
wKP1(:,0) = wKP1(:,ny) |
216 |
C /---------------------------------------------------------------\ |
217 |
C | uK <- U(K) | |
218 |
C \---------------------------------------------------------------/ |
219 |
uK(1:nx,1:ny) = U(_I3(K,:,:)) |
220 |
uK(nx+1,1:ny) = uK(1,1:ny) |
221 |
uK(:,0) = uK(:,ny) |
222 |
uKNM1(1:nx,1:ny) = UNM1(_I3(K,:,:)) |
223 |
uKNM1(nx+1,1:ny) = uKNM1(1,1:ny) |
224 |
uKNM1(:,0) = uKNM1(:,ny) |
225 |
uMskK(1:nx,1:ny) = UMASK(_I3(K,:,:)) |
226 |
uMskK(nx+1,1:ny) = uMskK(1,1:ny) |
227 |
uMskK(:,0) = uMskK(:,ny) |
228 |
C /---------------------------------------------------------------\ |
229 |
C | v[KM1,K,KP1] = V[(K-1),(K),(K+1)] | |
230 |
C \---------------------------------------------------------------/ |
231 |
vKP1 = 0. |
232 |
vKM1 = 0. |
233 |
vK(1:nx,1:ny) = V(_I3(K,:,:)) |
234 |
vK(0,1:ny) = vK(nx,1:ny) |
235 |
vK(:,ny+1) = vK(:,1) |
236 |
vKNM1(1:nx,1:ny) = VNM1(_I3(K,:,:)) |
237 |
vKNM1(0,1:ny) = vKNM1(nx,1:ny) |
238 |
vKNM1(:,ny+1) = vKNM1(:,1) |
239 |
vMskK(1:nx,1:ny) = VMASK(_I3(K,:,:)) |
240 |
vMskK(0,1:ny) = vMskK(nx,1:ny) |
241 |
vMskK(:,ny+1) = vMskK(:,1) |
242 |
IF ( .NOT. TOP_LAYER ) vKM1(1:nx,1:ny) = V(_I3(K-1,:,:)) |
243 |
IF ( .NOT. BOTTOM_LAYER ) vKP1(1:nx,1:ny) = V(_I3(K+1,:,:)) |
244 |
C /---------------------------------------------------------------\ |
245 |
C | Gu Section | |
246 |
C |===============================================================| |
247 |
C | In this section Gu = GuDiffusion | |
248 |
C | + GuAdvection | |
249 |
C | + GuMetric | |
250 |
C | + GuXZCoriolis | |
251 |
C | + GuForcing | |
252 |
C \---------------------------------------------------------------/ |
253 |
C /---------------------------------------------------------------\ |
254 |
C | Initialise local terms | |
255 |
C |===============================================================| |
256 |
C | rVk <- 1/Vu | |
257 |
C \---------------------------------------------------------------/ |
258 |
rVk = rUvol(_I3(K,:,:)) |
259 |
fZK = 0. |
260 |
fZKP1 = 0. |
261 |
C /---------------------------------------------------------------\ |
262 |
C | XY GuDiffusion + XY GuAdvection | |
263 |
C |===============================================================| |
264 |
C | Note: XY GuDiffusion = Laplacian + Biharmonic in XY-plane. | |
265 |
C | o Laplacian ( -a2MomXY.del^2(u) ) | |
266 |
C | o Biharmonic ( +a4MomXY.del^4(u) ) | |
267 |
C \---------------------------------------------------------------/ |
268 |
DO J = 1, ny |
269 |
DO I = 1, nx |
270 |
C /-------------------------------------------------------------\ |
271 |
C | fX <- 1/DXu*XAu*Ddx{U} | |
272 |
C | fY <- 1/DYu*YAu*Ddy{U} | |
273 |
C \-------------------------------------------------------------/ |
274 |
fX(I,J) = uGradxOp(I,J)*(uK(I+1,J)-uK(I,J)) |
275 |
fY(I,J) = uGradyOp(I,J)*(uK(I,J)-uK(I,J-1)) |
276 |
ENDDO |
277 |
ENDDO |
278 |
fX(0 ,1:ny) = fX(nx,1:ny) |
279 |
fY(1:nx,ny+1) = fY(1:nx,1) |
280 |
C WRITE(0,*) ' CALC_G fY ' |
281 |
C CALL PLOT_FIELD( fY(1:nx,1:ny), nX, nY ) |
282 |
DO J=1,ny |
283 |
DO I=1,nx |
284 |
C /-------------------------------------------------------------\ |
285 |
C | tmp <- 1/Vu*(Dx{fX}+Dy{fY}) | |
286 |
C \-------------------------------------------------------------/ |
287 |
tmp(I,J)=rVk(I,J)*fX(I,J)-rVk(I,J)*fX(I-1,J) |
288 |
& +rVk(I,J)*fY(I,J+1)-rVk(I,J)*fY(I,J) |
289 |
ENDDO |
290 |
ENDDO |
291 |
tmp(nx+1,1:ny) = tmp(1,1:ny) |
292 |
tmp(:,0) = tmp(:,ny) |
293 |
DO J = 1,ny |
294 |
DO I = 1,nx |
295 |
C /-------------------------------------------------------------\ |
296 |
C | Diffusive and advective flux at eastern face | |
297 |
C |=============================================================| |
298 |
C | fX <- -k2*fX + ( Laplacian,<Line 1> ) | |
299 |
C | k4/DXu*XAu*Ddx{tmp} + ( Biharmonic,<Line 2> ) | |
300 |
C | Bx{XA*U}*Bx{U} ( Advection,<Line 3-4> ) | |
301 |
C \-------------------------------------------------------------/ |
302 |
fX(I,J) = 0. |
303 |
_LPM(& -_DB1*a2MomXY*fX(I,J) ) |
304 |
_BHM(& +a4MomXY*uGradxOp(I,J)*(tmp(I+1,J)-tmp(I,J)) ) |
305 |
_ADM(& +0.25D0*(uK(I,J)*XAK(I,J)+uK(I+1,J)*XAK(I+1,J))* ) |
306 |
_ADM(& (uK(I,J)+uK(I+1,J)) ) |
307 |
C /-------------------------------------------------------------\ |
308 |
C | Diffusive and advective flux at southern face | |
309 |
C |=============================================================| |
310 |
C | fY <- -k2*fY + ( Laplacian,<Line 1> ) | |
311 |
C | k4/DYu*YAu*Ddy{tmp} + ( Biharmonic,<Line 2> ) | |
312 |
C | Bx{YA*V}*By{U} ( Advection,<Line 3-4> ) | |
313 |
C \-------------------------------------------------------------/ |
314 |
fY(I,J) = 0. |
315 |
_LPM(& -_DB1*a2MomXY*fY(I,J) ) |
316 |
_BHM(& +a4MomXY*uGradyOp(I,J)*(tmp(I,J)-tmp(I,J-1)) ) |
317 |
_ADM(& +0.25D0*(uK(I,J)+uK(I,J-1))* ) |
318 |
_ADM(& (vK(I,J)*YAK(I,J)+vK(I-1,J)*YAK(I-1,J)) ) |
319 |
_ADM(& *uMskK(I,J)*uMskK(I,J-1) ) |
320 |
C Note!!!! The last line above "*uMskK(I,J)*uMskK(I,J-1)" is the |
321 |
C Note!!!! boundary condition in the standard v010 CM5 code. The |
322 |
C Note!!!! FV paper used a different bc in which the last line is |
323 |
C Note!!!! omitted. |
324 |
ENDDO |
325 |
ENDDO |
326 |
|
327 |
C /---------------------------------------------------------------\ |
328 |
C | Z GuDiffusion + Z GuAdvection | |
329 |
C |===============================================================| |
330 |
C \---------------------------------------------------------------/ |
331 |
IF ( .NOT. TOP_LAYER ) THEN |
332 |
DO J = 1, ny |
333 |
DO I = 1, nx |
334 |
C /------------------------------------------------------------\ |
335 |
C | Diffusive and advective flux at upper face | |
336 |
C |============================================================| |
337 |
C | fZK <- -k2Z/DZu*ZAu*Ddz{U} + ( Laplacian,<Line 1> ) | |
338 |
C | Bx{ZA*-w}*Bz{U} ( Advection,<Line 2-3> ) | |
339 |
C \------------------------------------------------------------/ |
340 |
fZK(I,J) = 0. |
341 |
_LPM(& -_DB1*a2MomZ*uGradzOpK(I,J)*( U(_I3(K-1,I,J))-uK(I,J) ) ) |
342 |
_ADM(& +0.25D0*(uK(I,J)+U(_I3(K-1,I,J)))* ) |
343 |
_ADM(& (-wK(I,J)*ZAK(I,J)-wK(I-1,J)*ZAK(I-1,J)) ) |
344 |
_ADM(& *uMskK(I,J)*UMASK(_I3(K-1,I,J)) ) |
345 |
C Note!!!! The last line above "*uMskK(I,J)*UMASK(_I3(K-1,I,J))" is |
346 |
C Note!!!! boundary condition in the standard v010 CM5 code. The |
347 |
C Note!!!! FV paper used a different bc in which the last line is |
348 |
C Note!!!! omitted. |
349 |
ENDDO |
350 |
ENDDO |
351 |
ELSE |
352 |
DO J = 1, ny |
353 |
DO I = 1, nx |
354 |
C /------------------------------------------------------------\ |
355 |
C | Advective flux at top layer | |
356 |
C |============================================================| |
357 |
C | o Only active when using implcit free surface. | |
358 |
C \------------------------------------------------------------/ |
359 |
fZK(I,J) = 0. |
360 |
_ADM(& +freeSurfFac*0.25D0*(uK(I,J)+uK(I,J) )* ) |
361 |
_ADM(& (-wK(I,J)*ZAK(I,J)-wK(I-1,J)*ZAK(I-1,J)) ) |
362 |
ENDDO |
363 |
ENDDO |
364 |
ENDIF |
365 |
IF ( .NOT. BOTTOM_LAYER ) THEN |
366 |
DO J = 1, ny |
367 |
DO I = 1, nx |
368 |
C /------------------------------------------------------------\ |
369 |
C | Diffusive and advective flux at lower face | |
370 |
C |============================================================| |
371 |
C | fZKP1 <- -k2Z/DZu*ZAu*Ddz{U} + ( Laplacian,<Line 1> ) | |
372 |
C | Bx{ZA*-w}*Bz{U} ( Advection,<Line 2-3> ) | |
373 |
C \------------------------------------------------------------/ |
374 |
fZKP1(I,J) = 0. |
375 |
_LPM(& -_DB1*a2MomZ*uGradzOpKP1(I,J)*(uK(I,J)-U(_I3(K+1,I,J))) ) |
376 |
_ADM(& +0.25D0*(uK(I,J)+U(_I3(K+1,I,J)))* ) |
377 |
_ADM(& (-wKP1(I,J)*ZAKP1(I,J)-wKP1(I-1,J)*ZAKP1(I-1,J)) ) |
378 |
_ADM(& *uMskK(I,J)*UMASK(_I3(K+1,I,J)) ) |
379 |
C Note!!!! The last line above "*uMskK(I,J)*UMASK(_I3(K-1,I,J))" is |
380 |
C Note!!!! boundary condition in the standard v010 CM5 code. The |
381 |
C Note!!!! FV paper used a different bc in which the last line is |
382 |
C Note!!!! omitted. |
383 |
ENDDO |
384 |
ENDDO |
385 |
ELSE |
386 |
DO J = 1, ny |
387 |
DO I = 1, nx |
388 |
fZKP1(I,J) = 0. |
389 |
_LPM(& -_DB1*a2MomZ*rDzAtP(K)*2.D0*uK(I,J) ) |
390 |
_LPM(& *0.5*( ZAk(I,J)+ZAk(I-1,J) ) ) |
391 |
ENDDO |
392 |
ENDDO |
393 |
ENDIF |
394 |
C /---------------------------------------------------------------\ |
395 |
C | gU <- GuDiffusion + GuAdvection = -div(fX,fY,fZ) | |
396 |
C \---------------------------------------------------------------/ |
397 |
fX(0,1:ny) = fX(nx,1:ny) |
398 |
fY(:,ny+1) = fY(:,1) |
399 |
DO J = 1, ny |
400 |
DO I = 1, nx |
401 |
C /-------------------------------------------------------------\ |
402 |
C | o Minus divergence of the lateral fluxes | |
403 |
C | gU <- -1/Vu*(Ddx{fX}+Ddy{fY}) | |
404 |
C \-------------------------------------------------------------/ |
405 |
gU(_I3(K,I,J))=-fX(I,J)*rVk(I,J)+fX(I-1,J)*rVk(I,J) |
406 |
& -fY(I,J+1)*rVk(I,J)+fY(I,J)*rVk(I,J) |
407 |
ENDDO |
408 |
ENDDO |
409 |
DO J = 1, ny |
410 |
DO I = 1, nx |
411 |
C /------------------------------------------------------------\ |
412 |
C | o Minus flux at upper face | |
413 |
C | gU <- gU-1/Vu*fZK | |
414 |
C \------------------------------------------------------------/ |
415 |
gU(_I3(K,I,J))=gU(_I3(K,I,J))-fZK(I,J)*rVk(I,J) |
416 |
ENDDO |
417 |
ENDDO |
418 |
CcnhDebugStarts |
419 |
C Changed to make noslip bottom boundary. |
420 |
Cdbg IF ( .NOT. BOTTOM_LAYER ) THEN |
421 |
CcnhDebugEnds |
422 |
DO J = 1, ny |
423 |
DO I = 1, nx |
424 |
C /------------------------------------------------------------\ |
425 |
C | o Plus flux at lower face | |
426 |
C | gU <- gU+1/Vu*fZP1 | |
427 |
C \------------------------------------------------------------/ |
428 |
gU(_I3(K,I,J))=gU(_I3(K,I,J))+fZKP1(I,J)*rVk(I,J) |
429 |
ENDDO |
430 |
ENDDO |
431 |
CcnhDebugStarts |
432 |
C Changed to make noslip bottom boundary. |
433 |
Cdbg ENDIF |
434 |
CcnhDebugEnds |
435 |
|
436 |
C /---------------------------------------------------------------\ |
437 |
C | gU <- gU + GuMetric + GuXZCoriolis | |
438 |
C \---------------------------------------------------------------/ |
439 |
DO J = 1, ny |
440 |
DO I = 1, nx |
441 |
# ifdef _SPHERICAL_POLAR_METRIC_TERMS |
442 |
C /-------------------------------------------------------------\ |
443 |
C | o Spherical polar metric terms for U | |
444 |
C |=============================================================| |
445 |
C | gU <- gU - U/aEarth*BzBx{-w/G/RONIL} - U*tan(lat)*ByBx{V} | |
446 |
C \-------------------------------------------------------------/ |
447 |
gU(_I3(K,I,J)) = gU(_I3(K,I,J)) |
448 |
_SPM(& -uMask(_I3(K,I,J))*uK(I,J)/rZero* ) |
449 |
_SPM(&0.25D0*(wK(I-1,J)+wK(I,J)+wKP1(I-1,J)+wKP1(I,J))/(-G*RONIL)) |
450 |
_SPM(& +uMask(_I3(K,I,J))*uK(I,J)*uTphiOP(_I3(K,I,J))/rZero ) |
451 |
_SPM(& *0.25D0*(vK(I,J)+vK(I-1,J)+vK(I,J+1)+vK(I-1,J+1)) ) |
452 |
# endif |
453 |
C /-------------------------------------------------------------\ |
454 |
C | o XZ-plane coriolis terms for U | |
455 |
C |=============================================================| |
456 |
C | gU <- gU - 2*omega*cos(lat)*BzBx{-w/G/RONIL} | |
457 |
C \-------------------------------------------------------------/ |
458 |
gU(_I3(K,I,J)) = gU(_I3(K,I,J)) |
459 |
_XZC(& -fCorW(_I3(K,I,J))/ ) |
460 |
_XZC(&G/(-RONIL)*0.25D0*(wK(I-1,J)+wK(I,J)+wKP1(I-1,J)+wKP1(I,J))) |
461 |
ENDDO |
462 |
ENDDO |
463 |
|
464 |
#ifdef _MOMENTUM_FORCING |
465 |
C /---------------------------------------------------------------\ |
466 |
C | gU <- gU + GuForcing | |
467 |
C \---------------------------------------------------------------/ |
468 |
IF ( TOP_LAYER ) THEN |
469 |
C /--------------------------------------------------------------\ |
470 |
C | o Add a drag term to the top-layer of U | |
471 |
C | u* is the velocity of the disk relative to the surface. | |
472 |
C | lambda is the inverse of the time scale with which the | |
473 |
C | flow equilibrates withthe rotating disk. | |
474 |
C | to get lambda we use the following | |
475 |
C | du/dt = stress/rho/dz | |
476 |
C | stress = rho.Kz.DU/DZ | |
477 |
C | DZ = (2Kz/f)^1/2 (Ekman layer thiickness ) | |
478 |
C | DU = (U-u*) | |
479 |
C | F=1e-4, Kz=1e-3, dz of 25m gives an equilibration time | |
480 |
C | scale of 1.3 days. | |
481 |
C | gU(K=1) <- gU(K=1) - lambda(U-u*) | |
482 |
C \--------------------------------------------------------------/ |
483 |
gU(_I3(K,:,:)) = gU(_I3(K,:,:)) |
484 |
& + fu |
485 |
ENDIF |
486 |
#endif |
487 |
|
488 |
C ********************************************************************* |
489 |
C ****Gv SECTION******************************************************* |
490 |
C ********************************************************************* |
491 |
C In this section Gv = GvDiffusion |
492 |
C + GvAdvection |
493 |
C + GvMetric |
494 |
C + GvHorizontalCoriolis |
495 |
C + GvForcing |
496 |
C === Initialise local operators ====================================== |
497 |
C rVk: 1/VVolume |
498 |
fZK = 0. |
499 |
fZKP1 = 0. |
500 |
rVk = rVvol(_I3(K,:,:)) |
501 |
C === Lateral diffusive flux of V ===================================== |
502 |
C o Laplacian ( -a2MomXY.del^2(v) ) |
503 |
C o Biharmonic ( +a4MomXY.del^4(v) ) |
504 |
DO J = 1, ny |
505 |
DO I = 1, nx |
506 |
C fX: VAreaX*dx[V] |
507 |
C fY: VAreaY*dy[V] |
508 |
fX(I,J) = vGradxOp(I,J)*(vK(I,J)-vK(I-1,J)) |
509 |
fY(I,J) = vGradyOp(I,J)*(vK(I,J+1)-vK(I,J)) |
510 |
ENDDO |
511 |
ENDDO |
512 |
fX(nx+1,1:ny) = fX(1 ,1:ny) |
513 |
fY(1:nx,0 ) = fY(1:nx,ny) |
514 |
DO J=1,ny |
515 |
DO I=1,nx |
516 |
C o Del^2 V ( for use in biharmonic mixing ). |
517 |
C tmp: div(fX,fY) = 1/VVolume*(dx[fX]+dy[fy]) |
518 |
tmp(I,J)= fX(I+1,J)*rVk(I,J)-fX(I,J)*rVk(I,J) |
519 |
& +fY(I,J)*rVk(I,J)-fY(I,J-1)*rVk(I,J) |
520 |
ENDDO |
521 |
ENDDO |
522 |
tmp(0,1:ny) = tmp(nx,1:ny) |
523 |
tmp(:,ny+1) = tmp(:,1) |
524 |
DO J = 1,ny |
525 |
DO I = 1,nx |
526 |
C o Laplacian<Line 1> + Biharmonic<Line 2> + Advection<Line 3,4>. |
527 |
C fX: -K2xyMomentum*fX + K4xyMomentum*VAreaX*dx[tmp] + By[Ax*U]*Bx[V] |
528 |
C fY: -K2xyMomentum*fY + K4xyMomentum*VAreaY*dy[tmp] + By[Ay*V]*By[V] |
529 |
fX(I,J) = 0. |
530 |
_LPM(& -_DB1*a2MomXY*fX(I,J) ) |
531 |
_BHM(& +a4MomXY*vGradxOp(I,J)*(tmp(I,J)-tmp(I-1,J)) ) |
532 |
_ADM(& +0.25D0*(uK(I,J-1)*XAK(I,J-1)+uK(I,J)*XAK(I,J))* ) |
533 |
_ADM(& (vK(I-1,J)+vK(I,J)) ) |
534 |
_ADM(& *vMskK(I,J)*vMskK(I-1,J) ) |
535 |
C Note!!!! The last line above "*vMskK(I,J)*vMskK(I-1,J)" is |
536 |
C Note!!!! boundary condition in the standard v010 CM5 code. The |
537 |
C Note!!!! FV paper used a different bc in which the last line is |
538 |
C Note!!!! omitted. |
539 |
fY(I,J) = 0. |
540 |
_LPM(& -_DB1*a2MomXY*fY(I,J) ) |
541 |
_BHM(& +a4MomXY*vGradyOp(I,J)*(tmp(I,J+1)-tmp(I,J)) ) |
542 |
_ADM(& +0.25D0*(vK(I,J+1)*YAK(I,J+1)+vK(I,J)*YAK(I,J))* ) |
543 |
_ADM(& (vK(I,J+1)+vK(I,J)) ) |
544 |
C Note!!!! The last line above "*vMskK(I,J)*vMskK(I-1,J)" is |
545 |
ENDDO |
546 |
ENDDO |
547 |
C === Vertical diffusive flux of V ==================================== |
548 |
C === Vertical advective flux of V ==================================== |
549 |
C o Laplacian Diffusion<Line 1> + Advection<Line 2,3> |
550 |
C fZK: -KzMomentum*VAreaZ*dz[V] + Bz[V]*By[Az*W] |
551 |
IF ( .NOT. TOP_LAYER ) THEN |
552 |
c WRITE(0,*) ' fZK subsurface layer ' |
553 |
DO J = 1, ny |
554 |
DO I = 1, nx |
555 |
fZK(I,J) = 0. |
556 |
_LPM(& -_DB1*a2MomZ*vGradzOpK(I,J)*(V(_I3(K-1,I,J))-vK(I,J) ) ) |
557 |
_ADM(& +0.25D0*(vK(I,J)+V(_I3(K-1,I,J)))* ) |
558 |
_ADM(& (-wK(I,J)*ZAK(I,J)-wK(I,J-1)*ZAK(I,J-1)) ) |
559 |
_ADM(& *vMskK(I,J)*vMask(_I3(K-1,I,J))*_wdvdz ) |
560 |
C Note!!!! The last line above "*vMskK(I,J)*vMask(_I3(K-1,I,J))" is |
561 |
C Note!!!! boundary condition in the standard v010 CM5 code. The |
562 |
C Note!!!! FV paper used a different bc in which the last line is |
563 |
C Note!!!! omitted. |
564 |
ENDDO |
565 |
ENDDO |
566 |
ELSE |
567 |
c WRITE(0,*) ' fZK top layer ' |
568 |
DO J = 1, ny |
569 |
DO I = 1, nx |
570 |
fZK(I,J) = 0.D0 |
571 |
_ADM(& +freeSurfFac*0.25D0*(vK(I,J)+vK(I,J)) ) |
572 |
_ADM(& *(-wK(I,J)*ZAK(I,J)-wK(I,J-1)*ZAK(I,J-1))*_wdvdz ) |
573 |
ENDDO |
574 |
ENDDO |
575 |
ENDIF |
576 |
IF ( .NOT. BOTTOM_LAYER ) THEN |
577 |
DO J = 1, ny |
578 |
DO I = 1, nx |
579 |
fZKP1(I,J) = 0. |
580 |
_LPM(& -_DB1*a2MomZ*vGradzOpKP1(I,J)*( vK(I,J)-V(_I3(K+1,I,J))) ) |
581 |
_ADM(& +0.25D0*(vK(I,J)+V(_I3(K+1,I,J)))* ) |
582 |
_ADM(& (-wKP1(I,J)*ZAKP1(I,J)-wKP1(I,J-1)*ZAKP1(I,J-1)) ) |
583 |
_ADM(& *vMskK(I,J)*vMask(_I3(K+1,I,J))*_wdvdz ) |
584 |
C Note!!!! The last line above "*vMskK(I,J)*vMask(_I3(K+1,I,J))" is |
585 |
C Note!!!! boundary condition in the standard v010 CM5 code. The |
586 |
C Note!!!! FV paper used a different bc in which the last line is |
587 |
C Note!!!! omitted. |
588 |
ENDDO |
589 |
ENDDO |
590 |
ELSE |
591 |
DO J = 1, ny |
592 |
DO I = 1, nx |
593 |
fZKP1(I,J) = 0. |
594 |
_LPM(& -_DB1*a2MomZ*rDzAtP(K)*2.*vK(I,J) ) |
595 |
_LPM(& *0.5*( ZAk(I,J)+ZAk(I,J-1) ) ) |
596 |
|
597 |
ENDDO |
598 |
ENDDO |
599 |
ENDIF |
600 |
c WRITE(0,*) 'K, MAXVAL(fZK) = ', K, MAXVAL(fZK) |
601 |
c WRITE(0,*) 'K, MAXVAL(fZKP1) = ', K, MAXVAL(fZKP1) |
602 |
C === Combine advective and diffusive fluxes of V ===================== |
603 |
C o V tendency is minus divergence of diffusive and advective flux. |
604 |
C gV: -div(fX,fY,fZ) = -1/VVolume*(dx[fX]+dy[fY]+fZK-fZKP1) |
605 |
fX(nx+1,1:ny) = fX(1,1:ny) |
606 |
fY(:,0) = fY(:,ny) |
607 |
DO I = 1, nx |
608 |
DO J = 1, ny |
609 |
gV(_I3(K,I,J))=-fX(I+1,J)*rVK(I,J)+fX(I,J)*rVK(I,J) |
610 |
& -fY(I,J)*rVK(I,J)+fY(I,J-1)*rVK(I,J) |
611 |
ENDDO |
612 |
ENDDO |
613 |
DO J = 1, ny |
614 |
DO I = 1, nx |
615 |
gV(_I3(K,I,J))=gV(_I3(K,I,J))-fZK(I,J)*rVK(I,J) |
616 |
ENDDO |
617 |
ENDDO |
618 |
CcnhDebugStarts |
619 |
C Changed to provide noslip bottom boundary |
620 |
Cdbg IF ( .NOT. BOTTOM_LAYER ) THEN |
621 |
CcnhDebugEnds |
622 |
DO J = 1, ny |
623 |
DO I = 1, nx |
624 |
gV(_I3(K,I,J))=gV(_I3(K,I,J))+fZKP1(I,J)*rVK(I,J) |
625 |
ENDDO |
626 |
ENDDO |
627 |
CcnhDebugStarts |
628 |
Cdbg ENDIF |
629 |
CcnhDebugEnds |
630 |
C === Metric terms due to curvilinear coordinates ===================== |
631 |
DO J = 1, ny |
632 |
DO I = 1, nx |
633 |
C o Spherical polar metric terms. |
634 |
C gV: gV - V/radius*ByBz[w]+tan(latitude)*BxBy[u] |
635 |
gV(_I3(K,I,J)) = gV(_I3(K,I,J)) |
636 |
_SPM(& -vMask(_I3(K,I,J))*vK(I,J)/rZero* ) |
637 |
_SPM(& 0.25D0*(wK(I,J-1)+wK(I,J)+wKP1(I,J-1)+wKP1(I,J))/(-G*RONIL) ) |
638 |
_SPM(& -vMask(_I3(K,I,J))*vTphiOP(_I3(K,I,J))/rZero*0.25D0*0.25D0 ) |
639 |
_SPM(& *(uK(I,J)+uK(I+1,J)+uK(I,J-1)+uK(I+1,J-1)) ) |
640 |
_SPM(& *(uK(I,J)+uK(I+1,J)+uK(I,J-1)+uK(I+1,J-1)) ) |
641 |
ENDDO |
642 |
ENDDO |
643 |
|
644 |
#ifdef _MOMENTUM_FORCING |
645 |
IF ( TOP_LAYER ) THEN |
646 |
gV(_I3(K,:,:)) = gV(_I3(K,:,:)) |
647 |
& + fv |
648 |
ENDIF |
649 |
#endif |
650 |
|
651 |
C Hydrostatic pressure term |
652 |
DO J=1,Ny |
653 |
DO I=1,Nx |
654 |
tmp(I,J)=G*PH(_I3(K,I,J)) |
655 |
ENDDO |
656 |
ENDDO |
657 |
DO J=1,Ny |
658 |
tmp(0,J)=tmp(Nx,J) |
659 |
ENDDO |
660 |
DO I=1,Nx |
661 |
tmp(I,0)=tmp(I,Ny) |
662 |
ENDDO |
663 |
DO J = 1, ny |
664 |
DO I = 1, nx |
665 |
gU(_I3(K,I,J)) = gU(_I3(K,I,J)) -pDdxOp(_I3(K,I,J))*(tmp(I,J)-tmp(I-1,J)) |
666 |
gV(_I3(K,I,J)) = gV(_I3(K,I,J)) -pDdyOp(_I3(K,I,J))*(tmp(I,J)-tmp(I,J-1)) |
667 |
ENDDO |
668 |
ENDDO |
669 |
|
670 |
|
671 |
#ifdef _XY_CORIOLIS |
672 |
C /---------------------------------------------------------------\ |
673 |
C | Cd Scheme | |
674 |
C |===============================================================| |
675 |
C | o Gu = Gu + GuXYCoriolis | |
676 |
C | o Gv = Gv + GvXYCoriolis | |
677 |
C \---------------------------------------------------------------/ |
678 |
C /---------------------------------------------------------------\ |
679 |
C | tmp2 <- P(n) = P extrapolated to time level n. | |
680 |
C \---------------------------------------------------------------/ |
681 |
tmp2(1:nx,1:ny) = (1.5+abEps)*(pS+pH(_I3(K,:,:)) ) |
682 |
& -(0.5+abEps)*(psNM1+phNM1(_I3(K,:,:)) ) |
683 |
tmp2(1:nx,1:ny) = (1.5+abEps)*(pS ) |
684 |
& -(0.5+abEps)*(psNM1 ) |
685 |
tmp2(0,1:ny) = tmp2(nx,1:ny) |
686 |
DO J = 1,ny |
687 |
DO I = 1,nx |
688 |
C /-------------------------------------------------------------\ |
689 |
C | tmp <- -G*Ddx{P(n)} + gU | |
690 |
C \-------------------------------------------------------------/ |
691 |
tmp(I,J) = |
692 |
& -G*pDdxOp(_I3(K,I,J))*(tmp2(I,J)-tmp2(I-1,J))+ |
693 |
& gU(_I3(K,I,J))*uMask(_I3(K,I,J)) |
694 |
ENDDO |
695 |
ENDDO |
696 |
tmp(nx+1,1:ny) = tmp(1,1:ny) |
697 |
tmp(:,0 ) = tmp(:,ny ) |
698 |
DO J = 1,ny |
699 |
DO I = 1,nx |
700 |
C /-------------------------------------------------------------\ |
701 |
C | uAja <- uAja + DT*BxBy{-d/dx P(n) + gU}+f*V | |
702 |
C \-------------------------------------------------------------/ |
703 |
uAja(_I3(K,I,J)) = uAja(_I3(K,I,J)) |
704 |
& +1.*DelT*(0.*tmp(I,J)+ |
705 |
& 1.*0.25D0*(tmp(I,J)+tmp(I+1,J)+tmp(I,J-1)+tmp(I+1,J-1)) |
706 |
& +fCorV(_I3(K,I,J))*V(_I3(K,I,J)) ) |
707 |
C uAja: uAja - lambda(uAja-BxBy[U(n+0.5)] ( note rAja = (1-lambda) ) |
708 |
uAja(_I3(K,I,J)) = rAja*uAja(_I3(K,I,J)) |
709 |
& +(1.-rAja)*( (1.5+abEps)*0.25D0*(uK(I,J)+uK(I+1,J)+uK(I,J-1)+uK(I+1,J-1)) |
710 |
& -(0.5+abEps)*0.25D0*(uKNM1(I,J)+uKNM1(I+1,J)+uKNM1(I,J-1)+uKNM1(I+1,J-1)) ) |
711 |
ENDDO |
712 |
ENDDO |
713 |
CcnhDebugStarts |
714 |
C WRITE(0,*) ' CG2D MAXVAL(uAja) ', MAXVAL(uAja) |
715 |
CcnhDebugEnds |
716 |
C tmp2: P(n) = P extrapolated to time level n. |
717 |
tmp2(1:nx,1:ny) = (1.5+abEps)*(pS+pH(_I3(K,:,:)) ) |
718 |
& -(0.5+abEps)*(psNM1+phNM1(_I3(K,:,:)) ) |
719 |
tmp2(1:nx,1:ny) = (1.5+abEps)*(pS ) |
720 |
& -(0.5+abEps)*(psNM1 ) |
721 |
tmp2(1:nx,0) = tmp2(1:nx,ny) |
722 |
DO J = 1,ny |
723 |
DO I = 1,nx |
724 |
C tmp: -d/dy P(n) + gV |
725 |
tmp(I,J) = -G*pDdyOp(_I3(K,I,J))*(tmp2(I,J)-tmp2(I,J-1))+ |
726 |
& gV(_I3(K,I,J))*vMask(_I3(K,I,J)) |
727 |
ENDDO |
728 |
ENDDO |
729 |
tmp(0,1:ny) = tmp(nx,1:ny) |
730 |
tmp(:,ny+1) = tmp(: ,1 ) |
731 |
DO J = 1,ny |
732 |
DO I = 1,nx |
733 |
C vAja: vAja + DT*(BxBy[-d/dy P(n) + gV]-f*U) |
734 |
vAja(_I3(K,I,J)) = vAja(_I3(K,I,J)) |
735 |
& +DelT*(0.25D0*(tmp(I-1,J)+tmp(I,J)+tmp(I-1,J+1)+tmp(I,J+1)) |
736 |
& -fCorU(_I3(K,I,J))*U(_I3(K,I,J)) ) |
737 |
C vAja: vAja - lambda(vAja-BxBy[V(n+0.5)] ( note rAja = (1-lambda) ) |
738 |
vAja(_I3(K,I,J)) = rAja*vAja(_I3(K,I,J)) |
739 |
& +(1.-rAja)*( (1.5+abEps)*0.25D0*(vK(I,J)+vK(I-1,J)+vK(I,J+1)+vK(I-1,J+1)) |
740 |
& -(0.5+abEps)*0.25D0*(vKNM1(I,J)+vKNM1(I-1,J)+vKNM1(I,J+1)+vKNM1(I-1,J+1)) ) |
741 |
ENDDO |
742 |
ENDDO |
743 |
#endif |
744 |
|
745 |
C Add coriolis term from Cd grid into Gu and Gv. |
746 |
C Note: AB2 does not apply to Cd variables. |
747 |
gU(_I3(K,:,:)) = gU(_I3(K,:,:))*uMask(_I3(K,:,:)) |
748 |
tmp(1:nx,1:ny) = gUNM1(_I3(K,:,:)) |
749 |
gUNM1(_I3(K,:,:)) = gU(_I3(K,:,:)) |
750 |
gU(_I3(K,:,:)) = |
751 |
& (1.5+abEps)*gU(_I3(K,:,:))-(0.5+abEps)*tmp(1:nx,1:ny) |
752 |
& +fCorU(_I3(K,:,:))*vAja(_I3(K,:,:))*uMask(_I3(K,:,:)) |
753 |
IF ( BOTTOM_LAYER ) THEN |
754 |
_D(( ' S/R G_CALC: MAXVAL(GU)', MAXVAL(GU) )) |
755 |
_D(( ' MAXVAL(GU) @', MAXLOC(GU) )) |
756 |
ENDIF |
757 |
|
758 |
gV(_I3(K,:,:)) = gV(_I3(K,:,:))*vMask(_I3(K,:,:)) |
759 |
tmp(1:nx,1:ny) = gVNM1(_I3(K,:,:)) |
760 |
gVNM1(_I3(K,:,:)) = gV(_I3(K,:,:)) |
761 |
gV(_I3(K,:,:)) = |
762 |
& (1.5+abEps)*gV(_I3(K,:,:))-(0.5+abEps)*tmp(1:nx,1:ny) |
763 |
& -fCorV(_I3(K,:,:))*uAja(_I3(K,:,:))*vMask(_I3(K,:,:)) |
764 |
_D(( ' S/R G_CALC: MAXVAL(GV)', MAXVAL(GV) )) |
765 |
C ********************************************************************* |
766 |
C *** FILTER SECTION ************************************************** |
767 |
C ********************************************************************* |
768 |
C Apply filter to control energy build up at high-latitudes. |
769 |
|
770 |
C ********************************************************************* |
771 |
C *** RHS SECTION ***************************************************** |
772 |
C ********************************************************************* |
773 |
C Form right-hand side of elliptic eqaution |
774 |
C RHS2d: dx[Ax*Gu]+dy[Ay*Gv] |
775 |
C +1/delt*(dx[Ax*U]+dy[Ay*V]) |
776 |
tmp(1:nx,1:ny) = |
777 |
& gU(_I3(K,:,:))/RONIL/G |
778 |
& +rDelt*Uk(1:nx,1:ny)/RONIL/G |
779 |
tmp(nx+1,1:ny) = tmp(1,1:ny) |
780 |
DO J=1,ny |
781 |
DO I=1,nx |
782 |
RHS2d(I,J)=RHS2d(I,J)+XAk(I+1,J)*tmp(I+1,J)-XAk(I,J)*tmp(I,J) |
783 |
ENDDO |
784 |
ENDDO |
785 |
tmp(1:nx,1:ny) = |
786 |
& gV(_I3(K,:,:))/RONIL/G |
787 |
& +rDelt*Vk(1:nx,1:ny)/RONIL/G |
788 |
tmp(1:nx,ny+1) = tmp(1:nx,1) |
789 |
DO J=1,ny |
790 |
DO I=1,nx |
791 |
RHS2d(I,J)=RHS2d(I,J)+YAk(I,J+1)*tmp(I,J+1)-YAk(I,J)*tmp(I,J) |
792 |
ENDDO |
793 |
ENDDO |
794 |
CcnhDebugStarts |
795 |
IF ( BOTTOM_LAYER ) THEN |
796 |
C DO I=1,Nx |
797 |
C WRITE(0,*) ' CG2D RHS2d(I,30) ', rhs2d(I,30), I |
798 |
C ENDDO |
799 |
C DO J=1,Ny |
800 |
C WRITE(0,*) ' CG2D MAXVAL(RHS2d) ', |
801 |
C & MAXVAL(RHS2d(1:Nx,J)), J |
802 |
C WRITE(0,*) ' CG2D MINVAL(RHS2d) ', |
803 |
C & MINVAL(RHS2d(1:Nx,J)), J |
804 |
C & MAXVAL(ABS(fCorV(_I3(:,I,:))*uAja(_I3(:,I,:))*vMask(_I3(:,I,:)))) |
805 |
C ENDDO |
806 |
C CALL PLOT_FIELD( RHS2d(1:nx,1:ny), nX, nY ) |
807 |
C STOP |
808 |
ENDIF |
809 |
CcnhDebugEnds |
810 |
DO J=1,ny |
811 |
DO I=1,nx |
812 |
divH(I,J)= divH(I,J)+ |
813 |
& XAk(I+1,J)*uK(I+1,J)-XAk(I,J)*uK(I,J) |
814 |
& +YAk(I,J+1)*vK(I,J+1)-YAk(I,J)*vK(I,J) |
815 |
ENDDO |
816 |
ENDDO |
817 |
|
818 |
RETURN |
819 |
END |