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

Last change on this file since 1036 was 1036, checked in by raasch, 9 years ago

code has been put under the GNU General Public License (v3)

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