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

Last change on this file since 1320 was 1320, checked in by raasch, 11 years ago

ONLY-attribute added to USE-statements,
kind-parameters added to all INTEGER and REAL declaration statements,
kinds are defined in new module kinds,
old module precision_kind is removed,
revision history before 2012 removed,
comment fields (!:) to be used for variable explanations added to all variable declaration statements

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