1 |
cnh |
1.2 |
C $Id: update_t.F,v 1.1 1998/05/25 20:21:06 cnh Exp $ |
2 |
cnh |
1.1 |
#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 |
cnh |
1.2 |
C IF ( TOP_LAYER ) THEN |
338 |
|
|
C DO J =1,Ny |
339 |
|
|
C DO I=1, Nx |
340 |
cnh |
1.1 |
Cdbg gT(_I3(K,I,J)) = gT(_I3(K,I,J)) + HEAT(I,J) |
341 |
cnh |
1.2 |
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 |
cnh |
1.1 |
|
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 |