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
RevLine 
[1682]1!> @file sum_up_3d_data.f90
[2000]2!------------------------------------------------------------------------------!
[2696]3! This file is part of the PALM model system.
[1036]4!
[2000]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.
[1036]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!
[2718]17! Copyright 1997-2018 Leibniz Universitaet Hannover
[2000]18!------------------------------------------------------------------------------!
[1036]19!
[484]20! Current revisions:
[1]21! -----------------
[1360]22!
[2233]23!
[1321]24! Former revisions:
25! -----------------
26! $Id: sum_up_3d_data.f90 2817 2018-02-19 16:32:21Z knoop $
[2817]27! Preliminary gust module interface implemented
28!
29! 2798 2018-02-09 17:16:39Z suehring
[2798]30! Consider also default-type surfaces for surface temperature output.
31!
32! 2797 2018-02-08 13:24:35Z suehring
[2797]33! Enable output of ground-heat flux also at urban surfaces.
34!
35! 2790 2018-02-06 11:57:19Z suehring
[2790]36! Bugfix in summation of surface sensible and latent heat flux
37!
38! 2766 2018-01-22 17:17:47Z kanani
[2766]39! Removed preprocessor directive __chem
40!
41! 2743 2018-01-12 16:03:39Z suehring
[2743]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
[2742]45! Enable output of surface temperature
46!
47! 2735 2018-01-11 12:01:27Z suehring
[2735]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
[2716]51! Corrected "Former revisions" section
52!
53! 2696 2017-12-14 17:12:51Z kanani
54! - Change in file header (GPL part)
[2696]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
[2292]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
[1321]67!
[2233]68! 2232 2017-05-30 17:47:52Z suehring
69! Adjustments to new surface concept
70!
[2032]71! 2031 2016-10-21 15:11:58Z knoop
72! renamed variable rho to rho_ocean and rho_av to rho_ocean_av
73!
[2025]74! 2024 2016-10-12 16:42:37Z kanani
75! Added missing CASE for ssws*
76!
[2012]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!
[2008]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!
[2001]87! 2000 2016-08-20 18:09:15Z knoop
88! Forced header and separation lines into 80 columns
89!
[1993]90! 1992 2016-08-12 15:14:59Z suehring
91! Bugfix in summation of passive scalar
92!
[1977]93! 1976 2016-07-27 13:28:04Z maronga
94! Radiation actions are now done directly in the respective module
95!
[1973]96! 1972 2016-07-26 07:52:02Z maronga
97! Land surface actions are now done directly in the respective module
98!
[1961]99! 1960 2016-07-12 16:34:24Z suehring
100! Scalar surface flux added
101!
[1950]102! 1949 2016-06-17 07:19:16Z maronga
103! Bugfix: calculation of lai_av, c_veg_av and c_liq_av.
104!
[1851]105! 1849 2016-04-08 11:33:18Z hoffmann
106! precipitation_rate moved to arrays_3d
[1852]107!
[1789]108! 1788 2016-03-10 11:01:04Z maronga
109! Added z0q and z0q_av
110!
[1694]111! 1693 2015-10-27 08:35:45Z maronga
112! Last revision text corrected
113!
[1692]114! 1691 2015-10-26 16:17:44Z maronga
115! Added output of Obukhov length and radiative heating rates for RRTMG.
[1693]116! Corrected output of liquid water path.
[1692]117!
[1683]118! 1682 2015-10-07 23:56:08Z knoop
119! Code annotations made doxygen readable
120!
[1586]121! 1585 2015-04-30 07:05:52Z maronga
122! Adapted for RRTMG
123!
[1556]124! 1555 2015-03-04 17:44:27Z maronga
125! Added output of r_a and r_s
126!
[1552]127! 1551 2015-03-03 14:18:16Z maronga
128! Added support for land surface model and radiation model data.
129!
[1360]130! 1359 2014-04-11 17:15:14Z hoffmann
131! New particle structure integrated.
132!
[1354]133! 1353 2014-04-08 15:21:23Z heinze
134! REAL constants provided with KIND-attribute
135!
[1321]136! 1320 2014-03-20 08:40:49Z raasch
[1320]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
[1]144!
[1319]145! 1318 2014-03-17 13:35:16Z raasch
146! barrier argument removed from cpu_log,
147! module interfaces removed
148!
[1116]149! 1115 2013-03-26 18:16:16Z hoffmann
150! ql is calculated by calc_liquid_water_content
151!
[1054]152! 1053 2012-11-13 17:11:03Z hoffmann
153! +nr, prr, qr
154!
[1037]155! 1036 2012-10-22 13:43:42Z raasch
156! code put under GPL (PALM 3.9)
157!
[1008]158! 1007 2012-09-19 14:30:36Z franke
159! Bugfix in calculation of ql_vp
160!
[979]161! 978 2012-08-09 08:28:32Z fricke
162! +z0h*
163!
[1]164! Revision 1.1  2006/02/23 12:55:23  raasch
165! Initial revision
166!
167!
168! Description:
169! ------------
[1682]170!> Sum-up the values of 3d-arrays. The real averaging is later done in routine
171!> average_3d_data.
[1]172!------------------------------------------------------------------------------!
[1682]173 SUBROUTINE sum_up_3d_data
174 
[1]175
[1320]176    USE arrays_3d,                                                             &
[2743]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
[1]181
[1320]182    USE averaging,                                                             &
[2797]183        ONLY:  diss_av, e_av, ghf_av, kh_av, km_av, lpt_av, lwp_av, nc_av,     &
184               nr_av,                                                          &
[2696]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, &
[2735]187               qv_av, r_a_av, rho_ocean_av, s_av, sa_av, shf_av, ssws_av,      &
[2742]188               ts_av, tsurf_av, u_av, us_av, v_av, vpt_av, w_av, z0_av, z0h_av,&
189               z0q_av
[2696]190    USE chemistry_model_mod,                                                   &
191        ONLY:  chem_3d_data_averaging, chem_integrate, chem_species, nspec                                   
[1320]192
193    USE cloud_parameters,                                                      &
[2743]194        ONLY:  cp, l_d_cp, l_v, pt_d_t
[1320]195
196    USE control_parameters,                                                    &
[2696]197        ONLY:  air_chemistry, average_count_3d, cloud_physics, doav, doav_n,   &
198               land_surface, rho_surface, urban_surface, uv_exposure,          &
199               varnamelength
[1320]200
201    USE cpulog,                                                                &
202        ONLY:  cpu_log, log_point
203
[2817]204    USE gust_mod,                                                              &
205        ONLY:  gust_3d_data_averaging, gust_module_enabled
206
[1320]207    USE indices,                                                               &
208        ONLY:  nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt 
209
210    USE kinds
211
[1551]212    USE land_surface_model_mod,                                                &
[2232]213        ONLY:  lsm_3d_data_averaging
[1551]214
[1320]215    USE particle_attributes,                                                   &
[1359]216        ONLY:  grid_particles, number_of_particles, particles, prt_count
[1320]217
[1551]218    USE radiation_model_mod,                                                   &
[1976]219        ONLY:  radiation, radiation_3d_data_averaging
[1551]220
[2232]221    USE surface_mod,                                                           &
222        ONLY:  surf_def_h, surf_lsm_h, surf_usm_h
223
[2696]224    USE turbulence_closure_mod,                                                &
225        ONLY:  tcm_3d_data_averaging
226
[2007]227    USE urban_surface_mod,                                                     &
[2011]228        ONLY:  usm_average_3d_data
[1691]229
[2696]230    USE uv_exposure_model_mod,                                                &
231        ONLY:  uvem_3d_data_averaging
[2007]232
[2696]233
[1]234    IMPLICIT NONE
235
[2232]236    INTEGER(iwp) ::  i   !< grid index x direction
[2007]237    INTEGER(iwp) ::  ii  !< running index
[2232]238    INTEGER(iwp) ::  j   !< grid index y direction
239    INTEGER(iwp) ::  k   !< grid index x direction
240    INTEGER(iwp) ::  m   !< running index surface type
[1682]241    INTEGER(iwp) ::  n   !<
[1]242
[1682]243    REAL(wp)     ::  mean_r !<
244    REAL(wp)     ::  s_r2   !<
245    REAL(wp)     ::  s_r3   !<
[1]246
[2011]247    CHARACTER (LEN=varnamelength) ::  trimvar  !< TRIM of output-variable string
[2007]248
249
[1]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
[2007]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) )
[2011]264          IF ( urban_surface  .AND.  trimvar(1:4) == 'usm_' )  THEN
[2007]265             trimvar = 'usm_output'
266          ENDIF
267       
268          SELECT CASE ( trimvar )
[1]269
[2797]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
[1]276             CASE ( 'e' )
277                IF ( .NOT. ALLOCATED( e_av ) )  THEN
[667]278                   ALLOCATE( e_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1]279                ENDIF
[1353]280                e_av = 0.0_wp
[1]281
[771]282             CASE ( 'lpt' )
283                IF ( .NOT. ALLOCATED( lpt_av ) )  THEN
284                   ALLOCATE( lpt_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
285                ENDIF
[1353]286                lpt_av = 0.0_wp
[771]287
[1]288             CASE ( 'lwp*' )
289                IF ( .NOT. ALLOCATED( lwp_av ) )  THEN
[667]290                   ALLOCATE( lwp_av(nysg:nyng,nxlg:nxrg) )
[1]291                ENDIF
[1353]292                lwp_av = 0.0_wp
[1]293
[2292]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
[1053]300             CASE ( 'nr' )
301                IF ( .NOT. ALLOCATED( nr_av ) )  THEN
302                   ALLOCATE( nr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
303                ENDIF
[1353]304                nr_av = 0.0_wp
[1053]305
[1691]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
[1]312             CASE ( 'p' )
313                IF ( .NOT. ALLOCATED( p_av ) )  THEN
[667]314                   ALLOCATE( p_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1]315                ENDIF
[1353]316                p_av = 0.0_wp
[1]317
318             CASE ( 'pc' )
319                IF ( .NOT. ALLOCATED( pc_av ) )  THEN
[667]320                   ALLOCATE( pc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1]321                ENDIF
[1353]322                pc_av = 0.0_wp
[1]323
324             CASE ( 'pr' )
325                IF ( .NOT. ALLOCATED( pr_av ) )  THEN
[667]326                   ALLOCATE( pr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1]327                ENDIF
[1353]328                pr_av = 0.0_wp
[1]329
[1053]330             CASE ( 'prr' )
331                IF ( .NOT. ALLOCATED( prr_av ) )  THEN
332                   ALLOCATE( prr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
333                ENDIF
[1353]334                prr_av = 0.0_wp
[1053]335
[72]336             CASE ( 'prr*' )
337                IF ( .NOT. ALLOCATED( precipitation_rate_av ) )  THEN
[667]338                   ALLOCATE( precipitation_rate_av(nysg:nyng,nxlg:nxrg) )
[72]339                ENDIF
[1353]340                precipitation_rate_av = 0.0_wp
[72]341
[1]342             CASE ( 'pt' )
343                IF ( .NOT. ALLOCATED( pt_av ) )  THEN
[667]344                   ALLOCATE( pt_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1]345                ENDIF
[1353]346                pt_av = 0.0_wp
[1]347
348             CASE ( 'q' )
349                IF ( .NOT. ALLOCATED( q_av ) )  THEN
[667]350                   ALLOCATE( q_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1]351                ENDIF
[1353]352                q_av = 0.0_wp
[1]353
[1115]354             CASE ( 'qc' )
355                IF ( .NOT. ALLOCATED( qc_av ) )  THEN
356                   ALLOCATE( qc_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
357                ENDIF
[1353]358                qc_av = 0.0_wp
[1115]359
[1]360             CASE ( 'ql' )
361                IF ( .NOT. ALLOCATED( ql_av ) )  THEN
[667]362                   ALLOCATE( ql_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1]363                ENDIF
[1353]364                ql_av = 0.0_wp
[1]365
366             CASE ( 'ql_c' )
367                IF ( .NOT. ALLOCATED( ql_c_av ) )  THEN
[667]368                   ALLOCATE( ql_c_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1]369                ENDIF
[1353]370                ql_c_av = 0.0_wp
[1]371
372             CASE ( 'ql_v' )
373                IF ( .NOT. ALLOCATED( ql_v_av ) )  THEN
[667]374                   ALLOCATE( ql_v_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1]375                ENDIF
[1353]376                ql_v_av = 0.0_wp
[1]377
378             CASE ( 'ql_vp' )
379                IF ( .NOT. ALLOCATED( ql_vp_av ) )  THEN
[667]380                   ALLOCATE( ql_vp_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1]381                ENDIF
[1353]382                ql_vp_av = 0.0_wp
[1]383
[1053]384             CASE ( 'qr' )
385                IF ( .NOT. ALLOCATED( qr_av ) )  THEN
386                   ALLOCATE( qr_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
387                ENDIF
[1353]388                qr_av = 0.0_wp
[1053]389
[354]390             CASE ( 'qsws*' )
391                IF ( .NOT. ALLOCATED( qsws_av ) )  THEN
[667]392                   ALLOCATE( qsws_av(nysg:nyng,nxlg:nxrg) )
[354]393                ENDIF
[1353]394                qsws_av = 0.0_wp
[354]395
[1]396             CASE ( 'qv' )
397                IF ( .NOT. ALLOCATED( qv_av ) )  THEN
[667]398                   ALLOCATE( qv_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1]399                ENDIF
[1353]400                qv_av = 0.0_wp
[1]401
[2735]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
[2031]408             CASE ( 'rho_ocean' )
409                IF ( .NOT. ALLOCATED( rho_ocean_av ) )  THEN
410                   ALLOCATE( rho_ocean_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[96]411                ENDIF
[2031]412                rho_ocean_av = 0.0_wp
[96]413
[1]414             CASE ( 's' )
415                IF ( .NOT. ALLOCATED( s_av ) )  THEN
[667]416                   ALLOCATE( s_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1]417                ENDIF
[1353]418                s_av = 0.0_wp
[1]419
[96]420             CASE ( 'sa' )
421                IF ( .NOT. ALLOCATED( sa_av ) )  THEN
[667]422                   ALLOCATE( sa_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[96]423                ENDIF
[1353]424                sa_av = 0.0_wp
[96]425
[354]426             CASE ( 'shf*' )
427                IF ( .NOT. ALLOCATED( shf_av ) )  THEN
[667]428                   ALLOCATE( shf_av(nysg:nyng,nxlg:nxrg) )
[354]429                ENDIF
[1353]430                shf_av = 0.0_wp
[2024]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               
[354]437
[1]438             CASE ( 't*' )
439                IF ( .NOT. ALLOCATED( ts_av ) )  THEN
[667]440                   ALLOCATE( ts_av(nysg:nyng,nxlg:nxrg) )
[1]441                ENDIF
[1353]442                ts_av = 0.0_wp
[1]443
[2742]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
[1]450             CASE ( 'u' )
451                IF ( .NOT. ALLOCATED( u_av ) )  THEN
[667]452                   ALLOCATE( u_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1]453                ENDIF
[1353]454                u_av = 0.0_wp
[1]455
456             CASE ( 'u*' )
457                IF ( .NOT. ALLOCATED( us_av ) )  THEN
[667]458                   ALLOCATE( us_av(nysg:nyng,nxlg:nxrg) )
[1]459                ENDIF
[1353]460                us_av = 0.0_wp
[1]461
462             CASE ( 'v' )
463                IF ( .NOT. ALLOCATED( v_av ) )  THEN
[667]464                   ALLOCATE( v_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1]465                ENDIF
[1353]466                v_av = 0.0_wp
[1]467
468             CASE ( 'vpt' )
469                IF ( .NOT. ALLOCATED( vpt_av ) )  THEN
[667]470                   ALLOCATE( vpt_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1]471                ENDIF
[1353]472                vpt_av = 0.0_wp
[1]473
474             CASE ( 'w' )
475                IF ( .NOT. ALLOCATED( w_av ) )  THEN
[667]476                   ALLOCATE( w_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
[1]477                ENDIF
[1353]478                w_av = 0.0_wp
[1]479
[72]480             CASE ( 'z0*' )
481                IF ( .NOT. ALLOCATED( z0_av ) )  THEN
[667]482                   ALLOCATE( z0_av(nysg:nyng,nxlg:nxrg) )
[72]483                ENDIF
[1353]484                z0_av = 0.0_wp
[72]485
[978]486             CASE ( 'z0h*' )
487                IF ( .NOT. ALLOCATED( z0h_av ) )  THEN
488                   ALLOCATE( z0h_av(nysg:nyng,nxlg:nxrg) )
489                ENDIF
[1353]490                z0h_av = 0.0_wp
[978]491
[1788]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
[2007]497!             
498!--          Block of urban surface model outputs
499             CASE ( 'usm_output' )
[1788]500
[2007]501                CALL usm_average_3d_data( 'allocate', doav(ii) )
502             
503
[1]504             CASE DEFAULT
[1972]505
[1]506!
[2696]507!--             Turbulence closure module
508                CALL tcm_3d_data_averaging( 'allocate', doav(ii) )
509
510!
[1972]511!--             Land surface quantity
512                IF ( land_surface )  THEN
513                   CALL lsm_3d_data_averaging( 'allocate', doav(ii) )
514                ENDIF
515
516!
[1976]517!--             Radiation quantity
518                IF ( radiation )  THEN
519                   CALL radiation_3d_data_averaging( 'allocate', doav(ii) )
520                ENDIF
521
522!
[2817]523!--             Gust module quantities
524                IF ( gust_module_enabled )  THEN
525                   CALL gust_3d_data_averaging( 'allocate', doav(ii) )
526                ENDIF
527
528!
[2696]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!
[1]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!
[2007]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) )
[2011]559          IF ( urban_surface  .AND.  trimvar(1:4) == 'usm_' )  THEN
[2007]560             trimvar = 'usm_output'
561          ENDIF
562!
[1]563!--    Store the array chosen on the temporary array.
[2007]564       SELECT CASE ( trimvar )
[1]565
[2797]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
[1]584          CASE ( 'e' )
[667]585             DO  i = nxlg, nxrg
586                DO  j = nysg, nyng
[1]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
[771]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
[1]602          CASE ( 'lwp*' )
[667]603             DO  i = nxlg, nxrg
604                DO  j = nysg, nyng
[1691]605                   lwp_av(j,i) = lwp_av(j,i) + SUM( ql(nzb:nzt,j,i)            &
606                                               * dzw(1:nzt+1) ) * rho_surface
[1]607                ENDDO
608             ENDDO
609
[2292]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
[1053]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
[1691]628          CASE ( 'ol*' )
[2232]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)
[1691]633             ENDDO
[2232]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
[1691]644
[1]645          CASE ( 'p' )
[667]646             DO  i = nxlg, nxrg
647                DO  j = nysg, nyng
[1]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
[1359]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
[1353]671                      s_r3 = 0.0_wp
[1359]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
[1]680                      ENDDO
[1359]681
682                      IF ( s_r2 > 0.0_wp )  THEN
683                         mean_r = s_r3 / s_r2
[1]684                      ELSE
[1353]685                         mean_r = 0.0_wp
[1]686                      ENDIF
687                      pr_av(k,j,i) = pr_av(k,j,i) + mean_r
688                   ENDDO
689                ENDDO
690             ENDDO
691
[1359]692
[72]693          CASE ( 'pr*' )
[667]694             DO  i = nxlg, nxrg
695                DO  j = nysg, nyng
[72]696                   precipitation_rate_av(j,i) = precipitation_rate_av(j,i) + &
697                                                precipitation_rate(j,i)
698                ENDDO
699             ENDDO
700
[1]701          CASE ( 'pt' )
702             IF ( .NOT. cloud_physics ) THEN
[667]703             DO  i = nxlg, nxrg
704                DO  j = nysg, nyng
705                   DO  k = nzb, nzt+1
[1]706                         pt_av(k,j,i) = pt_av(k,j,i) + pt(k,j,i)
707                      ENDDO
708                   ENDDO
709                ENDDO
710             ELSE
[667]711             DO  i = nxlg, nxrg
712                DO  j = nysg, nyng
713                   DO  k = nzb, nzt+1
[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' )
[667]722             DO  i = nxlg, nxrg
723                DO  j = nysg, nyng
[1]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
[402]729
[1115]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
[1]739          CASE ( 'ql' )
[667]740             DO  i = nxlg, nxrg
741                DO  j = nysg, nyng
[1]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' )
[667]749             DO  i = nxlg, nxrg
750                DO  j = nysg, nyng
[1]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' )
[667]758             DO  i = nxlg, nxrg
759                DO  j = nysg, nyng
[1]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' )
[1007]767             DO  i = nxl, nxr
768                DO  j = nys, nyn
[1]769                   DO  k = nzb, nzt+1
[1359]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
[1007]779                      ENDDO
[1]780                   ENDDO
781                ENDDO
782             ENDDO
783
[1053]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
[402]793          CASE ( 'qsws*' )
[2743]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.
[2232]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)
[2790]801                k = surf_def_h(0)%k(m)
[2743]802                qsws_av(j,i) = qsws_av(j,i) + surf_def_h(0)%qsws(m) *          &
803                                              waterflux_output_conversion(k)
[402]804             ENDDO
[2232]805             DO  m = 1, surf_lsm_h%ns
806                i = surf_lsm_h%i(m)
807                j = surf_lsm_h%j(m)
[2743]808                qsws_av(j,i) = qsws_av(j,i) + surf_lsm_h%qsws(m) * l_v
[2232]809             ENDDO
810             DO  m = 1, surf_usm_h%ns
811                i = surf_usm_h%i(m)
812                j = surf_usm_h%j(m)
[2743]813                qsws_av(j,i) = qsws_av(j,i) + surf_usm_h%qsws(m) * l_v
[2232]814             ENDDO
[402]815
[1]816          CASE ( 'qv' )
[667]817             DO  i = nxlg, nxrg
818                DO  j = nysg, nyng
[1]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
[2735]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
[2031]845          CASE ( 'rho_ocean' )
[667]846             DO  i = nxlg, nxrg
847                DO  j = nysg, nyng
[96]848                   DO  k = nzb, nzt+1
[2031]849                      rho_ocean_av(k,j,i) = rho_ocean_av(k,j,i) + rho_ocean(k,j,i)
[96]850                   ENDDO
851                ENDDO
852             ENDDO
[402]853
[1]854          CASE ( 's' )
[667]855             DO  i = nxlg, nxrg
856                DO  j = nysg, nyng
[1]857                   DO  k = nzb, nzt+1
[1992]858                      s_av(k,j,i) = s_av(k,j,i) + s(k,j,i)
[1]859                   ENDDO
860                ENDDO
861             ENDDO
[402]862
[96]863          CASE ( 'sa' )
[667]864             DO  i = nxlg, nxrg
865                DO  j = nysg, nyng
[96]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
[402]871
872          CASE ( 'shf*' )
[2743]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.
[2232]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)
[2790]880                k = surf_def_h(0)%k(m)
[2743]881                shf_av(j,i) = shf_av(j,i) + surf_def_h(0)%shf(m)  *            &
882                                            heatflux_output_conversion(k)
[402]883             ENDDO
[2232]884             DO  m = 1, surf_lsm_h%ns
885                i = surf_lsm_h%i(m)
886                j = surf_lsm_h%j(m)
[2743]887                shf_av(j,i) = shf_av(j,i) + surf_lsm_h%shf(m) * cp
[2232]888             ENDDO
889             DO  m = 1, surf_usm_h%ns
890                i = surf_usm_h%i(m)
891                j = surf_usm_h%j(m)
[2743]892                shf_av(j,i) = shf_av(j,i) + surf_usm_h%shf(m) * cp
[2232]893             ENDDO
[402]894
[1960]895          CASE ( 'ssws*' )
[2232]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)
[1960]900             ENDDO
[2232]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
[1960]911
[1]912          CASE ( 't*' )
[2232]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)
[1]917             ENDDO
[2232]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
[1]928
[2742]929          CASE ( 'tsurf*' )
[2798]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
[2742]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
[1]948          CASE ( 'u' )
[667]949             DO  i = nxlg, nxrg
950                DO  j = nysg, nyng
[1]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*' )
[2232]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)
[1]962             ENDDO
[2232]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
[1]973
974          CASE ( 'v' )
[667]975             DO  i = nxlg, nxrg
976                DO  j = nysg, nyng
[1]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' )
[667]984             DO  i = nxlg, nxrg
985                DO  j = nysg, nyng
[1]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' )
[667]993             DO  i = nxlg, nxrg
994                DO  j = nysg, nyng
[1]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
[72]1001          CASE ( 'z0*' )
[2232]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)
[72]1006             ENDDO
[2232]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
[72]1017
[978]1018          CASE ( 'z0h*' )
[2232]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)
[978]1023             ENDDO
[2232]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
[978]1034
[1788]1035          CASE ( 'z0q*' )
[2232]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)
[1788]1040             ENDDO
[2232]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
[2007]1051!             
[2696]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.
[2007]1056          CASE ( 'usm_output' )
[2696]1057             CALL usm_average_3d_data( 'allocate', doav(ii) )
[2007]1058             CALL usm_average_3d_data( 'sum', doav(ii) )
[1788]1059
[1]1060          CASE DEFAULT
1061!
[2696]1062!--          Turbulence closure module
1063             CALL tcm_3d_data_averaging( 'sum', doav(ii) )
1064
1065!
[1972]1066!--          Land surface quantity
1067             IF ( land_surface )  THEN
1068                CALL lsm_3d_data_averaging( 'sum', doav(ii) )
1069             ENDIF
1070
1071!
[1976]1072!--          Radiation quantity
1073             IF ( radiation )  THEN
1074                CALL radiation_3d_data_averaging( 'sum', doav(ii) )
1075             ENDIF
1076
1077!
[2817]1078!--          Gust module quantities
1079             IF ( gust_module_enabled )  THEN
1080                CALL gust_3d_data_averaging( 'sum', doav(ii) )
1081             ENDIF
1082
1083!
[2696]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!
[1]1096!--          User-defined quantity
1097             CALL user_3d_data_averaging( 'sum', doav(ii) )
1098
1099       END SELECT
1100
1101    ENDDO
1102
[1318]1103    CALL cpu_log( log_point(34), 'sum_up_3d_data', 'stop' )
[1]1104
1105
1106 END SUBROUTINE sum_up_3d_data
Note: See TracBrowser for help on using the repository browser.