1 |
C $Id: update_t.F,v 1.1 1998/05/25 20:21:06 cnh Exp $ |
2 |
#include "CPP_OPTIONS.h" |
3 |
#include "CPP_MACROS.h" |
4 |
|
5 |
#define _DB0 0. |
6 |
|
7 |
C=================================================================== |
8 |
C Procedure name: UPDATE_T | |
9 |
C Function: Step forward the temperature field. | |
10 |
C Comments: | |
11 |
C=================================================================== |
12 |
CStartofinterface |
13 |
SUBROUTINE UPDATE_T( |
14 |
& U, V, W, K, T, GT ) |
15 |
IMPLICIT NONE |
16 |
C============== Global data ========================================== |
17 |
#include "SIZE.h" |
18 |
#include "AJAINF.h" |
19 |
#include "GRID.h" |
20 |
#include "PARAMS.h" |
21 |
#include "OPERATORS.h" |
22 |
#include "OLDG.h" |
23 |
#include "MASKS.h" |
24 |
#include "FORCING.h" |
25 |
C============= Routine arguments ===================================== |
26 |
C /--------------------------------------------------------------\ |
27 |
C | U, V, W - X,Y,Z Velocity ( m/s, m/s, Pa/s ). | |
28 |
C | K - Working level. | |
29 |
C | T - Potential temperature (oC). | |
30 |
C \--------------------------------------------------------------/ |
31 |
REAL U (Nx,Ny,Nz) |
32 |
REAL V (Nx,Ny,Nz) |
33 |
REAL W (Nx,Ny,Nz) |
34 |
INTEGER K |
35 |
REAL T (Nx,Ny,Nz) |
36 |
REAL gt (Nx,Ny,Nz) |
37 |
CEndofinterface |
38 |
C============= Local variables ====================================== |
39 |
C /--------------------------------------------------------------\ |
40 |
C |Slice Notation | |
41 |
C |phiK - XY slice of variable "phi" at level K. | |
42 |
C |phiKP1 - XY slice of variable "phi" at level K+1. | |
43 |
C |phiKM1 - XY slice of variable "phi" at level K-1. | |
44 |
C |phiNM1 - Variable "phi" at time level N-1. | |
45 |
C |phiU - Varaible "phi" interpolated to U grid. | |
46 |
C |phiV - Varaible "phi" interpolated to V grid. | |
47 |
C \--------------------------------------------------------------/ |
48 |
C /--------------------------------------------------------------\ |
49 |
C |Geometry terms | |
50 |
C | o xa,ya,za - X, Y, Z face areas ( m.Pa, m.Pa, m**2 ). | |
51 |
C | o rV - 1/Volume ( m**2.Pa ) | |
52 |
C | o ?Grad?Op - Generic gradient operator A[XYZ][t].d/d[xyz] | |
53 |
C | tGradxOp: AXt/DXt, tGradyOp: AYt/DYt etc... | |
54 |
C \--------------------------------------------------------------/ |
55 |
REAL xaK (Nx+1,0:Ny) |
56 |
REAL yaK (0:Nx,Ny+1) |
57 |
REAL tGradxOp (Nx ,Ny ) |
58 |
REAL tGradyOp (Nx ,Ny ) |
59 |
REAL tGradzOpK (Nx ,Ny ) |
60 |
REAL tGradzOpKP1(Nx ,Ny ) |
61 |
REAL zaK (0:Nx,0:Ny) |
62 |
REAL zaKP1 (Nx,Ny) |
63 |
REAL rVk (Nx ,Ny ) |
64 |
C /--------------------------------------------------------------\ |
65 |
C |State variables | |
66 |
C | o u, v, w - X,Y,Z velocity slices. | |
67 |
C | o t - Potential temperature. | |
68 |
C | o ph - Hydrostatic pressure. | |
69 |
C \--------------------------------------------------------------/ |
70 |
REAL uK (Nx+1,0:Ny) |
71 |
REAL wK (0:Nx,0:Ny) |
72 |
REAL wkP1 (0:Nx,0:Ny) |
73 |
REAL vK (0:Nx,Ny+1) |
74 |
REAL tK (0:Nx+1,0:Ny+1) |
75 |
REAL phK (60) |
76 |
C /--------------------------------------------------------------\ |
77 |
C |Work space arrays | |
78 |
C \--------------------------------------------------------------/ |
79 |
REAL tmp (0:Nx+1,0:Ny+1) |
80 |
REAL tmp2 (0:Nx+1,0:Ny+1) |
81 |
C /--------------------------------------------------------------\ |
82 |
C |Arrays for accumulating fluxes | |
83 |
C | o fX, fY, fZ - Flux in X, Y and Z. | |
84 |
C \--------------------------------------------------------------/ |
85 |
REAL fX (0:Nx+1,Ny) |
86 |
REAL fY (Nx,0:Ny+1) |
87 |
REAL fZK (Nx,Ny) |
88 |
REAL fZKP1 (Nx,Ny) |
89 |
C Loop counters: |
90 |
INTEGER I, J |
91 |
C Flags: |
92 |
LOGICAL TOP_LAYER |
93 |
LOGICAL BOTTOM_LAYER |
94 |
C === Set up controlling flags ======================================== |
95 |
TOP_LAYER = K .EQ. 1 |
96 |
BOTTOM_LAYER = K .EQ. Nz |
97 |
C _D(( ' K, Nz ', K, Nz )) |
98 |
C _D(( ' BOTTOM_LAYER ', BOTTOM_LAYER )) |
99 |
C === Set up the slices =============================================== |
100 |
C /---------------------------------------------------------------\ |
101 |
C | xaK <- xa(K) | |
102 |
C | yaK <- ya(K) | |
103 |
C | za[K,KP1] <- za[(K),(K+1)] | |
104 |
C \---------------------------------------------------------------/ |
105 |
DO J=1,Ny |
106 |
DO I=1,Nx |
107 |
xaK (I,J) = xa (_I3(K,I,J)) |
108 |
yaK (I,J) = ya (_I3(K,I,J)) |
109 |
zaK (I,J) = za (_I3(K,I,J)) |
110 |
ENDDO |
111 |
ENDDO |
112 |
IF ( .NOT. BOTTOM_LAYER ) THEN |
113 |
DO J=1,Ny |
114 |
DO I=1,Nx |
115 |
zaKP1(I,J) = za(_I3(K+1,I,J)) |
116 |
ENDDO |
117 |
ENDDO |
118 |
ENDIF |
119 |
fZK = 0. |
120 |
fZKP1 = 0. |
121 |
C /---------------------------------------------------------------\ |
122 |
C | tGradXOp <- 1/DX*xa | |
123 |
C | tGradYOp <- 1/DY*ya | |
124 |
C | tGradZOp[K,KP1] <- 1/DZ*ZA[K,K+1] | |
125 |
C \---------------------------------------------------------------/ |
126 |
DO J = 1, Ny |
127 |
DO I = 1, Nx |
128 |
tGradxOp (I,J) = xaK(I,J)*pDdxOp(_I3(K,I,J)) |
129 |
tGradyOp (I,J) = yaK(I,J)*pDdyOp(_I3(K,I,J)) |
130 |
tGradzOpK(I,J) = zaK(I,J)*pDdzOp(_I3(K,I,J)) |
131 |
ENDDO |
132 |
ENDDO |
133 |
IF ( .NOT. BOTTOM_LAYER ) THEN |
134 |
DO J = 1, Ny |
135 |
DO I = 1, Nx |
136 |
tGradzOpKP1(I,J) = zaKP1(I,J)*pDdzOp(I,J,K+1) |
137 |
ENDDO |
138 |
ENDDO |
139 |
ENDIF |
140 |
C /---------------------------------------------------------------\ |
141 |
C | w[K,KP1] <- W[(K),(K+1)] | |
142 |
C \---------------------------------------------------------------/ |
143 |
wKP1 = 0. |
144 |
DO J=1,Ny |
145 |
DO I=1,Nx |
146 |
wK(I,J) = W(_I3(K,I,J)) |
147 |
ENDDO |
148 |
ENDDO |
149 |
IF ( .NOT. BOTTOM_LAYER ) THEN |
150 |
DO J=1,Ny |
151 |
DO I=1,Nx |
152 |
wKP1(I,J) = W(_I3(K+1,I,J)) |
153 |
ENDDO |
154 |
ENDDO |
155 |
ENDIF |
156 |
C /---------------------------------------------------------------\ |
157 |
C | uK <- U(K) | |
158 |
C \---------------------------------------------------------------/ |
159 |
DO J=1,Ny |
160 |
DO I=1,Nx |
161 |
uK(I,J) = U(_I3(K,I,J)) |
162 |
ENDDO |
163 |
ENDDO |
164 |
C /---------------------------------------------------------------\ |
165 |
C | vK <- V(K) | |
166 |
C \---------------------------------------------------------------/ |
167 |
DO J=1,Ny |
168 |
DO I=1,Nx |
169 |
vK(I,J) = V(_I3(K,I,J)) |
170 |
ENDDO |
171 |
ENDDO |
172 |
C /---------------------------------------------------------------\ |
173 |
C | tK <- T(K) | |
174 |
C \---------------------------------------------------------------/ |
175 |
DO J=1,Ny |
176 |
DO I=1,Nx |
177 |
tK(I,J) = T(_I3(K,I,J)) |
178 |
ENDDO |
179 |
ENDDO |
180 |
DO J=1,Ny |
181 |
tK(0,J) = tK(Nx,J) |
182 |
tK(Nx+1,J) = tK(1 ,J) |
183 |
ENDDO |
184 |
DO I=0,Nx+1 |
185 |
tK(I,0) = tK(I , Ny) |
186 |
tK(I,Ny+1) = tK(I , 1) |
187 |
ENDDO |
188 |
C /---------------------------------------------------------------\ |
189 |
C | Gt Section | |
190 |
C |===============================================================| |
191 |
C | In this section Gt = GtDiffusion | |
192 |
C | + GtAdvection | |
193 |
C | + GtForcing | |
194 |
C \---------------------------------------------------------------/ |
195 |
C /---------------------------------------------------------------\ |
196 |
C | Initialise local terms | |
197 |
C | rVk <- 1/V | |
198 |
C \---------------------------------------------------------------/ |
199 |
DO J=1,Ny |
200 |
DO I=1,Nx |
201 |
rVk(I,J) = rPvol(_I3(K,I,J)) |
202 |
ENDDO |
203 |
ENDDO |
204 |
fZK = 0. |
205 |
fZKP1 = 0. |
206 |
C /---------------------------------------------------------------\ |
207 |
C | XY GtDiffusion + XY GtAdvection | |
208 |
C | Note: XY GtDiffusion = Laplacian + Biharmonic | |
209 |
C | o Laplacian ( -a2TempXY.del^2(u) ) | |
210 |
C | o Biharmonic ( +a4TempXY.del^4(u) ) | |
211 |
C \---------------------------------------------------------------/ |
212 |
DO J = 1, Ny |
213 |
DO I = 1, Nx |
214 |
C /-------------------------------------------------------------\ |
215 |
C | fX <- 1/DX*XA*Ddx{T} | |
216 |
C | fY <- 1/DY*YA*Ddy{T} | |
217 |
C \-------------------------------------------------------------/ |
218 |
fX(I,J) = tGradxOp(I,J)*(tK(I,J)-tK(I-1,J)) |
219 |
fY(I,J) = tGradyOp(I,J)*(tK(I,J)-tK(I,J-1)) |
220 |
ENDDO |
221 |
ENDDO |
222 |
DO J=1,Ny |
223 |
fX(Nx+1,J) = fX(1 ,J) |
224 |
ENDDO |
225 |
DO I=1,Nx |
226 |
fY(I,Ny+1) = fY(I, 1) |
227 |
ENDDO |
228 |
DO J=1,Ny |
229 |
DO I=1,Nx |
230 |
C /-------------------------------------------------------------\ |
231 |
C | tmp <- 1/V*(Dx{fX}+Dy{fY}) | |
232 |
C \-------------------------------------------------------------/ |
233 |
tmp(I,J)= rVk(I,J)*fX(I+1,J)-rVk(I,J)*fX(I,J) |
234 |
& +rVk(I,J)*fY(I,J+1)-rVk(I,J)*fY(I,J) |
235 |
ENDDO |
236 |
ENDDO |
237 |
tmp(0,1:Ny) = tmp(Nx,1:Ny) |
238 |
tmp(1:Nx,0) = tmp(1:Nx,Ny) |
239 |
DO J = 1,Ny |
240 |
DO I = 1,Nx |
241 |
C /-------------------------------------------------------------\ |
242 |
C | fX <- -k2*fX + ( Laplacian,<Line 1> ) | |
243 |
C | k4/DX*XA*Ddx{tmp} + ( Biharmonic,<Line 2> ) | |
244 |
C | XA*UBx{T} ( Advection,<Line 3> ) | |
245 |
C \-------------------------------------------------------------/ |
246 |
fX(I,J) = 0. |
247 |
_LPT(& -a2TempXY*fX(I,J) ) |
248 |
_BHT(& +a4TempXY*tGradxOp(I,J)*(tmp(I,J)-tmp(I-1,J)) ) |
249 |
_ADT(& +0.5*uK(I,J)*XAK(I,J)*( tK(I-1,J)+tK(I,J) ) ) |
250 |
C /-------------------------------------------------------------\ |
251 |
C | fY <- -k2*fY + ( Laplacian,<Line 1> ) | |
252 |
C | k4/DY*YA*Ddy{tmp} + ( Biharmonic,<Line 2> ) | |
253 |
C | YA*V*By{T} ( Advection,<Line 3> ) | |
254 |
C \-------------------------------------------------------------/ |
255 |
fY(I,J) = 0. |
256 |
_LPT(& -a2TempXY*fY(I,J) ) |
257 |
_BHT(& +a4TempXY*tGradyOp(I,J)*(tmp(I,J)-tmp(I,J-1)) ) |
258 |
_ADT(& +0.5*vK(I,J)*yaK(I,J)*( tK(I,J-1)+tK(I,J) ) ) |
259 |
ENDDO |
260 |
ENDDO |
261 |
C /---------------------------------------------------------------\ |
262 |
C | Z GtDiffusion + Z GtAdvection | |
263 |
C \---------------------------------------------------------------/ |
264 |
IF ( .NOT. TOP_LAYER ) THEN |
265 |
DO J = 1, Ny |
266 |
DO I = 1, Nx |
267 |
C /------------------------------------------------------------\ |
268 |
C | fZK <- -k2Z/DZ*ZA*Ddz{T} + ( Laplacian,<Line 1> ) | |
269 |
C | ZA*-w*Bz{T} ( Advection,<Line 2> ) | |
270 |
C \------------------------------------------------------------/ |
271 |
fZK(I,J) = 0. |
272 |
_LPT(& -a2TempZ*tGradzOpK(I,J)*(T(_I3(K-1,I,J))-tK(I,J) ) ) |
273 |
_ADT(& +0.5*(tK(I,J)+T(_I3(K-1,I,J)))*( -wK(I,J)*ZAK(I,J) ) ) |
274 |
ENDDO |
275 |
ENDDO |
276 |
ELSE |
277 |
DO J = 1, Ny |
278 |
DO I = 1, Nx |
279 |
C /------------------------------------------------------------\ |
280 |
C | fZK <- ZA*-w*Bz{T} ( Advection,<Line 1> ) | |
281 |
C \------------------------------------------------------------/ |
282 |
fZK(I,J) = 0. |
283 |
_ADT(& +0.5*freeSurfFac*(tK(I,J)+tK(I,J))*(-wK(I,J)*ZAK(I,J) ) ) |
284 |
ENDDO |
285 |
ENDDO |
286 |
ENDIF |
287 |
IF ( .NOT. BOTTOM_LAYER ) THEN |
288 |
DO J = 1, Ny |
289 |
DO I = 1, Nx |
290 |
C /------------------------------------------------------------\ |
291 |
C | fZKP1 <- -k2Z/DZ*ZA*Ddz{T} + ( Laplacian,<Line 1> ) | |
292 |
C | ZA*-w*Bz{T} ( Advection,<Line 2> ) | |
293 |
C \------------------------------------------------------------/ |
294 |
fZKP1(I,J) = 0. |
295 |
_LPT(& -a2TempZ*tGradzOpKP1(I,J)*( tK(I,J)-T(_I3(K+1,I,J)) ) ) |
296 |
_ADT(& +0.5*(tK(I,J)+T(_I3(K+1,I,J)))*(-wKP1(I,J))*zaKP1(I,J) ) |
297 |
ENDDO |
298 |
ENDDO |
299 |
ENDIF |
300 |
C /---------------------------------------------------------------\ |
301 |
C | gT <- GtDiffusion + GtAdvection = -div(fX,fY,fZ) | |
302 |
C \---------------------------------------------------------------/ |
303 |
fX(Nx+1,1:Ny) = fX(1,1:Ny) |
304 |
fY(:,Ny+1) = fY(:,1) |
305 |
DO J = 1, Ny |
306 |
DO I = 1, Nx |
307 |
C /-------------------------------------------------------------\ |
308 |
C | gT <- -1/V*(Ddx{fX}+Ddy{fY}) | |
309 |
C \-------------------------------------------------------------/ |
310 |
gt(_I3(K,I,J))=-fX(I+1,J)*rVk(I,J)+fX(I,J)*rVk(I,J) |
311 |
& -fY(I,J+1)*rVk(I,J)+fY(I,J)*rVk(I,J) |
312 |
ENDDO |
313 |
ENDDO |
314 |
DO J = 1, Ny |
315 |
DO I = 1, Nx |
316 |
C /------------------------------------------------------------\ |
317 |
C | gT <- gT-1/Vu*fZK | |
318 |
C \------------------------------------------------------------/ |
319 |
gt(_I3(K,I,J))=gt(_I3(K,I,J))-fZK(I,J)*rVk(I,J) |
320 |
ENDDO |
321 |
ENDDO |
322 |
IF ( .NOT. BOTTOM_LAYER ) THEN |
323 |
DO J = 1, Ny |
324 |
DO I = 1, Nx |
325 |
C /------------------------------------------------------------\ |
326 |
C | gT <- gT+1/Vu*fZP1 | |
327 |
C \------------------------------------------------------------/ |
328 |
gt(_I3(K,I,J))=gt(_I3(K,I,J))+fZKP1(I,J)*rVk(I,J) |
329 |
ENDDO |
330 |
ENDDO |
331 |
ENDIF |
332 |
C _D(( ' S/R UPDATE_T (GtDiff+GtAdv): MAXVAL(Gt )', MAXVAL(Gt))) |
333 |
C _D(( ' S/R UPDATE_T (GtDiff+GtAdv): MINVAL(Gt )', MINVAL(Gt))) |
334 |
C /---------------------------------------------------------------\ |
335 |
C | gT <- gT + GtForcing | |
336 |
C \---------------------------------------------------------------/ |
337 |
C IF ( TOP_LAYER ) THEN |
338 |
C DO J =1,Ny |
339 |
C DO I=1, Nx |
340 |
Cdbg gT(_I3(K,I,J)) = gT(_I3(K,I,J)) + HEAT(I,J) |
341 |
C gT(_I3(K,I,J)) = gT(_I3(K,I,J)) - |
342 |
C & 1./(30.*oneDay)*( T(_I3(K,I,J)) - Tsurf(I,J,K) ) |
343 |
C ENDDO |
344 |
C ENDDO |
345 |
C ENDIF |
346 |
|
347 |
C Final time derivative from AB2. |
348 |
gt(_I3(K,1:Nx,1:Ny)) = Gt(_I3(K,1:Nx,1:Ny))*pMask(_I3(K,:,:)) |
349 |
tmp = 0. |
350 |
tmp2 = 0. |
351 |
tmp (1:Nx,1:Ny) = gtNM1(_I3(K,:,:)) |
352 |
tmp2(1:Nx,1:Ny) = gt (_I3(K,:,:)) |
353 |
|
354 |
C _D(( ' S/R UPDATE_T (Gt Section): MAXVAL(GtNM1)', MAXVAL(GtNm1) )) |
355 |
C _D(( ' S/R UPDATE_T (Gt Section): MAXVAL(Gt)', MAXVAL(Gt) )) |
356 |
|
357 |
gtNM1(_I3(K,:,:)) = gt(_I3(K,:,:)) |
358 |
tmp2=(1.5+abEps)*tmp2-(0.5+abEps)*tmp |
359 |
gt(_I3(K,1:Nx,1:Ny))=tmp2(1:Nx,1:Ny) |
360 |
C |
361 |
RETURN |
362 |
END |