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

Last change on this file since 790 was 790, checked in by raasch, 12 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
RevLine 
[1]1 SUBROUTINE sum_up_3d_data
2
3!------------------------------------------------------------------------------!
[484]4! Current revisions:
[1]5! -----------------
[790]6! bugfix: calculation of 'pr' must depend on the particle weighting factor
[1]7!
8! Former revisions:
9! -----------------
[3]10! $Id: sum_up_3d_data.f90 790 2011-11-29 03:11:20Z raasch $
[77]11!
[772]12! 771 2011-10-27 10:56:21Z heinze
13! +lpt_av
14!
[668]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!
[482]18! 402 2009-10-21 11:59:41Z maronga
19! Bugfix in calculation of shf*_av, qsws*_av
20!
[392]21! 2009-08-25 08:35:52Z maronga
22! +shf*, qsws*
23!
[98]24! 96 2007-06-04 08:07:41Z raasch
25! +sum-up of density and salinity
26!
[77]27! 72 2007-03-19 08:20:46Z raasch
28! +sum-up of precipitation rate and roughness length (prr*, z0*)
29!
[3]30! RCS Log replace by Id keyword, revision history cleaned up
31!
[1]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
[667]73                   ALLOCATE( e_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1]74                ENDIF
75                e_av = 0.0
76
[771]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
[1]83             CASE ( 'lwp*' )
84                IF ( .NOT. ALLOCATED( lwp_av ) )  THEN
[667]85                   ALLOCATE( lwp_av(nysg:nyng,nxlg:nxrg) )
[1]86                ENDIF
87                lwp_av = 0.0
88
89             CASE ( 'p' )
90                IF ( .NOT. ALLOCATED( p_av ) )  THEN
[667]91                   ALLOCATE( p_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1]92                ENDIF
93                p_av = 0.0
94
95             CASE ( 'pc' )
96                IF ( .NOT. ALLOCATED( pc_av ) )  THEN
[667]97                   ALLOCATE( pc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1]98                ENDIF
99                pc_av = 0.0
100
101             CASE ( 'pr' )
102                IF ( .NOT. ALLOCATED( pr_av ) )  THEN
[667]103                   ALLOCATE( pr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1]104                ENDIF
105                pr_av = 0.0
106
[72]107             CASE ( 'prr*' )
108                IF ( .NOT. ALLOCATED( precipitation_rate_av ) )  THEN
[667]109                   ALLOCATE( precipitation_rate_av(nysg:nyng,nxlg:nxrg) )
[72]110                ENDIF
111                precipitation_rate_av = 0.0
112
[1]113             CASE ( 'pt' )
114                IF ( .NOT. ALLOCATED( pt_av ) )  THEN
[667]115                   ALLOCATE( pt_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1]116                ENDIF
117                pt_av = 0.0
118
119             CASE ( 'q' )
120                IF ( .NOT. ALLOCATED( q_av ) )  THEN
[667]121                   ALLOCATE( q_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1]122                ENDIF
123                q_av = 0.0
124
125             CASE ( 'ql' )
126                IF ( .NOT. ALLOCATED( ql_av ) )  THEN
[667]127                   ALLOCATE( ql_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1]128                ENDIF
129                ql_av = 0.0
130
131             CASE ( 'ql_c' )
132                IF ( .NOT. ALLOCATED( ql_c_av ) )  THEN
[667]133                   ALLOCATE( ql_c_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1]134                ENDIF
135                ql_c_av = 0.0
136
137             CASE ( 'ql_v' )
138                IF ( .NOT. ALLOCATED( ql_v_av ) )  THEN
[667]139                   ALLOCATE( ql_v_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1]140                ENDIF
141                ql_v_av = 0.0
142
143             CASE ( 'ql_vp' )
144                IF ( .NOT. ALLOCATED( ql_vp_av ) )  THEN
[667]145                   ALLOCATE( ql_vp_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1]146                ENDIF
147                ql_vp_av = 0.0
148
[354]149             CASE ( 'qsws*' )
150                IF ( .NOT. ALLOCATED( qsws_av ) )  THEN
[667]151                   ALLOCATE( qsws_av(nysg:nyng,nxlg:nxrg) )
[354]152                ENDIF
153                qsws_av = 0.0
154
[1]155             CASE ( 'qv' )
156                IF ( .NOT. ALLOCATED( qv_av ) )  THEN
[667]157                   ALLOCATE( qv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1]158                ENDIF
159                qv_av = 0.0
160
[96]161             CASE ( 'rho' )
162                IF ( .NOT. ALLOCATED( rho_av ) )  THEN
[667]163                   ALLOCATE( rho_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[96]164                ENDIF
165                rho_av = 0.0
166
[1]167             CASE ( 's' )
168                IF ( .NOT. ALLOCATED( s_av ) )  THEN
[667]169                   ALLOCATE( s_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1]170                ENDIF
171                s_av = 0.0
172
[96]173             CASE ( 'sa' )
174                IF ( .NOT. ALLOCATED( sa_av ) )  THEN
[667]175                   ALLOCATE( sa_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[96]176                ENDIF
177                sa_av = 0.0
178
[354]179             CASE ( 'shf*' )
180                IF ( .NOT. ALLOCATED( shf_av ) )  THEN
[667]181                   ALLOCATE( shf_av(nysg:nyng,nxlg:nxrg) )
[354]182                ENDIF
183                shf_av = 0.0
184
[1]185             CASE ( 't*' )
186                IF ( .NOT. ALLOCATED( ts_av ) )  THEN
[667]187                   ALLOCATE( ts_av(nysg:nyng,nxlg:nxrg) )
[1]188                ENDIF
189                ts_av = 0.0
190
191             CASE ( 'u' )
192                IF ( .NOT. ALLOCATED( u_av ) )  THEN
[667]193                   ALLOCATE( u_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1]194                ENDIF
195                u_av = 0.0
196
197             CASE ( 'u*' )
198                IF ( .NOT. ALLOCATED( us_av ) )  THEN
[667]199                   ALLOCATE( us_av(nysg:nyng,nxlg:nxrg) )
[1]200                ENDIF
201                us_av = 0.0
202
203             CASE ( 'v' )
204                IF ( .NOT. ALLOCATED( v_av ) )  THEN
[667]205                   ALLOCATE( v_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1]206                ENDIF
207                v_av = 0.0
208
209             CASE ( 'vpt' )
210                IF ( .NOT. ALLOCATED( vpt_av ) )  THEN
[667]211                   ALLOCATE( vpt_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1]212                ENDIF
213                vpt_av = 0.0
214
215             CASE ( 'w' )
216                IF ( .NOT. ALLOCATED( w_av ) )  THEN
[667]217                   ALLOCATE( w_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1]218                ENDIF
219                w_av = 0.0
220
[72]221             CASE ( 'z0*' )
222                IF ( .NOT. ALLOCATED( z0_av ) )  THEN
[667]223                   ALLOCATE( z0_av(nysg:nyng,nxlg:nxrg) )
[72]224                ENDIF
225                z0_av = 0.0
226
[1]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' )
[667]247             DO  i = nxlg, nxrg
248                DO  j = nysg, nyng
[1]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
[771]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
[1]264          CASE ( 'lwp*' )
[667]265             DO  i = nxlg, nxrg
266                DO  j = nysg, nyng
[1]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' )
[667]273             DO  i = nxlg, nxrg
274                DO  j = nysg, nyng
[1]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
[790]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
[1]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
[72]313          CASE ( 'pr*' )
[667]314             DO  i = nxlg, nxrg
315                DO  j = nysg, nyng
[72]316                   precipitation_rate_av(j,i) = precipitation_rate_av(j,i) + &
317                                                precipitation_rate(j,i)
318                ENDDO
319             ENDDO
320
[1]321          CASE ( 'pt' )
322             IF ( .NOT. cloud_physics ) THEN
[667]323             DO  i = nxlg, nxrg
324                DO  j = nysg, nyng
325                   DO  k = nzb, nzt+1
[1]326                         pt_av(k,j,i) = pt_av(k,j,i) + pt(k,j,i)
327                      ENDDO
328                   ENDDO
329                ENDDO
330             ELSE
[667]331             DO  i = nxlg, nxrg
332                DO  j = nysg, nyng
333                   DO  k = nzb, nzt+1
[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' )
[667]342             DO  i = nxlg, nxrg
343                DO  j = nysg, nyng
[1]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
[402]349
[1]350          CASE ( 'ql' )
[667]351             DO  i = nxlg, nxrg
352                DO  j = nysg, nyng
[1]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' )
[667]360             DO  i = nxlg, nxrg
361                DO  j = nysg, nyng
[1]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' )
[667]369             DO  i = nxlg, nxrg
370                DO  j = nysg, nyng
[1]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' )
[667]378             DO  i = nxlg, nxrg
379                DO  j = nysg, nyng
[1]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
[402]386          CASE ( 'qsws*' )
[667]387             DO  i = nxlg, nxrg
388                DO  j = nysg, nyng
[402]389                   qsws_av(j,i) = qsws_av(j,i) + qsws(j,i)
390                ENDDO
391             ENDDO
392
[1]393          CASE ( 'qv' )
[667]394             DO  i = nxlg, nxrg
395                DO  j = nysg, nyng
[1]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
[96]402          CASE ( 'rho' )
[667]403             DO  i = nxlg, nxrg
404                DO  j = nysg, nyng
[96]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
[402]410
[1]411          CASE ( 's' )
[667]412             DO  i = nxlg, nxrg
413                DO  j = nysg, nyng
[1]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
[402]419
[96]420          CASE ( 'sa' )
[667]421             DO  i = nxlg, nxrg
422                DO  j = nysg, nyng
[96]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
[402]428
429          CASE ( 'shf*' )
[667]430             DO  i = nxlg, nxrg
431                DO  j = nysg, nyng
[402]432                   shf_av(j,i) = shf_av(j,i) + shf(j,i)
433                ENDDO
434             ENDDO
435
[1]436          CASE ( 't*' )
[667]437             DO  i = nxlg, nxrg
438                DO  j = nysg, nyng
[1]439                   ts_av(j,i) = ts_av(j,i) + ts(j,i)
440                ENDDO
441             ENDDO
442
443          CASE ( 'u' )
[667]444             DO  i = nxlg, nxrg
445                DO  j = nysg, nyng
[1]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*' )
[667]453             DO  i = nxlg, nxrg
454                DO  j = nysg, nyng
[1]455                   us_av(j,i) = us_av(j,i) + us(j,i)
456                ENDDO
457             ENDDO
458
459          CASE ( 'v' )
[667]460             DO  i = nxlg, nxrg
461                DO  j = nysg, nyng
[1]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' )
[667]469             DO  i = nxlg, nxrg
470                DO  j = nysg, nyng
[1]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' )
[667]478             DO  i = nxlg, nxrg
479                DO  j = nysg, nyng
[1]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
[72]486          CASE ( 'z0*' )
[667]487             DO  i = nxlg, nxrg
488                DO  j = nysg, nyng
[72]489                   z0_av(j,i) = z0_av(j,i) + z0(j,i)
490                ENDDO
491             ENDDO
492
[1]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.