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

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

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


Revision 1.10 - (hide 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 molod 1.10 C $Header: /u/gcmpack/MITgcm/pkg/fizhi/fizhi_step_diag.F,v 1.9 2004/08/10 15:13:47 molod Exp $
2 molod 1.1 C $Name: $
3 molod 1.4
4 molod 1.6 #include "FIZHI_OPTIONS.h"
5 molod 1.2 subroutine fizhi_step_diag(myThid,p,uphy,vphy,thphy,sphy,qq,pk,dp,
6 molod 1.1 . radswt,radswg,swgclr,osr,osrclr,st4,dst4,tgz,tg0,radlwg,lwgclr,
7     . turbu,turbv,turbt,turbq,moistu,moistv,moistt,moistq,
8 molod 1.5 . lwdt,swdt,lwdtclr,swdtclr,dlwdtg,
9     . im1,im2,jm1,jm2,Nrphys,Nbi,Nbj,bi,bj,ntracer)
10 molod 1.1 C***********************************************************************
11     implicit none
12    
13 molod 1.3 #ifdef ALLOW_DIAGNOSTICS
14 molod 1.4 #include "SIZE.h"
15     #include "diagnostics_SIZE.h"
16 molod 1.1 #include "diagnostics.h"
17 molod 1.3 #endif
18 molod 1.1
19 molod 1.5 integer myThid,im1,im2,jm1,jm2,Nrphys,Nbi,Nbj,bi,bj,ntracer
20 molod 1.6 _RL p(im2,jm2,Nbi,Nbj)
21 molod 1.7 _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 molod 1.6 _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 molod 1.7 _RL turbu(im2,jm2,Nrphys,Nbi,Nbj)
34     _RL turbv(im2,jm2,Nrphys,Nbi,Nbj)
35 molod 1.6 _RL turbt(im2,jm2,Nrphys,Nbi,Nbj)
36     _RL turbq(im2,jm2,Nrphys,ntracer,Nbi,Nbj)
37 molod 1.7 _RL moistu(im2,jm2,Nrphys,Nbi,Nbj)
38     _RL moistv(im2,jm2,Nrphys,Nbi,Nbj)
39 molod 1.6 _RL moistt(im2,jm2,Nrphys,Nbi,Nbj)
40     _RL moistq(im2,jm2,Nrphys,ntracer,Nbi,Nbj)
41 molod 1.7 _RL lwdt(im2,jm2,Nrphys,Nbi,Nbj)
42     _RL swdt(im2,jm2,Nrphys,Nbi,Nbj)
43 molod 1.6 _RL lwdtclr(im2,jm2,Nrphys,Nbi,Nbj)
44     _RL swdtclr(im2,jm2,Nrphys,Nbi,Nbj)
45     _RL dlwdtg(im2,jm2,Nrphys,Nbi,Nbj)
46 molod 1.1
47 molod 1.2 integer i,j,L
48 molod 1.6 _RL pinv(im2,jm2), qbar(im2,jm2)
49 molod 1.1
50     C **********************************************************************
51    
52 molod 1.8 #ifdef ALLOW_DIAGNOSTICS
53 molod 1.1 do j=jm1,jm2
54     do i=im1,im2
55 molod 1.2 pinv(i,j) = 1.0 / p(i,j,bi,bj)
56 molod 1.1 enddo
57     enddo
58 molod 1.2
59 molod 1.1 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 molod 1.2 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 molod 1.1
70     c Net Solar Radiation at the Ground (W/m**2)
71     c ------------------------------------------
72 molod 1.2 if (iradswg.ne.0) then
73 molod 1.1 do j=jm1,jm2
74     do i=im1,im2
75 molod 1.2 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 molod 1.1
81     c Net Clear Sky Solar Radiation at the Ground (W/m**2)
82     c ----------------------------------------------------
83 molod 1.2 if (iswgclr.ne.0) then
84 molod 1.1 do j=jm1,jm2
85     do i=im1,im2
86 molod 1.2 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 molod 1.1
92 molod 1.2 c Outgoing Solar Radiation at top (W/m**2)
93 molod 1.1 c -----------------------------------------
94 molod 1.2 if (iosr.ne.0) then
95 molod 1.1 do j=jm1,jm2
96     do i=im1,im2
97 molod 1.2 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 molod 1.1
103 molod 1.2 c Outgoing Clear Sky Solar Radiation at top (W/m**2)
104 molod 1.1 c ---------------------------------------------------
105 molod 1.2 if (iosrclr.ne.0) then
106 molod 1.1 do j=jm1,jm2
107     do i=im1,im2
108 molod 1.2 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 molod 1.1
114     c Upward Longwave Flux at the Ground (W/m**2)
115     c -------------------------------------------
116 molod 1.2 if (ilwgup.ne.0) then
117 molod 1.1 do j=jm1,jm2
118     do i=im1,im2
119 molod 1.2 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 molod 1.1
125     c Net Longwave Flux at the Ground (W/m**2)
126     c ----------------------------------------
127 molod 1.2 if (iradlwg.ne.0) then
128 molod 1.1 do j=jm1,jm2
129     do i=im1,im2
130 molod 1.2 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 molod 1.1
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 molod 1.2 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 molod 1.1 enddo
146     enddo
147     endif
148 molod 1.2
149 molod 1.9 if( (bi.eq.1) .and. (bj.eq.1) ) then
150 molod 1.2 nradswt = nradswt + 1
151     nradswg = nradswg + 1
152     nswgclr = nswgclr + 1
153 molod 1.1 nosr = nosr + 1
154     nosrclr = nosrclr + 1
155     nradlwg = nradlwg + 1
156     nlwgclr = nlwgclr + 1
157     nlwgup = nlwgup + 1
158 molod 1.9 endif
159 molod 1.1
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 molod 1.2 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 molod 1.1 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 molod 1.2 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 molod 1.1 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 molod 1.2 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 molod 1.1 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 molod 1.2 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 molod 1.1 enddo
209     enddo
210     endif
211    
212 molod 1.2 c Longwave Heating (deg/day)
213     c --------------------------
214     if (iradlw.ne.0) then
215 molod 1.1 do j=jm1,jm2
216     do i=im1,im2
217 molod 1.2 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 molod 1.1 enddo
222     enddo
223     endif
224 molod 1.10 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 molod 1.1
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 molod 1.2 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 molod 1.1
255     c Solar Radiative Heating (deg/day)
256     c ---------------------------------
257 molod 1.2 if (iradsw.ne.0) then
258 molod 1.1 do j=jm1,jm2
259     do i=im1,im2
260 molod 1.2 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 molod 1.1
267     c Clear Sky Solar Radiative Heating (deg/day)
268     c -------------------------------------------
269 molod 1.2 if (iswclr.ne.0) then
270 molod 1.1 do j=jm1,jm2
271     do i=im1,im2
272 molod 1.2 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 molod 1.5 . pk(i,j,l,bi,bj)*pinv(i,j)*86400
275 molod 1.1 enddo
276     enddo
277     endif
278 molod 1.2
279 molod 1.1 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 molod 1.2 qdiag(i,j,iuwnd+L-1,bi,bj) = qdiag(i,j,iuwnd+L-1,bi,bj) +
285     . uphy(i,j,L,bi,bj)
286 molod 1.1 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 molod 1.2 qdiag(i,j,ivwnd+L-1,bi,bj) = qdiag(i,j,ivwnd+L-1,bi,bj) +
296     . vphy(i,j,L,bi,bj)
297 molod 1.1 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 molod 1.2 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 molod 1.1 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 molod 1.2 qdiag(i,j,itke+L-1,bi,bj) = qdiag(i,j,itke+L-1,bi,bj) + qq(i,j,L)
318 molod 1.1 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 molod 1.2 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 molod 1.1 enddo
330     enddo
331     endif
332    
333 molod 1.2 enddo
334 molod 1.1
335 molod 1.9 if( (bi.eq.1) .and. (bj.eq.1) ) then
336 molod 1.1 ndiabu = ndiabu + 1
337     ndiabv = ndiabv + 1
338     ndiabt = ndiabt + 1
339 molod 1.2 ndiabq = ndiabq + 1
340 molod 1.1 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 molod 1.10
350     nudiag3 = nudiag3 + 1
351     nudiag2 = nudiag2 + 1
352 molod 1.9 endif
353 molod 1.1
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 molod 1.2 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 molod 1.1 enddo
370     enddo
371     enddo
372     do j=jm1,jm2
373     do i=im1,im2
374 molod 1.2 qdiag(i,j,ivdtmoist,bi,bj) = qdiag(i,j,ivdtmoist,bi,bj) +
375 molod 1.5 . qbar(i,j)*pinv(i,j)*pinv(i,j)*86400
376 molod 1.1 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 molod 1.2 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 molod 1.1 enddo
394     enddo
395     enddo
396     do j=jm1,jm2
397     do i=im1,im2
398 molod 1.2 qdiag(i,j,ivdtturb,bi,bj) = qdiag(i,j,ivdtturb,bi,bj) +
399 molod 1.5 . qbar(i,j)*pinv(i,j)*pinv(i,j)*86400
400 molod 1.1 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 molod 1.2 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 molod 1.1 enddo
419     enddo
420     enddo
421     do j=jm1,jm2
422     do i=im1,im2
423 molod 1.2 qdiag(i,j,ivdtradlw,bi,bj) = qdiag(i,j,ivdtradlw,bi,bj) +
424 molod 1.5 . qbar(i,j)*pinv(i,j)*pinv(i,j)*86400
425 molod 1.1 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 molod 1.2 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 molod 1.1 enddo
443     enddo
444     enddo
445     do j=jm1,jm2
446     do i=im1,im2
447 molod 1.2 qdiag(i,j,ivdtradsw,bi,bj) = qdiag(i,j,ivdtradsw,bi,bj) +
448 molod 1.5 . qbar(i,j)*radswt(i,j,bi,bj)*pinv(i,j)*pinv(i,j)*86400
449 molod 1.1 enddo
450     enddo
451     endif
452    
453 molod 1.9 if( (bi.eq.1) .and. (bj.eq.1) ) then
454 molod 1.1 nvdtmoist = nvdtmoist + 1
455     nvdtturb = nvdtturb + 1
456     nvdtradlw = nvdtradlw + 1
457     nvdtradsw = nvdtradsw + 1
458 molod 1.9 endif
459 molod 1.1
460 molod 1.8 #endif
461 molod 1.1 return
462     end

  ViewVC Help
Powered by ViewVC 1.1.22