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

Last change on this file since 2876 was 2817, checked in by knoop, 6 years ago

Preliminary gust module interface implemented

  • Property svn:keywords set to Id
File size: 37.2 KB
Line 
1!> @file sum_up_3d_data.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! 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-2018 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: sum_up_3d_data.f90 2817 2018-02-19 16:32:21Z knoop $
27! Preliminary gust module interface implemented
28!
29! 2798 2018-02-09 17:16:39Z suehring
30! Consider also default-type surfaces for surface temperature output.
31!
32! 2797 2018-02-08 13:24:35Z suehring
33! Enable output of ground-heat flux also at urban surfaces.
34!
35! 2790 2018-02-06 11:57:19Z suehring
36! Bugfix in summation of surface sensible and latent heat flux
37!
38! 2766 2018-01-22 17:17:47Z kanani
39! Removed preprocessor directive __chem
40!
41! 2743 2018-01-12 16:03:39Z suehring
42! In case of natural- and urban-type surfaces output surfaces fluxes in W/m2.
43!
44! 2742 2018-01-12 14:59:47Z suehring
45! Enable output of surface temperature
46!
47! 2735 2018-01-11 12:01:27Z suehring
48! output of r_a moved from land-surface to consider also urban-type surfaces
49!
50! 2718 2018-01-02 08:49:38Z maronga
51! Corrected "Former revisions" section
52!
53! 2696 2017-12-14 17:12:51Z kanani
54! - Change in file header (GPL part)
55! - Implementation of uv exposure model (FK)
56! - output of diss_av, kh_av, km_av (turbulence_closure_mod) (TG)
57! - Implementation of chemistry module (FK)
58! - Workaround for sum-up usm arrays in case of restart runs, to avoid program
59!   crash (MS)
60!
61! 2292 2017-06-20 09:51:42Z schwenkel
62! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
63! includes two more prognostic equations for cloud drop concentration (nc) 
64! and cloud water content (qc).
65!
66! 2233 2017-05-30 18:08:54Z suehring
67!
68! 2232 2017-05-30 17:47:52Z suehring
69! Adjustments to new surface concept
70!
71! 2031 2016-10-21 15:11:58Z knoop
72! renamed variable rho to rho_ocean and rho_av to rho_ocean_av
73!
74! 2024 2016-10-12 16:42:37Z kanani
75! Added missing CASE for ssws*
76!
77! 2011 2016-09-19 17:29:57Z kanani
78! Flag urban_surface is now defined in module control_parameters,
79! changed prefix for urban surface model output to "usm_",
80! introduced control parameter varnamelength for LEN of trimvar.
81!
82! 2007 2016-08-24 15:47:17Z kanani
83! Added support for new urban surface model (temporary modifications of
84! SELECT CASE ( ) necessary, see variable trimvar),
85! added comments in variable declaration section
86!
87! 2000 2016-08-20 18:09:15Z knoop
88! Forced header and separation lines into 80 columns
89!
90! 1992 2016-08-12 15:14:59Z suehring
91! Bugfix in summation of passive scalar
92!
93! 1976 2016-07-27 13:28:04Z maronga
94! Radiation actions are now done directly in the respective module
95!
96! 1972 2016-07-26 07:52:02Z maronga
97! Land surface actions are now done directly in the respective module
98!
99! 1960 2016-07-12 16:34:24Z suehring
100! Scalar surface flux added
101!
102! 1949 2016-06-17 07:19:16Z maronga
103! Bugfix: calculation of lai_av, c_veg_av and c_liq_av.
104!
105! 1849 2016-04-08 11:33:18Z hoffmann
106! precipitation_rate moved to arrays_3d
107!
108! 1788 2016-03-10 11:01:04Z maronga
109! Added z0q and z0q_av
110!
111! 1693 2015-10-27 08:35:45Z maronga
112! Last revision text corrected
113!
114! 1691 2015-10-26 16:17:44Z maronga
115! Added output of Obukhov length and radiative heating rates for RRTMG.
116! Corrected output of liquid water path.
117!
118! 1682 2015-10-07 23:56:08Z knoop
119! Code annotations made doxygen readable
120!
121! 1585 2015-04-30 07:05:52Z maronga
122! Adapted for RRTMG
123!
124! 1555 2015-03-04 17:44:27Z maronga
125! Added output of r_a and r_s
126!
127! 1551 2015-03-03 14:18:16Z maronga
128! Added support for land surface model and radiation model data.
129!
130! 1359 2014-04-11 17:15:14Z hoffmann
131! New particle structure integrated.
132!
133! 1353 2014-04-08 15:21:23Z heinze
134! REAL constants provided with KIND-attribute
135!
136! 1320 2014-03-20 08:40:49Z raasch
137! ONLY-attribute added to USE-statements,
138! kind-parameters added to all INTEGER and REAL declaration statements,
139! kinds are defined in new module kinds,
140! old module precision_kind is removed,
141! revision history before 2012 removed,
142! comment fields (!:) to be used for variable explanations added to
143! all variable declaration statements
144!
145! 1318 2014-03-17 13:35:16Z raasch
146! barrier argument removed from cpu_log,
147! module interfaces removed
148!
149! 1115 2013-03-26 18:16:16Z hoffmann
150! ql is calculated by calc_liquid_water_content
151!
152! 1053 2012-11-13 17:11:03Z hoffmann
153! +nr, prr, qr
154!
155! 1036 2012-10-22 13:43:42Z raasch
156! code put under GPL (PALM 3.9)
157!
158! 1007 2012-09-19 14:30:36Z franke
159! Bugfix in calculation of ql_vp
160!
161! 978 2012-08-09 08:28:32Z fricke
162! +z0h*
163!
164! Revision 1.1  2006/02/23 12:55:23  raasch
165! Initial revision
166!
167!
168! Description:
169! ------------
170!> Sum-up the values of 3d-arrays. The real averaging is later done in routine
171!> average_3d_data.
172!------------------------------------------------------------------------------!
173 SUBROUTINE sum_up_3d_data
174 
175
176    USE arrays_3d,                                                             &
177        ONLY:  dzw, e, heatflux_output_conversion, nc, nr, p, pt,              &
178               precipitation_rate, q, qc, ql, ql_c,                            &
179               ql_v, qr, rho_ocean, s, sa, u, v, vpt, w,                       &
180               waterflux_output_conversion
181
182    USE averaging,                                                             &
183        ONLY:  diss_av, e_av, ghf_av, kh_av, km_av, lpt_av, lwp_av, nc_av,     &
184               nr_av,                                                          &
185               ol_av, p_av, pc_av, pr_av, prr_av, precipitation_rate_av, pt_av,&
186               q_av, qc_av, ql_av, ql_c_av, ql_v_av, ql_vp_av, qr_av, qsws_av, &
187               qv_av, r_a_av, rho_ocean_av, s_av, sa_av, shf_av, ssws_av,      &
188               ts_av, tsurf_av, u_av, us_av, v_av, vpt_av, w_av, z0_av, z0h_av,&
189               z0q_av
190    USE chemistry_model_mod,                                                   &
191        ONLY:  chem_3d_data_averaging, chem_integrate, chem_species, nspec                                   
192
193    USE cloud_parameters,                                                      &
194        ONLY:  cp, l_d_cp, l_v, pt_d_t
195
196    USE control_parameters,                                                    &
197        ONLY:  air_chemistry, average_count_3d, cloud_physics, doav, doav_n,   &
198               land_surface, rho_surface, urban_surface, uv_exposure,          &
199               varnamelength
200
201    USE cpulog,                                                                &
202        ONLY:  cpu_log, log_point
203
204    USE gust_mod,                                                              &
205        ONLY:  gust_3d_data_averaging, gust_module_enabled
206
207    USE indices,                                                               &
208        ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt 
209
210    USE kinds
211
212    USE land_surface_model_mod,                                                &
213        ONLY:  lsm_3d_data_averaging
214
215    USE particle_attributes,                                                   &
216        ONLY:  grid_particles, number_of_particles, particles, prt_count
217
218    USE radiation_model_mod,                                                   &
219        ONLY:  radiation, radiation_3d_data_averaging
220
221    USE surface_mod,                                                           &
222        ONLY:  surf_def_h, surf_lsm_h, surf_usm_h
223
224    USE turbulence_closure_mod,                                                &
225        ONLY:  tcm_3d_data_averaging
226
227    USE urban_surface_mod,                                                     &
228        ONLY:  usm_average_3d_data
229
230    USE uv_exposure_model_mod,                                                &
231        ONLY:  uvem_3d_data_averaging
232
233
234    IMPLICIT NONE
235
236    INTEGER(iwp) ::  i   !< grid index x direction
237    INTEGER(iwp) ::  ii  !< running index
238    INTEGER(iwp) ::  j   !< grid index y direction
239    INTEGER(iwp) ::  k   !< grid index x direction
240    INTEGER(iwp) ::  m   !< running index surface type
241    INTEGER(iwp) ::  n   !<
242
243    REAL(wp)     ::  mean_r !<
244    REAL(wp)     ::  s_r2   !<
245    REAL(wp)     ::  s_r3   !<
246
247    CHARACTER (LEN=varnamelength) ::  trimvar  !< TRIM of output-variable string
248
249
250    CALL cpu_log (log_point(34),'sum_up_3d_data','start')
251
252!
253!-- Allocate and initialize the summation arrays if called for the very first
254!-- time or the first time after average_3d_data has been called
255!-- (some or all of the arrays may have been already allocated
256!-- in read_3d_binary)
257    IF ( average_count_3d == 0 )  THEN
258
259       DO  ii = 1, doav_n
260!
261!--       Temporary solution to account for data output within the new urban
262!--       surface model (urban_surface_mod.f90), see also SELECT CASE ( trimvar )
263          trimvar = TRIM( doav(ii) )
264          IF ( urban_surface  .AND.  trimvar(1:4) == 'usm_' )  THEN
265             trimvar = 'usm_output'
266          ENDIF
267       
268          SELECT CASE ( trimvar )
269
270             CASE ( 'ghf*' )
271                IF ( .NOT. ALLOCATED( ghf_av ) )  THEN
272                   ALLOCATE( ghf_av(nysg:nyng,nxlg:nxrg) )
273                ENDIF
274                ghf_av = 0.0_wp
275
276             CASE ( 'e' )
277                IF ( .NOT. ALLOCATED( e_av ) )  THEN
278                   ALLOCATE( e_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
279                ENDIF
280                e_av = 0.0_wp
281
282             CASE ( 'lpt' )
283                IF ( .NOT. ALLOCATED( lpt_av ) )  THEN
284                   ALLOCATE( lpt_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
285                ENDIF
286                lpt_av = 0.0_wp
287
288             CASE ( 'lwp*' )
289                IF ( .NOT. ALLOCATED( lwp_av ) )  THEN
290                   ALLOCATE( lwp_av(nysg:nyng,nxlg:nxrg) )
291                ENDIF
292                lwp_av = 0.0_wp
293
294             CASE ( 'nc' )
295                IF ( .NOT. ALLOCATED( nc_av ) )  THEN
296                   ALLOCATE( nc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
297                ENDIF
298                nc_av = 0.0_wp
299
300             CASE ( 'nr' )
301                IF ( .NOT. ALLOCATED( nr_av ) )  THEN
302                   ALLOCATE( nr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
303                ENDIF
304                nr_av = 0.0_wp
305
306             CASE ( 'ol*' )
307                IF ( .NOT. ALLOCATED( ol_av ) )  THEN
308                   ALLOCATE( ol_av(nysg:nyng,nxlg:nxrg) )
309                ENDIF
310                ol_av = 0.0_wp
311
312             CASE ( 'p' )
313                IF ( .NOT. ALLOCATED( p_av ) )  THEN
314                   ALLOCATE( p_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
315                ENDIF
316                p_av = 0.0_wp
317
318             CASE ( 'pc' )
319                IF ( .NOT. ALLOCATED( pc_av ) )  THEN
320                   ALLOCATE( pc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
321                ENDIF
322                pc_av = 0.0_wp
323
324             CASE ( 'pr' )
325                IF ( .NOT. ALLOCATED( pr_av ) )  THEN
326                   ALLOCATE( pr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
327                ENDIF
328                pr_av = 0.0_wp
329
330             CASE ( 'prr' )
331                IF ( .NOT. ALLOCATED( prr_av ) )  THEN
332                   ALLOCATE( prr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
333                ENDIF
334                prr_av = 0.0_wp
335
336             CASE ( 'prr*' )
337                IF ( .NOT. ALLOCATED( precipitation_rate_av ) )  THEN
338                   ALLOCATE( precipitation_rate_av(nysg:nyng,nxlg:nxrg) )
339                ENDIF
340                precipitation_rate_av = 0.0_wp
341
342             CASE ( 'pt' )
343                IF ( .NOT. ALLOCATED( pt_av ) )  THEN
344                   ALLOCATE( pt_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
345                ENDIF
346                pt_av = 0.0_wp
347
348             CASE ( 'q' )
349                IF ( .NOT. ALLOCATED( q_av ) )  THEN
350                   ALLOCATE( q_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
351                ENDIF
352                q_av = 0.0_wp
353
354             CASE ( 'qc' )
355                IF ( .NOT. ALLOCATED( qc_av ) )  THEN
356                   ALLOCATE( qc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
357                ENDIF
358                qc_av = 0.0_wp
359
360             CASE ( 'ql' )
361                IF ( .NOT. ALLOCATED( ql_av ) )  THEN
362                   ALLOCATE( ql_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
363                ENDIF
364                ql_av = 0.0_wp
365
366             CASE ( 'ql_c' )
367                IF ( .NOT. ALLOCATED( ql_c_av ) )  THEN
368                   ALLOCATE( ql_c_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
369                ENDIF
370                ql_c_av = 0.0_wp
371
372             CASE ( 'ql_v' )
373                IF ( .NOT. ALLOCATED( ql_v_av ) )  THEN
374                   ALLOCATE( ql_v_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
375                ENDIF
376                ql_v_av = 0.0_wp
377
378             CASE ( 'ql_vp' )
379                IF ( .NOT. ALLOCATED( ql_vp_av ) )  THEN
380                   ALLOCATE( ql_vp_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
381                ENDIF
382                ql_vp_av = 0.0_wp
383
384             CASE ( 'qr' )
385                IF ( .NOT. ALLOCATED( qr_av ) )  THEN
386                   ALLOCATE( qr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
387                ENDIF
388                qr_av = 0.0_wp
389
390             CASE ( 'qsws*' )
391                IF ( .NOT. ALLOCATED( qsws_av ) )  THEN
392                   ALLOCATE( qsws_av(nysg:nyng,nxlg:nxrg) )
393                ENDIF
394                qsws_av = 0.0_wp
395
396             CASE ( 'qv' )
397                IF ( .NOT. ALLOCATED( qv_av ) )  THEN
398                   ALLOCATE( qv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
399                ENDIF
400                qv_av = 0.0_wp
401
402             CASE ( 'r_a*' )
403                IF ( .NOT. ALLOCATED( r_a_av ) )  THEN
404                   ALLOCATE( r_a_av(nysg:nyng,nxlg:nxrg) )
405                ENDIF
406                r_a_av = 0.0_wp
407
408             CASE ( 'rho_ocean' )
409                IF ( .NOT. ALLOCATED( rho_ocean_av ) )  THEN
410                   ALLOCATE( rho_ocean_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
411                ENDIF
412                rho_ocean_av = 0.0_wp
413
414             CASE ( 's' )
415                IF ( .NOT. ALLOCATED( s_av ) )  THEN
416                   ALLOCATE( s_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
417                ENDIF
418                s_av = 0.0_wp
419
420             CASE ( 'sa' )
421                IF ( .NOT. ALLOCATED( sa_av ) )  THEN
422                   ALLOCATE( sa_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
423                ENDIF
424                sa_av = 0.0_wp
425
426             CASE ( 'shf*' )
427                IF ( .NOT. ALLOCATED( shf_av ) )  THEN
428                   ALLOCATE( shf_av(nysg:nyng,nxlg:nxrg) )
429                ENDIF
430                shf_av = 0.0_wp
431               
432             CASE ( 'ssws*' )
433                IF ( .NOT. ALLOCATED( ssws_av ) )  THEN
434                   ALLOCATE( ssws_av(nysg:nyng,nxlg:nxrg) )
435                ENDIF
436                ssws_av = 0.0_wp               
437
438             CASE ( 't*' )
439                IF ( .NOT. ALLOCATED( ts_av ) )  THEN
440                   ALLOCATE( ts_av(nysg:nyng,nxlg:nxrg) )
441                ENDIF
442                ts_av = 0.0_wp
443
444             CASE ( 'tsurf*' )
445                IF ( .NOT. ALLOCATED( tsurf_av ) )  THEN
446                   ALLOCATE( tsurf_av(nysg:nyng,nxlg:nxrg) )
447                ENDIF
448                tsurf_av = 0.0_wp
449
450             CASE ( 'u' )
451                IF ( .NOT. ALLOCATED( u_av ) )  THEN
452                   ALLOCATE( u_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
453                ENDIF
454                u_av = 0.0_wp
455
456             CASE ( 'u*' )
457                IF ( .NOT. ALLOCATED( us_av ) )  THEN
458                   ALLOCATE( us_av(nysg:nyng,nxlg:nxrg) )
459                ENDIF
460                us_av = 0.0_wp
461
462             CASE ( 'v' )
463                IF ( .NOT. ALLOCATED( v_av ) )  THEN
464                   ALLOCATE( v_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
465                ENDIF
466                v_av = 0.0_wp
467
468             CASE ( 'vpt' )
469                IF ( .NOT. ALLOCATED( vpt_av ) )  THEN
470                   ALLOCATE( vpt_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
471                ENDIF
472                vpt_av = 0.0_wp
473
474             CASE ( 'w' )
475                IF ( .NOT. ALLOCATED( w_av ) )  THEN
476                   ALLOCATE( w_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
477                ENDIF
478                w_av = 0.0_wp
479
480             CASE ( 'z0*' )
481                IF ( .NOT. ALLOCATED( z0_av ) )  THEN
482                   ALLOCATE( z0_av(nysg:nyng,nxlg:nxrg) )
483                ENDIF
484                z0_av = 0.0_wp
485
486             CASE ( 'z0h*' )
487                IF ( .NOT. ALLOCATED( z0h_av ) )  THEN
488                   ALLOCATE( z0h_av(nysg:nyng,nxlg:nxrg) )
489                ENDIF
490                z0h_av = 0.0_wp
491
492             CASE ( 'z0q*' )
493                IF ( .NOT. ALLOCATED( z0q_av ) )  THEN
494                   ALLOCATE( z0q_av(nysg:nyng,nxlg:nxrg) )
495                ENDIF
496                z0q_av = 0.0_wp
497!             
498!--          Block of urban surface model outputs
499             CASE ( 'usm_output' )
500
501                CALL usm_average_3d_data( 'allocate', doav(ii) )
502             
503
504             CASE DEFAULT
505
506!
507!--             Turbulence closure module
508                CALL tcm_3d_data_averaging( 'allocate', doav(ii) )
509
510!
511!--             Land surface quantity
512                IF ( land_surface )  THEN
513                   CALL lsm_3d_data_averaging( 'allocate', doav(ii) )
514                ENDIF
515
516!
517!--             Radiation quantity
518                IF ( radiation )  THEN
519                   CALL radiation_3d_data_averaging( 'allocate', doav(ii) )
520                ENDIF
521
522!
523!--             Gust module quantities
524                IF ( gust_module_enabled )  THEN
525                   CALL gust_3d_data_averaging( 'allocate', doav(ii) )
526                ENDIF
527
528!
529!--             Chemical quantity                                           
530#if defined( __chem )               
531                IF ( air_chemistry  .AND.  trimvar(1:3) == 'kc_')  THEN
532                   CALL chem_3d_data_averaging( 'allocate', doav(ii) )
533                ENDIF
534#endif
535
536!
537!--             UV exposure quantity
538                IF ( uv_exposure  .AND.  trimvar(1:5) == 'uvem_')  THEN
539                   CALL uvem_3d_data_averaging( 'allocate', doav(ii) )
540                ENDIF
541
542!
543!--             User-defined quantity
544                CALL user_3d_data_averaging( 'allocate', doav(ii) )
545
546          END SELECT
547
548       ENDDO
549
550    ENDIF
551
552!
553!-- Loop of all variables to be averaged.
554    DO  ii = 1, doav_n
555!
556!--       Temporary solution to account for data output within the new urban
557!--       surface model (urban_surface_mod.f90), see also SELECT CASE ( trimvar )
558          trimvar = TRIM( doav(ii) )
559          IF ( urban_surface  .AND.  trimvar(1:4) == 'usm_' )  THEN
560             trimvar = 'usm_output'
561          ENDIF
562!
563!--    Store the array chosen on the temporary array.
564       SELECT CASE ( trimvar )
565
566          CASE ( 'ghf*' )
567             DO  m = 1, surf_lsm_h%ns
568                i   = surf_lsm_h%i(m)           
569                j   = surf_lsm_h%j(m)
570                ghf_av(j,i) = ghf_av(j,i) + surf_lsm_h%ghf(m)
571             ENDDO
572
573             DO  m = 1, surf_usm_h%ns
574                i   = surf_usm_h%i(m)           
575                j   = surf_usm_h%j(m)
576                ghf_av(j,i) = ghf_av(j,i) + surf_usm_h%frac(0,m)     *          &
577                                            surf_usm_h%wghf_eb(m)        +      &
578                                            surf_usm_h%frac(1,m)     *          &
579                                            surf_usm_h%wghf_eb_green(m)  +      &
580                                            surf_usm_h%frac(2,m)     *          &
581                                            surf_usm_h%wghf_eb_window(m)
582             ENDDO
583
584          CASE ( 'e' )
585             DO  i = nxlg, nxrg
586                DO  j = nysg, nyng
587                   DO  k = nzb, nzt+1
588                      e_av(k,j,i) = e_av(k,j,i) + e(k,j,i)
589                   ENDDO
590                ENDDO
591             ENDDO
592
593          CASE ( 'lpt' )
594             DO  i = nxlg, nxrg
595                DO  j = nysg, nyng
596                   DO  k = nzb, nzt+1
597                      lpt_av(k,j,i) = lpt_av(k,j,i) + pt(k,j,i)
598                   ENDDO
599                ENDDO
600             ENDDO
601
602          CASE ( 'lwp*' )
603             DO  i = nxlg, nxrg
604                DO  j = nysg, nyng
605                   lwp_av(j,i) = lwp_av(j,i) + SUM( ql(nzb:nzt,j,i)            &
606                                               * dzw(1:nzt+1) ) * rho_surface
607                ENDDO
608             ENDDO
609
610          CASE ( 'nc' )
611             DO  i = nxlg, nxrg
612                DO  j = nysg, nyng
613                   DO  k = nzb, nzt+1
614                      nc_av(k,j,i) = nc_av(k,j,i) + nc(k,j,i)
615                   ENDDO
616                ENDDO
617             ENDDO
618
619          CASE ( 'nr' )
620             DO  i = nxlg, nxrg
621                DO  j = nysg, nyng
622                   DO  k = nzb, nzt+1
623                      nr_av(k,j,i) = nr_av(k,j,i) + nr(k,j,i)
624                   ENDDO
625                ENDDO
626             ENDDO
627
628          CASE ( 'ol*' )
629             DO  m = 1, surf_def_h(0)%ns
630                i = surf_def_h(0)%i(m)
631                j = surf_def_h(0)%j(m)
632                ol_av(j,i) = ol_av(j,i) + surf_def_h(0)%ol(m)
633             ENDDO
634             DO  m = 1, surf_lsm_h%ns
635                i = surf_lsm_h%i(m)
636                j = surf_lsm_h%j(m)
637                ol_av(j,i) = ol_av(j,i) + surf_lsm_h%ol(m)
638             ENDDO
639             DO  m = 1, surf_usm_h%ns
640                i = surf_usm_h%i(m)
641                j = surf_usm_h%j(m)
642                ol_av(j,i) = ol_av(j,i) + surf_usm_h%ol(m)
643             ENDDO
644
645          CASE ( 'p' )
646             DO  i = nxlg, nxrg
647                DO  j = nysg, nyng
648                   DO  k = nzb, nzt+1
649                      p_av(k,j,i) = p_av(k,j,i) + p(k,j,i)
650                   ENDDO
651                ENDDO
652             ENDDO
653
654          CASE ( 'pc' )
655             DO  i = nxl, nxr
656                DO  j = nys, nyn
657                   DO  k = nzb, nzt+1
658                      pc_av(k,j,i) = pc_av(k,j,i) + prt_count(k,j,i)
659                   ENDDO
660                ENDDO
661             ENDDO
662
663          CASE ( 'pr' )
664             DO  i = nxl, nxr
665                DO  j = nys, nyn
666                   DO  k = nzb, nzt+1
667                      number_of_particles = prt_count(k,j,i)
668                      IF ( number_of_particles <= 0 )  CYCLE
669                      particles => grid_particles(k,j,i)%particles(1:number_of_particles)
670                      s_r2 = 0.0_wp
671                      s_r3 = 0.0_wp
672
673                      DO  n = 1, number_of_particles
674                         IF ( particles(n)%particle_mask )  THEN
675                            s_r2 = s_r2 + particles(n)%radius**2 * &
676                                particles(n)%weight_factor
677                            s_r3 = s_r3 + particles(n)%radius**3 * &
678                                particles(n)%weight_factor
679                         ENDIF
680                      ENDDO
681
682                      IF ( s_r2 > 0.0_wp )  THEN
683                         mean_r = s_r3 / s_r2
684                      ELSE
685                         mean_r = 0.0_wp
686                      ENDIF
687                      pr_av(k,j,i) = pr_av(k,j,i) + mean_r
688                   ENDDO
689                ENDDO
690             ENDDO
691
692
693          CASE ( 'pr*' )
694             DO  i = nxlg, nxrg
695                DO  j = nysg, nyng
696                   precipitation_rate_av(j,i) = precipitation_rate_av(j,i) + &
697                                                precipitation_rate(j,i)
698                ENDDO
699             ENDDO
700
701          CASE ( 'pt' )
702             IF ( .NOT. cloud_physics ) THEN
703             DO  i = nxlg, nxrg
704                DO  j = nysg, nyng
705                   DO  k = nzb, nzt+1
706                         pt_av(k,j,i) = pt_av(k,j,i) + pt(k,j,i)
707                      ENDDO
708                   ENDDO
709                ENDDO
710             ELSE
711             DO  i = nxlg, nxrg
712                DO  j = nysg, nyng
713                   DO  k = nzb, nzt+1
714                         pt_av(k,j,i) = pt_av(k,j,i) + pt(k,j,i) + l_d_cp * &
715                                                       pt_d_t(k) * ql(k,j,i)
716                      ENDDO
717                   ENDDO
718                ENDDO
719             ENDIF
720
721          CASE ( 'q' )
722             DO  i = nxlg, nxrg
723                DO  j = nysg, nyng
724                   DO  k = nzb, nzt+1
725                      q_av(k,j,i) = q_av(k,j,i) + q(k,j,i)
726                   ENDDO
727                ENDDO
728             ENDDO
729
730          CASE ( 'qc' )
731             DO  i = nxlg, nxrg
732                DO  j = nysg, nyng
733                   DO  k = nzb, nzt+1
734                      qc_av(k,j,i) = qc_av(k,j,i) + qc(k,j,i)
735                   ENDDO
736                ENDDO
737             ENDDO
738
739          CASE ( 'ql' )
740             DO  i = nxlg, nxrg
741                DO  j = nysg, nyng
742                   DO  k = nzb, nzt+1
743                      ql_av(k,j,i) = ql_av(k,j,i) + ql(k,j,i)
744                   ENDDO
745                ENDDO
746             ENDDO
747
748          CASE ( 'ql_c' )
749             DO  i = nxlg, nxrg
750                DO  j = nysg, nyng
751                   DO  k = nzb, nzt+1
752                      ql_c_av(k,j,i) = ql_c_av(k,j,i) + ql_c(k,j,i)
753                   ENDDO
754                ENDDO
755             ENDDO
756
757          CASE ( 'ql_v' )
758             DO  i = nxlg, nxrg
759                DO  j = nysg, nyng
760                   DO  k = nzb, nzt+1
761                      ql_v_av(k,j,i) = ql_v_av(k,j,i) + ql_v(k,j,i)
762                   ENDDO
763                ENDDO
764             ENDDO
765
766          CASE ( 'ql_vp' )
767             DO  i = nxl, nxr
768                DO  j = nys, nyn
769                   DO  k = nzb, nzt+1
770                      number_of_particles = prt_count(k,j,i)
771                      IF ( number_of_particles <= 0 )  CYCLE
772                      particles => grid_particles(k,j,i)%particles(1:number_of_particles)
773                      DO  n = 1, number_of_particles
774                         IF ( particles(n)%particle_mask )  THEN
775                            ql_vp_av(k,j,i) = ql_vp_av(k,j,i) + &
776                                              particles(n)%weight_factor / &
777                                              number_of_particles
778                         ENDIF
779                      ENDDO
780                   ENDDO
781                ENDDO
782             ENDDO
783
784          CASE ( 'qr' )
785             DO  i = nxlg, nxrg
786                DO  j = nysg, nyng
787                   DO  k = nzb, nzt+1
788                      qr_av(k,j,i) = qr_av(k,j,i) + qr(k,j,i)
789                   ENDDO
790                ENDDO
791             ENDDO
792
793          CASE ( 'qsws*' )
794!
795!--          In case of default surfaces, clean-up flux by density.
796!--          In case of land- and urban-surfaces, convert fluxes into
797!--          dynamic units.
798             DO  m = 1, surf_def_h(0)%ns
799                i = surf_def_h(0)%i(m)
800                j = surf_def_h(0)%j(m)
801                k = surf_def_h(0)%k(m)
802                qsws_av(j,i) = qsws_av(j,i) + surf_def_h(0)%qsws(m) *          &
803                                              waterflux_output_conversion(k)
804             ENDDO
805             DO  m = 1, surf_lsm_h%ns
806                i = surf_lsm_h%i(m)
807                j = surf_lsm_h%j(m)
808                qsws_av(j,i) = qsws_av(j,i) + surf_lsm_h%qsws(m) * l_v
809             ENDDO
810             DO  m = 1, surf_usm_h%ns
811                i = surf_usm_h%i(m)
812                j = surf_usm_h%j(m)
813                qsws_av(j,i) = qsws_av(j,i) + surf_usm_h%qsws(m) * l_v
814             ENDDO
815
816          CASE ( 'qv' )
817             DO  i = nxlg, nxrg
818                DO  j = nysg, nyng
819                   DO  k = nzb, nzt+1
820                      qv_av(k,j,i) = qv_av(k,j,i) + q(k,j,i) - ql(k,j,i)
821                   ENDDO
822                ENDDO
823             ENDDO
824
825          CASE ( 'r_a*' )
826             DO  m = 1, surf_lsm_h%ns
827                i   = surf_lsm_h%i(m)           
828                j   = surf_lsm_h%j(m)
829                r_a_av(j,i) = r_a_av(j,i) + surf_lsm_h%r_a(m)
830             ENDDO
831!
832!--          Please note, resistance is also applied at urban-type surfaces,
833!--          and is output only as a single variable. Here, tile approach is
834!--          already implemented, so for each surface fraction resistance
835!--          need to be summed-up.
836             DO  m = 1, surf_usm_h%ns
837                i   = surf_usm_h%i(m)           
838                j   = surf_usm_h%j(m)
839                r_a_av(j,i) = r_a_av(j,i) +                                    &
840                           ( surf_usm_h%frac(0,m) * surf_usm_h%r_a(m)       +  & 
841                             surf_usm_h%frac(1,m) * surf_usm_h%r_a_green(m) +  & 
842                             surf_usm_h%frac(2,m) * surf_usm_h%r_a_window(m) )
843             ENDDO
844
845          CASE ( 'rho_ocean' )
846             DO  i = nxlg, nxrg
847                DO  j = nysg, nyng
848                   DO  k = nzb, nzt+1
849                      rho_ocean_av(k,j,i) = rho_ocean_av(k,j,i) + rho_ocean(k,j,i)
850                   ENDDO
851                ENDDO
852             ENDDO
853
854          CASE ( 's' )
855             DO  i = nxlg, nxrg
856                DO  j = nysg, nyng
857                   DO  k = nzb, nzt+1
858                      s_av(k,j,i) = s_av(k,j,i) + s(k,j,i)
859                   ENDDO
860                ENDDO
861             ENDDO
862
863          CASE ( 'sa' )
864             DO  i = nxlg, nxrg
865                DO  j = nysg, nyng
866                   DO  k = nzb, nzt+1
867                      sa_av(k,j,i) = sa_av(k,j,i) + sa(k,j,i)
868                   ENDDO
869                ENDDO
870             ENDDO
871
872          CASE ( 'shf*' )
873!
874!--          In case of default surfaces, clean-up flux by density.
875!--          In case of land- and urban-surfaces, convert fluxes into
876!--          dynamic units.
877             DO  m = 1, surf_def_h(0)%ns
878                i = surf_def_h(0)%i(m)
879                j = surf_def_h(0)%j(m)
880                k = surf_def_h(0)%k(m)
881                shf_av(j,i) = shf_av(j,i) + surf_def_h(0)%shf(m)  *            &
882                                            heatflux_output_conversion(k)
883             ENDDO
884             DO  m = 1, surf_lsm_h%ns
885                i = surf_lsm_h%i(m)
886                j = surf_lsm_h%j(m)
887                shf_av(j,i) = shf_av(j,i) + surf_lsm_h%shf(m) * cp
888             ENDDO
889             DO  m = 1, surf_usm_h%ns
890                i = surf_usm_h%i(m)
891                j = surf_usm_h%j(m)
892                shf_av(j,i) = shf_av(j,i) + surf_usm_h%shf(m) * cp
893             ENDDO
894
895          CASE ( 'ssws*' )
896             DO  m = 1, surf_def_h(0)%ns
897                i = surf_def_h(0)%i(m)
898                j = surf_def_h(0)%j(m)
899                ssws_av(j,i) = ssws_av(j,i) + surf_def_h(0)%ssws(m)
900             ENDDO
901             DO  m = 1, surf_lsm_h%ns
902                i = surf_lsm_h%i(m)
903                j = surf_lsm_h%j(m)
904                ssws_av(j,i) = ssws_av(j,i) + surf_lsm_h%ssws(m)
905             ENDDO
906             DO  m = 1, surf_usm_h%ns
907                i = surf_usm_h%i(m)
908                j = surf_usm_h%j(m)
909                ssws_av(j,i) = ssws_av(j,i) + surf_usm_h%ssws(m)
910             ENDDO
911
912          CASE ( 't*' )
913             DO  m = 1, surf_def_h(0)%ns
914                i = surf_def_h(0)%i(m)
915                j = surf_def_h(0)%j(m)
916                ts_av(j,i) = ts_av(j,i) + surf_def_h(0)%ts(m)
917             ENDDO
918             DO  m = 1, surf_lsm_h%ns
919                i = surf_lsm_h%i(m)
920                j = surf_lsm_h%j(m)
921                ts_av(j,i) = ts_av(j,i) + surf_lsm_h%ts(m)
922             ENDDO
923             DO  m = 1, surf_usm_h%ns
924                i = surf_usm_h%i(m)
925                j = surf_usm_h%j(m)
926                ts_av(j,i) = ts_av(j,i) + surf_usm_h%ts(m)
927             ENDDO
928
929          CASE ( 'tsurf*' )
930             DO  m = 1, surf_def_h(0)%ns
931                i   = surf_def_h(0)%i(m)           
932                j   = surf_def_h(0)%j(m)
933                tsurf_av(j,i) = tsurf_av(j,i) + surf_def_h(0)%pt_surface(m)
934             ENDDO
935
936             DO  m = 1, surf_lsm_h%ns
937                i   = surf_lsm_h%i(m)           
938                j   = surf_lsm_h%j(m)
939                tsurf_av(j,i) = tsurf_av(j,i) + surf_lsm_h%pt_surface(m)
940             ENDDO
941
942             DO  m = 1, surf_usm_h%ns
943                i   = surf_usm_h%i(m)           
944                j   = surf_usm_h%j(m)
945                tsurf_av(j,i) = tsurf_av(j,i) + surf_usm_h%pt_surface(m)
946             ENDDO
947
948          CASE ( 'u' )
949             DO  i = nxlg, nxrg
950                DO  j = nysg, nyng
951                   DO  k = nzb, nzt+1
952                      u_av(k,j,i) = u_av(k,j,i) + u(k,j,i)
953                   ENDDO
954                ENDDO
955             ENDDO
956
957          CASE ( 'u*' )
958             DO  m = 1, surf_def_h(0)%ns
959                i = surf_def_h(0)%i(m)
960                j = surf_def_h(0)%j(m)
961                us_av(j,i) = us_av(j,i) + surf_def_h(0)%us(m)
962             ENDDO
963             DO  m = 1, surf_lsm_h%ns
964                i = surf_lsm_h%i(m)
965                j = surf_lsm_h%j(m)
966                us_av(j,i) = us_av(j,i) + surf_lsm_h%us(m)
967             ENDDO
968             DO  m = 1, surf_usm_h%ns
969                i = surf_usm_h%i(m)
970                j = surf_usm_h%j(m)
971                us_av(j,i) = us_av(j,i) + surf_usm_h%us(m)
972             ENDDO
973
974          CASE ( 'v' )
975             DO  i = nxlg, nxrg
976                DO  j = nysg, nyng
977                   DO  k = nzb, nzt+1
978                      v_av(k,j,i) = v_av(k,j,i) + v(k,j,i)
979                   ENDDO
980                ENDDO
981             ENDDO
982
983          CASE ( 'vpt' )
984             DO  i = nxlg, nxrg
985                DO  j = nysg, nyng
986                   DO  k = nzb, nzt+1
987                      vpt_av(k,j,i) = vpt_av(k,j,i) + vpt(k,j,i)
988                   ENDDO
989                ENDDO
990             ENDDO
991
992          CASE ( 'w' )
993             DO  i = nxlg, nxrg
994                DO  j = nysg, nyng
995                   DO  k = nzb, nzt+1
996                      w_av(k,j,i) = w_av(k,j,i) + w(k,j,i)
997                   ENDDO
998                ENDDO
999             ENDDO
1000
1001          CASE ( 'z0*' )
1002             DO  m = 1, surf_def_h(0)%ns
1003                i = surf_def_h(0)%i(m)
1004                j = surf_def_h(0)%j(m)
1005                z0_av(j,i) = z0_av(j,i) + surf_def_h(0)%z0(m)
1006             ENDDO
1007             DO  m = 1, surf_lsm_h%ns
1008                i = surf_lsm_h%i(m)
1009                j = surf_lsm_h%j(m)
1010                z0_av(j,i) = z0_av(j,i) + surf_lsm_h%z0(m)
1011             ENDDO
1012             DO  m = 1, surf_usm_h%ns
1013                i = surf_usm_h%i(m)
1014                j = surf_usm_h%j(m)
1015                z0_av(j,i) = z0_av(j,i) + surf_usm_h%z0(m)
1016             ENDDO
1017
1018          CASE ( 'z0h*' )
1019             DO  m = 1, surf_def_h(0)%ns
1020                i = surf_def_h(0)%i(m)
1021                j = surf_def_h(0)%j(m)
1022                z0h_av(j,i) = z0h_av(j,i) + surf_def_h(0)%z0h(m)
1023             ENDDO
1024             DO  m = 1, surf_lsm_h%ns
1025                i = surf_lsm_h%i(m)
1026                j = surf_lsm_h%j(m)
1027                z0h_av(j,i) = z0h_av(j,i) + surf_lsm_h%z0h(m)
1028             ENDDO
1029             DO  m = 1, surf_usm_h%ns
1030                i = surf_usm_h%i(m)
1031                j = surf_usm_h%j(m)
1032                z0h_av(j,i) = z0h_av(j,i) + surf_usm_h%z0h(m)
1033             ENDDO
1034
1035          CASE ( 'z0q*' )
1036             DO  m = 1, surf_def_h(0)%ns
1037                i = surf_def_h(0)%i(m)
1038                j = surf_def_h(0)%j(m)
1039                z0q_av(j,i) = z0q_av(j,i) + surf_def_h(0)%z0q(m)
1040             ENDDO
1041             DO  m = 1, surf_lsm_h%ns
1042                i = surf_lsm_h%i(m)
1043                j = surf_lsm_h%j(m)
1044                z0q_av(j,i) = z0q_av(j,i) + surf_lsm_h%z0q(m)
1045             ENDDO
1046             DO  m = 1, surf_usm_h%ns
1047                i = surf_usm_h%i(m)
1048                j = surf_usm_h%j(m)
1049                z0q_av(j,i) = z0q_av(j,i) + surf_usm_h%z0q(m)
1050             ENDDO
1051!             
1052!--       Block of urban surface model outputs.
1053!--       In case of urban surface variables it should be always checked
1054!--       if respective arrays are allocated, at least in case of a restart
1055!--       run, as usm arrays are not read from file at the moment.
1056          CASE ( 'usm_output' )
1057             CALL usm_average_3d_data( 'allocate', doav(ii) )
1058             CALL usm_average_3d_data( 'sum', doav(ii) )
1059
1060          CASE DEFAULT
1061!
1062!--          Turbulence closure module
1063             CALL tcm_3d_data_averaging( 'sum', doav(ii) )
1064
1065!
1066!--          Land surface quantity
1067             IF ( land_surface )  THEN
1068                CALL lsm_3d_data_averaging( 'sum', doav(ii) )
1069             ENDIF
1070
1071!
1072!--          Radiation quantity
1073             IF ( radiation )  THEN
1074                CALL radiation_3d_data_averaging( 'sum', doav(ii) )
1075             ENDIF
1076
1077!
1078!--          Gust module quantities
1079             IF ( gust_module_enabled )  THEN
1080                CALL gust_3d_data_averaging( 'sum', doav(ii) )
1081             ENDIF
1082
1083!
1084!--          Chemical quantity
1085             IF ( air_chemistry  .AND.  trimvar(1:3) == 'kc_')  THEN
1086                CALL chem_3d_data_averaging( 'sum',doav(ii) )
1087             ENDIF
1088
1089!
1090!--          UV exposure quantity
1091             IF ( uv_exposure )  THEN
1092                CALL uvem_3d_data_averaging( 'sum', doav(ii) )
1093             ENDIF
1094
1095!
1096!--          User-defined quantity
1097             CALL user_3d_data_averaging( 'sum', doav(ii) )
1098
1099       END SELECT
1100
1101    ENDDO
1102
1103    CALL cpu_log( log_point(34), 'sum_up_3d_data', 'stop' )
1104
1105
1106 END SUBROUTINE sum_up_3d_data
Note: See TracBrowser for help on using the repository browser.