source: palm/trunk/SOURCE/sum_up_3d_data.f90 @ 402

Last change on this file since 402 was 402, checked in by maronga, 12 years ago

bugfix: calculation of time-averaged surface-heatfluxes

  • Property svn:keywords set to Id
File size: 14.1 KB
Line 
1 SUBROUTINE sum_up_3d_data
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6! Bugfix in calculation of shf*_av, qsws*_av
7!
8! Former revisions:
9! -----------------
10! $Id: sum_up_3d_data.f90 402 2009-10-21 11:59:41Z maronga $
11!
12! 2009-08-25 08:35:52Z maronga
13! +shf*, qsws*
14!
15! 96 2007-06-04 08:07:41Z raasch
16! +sum-up of density and salinity
17!
18! 72 2007-03-19 08:20:46Z raasch
19! +sum-up of precipitation rate and roughness length (prr*, z0*)
20!
21! RCS Log replace by Id keyword, revision history cleaned up
22!
23! Revision 1.1  2006/02/23 12:55:23  raasch
24! Initial revision
25!
26!
27! Description:
28! ------------
29! Sum-up the values of 3d-arrays. The real averaging is later done in routine
30! average_3d_data.
31!------------------------------------------------------------------------------!
32
33    USE arrays_3d
34    USE averaging
35    USE cloud_parameters
36    USE control_parameters
37    USE cpulog
38    USE indices
39    USE interfaces
40    USE particle_attributes
41
42    IMPLICIT NONE
43
44    INTEGER ::  i, ii, j, k, n, psi
45
46    REAL    ::  mean_r, s_r3, s_r4
47
48
49    CALL cpu_log (log_point(34),'sum_up_3d_data','start')
50
51!
52!-- Allocate and initialize the summation arrays if called for the very first
53!-- time or the first time after average_3d_data has been called
54!-- (some or all of the arrays may have been already allocated
55!-- in read_3d_binary)
56    IF ( average_count_3d == 0 )  THEN
57
58       DO  ii = 1, doav_n
59
60          SELECT CASE ( TRIM( doav(ii) ) )
61
62             CASE ( 'e' )
63                IF ( .NOT. ALLOCATED( e_av ) )  THEN
64                   ALLOCATE( e_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
65                ENDIF
66                e_av = 0.0
67
68             CASE ( 'lwp*' )
69                IF ( .NOT. ALLOCATED( lwp_av ) )  THEN
70                   ALLOCATE( lwp_av(nys-1:nyn+1,nxl-1:nxr+1) )
71                ENDIF
72                lwp_av = 0.0
73
74             CASE ( 'p' )
75                IF ( .NOT. ALLOCATED( p_av ) )  THEN
76                   ALLOCATE( p_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
77                ENDIF
78                p_av = 0.0
79
80             CASE ( 'pc' )
81                IF ( .NOT. ALLOCATED( pc_av ) )  THEN
82                   ALLOCATE( pc_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
83                ENDIF
84                pc_av = 0.0
85
86             CASE ( 'pr' )
87                IF ( .NOT. ALLOCATED( pr_av ) )  THEN
88                   ALLOCATE( pr_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
89                ENDIF
90                pr_av = 0.0
91
92             CASE ( 'prr*' )
93                IF ( .NOT. ALLOCATED( precipitation_rate_av ) )  THEN
94                   ALLOCATE( precipitation_rate_av(nys-1:nyn+1,nxl-1:nxr+1) )
95                ENDIF
96                precipitation_rate_av = 0.0
97
98             CASE ( 'pt' )
99                IF ( .NOT. ALLOCATED( pt_av ) )  THEN
100                   ALLOCATE( pt_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
101                ENDIF
102                pt_av = 0.0
103
104             CASE ( 'q' )
105                IF ( .NOT. ALLOCATED( q_av ) )  THEN
106                   ALLOCATE( q_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
107                ENDIF
108                q_av = 0.0
109
110             CASE ( 'ql' )
111                IF ( .NOT. ALLOCATED( ql_av ) )  THEN
112                   ALLOCATE( ql_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
113                ENDIF
114                ql_av = 0.0
115
116             CASE ( 'ql_c' )
117                IF ( .NOT. ALLOCATED( ql_c_av ) )  THEN
118                   ALLOCATE( ql_c_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
119                ENDIF
120                ql_c_av = 0.0
121
122             CASE ( 'ql_v' )
123                IF ( .NOT. ALLOCATED( ql_v_av ) )  THEN
124                   ALLOCATE( ql_v_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
125                ENDIF
126                ql_v_av = 0.0
127
128             CASE ( 'ql_vp' )
129                IF ( .NOT. ALLOCATED( ql_vp_av ) )  THEN
130                   ALLOCATE( ql_vp_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
131                ENDIF
132                ql_vp_av = 0.0
133
134             CASE ( 'qsws*' )
135                IF ( .NOT. ALLOCATED( qsws_av ) )  THEN
136                   ALLOCATE( qsws_av(nys-1:nyn+1,nxl-1:nxr+1) )
137                ENDIF
138                qsws_av = 0.0
139
140             CASE ( 'qv' )
141                IF ( .NOT. ALLOCATED( qv_av ) )  THEN
142                   ALLOCATE( qv_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
143                ENDIF
144                qv_av = 0.0
145
146             CASE ( 'rho' )
147                IF ( .NOT. ALLOCATED( rho_av ) )  THEN
148                   ALLOCATE( rho_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
149                ENDIF
150                rho_av = 0.0
151
152             CASE ( 's' )
153                IF ( .NOT. ALLOCATED( s_av ) )  THEN
154                   ALLOCATE( s_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
155                ENDIF
156                s_av = 0.0
157
158             CASE ( 'sa' )
159                IF ( .NOT. ALLOCATED( sa_av ) )  THEN
160                   ALLOCATE( sa_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
161                ENDIF
162                sa_av = 0.0
163
164             CASE ( 'shf*' )
165                IF ( .NOT. ALLOCATED( shf_av ) )  THEN
166                   ALLOCATE( shf_av(nys-1:nyn+1,nxl-1:nxr+1) )
167                ENDIF
168                shf_av = 0.0
169
170             CASE ( 't*' )
171                IF ( .NOT. ALLOCATED( ts_av ) )  THEN
172                   ALLOCATE( ts_av(nys-1:nyn+1,nxl-1:nxr+1) )
173                ENDIF
174                ts_av = 0.0
175
176             CASE ( 'u' )
177                IF ( .NOT. ALLOCATED( u_av ) )  THEN
178                   ALLOCATE( u_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
179                ENDIF
180                u_av = 0.0
181
182             CASE ( 'u*' )
183                IF ( .NOT. ALLOCATED( us_av ) )  THEN
184                   ALLOCATE( us_av(nys-1:nyn+1,nxl-1:nxr+1) )
185                ENDIF
186                us_av = 0.0
187
188             CASE ( 'v' )
189                IF ( .NOT. ALLOCATED( v_av ) )  THEN
190                   ALLOCATE( v_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
191                ENDIF
192                v_av = 0.0
193
194             CASE ( 'vpt' )
195                IF ( .NOT. ALLOCATED( vpt_av ) )  THEN
196                   ALLOCATE( vpt_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
197                ENDIF
198                vpt_av = 0.0
199
200             CASE ( 'w' )
201                IF ( .NOT. ALLOCATED( w_av ) )  THEN
202                   ALLOCATE( w_av(nzb:nzt+1,nys-1:nyn+1,nxl-1:nxr+1) )
203                ENDIF
204                w_av = 0.0
205
206             CASE ( 'z0*' )
207                IF ( .NOT. ALLOCATED( z0_av ) )  THEN
208                   ALLOCATE( z0_av(nys-1:nyn+1,nxl-1:nxr+1) )
209                ENDIF
210                z0_av = 0.0
211
212             CASE DEFAULT
213!
214!--             User-defined quantity
215                CALL user_3d_data_averaging( 'allocate', doav(ii) )
216
217          END SELECT
218
219       ENDDO
220
221    ENDIF
222
223!
224!-- Loop of all variables to be averaged.
225    DO  ii = 1, doav_n
226
227!
228!--    Store the array chosen on the temporary array.
229       SELECT CASE ( TRIM( doav(ii) ) )
230
231          CASE ( 'e' )
232             DO  i = nxl-1, nxr+1
233                DO  j = nys-1, nyn+1
234                   DO  k = nzb, nzt+1
235                      e_av(k,j,i) = e_av(k,j,i) + e(k,j,i)
236                   ENDDO
237                ENDDO
238             ENDDO
239
240          CASE ( 'lwp*' )
241             DO  i = nxl-1, nxr+1
242                DO  j = nys-1, nyn+1
243                   lwp_av(j,i) = lwp_av(j,i) + SUM( ql(nzb:nzt,j,i) * &
244                                                    dzw(1:nzt+1) )
245                ENDDO
246             ENDDO
247
248          CASE ( 'p' )
249             DO  i = nxl-1, nxr+1
250                DO  j = nys-1, nyn+1
251                   DO  k = nzb, nzt+1
252                      p_av(k,j,i) = p_av(k,j,i) + p(k,j,i)
253                   ENDDO
254                ENDDO
255             ENDDO
256
257          CASE ( 'pc' )
258             DO  i = nxl, nxr
259                DO  j = nys, nyn
260                   DO  k = nzb, nzt+1
261                      pc_av(k,j,i) = pc_av(k,j,i) + prt_count(k,j,i)
262                   ENDDO
263                ENDDO
264             ENDDO
265
266          CASE ( 'pr' )
267             DO  i = nxl, nxr
268                DO  j = nys, nyn
269                   DO  k = nzb, nzt+1
270                      psi = prt_start_index(k,j,i)
271                      s_r3 = 0.0
272                      s_r4 = 0.0
273                      DO  n = psi, psi+prt_count(k,j,i)-1
274                         s_r3 = s_r3 + particles(n)%radius**3
275                         s_r4 = s_r4 + particles(n)%radius**4
276                      ENDDO
277                      IF ( s_r3 /= 0.0 )  THEN
278                         mean_r = s_r4 / s_r3
279                      ELSE
280                         mean_r = 0.0
281                      ENDIF
282                      pr_av(k,j,i) = pr_av(k,j,i) + mean_r
283                   ENDDO
284                ENDDO
285             ENDDO
286
287          CASE ( 'pr*' )
288             DO  i = nxl-1, nxr+1
289                DO  j = nys-1, nyn+1
290                   precipitation_rate_av(j,i) = precipitation_rate_av(j,i) + &
291                                                precipitation_rate(j,i)
292                ENDDO
293             ENDDO
294
295          CASE ( 'pt' )
296             IF ( .NOT. cloud_physics ) THEN
297                DO  i = nxl-1, nxr+1
298                   DO  j = nys-1, nyn+1
299                      DO  k = nzb, nzt+1
300                         pt_av(k,j,i) = pt_av(k,j,i) + pt(k,j,i)
301                      ENDDO
302                   ENDDO
303                ENDDO
304             ELSE
305                DO  i = nxl-1, nxr+1
306                   DO  j = nys-1, nyn+1
307                      DO  k = nzb, nzt+1
308                         pt_av(k,j,i) = pt_av(k,j,i) + pt(k,j,i) + l_d_cp * &
309                                                       pt_d_t(k) * ql(k,j,i)
310                      ENDDO
311                   ENDDO
312                ENDDO
313             ENDIF
314
315          CASE ( 'q' )
316             DO  i = nxl-1, nxr+1
317                DO  j = nys-1, nyn+1
318                   DO  k = nzb, nzt+1
319                      q_av(k,j,i) = q_av(k,j,i) + q(k,j,i)
320                   ENDDO
321                ENDDO
322             ENDDO
323
324          CASE ( 'ql' )
325             DO  i = nxl-1, nxr+1
326                DO  j = nys-1, nyn+1
327                   DO  k = nzb, nzt+1
328                      ql_av(k,j,i) = ql_av(k,j,i) + ql(k,j,i)
329                   ENDDO
330                ENDDO
331             ENDDO
332
333          CASE ( 'ql_c' )
334             DO  i = nxl-1, nxr+1
335                DO  j = nys-1, nyn+1
336                   DO  k = nzb, nzt+1
337                      ql_c_av(k,j,i) = ql_c_av(k,j,i) + ql_c(k,j,i)
338                   ENDDO
339                ENDDO
340             ENDDO
341
342          CASE ( 'ql_v' )
343             DO  i = nxl-1, nxr+1
344                DO  j = nys-1, nyn+1
345                   DO  k = nzb, nzt+1
346                      ql_v_av(k,j,i) = ql_v_av(k,j,i) + ql_v(k,j,i)
347                   ENDDO
348                ENDDO
349             ENDDO
350
351          CASE ( 'ql_vp' )
352             DO  i = nxl-1, nxr+1
353                DO  j = nys-1, nyn+1
354                   DO  k = nzb, nzt+1
355                      ql_vp_av(k,j,i) = ql_vp_av(k,j,i) + ql_vp(k,j,i)
356                   ENDDO
357                ENDDO
358             ENDDO
359
360          CASE ( 'qsws*' )
361             DO  i = nxl-1, nxr+1
362                DO  j = nys-1, nyn+1
363                   qsws_av(j,i) = qsws_av(j,i) + qsws(j,i)
364                ENDDO
365             ENDDO
366
367          CASE ( 'qv' )
368             DO  i = nxl-1, nxr+1
369                DO  j = nys-1, nyn+1
370                   DO  k = nzb, nzt+1
371                      qv_av(k,j,i) = qv_av(k,j,i) + q(k,j,i) - ql(k,j,i)
372                   ENDDO
373                ENDDO
374             ENDDO
375
376          CASE ( 'rho' )
377             DO  i = nxl-1, nxr+1
378                DO  j = nys-1, nyn+1
379                   DO  k = nzb, nzt+1
380                      rho_av(k,j,i) = rho_av(k,j,i) + rho(k,j,i)
381                   ENDDO
382                ENDDO
383             ENDDO
384
385          CASE ( 's' )
386             DO  i = nxl-1, nxr+1
387                DO  j = nys-1, nyn+1
388                   DO  k = nzb, nzt+1
389                      s_av(k,j,i) = s_av(k,j,i) + q(k,j,i)
390                   ENDDO
391                ENDDO
392             ENDDO
393
394          CASE ( 'sa' )
395             DO  i = nxl-1, nxr+1
396                DO  j = nys-1, nyn+1
397                   DO  k = nzb, nzt+1
398                      sa_av(k,j,i) = sa_av(k,j,i) + sa(k,j,i)
399                   ENDDO
400                ENDDO
401             ENDDO
402
403          CASE ( 'shf*' )
404             DO  i = nxl-1, nxr+1
405                DO  j = nys-1, nyn+1
406                   shf_av(j,i) = shf_av(j,i) + shf(j,i)
407                ENDDO
408             ENDDO
409
410          CASE ( 't*' )
411             DO  i = nxl-1, nxr+1
412                DO  j = nys-1, nyn+1
413                   ts_av(j,i) = ts_av(j,i) + ts(j,i)
414                ENDDO
415             ENDDO
416
417          CASE ( 'u' )
418             DO  i = nxl-1, nxr+1
419                DO  j = nys-1, nyn+1
420                   DO  k = nzb, nzt+1
421                      u_av(k,j,i) = u_av(k,j,i) + u(k,j,i)
422                   ENDDO
423                ENDDO
424             ENDDO
425
426          CASE ( 'u*' )
427             DO  i = nxl-1, nxr+1
428                DO  j = nys-1, nyn+1
429                   us_av(j,i) = us_av(j,i) + us(j,i)
430                ENDDO
431             ENDDO
432
433          CASE ( 'v' )
434             DO  i = nxl-1, nxr+1
435                DO  j = nys-1, nyn+1
436                   DO  k = nzb, nzt+1
437                      v_av(k,j,i) = v_av(k,j,i) + v(k,j,i)
438                   ENDDO
439                ENDDO
440             ENDDO
441
442          CASE ( 'vpt' )
443             DO  i = nxl-1, nxr+1
444                DO  j = nys-1, nyn+1
445                   DO  k = nzb, nzt+1
446                      vpt_av(k,j,i) = vpt_av(k,j,i) + vpt(k,j,i)
447                   ENDDO
448                ENDDO
449             ENDDO
450
451          CASE ( 'w' )
452             DO  i = nxl-1, nxr+1
453                DO  j = nys-1, nyn+1
454                   DO  k = nzb, nzt+1
455                      w_av(k,j,i) = w_av(k,j,i) + w(k,j,i)
456                   ENDDO
457                ENDDO
458             ENDDO
459
460          CASE ( 'z0*' )
461             DO  i = nxl-1, nxr+1
462                DO  j = nys-1, nyn+1
463                   z0_av(j,i) = z0_av(j,i) + z0(j,i)
464                ENDDO
465             ENDDO
466
467          CASE DEFAULT
468!
469!--          User-defined quantity
470             CALL user_3d_data_averaging( 'sum', doav(ii) )
471
472       END SELECT
473
474    ENDDO
475
476    CALL cpu_log (log_point(34),'sum_up_3d_data','stop','nobarrier')
477
478
479 END SUBROUTINE sum_up_3d_data
Note: See TracBrowser for help on using the repository browser.