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

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

last commit documented

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