Changeset 4646


Ignore:
Timestamp:
Aug 24, 2020 4:02:40 PM (2 months ago)
Author:
raasch
Message:

files re-formatted to follow the PALM coding standard

Location:
palm/trunk/SOURCE
Files:
7 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/fft_xy_mod.f90

    r4370 r4646  
    11!> @file fft_xy_mod.f90
    2 !------------------------------------------------------------------------------!
     2!--------------------------------------------------------------------------------------------------!
    33! This file is part of the PALM model system.
    44!
    5 ! PALM is free software: you can redistribute it and/or modify it under the
    6 ! terms of the GNU General Public License as published by the Free Software
    7 ! Foundation, either version 3 of the License, or (at your option) any later
    8 ! version.
    9 !
    10 ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
    11 ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
    12 ! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
    13 !
    14 ! You should have received a copy of the GNU General Public License along with
    15 ! PALM. If not, see <http://www.gnu.org/licenses/>.
     5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
     6! Public License as published by the Free Software Foundation, either version 3 of the License, or
     7! (at your option) any later version.
     8!
     9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
     10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
     11! Public License for more details.
     12!
     13! You should have received a copy of the GNU General Public License along with PALM. If not, see
     14! <http://www.gnu.org/licenses/>.
    1615!
    1716! Copyright 1997-2020 Leibniz Universitaet Hannover
    18 !------------------------------------------------------------------------------!
     17!--------------------------------------------------------------------------------------------------!
    1918!
    2019! Current revisions:
    2120! -----------------
    22 ! 
    23 ! 
     21!
     22!
    2423! Former revisions:
    2524! -----------------
    2625! $Id$
     26! file re-formatted to follow the PALM coding standard
     27!
     28! 4370 2020-01-10 14:00:44Z raasch
    2729! bugfix for Temperton-fft usage on GPU
    28 ! 
     30!
    2931! 4366 2020-01-09 08:12:43Z raasch
    3032! Vectorized Temperton-fft added
    31 ! 
     33!
    3234! 4360 2020-01-07 11:25:50Z suehring
    3335! Corrected "Former revisions" section
    34 ! 
     36!
    3537! 4069 2019-07-01 14:05:51Z Giersch
    3638! Code added to avoid compiler warnings
    37 ! 
     39!
    3840! 3655 2019-01-07 16:51:22Z knoop
    3941! OpenACC port for SPEC
     
    5052!------------------------------------------------------------------------------!
    5153 MODULE fft_xy
    52  
    53 
    54     USE control_parameters,                                                    &
     54
     55
     56    USE control_parameters,                                                                        &
    5557        ONLY:  fft_method, loop_optimization, message_string
    56        
     58
    5759    USE cuda_fft_interfaces
    58        
    59     USE indices,                                                               &
     60
     61    USE indices,                                                                                   &
    6062        ONLY:  nx, ny, nz
    6163
     
    6769
    6870    USE kinds
    69    
    70     USE singleton,                                                             &
     71
     72    USE singleton,                                                                                 &
    7173        ONLY: fftn
    72    
     74
    7375    USE temperton_fft
    74    
    75     USE transpose_indices,                                                     &
     76
     77    USE transpose_indices,                                                                         &
    7678        ONLY:  nxl_y, nxr_y, nyn_x, nys_x, nzb_x, nzb_y, nzt_x, nzt_y
    7779
     
    7981
    8082    PRIVATE
    81     PUBLIC fft_x, fft_x_1d, fft_y, fft_y_1d, fft_init, fft_x_m, fft_y_m, f_vec_x, temperton_fft_vec
     83    PUBLIC fft_init, f_vec_x, fft_x, fft_x_1d, fft_x_m, fft_y, fft_y_1d, fft_y_m, temperton_fft_vec
    8284
    8385    INTEGER(iwp), DIMENSION(:), ALLOCATABLE, SAVE ::  ifax_x  !<
     
    9193    REAL(wp), SAVE ::  sqr_dnx  !<
    9294    REAL(wp), SAVE ::  sqr_dny  !<
    93    
    94     REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE ::  trigs_x  !< 
     95
     96    REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE ::  trigs_x  !<
    9597    REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE ::  trigs_y  !<
    9698
     
    98100
    99101#if defined( __ibm )
    100     INTEGER(iwp), PARAMETER ::  nau1 = 20000  !< 
     102    INTEGER(iwp), PARAMETER ::  nau1 = 20000  !<
    101103    INTEGER(iwp), PARAMETER ::  nau2 = 22000  !<
    102104!
    103 !-- The following working arrays contain tables and have to be "save" and
    104 !-- shared in OpenMP sense
    105     REAL(wp), DIMENSION(nau1), SAVE ::  aux1  !<
     105!-- The following working arrays contain tables and have to be "save" and shared in OpenMP sense
     106    REAL(wp), DIMENSION(nau1), SAVE ::  aux1  !<
    106107    REAL(wp), DIMENSION(nau1), SAVE ::  auy1  !<
    107     REAL(wp), DIMENSION(nau1), SAVE ::  aux3  !< 
     108    REAL(wp), DIMENSION(nau1), SAVE ::  aux3  !<
    108109    REAL(wp), DIMENSION(nau1), SAVE ::  auy3  !<
    109    
     110
    110111#elif defined( __nec_fft )
    111112    INTEGER(iwp), SAVE ::  nz1  !<
    112    
     113
    113114    REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE ::  trig_xb  !<
    114     REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE ::  trig_xf  !< 
     115    REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE ::  trig_xf  !<
    115116    REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE ::  trig_yb  !<
    116117    REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE ::  trig_yf  !<
    117    
     118
    118119#elif defined( __cuda_fft )
    119120    INTEGER(C_INT), SAVE ::  plan_xf  !<
     
    126127#if defined( __fftw )
    127128    INCLUDE  'fftw3.f03'
     129    COMPLEX(KIND=C_DOUBLE_COMPLEX), DIMENSION(:), ALLOCATABLE, SAVE ::  x_out  !<
     130    COMPLEX(KIND=C_DOUBLE_COMPLEX), DIMENSION(:), ALLOCATABLE, SAVE ::  y_out  !<
     131
    128132    INTEGER(KIND=C_INT) ::  nx_c  !<
    129133    INTEGER(KIND=C_INT) ::  ny_c  !<
    130    
    131     COMPLEX(KIND=C_DOUBLE_COMPLEX), DIMENSION(:), ALLOCATABLE, SAVE ::  x_out  !<
    132     COMPLEX(KIND=C_DOUBLE_COMPLEX), DIMENSION(:), ALLOCATABLE, SAVE ::         &
    133        y_out  !<
    134    
    135     REAL(KIND=C_DOUBLE), DIMENSION(:), ALLOCATABLE, SAVE ::                    &
    136        x_in   !<
    137     REAL(KIND=C_DOUBLE), DIMENSION(:), ALLOCATABLE, SAVE ::                    &
    138        y_in   !<
     134
     135    REAL(KIND=C_DOUBLE), DIMENSION(:), ALLOCATABLE, SAVE ::  x_in   !<
     136    REAL(KIND=C_DOUBLE), DIMENSION(:), ALLOCATABLE, SAVE ::  y_in   !<
    139137    !$OMP THREADPRIVATE( x_out, y_out, x_in, y_in )
    140    
    141    
     138
     139
    142140    TYPE(C_PTR), SAVE ::  plan_xf, plan_xi, plan_yf, plan_yi
    143141#endif
     
    176174
    177175
    178 !------------------------------------------------------------------------------!
     176!--------------------------------------------------------------------------------------------------!
    179177! Description:
    180178! ------------
    181179!> @todo Missing subroutine description.
    182 !------------------------------------------------------------------------------!
     180!--------------------------------------------------------------------------------------------------!
    183181    SUBROUTINE fft_init
    184182
     
    192190!--    in OpenMP sense
    193191#if defined( __ibm )
     192       REAL(wp), DIMENSION(nau2)   ::  aux2   !<
     193       REAL(wp), DIMENSION(nau2)   ::  auy2   !<
     194       REAL(wp), DIMENSION(nau2)   ::  aux4   !<
     195       REAL(wp), DIMENSION(nau2)   ::  auy4   !<
    194196       REAL(wp), DIMENSION(0:nx+2) ::  workx  !<
    195197       REAL(wp), DIMENSION(0:ny+2) ::  worky  !<
    196        REAL(wp), DIMENSION(nau2)   ::  aux2   !<
    197        REAL(wp), DIMENSION(nau2)   ::  auy2   !<
    198        REAL(wp), DIMENSION(nau2)   ::  aux4   !<
    199        REAL(wp), DIMENSION(nau2)   ::  auy4   !<
    200198#elif defined( __nec_fft )
    201199       REAL(wp), DIMENSION(0:nx+3,nz+1)   ::  work_x  !<
     
    203201       REAL(wp), DIMENSION(6*(nx+3),nz+1) ::  workx   !<
    204202       REAL(wp), DIMENSION(6*(ny+3),nz+1) ::  worky   !<
    205 #endif 
     203#endif
    206204
    207205!
     
    233231!
    234232!--       Initialize tables for fft along x
    235           CALL DRCFT( 1, workx, 1, workx, 1, nx+1, 1,  1, sqr_dnx, aux1, nau1, &
    236                       aux2, nau2 )
    237           CALL DCRFT( 1, workx, 1, workx, 1, nx+1, 1, -1, sqr_dnx, aux3, nau1, &
    238                       aux4, nau2 )
     233          CALL DRCFT( 1, workx, 1, workx, 1, nx+1, 1,  1, sqr_dnx, aux1, nau1, aux2, nau2 )
     234          CALL DCRFT( 1, workx, 1, workx, 1, nx+1, 1, -1, sqr_dnx, aux3, nau1, aux4, nau2 )
    239235!
    240236!--       Initialize tables for fft along y
    241           CALL DRCFT( 1, worky, 1, worky, 1, ny+1, 1,  1, sqr_dny, auy1, nau1, &
    242                       auy2, nau2 )
    243           CALL DCRFT( 1, worky, 1, worky, 1, ny+1, 1, -1, sqr_dny, auy3, nau1, &
    244                       auy4, nau2 )
     237          CALL DRCFT( 1, worky, 1, worky, 1, ny+1, 1,  1, sqr_dny, auy1, nau1, auy2, nau2 )
     238          CALL DCRFT( 1, worky, 1, worky, 1, ny+1, 1, -1, sqr_dny, auy3, nau1, auy4, nau2 )
    245239#elif defined( __nec_fft )
    246           message_string = 'fft method "' // TRIM( fft_method) // &
    247                            '" currently does not work on NEC'
     240          message_string = 'fft method "' // TRIM( fft_method) // '" currently does not work on NEC'
    248241          CALL message( 'fft_init', 'PA0187', 1, 2, 0, 6, 0 )
    249242
    250           ALLOCATE( trig_xb(2*(nx+1)), trig_xf(2*(nx+1)),                      &
    251                     trig_yb(2*(ny+1)), trig_yf(2*(ny+1)) )
     243          ALLOCATE( trig_xb(2*(nx+1)), trig_xf(2*(nx+1)), trig_yb(2*(ny+1)), trig_yf(2*(ny+1)) )
    252244
    253245          work_x = 0.0_wp
     
    260252          CALL DZFFT( 0, nx+1, sqr_dnx, work_x, work_x, trig_xf, workx, 0 )
    261253          CALL ZDFFT( 0, nx+1, sqr_dnx, work_x, work_x, trig_xb, workx, 0 )
    262           CALL DZFFTM( 0, nx+1, nz1, sqr_dnx, work_x, nx+4, work_x, nx+4,      &
    263                        trig_xf, workx, 0 )
    264           CALL ZDFFTM( 0, nx+1, nz1, sqr_dnx, work_x, nx+4, work_x, nx+4,      &
    265                        trig_xb, workx, 0 )
     254          CALL DZFFTM( 0, nx+1, nz1, sqr_dnx, work_x, nx+4, work_x, nx+4, trig_xf, workx, 0 )
     255          CALL ZDFFTM( 0, nx+1, nz1, sqr_dnx, work_x, nx+4, work_x, nx+4, trig_xb, workx, 0 )
    266256!
    267257!--       Initialize tables for fft along y (non-vector and vector case (M))
    268258          CALL DZFFT( 0, ny+1, sqr_dny, work_y, work_y, trig_yf, worky, 0 )
    269259          CALL ZDFFT( 0, ny+1, sqr_dny, work_y, work_y, trig_yb, worky, 0 )
    270           CALL DZFFTM( 0, ny+1, nz1, sqr_dny, work_y, ny+4, work_y, ny+4,      &
    271                        trig_yf, worky, 0 )
    272           CALL ZDFFTM( 0, ny+1, nz1, sqr_dny, work_y, ny+4, work_y, ny+4,      &
    273                        trig_yb, worky, 0 )
     260          CALL DZFFTM( 0, ny+1, nz1, sqr_dny, work_y, ny+4, work_y, ny+4, trig_yf, worky, 0 )
     261          CALL ZDFFTM( 0, ny+1, nz1, sqr_dny, work_y, ny+4, work_y, ny+4, trig_yb, worky, 0 )
    274262#elif defined( __cuda_fft )
    275263          CALL CUFFTPLAN1D( plan_xf, nx+1, CUFFT_D2Z, (nyn_x-nys_x+1) * (nzt_x-nzb_x+1) )
     
    303291          ny_c = ny+1
    304292          !$OMP PARALLEL
    305           ALLOCATE( x_in(0:nx+2), y_in(0:ny+2), x_out(0:(nx+1)/2),             &
    306                     y_out(0:(ny+1)/2) )
     293          ALLOCATE( x_in(0:nx+2), y_in(0:ny+2), x_out(0:(nx+1)/2), y_out(0:(ny+1)/2) )
    307294          !$OMP END PARALLEL
    308295          plan_xf = FFTW_PLAN_DFT_R2C_1D( nx_c, x_in, x_out, FFTW_ESTIMATE )
     
    321308       ELSE
    322309
    323           message_string = 'fft method "' // TRIM( fft_method) // &
    324                            '" not available'
     310          message_string = 'fft method "' // TRIM( fft_method) // '" not available'
    325311          CALL message( 'fft_init', 'PA0189', 1, 2, 0, 6, 0 )
    326312       ENDIF
     
    329315
    330316
    331 !------------------------------------------------------------------------------!
     317!--------------------------------------------------------------------------------------------------!
    332318! Description:
    333319! ------------
    334 !> Fourier-transformation along x-direction.                 
     320!> Fourier-transformation along x-direction.
    335321!> Version for 2D-decomposition.
    336 !> It uses internal algorithms (Singleton or Temperton) or      
    337 !> system-specific routines, if they are available           
    338 !------------------------------------------------------------------------------!
    339  
     322!> It uses internal algorithms (Singleton or Temperton) or system-specific routines, if they are
     323!> available.
     324!--------------------------------------------------------------------------------------------------!
     325
    340326    SUBROUTINE fft_x( ar, direction, ar_2d, ar_inv )
    341327
     
    344330
    345331       CHARACTER (LEN=*) ::  direction  !<
    346        
     332
    347333       COMPLEX(wp), DIMENSION(:), ALLOCATABLE ::  cwork  !<
    348334
    349        INTEGER(iwp) ::  i          !< 
     335       INTEGER(iwp) ::  i          !<
    350336       INTEGER(iwp) ::  ishape(1)  !<
    351337       INTEGER(iwp) ::  j          !<
     
    354340
    355341       LOGICAL ::  forward_fft !<
    356        
     342
    357343       REAL(wp), DIMENSION(0:nx+2) ::  work   !<
    358344       REAL(wp), DIMENSION(nx+2)   ::  work1  !<
    359        
     345
    360346       REAL(wp), DIMENSION(:,:), ALLOCATABLE           ::  work_vec  !<
    361347       REAL(wp), DIMENSION(0:nx,nys_x:nyn_x), OPTIONAL ::  ar_2d     !<
    362348
     349       REAL(wp), DIMENSION(0:nx,nys_x:nyn_x,nzb_x:nzt_x)           ::  ar       !<
    363350       REAL(wp), DIMENSION(nys_x:nyn_x,nzb_x:nzt_x,0:nx), OPTIONAL ::  ar_inv   !<
    364        REAL(wp), DIMENSION(0:nx,nys_x:nyn_x,nzb_x:nzt_x)           ::  ar       !<
    365351
    366352#if defined( __ibm )
    367        REAL(wp), DIMENSION(nau2) ::  aux2  !< 
     353       REAL(wp), DIMENSION(nau2) ::  aux2  !<
    368354       REAL(wp), DIMENSION(nau2) ::  aux4  !<
    369355#elif defined( __nec_fft )
     
    387373
    388374!
    389 !--       Performing the fft with singleton's software works on every system,
    390 !--       since it is part of the model
     375!--       Performing the fft with singleton's software works on every system, since it is part of
     376!--       the model.
    391377          ALLOCATE( cwork(0:nx) )
    392      
    393           IF ( forward_fft )   then
     378
     379          IF ( forward_fft )  THEN
    394380
    395381             !$OMP PARALLEL PRIVATE ( cwork, i, ishape, j, k )
     
    425411                   cwork(0) = CMPLX( ar(0,j,k), 0.0_wp, KIND=wp )
    426412                   DO  i = 1, (nx+1)/2 - 1
    427                       cwork(i)      = CMPLX( ar(i,j,k), -ar(nx+1-i,j,k),       &
    428                                              KIND=wp )
    429                       cwork(nx+1-i) = CMPLX( ar(i,j,k),  ar(nx+1-i,j,k),       &
    430                                              KIND=wp )
     413                      cwork(i)      = CMPLX( ar(i,j,k), -ar(nx+1-i,j,k), KIND=wp )
     414                      cwork(nx+1-i) = CMPLX( ar(i,j,k),  ar(nx+1-i,j,k), KIND=wp )
    431415                   ENDDO
    432416                   cwork((nx+1)/2) = CMPLX( ar((nx+1)/2,j,k), 0.0_wp, KIND=wp )
     
    450434
    451435!
    452 !--       Performing the fft with Temperton's software works on every system,
    453 !--       since it is part of the model
     436!--       Performing the fft with Temperton's software works on every system, since it is part of
     437!--       the model.
    454438          IF ( forward_fft )  THEN
    455439
     
    633617                      x_out(0) = CMPLX( ar_2d(0,j), 0.0_wp, KIND=wp )
    634618                      DO  i = 1, (nx+1)/2 - 1
    635                          x_out(i) = CMPLX( ar_2d(i,j), ar_2d(nx+1-i,j),        &
    636                                            KIND=wp )
    637                       ENDDO
    638                       x_out((nx+1)/2) = CMPLX( ar_2d((nx+1)/2,j), 0.0_wp,      &
    639                                                KIND=wp )
     619                         x_out(i) = CMPLX( ar_2d(i,j), ar_2d(nx+1-i,j), KIND=wp )
     620                      ENDDO
     621                      x_out((nx+1)/2) = CMPLX( ar_2d((nx+1)/2,j), 0.0_wp, KIND=wp )
    640622
    641623                   ELSE
     
    645627                         x_out(i) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k), KIND=wp )
    646628                      ENDDO
    647                       x_out((nx+1)/2) = CMPLX( ar((nx+1)/2,j,k), 0.0_wp,       &
    648                                                KIND=wp )
     629                      x_out((nx+1)/2) = CMPLX( ar((nx+1)/2,j,k), 0.0_wp, KIND=wp )
    649630
    650631                   ENDIF
     
    670651                DO  j = nys_x, nyn_x
    671652
    672                    CALL DRCFT( 0, ar, 1, work, 1, nx+1, 1, 1, sqr_dnx, aux1,   &
    673                                nau1, aux2, nau2 )
     653                   CALL DRCFT( 0, ar, 1, work, 1, nx+1, 1, 1, sqr_dnx, aux1, nau1, aux2, nau2 )
    674654
    675655                   DO  i = 0, (nx+1)/2
     
    700680                   work(nx+2) = 0.0_wp
    701681
    702                    CALL DCRFT( 0, work, 1, work, 1, nx+1, 1, -1, sqr_dnx,      &
    703                                aux3, nau1, aux4, nau2 )
     682                   CALL DCRFT( 0, work, 1, work, 1, nx+1, 1, -1, sqr_dnx, aux3, nau1, aux4, nau2 )
    704683
    705684                   DO  i = 0, nx
     
    725704
    726705                   CALL DZFFT( 1, nx+1, sqr_dnx, work, work, trig_xf, work2, 0 )
    727      
     706
    728707                   DO  i = 0, (nx+1)/2
    729708                      ar(i,j,k) = work(2*i)
     
    797776
    798777                   DO  i = 1, (nx+1)/2 - 1
    799                       ar_tmp(i,j,k) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k),        &
    800                                              KIND=wp )
    801                    ENDDO
    802                    ar_tmp((nx+1)/2,j,k) = CMPLX( ar((nx+1)/2,j,k), 0.0_wp,     &
    803                                                  KIND=wp )
     778                      ar_tmp(i,j,k) = CMPLX( ar(i,j,k), ar(nx+1-i,j,k), KIND=wp )
     779                   ENDDO
     780                   ar_tmp((nx+1)/2,j,k) = CMPLX( ar((nx+1)/2,j,k), 0.0_wp, KIND=wp )
    804781
    805782                ENDDO
     
    818795    END SUBROUTINE fft_x
    819796
    820 !------------------------------------------------------------------------------!
     797!--------------------------------------------------------------------------------------------------!
    821798! Description:
    822799! ------------
    823800!> Fourier-transformation along x-direction.
    824801!> Version for 1D-decomposition.
    825 !> It uses internal algorithms (Singleton or Temperton) or
    826 !> system-specific routines, if they are available
    827 !------------------------------------------------------------------------------!
    828  
     802!> It uses internal algorithms (Singleton or Temperton) or system-specific routines, if they are
     803!> available.
     804!--------------------------------------------------------------------------------------------------!
     805
    829806    SUBROUTINE fft_x_1d( ar, direction )
    830807
     
    833810
    834811       CHARACTER (LEN=*) ::  direction  !<
    835        
     812
    836813       INTEGER(iwp) ::  i               !<
    837814       INTEGER(iwp) ::  ishape(1)       !<
     
    842819       REAL(wp), DIMENSION(0:nx+2) ::  work   !<
    843820       REAL(wp), DIMENSION(nx+2)   ::  work1  !<
    844        
     821
    845822       COMPLEX(wp), DIMENSION(:), ALLOCATABLE ::  cwork  !<
    846        
     823
    847824#if defined( __ibm )
    848825       REAL(wp), DIMENSION(nau2) ::  aux2       !<
     
    861838
    862839!
    863 !--       Performing the fft with singleton's software works on every system,
    864 !--       since it is part of the model
     840!--       Performing the fft with singleton's software works on every system, since it is part of
     841!--       the model.
    865842          ALLOCATE( cwork(0:nx) )
    866      
    867           IF ( forward_fft )   then
     843
     844          IF ( forward_fft )  THEN
    868845
    869846             DO  i = 0, nx
     
    902879
    903880!
    904 !--       Performing the fft with Temperton's software works on every system,
    905 !--       since it is part of the model
     881!--       Performing the fft with Temperton's software works on every system, since it is part of
     882!--       the model.
    906883          IF ( forward_fft )  THEN
    907884
     
    966943          IF ( forward_fft )  THEN
    967944
    968              CALL DRCFT( 0, ar, 1, work, 1, nx+1, 1, 1, sqr_dnx, aux1, nau1,   &
    969                          aux2, nau2 )
     945             CALL DRCFT( 0, ar, 1, work, 1, nx+1, 1, 1, sqr_dnx, aux1, nau1, aux2, nau2 )
    970946
    971947             DO  i = 0, (nx+1)/2
     
    987963             work(nx+2) = 0.0_wp
    988964
    989              CALL DCRFT( 0, work, 1, work, 1, nx+1, 1, -1, sqr_dnx, aux3, nau1, &
    990                          aux4, nau2 )
     965             CALL DCRFT( 0, work, 1, work, 1, nx+1, 1, -1, sqr_dnx, aux3, nau1, aux4, nau2 )
    991966
    992967             DO  i = 0, nx
     
    1001976
    1002977             CALL DZFFT( 1, nx+1, sqr_dnx, work, work, trig_xf, work2, 0 )
    1003      
     978
    1004979             DO  i = 0, (nx+1)/2
    1005980                ar(i) = work(2*i)
     
    10311006    END SUBROUTINE fft_x_1d
    10321007
    1033 !------------------------------------------------------------------------------!
     1008!--------------------------------------------------------------------------------------------------!
    10341009! Description:
    10351010! ------------
    10361011!> Fourier-transformation along y-direction.
    10371012!> Version for 2D-decomposition.
    1038 !> It uses internal algorithms (Singleton or Temperton) or
    1039 !> system-specific routines, if they are available.
    1040 !> 
     1013!> It uses internal algorithms (Singleton or Temperton) or system-specific routines, if they are
     1014!> available.
     1015!>
    10411016!> direction:  'forward' or 'backward'
    1042 !> ar, ar_tr:  3D data arrays 
     1017!> ar, ar_tr:  3D data arrays
    10431018!>             forward:   ar: before  ar_tr: after transformation
    10441019!>             backward:  ar_tr: before  ar: after transfosition
    1045 !> 
     1020!>
    10461021!> In case of non-overlapping transposition/transformation:
    1047 !> nxl_y_bound = nxl_y_l = nxl_y 
    1048 !> nxr_y_bound = nxr_y_l = nxr_y 
    1049 !> 
     1022!> nxl_y_bound = nxl_y_l = nxl_y
     1023!> nxr_y_bound = nxr_y_l = nxr_y
     1024!>
    10501025!> In case of overlapping transposition/transformation
    1051 !> - nxl_y_bound  and  nxr_y_bound have the original values of
    1052 !>   nxl_y, nxr_y.  ar_tr is dimensioned using these values.
    1053 !> - nxl_y_l = nxr_y_r.  ar is dimensioned with these values, so that
    1054 !>   transformation is carried out for a 2D-plane only.
    1055 !------------------------------------------------------------------------------!
    1056  
    1057     SUBROUTINE fft_y( ar, direction, ar_tr, nxl_y_bound, nxr_y_bound, nxl_y_l, &
    1058                       nxr_y_l, ar_inv )
     1026!> - nxl_y_bound  and  nxr_y_bound have the original values of nxl_y, nxr_y.  ar_tr is dimensioned
     1027!>   using these values.
     1028!> - nxl_y_l = nxr_y_r.  ar is dimensioned with these values, so that transformation is carried out
     1029!>   for a 2D-plane only.
     1030!--------------------------------------------------------------------------------------------------!
     1031
     1032    SUBROUTINE fft_y( ar, direction, ar_tr, nxl_y_bound, nxr_y_bound, nxl_y_l, nxr_y_l, ar_inv )
    10591033
    10601034
     
    10621036
    10631037       CHARACTER (LEN=*) ::  direction  !<
    1064        
     1038
    10651039       INTEGER(iwp) ::  i            !<
    1066        INTEGER(iwp) ::  j            !< 
     1040       INTEGER(iwp) ::  j            !<
    10671041       INTEGER(iwp) ::  jshape(1)    !<
    10681042       INTEGER(iwp) ::  k            !<
     
    10771051       REAL(wp), DIMENSION(0:ny+2) ::  work   !<
    10781052       REAL(wp), DIMENSION(ny+2)   ::  work1  !<
    1079        
     1053
    10801054       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  f_vec_y
    10811055       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  work_vec
     
    10861060
    10871061       COMPLEX(wp), DIMENSION(:), ALLOCATABLE ::  cwork  !<
    1088        
     1062
    10891063#if defined( __ibm )
    10901064       REAL(wp), DIMENSION(nau2) ::  auy2  !<
     
    10931067       REAL(wp), DIMENSION(6*(ny+1)) ::  work2  !<
    10941068#elif defined( __cuda_fft )
    1095        COMPLEX(dp), DIMENSION(0:(ny+1)/2,nxl_y:nxr_y,nzb_y:nzt_y) ::           &
    1096           ar_tmp  !<
     1069       COMPLEX(dp), DIMENSION(0:(ny+1)/2,nxl_y:nxr_y,nzb_y:nzt_y) ::  ar_tmp  !<
    10971070       !$ACC DECLARE CREATE(ar_tmp)
    10981071#endif
     
    11081081
    11091082!
    1110 !--       Performing the fft with singleton's software works on every system,
    1111 !--       since it is part of the model
     1083!--       Performing the fft with singleton's software works on every system, since it is part of
     1084!--       the model.
    11121085          ALLOCATE( cwork(0:ny) )
    11131086
    1114           IF ( forward_fft )   then
     1087          IF ( forward_fft )  THEN
    11151088
    11161089             !$OMP PARALLEL PRIVATE ( cwork, i, jshape, j, k )
     
    11461119                   cwork(0) = CMPLX( ar_tr(0,i,k), 0.0_wp, KIND=wp )
    11471120                   DO  j = 1, (ny+1)/2 - 1
    1148                       cwork(j)      = CMPLX( ar_tr(j,i,k), -ar_tr(ny+1-j,i,k), &
    1149                                              KIND=wp )
    1150                       cwork(ny+1-j) = CMPLX( ar_tr(j,i,k),  ar_tr(ny+1-j,i,k), &
    1151                                              KIND=wp )
    1152                    ENDDO
    1153                    cwork((ny+1)/2) = CMPLX( ar_tr((ny+1)/2,i,k), 0.0_wp,       &
    1154                                             KIND=wp )
     1121                      cwork(j)      = CMPLX( ar_tr(j,i,k), -ar_tr(ny+1-j,i,k), KIND=wp )
     1122                      cwork(ny+1-j) = CMPLX( ar_tr(j,i,k),  ar_tr(ny+1-j,i,k), KIND=wp )
     1123                   ENDDO
     1124                   cwork((ny+1)/2) = CMPLX( ar_tr((ny+1)/2,i,k), 0.0_wp, KIND=wp )
    11551125
    11561126                   jshape = SHAPE( cwork )
     
    11721142
    11731143!
    1174 !--       Performing the fft with Temperton's software works on every system,
    1175 !--       since it is part of the model
     1144!--       Performing the fft with Temperton's software works on every system, since it is part of
     1145!--       the model.
    11761146          IF ( forward_fft )  THEN
    11771147
     
    13611331                   y_out(0) = CMPLX( ar_tr(0,i,k), 0.0_wp, KIND=wp )
    13621332                   DO  j = 1, (ny+1)/2 - 1
    1363                       y_out(j) = CMPLX( ar_tr(j,i,k), ar_tr(ny+1-j,i,k),       &
    1364                                         KIND=wp )
    1365                    ENDDO
    1366                    y_out((ny+1)/2) = CMPLX( ar_tr((ny+1)/2,i,k), 0.0_wp,       &
    1367                                             KIND=wp )
     1333                      y_out(j) = CMPLX( ar_tr(j,i,k), ar_tr(ny+1-j,i,k), KIND=wp )
     1334                   ENDDO
     1335                   y_out((ny+1)/2) = CMPLX( ar_tr((ny+1)/2,i,k), 0.0_wp, KIND=wp )
    13681336
    13691337                   CALL FFTW_EXECUTE_DFT_C2R( plan_yi, y_out, y_in )
     
    13871355                DO  i = nxl_y_l, nxr_y_l
    13881356
    1389                    CALL DRCFT( 0, ar, 1, work, 1, ny+1, 1, 1, sqr_dny, auy1,   &
    1390                                nau1, auy2, nau2 )
     1357                   CALL DRCFT( 0, ar, 1, work, 1, ny+1, 1, 1, sqr_dny, auy1, nau1, auy2, nau2 )
    13911358
    13921359                   DO  j = 0, (ny+1)/2
     
    14171384                   work(ny+2) = 0.0_wp
    14181385
    1419                    CALL DCRFT( 0, work, 1, work, 1, ny+1, 1, -1, sqr_dny,      &
    1420                                auy3, nau1, auy4, nau2 )
     1386                   CALL DCRFT( 0, work, 1, work, 1, ny+1, 1, -1, sqr_dny, auy3, nau1, auy4, nau2 )
    14211387
    14221388                   DO  j = 0, ny
     
    14911457
    14921458                   DO  j = 0, (ny+1)/2
    1493                       ar(j,i,k)      = REAL( ar_tmp(j,i,k), KIND=wp )  * dny
     1459                      ar(j,i,k)      = REAL( ar_tmp(j,i,k), KIND=wp ) * dny
    14941460                   ENDDO
    14951461
     
    15111477
    15121478                   DO  j = 1, (ny+1)/2 - 1
    1513                       ar_tmp(j,i,k) = CMPLX( ar(j,i,k), ar(ny+1-j,i,k),        &
    1514                                              KIND=wp )
    1515                    ENDDO
    1516                    ar_tmp((ny+1)/2,i,k) = CMPLX( ar((ny+1)/2,i,k), 0.0_wp,     &
    1517                                                  KIND=wp )
     1479                      ar_tmp(j,i,k) = CMPLX( ar(j,i,k), ar(ny+1-j,i,k), KIND=wp )
     1480                   ENDDO
     1481                   ar_tmp((ny+1)/2,i,k) = CMPLX( ar((ny+1)/2,i,k), 0.0_wp, KIND=wp )
    15181482
    15191483                ENDDO
     
    15321496    END SUBROUTINE fft_y
    15331497
    1534 !------------------------------------------------------------------------------!
     1498!--------------------------------------------------------------------------------------------------!
    15351499! Description:
    15361500! ------------
    15371501!> Fourier-transformation along y-direction.
    15381502!> Version for 1D-decomposition.
    1539 !> It uses internal algorithms (Singleton or Temperton) or
    1540 !> system-specific routines, if they are available.
    1541 !------------------------------------------------------------------------------!
    1542  
     1503!> It uses internal algorithms (Singleton or Temperton) or system-specific routines, if they are
     1504!> available.
     1505!--------------------------------------------------------------------------------------------------!
     1506
    15431507    SUBROUTINE fft_y_1d( ar, direction )
    15441508
     
    15471511
    15481512       CHARACTER (LEN=*) ::  direction
    1549        
     1513
    15501514       INTEGER(iwp) ::  j          !<
    15511515       INTEGER(iwp) ::  jshape(1)  !<
     
    15561520       REAL(wp), DIMENSION(0:ny+2)  ::  work   !<
    15571521       REAL(wp), DIMENSION(ny+2)    ::  work1  !<
    1558        
     1522
    15591523       COMPLEX(wp), DIMENSION(:), ALLOCATABLE ::  cwork  !<
    1560        
     1524
    15611525#if defined( __ibm )
    15621526       REAL(wp), DIMENSION(nau2) ::  auy2  !<
     
    15751539
    15761540!
    1577 !--       Performing the fft with singleton's software works on every system,
    1578 !--       since it is part of the model
     1541!--       Performing the fft with singleton's software works on every system, since it is part of
     1542!--       the model.
    15791543          ALLOCATE( cwork(0:ny) )
    15801544
     
    16181582
    16191583!
    1620 !--       Performing the fft with Temperton's software works on every system,
    1621 !--       since it is part of the model
     1584!--       Performing the fft with Temperton's software works on every system, since it is part of
     1585!--       the model.
    16221586          IF ( forward_fft )  THEN
    16231587
     
    16821646          IF ( forward_fft )  THEN
    16831647
    1684              CALL DRCFT( 0, ar, 1, work, 1, ny+1, 1, 1, sqr_dny, auy1, nau1,   &
    1685                          auy2, nau2 )
     1648             CALL DRCFT( 0, ar, 1, work, 1, ny+1, 1, 1, sqr_dny, auy1, nau1, auy2, nau2 )
    16861649
    16871650             DO  j = 0, (ny+1)/2
     
    17031666             work(ny+2) = 0.0_wp
    17041667
    1705              CALL DCRFT( 0, work, 1, work, 1, ny+1, 1, -1, sqr_dny, auy3,      &
    1706                          nau1, auy4, nau2 )
     1668             CALL DCRFT( 0, work, 1, work, 1, ny+1, 1, -1, sqr_dny, auy3, nau1, auy4, nau2 )
    17071669
    17081670             DO  j = 0, ny
     
    17471709    END SUBROUTINE fft_y_1d
    17481710
    1749 !------------------------------------------------------------------------------!
     1711!--------------------------------------------------------------------------------------------------!
    17501712! Description:
    17511713! ------------
    17521714!> Fourier-transformation along x-direction.
    1753 !> Version for 1d domain decomposition
     1715!> Version for 1d domain decomposition,
    17541716!> using multiple 1D FFT from Math Keisan on NEC or Temperton-algorithm
    1755 !> (no singleton-algorithm on NEC because it does not vectorize)
    1756 !------------------------------------------------------------------------------!
    1757  
     1717!> (no singleton-algorithm on NEC because it does not vectorize).
     1718!--------------------------------------------------------------------------------------------------!
     1719
    17581720    SUBROUTINE fft_x_m( ar, direction )
    17591721
     
    17621724
    17631725       CHARACTER (LEN=*) ::  direction  !<
    1764        
     1726
    17651727       INTEGER(iwp) ::  i     !<
    17661728       INTEGER(iwp) ::  k     !<
     
    17731735       REAL(wp), DIMENSION(0:nx+3,nz+1)   ::  ai     !<
    17741736       REAL(wp), DIMENSION(6*(nx+4),nz+1) ::  work1  !<
    1775        
     1737
    17761738#if defined( __nec_fft )
    17771739       COMPLEX(wp), DIMENSION(:,:), ALLOCATABLE ::  work
     
    18271789
    18281790!
    1829 !--          Tables are initialized once more. This call should not be
    1830 !--          necessary, but otherwise program aborts in asymmetric case
    1831              CALL DZFFTM( 0, nx+1, nz1, sqr_dnx, work, nx+4, work, nx+4,       &
    1832                           trig_xf, work1, 0 )
     1791!--          Tables are initialized once more. This call should not be necessary, but otherwise
     1792!--          program aborts in asymmetric case.
     1793             CALL DZFFTM( 0, nx+1, nz1, sqr_dnx, work, nx+4, work, nx+4, trig_xf, work1, 0 )
    18331794
    18341795             ai(0:nx,1:nz) = ar(0:nx,1:nz)
     
    18371798             ENDIF
    18381799
    1839              CALL DZFFTM( 1, nx+1, nz1, sqr_dnx, ai, siza, work, sizw,         &
    1840                           trig_xf, work1, 0 )
     1800             CALL DZFFTM( 1, nx+1, nz1, sqr_dnx, ai, siza, work, sizw, trig_xf, work1, 0 )
    18411801
    18421802             DO  k = 1, nz
     
    18541814!--          Tables are initialized once more. This call should not be
    18551815!--          necessary, but otherwise program aborts in asymmetric case
    1856              CALL ZDFFTM( 0, nx+1, nz1, sqr_dnx, work, nx+4, work, nx+4,       &
    1857                           trig_xb, work1, 0 )
     1816             CALL ZDFFTM( 0, nx+1, nz1, sqr_dnx, work, nx+4, work, nx+4, trig_xb, work1, 0 )
    18581817
    18591818             IF ( nz1 > nz )  THEN
     
    18681827             ENDDO
    18691828
    1870              CALL ZDFFTM( -1, nx+1, nz1, sqr_dnx, work, sizw, ai, siza, &
    1871                           trig_xb, work1, 0 )
     1829             CALL ZDFFTM( -1, nx+1, nz1, sqr_dnx, work, sizw, ai, siza, trig_xb, work1, 0 )
    18721830
    18731831             ar(0:nx,1:nz) = ai(0:nx,1:nz)
     
    18821840    END SUBROUTINE fft_x_m
    18831841
    1884 !------------------------------------------------------------------------------!
     1842!--------------------------------------------------------------------------------------------------!
    18851843! Description:
    18861844! ------------
    18871845!> Fourier-transformation along y-direction.
    1888 !> Version for 1d domain decomposition
     1846!> Version for 1d domain decomposition,
    18891847!> using multiple 1D FFT from Math Keisan on NEC or Temperton-algorithm
    1890 !> (no singleton-algorithm on NEC because it does not vectorize)
    1891 !------------------------------------------------------------------------------!
    1892  
     1848!> (no singleton-algorithm on NEC because it does not vectorize).
     1849!--------------------------------------------------------------------------------------------------!
     1850
    18931851    SUBROUTINE fft_y_m( ar, ny1, direction )
    18941852
     
    18971855
    18981856       CHARACTER (LEN=*) ::  direction  !<
    1899        
    1900        INTEGER(iwp) ::  j     !< 
     1857
     1858       INTEGER(iwp) ::  j     !<
    19011859       INTEGER(iwp) ::  k     !<
    19021860       INTEGER(iwp) ::  ny1   !<
     
    19641922
    19651923!
    1966 !--          Tables are initialized once more. This call should not be
    1967 !--          necessary, but otherwise program aborts in asymmetric case
    1968              CALL DZFFTM( 0, ny+1, nz1, sqr_dny, work, ny+4, work, ny+4, &
    1969                           trig_yf, work1, 0 )
     1924!--          Tables are initialized once more. This call should not be necessary, but otherwise
     1925!--          program aborts in asymmetric case.
     1926             CALL DZFFTM( 0, ny+1, nz1, sqr_dny, work, ny+4, work, ny+4, trig_yf, work1, 0 )
    19701927
    19711928             ai(0:ny,1:nz) = ar(0:ny,1:nz)
     
    19741931             ENDIF
    19751932
    1976              CALL DZFFTM( 1, ny+1, nz1, sqr_dny, ai, siza, work, sizw, &
    1977                           trig_yf, work1, 0 )
     1933             CALL DZFFTM( 1, ny+1, nz1, sqr_dny, ai, siza, work, sizw, trig_yf, work1, 0 )
    19781934
    19791935             DO  k = 1, nz
     
    19891945
    19901946!
    1991 !--          Tables are initialized once more. This call should not be
    1992 !--          necessary, but otherwise program aborts in asymmetric case
    1993              CALL ZDFFTM( 0, ny+1, nz1, sqr_dny, work, ny+4, work, ny+4, &
    1994                           trig_yb, work1, 0 )
     1947!--          Tables are initialized once more. This call should not be necessary, but otherwise
     1948!--          program aborts in asymmetric case.
     1949             CALL ZDFFTM( 0, ny+1, nz1, sqr_dny, work, ny+4, work, ny+4, trig_yb, work1, 0 )
    19951950
    19961951             IF ( nz1 > nz )  THEN
     
    20051960             ENDDO
    20061961
    2007              CALL ZDFFTM( -1, ny+1, nz1, sqr_dny, work, sizw, ai, siza, &
    2008                           trig_yb, work1, 0 )
     1962             CALL ZDFFTM( -1, ny+1, nz1, sqr_dny, work, sizw, ai, siza, trig_yb, work1, 0 )
    20091963
    20101964             ar(0:ny,1:nz) = ai(0:ny,1:nz)
  • palm/trunk/SOURCE/flow_statistics.f90

    r4581 r4646  
    11!> @file flow_statistics.f90
    2 !------------------------------------------------------------------------------!
     2!--------------------------------------------------------------------------------------------------!
    33! This file is part of the PALM model system.
    44!
    5 ! PALM is free software: you can redistribute it and/or modify it under the
    6 ! terms of the GNU General Public License as published by the Free Software
    7 ! Foundation, either version 3 of the License, or (at your option) any later
    8 ! version.
    9 !
    10 ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
    11 ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
    12 ! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
    13 !
    14 ! You should have received a copy of the GNU General Public License along with
    15 ! PALM. If not, see <http://www.gnu.org/licenses/>.
     5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
     6! Public License as published by the Free Software Foundation, either version 3 of the License, or
     7! (at your option) any later version.
     8!
     9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
     10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
     11! Public License for more details.
     12!
     13! You should have received a copy of the GNU General Public License along with PALM. If not, see
     14! <http://www.gnu.org/licenses/>.
    1615!
    1716! Copyright 1997-2020 Leibniz Universitaet Hannover
    18 !------------------------------------------------------------------------------!
     17!--------------------------------------------------------------------------------------------------!
    1918!
    2019! Current revisions:
    2120! ------------------
    22 ! 
    23 ! 
     21!
     22!
    2423! Former revisions:
    2524! -----------------
    2625! $Id$
     26! file re-formatted to follow the PALM coding standard
     27!
     28! 4581 2020-06-29 08:49:58Z suehring
    2729! Formatting adjustment
    28 ! 
     30!
    2931! 4551 2020-06-02 10:22:25Z suehring
    3032! Bugfix in summation for statistical regions
    31 ! 
     33!
    3234! 4521 2020-05-06 11:39:49Z schwenkel
    3335! Rename variable
    34 ! 
     36!
    3537! 4502 2020-04-17 16:14:16Z schwenkel
    3638! Implementation of ice microphysics
    37 ! 
     39!
    3840! 4472 2020-03-24 12:21:00Z Giersch
    3941! Calculations of the Kolmogorov lengt scale eta implemented
     
    4345!
    4446! 4463 2020-03-17 09:27:36Z Giersch
    45 ! Calculate horizontally averaged profiles of all velocity components at the
    46 ! same place
     47! Calculate horizontally averaged profiles of all velocity components at the same place
    4748!
    4849! 4444 2020-03-05 15:59:50Z raasch
     
    5051!
    5152! 4442 2020-03-04 19:21:13Z suehring
    52 ! Change order of dimension in surface array %frac to allow for better
    53 ! vectorization.
     53! Change order of dimension in surface array %frac to allow for better vectorization.
    5454!
    5555! 4441 2020-03-04 19:20:35Z suehring
    56 ! Introduction of wall_flags_total_0, which currently sets bits based on static
    57 ! topography information used in wall_flags_static_0
     56! Introduction of wall_flags_total_0, which currently sets bits based on static topography
     57! information used in wall_flags_static_0
    5858!
    5959! 4329 2019-12-10 15:46:36Z motisi
     
    6767!
    6868! 4039 2019-06-18 10:32:41Z suehring
    69 ! Correct conversion to kinematic scalar fluxes in case of pw-scheme and
    70 ! statistic regions
     69! Correct conversion to kinematic scalar fluxes in case of pw-scheme and statistic regions
    7170!
    7271! 3828 2019-03-27 19:36:23Z raasch
     
    8281! Description:
    8382! ------------
    84 !> Compute average profiles and further average flow quantities for the different
    85 !> user-defined (sub-)regions. The region indexed 0 is the total model domain.
     83!> Compute average profiles and further average flow quantities for the different user-defined
     84!> (sub-)regions. The region indexed 0 is the total model domain.
    8685!>
    87 !> @note For simplicity, nzb_s_inner and nzb_diff_s_inner are being used as a
    88 !>       lower vertical index for k-loops for all variables, although strictly
    89 !>       speaking the k-loops would have to be split up according to the staggered
    90 !>       grid. However, this implies no error since staggered velocity components
    91 !>       are zero at the walls and inside buildings.
    92 !------------------------------------------------------------------------------!
     86!> @note For simplicity, nzb_s_inner and nzb_diff_s_inner are used as a lower vertical index for
     87!>       k-loops for all variables, although strictly speaking the k-loops would have to be split
     88!>       up according to the staggered grid. However, this implies no error since staggered velocity
     89!>       components are zero at the walls and inside buildings.
     90!--------------------------------------------------------------------------------------------------!
    9391 SUBROUTINE flow_statistics
    9492
    9593
    96     USE arrays_3d,                                                             &
    97         ONLY:  ddzu, ddzw, e, heatflux_output_conversion, hyp, km, kh,         &
    98                momentumflux_output_conversion, nc, ni, nr, p, prho, prr, pt, q,&
    99                qc, qi, ql, qr, rho_air, rho_air_zw, rho_ocean, s,              &
    100                sa, u, ug, v, vg, vpt, w, w_subs, waterflux_output_conversion,  &
    101                zw, d_exner
    102 
    103     USE basic_constants_and_equations_mod,                                     &
     94    USE arrays_3d,                                                                                 &
     95        ONLY:  ddzu, ddzw, d_exner, e, heatflux_output_conversion, hyp, km, kh,                    &
     96               momentumflux_output_conversion, nc, ni, nr, p, prho, prr, pt, q, qc, qi, ql, qr,    &
     97               rho_air, rho_air_zw, rho_ocean, s, sa, u, ug, v, vg, vpt, w, w_subs,                &
     98               waterflux_output_conversion, zw
     99
     100    USE basic_constants_and_equations_mod,                                                         &
    104101        ONLY:  g, lv_d_cp
    105102
    106     USE bulk_cloud_model_mod,                                                  &
    107         ONLY: bulk_cloud_model, microphysics_morrison, microphysics_seifert,   &
    108               microphysics_ice_phase
    109 
    110     USE chem_modules,                                                          &
     103    USE bulk_cloud_model_mod,                                                                      &
     104        ONLY: bulk_cloud_model, microphysics_morrison, microphysics_seifert, microphysics_ice_phase
     105
     106    USE chem_modules,                                                                              &
    111107        ONLY:  max_pr_cs
    112108
    113     USE control_parameters,                                                    &
    114         ONLY:   air_chemistry, average_count_pr, cloud_droplets, do_sum,       &
    115                 dt_3d, humidity, initializing_actions, kolmogorov_length_scale,&
    116                 land_surface, large_scale_forcing, large_scale_subsidence,     &
    117                 max_pr_user, message_string, neutral, ocean_mode,              &
    118                 passive_scalar, simulated_time, simulated_time_at_begin,       &
    119                 use_subsidence_tendencies, use_surface_fluxes, use_top_fluxes, &
    120                 ws_scheme_mom, ws_scheme_sca, salsa, max_pr_salsa
    121 
    122     USE cpulog,                                                                &
     109    USE control_parameters,                                                                        &
     110        ONLY:   air_chemistry, average_count_pr, cloud_droplets, do_sum, dt_3d, humidity,          &
     111                initializing_actions, kolmogorov_length_scale, land_surface, large_scale_forcing,  &
     112                large_scale_subsidence, max_pr_salsa, max_pr_user, message_string, neutral,        &
     113                ocean_mode, passive_scalar, salsa, simulated_time, simulated_time_at_begin,        &
     114                use_subsidence_tendencies, use_surface_fluxes, use_top_fluxes, ws_scheme_mom,      &
     115                ws_scheme_sca
     116
     117    USE cpulog,                                                                                    &
    123118        ONLY:   cpu_log, log_point
    124119
    125     USE grid_variables,                                                        &
     120    USE grid_variables,                                                                            &
    126121        ONLY:   ddx, ddy
    127122
    128     USE indices,                                                               &
    129         ONLY:   ngp_2dh, ngp_2dh_s_inner, ngp_3d, ngp_3d_inner, nxl, nxr, nyn, &
    130                 nys, nzb, nzt, topo_min_level, wall_flags_total_0
     123    USE indices,                                                                                   &
     124        ONLY:   ngp_2dh, ngp_2dh_s_inner, ngp_3d, ngp_3d_inner, nxl, nxr, nyn, nys, nzb, nzt,      &
     125                topo_min_level, wall_flags_total_0
    131126
    132127#if defined( __parallel )
    133     USE indices,                                                               &
     128    USE indices,                                                                                   &
    134129        ONLY:  ngp_sums, ngp_sums_ls
    135130#endif
     
    137132    USE kinds
    138133
    139     USE land_surface_model_mod,                                                &
     134    USE land_surface_model_mod,                                                                    &
    140135        ONLY:   m_soil_h, nzb_soil, nzt_soil, t_soil_h
    141136
    142     USE lsf_nudging_mod,                                                       &
     137    USE lsf_nudging_mod,                                                                           &
    143138        ONLY:   td_lsa_lpt, td_lsa_q, td_sub_lpt, td_sub_q, time_vert
    144139
    145     USE module_interface,                                                      &
     140    USE module_interface,                                                                          &
    146141        ONLY:  module_interface_statistics
    147142
    148     USE netcdf_interface,                                                      &
     143    USE netcdf_interface,                                                                          &
    149144        ONLY:  dots_rad, dots_soil, dots_max
    150145
    151146    USE pegrid
    152147
    153     USE radiation_model_mod,                                                   &
    154         ONLY:  radiation, radiation_scheme,                                    &
    155                rad_lw_in, rad_lw_out, rad_lw_cs_hr, rad_lw_hr,                 &
     148    USE radiation_model_mod,                                                                       &
     149        ONLY:  radiation, radiation_scheme,                                                        &
     150               rad_lw_in, rad_lw_out, rad_lw_cs_hr, rad_lw_hr,                                     &
    156151               rad_sw_in, rad_sw_out, rad_sw_cs_hr, rad_sw_hr
    157152
    158153    USE statistics
    159154
    160     USE surface_mod,                                                           &
     155    USE surface_mod,                                                                               &
    161156        ONLY :  surf_def_h, surf_lsm_h, surf_usm_h
    162157
     
    215210    REAL(wp) ::  sums_ll(nzb:nzt+1,2)  !<
    216211
     212
    217213    CALL cpu_log( log_point(10), 'flow_statistics', 'start' )
    218214
    219215
    220216!
    221 !-- To be on the safe side, check whether flow_statistics has already been
    222 !-- called once after the current time step
     217!-- To be on the safe side, check whether flow_statistics has already been called once after the
     218!-- current time step.
    223219    IF ( flow_statistics_called )  THEN
    224220
    225        message_string = 'flow_statistics is called two times within one ' // &
    226                         'timestep'
     221       message_string = 'flow_statistics is called two times within one ' // 'timestep'
    227222       CALL message( 'flow_statistics', 'PA0190', 1, 2, 0, 6, 0 )
    228223
     
    243238
    244239!
    245 !--    Store sums that have been computed in other subroutines in summation
    246 !--    array
     240!--    Store sums that have been computed in other subroutines in summation array
    247241       sums_l(:,11,:) = sums_l_l(:,sr,:)      ! mixing length from diffusivities
    248242!--    WARNING: next line still has to be adjusted for OpenMP
    249        sums_l(:,21,0) = sums_wsts_bc_l(:,sr) *                                 &
     243       sums_l(:,21,0) = sums_wsts_bc_l(:,sr) *                                                     &
    250244                        heatflux_output_conversion  ! heat flux from advec_s_bc
    251245       sums_l(nzb+9,pr_palm,0)  = sums_divold_l(sr)  ! old divergence from pres
     
    253247
    254248!
    255 !--    When calcuating horizontally-averaged total (resolved- plus subgrid-
    256 !--    scale) vertical fluxes and velocity variances by using commonly-
    257 !--    applied Reynolds-based methods ( e.g. <w'pt'> = (w-<w>)*(pt-<pt>) )
    258 !--    in combination with the 5th order advection scheme, pronounced
    259 !--    artificial kinks could be observed in the vertical profiles near the
    260 !--    surface. Please note: these kinks were not related to the model truth,
    261 !--    i.e. these kinks are just related to an evaluation problem.
    262 !--    In order avoid these kinks, vertical fluxes and horizontal as well
    263 !--    vertical velocity variances are calculated directly within the advection
    264 !--    routines, according to the numerical discretization, to evaluate the
    265 !--    statistical quantities as they will appear within the prognostic
    266 !--    equations.
    267 !--    Copy the turbulent quantities, evaluated in the advection routines to
    268 !--    the local array sums_l() for further computations.
     249!--    When calcuating horizontally-averaged total (resolved- plus subgrid-scale) vertical fluxes
     250!--    and velocity variances by using commonly-applied Reynolds-based methods
     251!--    ( e.g. <w'pt'> = (w-<w>)*(pt-<pt>) ) in combination with the 5th order advection scheme,
     252!--    pronounced artificial kinks could be observed in the vertical profiles near the surface.
     253!--    Please note: these kinks were not related to the model truth, i.e. these kinks are just
     254!--    related to an evaluation problem.
     255!--    In order avoid these kinks, vertical fluxes and horizontal as well vertical velocity
     256!--    variances are calculated directly within the advection routines, according to the numerical
     257!--    discretization, to evaluate the statistical quantities as they will appear within the
     258!--    prognostic equations.
     259!--    Copy the turbulent quantities, evaluated in the advection routines to the local array
     260!--    sums_l() for further computations.
    269261       IF ( ws_scheme_mom .AND. sr == 0 )  THEN
    270262
    271263!
    272 !--       According to the Neumann bc for the horizontal velocity components,
    273 !--       the corresponding fluxes has to satisfiy the same bc.
     264!--       According to the Neumann bc for the horizontal velocity components, the corresponding
     265!--       fluxes has to satisfiy the same bc.
    274266          IF ( ocean_mode )  THEN
    275267             sums_us2_ws_l(nzt+1,:) = sums_us2_ws_l(nzt,:)
     
    280272!
    281273!--          Swap the turbulent quantities evaluated in advec_ws.
    282              sums_l(:,13,i) = sums_wsus_ws_l(:,i)                              &
    283                               * momentumflux_output_conversion ! w*u*
    284              sums_l(:,15,i) = sums_wsvs_ws_l(:,i)                              &
    285                               * momentumflux_output_conversion ! w*v*
    286              sums_l(:,30,i) = sums_us2_ws_l(:,i)        ! u*2
    287              sums_l(:,31,i) = sums_vs2_ws_l(:,i)        ! v*2
    288              sums_l(:,32,i) = sums_ws2_ws_l(:,i)        ! w*2
    289              sums_l(:,34,i) = sums_l(:,34,i) + 0.5_wp *                        &
    290                               ( sums_us2_ws_l(:,i) + sums_vs2_ws_l(:,i) +      &
    291                                 sums_ws2_ws_l(:,i) )    ! e*
     274             sums_l(:,13,i) = sums_wsus_ws_l(:,i) * momentumflux_output_conversion ! w*u*
     275             sums_l(:,15,i) = sums_wsvs_ws_l(:,i) * momentumflux_output_conversion ! w*v*
     276             sums_l(:,30,i) = sums_us2_ws_l(:,i)                                   ! u*2
     277             sums_l(:,31,i) = sums_vs2_ws_l(:,i)                                   ! v*2
     278             sums_l(:,32,i) = sums_ws2_ws_l(:,i)                                   ! w*2
     279             sums_l(:,34,i) = sums_l(:,34,i) + 0.5_wp *                                            &
     280                              ( sums_us2_ws_l(:,i) + sums_vs2_ws_l(:,i) + sums_ws2_ws_l(:,i) )  ! e*
    292281          ENDDO
    293282
     
    297286
    298287          DO  i = 0, threads_per_task-1
    299              sums_l(:,17,i)                        = sums_wspts_ws_l(:,i)      &
    300                                            * heatflux_output_conversion  ! w*pt*
    301              IF ( ocean_mode     ) sums_l(:,66,i)  = sums_wssas_ws_l(:,i) ! w*sa*
    302              IF ( humidity       ) sums_l(:,49,i)  = sums_wsqs_ws_l(:,i)       &
    303                                            * waterflux_output_conversion  ! w*q*
    304              IF ( passive_scalar ) sums_l(:,114,i) = sums_wsss_ws_l(:,i)  ! w*s*
     288             sums_l(:,17,i)                        = sums_wspts_ws_l(:,i)                          &
     289                                                     * heatflux_output_conversion   ! w*pt*
     290             IF ( ocean_mode     ) sums_l(:,66,i)  = sums_wssas_ws_l(:,i)           ! w*sa*
     291             IF ( humidity       ) sums_l(:,49,i)  = sums_wsqs_ws_l(:,i)                           &
     292                                                     * waterflux_output_conversion  ! w*q*
     293             IF ( passive_scalar ) sums_l(:,114,i) = sums_wsss_ws_l(:,i)            ! w*s*
    305294          ENDDO
    306295
     
    308297!
    309298!--    Horizontally averaged profiles of horizontal velocities and temperature.
    310 !--    They must have been computed before, because they are already required
    311 !--    for other horizontal averages.
     299!--    They must have been computed before, because they are already required for other horizontal
     300!--    averages.
    312301       tn = 0
    313302       !$OMP PARALLEL PRIVATE( i, j, k, tn, flag )
     
    321310                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 22 ) )
    322311                !$ACC ATOMIC
    323                 sums_l(k,1,tn)  = sums_l(k,1,tn)  + u(k,j,i)  * rmask(j,i,sr)  &
    324                                                               * flag
    325                 !$ACC ATOMIC
    326                 sums_l(k,2,tn)  = sums_l(k,2,tn)  + v(k,j,i)  * rmask(j,i,sr)  &
    327                                                               * flag
    328                 !$ACC ATOMIC
    329                 sums_l(k,4,tn)  = sums_l(k,4,tn)  + pt(k,j,i) * rmask(j,i,sr)  &
    330                                                               * flag
     312                sums_l(k,1,tn)  = sums_l(k,1,tn)  + u(k,j,i)  * rmask(j,i,sr) * flag
     313                !$ACC ATOMIC
     314                sums_l(k,2,tn)  = sums_l(k,2,tn)  + v(k,j,i)  * rmask(j,i,sr) * flag
     315                !$ACC ATOMIC
     316                sums_l(k,4,tn)  = sums_l(k,4,tn)  + pt(k,j,i) * rmask(j,i,sr) * flag
    331317             ENDDO
    332318          ENDDO
     
    341327             DO  j =  nys, nyn
    342328                DO  k = nzb, nzt+1
    343                    sums_l(k,23,tn)  = sums_l(k,23,tn) + sa(k,j,i)              &
    344                                     * rmask(j,i,sr)                            &
    345                                     * MERGE( 1.0_wp, 0.0_wp,                   &
    346                                              BTEST( wall_flags_total_0(k,j,i), 22 ) )
     329                   sums_l(k,23,tn)  = sums_l(k,23,tn) + sa(k,j,i)                                  &
     330                                      * rmask(j,i,sr)                                              &
     331                                      * MERGE( 1.0_wp, 0.0_wp,                                     &
     332                                               BTEST( wall_flags_total_0(k,j,i), 22 ) )
    347333                ENDDO
    348334             ENDDO
     
    351337
    352338!
    353 !--    Horizontally averaged profiles of virtual potential temperature,
    354 !--    total water content, water vapor mixing ratio and liquid water potential
    355 !--    temperature
     339!--    Horizontally averaged profiles of virtual potential temperature, total water content, water
     340!--    vapor mixing ratio and liquid water potential temperature
    356341       IF ( humidity )  THEN
    357342          !$OMP DO
     
    360345                DO  k = nzb, nzt+1
    361346                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 22 ) )
    362                    sums_l(k,44,tn)  = sums_l(k,44,tn) +                        &
    363                                       vpt(k,j,i) * rmask(j,i,sr) * flag
    364                    sums_l(k,41,tn)  = sums_l(k,41,tn) +                        &
    365                                       q(k,j,i) * rmask(j,i,sr)   * flag
     347                   sums_l(k,44,tn)  = sums_l(k,44,tn) + vpt(k,j,i) * rmask(j,i,sr) * flag
     348                   sums_l(k,41,tn)  = sums_l(k,41,tn) + q(k,j,i) * rmask(j,i,sr)   * flag
    366349                ENDDO
    367350             ENDDO
     
    374357                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 22 ) )
    375358                      sums_l(k,42,tn) = sums_l(k,42,tn) +                      &
    376                                       ( q(k,j,i) - ql(k,j,i) ) * rmask(j,i,sr) &
    377                                                                * flag
    378                       sums_l(k,43,tn) = sums_l(k,43,tn) + (                    &
    379                                       pt(k,j,i) + lv_d_cp * d_exner(k) * ql(k,j,i) &
    380                                                           ) * rmask(j,i,sr)    &
    381                                                             * flag
     359                                        ( q(k,j,i) - ql(k,j,i) ) * rmask(j,i,sr) * flag
     360                      sums_l(k,43,tn) = sums_l(k,43,tn) + (                                        &
     361                                           pt(k,j,i) + lv_d_cp * d_exner(k) * ql(k,j,i)            &
     362                                                          ) * rmask(j,i,sr) * flag
    382363                   ENDDO
    383364                ENDDO
     
    393374             DO  j =  nys, nyn
    394375                DO  k = nzb, nzt+1
    395                    sums_l(k,115,tn)  = sums_l(k,115,tn) + s(k,j,i)             &
    396                                     * rmask(j,i,sr)                            &
    397                                     * MERGE( 1.0_wp, 0.0_wp,                   &
    398                                              BTEST( wall_flags_total_0(k,j,i), 22 ) )
     376                   sums_l(k,115,tn)  = sums_l(k,115,tn) + s(k,j,i)                                 &
     377                                       * rmask(j,i,sr)                                             &
     378                                       * MERGE( 1.0_wp, 0.0_wp,                                    &
     379                                                BTEST( wall_flags_total_0(k,j,i), 22 ) )
    399380                ENDDO
    400381             ENDDO
     
    430411!--    Compute total sum from local sums
    431412       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    432        CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, MPI_REAL, &
    433                            MPI_SUM, comm2d, ierr )
     413       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), nzt+2-nzb, MPI_REAL, MPI_SUM, comm2d,    &
     414                           ierr )
    434415       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    435        CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, MPI_REAL, &
    436                            MPI_SUM, comm2d, ierr )
     416       CALL MPI_ALLREDUCE( sums_l(nzb,2,0), sums(nzb,2), nzt+2-nzb, MPI_REAL, MPI_SUM, comm2d,    &
     417                           ierr )
    437418       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    438        CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, MPI_REAL, &
    439                            MPI_SUM, comm2d, ierr )
     419       CALL MPI_ALLREDUCE( sums_l(nzb,4,0), sums(nzb,4), nzt+2-nzb, MPI_REAL, MPI_SUM, comm2d,    &
     420                           ierr )
    440421       IF ( ocean_mode )  THEN
    441422          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    442           CALL MPI_ALLREDUCE( sums_l(nzb,23,0), sums(nzb,23), nzt+2-nzb,       &
    443                               MPI_REAL, MPI_SUM, comm2d, ierr )
     423          CALL MPI_ALLREDUCE( sums_l(nzb,23,0), sums(nzb,23), nzt+2-nzb, MPI_REAL, MPI_SUM, comm2d,&
     424                              ierr )
    444425       ENDIF
    445426       IF ( humidity ) THEN
    446427          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    447           CALL MPI_ALLREDUCE( sums_l(nzb,44,0), sums(nzb,44), nzt+2-nzb,       &
    448                               MPI_REAL, MPI_SUM, comm2d, ierr )
     428          CALL MPI_ALLREDUCE( sums_l(nzb,44,0), sums(nzb,44), nzt+2-nzb, MPI_REAL, MPI_SUM, comm2d,&
     429                              ierr )
    449430          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    450           CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb,       &
    451                               MPI_REAL, MPI_SUM, comm2d, ierr )
     431          CALL MPI_ALLREDUCE( sums_l(nzb,41,0), sums(nzb,41), nzt+2-nzb, MPI_REAL, MPI_SUM, comm2d,&
     432                              ierr )
    452433          IF ( bulk_cloud_model ) THEN
    453434             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    454              CALL MPI_ALLREDUCE( sums_l(nzb,42,0), sums(nzb,42), nzt+2-nzb,    &
    455                                  MPI_REAL, MPI_SUM, comm2d, ierr )
     435             CALL MPI_ALLREDUCE( sums_l(nzb,42,0), sums(nzb,42), nzt+2-nzb, MPI_REAL, MPI_SUM,     &
     436                                 comm2d, ierr )
    456437             IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    457              CALL MPI_ALLREDUCE( sums_l(nzb,43,0), sums(nzb,43), nzt+2-nzb,    &
    458                                  MPI_REAL, MPI_SUM, comm2d, ierr )
     438             CALL MPI_ALLREDUCE( sums_l(nzb,43,0), sums(nzb,43), nzt+2-nzb, MPI_REAL, MPI_SUM,     &
     439                                 comm2d, ierr )
    459440          ENDIF
    460441       ENDIF
     
    462443       IF ( passive_scalar )  THEN
    463444          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    464           CALL MPI_ALLREDUCE( sums_l(nzb,115,0), sums(nzb,115), nzt+2-nzb,       &
    465                               MPI_REAL, MPI_SUM, comm2d, ierr )
     445          CALL MPI_ALLREDUCE( sums_l(nzb,115,0), sums(nzb,115), nzt+2-nzb, MPI_REAL, MPI_SUM,      &
     446                              comm2d, ierr )
    466447       ENDIF
    467448#else
     
    482463
    483464!
    484 !--    Final values are obtained by division by the total number of grid points
    485 !--    used for summation. After that store profiles.
     465!--    Final values are obtained by division by the total number of grid points used for summation.
     466!--    After that store profiles.
    486467       sums(:,1) = sums(:,1) / ngp_2dh(sr)
    487468       sums(:,2) = sums(:,2) / ngp_2dh(sr)
     
    517498!
    518499!--    Passive scalar
    519        IF ( passive_scalar )  hom(:,1,115,sr) = sums(:,115) /                  &
    520             ngp_2dh_s_inner(:,sr)                    ! s
    521 
    522 !
    523 !--    Horizontally averaged profiles of the remaining prognostic variables,
    524 !--    variances, the total and the perturbation energy (single values in last
    525 !--    column of sums_l) and some diagnostic quantities.
    526 !--    NOTE: for simplicity, nzb_s_inner is used below, although strictly
    527 !--    ----  speaking the following k-loop would have to be split up and
    528 !--          rearranged according to the staggered grid.
    529 !--          However, this implies no error since staggered velocity components
    530 !--          are zero at the walls and inside buildings.
     500       IF ( passive_scalar )  hom(:,1,115,sr) = sums(:,115) / ngp_2dh_s_inner(:,sr)  ! s
     501
     502!
     503!--    Horizontally averaged profiles of the remaining prognostic variables, variances, the total
     504!--    and the perturbation energy (single values in last column of sums_l) and some diagnostic
     505!--    quantities.
     506!--    NOTE: for simplicity, nzb_s_inner is used below, although strictly speaking the following
     507!--    ----  k-loop would have to be split up and rearranged according to the staggered grid.
     508!--          However, this implies no error since staggered velocity components are zero at the
     509!--          walls and inside buildings.
    531510       tn = 0
    532511       !$OMP PARALLEL PRIVATE( i, j, k, pts, sums_ll,                          &
     
    554533!--             Prognostic and diagnostic variables
    555534                !$ACC ATOMIC
    556                 sums_l(k,3,tn)  = sums_l(k,3,tn)  + w(k,j,i)  * rmask(j,i,sr)  &
    557                                                               * flag
    558                 !$ACC ATOMIC
    559                 sums_l(k,8,tn)  = sums_l(k,8,tn)  + e(k,j,i)  * rmask(j,i,sr)  &
    560                                                               * flag
    561                 !$ACC ATOMIC
    562                 sums_l(k,9,tn)  = sums_l(k,9,tn)  + km(k,j,i) * rmask(j,i,sr)  &
    563                                                               * flag
    564                 !$ACC ATOMIC
    565                 sums_l(k,10,tn) = sums_l(k,10,tn) + kh(k,j,i) * rmask(j,i,sr)  &
    566                                                               * flag
    567                 !$ACC ATOMIC
    568                 sums_l(k,40,tn) = sums_l(k,40,tn) + ( p(k,j,i)                 &
    569                                          / momentumflux_output_conversion(k) ) &
    570                                                               * flag
     535                sums_l(k,3,tn)  = sums_l(k,3,tn)  + w(k,j,i)  * rmask(j,i,sr) * flag
     536                !$ACC ATOMIC
     537                sums_l(k,8,tn)  = sums_l(k,8,tn)  + e(k,j,i)  * rmask(j,i,sr) * flag
     538                !$ACC ATOMIC
     539                sums_l(k,9,tn)  = sums_l(k,9,tn)  + km(k,j,i) * rmask(j,i,sr) * flag
     540                !$ACC ATOMIC
     541                sums_l(k,10,tn) = sums_l(k,10,tn) + kh(k,j,i) * rmask(j,i,sr) * flag
     542                !$ACC ATOMIC
     543                sums_l(k,40,tn) = sums_l(k,40,tn) + ( p(k,j,i)                                     &
     544                                         / momentumflux_output_conversion(k) ) * flag
    571545
    572546                !$ACC ATOMIC
    573547                sums_l(k,33,tn) = sums_l(k,33,tn) + &
    574                                   ( pt(k,j,i)-hom(k,1,4,sr) )**2 * rmask(j,i,sr)&
    575                                                                  * flag
     548                                  ( pt(k,j,i)-hom(k,1,4,sr) )**2 * rmask(j,i,sr) * flag
    576549#ifndef _OPENACC
    577550                IF ( humidity )  THEN
    578                    sums_l(k,70,tn) = sums_l(k,70,tn) + &
    579                                   ( q(k,j,i)-hom(k,1,41,sr) )**2 * rmask(j,i,sr)&
    580                                                                  * flag
     551                   sums_l(k,70,tn) = sums_l(k,70,tn) +                                             &
     552                                     ( q(k,j,i)-hom(k,1,41,sr) )**2 * rmask(j,i,sr) * flag
    581553                ENDIF
    582554                IF ( passive_scalar )  THEN
    583                    sums_l(k,116,tn) = sums_l(k,116,tn) + &
    584                                   ( s(k,j,i)-hom(k,1,115,sr) )**2 * rmask(j,i,sr)&
    585                                                                   * flag
     555                   sums_l(k,116,tn) = sums_l(k,116,tn) +                                           &
     556                                      ( s(k,j,i)-hom(k,1,115,sr) )**2 * rmask(j,i,sr) * flag
    586557                ENDIF
    587558#endif
     
    590561!--             (Computation of the skewness of w further below)
    591562                !$ACC ATOMIC
    592                 sums_l(k,38,tn) = sums_l(k,38,tn) + w(k,j,i)**3 * rmask(j,i,sr) &
    593                                                                 * flag
    594 
    595                 sums_l_etot  = sums_l_etot + &
    596                                         0.5_wp * ( u(k,j,i)**2 + v(k,j,i)**2 +  &
    597                                         w(k,j,i)**2 )            * rmask(j,i,sr)&
    598                                                                  * flag
    599 
    600 !
    601 !--             Computation of the Kolmogorov length scale. Calculation is based
    602 !--             on gradients of the deviations from the horizontal mean.
     563                sums_l(k,38,tn) = sums_l(k,38,tn) + w(k,j,i)**3 * rmask(j,i,sr) * flag
     564
     565                sums_l_etot  = sums_l_etot + 0.5_wp * ( u(k,j,i)**2 + v(k,j,i)**2 + w(k,j,i)**2 )  &
     566                                             * rmask(j,i,sr) * flag
     567
     568!
     569!--             Computation of the Kolmogorov length scale. Calculation is based on gradients of the
     570!--             deviations from the horizontal mean.
    603571!--             Kolmogorov scale at the boundaries (k=0/z=0m and k=nzt+1) is set to zero.
    604572                IF ( kolmogorov_length_scale .AND. k /= nzb .AND. k /= nzt+1) THEN
     
    607575!
    608576!--                Calculate components of the fluctuating rate-of-strain tensor
    609 !--                (0.5*(del u'_i/del x_j + del u'_j/del x_i)) defined in the
    610 !--                center of each grid box.
    611                    du_dx = ( ( u(k,j,i+1) - hom(k,1,1,sr) ) -                  &
     577!--                (0.5*(del u'_i/del x_j + del u'_j/del x_i)) defined in the center of each grid
     578!--                box.
     579                   du_dx = ( ( u(k,j,i+1) - hom(k,1,1,sr) ) -                                      &
    612580                             ( u(k,j,i) - hom(k,1,1,sr) ) ) * ddx
    613                    du_dy = 0.25_wp * ddy *                                     &
    614                            ( ( u(k,j+1,i) - hom(k,1,1,sr) ) -                  &
    615                              ( u(k,j-1,i) - hom(k,1,1,sr) ) +                  &
    616                              ( u(k,j+1,i+1) - hom(k,1,1,sr) ) -                &
     581                   du_dy = 0.25_wp * ddy *                                                         &
     582                           ( ( u(k,j+1,i) - hom(k,1,1,sr) ) -                                      &
     583                             ( u(k,j-1,i) - hom(k,1,1,sr) ) +                                      &
     584                             ( u(k,j+1,i+1) - hom(k,1,1,sr) ) -                                    &
    617585                             ( u(k,j-1,i+1) - hom(k,1,1,sr) ) )
    618                    du_dz = 0.25_wp * ( ( ( u(k+1,j,i) - hom(k+1,1,1,sr) ) -    &
    619                                          ( u(k,j,i) - hom(k,1,1,sr) ) ) *      &
    620                                         ddzu(k+1) +                            &
    621                                        ( ( u(k,j,i) - hom(k,1,1,sr) ) -        &
    622                                          ( u(k-1,j,i) - hom(k-1,1,1,sr) ) )*   &
    623                                         ddzu(k) +                              &
    624                                        ( ( u(k+1,j,i+1) - hom(k+1,1,1,sr) )-   &
    625                                          ( u(k,j,i+1) - hom(k,1,1,sr) ) ) *    &
    626                                         ddzu(k+1) +                            &
    627                                        ( ( u(k,j,i+1) - hom(k,1,1,sr) ) -      &
    628                                          ( u(k-1,j,i+1) - hom(k-1,1,1,sr) ) ) *&
    629                                         ddzu(k) )
    630 
    631                    dv_dx = 0.25_wp * ddx *                                     &
    632                            ( ( v(k,j,i+1) - hom(k,1,2,sr) ) -                  &
    633                              ( v(k,j,i-1) - hom(k,1,2,sr) ) +                  &
    634                              ( v(k,j+1,i+1) - hom(k,1,2,sr) ) -                &
     586                   du_dz = 0.25_wp * ( ( ( u(k+1,j,i) - hom(k+1,1,1,sr) ) -                        &
     587                                         ( u(k,j,i) - hom(k,1,1,sr) ) ) *                          &
     588                                       ddzu(k+1) +                                                 &
     589                                       ( ( u(k,j,i) - hom(k,1,1,sr) ) -                            &
     590                                         ( u(k-1,j,i) - hom(k-1,1,1,sr) ) ) *                      &
     591                                       ddzu(k) +                                                   &
     592                                       ( ( u(k+1,j,i+1) - hom(k+1,1,1,sr) ) -                      &
     593                                         ( u(k,j,i+1) - hom(k,1,1,sr) ) ) *                        &
     594                                       ddzu(k+1) +                                                 &
     595                                       ( ( u(k,j,i+1) - hom(k,1,1,sr) ) -                          &
     596                                         ( u(k-1,j,i+1) - hom(k-1,1,1,sr) ) ) *                    &
     597                                       ddzu(k) )
     598
     599                   dv_dx = 0.25_wp * ddx *                                                         &
     600                           ( ( v(k,j,i+1) - hom(k,1,2,sr) ) -                                      &
     601                             ( v(k,j,i-1) - hom(k,1,2,sr) ) +                                      &
     602                             ( v(k,j+1,i+1) - hom(k,1,2,sr) ) -                                    &
    635603                             ( v(k,j+1,i-1) - hom(k,1,2,sr) ) )
    636                    dv_dy = ( ( v(k,j+1,i) - hom(k,1,2,sr) ) -                  &
    637                              ( v(k,j,i) - hom(k,1,2,sr) ) ) * ddy
    638                    dv_dz = 0.25_wp * ( ( ( v(k+1,j,i) - hom(k+1,1,2,sr) ) -    &
    639                                          ( v(k,j,i) - hom(k,1,2,sr) ) ) *      &
    640                                         ddzu(k+1) +                            &
    641                                        ( ( v(k,j,i) - hom(k,1,2,sr) ) -        &
    642                                          ( v(k-1,j,i) - hom(k-1,1,2,sr) ) ) *  &
    643                                         ddzu(k) +                              &
    644                                        ( ( v(k+1,j+1,i) - hom(k+1,1,2,sr) ) -  &
    645                                          ( v(k,j+1,i) - hom(k,1,2,sr) ) ) *    &
    646                                         ddzu(k+1) +                            &
    647                                        ( ( v(k,j+1,i) - hom(k,1,2,sr) ) -      &
    648                                          ( v(k-1,j+1,i) - hom(k-1,1,2,sr) ) ) *&
    649                                         ddzu(k) )
    650 
    651                    dw_dx = 0.25_wp * ddx * ( w(k,j,i+1) - w(k,j,i-1) +         &
    652                                              w(k-1,j,i+1) - w(k-1,j,i-1) )
    653                    dw_dy = 0.25_wp * ddy * ( w(k,j+1,i) - w(k,j-1,i) +         &
    654                                              w(k-1,j+1,i) - w(k-1,j-1,i) )
     604                   dv_dy = ( ( v(k,j+1,i) - hom(k,1,2,sr) ) - ( v(k,j,i) - hom(k,1,2,sr) ) ) * ddy
     605                   dv_dz = 0.25_wp * ( ( ( v(k+1,j,i) - hom(k+1,1,2,sr) ) -                        &
     606                                         ( v(k,j,i) - hom(k,1,2,sr) ) ) *                          &
     607                                       ddzu(k+1) +                                                 &
     608                                       ( ( v(k,j,i) - hom(k,1,2,sr) ) -                            &
     609                                         ( v(k-1,j,i) - hom(k-1,1,2,sr) ) ) *                      &
     610                                       ddzu(k) +                                                   &
     611                                       ( ( v(k+1,j+1,i) - hom(k+1,1,2,sr) ) -                      &
     612                                         ( v(k,j+1,i) - hom(k,1,2,sr) ) ) *                        &
     613                                       ddzu(k+1) +                                                 &
     614                                       ( ( v(k,j+1,i) - hom(k,1,2,sr) ) -                          &
     615                                         ( v(k-1,j+1,i) - hom(k-1,1,2,sr) ) ) *                    &
     616                                       ddzu(k) )
     617
     618                   dw_dx = 0.25_wp * ddx * ( w(k,j,i+1) - w(k,j,i-1) + w(k-1,j,i+1) - w(k-1,j,i-1) )
     619                   dw_dy = 0.25_wp * ddy * ( w(k,j+1,i) - w(k,j-1,i) + w(k-1,j+1,i) - w(k-1,j-1,i) )
    655620                   dw_dz = ( w(k,j,i) - w(k-1,j,i) ) * ddzw(k)
    656621
     
    667632                   s33 = 0.5_wp * ( dw_dz + dw_dz )
    668633
    669 !--                Calculate 3D instantaneous energy dissipation rate after
    670 !--                Pope (2000): Turbulent flows, p.259. It is defined in the center
    671 !--                of each grid volume.
    672                    dissipation = 2.0_wp * km(k,j,i) *                          &
    673                                 ( s11*s11 + s21*s21 + s31*s31 +                &
    674                                  s12*s12 + s22*s22 + s32*s32 +                 &
    675                                 s13*s13 + s23*s23 + s33*s33 )
     634!--                Calculate 3D instantaneous energy dissipation rate following Pope (2000):
     635!--                Turbulent flows, p.259. It is defined in the center of each grid volume.
     636                   dissipation = 2.0_wp * km(k,j,i) *                                              &
     637                                ( s11*s11 + s21*s21 + s31*s31 +                                    &
     638                                  s12*s12 + s22*s22 + s32*s32 +                                    &
     639                                  s13*s13 + s23*s23 + s33*s33 )
    676640                   eta         = ( km(k,j,i)**3.0_wp / ( dissipation+1.0E-12 ) )**(1.0_wp/4.0_wp)
    677641
    678642                   !$ACC ATOMIC
    679                    sums_l(k,121,tn) = sums_l(k,121,tn) + eta * rmask(j,i,sr)   &
    680                                                                 * flag
     643                   sums_l(k,121,tn) = sums_l(k,121,tn) + eta * rmask(j,i,sr) * flag
    681644
    682645
     
    685648             ENDDO !k-loop
    686649!
    687 !--          Total and perturbation energy for the total domain (being
    688 !--          collected in the last column of sums_l). Summation of these
    689 !--          quantities is seperated from the previous loop in order to
    690 !--          allow vectorization of that loop.
     650!--          Total and perturbation energy for the total domain (being collected in the last column
     651!--          of sums_l). Summation of these quantities is seperated from the previous loop in order
     652!--          to allow vectorization of that loop.
    691653             !$ACC ATOMIC
    692654             sums_l(nzb+4,pr_palm,tn) = sums_l(nzb+4,pr_palm,tn) + sums_l_etot
    693655!
    694656!--          2D-arrays (being collected in the last column of sums_l)
    695              IF ( surf_def_h(0)%end_index(j,i) >=                              &
    696                   surf_def_h(0)%start_index(j,i) )  THEN
     657             IF ( surf_def_h(0)%end_index(j,i) >= surf_def_h(0)%start_index(j,i) )  THEN
    697658                m = surf_def_h(0)%start_index(j,i)
    698659                !$ACC ATOMIC
    699                 sums_l(nzb,pr_palm,tn)   = sums_l(nzb,pr_palm,tn) +            &
    700                                         surf_def_h(0)%us(m)   * rmask(j,i,sr)
    701                 !$ACC ATOMIC
    702                 sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) +          &
    703                                         surf_def_h(0)%usws(m) * rmask(j,i,sr)
    704                 !$ACC ATOMIC
    705                 sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) +          &
    706                                         surf_def_h(0)%vsws(m) * rmask(j,i,sr)
    707                 !$ACC ATOMIC
    708                 sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) +          &
    709                                         surf_def_h(0)%ts(m)   * rmask(j,i,sr)
     660                sums_l(nzb,pr_palm,tn)   = sums_l(nzb,pr_palm,tn) +                                &
     661                                           surf_def_h(0)%us(m)   * rmask(j,i,sr)
     662                !$ACC ATOMIC
     663                sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) +                              &
     664                                           surf_def_h(0)%usws(m) * rmask(j,i,sr)
     665                !$ACC ATOMIC
     666                sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) +                              &
     667                                           surf_def_h(0)%vsws(m) * rmask(j,i,sr)
     668                !$ACC ATOMIC
     669                sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) +                              &
     670                                           surf_def_h(0)%ts(m)   * rmask(j,i,sr)
    710671#ifndef _OPENACC
    711672                IF ( humidity )  THEN
    712                    sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) +     &
    713                                             surf_def_h(0)%qs(m)   * rmask(j,i,sr)
     673                   sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) +                         &
     674                                               surf_def_h(0)%qs(m)   * rmask(j,i,sr)
    714675                ENDIF
    715676                IF ( passive_scalar )  THEN
    716                    sums_l(nzb+13,pr_palm,tn) = sums_l(nzb+13,pr_palm,tn) +     &
    717                                             surf_def_h(0)%ss(m)   * rmask(j,i,sr)
     677                   sums_l(nzb+13,pr_palm,tn) = sums_l(nzb+13,pr_palm,tn) +                         &
     678                                               surf_def_h(0)%ss(m)   * rmask(j,i,sr)
    718679                ENDIF
    719680#endif
     
    721682!--             Summation of surface temperature.
    722683                !$ACC ATOMIC
    723                 sums_l(nzb+14,pr_palm,tn) = sums_l(nzb+14,pr_palm,tn)   +      &
    724                                             surf_def_h(0)%pt_surface(m) *      &
    725                                             rmask(j,i,sr)
     684                sums_l(nzb+14,pr_palm,tn) = sums_l(nzb+14,pr_palm,tn)   +                          &
     685                                            surf_def_h(0)%pt_surface(m) * rmask(j,i,sr)
    726686             ENDIF
    727687             IF ( surf_lsm_h%end_index(j,i) >= surf_lsm_h%start_index(j,i) )  THEN
    728688                m = surf_lsm_h%start_index(j,i)
    729689                !$ACC ATOMIC
    730                 sums_l(nzb,pr_palm,tn)   = sums_l(nzb,pr_palm,tn) +            &
    731                                         surf_lsm_h%us(m)   * rmask(j,i,sr)
    732                 !$ACC ATOMIC
    733                 sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) +          &
    734                                         surf_lsm_h%usws(m) * rmask(j,i,sr)
    735                 !$ACC ATOMIC
    736                 sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) +          &
    737                                         surf_lsm_h%vsws(m) * rmask(j,i,sr)
    738                 !$ACC ATOMIC
    739                 sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) +          &
    740                                         surf_lsm_h%ts(m)   * rmask(j,i,sr)
     690                sums_l(nzb,pr_palm,tn)   = sums_l(nzb,pr_palm,tn) +                                &
     691                                           surf_lsm_h%us(m)   * rmask(j,i,sr)
     692                !$ACC ATOMIC
     693                sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) +                              &
     694                                           surf_lsm_h%usws(m) * rmask(j,i,sr)
     695                !$ACC ATOMIC
     696                sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) +                              &
     697                                           surf_lsm_h%vsws(m) * rmask(j,i,sr)
     698                !$ACC ATOMIC
     699                sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) +                              &
     700                                           surf_lsm_h%ts(m)   * rmask(j,i,sr)
    741701#ifndef _OPENACC
    742702                IF ( humidity )  THEN
    743                    sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) +     &
    744                                             surf_lsm_h%qs(m)   * rmask(j,i,sr)
     703                   sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) +                         &
     704                                               surf_lsm_h%qs(m)   * rmask(j,i,sr)
    745705                ENDIF
    746706                IF ( passive_scalar )  THEN
    747                    sums_l(nzb+13,pr_palm,tn) = sums_l(nzb+13,pr_palm,tn) +     &
    748                                             surf_lsm_h%ss(m)   * rmask(j,i,sr)
     707                   sums_l(nzb+13,pr_palm,tn) = sums_l(nzb+13,pr_palm,tn) +                         &
     708                                               surf_lsm_h%ss(m)   * rmask(j,i,sr)
    749709                ENDIF
    750710#endif
     
    752712!--             Summation of surface temperature.
    753713                !$ACC ATOMIC
    754                 sums_l(nzb+14,pr_palm,tn) = sums_l(nzb+14,pr_palm,tn)   +      &
    755                                             surf_lsm_h%pt_surface(m)    *      &
    756                                             rmask(j,i,sr)
     714                sums_l(nzb+14,pr_palm,tn) = sums_l(nzb+14,pr_palm,tn) +                            &
     715                                            surf_lsm_h%pt_surface(m) * rmask(j,i,sr)
    757716             ENDIF
    758717             IF ( surf_usm_h%end_index(j,i) >= surf_usm_h%start_index(j,i) )  THEN
    759718                m = surf_usm_h%start_index(j,i)
    760719                !$ACC ATOMIC
    761                 sums_l(nzb,pr_palm,tn)   = sums_l(nzb,pr_palm,tn) +            &
    762                                         surf_usm_h%us(m)   * rmask(j,i,sr)
    763                 !$ACC ATOMIC
    764                 sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) +          &
    765                                         surf_usm_h%usws(m) * rmask(j,i,sr)
    766                 !$ACC ATOMIC
    767                 sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) +          &
    768                                         surf_usm_h%vsws(m) * rmask(j,i,sr)
    769                 !$ACC ATOMIC
    770                 sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) +          &
    771                                         surf_usm_h%ts(m)   * rmask(j,i,sr)
     720                sums_l(nzb,pr_palm,tn)   = sums_l(nzb,pr_palm,tn) +                                &
     721                                           surf_usm_h%us(m)   * rmask(j,i,sr)
     722                !$ACC ATOMIC
     723                sums_l(nzb+1,pr_palm,tn) = sums_l(nzb+1,pr_palm,tn) +                              &
     724                                           surf_usm_h%usws(m) * rmask(j,i,sr)
     725                !$ACC ATOMIC
     726                sums_l(nzb+2,pr_palm,tn) = sums_l(nzb+2,pr_palm,tn) +                              &
     727                                           surf_usm_h%vsws(m) * rmask(j,i,sr)
     728                !$ACC ATOMIC
     729                sums_l(nzb+3,pr_palm,tn) = sums_l(nzb+3,pr_palm,tn) +                              &
     730                                           surf_usm_h%ts(m)   * rmask(j,i,sr)
    772731#ifndef _OPENACC
    773732                IF ( humidity )  THEN
    774                    sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) +     &
    775                                             surf_usm_h%qs(m)   * rmask(j,i,sr)
     733                   sums_l(nzb+12,pr_palm,tn) = sums_l(nzb+12,pr_palm,tn) +                         &
     734                                               surf_usm_h%qs(m)   * rmask(j,i,sr)
    776735                ENDIF
    777736                IF ( passive_scalar )  THEN
    778                    sums_l(nzb+13,pr_palm,tn) = sums_l(nzb+13,pr_palm,tn) +     &
    779                                             surf_usm_h%ss(m)  * rmask(j,i,sr)
     737                   sums_l(nzb+13,pr_palm,tn) = sums_l(nzb+13,pr_palm,tn) +                         &
     738                                               surf_usm_h%ss(m) * rmask(j,i,sr)
    780739                ENDIF
    781740#endif
     
    783742!--             Summation of surface temperature.
    784743                !$ACC ATOMIC
    785                 sums_l(nzb+14,pr_palm,tn) = sums_l(nzb+14,pr_palm,tn)   +      &
    786                                             surf_usm_h%pt_surface(m)    *      &
    787                                             rmask(j,i,sr)
     744                sums_l(nzb+14,pr_palm,tn) = sums_l(nzb+14,pr_palm,tn) +                            &
     745                                            surf_usm_h%pt_surface(m)  * rmask(j,i,sr)
    788746             ENDIF
    789747          ENDDO !j-loop
     
    798756!--    Computation of statistics when ws-scheme is not used. Else these
    799757!--    quantities are evaluated in the advection routines.
    800        IF ( .NOT. ws_scheme_mom .OR. sr /= 0 .OR. simulated_time == 0.0_wp )   &
    801        THEN
     758       IF ( .NOT. ws_scheme_mom .OR. sr /= 0 .OR. simulated_time == 0.0_wp )  THEN
    802759          !$OMP DO
    803760          DO  i = nxl, nxr
     
    812769                   vst2 = ( v(k,j,i) - hom(k,1,2,sr) )**2
    813770
    814                    sums_l(k,30,tn) = sums_l(k,30,tn) + ust2 * rmask(j,i,sr)    &
    815                                                             * flag
    816                    sums_l(k,31,tn) = sums_l(k,31,tn) + vst2 * rmask(j,i,sr)    &
    817                                                             * flag
    818                    sums_l(k,32,tn) = sums_l(k,32,tn) + w2   * rmask(j,i,sr)    &
    819                                                             * flag
     771                   sums_l(k,30,tn) = sums_l(k,30,tn) + ust2 * rmask(j,i,sr) * flag
     772                   sums_l(k,31,tn) = sums_l(k,31,tn) + vst2 * rmask(j,i,sr) * flag
     773                   sums_l(k,32,tn) = sums_l(k,32,tn) + w2   * rmask(j,i,sr) * flag
    820774!
    821775!--                Perturbation energy
    822 
    823                    sums_l(k,34,tn) = sums_l(k,34,tn) + 0.5_wp *                &
    824                                   ( ust2 + vst2 + w2 )      * rmask(j,i,sr)    &
    825                                                             * flag
     776                   sums_l(k,34,tn) = sums_l(k,34,tn) +                                             &
     777                                     0.5_wp * ( ust2 + vst2 + w2 ) * rmask(j,i,sr) * flag
    826778                ENDDO
    827779             ENDDO
     
    829781       ENDIF
    830782!
    831 !--    Computaion of domain-averaged perturbation energy. Please note,
    832 !--    to prevent that perturbation energy is larger (even if only slightly)
    833 !--    than the total kinetic energy, calculation is based on deviations from
    834 !--    the horizontal mean, instead of spatial descretization of the advection
     783!--    Computaion of domain-averaged perturbation energy. Please note, to prevent that perturbation
     784!--    energy is larger (even if only slightly) than the total kinetic energy, calculation is based
     785!--    on deviations from the horizontal mean, instead of spatial descretization of the advection
    835786!--    term.
    836787       !$OMP DO
     
    849800
    850801                !$ACC ATOMIC
    851                 sums_l(nzb+5,pr_palm,tn) = sums_l(nzb+5,pr_palm,tn)            &
    852                                  + 0.5_wp * ( ust2 + vst2 + w2 )               &
    853                                  * rmask(j,i,sr)                               &
    854                                  * flag
     802                sums_l(nzb+5,pr_palm,tn) = sums_l(nzb+5,pr_palm,tn)                                &
     803                                           + 0.5_wp * ( ust2 + vst2 + w2 ) * rmask(j,i,sr) * flag
    855804
    856805             ENDDO
     
    861810!
    862811!--    Horizontally averaged profiles of the vertical fluxes
    863 
    864812       !$OMP DO
    865813       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, k, l, m) &
     
    873821          DO  j = nys, nyn
    874822!
    875 !--          Subgridscale fluxes (without Prandtl layer from k=nzb,
    876 !--          oterwise from k=nzb+1)
    877 !--          NOTE: for simplicity, nzb_diff_s_inner is used below, although
    878 !--          ----  strictly speaking the following k-loop would have to be
    879 !--                split up according to the staggered grid.
    880 !--                However, this implies no error since staggered velocity
    881 !--                components are zero at the walls and inside buildings.
    882 !--          Flag 23 is used to mask surface fluxes as well as model-top fluxes,
    883 !--          which are added further below.
     823!--          Subgridscale fluxes (without Prandtl layer from k=nzb, oterwise from k=nzb+1)
     824!--          NOTE: for simplicity, nzb_diff_s_inner is used below, although strictly speaking the
     825!--          ----  following k-loop would have to be split up according to the staggered grid.
     826!--                However, this implies no error since staggered velocity components are zero at
     827!--                the walls and inside buildings.
     828!--          Flag 23 is used to mask surface fluxes as well as model-top fluxes, which are added
     829!--          further below.
    884830             DO  k = nzb, nzt
    885                 flag = MERGE( 1.0_wp, 0.0_wp,                                  &
    886                               BTEST( wall_flags_total_0(k,j,i), 23 ) ) *       &
    887                        MERGE( 1.0_wp, 0.0_wp,                                  &
    888                               BTEST( wall_flags_total_0(k,j,i), 9  ) )
     831                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 23 ) ) *           &
     832                       MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 9  ) )
    889833!
    890834!--             Momentum flux w"u"
    891835                !$ACC ATOMIC
    892                 sums_l(k,12,tn) = sums_l(k,12,tn) - 0.25_wp * (                &
    893                                km(k,j,i)+km(k+1,j,i)+km(k,j,i-1)+km(k+1,j,i-1) &
    894                                                            ) * (               &
    895                                    ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)     &
    896                                  + ( w(k,j,i)   - w(k,j,i-1) ) * ddx           &
    897                                                            ) * rmask(j,i,sr)   &
    898                                          * rho_air_zw(k)                       &
    899                                          * momentumflux_output_conversion(k)   &
    900                                          * flag
     836                sums_l(k,12,tn) = sums_l(k,12,tn) - 0.25_wp * (                                    &
     837                                  km(k,j,i) + km(k+1,j,i) + km(k,j,i-1) + km(k+1,j,i-1)            &
     838                                                              ) * (                                &
     839                                      ( u(k+1,j,i) - u(k,j,i)   ) * ddzu(k+1)                      &
     840                                    + ( w(k,j,i)   - w(k,j,i-1) ) * ddx                            &
     841                                                                  ) * rmask(j,i,sr)                &
     842                                            * rho_air_zw(k)                                        &
     843                                            * momentumflux_output_conversion(k)                    &
     844                                            * flag
    901845!
    902846!--             Momentum flux w"v"
    903847                !$ACC ATOMIC
    904                 sums_l(k,14,tn) = sums_l(k,14,tn) - 0.25_wp * (                &
    905                                km(k,j,i)+km(k+1,j,i)+km(k,j-1,i)+km(k+1,j-1,i) &
    906                                                            ) * (               &
    907                                    ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)     &
    908                                  + ( w(k,j,i)   - w(k,j-1,i) ) * ddy           &
    909                                                            ) * rmask(j,i,sr)   &
    910                                          * rho_air_zw(k)                       &
    911                                          * momentumflux_output_conversion(k)   &
    912                                          * flag
     848                sums_l(k,14,tn) = sums_l(k,14,tn) - 0.25_wp * (                                    &
     849                                  km(k,j,i) + km(k+1,j,i) + km(k,j-1,i) + km(k+1,j-1,i)            &
     850                                                              ) * (                                &
     851                                      ( v(k+1,j,i) - v(k,j,i)   ) * ddzu(k+1)                      &
     852                                    + ( w(k,j,i)   - w(k,j-1,i) ) * ddy                            &
     853                                                                  ) * rmask(j,i,sr)                &
     854                                            * rho_air_zw(k)                                        &
     855                                            * momentumflux_output_conversion(k)                    &
     856                                            * flag
    913857!
    914858!--             Heat flux w"pt"
    915859                !$ACC ATOMIC
    916                 sums_l(k,16,tn) = sums_l(k,16,tn)                              &
    917                                          - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
    918                                                * ( pt(k+1,j,i) - pt(k,j,i) )   &
    919                                                * rho_air_zw(k)                 &
    920                                                * heatflux_output_conversion(k) &
    921                                                * ddzu(k+1) * rmask(j,i,sr)     &
     860                sums_l(k,16,tn) = sums_l(k,16,tn)                                                  &
     861                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )                    &
     862                                               * ( pt(k+1,j,i) - pt(k,j,i) )                       &
     863                                               * rho_air_zw(k)                                     &
     864                                               * heatflux_output_conversion(k)                     &
     865                                               * ddzu(k+1) * rmask(j,i,sr)                         &
    922866                                               * flag
    923867
     
    926870#ifndef _OPENACC
    927871                IF ( ocean_mode )  THEN
    928                    sums_l(k,65,tn) = sums_l(k,65,tn)                           &
    929                                          - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
    930                                                * ( sa(k+1,j,i) - sa(k,j,i) )   &
    931                                                * ddzu(k+1) * rmask(j,i,sr)     &
     872                   sums_l(k,65,tn) = sums_l(k,65,tn)                                               &
     873                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )                    &
     874                                               * ( sa(k+1,j,i) - sa(k,j,i) )                       &
     875                                               * ddzu(k+1) * rmask(j,i,sr)                         &
    932876                                               * flag
    933877                ENDIF
     
    936880!--             Buoyancy flux, water flux (humidity flux) w"q"
    937881                IF ( humidity ) THEN
    938                    sums_l(k,45,tn) = sums_l(k,45,tn)                           &
    939                                          - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
    940                                                * ( vpt(k+1,j,i) - vpt(k,j,i) ) &
    941                                                * rho_air_zw(k)                 &
    942                                                * heatflux_output_conversion(k) &
     882                   sums_l(k,45,tn) = sums_l(k,45,tn)                                               &
     883                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )                    &
     884                                               * ( vpt(k+1,j,i) - vpt(k,j,i) )                     &
     885                                               * rho_air_zw(k)                                     &
     886                                               * heatflux_output_conversion(k)                     &
    943887                                               * ddzu(k+1) * rmask(j,i,sr) * flag
    944                    sums_l(k,48,tn) = sums_l(k,48,tn)                           &
    945                                          - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
    946                                                * ( q(k+1,j,i) - q(k,j,i) )     &
    947                                                * rho_air_zw(k)                 &
    948                                                * waterflux_output_conversion(k)&
     888                   sums_l(k,48,tn) = sums_l(k,48,tn)                                               &
     889                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )                    &
     890                                               * ( q(k+1,j,i) - q(k,j,i) )                         &
     891                                               * rho_air_zw(k)                                     &
     892                                               * waterflux_output_conversion(k)                    &
    949893                                               * ddzu(k+1) * rmask(j,i,sr) * flag
    950894
    951895                   IF ( bulk_cloud_model ) THEN
    952                       sums_l(k,51,tn) = sums_l(k,51,tn)                        &
    953                                          - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
    954                                                * ( ( q(k+1,j,i) - ql(k+1,j,i) )&
    955                                                 - ( q(k,j,i) - ql(k,j,i) ) )   &
    956                                                * rho_air_zw(k)                 &
    957                                                * waterflux_output_conversion(k)&
     896                      sums_l(k,51,tn) = sums_l(k,51,tn)                                            &
     897                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )                    &
     898                                               * ( ( q(k+1,j,i) - ql(k+1,j,i) )                    &
     899                                                - ( q(k,j,i) - ql(k,j,i) ) )                       &
     900                                               * rho_air_zw(k)                                     &
     901                                               * waterflux_output_conversion(k)                    &
    958902                                               * ddzu(k+1) * rmask(j,i,sr) * flag
    959903                   ENDIF
     
    963907!--             Passive scalar flux
    964908                IF ( passive_scalar )  THEN
    965                    sums_l(k,117,tn) = sums_l(k,117,tn)                         &
    966                                          - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )&
    967                                                   * ( s(k+1,j,i) - s(k,j,i) )  &
    968                                                   * ddzu(k+1) * rmask(j,i,sr)  &
     909                   sums_l(k,117,tn) = sums_l(k,117,tn)                                             &
     910                                         - 0.5_wp * ( kh(k,j,i) + kh(k+1,j,i) )                    &
     911                                                  * ( s(k+1,j,i) - s(k,j,i) )                      &
     912                                                  * ddzu(k+1) * rmask(j,i,sr)                      &
    969913                                                              * flag
    970914                ENDIF
     
    983927                   ki = -1 + l
    984928                   IF ( surf_def_h(l)%ns >= 1 )  THEN
    985                       DO  m = surf_def_h(l)%start_index(j,i),                  &
     929                      DO  m = surf_def_h(l)%start_index(j,i),                                      &
    986930                              surf_def_h(l)%end_index(j,i)
    987931                         k = surf_def_h(l)%k(m)
    988932
    989933                         !$ACC ATOMIC
    990                          sums_l(k+ki,12,tn) = sums_l(k+ki,12,tn) + &
    991                                     momentumflux_output_conversion(k+ki) * &
    992                                     surf_def_h(l)%usws(m) * rmask(j,i,sr)     ! w"u"
     934                         sums_l(k+ki,12,tn) = sums_l(k+ki,12,tn) +                                 &
     935                                              momentumflux_output_conversion(k+ki) *              &
     936                                              surf_def_h(l)%usws(m) * rmask(j,i,sr)     ! w"u"
    993937                         !$ACC ATOMIC
    994                          sums_l(k+ki,14,tn) = sums_l(k+ki,14,tn) + &
    995                                     momentumflux_output_conversion(k+ki) * &
    996                                     surf_def_h(l)%vsws(m) * rmask(j,i,sr)     ! w"v"
     938                         sums_l(k+ki,14,tn) = sums_l(k+ki,14,tn) +                                 &
     939                                              momentumflux_output_conversion(k+ki) *              &
     940                                              surf_def_h(l)%vsws(m) * rmask(j,i,sr)     ! w"v"
    997941                         !$ACC ATOMIC
    998                          sums_l(k+ki,16,tn) = sums_l(k+ki,16,tn) + &
    999                                     heatflux_output_conversion(k+ki) * &
    1000                                     surf_def_h(l)%shf(m)  * rmask(j,i,sr)     ! w"pt"
     942                         sums_l(k+ki,16,tn) = sums_l(k+ki,16,tn) +                                 &
     943                                              heatflux_output_conversion(k+ki) *                  &
     944                                              surf_def_h(l)%shf(m)  * rmask(j,i,sr)     ! w"pt"
    1001945#if 0
    1002                          sums_l(k+ki,58,tn) = sums_l(k+ki,58,tn) + &
    1003                                     0.0_wp * rmask(j,i,sr)        ! u"pt"
    1004                          sums_l(k+ki,61,tn) = sums_l(k+ki,61,tn) + &
    1005                                     0.0_wp * rmask(j,i,sr)        ! v"pt"
     946                         sums_l(k+ki,58,tn) = sums_l(k+ki,58,tn) +                                 &
     947                                              0.0_wp * rmask(j,i,sr)                    ! u"pt"
     948                         sums_l(k+ki,61,tn) = sums_l(k+ki,61,tn) +                                 &
     949                                              0.0_wp * rmask(j,i,sr)                    ! v"pt"
    1006950#endif
    1007951#ifndef _OPENACC
    1008952                         IF ( ocean_mode )  THEN
    1009                             sums_l(k+ki,65,tn) = sums_l(k+ki,65,tn) + &
    1010                                        surf_def_h(l)%sasws(m) * rmask(j,i,sr)  ! w"sa"
     953                            sums_l(k+ki,65,tn) = sums_l(k+ki,65,tn) +                              &
     954                                                 surf_def_h(l)%sasws(m) * rmask(j,i,sr)  ! w"sa"
    1011955                         ENDIF
    1012956                         IF ( humidity )  THEN
    1013                             sums_l(k+ki,48,tn) = sums_l(k+ki,48,tn) +                     &
    1014                                        waterflux_output_conversion(k+ki) *      &
    1015                                        surf_def_h(l)%qsws(m) * rmask(j,i,sr)  ! w"q" (w"qv")
    1016                             sums_l(k+ki,45,tn) = sums_l(k+ki,45,tn) + (                   &
    1017                                        ( 1.0_wp + 0.61_wp * q(k+ki,j,i) ) *     &
    1018                                        surf_def_h(l)%shf(m) + 0.61_wp * pt(k+ki,j,i) *      &
    1019                                                   surf_def_h(l)%qsws(m) )                  &
    1020                                        * heatflux_output_conversion(k+ki)
     957                            sums_l(k+ki,48,tn) = sums_l(k+ki,48,tn) +                              &
     958                                                 waterflux_output_conversion(k+ki) *               &
     959                                                 surf_def_h(l)%qsws(m) * rmask(j,i,sr)  ! w"q" (w"qv")
     960                            sums_l(k+ki,45,tn) = sums_l(k+ki,45,tn) +  (                           &
     961                                                 ( 1.0_wp + 0.61_wp * q(k+ki,j,i) ) *              &
     962                                                 surf_def_h(l)%shf(m) + 0.61_wp * pt(k+ki,j,i) *   &
     963                                                 surf_def_h(l)%qsws(m) )                           &
     964                                                 * heatflux_output_conversion(k+ki)
    1021965                            IF ( cloud_droplets )  THEN
    1022                                sums_l(k+ki,45,tn) = sums_l(k+ki,45,tn) + (                &
    1023                                          ( 1.0_wp + 0.61_wp * q(k+ki,j,i) -     &
    1024                                            ql(k+ki,j,i) ) * surf_def_h(l)%shf(m) +          &
    1025                                            0.61_wp * pt(k+ki,j,i) * surf_def_h(l)%qsws(m) ) &
    1026                                           * heatflux_output_conversion(k+ki)
     966                               sums_l(k+ki,45,tn) = sums_l(k+ki,45,tn) +      (                    &
     967                                                    ( 1.0_wp + 0.61_wp * q(k+ki,j,i) -             &
     968                                                      ql(k+ki,j,i) ) * surf_def_h(l)%shf(m) +      &
     969                                                      0.61_wp * pt(k+ki,j,i)                       &
     970                                                      * surf_def_h(l)%qsws(m) )                    &
     971                                                    * heatflux_output_conversion(k+ki)
    1027972                            ENDIF
    1028973                            IF ( bulk_cloud_model )  THEN
    1029974!
    1030975!--                            Formula does not work if ql(k+ki) /= 0.0
    1031                                sums_l(k+ki,51,tn) = sums_l(k+ki,51,tn) +                  &
    1032                                           waterflux_output_conversion(k+ki) *   &
    1033                                           surf_def_h(l)%qsws(m) * rmask(j,i,sr) ! w"q" (w"qv")
     976                               sums_l(k+ki,51,tn) = sums_l(k+ki,51,tn) +                           &
     977                                                    waterflux_output_conversion(k+ki) *            &
     978                                                    surf_def_h(l)%qsws(m) * rmask(j,i,sr) ! w"q" (w"qv")
    1034979                            ENDIF
    1035980                         ENDIF
    1036981                         IF ( passive_scalar )  THEN
    1037                             sums_l(k+ki,117,tn) = sums_l(k+ki,117,tn) +                     &
    1038                                         surf_def_h(l)%ssws(m) * rmask(j,i,sr) ! w"s"
     982                            sums_l(k+ki,117,tn) = sums_l(k+ki,117,tn) +                            &
     983                                                  surf_def_h(l)%ssws(m) * rmask(j,i,sr) ! w"s"
    1039984                         ENDIF
    1040985#endif
     
    1044989                   ENDIF
    1045990                ENDDO
    1046                 IF ( surf_lsm_h%end_index(j,i) >=                              &
    1047                      surf_lsm_h%start_index(j,i) )  THEN
     991                IF ( surf_lsm_h%end_index(j,i) >= surf_lsm_h%start_index(j,i) )  THEN
    1048992                   m = surf_lsm_h%start_index(j,i)
    1049993                   !$ACC ATOMIC
    1050                    sums_l(nzb,12,tn) = sums_l(nzb,12,tn) + &
    1051                                     momentumflux_output_conversion(nzb) * &
    1052                                     surf_lsm_h%usws(m) * rmask(j,i,sr)     ! w"u"
     994                   sums_l(nzb,12,tn) = sums_l(nzb,12,tn) +                                         &
     995                                       momentumflux_output_conversion(nzb) *                      &
     996                                       surf_lsm_h%usws(m) * rmask(j,i,sr)     ! w"u"
    1053997                   !$ACC ATOMIC
    1054                    sums_l(nzb,14,tn) = sums_l(nzb,14,tn) + &
    1055                                     momentumflux_output_conversion(nzb) * &
    1056                                     surf_lsm_h%vsws(m) * rmask(j,i,sr)     ! w"v"
     998                   sums_l(nzb,14,tn) = sums_l(nzb,14,tn) +                                         &
     999                                       momentumflux_output_conversion(nzb) *                      &
     1000                                       surf_lsm_h%vsws(m) * rmask(j,i,sr)     ! w"v"
    10571001                   !$ACC ATOMIC
    1058                    sums_l(nzb,16,tn) = sums_l(nzb,16,tn) + &
    1059                                     heatflux_output_conversion(nzb) * &
    1060                                     surf_lsm_h%shf(m)  * rmask(j,i,sr)     ! w"pt"
     1002                   sums_l(nzb,16,tn) = sums_l(nzb,16,tn) +                                         &
     1003                                       heatflux_output_conversion(nzb) *                          &
     1004                                       surf_lsm_h%shf(m)  * rmask(j,i,sr)     ! w"pt"
    10611005#if 0
    1062                    sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + &
    1063                                     0.0_wp * rmask(j,i,sr)        ! u"pt"
    1064                    sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + &
    1065                                     0.0_wp * rmask(j,i,sr)        ! v"pt"
     1006                   sums_l(nzb,58,tn) = sums_l(nzb,58,tn) +                                         &
     1007                                       0.0_wp * rmask(j,i,sr)        ! u"pt"
     1008                   sums_l(nzb,61,tn) = sums_l(nzb,61,tn) +                                         &
     1009                                       0.0_wp * rmask(j,i,sr)        ! v"pt"
    10661010#endif
    10671011#ifndef _OPENACC
    10681012                   IF ( ocean_mode )  THEN
    1069                       sums_l(nzb,65,tn) = sums_l(nzb,65,tn) + &
    1070                                        surf_lsm_h%sasws(m) * rmask(j,i,sr)  ! w"sa"
     1013                      sums_l(nzb,65,tn) = sums_l(nzb,65,tn) +                                      &
     1014                                          surf_lsm_h%sasws(m) * rmask(j,i,sr)  ! w"sa"
    10711015                   ENDIF
    10721016                   IF ( humidity )  THEN
    1073                       sums_l(nzb,48,tn) = sums_l(nzb,48,tn) +                     &
    1074                                        waterflux_output_conversion(nzb) *      &
    1075                                        surf_lsm_h%qsws(m) * rmask(j,i,sr)  ! w"q" (w"qv")
    1076                       sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                   &
    1077                                        ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) *     &
    1078                                        surf_lsm_h%shf(m) + 0.61_wp * pt(nzb,j,i) *      &
    1079                                                   surf_lsm_h%qsws(m) )                  &
    1080                                        * heatflux_output_conversion(nzb)
     1017                      sums_l(nzb,48,tn) = sums_l(nzb,48,tn) +                                      &
     1018                                          waterflux_output_conversion(nzb) *                       &
     1019                                          surf_lsm_h%qsws(m) * rmask(j,i,sr)  ! w"q" (w"qv")
     1020                      sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                                    &
     1021                                          ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) * surf_lsm_h%shf(m) +  &
     1022                                            0.61_wp * pt(nzb,j,i) * surf_lsm_h%qsws(m) )           &
     1023                                          * heatflux_output_conversion(nzb)
    10811024                      IF ( cloud_droplets )  THEN
    1082                          sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                &
    1083                                          ( 1.0_wp + 0.61_wp * q(nzb,j,i) -     &
    1084                                            ql(nzb,j,i) ) * surf_lsm_h%shf(m) +          &
    1085                                            0.61_wp * pt(nzb,j,i) * surf_lsm_h%qsws(m) ) &
    1086                                           * heatflux_output_conversion(nzb)
     1025                         sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                                 &
     1026                                             ( 1.0_wp + 0.61_wp * q(nzb,j,i) -                     &
     1027                                               ql(nzb,j,i) ) * surf_lsm_h%shf(m) +                 &
     1028                                               0.61_wp * pt(nzb,j,i) * surf_lsm_h%qsws(m) )        &
     1029                                             * heatflux_output_conversion(nzb)
    10871030                      ENDIF
    10881031                      IF ( bulk_cloud_model )  THEN
    10891032!
    10901033!--                      Formula does not work if ql(nzb) /= 0.0
    1091                          sums_l(nzb,51,tn) = sums_l(nzb,51,tn) +                  &
    1092                                           waterflux_output_conversion(nzb) *   &
    1093                                           surf_lsm_h%qsws(m) * rmask(j,i,sr) ! w"q" (w"qv")
     1034                         sums_l(nzb,51,tn) = sums_l(nzb,51,tn) +                                   &
     1035                                             waterflux_output_conversion(nzb) *                    &
     1036                                             surf_lsm_h%qsws(m) * rmask(j,i,sr) ! w"q" (w"qv")
    10941037                      ENDIF
    10951038                   ENDIF
    10961039                   IF ( passive_scalar )  THEN
    1097                       sums_l(nzb,117,tn) = sums_l(nzb,117,tn) +                     &
    1098                                         surf_lsm_h%ssws(m) * rmask(j,i,sr) ! w"s"
     1040                      sums_l(nzb,117,tn) = sums_l(nzb,117,tn) +                                    &
     1041                                           surf_lsm_h%ssws(m) * rmask(j,i,sr) ! w"s"
    10991042                   ENDIF
    11001043#endif
    11011044
    11021045                ENDIF
    1103                 IF ( surf_usm_h%end_index(j,i) >=                              &
    1104                      surf_usm_h%start_index(j,i) )  THEN
     1046                IF ( surf_usm_h%end_index(j,i) >= surf_usm_h%start_index(j,i) )  THEN
    11051047                   m = surf_usm_h%start_index(j,i)
    11061048                   !$ACC ATOMIC
    1107                    sums_l(nzb,12,tn) = sums_l(nzb,12,tn) + &
    1108                                     momentumflux_output_conversion(nzb) * &
    1109                                     surf_usm_h%usws(m) * rmask(j,i,sr)     ! w"u"
     1049                   sums_l(nzb,12,tn) = sums_l(nzb,12,tn) +                                         &
     1050                                       momentumflux_output_conversion(nzb) *                      &
     1051                                       surf_usm_h%usws(m) * rmask(j,i,sr)                ! w"u"
    11101052                   !$ACC ATOMIC
    1111                    sums_l(nzb,14,tn) = sums_l(nzb,14,tn) + &
    1112                                     momentumflux_output_conversion(nzb) * &
    1113                                     surf_usm_h%vsws(m) * rmask(j,i,sr)     ! w"v"
     1053                   sums_l(nzb,14,tn) = sums_l(nzb,14,tn) +                                         &
     1054                                       momentumflux_output_conversion(nzb) *                      &
     1055                                       surf_usm_h%vsws(m) * rmask(j,i,sr)                ! w"v"
    11141056                   !$ACC ATOMIC
    1115                    sums_l(nzb,16,tn) = sums_l(nzb,16,tn) + &
    1116                                     heatflux_output_conversion(nzb) * &
    1117                                     surf_usm_h%shf(m)  * rmask(j,i,sr)     ! w"pt"
     1057                   sums_l(nzb,16,tn) = sums_l(nzb,16,tn) +                                         &
     1058                                       heatflux_output_conversion(nzb) *                          &
     1059                                       surf_usm_h%shf(m)  * rmask(j,i,sr)                ! w"pt"
    11181060#if 0
    1119                    sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + &
    1120                                     0.0_wp * rmask(j,i,sr)        ! u"pt"
    1121                    sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + &
    1122                                     0.0_wp * rmask(j,i,sr)        ! v"pt"
     1061                   sums_l(nzb,58,tn) = sums_l(nzb,58,tn) + 0.0_wp * rmask(j,i,sr)        ! u"pt"
     1062                   sums_l(nzb,61,tn) = sums_l(nzb,61,tn) + 0.0_wp * rmask(j,i,sr)        ! v"pt"
    11231063#endif
    11241064#ifndef _OPENACC
    11251065                   IF ( ocean_mode )  THEN
    1126                       sums_l(nzb,65,tn) = sums_l(nzb,65,tn) + &
    1127                                        surf_usm_h%sasws(m) * rmask(j,i,sr)  ! w"sa"
     1066                      sums_l(nzb,65,tn) = sums_l(nzb,65,tn) +                                      &
     1067                                          surf_usm_h%sasws(m) * rmask(j,i,sr)            ! w"sa"
    11281068                   ENDIF
    11291069                   IF ( humidity )  THEN
    1130                       sums_l(nzb,48,tn) = sums_l(nzb,48,tn) +                     &
    1131                                        waterflux_output_conversion(nzb) *      &
    1132                                        surf_usm_h%qsws(m) * rmask(j,i,sr)  ! w"q" (w"qv")
    1133                       sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                   &
    1134                                        ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) *     &
    1135                                        surf_usm_h%shf(m) + 0.61_wp * pt(nzb,j,i) *      &
    1136                                                   surf_usm_h%qsws(m) )                  &
    1137                                        * heatflux_output_conversion(nzb)
     1070                      sums_l(nzb,48,tn) = sums_l(nzb,48,tn) +                                      &
     1071                                          waterflux_output_conversion(nzb) *                       &
     1072                                          surf_usm_h%qsws(m) * rmask(j,i,sr)             ! w"q" (w"qv")
     1073                      sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                                    &
     1074                                          ( 1.0_wp + 0.61_wp * q(nzb,j,i) ) *                      &
     1075                                          surf_usm_h%shf(m) + 0.61_wp * pt(nzb,j,i) *              &
     1076                                          surf_usm_h%qsws(m)  )                                    &
     1077                                          * heatflux_output_conversion(nzb)
    11381078                      IF ( cloud_droplets )  THEN
    1139                          sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                &
    1140                                          ( 1.0_wp + 0.61_wp * q(nzb,j,i) -     &
    1141                                            ql(nzb,j,i) ) * surf_usm_h%shf(m) +          &
    1142                                            0.61_wp * pt(nzb,j,i) * surf_usm_h%qsws(m) ) &
    1143                                           * heatflux_output_conversion(nzb)
     1079                         sums_l(nzb,45,tn) = sums_l(nzb,45,tn) + (                                 &
     1080                                             ( 1.0_wp + 0.61_wp * q(nzb,j,i) -                     &
     1081                                              ql(nzb,j,i) ) * surf_usm_h%shf(m) +                  &
     1082                                              0.61_wp * pt(nzb,j,i) * surf_usm_h%qsws(m) )        &
     1083                                             * heatflux_output_conversion(nzb)
    11441084                      ENDIF
    11451085                      IF ( bulk_cloud_model )  THEN
    11461086!
    11471087!--                      Formula does not work if ql(nzb) /= 0.0
    1148                          sums_l(nzb,51,tn) = sums_l(nzb,51,tn) +                  &
    1149                                           waterflux_output_conversion(nzb) *   &
    1150                                           surf_usm_h%qsws(m) * rmask(j,i,sr) ! w"q" (w"qv")
     1088                         sums_l(nzb,51,tn) = sums_l(nzb,51,tn) +                                   &
     1089                                             waterflux_output_conversion(nzb) *                    &
     1090                                             surf_usm_h%qsws(m) * rmask(j,i,sr)          ! w"q" (w"qv")
    11511091                      ENDIF
    11521092                   ENDIF
    11531093                   IF ( passive_scalar )  THEN
    1154                       sums_l(nzb,117,tn) = sums_l(nzb,117,tn) +                     &
    1155                                         surf_usm_h%ssws(m) * rmask(j,i,sr) ! w"s"
     1094                      sums_l(nzb,117,tn) = sums_l(nzb,117,tn) +                                    &
     1095                                           surf_usm_h%ssws(m) * rmask(j,i,sr)            ! w"s"
    11561096                   ENDIF
    11571097#endif
     
    11631103#ifndef _OPENACC
    11641104             IF ( .NOT. neutral )  THEN
    1165                 IF ( surf_def_h(0)%end_index(j,i) >=                           &
    1166                      surf_def_h(0)%start_index(j,i) )  THEN
     1105                IF ( surf_def_h(0)%end_index(j,i) >= surf_def_h(0)%start_index(j,i) )  THEN
    11671106                   m = surf_def_h(0)%start_index(j,i)
    1168                    sums_l(nzb,112,tn) = sums_l(nzb,112,tn) +                   &
    1169                                         surf_def_h(0)%ol(m)  * rmask(j,i,sr) ! L
    1170                 ENDIF
    1171                 IF ( surf_lsm_h%end_index(j,i) >=                              &
    1172                      surf_lsm_h%start_index(j,i) )  THEN
     1107                   sums_l(nzb,112,tn) = sums_l(nzb,112,tn) + surf_def_h(0)%ol(m) * rmask(j,i,sr) ! L
     1108                ENDIF
     1109                IF ( surf_lsm_h%end_index(j,i) >= surf_lsm_h%start_index(j,i) )  THEN
    11731110                   m = surf_lsm_h%start_index(j,i)
    1174                    sums_l(nzb,112,tn) = sums_l(nzb,112,tn) +                   &
    1175                                         surf_lsm_h%ol(m)  * rmask(j,i,sr) ! L
    1176                 ENDIF
    1177                 IF ( surf_usm_h%end_index(j,i) >=                              &
    1178                      surf_usm_h%start_index(j,i) )  THEN
     1111                   sums_l(nzb,112,tn) = sums_l(nzb,112,tn) + surf_lsm_h%ol(m)    * rmask(j,i,sr) ! L
     1112                ENDIF
     1113                IF ( surf_usm_h%end_index(j,i) >= surf_usm_h%start_index(j,i) )  THEN
    11791114                   m = surf_usm_h%start_index(j,i)
    1180                    sums_l(nzb,112,tn) = sums_l(nzb,112,tn) +                   &
    1181                                         surf_usm_h%ol(m)  * rmask(j,i,sr) ! L
     1115                   sums_l(nzb,112,tn) = sums_l(nzb,112,tn) + surf_usm_h%ol(m)    * rmask(j,i,sr) ! L
    11821116                ENDIF
    11831117             ENDIF
    11841118
    11851119             IF ( radiation )  THEN
    1186                 IF ( surf_def_h(0)%end_index(j,i) >=                           &
    1187                      surf_def_h(0)%start_index(j,i) )  THEN
     1120                IF ( surf_def_h(0)%end_index(j,i) >= surf_def_h(0)%start_index(j,i) )  THEN
    11881121                   m = surf_def_h(0)%start_index(j,i)
    1189                    sums_l(nzb,99,tn)  = sums_l(nzb,99,tn)  +                   &
    1190                                         surf_def_h(0)%rad_net(m) * rmask(j,i,sr)
    1191                    sums_l(nzb,100,tn) = sums_l(nzb,100,tn)  +                  &
    1192                                         surf_def_h(0)%rad_lw_in(m) * rmask(j,i,sr)
    1193                    sums_l(nzb,101,tn) = sums_l(nzb,101,tn)  +                  &
     1122                   sums_l(nzb,99,tn)  = sums_l(nzb,99,tn)   +                                      &
     1123                                        surf_def_h(0)%rad_net(m)    * rmask(j,i,sr)
     1124                   sums_l(nzb,100,tn) = sums_l(nzb,100,tn)  +                                      &
     1125                                        surf_def_h(0)%rad_lw_in(m)  * rmask(j,i,sr)
     1126                   sums_l(nzb,101,tn) = sums_l(nzb,101,tn)  +                                      &
    11941127                                        surf_def_h(0)%rad_lw_out(m) * rmask(j,i,sr)
    1195                    sums_l(nzb,102,tn) = sums_l(nzb,102,tn)  +                  &
    1196                                         surf_def_h(0)%rad_sw_in(m) * rmask(j,i,sr)
    1197                    sums_l(nzb,103,tn) = sums_l(nzb,103,tn)  +                  &
     1128                   sums_l(nzb,102,tn) = sums_l(nzb,102,tn)  +                                      &
     1129                                        surf_def_h(0)%rad_sw_in(m)  * rmask(j,i,sr)
     1130                   sums_l(nzb,103,tn) = sums_l(nzb,103,tn)  +                                      &
    11981131                                        surf_def_h(0)%rad_sw_out(m) * rmask(j,i,sr)
    11991132                ENDIF
    1200                 IF ( surf_lsm_h%end_index(j,i) >=                              &
    1201                      surf_lsm_h%start_index(j,i) )  THEN
     1133                IF ( surf_lsm_h%end_index(j,i) >= surf_lsm_h%start_index(j,i) )  THEN
    12021134                   m = surf_lsm_h%start_index(j,i)
    1203                    sums_l(nzb,99,tn)  = sums_l(nzb,99,tn)  +                   &
    1204                                         surf_lsm_h%rad_net(m) * rmask(j,i,sr)
    1205                    sums_l(nzb,100,tn) = sums_l(nzb,100,tn)  +                  &
    1206                                         surf_lsm_h%rad_lw_in(m) * rmask(j,i,sr)
    1207                    sums_l(nzb,101,tn) = sums_l(nzb,101,tn)  +                  &
     1135                   sums_l(nzb,99,tn)  = sums_l(nzb,99,tn)   +                                      &
     1136                                        surf_lsm_h%rad_net(m)    * rmask(j,i,sr)
     1137                   sums_l(nzb,100,tn) = sums_l(nzb,100,tn)  +                                      &
     1138                                        surf_lsm_h%rad_lw_in(m)  * rmask(j,i,sr)
     1139                   sums_l(nzb,101,tn) = sums_l(nzb,101,tn)  +                                      &
    12081140                                        surf_lsm_h%rad_lw_out(m) * rmask(j,i,sr)
    1209                    sums_l(nzb,102,tn) = sums_l(nzb,102,tn)  +                  &
    1210                                         surf_lsm_h%rad_sw_in(m) * rmask(j,i,sr)
    1211                    sums_l(nzb,103,tn) = sums_l(nzb,103,tn)  +                  &
     1141                   sums_l(nzb,102,tn) = sums_l(nzb,102,tn)  +                                      &
     1142                                        surf_lsm_h%rad_sw_in(m)  * rmask(j,i,sr)
     1143                   sums_l(nzb,103,tn) = sums_l(nzb,103,tn)  +                                      &
    12121144                                        surf_lsm_h%rad_sw_out(m) * rmask(j,i,sr)
    12131145                ENDIF
    1214                 IF ( surf_usm_h%end_index(j,i) >=                              &
    1215                      surf_usm_h%start_index(j,i) )  THEN
     1146                IF ( surf_usm_h%end_index(j,i) >= surf_usm_h%start_index(j,i) )  THEN
    12161147                   m = surf_usm_h%start_index(j,i)
    1217                    sums_l(nzb,99,tn)  = sums_l(nzb,99,tn)  +                   &
    1218                                         surf_usm_h%rad_net(m) * rmask(j,i,sr)
    1219                    sums_l(nzb,100,tn) = sums_l(nzb,100,tn)  +                  &
    1220                                         surf_usm_h%rad_lw_in(m) * rmask(j,i,sr)
    1221                    sums_l(nzb,101,tn) = sums_l(nzb,101,tn)  +                  &
     1148                   sums_l(nzb,99,tn)  = sums_l(nzb,99,tn)   +                                      &
     1149                                        surf_usm_h%rad_net(m)    * rmask(j,i,sr)
     1150                   sums_l(nzb,100,tn) = sums_l(nzb,100,tn)  +                                      &
     1151                                        surf_usm_h%rad_lw_in(m)  * rmask(j,i,sr)
     1152                   sums_l(nzb,101,tn) = sums_l(nzb,101,tn)  +                                      &
    12221153                                        surf_usm_h%rad_lw_out(m) * rmask(j,i,sr)
    1223                    sums_l(nzb,102,tn) = sums_l(nzb,102,tn)  +                  &
    1224                                         surf_usm_h%rad_sw_in(m) * rmask(j,i,sr)
    1225                    sums_l(nzb,103,tn) = sums_l(nzb,103,tn)  +                  &
     1154                   sums_l(nzb,102,tn) = sums_l(nzb,102,tn)  +                                      &
     1155                                        surf_usm_h%rad_sw_in(m)  * rmask(j,i,sr)
     1156                   sums_l(nzb,103,tn) = sums_l(nzb,103,tn)  +                                      &
    12261157                                        surf_usm_h%rad_sw_out(m) * rmask(j,i,sr)
    12271158                ENDIF
     
    12301161                IF ( radiation_scheme == 'rrtmg' )  THEN
    12311162
    1232                    IF ( surf_def_h(0)%end_index(j,i) >=                        &
    1233                         surf_def_h(0)%start_index(j,i) )  THEN
     1163                   IF ( surf_def_h(0)%end_index(j,i) >= surf_def_h(0)%start_index(j,i) )  THEN
    12341164                      m = surf_def_h(0)%start_index(j,i)
    1235                       sums_l(nzb,108,tn)  = sums_l(nzb,108,tn)  +              &
    1236                                    surf_def_h(0)%rrtm_aldif(m,0) * rmask(j,i,sr)
    1237                       sums_l(nzb,109,tn) = sums_l(nzb,109,tn)  +               &
    1238                                    surf_def_h(0)%rrtm_aldir(m,0) * rmask(j,i,sr)
    1239                       sums_l(nzb,110,tn) = sums_l(nzb,110,tn)  +               &
    1240                                    surf_def_h(0)%rrtm_asdif(m,0) * rmask(j,i,sr)
    1241                       sums_l(nzb,111,tn) = sums_l(nzb,111,tn)  +               &
    1242                                    surf_def_h(0)%rrtm_asdir(m,0) * rmask(j,i,sr)
     1165                      sums_l(nzb,108,tn)  = sums_l(nzb,108,tn)  +                                  &
     1166                                            surf_def_h(0)%rrtm_aldif(m,0) * rmask(j,i,sr)
     1167                      sums_l(nzb,109,tn) = sums_l(nzb,109,tn)  +                                   &
     1168                                           surf_def_h(0)%rrtm_aldir(m,0) * rmask(j,i,sr)
     1169                      sums_l(nzb,110,tn) = sums_l(nzb,110,tn)  +                                   &
     1170                                           surf_def_h(0)%rrtm_asdif(m,0) * rmask(j,i,sr)
     1171                      sums_l(nzb,111,tn) = sums_l(nzb,111,tn)  +                                   &
     1172                                           surf_def_h(0)%rrtm_asdir(m,0) * rmask(j,i,sr)
    12431173                   ENDIF
    1244                    IF ( surf_lsm_h%end_index(j,i) >=                           &
    1245                         surf_lsm_h%start_index(j,i) )  THEN
     1174                   IF ( surf_lsm_h%end_index(j,i) >= surf_lsm_h%start_index(j,i) )  THEN
    12461175                      m = surf_lsm_h%start_index(j,i)
    1247                       sums_l(nzb,108,tn)  = sums_l(nzb,108,tn)  +              &
    1248                                SUM( surf_lsm_h%frac(m,:) *                     &
    1249                                     surf_lsm_h%rrtm_aldif(m,:) ) * rmask(j,i,sr)
    1250                       sums_l(nzb,109,tn) = sums_l(nzb,109,tn)  +               &
    1251                                SUM( surf_lsm_h%frac(m,:) *                     &
    1252                                     surf_lsm_h%rrtm_aldir(m,:) ) * rmask(j,i,sr)
    1253                       sums_l(nzb,110,tn) = sums_l(nzb,110,tn)  +               &
    1254                                SUM( surf_lsm_h%frac(m,:) *                     &
    1255                                     surf_lsm_h%rrtm_asdif(m,:) ) * rmask(j,i,sr)
    1256                       sums_l(nzb,111,tn) = sums_l(nzb,111,tn)  +               &
    1257                                SUM( surf_lsm_h%frac(m,:) *                     &
    1258                                     surf_lsm_h%rrtm_asdir(m,:) ) * rmask(j,i,sr)
     1176                      sums_l(nzb,108,tn)  = sums_l(nzb,108,tn)  +                                  &
     1177                                            SUM( surf_lsm_h%frac(m,:) *                            &
     1178                                                 surf_lsm_h%rrtm_aldif(m,:) ) * rmask(j,i,sr)
     1179                      sums_l(nzb,109,tn) = sums_l(nzb,109,tn)  +                                   &
     1180                                           SUM( surf_lsm_h%frac(m,:) *                             &
     1181                                                surf_lsm_h%rrtm_aldir(m,:) ) * rmask(j,i,sr)
     1182                      sums_l(nzb,110,tn) = sums_l(nzb,110,tn)  +                                   &
     1183                                           SUM( surf_lsm_h%frac(m,:) *                             &
     1184                                                surf_lsm_h%rrtm_asdif(m,:) ) * rmask(j,i,sr)
     1185                      sums_l(nzb,111,tn) = sums_l(nzb,111,tn)  +                                   &
     1186                                           SUM( surf_lsm_h%frac(m,:) *                             &
     1187                                                surf_lsm_h%rrtm_asdir(m,:) ) * rmask(j,i,sr)
    12591188                   ENDIF
    1260                    IF ( surf_usm_h%end_index(j,i) >=                           &
    1261                         surf_usm_h%start_index(j,i) )  THEN
     1189                   IF ( surf_usm_h%end_index(j,i) >= surf_usm_h%start_index(j,i) )  THEN
    12621190                      m = surf_usm_h%start_index(j,i)
    1263                       sums_l(nzb,108,tn)  = sums_l(nzb,108,tn)  +              &
    1264                                SUM( surf_usm_h%frac(m,:) *                     &
    1265                                     surf_usm_h%rrtm_aldif(m,:) ) * rmask(j,i,sr)
    1266                       sums_l(nzb,109,tn) = sums_l(nzb,109,tn)  +               &
    1267                                SUM( surf_usm_h%frac(m,:) *                     &
    1268                                     surf_usm_h%rrtm_aldir(m,:) ) * rmask(j,i,sr)
    1269                       sums_l(nzb,110,tn) = sums_l(nzb,110,tn)  +               &
    1270                                SUM( surf_usm_h%frac(m,:) *                     &
    1271                                     surf_usm_h%rrtm_asdif(m,:) ) * rmask(j,i,sr)
    1272                       sums_l(nzb,111,tn) = sums_l(nzb,111,tn)  +               &
    1273                                SUM( surf_usm_h%frac(m,:) *                     &
    1274                                     surf_usm_h%rrtm_asdir(m,:) ) * rmask(j,i,sr)
     1191                      sums_l(nzb,108,tn)  = sums_l(nzb,108,tn)  +                                  &
     1192                                            SUM( surf_usm_h%frac(m,:) *                            &
     1193                                                 surf_usm_h%rrtm_aldif(m,:) ) * rmask(j,i,sr)
     1194                      sums_l(nzb,109,tn) = sums_l(nzb,109,tn)  +                                   &
     1195                                           SUM( surf_usm_h%frac(m,:) *                             &
     1196                                                surf_usm_h%rrtm_aldir(m,:) ) * rmask(j,i,sr)
     1197                      sums_l(nzb,110,tn) = sums_l(nzb,110,tn)  +                                   &
     1198                                           SUM( surf_usm_h%frac(m,:) *                             &
     1199                                                surf_usm_h%rrtm_asdif(m,:) ) * rmask(j,i,sr)
     1200                      sums_l(nzb,111,tn) = sums_l(nzb,111,tn)  +                                   &
     1201                                           SUM( surf_usm_h%frac(m,:) *                             &
     1202                                                surf_usm_h%rrtm_asdir(m,:) ) * rmask(j,i,sr)
    12751203                   ENDIF
    12761204
     
    12841212                m = surf_def_h(2)%start_index(j,i)
    12851213                !$ACC ATOMIC
    1286                 sums_l(nzt,12,tn) = sums_l(nzt,12,tn) + &
    1287                                     momentumflux_output_conversion(nzt) * &
     1214                sums_l(nzt,12,tn) = sums_l(nzt,12,tn) +                                            &
     1215                                    momentumflux_output_conversion(nzt) *                          &
    12881216                                    surf_def_h(2)%usws(m) * rmask(j,i,sr)    ! w"u"
    12891217                !$ACC ATOMIC
    1290                 sums_l(nzt+1,12,tn) = sums_l(nzt+1,12,tn) + &
    1291                                     momentumflux_output_conversion(nzt+1) * &
    1292                                     surf_def_h(2)%usws(m) * rmask(j,i,sr)    ! w"u"
    1293                 !$ACC ATOMIC
    1294                 sums_l(nzt,14,tn) = sums_l(nzt,14,tn) + &
    1295                                     momentumflux_output_conversion(nzt) * &
     1218                sums_l(nzt+1,12,tn) = sums_l(nzt+1,12,tn) +                                        &
     1219                                      momentumflux_output_conversion(nzt+1) *                      &
     1220                                      surf_def_h(2)%usws(m) * rmask(j,i,sr)  ! w"u"
     1221                !$ACC ATOMIC
     1222                sums_l(nzt,14,tn) = sums_l(nzt,14,tn) +                                            &
     1223                                    momentumflux_output_conversion(nzt) *                          &
    12961224                                    surf_def_h(2)%vsws(m) * rmask(j,i,sr)    ! w"v"
    12971225                !$ACC ATOMIC
    1298                 sums_l(nzt+1,14,tn) = sums_l(nzt+1,14,tn) + &
    1299                                     momentumflux_output_conversion(nzt+1) * &
    1300                                     surf_def_h(2)%vsws(m) * rmask(j,i,sr)    ! w"v"
    1301                 !$ACC ATOMIC
    1302                 sums_l(nzt,16,tn) = sums_l(nzt,16,tn) + &
    1303                                     heatflux_output_conversion(nzt) * &
    1304                                     surf_def_h(2)%shf(m)  * rmask(j,i,sr)   ! w"pt"
    1305                 !$ACC ATOMIC
    1306                 sums_l(nzt+1,16,tn) = sums_l(nzt+1,16,tn) + &
    1307                                     heatflux_output_conversion(nzt+1) * &
    1308                                     surf_def_h(2)%shf(m)  * rmask(j,i,sr)   ! w"pt"
     1226                sums_l(nzt+1,14,tn) = sums_l(nzt+1,14,tn) +                                        &
     1227                                      momentumflux_output_conversion(nzt+1) *                      &
     1228                                      surf_def_h(2)%vsws(m) * rmask(j,i,sr)  ! w"v"
     1229                !$ACC ATOMIC
     1230                sums_l(nzt,16,tn) = sums_l(nzt,16,tn) +                                            &
     1231                                    heatflux_output_conversion(nzt) *                              &
     1232                                    surf_def_h(2)%shf(m)  * rmask(j,i,sr)    ! w"pt"
     1233                !$ACC ATOMIC
     1234                sums_l(nzt+1,16,tn) = sums_l(nzt+1,16,tn) +                                        &
     1235                                      heatflux_output_conversion(nzt+1) *                          &
     1236                                      surf_def_h(2)%shf(m)  * rmask(j,i,sr)  ! w"pt"
    13091237#if 0
    1310                 sums_l(nzt:nzt+1,58,tn) = sums_l(nzt:nzt+1,58,tn) + &
    1311                                     0.0_wp * rmask(j,i,sr)        ! u"pt"
    1312                 sums_l(nzt:nzt+1,61,tn) = sums_l(nzt:nzt+1,61,tn) + &
    1313                                     0.0_wp * rmask(j,i,sr)        ! v"pt"
     1238                sums_l(nzt:nzt+1,58,tn) = sums_l(nzt:nzt+1,58,tn) +                                &
     1239                                          0.0_wp * rmask(j,i,sr)             ! u"pt"
     1240                sums_l(nzt:nzt+1,61,tn) = sums_l(nzt:nzt+1,61,tn) +                                &
     1241                                          0.0_wp * rmask(j,i,sr)             ! v"pt"
    13141242#endif
    13151243#ifndef _OPENACC
     
    13191247                ENDIF
    13201248                IF ( humidity )  THEN
    1321                    sums_l(nzt,48,tn) = sums_l(nzt,48,tn) +                     &
    1322                                        waterflux_output_conversion(nzt) *      &
     1249                   sums_l(nzt,48,tn) = sums_l(nzt,48,tn) +                                         &
     1250                                       waterflux_output_conversion(nzt) *                          &
    13231251                                       surf_def_h(2)%qsws(m) * rmask(j,i,sr) ! w"q" (w"qv")
    1324                    sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (                   &
    1325                                        ( 1.0_wp + 0.61_wp * q(nzt,j,i) ) *     &
    1326                                        surf_def_h(2)%shf(m) +                  &
    1327                                        0.61_wp * pt(nzt,j,i) *    &
    1328                                        surf_def_h(2)%qsws(m) )      &
     1252                   sums_l(nzt,45,tn) = sums_l(nzt,45,tn) +   (                                     &
     1253                                       ( 1.0_wp + 0.61_wp * q(nzt,j,i) ) *                         &
     1254                                       surf_def_h(2)%shf(m) +                                      &
     1255                                       0.61_wp * pt(nzt,j,i) *                                     &
     1256                                       surf_def_h(2)%qsws(m) )                                     &
    13291257                                       * heatflux_output_conversion(nzt)
    13301258                   IF ( cloud_droplets )  THEN
    1331                       sums_l(nzt,45,tn) = sums_l(nzt,45,tn) + (                &
    1332                                           ( 1.0_wp + 0.61_wp * q(nzt,j,i) -    &
    1333                                             ql(nzt,j,i) ) *                    &
    1334                                             surf_def_h(2)%shf(m) +             &
    1335                                            0.61_wp * pt(nzt,j,i) *             &
    1336                                            surf_def_h(2)%qsws(m) )&
     1259                      sums_l(nzt,45,tn) = sums_l(nzt,45,tn) +    (                                 &
     1260                                          ( 1.0_wp + 0.61_wp * q(nzt,j,i) -                        &
     1261                                            ql(nzt,j,i) ) *                                        &
     1262                                            surf_def_h(2)%shf(m) +                                 &
     1263                                           0.61_wp * pt(nzt,j,i) *                                 &
     1264                                           surf_def_h(2)%qsws(m) )                                 &
    13371265                                           * heatflux_output_conversion(nzt)
    13381266                   ENDIF
     
    13401268!
    13411269!--                   Formula does not work if ql(nzb) /= 0.0
    1342                       sums_l(nzt,51,tn) = sums_l(nzt,51,tn) + &   ! w"q" (w"qv")
    1343                                           waterflux_output_conversion(nzt) *   &
     1270                      sums_l(nzt,51,tn) = sums_l(nzt,51,tn) +              &  ! w"q" (w"qv")
     1271                                          waterflux_output_conversion(nzt) *                       &
    13441272                                          surf_def_h(2)%qsws(m) * rmask(j,i,sr)
    13451273                   ENDIF
    13461274                ENDIF
    13471275                IF ( passive_scalar )  THEN
    1348                    sums_l(nzt,117,tn) = sums_l(nzt,117,tn) + &
     1276                   sums_l(nzt,117,tn) = sums_l(nzt,117,tn) +                                       &
    13491277                                        surf_def_h(2)%ssws(m) * rmask(j,i,sr) ! w"s"
    13501278                ENDIF
     
    13541282!
    13551283!--          Resolved fluxes (can be computed for all horizontal points)
    1356 !--          NOTE: for simplicity, nzb_s_inner is used below, although strictly
    1357 !--          ----  speaking the following k-loop would have to be split up and
    1358 !--                rearranged according to the staggered grid.
     1284!--          NOTE: for simplicity, nzb_s_inner is used below, although strictly speaking the
     1285!--          ----  following k-loop would have to be split up and rearranged according to the
     1286!--                staggered grid.
    13591287             DO  k = nzb, nzt
    13601288                flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 22 ) )
    1361                 ust = 0.5_wp * ( u(k,j,i)   - hom(k,1,1,sr) +                  &
     1289                ust = 0.5_wp * ( u(k,j,i)   - hom(k,1,1,sr) +                                      &
    13621290                                 u(k+1,j,i) - hom(k+1,1,1,sr) )
    1363                 vst = 0.5_wp * ( v(k,j,i)   - hom(k,1,2,sr) +                  &
     1291                vst = 0.5_wp * ( v(k,j,i)   - hom(k,1,2,sr) +                                      &
    13641292                                 v(k+1,j,i) - hom(k+1,1,2,sr) )
    1365                 pts = 0.5_wp * ( pt(k,j,i)   - hom(k,1,4,sr) +                 &
     1293                pts = 0.5_wp * ( pt(k,j,i)   - hom(k,1,4,sr) +                                     &
    13661294                                 pt(k+1,j,i) - hom(k+1,1,4,sr) )
    1367 
     1295!
    13681296!--             Higher moments
    13691297                !$ACC ATOMIC
    1370                 sums_l(k,35,tn) = sums_l(k,35,tn) + pts * w(k,j,i)**2 *        &
    1371                                                     rmask(j,i,sr) * flag
    1372                 !$ACC ATOMIC
    1373                 sums_l(k,36,tn) = sums_l(k,36,tn) + pts**2 * w(k,j,i) *        &
    1374                                                     rmask(j,i,sr) * flag
    1375 
    1376 !
    1377 !--             Salinity flux and density (density does not belong to here,
    1378 !--             but so far there is no other suitable place to calculate)
     1298                sums_l(k,35,tn) = sums_l(k,35,tn) + pts * w(k,j,i)**2 * rmask(j,i,sr) * flag
     1299                !$ACC ATOMIC
     1300                sums_l(k,36,tn) = sums_l(k,36,tn) + pts**2 * w(k,j,i) * rmask(j,i,sr) * flag
     1301
     1302!
     1303!--             Salinity flux and density (density does not belong to here, but so far there is no
     1304!--             other suitable place to calculate)
    13791305#ifndef _OPENACC
    13801306                IF ( ocean_mode )  THEN
    13811307                   IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
    1382                       pts = 0.5_wp * ( sa(k,j,i)   - hom(k,1,23,sr) +          &
     1308                      pts = 0.5_wp * ( sa(k,j,i)   - hom(k,1,23,sr) +                              &
    13831309                                       sa(k+1,j,i) - hom(k+1,1,23,sr) )
    1384                       sums_l(k,66,tn) = sums_l(k,66,tn) + pts * w(k,j,i) *     &
     1310                      sums_l(k,66,tn) = sums_l(k,66,tn) + pts * w(k,j,i) *                         &
    13851311                                        rmask(j,i,sr) * flag
    13861312                   ENDIF
    1387                    sums_l(k,64,tn) = sums_l(k,64,tn) + rho_ocean(k,j,i) *      &
    1388                                                        rmask(j,i,sr) * flag
    1389                    sums_l(k,71,tn) = sums_l(k,71,tn) + prho(k,j,i) *           &
    1390                                                        rmask(j,i,sr) * flag
    1391                 ENDIF
    1392 
    1393 !
    1394 !--             Buoyancy flux, water flux, humidity flux, liquid water
    1395 !--             content, rain drop concentration and rain water content
     1313                   sums_l(k,64,tn) = sums_l(k,64,tn) + rho_ocean(k,j,i) * rmask(j,i,sr) * flag
     1314                   sums_l(k,71,tn) = sums_l(k,71,tn) + prho(k,j,i)      * rmask(j,i,sr) * flag
     1315                ENDIF
     1316
     1317!
     1318!--             Buoyancy flux, water flux, humidity flux, liquid water content, rain drop
     1319!--             concentration and rain water content
    13961320                IF ( humidity )  THEN
    1397                    IF ( bulk_cloud_model .OR. cloud_droplets )  THEN
    1398                       pts = 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +         &
    1399                                     vpt(k+1,j,i) - hom(k+1,1,44,sr) )
    1400                       sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) *     &
    1401                                                rho_air_zw(k) *                 &
    1402                                                heatflux_output_conversion(k) * &
     1321                   IF ( bulk_cloud_model  .OR. cloud_droplets )  THEN
     1322                      pts = 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +                             &
     1323                                       vpt(k+1,j,i) - hom(k+1,1,44,sr) )
     1324                      sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) *                         &
     1325                                                          rho_air_zw(k) *                          &
     1326                                                          heatflux_output_conversion(k) *          &
    14031327                                                          rmask(j,i,sr) * flag
    1404                       sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) * rmask(j,i,sr) &
    1405                                                                     * flag
     1328                      sums_l(k,54,tn) = sums_l(k,54,tn) + ql(k,j,i) * rmask(j,i,sr) * flag
    14061329
    14071330                      IF ( .NOT. cloud_droplets )  THEN
    1408                          pts = 0.5_wp *                                        &
    1409                               ( ( q(k,j,i) - ql(k,j,i) ) -                     &
    1410                               hom(k,1,42,sr) +                                 &
    1411                               ( q(k+1,j,i) - ql(k+1,j,i) ) -                   &
    1412                               hom(k+1,1,42,sr) )
    1413                          sums_l(k,52,tn) = sums_l(k,52,tn) + pts * w(k,j,i) *  &
    1414                                              rho_air_zw(k) *                   &
    1415                                              waterflux_output_conversion(k) *  &
    1416                                                              rmask(j,i,sr)  *  &
    1417                                                              flag
    1418                          sums_l(k,75,tn) = sums_l(k,75,tn) + qc(k,j,i) *       &
    1419                                                              rmask(j,i,sr) *   &
    1420                                                              flag
    1421                          sums_l(k,76,tn) = sums_l(k,76,tn) + prr(k,j,i) *      &
    1422                                                              rmask(j,i,sr) *   &
    1423                                                              flag
     1331                         pts = 0.5_wp *                                                            &
     1332                               ( ( q(k,j,i) - ql(k,j,i) ) -                                        &
     1333                               hom(k,1,42,sr) +                                                    &
     1334                               ( q(k+1,j,i) - ql(k+1,j,i) ) -                                      &
     1335                               hom(k+1,1,42,sr) )
     1336                         sums_l(k,52,tn) = sums_l(k,52,tn) + pts * w(k,j,i) *                      &
     1337                                             rho_air_zw(k) *                                       &
     1338                                             waterflux_output_conversion(k) *                      &
     1339                                             rmask(j,i,sr) * flag
     1340                         sums_l(k,75,tn) = sums_l(k,75,tn) + qc(k,j,i)  * rmask(j,i,sr) * flag
     1341                         sums_l(k,76,tn) = sums_l(k,76,tn) + prr(k,j,i) * rmask(j,i,sr) * flag
    14241342                         IF ( microphysics_morrison )  THEN
    1425                             sums_l(k,123,tn) = sums_l(k,123,tn) + nc(k,j,i) *  &
    1426                                                                 rmask(j,i,sr) *&
    1427                                                                 flag
     1343                            sums_l(k,123,tn) = sums_l(k,123,tn) + nc(k,j,i) * rmask(j,i,sr) * flag
    14281344                         ENDIF
    14291345                         IF ( microphysics_ice_phase )  THEN
    1430                             sums_l(k,124,tn) = sums_l(k,124,tn) + ni(k,j,i) *  &
    1431                                                                 rmask(j,i,sr) *&
    1432                                                                 flag
    1433                             sums_l(k,125,tn) = sums_l(k,125,tn) + qi(k,j,i) *  &
    1434                                                                 rmask(j,i,sr) *&
    1435                                                                 flag
     1346                            sums_l(k,124,tn) = sums_l(k,124,tn) + ni(k,j,i) * rmask(j,i,sr) * flag
     1347                            sums_l(k,125,tn) = sums_l(k,125,tn) + qi(k,j,i) * rmask(j,i,sr) * flag
    14361348                         ENDIF
    14371349
    14381350                         IF ( microphysics_seifert )  THEN
    1439                             sums_l(k,73,tn) = sums_l(k,73,tn) + nr(k,j,i) *    &
    1440                                                                 rmask(j,i,sr) *&
    1441                                                                 flag
    1442                             sums_l(k,74,tn) = sums_l(k,74,tn) + qr(k,j,i) *    &
    1443                                                                 rmask(j,i,sr) *&
    1444                                                                 flag
     1351                            sums_l(k,73,tn) = sums_l(k,73,tn) + nr(k,j,i) * rmask(j,i,sr) * flag
     1352                            sums_l(k,74,tn) = sums_l(k,74,tn) + qr(k,j,i) * rmask(j,i,sr) * flag
    14451353                         ENDIF
    14461354                      ENDIF
    14471355
    14481356                   ELSE
    1449                       IF( .NOT. ws_scheme_sca .OR. sr /= 0 )  THEN
    1450                          pts = 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +      &
     1357                      IF( .NOT. ws_scheme_sca  .OR. sr /= 0 )  THEN
     1358                         pts = 0.5_wp * ( vpt(k,j,i)   - hom(k,1,44,sr) +                          &
    14511359                                          vpt(k+1,j,i) - hom(k+1,1,44,sr) )
    1452                          sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) *  &
    1453                                               rho_air_zw(k) *                  &
    1454                                               heatflux_output_conversion(k) *  &
    1455                                                              rmask(j,i,sr)  *  &
    1456                                                              flag
    1457                       ELSE IF ( ws_scheme_sca .AND. sr == 0 )  THEN
    1458                          sums_l(k,46,tn) = ( ( 1.0_wp + 0.61_wp *              &
    1459                                                hom(k,1,41,sr) ) *              &
    1460                                              sums_l(k,17,tn) +                 &
    1461                                              0.61_wp * hom(k,1,4,sr) *         &
    1462                                              sums_l(k,49,tn)                   &
    1463                                            ) * heatflux_output_conversion(k) * &
    1464                                                flag
     1360                         sums_l(k,46,tn) = sums_l(k,46,tn) + pts * w(k,j,i) *                      &
     1361                                                             rho_air_zw(k)  *                      &
     1362                                                             heatflux_output_conversion(k) *       &
     1363                                                             rmask(j,i,sr)  * flag
     1364                      ELSE IF ( ws_scheme_sca  .AND.  sr == 0 )  THEN
     1365                         sums_l(k,46,tn) = ( ( 1.0_wp + 0.61_wp *                                  &
     1366                                               hom(k,1,41,sr) ) *                                  &
     1367                                             sums_l(k,17,tn) +                                     &
     1368                                             0.61_wp * hom(k,1,4,sr) *                             &
     1369                                             sums_l(k,49,tn)                                       &
     1370                                           ) * heatflux_output_conversion(k) * flag
    14651371                      END IF
    14661372                   END IF
     
    14681374!
    14691375!--             Passive scalar flux
    1470                 IF ( passive_scalar .AND. ( .NOT. ws_scheme_sca                &
    1471                      .OR. sr /= 0 ) )  THEN
    1472                    pts = 0.5_wp * ( s(k,j,i)   - hom(k,1,115,sr) +             &
     1376                IF ( passive_scalar  .AND.  ( .NOT. ws_scheme_sca .OR. sr /= 0 ) )  THEN
     1377                   pts = 0.5_wp * ( s(k,j,i)   - hom(k,1,115,sr) +                                 &
    14731378                                    s(k+1,j,i) - hom(k+1,1,115,sr) )
    1474                    sums_l(k,114,tn) = sums_l(k,114,tn) + pts * w(k,j,i) *      &
    1475                                                        rmask(j,i,sr) * flag
     1379                   sums_l(k,114,tn) = sums_l(k,114,tn) + pts * w(k,j,i) * rmask(j,i,sr) * flag
    14761380                ENDIF
    14771381#endif
     
    14811385!--             has to be adjusted
    14821386                !$ACC ATOMIC
    1483                 sums_l(k,37,tn) = sums_l(k,37,tn) + w(k,j,i) * 0.5_wp *        &
    1484                                              ( ust**2 + vst**2 + w(k,j,i)**2 ) &
    1485                                            * rho_air_zw(k)                     &
    1486                                            * momentumflux_output_conversion(k) &
    1487                                            * rmask(j,i,sr) * flag
     1387                sums_l(k,37,tn) = sums_l(k,37,tn) + w(k,j,i) * 0.5_wp *                            &
     1388                                                    ( ust**2 + vst**2 + w(k,j,i)**2 )              &
     1389                                                    * rho_air_zw(k)                                &
     1390                                                    * momentumflux_output_conversion(k)            &
     1391                                                    * rmask(j,i,sr) * flag
    14881392             ENDDO
    14891393          ENDDO
     
    15051409             j = surf_lsm_h%j(m)
    15061410
    1507              IF ( i >= nxl  .AND.  i <= nxr  .AND.                             &
    1508                   j >= nys  .AND.  j <= nyn )  THEN
     1411             IF ( i >= nxl  .AND.  i <= nxr  .AND.  j >= nys  .AND.  j <= nyn )  THEN
    15091412                sums_l(nzb,93,tn)  = sums_l(nzb,93,tn) + surf_lsm_h%ghf(m)       * rmask(j,i,sr)
    15101413                sums_l(nzb,94,tn)  = sums_l(nzb,94,tn) + surf_lsm_h%qsws_liq(m)  * rmask(j,i,sr)
     
    15261429             j = surf_lsm_h%j(m)
    15271430
    1528              IF ( i >= nxl  .AND.  i <= nxr  .AND.                             &
    1529                   j >= nys  .AND.  j <= nyn )  THEN
     1431             IF ( i >= nxl  .AND.  i <= nxr  .AND.  j >= nys  .AND.  j <= nyn )  THEN
    15301432
    15311433                DO  k = nzb_soil, nzt_soil
    1532                    sums_l(k,89,tn)  = sums_l(k,89,tn)  + t_soil_h%var_2d(k,m)  &
    1533                                       * rmask(j,i,sr)
    1534                    sums_l(k,91,tn)  = sums_l(k,91,tn)  + m_soil_h%var_2d(k,m)  &
    1535                                       * rmask(j,i,sr)
     1434                   sums_l(k,89,tn)  = sums_l(k,89,tn)  + t_soil_h%var_2d(k,m) * rmask(j,i,sr)
     1435                   sums_l(k,91,tn)  = sums_l(k,91,tn)  + m_soil_h%var_2d(k,m) * rmask(j,i,sr)
    15361436                ENDDO
    15371437             ENDIF
     
    15401440       ENDIF
    15411441!
    1542 !--    For speed optimization fluxes which have been computed in part directly
    1543 !--    inside the WS advection routines are treated seperatly
     1442!--    For speed optimization fluxes which have been computed in part directly inside the WS
     1443!--    advection routines are treated seperatly.
    15441444!--    Momentum fluxes first:
    15451445
     
    15531453                DO  k = nzb, nzt
    15541454!
    1555 !--                Flag 23 is used to mask surface fluxes as well as model-top
    1556 !--                fluxes, which are added further below.
    1557                    flag = MERGE( 1.0_wp, 0.0_wp,                               &
    1558                                  BTEST( wall_flags_total_0(k,j,i), 23 ) ) *    &
    1559                           MERGE( 1.0_wp, 0.0_wp,                               &
    1560                                  BTEST( wall_flags_total_0(k,j,i), 9  ) )
    1561 
    1562                    ust = 0.5_wp * ( u(k,j,i)   - hom(k,1,1,sr) +               &
     1455!--                Flag 23 is used to mask surface fluxes as well as model-top fluxes, which are
     1456!--                added further below.
     1457                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 23 ) ) *        &
     1458                          MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 9  ) )
     1459
     1460                   ust = 0.5_wp * ( u(k,j,i)   - hom(k,1,1,sr) +                                   &
    15631461                                    u(k+1,j,i) - hom(k+1,1,1,sr) )
    1564                    vst = 0.5_wp * ( v(k,j,i)   - hom(k,1,2,sr) +               &
     1462                   vst = 0.5_wp * ( v(k,j,i)   - hom(k,1,2,sr) +                                   &
    15651463                                    v(k+1,j,i) - hom(k+1,1,2,sr) )
    15661464!
    15671465!--                Momentum flux w*u*
    1568                    sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5_wp *                &
    1569                                                      ( w(k,j,i-1) + w(k,j,i) ) &
    1570                                            * rho_air_zw(k)                     &
    1571                                            * momentumflux_output_conversion(k) &
    1572                                                      * ust * rmask(j,i,sr)     &
     1466                   sums_l(k,13,tn) = sums_l(k,13,tn) + 0.5_wp *                                    &
     1467                                                     ( w(k,j,i-1) + w(k,j,i) )                     &
     1468                                                     * rho_air_zw(k)                               &
     1469                                                     * momentumflux_output_conversion(k)          &
     1470                                                     * ust * rmask(j,i,sr)                         &
    15731471                                                           * flag
    15741472!
    15751473!--                Momentum flux w*v*
    1576                    sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5_wp *                &
    1577                                                      ( w(k,j-1,i) + w(k,j,i) ) &
    1578                                            * rho_air_zw(k)                     &
    1579                                            * momentumflux_output_conversion(k) &
    1580                                                      * vst * rmask(j,i,sr)     &
     1474                   sums_l(k,15,tn) = sums_l(k,15,tn) + 0.5_wp * ( w(k,j-1,i) + w(k,j,i) )          &
     1475                                                     * rho_air_zw(k)                               &
     1476                                                     * momentumflux_output_conversion(k)           &
     1477                                                     * vst * rmask(j,i,sr)                         &
    15811478                                                           * flag
    15821479                ENDDO
     
    15901487             DO  j = nys, nyn
    15911488                DO  k = nzb, nzt
    1592                    flag = MERGE( 1.0_wp, 0.0_wp,                               &
    1593                                  BTEST( wall_flags_total_0(k,j,i), 23 ) ) *    &
    1594                           MERGE( 1.0_wp, 0.0_wp,                               &
    1595                                  BTEST( wall_flags_total_0(k,j,i), 9  ) )
     1489                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 23 ) ) *        &
     1490                          MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 9  ) )
    15961491!
    15971492!--                Vertical heat flux
    1598                    sums_l(k,17,tn) = sums_l(k,17,tn) + 0.5_wp *                &
    1599                            ( pt(k,j,i)   - hom(k,1,4,sr) +                     &
    1600                              pt(k+1,j,i) - hom(k+1,1,4,sr) )                   &
    1601                            * rho_air_zw(k)                                     &
    1602                            * heatflux_output_conversion(k)                     &
    1603                            * w(k,j,i) * rmask(j,i,sr) * flag
     1493                   sums_l(k,17,tn) = sums_l(k,17,tn) + 0.5_wp *                                    &
     1494                                     ( pt(k,j,i)   - hom(k,1,4,sr) +                               &
     1495                                       pt(k+1,j,i) - hom(k+1,1,4,sr) )                             &
     1496                                     * rho_air_zw(k)                                               &
     1497                                     * heatflux_output_conversion(k)                               &
     1498                                     * w(k,j,i) * rmask(j,i,sr) * flag
    16041499                   IF ( humidity )  THEN
    1605                       pts = 0.5_wp * ( q(k,j,i)   - hom(k,1,41,sr) +           &
    1606                                       q(k+1,j,i) - hom(k+1,1,41,sr) )
    1607                       sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) *     &
    1608                                        rho_air_zw(k) *                         &
    1609                                        waterflux_output_conversion(k) *        &
    1610                                        rmask(j,i,sr) * flag
     1500                      pts = 0.5_wp * ( q(k,j,i)   - hom(k,1,41,sr) +                               &
     1501                                       q(k+1,j,i) - hom(k+1,1,41,sr) )
     1502                      sums_l(k,49,tn) = sums_l(k,49,tn) + pts * w(k,j,i) *                         &
     1503                                                          rho_air_zw(k) *                         &
     1504                                                          waterflux_output_conversion(k) *         &
     1505                                                          rmask(j,i,sr) * flag
    16111506                   ENDIF
    16121507                   IF ( passive_scalar )  THEN
    1613                       pts = 0.5_wp * ( s(k,j,i)   - hom(k,1,115,sr) +           &
    1614                                       s(k+1,j,i) - hom(k+1,1,115,sr) )
    1615                       sums_l(k,114,tn) = sums_l(k,114,tn) + pts * w(k,j,i) *    &
    1616                                         rmask(j,i,sr) * flag
     1508                      pts = 0.5_wp * ( s(k,j,i)   - hom(k,1,115,sr) +                              &
     1509                                       s(k+1,j,i) - hom(k+1,1,115,sr) )
     1510                      sums_l(k,114,tn) = sums_l(k,114,tn) + pts * w(k,j,i) *  rmask(j,i,sr) * flag
    16171511                   ENDIF
    16181512                ENDDO
     
    16301524
    16311525!
    1632 !--    Divergence of vertical flux of resolved scale energy and pressure
    1633 !--    fluctuations as well as flux of pressure fluctuation itself (68).
     1526!--    Divergence of vertical flux of resolved scale energy and pressure fluctuations as well as
     1527!--    flux of pressure fluctuation itself (68).
    16341528!--    First calculate the products, then the divergence.
    16351529!--    Calculation is time consuming. Do it only, if profiles shall be plotted.
    1636        IF ( hom(nzb+1,2,55,0) /= 0.0_wp  .OR.  hom(nzb+1,2,68,0) /= 0.0_wp )   &
    1637        THEN
     1530       IF ( hom(nzb+1,2,55,0) /= 0.0_wp  .OR.  hom(nzb+1,2,68,0) /= 0.0_wp )  THEN
    16381531          sums_ll = 0.0_wp  ! local array
    16391532
     
    16441537                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
    16451538
    1646                    sums_ll(k,1) = sums_ll(k,1) + 0.5_wp * w(k,j,i) * (         &
    1647                   ( 0.25_wp * ( u(k,j,i)+u(k+1,j,i)+u(k,j,i+1)+u(k+1,j,i+1) )  &
    1648                             - 0.5_wp * ( hom(k,1,1,sr) + hom(k+1,1,1,sr) ) )**2&
    1649                 + ( 0.25_wp * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i) )  &
    1650                             - 0.5_wp * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) ) )**2&
    1651                 + w(k,j,i)**2                                        ) * flag * rmask(j,i,sr)
    1652 
    1653                    sums_ll(k,2) = sums_ll(k,2) + 0.5_wp * w(k,j,i)             &
    1654                                        * ( ( p(k,j,i) + p(k+1,j,i) )           &
    1655                                          / momentumflux_output_conversion(k) ) &
     1539                   sums_ll(k,1) = sums_ll(k,1) + 0.5_wp * w(k,j,i) * (                             &
     1540                                    ( 0.25_wp * ( u(k,j,i)+u(k+1,j,i)+u(k,j,i+1)+u(k+1,j,i+1) )    &
     1541                                              - 0.5_wp * ( hom(k,1,1,sr) + hom(k+1,1,1,sr) ) )**2  &
     1542                                  + ( 0.25_wp * ( v(k,j,i)+v(k+1,j,i)+v(k,j+1,i)+v(k+1,j+1,i) )    &
     1543                                              - 0.5_wp * ( hom(k,1,2,sr) + hom(k+1,1,2,sr) ) )**2  &
     1544                                  + w(k,j,i)**2                      ) * flag * rmask(j,i,sr)
     1545
     1546                   sums_ll(k,2) = sums_ll(k,2) + 0.5_wp * w(k,j,i)                                 &
     1547                                       * ( ( p(k,j,i) + p(k+1,j,i) )                               &
     1548                                         / momentumflux_output_conversion(k) )                     &
    16561549                                       * flag * rmask(j,i,sr)
    16571550
     
    16771570!
    16781571!--    Divergence of vertical flux of SGS TKE and the flux itself (69)
    1679        IF ( hom(nzb+1,2,57,0) /= 0.0_wp  .OR.  hom(nzb+1,2,69,0) /= 0.0_wp )   &
    1680        THEN
     1572       IF ( hom(nzb+1,2,57,0) /= 0.0_wp  .OR.  hom(nzb+1,2,69,0) /= 0.0_wp )  THEN
    16811573          !$OMP DO
    16821574          DO  i = nxl, nxr
     
    16861578                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
    16871579
    1688                    sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5_wp * (              &
    1689                    (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
    1690                  - (km(k-1,j,i)+km(k,j,i)) * (e(k,j,i)-e(k-1,j,i)) * ddzu(k)   &
    1691                                                                 ) * ddzw(k)    &
     1580                   sums_l(k,57,tn) = sums_l(k,57,tn) - 0.5_wp * (                                  &
     1581                                       (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
     1582                                     - (km(k-1,j,i)+km(k,j,i)) * (e(k,j,i)-e(k-1,j,i)) * ddzu(k)   &
     1583                                                                ) * ddzw(k)                        &
    16921584                                                                  * flag * rmask(j,i,sr)
    16931585
    1694                    sums_l(k,69,tn) = sums_l(k,69,tn) - 0.5_wp * (              &
    1695                    (km(k,j,i)+km(k+1,j,i)) * (e(k+1,j,i)-e(k,j,i)) * ddzu(k+1) &
    1696                                                                 )  * flag * rmask(j,i,sr)
     1586                   sums_l(k,69,tn) = sums_l(k,69,tn) - 0.5_wp * (                                  &
     1587                                        ( km(k,j,i) + km(k+1,j,i) ) *                              &
     1588                                        ( e(k+1,j,i) - e(k,j,i) ) * ddzu(k+1)                      &
     1589                                                                ) * flag * rmask(j,i,sr)
    16971590
    16981591                ENDDO
     
    17161609!
    17171610!--                Subgrid horizontal heat fluxes u"pt", v"pt"
    1718                    sums_l(k,58,tn) = sums_l(k,58,tn) - 0.5_wp *                &
    1719                                                    ( kh(k,j,i) + kh(k,j,i-1) ) &
    1720                                                  * ( pt(k,j,i-1) - pt(k,j,i) ) &
    1721                                                * rho_air_zw(k)                 &
    1722                                                * heatflux_output_conversion(k) &
    1723                                                  * ddx * rmask(j,i,sr) * flag
    1724                    sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5_wp *                &
    1725                                                    ( kh(k,j,i) + kh(k,j-1,i) ) &
    1726                                                  * ( pt(k,j-1,i) - pt(k,j,i) ) &
    1727                                                * rho_air_zw(k)                 &
    1728                                                * heatflux_output_conversion(k) &
    1729                                                  * ddy * rmask(j,i,sr) * flag
     1611                   sums_l(k,58,tn) = sums_l(k,58,tn) - 0.5_wp *                                    &
     1612                                                        ( kh(k,j,i) + kh(k,j,i-1) )                &
     1613                                                      * ( pt(k,j,i-1) - pt(k,j,i) )                &
     1614                                                      * rho_air_zw(k)                              &
     1615                                                      * heatflux_output_conversion(k)              &
     1616                                                      * ddx * rmask(j,i,sr) * flag
     1617                   sums_l(k,61,tn) = sums_l(k,61,tn) - 0.5_wp *                                    &
     1618                                                        ( kh(k,j,i) + kh(k,j-1,i) )                &
     1619                                                      * ( pt(k,j-1,i) - pt(k,j,i) )                &
     1620                                                      * rho_air_zw(k)                              &
     1621                                                      * heatflux_output_conversion(k)              &
     1622                                                      * ddy * rmask(j,i,sr) * flag
    17301623!
    17311624!--                Resolved horizontal heat fluxes u*pt*, v*pt*
    1732                    sums_l(k,59,tn) = sums_l(k,59,tn) +                         &
    1733                                                   ( u(k,j,i) - hom(k,1,1,sr) ) &
    1734                                     * 0.5_wp * ( pt(k,j,i-1) - hom(k,1,4,sr) + &
    1735                                                  pt(k,j,i)   - hom(k,1,4,sr) ) &
    1736                                                * heatflux_output_conversion(k) &
    1737                                                * flag
    1738                    pts = 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) +              &
     1625                   sums_l(k,59,tn) = sums_l(k,59,tn) +              ( u(k,j,i) - hom(k,1,1,sr) )   &
     1626                                                      * 0.5_wp * ( pt(k,j,i-1) - hom(k,1,4,sr) +   &
     1627                                                                   pt(k,j,i)   - hom(k,1,4,sr) )   &
     1628                                                                 * heatflux_output_conversion(k)   &
     1629                                                                 * flag
     1630                   pts = 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) +                                  &
    17391631                                    pt(k,j,i)   - hom(k,1,4,sr) )
    1740                    sums_l(k,62,tn) = sums_l(k,62,tn) +                         &
    1741                                                   ( v(k,j,i) - hom(k,1,2,sr) ) &
    1742                                     * 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) + &
    1743                                                  pt(k,j,i)   - hom(k,1,4,sr) ) &
    1744                                                * heatflux_output_conversion(k) &
    1745                                                * flag
     1632                   sums_l(k,62,tn) = sums_l(k,62,tn) +              ( v(k,j,i) - hom(k,1,2,sr) )   &
     1633                                                      * 0.5_wp * ( pt(k,j-1,i) - hom(k,1,4,sr) +   &
     1634                                                                   pt(k,j,i)   - hom(k,1,4,sr) )   &
     1635                                                                 * heatflux_output_conversion(k)   &
     1636                                                                 * flag
    17461637                ENDDO
    17471638             ENDDO
     
    17731664          ENDIF
    17741665
    1775           fac = ( simulated_time - dt_3d - time_vert(nt) )                     &
    1776                 / ( time_vert(nt+1)-time_vert(nt) )
     1666          fac = ( simulated_time - dt_3d - time_vert(nt) ) / ( time_vert(nt+1)-time_vert(nt) )
    17771667
    17781668
    17791669          DO  k = nzb, nzt
    1780              sums_ls_l(k,0) = td_lsa_lpt(k,nt)                                 &
    1781                               + fac * ( td_lsa_lpt(k,nt+1) - td_lsa_lpt(k,nt) )
    1782              sums_ls_l(k,1) = td_lsa_q(k,nt)                                   &
    1783                               + fac * ( td_lsa_q(k,nt+1) - td_lsa_q(k,nt) )
     1670             sums_ls_l(k,0) = td_lsa_lpt(k,nt) + fac * ( td_lsa_lpt(k,nt+1) - td_lsa_lpt(k,nt) )
     1671             sums_ls_l(k,1) = td_lsa_q(k,nt)   + fac * ( td_lsa_q(k,nt+1)   - td_lsa_q(k,nt) )
    17841672          ENDDO
    17851673
     
    17871675          sums_ls_l(nzt+1,1) = sums_ls_l(nzt,1)
    17881676
    1789           IF ( large_scale_subsidence .AND. use_subsidence_tendencies )  THEN
     1677          IF ( large_scale_subsidence  .AND. use_subsidence_tendencies )  THEN
    17901678
    17911679             DO  k = nzb, nzt
    1792                 sums_ls_l(k,2) = td_sub_lpt(k,nt) + fac *                      &
    1793                                  ( td_sub_lpt(k,nt+1) - td_sub_lpt(k,nt) )
    1794                 sums_ls_l(k,3) = td_sub_q(k,nt) + fac *                        &
    1795                                  ( td_sub_q(k,nt+1) - td_sub_q(k,nt) )
     1680                sums_ls_l(k,2) = td_sub_lpt(k,nt) + fac * ( td_sub_lpt(k,nt+1) - td_sub_lpt(k,nt) )
     1681                sums_ls_l(k,3) = td_sub_q(k,nt)   + fac * ( td_sub_q(k,nt+1)   - td_sub_q(k,nt) )
    17961682             ENDDO
    17971683
     
    18061692       !$OMP PARALLEL PRIVATE( i, j, k, tn )
    18071693       !$ tn = omp_get_thread_num()
    1808        IF ( radiation .AND. radiation_scheme == 'rrtmg' )  THEN
     1694       IF ( radiation  .AND. radiation_scheme == 'rrtmg' )  THEN
    18091695          !$OMP DO
    18101696          DO  i = nxl, nxr
     
    18131699                   flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
    18141700
    1815                    sums_l(k,100,tn)  = sums_l(k,100,tn)  + rad_lw_in(k,j,i)    &
    1816                                        * rmask(j,i,sr) * flag
    1817                    sums_l(k,101,tn)  = sums_l(k,101,tn)  + rad_lw_out(k,j,i)   &
    1818                                        * rmask(j,i,sr) * flag
    1819                    sums_l(k,102,tn)  = sums_l(k,102,tn)  + rad_sw_in(k,j,i)    &
    1820                                        * rmask(j,i,sr) * flag
    1821                    sums_l(k,103,tn)  = sums_l(k,103,tn)  + rad_sw_out(k,j,i)   &
    1822                                        * rmask(j,i,sr) * flag
    1823                    sums_l(k,104,tn)  = sums_l(k,104,tn)  + rad_lw_cs_hr(k,j,i) &
    1824                                        * rmask(j,i,sr) * flag
    1825                    sums_l(k,105,tn)  = sums_l(k,105,tn)  + rad_lw_hr(k,j,i)    &
    1826                                        * rmask(j,i,sr) * flag
    1827                    sums_l(k,106,tn)  = sums_l(k,106,tn)  + rad_sw_cs_hr(k,j,i) &
    1828                                        * rmask(j,i,sr) * flag
    1829                    sums_l(k,107,tn)  = sums_l(k,107,tn)  + rad_sw_hr(k,j,i)    &
    1830                                        * rmask(j,i,sr) * flag
     1701                   sums_l(k,100,tn)  = sums_l(k,100,tn) + rad_lw_in(k,j,i)    * rmask(j,i,sr) * flag
     1702                   sums_l(k,101,tn)  = sums_l(k,101,tn) + rad_lw_out(k,j,i)   * rmask(j,i,sr) * flag
     1703                   sums_l(k,102,tn)  = sums_l(k,102,tn) + rad_sw_in(k,j,i)    * rmask(j,i,sr) * flag
     1704                   sums_l(k,103,tn)  = sums_l(k,103,tn) + rad_sw_out(k,j,i)   * rmask(j,i,sr) * flag
     1705                   sums_l(k,104,tn)  = sums_l(k,104,tn) + rad_lw_cs_hr(k,j,i) * rmask(j,i,sr) * flag
     1706                   sums_l(k,105,tn)  = sums_l(k,105,tn) + rad_lw_hr(k,j,i)    * rmask(j,i,sr) * flag
     1707                   sums_l(k,106,tn)  = sums_l(k,106,tn) + rad_sw_cs_hr(k,j,i) * rmask(j,i,sr) * flag
     1708                   sums_l(k,107,tn)  = sums_l(k,107,tn) + rad_sw_hr(k,j,i)    * rmask(j,i,sr) * flag
    18311709                ENDDO
    18321710             ENDDO
     
    18481726                                      sums_l(:,45:pr_palm,i)
    18491727             IF ( max_pr_user > 0 )  THEN
    1850                 sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) = &
    1851                                    sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) + &
    1852                                    sums_l(:,pr_palm+1:pr_palm+max_pr_user,i)
     1728                sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) =                                        &
     1729                                                       sums_l(:,pr_palm+1:pr_palm+max_pr_user,0) + &
     1730                                                       sums_l(:,pr_palm+1:pr_palm+max_pr_user,i)
    18531731             ENDIF
    18541732
     
    18771755!--    Compute total sum from local sums
    18781756       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    1879        CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), ngp_sums, MPI_REAL,   &
    1880                            MPI_SUM, comm2d, ierr )
     1757       CALL MPI_ALLREDUCE( sums_l(nzb,1,0), sums(nzb,1), ngp_sums, MPI_REAL, MPI_SUM, comm2d, ierr )
    18811758       IF ( large_scale_forcing )  THEN
    1882           CALL MPI_ALLREDUCE( sums_ls_l(nzb,2), sums(nzb,83), ngp_sums_ls,     &
    1883                               MPI_REAL, MPI_SUM, comm2d, ierr )
     1759          CALL MPI_ALLREDUCE( sums_ls_l(nzb,2), sums(nzb,83), ngp_sums_ls, MPI_REAL, MPI_SUM,      &
     1760                              comm2d, ierr )
    18841761       ENDIF
    18851762
     
    18871764          IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    18881765          DO  i = 1, max_pr_cs
    1889              CALL MPI_ALLREDUCE( sums_l(nzb,pr_palm+max_pr_user+i,0),          &
    1890                                  sums(nzb,pr_palm+max_pr_user+i),              &
     1766             CALL MPI_ALLREDUCE( sums_l(nzb,pr_palm+max_pr_user+i,0),                              &
     1767                                 sums(nzb,pr_palm+max_pr_user+i),                                  &
    18911768                                 nzt+2-nzb, MPI_REAL, MPI_SUM, comm2d, ierr )
    18921769          ENDDO
     
    19101787
    19111788!
    1912 !--    Final values are obtained by division by the total number of grid points
    1913 !--    used for summation. After that store profiles.
    1914 !--    Check, if statistical regions do contain at least one grid point at the
    1915 !--    respective k-level, otherwise division by zero will lead to undefined
    1916 !--    values, which may cause e.g. problems with NetCDF output
     1789!--    Final values are obtained by division by the total number of grid points used for summation.
     1790!--    After that store profiles.
     1791!--    Check, if statistical regions do contain at least one grid point at the respective k-level,
     1792!--    otherwise division by zero will lead to undefined values, which may cause e.g. problems with
     1793!--    NetCDF output.
    19171794!--    Profiles:
    19181795       DO  k = nzb, nzt+1
     
    19421819
    19431820!--    u* and so on
    1944 !--    As sums(nzb:nzb+3,pr_palm) are full 2D arrays (us, usws, vsws, ts) whose
    1945 !--    size is always ( nx + 1 ) * ( ny + 1 ), defined at the first grid layer
    1946 !--    above the topography, they are being divided by ngp_2dh(sr)
    1947        sums(nzb:nzb+3,pr_palm)    = sums(nzb:nzb+3,pr_palm)    / &
    1948                                     ngp_2dh(sr)
    1949        sums(nzb+12,pr_palm)       = sums(nzb+12,pr_palm)       / &    ! qs
    1950                                     ngp_2dh(sr)
    1951        sums(nzb+13,pr_palm)       = sums(nzb+13,pr_palm)       / &    ! ss
    1952                                     ngp_2dh(sr)
    1953        sums(nzb+14,pr_palm)       = sums(nzb+14,pr_palm)       / &    ! surface temperature
    1954                                     ngp_2dh(sr)
     1821!--    As sums(nzb:nzb+3,pr_palm) are full 2D arrays (us, usws, vsws, ts) whose size is always
     1822!--    ( nx + 1 ) * ( ny + 1 ), defined at the first grid layer above the topography, they are
     1823!--    divided by ngp_2dh(sr)
     1824       sums(nzb:nzb+3,pr_palm)    = sums(nzb:nzb+3,pr_palm)  /  ngp_2dh(sr)
     1825       sums(nzb+12,pr_palm)       = sums(nzb+12,pr_palm)     /  ngp_2dh(sr)    ! qs
     1826       sums(nzb+13,pr_palm)       = sums(nzb+13,pr_palm)     /  ngp_2dh(sr)    ! ss
     1827       sums(nzb+14,pr_palm)       = sums(nzb+14,pr_palm)     /  ngp_2dh(sr)    ! surface temperature
     1828
    19551829!--    eges, e*
    1956        sums(nzb+4:nzb+5,pr_palm)  = sums(nzb+4:nzb+5,pr_palm)  / &
    1957                                     ngp_3d(sr)
     1830       sums(nzb+4:nzb+5,pr_palm)  = sums(nzb+4:nzb+5,pr_palm)  /  ngp_3d(sr)
    19581831!--    Old and new divergence
    1959        sums(nzb+9:nzb+10,pr_palm) = sums(nzb+9:nzb+10,pr_palm) / &
    1960                                     ngp_3d_inner(sr)
     1832       sums(nzb+9:nzb+10,pr_palm) = sums(nzb+9:nzb+10,pr_palm) /  ngp_3d_inner(sr)
    19611833
    19621834!--    User-defined profiles
    19631835       IF ( max_pr_user > 0 )  THEN
    19641836          DO  k = nzb, nzt+1
    1965              sums(k,pr_palm+1:pr_palm+max_pr_user) = &
    1966                                     sums(k,pr_palm+1:pr_palm+max_pr_user) / &
    1967                                     ngp_2dh_s_inner(k,sr)
    1968           ENDDO
    1969        ENDIF
    1970 
    1971        IF ( air_chemistry ) THEN
     1837             sums(k,pr_palm+1:pr_palm+max_pr_user) =  sums(k,pr_palm+1:pr_palm+max_pr_user) /      &
     1838                                                      ngp_2dh_s_inner(k,sr)
     1839          ENDDO
     1840       ENDIF
     1841
     1842       IF ( air_chemistry )  THEN
    19721843          IF ( max_pr_cs > 0 )  THEN
    19731844             DO k = nzb, nzt+1
    1974                 sums(k, pr_palm+1:pr_palm+max_pr_user+max_pr_cs) = &
    1975                                  sums(k, pr_palm+1:pr_palm+max_pr_user+max_pr_cs) / &
    1976                                  ngp_2dh_s_inner(k,sr)
     1845                sums(k, pr_palm+1:pr_palm+max_pr_user+max_pr_cs) =                                 &
     1846                                                sums(k, pr_palm+1:pr_palm+max_pr_user+max_pr_cs) / &
     1847                                                ngp_2dh_s_inner(k,sr)
    19771848             ENDDO
    19781849          ENDIF
    19791850       ENDIF
    19801851
    1981        IF ( salsa ) THEN
     1852       IF ( salsa )  THEN
    19821853          IF ( max_pr_salsa > 0 )  THEN
    19831854             DO k = nzb, nzt+1
     
    20201891       hom(:,1,37,sr) = sums(:,37)     ! w*e*
    20211892       hom(:,1,38,sr) = sums(:,38)     ! w*3
    2022        hom(:,1,39,sr) = sums(:,38) / ( abs( sums(:,32) ) + 1E-20_wp )**1.5_wp   ! Sw
     1893       hom(:,1,39,sr) = sums(:,38) / ( ABS( sums(:,32) ) + 1E-20_wp )**1.5_wp   ! Sw
    20231894       hom(:,1,40,sr) = sums(:,40)     ! p
    20241895       hom(:,1,45,sr) = sums(:,45)     ! w"vpt"
     
    21241995       hom(:,1,120,sr) = rho_air_zw    ! rho_air_zw in Kg/m^3
    21251996
    2126        IF ( kolmogorov_length_scale ) THEN
     1997       IF ( kolmogorov_length_scale )  THEN
    21271998          hom(:,1,121,sr) = sums(:,121) * 1E3_wp  ! eta in mm
    21281999       ENDIF
     
    21522023!
    21532024!--    Determine the boundary layer height using two different schemes.
    2154 !--    First scheme: Starting from the Earth's (Ocean's) surface, look for the
    2155 !--    first relative minimum (maximum) of the total heat flux.
    2156 !--    The corresponding height is assumed as the boundary layer height, if it
    2157 !--    is less than 1.5 times the height where the heat flux becomes negative
    2158 !--    (positive) for the first time. Attention: the resolved vertical sensible
    2159 !--    heat flux (hom(:,1,17,sr) = w*pt*) is not known at the beginning because
    2160 !--    the calculation happens in advec_s_ws which is called after
    2161 !--    flow_statistics. Therefore z_i is directly taken from restart data at
    2162 !--    the beginning of restart runs.
    2163        IF ( TRIM( initializing_actions ) /= 'read_restart_data' .OR.           &
     2025!--    First scheme: Starting from the Earth's (Ocean's) surface, look for the first relative
     2026!--    minimum (maximum) of the total heat flux.
     2027!--    The corresponding height is assumed as the boundary layer height, if it is less than 1.5
     2028!--    times the height where the heat flux becomes negative (positive) for the first time.
     2029!--    Attention: the resolved vertical sensible heat flux (hom(:,1,17,sr) = w*pt*) is not known at
     2030!--    the beginning because the calculation happens in advec_s_ws which is called after
     2031!--    flow_statistics. Therefore z_i is directly taken from restart data at the beginning of
     2032!--    restart runs.
     2033       IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .OR.                              &
    21642034            simulated_time_at_begin /= simulated_time ) THEN
    21652035
     
    21732043                   height = zw(k)
    21742044                ENDIF
    2175                 IF ( hom(k,1,18,sr) < -1.0E-8_wp  .AND.                        &
    2176                      hom(k-1,1,18,sr) > hom(k,1,18,sr) )  THEN
     2045                IF ( hom(k,1,18,sr) < -1.0E-8_wp  .AND.  hom(k-1,1,18,sr) > hom(k,1,18,sr) )  THEN
    21772046                   IF ( zw(k) < 1.5_wp * height )  THEN
    21782047                      z_i(1) = zw(k)
     
    21892058                   height = zw(k)
    21902059                ENDIF
    2191                 IF ( hom(k,1,18,sr) < -1.0E-8_wp  .AND.                        &
    2192                      hom(k+1,1,18,sr) > hom(k,1,18,sr) )  THEN
     2060                IF ( hom(k,1,18,sr) < -1.0E-8_wp  .AND.  hom(k+1,1,18,sr) > hom(k,1,18,sr) )  THEN
    21932061                   IF ( zw(k) < 1.5_wp * height )  THEN
    21942062                      z_i(1) = zw(k)
     
    22022070
    22032071!
    2204 !--       Second scheme: Gradient scheme from Sullivan et al. (1998), modified
    2205 !--       by Uhlenbrock(2006). The boundary layer height is the height with the
    2206 !--       maximal local temperature gradient: starting from the second (the
    2207 !--       last but one) vertical gridpoint, the local gradient must be at least
    2208 !--       0.2K/100m and greater than the next four gradients.
     2072!--       Second scheme: Gradient scheme from Sullivan et al. (1998), modified by Uhlenbrock(2006).
     2073!--       The boundary layer height is the height with the maximal local temperature gradient:
     2074!--       starting from the second (the last but one) vertical gridpoint, the local gradient must be
     2075!--       at least 0.2K/100m and greater than the next four gradients.
    22092076!--       WARNING: The threshold value of 0.2K/100m must be adjusted for the
    22102077!--       ocean case!
     
    22172084          IF ( ocean_mode )  THEN
    22182085             DO  k = nzt+1, nzb+5, -1
    2219                 IF ( dptdz(k) > dptdz_threshold  .AND.                         &
    2220                      dptdz(k) > dptdz(k-1)  .AND.  dptdz(k) > dptdz(k-2)  .AND.&
     2086                IF ( dptdz(k) > dptdz_threshold  .AND.                                             &
     2087                     dptdz(k) > dptdz(k-1)  .AND.  dptdz(k) > dptdz(k-2)  .AND.                    &
    22212088                     dptdz(k) > dptdz(k-3)  .AND.  dptdz(k) > dptdz(k-4) )  THEN
    22222089                   z_i(2) = zw(k-1)
     
    22262093          ELSE
    22272094             DO  k = nzb+1, nzt-3
    2228                 IF ( dptdz(k) > dptdz_threshold  .AND.                         &
    2229                      dptdz(k) > dptdz(k+1)  .AND.  dptdz(k) > dptdz(k+2)  .AND.&
     2095                IF ( dptdz(k) > dptdz_threshold  .AND.                                             &
     2096                     dptdz(k) > dptdz(k+1)  .AND.  dptdz(k) > dptdz(k+2)  .AND.                    &
    22302097                     dptdz(k) > dptdz(k+3)  .AND.  dptdz(k) > dptdz(k+4) )  THEN
    22312098                   z_i(2) = zw(k-1)
     
    22412108
    22422109!
    2243 !--    Determine vertical index which is nearest to the mean surface level
    2244 !--    height of the respective statistic region
     2110!--    Determine vertical index which is nearest to the mean surface level height of the respective
     2111!--    statistic region
    22452112       DO  k = nzb, nzt
    22462113          IF ( zw(k) >= mean_surface_level_height(sr) )  THEN
     
    22512118
    22522119!
    2253 !--    Computation of both the characteristic vertical velocity and
    2254 !--    the characteristic convective boundary layer temperature.
    2255 !--    The inversion height entering into the equation is defined with respect
    2256 !--    to the mean surface level height of the respective statistic region.
    2257 !--    The horizontal average at surface level index + 1 is input for the
    2258 !--    average temperature.
    2259        IF ( hom(k_surface_level,1,18,sr) > 1.0E-8_wp  .AND.  z_i(1) /= 0.0_wp )&
    2260        THEN
    2261           hom(nzb+8,1,pr_palm,sr) =                                            &
    2262              ( g / hom(k_surface_level+1,1,4,sr) *                             &
    2263              ( hom(k_surface_level,1,18,sr) /                                  &
    2264              ( heatflux_output_conversion(nzb) * rho_air(nzb) ) )              &
    2265              * ABS( z_i(1) - mean_surface_level_height(sr) ) )**0.333333333_wp
     2120!--    Computation of both the characteristic vertical velocity and the characteristic convective
     2121!--    boundary layer temperature.
     2122!--    The inversion height entering into the equation is defined with respect to the mean surface
     2123!--    level height of the respective statistic region.
     2124!--    The horizontal average at surface level index + 1 is input for the average temperature.
     2125       IF ( hom(k_surface_level,1,18,sr) > 1.0E-8_wp  .AND.  z_i(1) /= 0.0_wp )  THEN
     2126          hom(nzb+8,1,pr_palm,sr) =                                                                &
     2127                                   ( g / hom(k_surface_level+1,1,4,sr) *                           &
     2128                                   ( hom(k_surface_level,1,18,sr) /                                &
     2129                                   ( heatflux_output_conversion(nzb) * rho_air(nzb) ) )            &
     2130                                   * ABS( z_i(1) - mean_surface_level_height(sr) ) )**0.333333333_wp
    22662131       ELSE
    22672132          hom(nzb+8,1,pr_palm,sr)  = 0.0_wp
     
    22692134
    22702135!
    2271 !--    Collect the time series quantities. Please note, timeseries quantities
    2272 !--    which are collected from horizontally averaged profiles, e.g. wpt
    2273 !--    or pt(zp), are treated specially. In case of elevated model surfaces,
    2274 !--    index nzb+1 might be within topography and data will be zero. Therefore,
    2275 !--    take value for the first atmosphere index, which is topo_min_level+1.
     2136!--    Collect the time series quantities. Please note, timeseries quantities which are collected
     2137!--    from horizontally averaged profiles, e.g. wpt or pt(zp), are treated specially. In case of
     2138!--    elevated model surfaces, index nzb+1 might be within topography and data will be zero.
     2139!--    Therefore, take value for the first atmosphere index, which is topo_min_level+1.
    22762140       ts_value(1,sr) = hom(nzb+4,1,pr_palm,sr)        ! E
    22772141       ts_value(2,sr) = hom(nzb+5,1,pr_palm,sr)        ! E*
  • palm/trunk/SOURCE/global_min_max.f90

    r4429 r4646  
    1 !> @file global_min_max.f90
    2 !------------------------------------------------------------------------------!
     1!--------------------------------------------------------------------------------------------------!
    32! This file is part of the PALM model system.
    43!
    5 ! PALM is free software: you can redistribute it and/or modify it under the
    6 ! terms of the GNU General Public License as published by the Free Software
    7 ! Foundation, either version 3 of the License, or (at your option) any later
    8 ! version.
    9 !
    10 ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
    11 ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
    12 ! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
    13 !
    14 ! You should have received a copy of the GNU General Public License along with
    15 ! PALM. If not, see <http://www.gnu.org/licenses/>.
     4! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
     5! Public License as published by the Free Software Foundation, either version 3 of the License, or
     6! (at your option) any later version.
     7!
     8! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
     9! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
     10! Public License for more details.
     11!
     12! You should have received a copy of the GNU General Public License along with PALM. If not, see
     13! <http://www.gnu.org/licenses/>.
    1614!
    1715! Copyright 1997-2020 Leibniz Universitaet Hannover
    18 !------------------------------------------------------------------------------!
     16!--------------------------------------------------------------------------------------------------!
    1917!
    2018! Current revisions:
    2119! ------------------
    22 ! 
    23 ! 
     20!
     21!
    2422! Former revisions:
    2523! -----------------
    2624! $Id$
     25! file re-formatted to follow the PALM coding standard
     26!
     27! 4429 2020-02-27 15:24:30Z raasch
    2728! bugfix: cpp-directives added for serial mode
    28 ! 
     29!
    2930! 4360 2020-01-07 11:25:50Z suehring
    3031! OpenACC support added
    31 ! 
     32!
    3233! 4182 2019-08-22 15:20:23Z scharf
    3334! Corrected "Former revisions" section
    34 ! 
     35!
    3536! 3655 2019-01-07 16:51:22Z knoop
    3637! Corrected "Former revisions" section
     
    4344! ------------
    4445!> Determine the array minimum/maximum and the corresponding indices.
    45 !------------------------------------------------------------------------------!
    46  SUBROUTINE global_min_max( i1, i2, j1, j2, k1, k2, ar, mode, offset, value, &
    47                             value_ijk, value1, value1_ijk )
    48  
    49 
    50     USE indices,                                                               &
     46!--------------------------------------------------------------------------------------------------!
     47 SUBROUTINE global_min_max( i1, i2, j1, j2, k1, k2, ar, mode, offset, value, value_ijk, value1,    &
     48                            value1_ijk )
     49
     50
     51    USE indices,                                                                                   &
    5152        ONLY:  nbgp, ny, nx
    52        
     53
    5354    USE kinds
    54    
     55
    5556    USE pegrid
    5657
     
    7273    INTEGER(iwp) ::  k1             !<
    7374    INTEGER(iwp) ::  k2             !<
    74     INTEGER(iwp) ::  fmax_ijk(3)    !<
    75     INTEGER(iwp) ::  fmax_ijk_l(3)  !<
    76     INTEGER(iwp) ::  fmin_ijk(3)    !<
    77     INTEGER(iwp) ::  fmin_ijk_l(3)  !<
    7875    INTEGER(iwp) ::  value_ijk(3)   !<
    79    
    80     INTEGER(iwp), OPTIONAL ::  value1_ijk(3)  !<
    81    
    82     REAL(wp) ::  offset                 !<
    83     REAL(wp) ::  value                  !<
     76
     77    INTEGER(iwp), DIMENSION(3) ::  fmax_ijk    !<
     78    INTEGER(iwp), DIMENSION(3) ::  fmax_ijk_l  !<
     79    INTEGER(iwp), DIMENSION(3) ::  fmin_ijk    !<
     80    INTEGER(iwp), DIMENSION(3) ::  fmin_ijk_l  !<
     81
     82    INTEGER(iwp), DIMENSION(3), OPTIONAL ::  value1_ijk  !<
     83
     84    REAL(wp) ::  offset            !<
     85    REAL(wp) ::  value             !<
     86    REAL(wp), OPTIONAL ::  value1  !<
     87
    8488    REAL(wp) ::  ar(i1:i2,j1:j2,k1:k2)  !<
    85    
     89
    8690#if defined( __ibm )
    8791    REAL(sp) ::  fmax(2)    !<
     
    8993    REAL(sp) ::  fmin(2)    !<
    9094    REAL(sp) ::  fmin_l(2)  !<
    91              ! on 32bit-machines MPI_2REAL must not be replaced 
     95             ! on 32bit-machines MPI_2REAL must not be replaced
    9296             ! by MPI_2DOUBLE_PRECISION
    9397#else
    94     REAL(wp) ::  fmax(2)    !<
    95     REAL(wp) ::  fmax_l(2)  !<
    96     REAL(wp) ::  fmin(2)    !<
    97     REAL(wp) ::  fmin_l(2)  !<
    98 #endif
     98    REAL(wp), DIMENSION(2) ::  fmax    !<
     99    REAL(wp), DIMENSION(2) ::  fmax_l  !<
     100    REAL(wp), DIMENSION(2) ::  fmin    !<
     101    REAL(wp), DIMENSION(2) ::  fmin_l  !<
     102#endif
     103
    99104#if defined( _OPENACC )
     105    INTEGER(iwp) ::  count_eq   !< counter for locations of maximum
    100106    REAL(wp)     ::  red        !< scalar for reduction with OpenACC
    101     INTEGER(iwp) ::  count_eq   !< counter for locations of maximum
    102 #endif
    103     REAL(wp), OPTIONAL ::  value1  !<
     107#endif
    104108
    105109
     
    119123       fmin_l(2)  = myid
    120124       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    121        CALL MPI_ALLREDUCE( fmin_l, fmin, 1, MPI_2REAL, MPI_MINLOC, comm2d, &
    122                            ierr )
     125       CALL MPI_ALLREDUCE( fmin_l, fmin, 1, MPI_2REAL, MPI_MINLOC, comm2d, ierr )
    123126
    124127!
     
    127130       IF ( id_fmin /= 0 )  THEN
    128131          IF ( myid == 0 )  THEN
    129              CALL MPI_RECV( fmin_ijk, 3, MPI_INTEGER, id_fmin, 0, comm2d, &
    130                             status, ierr )
     132             CALL MPI_RECV( fmin_ijk, 3, MPI_INTEGER, id_fmin, 0, comm2d, status, ierr )
    131133          ELSEIF ( myid == id_fmin )  THEN
    132134             CALL MPI_SEND( fmin_ijk_l, 3, MPI_INTEGER, 0, 0, comm2d, ierr )
     
    160162       fmax_l(2)  = myid
    161163       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    162        CALL MPI_ALLREDUCE( fmax_l, fmax, 1, MPI_2REAL, MPI_MAXLOC, comm2d, &
    163                            ierr )
     164       CALL MPI_ALLREDUCE( fmax_l, fmax, 1, MPI_2REAL, MPI_MAXLOC, comm2d, ierr )
    164165
    165166!
     
    168169       IF ( id_fmax /= 0 )  THEN
    169170          IF ( myid == 0 )  THEN
    170              CALL MPI_RECV( fmax_ijk, 3, MPI_INTEGER, id_fmax, 0, comm2d, &
    171                             status, ierr )
     171             CALL MPI_RECV( fmax_ijk, 3, MPI_INTEGER, id_fmax, 0, comm2d, status, ierr )
    172172          ELSEIF ( myid == id_fmax )  THEN
    173173             CALL MPI_SEND( fmax_ijk_l, 3, MPI_INTEGER, 0, 0, comm2d, ierr )
     
    177177       ENDIF
    178178!
    179 !--    send the indices of the just determined array maximum to other PEs
     179!--    Send the indices of the just determined array maximum to other PEs
    180180       CALL MPI_BCAST( fmax_ijk, 3, MPI_INTEGER, 0, comm2d, ierr )
    181181#else
     
    226226       IF ( count_eq == 1 ) THEN
    227227!
    228 !--       We found a single maximum element and correctly got its position. Transfer its
    229 !--       value to handle the negative case correctly.
     228!--       We found a single maximum element and correctly got its position. Transfer its value to
     229!--       handle the negative case correctly.
    230230          !$ACC UPDATE HOST(ar(fmax_ijk_l(1):fmax_ijk_l(1),fmax_ijk_l(2),fmax_ijk_l(3)))
    231231       ELSE
     
    257257#if defined( _OPENACC )
    258258!
    259 !--       Close ELSE case from above
     259!--    Close ELSE case from above
    260260       ENDIF
    261261#endif
     
    263263!
    264264!--    Set a flag in case that the determined value is negative.
    265 !--    A constant offset has to be subtracted in order to handle the special
    266 !--    case i=0 correctly
     265!--    A constant offset has to be subtracted in order to handle the special case i=0 correctly.
    267266       IF ( ar(fmax_ijk_l(1),fmax_ijk_l(2),fmax_ijk_l(3)) < 0.0_wp )  THEN
    268267          fmax_ijk_l(1) = -fmax_ijk_l(1) - 10
     
    272271       fmax_l(2)  = myid
    273272       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    274        CALL MPI_ALLREDUCE( fmax_l, fmax, 1, MPI_2REAL, MPI_MAXLOC, comm2d, &
    275                            ierr )
     273       CALL MPI_ALLREDUCE( fmax_l, fmax, 1, MPI_2REAL, MPI_MAXLOC, comm2d, ierr )
    276274
    277275!
     
    280278       IF ( id_fmax /= 0 )  THEN
    281279          IF ( myid == 0 )  THEN
    282              CALL MPI_RECV( fmax_ijk, 3, MPI_INTEGER, id_fmax, 0, comm2d, &
    283                             status, ierr )
     280             CALL MPI_RECV( fmax_ijk, 3, MPI_INTEGER, id_fmax, 0, comm2d, status, ierr )
    284281          ELSEIF ( myid == id_fmax )  THEN
    285282             CALL MPI_SEND( fmax_ijk_l, 3, MPI_INTEGER, 0, 0, comm2d, ierr )
     
    311308          DO  j = j1, j2
    312309!
    313 !--          Attention: the lowest gridpoint is excluded here, because there
    314 !--          ---------  is no advection at nzb=0 and mode 'absoff' is only
    315 !--                     used for calculating u,v extrema for CFL-criteria
     310!--          Attention: the lowest gridpoint is excluded here, because there is no advection at
     311!--          ---------- nzb=0 and mode 'absoff' is only used for calculating u,v extrema for
     312!--                     CFL-criteria.
    316313             DO  i = i1+1, i2
    317314                IF ( ABS( ar(i,j,k) - offset ) > fmax_l(1) )  THEN
     
    327324!
    328325!--    Set a flag in case that the determined value is negative.
    329 !--    A constant offset has to be subtracted in order to handle the special
    330 !--    case i=0 correctly
     326!--    A constant offset has to be subtracted in order to handle the special case i=0 correctly.
    331327       IF ( ar(fmax_ijk_l(1),fmax_ijk_l(2),fmax_ijk_l(3)) < 0.0_wp )  THEN
    332328          fmax_ijk_l(1) = -fmax_ijk_l(1) - 10
     
    336332       fmax_l(2)  = myid
    337333       IF ( collective_wait )  CALL MPI_BARRIER( comm2d, ierr )
    338        CALL MPI_ALLREDUCE( fmax_l, fmax, 1, MPI_2REAL, MPI_MAXLOC, comm2d, &
    339                            ierr )
     334       CALL MPI_ALLREDUCE( fmax_l, fmax, 1, MPI_2REAL, MPI_MAXLOC, comm2d, ierr )
    340335
    341336!
     
    344339       IF ( id_fmax /= 0 )  THEN
    345340          IF ( myid == 0 )  THEN
    346              CALL MPI_RECV( fmax_ijk, 3, MPI_INTEGER, id_fmax, 0, comm2d, &
    347                             status, ierr )
     341             CALL MPI_RECV( fmax_ijk, 3, MPI_INTEGER, id_fmax, 0, comm2d, status, ierr )
    348342          ELSEIF ( myid == id_fmax )  THEN
    349343             CALL MPI_SEND( fmax_ijk_l, 3, MPI_INTEGER, 0, 0, comm2d, ierr )
  • palm/trunk/SOURCE/gust_mod.f90

    r4535 r4646  
    11!> @file gust_mod.f90
    2 !------------------------------------------------------------------------------!
    3 ! This file is part of PALM.
    4 !
    5 ! PALM is free software: you can redistribute it and/or modify it under the
    6 ! terms of the GNU General Public License as published by the Free Software
    7 ! Foundation, either version 3 of the License, or (at your option) any later
    8 ! version.
    9 !
    10 ! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
    11 ! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
    12 ! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
    13 !
    14 ! You should have received a copy of the GNU General Public License along with
    15 ! PALM. If not, see <http://www.gnu.org/licenses/>.
     2!--------------------------------------------------------------------------------------------------!
     3! This file is part of the PALM model system.
     4!
     5! PALM is free software: you can redistribute it and/or modify it under the terms of the GNU General
     6! Public License as published by the Free Software Foundation, either version 3 of the License, or
     7! (at your option) any later version.
     8!
     9! PALM is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the
     10! implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General
     11! Public License for more details.
     12!
     13! You should have received a copy of the GNU General Public License along with PALM. If not, see
     14! <http://www.gnu.org/licenses/>.
    1615!
    1716! Copyright 1997-2020 Leibniz Universitaet Hannover
    18 !------------------------------------------------------------------------------!
     17!--------------------------------------------------------------------------------------------------!
    1918!
    2019! Current revisions:
    2120! -----------------
    22 ! 
    23 ! 
     21!
     22!
    2423! Former revisions:
    2524! -----------------
    2625! $Id$
     26! file re-formatted to follow the PALM coding standard
     27!
     28! 4535 2020-05-15 12:07:23Z raasch
    2729! bugfix for restart data format query
    28 ! 
     30!
    2931! 4517 2020-05-03 14:29:30Z raasch
    3032! added restart with MPI-IO for reading local arrays
    31 ! 
     33!
    3234! 4495 2020-04-13 20:11:20Z raasch
    3335! restart data handling with MPI-IO added
    34 ! 
     36!
    3537! 4360 2020-01-07 11:25:50Z suehring
    36 ! CASE statement for dummy variable u2_av in gust_rrd_local changed to avoid
    37 ! unintended interdependencies with user-defined variables
    38 ! 
     38! CASE statement for dummy variable u2_av in gust_rrd_local changed to avoid unintended
     39! interdependencies with user-defined variables
     40!
    3941! 3837 2019-03-28 16:55:58Z knoop
    4042! unused variable for file index removed from rrd-subroutines parameter list
    41 ! 
     43!
    4244! 3725 2019-02-07 10:11:02Z raasch
    4345! dummy statement modified to avoid compiler warnings about unused variables
    44 ! 
     46!
    4547! 3685 2019-01-21 01:02:11Z knoop
    4648! Some interface calls moved to module_interface + cleanup
    47 ! 
     49!
    4850! 3665 2019-01-10 08:28:24Z raasch
    4951! dummy statements added to avoid compiler warnings about unused variables
    50 ! 
     52!
    5153! 3655 2019-01-07 16:51:22Z knoop
    5254! Bugfix: domain bounds of local_pf corrected
    53 ! 
    54 ! 
     55!
     56!
    5557! Interfaces concerning data output updated
    56 ! 
    57 ! 
     58!
     59!
    5860! renamed gust_par to gust_parameters
    59 ! 
    60 ! 
     61!
     62!
    6163! Initial interface definition
    6264!
    63 ! 
     65!
    6466! Description:
    6567! ------------
     
    6769!>
    6870!> @todo This is just a dummy module. The actual module ist not released yet.
    69 !------------------------------------------------------------------------------!
     71!--------------------------------------------------------------------------------------------------!
    7072 MODULE gust_mod
    7173
     
    7375        ONLY:  restart_data_format_output
    7476
    75     USE indices,                                                               &
     77    USE indices,                                                                                   &
    7678        ONLY:  nxl, nxlg, nxr, nxrg, nys, nysg, nyn, nyng, nzb, nzt
    7779
     
    9496!
    9597!-- Public functions
    96     PUBLIC &
    97        gust_parin, &
    98        gust_check_parameters, &
    99        gust_check_data_output_pr, &
    100        gust_check_data_output, &
    101        gust_init_arrays, &
    102        gust_init, &
    103        gust_define_netcdf_grid, &
    104        gust_header, &
    105        gust_actions, &
    106        gust_prognostic_equations, &
    107        gust_swap_timelevel, &
    108        gust_3d_data_averaging, &
    109        gust_data_output_2d, &
    110        gust_data_output_3d, &
    111        gust_statistics, &
    112        gust_rrd_global, &
    113        gust_wrd_global, &
    114        gust_rrd_local, &
     98    PUBLIC                                                                                         &
     99       gust_parin,                                                                                 &
     100       gust_check_parameters,                                                                      &
     101       gust_check_data_output_pr,                                                                  &
     102       gust_check_data_output,                                                                     &
     103       gust_init_arrays,                                                                           &
     104       gust_init,                                                                                  &
     105       gust_define_netcdf_grid,                                                                    &
     106       gust_header,                                                                                &
     107       gust_actions,                                                                               &
     108       gust_prognostic_equations,                                                                  &
     109       gust_swap_timelevel,                                                                        &
     110       gust_3d_data_averaging,                                                                     &
     111       gust_data_output_2d,                                                                        &
     112       gust_data_output_3d,                                                                        &
     113       gust_statistics,                                                                            &
     114       gust_rrd_global,                                                                            &
     115       gust_wrd_global,                                                                            &
     116       gust_rrd_local,                                                                             &
    115117       gust_wrd_local
    116118!
    117119!-- Public parameters, constants and initial values
    118     PUBLIC &
     120    PUBLIC                                                                                         &
    119121       gust_module_enabled
    120122
     
    203205
    204206
    205 !------------------------------------------------------------------------------!
     207!--------------------------------------------------------------------------------------------------!
    206208! Description:
    207209! ------------
    208210!> Parin for &gust_parameters for gust module
    209 !------------------------------------------------------------------------------!
     211!--------------------------------------------------------------------------------------------------!
    210212    SUBROUTINE gust_parin
    211213
     
    215217       CHARACTER (LEN=80)  ::  line  !< dummy string that contains the current line of the parameter file
    216218
    217        NAMELIST /gust_parameters/  &
     219       NAMELIST /gust_parameters/                                                                  &
    218220          gust_module_enabled
    219221
     
    240242
    241243
    242 !------------------------------------------------------------------------------!
     244!--------------------------------------------------------------------------------------------------!
    243245! Description:
    244246! ------------
    245247!> Check parameters routine for gust module
    246 !------------------------------------------------------------------------------!
     248!--------------------------------------------------------------------------------------------------!
    247249    SUBROUTINE gust_check_parameters
    248250
     
    254256
    255257
    256 !------------------------------------------------------------------------------!
     258!--------------------------------------------------------------------------------------------------!
    257259! Description:
    258260! ------------
    259261!> Check data output of profiles for gust module
    260 !------------------------------------------------------------------------------!
     262!--------------------------------------------------------------------------------------------------!
    261263    SUBROUTINE gust_check_data_output_pr( variable, var_count, unit, dopr_unit )
    262264
     
    264266       IMPLICIT NONE
    265267
     268       CHARACTER (LEN=*) ::  dopr_unit !< local value of dopr_unit
    266269       CHARACTER (LEN=*) ::  unit      !<
    267270       CHARACTER (LEN=*) ::  variable  !<
    268        CHARACTER (LEN=*) ::  dopr_unit !< local value of dopr_unit
    269271
    270272       INTEGER(iwp) ::  var_count      !<
     
    276278    END SUBROUTINE gust_check_data_output_pr
    277279
    278 !------------------------------------------------------------------------------!
     280!--------------------------------------------------------------------------------------------------!
    279281! Description:
    280282! ------------
    281283!> Check data output for gust module
    282 !------------------------------------------------------------------------------!
     284!--------------------------------------------------------------------------------------------------!
    283285    SUBROUTINE gust_check_data_output( var, unit )
    284286
     
    296298
    297299
    298 !------------------------------------------------------------------------------!
     300!--------------------------------------------------------------------------------------------------!
    299301! Description:
    300302! ------------
    301303!> Allocate gust module arrays and define pointers
    302 !------------------------------------------------------------------------------!
     304!--------------------------------------------------------------------------------------------------!
    303305    SUBROUTINE gust_init_arrays
    304306
     
    310312
    311313
    312 !------------------------------------------------------------------------------!
     314!--------------------------------------------------------------------------------------------------!
    313315! Description:
    314316! ------------
    315317!> Initialization of the gust module
    316 !------------------------------------------------------------------------------!
     318!--------------------------------------------------------------------------------------------------!
    317319    SUBROUTINE gust_init
    318320
     
    324326
    325327
    326 !------------------------------------------------------------------------------!
     328!--------------------------------------------------------------------------------------------------!
    327329!
    328330! Description:
     
    330332!> Subroutine defining appropriate grid for netcdf variables.
    331333!> It is called out from subroutine netcdf.
    332 !------------------------------------------------------------------------------!
     334!--------------------------------------------------------------------------------------------------!
    333335    SUBROUTINE gust_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
    334336
     
    336338       IMPLICIT NONE
    337339
    338        CHARACTER (LEN=*), INTENT(IN)  ::  var         !<
    339        LOGICAL, INTENT(IN)           ::  found       !<
    340340       CHARACTER (LEN=*), INTENT(IN) ::  grid_x      !<
    341341       CHARACTER (LEN=*), INTENT(IN) ::  grid_y      !<
    342342       CHARACTER (LEN=*), INTENT(IN) ::  grid_z      !<
     343       CHARACTER (LEN=*), INTENT(IN) ::  var         !<
     344
     345       LOGICAL, INTENT(IN)           ::  found       !<
    343346
    344347!
     
    349352
    350353
    351 !------------------------------------------------------------------------------!
     354!--------------------------------------------------------------------------------------------------!
    352355! Description:
    353356! ------------
    354357!> Header output for gust module
    355 !------------------------------------------------------------------------------!
     358!--------------------------------------------------------------------------------------------------!
    356359    SUBROUTINE gust_header ( io )
    357360
     
    368371
    369372
    370 !------------------------------------------------------------------------------!
     373!--------------------------------------------------------------------------------------------------!
    371374! Description:
    372375! ------------
    373376!> Call for all grid points
    374 !------------------------------------------------------------------------------!
     377!--------------------------------------------------------------------------------------------------!
    375378    SUBROUTINE gust_actions( location )
    376379
     
    387390
    388391
    389 !------------------------------------------------------------------------------!
     392!--------------------------------------------------------------------------------------------------!
    390393! Description:
    391394! ------------
    392395!> Call for grid point i,j
    393 !------------------------------------------------------------------------------!
     396!--------------------------------------------------------------------------------------------------!
    394397    SUBROUTINE gust_actions_ij( i, j, location )
    395398
     
    409412
    410413
    411 !------------------------------------------------------------------------------!
     414!--------------------------------------------------------------------------------------------------!
    412415! Description:
    413416! ------------
    414417!> Call for all grid points
    415 !------------------------------------------------------------------------------!
     418!--------------------------------------------------------------------------------------------------!
    416419    SUBROUTINE gust_prognostic_equations()
    417420
     
    423426
    424427
    425 !------------------------------------------------------------------------------!
     428!--------------------------------------------------------------------------------------------------!
    426429! Description:
    427430! ------------
    428431!> Call for grid point i,j
    429 !------------------------------------------------------------------------------!
     432!--------------------------------------------------------------------------------------------------!
    430433    SUBROUTINE gust_prognostic_equations_ij( i, j, i_omp_start, tn )
    431434
    432435
    433436       INTEGER(iwp), INTENT(IN) ::  i            !< grid index in x-direction
     437       INTEGER(iwp), INTENT(IN) ::  i_omp_start  !< first loop index of i-loop in prognostic_equations
    434438       INTEGER(iwp), INTENT(IN) ::  j            !< grid index in y-direction
    435        INTEGER(iwp), INTENT(IN) ::  i_omp_start  !< first loop index of i-loop in prognostic_equations
    436439       INTEGER(iwp), INTENT(IN) ::  tn           !< task number of openmp task
    437440
     
    443446
    444447
    445 !------------------------------------------------------------------------------!
     448!--------------------------------------------------------------------------------------------------!
    446449! Description:
    447450! ------------
    448451!> Swapping of timelevels
    449 !------------------------------------------------------------------------------!
     452!--------------------------------------------------------------------------------------------------!
    450453    SUBROUTINE gust_swap_timelevel ( mod_count )
    451454
     
    462465
    463466
    464 !------------------------------------------------------------------------------!
     467!--------------------------------------------------------------------------------------------------!
    465468!
    466469! Description:
    467470! ------------
    468471!> Subroutine for averaging 3D data
    469 !------------------------------------------------------------------------------!
     472!--------------------------------------------------------------------------------------------------!
    470473    SUBROUTINE gust_3d_data_averaging( mode, variable )
    471474
     
    474477
    475478       CHARACTER (LEN=*) ::  mode    !<
    476        CHARACTER (LEN=*) :: variable !<
     479       CHARACTER (LEN=*) ::  variable !<
    477480
    478481!
     
    482485    END SUBROUTINE gust_3d_data_averaging
    483486
    484 !------------------------------------------------------------------------------!
     487!--------------------------------------------------------------------------------------------------!
    485488!
    486489! Description:
    487490! ------------
    488491!> Subroutine defining 2D output variables
    489 !------------------------------------------------------------------------------!
    490     SUBROUTINE gust_data_output_2d( av, variable, found, grid, mode, local_pf, &
    491                                     two_d, nzb_do, nzt_do, fill_value )
     492!--------------------------------------------------------------------------------------------------!
     493    SUBROUTINE gust_data_output_2d( av, variable, found, grid, mode, local_pf, two_d, nzb_do,      &
     494                                    nzt_do, fill_value )
    492495
    493496
     
    495498
    496499       CHARACTER (LEN=*), INTENT(INOUT) ::  grid       !< name of vertical grid
    497        CHARACTER (LEN=*), INTENT(IN) ::  mode       !< either 'xy', 'xz' or 'yz'
    498        CHARACTER (LEN=*), INTENT(IN) ::  variable   !< name of variable
     500       CHARACTER (LEN=*), INTENT(IN)    ::  mode       !< either 'xy', 'xz' or 'yz'
     501       CHARACTER (LEN=*), INTENT(IN)    ::  variable   !< name of variable
    499502
    500503       INTEGER(iwp), INTENT(IN) ::  av        !< flag for (non-)average output
     
    508511
    509512       REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do), INTENT(INOUT) ::  local_pf !< local
    510           !< array to which output data is resorted to
     513                                                                                      !< array to which output data is resorted to
    511514
    512515!
     
    519522
    520523
    521 !------------------------------------------------------------------------------!
     524!--------------------------------------------------------------------------------------------------!
    522525!
    523526! Description:
    524527! ------------
    525528!> Subroutine defining 3D output variables
    526 !------------------------------------------------------------------------------!
     529!--------------------------------------------------------------------------------------------------!
    527530    SUBROUTINE gust_data_output_3d( av, variable, found, local_pf, fill_value, nzb_do, nzt_do )
    528531
     
    550553
    551554
    552 !------------------------------------------------------------------------------!
     555!--------------------------------------------------------------------------------------------------!
    553556! Description:
    554557! ------------
    555558!> This routine computes profile and timeseries data for the gust module.
    556 !------------------------------------------------------------------------------!
     559!--------------------------------------------------------------------------------------------------!
    557560    SUBROUTINE gust_statistics( mode, sr, tn, dots_max )
    558561
     
    573576
    574577
    575 !------------------------------------------------------------------------------!
     578!--------------------------------------------------------------------------------------------------!
    576579! Description:
    577580! ------------
    578581!> Read module-specific global restart data (Fortran binary format).
    579 !------------------------------------------------------------------------------!
     582!--------------------------------------------------------------------------------------------------!
    580583    SUBROUTINE gust_rrd_global_ftn( found )
    581584
    582585
    583        USE control_parameters,                                                 &
     586       USE control_parameters,                                                                     &
    584587           ONLY: length, restart_string
    585588
     
    608611
    609612
    610 !------------------------------------------------------------------------------!
     613!--------------------------------------------------------------------------------------------------!
    611614! Description:
    612615! ------------
    613616!> Read module-specific global restart data (MPI-IO).
    614 !------------------------------------------------------------------------------!
     617!--------------------------------------------------------------------------------------------------!
    615618    SUBROUTINE gust_rrd_global_mpi
    616619
     
    622625
    623626
    624 !------------------------------------------------------------------------------!
     627!--------------------------------------------------------------------------------------------------!
    625628! Description:
    626629! ------------
    627630!> Read module-specific local restart data arrays (Fortran binary format).
    628 !------------------------------------------------------------------------------!
    629     SUBROUTINE gust_rrd_local_ftn( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc,         &
    630                                    nxr_on_file, nynf, nync, nyn_on_file, nysf,     &
    631                                    nysc, nys_on_file, tmp_2d, tmp_3d, found )
     631!--------------------------------------------------------------------------------------------------!
     632    SUBROUTINE gust_rrd_local_ftn( k, nxlf, nxlc, nxl_on_file, nxrf, nxrc, nxr_on_file, nynf, nync,&
     633                                   nyn_on_file, nysf, nysc, nys_on_file, tmp_2d, tmp_3d, found )
    632634
    633635
     
    696698
    697699
    698 !------------------------------------------------------------------------------!
     700!--------------------------------------------------------------------------------------------------!
    699701! Description:
    700702! ------------
    701703!> Read module-specific local restart data arrays (MPI-IO).
    702 !------------------------------------------------------------------------------!
     704!--------------------------------------------------------------------------------------------------!
    703705    SUBROUTINE gust_rrd_local_mpi
    704706