source: palm/trunk/SOURCE/flow_statistics.f90 @ 2320

Last change on this file since 2320 was 2320, checked in by suehring, 4 years ago

large-scale forcing and nudging modularized

  • Property svn:keywords set to Id
File size: 95.0 KB
RevLine 
[1682]1!> @file flow_statistics.f90
[2000]2!------------------------------------------------------------------------------!
[1036]3! This file is part of PALM.
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!
[2101]17! Copyright 1997-2017 Leibniz Universitaet Hannover
[2000]18!------------------------------------------------------------------------------!
[1036]19!
[254]20! Current revisions:
[1]21! -----------------
[1961]22!
[2233]23!
[1739]24! Former revisions:
25! -----------------
26! $Id: flow_statistics.f90 2320 2017-07-21 12:47:43Z suehring $
[2320]27! Modularize large-scale forcing and nudging
28!
29! 2296 2017-06-28 07:53:56Z maronga
[2296]30! Enabled output of radiation quantities for radiation_scheme = 'constant'
31!
32! 2292 2017-06-20 09:51:42Z schwenkel
[2292]33! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
34! includes two more prognostic equations for cloud drop concentration (nc) 
35! and cloud water content (qc).
36!
37! 2270 2017-06-09 12:18:47Z maronga
[2270]38! Revised numbering (removed 2 timeseries)
39!
40! 2252 2017-06-07 09:35:37Z knoop
[2252]41! perturbation pressure now depending on flux_output_mode
42!
43! 2233 2017-05-30 18:08:54Z suehring
[1739]44!
[2233]45! 2232 2017-05-30 17:47:52Z suehring
46! Adjustments to new topography and surface concept
47!
[2119]48! 2118 2017-01-17 16:38:49Z raasch
49! OpenACC version of subroutine removed
50!
[2074]51! 2073 2016-11-30 14:34:05Z raasch
52! openmp bugfix: large scale forcing calculations cannot be executed thread
53! parallel
54!
[2038]55! 2037 2016-10-26 11:15:40Z knoop
56! Anelastic approximation implemented
57!
[2032]58! 2031 2016-10-21 15:11:58Z knoop
59! renamed variable rho to rho_ocean
60!
[2027]61! 2026 2016-10-18 10:27:02Z suehring
62! Bugfix, enable output of s*2.
63! Change, calculation of domain-averaged perturbation energy.
64! Some formatting adjustments.
65!
[2001]66! 2000 2016-08-20 18:09:15Z knoop
67! Forced header and separation lines into 80 columns
68!
[1977]69! 1976 2016-07-27 13:28:04Z maronga
70! Removed some unneeded __rrtmg preprocessor directives
71!
[1961]72! 1960 2016-07-12 16:34:24Z suehring
73! Separate humidity and passive scalar
74!
[1919]75! 1918 2016-05-27 14:35:57Z raasch
76! in case of Wicker-Skamarock scheme, calculate disturbance kinetic energy here,
77! if flow_statistics is called before the first initial time step
78!
[1854]79! 1853 2016-04-11 09:00:35Z maronga
80! Adjusted for use with radiation_scheme = constant
81!
[1851]82! 1849 2016-04-08 11:33:18Z hoffmann
83! prr moved to arrays_3d
84!
[1823]85! 1822 2016-04-07 07:49:42Z hoffmann
86! Output of bulk microphysics simplified.
87!
[1816]88! 1815 2016-04-06 13:49:59Z raasch
89! cpp-directives for intel openmp bug removed
90!
[1784]91! 1783 2016-03-06 18:36:17Z raasch
92! +module netcdf_interface
93!
[1748]94! 1747 2016-02-08 12:25:53Z raasch
95! small bugfixes for accelerator version
96!
[1739]97! 1738 2015-12-18 13:56:05Z raasch
[1738]98! bugfixes for calculations in statistical regions which do not contain grid
99! points in the lowest vertical levels, mean surface level height considered
100! in the calculation of the characteristic vertical velocity,
101! old upstream parts removed
[1383]102!
[1710]103! 1709 2015-11-04 14:47:01Z maronga
104! Updated output of Obukhov length
105!
[1692]106! 1691 2015-10-26 16:17:44Z maronga
107! Revised calculation of Obukhov length. Added output of radiative heating >
108! rates for RRTMG.
109!
[1683]110! 1682 2015-10-07 23:56:08Z knoop
111! Code annotations made doxygen readable
112!
[1659]113! 1658 2015-09-18 10:52:53Z raasch
114! bugfix: temporary reduction variables in the openacc branch are now
115! initialized to zero
116!
[1655]117! 1654 2015-09-17 09:20:17Z raasch
118! FORTRAN bugfix of r1652
119!
[1653]120! 1652 2015-09-17 08:12:24Z raasch
121! bugfix in calculation of energy production by turbulent transport of TKE
122!
[1594]123! 1593 2015-05-16 13:58:02Z raasch
124! FORTRAN errors removed from openacc branch
125!
[1586]126! 1585 2015-04-30 07:05:52Z maronga
127! Added output of timeseries and profiles for RRTMG
128!
[1572]129! 1571 2015-03-12 16:12:49Z maronga
130! Bugfix: output of rad_net and rad_sw_in
131!
[1568]132! 1567 2015-03-10 17:57:55Z suehring
133! Reverse modifications made for monotonic limiter.
134!
[1558]135! 1557 2015-03-05 16:43:04Z suehring
136! Adjustments for monotonic limiter
137!
[1556]138! 1555 2015-03-04 17:44:27Z maronga
139! Added output of r_a and r_s.
140!
[1552]141! 1551 2015-03-03 14:18:16Z maronga
142! Added suppport for land surface model and radiation model output.
143!
[1499]144! 1498 2014-12-03 14:09:51Z suehring
145! Comments added
146!
[1483]147! 1482 2014-10-18 12:34:45Z raasch
148! missing ngp_sums_ls added in accelerator version
149!
[1451]150! 1450 2014-08-21 07:31:51Z heinze
151! bugfix: calculate fac only for simulated_time >= 0.0
152!
[1397]153! 1396 2014-05-06 13:37:41Z raasch
154! bugfix: "copyin" replaced by "update device" in openacc-branch
155!
[1387]156! 1386 2014-05-05 13:55:30Z boeske
157! bugfix: simulated time before the last timestep is needed for the correct
158! calculation of the profiles of large scale forcing tendencies
159!
[1383]160! 1382 2014-04-30 12:15:41Z boeske
[1382]161! Renamed variables which store large scale forcing tendencies
162! pt_lsa -> td_lsa_lpt, pt_subs -> td_sub_lpt,
163! q_lsa  -> td_lsa_q,   q_subs  -> td_sub_q,
164! added Neumann boundary conditions for profile data output of large scale
165! advection and subsidence terms at nzt+1
[1354]166!
[1375]167! 1374 2014-04-25 12:55:07Z raasch
168! bugfix: syntax errors removed from openacc-branch
169! missing variables added to ONLY-lists
170!
[1366]171! 1365 2014-04-22 15:03:56Z boeske
172! Output of large scale advection, large scale subsidence and nudging tendencies
173! +sums_ls_l, ngp_sums_ls, use_subsidence_tendencies
174!
[1354]175! 1353 2014-04-08 15:21:23Z heinze
176! REAL constants provided with KIND-attribute
177!
[1323]178! 1322 2014-03-20 16:38:49Z raasch
179! REAL constants defined as wp-kind
180!
[1321]181! 1320 2014-03-20 08:40:49Z raasch
[1320]182! ONLY-attribute added to USE-statements,
183! kind-parameters added to all INTEGER and REAL declaration statements,
184! kinds are defined in new module kinds,
185! revision history before 2012 removed,
186! comment fields (!:) to be used for variable explanations added to
187! all variable declaration statements
[1008]188!
[1300]189! 1299 2014-03-06 13:15:21Z heinze
190! Output of large scale vertical velocity w_subs
191!
[1258]192! 1257 2013-11-08 15:18:40Z raasch
193! openacc "end parallel" replaced by "end parallel loop"
194!
[1242]195! 1241 2013-10-30 11:36:58Z heinze
196! Output of ug and vg
197!
[1222]198! 1221 2013-09-10 08:59:13Z raasch
199! ported for openACC in separate #else branch
200!
[1182]201! 1179 2013-06-14 05:57:58Z raasch
202! comment for profile 77 added
203!
[1116]204! 1115 2013-03-26 18:16:16Z hoffmann
205! ql is calculated by calc_liquid_water_content
206!
[1112]207! 1111 2013-03-08 23:54:10Z raasch
208! openACC directive added
209!
[1054]210! 1053 2012-11-13 17:11:03Z hoffmann
[1112]211! additions for two-moment cloud physics scheme:
[1054]212! +nr, qr, qc, prr
213!
[1037]214! 1036 2012-10-22 13:43:42Z raasch
215! code put under GPL (PALM 3.9)
216!
[1008]217! 1007 2012-09-19 14:30:36Z franke
[1007]218! Calculation of buoyancy flux for humidity in case of WS-scheme is now using
219! turbulent fluxes of WS-scheme
220! Bugfix: Calculation of subgridscale buoyancy flux for humidity and cloud
221! droplets at nzb and nzt added
[700]222!
[802]223! 801 2012-01-10 17:30:36Z suehring
224! Calculation of turbulent fluxes in advec_ws is now thread-safe.
225!
[1]226! Revision 1.1  1997/08/11 06:15:17  raasch
227! Initial revision
228!
229!
230! Description:
231! ------------
[1682]232!> Compute average profiles and further average flow quantities for the different
233!> user-defined (sub-)regions. The region indexed 0 is the total model domain.
234!>
235!> @note For simplicity, nzb_s_inner and nzb_diff_s_inner are being used as a
236!>       lower vertical index for k-loops for all variables, although strictly
237!>       speaking the k-loops would have to be split up according to the staggered
238!>       grid. However, this implies no error since staggered velocity components
239!>       are zero at the walls and inside buildings.
[1]240!------------------------------------------------------------------------------!
[1682]241 SUBROUTINE flow_statistics
242 
[1]243
[1320]244    USE arrays_3d,                                                             &
[2037]245        ONLY:  ddzu, ddzw, e, heatflux_output_conversion, hyp, km, kh,         &
[2292]246               momentumflux_output_conversion, nc, nr, p, prho, prr, pt, q,    &
[2232]247               qc, ql, qr, rho_air, rho_air_zw, rho_ocean, s,                  &
[2320]248               sa, u, ug, v, vg, vpt, w, w_subs, waterflux_output_conversion, zw
[1320]249       
250    USE cloud_parameters,                                                      &
[1849]251        ONLY:   l_d_cp, pt_d_t
[1320]252       
253    USE control_parameters,                                                    &
[1551]254        ONLY:   average_count_pr, cloud_droplets, cloud_physics, do_sum,       &
[2232]255                dt_3d, g, humidity, kappa, land_surface, large_scale_forcing,  &
[1691]256                large_scale_subsidence, max_pr_user, message_string, neutral,  &
[2292]257                microphysics_morrison, microphysics_seifert, ocean,            &
258                passive_scalar, simulated_time, use_subsidence_tendencies,     &
259                use_surface_fluxes, use_top_fluxes, ws_scheme_mom,             &
260                ws_scheme_sca
[1320]261       
262    USE cpulog,                                                                &
[1551]263        ONLY:   cpu_log, log_point
[1320]264       
265    USE grid_variables,                                                        &
[1551]266        ONLY:   ddx, ddy
[1320]267       
268    USE indices,                                                               &
[1551]269        ONLY:   ngp_2dh, ngp_2dh_s_inner, ngp_3d, ngp_3d_inner, ngp_sums,      &
[2232]270                ngp_sums_ls, nxl, nxr, nyn, nys, nzb, nzt, wall_flags_0
[1320]271       
272    USE kinds
273   
[1551]274    USE land_surface_model_mod,                                                &
[2232]275        ONLY:   m_soil_h, nzb_soil, nzt_soil, t_soil_h
[1551]276
[2320]277    USE lsf_nudging_mod,                                                       &
278        ONLY:   td_lsa_lpt, td_lsa_q, td_sub_lpt, td_sub_q, time_vert
279
[1783]280    USE netcdf_interface,                                                      &
281        ONLY:  dots_rad, dots_soil
282
[1]283    USE pegrid
[1551]284
285    USE radiation_model_mod,                                                   &
[1783]286        ONLY:  radiation, radiation_scheme, rad_net,                           &
[1691]287               rad_lw_in, rad_lw_out, rad_lw_cs_hr, rad_lw_hr,                 &
288               rad_sw_in, rad_sw_out, rad_sw_cs_hr, rad_sw_hr
[1585]289
290#if defined ( __rrtmg )
291    USE radiation_model_mod,                                                   &
292        ONLY:  rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir
293#endif
294 
[1]295    USE statistics
296
[2232]297       USE surface_mod,                                                        &
298          ONLY :  surf_def_h, surf_lsm_h, surf_usm_h
[1691]299
[2232]300
[1]301    IMPLICIT NONE
302
[1682]303    INTEGER(iwp) ::  i                   !<
304    INTEGER(iwp) ::  j                   !<
305    INTEGER(iwp) ::  k                   !<
[2232]306    INTEGER(iwp) ::  ki                  !<
[1738]307    INTEGER(iwp) ::  k_surface_level     !<
[2232]308    INTEGER(iwp) ::  m                   !< loop variable over all horizontal wall elements
309    INTEGER(iwp) ::  l                   !< loop variable over surface facing -- up- or downward-facing
[1682]310    INTEGER(iwp) ::  nt                  !<
311    INTEGER(iwp) ::  omp_get_thread_num  !<
312    INTEGER(iwp) ::  sr                  !<
313    INTEGER(iwp) ::  tn                  !<
[2232]314
[1682]315    LOGICAL ::  first  !<
[1320]316   
[1682]317    REAL(wp) ::  dptdz_threshold  !<
318    REAL(wp) ::  fac              !<
[2232]319    REAL(wp) ::  flag             !<
[1682]320    REAL(wp) ::  height           !<
321    REAL(wp) ::  pts              !<
322    REAL(wp) ::  sums_l_eper      !<
323    REAL(wp) ::  sums_l_etot      !<
324    REAL(wp) ::  ust              !<
325    REAL(wp) ::  ust2             !<
326    REAL(wp) ::  u2               !<
327    REAL(wp) ::  vst              !<
328    REAL(wp) ::  vst2             !<
329    REAL(wp) ::  v2               !<
330    REAL(wp) ::  w2               !<
331    REAL(wp) ::  z_i(2)           !<
[1320]332   
[1682]333    REAL(wp) ::  dptdz(nzb+1:nzt+1)    !<
334    REAL(wp) ::  sums_ll(nzb:nzt+1,2)  !<
[1]335
336    CALL cpu_log( log_point(10), 'flow_statistics', 'start' )
337
[1221]338
[1]339!
340!-- To be on the safe side, check whether flow_statistics has already been
341!-- called once after the current time step
342    IF ( flow_statistics_called )  THEN
[254]343
[274]344       message_string = 'flow_statistics is called two times within one ' // &
345                        'timestep'
[254]346       CALL message( 'flow_statistics', 'PA0190', 1, 2, 0, 6, 0 )
[1007]347
[1]348    ENDIF
349
350!
351!-- Compute statistics for each (sub-)region
352    DO  sr = 0, statistic_regions
353
354!
355!--    Initialize (local) summation array
[1353]356       sums_l = 0.0_wp
[1]357
358!
359!--    Store sums that have been computed in other subroutines in summation
360!--    array
361       sums_l(:,11,:) = sums_l_l(:,sr,:)      ! mixing length from diffusivities
362!--    WARNING: next line still has to be adjusted for OpenMP
[2037]363       sums_l(:,21,0) = sums_wsts_bc_l(:,sr) *                                 &
364                        heatflux_output_conversion  ! heat flux from advec_s_bc
[87]365       sums_l(nzb+9,pr_palm,0)  = sums_divold_l(sr)  ! old divergence from pres
366       sums_l(nzb+10,pr_palm,0) = sums_divnew_l(sr)  ! new divergence from pres
[1]367
[667]368!
[1498]369!--    When calcuating horizontally-averaged total (resolved- plus subgrid-
370!--    scale) vertical fluxes and velocity variances by using commonly-
371!--    applied Reynolds-based methods ( e.g. <w'pt'> = (w-<w>)*(pt-<pt>) )
372!--    in combination with the 5th order advection scheme, pronounced
373!--    artificial kinks could be observed in the vertical profiles near the
374!--    surface. Please note: these kinks were not related to the model truth,
375!--    i.e. these kinks are just related to an evaluation problem.   
376!--    In order avoid these kinks, vertical fluxes and horizontal as well
377!--    vertical velocity variances are calculated directly within the advection
378!--    routines, according to the numerical discretization, to evaluate the
379!--    statistical quantities as they will appear within the prognostic
380!--    equations.
[667]381!--    Copy the turbulent quantities, evaluated in the advection routines to
[1498]382!--    the local array sums_l() for further computations.
[743]383       IF ( ws_scheme_mom .AND. sr == 0 )  THEN
[696]384
[1007]385!
[673]386!--       According to the Neumann bc for the horizontal velocity components,
387!--       the corresponding fluxes has to satisfiy the same bc.
388          IF ( ocean )  THEN
[801]389             sums_us2_ws_l(nzt+1,:) = sums_us2_ws_l(nzt,:)
[1007]390             sums_vs2_ws_l(nzt+1,:) = sums_vs2_ws_l(nzt,:)
[673]391          ENDIF
[696]392
393          DO  i = 0, threads_per_task-1
[1007]394!
[696]395!--          Swap the turbulent quantities evaluated in advec_ws.
[2037]396             sums_l(:,13,i) = sums_wsus_ws_l(:,i)                              &
397                              * momentumflux_output_conversion ! w*u*
398             sums_l(:,15,i) = sums_wsvs_ws_l(:,i)                              &
399                              * momentumflux_output_conversion ! w*v*
[801]400             sums_l(:,30,i) = sums_us2_ws_l(:,i)        ! u*2
401             sums_l(:,31,i) = sums_vs2_ws_l(:,i)        ! v*2
402             sums_l(:,32,i) = sums_ws2_ws_l(:,i)        ! w*2
[1353]403             sums_l(:,34,i) = sums_l(:,34,i) + 0.5_wp *                        & 
[1320]404                              ( sums_us2_ws_l(:,i) + sums_vs2_ws_l(:,i) +      &
[801]405                                sums_ws2_ws_l(:,i) )    ! e*
[667]406          ENDDO
[696]407
[667]408       ENDIF
[696]409
[1567]410       IF ( ws_scheme_sca .AND. sr == 0 )  THEN
[696]411
412          DO  i = 0, threads_per_task-1
[2037]413             sums_l(:,17,i)                        = sums_wspts_ws_l(:,i)      &
414                                           * heatflux_output_conversion  ! w*pt*
[1960]415             IF ( ocean          ) sums_l(:,66,i)  = sums_wssas_ws_l(:,i) ! w*sa*
[2037]416             IF ( humidity       ) sums_l(:,49,i)  = sums_wsqs_ws_l(:,i)       &
417                                           * waterflux_output_conversion  ! w*q*
[2270]418             IF ( passive_scalar ) sums_l(:,114,i) = sums_wsss_ws_l(:,i)  ! w*s*
[696]419          ENDDO
420
[667]421       ENDIF
[305]422!
[1]423!--    Horizontally averaged profiles of horizontal velocities and temperature.
424!--    They must have been computed before, because they are already required
425!--    for other horizontal averages.
426       tn = 0
[2232]427       !$OMP PARALLEL PRIVATE( i, j, k, tn, flag )
428       !$ tn = omp_get_thread_num()
[1]429       !$OMP DO
430       DO  i = nxl, nxr
431          DO  j =  nys, nyn
[2232]432             DO  k = nzb, nzt+1
433                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) )
434                sums_l(k,1,tn)  = sums_l(k,1,tn)  + u(k,j,i)  * rmask(j,i,sr)  &
435                                                              * flag
436                sums_l(k,2,tn)  = sums_l(k,2,tn)  + v(k,j,i)  * rmask(j,i,sr)  &
437                                                              * flag
438                sums_l(k,4,tn)  = sums_l(k,4,tn)  + pt(k,j,i) * rmask(j,i,sr)  &
439                                                              * flag
[1]440             ENDDO
441          ENDDO
442       ENDDO
443
444!
[96]445!--    Horizontally averaged profile of salinity
446       IF ( ocean )  THEN
447          !$OMP DO
448          DO  i = nxl, nxr
449             DO  j =  nys, nyn
[2232]450                DO  k = nzb, nzt+1
451                   sums_l(k,23,tn)  = sums_l(k,23,tn) + sa(k,j,i)              &
452                                    * rmask(j,i,sr)                            &
453                                    * MERGE( 1.0_wp, 0.0_wp,                   &
454                                             BTEST( wall_flags_0(k,j,i), 22 ) )
[96]455                ENDDO
456             ENDDO
457          ENDDO
458       ENDIF
459
460!
[1]461!--    Horizontally averaged profiles of virtual potential temperature,
462!--    total water content, specific humidity and liquid water potential
463!--    temperature
[75]464       IF ( humidity )  THEN
[1]465          !$OMP DO
466          DO  i = nxl, nxr
467             DO  j =  nys, nyn
[2232]468                DO  k = nzb, nzt+1
469                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) )
470                   sums_l(k,44,tn)  = sums_l(k,44,tn) +                        &
471                                      vpt(k,j,i) * rmask(j,i,sr) * flag
472                   sums_l(k,41,tn)  = sums_l(k,41,tn) +                        &
473                                      q(k,j,i) * rmask(j,i,sr)   * flag
[1]474                ENDDO
475             ENDDO
476          ENDDO
477          IF ( cloud_physics )  THEN
478             !$OMP DO
479             DO  i = nxl, nxr
480                DO  j =  nys, nyn
[2232]481                   DO  k = nzb, nzt+1
482                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) )
483                      sums_l(k,42,tn) = sums_l(k,42,tn) +                      &
484                                      ( q(k,j,i) - ql(k,j,i) ) * rmask(j,i,sr) &
485                                                               * flag
486                      sums_l(k,43,tn) = sums_l(k,43,tn) + (                    &
[1]487                                      pt(k,j,i) + l_d_cp*pt_d_t(k) * ql(k,j,i) &
[2232]488                                                          ) * rmask(j,i,sr)    &
489                                                            * flag
[1]490                   ENDDO
491                ENDDO
492             ENDDO
493          ENDIF
494       ENDIF
495
496!
497!--    Horizontally averaged profiles of passive scalar
498       IF ( passive_scalar )  THEN
499          !$OMP DO
500          DO  i = nxl, nxr
501             DO  j =  nys, nyn
[2232]502                DO  k = nzb, nzt+1
[2270]503                   sums_l(k,115,tn)  = sums_l(k,115,tn) + s(k,j,i)             &
[2232]504                                    * rmask(j,i,sr)                            &
505                                    * MERGE( 1.0_wp, 0.0_wp,                   &
506                                             BTEST( wall_flags_0(k,j,i), 22 ) )
[1]507                ENDDO
508             ENDDO
509          ENDDO
510       ENDIF
511       !$OMP END PARALLEL
512!
513!--    Summation of thread sums
514       IF ( threads_per_task > 1 )  THEN
515          DO  i = 1, threads_per_task-1
516             sums_l(:,1,0) = sums_l(:,1,0) + sums_l(:,1,i)
517             sums_l(:,2,0) = sums_l(:,2,0) + sums_l(:,2,i)
518             sums_l(:,4,0) = sums_l(:,4,0) + sums_l(:,4,i)
[96]519             IF ( ocean )  THEN
520                sums_l(:,23,0) = sums_l(:,23,0) + sums_l(:,23,i)
521             ENDIF
[75]522             IF ( humidity )  THEN
[1]523                sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i)
524                sums_l(:,44,0) = sums_l(:,44,0) + sums_l(:,44,i)
525                IF ( cloud_physics )  THEN
526                   sums_l(:,42,0) = sums_l(:,42,0) + sums_l(:,42,i)
527                   sums_l(:,43,0) = sums_l(:,43,0) + sums_l(:,43,i)
528                ENDIF
529             ENDIF
530             IF ( passive_scalar )  THEN
[2270]531                sums_l(:,115,0) = sums_l(:,115,0) + sums_l(:,115,i)
[1]532             ENDIF
533          ENDDO
534       ENDIF
535
536#if defined( __parallel )
537!
538!--    Compute total sum from local sums
[622]539       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
[1320]540       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, MPI_REAL,  &
[1]541                           MPI_SUM, comm2d, ierr )
[622]542       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
[1320]543       CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, MPI_REAL,  &
[1]544                           MPI_SUM, comm2d, ierr )
[622]545       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
[1320]546       CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, MPI_REAL,  &
[1]547                           MPI_SUM, comm2d, ierr )
[96]548       IF ( ocean )  THEN
[622]549          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
[1320]550          CALL MPI_ALLREDUCE( sums_l(nzb,23,0), sums(nzb,23), nzt+2-nzb,       &
[96]551                              MPI_REAL, MPI_SUM, comm2d, ierr )
552       ENDIF
[75]553       IF ( humidity ) THEN
[622]554          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
[1320]555          CALL MPI_ALLREDUCE( sums_l(nzb,44,0), sums(nzb,44), nzt+2-nzb,       &
[1]556                              MPI_REAL, MPI_SUM, comm2d, ierr )
[622]557          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
[1320]558          CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb,       &
[1]559                              MPI_REAL, MPI_SUM, comm2d, ierr )
560          IF ( cloud_physics ) THEN
[622]561             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
[1320]562             CALL MPI_ALLREDUCE( sums_l(nzb,42,0), sums(nzb,42), nzt+2-nzb,    &
[1]563                                 MPI_REAL, MPI_SUM, comm2d, ierr )
[622]564             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
[1320]565             CALL MPI_ALLREDUCE( sums_l(nzb,43,0), sums(nzb,43), nzt+2-nzb,    &
[1]566                                 MPI_REAL, MPI_SUM, comm2d, ierr )
567          ENDIF
568       ENDIF
569
570       IF ( passive_scalar )  THEN
[622]571          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
[2270]572          CALL MPI_ALLREDUCE( sums_l(nzb,115,0), sums(nzb,115), nzt+2-nzb,       &
[1]573                              MPI_REAL, MPI_SUM, comm2d, ierr )
574       ENDIF
575#else
576       sums(:,1) = sums_l(:,1,0)
577       sums(:,2) = sums_l(:,2,0)
578       sums(:,4) = sums_l(:,4,0)
[96]579       IF ( ocean )  sums(:,23) = sums_l(:,23,0)
[75]580       IF ( humidity ) THEN
[1]581          sums(:,44) = sums_l(:,44,0)
582          sums(:,41) = sums_l(:,41,0)
583          IF ( cloud_physics ) THEN
584             sums(:,42) = sums_l(:,42,0)
585             sums(:,43) = sums_l(:,43,0)
586          ENDIF
587       ENDIF
[2270]588       IF ( passive_scalar )  sums(:,115) = sums_l(:,115,0)
[1]589#endif
590
591!
592!--    Final values are obtained by division by the total number of grid points
593!--    used for summation. After that store profiles.
[132]594       sums(:,1) = sums(:,1) / ngp_2dh(sr)
595       sums(:,2) = sums(:,2) / ngp_2dh(sr)
596       sums(:,4) = sums(:,4) / ngp_2dh_s_inner(:,sr)
[1]597       hom(:,1,1,sr) = sums(:,1)             ! u
598       hom(:,1,2,sr) = sums(:,2)             ! v
599       hom(:,1,4,sr) = sums(:,4)             ! pt
600
[667]601
[1]602!
[96]603!--    Salinity
604       IF ( ocean )  THEN
[132]605          sums(:,23) = sums(:,23) / ngp_2dh_s_inner(:,sr)
[96]606          hom(:,1,23,sr) = sums(:,23)             ! sa
607       ENDIF
608
609!
[1]610!--    Humidity and cloud parameters
[75]611       IF ( humidity ) THEN
[132]612          sums(:,44) = sums(:,44) / ngp_2dh_s_inner(:,sr)
613          sums(:,41) = sums(:,41) / ngp_2dh_s_inner(:,sr)
[1]614          hom(:,1,44,sr) = sums(:,44)             ! vpt
615          hom(:,1,41,sr) = sums(:,41)             ! qv (q)
616          IF ( cloud_physics ) THEN
[132]617             sums(:,42) = sums(:,42) / ngp_2dh_s_inner(:,sr)
618             sums(:,43) = sums(:,43) / ngp_2dh_s_inner(:,sr)
[1]619             hom(:,1,42,sr) = sums(:,42)             ! qv
620             hom(:,1,43,sr) = sums(:,43)             ! pt
621          ENDIF
622       ENDIF
623
624!
625!--    Passive scalar
[2270]626       IF ( passive_scalar )  hom(:,1,115,sr) = sums(:,115) /                  &
[1960]627            ngp_2dh_s_inner(:,sr)                    ! s
[1]628
629!
630!--    Horizontally averaged profiles of the remaining prognostic variables,
631!--    variances, the total and the perturbation energy (single values in last
632!--    column of sums_l) and some diagnostic quantities.
[132]633!--    NOTE: for simplicity, nzb_s_inner is used below, although strictly
[1]634!--    ----  speaking the following k-loop would have to be split up and
635!--          rearranged according to the staggered grid.
[132]636!--          However, this implies no error since staggered velocity components
637!--          are zero at the walls and inside buildings.
[1]638       tn = 0
[1815]639       !$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper,             &
640       !$OMP                   sums_l_etot, tn, ust, ust2, u2, vst, vst2, v2,  &
[2232]641       !$OMP                   w2, flag, m, ki, l )
642       !$ tn = omp_get_thread_num()
[1]643       !$OMP DO
644       DO  i = nxl, nxr
645          DO  j =  nys, nyn
[1353]646             sums_l_etot = 0.0_wp
[2232]647             DO  k = nzb, nzt+1
648                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) )
[1]649!
650!--             Prognostic and diagnostic variables
[2232]651                sums_l(k,3,tn)  = sums_l(k,3,tn)  + w(k,j,i)  * rmask(j,i,sr)  &
652                                                              * flag
653                sums_l(k,8,tn)  = sums_l(k,8,tn)  + e(k,j,i)  * rmask(j,i,sr)  &
654                                                              * flag
655                sums_l(k,9,tn)  = sums_l(k,9,tn)  + km(k,j,i) * rmask(j,i,sr)  &
656                                                              * flag
657                sums_l(k,10,tn) = sums_l(k,10,tn) + kh(k,j,i) * rmask(j,i,sr)  &
658                                                              * flag
[2252]659                sums_l(k,40,tn) = sums_l(k,40,tn) + ( p(k,j,i)                 &
660                                         / momentumflux_output_conversion(k) ) &
661                                                              * flag
[1]662
663                sums_l(k,33,tn) = sums_l(k,33,tn) + &
[2232]664                                  ( pt(k,j,i)-hom(k,1,4,sr) )**2 * rmask(j,i,sr)&
665                                                                 * flag
[624]666
667                IF ( humidity )  THEN
668                   sums_l(k,70,tn) = sums_l(k,70,tn) + &
[2232]669                                  ( q(k,j,i)-hom(k,1,41,sr) )**2 * rmask(j,i,sr)&
670                                                                 * flag
[624]671                ENDIF
[1960]672                IF ( passive_scalar )  THEN
[2270]673                   sums_l(k,116,tn) = sums_l(k,116,tn) + &
674                                  ( s(k,j,i)-hom(k,1,115,sr) )**2 * rmask(j,i,sr)&
[2232]675                                                                  * flag
[1960]676                ENDIF
[699]677!
678!--             Higher moments
679!--             (Computation of the skewness of w further below)
[2232]680                sums_l(k,38,tn) = sums_l(k,38,tn) + w(k,j,i)**3 * rmask(j,i,sr) &
681                                                                * flag
[667]682
[1]683                sums_l_etot  = sums_l_etot + &
[2232]684                                        0.5_wp * ( u(k,j,i)**2 + v(k,j,i)**2 +  &
685                                        w(k,j,i)**2 )            * rmask(j,i,sr)&
686                                                                 * flag
[1]687             ENDDO
688!
689!--          Total and perturbation energy for the total domain (being
690!--          collected in the last column of sums_l). Summation of these
691!--          quantities is seperated from the previous loop in order to
692!--          allow vectorization of that loop.
[87]693             sums_l(nzb+4,pr_palm,tn) = sums_l(nzb+4,pr_palm,tn) + sums_l_etot
[1]694!
695!--          2D-arrays (being collected in the last column of sums_l)
[2232]696             IF ( surf_def_h(0)%ns >= 1 )  THEN
697                m = surf_def_h(0)%start_index(j,i)
698                sums_l(nzb,pr_palm,tn)   = sums_l(nzb,pr_palm,tn) +               &
699                                        surf_def_h(0)%us(m)   * rmask(j,i,sr)
700                sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) +             &
701                                        surf_def_h(0)%usws(m) * rmask(j,i,sr)
702                sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) +             &
703                                        surf_def_h(0)%vsws(m) * rmask(j,i,sr)
704                sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) +             &
705                                        surf_def_h(0)%ts(m)   * rmask(j,i,sr)
706                IF ( humidity )  THEN
707                   sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) +        &
708                                            surf_def_h(0)%qs(m)   * rmask(j,i,sr)
709                ENDIF
710                IF ( passive_scalar )  THEN
711                   sums_l(nzb+13,pr_palm,tn) = sums_l(nzb+13,pr_palm,tn) +        &
712                                            surf_def_h(0)%ss(m)   * rmask(j,i,sr)
713                ENDIF
[197]714             ENDIF
[2232]715             IF ( surf_lsm_h%ns >= 1 )  THEN
716                m = surf_lsm_h%start_index(j,i)
717                sums_l(nzb,pr_palm,tn)   = sums_l(nzb,pr_palm,tn) +               &
718                                        surf_lsm_h%us(m)   * rmask(j,i,sr)
719                sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) +             &
720                                        surf_lsm_h%usws(m) * rmask(j,i,sr)
721                sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) +             &
722                                        surf_lsm_h%vsws(m) * rmask(j,i,sr)
723                sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) +             &
724                                        surf_lsm_h%ts(m)   * rmask(j,i,sr)
725                IF ( humidity )  THEN
726                   sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) +        &
727                                            surf_lsm_h%qs(m)   * rmask(j,i,sr)
728                ENDIF
729                IF ( passive_scalar )  THEN
730                   sums_l(nzb+13,pr_palm,tn) = sums_l(nzb+13,pr_palm,tn) +        &
731                                            surf_lsm_h%ss(m)   * rmask(j,i,sr)
732                ENDIF
[1960]733             ENDIF
[2232]734             IF ( surf_usm_h%ns >= 1 )  THEN
735                m = surf_lsm_h%start_index(j,i)
736                sums_l(nzb,pr_palm,tn)   = sums_l(nzb,pr_palm,tn) +               &
737                                        surf_usm_h%us(m)   * rmask(j,i,sr)
738                sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) +             &
739                                        surf_usm_h%usws(m) * rmask(j,i,sr)
740                sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) +             &
741                                        surf_usm_h%vsws(m) * rmask(j,i,sr)
742                sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) +             &
743                                        surf_usm_h%ts(m)   * rmask(j,i,sr)
744                IF ( humidity )  THEN
745                   sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) +        &
746                                            surf_usm_h%qs(m)   * rmask(j,i,sr)
747                ENDIF
748                IF ( passive_scalar )  THEN
749                   sums_l(nzb+13,pr_palm,tn) = sums_l(nzb+13,pr_palm,tn) +        &
750                                            surf_usm_h%ss(m)   * rmask(j,i,sr)
751                ENDIF
752             ENDIF
[1]753          ENDDO
754       ENDDO
755
756!
[667]757!--    Computation of statistics when ws-scheme is not used. Else these
758!--    quantities are evaluated in the advection routines.
[1918]759       IF ( .NOT. ws_scheme_mom .OR. sr /= 0 .OR. simulated_time == 0.0_wp )   &
760       THEN
[667]761          !$OMP DO
762          DO  i = nxl, nxr
763             DO  j =  nys, nyn
[2232]764                DO  k = nzb, nzt+1
765                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) )
766
[667]767                   u2   = u(k,j,i)**2
768                   v2   = v(k,j,i)**2
769                   w2   = w(k,j,i)**2
770                   ust2 = ( u(k,j,i) - hom(k,1,1,sr) )**2
771                   vst2 = ( v(k,j,i) - hom(k,1,2,sr) )**2
772
[2232]773                   sums_l(k,30,tn) = sums_l(k,30,tn) + ust2 * rmask(j,i,sr)    &
774                                                            * flag
775                   sums_l(k,31,tn) = sums_l(k,31,tn) + vst2 * rmask(j,i,sr)    &
776                                                            * flag
777                   sums_l(k,32,tn) = sums_l(k,32,tn) + w2   * rmask(j,i,sr)    &
778                                                            * flag
[667]779!
[2026]780!--                Perturbation energy
[667]781
[1353]782                   sums_l(k,34,tn) = sums_l(k,34,tn) + 0.5_wp *                &
[2232]783                                  ( ust2 + vst2 + w2 )      * rmask(j,i,sr)    &
784                                                            * flag
[667]785                ENDDO
786             ENDDO
787          ENDDO
788       ENDIF
[2026]789!
790!--    Computaion of domain-averaged perturbation energy. Please note,
791!--    to prevent that perturbation energy is larger (even if only slightly)
792!--    than the total kinetic energy, calculation is based on deviations from
793!--    the horizontal mean, instead of spatial descretization of the advection
794!--    term.
795       !$OMP DO
796       DO  i = nxl, nxr
797          DO  j =  nys, nyn
[2232]798             DO  k = nzb, nzt+1
799                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) )
800
[2026]801                w2   = w(k,j,i)**2
802                ust2 = ( u(k,j,i) - hom(k,1,1,sr) )**2
803                vst2 = ( v(k,j,i) - hom(k,1,2,sr) )**2
804                w2   = w(k,j,i)**2
[1241]805
[2026]806                sums_l(nzb+5,pr_palm,tn) = sums_l(nzb+5,pr_palm,tn)            &
[2232]807                                 + 0.5_wp * ( ust2 + vst2 + w2 )               &
808                                 * rmask(j,i,sr)                               &
809                                 * flag
[2026]810             ENDDO
811          ENDDO
812       ENDDO
813
[667]814!
[1]815!--    Horizontally averaged profiles of the vertical fluxes
[667]816
[1]817       !$OMP DO
818       DO  i = nxl, nxr
819          DO  j = nys, nyn
820!
821!--          Subgridscale fluxes (without Prandtl layer from k=nzb,
822!--          oterwise from k=nzb+1)
[132]823!--          NOTE: for simplicity, nzb_diff_s_inner is used below, although
[1]824!--          ----  strictly speaking the following k-loop would have to be
825!--                split up according to the staggered grid.
[132]826!--                However, this implies no error since staggered velocity
827!--                components are zero at the walls and inside buildings.
[2232]828!--          Flag 23 is used to mask surface fluxes as well as model-top fluxes,
829!--          which are added further below.
830             DO  k = nzb, nzt
831                flag = MERGE( 1.0_wp, 0.0_wp,                                  &
832                              BTEST( wall_flags_0(k,j,i), 23 ) ) *             &
833                       MERGE( 1.0_wp, 0.0_wp,                                  &
834                              BTEST( wall_flags_0(k,j,i), 9  ) )
[1]835!
836!--             Momentum flux w"u"
[1353]837                sums_l(k,12,tn) = sums_l(k,12,tn) - 0.25_wp * (                &
[1]838                               km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) &
839                                                           ) * (               &
840                                   ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)     &
841                                 + ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
[2037]842                                                           ) * rmask(j,i,sr)   &
843                                         * rho_air_zw(k)                       &
[2232]844                                         * momentumflux_output_conversion(k)   &
845                                         * flag
[1]846!
847!--             Momentum flux w"v"
[1353]848                sums_l(k,14,tn) = sums_l(k,14,tn) - 0.25_wp * (                &
[1]849                               km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) &
850                                                           ) * (               &
851                                   ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
852                                 + ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
[2037]853                                                           ) * rmask(j,i,sr)   &
854                                         * rho_air_zw(k)                       &
[2232]855                                         * momentumflux_output_conversion(k)   &
856                                         * flag
[1]857!
858!--             Heat flux w"pt"
859                sums_l(k,16,tn) = sums_l(k,16,tn)                              &
[1353]860                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
[1]861                                               * ( pt(k+1,j,i) - pt(k,j,i) )   &
[2037]862                                               * rho_air_zw(k)                 &
863                                               * heatflux_output_conversion(k) &
[2232]864                                               * ddzu(k+1) * rmask(j,i,sr)     &
865                                               * flag
[1]866
867
868!
[96]869!--             Salinity flux w"sa"
870                IF ( ocean )  THEN
871                   sums_l(k,65,tn) = sums_l(k,65,tn)                           &
[1353]872                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
[96]873                                               * ( sa(k+1,j,i) - sa(k,j,i) )   &
[2232]874                                               * ddzu(k+1) * rmask(j,i,sr)     &
875                                               * flag
[96]876                ENDIF
877
878!
[1]879!--             Buoyancy flux, water flux (humidity flux) w"q"
[75]880                IF ( humidity ) THEN
[1]881                   sums_l(k,45,tn) = sums_l(k,45,tn)                           &
[1353]882                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
[1]883                                               * ( vpt(k+1,j,i) - vpt(k,j,i) ) &
[2037]884                                               * rho_air_zw(k)                 &
885                                               * heatflux_output_conversion(k) &
[2232]886                                               * ddzu(k+1) * rmask(j,i,sr) * flag
[1]887                   sums_l(k,48,tn) = sums_l(k,48,tn)                           &
[1353]888                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
[1]889                                               * ( q(k+1,j,i) - q(k,j,i) )     &
[2037]890                                               * rho_air_zw(k)                 &
891                                               * waterflux_output_conversion(k)&
[2232]892                                               * ddzu(k+1) * rmask(j,i,sr) * flag
[1007]893
[1]894                   IF ( cloud_physics ) THEN
895                      sums_l(k,51,tn) = sums_l(k,51,tn)                        &
[1353]896                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
[1]897                                               * ( ( q(k+1,j,i) - ql(k+1,j,i) )&
898                                                - ( q(k,j,i) - ql(k,j,i) ) )   &
[2037]899                                               * rho_air_zw(k)                 &
900                                               * waterflux_output_conversion(k)&
[2232]901                                               * ddzu(k+1) * rmask(j,i,sr) * flag
[1]902                   ENDIF
903                ENDIF
904
905!
906!--             Passive scalar flux
907                IF ( passive_scalar )  THEN
[2270]908                   sums_l(k,117,tn) = sums_l(k,117,tn)                         &
[1353]909                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
[2026]910                                                  * ( s(k+1,j,i) - s(k,j,i) )  &
[2232]911                                                  * ddzu(k+1) * rmask(j,i,sr)  &
912                                                              * flag
[1]913                ENDIF
914
915             ENDDO
916
917!
918!--          Subgridscale fluxes in the Prandtl layer
919             IF ( use_surface_fluxes )  THEN
[2232]920                DO  l = 0, 1
921                   ki = MERGE( -1, 0, l == 0 )
922                   IF ( surf_def_h(l)%ns >= 1 )  THEN
923                      DO  m = surf_def_h(l)%start_index(j,i), surf_def_h(l)%end_index(j,i)
924                         k = surf_def_h(l)%k(m)
925
926                         sums_l(k+ki,12,tn) = sums_l(k+ki,12,tn) + &
927                                    momentumflux_output_conversion(k+ki) * &
928                                    surf_def_h(l)%usws(m) * rmask(j,i,sr)     ! w"u"
929                         sums_l(k+ki,14,tn) = sums_l(k+ki,14,tn) + &
930                                    momentumflux_output_conversion(k+ki) * &
931                                    surf_def_h(l)%vsws(m) * rmask(j,i,sr)     ! w"v"
932                         sums_l(k+ki,16,tn) = sums_l(k+ki,16,tn) + &
933                                    heatflux_output_conversion(k+ki) * &
934                                    surf_def_h(l)%shf(m)  * rmask(j,i,sr)     ! w"pt"
935                         sums_l(k+ki,58,tn) = sums_l(k+ki,58,tn) + &
936                                    0.0_wp * rmask(j,i,sr)        ! u"pt"
937                         sums_l(k+ki,61,tn) = sums_l(k+ki,61,tn) + &
938                                    0.0_wp * rmask(j,i,sr)        ! v"pt"
939                         IF ( ocean )  THEN
940                            sums_l(k+ki,65,tn) = sums_l(k+ki,65,tn) + &
941                                       surf_def_h(l)%sasws(m) * rmask(j,i,sr)  ! w"sa"
942                         ENDIF
943                         IF ( humidity )  THEN
944                            sums_l(k+ki,48,tn) = sums_l(k+ki,48,tn) +                     &
945                                       waterflux_output_conversion(k+ki) *      &
946                                       surf_def_h(l)%qsws(m) * rmask(j,i,sr)  ! w"q" (w"qv")
947                            sums_l(k+ki,45,tn) = sums_l(k+ki,45,tn) + (                   &
948                                       ( 1.0_wp + 0.61_wp * q(k+ki,j,i) ) *     &
949                                       surf_def_h(l)%shf(m) + 0.61_wp * pt(k+ki,j,i) *      &
950                                                  surf_def_h(l)%qsws(m) )                  &
951                                       * heatflux_output_conversion(k+ki)
952                            IF ( cloud_droplets )  THEN
953                               sums_l(k+ki,45,tn) = sums_l(k+ki,45,tn) + (                &
954                                         ( 1.0_wp + 0.61_wp * q(k+ki,j,i) -     &
955                                           ql(k+ki,j,i) ) * surf_def_h(l)%shf(m) +          &
956                                           0.61_wp * pt(k+ki,j,i) * surf_def_h(l)%qsws(m) ) &
957                                          * heatflux_output_conversion(k+ki)
958                            ENDIF
959                            IF ( cloud_physics )  THEN
960!
961!--                            Formula does not work if ql(k+ki) /= 0.0
962                               sums_l(k+ki,51,tn) = sums_l(k+ki,51,tn) +                  &
963                                          waterflux_output_conversion(k+ki) *   &
964                                          surf_def_h(l)%qsws(m) * rmask(j,i,sr) ! w"q" (w"qv")
965                            ENDIF
966                         ENDIF
967                         IF ( passive_scalar )  THEN
[2270]968                            sums_l(k+ki,117,tn) = sums_l(k+ki,117,tn) +                     &
[2232]969                                        surf_def_h(l)%ssws(m) * rmask(j,i,sr) ! w"s"
970                         ENDIF
971
972                      ENDDO
973
974                   ENDIF
975                ENDDO
976                IF ( surf_lsm_h%ns >= 1 )  THEN
977                   m = surf_lsm_h%start_index(j,i)
978                   sums_l(nzb,12,tn) = sums_l(nzb,12,tn) + &
[2037]979                                    momentumflux_output_conversion(nzb) * &
[2232]980                                    surf_lsm_h%usws(m) * rmask(j,i,sr)     ! w"u"
981                   sums_l(nzb,14,tn) = sums_l(nzb,14,tn) + &
[2037]982                                    momentumflux_output_conversion(nzb) * &
[2232]983                                    surf_lsm_h%vsws(m) * rmask(j,i,sr)     ! w"v"
984                   sums_l(nzb,16,tn) = sums_l(nzb,16,tn) + &
[2037]985                                    heatflux_output_conversion(nzb) * &
[2232]986                                    surf_lsm_h%shf(m)  * rmask(j,i,sr)     ! w"pt"
987                   sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + &
[1353]988                                    0.0_wp * rmask(j,i,sr)        ! u"pt"
[2232]989                   sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + &
[1353]990                                    0.0_wp * rmask(j,i,sr)        ! v"pt"
[2232]991                   IF ( ocean )  THEN
992                      sums_l(nzb,65,tn) = sums_l(nzb,65,tn) + &
993                                       surf_lsm_h%sasws(m) * rmask(j,i,sr)  ! w"sa"
994                   ENDIF
995                   IF ( humidity )  THEN
996                      sums_l(nzb,48,tn) = sums_l(nzb,48,tn) +                     &
997                                       waterflux_output_conversion(nzb) *      &
998                                       surf_lsm_h%qsws(m) * rmask(j,i,sr)  ! w"q" (w"qv")
999                      sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                   &
1000                                       ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) *     &
1001                                       surf_lsm_h%shf(m) + 0.61_wp * pt(nzb,j,i) *      &
1002                                                  surf_lsm_h%qsws(m) )                  &
1003                                       * heatflux_output_conversion(nzb)
1004                      IF ( cloud_droplets )  THEN
1005                         sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                &
1006                                         ( 1.0_wp + 0.61_wp * q(nzb,j,i) -     &
1007                                           ql(nzb,j,i) ) * surf_lsm_h%shf(m) +          &
1008                                           0.61_wp * pt(nzb,j,i) * surf_lsm_h%qsws(m) ) &
1009                                          * heatflux_output_conversion(nzb)
1010                      ENDIF
1011                      IF ( cloud_physics )  THEN
1012!
1013!--                      Formula does not work if ql(nzb) /= 0.0
1014                         sums_l(nzb,51,tn) = sums_l(nzb,51,tn) +                  &
1015                                          waterflux_output_conversion(nzb) *   &
1016                                          surf_lsm_h%qsws(m) * rmask(j,i,sr) ! w"q" (w"qv")
1017                      ENDIF
1018                   ENDIF
1019                   IF ( passive_scalar )  THEN
[2270]1020                      sums_l(nzb,117,tn) = sums_l(nzb,117,tn) +                     &
[2232]1021                                        surf_lsm_h%ssws(m) * rmask(j,i,sr) ! w"s"
1022                   ENDIF
1023
1024
[96]1025                ENDIF
[2232]1026                IF ( surf_usm_h%ns >= 1 )  THEN
1027                   m = surf_usm_h%start_index(j,i)
1028                   sums_l(nzb,12,tn) = sums_l(nzb,12,tn) + &
1029                                    momentumflux_output_conversion(nzb) * &
1030                                    surf_usm_h%usws(m) * rmask(j,i,sr)     ! w"u"
1031                   sums_l(nzb,14,tn) = sums_l(nzb,14,tn) + &
1032                                    momentumflux_output_conversion(nzb) * &
1033                                    surf_usm_h%vsws(m) * rmask(j,i,sr)     ! w"v"
1034                   sums_l(nzb,16,tn) = sums_l(nzb,16,tn) + &
1035                                    heatflux_output_conversion(nzb) * &
1036                                    surf_usm_h%shf(m)  * rmask(j,i,sr)     ! w"pt"
1037                   sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + &
1038                                    0.0_wp * rmask(j,i,sr)        ! u"pt"
1039                   sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + &
1040                                    0.0_wp * rmask(j,i,sr)        ! v"pt"
1041                   IF ( ocean )  THEN
1042                      sums_l(nzb,65,tn) = sums_l(nzb,65,tn) + &
1043                                       surf_usm_h%sasws(m) * rmask(j,i,sr)  ! w"sa"
1044                   ENDIF
1045                   IF ( humidity )  THEN
1046                      sums_l(nzb,48,tn) = sums_l(nzb,48,tn) +                     &
[2037]1047                                       waterflux_output_conversion(nzb) *      &
[2232]1048                                       surf_usm_h%qsws(m) * rmask(j,i,sr)  ! w"q" (w"qv")
1049                      sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                   &
[1353]1050                                       ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) *     &
[2232]1051                                       surf_usm_h%shf(m) + 0.61_wp * pt(nzb,j,i) *      &
1052                                                  surf_usm_h%qsws(m) )                  &
[2037]1053                                       * heatflux_output_conversion(nzb)
[2232]1054                      IF ( cloud_droplets )  THEN
1055                         sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                &
[1353]1056                                         ( 1.0_wp + 0.61_wp * q(nzb,j,i) -     &
[2232]1057                                           ql(nzb,j,i) ) * surf_usm_h%shf(m) +          &
1058                                           0.61_wp * pt(nzb,j,i) * surf_usm_h%qsws(m) ) &
[2037]1059                                          * heatflux_output_conversion(nzb)
[2232]1060                      ENDIF
1061                      IF ( cloud_physics )  THEN
[1]1062!
[2232]1063!--                      Formula does not work if ql(nzb) /= 0.0
1064                         sums_l(nzb,51,tn) = sums_l(nzb,51,tn) +                  &
[2037]1065                                          waterflux_output_conversion(nzb) *   &
[2232]1066                                          surf_usm_h%qsws(m) * rmask(j,i,sr) ! w"q" (w"qv")
1067                      ENDIF
[1]1068                   ENDIF
[2232]1069                   IF ( passive_scalar )  THEN
[2270]1070                      sums_l(nzb,117,tn) = sums_l(nzb,117,tn) +                     &
[2232]1071                                        surf_usm_h%ssws(m) * rmask(j,i,sr) ! w"s"
1072                   ENDIF
1073
1074
[1]1075                ENDIF
[2232]1076
[1]1077             ENDIF
1078
[1691]1079             IF ( .NOT. neutral )  THEN
[2232]1080                IF ( surf_def_h(0)%ns >= 1 )  THEN
1081                   m = surf_def_h(0)%start_index(j,i)
[2270]1082                   sums_l(nzb,112,tn) = sums_l(nzb,112,tn) +                      &
[2232]1083                                        surf_def_h(0)%ol(m)  * rmask(j,i,sr) ! L
1084                ENDIF
1085                IF ( surf_lsm_h%ns >= 1 )  THEN
1086                   m = surf_lsm_h%start_index(j,i)
[2270]1087                   sums_l(nzb,112,tn) = sums_l(nzb,112,tn) +                      &
[2232]1088                                        surf_lsm_h%ol(m)  * rmask(j,i,sr) ! L
1089                ENDIF
1090                IF ( surf_usm_h%ns >= 1 )  THEN
1091                   m = surf_usm_h%start_index(j,i)
[2270]1092                   sums_l(nzb,112,tn) = sums_l(nzb,112,tn) +                      &
[2232]1093                                        surf_usm_h%ol(m)  * rmask(j,i,sr) ! L
1094                ENDIF
[1691]1095             ENDIF
1096
[2296]1097             IF ( radiation )  THEN
[2270]1098                sums_l(nzb,99,tn)  = sums_l(nzb,99,tn)  + rad_net(j,i)
1099                sums_l(nzb,100,tn)  = sums_l(nzb,100,tn)  + rad_lw_in(nzb,j,i)
1100                sums_l(nzb,101,tn)  = sums_l(nzb,101,tn)  + rad_lw_out(nzb,j,i)
1101                sums_l(nzb,102,tn)  = sums_l(nzb,102,tn)  + rad_sw_in(nzb,j,i)
1102                sums_l(nzb,103,tn)  = sums_l(nzb,103,tn)  + rad_sw_out(nzb,j,i)
[1585]1103
1104#if defined ( __rrtmg )
1105                IF ( radiation_scheme == 'rrtmg' )  THEN
[2270]1106                   sums_l(nzb,108,tn)  = sums_l(nzb,108,tn)  + rrtm_aldif(0,j,i)
1107                   sums_l(nzb,109,tn)  = sums_l(nzb,109,tn)  + rrtm_aldir(0,j,i)
1108                   sums_l(nzb,110,tn)  = sums_l(nzb,110,tn)  + rrtm_asdif(0,j,i)
1109                   sums_l(nzb,111,tn)  = sums_l(nzb,111,tn)  + rrtm_asdir(0,j,i)
[1585]1110                ENDIF
1111#endif
[1551]1112             ENDIF
[1]1113!
[19]1114!--          Subgridscale fluxes at the top surface
1115             IF ( use_top_fluxes )  THEN
[2232]1116                m = surf_def_h(2)%start_index(j,i)
[550]1117                sums_l(nzt:nzt+1,12,tn) = sums_l(nzt:nzt+1,12,tn) + &
[2037]1118                                    momentumflux_output_conversion(nzt:nzt+1) * &
[2232]1119                                    surf_def_h(2)%usws(m) * rmask(j,i,sr)    ! w"u"
[550]1120                sums_l(nzt:nzt+1,14,tn) = sums_l(nzt:nzt+1,14,tn) + &
[2037]1121                                    momentumflux_output_conversion(nzt:nzt+1) * &
[2232]1122                                    surf_def_h(2)%vsws(m) * rmask(j,i,sr)    ! w"v"
[550]1123                sums_l(nzt:nzt+1,16,tn) = sums_l(nzt:nzt+1,16,tn) + &
[2037]1124                                    heatflux_output_conversion(nzt:nzt+1) * &
[2232]1125                                    surf_def_h(2)%shf(m)  * rmask(j,i,sr)   ! w"pt"
[550]1126                sums_l(nzt:nzt+1,58,tn) = sums_l(nzt:nzt+1,58,tn) + &
[1353]1127                                    0.0_wp * rmask(j,i,sr)        ! u"pt"
[550]1128                sums_l(nzt:nzt+1,61,tn) = sums_l(nzt:nzt+1,61,tn) + &
[1353]1129                                    0.0_wp * rmask(j,i,sr)        ! v"pt"
[550]1130
[96]1131                IF ( ocean )  THEN
1132                   sums_l(nzt,65,tn) = sums_l(nzt,65,tn) + &
[2232]1133                                       surf_def_h(2)%sasws(m) * rmask(j,i,sr)  ! w"sa"
[96]1134                ENDIF
[75]1135                IF ( humidity )  THEN
[1353]1136                   sums_l(nzt,48,tn) = sums_l(nzt,48,tn) +                     &
[2037]1137                                       waterflux_output_conversion(nzt) *      &
[2232]1138                                       surf_def_h(2)%qsws(m) * rmask(j,i,sr) ! w"q" (w"qv")
[1353]1139                   sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (                   &
1140                                       ( 1.0_wp + 0.61_wp * q(nzt,j,i) ) *     &
[2232]1141                                       surf_def_h(2)%shf(m) +                  &
1142                                       0.61_wp * pt(nzt,j,i) *    &
1143                                       surf_def_h(2)%qsws(m) )      &
[2037]1144                                       * heatflux_output_conversion(nzt)
[1007]1145                   IF ( cloud_droplets )  THEN
[1353]1146                      sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (                &
1147                                          ( 1.0_wp + 0.61_wp * q(nzt,j,i) -    &
[2232]1148                                            ql(nzt,j,i) ) *                    &
1149                                            surf_def_h(2)%shf(m) +             &
1150                                           0.61_wp * pt(nzt,j,i) *             &
1151                                           surf_def_h(2)%qsws(m) )&
[2037]1152                                           * heatflux_output_conversion(nzt)
[1007]1153                   ENDIF
[19]1154                   IF ( cloud_physics )  THEN
1155!
1156!--                   Formula does not work if ql(nzb) /= 0.0
1157                      sums_l(nzt,51,tn) = sums_l(nzt,51,tn) + &   ! w"q" (w"qv")
[2037]1158                                          waterflux_output_conversion(nzt) *   &
[2232]1159                                          surf_def_h(2)%qsws(m) * rmask(j,i,sr)
[19]1160                   ENDIF
1161                ENDIF
1162                IF ( passive_scalar )  THEN
[2270]1163                   sums_l(nzt,117,tn) = sums_l(nzt,117,tn) + &
[2232]1164                                        surf_def_h(2)%ssws(m) * rmask(j,i,sr) ! w"s"
[19]1165                ENDIF
1166             ENDIF
1167
1168!
[1]1169!--          Resolved fluxes (can be computed for all horizontal points)
[132]1170!--          NOTE: for simplicity, nzb_s_inner is used below, although strictly
[1]1171!--          ----  speaking the following k-loop would have to be split up and
1172!--                rearranged according to the staggered grid.
[2232]1173             DO  k = nzb, nzt
1174                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 22 ) )
[1353]1175                ust = 0.5_wp * ( u(k,j,i)   - hom(k,1,1,sr) +                  &
1176                                 u(k+1,j,i) - hom(k+1,1,1,sr) )
1177                vst = 0.5_wp * ( v(k,j,i)   - hom(k,1,2,sr) +                  &
1178                                 v(k+1,j,i) - hom(k+1,1,2,sr) )
1179                pts = 0.5_wp * ( pt(k,j,i)   - hom(k,1,4,sr) +                 &
1180                                 pt(k+1,j,i) - hom(k+1,1,4,sr) )
[667]1181
[1]1182!--             Higher moments
[1353]1183                sums_l(k,35,tn) = sums_l(k,35,tn) + pts * w(k,j,i)**2 *        &
[2232]1184                                                    rmask(j,i,sr) * flag
[1353]1185                sums_l(k,36,tn) = sums_l(k,36,tn) + pts**2 * w(k,j,i) *        &
[2232]1186                                                    rmask(j,i,sr) * flag
[1]1187
1188!
[96]1189!--             Salinity flux and density (density does not belong to here,
[97]1190!--             but so far there is no other suitable place to calculate)
[96]1191                IF ( ocean )  THEN
[1567]1192                   IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
[1353]1193                      pts = 0.5_wp * ( sa(k,j,i)   - hom(k,1,23,sr) +          &
1194                                       sa(k+1,j,i) - hom(k+1,1,23,sr) )
1195                      sums_l(k,66,tn) = sums_l(k,66,tn) + pts * w(k,j,i) *     &
[2232]1196                                        rmask(j,i,sr) * flag
[667]1197                   ENDIF
[2232]1198                   sums_l(k,64,tn) = sums_l(k,64,tn) + rho_ocean(k,j,i) *      &
1199                                                       rmask(j,i,sr) * flag
[1353]1200                   sums_l(k,71,tn) = sums_l(k,71,tn) + prho(k,j,i) *           &
[2232]1201                                                       rmask(j,i,sr) * flag
[96]1202                ENDIF
1203
1204!
[1053]1205!--             Buoyancy flux, water flux, humidity flux, liquid water
1206!--             content, rain drop concentration and rain water content
[75]1207                IF ( humidity )  THEN
[1007]1208                   IF ( cloud_physics .OR. cloud_droplets )  THEN
[1353]1209                      pts = 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +         &
[1007]1210                                    vpt(k+1,j,i) - hom(k+1,1,44,sr) )
[1353]1211                      sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) *     &
[2037]1212                                               heatflux_output_conversion(k) * &
[2232]1213                                                          rmask(j,i,sr) * flag
1214                      sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) * rmask(j,i,sr) &
1215                                                                    * flag
[1822]1216
[1053]1217                      IF ( .NOT. cloud_droplets )  THEN
[1353]1218                         pts = 0.5_wp *                                        &
[1115]1219                              ( ( q(k,j,i) - ql(k,j,i) ) -                     &
1220                              hom(k,1,42,sr) +                                 &
1221                              ( q(k+1,j,i) - ql(k+1,j,i) ) -                   &
[1053]1222                              hom(k+1,1,42,sr) )
[1115]1223                         sums_l(k,52,tn) = sums_l(k,52,tn) + pts * w(k,j,i) *  &
[2037]1224                                             waterflux_output_conversion(k) *  &
[2232]1225                                                             rmask(j,i,sr)  *  &
1226                                                             flag
[1822]1227                         sums_l(k,75,tn) = sums_l(k,75,tn) + qc(k,j,i) *       &
[2232]1228                                                             rmask(j,i,sr) *   &
1229                                                             flag
[1822]1230                         sums_l(k,76,tn) = sums_l(k,76,tn) + prr(k,j,i) *      &
[2232]1231                                                             rmask(j,i,sr) *   &
1232                                                             flag
[2292]1233                         IF ( microphysics_morrison )  THEN
1234                            sums_l(k,123,tn) = sums_l(k,123,tn) + nc(k,j,i) *  &
1235                                                                rmask(j,i,sr) *&
1236                                                                flag
1237                         ENDIF
[1822]1238                         IF ( microphysics_seifert )  THEN
1239                            sums_l(k,73,tn) = sums_l(k,73,tn) + nr(k,j,i) *    &
[2232]1240                                                                rmask(j,i,sr) *&
1241                                                                flag
[1822]1242                            sums_l(k,74,tn) = sums_l(k,74,tn) + qr(k,j,i) *    &
[2232]1243                                                                rmask(j,i,sr) *&
1244                                                                flag
[1053]1245                         ENDIF
1246                      ENDIF
[1822]1247
[1007]1248                   ELSE
[1567]1249                      IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
[1353]1250                         pts = 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +      &
1251                                          vpt(k+1,j,i) - hom(k+1,1,44,sr) )
1252                         sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) *  &
[2037]1253                                              heatflux_output_conversion(k) *  &
[2232]1254                                                             rmask(j,i,sr)  *  &
1255                                                             flag
[1567]1256                      ELSE IF ( ws_scheme_sca .AND. sr == 0 )  THEN
[2037]1257                         sums_l(k,46,tn) = ( ( 1.0_wp + 0.61_wp *              & 
1258                                               hom(k,1,41,sr) ) *              &
1259                                             sums_l(k,17,tn) +                 &
1260                                             0.61_wp * hom(k,1,4,sr) *         &
1261                                             sums_l(k,49,tn)                   &
[2232]1262                                           ) * heatflux_output_conversion(k) * &
1263                                               flag
[1007]1264                      END IF
1265                   END IF
[1]1266                ENDIF
1267!
1268!--             Passive scalar flux
[1353]1269                IF ( passive_scalar .AND. ( .NOT. ws_scheme_sca                &
[1567]1270                     .OR. sr /= 0 ) )  THEN
[2270]1271                   pts = 0.5_wp * ( s(k,j,i)   - hom(k,1,115,sr) +             &
1272                                    s(k+1,j,i) - hom(k+1,1,115,sr) )
1273                   sums_l(k,114,tn) = sums_l(k,114,tn) + pts * w(k,j,i) *      &
[2232]1274                                                       rmask(j,i,sr) * flag
[1]1275                ENDIF
1276
1277!
1278!--             Energy flux w*e*
[667]1279!--             has to be adjusted
[1353]1280                sums_l(k,37,tn) = sums_l(k,37,tn) + w(k,j,i) * 0.5_wp *        &
1281                                             ( ust**2 + vst**2 + w(k,j,i)**2 ) &
[2037]1282                                           * momentumflux_output_conversion(k) &
[2232]1283                                           * rmask(j,i,sr) * flag
[1]1284             ENDDO
1285          ENDDO
1286       ENDDO
[2232]1287       !$OMP END PARALLEL
[709]1288!
[2232]1289!--    Treat land-surface quantities according to new wall model structure.
1290       IF ( land_surface )  THEN
1291          tn = 0
1292          !$OMP PARALLEL PRIVATE( i, j, m, tn )
1293          !$ tn = omp_get_thread_num()
1294          !$OMP DO
1295          DO  m = 1, surf_lsm_h%ns
1296             i = surf_lsm_h%i(m)
1297             j = surf_lsm_h%j(m)
1298       
1299             IF ( i >= nxl  .AND.  i <= nxr  .AND.                             &
1300                  j >= nys  .AND.  j <= nyn )  THEN
[2270]1301                sums_l(nzb,93,tn)  = sums_l(nzb,93,tn) + surf_lsm_h%ghf(m)
1302                sums_l(nzb,94,tn)  = sums_l(nzb,94,tn) + surf_lsm_h%qsws_liq(m)
1303                sums_l(nzb,95,tn)  = sums_l(nzb,95,tn) + surf_lsm_h%qsws_soil(m)
1304                sums_l(nzb,96,tn)  = sums_l(nzb,96,tn) + surf_lsm_h%qsws_veg(m)
1305                sums_l(nzb,97,tn)  = sums_l(nzb,97,tn) + surf_lsm_h%r_a(m)
1306                sums_l(nzb,98,tn) = sums_l(nzb,98,tn)+ surf_lsm_h%r_s(m)
[2232]1307             ENDIF
1308          ENDDO
1309          !$OMP END PARALLEL
1310
1311          tn = 0
1312          !$OMP PARALLEL PRIVATE( i, j, k, m, tn )
1313          !$ tn = omp_get_thread_num()
1314          !$OMP DO
1315          DO  m = 1, surf_lsm_h%ns
1316
1317             i = surf_lsm_h%i(m)           
1318             j = surf_lsm_h%j(m)
1319
1320             IF ( i >= nxl  .AND.  i <= nxr  .AND.                             &
1321                  j >= nys  .AND.  j <= nyn )  THEN
1322
1323                DO  k = nzb_soil, nzt_soil
1324                   sums_l(k,89,tn)  = sums_l(k,89,tn)  + t_soil_h%var_2d(k,m)  &
1325                                      * rmask(j,i,sr)
1326                   sums_l(k,91,tn)  = sums_l(k,91,tn)  + m_soil_h%var_2d(k,m)  &
1327                                      * rmask(j,i,sr)
1328                ENDDO
1329             ENDIF
1330          ENDDO
1331          !$OMP END PARALLEL
1332       ENDIF
1333!
[709]1334!--    For speed optimization fluxes which have been computed in part directly
1335!--    inside the WS advection routines are treated seperatly
1336!--    Momentum fluxes first:
[2232]1337
1338       tn = 0
1339       !$OMP PARALLEL PRIVATE( i, j, k, tn, flag )
1340       !$ tn = omp_get_thread_num()
[743]1341       IF ( .NOT. ws_scheme_mom .OR. sr /= 0  )  THEN
[2232]1342          !$OMP DO
1343          DO  i = nxl, nxr
1344             DO  j = nys, nyn
1345                DO  k = nzb, nzt
[1007]1346!
[2232]1347!--                Flag 23 is used to mask surface fluxes as well as model-top
1348!--                fluxes, which are added further below.
1349                   flag = MERGE( 1.0_wp, 0.0_wp,                               &
1350                                 BTEST( wall_flags_0(k,j,i), 23 ) ) *          &
1351                          MERGE( 1.0_wp, 0.0_wp,                               &
1352                                 BTEST( wall_flags_0(k,j,i), 9  ) )
1353
1354                   ust = 0.5_wp * ( u(k,j,i)   - hom(k,1,1,sr) +               &
1355                                    u(k+1,j,i) - hom(k+1,1,1,sr) )
1356                   vst = 0.5_wp * ( v(k,j,i)   - hom(k,1,2,sr) +               &
1357                                    v(k+1,j,i) - hom(k+1,1,2,sr) )
[667]1358!
[2232]1359!--                Momentum flux w*u*
1360                   sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5_wp *                &
1361                                                     ( w(k,j,i-1) + w(k,j,i) ) &
1362                                           * momentumflux_output_conversion(k) &
1363                                                     * ust * rmask(j,i,sr)     &
1364                                                           * flag
1365!
1366!--                Momentum flux w*v*
1367                   sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5_wp *                &
1368                                                     ( w(k,j-1,i) + w(k,j,i) ) &
1369                                           * momentumflux_output_conversion(k) &
1370                                                     * vst * rmask(j,i,sr)     &
1371                                                           * flag
1372                ENDDO
1373             ENDDO
1374          ENDDO
[1]1375
[667]1376       ENDIF
[1567]1377       IF ( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
[2232]1378          !$OMP DO
1379          DO  i = nxl, nxr
1380             DO  j = nys, nyn
1381                DO  k = nzb, nzt
1382                   flag = MERGE( 1.0_wp, 0.0_wp,                               &
1383                                 BTEST( wall_flags_0(k,j,i), 23 ) ) *          &
1384                          MERGE( 1.0_wp, 0.0_wp,                               &
1385                                 BTEST( wall_flags_0(k,j,i), 9  ) )
[709]1386!
[2232]1387!--                Vertical heat flux
1388                   sums_l(k,17,tn) = sums_l(k,17,tn) + 0.5_wp *                &
[1353]1389                           ( pt(k,j,i)   - hom(k,1,4,sr) +                     &
1390                             pt(k+1,j,i) - hom(k+1,1,4,sr) )                   &
[2037]1391                           * heatflux_output_conversion(k)                     &
[2232]1392                           * w(k,j,i) * rmask(j,i,sr) * flag
1393                   IF ( humidity )  THEN
1394                      pts = 0.5_wp * ( q(k,j,i)   - hom(k,1,41,sr) +           &
[1353]1395                                      q(k+1,j,i) - hom(k+1,1,41,sr) )
[2232]1396                      sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) *     &
[2037]1397                                       waterflux_output_conversion(k) *        &
[2232]1398                                       rmask(j,i,sr) * flag
1399                   ENDIF
1400                   IF ( passive_scalar )  THEN
[2270]1401                      pts = 0.5_wp * ( s(k,j,i)   - hom(k,1,115,sr) +           &
1402                                      s(k+1,j,i) - hom(k+1,1,115,sr) )
1403                      sums_l(k,114,tn) = sums_l(k,114,tn) + pts * w(k,j,i) *    &
[2232]1404                                        rmask(j,i,sr) * flag
1405                   ENDIF
1406                ENDDO
1407             ENDDO
1408          ENDDO
[667]1409
1410       ENDIF
1411
[1]1412!
[97]1413!--    Density at top follows Neumann condition
[388]1414       IF ( ocean )  THEN
1415          sums_l(nzt+1,64,tn) = sums_l(nzt,64,tn)
1416          sums_l(nzt+1,71,tn) = sums_l(nzt,71,tn)
1417       ENDIF
[97]1418
1419!
[1]1420!--    Divergence of vertical flux of resolved scale energy and pressure
[106]1421!--    fluctuations as well as flux of pressure fluctuation itself (68).
1422!--    First calculate the products, then the divergence.
[1]1423!--    Calculation is time consuming. Do it only, if profiles shall be plotted.
[1691]1424       IF ( hom(nzb+1,2,55,0) /= 0.0_wp  .OR.  hom(nzb+1,2,68,0) /= 0.0_wp )   &
1425       THEN
[1353]1426          sums_ll = 0.0_wp  ! local array
[1]1427
1428          !$OMP DO
1429          DO  i = nxl, nxr
1430             DO  j = nys, nyn
[2232]1431                DO  k = nzb+1, nzt
1432                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
[1]1433
[1353]1434                   sums_ll(k,1) = sums_ll(k,1) + 0.5_wp * w(k,j,i) * (         &
[1652]1435                  ( 0.25_wp * ( u(k,j,i)+u(k+1,j,i)+u(k,j,i+1)+u(k+1,j,i+1) )  &
1436                            - 0.5_wp * ( hom(k,1,1,sr) + hom(k+1,1,1,sr) ) )**2&
1437                + ( 0.25_wp * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i) )  &
[1654]1438                            - 0.5_wp * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) ) )**2&
[2232]1439                + w(k,j,i)**2                                        ) * flag
[1]1440
[1353]1441                   sums_ll(k,2) = sums_ll(k,2) + 0.5_wp * w(k,j,i)             &
[2252]1442                                       * ( ( p(k,j,i) + p(k+1,j,i) )           &
1443                                         / momentumflux_output_conversion(k) ) &
1444                                       * flag
[1]1445
1446                ENDDO
1447             ENDDO
1448          ENDDO
[1353]1449          sums_ll(0,1)     = 0.0_wp    ! because w is zero at the bottom
1450          sums_ll(nzt+1,1) = 0.0_wp
1451          sums_ll(0,2)     = 0.0_wp
1452          sums_ll(nzt+1,2) = 0.0_wp
[1]1453
[678]1454          DO  k = nzb+1, nzt
[1]1455             sums_l(k,55,tn) = ( sums_ll(k,1) - sums_ll(k-1,1) ) * ddzw(k)
1456             sums_l(k,56,tn) = ( sums_ll(k,2) - sums_ll(k-1,2) ) * ddzw(k)
[106]1457             sums_l(k,68,tn) = sums_ll(k,2)
[1]1458          ENDDO
1459          sums_l(nzb,55,tn) = sums_l(nzb+1,55,tn)
1460          sums_l(nzb,56,tn) = sums_l(nzb+1,56,tn)
[1353]1461          sums_l(nzb,68,tn) = 0.0_wp    ! because w* = 0 at nzb
[1]1462
1463       ENDIF
1464
1465!
[106]1466!--    Divergence of vertical flux of SGS TKE and the flux itself (69)
[1691]1467       IF ( hom(nzb+1,2,57,0) /= 0.0_wp  .OR.  hom(nzb+1,2,69,0) /= 0.0_wp )   &
1468       THEN
[1]1469          !$OMP DO
1470          DO  i = nxl, nxr
1471             DO  j = nys, nyn
[2232]1472                DO  k = nzb+1, nzt
[1]1473
[2232]1474                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1475
[1353]1476                   sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5_wp * (              &
[1]1477                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
1478                 - (km(k-1,j,i)+km(k,j,i)) * (e(k,j,i)-e(k-1,j,i)) * ddzu(k)   &
[2232]1479                                                                ) * ddzw(k)    &
1480                                                                  * flag
[1]1481
[1353]1482                   sums_l(k,69,tn) = sums_l(k,69,tn) - 0.5_wp * (              &
[106]1483                   (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
[2232]1484                                                                )  * flag
[106]1485
[1]1486                ENDDO
1487             ENDDO
1488          ENDDO
1489          sums_l(nzb,57,tn) = sums_l(nzb+1,57,tn)
[106]1490          sums_l(nzb,69,tn) = sums_l(nzb+1,69,tn)
[1]1491
1492       ENDIF
1493
1494!
1495!--    Horizontal heat fluxes (subgrid, resolved, total).
1496!--    Do it only, if profiles shall be plotted.
[1353]1497       IF ( hom(nzb+1,2,58,0) /= 0.0_wp ) THEN
[1]1498
1499          !$OMP DO
1500          DO  i = nxl, nxr
1501             DO  j = nys, nyn
[2232]1502                DO  k = nzb+1, nzt
1503                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
[1]1504!
1505!--                Subgrid horizontal heat fluxes u"pt", v"pt"
[1353]1506                   sums_l(k,58,tn) = sums_l(k,58,tn) - 0.5_wp *                &
[1]1507                                                   ( kh(k,j,i) + kh(k,j,i-1) ) &
1508                                                 * ( pt(k,j,i-1) - pt(k,j,i) ) &
[2037]1509                                               * rho_air_zw(k)                 &
1510                                               * heatflux_output_conversion(k) &
[2232]1511                                                 * ddx * rmask(j,i,sr) * flag
[1353]1512                   sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5_wp *                &
[1]1513                                                   ( kh(k,j,i) + kh(k,j-1,i) ) &
1514                                                 * ( pt(k,j-1,i) - pt(k,j,i) ) &
[2037]1515                                               * rho_air_zw(k)                 &
1516                                               * heatflux_output_conversion(k) &
[2232]1517                                                 * ddy * rmask(j,i,sr) * flag
[1]1518!
1519!--                Resolved horizontal heat fluxes u*pt*, v*pt*
1520                   sums_l(k,59,tn) = sums_l(k,59,tn) +                         &
1521                                                  ( u(k,j,i) - hom(k,1,1,sr) ) &
[1353]1522                                    * 0.5_wp * ( pt(k,j,i-1) - hom(k,1,4,sr) + &
[2037]1523                                                 pt(k,j,i)   - hom(k,1,4,sr) ) &
[2232]1524                                               * heatflux_output_conversion(k) &
1525                                               * flag
[1353]1526                   pts = 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) +              &
1527                                    pt(k,j,i)   - hom(k,1,4,sr) )
[1]1528                   sums_l(k,62,tn) = sums_l(k,62,tn) +                         &
1529                                                  ( v(k,j,i) - hom(k,1,2,sr) ) &
[1353]1530                                    * 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
[2037]1531                                                 pt(k,j,i)   - hom(k,1,4,sr) ) &
[2232]1532                                               * heatflux_output_conversion(k) &
1533                                               * flag
[1]1534                ENDDO
1535             ENDDO
1536          ENDDO
1537!
1538!--       Fluxes at the surface must be zero (e.g. due to the Prandtl-layer)
[1353]1539          sums_l(nzb,58,tn) = 0.0_wp
1540          sums_l(nzb,59,tn) = 0.0_wp
1541          sums_l(nzb,60,tn) = 0.0_wp
1542          sums_l(nzb,61,tn) = 0.0_wp
1543          sums_l(nzb,62,tn) = 0.0_wp
1544          sums_l(nzb,63,tn) = 0.0_wp
[1]1545
1546       ENDIF
[2073]1547       !$OMP END PARALLEL
[87]1548
1549!
[1365]1550!--    Collect current large scale advection and subsidence tendencies for
1551!--    data output
[1691]1552       IF ( large_scale_forcing  .AND.  ( simulated_time > 0.0_wp ) )  THEN
[1365]1553!
1554!--       Interpolation in time of LSF_DATA
1555          nt = 1
[1386]1556          DO WHILE ( simulated_time - dt_3d > time_vert(nt) )
[1365]1557             nt = nt + 1
1558          ENDDO
[1386]1559          IF ( simulated_time - dt_3d /= time_vert(nt) )  THEN
[1365]1560            nt = nt - 1
1561          ENDIF
1562
[1386]1563          fac = ( simulated_time - dt_3d - time_vert(nt) )                     &
[1365]1564                / ( time_vert(nt+1)-time_vert(nt) )
1565
1566
1567          DO  k = nzb, nzt
[1382]1568             sums_ls_l(k,0) = td_lsa_lpt(k,nt)                                 &
1569                              + fac * ( td_lsa_lpt(k,nt+1) - td_lsa_lpt(k,nt) )
1570             sums_ls_l(k,1) = td_lsa_q(k,nt)                                   &
1571                              + fac * ( td_lsa_q(k,nt+1) - td_lsa_q(k,nt) )
[1365]1572          ENDDO
1573
[1382]1574          sums_ls_l(nzt+1,0) = sums_ls_l(nzt,0)
1575          sums_ls_l(nzt+1,1) = sums_ls_l(nzt,1)
1576
[1365]1577          IF ( large_scale_subsidence .AND. use_subsidence_tendencies )  THEN
1578
1579             DO  k = nzb, nzt
[1382]1580                sums_ls_l(k,2) = td_sub_lpt(k,nt) + fac *                      &
1581                                 ( td_sub_lpt(k,nt+1) - td_sub_lpt(k,nt) )
1582                sums_ls_l(k,3) = td_sub_q(k,nt) + fac *                        &
1583                                 ( td_sub_q(k,nt+1) - td_sub_q(k,nt) )
[1365]1584             ENDDO
1585
[1382]1586             sums_ls_l(nzt+1,2) = sums_ls_l(nzt,2)
1587             sums_ls_l(nzt+1,3) = sums_ls_l(nzt,3)
1588
[1365]1589          ENDIF
1590
1591       ENDIF
1592
[2232]1593       tn = 0
[2073]1594       !$OMP PARALLEL PRIVATE( i, j, k, tn )
[2232]1595       !$ tn = omp_get_thread_num()       
[1585]1596       IF ( radiation .AND. radiation_scheme == 'rrtmg' )  THEN
1597          !$OMP DO
1598          DO  i = nxl, nxr
1599             DO  j =  nys, nyn
[2232]1600                DO  k = nzb+1, nzt+1
1601                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
1602
[2270]1603                   sums_l(k,100,tn)  = sums_l(k,100,tn)  + rad_lw_in(k,j,i)    &
[2232]1604                                       * rmask(j,i,sr) * flag
[2270]1605                   sums_l(k,101,tn)  = sums_l(k,101,tn)  + rad_lw_out(k,j,i)   &
[2232]1606                                       * rmask(j,i,sr) * flag
[2270]1607                   sums_l(k,102,tn)  = sums_l(k,102,tn)  + rad_sw_in(k,j,i)    &
[2232]1608                                       * rmask(j,i,sr) * flag
[2270]1609                   sums_l(k,103,tn)  = sums_l(k,103,tn)  + rad_sw_out(k,j,i)   &
[2232]1610                                       * rmask(j,i,sr) * flag
[2270]1611                   sums_l(k,104,tn)  = sums_l(k,104,tn)  + rad_lw_cs_hr(k,j,i) &
[2232]1612                                       * rmask(j,i,sr) * flag
[2270]1613                   sums_l(k,105,tn)  = sums_l(k,105,tn)  + rad_lw_hr(k,j,i)    &
[2232]1614                                       * rmask(j,i,sr) * flag
[2270]1615                   sums_l(k,106,tn)  = sums_l(k,106,tn)  + rad_sw_cs_hr(k,j,i) &
[2232]1616                                       * rmask(j,i,sr) * flag
[2270]1617                   sums_l(k,107,tn)  = sums_l(k,107,tn)  + rad_sw_hr(k,j,i)    &
[2232]1618                                       * rmask(j,i,sr) * flag
[1585]1619                ENDDO
1620             ENDDO
1621          ENDDO
1622       ENDIF
[1365]1623!
[87]1624!--    Calculate the user-defined profiles
1625       CALL user_statistics( 'profiles', sr, tn )
[1]1626       !$OMP END PARALLEL
1627
1628!
1629!--    Summation of thread sums
1630       IF ( threads_per_task > 1 )  THEN
1631          DO  i = 1, threads_per_task-1
1632             sums_l(:,3,0)          = sums_l(:,3,0) + sums_l(:,3,i)
1633             sums_l(:,4:40,0)       = sums_l(:,4:40,0) + sums_l(:,4:40,i)
[87]1634             sums_l(:,45:pr_palm,0) = sums_l(:,45:pr_palm,0) + &
1635                                      sums_l(:,45:pr_palm,i)
1636             IF ( max_pr_user > 0 )  THEN
1637                sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) = &
1638                                   sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) + &
1639                                   sums_l(:,pr_palm+1:pr_palm+max_pr_user,i)
1640             ENDIF
[1]1641          ENDDO
1642       ENDIF
1643
1644#if defined( __parallel )
[667]1645
[1]1646!
1647!--    Compute total sum from local sums
[622]1648       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
[1365]1649       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), ngp_sums, MPI_REAL,   &
[1]1650                           MPI_SUM, comm2d, ierr )
[1365]1651       IF ( large_scale_forcing )  THEN
1652          CALL MPI_ALLREDUCE( sums_ls_l(nzb,2), sums(nzb,83), ngp_sums_ls,     &
1653                              MPI_REAL, MPI_SUM, comm2d, ierr )
1654       ENDIF
[1]1655#else
1656       sums = sums_l(:,:,0)
[1365]1657       IF ( large_scale_forcing )  THEN
1658          sums(:,81:88) = sums_ls_l
1659       ENDIF
[1]1660#endif
1661
1662!
1663!--    Final values are obtained by division by the total number of grid points
1664!--    used for summation. After that store profiles.
[1738]1665!--    Check, if statistical regions do contain at least one grid point at the
1666!--    respective k-level, otherwise division by zero will lead to undefined
1667!--    values, which may cause e.g. problems with NetCDF output
[1]1668!--    Profiles:
1669       DO  k = nzb, nzt+1
[1738]1670          sums(k,3)             = sums(k,3)             / ngp_2dh(sr)
1671          sums(k,12:22)         = sums(k,12:22)         / ngp_2dh(sr)
1672          sums(k,30:32)         = sums(k,30:32)         / ngp_2dh(sr)
1673          sums(k,35:39)         = sums(k,35:39)         / ngp_2dh(sr)
1674          sums(k,45:53)         = sums(k,45:53)         / ngp_2dh(sr)
1675          sums(k,55:63)         = sums(k,55:63)         / ngp_2dh(sr)
1676          sums(k,81:88)         = sums(k,81:88)         / ngp_2dh(sr)
[2270]1677          sums(k,89:112)        = sums(k,89:112)        / ngp_2dh(sr)
1678          sums(k,114)           = sums(k,114)           / ngp_2dh(sr)
1679          sums(k,117)           = sums(k,117)           / ngp_2dh(sr)
[1738]1680          IF ( ngp_2dh_s_inner(k,sr) /= 0 )  THEN
1681             sums(k,8:11)          = sums(k,8:11)          / ngp_2dh_s_inner(k,sr)
1682             sums(k,23:29)         = sums(k,23:29)         / ngp_2dh_s_inner(k,sr)
1683             sums(k,33:34)         = sums(k,33:34)         / ngp_2dh_s_inner(k,sr)
1684             sums(k,40)            = sums(k,40)            / ngp_2dh_s_inner(k,sr)
1685             sums(k,54)            = sums(k,54)            / ngp_2dh_s_inner(k,sr)
1686             sums(k,64)            = sums(k,64)            / ngp_2dh_s_inner(k,sr)
1687             sums(k,70:80)         = sums(k,70:80)         / ngp_2dh_s_inner(k,sr)
[2270]1688             sums(k,116)           = sums(k,116)           / ngp_2dh_s_inner(k,sr)
1689             sums(k,118:pr_palm-2) = sums(k,118:pr_palm-2) / ngp_2dh_s_inner(k,sr)
[1738]1690          ENDIF
[1]1691       ENDDO
[667]1692
[1]1693!--    u* and so on
[87]1694!--    As sums(nzb:nzb+3,pr_palm) are full 2D arrays (us, usws, vsws, ts) whose
[1]1695!--    size is always ( nx + 1 ) * ( ny + 1 ), defined at the first grid layer
1696!--    above the topography, they are being divided by ngp_2dh(sr)
[87]1697       sums(nzb:nzb+3,pr_palm)    = sums(nzb:nzb+3,pr_palm)    / &
[1]1698                                    ngp_2dh(sr)
[197]1699       sums(nzb+12,pr_palm)       = sums(nzb+12,pr_palm)       / &    ! qs
1700                                    ngp_2dh(sr)
[1960]1701       sums(nzb+13,pr_palm)       = sums(nzb+13,pr_palm)       / &    ! ss
1702                                    ngp_2dh(sr)
[1]1703!--    eges, e*
[87]1704       sums(nzb+4:nzb+5,pr_palm)  = sums(nzb+4:nzb+5,pr_palm)  / &
[132]1705                                    ngp_3d(sr)
[1]1706!--    Old and new divergence
[87]1707       sums(nzb+9:nzb+10,pr_palm) = sums(nzb+9:nzb+10,pr_palm) / &
[1]1708                                    ngp_3d_inner(sr)
1709
[87]1710!--    User-defined profiles
1711       IF ( max_pr_user > 0 )  THEN
1712          DO  k = nzb, nzt+1
1713             sums(k,pr_palm+1:pr_palm+max_pr_user) = &
1714                                    sums(k,pr_palm+1:pr_palm+max_pr_user) / &
[132]1715                                    ngp_2dh_s_inner(k,sr)
[87]1716          ENDDO
1717       ENDIF
[1007]1718
[1]1719!
1720!--    Collect horizontal average in hom.
1721!--    Compute deduced averages (e.g. total heat flux)
1722       hom(:,1,3,sr)  = sums(:,3)      ! w
1723       hom(:,1,8,sr)  = sums(:,8)      ! e     profiles 5-7 are initial profiles
1724       hom(:,1,9,sr)  = sums(:,9)      ! km
1725       hom(:,1,10,sr) = sums(:,10)     ! kh
1726       hom(:,1,11,sr) = sums(:,11)     ! l
1727       hom(:,1,12,sr) = sums(:,12)     ! w"u"
1728       hom(:,1,13,sr) = sums(:,13)     ! w*u*
1729       hom(:,1,14,sr) = sums(:,14)     ! w"v"
1730       hom(:,1,15,sr) = sums(:,15)     ! w*v*
1731       hom(:,1,16,sr) = sums(:,16)     ! w"pt"
1732       hom(:,1,17,sr) = sums(:,17)     ! w*pt*
1733       hom(:,1,18,sr) = sums(:,16) + sums(:,17)    ! wpt
1734       hom(:,1,19,sr) = sums(:,12) + sums(:,13)    ! wu
1735       hom(:,1,20,sr) = sums(:,14) + sums(:,15)    ! wv
1736       hom(:,1,21,sr) = sums(:,21)     ! w*pt*BC
1737       hom(:,1,22,sr) = sums(:,16) + sums(:,21)    ! wptBC
[96]1738                                       ! profile 24 is initial profile (sa)
1739                                       ! profiles 25-29 left empty for initial
[1]1740                                       ! profiles
1741       hom(:,1,30,sr) = sums(:,30)     ! u*2
1742       hom(:,1,31,sr) = sums(:,31)     ! v*2
1743       hom(:,1,32,sr) = sums(:,32)     ! w*2
1744       hom(:,1,33,sr) = sums(:,33)     ! pt*2
1745       hom(:,1,34,sr) = sums(:,34)     ! e*
1746       hom(:,1,35,sr) = sums(:,35)     ! w*2pt*
1747       hom(:,1,36,sr) = sums(:,36)     ! w*pt*2
1748       hom(:,1,37,sr) = sums(:,37)     ! w*e*
1749       hom(:,1,38,sr) = sums(:,38)     ! w*3
[1353]1750       hom(:,1,39,sr) = sums(:,38) / ( abs( sums(:,32) ) + 1E-20_wp )**1.5_wp   ! Sw
[1]1751       hom(:,1,40,sr) = sums(:,40)     ! p
[531]1752       hom(:,1,45,sr) = sums(:,45)     ! w"vpt"
[1]1753       hom(:,1,46,sr) = sums(:,46)     ! w*vpt*       
1754       hom(:,1,47,sr) = sums(:,45) + sums(:,46)    ! wvpt
1755       hom(:,1,48,sr) = sums(:,48)     ! w"q" (w"qv")
1756       hom(:,1,49,sr) = sums(:,49)     ! w*q* (w*qv*)
1757       hom(:,1,50,sr) = sums(:,48) + sums(:,49)    ! wq (wqv)
1758       hom(:,1,51,sr) = sums(:,51)     ! w"qv"
1759       hom(:,1,52,sr) = sums(:,52)     ! w*qv*       
1760       hom(:,1,53,sr) = sums(:,52) + sums(:,51)    ! wq (wqv)
1761       hom(:,1,54,sr) = sums(:,54)     ! ql
1762       hom(:,1,55,sr) = sums(:,55)     ! w*u*u*/dz
1763       hom(:,1,56,sr) = sums(:,56)     ! w*p*/dz
[2031]1764       hom(:,1,57,sr) = sums(:,57)     ! ( w"e + w"p"/rho_ocean )/dz
[1]1765       hom(:,1,58,sr) = sums(:,58)     ! u"pt"
1766       hom(:,1,59,sr) = sums(:,59)     ! u*pt*
1767       hom(:,1,60,sr) = sums(:,58) + sums(:,59)    ! upt_t
1768       hom(:,1,61,sr) = sums(:,61)     ! v"pt"
1769       hom(:,1,62,sr) = sums(:,62)     ! v*pt*
1770       hom(:,1,63,sr) = sums(:,61) + sums(:,62)    ! vpt_t
[2031]1771       hom(:,1,64,sr) = sums(:,64)     ! rho_ocean
[96]1772       hom(:,1,65,sr) = sums(:,65)     ! w"sa"
1773       hom(:,1,66,sr) = sums(:,66)     ! w*sa*
1774       hom(:,1,67,sr) = sums(:,65) + sums(:,66)    ! wsa
[106]1775       hom(:,1,68,sr) = sums(:,68)     ! w*p*
[2031]1776       hom(:,1,69,sr) = sums(:,69)     ! w"e + w"p"/rho_ocean
[197]1777       hom(:,1,70,sr) = sums(:,70)     ! q*2
[388]1778       hom(:,1,71,sr) = sums(:,71)     ! prho
[2252]1779       hom(:,1,72,sr) = hyp * 1E-2_wp  ! hyp in hPa
[2292]1780       hom(:,1,123,sr) = sums(:,123)   ! nc
[1053]1781       hom(:,1,73,sr) = sums(:,73)     ! nr
1782       hom(:,1,74,sr) = sums(:,74)     ! qr
1783       hom(:,1,75,sr) = sums(:,75)     ! qc
1784       hom(:,1,76,sr) = sums(:,76)     ! prr (precipitation rate)
[1179]1785                                       ! 77 is initial density profile
[1241]1786       hom(:,1,78,sr) = ug             ! ug
1787       hom(:,1,79,sr) = vg             ! vg
[1299]1788       hom(:,1,80,sr) = w_subs         ! w_subs
[1]1789
[1365]1790       IF ( large_scale_forcing )  THEN
[1382]1791          hom(:,1,81,sr) = sums_ls_l(:,0)          ! td_lsa_lpt
1792          hom(:,1,82,sr) = sums_ls_l(:,1)          ! td_lsa_q
[1365]1793          IF ( use_subsidence_tendencies )  THEN
[1382]1794             hom(:,1,83,sr) = sums_ls_l(:,2)       ! td_sub_lpt
1795             hom(:,1,84,sr) = sums_ls_l(:,3)       ! td_sub_q
[1365]1796          ELSE
[1382]1797             hom(:,1,83,sr) = sums(:,83)           ! td_sub_lpt
1798             hom(:,1,84,sr) = sums(:,84)           ! td_sub_q
[1365]1799          ENDIF
[1382]1800          hom(:,1,85,sr) = sums(:,85)              ! td_nud_lpt
1801          hom(:,1,86,sr) = sums(:,86)              ! td_nud_q
1802          hom(:,1,87,sr) = sums(:,87)              ! td_nud_u
1803          hom(:,1,88,sr) = sums(:,88)              ! td_nud_v
[1365]1804       ENDIF
1805
[1551]1806       IF ( land_surface )  THEN
1807          hom(:,1,89,sr) = sums(:,89)              ! t_soil
1808                                                   ! 90 is initial t_soil profile
1809          hom(:,1,91,sr) = sums(:,91)              ! m_soil
1810                                                   ! 92 is initial m_soil profile
[2270]1811          hom(:,1,93,sr)  = sums(:,93)             ! ghf
1812          hom(:,1,94,sr)  = sums(:,94)             ! qsws_liq
1813          hom(:,1,95,sr)  = sums(:,95)             ! qsws_soil
1814          hom(:,1,96,sr)  = sums(:,96)             ! qsws_veg
1815          hom(:,1,97,sr)  = sums(:,97)             ! r_a
1816          hom(:,1,98,sr) = sums(:,98)              ! r_s
[1555]1817
[1551]1818       ENDIF
1819
1820       IF ( radiation )  THEN
[2270]1821          hom(:,1,99,sr) = sums(:,99)            ! rad_net
1822          hom(:,1,100,sr) = sums(:,100)            ! rad_lw_in
1823          hom(:,1,101,sr) = sums(:,101)            ! rad_lw_out
1824          hom(:,1,102,sr) = sums(:,102)            ! rad_sw_in
1825          hom(:,1,103,sr) = sums(:,103)            ! rad_sw_out
[1585]1826
[1691]1827          IF ( radiation_scheme == 'rrtmg' )  THEN
[2270]1828             hom(:,1,104,sr) = sums(:,104)            ! rad_lw_cs_hr
1829             hom(:,1,105,sr) = sums(:,105)            ! rad_lw_hr
1830             hom(:,1,106,sr) = sums(:,106)            ! rad_sw_cs_hr
1831             hom(:,1,107,sr) = sums(:,107)            ! rad_sw_hr
[1691]1832
[2270]1833             hom(:,1,108,sr) = sums(:,108)            ! rrtm_aldif
1834             hom(:,1,109,sr) = sums(:,109)            ! rrtm_aldir
1835             hom(:,1,110,sr) = sums(:,110)            ! rrtm_asdif
1836             hom(:,1,111,sr) = sums(:,111)            ! rrtm_asdir
[1585]1837          ENDIF
[1551]1838       ENDIF
1839
[2270]1840       hom(:,1,112,sr) = sums(:,112)            !: L
[1691]1841
[1960]1842       IF ( passive_scalar )  THEN
[2270]1843          hom(:,1,117,sr) = sums(:,117)     ! w"s"
1844          hom(:,1,114,sr) = sums(:,114)     ! w*s*
1845          hom(:,1,118,sr) = sums(:,117) + sums(:,114)    ! ws
1846          hom(:,1,116,sr) = sums(:,116)     ! s*2
[1960]1847       ENDIF
1848
[2270]1849       hom(:,1,119,sr) = rho_air       ! rho_air in Kg/m^3
1850       hom(:,1,120,sr) = rho_air_zw    ! rho_air_zw in Kg/m^3
[2037]1851
[667]1852       hom(:,1,pr_palm,sr) =   sums(:,pr_palm)
[1]1853                                       ! u*, w'u', w'v', t* (in last profile)
1854
[87]1855       IF ( max_pr_user > 0 )  THEN    ! user-defined profiles
1856          hom(:,1,pr_palm+1:pr_palm+max_pr_user,sr) = &
1857                               sums(:,pr_palm+1:pr_palm+max_pr_user)
1858       ENDIF
1859
[1]1860!
1861!--    Determine the boundary layer height using two different schemes.
[94]1862!--    First scheme: Starting from the Earth's (Ocean's) surface, look for the
1863!--    first relative minimum (maximum) of the total heat flux.
1864!--    The corresponding height is assumed as the boundary layer height, if it
1865!--    is less than 1.5 times the height where the heat flux becomes negative
1866!--    (positive) for the first time.
[1353]1867       z_i(1) = 0.0_wp
[1]1868       first = .TRUE.
[667]1869
[97]1870       IF ( ocean )  THEN
1871          DO  k = nzt, nzb+1, -1
[1738]1872             IF ( first  .AND.  hom(k,1,18,sr) < -1.0E-8_wp )  THEN
[97]1873                first = .FALSE.
1874                height = zw(k)
1875             ENDIF
[1738]1876             IF ( hom(k,1,18,sr) < -1.0E-8_wp  .AND.                           &
[97]1877                  hom(k-1,1,18,sr) > hom(k,1,18,sr) )  THEN
[1353]1878                IF ( zw(k) < 1.5_wp * height )  THEN
[97]1879                   z_i(1) = zw(k)
1880                ELSE
1881                   z_i(1) = height
1882                ENDIF
1883                EXIT
1884             ENDIF
1885          ENDDO
1886       ELSE
[94]1887          DO  k = nzb, nzt-1
[1738]1888             IF ( first  .AND.  hom(k,1,18,sr) < -1.0E-8_wp )  THEN
[94]1889                first = .FALSE.
1890                height = zw(k)
[1]1891             ENDIF
[1738]1892             IF ( hom(k,1,18,sr) < -1.0E-8_wp  .AND.                           &
[94]1893                  hom(k+1,1,18,sr) > hom(k,1,18,sr) )  THEN
[1353]1894                IF ( zw(k) < 1.5_wp * height )  THEN
[94]1895                   z_i(1) = zw(k)
1896                ELSE
1897                   z_i(1) = height
1898                ENDIF
1899                EXIT
1900             ENDIF
1901          ENDDO
[97]1902       ENDIF
[1]1903
1904!
[291]1905!--    Second scheme: Gradient scheme from Sullivan et al. (1998), modified
1906!--    by Uhlenbrock(2006). The boundary layer height is the height with the
1907!--    maximal local temperature gradient: starting from the second (the last
1908!--    but one) vertical gridpoint, the local gradient must be at least
1909!--    0.2K/100m and greater than the next four gradients.
1910!--    WARNING: The threshold value of 0.2K/100m must be adjusted for the
1911!--             ocean case!
[1353]1912       z_i(2) = 0.0_wp
[291]1913       DO  k = nzb+1, nzt+1
1914          dptdz(k) = ( hom(k,1,4,sr) - hom(k-1,1,4,sr) ) * ddzu(k)
1915       ENDDO
[1322]1916       dptdz_threshold = 0.2_wp / 100.0_wp
[291]1917
[97]1918       IF ( ocean )  THEN
[291]1919          DO  k = nzt+1, nzb+5, -1
1920             IF ( dptdz(k) > dptdz_threshold  .AND.                           &
1921                  dptdz(k) > dptdz(k-1)  .AND.  dptdz(k) > dptdz(k-2)  .AND.  &
1922                  dptdz(k) > dptdz(k-3)  .AND.  dptdz(k) > dptdz(k-4) )  THEN
1923                z_i(2) = zw(k-1)
[97]1924                EXIT
1925             ENDIF
1926          ENDDO
1927       ELSE
[291]1928          DO  k = nzb+1, nzt-3
1929             IF ( dptdz(k) > dptdz_threshold  .AND.                           &
1930                  dptdz(k) > dptdz(k+1)  .AND.  dptdz(k) > dptdz(k+2)  .AND.  &
1931                  dptdz(k) > dptdz(k+3)  .AND.  dptdz(k) > dptdz(k+4) )  THEN
1932                z_i(2) = zw(k-1)
[97]1933                EXIT
1934             ENDIF
1935          ENDDO
1936       ENDIF
[1]1937
[87]1938       hom(nzb+6,1,pr_palm,sr) = z_i(1)
1939       hom(nzb+7,1,pr_palm,sr) = z_i(2)
[1]1940
1941!
[1738]1942!--    Determine vertical index which is nearest to the mean surface level
1943!--    height of the respective statistic region
1944       DO  k = nzb, nzt
1945          IF ( zw(k) >= mean_surface_level_height(sr) )  THEN
1946             k_surface_level = k
1947             EXIT
1948          ENDIF
1949       ENDDO
1950!
[1]1951!--    Computation of both the characteristic vertical velocity and
1952!--    the characteristic convective boundary layer temperature.
[1738]1953!--    The inversion height entering into the equation is defined with respect
1954!--    to the mean surface level height of the respective statistic region.
1955!--    The horizontal average at surface level index + 1 is input for the
1956!--    average temperature.
1957       IF ( hom(k_surface_level,1,18,sr) > 1.0E-8_wp  .AND.  z_i(1) /= 0.0_wp )&
1958       THEN
[2252]1959          hom(nzb+8,1,pr_palm,sr) =                                            &
[2037]1960             ( g / hom(k_surface_level+1,1,4,sr) *                             &
[2252]1961             ( hom(k_surface_level,1,18,sr) /                                  &
1962             ( heatflux_output_conversion(nzb) * rho_air(nzb) ) )              &
[1738]1963             * ABS( z_i(1) - mean_surface_level_height(sr) ) )**0.333333333_wp
[1]1964       ELSE
[1353]1965          hom(nzb+8,1,pr_palm,sr)  = 0.0_wp
[1]1966       ENDIF
1967
[48]1968!
1969!--    Collect the time series quantities
[87]1970       ts_value(1,sr) = hom(nzb+4,1,pr_palm,sr)     ! E
1971       ts_value(2,sr) = hom(nzb+5,1,pr_palm,sr)     ! E*
[48]1972       ts_value(3,sr) = dt_3d
[87]1973       ts_value(4,sr) = hom(nzb,1,pr_palm,sr)       ! u*
1974       ts_value(5,sr) = hom(nzb+3,1,pr_palm,sr)     ! th*
[48]1975       ts_value(6,sr) = u_max
1976       ts_value(7,sr) = v_max
1977       ts_value(8,sr) = w_max
[87]1978       ts_value(9,sr) = hom(nzb+10,1,pr_palm,sr)    ! new divergence
1979       ts_value(10,sr) = hom(nzb+9,1,pr_palm,sr)    ! old Divergence
1980       ts_value(11,sr) = hom(nzb+6,1,pr_palm,sr)    ! z_i(1)
1981       ts_value(12,sr) = hom(nzb+7,1,pr_palm,sr)    ! z_i(2)
1982       ts_value(13,sr) = hom(nzb+8,1,pr_palm,sr)    ! w*
[48]1983       ts_value(14,sr) = hom(nzb,1,16,sr)           ! w'pt'   at k=0
1984       ts_value(15,sr) = hom(nzb+1,1,16,sr)         ! w'pt'   at k=1
1985       ts_value(16,sr) = hom(nzb+1,1,18,sr)         ! wpt     at k=1
1986       ts_value(17,sr) = hom(nzb,1,4,sr)            ! pt(0)
1987       ts_value(18,sr) = hom(nzb+1,1,4,sr)          ! pt(zp)
[197]1988       ts_value(19,sr) = hom(nzb+1,1,pr_palm,sr)    ! u'w'    at k=0
1989       ts_value(20,sr) = hom(nzb+2,1,pr_palm,sr)    ! v'w'    at k=0
[343]1990       ts_value(21,sr) = hom(nzb,1,48,sr)           ! w"q"    at k=0
[1709]1991
1992       IF ( .NOT. neutral )  THEN
[2270]1993          ts_value(22,sr) = hom(nzb,1,112,sr)          ! L
[48]1994       ELSE
[1709]1995          ts_value(22,sr) = 1.0E10_wp
[48]1996       ENDIF
[1]1997
[343]1998       ts_value(23,sr) = hom(nzb+12,1,pr_palm,sr)   ! q*
[1551]1999
[1960]2000       IF ( passive_scalar )  THEN
[2270]2001          ts_value(24,sr) = hom(nzb+13,1,117,sr)       ! w"s" ( to do ! )
[1960]2002          ts_value(25,sr) = hom(nzb+13,1,pr_palm,sr)   ! s*
2003       ENDIF
2004
[1]2005!
[1551]2006!--    Collect land surface model timeseries
2007       IF ( land_surface )  THEN
[2270]2008          ts_value(dots_soil  ,sr) = hom(nzb,1,93,sr)           ! ghf
2009          ts_value(dots_soil+1,sr) = hom(nzb,1,94,sr)           ! qsws_liq
2010          ts_value(dots_soil+2,sr) = hom(nzb,1,95,sr)           ! qsws_soil
2011          ts_value(dots_soil+3,sr) = hom(nzb,1,96,sr)           ! qsws_veg
2012          ts_value(dots_soil+4,sr) = hom(nzb,1,97,sr)           ! r_a
2013          ts_value(dots_soil+5,sr) = hom(nzb,1,98,sr)           ! r_s
[1551]2014       ENDIF
2015!
2016!--    Collect radiation model timeseries
2017       IF ( radiation )  THEN
[2270]2018          ts_value(dots_rad,sr)   = hom(nzb,1,99,sr)           ! rad_net
2019          ts_value(dots_rad+1,sr) = hom(nzb,1,100,sr)          ! rad_lw_in
2020          ts_value(dots_rad+2,sr) = hom(nzb,1,101,sr)          ! rad_lw_out
2021          ts_value(dots_rad+3,sr) = hom(nzb,1,102,sr)          ! rad_sw_in
2022          ts_value(dots_rad+4,sr) = hom(nzb,1,103,sr)          ! rad_sw_out
[1585]2023
2024          IF ( radiation_scheme == 'rrtmg' )  THEN
[2270]2025             ts_value(dots_rad+5,sr) = hom(nzb,1,108,sr)          ! rrtm_aldif
2026             ts_value(dots_rad+6,sr) = hom(nzb,1,109,sr)          ! rrtm_aldir
2027             ts_value(dots_rad+7,sr) = hom(nzb,1,110,sr)          ! rrtm_asdif
2028             ts_value(dots_rad+8,sr) = hom(nzb,1,111,sr)          ! rrtm_asdir
[1585]2029          ENDIF
2030
[1551]2031       ENDIF
2032
2033!
[48]2034!--    Calculate additional statistics provided by the user interface
[87]2035       CALL user_statistics( 'time_series', sr, 0 )
[1]2036
[48]2037    ENDDO    ! loop of the subregions
2038
[1]2039!
[1918]2040!-- If required, sum up horizontal averages for subsequent time averaging.
2041!-- Do not sum, if flow statistics is called before the first initial time step.
2042    IF ( do_sum  .AND.  simulated_time /= 0.0_wp )  THEN
[1353]2043       IF ( average_count_pr == 0 )  hom_sum = 0.0_wp
[1]2044       hom_sum = hom_sum + hom(:,1,:,:)
2045       average_count_pr = average_count_pr + 1
2046       do_sum = .FALSE.
2047    ENDIF
2048
2049!
2050!-- Set flag for other UPs (e.g. output routines, but also buoyancy).
2051!-- This flag is reset after each time step in time_integration.
2052    flow_statistics_called = .TRUE.
2053
2054    CALL cpu_log( log_point(10), 'flow_statistics', 'stop' )
2055
2056
2057 END SUBROUTINE flow_statistics
Note: See TracBrowser for help on using the repository browser.