/[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.4 - (hide annotations) (download)
Tue Jul 13 23:44:43 2004 UTC (19 years, 11 months ago) by molod
Branch: MAIN
Changes since 1.3: +5 -2 lines
Debugging

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

  ViewVC Help
Powered by ViewVC 1.1.22