[1682]  1  !> @file flow_statistics.f90 

[1036]  2  !! 

 3  ! This file is part of PALM. 

 4  ! 

 5  ! PALM is free software: you can redistribute it and/or modify it under the terms 

 6  ! of the GNU General Public License as published by the Free Software Foundation, 

 7  ! either version 3 of the License, or (at your option) any later version. 

 8  ! 

 9  ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY 

 10  ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR 

 11  ! A PARTICULAR PURPOSE. See the GNU General Public License for more details. 

 12  ! 

 13  ! You should have received a copy of the GNU General Public License along with 

 14  ! PALM. If not, see <http://www.gnu.org/licenses/>. 

 15  ! 

[1691]  16  ! Copyright 19972015 Leibniz Universitaet Hannover 

[1036]  17  !! 

 18  ! 

[254]  19  ! Current revisions: 

[1]  20  !  

[1739]  21  ! 

[1784]  22  ! 

[1739]  23  ! Former revisions: 

 24  !  

 25  ! $Id: flow_statistics.f90 1784 20160306 19:14:40Z raasch $ 

 26  ! 

[1784]  27  ! 1783 20160306 18:36:17Z raasch 

 28  ! +module netcdf_interface 

 29  ! 

[1748]  30  ! 1747 20160208 12:25:53Z raasch 

 31  ! small bugfixes for accelerator version 

 32  ! 

[1739]  33  ! 1738 20151218 13:56:05Z raasch 

[1738]  34  ! bugfixes for calculations in statistical regions which do not contain grid 

 35  ! points in the lowest vertical levels, mean surface level height considered 

 36  ! in the calculation of the characteristic vertical velocity, 

 37  ! old upstream parts removed 

[1383]  38  ! 

[1710]  39  ! 1709 20151104 14:47:01Z maronga 

 40  ! Updated output of Obukhov length 

 41  ! 

[1692]  42  ! 1691 20151026 16:17:44Z maronga 

 43  ! Revised calculation of Obukhov length. Added output of radiative heating > 

 44  ! rates for RRTMG. 

 45  ! 

[1683]  46  ! 1682 20151007 23:56:08Z knoop 

 47  ! Code annotations made doxygen readable 

 48  ! 

[1659]  49  ! 1658 20150918 10:52:53Z raasch 

 50  ! bugfix: temporary reduction variables in the openacc branch are now 

 51  ! initialized to zero 

 52  ! 

[1655]  53  ! 1654 20150917 09:20:17Z raasch 

 54  ! FORTRAN bugfix of r1652 

 55  ! 

[1653]  56  ! 1652 20150917 08:12:24Z raasch 

 57  ! bugfix in calculation of energy production by turbulent transport of TKE 

 58  ! 

[1594]  59  ! 1593 20150516 13:58:02Z raasch 

 60  ! FORTRAN errors removed from openacc branch 

 61  ! 

[1586]  62  ! 1585 20150430 07:05:52Z maronga 

 63  ! Added output of timeseries and profiles for RRTMG 

 64  ! 

[1572]  65  ! 1571 20150312 16:12:49Z maronga 

 66  ! Bugfix: output of rad_net and rad_sw_in 

 67  ! 

[1568]  68  ! 1567 20150310 17:57:55Z suehring 

 69  ! Reverse modifications made for monotonic limiter. 

 70  ! 

[1558]  71  ! 1557 20150305 16:43:04Z suehring 

 72  ! Adjustments for monotonic limiter 

 73  ! 

[1556]  74  ! 1555 20150304 17:44:27Z maronga 

 75  ! Added output of r_a and r_s. 

 76  ! 

[1552]  77  ! 1551 20150303 14:18:16Z maronga 

 78  ! Added suppport for land surface model and radiation model output. 

 79  ! 

[1499]  80  ! 1498 20141203 14:09:51Z suehring 

 81  ! Comments added 

 82  ! 

[1483]  83  ! 1482 20141018 12:34:45Z raasch 

 84  ! missing ngp_sums_ls added in accelerator version 

 85  ! 

[1451]  86  ! 1450 20140821 07:31:51Z heinze 

 87  ! bugfix: calculate fac only for simulated_time >= 0.0 

 88  ! 

[1397]  89  ! 1396 20140506 13:37:41Z raasch 

 90  ! bugfix: "copyin" replaced by "update device" in openaccbranch 

 91  ! 

[1387]  92  ! 1386 20140505 13:55:30Z boeske 

 93  ! bugfix: simulated time before the last timestep is needed for the correct 

 94  ! calculation of the profiles of large scale forcing tendencies 

 95  ! 

[1383]  96  ! 1382 20140430 12:15:41Z boeske 

[1382]  97  ! Renamed variables which store large scale forcing tendencies 

 98  ! pt_lsa > td_lsa_lpt, pt_subs > td_sub_lpt, 

 99  ! q_lsa > td_lsa_q, q_subs > td_sub_q, 

 100  ! added Neumann boundary conditions for profile data output of large scale 

 101  ! advection and subsidence terms at nzt+1 

[1354]  102  ! 

[1375]  103  ! 1374 20140425 12:55:07Z raasch 

 104  ! bugfix: syntax errors removed from openaccbranch 

 105  ! missing variables added to ONLYlists 

 106  ! 

[1366]  107  ! 1365 20140422 15:03:56Z boeske 

 108  ! Output of large scale advection, large scale subsidence and nudging tendencies 

 109  ! +sums_ls_l, ngp_sums_ls, use_subsidence_tendencies 

 110  ! 

[1354]  111  ! 1353 20140408 15:21:23Z heinze 

 112  ! REAL constants provided with KINDattribute 

 113  ! 

[1323]  114  ! 1322 20140320 16:38:49Z raasch 

 115  ! REAL constants defined as wpkind 

 116  ! 

[1321]  117  ! 1320 20140320 08:40:49Z raasch 

[1320]  118  ! ONLYattribute added to USEstatements, 

 119  ! kindparameters added to all INTEGER and REAL declaration statements, 

 120  ! kinds are defined in new module kinds, 

 121  ! revision history before 2012 removed, 

 122  ! comment fields (!:) to be used for variable explanations added to 

 123  ! all variable declaration statements 

[1008]  124  ! 

[1300]  125  ! 1299 20140306 13:15:21Z heinze 

 126  ! Output of large scale vertical velocity w_subs 

 127  ! 

[1258]  128  ! 1257 20131108 15:18:40Z raasch 

 129  ! openacc "end parallel" replaced by "end parallel loop" 

 130  ! 

[1242]  131  ! 1241 20131030 11:36:58Z heinze 

 132  ! Output of ug and vg 

 133  ! 

[1222]  134  ! 1221 20130910 08:59:13Z raasch 

 135  ! ported for openACC in separate #else branch 

 136  ! 

[1182]  137  ! 1179 20130614 05:57:58Z raasch 

 138  ! comment for profile 77 added 

 139  ! 

[1116]  140  ! 1115 20130326 18:16:16Z hoffmann 

 141  ! ql is calculated by calc_liquid_water_content 

 142  ! 

[1112]  143  ! 1111 20130308 23:54:10Z raasch 

 144  ! openACC directive added 

 145  ! 

[1054]  146  ! 1053 20121113 17:11:03Z hoffmann 

[1112]  147  ! additions for twomoment cloud physics scheme: 

[1054]  148  ! +nr, qr, qc, prr 

 149  ! 

[1037]  150  ! 1036 20121022 13:43:42Z raasch 

 151  ! code put under GPL (PALM 3.9) 

 152  ! 

[1008]  153  ! 1007 20120919 14:30:36Z franke 

[1007]  154  ! Calculation of buoyancy flux for humidity in case of WSscheme is now using 

 155  ! turbulent fluxes of WSscheme 

 156  ! Bugfix: Calculation of subgridscale buoyancy flux for humidity and cloud 

 157  ! droplets at nzb and nzt added 

[700]  158  ! 

[802]  159  ! 801 20120110 17:30:36Z suehring 

 160  ! Calculation of turbulent fluxes in advec_ws is now threadsafe. 

 161  ! 

[1]  162  ! Revision 1.1 1997/08/11 06:15:17 raasch 

 163  ! Initial revision 

 164  ! 

 165  ! 

 166  ! Description: 

 167  !  

[1682]  168  !> Compute average profiles and further average flow quantities for the different 

 169  !> userdefined (sub)regions. The region indexed 0 is the total model domain. 

 170  !> 

 171  !> @note For simplicity, nzb_s_inner and nzb_diff_s_inner are being used as a 

 172  !> lower vertical index for kloops for all variables, although strictly 

 173  !> speaking the kloops would have to be split up according to the staggered 

 174  !> grid. However, this implies no error since staggered velocity components 

 175  !> are zero at the walls and inside buildings. 

[1]  176  !! 

[1682]  177  #if ! defined( __openacc ) 

 178  SUBROUTINE flow_statistics 

 179  

[1]  180  

[1320]  181  USE arrays_3d, & 

[1691]  182  ONLY: ddzu, ddzw, e, hyp, km, kh, nr, ol, p, prho, pt, q, qc, ql, qr, & 

 183  qs, qsws, qswst, rho, sa, saswsb, saswst, shf, td_lsa_lpt, & 

 184  td_lsa_q, td_sub_lpt, td_sub_q, time_vert, ts, tswst, u, ug, us,& 

 185  usws, uswst, vsws, v, vg, vpt, vswst, w, w_subs, zw 

[1320]  186  

 187  USE cloud_parameters, & 

[1551]  188  ONLY: l_d_cp, prr, pt_d_t 

[1320]  189  

 190  USE control_parameters, & 

[1551]  191  ONLY: average_count_pr, cloud_droplets, cloud_physics, do_sum, & 

[1365]  192  dt_3d, g, humidity, icloud_scheme, kappa, large_scale_forcing, & 

[1691]  193  large_scale_subsidence, max_pr_user, message_string, neutral, & 

 194  ocean, passive_scalar, precipitation, simulated_time, & 

[1365]  195  use_subsidence_tendencies, use_surface_fluxes, use_top_fluxes, & 

 196  ws_scheme_mom, ws_scheme_sca 

[1320]  197  

 198  USE cpulog, & 

[1551]  199  ONLY: cpu_log, log_point 

[1320]  200  

 201  USE grid_variables, & 

[1551]  202  ONLY: ddx, ddy 

[1320]  203  

 204  USE indices, & 

[1551]  205  ONLY: ngp_2dh, ngp_2dh_s_inner, ngp_3d, ngp_3d_inner, ngp_sums, & 

[1365]  206  ngp_sums_ls, nxl, nxr, nyn, nys, nzb, nzb_diff_s_inner, & 

 207  nzb_s_inner, nzt, nzt_diff 

[1320]  208  

 209  USE kinds 

 210  

[1551]  211  USE land_surface_model_mod, & 

[1783]  212  ONLY: ghf_eb, land_surface, m_soil, nzb_soil, nzt_soil, & 

[1555]  213  qsws_eb, qsws_liq_eb, qsws_soil_eb, qsws_veg_eb, r_a, r_s, & 

 214  shf_eb, t_soil 

[1551]  215  

[1783]  216  USE netcdf_interface, & 

 217  ONLY: dots_rad, dots_soil 

 218  

[1]  219  USE pegrid 

[1551]  220  

 221  USE radiation_model_mod, & 

[1783]  222  ONLY: radiation, radiation_scheme, rad_net, & 

[1691]  223  rad_lw_in, rad_lw_out, rad_lw_cs_hr, rad_lw_hr, & 

 224  rad_sw_in, rad_sw_out, rad_sw_cs_hr, rad_sw_hr 

[1585]  225  

 226  #if defined ( __rrtmg ) 

 227  USE radiation_model_mod, & 

 228  ONLY: rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir 

 229  #endif 

 230  

[1]  231  USE statistics 

 232  

[1691]  233  

[1]  234  IMPLICIT NONE 

 235  

[1682]  236  INTEGER(iwp) :: i !< 

 237  INTEGER(iwp) :: j !< 

 238  INTEGER(iwp) :: k !< 

[1738]  239  INTEGER(iwp) :: k_surface_level !< 

[1682]  240  INTEGER(iwp) :: nt !< 

 241  INTEGER(iwp) :: omp_get_thread_num !< 

 242  INTEGER(iwp) :: sr !< 

 243  INTEGER(iwp) :: tn !< 

[1320]  244  

[1682]  245  LOGICAL :: first !< 

[1320]  246  

[1682]  247  REAL(wp) :: dptdz_threshold !< 

 248  REAL(wp) :: fac !< 

 249  REAL(wp) :: height !< 

 250  REAL(wp) :: pts !< 

 251  REAL(wp) :: sums_l_eper !< 

 252  REAL(wp) :: sums_l_etot !< 

 253  REAL(wp) :: ust !< 

 254  REAL(wp) :: ust2 !< 

 255  REAL(wp) :: u2 !< 

 256  REAL(wp) :: vst !< 

 257  REAL(wp) :: vst2 !< 

 258  REAL(wp) :: v2 !< 

 259  REAL(wp) :: w2 !< 

 260  REAL(wp) :: z_i(2) !< 

[1320]  261  

[1682]  262  REAL(wp) :: dptdz(nzb+1:nzt+1) !< 

 263  REAL(wp) :: sums_ll(nzb:nzt+1,2) !< 

[1]  264  

 265  CALL cpu_log( log_point(10), 'flow_statistics', 'start' ) 

 266  

[1691]  267  !$acc update host( km, kh, e, ol, pt, qs, qsws, shf, ts, u, usws, v, vsws, w ) 

[1221]  268  

[1]  269  ! 

 270  ! To be on the safe side, check whether flow_statistics has already been 

 271  ! called once after the current time step 

 272  IF ( flow_statistics_called ) THEN 

[254]  273  

[274]  274  message_string = 'flow_statistics is called two times within one ' // & 

 275  'timestep' 

[254]  276  CALL message( 'flow_statistics', 'PA0190', 1, 2, 0, 6, 0 ) 

[1007]  277  

[1]  278  ENDIF 

 279  

 280  ! 

 281  ! Compute statistics for each (sub)region 

 282  DO sr = 0, statistic_regions 

 283  

 284  ! 

 285  ! Initialize (local) summation array 

[1353]  286  sums_l = 0.0_wp 

[1]  287  

 288  ! 

 289  ! Store sums that have been computed in other subroutines in summation 

 290  ! array 

 291  sums_l(:,11,:) = sums_l_l(:,sr,:) ! mixing length from diffusivities 

 292  ! WARNING: next line still has to be adjusted for OpenMP 

 293  sums_l(:,21,0) = sums_wsts_bc_l(:,sr) ! heat flux from advec_s_bc 

[87]  294  sums_l(nzb+9,pr_palm,0) = sums_divold_l(sr) ! old divergence from pres 

 295  sums_l(nzb+10,pr_palm,0) = sums_divnew_l(sr) ! new divergence from pres 

[1]  296  

[667]  297  ! 

[1498]  298  ! When calcuating horizontallyaveraged total (resolved plus subgrid 

 299  ! scale) vertical fluxes and velocity variances by using commonly 

 300  ! applied Reynoldsbased methods ( e.g. <w'pt'> = (w<w>)*(pt<pt>) ) 

 301  ! in combination with the 5th order advection scheme, pronounced 

 302  ! artificial kinks could be observed in the vertical profiles near the 

 303  ! surface. Please note: these kinks were not related to the model truth, 

 304  ! i.e. these kinks are just related to an evaluation problem. 

 305  ! In order avoid these kinks, vertical fluxes and horizontal as well 

 306  ! vertical velocity variances are calculated directly within the advection 

 307  ! routines, according to the numerical discretization, to evaluate the 

 308  ! statistical quantities as they will appear within the prognostic 

 309  ! equations. 

[667]  310  ! Copy the turbulent quantities, evaluated in the advection routines to 

[1498]  311  ! the local array sums_l() for further computations. 

[743]  312  IF ( ws_scheme_mom .AND. sr == 0 ) THEN 

[696]  313  

[1007]  314  ! 

[673]  315  ! According to the Neumann bc for the horizontal velocity components, 

 316  ! the corresponding fluxes has to satisfiy the same bc. 

 317  IF ( ocean ) THEN 

[801]  318  sums_us2_ws_l(nzt+1,:) = sums_us2_ws_l(nzt,:) 

[1007]  319  sums_vs2_ws_l(nzt+1,:) = sums_vs2_ws_l(nzt,:) 

[673]  320  ENDIF 

[696]  321  

 322  DO i = 0, threads_per_task1 

[1007]  323  ! 

[696]  324  ! Swap the turbulent quantities evaluated in advec_ws. 

[801]  325  sums_l(:,13,i) = sums_wsus_ws_l(:,i) ! w*u* 

 326  sums_l(:,15,i) = sums_wsvs_ws_l(:,i) ! w*v* 

 327  sums_l(:,30,i) = sums_us2_ws_l(:,i) ! u*2 

 328  sums_l(:,31,i) = sums_vs2_ws_l(:,i) ! v*2 

 329  sums_l(:,32,i) = sums_ws2_ws_l(:,i) ! w*2 

[1353]  330  sums_l(:,34,i) = sums_l(:,34,i) + 0.5_wp * & 

[1320]  331  ( sums_us2_ws_l(:,i) + sums_vs2_ws_l(:,i) + & 

[801]  332  sums_ws2_ws_l(:,i) ) ! e* 

[696]  333  DO k = nzb, nzt 

[1353]  334  sums_l(nzb+5,pr_palm,i) = sums_l(nzb+5,pr_palm,i) + 0.5_wp * ( & 

[1320]  335  sums_us2_ws_l(k,i) + & 

 336  sums_vs2_ws_l(k,i) + & 

[801]  337  sums_ws2_ws_l(k,i) ) 

[696]  338  ENDDO 

[667]  339  ENDDO 

[696]  340  

[667]  341  ENDIF 

[696]  342  

[1567]  343  IF ( ws_scheme_sca .AND. sr == 0 ) THEN 

[696]  344  

 345  DO i = 0, threads_per_task1 

[801]  346  sums_l(:,17,i) = sums_wspts_ws_l(:,i) ! w*pt* from advec_s_ws 

 347  IF ( ocean ) sums_l(:,66,i) = sums_wssas_ws_l(:,i) ! w*sa* 

[696]  348  IF ( humidity .OR. passive_scalar ) sums_l(:,49,i) = & 

[801]  349  sums_wsqs_ws_l(:,i) !w*q* 

[696]  350  ENDDO 

 351  

[667]  352  ENDIF 

[305]  353  ! 

[1]  354  ! Horizontally averaged profiles of horizontal velocities and temperature. 

 355  ! They must have been computed before, because they are already required 

 356  ! for other horizontal averages. 

 357  tn = 0 

[667]  358  

[1]  359  !$OMP PARALLEL PRIVATE( i, j, k, tn ) 

[82]  360  #if defined( __intel_openmp_bug ) 

[1]  361  tn = omp_get_thread_num() 

 362  #else 

 363  !$ tn = omp_get_thread_num() 

 364  #endif 

 365  

 366  !$OMP DO 

 367  DO i = nxl, nxr 

 368  DO j = nys, nyn 

[132]  369  DO k = nzb_s_inner(j,i), nzt+1 

[1]  370  sums_l(k,1,tn) = sums_l(k,1,tn) + u(k,j,i) * rmask(j,i,sr) 

 371  sums_l(k,2,tn) = sums_l(k,2,tn) + v(k,j,i) * rmask(j,i,sr) 

 372  sums_l(k,4,tn) = sums_l(k,4,tn) + pt(k,j,i) * rmask(j,i,sr) 

 373  ENDDO 

 374  ENDDO 

 375  ENDDO 

 376  

 377  ! 

[96]  378  ! Horizontally averaged profile of salinity 

 379  IF ( ocean ) THEN 

 380  !$OMP DO 

 381  DO i = nxl, nxr 

 382  DO j = nys, nyn 

[132]  383  DO k = nzb_s_inner(j,i), nzt+1 

[96]  384  sums_l(k,23,tn) = sums_l(k,23,tn) + & 

 385  sa(k,j,i) * rmask(j,i,sr) 

 386  ENDDO 

 387  ENDDO 

 388  ENDDO 

 389  ENDIF 

 390  

 391  ! 

[1]  392  ! Horizontally averaged profiles of virtual potential temperature, 

 393  ! total water content, specific humidity and liquid water potential 

 394  ! temperature 

[75]  395  IF ( humidity ) THEN 

[1]  396  !$OMP DO 

 397  DO i = nxl, nxr 

 398  DO j = nys, nyn 

[132]  399  DO k = nzb_s_inner(j,i), nzt+1 

[1]  400  sums_l(k,44,tn) = sums_l(k,44,tn) + & 

 401  vpt(k,j,i) * rmask(j,i,sr) 

 402  sums_l(k,41,tn) = sums_l(k,41,tn) + & 

 403  q(k,j,i) * rmask(j,i,sr) 

 404  ENDDO 

 405  ENDDO 

 406  ENDDO 

 407  IF ( cloud_physics ) THEN 

 408  !$OMP DO 

 409  DO i = nxl, nxr 

 410  DO j = nys, nyn 

[132]  411  DO k = nzb_s_inner(j,i), nzt+1 

[1]  412  sums_l(k,42,tn) = sums_l(k,42,tn) + & 

 413  ( q(k,j,i)  ql(k,j,i) ) * rmask(j,i,sr) 

 414  sums_l(k,43,tn) = sums_l(k,43,tn) + ( & 

 415  pt(k,j,i) + l_d_cp*pt_d_t(k) * ql(k,j,i) & 

 416  ) * rmask(j,i,sr) 

 417  ENDDO 

 418  ENDDO 

 419  ENDDO 

 420  ENDIF 

 421  ENDIF 

 422  

 423  ! 

 424  ! Horizontally averaged profiles of passive scalar 

 425  IF ( passive_scalar ) THEN 

 426  !$OMP DO 

 427  DO i = nxl, nxr 

 428  DO j = nys, nyn 

[132]  429  DO k = nzb_s_inner(j,i), nzt+1 

[1]  430  sums_l(k,41,tn) = sums_l(k,41,tn) + q(k,j,i) * rmask(j,i,sr) 

 431  ENDDO 

 432  ENDDO 

 433  ENDDO 

 434  ENDIF 

 435  !$OMP END PARALLEL 

 436  ! 

 437  ! Summation of thread sums 

 438  IF ( threads_per_task > 1 ) THEN 

 439  DO i = 1, threads_per_task1 

 440  sums_l(:,1,0) = sums_l(:,1,0) + sums_l(:,1,i) 

 441  sums_l(:,2,0) = sums_l(:,2,0) + sums_l(:,2,i) 

 442  sums_l(:,4,0) = sums_l(:,4,0) + sums_l(:,4,i) 

[96]  443  IF ( ocean ) THEN 

 444  sums_l(:,23,0) = sums_l(:,23,0) + sums_l(:,23,i) 

 445  ENDIF 

[75]  446  IF ( humidity ) THEN 

[1]  447  sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i) 

 448  sums_l(:,44,0) = sums_l(:,44,0) + sums_l(:,44,i) 

 449  IF ( cloud_physics ) THEN 

 450  sums_l(:,42,0) = sums_l(:,42,0) + sums_l(:,42,i) 

 451  sums_l(:,43,0) = sums_l(:,43,0) + sums_l(:,43,i) 

 452  ENDIF 

 453  ENDIF 

 454  IF ( passive_scalar ) THEN 

 455  sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i) 

 456  ENDIF 

 457  ENDDO 

 458  ENDIF 

 459  

 460  #if defined( __parallel ) 

 461  ! 

 462  ! Compute total sum from local sums 

[622]  463  IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 

[1320]  464  CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2nzb, MPI_REAL, & 

[1]  465  MPI_SUM, comm2d, ierr ) 

[622]  466  IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 

[1320]  467  CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2nzb, MPI_REAL, & 

[1]  468  MPI_SUM, comm2d, ierr ) 

[622]  469  IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 

[1320]  470  CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2nzb, MPI_REAL, & 

[1]  471  MPI_SUM, comm2d, ierr ) 

[96]  472  IF ( ocean ) THEN 

[622]  473  IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 

[1320]  474  CALL MPI_ALLREDUCE( sums_l(nzb,23,0), sums(nzb,23), nzt+2nzb, & 

[96]  475  MPI_REAL, MPI_SUM, comm2d, ierr ) 

 476  ENDIF 

[75]  477  IF ( humidity ) THEN 

[622]  478  IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 

[1320]  479  CALL MPI_ALLREDUCE( sums_l(nzb,44,0), sums(nzb,44), nzt+2nzb, & 

[1]  480  MPI_REAL, MPI_SUM, comm2d, ierr ) 

[622]  481  IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 

[1320]  482  CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2nzb, & 

[1]  483  MPI_REAL, MPI_SUM, comm2d, ierr ) 

 484  IF ( cloud_physics ) THEN 

[622]  485  IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 

[1320]  486  CALL MPI_ALLREDUCE( sums_l(nzb,42,0), sums(nzb,42), nzt+2nzb, & 

[1]  487  MPI_REAL, MPI_SUM, comm2d, ierr ) 

[622]  488  IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 

[1320]  489  CALL MPI_ALLREDUCE( sums_l(nzb,43,0), sums(nzb,43), nzt+2nzb, & 

[1]  490  MPI_REAL, MPI_SUM, comm2d, ierr ) 

 491  ENDIF 

 492  ENDIF 

 493  

 494  IF ( passive_scalar ) THEN 

[622]  495  IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 

[1320]  496  CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2nzb, & 

[1]  497  MPI_REAL, MPI_SUM, comm2d, ierr ) 

 498  ENDIF 

 499  #else 

 500  sums(:,1) = sums_l(:,1,0) 

 501  sums(:,2) = sums_l(:,2,0) 

 502  sums(:,4) = sums_l(:,4,0) 

[96]  503  IF ( ocean ) sums(:,23) = sums_l(:,23,0) 

[75]  504  IF ( humidity ) THEN 

[1]  505  sums(:,44) = sums_l(:,44,0) 

 506  sums(:,41) = sums_l(:,41,0) 

 507  IF ( cloud_physics ) THEN 

 508  sums(:,42) = sums_l(:,42,0) 

 509  sums(:,43) = sums_l(:,43,0) 

 510  ENDIF 

 511  ENDIF 

 512  IF ( passive_scalar ) sums(:,41) = sums_l(:,41,0) 

 513  #endif 

 514  

 515  ! 

 516  ! Final values are obtained by division by the total number of grid points 

 517  ! used for summation. After that store profiles. 

[132]  518  sums(:,1) = sums(:,1) / ngp_2dh(sr) 

 519  sums(:,2) = sums(:,2) / ngp_2dh(sr) 

 520  sums(:,4) = sums(:,4) / ngp_2dh_s_inner(:,sr) 

[1]  521  hom(:,1,1,sr) = sums(:,1) ! u 

 522  hom(:,1,2,sr) = sums(:,2) ! v 

 523  hom(:,1,4,sr) = sums(:,4) ! pt 

 524  

[667]  525  

[1]  526  ! 

[96]  527  ! Salinity 

 528  IF ( ocean ) THEN 

[132]  529  sums(:,23) = sums(:,23) / ngp_2dh_s_inner(:,sr) 

[96]  530  hom(:,1,23,sr) = sums(:,23) ! sa 

 531  ENDIF 

 532  

 533  ! 

[1]  534  ! Humidity and cloud parameters 

[75]  535  IF ( humidity ) THEN 

[132]  536  sums(:,44) = sums(:,44) / ngp_2dh_s_inner(:,sr) 

 537  sums(:,41) = sums(:,41) / ngp_2dh_s_inner(:,sr) 

[1]  538  hom(:,1,44,sr) = sums(:,44) ! vpt 

 539  hom(:,1,41,sr) = sums(:,41) ! qv (q) 

 540  IF ( cloud_physics ) THEN 

[132]  541  sums(:,42) = sums(:,42) / ngp_2dh_s_inner(:,sr) 

 542  sums(:,43) = sums(:,43) / ngp_2dh_s_inner(:,sr) 

[1]  543  hom(:,1,42,sr) = sums(:,42) ! qv 

 544  hom(:,1,43,sr) = sums(:,43) ! pt 

 545  ENDIF 

 546  ENDIF 

 547  

 548  ! 

 549  ! Passive scalar 

[1320]  550  IF ( passive_scalar ) hom(:,1,41,sr) = sums(:,41) / & 

[132]  551  ngp_2dh_s_inner(:,sr) ! s (q) 

[1]  552  

 553  ! 

 554  ! Horizontally averaged profiles of the remaining prognostic variables, 

 555  ! variances, the total and the perturbation energy (single values in last 

 556  ! column of sums_l) and some diagnostic quantities. 

[132]  557  ! NOTE: for simplicity, nzb_s_inner is used below, although strictly 

[1]  558  !  speaking the following kloop would have to be split up and 

 559  ! rearranged according to the staggered grid. 

[132]  560  ! However, this implies no error since staggered velocity components 

 561  ! are zero at the walls and inside buildings. 

[1]  562  tn = 0 

[82]  563  #if defined( __intel_openmp_bug ) 

[1]  564  !$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper, sums_l_etot, & 

 565  !$OMP tn, ust, ust2, u2, vst, vst2, v2, w2 ) 

 566  tn = omp_get_thread_num() 

 567  #else 

 568  !$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper, sums_l_etot, tn, ust, ust2, u2, vst, vst2, v2, w2 ) 

 569  !$ tn = omp_get_thread_num() 

 570  #endif 

 571  !$OMP DO 

 572  DO i = nxl, nxr 

 573  DO j = nys, nyn 

[1353]  574  sums_l_etot = 0.0_wp 

[132]  575  DO k = nzb_s_inner(j,i), nzt+1 

[1]  576  ! 

 577  ! Prognostic and diagnostic variables 

 578  sums_l(k,3,tn) = sums_l(k,3,tn) + w(k,j,i) * rmask(j,i,sr) 

 579  sums_l(k,8,tn) = sums_l(k,8,tn) + e(k,j,i) * rmask(j,i,sr) 

 580  sums_l(k,9,tn) = sums_l(k,9,tn) + km(k,j,i) * rmask(j,i,sr) 

 581  sums_l(k,10,tn) = sums_l(k,10,tn) + kh(k,j,i) * rmask(j,i,sr) 

 582  sums_l(k,40,tn) = sums_l(k,40,tn) + p(k,j,i) 

 583  

 584  sums_l(k,33,tn) = sums_l(k,33,tn) + & 

 585  ( pt(k,j,i)hom(k,1,4,sr) )**2 * rmask(j,i,sr) 

[624]  586  

 587  IF ( humidity ) THEN 

 588  sums_l(k,70,tn) = sums_l(k,70,tn) + & 

 589  ( q(k,j,i)hom(k,1,41,sr) )**2 * rmask(j,i,sr) 

 590  ENDIF 

[1007]  591  

[699]  592  ! 

 593  ! Higher moments 

 594  ! (Computation of the skewness of w further below) 

 595  sums_l(k,38,tn) = sums_l(k,38,tn) + w(k,j,i)**3 * rmask(j,i,sr) 

[667]  596  

[1]  597  sums_l_etot = sums_l_etot + & 

[1353]  598  0.5_wp * ( u(k,j,i)**2 + v(k,j,i)**2 + & 

[667]  599  w(k,j,i)**2 ) * rmask(j,i,sr) 

[1]  600  ENDDO 

 601  ! 

 602  ! Total and perturbation energy for the total domain (being 

 603  ! collected in the last column of sums_l). Summation of these 

 604  ! quantities is seperated from the previous loop in order to 

 605  ! allow vectorization of that loop. 

[87]  606  sums_l(nzb+4,pr_palm,tn) = sums_l(nzb+4,pr_palm,tn) + sums_l_etot 

[1]  607  ! 

 608  ! 2Darrays (being collected in the last column of sums_l) 

[1320]  609  sums_l(nzb,pr_palm,tn) = sums_l(nzb,pr_palm,tn) + & 

[1]  610  us(j,i) * rmask(j,i,sr) 

[1320]  611  sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) + & 

[1]  612  usws(j,i) * rmask(j,i,sr) 

[1320]  613  sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) + & 

[1]  614  vsws(j,i) * rmask(j,i,sr) 

[1320]  615  sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) + & 

[1]  616  ts(j,i) * rmask(j,i,sr) 

[197]  617  IF ( humidity ) THEN 

[1320]  618  sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) + & 

[197]  619  qs(j,i) * rmask(j,i,sr) 

 620  ENDIF 

[1]  621  ENDDO 

 622  ENDDO 

 623  

 624  ! 

[667]  625  ! Computation of statistics when wsscheme is not used. Else these 

 626  ! quantities are evaluated in the advection routines. 

[743]  627  IF ( .NOT. ws_scheme_mom .OR. sr /= 0 ) THEN 

[667]  628  !$OMP DO 

 629  DO i = nxl, nxr 

 630  DO j = nys, nyn 

[1353]  631  sums_l_eper = 0.0_wp 

[667]  632  DO k = nzb_s_inner(j,i), nzt+1 

 633  u2 = u(k,j,i)**2 

 634  v2 = v(k,j,i)**2 

 635  w2 = w(k,j,i)**2 

 636  ust2 = ( u(k,j,i)  hom(k,1,1,sr) )**2 

 637  vst2 = ( v(k,j,i)  hom(k,1,2,sr) )**2 

 638  

 639  sums_l(k,30,tn) = sums_l(k,30,tn) + ust2 * rmask(j,i,sr) 

 640  sums_l(k,31,tn) = sums_l(k,31,tn) + vst2 * rmask(j,i,sr) 

 641  sums_l(k,32,tn) = sums_l(k,32,tn) + w2 * rmask(j,i,sr) 

 642  ! 

 643  ! Perturbation energy 

 644  

[1353]  645  sums_l(k,34,tn) = sums_l(k,34,tn) + 0.5_wp * & 

[667]  646  ( ust2 + vst2 + w2 ) * rmask(j,i,sr) 

[1353]  647  sums_l_eper = sums_l_eper + & 

 648  0.5_wp * ( ust2+vst2+w2 ) * rmask(j,i,sr) 

[667]  649  

 650  ENDDO 

[1353]  651  sums_l(nzb+5,pr_palm,tn) = sums_l(nzb+5,pr_palm,tn) & 

[667]  652  + sums_l_eper 

 653  ENDDO 

 654  ENDDO 

 655  ENDIF 

[1241]  656  

[667]  657  ! 

[1]  658  ! Horizontally averaged profiles of the vertical fluxes 

[667]  659  

[1]  660  !$OMP DO 

 661  DO i = nxl, nxr 

 662  DO j = nys, nyn 

 663  ! 

 664  ! Subgridscale fluxes (without Prandtl layer from k=nzb, 

 665  ! oterwise from k=nzb+1) 

[132]  666  ! NOTE: for simplicity, nzb_diff_s_inner is used below, although 

[1]  667  !  strictly speaking the following kloop would have to be 

 668  ! split up according to the staggered grid. 

[132]  669  ! However, this implies no error since staggered velocity 

 670  ! components are zero at the walls and inside buildings. 

 671  

 672  DO k = nzb_diff_s_inner(j,i)1, nzt_diff 

[1]  673  ! 

 674  ! Momentum flux w"u" 

[1353]  675  sums_l(k,12,tn) = sums_l(k,12,tn)  0.25_wp * ( & 

[1]  676  km(k,j,i)+km(k+1,j,i)+km(k,j,i1)+km(k+1,j,i1) & 

 677  ) * ( & 

 678  ( u(k+1,j,i)  u(k,j,i) ) * ddzu(k+1) & 

 679  + ( w(k,j,i)  w(k,j,i1) ) * ddx & 

 680  ) * rmask(j,i,sr) 

 681  ! 

 682  ! Momentum flux w"v" 

[1353]  683  sums_l(k,14,tn) = sums_l(k,14,tn)  0.25_wp * ( & 

[1]  684  km(k,j,i)+km(k+1,j,i)+km(k,j1,i)+km(k+1,j1,i) & 

 685  ) * ( & 

 686  ( v(k+1,j,i)  v(k,j,i) ) * ddzu(k+1) & 

 687  + ( w(k,j,i)  w(k,j1,i) ) * ddy & 

 688  ) * rmask(j,i,sr) 

 689  ! 

 690  ! Heat flux w"pt" 

 691  sums_l(k,16,tn) = sums_l(k,16,tn) & 

[1353]  692   0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )& 

[1]  693  * ( pt(k+1,j,i)  pt(k,j,i) ) & 

 694  * ddzu(k+1) * rmask(j,i,sr) 

 695  

 696  

 697  ! 

[96]  698  ! Salinity flux w"sa" 

 699  IF ( ocean ) THEN 

 700  sums_l(k,65,tn) = sums_l(k,65,tn) & 

[1353]  701   0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )& 

[96]  702  * ( sa(k+1,j,i)  sa(k,j,i) ) & 

 703  * ddzu(k+1) * rmask(j,i,sr) 

 704  ENDIF 

 705  

 706  ! 

[1]  707  ! Buoyancy flux, water flux (humidity flux) w"q" 

[75]  708  IF ( humidity ) THEN 

[1]  709  sums_l(k,45,tn) = sums_l(k,45,tn) & 

[1353]  710   0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )& 

[1]  711  * ( vpt(k+1,j,i)  vpt(k,j,i) ) & 

 712  * ddzu(k+1) * rmask(j,i,sr) 

 713  sums_l(k,48,tn) = sums_l(k,48,tn) & 

[1353]  714   0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )& 

[1]  715  * ( q(k+1,j,i)  q(k,j,i) ) & 

 716  * ddzu(k+1) * rmask(j,i,sr) 

[1007]  717  

[1]  718  IF ( cloud_physics ) THEN 

 719  sums_l(k,51,tn) = sums_l(k,51,tn) & 

[1353]  720   0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )& 

[1]  721  * ( ( q(k+1,j,i)  ql(k+1,j,i) )& 

 722   ( q(k,j,i)  ql(k,j,i) ) ) & 

 723  * ddzu(k+1) * rmask(j,i,sr) 

 724  ENDIF 

 725  ENDIF 

 726  

 727  ! 

 728  ! Passive scalar flux 

 729  IF ( passive_scalar ) THEN 

 730  sums_l(k,48,tn) = sums_l(k,48,tn) & 

[1353]  731   0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )& 

[1]  732  * ( q(k+1,j,i)  q(k,j,i) ) & 

 733  * ddzu(k+1) * rmask(j,i,sr) 

 734  ENDIF 

 735  

 736  ENDDO 

 737  

 738  ! 

 739  ! Subgridscale fluxes in the Prandtl layer 

 740  IF ( use_surface_fluxes ) THEN 

 741  sums_l(nzb,12,tn) = sums_l(nzb,12,tn) + & 

 742  usws(j,i) * rmask(j,i,sr) ! w"u" 

 743  sums_l(nzb,14,tn) = sums_l(nzb,14,tn) + & 

 744  vsws(j,i) * rmask(j,i,sr) ! w"v" 

 745  sums_l(nzb,16,tn) = sums_l(nzb,16,tn) + & 

 746  shf(j,i) * rmask(j,i,sr) ! w"pt" 

 747  sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + & 

[1353]  748  0.0_wp * rmask(j,i,sr) ! u"pt" 

[1]  749  sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + & 

[1353]  750  0.0_wp * rmask(j,i,sr) ! v"pt" 

[96]  751  IF ( ocean ) THEN 

 752  sums_l(nzb,65,tn) = sums_l(nzb,65,tn) + & 

 753  saswsb(j,i) * rmask(j,i,sr) ! w"sa" 

 754  ENDIF 

[75]  755  IF ( humidity ) THEN 

[1353]  756  sums_l(nzb,48,tn) = sums_l(nzb,48,tn) + & 

[1]  757  qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv") 

[1353]  758  sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + ( & 

 759  ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) * & 

 760  shf(j,i) + 0.61_wp * pt(nzb,j,i) * & 

[1007]  761  qsws(j,i) ) 

 762  IF ( cloud_droplets ) THEN 

[1353]  763  sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + ( & 

 764  ( 1.0_wp + 0.61_wp * q(nzb,j,i)  & 

 765  ql(nzb,j,i) ) * shf(j,i) + & 

 766  0.61_wp * pt(nzb,j,i) * qsws(j,i) ) 

[1007]  767  ENDIF 

[1]  768  IF ( cloud_physics ) THEN 

 769  ! 

 770  ! Formula does not work if ql(nzb) /= 0.0 

[1691]  771  sums_l(nzb,51,tn) = sums_l(nzb,51,tn) + & 

 772  qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv") 

[1]  773  ENDIF 

 774  ENDIF 

 775  IF ( passive_scalar ) THEN 

[1691]  776  sums_l(nzb,48,tn) = sums_l(nzb,48,tn) + & 

 777  qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv") 

[1]  778  ENDIF 

 779  ENDIF 

 780  

[1691]  781  IF ( .NOT. neutral ) THEN 

 782  sums_l(nzb,114,tn) = sums_l(nzb,114,tn) + & 

 783  ol(j,i) * rmask(j,i,sr) ! L 

 784  ENDIF 

 785  

 786  

[1551]  787  IF ( land_surface ) THEN 

[1555]  788  sums_l(nzb,93,tn) = sums_l(nzb,93,tn) + ghf_eb(j,i) 

 789  sums_l(nzb,94,tn) = sums_l(nzb,94,tn) + shf_eb(j,i) 

 790  sums_l(nzb,95,tn) = sums_l(nzb,95,tn) + qsws_eb(j,i) 

 791  sums_l(nzb,96,tn) = sums_l(nzb,96,tn) + qsws_liq_eb(j,i) 

 792  sums_l(nzb,97,tn) = sums_l(nzb,97,tn) + qsws_soil_eb(j,i) 

 793  sums_l(nzb,98,tn) = sums_l(nzb,98,tn) + qsws_veg_eb(j,i) 

 794  sums_l(nzb,99,tn) = sums_l(nzb,99,tn) + r_a(j,i) 

 795  sums_l(nzb,100,tn) = sums_l(nzb,100,tn)+ r_s(j,i) 

[1551]  796  ENDIF 

 797  

 798  IF ( radiation ) THEN 

[1555]  799  sums_l(nzb,101,tn) = sums_l(nzb,101,tn) + rad_net(j,i) 

[1585]  800  sums_l(nzb,102,tn) = sums_l(nzb,102,tn) + rad_lw_in(nzb,j,i) 

 801  sums_l(nzb,103,tn) = sums_l(nzb,103,tn) + rad_lw_out(nzb,j,i) 

 802  sums_l(nzb,104,tn) = sums_l(nzb,104,tn) + rad_sw_in(nzb,j,i) 

 803  sums_l(nzb,105,tn) = sums_l(nzb,105,tn) + rad_sw_out(nzb,j,i) 

 804  

 805  #if defined ( __rrtmg ) 

 806  IF ( radiation_scheme == 'rrtmg' ) THEN 

[1691]  807  sums_l(nzb,110,tn) = sums_l(nzb,110,tn) + rrtm_aldif(0,j,i) 

 808  sums_l(nzb,111,tn) = sums_l(nzb,111,tn) + rrtm_aldir(0,j,i) 

 809  sums_l(nzb,112,tn) = sums_l(nzb,112,tn) + rrtm_asdif(0,j,i) 

 810  sums_l(nzb,113,tn) = sums_l(nzb,113,tn) + rrtm_asdir(0,j,i) 

[1585]  811  ENDIF 

 812  #endif 

 813  

[1551]  814  ENDIF 

 815  

[1]  816  ! 

[19]  817  ! Subgridscale fluxes at the top surface 

 818  IF ( use_top_fluxes ) THEN 

[550]  819  sums_l(nzt:nzt+1,12,tn) = sums_l(nzt:nzt+1,12,tn) + & 

[102]  820  uswst(j,i) * rmask(j,i,sr) ! w"u" 

[550]  821  sums_l(nzt:nzt+1,14,tn) = sums_l(nzt:nzt+1,14,tn) + & 

[102]  822  vswst(j,i) * rmask(j,i,sr) ! w"v" 

[550]  823  sums_l(nzt:nzt+1,16,tn) = sums_l(nzt:nzt+1,16,tn) + & 

[19]  824  tswst(j,i) * rmask(j,i,sr) ! w"pt" 

[550]  825  sums_l(nzt:nzt+1,58,tn) = sums_l(nzt:nzt+1,58,tn) + & 

[1353]  826  0.0_wp * rmask(j,i,sr) ! u"pt" 

[550]  827  sums_l(nzt:nzt+1,61,tn) = sums_l(nzt:nzt+1,61,tn) + & 

[1353]  828  0.0_wp * rmask(j,i,sr) ! v"pt" 

[550]  829  

[96]  830  IF ( ocean ) THEN 

 831  sums_l(nzt,65,tn) = sums_l(nzt,65,tn) + & 

 832  saswst(j,i) * rmask(j,i,sr) ! w"sa" 

 833  ENDIF 

[75]  834  IF ( humidity ) THEN 

[1353]  835  sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + & 

[388]  836  qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv") 

[1353]  837  sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + ( & 

 838  ( 1.0_wp + 0.61_wp * q(nzt,j,i) ) * & 

 839  tswst(j,i) + 0.61_wp * pt(nzt,j,i) * & 

 840  qswst(j,i) ) 

[1007]  841  IF ( cloud_droplets ) THEN 

[1353]  842  sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + ( & 

 843  ( 1.0_wp + 0.61_wp * q(nzt,j,i)  & 

 844  ql(nzt,j,i) ) * tswst(j,i) + & 

 845  0.61_wp * pt(nzt,j,i) * qswst(j,i) ) 

[1007]  846  ENDIF 

[19]  847  IF ( cloud_physics ) THEN 

 848  ! 

 849  ! Formula does not work if ql(nzb) /= 0.0 

 850  sums_l(nzt,51,tn) = sums_l(nzt,51,tn) + & ! w"q" (w"qv") 

 851  qswst(j,i) * rmask(j,i,sr) 

 852  ENDIF 

 853  ENDIF 

 854  IF ( passive_scalar ) THEN 

 855  sums_l(nzt,48,tn) = sums_l(nzt,48,tn) + & 

[388]  856  qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv") 

[19]  857  ENDIF 

 858  ENDIF 

 859  

 860  ! 

[1]  861  ! Resolved fluxes (can be computed for all horizontal points) 

[132]  862  ! NOTE: for simplicity, nzb_s_inner is used below, although strictly 

[1]  863  !  speaking the following kloop would have to be split up and 

 864  ! rearranged according to the staggered grid. 

[132]  865  DO k = nzb_s_inner(j,i), nzt 

[1353]  866  ust = 0.5_wp * ( u(k,j,i)  hom(k,1,1,sr) + & 

 867  u(k+1,j,i)  hom(k+1,1,1,sr) ) 

 868  vst = 0.5_wp * ( v(k,j,i)  hom(k,1,2,sr) + & 

 869  v(k+1,j,i)  hom(k+1,1,2,sr) ) 

 870  pts = 0.5_wp * ( pt(k,j,i)  hom(k,1,4,sr) + & 

 871  pt(k+1,j,i)  hom(k+1,1,4,sr) ) 

[667]  872  

[1]  873  ! Higher moments 

[1353]  874  sums_l(k,35,tn) = sums_l(k,35,tn) + pts * w(k,j,i)**2 * & 

[1]  875  rmask(j,i,sr) 

[1353]  876  sums_l(k,36,tn) = sums_l(k,36,tn) + pts**2 * w(k,j,i) * & 

[1]  877  rmask(j,i,sr) 

 878  

 879  ! 

[96]  880  ! Salinity flux and density (density does not belong to here, 

[97]  881  ! but so far there is no other suitable place to calculate) 

[96]  882  IF ( ocean ) THEN 

[1567]  883  IF( .NOT. ws_scheme_sca .OR. sr /= 0 ) THEN 

[1353]  884  pts = 0.5_wp * ( sa(k,j,i)  hom(k,1,23,sr) + & 

 885  sa(k+1,j,i)  hom(k+1,1,23,sr) ) 

 886  sums_l(k,66,tn) = sums_l(k,66,tn) + pts * w(k,j,i) * & 

 887  rmask(j,i,sr) 

[667]  888  ENDIF 

[1353]  889  sums_l(k,64,tn) = sums_l(k,64,tn) + rho(k,j,i) * & 

[96]  890  rmask(j,i,sr) 

[1353]  891  sums_l(k,71,tn) = sums_l(k,71,tn) + prho(k,j,i) * & 

[388]  892  rmask(j,i,sr) 

[96]  893  ENDIF 

 894  

 895  ! 

[1053]  896  ! Buoyancy flux, water flux, humidity flux, liquid water 

 897  ! content, rain drop concentration and rain water content 

[75]  898  IF ( humidity ) THEN 

[1007]  899  IF ( cloud_physics .OR. cloud_droplets ) THEN 

[1353]  900  pts = 0.5_wp * ( vpt(k,j,i)  hom(k,1,44,sr) + & 

[1007]  901  vpt(k+1,j,i)  hom(k+1,1,44,sr) ) 

[1353]  902  sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * & 

[1]  903  rmask(j,i,sr) 

[1053]  904  IF ( .NOT. cloud_droplets ) THEN 

[1353]  905  pts = 0.5_wp * & 

[1115]  906  ( ( q(k,j,i)  ql(k,j,i) )  & 

 907  hom(k,1,42,sr) + & 

 908  ( q(k+1,j,i)  ql(k+1,j,i) )  & 

[1053]  909  hom(k+1,1,42,sr) ) 

[1115]  910  sums_l(k,52,tn) = sums_l(k,52,tn) + pts * w(k,j,i) * & 

[1053]  911  rmask(j,i,sr) 

 912  IF ( icloud_scheme == 0 ) THEN 

[1115]  913  sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) * & 

[1053]  914  rmask(j,i,sr) 

[1115]  915  sums_l(k,75,tn) = sums_l(k,75,tn) + qc(k,j,i) * & 

[1053]  916  rmask(j,i,sr) 

[1115]  917  IF ( precipitation ) THEN 

 918  sums_l(k,73,tn) = sums_l(k,73,tn) + nr(k,j,i) * & 

 919  rmask(j,i,sr) 

 920  sums_l(k,74,tn) = sums_l(k,74,tn) + qr(k,j,i) * & 

 921  rmask(j,i,sr) 

 922  sums_l(k,76,tn) = sums_l(k,76,tn) + prr(k,j,i) *& 

 923  rmask(j,i,sr) 

 924  ENDIF 

[1053]  925  ELSE 

[1115]  926  sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) * & 

[1053]  927  rmask(j,i,sr) 

 928  ENDIF 

 929  ELSE 

[1115]  930  sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) * & 

[1053]  931  rmask(j,i,sr) 

 932  ENDIF 

[1007]  933  ELSE 

[1567]  934  IF( .NOT. ws_scheme_sca .OR. sr /= 0 ) THEN 

[1353]  935  pts = 0.5_wp * ( vpt(k,j,i)  hom(k,1,44,sr) + & 

 936  vpt(k+1,j,i)  hom(k+1,1,44,sr) ) 

 937  sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) * & 

[1007]  938  rmask(j,i,sr) 

[1567]  939  ELSE IF ( ws_scheme_sca .AND. sr == 0 ) THEN 

[1353]  940  sums_l(k,46,tn) = ( 1.0_wp + 0.61_wp * & 

 941  hom(k,1,41,sr) ) * & 

 942  sums_l(k,17,tn) + & 

 943  0.61_wp * hom(k,1,4,sr) * & 

 944  sums_l(k,49,tn) 

[1007]  945  END IF 

 946  END IF 

[1]  947  ENDIF 

 948  ! 

 949  ! Passive scalar flux 

[1353]  950  IF ( passive_scalar .AND. ( .NOT. ws_scheme_sca & 

[1567]  951  .OR. sr /= 0 ) ) THEN 

[1353]  952  pts = 0.5_wp * ( q(k,j,i)  hom(k,1,41,sr) + & 

 953  q(k+1,j,i)  hom(k+1,1,41,sr) ) 

 954  sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) * & 

[1]  955  rmask(j,i,sr) 

 956  ENDIF 

 957  

 958  ! 

 959  ! Energy flux w*e* 

[667]  960  ! has to be adjusted 

[1353]  961  sums_l(k,37,tn) = sums_l(k,37,tn) + w(k,j,i) * 0.5_wp * & 

 962  ( ust**2 + vst**2 + w(k,j,i)**2 ) & 

[667]  963  * rmask(j,i,sr) 

[1]  964  ENDDO 

 965  ENDDO 

 966  ENDDO 

[709]  967  ! 

 968  ! For speed optimization fluxes which have been computed in part directly 

 969  ! inside the WS advection routines are treated seperatly 

 970  ! Momentum fluxes first: 

[743]  971  IF ( .NOT. ws_scheme_mom .OR. sr /= 0 ) THEN 

[667]  972  !$OMP DO 

 973  DO i = nxl, nxr 

 974  DO j = nys, nyn 

 975  DO k = nzb_diff_s_inner(j,i)1, nzt_diff 

[1353]  976  ust = 0.5_wp * ( u(k,j,i)  hom(k,1,1,sr) + & 

 977  u(k+1,j,i)  hom(k+1,1,1,sr) ) 

 978  vst = 0.5_wp * ( v(k,j,i)  hom(k,1,2,sr) + & 

 979  v(k+1,j,i)  hom(k+1,1,2,sr) ) 

[1007]  980  ! 

[667]  981  ! Momentum flux w*u* 

[1353]  982  sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5_wp * & 

 983  ( w(k,j,i1) + w(k,j,i) ) & 

[667]  984  * ust * rmask(j,i,sr) 

 985  ! 

 986  ! Momentum flux w*v* 

[1353]  987  sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5_wp * & 

 988  ( w(k,j1,i) + w(k,j,i) ) & 

[667]  989  * vst * rmask(j,i,sr) 

 990  ENDDO 

 991  ENDDO 

 992  ENDDO 

[1]  993  

[667]  994  ENDIF 

[1567]  995  IF ( .NOT. ws_scheme_sca .OR. sr /= 0 ) THEN 

[667]  996  !$OMP DO 

 997  DO i = nxl, nxr 

 998  DO j = nys, nyn 

[709]  999  DO k = nzb_diff_s_inner(j,i)1, nzt_diff 

 1000  ! 

 1001  ! Vertical heat flux 

[1353]  1002  sums_l(k,17,tn) = sums_l(k,17,tn) + 0.5_wp * & 

 1003  ( pt(k,j,i)  hom(k,1,4,sr) + & 

 1004  pt(k+1,j,i)  hom(k+1,1,4,sr) ) & 

[667]  1005  * w(k,j,i) * rmask(j,i,sr) 

 1006  IF ( humidity ) THEN 

[1353]  1007  pts = 0.5_wp * ( q(k,j,i)  hom(k,1,41,sr) + & 

 1008  q(k+1,j,i)  hom(k+1,1,41,sr) ) 

 1009  sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) * & 

 1010  rmask(j,i,sr) 

[667]  1011  ENDIF 

 1012  ENDDO 

 1013  ENDDO 

 1014  ENDDO 

 1015  

 1016  ENDIF 

 1017  

[1]  1018  ! 

[97]  1019  ! Density at top follows Neumann condition 

[388]  1020  IF ( ocean ) THEN 

 1021  sums_l(nzt+1,64,tn) = sums_l(nzt,64,tn) 

 1022  sums_l(nzt+1,71,tn) = sums_l(nzt,71,tn) 

 1023  ENDIF 

[97]  1024  

 1025  ! 

[1]  1026  ! Divergence of vertical flux of resolved scale energy and pressure 

[106]  1027  ! fluctuations as well as flux of pressure fluctuation itself (68). 

 1028  ! First calculate the products, then the divergence. 

[1]  1029  ! Calculation is time consuming. Do it only, if profiles shall be plotted. 

[1691]  1030  IF ( hom(nzb+1,2,55,0) /= 0.0_wp .OR. hom(nzb+1,2,68,0) /= 0.0_wp ) & 

 1031  THEN 

[1353]  1032  sums_ll = 0.0_wp ! local array 

[1]  1033  

 1034  !$OMP DO 

 1035  DO i = nxl, nxr 

 1036  DO j = nys, nyn 

[132]  1037  DO k = nzb_s_inner(j,i)+1, nzt 

[1]  1038  

[1353]  1039  sums_ll(k,1) = sums_ll(k,1) + 0.5_wp * w(k,j,i) * ( & 

[1652]  1040  ( 0.25_wp * ( u(k,j,i)+u(k+1,j,i)+u(k,j,i+1)+u(k+1,j,i+1) ) & 

 1041   0.5_wp * ( hom(k,1,1,sr) + hom(k+1,1,1,sr) ) )**2& 

 1042  + ( 0.25_wp * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i) ) & 

[1654]  1043   0.5_wp * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) ) )**2& 

[1353]  1044  + w(k,j,i)**2 ) 

[1]  1045  

[1353]  1046  sums_ll(k,2) = sums_ll(k,2) + 0.5_wp * w(k,j,i) & 

[1]  1047  * ( p(k,j,i) + p(k+1,j,i) ) 

 1048  

 1049  ENDDO 

 1050  ENDDO 

 1051  ENDDO 

[1353]  1052  sums_ll(0,1) = 0.0_wp ! because w is zero at the bottom 

 1053  sums_ll(nzt+1,1) = 0.0_wp 

 1054  sums_ll(0,2) = 0.0_wp 

 1055  sums_ll(nzt+1,2) = 0.0_wp 

[1]  1056  

[678]  1057  DO k = nzb+1, nzt 

[1]  1058  sums_l(k,55,tn) = ( sums_ll(k,1)  sums_ll(k1,1) ) * ddzw(k) 

 1059  sums_l(k,56,tn) = ( sums_ll(k,2)  sums_ll(k1,2) ) * ddzw(k) 

[106]  1060  sums_l(k,68,tn) = sums_ll(k,2) 

[1]  1061  ENDDO 

 1062  sums_l(nzb,55,tn) = sums_l(nzb+1,55,tn) 

 1063  sums_l(nzb,56,tn) = sums_l(nzb+1,56,tn) 

[1353]  1064  sums_l(nzb,68,tn) = 0.0_wp ! because w* = 0 at nzb 

[1]  1065  

 1066  ENDIF 

 1067  

 1068  ! 

[106]  1069  ! Divergence of vertical flux of SGS TKE and the flux itself (69) 

[1691]  1070  IF ( hom(nzb+1,2,57,0) /= 0.0_wp .OR. hom(nzb+1,2,69,0) /= 0.0_wp ) & 

 1071  THEN 

[1]  1072  !$OMP DO 

 1073  DO i = nxl, nxr 

 1074  DO j = nys, nyn 

[132]  1075  DO k = nzb_s_inner(j,i)+1, nzt 

[1]  1076  

[1353]  1077  sums_l(k,57,tn) = sums_l(k,57,tn)  0.5_wp * ( & 

[1]  1078  (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)e(k,j,i)) * ddzu(k+1) & 

 1079   (km(k1,j,i)+km(k,j,i)) * (e(k,j,i)e(k1,j,i)) * ddzu(k) & 

[1353]  1080  ) * ddzw(k) 

[1]  1081  

[1353]  1082  sums_l(k,69,tn) = sums_l(k,69,tn)  0.5_wp * ( & 

[106]  1083  (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)e(k,j,i)) * ddzu(k+1) & 

[1353]  1084  ) 

[106]  1085  

[1]  1086  ENDDO 

 1087  ENDDO 

 1088  ENDDO 

 1089  sums_l(nzb,57,tn) = sums_l(nzb+1,57,tn) 

[106]  1090  sums_l(nzb,69,tn) = sums_l(nzb+1,69,tn) 

[1]  1091  

 1092  ENDIF 

 1093  

 1094  ! 

 1095  ! Horizontal heat fluxes (subgrid, resolved, total). 

 1096  ! Do it only, if profiles shall be plotted. 

[1353]  1097  IF ( hom(nzb+1,2,58,0) /= 0.0_wp ) THEN 

[1]  1098  

 1099  !$OMP DO 

 1100  DO i = nxl, nxr 

 1101  DO j = nys, nyn 

[132]  1102  DO k = nzb_s_inner(j,i)+1, nzt 

[1]  1103  ! 

 1104  ! Subgrid horizontal heat fluxes u"pt", v"pt" 

[1353]  1105  sums_l(k,58,tn) = sums_l(k,58,tn)  0.5_wp * & 

[1]  1106  ( kh(k,j,i) + kh(k,j,i1) ) & 

 1107  * ( pt(k,j,i1)  pt(k,j,i) ) & 

 1108  * ddx * rmask(j,i,sr) 

[1353]  1109  sums_l(k,61,tn) = sums_l(k,61,tn)  0.5_wp * & 

[1]  1110  ( kh(k,j,i) + kh(k,j1,i) ) & 

 1111  * ( pt(k,j1,i)  pt(k,j,i) ) & 

 1112  * ddy * rmask(j,i,sr) 

 1113  ! 

 1114  ! Resolved horizontal heat fluxes u*pt*, v*pt* 

 1115  sums_l(k,59,tn) = sums_l(k,59,tn) + & 

 1116  ( u(k,j,i)  hom(k,1,1,sr) ) & 

[1353]  1117  * 0.5_wp * ( pt(k,j,i1)  hom(k,1,4,sr) + & 

[1]  1118  pt(k,j,i)  hom(k,1,4,sr) ) 

[1353]  1119  pts = 0.5_wp * ( pt(k,j1,i)  hom(k,1,4,sr) + & 

 1120  pt(k,j,i)  hom(k,1,4,sr) ) 

[1]  1121  sums_l(k,62,tn) = sums_l(k,62,tn) + & 

 1122  ( v(k,j,i)  hom(k,1,2,sr) ) & 

[1353]  1123  * 0.5_wp * ( pt(k,j1,i)  hom(k,1,4,sr) + & 

[1]  1124  pt(k,j,i)  hom(k,1,4,sr) ) 

 1125  ENDDO 

 1126  ENDDO 

 1127  ENDDO 

 1128  ! 

 1129  ! Fluxes at the surface must be zero (e.g. due to the Prandtllayer) 

[1353]  1130  sums_l(nzb,58,tn) = 0.0_wp 

 1131  sums_l(nzb,59,tn) = 0.0_wp 

 1132  sums_l(nzb,60,tn) = 0.0_wp 

 1133  sums_l(nzb,61,tn) = 0.0_wp 

 1134  sums_l(nzb,62,tn) = 0.0_wp 

 1135  sums_l(nzb,63,tn) = 0.0_wp 

[1]  1136  

 1137  ENDIF 

[87]  1138  

 1139  ! 

[1365]  1140  ! Collect current large scale advection and subsidence tendencies for 

 1141  ! data output 

[1691]  1142  IF ( large_scale_forcing .AND. ( simulated_time > 0.0_wp ) ) THEN 

[1365]  1143  ! 

 1144  ! Interpolation in time of LSF_DATA 

 1145  nt = 1 

[1386]  1146  DO WHILE ( simulated_time  dt_3d > time_vert(nt) ) 

[1365]  1147  nt = nt + 1 

 1148  ENDDO 

[1386]  1149  IF ( simulated_time  dt_3d /= time_vert(nt) ) THEN 

[1365]  1150  nt = nt  1 

 1151  ENDIF 

 1152  

[1386]  1153  fac = ( simulated_time  dt_3d  time_vert(nt) ) & 

[1365]  1154  / ( time_vert(nt+1)time_vert(nt) ) 

 1155  

 1156  

 1157  DO k = nzb, nzt 

[1382]  1158  sums_ls_l(k,0) = td_lsa_lpt(k,nt) & 

 1159  + fac * ( td_lsa_lpt(k,nt+1)  td_lsa_lpt(k,nt) ) 

 1160  sums_ls_l(k,1) = td_lsa_q(k,nt) & 

 1161  + fac * ( td_lsa_q(k,nt+1)  td_lsa_q(k,nt) ) 

[1365]  1162  ENDDO 

 1163  

[1382]  1164  sums_ls_l(nzt+1,0) = sums_ls_l(nzt,0) 

 1165  sums_ls_l(nzt+1,1) = sums_ls_l(nzt,1) 

 1166  

[1365]  1167  IF ( large_scale_subsidence .AND. use_subsidence_tendencies ) THEN 

 1168  

 1169  DO k = nzb, nzt 

[1382]  1170  sums_ls_l(k,2) = td_sub_lpt(k,nt) + fac * & 

 1171  ( td_sub_lpt(k,nt+1)  td_sub_lpt(k,nt) ) 

 1172  sums_ls_l(k,3) = td_sub_q(k,nt) + fac * & 

 1173  ( td_sub_q(k,nt+1)  td_sub_q(k,nt) ) 

[1365]  1174  ENDDO 

 1175  

[1382]  1176  sums_ls_l(nzt+1,2) = sums_ls_l(nzt,2) 

 1177  sums_ls_l(nzt+1,3) = sums_ls_l(nzt,3) 

 1178  

[1365]  1179  ENDIF 

 1180  

 1181  ENDIF 

 1182  

[1551]  1183  

 1184  IF ( land_surface ) THEN 

 1185  !$OMP DO 

 1186  DO i = nxl, nxr 

 1187  DO j = nys, nyn 

 1188  DO k = nzb_soil, nzt_soil 

[1691]  1189  sums_l(k,89,tn) = sums_l(k,89,tn) + t_soil(k,j,i) & 

 1190  * rmask(j,i,sr) 

 1191  sums_l(k,91,tn) = sums_l(k,91,tn) + m_soil(k,j,i) & 

 1192  * rmask(j,i,sr) 

[1551]  1193  ENDDO 

 1194  ENDDO 

 1195  ENDDO 

 1196  ENDIF 

 1197  

[1585]  1198  IF ( radiation .AND. radiation_scheme == 'rrtmg' ) THEN 

 1199  !$OMP DO 

 1200  DO i = nxl, nxr 

 1201  DO j = nys, nyn 

 1202  DO k = nzb_s_inner(j,i)+1, nzt+1 

[1691]  1203  sums_l(k,102,tn) = sums_l(k,102,tn) + rad_lw_in(k,j,i) & 

 1204  * rmask(j,i,sr) 

 1205  sums_l(k,103,tn) = sums_l(k,103,tn) + rad_lw_out(k,j,i) & 

 1206  * rmask(j,i,sr) 

 1207  sums_l(k,104,tn) = sums_l(k,104,tn) + rad_sw_in(k,j,i) & 

 1208  * rmask(j,i,sr) 

 1209  sums_l(k,105,tn) = sums_l(k,105,tn) + rad_sw_out(k,j,i) & 

 1210  * rmask(j,i,sr) 

 1211  sums_l(k,106,tn) = sums_l(k,106,tn) + rad_lw_cs_hr(k,j,i) & 

 1212  * rmask(j,i,sr) 

 1213  sums_l(k,107,tn) = sums_l(k,107,tn) + rad_lw_hr(k,j,i) & 

 1214  * rmask(j,i,sr) 

 1215  sums_l(k,108,tn) = sums_l(k,108,tn) + rad_sw_cs_hr(k,j,i) & 

 1216  * rmask(j,i,sr) 

[1701]  1217  sums_l(k,109,tn) = sums_l(k,109,tn) + rad_sw_hr(k,j,i) & 

[1691]  1218  * rmask(j,i,sr) 

[1585]  1219  ENDDO 

 1220  ENDDO 

 1221  ENDDO 

 1222  ENDIF 

[1365]  1223  ! 

[87]  1224  ! Calculate the userdefined profiles 

 1225  CALL user_statistics( 'profiles', sr, tn ) 

[1]  1226  !$OMP END PARALLEL 

 1227  

 1228  ! 

 1229  ! Summation of thread sums 

 1230  IF ( threads_per_task > 1 ) THEN 

 1231  DO i = 1, threads_per_task1 

 1232  sums_l(:,3,0) = sums_l(:,3,0) + sums_l(:,3,i) 

 1233  sums_l(:,4:40,0) = sums_l(:,4:40,0) + sums_l(:,4:40,i) 

[87]  1234  sums_l(:,45:pr_palm,0) = sums_l(:,45:pr_palm,0) + & 

 1235  sums_l(:,45:pr_palm,i) 

 1236  IF ( max_pr_user > 0 ) THEN 

 1237  sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) = & 

 1238  sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) + & 

 1239  sums_l(:,pr_palm+1:pr_palm+max_pr_user,i) 

 1240  ENDIF 

[1]  1241  ENDDO 

 1242  ENDIF 

 1243  

 1244  #if defined( __parallel ) 

[667]  1245  

[1]  1246  ! 

 1247  ! Compute total sum from local sums 

[622]  1248  IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 

[1365]  1249  CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), ngp_sums, MPI_REAL, & 

[1]  1250  MPI_SUM, comm2d, ierr ) 

[1365]  1251  IF ( large_scale_forcing ) THEN 

 1252  CALL MPI_ALLREDUCE( sums_ls_l(nzb,2), sums(nzb,83), ngp_sums_ls, & 

 1253  MPI_REAL, MPI_SUM, comm2d, ierr ) 

 1254  ENDIF 

[1]  1255  #else 

 1256  sums = sums_l(:,:,0) 

[1365]  1257  IF ( large_scale_forcing ) THEN 

 1258  sums(:,81:88) = sums_ls_l 

 1259  ENDIF 

[1]  1260  #endif 

 1261  

 1262  ! 

 1263  ! Final values are obtained by division by the total number of grid points 

 1264  ! used for summation. After that store profiles. 

[1738]  1265  ! Check, if statistical regions do contain at least one grid point at the 

 1266  ! respective klevel, otherwise division by zero will lead to undefined 

 1267  ! values, which may cause e.g. problems with NetCDF output 

[1]  1268  ! Profiles: 

 1269  DO k = nzb, nzt+1 

[1738]  1270  sums(k,3) = sums(k,3) / ngp_2dh(sr) 

 1271  sums(k,12:22) = sums(k,12:22) / ngp_2dh(sr) 

 1272  sums(k,30:32) = sums(k,30:32) / ngp_2dh(sr) 

 1273  sums(k,35:39) = sums(k,35:39) / ngp_2dh(sr) 

 1274  sums(k,45:53) = sums(k,45:53) / ngp_2dh(sr) 

 1275  sums(k,55:63) = sums(k,55:63) / ngp_2dh(sr) 

 1276  sums(k,81:88) = sums(k,81:88) / ngp_2dh(sr) 

 1277  sums(k,89:114) = sums(k,89:114) / ngp_2dh(sr) 

 1278  IF ( ngp_2dh_s_inner(k,sr) /= 0 ) THEN 

 1279  sums(k,8:11) = sums(k,8:11) / ngp_2dh_s_inner(k,sr) 

 1280  sums(k,23:29) = sums(k,23:29) / ngp_2dh_s_inner(k,sr) 

 1281  sums(k,33:34) = sums(k,33:34) / ngp_2dh_s_inner(k,sr) 

 1282  sums(k,40) = sums(k,40) / ngp_2dh_s_inner(k,sr) 

 1283  sums(k,54) = sums(k,54) / ngp_2dh_s_inner(k,sr) 

 1284  sums(k,64) = sums(k,64) / ngp_2dh_s_inner(k,sr) 

 1285  sums(k,70:80) = sums(k,70:80) / ngp_2dh_s_inner(k,sr) 

 1286  sums(k,115:pr_palm2) = sums(k,115:pr_palm2) / ngp_2dh_s_inner(k,sr) 

 1287  ENDIF 

[1]  1288  ENDDO 

[667]  1289  

[1]  1290  ! u* and so on 

[87]  1291  ! As sums(nzb:nzb+3,pr_palm) are full 2D arrays (us, usws, vsws, ts) whose 

[1]  1292  ! size is always ( nx + 1 ) * ( ny + 1 ), defined at the first grid layer 

 1293  ! above the topography, they are being divided by ngp_2dh(sr) 

[87]  1294  sums(nzb:nzb+3,pr_palm) = sums(nzb:nzb+3,pr_palm) / & 

[1]  1295  ngp_2dh(sr) 

[197]  1296  sums(nzb+12,pr_palm) = sums(nzb+12,pr_palm) / & ! qs 

 1297  ngp_2dh(sr) 

[1]  1298  ! eges, e* 

[87]  1299  sums(nzb+4:nzb+5,pr_palm) = sums(nzb+4:nzb+5,pr_palm) / & 

[132]  1300  ngp_3d(sr) 

[1]  1301  ! Old and new divergence 

[87]  1302  sums(nzb+9:nzb+10,pr_palm) = sums(nzb+9:nzb+10,pr_palm) / & 

[1]  1303  ngp_3d_inner(sr) 

 1304  

[87]  1305  ! Userdefined profiles 

 1306  IF ( max_pr_user > 0 ) THEN 

 1307  DO k = nzb, nzt+1 

 1308  sums(k,pr_palm+1:pr_palm+max_pr_user) = & 

 1309  sums(k,pr_palm+1:pr_palm+max_pr_user) / & 

[132]  1310  ngp_2dh_s_inner(k,sr) 

[87]  1311  ENDDO 

 1312  ENDIF 

[1007]  1313  

[1]  1314  ! 

 1315  ! Collect horizontal average in hom. 

 1316  ! Compute deduced averages (e.g. total heat flux) 

 1317  hom(:,1,3,sr) = sums(:,3) ! w 

 1318  hom(:,1,8,sr) = sums(:,8) ! e profiles 57 are initial profiles 

 1319  hom(:,1,9,sr) = sums(:,9) ! km 

 1320  hom(:,1,10,sr) = sums(:,10) ! kh 

 1321  hom(:,1,11,sr) = sums(:,11) ! l 

 1322  hom(:,1,12,sr) = sums(:,12) ! w"u" 

 1323  hom(:,1,13,sr) = sums(:,13) ! w*u* 

 1324  hom(:,1,14,sr) = sums(:,14) ! w"v" 

 1325  hom(:,1,15,sr) = sums(:,15) ! w*v* 

 1326  hom(:,1,16,sr) = sums(:,16) ! w"pt" 

 1327  hom(:,1,17,sr) = sums(:,17) ! w*pt* 

 1328  hom(:,1,18,sr) = sums(:,16) + sums(:,17) ! wpt 

 1329  hom(:,1,19,sr) = sums(:,12) + sums(:,13) ! wu 

 1330  hom(:,1,20,sr) = sums(:,14) + sums(:,15) ! wv 

 1331  hom(:,1,21,sr) = sums(:,21) ! w*pt*BC 

 1332  hom(:,1,22,sr) = sums(:,16) + sums(:,21) ! wptBC 

[96]  1333  ! profile 24 is initial profile (sa) 

 1334  ! profiles 2529 left empty for initial 

[1]  1335  ! profiles 

 1336  hom(:,1,30,sr) = sums(:,30) ! u*2 

 1337  hom(:,1,31,sr) = sums(:,31) ! v*2 

 1338  hom(:,1,32,sr) = sums(:,32) ! w*2 

 1339  hom(:,1,33,sr) = sums(:,33) ! pt*2 

 1340  hom(:,1,34,sr) = sums(:,34) ! e* 

 1341  hom(:,1,35,sr) = sums(:,35) ! w*2pt* 

 1342  hom(:,1,36,sr) = sums(:,36) ! w*pt*2 

 1343  hom(:,1,37,sr) = sums(:,37) ! w*e* 

 1344  hom(:,1,38,sr) = sums(:,38) ! w*3 

[1353]  1345  hom(:,1,39,sr) = sums(:,38) / ( abs( sums(:,32) ) + 1E20_wp )**1.5_wp ! Sw 

[1]  1346  hom(:,1,40,sr) = sums(:,40) ! p 

[531]  1347  hom(:,1,45,sr) = sums(:,45) ! w"vpt" 

[1]  1348  hom(:,1,46,sr) = sums(:,46) ! w*vpt* 

 1349  hom(:,1,47,sr) = sums(:,45) + sums(:,46) ! wvpt 

 1350  hom(:,1,48,sr) = sums(:,48) ! w"q" (w"qv") 

 1351  hom(:,1,49,sr) = sums(:,49) ! w*q* (w*qv*) 

 1352  hom(:,1,50,sr) = sums(:,48) + sums(:,49) ! wq (wqv) 

 1353  hom(:,1,51,sr) = sums(:,51) ! w"qv" 

 1354  hom(:,1,52,sr) = sums(:,52) ! w*qv* 

 1355  hom(:,1,53,sr) = sums(:,52) + sums(:,51) ! wq (wqv) 

 1356  hom(:,1,54,sr) = sums(:,54) ! ql 

 1357  hom(:,1,55,sr) = sums(:,55) ! w*u*u*/dz 

 1358  hom(:,1,56,sr) = sums(:,56) ! w*p*/dz 

[106]  1359  hom(:,1,57,sr) = sums(:,57) ! ( w"e + w"p"/rho )/dz 

[1]  1360  hom(:,1,58,sr) = sums(:,58) ! u"pt" 

 1361  hom(:,1,59,sr) = sums(:,59) ! u*pt* 

 1362  hom(:,1,60,sr) = sums(:,58) + sums(:,59) ! upt_t 

 1363  hom(:,1,61,sr) = sums(:,61) ! v"pt" 

 1364  hom(:,1,62,sr) = sums(:,62) ! v*pt* 

 1365  hom(:,1,63,sr) = sums(:,61) + sums(:,62) ! vpt_t 

[96]  1366  hom(:,1,64,sr) = sums(:,64) ! rho 

 1367  hom(:,1,65,sr) = sums(:,65) ! w"sa" 

 1368  hom(:,1,66,sr) = sums(:,66) ! w*sa* 

 1369  hom(:,1,67,sr) = sums(:,65) + sums(:,66) ! wsa 

[106]  1370  hom(:,1,68,sr) = sums(:,68) ! w*p* 

 1371  hom(:,1,69,sr) = sums(:,69) ! w"e + w"p"/rho 

[197]  1372  hom(:,1,70,sr) = sums(:,70) ! q*2 

[388]  1373  hom(:,1,71,sr) = sums(:,71) ! prho 

[1353]  1374  hom(:,1,72,sr) = hyp * 1E4_wp ! hyp in dbar 

[1053]  1375  hom(:,1,73,sr) = sums(:,73) ! nr 

 1376  hom(:,1,74,sr) = sums(:,74) ! qr 

 1377  hom(:,1,75,sr) = sums(:,75) ! qc 

 1378  hom(:,1,76,sr) = sums(:,76) ! prr (precipitation rate) 

[1179]  1379  ! 77 is initial density profile 

[1241]  1380  hom(:,1,78,sr) = ug ! ug 

 1381  hom(:,1,79,sr) = vg ! vg 

[1299]  1382  hom(:,1,80,sr) = w_subs ! w_subs 

[1]  1383  

[1365]  1384  IF ( large_scale_forcing ) THEN 

[1382]  1385  hom(:,1,81,sr) = sums_ls_l(:,0) ! td_lsa_lpt 

 1386  hom(:,1,82,sr) = sums_ls_l(:,1) ! td_lsa_q 

[1365]  1387  IF ( use_subsidence_tendencies ) THEN 

[1382]  1388  hom(:,1,83,sr) = sums_ls_l(:,2) ! td_sub_lpt 

 1389  hom(:,1,84,sr) = sums_ls_l(:,3) ! td_sub_q 

[1365]  1390  ELSE 

[1382]  1391  hom(:,1,83,sr) = sums(:,83) ! td_sub_lpt 

 1392  hom(:,1,84,sr) = sums(:,84) ! td_sub_q 

[1365]  1393  ENDIF 

[1382]  1394  hom(:,1,85,sr) = sums(:,85) ! td_nud_lpt 

 1395  hom(:,1,86,sr) = sums(:,86) ! td_nud_q 

 1396  hom(:,1,87,sr) = sums(:,87) ! td_nud_u 

 1397  hom(:,1,88,sr) = sums(:,88) ! td_nud_v 

[1365]  1398  ENDIF 

 1399  

[1551]  1400  IF ( land_surface ) THEN 

 1401  hom(:,1,89,sr) = sums(:,89) ! t_soil 

 1402  ! 90 is initial t_soil profile 

 1403  hom(:,1,91,sr) = sums(:,91) ! m_soil 

 1404  ! 92 is initial m_soil profile 

[1555]  1405  hom(:,1,93,sr) = sums(:,93) ! ghf_eb 

 1406  hom(:,1,94,sr) = sums(:,94) ! shf_eb 

 1407  hom(:,1,95,sr) = sums(:,95) ! qsws_eb 

 1408  hom(:,1,96,sr) = sums(:,96) ! qsws_liq_eb 

 1409  hom(:,1,97,sr) = sums(:,97) ! qsws_soil_eb 

 1410  hom(:,1,98,sr) = sums(:,98) ! qsws_veg_eb 

 1411  hom(:,1,99,sr) = sums(:,99) ! r_a 

 1412  hom(:,1,100,sr) = sums(:,100) ! r_s 

 1413  

[1551]  1414  ENDIF 

 1415  

 1416  IF ( radiation ) THEN 

[1585]  1417  hom(:,1,101,sr) = sums(:,101) ! rad_net 

 1418  hom(:,1,102,sr) = sums(:,102) ! rad_lw_in 

 1419  hom(:,1,103,sr) = sums(:,103) ! rad_lw_out 

 1420  hom(:,1,104,sr) = sums(:,104) ! rad_sw_in 

 1421  hom(:,1,105,sr) = sums(:,105) ! rad_sw_out 

 1422  

[1691]  1423  IF ( radiation_scheme == 'rrtmg' ) THEN 

[1701]  1424  #if defined ( __rrtmg ) 

[1691]  1425  hom(:,1,106,sr) = sums(:,106) ! rad_lw_cs_hr 

 1426  hom(:,1,107,sr) = sums(:,107) ! rad_lw_hr 

 1427  hom(:,1,108,sr) = sums(:,108) ! rad_sw_cs_hr 

 1428  hom(:,1,109,sr) = sums(:,109) ! rad_sw_hr 

 1429  

 1430  hom(:,1,110,sr) = sums(:,110) ! rrtm_aldif 

 1431  hom(:,1,111,sr) = sums(:,111) ! rrtm_aldir 

 1432  hom(:,1,112,sr) = sums(:,112) ! rrtm_asdif 

 1433  hom(:,1,113,sr) = sums(:,113) ! rrtm_asdir 

 1434  #endif 

[1585]  1435  ENDIF 

 1436  

[1691]  1437  

[1551]  1438  ENDIF 

 1439  

[1691]  1440  hom(:,1,114,sr) = sums(:,114) !: L 

 1441  

[667]  1442  hom(:,1,pr_palm,sr) = sums(:,pr_palm) 

[1]  1443  ! u*, w'u', w'v', t* (in last profile) 

 1444  

[87]  1445  IF ( max_pr_user > 0 ) THEN ! userdefined profiles 

 1446  hom(:,1,pr_palm+1:pr_palm+max_pr_user,sr) = & 

 1447  sums(:,pr_palm+1:pr_palm+max_pr_user) 

 1448  ENDIF 

 1449  

[1]  1450  ! 

 1451  ! Determine the boundary layer height using two different schemes. 

[94]  1452  ! First scheme: Starting from the Earth's (Ocean's) surface, look for the 

 1453  ! first relative minimum (maximum) of the total heat flux. 

 1454  ! The corresponding height is assumed as the boundary layer height, if it 

 1455  ! is less than 1.5 times the height where the heat flux becomes negative 

 1456  ! (positive) for the first time. 

[1353]  1457  z_i(1) = 0.0_wp 

[1]  1458  first = .TRUE. 

[667]  1459  

[97]  1460  IF ( ocean ) THEN 

 1461  DO k = nzt, nzb+1, 1 

[1738]  1462  IF ( first .AND. hom(k,1,18,sr) < 1.0E8_wp ) THEN 

[97]  1463  first = .FALSE. 

 1464  height = zw(k) 

 1465  ENDIF 

[1738]  1466  IF ( hom(k,1,18,sr) < 1.0E8_wp .AND. & 

[97]  1467  hom(k1,1,18,sr) > hom(k,1,18,sr) ) THEN 

[1353]  1468  IF ( zw(k) < 1.5_wp * height ) THEN 

[97]  1469  z_i(1) = zw(k) 

 1470  ELSE 

 1471  z_i(1) = height 

 1472  ENDIF 

 1473  EXIT 

 1474  ENDIF 

 1475  ENDDO 

 1476  ELSE 

[94]  1477  DO k = nzb, nzt1 

[1738]  1478  IF ( first .AND. hom(k,1,18,sr) < 1.0E8_wp ) THEN 

[94]  1479  first = .FALSE. 

 1480  height = zw(k) 

[1]  1481  ENDIF 

[1738]  1482  IF ( hom(k,1,18,sr) < 1.0E8_wp .AND. & 

[94]  1483  hom(k+1,1,18,sr) > hom(k,1,18,sr) ) THEN 

[1353]  1484  IF ( zw(k) < 1.5_wp * height ) THEN 

[94]  1485  z_i(1) = zw(k) 

 1486  ELSE 

 1487  z_i(1) = height 

 1488  ENDIF 

 1489  EXIT 

 1490  ENDIF 

 1491  ENDDO 

[97]  1492  ENDIF 

[1]  1493  

 1494  ! 

[291]  1495  ! Second scheme: Gradient scheme from Sullivan et al. (1998), modified 

 1496  ! by Uhlenbrock(2006). The boundary layer height is the height with the 

 1497  ! maximal local temperature gradient: starting from the second (the last 

 1498  ! but one) vertical gridpoint, the local gradient must be at least 

 1499  ! 0.2K/100m and greater than the next four gradients. 

 1500  ! WARNING: The threshold value of 0.2K/100m must be adjusted for the 

 1501  ! ocean case! 

[1353]  1502  z_i(2) = 0.0_wp 

[291]  1503  DO k = nzb+1, nzt+1 

 1504  dptdz(k) = ( hom(k,1,4,sr)  hom(k1,1,4,sr) ) * ddzu(k) 

 1505  ENDDO 

[1322]  1506  dptdz_threshold = 0.2_wp / 100.0_wp 

[291]  1507  

[97]  1508  IF ( ocean ) THEN 

[291]  1509  DO k = nzt+1, nzb+5, 1 

 1510  IF ( dptdz(k) > dptdz_threshold .AND. & 

 1511  dptdz(k) > dptdz(k1) .AND. dptdz(k) > dptdz(k2) .AND. & 

 1512  dptdz(k) > dptdz(k3) .AND. dptdz(k) > dptdz(k4) ) THEN 

 1513  z_i(2) = zw(k1) 

[97]  1514  EXIT 

 1515  ENDIF 

 1516  ENDDO 

 1517  ELSE 

[291]  1518  DO k = nzb+1, nzt3 

 1519  IF ( dptdz(k) > dptdz_threshold .AND. & 

 1520  dptdz(k) > dptdz(k+1) .AND. dptdz(k) > dptdz(k+2) .AND. & 

 1521  dptdz(k) > dptdz(k+3) .AND. dptdz(k) > dptdz(k+4) ) THEN 

 1522  z_i(2) = zw(k1) 

[97]  1523  EXIT 

 1524  ENDIF 

 1525  ENDDO 

 1526  ENDIF 

[1]  1527  

[87]  1528  hom(nzb+6,1,pr_palm,sr) = z_i(1) 

 1529  hom(nzb+7,1,pr_palm,sr) = z_i(2) 

[1]  1530  

 1531  ! 

[1738]  1532  ! Determine vertical index which is nearest to the mean surface level 

 1533  ! height of the respective statistic region 

 1534  DO k = nzb, nzt 

 1535  IF ( zw(k) >= mean_surface_level_height(sr) ) THEN 

 1536  k_surface_level = k 

 1537  EXIT 

 1538  ENDIF 

 1539  ENDDO 

 1540  ! 

[1]  1541  ! Computation of both the characteristic vertical velocity and 

 1542  ! the characteristic convective boundary layer temperature. 

[1738]  1543  ! The inversion height entering into the equation is defined with respect 

 1544  ! to the mean surface level height of the respective statistic region. 

 1545  ! The horizontal average at surface level index + 1 is input for the 

 1546  ! average temperature. 

 1547  IF ( hom(k_surface_level,1,18,sr) > 1.0E8_wp .AND. z_i(1) /= 0.0_wp )& 

 1548  THEN 

 1549  hom(nzb+8,1,pr_palm,sr) = & 

 1550  ( g / hom(k_surface_level+1,1,4,sr) * hom(k_surface_level,1,18,sr)& 

 1551  * ABS( z_i(1)  mean_surface_level_height(sr) ) )**0.333333333_wp 

[1]  1552  ELSE 

[1353]  1553  hom(nzb+8,1,pr_palm,sr) = 0.0_wp 

[1]  1554  ENDIF 

 1555  

[48]  1556  ! 

 1557  ! Collect the time series quantities 

[87]  1558  ts_value(1,sr) = hom(nzb+4,1,pr_palm,sr) ! E 

 1559  ts_value(2,sr) = hom(nzb+5,1,pr_palm,sr) ! E* 

[48]  1560  ts_value(3,sr) = dt_3d 

[87]  1561  ts_value(4,sr) = hom(nzb,1,pr_palm,sr) ! u* 

 1562  ts_value(5,sr) = hom(nzb+3,1,pr_palm,sr) ! th* 

[48]  1563  ts_value(6,sr) = u_max 

 1564  ts_value(7,sr) = v_max 

 1565  ts_value(8,sr) = w_max 

[87]  1566  ts_value(9,sr) = hom(nzb+10,1,pr_palm,sr) ! new divergence 

 1567  ts_value(10,sr) = hom(nzb+9,1,pr_palm,sr) ! old Divergence 

 1568  ts_value(11,sr) = hom(nzb+6,1,pr_palm,sr) ! z_i(1) 

 1569  ts_value(12,sr) = hom(nzb+7,1,pr_palm,sr) ! z_i(2) 

 1570  ts_value(13,sr) = hom(nzb+8,1,pr_palm,sr) ! w* 

[48]  1571  ts_value(14,sr) = hom(nzb,1,16,sr) ! w'pt' at k=0 

 1572  ts_value(15,sr) = hom(nzb+1,1,16,sr) ! w'pt' at k=1 

 1573  ts_value(16,sr) = hom(nzb+1,1,18,sr) ! wpt at k=1 

 1574  ts_value(17,sr) = hom(nzb,1,4,sr) ! pt(0) 

 1575  ts_value(18,sr) = hom(nzb+1,1,4,sr) ! pt(zp) 

[197]  1576  ts_value(19,sr) = hom(nzb+1,1,pr_palm,sr) ! u'w' at k=0 

 1577  ts_value(20,sr) = hom(nzb+2,1,pr_palm,sr) ! v'w' at k=0 

[343]  1578  ts_value(21,sr) = hom(nzb,1,48,sr) ! w"q" at k=0 

[1709]  1579  

 1580  IF ( .NOT. neutral ) THEN 

 1581  ts_value(22,sr) = hom(nzb,1,114,sr) ! L 

[48]  1582  ELSE 

[1709]  1583  ts_value(22,sr) = 1.0E10_wp 

[48]  1584  ENDIF 

[1]  1585  

[343]  1586  ts_value(23,sr) = hom(nzb+12,1,pr_palm,sr) ! q* 

[1551]  1587  

[1]  1588  ! 

[1551]  1589  ! Collect land surface model timeseries 

 1590  IF ( land_surface ) THEN 

 1591  ts_value(dots_soil ,sr) = hom(nzb,1,93,sr) ! ghf_eb 

 1592  ts_value(dots_soil+1,sr) = hom(nzb,1,94,sr) ! shf_eb 

 1593  ts_value(dots_soil+2,sr) = hom(nzb,1,95,sr) ! qsws_eb 

 1594  ts_value(dots_soil+3,sr) = hom(nzb,1,96,sr) ! qsws_liq_eb 

 1595  ts_value(dots_soil+4,sr) = hom(nzb,1,97,sr) ! qsws_soil_eb 

 1596  ts_value(dots_soil+5,sr) = hom(nzb,1,98,sr) ! qsws_veg_eb 

[1555]  1597  ts_value(dots_soil+6,sr) = hom(nzb,1,99,sr) ! r_a 

 1598  ts_value(dots_soil+7,sr) = hom(nzb,1,100,sr) ! r_s 

[1551]  1599  ENDIF 

 1600  ! 

 1601  ! Collect radiation model timeseries 

 1602  IF ( radiation ) THEN 

[1585]  1603  ts_value(dots_rad,sr) = hom(nzb,1,101,sr) ! rad_net 

 1604  ts_value(dots_rad+1,sr) = hom(nzb,1,102,sr) ! rad_lw_in 

 1605  ts_value(dots_rad+2,sr) = hom(nzb,1,103,sr) ! rad_lw_out 

[1701]  1606  ts_value(dots_rad+3,sr) = hom(nzb,1,104,sr) ! rad_sw_in 

 1607  ts_value(dots_rad+4,sr) = hom(nzb,1,105,sr) ! rad_sw_out 

[1585]  1608  

 1609  #if defined ( __rrtmg ) 

 1610  IF ( radiation_scheme == 'rrtmg' ) THEN 

[1691]  1611  ts_value(dots_rad+5,sr) = hom(nzb,1,110,sr) ! rrtm_aldif 

 1612  ts_value(dots_rad+6,sr) = hom(nzb,1,111,sr) ! rrtm_aldir 

 1613  ts_value(dots_rad+7,sr) = hom(nzb,1,112,sr) ! rrtm_asdif 

 1614  ts_value(dots_rad+8,sr) = hom(nzb,1,113,sr) ! rrtm_asdir 

[1585]  1615  ENDIF 

 1616  #endif 

 1617  

[1551]  1618  ENDIF 

 1619  

 1620  ! 

[48]  1621  ! Calculate additional statistics provided by the user interface 

[87]  1622  CALL user_statistics( 'time_series', sr, 0 ) 

[1]  1623  

[48]  1624  ENDDO ! loop of the subregions 

 1625  

[1]  1626  ! 

 1627  ! If required, sum up horizontal averages for subsequent time averaging 

 1628  IF ( do_sum ) THEN 

[1353]  1629  IF ( average_count_pr == 0 ) hom_sum = 0.0_wp 

[1]  1630  hom_sum = hom_sum + hom(:,1,:,:) 

 1631  average_count_pr = average_count_pr + 1 

 1632  do_sum = .FALSE. 

 1633  ENDIF 

 1634  

 1635  ! 

 1636  ! Set flag for other UPs (e.g. output routines, but also buoyancy). 

 1637  ! This flag is reset after each time step in time_integration. 

 1638  flow_statistics_called = .TRUE. 

 1639  

 1640  CALL cpu_log( log_point(10), 'flow_statistics', 'stop' ) 

 1641  

 1642  

 1643  END SUBROUTINE flow_statistics 

[1221]  1644  

 1645  

 1646  #else 

 1647  

 1648  

 1649  !! 

[1682]  1650  ! Description: 

 1651  !  

 1652  !> flow statistics  accelerator version 

[1221]  1653  !! 

 1654  SUBROUTINE flow_statistics 

 1655  

[1320]  1656  USE arrays_3d, & 

[1382]  1657  ONLY: ddzu, ddzw, e, hyp, km, kh, nr, p, prho, pt, q, qc, ql, qr, qs, & 

 1658  qsws, qswst, rho, sa, saswsb, saswst, shf, td_lsa_lpt, td_lsa_q,& 

 1659  td_sub_lpt, td_sub_q, time_vert, ts, tswst, u, ug, us, usws, & 

 1660  uswst, vsws, v, vg, vpt, vswst, w, w_subs, zw 

[1365]  1661  

[1320]  1662  

 1663  USE cloud_parameters, & 

 1664  ONLY: l_d_cp, prr, pt_d_t 

 1665  

 1666  USE control_parameters, & 

[1365]  1667  ONLY : average_count_pr, cloud_droplets, cloud_physics, do_sum, & 

 1668  dt_3d, g, humidity, icloud_scheme, kappa, large_scale_forcing, & 

[1747]  1669  large_scale_subsidence, max_pr_user, message_string, neutral, & 

 1670  ocean, passive_scalar, precipitation, simulated_time, & 

[1365]  1671  use_subsidence_tendencies, use_surface_fluxes, use_top_fluxes, & 

 1672  ws_scheme_mom, ws_scheme_sca 

[1320]  1673  

 1674  USE cpulog, & 

 1675  ONLY: cpu_log, log_point 

 1676  

 1677  USE grid_variables, & 

 1678  ONLY: ddx, ddy 

 1679  

 1680  USE indices, & 

[1482]  1681  ONLY: ngp_2dh, ngp_2dh_s_inner, ngp_3d, ngp_3d_inner, ngp_sums, & 

 1682  ngp_sums_ls, nxl, nxr, nyn, nys, nzb, nzb_diff_s_inner, & 

 1683  nzb_s_inner, nzt, nzt_diff, rflags_invers 

[1320]  1684  

 1685  USE kinds 

 1686  

[1593]  1687  USE land_surface_model_mod, & 

[1784]  1688  ONLY: ghf_eb, land_surface, m_soil, nzb_soil, nzt_soil, & 

[1593]  1689  qsws_eb, qsws_liq_eb, qsws_soil_eb, qsws_veg_eb, r_a, r_s, & 

 1690  shf_eb, t_soil 

 1691  

[1784]  1692  USE netcdf_interface, & 

 1693  ONLY: dots_rad, dots_soil 

 1694  

[1221]  1695  USE pegrid 

[1320]  1696  

[1593]  1697  USE radiation_model_mod, & 

[1784]  1698  ONLY: radiation, radiation_scheme, rad_net, & 

[1593]  1699  rad_lw_in, rad_lw_out, rad_sw_in, rad_sw_out 

 1700  

 1701  #if defined ( __rrtmg ) 

 1702  USE radiation_model_mod, & 

[1691]  1703  ONLY: rrtm_aldif, rrtm_aldir, rrtm_asdif, rrtm_asdir, rad_lw_cs_hr, & 

 1704  rad_lw_hr, rad_sw_cs_hr, rad_sw_hr 

[1593]  1705  #endif 

 1706  

[1221]  1707  USE statistics 

 1708  

 1709  IMPLICIT NONE 

 1710  

[1682]  1711  INTEGER(iwp) :: i !< 

 1712  INTEGER(iwp) :: j !< 

 1713  INTEGER(iwp) :: k !< 

[1738]  1714  INTEGER(iwp) :: k_surface_level !< 

[1682]  1715  INTEGER(iwp) :: nt !< 

 1716  INTEGER(iwp) :: omp_get_thread_num !< 

 1717  INTEGER(iwp) :: sr !< 

 1718  INTEGER(iwp) :: tn !< 

[1320]  1719  

[1682]  1720  LOGICAL :: first !< 

[1320]  1721  

[1682]  1722  REAL(wp) :: dptdz_threshold !< 

 1723  REAL(wp) :: fac !< 

 1724  REAL(wp) :: height !< 

 1725  REAL(wp) :: pts !< 

 1726  REAL(wp) :: sums_l_eper !< 

 1727  REAL(wp) :: sums_l_etot !< 

 1728  REAL(wp) :: s1 !< 

 1729  REAL(wp) :: s2 !< 

 1730  REAL(wp) :: s3 !< 

 1731  REAL(wp) :: s4 !< 

 1732  REAL(wp) :: s5 !< 

 1733  REAL(wp) :: s6 !< 

 1734  REAL(wp) :: s7 !< 

 1735  REAL(wp) :: ust !< 

 1736  REAL(wp) :: ust2 !< 

 1737  REAL(wp) :: u2 !< 

 1738  REAL(wp) :: vst !< 

 1739  REAL(wp) :: vst2 !< 

 1740  REAL(wp) :: v2 !< 

 1741  REAL(wp) :: w2 !< 

 1742  REAL(wp) :: z_i(2) !< 

[1221]  1743  

[1682]  1744  REAL(wp) :: dptdz(nzb+1:nzt+1) !< 

 1745  REAL(wp) :: sums_ll(nzb:nzt+1,2) !< 

[1320]  1746  

[1221]  1747  CALL cpu_log( log_point(10), 'flow_statistics', 'start' ) 

 1748  

 1749  ! 

 1750  ! To be on the safe side, check whether flow_statistics has already been 

 1751  ! called once after the current time step 

 1752  IF ( flow_statistics_called ) THEN 

 1753  

 1754  message_string = 'flow_statistics is called two times within one ' // & 

 1755  'timestep' 

 1756  CALL message( 'flow_statistics', 'PA0190', 1, 2, 0, 6, 0 ) 

 1757  

 1758  ENDIF 

 1759  

[1396]  1760  !$acc data create( sums, sums_l ) 

 1761  !$acc update device( hom ) 

[1221]  1762  

 1763  ! 

 1764  ! Compute statistics for each (sub)region 

 1765  DO sr = 0, statistic_regions 

 1766  

 1767  ! 

 1768  ! Initialize (local) summation array 

[1353]  1769  sums_l = 0.0_wp 

[1221]  1770  

 1771  ! 

 1772  ! Store sums that have been computed in other subroutines in summation 

 1773  ! array 

 1774  sums_l(:,11,:) = sums_l_l(:,sr,:) ! mixing length from diffusivities 

 1775  ! WARNING: next line still has to be adjusted for OpenMP 

 1776  sums_l(:,21,0) = sums_wsts_bc_l(:,sr) ! heat flux from advec_s_bc 

 1777  sums_l(nzb+9,pr_palm,0) = sums_divold_l(sr) ! old divergence from pres 

 1778  sums_l(nzb+10,pr_palm,0) = sums_divnew_l(sr) ! new divergence from pres 

 1779  

 1780  ! 

[1498]  1781  ! When calcuating horizontallyaveraged total (resolved plus subgrid 

 1782  ! scale) vertical fluxes and velocity variances by using commonly 

 1783  ! applied Reynoldsbased methods ( e.g. <w'pt'> = (w<w>)*(pt<pt>) ) 

 1784  ! in combination with the 5th order advection scheme, pronounced 

 1785  ! artificial kinks could be observed in the vertical profiles near the 

 1786  ! surface. Please note: these kinks were not related to the model truth, 

 1787  ! i.e. these kinks are just related to an evaluation problem. 

 1788  ! In order avoid these kinks, vertical fluxes and horizontal as well 

 1789  ! vertical velocity variances are calculated directly within the advection 

 1790  ! routines, according to the numerical discretization, to evaluate the 

 1791  ! statistical quantities as they will appear within the prognostic 

 1792  ! equations. 

[1221]  1793  ! Copy the turbulent quantities, evaluated in the advection routines to 

[1498]  1794  ! the local array sums_l() for further computations. 

[1221]  1795  IF ( ws_scheme_mom .AND. sr == 0 ) THEN 

 1796  

 1797  ! 

 1798  ! According to the Neumann bc for the horizontal velocity components, 

 1799  ! the corresponding fluxes has to satisfiy the same bc. 

 1800  IF ( ocean ) THEN 

 1801  sums_us2_ws_l(nzt+1,:) = sums_us2_ws_l(nzt,:) 

 1802  sums_vs2_ws_l(nzt+1,:) = sums_vs2_ws_l(nzt,:) 

 1803  ENDIF 

 1804  

 1805  DO i = 0, threads_per_task1 

 1806  ! 

 1807  ! Swap the turbulent quantities evaluated in advec_ws. 

 1808  sums_l(:,13,i) = sums_wsus_ws_l(:,i) ! w*u* 

 1809  sums_l(:,15,i) = sums_wsvs_ws_l(:,i) ! w*v* 

 1810  sums_l(:,30,i) = sums_us2_ws_l(:,i) ! u*2 

 1811  sums_l(:,31,i) = sums_vs2_ws_l(:,i) ! v*2 

 1812  sums_l(:,32,i) = sums_ws2_ws_l(:,i) ! w*2 

[1353]  1813  sums_l(:,34,i) = sums_l(:,34,i) + 0.5_wp * & 

 1814  ( sums_us2_ws_l(:,i) + sums_vs2_ws_l(:,i) + & 

[1221]  1815  sums_ws2_ws_l(:,i) ) ! e* 

 1816  DO k = nzb, nzt 

[1353]  1817  sums_l(nzb+5,pr_palm,i) = sums_l(nzb+5,pr_palm,i) + 0.5_wp * ( & 

 1818  sums_us2_ws_l(k,i) + & 

 1819  sums_vs2_ws_l(k,i) + & 

 1820  sums_ws2_ws_l(k,i) ) 

[1221]  1821  ENDDO 

 1822  ENDDO 

 1823  

 1824  ENDIF 

 1825  

[1567]  1826  IF ( ws_scheme_sca .AND. sr == 0 ) THEN 

[1221]  1827  

 1828  DO i = 0, threads_per_task1 

 1829  sums_l(:,17,i) = sums_wspts_ws_l(:,i) ! w*pt* from advec_s_ws 

 1830  IF ( ocean ) sums_l(:,66,i) = sums_wssas_ws_l(:,i) ! w*sa* 

 1831  IF ( humidity .OR. passive_scalar ) sums_l(:,49,i) = & 

 1832  sums_wsqs_ws_l(:,i) !w*q* 

 1833  ENDDO 

 1834  

 1835  ENDIF 

 1836  ! 

 1837  ! Horizontally averaged profiles of horizontal velocities and temperature. 

 1838  ! They must have been computed before, because they are already required 

 1839  ! for other horizontal averages. 

 1840  tn = 0 

 1841  

 1842  !$OMP PARALLEL PRIVATE( i, j, k, tn ) 

 1843  #if defined( __intel_openmp_bug ) 

 1844  tn = omp_get_thread_num() 

 1845  #else 

 1846  !$ tn = omp_get_thread_num() 

 1847  #endif 

 1848  

 1849  !$acc update device( sums_l ) 

 1850  

 1851  !$OMP DO 

 1852  !$acc parallel loop gang present( pt, rflags_invers, rmask, sums_l, u, v ) create( s1, s2, s3 ) 

 1853  DO k = nzb, nzt+1 

[1658]  1854  s1 = 0 

 1855  s2 = 0 

 1856  s3 = 0 

[1221]  1857  !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 ) 

 1858  DO i = nxl, nxr 

 1859  DO j = nys, nyn 

 1860  ! 

 1861  ! k+1 is used in rflags since rflags is set 0 at surface points 

 1862  s1 = s1 + u(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1) 

 1863  s2 = s2 + v(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1) 

 1864  s3 = s3 + pt(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1) 

 1865  ENDDO 

 1866  ENDDO 

 1867  sums_l(k,1,tn) = s1 

 1868  sums_l(k,2,tn) = s2 

 1869  sums_l(k,4,tn) = s3 

 1870  ENDDO 

[1257]  1871  !$acc end parallel loop 

[1221]  1872  

 1873  ! 

 1874  ! Horizontally averaged profile of salinity 

 1875  IF ( ocean ) THEN 

 1876  !$OMP DO 

 1877  !$acc parallel loop gang present( rflags_invers, rmask, sums_l, sa ) create( s1 ) 

 1878  DO k = nzb, nzt+1 

[1658]  1879  s1 = 0 

[1221]  1880  !$acc loop vector collapse( 2 ) reduction( +: s1 ) 

 1881  DO i = nxl, nxr 

 1882  DO j = nys, nyn 

 1883  s1 = s1 + sa(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1) 

 1884  ENDDO 

 1885  ENDDO 

 1886  sums_l(k,23,tn) = s1 

 1887  ENDDO 

[1257]  1888  !$acc end parallel loop 

[1221]  1889  ENDIF 

 1890  

 1891  ! 

 1892  ! Horizontally averaged profiles of virtual potential temperature, 

 1893  ! total water content, specific humidity and liquid water potential 

 1894  ! temperature 

 1895  IF ( humidity ) THEN 

 1896  

 1897  !$OMP DO 

 1898  !$acc parallel loop gang present( q, rflags_invers, rmask, sums_l, vpt ) create( s1, s2 ) 

 1899  DO k = nzb, nzt+1 

[1658]  1900  s1 = 0 

 1901  s2 = 0 

[1221]  1902  !$acc loop vector collapse( 2 ) reduction( +: s1, s2 ) 

 1903  DO i = nxl, nxr 

 1904  DO j = nys, nyn 

 1905  s1 = s1 + q(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1) 

 1906  s2 = s2 + vpt(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1) 

 1907  ENDDO 

 1908  ENDDO 

 1909  sums_l(k,41,tn) = s1 

 1910  sums_l(k,44,tn) = s2 

 1911  ENDDO 

[1257]  1912  !$acc end parallel loop 

[1221]  1913  

 1914  IF ( cloud_physics ) THEN 

 1915  !$OMP DO 

 1916  !$acc parallel loop gang present( pt, q, ql, rflags_invers, rmask, sums_l ) create( s1, s2 ) 

 1917  DO k = nzb, nzt+1 

[1658]  1918  s1 = 0 

 1919  s2 = 0 

[1221]  1920  !$acc loop vector collapse( 2 ) reduction( +: s1, s2 ) 

 1921  DO i = nxl, nxr 

 1922  DO j = nys, nyn 

 1923  s1 = s1 + ( q(k,j,i)  ql(k,j,i) ) * & 

 1924  rmask(j,i,sr) * rflags_invers(j,i,k+1) 

 1925  s2 = s2 + ( pt(k,j,i) + l_d_cp*pt_d_t(k) * ql(k,j,i) ) * & 

 1926  rmask(j,i,sr) * rflags_invers(j,i,k+1) 

 1927  ENDDO 

 1928  ENDDO 

 1929  sums_l(k,42,tn) = s1 

 1930  sums_l(k,43,tn) = s2 

 1931  ENDDO 

[1257]  1932  !$acc end parallel loop 

[1221]  1933  ENDIF 

 1934  ENDIF 

 1935  

 1936  ! 

 1937  ! Horizontally averaged profiles of passive scalar 

 1938  IF ( passive_scalar ) THEN 

 1939  !$OMP DO 

 1940  !$acc parallel loop gang present( q, rflags_invers, rmask, sums_l ) create( s1 ) 

 1941  DO k = nzb, nzt+1 

[1658]  1942  s1 = 0 

[1221]  1943  !$acc loop vector collapse( 2 ) reduction( +: s1 ) 

 1944  DO i = nxl, nxr 

 1945  DO j = nys, nyn 

 1946  s1 = s1 + q(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1) 

 1947  ENDDO 

 1948  ENDDO 

 1949  sums_l(k,41,tn) = s1 

 1950  ENDDO 

[1257]  1951  !$acc end parallel loop 

[1221]  1952  ENDIF 

 1953  !$OMP END PARALLEL 

 1954  

 1955  ! 

 1956  ! Summation of thread sums 

 1957  IF ( threads_per_task > 1 ) THEN 

 1958  DO i = 1, threads_per_task1 

 1959  !$acc parallel present( sums_l ) 

 1960  sums_l(:,1,0) = sums_l(:,1,0) + sums_l(:,1,i) 

 1961  sums_l(:,2,0) = sums_l(:,2,0) + sums_l(:,2,i) 

 1962  sums_l(:,4,0) = sums_l(:,4,0) + sums_l(:,4,i) 

 1963  !$acc end parallel 

 1964  IF ( ocean ) THEN 

 1965  !$acc parallel present( sums_l ) 

 1966  sums_l(:,23,0) = sums_l(:,23,0) + sums_l(:,23,i) 

 1967  !$acc end parallel 

 1968  ENDIF 

 1969  IF ( humidity ) THEN 

 1970  !$acc parallel present( sums_l ) 

 1971  sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i) 

 1972  sums_l(:,44,0) = sums_l(:,44,0) + sums_l(:,44,i) 

 1973  !$acc end parallel 

 1974  IF ( cloud_physics ) THEN 

 1975  !$acc parallel present( sums_l ) 

 1976  sums_l(:,42,0) = sums_l(:,42,0) + sums_l(:,42,i) 

 1977  sums_l(:,43,0) = sums_l(:,43,0) + sums_l(:,43,i) 

 1978  !$acc end parallel 

 1979  ENDIF 

 1980  ENDIF 

 1981  IF ( passive_scalar ) THEN 

 1982  !$acc parallel present( sums_l ) 

 1983  sums_l(:,41,0) = sums_l(:,41,0) + sums_l(:,41,i) 

 1984  !$acc end parallel 

 1985  ENDIF 

 1986  ENDDO 

 1987  ENDIF 

 1988  

 1989  #if defined( __parallel ) 

 1990  ! 

 1991  ! Compute total sum from local sums 

 1992  !$acc update host( sums_l ) 

 1993  IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 

[1353]  1994  CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2nzb, MPI_REAL, & 

[1221]  1995  MPI_SUM, comm2d, ierr ) 

 1996  IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 

[1353]  1997  CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2nzb, MPI_REAL, & 

[1221]  1998  MPI_SUM, comm2d, ierr ) 

 1999  IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 

[1353]  2000  CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2nzb, MPI_REAL, & 

[1221]  2001  MPI_SUM, comm2d, ierr ) 

 2002  IF ( ocean ) THEN 

 2003  IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 

[1353]  2004  CALL MPI_ALLREDUCE( sums_l(nzb,23,0), sums(nzb,23), nzt+2nzb, & 

[1221]  2005  MPI_REAL, MPI_SUM, comm2d, ierr ) 

 2006  ENDIF 

 2007  IF ( humidity ) THEN 

 2008  IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 

[1353]  2009  CALL MPI_ALLREDUCE( sums_l(nzb,44,0), sums(nzb,44), nzt+2nzb, & 

[1221]  2010  MPI_REAL, MPI_SUM, comm2d, ierr ) 

 2011  IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 

[1353]  2012  CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2nzb, & 

[1221]  2013  MPI_REAL, MPI_SUM, comm2d, ierr ) 

 2014  IF ( cloud_physics ) THEN 

 2015  IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 

[1353]  2016  CALL MPI_ALLREDUCE( sums_l(nzb,42,0), sums(nzb,42), nzt+2nzb, & 

[1221]  2017  MPI_REAL, MPI_SUM, comm2d, ierr ) 

 2018  IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 

[1353]  2019  CALL MPI_ALLREDUCE( sums_l(nzb,43,0), sums(nzb,43), nzt+2nzb, & 

[1221]  2020  MPI_REAL, MPI_SUM, comm2d, ierr ) 

 2021  ENDIF 

 2022  ENDIF 

 2023  

 2024  IF ( passive_scalar ) THEN 

 2025  IF ( collective_wait ) CALL MPI_BARRIER( comm2d, ierr ) 

[1353]  2026  CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2nzb, & 

[1221]  2027  MPI_REAL, MPI_SUM, comm2d, ierr ) 

 2028  ENDIF 

 2029  !$acc update device( sums ) 

 2030  #else 

 2031  !$acc parallel present( sums, sums_l ) 

 2032  sums(:,1) = sums_l(:,1,0) 

 2033  sums(:,2) = sums_l(:,2,0) 

 2034  sums(:,4) = sums_l(:,4,0) 

 2035  !$acc end parallel 

 2036  IF ( ocean ) THEN 

 2037  !$acc parallel present( sums, sums_l ) 

 2038  sums(:,23) = sums_l(:,23,0) 

 2039  !$acc end parallel 

 2040  ENDIF 

 2041  IF ( humidity ) THEN 

 2042  !$acc parallel present( sums, sums_l ) 

 2043  sums(:,44) = sums_l(:,44,0) 

 2044  sums(:,41) = sums_l(:,41,0) 

 2045  !$acc end parallel 

 2046  IF ( cloud_physics ) THEN 

 2047  !$acc parallel present( sums, sums_l ) 

 2048  sums(:,42) = sums_l(:,42,0) 

 2049  sums(:,43) = sums_l(:,43,0) 

 2050  !$acc end parallel 

 2051  ENDIF 

 2052  ENDIF 

 2053  IF ( passive_scalar ) THEN 

 2054  !$acc parallel present( sums, sums_l ) 

 2055  sums(:,41) = sums_l(:,41,0) 

 2056  !$acc end parallel 

 2057  ENDIF 

 2058  #endif 

 2059  

 2060  ! 

 2061  ! Final values are obtained by division by the total number of grid points 

 2062  ! used for summation. After that store profiles. 

 2063  !$acc parallel present( hom, ngp_2dh, ngp_2dh_s_inner, sums ) 

 2064  sums(:,1) = sums(:,1) / ngp_2dh(sr) 

 2065  sums(:,2) = sums(:,2) / ngp_2dh(sr) 

 2066  sums(:,4) = sums(:,4) / ngp_2dh_s_inner(:,sr) 

 2067  hom(:,1,1,sr) = sums(:,1) ! u 

 2068  hom(:,1,2,sr) = sums(:,2) ! v 

 2069  hom(:,1,4,sr) = sums(:,4) ! pt 

 2070  !$acc end parallel 

 2071  

 2072  ! 

 2073  ! Salinity 

 2074  IF ( ocean ) THEN 

 2075  !$acc parallel present( hom, ngp_2dh_s_inner, sums ) 

 2076  sums(:,23) = sums(:,23) / ngp_2dh_s_inner(:,sr) 

 2077  hom(:,1,23,sr) = sums(:,23) ! sa 

 2078  !$acc end parallel 

 2079  ENDIF 

 2080  

 2081  ! 

 2082  ! Humidity and cloud parameters 

 2083  IF ( humidity ) THEN 

 2084  !$acc parallel present( hom, ngp_2dh_s_inner, sums ) 

 2085  sums(:,44) = sums(:,44) / ngp_2dh_s_inner(:,sr) 

 2086  sums(:,41) = sums(:,41) / ngp_2dh_s_inner(:,sr) 

 2087  hom(:,1,44,sr) = sums(:,44) ! vpt 

 2088  hom(:,1,41,sr) = sums(:,41) ! qv (q) 

 2089  !$acc end parallel 

 2090  IF ( cloud_physics ) THEN 

 2091  !$acc parallel present( hom, ngp_2dh_s_inner, sums ) 

 2092  sums(:,42) = sums(:,42) / ngp_2dh_s_inner(:,sr) 

 2093  sums(:,43) = sums(:,43) / ngp_2dh_s_inner(:,sr) 

 2094  hom(:,1,42,sr) = sums(:,42) ! qv 

 2095  hom(:,1,43,sr) = sums(:,43) ! pt 

 2096  !$acc end parallel 

 2097  ENDIF 

 2098  ENDIF 

 2099  

 2100  ! 

 2101  ! Passive scalar 

 2102  IF ( passive_scalar ) THEN 

 2103  !$acc parallel present( hom, ngp_2dh_s_inner, sums ) 

 2104  sums(:,41) = sums(:,41) / ngp_2dh_s_inner(:,sr) 

 2105  hom(:,1,41,sr) = sums(:,41) ! s (q) 

 2106  !$acc end parallel 

 2107  ENDIF 

 2108  

 2109  ! 

 2110  ! Horizontally averaged profiles of the remaining prognostic variables, 

 2111  ! variances, the total and the perturbation energy (single values in last 

 2112  ! column of sums_l) and some diagnostic quantities. 

 2113  ! NOTE: for simplicity, nzb_s_inner is used below, although strictly 

 2114  !  speaking the following kloop would have to be split up and 

 2115  ! rearranged according to the staggered grid. 

 2116  ! However, this implies no error since staggered velocity components 

 2117  ! are zero at the walls and inside buildings. 

 2118  tn = 0 

 2119  #if defined( __intel_openmp_bug ) 

 2120  !$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper, sums_l_etot, & 

 2121  !$OMP tn, ust, ust2, u2, vst, vst2, v2, w2 ) 

 2122  tn = omp_get_thread_num() 

 2123  #else 

 2124  !$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll, sums_l_eper, sums_l_etot, tn, ust, ust2, u2, vst, vst2, v2, w2 ) 

 2125  !$ tn = omp_get_thread_num() 

 2126  #endif 

 2127  !$OMP DO 

 2128  !$acc parallel loop gang present( e, hom, kh, km, p, pt, w, rflags_invers, rmask, sums_l ) create( s1, s2, s3, s4, s5, s6, s7 ) 

 2129  DO k = nzb, nzt+1 

[1658]  2130  s1 = 0 

 2131  s2 = 0 

 2132  s3 = 0 

 2133  s4 = 0 

 2134  s5 = 0 

 2135  s6 = 0 

 2136  s7 = 0 

[1221]  2137  !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4, s5, s6, s7 ) 

 2138  DO i = nxl, nxr 

 2139  DO j = nys, nyn 

 2140  ! 

 2141  ! Prognostic and diagnostic variables 

 2142  s1 = s1 + w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1) 

 2143  s2 = s2 + e(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1) 

 2144  s3 = s3 + km(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1) 

 2145  s4 = s4 + kh(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1) 

 2146  s5 = s5 + p(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1) 

 2147  s6 = s6 + ( pt(k,j,i)hom(k,1,4,sr) )**2 * rmask(j,i,sr) * & 

 2148  rflags_invers(j,i,k+1) 

 2149  ! 

 2150  ! Higher moments 

 2151  ! (Computation of the skewness of w further below) 

 2152  s7 = s7 + w(k,j,i)**3 * rmask(j,i,sr) * rflags_invers(j,i,k+1) 

 2153  ENDDO 

 2154  ENDDO 

 2155  sums_l(k,3,tn) = s1 

 2156  sums_l(k,8,tn) = s2 

 2157  sums_l(k,9,tn) = s3 

 2158  sums_l(k,10,tn) = s4 

 2159  sums_l(k,40,tn) = s5 

 2160  sums_l(k,33,tn) = s6 

 2161  sums_l(k,38,tn) = s7 

 2162  ENDDO 

[1257]  2163  !$acc end parallel loop 

[1221]  2164  

 2165  IF ( humidity ) THEN 

 2166  !$OMP DO 

 2167  !$acc parallel loop gang present( hom, q, rflags_invers, rmask, sums_l ) create( s1 ) 

 2168  DO k = nzb, nzt+1 

[1658]  2169  s1 = 0 

[1221]  2170  !$acc loop vector collapse( 2 ) reduction( +: s1 ) 

 2171  DO i = nxl, nxr 

 2172  DO j = nys, nyn 

 2173  s1 = s1 + ( q(k,j,i)hom(k,1,41,sr) )**2 * rmask(j,i,sr) * & 

 2174  rflags_invers(j,i,k+1) 

 2175  ENDDO 

 2176  ENDDO 

 2177  sums_l(k,70,tn) = s1 

 2178  ENDDO 

[1257]  2179  !$acc end parallel loop 

[1221]  2180  ENDIF 

 2181  

 2182  ! 

 2183  ! Total and perturbation energy for the total domain (being 

 2184  ! collected in the last column of sums_l). 

[1658]  2185  s1 = 0 

[1221]  2186  !$OMP DO 

 2187  !$acc parallel loop collapse(3) present( rflags_invers, rmask, u, v, w ) reduction(+:s1) 

 2188  DO i = nxl, nxr 

 2189  DO j = nys, nyn 

 2190  DO k = nzb, nzt+1 

[1353]  2191  s1 = s1 + 0.5_wp * & 

 2192  ( u(k,j,i)**2 + v(k,j,i)**2 + w(k,j,i)**2 ) * & 

[1221]  2193  rmask(j,i,sr) * rflags_invers(j,i,k+1) 

 2194  ENDDO 

 2195  ENDDO 

 2196  ENDDO 

[1257]  2197  !$acc end parallel loop 

[1221]  2198  !$acc parallel present( sums_l ) 

 2199  sums_l(nzb+4,pr_palm,tn) = s1 

 2200  !$acc end parallel 

 2201  

 2202  !$OMP DO 

 2203  !$acc parallel present( rmask, sums_l, us, usws, vsws, ts ) create( s1, s2, s3, s4 ) 

[1658]  2204  s1 = 0 

 2205  s2 = 0 

 2206  s3 = 0 

 2207  s4 = 0 

[1221]  2208  !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4 ) 

 2209  DO i = nxl, nxr 

 2210  DO j = nys, nyn 

 2211  ! 

 2212  ! 2Darrays (being collected in the last column of sums_l) 

 2213  s1 = s1 + us(j,i) * rmask(j,i,sr) 

 2214  s2 = s2 + usws(j,i) * rmask(j,i,sr) 

 2215  s3 = s3 + vsws(j,i) * rmask(j,i,sr) 

 2216  s4 = s4 + ts(j,i) * rmask(j,i,sr) 

 2217  ENDDO 

 2218  ENDDO 

 2219  sums_l(nzb,pr_palm,tn) = s1 

 2220  sums_l(nzb+1,pr_palm,tn) = s2 

 2221  sums_l(nzb+2,pr_palm,tn) = s3 

 2222  sums_l(nzb+3,pr_palm,tn) = s4 

 2223  !$acc end parallel 

 2224  

 2225  IF ( humidity ) THEN 

 2226  !$acc parallel present( qs, rmask, sums_l ) create( s1 ) 

[1658]  2227  s1 = 0 

[1221]  2228  !$acc loop vector collapse( 2 ) reduction( +: s1 ) 

 2229  DO i = nxl, nxr 

 2230  DO j = nys, nyn 

 2231  s1 = s1 + qs(j,i) * rmask(j,i,sr) 

 2232  ENDDO 

 2233  ENDDO 

 2234  sums_l(nzb+12,pr_palm,tn) = s1 

 2235  !$acc end parallel 

 2236  ENDIF 

 2237  

 2238  ! 

 2239  ! Computation of statistics when wsscheme is not used. Else these 

 2240  ! quantities are evaluated in the advection routines. 

 2241  IF ( .NOT. ws_scheme_mom .OR. sr /= 0 ) THEN 

 2242  

 2243  !$OMP DO 

 2244  !$acc parallel loop gang present( u, v, w, rflags_invers, rmask, sums_l ) create( s1, s2, s3, s4, ust2, vst2, w2 ) 

 2245  DO k = nzb, nzt+1 

[1658]  2246  s1 = 0 

 2247  s2 = 0 

 2248  s3 = 0 

 2249  s4 = 0 

[1221]  2250  !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4 ) 

 2251  DO i = nxl, nxr 

 2252  DO j = nys, nyn 

 2253  ust2 = ( u(k,j,i)  hom(k,1,1,sr) )**2 

 2254  vst2 = ( v(k,j,i)  hom(k,1,2,sr) )**2 

 2255  w2 = w(k,j,i)**2 

 2256  

 2257  s1 = s1 + ust2 * rmask(j,i,sr) * rflags_invers(j,i,k+1) 

 2258  s2 = s2 + vst2 * rmask(j,i,sr) * rflags_invers(j,i,k+1) 

 2259  s3 = s3 + w2 * rmask(j,i,sr) * rflags_invers(j,i,k+1) 

 2260  ! 

 2261  ! Perturbation energy 

[1353]  2262  s4 = s4 + 0.5_wp * ( ust2 + vst2 + w2 ) * rmask(j,i,sr) * & 

[1221]  2263  rflags_invers(j,i,k+1) 

 2264  ENDDO 

 2265  ENDDO 

 2266  sums_l(k,30,tn) = s1 

 2267  sums_l(k,31,tn) = s2 

 2268  sums_l(k,32,tn) = s3 

 2269  sums_l(k,34,tn) = s4 

 2270  ENDDO 

[1257]  2271  !$acc end parallel loop 

[1221]  2272  ! 

 2273  ! Total perturbation TKE 

 2274  !$OMP DO 

 2275  !$acc parallel present( sums_l ) create( s1 ) 

[1658]  2276  s1 = 0 

[1221]  2277  !$acc loop reduction( +: s1 ) 

 2278  DO k = nzb, nzt+1 

 2279  s1 = s1 + sums_l(k,34,tn) 

 2280  ENDDO 

 2281  sums_l(nzb+5,pr_palm,tn) = s1 

 2282  !$acc end parallel 

 2283  

 2284  ENDIF 

 2285  

 2286  ! 

 2287  ! Horizontally averaged profiles of the vertical fluxes 

 2288  

 2289  ! 

 2290  ! Subgridscale fluxes. 

 2291  ! WARNING: If a Prandtllayer is used (k=nzb for flat terrain), the fluxes 

 2292  !  should be calculated there in a different way. This is done 

 2293  ! in the next loop further below, where results from this loop are 

 2294  ! overwritten. However, THIS WORKS IN CASE OF FLAT TERRAIN ONLY! 

 2295  ! The nonflat case still has to be handled. 

 2296  ! NOTE: for simplicity, nzb_s_inner is used below, although 

 2297  !  strictly speaking the following kloop would have to be 

 2298  ! split up according to the staggered grid. 

 2299  ! However, this implies no error since staggered velocity 

 2300  ! components are zero at the walls and inside buildings. 

 2301  !$OMP DO 

 2302  !$acc parallel loop gang present( ddzu, kh, km, pt, u, v, w, rflags_invers, rmask, sums_l ) create( s1, s2, s3 ) 

 2303  DO k = nzb, nzt_diff 

[1658]  2304  s1 = 0 

 2305  s2 = 0 

 2306  s3 = 0 

[1221]  2307  !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 ) 

 2308  DO i = nxl, nxr 

 2309  DO j = nys, nyn 

 2310  

 2311  ! 

 2312  ! Momentum flux w"u" 

[1353]  2313  s1 = s1  0.25_wp * ( & 

[1221]  2314  km(k,j,i)+km(k+1,j,i)+km(k,j,i1)+km(k+1,j,i1) & 

 2315  ) * ( & 

 2316  ( u(k+1,j,i)  u(k,j,i) ) * ddzu(k+1) & 

 2317  + ( w(k,j,i)  w(k,j,i1) ) * ddx & 

 2318  ) & 

 2319  * rmask(j,i,sr) * rflags_invers(j,i,k+1) 

 2320  ! 

 2321  ! Momentum flux w"v" 

[1353]  2322  s2 = s2  0.25_wp * ( & 

[1221]  2323  km(k,j,i)+km(k+1,j,i)+km(k,j1,i)+km(k+1,j1,i) & 

 2324  ) * ( & 

 2325  ( v(k+1,j,i)  v(k,j,i) ) * ddzu(k+1) & 

 2326  + ( w(k,j,i)  w(k,j1,i) ) * ddy & 

 2327  ) & 

 2328  * rmask(j,i,sr) * rflags_invers(j,i,k+1) 

 2329  ! 

 2330  ! Heat flux w"pt" 

[1353]  2331  s3 = s3  0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) ) & 

 2332  * ( pt(k+1,j,i)  pt(k,j,i) ) & 

 2333  * ddzu(k+1) * rmask(j,i,sr) & 

 2334  * rflags_invers(j,i,k+1) 

[1221]  2335  ENDDO 

 2336  ENDDO 

 2337  sums_l(k,12,tn) = s1 

 2338  sums_l(k,14,tn) = s2 

 2339  sums_l(k,16,tn) = s3 

 2340  ENDDO 

[1257]  2341  !$acc end parallel loop 

[1221]  2342  

 2343  ! 

 2344  ! Salinity flux w"sa" 

 2345  IF ( ocean ) THEN 

 2346  !$acc parallel loop gang present( ddzu, kh, sa, rflags_invers, rmask, sums_l ) create( s1 ) 

 2347  DO k = nzb, nzt_diff 

[1658]  2348  s1 = 0 

[1221]  2349  !$acc loop vector collapse( 2 ) reduction( +: s1 ) 

 2350  DO i = nxl, nxr 

 2351  DO j = nys, nyn 

[1353]  2352  s1 = s1  0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) ) & 

 2353  * ( sa(k+1,j,i)  sa(k,j,i) ) & 

 2354  * ddzu(k+1) * rmask(j,i,sr) & 

 2355  * rflags_invers(j,i,k+1) 

[1221]  2356  ENDDO 

 2357  ENDDO 

 2358  sums_l(k,65,tn) = s1 

 2359  ENDDO 

[1257]  2360  !$acc end parallel loop 

[1221]  2361  ENDIF 

 2362  

 2363  ! 

 2364  ! Buoyancy flux, water flux (humidity flux) w"q" 

 2365  IF ( humidity ) THEN 

 2366  

 2367  !$acc parallel loop gang present( ddzu, kh, q, vpt, rflags_invers, rmask, sums_l ) create( s1, s2 ) 

 2368  DO k = nzb, nzt_diff 

[1658]  2369  s1 = 0 

 2370  s2 = 0 

[1221]  2371  !$acc loop vector collapse( 2 ) reduction( +: s1, s2 ) 

 2372  DO i = nxl, nxr 

 2373  DO j = nys, nyn 

[1353]  2374  s1 = s1  0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) ) & 

 2375  * ( vpt(k+1,j,i)  vpt(k,j,i) ) & 

[1374]  2376  * ddzu(k+1) * rmask(j,i,sr) & 

[1353]  2377  * rflags_invers(j,i,k+1) 

 2378  s2 = s2  0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) ) & 

 2379  * ( q(k+1,j,i)  q(k,j,i) ) & 

 2380  * ddzu(k+1) * rmask(j,i,sr) & 

 2381  * rflags_invers(j,i,k+1) 

[1221]  2382  ENDDO 

 2383  ENDDO 

 2384  sums_l(k,45,tn) = s1 

 2385  sums_l(k,48,tn) = s2 

 2386  ENDDO 

[1257]  2387  !$acc end parallel loop 

[1221]  2388  

 2389  IF ( cloud_physics ) THEN 

 2390  

 2391  !$acc parallel loop gang present( ddzu, kh, q, ql, rflags_invers, rmask, sums_l ) create( s1 ) 

 2392  DO k = nzb, nzt_diff 

[1658]  2393  s1 = 0 

[1221]  2394  !$acc loop vector collapse( 2 ) reduction( +: s1 ) 

 2395  DO i = nxl, nxr 

 2396  DO j = nys, nyn 

[1353]  2397  s1 = s1  0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) ) & 

 2398  * ( ( q(k+1,j,i)  ql(k+1,j,i) ) & 

 2399   ( q(k,j,i)  ql(k,j,i) ) ) & 

 2400  * ddzu(k+1) * rmask(j,i,sr) & 

 2401  * rflags_invers(j,i,k+1) 

[1221]  2402  ENDDO 

 2403  ENDDO 

 2404  sums_l(k,51,tn) = s1 

 2405  ENDDO 

[1257]  2406  !$acc end parallel loop 

[1221]  2407  

 2408  ENDIF 

 2409  

 2410  ENDIF 

 2411  ! 

 2412  ! Passive scalar flux 

 2413  IF ( passive_scalar ) THEN 

 2414  

 2415  !$acc parallel loop gang present( ddzu, kh, q, rflags_invers, rmask, sums_l ) create( s1 ) 

 2416  DO k = nzb, nzt_diff 

[1658]  2417  s1 = 0 

[1221]  2418  !$acc loop vector collapse( 2 ) reduction( +: s1 ) 

 2419  DO i = nxl, nxr 

 2420  DO j = nys, nyn 

[1353]  2421  s1 = s1  0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) ) & 

 2422  * ( q(k+1,j,i)  q(k,j,i) ) & 

 2423  * ddzu(k+1) * rmask(j,i,sr) & 

 2424  * rflags_invers(j,i,k+1) 

[1221]  2425  ENDDO 

 2426  ENDDO 

 2427  sums_l(k,48,tn) = s1 

 2428  ENDDO 

[1257]  2429  !$acc end parallel loop 

[1221]  2430  

 2431  ENDIF 

 2432  

 2433  IF ( use_surface_fluxes ) THEN 

 2434  

 2435  !$OMP DO 

 2436  !$acc parallel present( rmask, shf, sums_l, usws, vsws ) create( s1, s2, s3, s4, s5 ) 

[1658]  2437  s1 = 0 

 2438  s2 = 0 

 2439  s3 = 0 

 2440  s4 = 0 

 2441  s5 = 0 

[1221]  2442  !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4, s5 ) 

 2443  DO i = nxl, nxr 

 2444  DO j = nys, nyn 

 2445  ! 

 2446  ! Subgridscale fluxes in the Prandtl layer 

 2447  s1 = s1 + usws(j,i) * rmask(j,i,sr) ! w"u" 

 2448  s2 = s2 + vsws(j,i) * rmask(j,i,sr) ! w"v" 

 2449  s3 = s3 + shf(j,i) * rmask(j,i,sr) ! w"pt" 

[1353]  2450  s4 = s4 + 0.0_wp * rmask(j,i,sr) ! u"pt" 

 2451  s5 = s5 + 0.0_wp * rmask(j,i,sr) ! v"pt" 

[1221]  2452  ENDDO 

 2453  ENDDO 

 2454  sums_l(nzb,12,tn) = s1 

 2455  sums_l(nzb,14,tn) = s2 

 2456  sums_l(nzb,16,tn) = s3 

 2457  sums_l(nzb,58,tn) = s4 

 2458  sums_l(nzb,61,tn) = s5 

 2459  !$acc end parallel 

 2460  

 2461  IF ( ocean ) THEN 

 2462  

 2463  !$OMP DO 

 2464  !$acc parallel present( rmask, saswsb, sums_l ) create( s1 ) 

[1658]  2465  s1 = 0 

[1221]  2466  !$acc loop vector collapse( 2 ) reduction( +: s1 ) 

 2467  DO i = nxl, nxr 

 2468  DO j = nys, nyn 

 2469  s1 = s1 + saswsb(j,i) * rmask(j,i,sr) ! w"sa" 

 2470  ENDDO 

 2471  ENDDO 

 2472  sums_l(nzb,65,tn) = s1 

 2473  !$acc end parallel 

 2474  

 2475  ENDIF 

 2476  

 2477  IF ( humidity ) THEN 

 2478  

 2479  !$OMP DO 

 2480  !$acc parallel present( pt, q, qsws, rmask, shf, sums_l ) create( s1, s2 ) 

[1658]  2481  s1 = 0 

 2482  s2 = 0 

[1221]  2483  !$acc loop vector collapse( 2 ) reduction( +: s1, s2 ) 

 2484  DO i = nxl, nxr 

 2485  DO j = nys, nyn 

 2486  s1 = s1 + qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv") 

[1353]  2487  s2 = s2 + ( ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) * shf(j,i) & 

 2488  + 0.61_wp * pt(nzb,j,i) * qsws(j,i) ) 

[1221]  2489  ENDDO 

 2490  ENDDO 

 2491  sums_l(nzb,48,tn) = s1 

 2492  sums_l(nzb,45,tn) = s2 

 2493  !$acc end parallel 

 2494  

 2495  IF ( cloud_droplets ) THEN 

 2496  

 2497  !$OMP DO 

 2498  !$acc parallel present( pt, q, ql, qsws, rmask, shf, sums_l ) create( s1 ) 

[1658]  2499  s1 = 0 

[1221]  2500  !$acc loop vector collapse( 2 ) reduction( +: s1 ) 

 2501  DO i = nxl, nxr 

 2502  DO j = nys, nyn 

[1353]  2503  s1 = s1 + ( ( 1.0_wp + & 

 2504  0.61_wp * q(nzb,j,i)  ql(nzb,j,i) ) * & 

 2505  shf(j,i) + 0.61_wp * pt(nzb,j,i) * qsws(j,i) ) 

[1221]  2506  ENDDO 

 2507  ENDDO 

 2508  sums_l(nzb,45,tn) = s1 

 2509  !$acc end parallel 

 2510  

 2511  ENDIF 

 2512  

 2513  IF ( cloud_physics ) THEN 

 2514  

 2515  !$OMP DO 

 2516  !$acc parallel present( qsws, rmask, sums_l ) create( s1 ) 

[1658]  2517  s1 = 0 

[1221]  2518  !$acc loop vector collapse( 2 ) reduction( +: s1 ) 

 2519  DO i = nxl, nxr 

 2520  DO j = nys, nyn 

 2521  ! 

 2522  ! Formula does not work if ql(nzb) /= 0.0 

 2523  s1 = s1 + qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv") 

 2524  ENDDO 

 2525  ENDDO 

 2526  sums_l(nzb,51,tn) = s1 

 2527  !$acc end parallel 

 2528  

 2529  ENDIF 

 2530  

 2531  ENDIF 

 2532  

 2533  IF ( passive_scalar ) THEN 

 2534  

 2535  !$OMP DO 

 2536  !$acc parallel present( qsws, rmask, sums_l ) create( s1 ) 

[1658]  2537  s1 = 0 

[1221]  2538  !$acc loop vector collapse( 2 ) reduction( +: s1 ) 

 2539  DO i = nxl, nxr 

 2540  DO j = nys, nyn 

 2541  s1 = s1 + qsws(j,i) * rmask(j,i,sr) ! w"q" (w"qv") 

 2542  ENDDO 

 2543  ENDDO 

 2544  sums_l(nzb,48,tn) = s1 

 2545  !$acc end parallel 

 2546  

 2547  ENDIF 

 2548  

 2549  ENDIF 

 2550  

 2551  ! 

 2552  ! Subgridscale fluxes at the top surface 

 2553  IF ( use_top_fluxes ) THEN 

 2554  

 2555  !$OMP DO 

 2556  !$acc parallel present( rmask, sums_l, tswst, uswst, vswst ) create( s1, s2, s3, s4, s5 ) 

[1658]  2557  s1 = 0 

 2558  s2 = 0 

 2559  s3 = 0 

 2560  s4 = 0 

 2561  s5 = 0 

[1221]  2562  !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3, s4, s5 ) 

 2563  DO i = nxl, nxr 

 2564  DO j = nys, nyn 

 2565  s1 = s1 + uswst(j,i) * rmask(j,i,sr) ! w"u" 

 2566  s2 = s2 + vswst(j,i) * rmask(j,i,sr) ! w"v" 

 2567  s3 = s3 + tswst(j,i) * rmask(j,i,sr) ! w"pt" 

[1353]  2568  s4 = s4 + 0.0_wp * rmask(j,i,sr) ! u"pt" 

 2569  s5 = s5 + 0.0_wp * rmask(j,i,sr) ! v"pt" 

[1221]  2570  ENDDO 

 2571  ENDDO 

 2572  sums_l(nzt:nzt+1,12,tn) = s1 

 2573  sums_l(nzt:nzt+1,14,tn) = s2 

 2574  sums_l(nzt:nzt+1,16,tn) = s3 

 2575  sums_l(nzt:nzt+1,58,tn) = s4 

 2576  sums_l(nzt:nzt+1,61,tn) = s5 

 2577  !$acc end parallel 

 2578  

 2579  IF ( ocean ) THEN 

 2580  

 2581  !$OMP DO 

 2582  !$acc parallel present( rmask, saswst, sums_l ) create( s1 ) 

[1658]  2583  s1 = 0 

[1221]  2584  !$acc loop vector collapse( 2 ) reduction( +: s1 ) 

 2585  DO i = nxl, nxr 

 2586  DO j = nys, nyn 

 2587  s1 = s1 + saswst(j,i) * rmask(j,i,sr) ! w"sa" 

 2588  ENDDO 

 2589  ENDDO 

 2590  sums_l(nzt,65,tn) = s1 

 2591  !$acc end parallel 

 2592  

 2593  ENDIF 

 2594  

 2595  IF ( humidity ) THEN 

 2596  

 2597  !$OMP DO 

 2598  !$acc parallel present( pt, q, qswst, rmask, tswst, sums_l ) create( s1, s2 ) 

[1658]  2599  s1 = 0 

 2600  s2 = 0 

[1221]  2601  !$acc loop vector collapse( 2 ) reduction( +: s1, s2 ) 

 2602  DO i = nxl, nxr 

 2603  DO j = nys, nyn 

 2604  s1 = s1 + qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv") 

[1353]  2605  s2 = s2 + ( ( 1.0_wp + 0.61_wp * q(nzt,j,i) ) * tswst(j,i) +& 

 2606  0.61_wp * pt(nzt,j,i) * qswst(j,i) ) 

[1221]  2607  ENDDO 

 2608  ENDDO 

 2609  sums_l(nzt,48,tn) = s1 

 2610  sums_l(nzt,45,tn) = s2 

 2611  !$acc end parallel 

 2612  

 2613  IF ( cloud_droplets ) THEN 

 2614  

 2615  !$OMP DO 

 2616  !$acc parallel present( pt, q, ql, qswst, rmask, tswst, sums_l ) create( s1 ) 

[1658]  2617  s1 = 0 

[1221]  2618  !$acc loop vector collapse( 2 ) reduction( +: s1 ) 

 2619  DO i = nxl, nxr 

 2620  DO j = nys, nyn 

[1353]  2621  s1 = s1 + ( ( 1.0_wp + & 

 2622  0.61_wp * q(nzt,j,i)  ql(nzt,j,i) ) * & 

 2623  tswst(j,i) + & 

 2624  0.61_wp * pt(nzt,j,i) * qswst(j,i) ) 

[1221]  2625  ENDDO 

 2626  ENDDO 

 2627  sums_l(nzt,45,tn) = s1 

 2628  !$acc end parallel 

 2629  

 2630  ENDIF 

 2631  

 2632  IF ( cloud_physics ) THEN 

 2633  

 2634  !$OMP DO 

 2635  !$acc parallel present( qswst, rmask, sums_l ) create( s1 ) 

[1658]  2636  s1 = 0 

[1221]  2637  !$acc loop vector collapse( 2 ) reduction( +: s1 ) 

 2638  DO i = nxl, nxr 

 2639  DO j = nys, nyn 

 2640  ! 

 2641  ! Formula does not work if ql(nzb) /= 0.0 

 2642  s1 = s1 + qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv") 

 2643  ENDDO 

 2644  ENDDO 

 2645  sums_l(nzt,51,tn) = s1 

 2646  !$acc end parallel 

 2647  

 2648  ENDIF 

 2649  

 2650  ENDIF 

 2651  

 2652  IF ( passive_scalar ) THEN 

 2653  

 2654  !$OMP DO 

 2655  !$acc parallel present( qswst, rmask, sums_l ) create( s1 ) 

[1658]  2656  s1 = 0 

[1221]  2657  !$acc loop vector collapse( 2 ) reduction( +: s1 ) 

 2658  DO i = nxl, nxr 

 2659  DO j = nys, nyn 

 2660  s1 = s1 + qswst(j,i) * rmask(j,i,sr) ! w"q" (w"qv") 

 2661  ENDDO 

 2662  ENDDO 

 2663  sums_l(nzt,48,tn) = s1 

 2664  !$acc end parallel 

 2665  

 2666  ENDIF 

 2667  

 2668  ENDIF 

 2669  

 2670  ! 

 2671  ! Resolved fluxes (can be computed for all horizontal points) 

 2672  ! NOTE: for simplicity, nzb_s_inner is used below, although strictly 

 2673  !  speaking the following kloop would have to be split up and 

 2674  ! rearranged according to the staggered grid. 

 2675  !$acc parallel loop gang present( hom, pt, rflags_invers, rmask, sums_l, u, v, w ) create( s1, s2, s3 ) 

 2676  DO k = nzb, nzt_diff 

[1658]  2677  s1 = 0 

 2678  s2 = 0 

 2679  s3 = 0 

[1221]  2680  !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 ) 

 2681  DO i = nxl, nxr 

 2682  DO j = nys, nyn 

[1353]  2683  ust = 0.5_wp * ( u(k,j,i)  hom(k,1,1,sr) + & 

 2684  u(k+1,j,i)  hom(k+1,1,1,sr) ) 

 2685  vst = 0.5_wp * ( v(k,j,i)  hom(k,1,2,sr) + & 

 2686  v(k+1,j,i)  hom(k+1,1,2,sr) ) 

 2687  pts = 0.5_wp * ( pt(k,j,i)  hom(k,1,4,sr) + & 

 2688  pt(k+1,j,i)  hom(k+1,1,4,sr) ) 

[1221]  2689  ! 

 2690  ! Higher moments 

 2691  s1 = s1 + pts * w(k,j,i)**2 * rmask(j,i,sr) * rflags_invers(j,i,k+1) 

 2692  s2 = s2 + pts**2 * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1) 

 2693  ! 

 2694  ! Energy flux w*e* (has to be adjusted?) 

[1353]  2695  s3 = s3 + w(k,j,i) * 0.5_wp * ( ust**2 + vst**2 + w(k,j,i)**2 )& 

[1221]  2696  * rmask(j,i,sr) * rflags_invers(j,i,k+1) 

 2697  ENDDO 

 2698  ENDDO 

 2699  sums_l(k,35,tn) = s1 

 2700  sums_l(k,36,tn) = s2 

 2701  sums_l(k,37,tn) = s3 

 2702  ENDDO 

[1257]  2703  !$acc end parallel loop 

[1221]  2704  

 2705  ! 

 2706  ! Salinity flux and density (density does not belong to here, 

 2707  ! but so far there is no other suitable place to calculate) 

 2708  IF ( ocean ) THEN 

 2709  

[1567]  2710  IF( .NOT. ws_scheme_sca .OR. sr /= 0 ) THEN 

[1221]  2711  

 2712  !$acc parallel loop gang present( hom, rflags_invers, rmask, sa, sums_l, w ) create( s1 ) 

 2713  DO k = nzb, nzt_diff 

[1658]  2714  s1 = 0 

[1221]  2715  !$acc loop vector collapse( 2 ) reduction( +: s1 ) 

 2716  DO i = nxl, nxr 

 2717  DO j = nys, nyn 

[1353]  2718  s1 = s1 + 0.5_wp * ( sa(k,j,i)  hom(k,1,23,sr) + & 

 2719  sa(k+1,j,i)  hom(k+1,1,23,sr) ) & 

 2720  * w(k,j,i) * rmask(j,i,sr) & 

 2721  * rflags_invers(j,i,k+1) 

[1221]  2722  ENDDO 

 2723  ENDDO 

 2724  sums_l(k,66,tn) = s1 

 2725  ENDDO 

[1257]  2726  !$acc end parallel loop 

[1221]  2727  

 2728  ENDIF 

 2729  

 2730  !$acc parallel loop gang present( rflags_invers, rho, prho, rmask, sums_l ) create( s1, s2 ) 

 2731  DO k = nzb, nzt_diff 

[1658]  2732  s1 = 0 

 2733  s2 = 0 

[1221]  2734  !$acc loop vector collapse( 2 ) reduction( +: s1, s2 ) 

 2735  DO i = nxl, nxr 

 2736  DO j = nys, nyn 

 2737  s1 = s1 + rho(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1) 

 2738  s2 = s2 + prho(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1) 

 2739  ENDDO 

 2740  ENDDO 

 2741  sums_l(k,64,tn) = s1 

 2742  sums_l(k,71,tn) = s2 

 2743  ENDDO 

[1257]  2744  !$acc end parallel loop 

[1221]  2745  

 2746  ENDIF 

 2747  

 2748  ! 

 2749  ! Buoyancy flux, water flux, humidity flux, liquid water 

 2750  ! content, rain drop concentration and rain water content 

 2751  IF ( humidity ) THEN 

 2752  

 2753  IF ( cloud_physics .OR. cloud_droplets ) THEN 

 2754  

 2755  !$acc parallel loop gang present( hom, rflags_invers, rmask, sums_l, vpt, w ) create( s1 ) 

 2756  DO k = nzb, nzt_diff 

[1658]  2757  s1 = 0 

[1221]  2758  !$acc loop vector collapse( 2 ) reduction( +: s1 ) 

 2759  DO i = nxl, nxr 

 2760  DO j = nys, nyn 

[1353]  2761  s1 = s1 + 0.5_wp * ( vpt(k,j,i)  hom(k,1,44,sr) + & 

 2762  vpt(k+1,j,i)  hom(k+1,1,44,sr) ) * & 

 2763  w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1) 

[1221]  2764  ENDDO 

 2765  ENDDO 

 2766  sums_l(k,46,tn) = s1 

 2767  ENDDO 

[1257]  2768  !$acc end parallel loop 

[1221]  2769  

 2770  IF ( .NOT. cloud_droplets ) THEN 

 2771  

 2772  !$acc parallel loop gang present( hom, q, ql, rflags_invers, rmask, sums_l, w ) create( s1 ) 

 2773  DO k = nzb, nzt_diff 

[1658]  2774  s1 = 0 

[1221]  2775  !$acc loop vector collapse( 2 ) reduction( +: s1 ) 

 2776  DO i = nxl, nxr 

 2777  DO j = nys, nyn 

[1353]  2778  s1 = s1 + 0.5_wp * ( ( q(k,j,i)  ql(k,j,i) )  hom(k,1,42,sr) + & 

 2779  ( q(k+1,j,i)  ql(k+1,j,i) )  hom(k+1,1,42,sr) ) & 

 2780  * w(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1) 

[1221]  2781  ENDDO 

 2782  ENDDO 

 2783  sums_l(k,52,tn) = s1 

 2784  ENDDO 

[1257]  2785  !$acc end parallel loop 

[1221]  2786  

 2787  IF ( icloud_scheme == 0 ) THEN 

 2788  

 2789  !$acc parallel loop gang present( qc, ql, rflags_invers, rmask, sums_l ) create( s1, s2 ) 

 2790  DO k = nzb, nzt_diff 

[1658]  2791  s1 = 0 

[1221]  2792  !$acc loop vector collapse( 2 ) reduction( +: s1, s2 ) 

 2793  DO i = nxl, nxr 

 2794  DO j = nys, nyn 

 2795  s1 = s1 + ql(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1) 

 2796  s2 = s2 + qc(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1) 

 2797  ENDDO 

 2798  ENDDO 

 2799  sums_l(k,54,tn) = s1 

 2800  sums_l(k,75,tn) = s2 

 2801  ENDDO 

[1257]  2802  !$acc end parallel loop 

[1221]  2803  

 2804  IF ( precipitation ) THEN 

 2805  

 2806  !$acc parallel loop gang present( nr, qr, prr, rflags_invers, rmask, sums_l ) create( s1, s2, s3 ) 

 2807  DO k = nzb, nzt_diff 

[1658]  2808  s1 = 0 

 2809  s2 = 0 

 2810  s3 = 0 

[1221]  2811  !$acc loop vector collapse( 2 ) reduction( +: s1, s2, s3 ) 

 2812  DO i = nxl, nxr 

 2813  DO j = nys, nyn 

 2814  s1 = s1 + nr(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1) 

 2815  s2 = s2 + qr(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1) 

 2816  s3 = s3 + prr(k,j,i) * rmask(j,i,sr) * rflags_invers(j,i,k+1) 

 2817  ENDDO 

 2818  ENDDO 

 2819  sums_l(k,73,tn) = s1 

 2820  sums_l(k,74,tn) = s2 

 2821  sums_l(k,76,tn) = s3 

 2822  ENDDO 


