/[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.10 - (show annotations) (download)
Thu Aug 12 15:21:22 2004 UTC (19 years, 10 months ago) by molod
Branch: MAIN
Changes since 1.9: +21 -1 lines
Debugging - Code for work-around for input vegetation dataset

1 C $Header: /u/gcmpack/MITgcm/pkg/fizhi/fizhi_step_diag.F,v 1.9 2004/08/10 15:13:47 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 if (iudiag2.ne.0) then
225 do j=jm1,jm2
226 do i=im1,im2
227 qdiag(i,j,iudiag2+l-1,bi,bj) = qdiag(i,j,iudiag2+l-1,bi,bj) +
228 . lwdt(i,j,l,bi,bj) * pk(i,j,l,bi,bj)*pinv(i,j)*86400
229 enddo
230 enddo
231 endif
232 if (iudiag3.ne.0) then
233 do j=jm1,jm2
234 do i=im1,im2
235 qdiag(i,j,iudiag3+l-1,bi,bj) = qdiag(i,j,iudiag3+l-1,bi,bj) +
236 . dlwdtg (i,j,L,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj))
237 . * pk(i,j,l,bi,bj)*pinv(i,j)*86400
238 enddo
239 enddo
240 endif
241
242 c Longwave Heating Clear-Sky (deg/day)
243 c ------------------------------------
244 if (ilwclr.ne.0) then
245 do j=jm1,jm2
246 do i=im1,im2
247 qdiag(i,j,ilwclr+l-1,bi,bj) = qdiag(i,j,ilwclr+l-1,bi,bj) +
248 . ( lwdtclr(i,j,l,bi,bj) +
249 . dlwdtg(i,j,L,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj)) )
250 . * pk(i,j,l,bi,bj)*pinv(i,j)*86400
251 enddo
252 enddo
253 endif
254
255 c Solar Radiative Heating (deg/day)
256 c ---------------------------------
257 if (iradsw.ne.0) then
258 do j=jm1,jm2
259 do i=im1,im2
260 qdiag(i,j,iradsw+l-1,bi,bj) = qdiag(i,j,iradsw+l-1,bi,bj) +
261 . + swdt(i,j,l,bi,bj)*radswt(i,j,bi,bj)*
262 . pk(i,j,l,bi,bj)*pinv(i,j)*86400
263 enddo
264 enddo
265 endif
266
267 c Clear Sky Solar Radiative Heating (deg/day)
268 c -------------------------------------------
269 if (iswclr.ne.0) then
270 do j=jm1,jm2
271 do i=im1,im2
272 qdiag(i,j,iswclr+l-1,bi,bj) = qdiag(i,j,iswclr+l-1,bi,bj) +
273 . swdtclr(i,j,l,bi,bj)*radswt(i,j,bi,bj)*
274 . pk(i,j,l,bi,bj)*pinv(i,j)*86400
275 enddo
276 enddo
277 endif
278
279 c Averaged U-Field (m/sec)
280 c ------------------------
281 if( iuwnd.ne.0 ) then
282 do j=jm1,jm2
283 do i=im1,im2
284 qdiag(i,j,iuwnd+L-1,bi,bj) = qdiag(i,j,iuwnd+L-1,bi,bj) +
285 . uphy(i,j,L,bi,bj)
286 enddo
287 enddo
288 endif
289
290 c Averaged V-Field (m/sec)
291 c ------------------------
292 if( ivwnd.ne.0 ) then
293 do j=jm1,jm2
294 do i=im1,im2
295 qdiag(i,j,ivwnd+L-1,bi,bj) = qdiag(i,j,ivwnd+L-1,bi,bj) +
296 . vphy(i,j,L,bi,bj)
297 enddo
298 enddo
299 endif
300
301 c Averaged T-Field (deg)
302 c ----------------------
303 if( itmpu.ne.0 ) then
304 do j=jm1,jm2
305 do i=im1,im2
306 qdiag(i,j,itmpu+L-1,bi,bj) = qdiag(i,j,itmpu+L-1,bi,bj) +
307 . thphy(i,j,L,bi,bj)*pk(i,j,L,bi,bj)
308 enddo
309 enddo
310 endif
311
312 c Averaged QQ-Field (m/sec)**2
313 c ----------------------------
314 if( itke.ne.0 ) then
315 do j=jm1,jm2
316 do i=im1,im2
317 qdiag(i,j,itke+L-1,bi,bj) = qdiag(i,j,itke+L-1,bi,bj) + qq(i,j,L)
318 enddo
319 enddo
320 endif
321
322 c Averaged Q-Field (g/kg)
323 c -----------------------
324 if( isphu.ne.0 ) then
325 do j=jm1,jm2
326 do i=im1,im2
327 qdiag(i,j,isphu+L-1,bi,bj) = qdiag(i,j,isphu+L-1,bi,bj) +
328 . sphy(i,j,L,bi,bj)*1000
329 enddo
330 enddo
331 endif
332
333 enddo
334
335 if( (bi.eq.1) .and. (bj.eq.1) ) then
336 ndiabu = ndiabu + 1
337 ndiabv = ndiabv + 1
338 ndiabt = ndiabt + 1
339 ndiabq = ndiabq + 1
340 nradlw = nradlw + 1
341 nlwclr = nlwclr + 1
342 nradsw = nradsw + 1
343 nswclr = nswclr + 1
344 nuwnd = nuwnd + 1
345 nvwnd = nvwnd + 1
346 ntmpu = ntmpu + 1
347 ntke = ntke + 1
348 nsphu = nsphu + 1
349
350 nudiag3 = nudiag3 + 1
351 nudiag2 = nudiag2 + 1
352 endif
353
354 C **********************************************************************
355
356 c Vertically Averaged Moist-T Increment (K/day)
357 c ---------------------------------------------
358 if( ivdtmoist.ne.0 ) then
359 do j=jm1,jm2
360 do i=im1,im2
361 qbar(i,j) = 0.0
362 enddo
363 enddo
364 do L=1,Nrphys
365 do j=jm1,jm2
366 do i=im1,im2
367 qbar(i,j) = qbar(i,j) +
368 . moistt(i,j,L,bi,bj)*pk(i,j,l,bi,bj)*dp(i,j,L,bi,bj)
369 enddo
370 enddo
371 enddo
372 do j=jm1,jm2
373 do i=im1,im2
374 qdiag(i,j,ivdtmoist,bi,bj) = qdiag(i,j,ivdtmoist,bi,bj) +
375 . qbar(i,j)*pinv(i,j)*pinv(i,j)*86400
376 enddo
377 enddo
378 endif
379
380 c Vertically Averaged Turb-T Increment (K/day)
381 c --------------------------------------------
382 if( ivdtturb.ne.0 ) then
383 do j=jm1,jm2
384 do i=im1,im2
385 qbar(i,j) = 0.0
386 enddo
387 enddo
388 do L=1,Nrphys
389 do j=jm1,jm2
390 do i=im1,im2
391 qbar(i,j) = qbar(i,j) +
392 . turbt(i,j,L,bi,bj)*pk(i,j,l,bi,bj)*dp(i,j,L,bi,bj)
393 enddo
394 enddo
395 enddo
396 do j=jm1,jm2
397 do i=im1,im2
398 qdiag(i,j,ivdtturb,bi,bj) = qdiag(i,j,ivdtturb,bi,bj) +
399 . qbar(i,j)*pinv(i,j)*pinv(i,j)*86400
400 enddo
401 enddo
402 endif
403
404 c Vertically Averaged RADLW Temperature Increment (K/day)
405 c -------------------------------------------------------
406 if( ivdtradlw.ne.0 ) then
407 do j=jm1,jm2
408 do i=im1,im2
409 qbar(i,j) = 0.0
410 enddo
411 enddo
412 do L=1,Nrphys
413 do j=jm1,jm2
414 do i=im1,im2
415 qbar(i,j) = qbar(i,j) + ( lwdt(i,j,L,bi,bj) +
416 . dlwdtg(i,j,L,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj)) )
417 . *pk(i,j,l,bi,bj)*dp(i,j,L,bi,bj)
418 enddo
419 enddo
420 enddo
421 do j=jm1,jm2
422 do i=im1,im2
423 qdiag(i,j,ivdtradlw,bi,bj) = qdiag(i,j,ivdtradlw,bi,bj) +
424 . qbar(i,j)*pinv(i,j)*pinv(i,j)*86400
425 enddo
426 enddo
427 endif
428
429 c Vertically Averaged RADSW Temperature Increment (K/day)
430 c -------------------------------------------------------
431 if( ivdtradsw.ne.0 ) then
432 do j=jm1,jm2
433 do i=im1,im2
434 qbar(i,j) = 0.0
435 enddo
436 enddo
437 do L=1,Nrphys
438 do j=jm1,jm2
439 do i=im1,im2
440 qbar(i,j) = qbar(i,j) +
441 . swdt(i,j,L,bi,bj)*pk(i,j,l,bi,bj)*dp(i,j,L,bi,bj)
442 enddo
443 enddo
444 enddo
445 do j=jm1,jm2
446 do i=im1,im2
447 qdiag(i,j,ivdtradsw,bi,bj) = qdiag(i,j,ivdtradsw,bi,bj) +
448 . qbar(i,j)*radswt(i,j,bi,bj)*pinv(i,j)*pinv(i,j)*86400
449 enddo
450 enddo
451 endif
452
453 if( (bi.eq.1) .and. (bj.eq.1) ) then
454 nvdtmoist = nvdtmoist + 1
455 nvdtturb = nvdtturb + 1
456 nvdtradlw = nvdtradlw + 1
457 nvdtradsw = nvdtradsw + 1
458 endif
459
460 #endif
461 return
462 end

  ViewVC Help
Powered by ViewVC 1.1.22