/[MITgcm]/MITgcm/pkg/fizhi/fizhi_step_diag.F
ViewVC logotype

Contents of /MITgcm/pkg/fizhi/fizhi_step_diag.F

Parent Directory Parent Directory | Revision Log Revision Log | View Revision Graph Revision Graph


Revision 1.9 - (show annotations) (download)
Tue Aug 10 15:13:47 2004 UTC (19 years, 10 months ago) by molod
Branch: MAIN
Changes since 1.8: +7 -1 lines
Debugging

1 C $Header: /u/gcmpack/MITgcm/pkg/fizhi/fizhi_step_diag.F,v 1.8 2004/07/30 18:55:50 molod Exp $
2 C $Name: $
3
4 #include "FIZHI_OPTIONS.h"
5 subroutine fizhi_step_diag(myThid,p,uphy,vphy,thphy,sphy,qq,pk,dp,
6 . radswt,radswg,swgclr,osr,osrclr,st4,dst4,tgz,tg0,radlwg,lwgclr,
7 . turbu,turbv,turbt,turbq,moistu,moistv,moistt,moistq,
8 . lwdt,swdt,lwdtclr,swdtclr,dlwdtg,
9 . im1,im2,jm1,jm2,Nrphys,Nbi,Nbj,bi,bj,ntracer)
10 C***********************************************************************
11 implicit none
12
13 #ifdef ALLOW_DIAGNOSTICS
14 #include "SIZE.h"
15 #include "diagnostics_SIZE.h"
16 #include "diagnostics.h"
17 #endif
18
19 integer myThid,im1,im2,jm1,jm2,Nrphys,Nbi,Nbj,bi,bj,ntracer
20 _RL p(im2,jm2,Nbi,Nbj)
21 _RL uphy(im2,jm2,Nrphys,Nbi,Nbj)
22 _RL vphy(im2,jm2,Nrphys,Nbi,Nbj)
23 _RL thphy(im2,jm2,Nrphys,Nbi,Nbj)
24 _RL sphy(im2,jm2,Nrphys,Nbi,Nbj)
25 _RL qq(im2,jm2,Nrphys),pk(im2,jm2,Nrphys,Nbi,Nbj)
26 _RL dp(im2,jm2,Nrphys,Nbi,Nbj)
27 _RL radswt(im2,jm2,Nbi,Nbj),radswg(im2,jm2,Nbi,Nbj)
28 _RL swgclr(im2,jm2,Nbi,Nbj),osr(im2,jm2,Nbi,Nbj)
29 _RL osrclr(im2,jm2,Nbi,Nbj),st4(im2,jm2,Nbi,Nbj)
30 _RL dst4(im2,jm2,Nbi,Nbj),tgz(im2,jm2,Nbi,Nbj)
31 _RL tg0(im2,jm2,Nbi,Nbj),radlwg(im2,jm2,Nbi,Nbj)
32 _RL lwgclr(im2,jm2,Nbi,Nbj)
33 _RL turbu(im2,jm2,Nrphys,Nbi,Nbj)
34 _RL turbv(im2,jm2,Nrphys,Nbi,Nbj)
35 _RL turbt(im2,jm2,Nrphys,Nbi,Nbj)
36 _RL turbq(im2,jm2,Nrphys,ntracer,Nbi,Nbj)
37 _RL moistu(im2,jm2,Nrphys,Nbi,Nbj)
38 _RL moistv(im2,jm2,Nrphys,Nbi,Nbj)
39 _RL moistt(im2,jm2,Nrphys,Nbi,Nbj)
40 _RL moistq(im2,jm2,Nrphys,ntracer,Nbi,Nbj)
41 _RL lwdt(im2,jm2,Nrphys,Nbi,Nbj)
42 _RL swdt(im2,jm2,Nrphys,Nbi,Nbj)
43 _RL lwdtclr(im2,jm2,Nrphys,Nbi,Nbj)
44 _RL swdtclr(im2,jm2,Nrphys,Nbi,Nbj)
45 _RL dlwdtg(im2,jm2,Nrphys,Nbi,Nbj)
46
47 integer i,j,L
48 _RL pinv(im2,jm2), qbar(im2,jm2)
49
50 C **********************************************************************
51
52 #ifdef ALLOW_DIAGNOSTICS
53 do j=jm1,jm2
54 do i=im1,im2
55 pinv(i,j) = 1.0 / p(i,j,bi,bj)
56 enddo
57 enddo
58
59 c Incident Solar Radiation (W/m**2)
60 c ---------------------------------
61 if (iradswt.ne.0) then
62 do j=jm1,jm2
63 do i=im1,im2
64 qdiag(i,j,iradswt,bi,bj)= qdiag(i,j,iradswt,bi,bj) +
65 . radswt(i,j,bi,bj)
66 enddo
67 enddo
68 endif
69
70 c Net Solar Radiation at the Ground (W/m**2)
71 c ------------------------------------------
72 if (iradswg.ne.0) then
73 do j=jm1,jm2
74 do i=im1,im2
75 qdiag(i,j,iradswg,bi,bj) = qdiag(i,j,iradswg,bi,bj) +
76 . radswg(i,j,bi,bj)*radswt(i,j,bi,bj)
77 enddo
78 enddo
79 endif
80
81 c Net Clear Sky Solar Radiation at the Ground (W/m**2)
82 c ----------------------------------------------------
83 if (iswgclr.ne.0) then
84 do j=jm1,jm2
85 do i=im1,im2
86 qdiag(i,j,iswgclr,bi,bj) = qdiag(i,j,iswgclr,bi,bj) +
87 . swgclr(i,j,bi,bj)*radswt(i,j,bi,bj)
88 enddo
89 enddo
90 endif
91
92 c Outgoing Solar Radiation at top (W/m**2)
93 c -----------------------------------------
94 if (iosr.ne.0) then
95 do j=jm1,jm2
96 do i=im1,im2
97 qdiag(i,j,iosr,bi,bj) = qdiag(i,j,iosr,bi,bj) +
98 . (1.0-osr(i,j,bi,bj))*radswt(i,j,bi,bj)
99 enddo
100 enddo
101 endif
102
103 c Outgoing Clear Sky Solar Radiation at top (W/m**2)
104 c ---------------------------------------------------
105 if (iosrclr.ne.0) then
106 do j=jm1,jm2
107 do i=im1,im2
108 qdiag(i,j,iosrclr,bi,bj) = qdiag(i,j,iosrclr,bi,bj) +
109 . (1.0-osrclr(i,j,bi,bj))*radswt(i,j,bi,bj)
110 enddo
111 enddo
112 endif
113
114 c Upward Longwave Flux at the Ground (W/m**2)
115 c -------------------------------------------
116 if (ilwgup.ne.0) then
117 do j=jm1,jm2
118 do i=im1,im2
119 qdiag(i,j,ilwgup,bi,bj) = qdiag(i,j,ilwgup,bi,bj) + st4(i,j,bi,bj)
120 . + dst4(i,j,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj))
121 enddo
122 enddo
123 endif
124
125 c Net Longwave Flux at the Ground (W/m**2)
126 c ----------------------------------------
127 if (iradlwg.ne.0) then
128 do j=jm1,jm2
129 do i=im1,im2
130 qdiag(i,j,iradlwg,bi,bj) = qdiag(i,j,iradlwg,bi,bj) +
131 . radlwg(i,j,bi,bj) +
132 . dst4(i,j,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj))
133 enddo
134 enddo
135 endif
136
137 c Net Longwave Flux at the Ground Clear Sky (W/m**2)
138 c --------------------------------------------------
139 if (ilwgclr.ne.0) then
140 do j=jm1,jm2
141 do i=im1,im2
142 qdiag(i,j,ilwgclr,bi,bj) = qdiag(i,j,ilwgclr,bi,bj) +
143 . lwgclr(i,j,bi,bj) +
144 . dst4(i,j,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj))
145 enddo
146 enddo
147 endif
148
149 if( (bi.eq.1) .and. (bj.eq.1) ) then
150 nradswt = nradswt + 1
151 nradswg = nradswg + 1
152 nswgclr = nswgclr + 1
153 nosr = nosr + 1
154 nosrclr = nosrclr + 1
155 nradlwg = nradlwg + 1
156 nlwgclr = nlwgclr + 1
157 nlwgup = nlwgup + 1
158 endif
159
160 C **********************************************************************
161 do L=1,Nrphys
162
163 c Total Diabatic U-Tendency (m/sec/day)
164 c -------------------------------------
165 if( idiabu.ne.0 ) then
166 do j=jm1,jm2
167 do i=im1,im2
168 qdiag(i,j,idiabu+L-1,bi,bj) = qdiag(i,j,idiabu+L-1,bi,bj)
169 . + ( moistu (i,j,L,bi,bj) + turbu(i,j,L,bi,bj) )*86400
170 enddo
171 enddo
172 endif
173
174 c Total Diabatic V-Tendency (m/sec/day)
175 c -------------------------------------
176 if( idiabv.ne.0 ) then
177 do j=jm1,jm2
178 do i=im1,im2
179 qdiag(i,j,idiabv+L-1,bi,bj) = qdiag(i,j,idiabv+L-1,bi,bj)
180 . + ( moistv (i,j,L,bi,bj) + turbv(i,j,L,bi,bj) )*86400
181 enddo
182 enddo
183 endif
184
185 c Total Diabatic T-Tendency (deg/day)
186 c -----------------------------------
187 if( idiabt.ne.0 ) then
188 do j=jm1,jm2
189 do i=im1,im2
190 qdiag(i,j,idiabt+L-1,bi,bj) = qdiag(i,j,idiabt+L-1,bi,bj) +
191 . ( turbt(i,j,L,bi,bj) + moistt(i,j,L,bi,bj) +
192 . lwdt(i,j,L,bi,bj) +
193 . dlwdtg(i,j,L,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj)) +
194 . swdt(i,j,L,bi,bj)*radswt(i,j,bi,bj) )
195 . * pk(i,j,L,bi,bj)*pinv(i,j)*86400
196 enddo
197 enddo
198 endif
199
200 c Total Diabatic Q-Tendency (g/kg/day)
201 c ------------------------------------
202 if( idiabq.ne.0 ) then
203 do j=jm1,jm2
204 do i=im1,im2
205 qdiag(i,j,idiabq+L-1,bi,bj) = qdiag(i,j,idiabq+L-1,bi,bj) +
206 . ( turbq(i,j,L,1,bi,bj) + moistq(i,j,L,1,bi,bj) ) *
207 . pinv(i,j)*86400*1000
208 enddo
209 enddo
210 endif
211
212 c Longwave Heating (deg/day)
213 c --------------------------
214 if (iradlw.ne.0) then
215 do j=jm1,jm2
216 do i=im1,im2
217 qdiag(i,j,iradlw+l-1,bi,bj) = qdiag(i,j,iradlw+l-1,bi,bj) +
218 . ( lwdt(i,j,l,bi,bj) +
219 . dlwdtg (i,j,L,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj)) )
220 . * pk(i,j,l,bi,bj)*pinv(i,j)*86400
221 enddo
222 enddo
223 endif
224
225 c Longwave Heating Clear-Sky (deg/day)
226 c ------------------------------------
227 if (ilwclr.ne.0) then
228 do j=jm1,jm2
229 do i=im1,im2
230 qdiag(i,j,ilwclr+l-1,bi,bj) = qdiag(i,j,ilwclr+l-1,bi,bj) +
231 . ( lwdtclr(i,j,l,bi,bj) +
232 . dlwdtg(i,j,L,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj)) )
233 . * pk(i,j,l,bi,bj)*pinv(i,j)*86400
234 enddo
235 enddo
236 endif
237
238 c Solar Radiative Heating (deg/day)
239 c ---------------------------------
240 if (iradsw.ne.0) then
241 do j=jm1,jm2
242 do i=im1,im2
243 qdiag(i,j,iradsw+l-1,bi,bj) = qdiag(i,j,iradsw+l-1,bi,bj) +
244 . + swdt(i,j,l,bi,bj)*radswt(i,j,bi,bj)*
245 . pk(i,j,l,bi,bj)*pinv(i,j)*86400
246 enddo
247 enddo
248 endif
249
250 c Clear Sky Solar Radiative Heating (deg/day)
251 c -------------------------------------------
252 if (iswclr.ne.0) then
253 do j=jm1,jm2
254 do i=im1,im2
255 qdiag(i,j,iswclr+l-1,bi,bj) = qdiag(i,j,iswclr+l-1,bi,bj) +
256 . swdtclr(i,j,l,bi,bj)*radswt(i,j,bi,bj)*
257 . pk(i,j,l,bi,bj)*pinv(i,j)*86400
258 enddo
259 enddo
260 endif
261
262 c Averaged U-Field (m/sec)
263 c ------------------------
264 if( iuwnd.ne.0 ) then
265 do j=jm1,jm2
266 do i=im1,im2
267 qdiag(i,j,iuwnd+L-1,bi,bj) = qdiag(i,j,iuwnd+L-1,bi,bj) +
268 . uphy(i,j,L,bi,bj)
269 enddo
270 enddo
271 endif
272
273 c Averaged V-Field (m/sec)
274 c ------------------------
275 if( ivwnd.ne.0 ) then
276 do j=jm1,jm2
277 do i=im1,im2
278 qdiag(i,j,ivwnd+L-1,bi,bj) = qdiag(i,j,ivwnd+L-1,bi,bj) +
279 . vphy(i,j,L,bi,bj)
280 enddo
281 enddo
282 endif
283
284 c Averaged T-Field (deg)
285 c ----------------------
286 if( itmpu.ne.0 ) then
287 do j=jm1,jm2
288 do i=im1,im2
289 qdiag(i,j,itmpu+L-1,bi,bj) = qdiag(i,j,itmpu+L-1,bi,bj) +
290 . thphy(i,j,L,bi,bj)*pk(i,j,L,bi,bj)
291 enddo
292 enddo
293 endif
294
295 c Averaged QQ-Field (m/sec)**2
296 c ----------------------------
297 if( itke.ne.0 ) then
298 do j=jm1,jm2
299 do i=im1,im2
300 qdiag(i,j,itke+L-1,bi,bj) = qdiag(i,j,itke+L-1,bi,bj) + qq(i,j,L)
301 enddo
302 enddo
303 endif
304
305 c Averaged Q-Field (g/kg)
306 c -----------------------
307 if( isphu.ne.0 ) then
308 do j=jm1,jm2
309 do i=im1,im2
310 qdiag(i,j,isphu+L-1,bi,bj) = qdiag(i,j,isphu+L-1,bi,bj) +
311 . sphy(i,j,L,bi,bj)*1000
312 enddo
313 enddo
314 endif
315
316 enddo
317
318 if( (bi.eq.1) .and. (bj.eq.1) ) then
319 ndiabu = ndiabu + 1
320 ndiabv = ndiabv + 1
321 ndiabt = ndiabt + 1
322 ndiabq = ndiabq + 1
323 nradlw = nradlw + 1
324 nlwclr = nlwclr + 1
325 nradsw = nradsw + 1
326 nswclr = nswclr + 1
327 nuwnd = nuwnd + 1
328 nvwnd = nvwnd + 1
329 ntmpu = ntmpu + 1
330 ntke = ntke + 1
331 nsphu = nsphu + 1
332 endif
333
334 C **********************************************************************
335
336 c Vertically Averaged Moist-T Increment (K/day)
337 c ---------------------------------------------
338 if( ivdtmoist.ne.0 ) then
339 do j=jm1,jm2
340 do i=im1,im2
341 qbar(i,j) = 0.0
342 enddo
343 enddo
344 do L=1,Nrphys
345 do j=jm1,jm2
346 do i=im1,im2
347 qbar(i,j) = qbar(i,j) +
348 . moistt(i,j,L,bi,bj)*pk(i,j,l,bi,bj)*dp(i,j,L,bi,bj)
349 enddo
350 enddo
351 enddo
352 do j=jm1,jm2
353 do i=im1,im2
354 qdiag(i,j,ivdtmoist,bi,bj) = qdiag(i,j,ivdtmoist,bi,bj) +
355 . qbar(i,j)*pinv(i,j)*pinv(i,j)*86400
356 enddo
357 enddo
358 endif
359
360 c Vertically Averaged Turb-T Increment (K/day)
361 c --------------------------------------------
362 if( ivdtturb.ne.0 ) then
363 do j=jm1,jm2
364 do i=im1,im2
365 qbar(i,j) = 0.0
366 enddo
367 enddo
368 do L=1,Nrphys
369 do j=jm1,jm2
370 do i=im1,im2
371 qbar(i,j) = qbar(i,j) +
372 . turbt(i,j,L,bi,bj)*pk(i,j,l,bi,bj)*dp(i,j,L,bi,bj)
373 enddo
374 enddo
375 enddo
376 do j=jm1,jm2
377 do i=im1,im2
378 qdiag(i,j,ivdtturb,bi,bj) = qdiag(i,j,ivdtturb,bi,bj) +
379 . qbar(i,j)*pinv(i,j)*pinv(i,j)*86400
380 enddo
381 enddo
382 endif
383
384 c Vertically Averaged RADLW Temperature Increment (K/day)
385 c -------------------------------------------------------
386 if( ivdtradlw.ne.0 ) then
387 do j=jm1,jm2
388 do i=im1,im2
389 qbar(i,j) = 0.0
390 enddo
391 enddo
392 do L=1,Nrphys
393 do j=jm1,jm2
394 do i=im1,im2
395 qbar(i,j) = qbar(i,j) + ( lwdt(i,j,L,bi,bj) +
396 . dlwdtg(i,j,L,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj)) )
397 . *pk(i,j,l,bi,bj)*dp(i,j,L,bi,bj)
398 enddo
399 enddo
400 enddo
401 do j=jm1,jm2
402 do i=im1,im2
403 qdiag(i,j,ivdtradlw,bi,bj) = qdiag(i,j,ivdtradlw,bi,bj) +
404 . qbar(i,j)*pinv(i,j)*pinv(i,j)*86400
405 enddo
406 enddo
407 endif
408
409 c Vertically Averaged RADSW Temperature Increment (K/day)
410 c -------------------------------------------------------
411 if( ivdtradsw.ne.0 ) then
412 do j=jm1,jm2
413 do i=im1,im2
414 qbar(i,j) = 0.0
415 enddo
416 enddo
417 do L=1,Nrphys
418 do j=jm1,jm2
419 do i=im1,im2
420 qbar(i,j) = qbar(i,j) +
421 . swdt(i,j,L,bi,bj)*pk(i,j,l,bi,bj)*dp(i,j,L,bi,bj)
422 enddo
423 enddo
424 enddo
425 do j=jm1,jm2
426 do i=im1,im2
427 qdiag(i,j,ivdtradsw,bi,bj) = qdiag(i,j,ivdtradsw,bi,bj) +
428 . qbar(i,j)*radswt(i,j,bi,bj)*pinv(i,j)*pinv(i,j)*86400
429 enddo
430 enddo
431 endif
432
433 if( (bi.eq.1) .and. (bj.eq.1) ) then
434 nvdtmoist = nvdtmoist + 1
435 nvdtturb = nvdtturb + 1
436 nvdtradlw = nvdtradlw + 1
437 nvdtradsw = nvdtradsw + 1
438 endif
439
440 #endif
441 return
442 end

  ViewVC Help
Powered by ViewVC 1.1.22