/[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.2 - (show 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 C $Header: $
2 C $Name: $
3
4 #include "CPP_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,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 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
34 integer i,j,L
35 real pinv(im2,jm2), qbar(im2,jm2)
36
37 C **********************************************************************
38
39 do j=jm1,jm2
40 do i=im1,im2
41 pinv(i,j) = 1.0 / p(i,j,bi,bj)
42 enddo
43 enddo
44
45 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 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
56 c Net Solar Radiation at the Ground (W/m**2)
57 c ------------------------------------------
58 if (iradswg.ne.0) then
59 do j=jm1,jm2
60 do i=im1,im2
61 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
67 c Net Clear Sky Solar Radiation at the Ground (W/m**2)
68 c ----------------------------------------------------
69 if (iswgclr.ne.0) then
70 do j=jm1,jm2
71 do i=im1,im2
72 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
78 c Outgoing Solar Radiation at top (W/m**2)
79 c -----------------------------------------
80 if (iosr.ne.0) then
81 do j=jm1,jm2
82 do i=im1,im2
83 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
89 c Outgoing Clear Sky Solar Radiation at top (W/m**2)
90 c ---------------------------------------------------
91 if (iosrclr.ne.0) then
92 do j=jm1,jm2
93 do i=im1,im2
94 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
100 c Upward Longwave Flux at the Ground (W/m**2)
101 c -------------------------------------------
102 if (ilwgup.ne.0) then
103 do j=jm1,jm2
104 do i=im1,im2
105 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
111 c Net Longwave Flux at the Ground (W/m**2)
112 c ----------------------------------------
113 if (iradlwg.ne.0) then
114 do j=jm1,jm2
115 do i=im1,im2
116 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
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 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 enddo
132 enddo
133 endif
134
135 nradswt = nradswt + 1
136 nradswg = nradswg + 1
137 nswgclr = nswgclr + 1
138 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 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 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 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 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 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 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 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 enddo
193 enddo
194 endif
195
196 c Longwave Heating (deg/day)
197 c --------------------------
198 if (iradlw.ne.0) then
199 do j=jm1,jm2
200 do i=im1,im2
201 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 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 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
222 c Solar Radiative Heating (deg/day)
223 c ---------------------------------
224 if (iradsw.ne.0) then
225 do j=jm1,jm2
226 do i=im1,im2
227 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
234 c Clear Sky Solar Radiative Heating (deg/day)
235 c -------------------------------------------
236 if (iswclr.ne.0) then
237 do j=jm1,jm2
238 do i=im1,im2
239 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 enddo
243 enddo
244 endif
245
246 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 qdiag(i,j,iuwnd+L-1,bi,bj) = qdiag(i,j,iuwnd+L-1,bi,bj) +
252 . uphy(i,j,L,bi,bj)
253 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 qdiag(i,j,ivwnd+L-1,bi,bj) = qdiag(i,j,ivwnd+L-1,bi,bj) +
263 . vphy(i,j,L,bi,bj)
264 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 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 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 qdiag(i,j,itke+L-1,bi,bj) = qdiag(i,j,itke+L-1,bi,bj) + qq(i,j,L)
285 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 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 enddo
297 enddo
298 endif
299
300 enddo
301
302 ndiabu = ndiabu + 1
303 ndiabv = ndiabv + 1
304 ndiabt = ndiabt + 1
305 ndiabq = ndiabq + 1
306 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 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 enddo
332 enddo
333 enddo
334 do j=jm1,jm2
335 do i=im1,im2
336 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 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 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 enddo
356 enddo
357 enddo
358 do j=jm1,jm2
359 do i=im1,im2
360 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 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 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 enddo
381 enddo
382 enddo
383 do j=jm1,jm2
384 do i=im1,im2
385 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 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 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 enddo
405 enddo
406 enddo
407 do j=jm1,jm2
408 do i=im1,im2
409 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 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