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

Last change on this file since 790 was 790, checked in by raasch, 10 years ago

Bugfix for output of mean particle radius + preliminary works for implementing the Wang collision kernel

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