/[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.12 - (show annotations) (download)
Fri Oct 22 02:24:04 2004 UTC (19 years, 8 months ago) by jmc
Branch: MAIN
CVS Tags: checkpoint55j_post, checkpoint56b_post, checkpoint56c_post, checkpoint57a_post, checkpoint57a_pre, checkpoint55i_post, checkpoint57, checkpoint56, checkpoint56a_post
Changes since 1.11: +2 -2 lines
fix a syntax error (with ifc)

1 C $Header: /u/gcmpack/MITgcm/pkg/fizhi/fizhi_step_diag.F,v 1.11 2004/08/29 19:43:43 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
225 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 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
238 c Solar Radiative Heating (deg/day)
239 c ---------------------------------
240 if (iradsw.ne.0) then
241 do j=jm1,jm2
242 do i=im1,im2
243 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
250 c Clear Sky Solar Radiative Heating (deg/day)
251 c -------------------------------------------
252 if (iswclr.ne.0) then
253 do j=jm1,jm2
254 do i=im1,im2
255 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 . pk(i,j,l,bi,bj)*pinv(i,j)*86400
258 enddo
259 enddo
260 endif
261
262 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 qdiag(i,j,iuwnd+L-1,bi,bj) = qdiag(i,j,iuwnd+L-1,bi,bj) +
268 . uphy(i,j,L,bi,bj)
269 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 qdiag(i,j,ivwnd+L-1,bi,bj) = qdiag(i,j,ivwnd+L-1,bi,bj) +
279 . vphy(i,j,L,bi,bj)
280 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 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 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 qdiag(i,j,itke+L-1,bi,bj) = qdiag(i,j,itke+L-1,bi,bj) + qq(i,j,L)
301 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 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 enddo
313 enddo
314 endif
315
316 enddo
317
318 if( (bi.eq.1) .and. (bj.eq.1) ) then
319
320 ndiabu = ndiabu + 1
321 ndiabv = ndiabv + 1
322 ndiabt = ndiabt + 1
323 ndiabq = ndiabq + 1
324 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
334 endif
335
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 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 enddo
352 enddo
353 enddo
354 do j=jm1,jm2
355 do i=im1,im2
356 qdiag(i,j,ivdtmoist,bi,bj) = qdiag(i,j,ivdtmoist,bi,bj) +
357 . qbar(i,j)*pinv(i,j)*pinv(i,j)*86400
358 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 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 enddo
376 enddo
377 enddo
378 do j=jm1,jm2
379 do i=im1,im2
380 qdiag(i,j,ivdtturb,bi,bj) = qdiag(i,j,ivdtturb,bi,bj) +
381 . qbar(i,j)*pinv(i,j)*pinv(i,j)*86400
382 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 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 enddo
401 enddo
402 enddo
403 do j=jm1,jm2
404 do i=im1,im2
405 qdiag(i,j,ivdtradlw,bi,bj) = qdiag(i,j,ivdtradlw,bi,bj) +
406 . qbar(i,j)*pinv(i,j)*pinv(i,j)*86400
407 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 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 enddo
425 enddo
426 enddo
427 do j=jm1,jm2
428 do i=im1,im2
429 qdiag(i,j,ivdtradsw,bi,bj) = qdiag(i,j,ivdtradsw,bi,bj) +
430 . qbar(i,j)*radswt(i,j,bi,bj)*pinv(i,j)*pinv(i,j)*86400
431 enddo
432 enddo
433 endif
434
435 if( (bi.eq.1) .and. (bj.eq.1) ) then
436 nvdtmoist = nvdtmoist + 1
437 nvdtturb = nvdtturb + 1
438 nvdtradlw = nvdtradlw + 1
439 nvdtradsw = nvdtradsw + 1
440 endif
441
442 #endif
443 return
444 end

  ViewVC Help
Powered by ViewVC 1.1.22