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

Last change on this file since 1357 was 1354, checked in by heinze, 11 years ago

last commit documented

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