/[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.11 - (hide annotations) (download)
Sun Aug 29 19:43:43 2004 UTC (19 years, 10 months ago) by molod
Branch: MAIN
CVS Tags: checkpoint55c_post, checkpoint54e_post, checkpoint55d_pre, checkpoint55h_post, checkpoint55b_post, checkpoint55, checkpoint54f_post, checkpoint55g_post, checkpoint55f_post, checkpoint55e_post, checkpoint55a_post, checkpoint55d_post
Changes since 1.10: +3 -21 lines
Remove user diagnostics calls from checked in version

1 molod 1.11 C $Header: /u/gcmpack/MITgcm/pkg/fizhi/fizhi_step_diag.F,v 1.10 2004/08/12 15:21:22 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.11
225 molod 1.1 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 molod 1.2 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 molod 1.1
238     c Solar Radiative Heating (deg/day)
239     c ---------------------------------
240 molod 1.2 if (iradsw.ne.0) then
241 molod 1.1 do j=jm1,jm2
242     do i=im1,im2
243 molod 1.2 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 molod 1.1
250     c Clear Sky Solar Radiative Heating (deg/day)
251     c -------------------------------------------
252 molod 1.2 if (iswclr.ne.0) then
253 molod 1.1 do j=jm1,jm2
254     do i=im1,im2
255 molod 1.2 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 molod 1.5 . pk(i,j,l,bi,bj)*pinv(i,j)*86400
258 molod 1.1 enddo
259     enddo
260     endif
261 molod 1.2
262 molod 1.1 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 molod 1.2 qdiag(i,j,iuwnd+L-1,bi,bj) = qdiag(i,j,iuwnd+L-1,bi,bj) +
268     . uphy(i,j,L,bi,bj)
269 molod 1.1 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 molod 1.2 qdiag(i,j,ivwnd+L-1,bi,bj) = qdiag(i,j,ivwnd+L-1,bi,bj) +
279     . vphy(i,j,L,bi,bj)
280 molod 1.1 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 molod 1.2 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 molod 1.1 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 molod 1.2 qdiag(i,j,itke+L-1,bi,bj) = qdiag(i,j,itke+L-1,bi,bj) + qq(i,j,L)
301 molod 1.1 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 molod 1.2 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 molod 1.1 enddo
313     enddo
314     endif
315    
316 molod 1.2 enddo
317 molod 1.1
318 molod 1.9 if( (bi.eq.1) .and. (bj.eq.1) ) then
319 molod 1.11
320 molod 1.1 ndiabu = ndiabu + 1
321     ndiabv = ndiabv + 1
322     ndiabt = ndiabt + 1
323 molod 1.2 ndiabq = ndiabq + 1
324 molod 1.1 nradlw = nradlw + 1
325     nlwclr = nlwclr + 1
326     nradsw = nradsw + 1
327     nswclr = nswclr + 1
328     nuwnd = nuwnd + 1
329     nvwnd = nvwnd + 1
330     ntmpu = ntmpu + 1
331     ntke = ntke + 1
332     nsphu = nsphu + 1
333 molod 1.10
334 molod 1.9 endif
335 molod 1.1
336     C **********************************************************************
337    
338     c Vertically Averaged Moist-T Increment (K/day)
339     c ---------------------------------------------
340     if( ivdtmoist.ne.0 ) then
341     do j=jm1,jm2
342     do i=im1,im2
343     qbar(i,j) = 0.0
344     enddo
345     enddo
346     do L=1,Nrphys
347     do j=jm1,jm2
348     do i=im1,im2
349 molod 1.2 qbar(i,j) = qbar(i,j) +
350     . moistt(i,j,L,bi,bj)*pk(i,j,l,bi,bj)*dp(i,j,L,bi,bj)
351 molod 1.1 enddo
352     enddo
353     enddo
354     do j=jm1,jm2
355     do i=im1,im2
356 molod 1.2 qdiag(i,j,ivdtmoist,bi,bj) = qdiag(i,j,ivdtmoist,bi,bj) +
357 molod 1.5 . qbar(i,j)*pinv(i,j)*pinv(i,j)*86400
358 molod 1.1 enddo
359     enddo
360     endif
361    
362     c Vertically Averaged Turb-T Increment (K/day)
363     c --------------------------------------------
364     if( ivdtturb.ne.0 ) then
365     do j=jm1,jm2
366     do i=im1,im2
367     qbar(i,j) = 0.0
368     enddo
369     enddo
370     do L=1,Nrphys
371     do j=jm1,jm2
372     do i=im1,im2
373 molod 1.2 qbar(i,j) = qbar(i,j) +
374     . turbt(i,j,L,bi,bj)*pk(i,j,l,bi,bj)*dp(i,j,L,bi,bj)
375 molod 1.1 enddo
376     enddo
377     enddo
378     do j=jm1,jm2
379     do i=im1,im2
380 molod 1.2 qdiag(i,j,ivdtturb,bi,bj) = qdiag(i,j,ivdtturb,bi,bj) +
381 molod 1.5 . qbar(i,j)*pinv(i,j)*pinv(i,j)*86400
382 molod 1.1 enddo
383     enddo
384     endif
385    
386     c Vertically Averaged RADLW Temperature Increment (K/day)
387     c -------------------------------------------------------
388     if( ivdtradlw.ne.0 ) then
389     do j=jm1,jm2
390     do i=im1,im2
391     qbar(i,j) = 0.0
392     enddo
393     enddo
394     do L=1,Nrphys
395     do j=jm1,jm2
396     do i=im1,im2
397 molod 1.2 qbar(i,j) = qbar(i,j) + ( lwdt(i,j,L,bi,bj) +
398     . dlwdtg(i,j,L,bi,bj)*(tgz(i,j,bi,bj)-tg0(i,j,bi,bj)) )
399     . *pk(i,j,l,bi,bj)*dp(i,j,L,bi,bj)
400 molod 1.1 enddo
401     enddo
402     enddo
403     do j=jm1,jm2
404     do i=im1,im2
405 molod 1.2 qdiag(i,j,ivdtradlw,bi,bj) = qdiag(i,j,ivdtradlw,bi,bj) +
406 molod 1.5 . qbar(i,j)*pinv(i,j)*pinv(i,j)*86400
407 molod 1.1 enddo
408     enddo
409     endif
410    
411     c Vertically Averaged RADSW Temperature Increment (K/day)
412     c -------------------------------------------------------
413     if( ivdtradsw.ne.0 ) then
414     do j=jm1,jm2
415     do i=im1,im2
416     qbar(i,j) = 0.0
417     enddo
418     enddo
419     do L=1,Nrphys
420     do j=jm1,jm2
421     do i=im1,im2
422 molod 1.2 qbar(i,j) = qbar(i,j) +
423     . swdt(i,j,L,bi,bj)*pk(i,j,l,bi,bj)*dp(i,j,L,bi,bj)
424 molod 1.1 enddo
425     enddo
426     enddo
427     do j=jm1,jm2
428     do i=im1,im2
429 molod 1.2 qdiag(i,j,ivdtradsw,bi,bj) = qdiag(i,j,ivdtradsw,bi,bj) +
430 molod 1.5 . qbar(i,j)*radswt(i,j,bi,bj)*pinv(i,j)*pinv(i,j)*86400
431 molod 1.1 enddo
432     enddo
433     endif
434    
435 molod 1.9 if( (bi.eq.1) .and. (bj.eq.1) ) then
436 molod 1.1 nvdtmoist = nvdtmoist + 1
437     nvdtturb = nvdtturb + 1
438     nvdtradlw = nvdtradlw + 1
439     nvdtradsw = nvdtradsw + 1
440 molod 1.9 endif
441 molod 1.1
442 molod 1.8 #endif
443 molod 1.1 return
444     end

  ViewVC Help
Powered by ViewVC 1.1.22