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

Last change on this file since 1093 was 1054, checked in by hoffmann, 12 years ago

last commit documented

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