/[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.2 - (hide annotations) (download)
Thu Jun 24 18:56:57 2004 UTC (20 years ago) by molod
Branch: MAIN
CVS Tags: checkpoint54a_pre, checkpoint54a_post, checkpoint54b_post, checkpoint54, checkpoint53g_post, checkpoint53f_post
Changes since 1.1: +141 -441 lines
Fizhi diagnostics that are computed every dynamics time step

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

  ViewVC Help
Powered by ViewVC 1.1.22