Home


About us
PALM
PALM-4U
Installation
Release notes
Documentation
Getting Help
Gallery
Developers area

Changeset 4510 for palm


Ignore:
Timestamp:
Apr 29, 2020 2:19:18 PM (15 months ago)
Author:
raasch
Message:

files re-formatted to follow the PALM coding standard

Location:
palm/trunk/SOURCE
Files:
3 edited

Legend:

Unmodified
Added
Removed

  • palm/trunk/SOURCE/tridia_solver_mod.f90

    r4360 r4510  
    11!> @file tridia_solver_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!--------------------------------------------------------------------------------------------------!
     18!
    1919!
    2020! Current revisions:
    21 ! ------------------
     21! -----------------
    2222!
    2323!
    … …  
    2525! -----------------
    2626! $Id$
     27! file re-formatted to follow the PALM coding standard
     28!
     29! 4360 2020-01-07 11:25:50Z suehring
    2730! Added missing OpenMP directives
    28 ! 
     31!
    2932! 4182 2019-08-22 15:20:23Z scharf
    3033! Corrected "Former revisions" section
    31 ! 
     34!
    3235! 3761 2019-02-25 15:31:42Z raasch
    3336! OpenACC modification to prevent compiler warning about unused variable
    34 ! 
     37!
    3538! 3690 2019-01-22 22:56:42Z knoop
    3639! OpenACC port for SPEC
    … …  
    3841! 1212 2013-08-15 08:46:27Z raasch
    3942! Initial revision.
    40 ! Routines have been moved to seperate module from former file poisfft to here.
    41 ! The tridiagonal matrix coefficients of array tri are calculated only once at
    42 ! the beginning, i.e. routine split is called within tridia_init.
    43 !
    44 !
    45 ! Description:
    46 ! ------------
    47 !> solves the linear system of equations:
     43! Routines have been moved to seperate module from former file poisfft to here. The tridiagonal
     44! matrix coefficients of array tri are calculated only once at the beginning, i.e. routine split is
     45! called within tridia_init.
     46!
     47!
     48! Description:
     49! ------------
     50!> Solves the linear system of equations:
    4851!>
    49 !> -(4 pi^2(i^2/(dx^2*nnx^2)+j^2/(dy^2*nny^2))+
    50 !>   1/(dzu(k)*dzw(k))+1/(dzu(k-1)*dzw(k)))*p(i,j,k)+
     52!> -(4 pi^2(i^2/(dx^2*nnx^2)+j^2/(dy^2*nny^2))+ 1/(dzu(k)*dzw(k))+1/(dzu(k-1)*dzw(k)))*p(i,j,k)+
    5153!> 1/(dzu(k)*dzw(k))*p(i,j,k+1)+1/(dzu(k-1)*dzw(k))*p(i,j,k-1)=d(i,j,k)
    5254!>
    5355!> by using the Thomas algorithm
    54 !------------------------------------------------------------------------------!
     56!--------------------------------------------------------------------------------------------------!
    5557
    5658#define __acc_fft_device ( defined( _OPENACC ) && ( defined ( __cuda_fft ) ) )
    5759
    5860 MODULE tridia_solver
    59  
    60 
    61     USE basic_constants_and_equations_mod,                                     &
     61
     62
     63    USE basic_constants_and_equations_mod,                                                         &
    6264        ONLY:  pi
    6365
    64     USE indices,                                                               &
    65         ONLY:  nx, ny, nz
     66    USE indices,                                                                                   &
     67        ONLY:  nx,                                                                                 &
     68               ny,                                                                                 &
     69               nz
    6670
    6771    USE kinds
    6872
    69     USE transpose_indices,                                                     &
    70         ONLY:  nxl_z, nyn_z, nxr_z, nys_z
    71 
    72     IMPLICIT NONE
    73 
    74     REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ddzuw !<
     73    USE transpose_indices,                                                                         &
     74        ONLY:  nxl_z,                                                                              &
     75               nyn_z,                                                                              &
     76               nxr_z,                                                                              &
     77               nys_z
     78
     79    IMPLICIT NONE
     80
     81    REAL(wp), DIMENSION(:,:), ALLOCATABLE ::  ddzuw  !<
    7582
    7683    PRIVATE
    … …  
    8996
    9097
    91 !------------------------------------------------------------------------------!
     98!--------------------------------------------------------------------------------------------------!
    9299! Description:
    93100! ------------
    94101!> @todo Missing subroutine description.
    95 !------------------------------------------------------------------------------!
    96     SUBROUTINE tridia_init
    97 
    98        USE arrays_3d,                                                          &
    99            ONLY:  ddzu_pres, ddzw, rho_air_zw
     102!--------------------------------------------------------------------------------------------------!
     103 SUBROUTINE tridia_init
     104
     105    USE arrays_3d,                                                                                 &
     106        ONLY:  ddzu_pres,                                                                          &
     107               ddzw,                                                                               &
     108               rho_air_zw
    100109
    101110#if defined( _OPENACC )
    102        USE arrays_3d,                                                          &
     111       USE arrays_3d,                                                                              &
    103112           ONLY:  tri
    104113#endif
    105114
    106        IMPLICIT NONE
    107 
    108        INTEGER(iwp) ::  k !<
    109 
    110        ALLOCATE( ddzuw(0:nz-1,3) )
    111 
    112        DO  k = 0, nz-1
    113           ddzuw(k,1) = ddzu_pres(k+1) * ddzw(k+1) * rho_air_zw(k)
    114           ddzuw(k,2) = ddzu_pres(k+2) * ddzw(k+1) * rho_air_zw(k+1)
    115           ddzuw(k,3) = -1.0_wp * &
    116                        ( ddzu_pres(k+2) * ddzw(k+1) * rho_air_zw(k+1) +        &
    117                          ddzu_pres(k+1) * ddzw(k+1) * rho_air_zw(k) )
    118        ENDDO
    119 !
    120 !--    Calculate constant coefficients of the tridiagonal matrix
    121        CALL maketri
    122        CALL split
    123 
    124 #if __acc_fft_device
    125        !$ACC ENTER DATA &
    126        !$ACC COPYIN(ddzuw(0:nz-1,1:3)) &
    127        !$ACC COPYIN(tri(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1,1:2))
    128 #endif
    129 
    130     END SUBROUTINE tridia_init
    131 
    132 
    133 !------------------------------------------------------------------------------!
    134 ! Description:
    135 ! ------------
    136 !> Computes the i- and j-dependent component of the matrix
    137 !> Provide the constant coefficients of the tridiagonal matrix for solution
    138 !> of the Poisson equation in Fourier space.
    139 !> The coefficients are computed following the method of
    140 !> Schmidt et al. (DFVLR-Mitteilung 84-15), which departs from Stephan
    141 !> Siano's original version by discretizing the Poisson equation,
    142 !> before it is Fourier-transformed.
    143 !------------------------------------------------------------------------------!
    144     SUBROUTINE maketri
    145 
    146 
    147           USE arrays_3d,                                                       &
    148               ONLY:  tric, rho_air
    149 
    150           USE control_parameters,                                              &
    151               ONLY:  ibc_p_b, ibc_p_t
    152 
    153           USE grid_variables,                                                  &
    154               ONLY:  dx, dy
    155 
    156 
    157           IMPLICIT NONE
    158 
    159           INTEGER(iwp) ::  i    !<
    160           INTEGER(iwp) ::  j    !<
    161           INTEGER(iwp) ::  k    !<
    162           INTEGER(iwp) ::  nnxh !<
    163           INTEGER(iwp) ::  nnyh !<
    164 
    165           REAL(wp)    ::  ll(nxl_z:nxr_z,nys_z:nyn_z) !<
    166 
    167 
    168           nnxh = ( nx + 1 ) / 2
    169           nnyh = ( ny + 1 ) / 2
    170 
    171           DO  j = nys_z, nyn_z
    172              DO  i = nxl_z, nxr_z
    173                 IF ( j >= 0  .AND.  j <= nnyh )  THEN
    174                    IF ( i >= 0  .AND.  i <= nnxh )  THEN
    175                       ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * i ) / &
    176                                             REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + &
    177                                 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) / &
    178                                             REAL( ny+1, KIND=wp ) ) ) / ( dy * dy )
    179                    ELSE
    180                       ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( nx+1-i ) ) / &
    181                                             REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + &
    182                                 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) / &
    183                                             REAL( ny+1, KIND=wp ) ) ) / ( dy * dy )
    184                    ENDIF
    185                 ELSE
    186                    IF ( i >= 0  .AND.  i <= nnxh )  THEN
    187                       ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * i ) / &
    188                                             REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + &
    189                                 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( ny+1-j ) ) / &
    190                                             REAL( ny+1, KIND=wp ) ) ) / ( dy * dy )
    191                    ELSE
    192                       ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( nx+1-i ) ) / &
    193                                             REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) + &
    194                                 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( ny+1-j ) ) / &
    195                                             REAL( ny+1, KIND=wp ) ) ) / ( dy * dy )
    196                    ENDIF
    197                 ENDIF
    198              ENDDO
    199           ENDDO
    200 
    201           DO  k = 0, nz-1
    202              DO  j = nys_z, nyn_z
    203                 DO  i = nxl_z, nxr_z
    204                    tric(i,j,k) = ddzuw(k,3) - ll(i,j) * rho_air(k+1)
    205                 ENDDO
    206              ENDDO
    207           ENDDO
    208 
    209           IF ( ibc_p_b == 1 )  THEN
    210              DO  j = nys_z, nyn_z
    211                 DO  i = nxl_z, nxr_z
    212                    tric(i,j,0) = tric(i,j,0) + ddzuw(0,1)
    213                 ENDDO
    214              ENDDO
    215           ENDIF
    216           IF ( ibc_p_t == 1 )  THEN
    217              DO  j = nys_z, nyn_z
    218                 DO  i = nxl_z, nxr_z
    219                    tric(i,j,nz-1) = tric(i,j,nz-1) + ddzuw(nz-1,2)
    220                 ENDDO
    221              ENDDO
    222           ENDIF
    223 
    224     END SUBROUTINE maketri
    225 
    226 
    227 !------------------------------------------------------------------------------!
    228 ! Description:
    229 ! ------------
    230 !> Substitution (Forward and Backward) (Thomas algorithm)
    231 !------------------------------------------------------------------------------!
    232     SUBROUTINE tridia_substi( ar )
    233 
    234 
    235           USE arrays_3d,                                                       &
    236               ONLY:  tri
    237 
    238           USE control_parameters,                                              &
    239               ONLY:  ibc_p_b, ibc_p_t
    240 
    241           IMPLICIT NONE
    242 
    243           INTEGER(iwp) ::  i !<
    244           INTEGER(iwp) ::  j !<
    245           INTEGER(iwp) ::  k !<
    246 
    247           REAL(wp)     ::  ar(nxl_z:nxr_z,nys_z:nyn_z,1:nz) !<
    248 
    249           REAL(wp), DIMENSION(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1)   ::  ar1 !<
    250 #if __acc_fft_device
    251           !$ACC DECLARE CREATE(ar1)
    252 #endif
    253 
    254           !$OMP PARALLEL PRIVATE(i,j,k)
    255 
    256 !
    257 !--       Forward substitution
    258 #if __acc_fft_device
    259           !$ACC PARALLEL PRESENT(ar, ar1, tri) PRIVATE(i,j,k)
    260 #endif
    261           DO  k = 0, nz - 1
    262 #if __acc_fft_device
    263              !$ACC LOOP COLLAPSE(2)
    264 #endif
    265              !$OMP DO
    266              DO  j = nys_z, nyn_z
    267                 DO  i = nxl_z, nxr_z
    268 
    269                    IF ( k == 0 )  THEN
    270                       ar1(i,j,k) = ar(i,j,k+1)
    271                    ELSE
    272                       ar1(i,j,k) = ar(i,j,k+1) - tri(i,j,k,2) * ar1(i,j,k-1)
    273                    ENDIF
    274 
    275                 ENDDO
    276              ENDDO
    277           ENDDO
    278 #if __acc_fft_device
    279           !$ACC END PARALLEL
    280 #endif
    281 
    282 !
    283 !--       Backward substitution
    284 !--       Note, the 1.0E-20 in the denominator is due to avoid divisions
    285 !--       by zero appearing if the pressure bc is set to neumann at the top of
    286 !--       the model domain.
    287 #if __acc_fft_device
    288           !$ACC PARALLEL PRESENT(ar, ar1, ddzuw, tri) PRIVATE(i,j,k)
    289 #endif
    290           DO  k = nz-1, 0, -1
    291 #if __acc_fft_device
    292              !$ACC LOOP COLLAPSE(2)
    293 #endif
    294              !$OMP DO
    295              DO  j = nys_z, nyn_z
    296                 DO  i = nxl_z, nxr_z
    297 
    298                    IF ( k == nz-1 )  THEN
    299                       ar(i,j,k+1) = ar1(i,j,k) / ( tri(i,j,k,1) + 1.0E-20_wp )
    300                    ELSE
    301                       ar(i,j,k+1) = ( ar1(i,j,k) - ddzuw(k,2) * ar(i,j,k+2) ) &
    302                               / tri(i,j,k,1)
    303                    ENDIF
    304                 ENDDO
    305              ENDDO
    306           ENDDO
    307 #if __acc_fft_device
    308           !$ACC END PARALLEL
    309 #endif
    310 
    311           !$OMP END PARALLEL
    312 
    313 !
    314 !--       Indices i=0, j=0 correspond to horizontally averaged pressure.
    315 !--       The respective values of ar should be zero at all k-levels if
    316 !--       acceleration of horizontally averaged vertical velocity is zero.
    317           IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1 )  THEN
    318              IF ( nys_z == 0  .AND.  nxl_z == 0 )  THEN
    319 #if __acc_fft_device
    320                 !$ACC PARALLEL LOOP PRESENT(ar)
    321 #endif
    322                 DO  k = 1, nz
    323                    ar(nxl_z,nys_z,k) = 0.0_wp
    324                 ENDDO
     115    IMPLICIT NONE
     116
     117    INTEGER(iwp) ::  k  !<
     118
     119    ALLOCATE( ddzuw(0:nz-1,3) )
     120
     121    DO  k = 0, nz-1
     122       ddzuw(k,1) = ddzu_pres(k+1) * ddzw(k+1) * rho_air_zw(k)
     123       ddzuw(k,2) = ddzu_pres(k+2) * ddzw(k+1) * rho_air_zw(k+1)
     124       ddzuw(k,3) = -1.0_wp * ( ddzu_pres(k+2) * ddzw(k+1) * rho_air_zw(k+1) +                     &
     125                                ddzu_pres(k+1) * ddzw(k+1) * rho_air_zw(k) )
     126    ENDDO
     127!
     128!-- Calculate constant coefficients of the tridiagonal matrix
     129    CALL maketri
     130    CALL split
     131
     132#if __acc_fft_device
     133    !$ACC ENTER DATA &
     134    !$ACC COPYIN(ddzuw(0:nz-1,1:3)) &
     135    !$ACC COPYIN(tri(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1,1:2))
     136#endif
     137
     138 END SUBROUTINE tridia_init
     139
     140
     141!--------------------------------------------------------------------------------------------------!
     142! Description:
     143! ------------
     144!> Computes the i- and j-dependent component of the matrix.
     145!> Provide the constant coefficients of the tridiagonal matrix for solution of the Poisson equation
     146!> in Fourier space. The coefficients are computed following the method of Schmidt et al.
     147!> (DFVLR-Mitteilung 84-15), which departs from Stephan Siano's original version by discretizing the
     148!> Poisson equation, before it is Fourier-transformed.
     149!--------------------------------------------------------------------------------------------------!
     150 SUBROUTINE maketri
     151
     152
     153    USE arrays_3d,                                                                                 &
     154        ONLY:  tric,                                                                               &
     155               rho_air
     156
     157    USE control_parameters,                                                                        &
     158        ONLY:  ibc_p_b,                                                                            &
     159               ibc_p_t
     160
     161    USE grid_variables,                                                                            &
     162        ONLY:  dx,                                                                                 &
     163               dy
     164
     165
     166    IMPLICIT NONE
     167
     168    INTEGER(iwp) ::  i     !<
     169    INTEGER(iwp) ::  j     !<
     170    INTEGER(iwp) ::  k     !<
     171    INTEGER(iwp) ::  nnxh  !<
     172    INTEGER(iwp) ::  nnyh  !<
     173
     174    REAL(wp)    ::  ll(nxl_z:nxr_z,nys_z:nyn_z)  !<
     175
     176
     177    nnxh = ( nx + 1 ) / 2
     178    nnyh = ( ny + 1 ) / 2
     179
     180    DO  j = nys_z, nyn_z
     181       DO  i = nxl_z, nxr_z
     182          IF ( j >= 0  .AND.  j <= nnyh )  THEN
     183             IF ( i >= 0  .AND.  i <= nnxh )  THEN
     184                ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * i ) /                           &
     185                                     REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) +                     &
     186                          2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) /                           &
     187                                     REAL( ny+1, KIND=wp ) ) ) / ( dy * dy )
     188             ELSE
     189                ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( nx+1-i ) ) /                  &
     190                                     REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) +                     &
     191                          2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) /                           &
     192                                     REAL( ny+1, KIND=wp ) ) ) / ( dy * dy )
     193             ENDIF
     194          ELSE
     195             IF ( i >= 0  .AND.  i <= nnxh )  THEN
     196                ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * i ) /                           &
     197                                     REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) +                     &
     198                          2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( ny+1-j ) ) /                  &
     199                                     REAL( ny+1, KIND=wp ) ) ) / ( dy * dy )
     200             ELSE
     201                ll(i,j) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( nx+1-i ) ) /                  &
     202                                     REAL( nx+1, KIND=wp ) ) ) / ( dx * dx ) +                     &
     203                          2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( ny+1-j ) ) /                  &
     204                                     REAL( ny+1, KIND=wp ) ) ) / ( dy * dy )
    325205             ENDIF
    326206          ENDIF
    327 
    328     END SUBROUTINE tridia_substi
    329 
    330 
    331 !------------------------------------------------------------------------------!
    332 ! Description:
    333 ! ------------
    334 !> Substitution (Forward and Backward) (Thomas algorithm)
    335 !------------------------------------------------------------------------------!
    336     SUBROUTINE tridia_substi_overlap( ar, jj )
    337 
    338 
    339           USE arrays_3d,                                                       &
    340               ONLY:  tri
    341 
    342           USE control_parameters,                                              &
    343               ONLY:  ibc_p_b, ibc_p_t
    344 
    345           IMPLICIT NONE
    346 
    347           INTEGER(iwp) ::  i  !<
    348           INTEGER(iwp) ::  j  !<
    349           INTEGER(iwp) ::  jj !<
    350           INTEGER(iwp) ::  k  !<
    351 
    352           REAL(wp)     ::  ar(nxl_z:nxr_z,nys_z:nyn_z,1:nz) !<
    353 
    354           REAL(wp), DIMENSION(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1) ::  ar1 !<
    355 
    356 !
    357 !--       Forward substitution
    358           DO  k = 0, nz - 1
    359              DO  j = nys_z, nyn_z
    360                 DO  i = nxl_z, nxr_z
    361 
    362                    IF ( k == 0 )  THEN
    363                       ar1(i,j,k) = ar(i,j,k+1)
    364                    ELSE
    365                       ar1(i,j,k) = ar(i,j,k+1) - tri(i,jj,k,2) * ar1(i,j,k-1)
    366                    ENDIF
    367 
    368                 ENDDO
    369              ENDDO
    370           ENDDO
    371 
    372 !
    373 !--       Backward substitution
    374 !--       Note, the 1.0E-20 in the denominator is due to avoid divisions
    375 !--       by zero appearing if the pressure bc is set to neumann at the top of
    376 !--       the model domain.
    377           DO  k = nz-1, 0, -1
    378              DO  j = nys_z, nyn_z
    379                 DO  i = nxl_z, nxr_z
    380 
    381                    IF ( k == nz-1 )  THEN
    382                       ar(i,j,k+1) = ar1(i,j,k) / ( tri(i,jj,k,1) + 1.0E-20_wp )
    383                    ELSE
    384                       ar(i,j,k+1) = ( ar1(i,j,k) - ddzuw(k,2) * ar(i,j,k+2) ) &
    385                               / tri(i,jj,k,1)
    386                    ENDIF
    387                 ENDDO
    388              ENDDO
    389           ENDDO
    390 
    391 !
    392 !--       Indices i=0, j=0 correspond to horizontally averaged pressure.
    393 !--       The respective values of ar should be zero at all k-levels if
    394 !--       acceleration of horizontally averaged vertical velocity is zero.
    395           IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1 )  THEN
    396              IF ( nys_z == 0  .AND.  nxl_z == 0 )  THEN
    397                 DO  k = 1, nz
    398                    ar(nxl_z,nys_z,k) = 0.0_wp
    399                 ENDDO
     207       ENDDO
     208    ENDDO
     209
     210    DO  k = 0, nz-1
     211       DO  j = nys_z, nyn_z
     212          DO  i = nxl_z, nxr_z
     213             tric(i,j,k) = ddzuw(k,3) - ll(i,j) * rho_air(k+1)
     214          ENDDO
     215       ENDDO
     216    ENDDO
     217
     218    IF ( ibc_p_b == 1 )  THEN
     219       DO  j = nys_z, nyn_z
     220          DO  i = nxl_z, nxr_z
     221             tric(i,j,0) = tric(i,j,0) + ddzuw(0,1)
     222          ENDDO
     223       ENDDO
     224    ENDIF
     225    IF ( ibc_p_t == 1 )  THEN
     226       DO  j = nys_z, nyn_z
     227          DO  i = nxl_z, nxr_z
     228             tric(i,j,nz-1) = tric(i,j,nz-1) + ddzuw(nz-1,2)
     229          ENDDO
     230       ENDDO
     231    ENDIF
     232
     233 END SUBROUTINE maketri
     234
     235
     236!--------------------------------------------------------------------------------------------------!
     237! Description:
     238! ------------
     239!> Substitution (Forward and Backward) (Thomas algorithm).
     240!--------------------------------------------------------------------------------------------------!
     241 SUBROUTINE tridia_substi( ar )
     242
     243
     244    USE arrays_3d,                                                                                 &
     245        ONLY:  tri
     246
     247    USE control_parameters,                                                                        &
     248        ONLY:  ibc_p_b,                                                                            &
     249               ibc_p_t
     250
     251    IMPLICIT NONE
     252
     253    INTEGER(iwp) ::  i !<
     254    INTEGER(iwp) ::  j !<
     255    INTEGER(iwp) ::  k !<
     256
     257    REAL(wp)     ::  ar(nxl_z:nxr_z,nys_z:nyn_z,1:nz)  !<
     258
     259    REAL(wp), DIMENSION(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1) ::  ar1  !<
     260#if __acc_fft_device
     261    !$ACC DECLARE CREATE(ar1)
     262#endif
     263
     264    !$OMP PARALLEL PRIVATE(i,j,k)
     265
     266!
     267!-- Forward substitution
     268#if __acc_fft_device
     269    !$ACC PARALLEL PRESENT(ar, ar1, tri) PRIVATE(i,j,k)
     270#endif
     271    DO  k = 0, nz - 1
     272#if __acc_fft_device
     273       !$ACC LOOP COLLAPSE(2)
     274#endif
     275       !$OMP DO
     276       DO  j = nys_z, nyn_z
     277          DO  i = nxl_z, nxr_z
     278
     279             IF ( k == 0 )  THEN
     280                ar1(i,j,k) = ar(i,j,k+1)
     281             ELSE
     282                ar1(i,j,k) = ar(i,j,k+1) - tri(i,j,k,2) * ar1(i,j,k-1)
    400283             ENDIF
    401           ENDIF
    402 
    403     END SUBROUTINE tridia_substi_overlap
    404 
    405 
    406 !------------------------------------------------------------------------------!
    407 ! Description:
    408 ! ------------
    409 !> Splitting of the tridiagonal matrix (Thomas algorithm)
    410 !------------------------------------------------------------------------------!
    411     SUBROUTINE split
    412 
    413 
    414           USE arrays_3d,                                                       &
    415               ONLY:  tri, tric
    416 
    417           IMPLICIT NONE
    418 
    419           INTEGER(iwp) ::  i !<
    420           INTEGER(iwp) ::  j !<
    421           INTEGER(iwp) ::  k !<
    422 !
    423 !--       Splitting
    424           DO  j = nys_z, nyn_z
    425              DO  i = nxl_z, nxr_z
    426                 tri(i,j,0,1) = tric(i,j,0)
    427              ENDDO
    428           ENDDO
    429 
    430           DO  k = 1, nz-1
    431              DO  j = nys_z, nyn_z
    432                 DO  i = nxl_z, nxr_z
    433                    tri(i,j,k,2) = ddzuw(k,1) / tri(i,j,k-1,1)
    434                    tri(i,j,k,1) = tric(i,j,k) - ddzuw(k-1,2) * tri(i,j,k,2)
    435                 ENDDO
    436              ENDDO
    437           ENDDO
    438 
    439     END SUBROUTINE split
    440 
    441 
    442 !------------------------------------------------------------------------------!
    443 ! Description:
    444 ! ------------
    445 !> Solves the linear system of equations for a 1d-decomposition along x (see
    446 !> tridia)
     284
     285          ENDDO
     286       ENDDO
     287    ENDDO
     288#if __acc_fft_device
     289    !$ACC END PARALLEL
     290#endif
     291
     292!
     293!-- Backward substitution
     294!-- Note, the 1.0E-20 in the denominator is due to avoid divisions by zero appearing if the
     295!-- pressure bc is set to neumann at the top of the model domain.
     296#if __acc_fft_device
     297    !$ACC PARALLEL PRESENT(ar, ar1, ddzuw, tri) PRIVATE(i,j,k)
     298#endif
     299    DO  k = nz-1, 0, -1
     300#if __acc_fft_device
     301       !$ACC LOOP COLLAPSE(2)
     302#endif
     303       !$OMP DO
     304       DO  j = nys_z, nyn_z
     305          DO  i = nxl_z, nxr_z
     306
     307             IF ( k == nz-1 )  THEN
     308                ar(i,j,k+1) = ar1(i,j,k) / ( tri(i,j,k,1) + 1.0E-20_wp )
     309             ELSE
     310                ar(i,j,k+1) = ( ar1(i,j,k) - ddzuw(k,2) * ar(i,j,k+2) ) / tri(i,j,k,1)
     311             ENDIF
     312          ENDDO
     313       ENDDO
     314    ENDDO
     315#if __acc_fft_device
     316    !$ACC END PARALLEL
     317#endif
     318
     319    !$OMP END PARALLEL
     320
     321!
     322!-- Indices i=0, j=0 correspond to horizontally averaged pressure. The respective values of ar
     323!-- should be zero at all k-levels if acceleration of horizontally averaged vertical velocity
     324!-- is zero.
     325    IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1 )  THEN
     326       IF ( nys_z == 0  .AND.  nxl_z == 0 )  THEN
     327#if __acc_fft_device
     328             !$ACC PARALLEL LOOP PRESENT(ar)
     329#endif
     330          DO  k = 1, nz
     331             ar(nxl_z,nys_z,k) = 0.0_wp
     332          ENDDO
     333       ENDIF
     334    ENDIF
     335
     336 END SUBROUTINE tridia_substi
     337
     338
     339!--------------------------------------------------------------------------------------------------!
     340! Description:
     341! ------------
     342!> Substitution (Forward and Backward) (Thomas algorithm).
     343!--------------------------------------------------------------------------------------------------!
     344 SUBROUTINE tridia_substi_overlap( ar, jj )
     345
     346
     347    USE arrays_3d,                                                                                 &
     348        ONLY:  tri
     349
     350    USE control_parameters,                                                                        &
     351        ONLY:  ibc_p_b,                                                                            &
     352               ibc_p_t
     353
     354    IMPLICIT NONE
     355
     356    INTEGER(iwp) ::  i  !<
     357    INTEGER(iwp) ::  j  !<
     358    INTEGER(iwp) ::  jj !<
     359    INTEGER(iwp) ::  k  !<
     360
     361    REAL(wp)     ::  ar(nxl_z:nxr_z,nys_z:nyn_z,1:nz)  !<
     362
     363    REAL(wp), DIMENSION(nxl_z:nxr_z,nys_z:nyn_z,0:nz-1) ::  ar1  !<
     364
     365!
     366!-- Forward substitution
     367    DO  k = 0, nz - 1
     368       DO  j = nys_z, nyn_z
     369          DO  i = nxl_z, nxr_z
     370
     371             IF ( k == 0 )  THEN
     372                ar1(i,j,k) = ar(i,j,k+1)
     373             ELSE
     374                ar1(i,j,k) = ar(i,j,k+1) - tri(i,jj,k,2) * ar1(i,j,k-1)
     375             ENDIF
     376
     377          ENDDO
     378       ENDDO
     379    ENDDO
     380
     381!
     382!-- Backward substitution
     383!-- Note, the 1.0E-20 in the denominator is due to avoid divisions by zero appearing if the
     384!-- pressure bc is set to neumann at the top of the model domain.
     385    DO  k = nz-1, 0, -1
     386       DO  j = nys_z, nyn_z
     387          DO  i = nxl_z, nxr_z
     388
     389             IF ( k == nz-1 )  THEN
     390                ar(i,j,k+1) = ar1(i,j,k) / ( tri(i,jj,k,1) + 1.0E-20_wp )
     391             ELSE
     392                ar(i,j,k+1) = ( ar1(i,j,k) - ddzuw(k,2) * ar(i,j,k+2) ) / tri(i,jj,k,1)
     393             ENDIF
     394          ENDDO
     395       ENDDO
     396    ENDDO
     397
     398!
     399!-- Indices i=0, j=0 correspond to horizontally averaged pressure. The respective values of ar
     400!-- should be zero at all k-levels if acceleration of horizontally averaged vertical velocity
     401!-- is zero.
     402    IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1 )  THEN
     403       IF ( nys_z == 0  .AND.  nxl_z == 0 )  THEN
     404          DO  k = 1, nz
     405             ar(nxl_z,nys_z,k) = 0.0_wp
     406          ENDDO
     407       ENDIF
     408    ENDIF
     409
     410 END SUBROUTINE tridia_substi_overlap
     411
     412
     413!--------------------------------------------------------------------------------------------------!
     414! Description:
     415! ------------
     416!> Splitting of the tridiagonal matrix (Thomas algorithm).
     417!--------------------------------------------------------------------------------------------------!
     418 SUBROUTINE split
     419
     420
     421    USE arrays_3d,                                                                                 &
     422        ONLY:  tri,                                                                                &
     423               tric
     424
     425    IMPLICIT NONE
     426
     427    INTEGER(iwp) ::  i  !<
     428    INTEGER(iwp) ::  j  !<
     429    INTEGER(iwp) ::  k  !<
     430!
     431!   Splitting
     432    DO  j = nys_z, nyn_z
     433       DO  i = nxl_z, nxr_z
     434          tri(i,j,0,1) = tric(i,j,0)
     435       ENDDO
     436    ENDDO
     437
     438    DO  k = 1, nz-1
     439       DO  j = nys_z, nyn_z
     440          DO  i = nxl_z, nxr_z
     441             tri(i,j,k,2) = ddzuw(k,1) / tri(i,j,k-1,1)
     442             tri(i,j,k,1) = tric(i,j,k) - ddzuw(k-1,2) * tri(i,j,k,2)
     443          ENDDO
     444       ENDDO
     445    ENDDO
     446
     447 END SUBROUTINE split
     448
     449
     450!--------------------------------------------------------------------------------------------------!
     451! Description:
     452! ------------
     453!> Solves the linear system of equations for a 1d-decomposition along x (see tridia).
    447454!>
    448 !> @attention when using the intel compilers older than 12.0, array tri must
    449 !>            be passed as an argument to the contained subroutines. Otherwise
    450 !>            addres faults will occur. This feature can be activated with
    451 !>            cpp-switch __intel11
    452 !>            On NEC, tri should not be passed (except for routine substi_1dd)
    453 !>            because this causes very bad performance.
    454 !------------------------------------------------------------------------------!
    455  
    456     SUBROUTINE tridia_1dd( ddx2, ddy2, nx, ny, j, ar, tri_for_1d )
    457 
    458 
    459        USE arrays_3d,                                                          &
    460            ONLY:  ddzu_pres, ddzw, rho_air, rho_air_zw
    461 
    462        USE control_parameters,                                                 &
    463            ONLY:  ibc_p_b, ibc_p_t
    464 
    465        IMPLICIT NONE
    466 
    467        INTEGER(iwp) ::  i                  !<
    468        INTEGER(iwp) ::  j                  !<
    469        INTEGER(iwp) ::  k                  !<
    470        INTEGER(iwp) ::  nnyh               !<
    471        INTEGER(iwp) ::  nx                 !<
    472        INTEGER(iwp) ::  ny                 !<
    473 
    474        REAL(wp)     ::  ddx2 !<
    475        REAL(wp)     ::  ddy2 !<
    476 
    477        REAL(wp), DIMENSION(0:nx,1:nz)     ::  ar         !<
    478        REAL(wp), DIMENSION(5,0:nx,0:nz-1) ::  tri_for_1d !<
    479 
    480 
    481        nnyh = ( ny + 1 ) / 2
    482 
    483 !
    484 !--    Define constant elements of the tridiagonal matrix.
    485 !--    The compiler on SX6 does loop exchange. If 0:nx is a high power of 2,
    486 !--    the exchanged loops create bank conflicts. The following directive
    487 !--    prohibits loop exchange and the loops perform much better.
     455!> @attention When using intel compilers older than 12.0, array tri must be passed as an argument to
     456!>            the contained subroutines. Otherwise address faults will occur. This feature can be
     457!>            activated with cpp-switch __intel11. On NEC, tri should not be passed
     458!>            (except for routine substi_1dd) because this causes very bad performance.
     459!--------------------------------------------------------------------------------------------------!
     460
     461 SUBROUTINE tridia_1dd( ddx2, ddy2, nx, ny, j, ar, tri_for_1d )
     462
     463
     464    USE arrays_3d,                                                                                 &
     465        ONLY:  ddzu_pres,                                                                          &
     466               ddzw,                                                                               &
     467               rho_air,                                                                            &
     468               rho_air_zw
     469
     470    USE control_parameters,                                                                        &
     471        ONLY:  ibc_p_b,                                                                            &
     472               ibc_p_t
     473
     474    IMPLICIT NONE
     475
     476    INTEGER(iwp) ::  i     !<
     477    INTEGER(iwp) ::  j     !<
     478    INTEGER(iwp) ::  k     !<
     479    INTEGER(iwp) ::  nnyh  !<
     480    INTEGER(iwp) ::  nx    !<
     481    INTEGER(iwp) ::  ny    !<
     482
     483    REAL(wp)     ::  ddx2  !<
     484    REAL(wp)     ::  ddy2  !<
     485
     486    REAL(wp), DIMENSION(0:nx,1:nz)     ::  ar          !<
     487    REAL(wp), DIMENSION(5,0:nx,0:nz-1) ::  tri_for_1d  !<
     488
     489
     490    nnyh = ( ny + 1 ) / 2
     491
     492!
     493!-- Define constant elements of the tridiagonal matrix. The compiler on SX6 does loop exchange.
     494!-- If 0:nx is a high power of 2, the exchanged loops create bank conflicts. The following directive
     495!-- prohibits loop exchange and the loops perform much better.
    488496!CDIR NOLOOPCHG
    489        DO  k = 0, nz-1
    490           DO  i = 0,nx
    491              tri_for_1d(2,i,k) = ddzu_pres(k+1) * ddzw(k+1) * rho_air_zw(k)
    492              tri_for_1d(3,i,k) = ddzu_pres(k+2) * ddzw(k+1) * rho_air_zw(k+1)
    493           ENDDO
    494        ENDDO
    495 
    496        IF ( j <= nnyh )  THEN
    497           CALL maketri_1dd( j )
     497    DO  k = 0, nz-1
     498       DO  i = 0,nx
     499          tri_for_1d(2,i,k) = ddzu_pres(k+1) * ddzw(k+1) * rho_air_zw(k)
     500          tri_for_1d(3,i,k) = ddzu_pres(k+2) * ddzw(k+1) * rho_air_zw(k+1)
     501       ENDDO
     502    ENDDO
     503
     504    IF ( j <= nnyh )  THEN
     505       CALL maketri_1dd( j )
     506    ELSE
     507       CALL maketri_1dd( ny+1-j )
     508    ENDIF
     509
     510    CALL split_1dd
     511    CALL substi_1dd( ar, tri_for_1d )
     512
     513 CONTAINS
     514
     515
     516!--------------------------------------------------------------------------------------------------!
     517! Description:
     518! ------------
     519!> Computes the i- and j-dependent component of the matrix.
     520!--------------------------------------------------------------------------------------------------!
     521 SUBROUTINE maketri_1dd( j )
     522
     523    IMPLICIT NONE
     524
     525    INTEGER(iwp) ::  i    !<
     526    INTEGER(iwp) ::  j    !<
     527    INTEGER(iwp) ::  k    !<
     528    INTEGER(iwp) ::  nnxh !<
     529
     530    REAL(wp)     ::  a !<
     531    REAL(wp)     ::  c !<
     532
     533    REAL(wp), DIMENSION(0:nx) ::  l !<
     534
     535
     536    nnxh = ( nx + 1 ) / 2
     537!
     538!-- Provide the tridiagonal matrix for solution of the Poisson equation in Fourier space.
     539!-- The coefficients are computed following the method of Schmidt et al. (DFVLR-Mitteilung 84-15),
     540!-- which departs from Stephan Siano's original version by discretizing the Poisson equation,
     541!-- before it is Fourier-transformed.
     542    DO  i = 0, nx
     543       IF ( i >= 0 .AND. i <= nnxh ) THEN
     544          l(i) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * i ) /                                    &
     545                            REAL( nx+1, KIND=wp ) ) ) * ddx2 +                                     &
     546                 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) /                                    &
     547                            REAL( ny+1, KIND=wp ) ) ) * ddy2
    498548       ELSE
    499           CALL maketri_1dd( ny+1-j )
     549          l(i) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( nx+1-i ) ) /                           &
     550                            REAL( nx+1, KIND=wp ) ) ) * ddx2 +                                     &
     551                 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) /                                    &
     552                            REAL( ny+1, KIND=wp ) ) ) * ddy2
    500553       ENDIF
    501 
    502        CALL split_1dd
    503        CALL substi_1dd( ar, tri_for_1d )
    504 
    505     CONTAINS
    506 
    507 
    508 !------------------------------------------------------------------------------!
    509 ! Description:
    510 ! ------------
    511 !> computes the i- and j-dependent component of the matrix
    512 !------------------------------------------------------------------------------!
    513        SUBROUTINE maketri_1dd( j )
    514 
    515           IMPLICIT NONE
    516 
    517           INTEGER(iwp) ::  i    !<
    518           INTEGER(iwp) ::  j    !<
    519           INTEGER(iwp) ::  k    !<
    520           INTEGER(iwp) ::  nnxh !<
    521 
    522           REAL(wp)     ::  a !<
    523           REAL(wp)     ::  c !<
    524 
    525           REAL(wp), DIMENSION(0:nx) ::  l !<
    526 
    527 
    528           nnxh = ( nx + 1 ) / 2
    529 !
    530 !--       Provide the tridiagonal matrix for solution of the Poisson equation in
    531 !--       Fourier space. The coefficients are computed following the method of
    532 !--       Schmidt et al. (DFVLR-Mitteilung 84-15), which departs from Stephan
    533 !--       Siano's original version by discretizing the Poisson equation,
    534 !--       before it is Fourier-transformed
    535           DO  i = 0, nx
    536              IF ( i >= 0 .AND. i <= nnxh ) THEN
    537                 l(i) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * i ) / &
    538                                    REAL( nx+1, KIND=wp ) ) ) * ddx2 + &
    539                        2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) / &
    540                                    REAL( ny+1, KIND=wp ) ) ) * ddy2
    541              ELSE
    542                 l(i) = 2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * ( nx+1-i ) ) / &
    543                                    REAL( nx+1, KIND=wp ) ) ) * ddx2 + &
    544                        2.0_wp * ( 1.0_wp - COS( ( 2.0_wp * pi * j ) / &
    545                                    REAL( ny+1, KIND=wp ) ) ) * ddy2
    546              ENDIF
    547           ENDDO
    548 
    549           DO  k = 0, nz-1
    550              DO  i = 0, nx
    551                 a = -1.0_wp * ddzu_pres(k+2) * ddzw(k+1) * rho_air_zw(k+1)
    552                 c = -1.0_wp * ddzu_pres(k+1) * ddzw(k+1) * rho_air_zw(k)
    553                 tri_for_1d(1,i,k) = a + c - l(i) * rho_air(k+1)
    554              ENDDO
    555           ENDDO
    556           IF ( ibc_p_b == 1 )  THEN
    557              DO  i = 0, nx
    558                 tri_for_1d(1,i,0) = tri_for_1d(1,i,0) + tri_for_1d(2,i,0)
    559              ENDDO
    560           ENDIF
    561           IF ( ibc_p_t == 1 )  THEN
    562              DO  i = 0, nx
    563                 tri_for_1d(1,i,nz-1) = tri_for_1d(1,i,nz-1) + tri_for_1d(3,i,nz-1)
    564              ENDDO
    565           ENDIF
    566 
    567        END SUBROUTINE maketri_1dd
    568 
    569 
    570 !------------------------------------------------------------------------------!
    571 ! Description:
    572 ! ------------
    573 !> Splitting of the tridiagonal matrix (Thomas algorithm)
    574 !------------------------------------------------------------------------------!
    575        SUBROUTINE split_1dd
    576 
    577           IMPLICIT NONE
    578 
    579           INTEGER(iwp) ::  i !<
    580           INTEGER(iwp) ::  k !<
    581 
    582 
    583 !
    584 !--       Splitting
    585           DO  i = 0, nx
    586              tri_for_1d(4,i,0) = tri_for_1d(1,i,0)
    587           ENDDO
    588           DO  k = 1, nz-1
    589              DO  i = 0, nx
    590                 tri_for_1d(5,i,k) = tri_for_1d(2,i,k) / tri_for_1d(4,i,k-1)
    591                 tri_for_1d(4,i,k) = tri_for_1d(1,i,k) - tri_for_1d(3,i,k-1) * tri_for_1d(5,i,k)
    592              ENDDO
    593           ENDDO
    594 
    595        END SUBROUTINE split_1dd
    596 
    597 
    598 !------------------------------------------------------------------------------!
    599 ! Description:
    600 ! ------------
    601 !> Substitution (Forward and Backward) (Thomas algorithm)
    602 !------------------------------------------------------------------------------!
    603        SUBROUTINE substi_1dd( ar, tri_for_1d )
    604 
    605 
    606           IMPLICIT NONE
    607 
    608           INTEGER(iwp) ::  i !<
    609           INTEGER(iwp) ::  k !<
    610 
    611           REAL(wp), DIMENSION(0:nx,nz)       ::  ar         !<
    612           REAL(wp), DIMENSION(0:nx,0:nz-1)   ::  ar1        !<
    613           REAL(wp), DIMENSION(5,0:nx,0:nz-1) ::  tri_for_1d !<
    614 
    615 !
    616 !--       Forward substitution
    617           DO  i = 0, nx
    618              ar1(i,0) = ar(i,1)
    619           ENDDO
    620           DO  k = 1, nz-1
    621              DO  i = 0, nx
    622                 ar1(i,k) = ar(i,k+1) - tri_for_1d(5,i,k) * ar1(i,k-1)
    623              ENDDO
    624           ENDDO
    625 
    626 !
    627 !--       Backward substitution
    628 !--       Note, the add of 1.0E-20 in the denominator is due to avoid divisions
    629 !--       by zero appearing if the pressure bc is set to neumann at the top of
    630 !--       the model domain.
    631           DO  i = 0, nx
    632              ar(i,nz) = ar1(i,nz-1) / ( tri_for_1d(4,i,nz-1) + 1.0E-20_wp )
    633           ENDDO
    634           DO  k = nz-2, 0, -1
    635              DO  i = 0, nx
    636                 ar(i,k+1) = ( ar1(i,k) - tri_for_1d(3,i,k) * ar(i,k+2) ) &
    637                             / tri_for_1d(4,i,k)
    638              ENDDO
    639           ENDDO
    640 
    641 !
    642 !--       Indices i=0, j=0 correspond to horizontally averaged pressure.
    643 !--       The respective values of ar should be zero at all k-levels if
    644 !--       acceleration of horizontally averaged vertical velocity is zero.
    645           IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1 )  THEN
    646              IF ( j == 0 )  THEN
    647                 DO  k = 1, nz
    648                    ar(0,k) = 0.0_wp
    649                 ENDDO
    650              ENDIF
    651           ENDIF
    652 
    653        END SUBROUTINE substi_1dd
    654 
    655     END SUBROUTINE tridia_1dd
     554    ENDDO
     555
     556    DO  k = 0, nz-1
     557       DO  i = 0, nx
     558          a = -1.0_wp * ddzu_pres(k+2) * ddzw(k+1) * rho_air_zw(k+1)
     559          c = -1.0_wp * ddzu_pres(k+1) * ddzw(k+1) * rho_air_zw(k)
     560          tri_for_1d(1,i,k) = a + c - l(i) * rho_air(k+1)
     561       ENDDO
     562    ENDDO
     563    IF ( ibc_p_b == 1 )  THEN
     564       DO  i = 0, nx
     565          tri_for_1d(1,i,0) = tri_for_1d(1,i,0) + tri_for_1d(2,i,0)
     566       ENDDO
     567    ENDIF
     568    IF ( ibc_p_t == 1 )  THEN
     569       DO  i = 0, nx
     570          tri_for_1d(1,i,nz-1) = tri_for_1d(1,i,nz-1) + tri_for_1d(3,i,nz-1)
     571       ENDDO
     572    ENDIF
     573
     574 END SUBROUTINE maketri_1dd
     575
     576
     577!--------------------------------------------------------------------------------------------------!
     578! Description:
     579! ------------
     580!> Splitting of the tridiagonal matrix (Thomas algorithm).
     581!--------------------------------------------------------------------------------------------------!
     582 SUBROUTINE split_1dd
     583
     584    IMPLICIT NONE
     585
     586    INTEGER(iwp) ::  i  !<
     587    INTEGER(iwp) ::  k  !<
     588
     589
     590!
     591!-- Splitting
     592    DO  i = 0, nx
     593       tri_for_1d(4,i,0) = tri_for_1d(1,i,0)
     594    ENDDO
     595    DO  k = 1, nz-1
     596       DO  i = 0, nx
     597          tri_for_1d(5,i,k) = tri_for_1d(2,i,k) / tri_for_1d(4,i,k-1)
     598          tri_for_1d(4,i,k) = tri_for_1d(1,i,k) - tri_for_1d(3,i,k-1) * tri_for_1d(5,i,k)
     599       ENDDO
     600    ENDDO
     601
     602 END SUBROUTINE split_1dd
     603
     604
     605!--------------------------------------------------------------------------------------------------!
     606! Description:
     607! ------------
     608!> Substitution (Forward and Backward) (Thomas algorithm).
     609!--------------------------------------------------------------------------------------------------!
     610 SUBROUTINE substi_1dd( ar, tri_for_1d )
     611
     612
     613    IMPLICIT NONE
     614
     615    INTEGER(iwp) ::  i  !<
     616    INTEGER(iwp) ::  k  !<
     617
     618    REAL(wp), DIMENSION(0:nx,nz)       ::  ar          !<
     619    REAL(wp), DIMENSION(0:nx,0:nz-1)   ::  ar1         !<
     620    REAL(wp), DIMENSION(5,0:nx,0:nz-1) ::  tri_for_1d  !<
     621
     622!
     623!-- Forward substitution
     624    DO  i = 0, nx
     625       ar1(i,0) = ar(i,1)
     626    ENDDO
     627    DO  k = 1, nz-1
     628       DO  i = 0, nx
     629          ar1(i,k) = ar(i,k+1) - tri_for_1d(5,i,k) * ar1(i,k-1)
     630       ENDDO
     631    ENDDO
     632
     633!
     634!-- Backward substitution
     635!-- Note, the add of 1.0E-20 in the denominator is due to avoid divisions by zero appearing if the
     636!-- pressure bc is set to neumann at the top of the model domain.
     637    DO  i = 0, nx
     638       ar(i,nz) = ar1(i,nz-1) / ( tri_for_1d(4,i,nz-1) + 1.0E-20_wp )
     639    ENDDO
     640    DO  k = nz-2, 0, -1
     641       DO  i = 0, nx
     642          ar(i,k+1) = ( ar1(i,k) - tri_for_1d(3,i,k) * ar(i,k+2) ) / tri_for_1d(4,i,k)
     643       ENDDO
     644    ENDDO
     645
     646!
     647!-- Indices i=0, j=0 correspond to horizontally averaged pressure. The respective values of ar
     648!-- should be zero at all k-levels if acceleration of horizontally averaged vertical velocity is
     649!-- zero.
     650    IF ( ibc_p_b == 1  .AND.  ibc_p_t == 1 )  THEN
     651       IF ( j == 0 )  THEN
     652          DO  k = 1, nz
     653             ar(0,k) = 0.0_wp
     654          ENDDO
     655       ENDIF
     656    ENDIF
     657
     658 END SUBROUTINE substi_1dd
     659
     660 END SUBROUTINE tridia_1dd
    656661
    657662

  • palm/trunk/SOURCE/turbulence_closure_mod.f90

    r4495 r4510  
    11!> @file turbulence_closure_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/>.
    16 !
    17 ! Copyright 2017-2020 Leibniz Universitaet Hannover
    18 !--------------------------------------------------------------------------------!
     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/>.
     15!
     16! Copyright 1997-2020 Leibniz Universitaet Hannover
     17!--------------------------------------------------------------------------------------------------!
     18!
    1919!
    2020! Current revisions:
    … …  
    2525! -----------------
    2626! $Id$
     27! file re-formatted to follow the PALM coding standard
     28!
     29! 4495 2020-04-13 20:11:20Z raasch
    2730! workaround for Intel14 compiler added
    28 ! 
     31!
    2932! 4486 2020-04-02 20:45:12Z maronga
    30 ! Bugfix: include topography in calculation of distance_to_wall (1.5-order-dai
    31 ! closure)
    32 !
     33! Bugfix: include topography in calculation of distance_to_wall (1.5-order-dai closure)
     34!
    3335! 4481 2020-03-31 18:55:54Z maronga
    34 ! - added new LES closure after Dai et al. (2020), which provides much better grid
    35 !   convergence in stable boundary layer runs. The implementation is experimental
    36 !   at the moment and should be used with special care. The new SGS closure can be
    37 !   switched on via turbulence_closure = '1.5-order-dai'
    38 ! - variable ml_wall_adjusted renamed to delta as it represents a grid size and
    39 !   not a mixing length (see Equations 14 and 18 in Maronga et al. 2015, GMD)
     36! - added new LES closure after Dai et al. (2020), which provides much better grid convergence in
     37!   stable boundary layer runs. The implementation is experimental at the moment and should be used
     38!   with special care. The new SGS closure can be switched on via
     39!   turbulence_closure = '1.5-order-dai'
     40! - variable ml_wall_adjusted renamed to delta as it represents a grid size and not a mixing length
     41!   (see Equations 14 and 18 in Maronga et al. 2015, GMD)
    4042! - nameing of turbulence closures revised:
    4143!        'Moeng_Wyngaard' to '1.5-order'
    … …  
    4446! - LOGICAL steering variable renamed:
    4547!        les_mw to les_default
    46 ! 
     48!
    4749! 4473 2020-03-25 21:04:07Z gronemeier
    4850! - rename l-grid to gridsize-geometric-mean
    … …  
    5456!          l to ml
    5557! - adjust some comments
    56 ! - corrected definition of wall-adjusted mixing length to include
    57 !   gridsize-geometric-mean
     58! - corrected definition of wall-adjusted mixing length to include gridsize-geometric-mean
    5859! - moved definition of wall_adjustment_factor to this module
    5960!
    … …  
    6869!
    6970! 4346 2019-12-18 11:55:56Z motisi
    70 ! Introduction of wall_flags_total_0, which currently sets bits based on static
    71 ! topography information used in wall_flags_static_0
     71! Introduction of wall_flags_total_0, which currently sets bits based on static topography
     72! information used in wall_flags_static_0
    7273!
    7374! 4329 2019-12-10 15:46:36Z motisi
    … …  
    8182!
    8283! 4170 2019-08-19 17:12:31Z gronemeier
    83 ! - add performance optimizations according to K. Ketelsen
    84 !   to diffusion_e and tcm_diffusivities_default
     84! - add performance optimizations according to K. Ketelsen to diffusion_e and
     85!   tcm_diffusivities_default
    8586! - bugfix in calculating l_wall for vertical walls
    8687! - bugfix in using l_wall in initialization (consider wall_adjustment_factor)
    … …  
    9192!
    9293! 4110 2019-07-22 17:05:21Z suehring
    93 ! pass integer flag array as well as boundary flags to WS scalar advection
    94 ! routine
     94! pass integer flag array as well as boundary flags to WS scalar advection routine
    9595!
    9696! 4109 2019-07-22 17:00:34Z suehring
    9797! - Modularize setting of boundary conditions for TKE and dissipation
    9898! - Neumann boundary condition for TKE at model top is set also in child domain
    99 ! - Revise setting of Neumann boundary conditions at non-cyclic lateral
    100 !   boundaries
    101 ! - Bugfix, set Neumann boundary condition for TKE at vertical wall instead of
    102 !   an implicit Dirichlet boundary condition which implied a sink of TKE
    103 !   at vertical walls
     99! - Revise setting of Neumann boundary conditions at non-cyclic lateral boundaries
     100! - Bugfix, set Neumann boundary condition for TKE at vertical wall instead of an implicit Dirichlet
     101!   boundary condition which implied a sink of TKE at vertical walls
    104102!
    105103! 4048 2019-06-21 21:00:21Z knoop
    … …  
    113111!
    114112! 3719 2019-02-06 13:10:18Z kanani
    115 ! Changed log_point to log_point_s, otherwise this overlaps with
    116 ! 'all progn.equations' cpu measurement.
     113! Changed log_point to log_point_s, otherwise this overlaps with 'all progn.equations' cpu
     114! measurement.
    117115!
    118116! 3684 2019-01-20 20:20:58Z knoop
    … …  
    136134!> @todo Check for random disturbances
    137135!> @note <Enter notes on the module>
    138 !-----------------------------------------------------------------------------!
     136!--------------------------------------------------------------------------------------------------!
    139137 MODULE turbulence_closure_mod
    140138
    141139
    142     USE arrays_3d,                                                            &
    143         ONLY:  diss, diss_1, diss_2, diss_3, diss_p, dzu, e, e_1, e_2, e_3,   &
    144                e_p, kh, km, mean_inflow_profiles, prho, pt, tdiss_m,          &
    145                te_m, tend, u, v, vpt, w
    146 
    147     USE basic_constants_and_equations_mod,                                    &
    148         ONLY:  g, kappa, lv_d_cp, lv_d_rd, rd_d_rv
    149 
    150     USE control_parameters,                                                   &
    151         ONLY:  bc_dirichlet_l,                                                &
    152                bc_dirichlet_n,                                                &
    153                bc_dirichlet_r,                                                &
    154                bc_dirichlet_s,                                                &
    155                bc_radiation_l,                                                &
    156                bc_radiation_n,                                                &
    157                bc_radiation_r,                                                &
    158                bc_radiation_s,                                                &
    159                child_domain,                                                  &
    160                constant_diffusion, dt_3d, e_init, humidity,                   &
    161                initializing_actions, intermediate_timestep_count,             &
    162                intermediate_timestep_count_max, km_constant,                  &
    163                les_dai, les_dynamic, les_default,                             &
    164                ocean_mode, plant_canopy, prandtl_number,                      &
    165                pt_reference, rans_mode, rans_tke_e, rans_tke_l,               &
    166                timestep_scheme, turbulence_closure,                           &
    167                turbulent_inflow, use_upstream_for_tke, vpt_reference,         &
    168                ws_scheme_sca, current_timestep_number
    169 
    170     USE advec_ws,                                                             &
     140    USE arrays_3d,                                                                                 &
     141        ONLY:  diss,                                                                               &
     142               diss_1,                                                                             &
     143               diss_2,                                                                             &
     144               diss_3,                                                                             &
     145               diss_p,                                                                             &
     146               dzu,                                                                                &
     147               e,                                                                                  &
     148               e_1,                                                                                &
     149               e_2,                                                                                &
     150               e_3,                                                                                &
     151               e_p,                                                                                &
     152               kh,                                                                                 &
     153               km,                                                                                 &
     154               mean_inflow_profiles,                                                               &
     155               prho,                                                                               &
     156               pt,                                                                                 &
     157               tdiss_m,                                                                            &
     158               te_m,                                                                               &
     159               tend,                                                                               &
     160               u,                                                                                  &
     161               v,                                                                                  &
     162               vpt,                                                                                &
     163               w
     164
     165    USE basic_constants_and_equations_mod,                                                         &
     166        ONLY:  g,                                                                                  &
     167               kappa,                                                                              &
     168               lv_d_cp,                                                                            &
     169               lv_d_rd,                                                                            &
     170               rd_d_rv
     171
     172    USE control_parameters,                                                                        &
     173        ONLY:  bc_dirichlet_l,                                                                     &
     174               bc_dirichlet_n,                                                                     &
     175               bc_dirichlet_r,                                                                     &
     176               bc_dirichlet_s,                                                                     &
     177               bc_radiation_l,                                                                     &
     178               bc_radiation_n,                                                                     &
     179               bc_radiation_r,                                                                     &
     180               bc_radiation_s,                                                                     &
     181               child_domain,                                                                       &
     182               constant_diffusion,                                                                 &
     183               current_timestep_number,                                                            &
     184               dt_3d,                                                                              &
     185               e_init,                                                                             &
     186               humidity,                                                                           &
     187               initializing_actions,                                                               &
     188               intermediate_timestep_count,                                                        &
     189               intermediate_timestep_count_max,                                                    &
     190               km_constant,                                                                        &
     191               les_dai,                                                                            &
     192               les_dynamic,                                                                        &
     193               les_default,                                                                        &
     194               ocean_mode,                                                                         &
     195               plant_canopy,                                                                       &
     196               prandtl_number,                                                                     &
     197               pt_reference,                                                                       &
     198               rans_mode,                                                                          &
     199               rans_tke_e,                                                                         &
     200               rans_tke_l,                                                                         &
     201               timestep_scheme,                                                                    &
     202               turbulence_closure,                                                                 &
     203               turbulent_inflow,                                                                   &
     204               use_upstream_for_tke,                                                               &
     205               vpt_reference,                                                                      &
     206               ws_scheme_sca
     207
     208
     209    USE advec_ws,                                                                                  &
    171210        ONLY:  advec_s_ws
    172211
    173     USE advec_s_bc_mod,                                                       &
     212    USE advec_s_bc_mod,                                                                            &
    174213        ONLY:  advec_s_bc
    175214
    176     USE advec_s_pw_mod,                                                       &
     215    USE advec_s_pw_mod,                                                                            &
    177216        ONLY:  advec_s_pw
    178217
    179     USE advec_s_up_mod,                                                       &
     218    USE advec_s_up_mod,                                                                            &
    180219        ONLY:  advec_s_up
    181220
    182     USE cpulog,                                                               &
     221    USE cpulog,                                                                                    &
    183222        ONLY:  cpu_log, log_point_s
    184223
    185     USE indices,                                                              &
    186         ONLY:  advc_flags_s,                                                  &
    187                nbgp, nxl, nxlg, nxr, nxrg, nyn, nyng, nys, nysg, nzb, nzt,    &
    188                topo_top_ind,                                                  &
     224    USE indices,                                                                                   &
     225        ONLY:  advc_flags_s,                                                                       &
     226               nbgp,                                                                               &
     227               nxl,                                                                                &
     228               nxlg,                                                                               &
     229               nxr,                                                                                &
     230               nxrg,                                                                               &
     231               nyn,                                                                                &
     232               nyng,                                                                               &
     233               nys,                                                                                &
     234               nysg,                                                                               &
     235               nzb,                                                                                &
     236               nzt,                                                                                &
     237               topo_top_ind,                                                                       &
    189238               wall_flags_total_0
    190239
    191240    USE kinds
    192241
    193     USE ocean_mod,                                                            &
     242    USE ocean_mod,                                                                                 &
    194243        ONLY:  prho_reference
    195244
    196245    USE pegrid
    197246
    198     USE plant_canopy_model_mod,                                               &
     247    USE plant_canopy_model_mod,                                                                    &
    199248        ONLY:  pcm_tendency
    200249
    201     USE statistics,                                                           &
    202         ONLY:  hom, hom_sum, statistic_regions
    203 
    204     USE surface_mod,                                                          &
    205         ONLY:  bc_h,                                                          &
    206                bc_v,                                                          &
    207                surf_def_h,                                                    &
    208                surf_def_v,                                                    &
    209                surf_lsm_h,                                                    &
    210                surf_lsm_v,                                                    &
    211                surf_usm_h,                                                    &
     250    USE statistics,                                                                                &
     251        ONLY:  hom,                                                                                &
     252               hom_sum,                                                                            &
     253               statistic_regions
     254
     255    USE surface_mod,                                                                               &
     256        ONLY:  bc_h,                                                                               &
     257               bc_v,                                                                               &
     258               surf_def_h,                                                                         &
     259               surf_def_v,                                                                         &
     260               surf_lsm_h,                                                                         &
     261               surf_lsm_v,                                                                         &
     262               surf_usm_h,                                                                         &
    212263               surf_usm_v
    213264
    … …  
    225276    REAL(wp) ::  wall_adjustment_factor = 1.8_wp  !< adjustment factor for mixing length
    226277
    227     REAL(wp), DIMENSION(0:4) :: rans_const_c = &        !< model constants for RANS mode (namelist param)
    228        (/ 0.55_wp, 1.44_wp, 1.92_wp, 1.44_wp, 0.0_wp /) !> default values fit for standard-tke-e closure
    229 
    230     REAL(wp), DIMENSION(2) :: rans_const_sigma = &      !< model constants for RANS mode, sigma values (namelist param)
    231        (/ 1.0_wp, 1.30_wp /)                            !> (sigma_e, sigma_diss)
    232 
    233     REAL(wp), ALLOCATABLE, DIMENSION(:)     ::  ml_blackadar             !< mixing length according to Blackadar
    234     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  delta                    !< grid size, possibly limited by wall adjustment factor
    235     REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  distance_to_wall         !< distance to the surface/wall
     278    REAL(wp), DIMENSION(0:4) ::  rans_const_c = &        !< model constants for RANS mode (namelist param)
     279       (/ 0.55_wp, 1.44_wp, 1.92_wp, 1.44_wp, 0.0_wp /)  !> default values fit for standard-tke-e closure
     280
     281    REAL(wp), DIMENSION(2) ::  rans_const_sigma = &      !< model constants for RANS mode, sigma values (namelist param)
     282       (/ 1.0_wp, 1.30_wp /)                             !> (sigma_e, sigma_diss)
     283
     284    REAL(wp), ALLOCATABLE, DIMENSION(:)     ::  ml_blackadar      !< mixing length according to Blackadar
     285    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  delta             !< grid size, possibly limited by wall adjustment factor
     286    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  distance_to_wall  !< distance to the surface/wall
    236287
    237288!
    238289!-- Public variables
    239     PUBLIC c_0, rans_const_c, rans_const_sigma
     290    PUBLIC c_0,                                                                                    &
     291           rans_const_c,                                                                           &
     292           rans_const_sigma
    240293
    241294    SAVE
    … …  
    244297!
    245298!-- Public subroutines
    246     PUBLIC                                                                     &
    247        tcm_boundary_conds,                                                     &
    248        tcm_check_parameters,                                                   &
    249        tcm_check_data_output,                                                  &
    250        tcm_define_netcdf_grid,                                                 &
    251        tcm_init_arrays,                                                        &
    252        tcm_init,                                                               &
    253        tcm_actions,                                                            &
    254        tcm_prognostic_equations,                                               &
    255        tcm_swap_timelevel,                                                     &
    256        tcm_3d_data_averaging,                                                  &
    257        tcm_data_output_2d,                                                     &
    258        tcm_data_output_3d,                                                     &
    259        tcm_diffusivities
     299    PUBLIC  tcm_actions,                                                                           &
     300            tcm_boundary_conds,                                                                    &
     301            tcm_check_parameters,                                                                  &
     302            tcm_check_data_output,                                                                 &
     303            tcm_data_output_2d,                                                                    &
     304            tcm_data_output_3d,                                                                    &
     305            tcm_define_netcdf_grid,                                                                &
     306            tcm_diffusivities,                                                                     &
     307            tcm_init_arrays,                                                                       &
     308            tcm_init,                                                                              &
     309            tcm_prognostic_equations,                                                              &
     310            tcm_swap_timelevel,                                                                    &
     311            tcm_3d_data_averaging
     312
     313
     314
    260315
    261316!
    … …  
    342397 CONTAINS
    343398
    344 !------------------------------------------------------------------------------!
     399!--------------------------------------------------------------------------------------------------!
    345400! Description:
    346401! ------------
    347402!> Check parameters routine for turbulence closure module.
    348 !------------------------------------------------------------------------------!
     403!--------------------------------------------------------------------------------------------------!
    349404 SUBROUTINE tcm_boundary_conds
    350405
    351     USE pmc_interface,                                                         &
     406    USE pmc_interface,                                                                             &
    352407        ONLY : rans_mode_parent
    353408
    … …  
    364419!
    365420!--    In LES mode, Neumann conditions with de/x_i=0 are assumed at solid walls.
    366 !--    Note, only TKE is prognostic in this case and dissipation is only
    367 !--    a diagnostic quantity.
     421!--    Note, only TKE is prognostic in this case and dissipation is only a diagnostic quantity.
    368422       IF ( .NOT. rans_mode )  THEN
    369423!
    … …  
    384438          DO  l = 0, 3
    385439!
    386 !--          Note concerning missing ACC directive for this loop: Even though
    387 !--          the data structure bc_v is present, it may not contain any
    388 !--          allocated arrays in the flat but also in a topography case,
    389 !--          leading to a runtime error. Therefore, omit ACC directives
    390 !--          for this loop, in contrast to the bc_h loop.
     440!--          Note concerning missing ACC directive for this loop: Even though the data structure
     441!--          bc_v is present, it may not contain any allocated arrays in the flat but also in a
     442!--          topography case, leading to a runtime error. Therefore, omit ACC directives for this
     443!--          loop, in contrast to the bc_h loop.
    391444             !$OMP PARALLEL DO PRIVATE( i, j, k )
    392445             DO  m = 1, bc_v(l)%ns
    … …  
    402455!
    403456!--       Use wall function within constant-flux layer
    404 !--       Note, grid points listed in bc_h are not included in any calculations in RANS mode and
    405 !--       are therefore not set here.
     457!--       Note, grid points listed in bc_h are not included in any calculations in RANS mode and are
     458!--       therefore not set here.
    406459!
    407460!--       Upward-facing surfaces
    … …  
    459512       ENDIF
    460513!
    461 !--    Set Neumann boundary condition for TKE at model top. Do this also
    462 !--    in case of a nested run.
     514!--    Set Neumann boundary condition for TKE at model top. Do this also in case of a nested run.
    463515       !$ACC KERNELS PRESENT(e_p)
    464516       e_p(nzt+1,:,:) = e_p(nzt,:,:)
    465517       !$ACC END KERNELS
    466518!
    467 !--    Nesting case: if parent operates in RANS mode and child in LES mode,
    468 !--    no TKE is transfered. This case, set Neumann conditions at lateral and
    469 !--    top child boundaries.
    470 !--    If not ( both either in RANS or in LES mode ), TKE boundary condition
    471 !--    is treated in the nesting.
     519!--    Nesting case: if parent operates in RANS mode and child in LES mode, no TKE is transfered.
     520!--    This case, set Neumann conditions at lateral and top child boundaries.
     521!--    If not ( both either in RANS or in LES mode ), TKE boundary condition is treated in the
     522!--    nesting.
    472523       If ( child_domain )  THEN
    473524          IF ( rans_mode_parent  .AND.  .NOT. rans_mode )  THEN
    … …  
    482533       ENDIF
    483534!
    484 !--    At in- and outflow boundaries also set Neumann boundary conditions
    485 !--    for the SGS-TKE. An exception is made for the child domain if
    486 !--    both parent and child operate in RANS mode. This case no
    487 !--    lateral Neumann boundary conditions will be set but Dirichlet
    488 !--    conditions will be set in the nesting.
    489        IF ( .NOT. child_domain  .AND.  .NOT. rans_mode_parent  .AND.           &
    490             .NOT. rans_mode )  THEN
     535!--    At in- and outflow boundaries also set Neumann boundary conditions for the SGS-TKE. An
     536!--    exception is made for the child domain if both parent and child operate in RANS mode. This
     537!--    case no lateral Neumann boundary conditions will be set but Dirichlet conditions will be set
     538!--    in the nesting.
     539       IF ( .NOT. child_domain  .AND.  .NOT. rans_mode_parent  .AND.  .NOT. rans_mode )  THEN
    491540          IF ( bc_dirichlet_s  .OR.  bc_radiation_s )  THEN
    492541             e_p(:,nys-1,:) = e_p(:,nys,:)
    … …  
    519568          j = surf_def_h(0)%j(m)
    520569          k = surf_def_h(0)%k(m)
    521           diss_p(k,j,i) = surf_def_h(0)%us(m)**3          &
    522                         / ( kappa * surf_def_h(0)%z_mo(m) )
     570          diss_p(k,j,i) = surf_def_h(0)%us(m)**3 / ( kappa * surf_def_h(0)%z_mo(m) )
    523571       ENDDO
    524572!
    … …  
    528576          j = surf_lsm_h%j(m)
    529577          k = surf_lsm_h%k(m)
    530           diss_p(k,j,i) = surf_lsm_h%us(m)**3          &
    531                         / ( kappa * surf_lsm_h%z_mo(m) )
     578          diss_p(k,j,i) = surf_lsm_h%us(m)**3 / ( kappa * surf_lsm_h%z_mo(m) )
    532579       ENDDO
    533580!
    … …  
    537584          j = surf_usm_h%j(m)
    538585          k = surf_usm_h%k(m)
    539           diss_p(k,j,i) = surf_usm_h%us(m)**3          &
    540                         / ( kappa * surf_usm_h%z_mo(m) )
     586          diss_p(k,j,i) = surf_usm_h%us(m)**3 / ( kappa * surf_usm_h%z_mo(m) )
    541587       ENDDO
    542588!
    … …  
    549595             j = surf_def_v(l)%j(m)
    550596             k = surf_def_v(l)%k(m)
    551              diss_p(k,j,i) = surf_def_v(l)%us(m)**3          &
    552                            / ( kappa * surf_def_v(l)%z_mo(m) )
     597             diss_p(k,j,i) = surf_def_v(l)%us(m)**3 / ( kappa * surf_def_v(l)%z_mo(m) )
    553598          ENDDO
    554599!
    … …  
    558603             j = surf_lsm_v(l)%j(m)
    559604             k = surf_lsm_v(l)%k(m)
    560              diss_p(k,j,i) = surf_lsm_v(l)%us(m)**3          &
    561                            / ( kappa * surf_lsm_v(l)%z_mo(m) )
     605             diss_p(k,j,i) = surf_lsm_v(l)%us(m)**3 / ( kappa * surf_lsm_v(l)%z_mo(m) )
    562606          ENDDO
    563607!
    … …  
    567611             j = surf_usm_v(l)%j(m)
    568612             k = surf_usm_v(l)%k(m)
    569              diss_p(k,j,i) = surf_usm_v(l)%us(m)**3          &
    570                            / ( kappa * surf_usm_v(l)%z_mo(m) )
     613             diss_p(k,j,i) = surf_usm_v(l)%us(m)**3 / ( kappa * surf_usm_v(l)%z_mo(m) )
    571614          ENDDO
    572615       ENDDO
    573616!
    574 !--    Limit change of diss to be between -90% and +100%. Also, set an absolute
    575 !--    minimum value
     617!--    Limit change of diss to be between -90% and +100%. Also, set an absolute minimum value
    576618       DO  i = nxl, nxr
    577619          DO  j = nys, nyn
    578620             DO  k = nzb, nzt+1
    579                 diss_p(k,j,i) = MAX( MIN( diss_p(k,j,i),          &
    580                                           2.0_wp * diss(k,j,i) ), &
    581                                      0.1_wp * diss(k,j,i),        &
    582                                      0.0001_wp )
     621                diss_p(k,j,i) = MAX( MIN( diss_p(k,j,i), 2.0_wp * diss(k,j,i) ),                   &
     622                                          0.1_wp * diss(k,j,i), 0.0001_wp )
    583623             ENDDO
    584624          ENDDO
    … …  
    591631 END SUBROUTINE tcm_boundary_conds
    592632
    593 !------------------------------------------------------------------------------!
     633!--------------------------------------------------------------------------------------------------!
    594634! Description:
    595635! ------------
    596636!> Check parameters routine for turbulence closure module.
    597 !------------------------------------------------------------------------------!
     637!--------------------------------------------------------------------------------------------------!
    598638 SUBROUTINE tcm_check_parameters
    599639
    600     USE control_parameters,                                                    &
    601         ONLY:  message_string, turbulent_inflow, turbulent_outflow
     640    USE control_parameters,                                                                        &
     641        ONLY:  message_string,                                                                     &
     642               turbulent_inflow,                                                                   &
     643               turbulent_outflow
    602644
    603645    IMPLICIT NONE
    … …  
    625667
    626668       CASE DEFAULT
    627           message_string = 'Unknown turbulence closure: ' //                &
    628                            TRIM( turbulence_closure )
     669          message_string = 'Unknown turbulence closure: ' // TRIM( turbulence_closure )
    629670          CALL message( 'tcm_check_parameters', 'PA0500', 1, 2, 0, 6, 0 )
    630671
    … …  
    641682       c_1 = rans_const_c(1)
    642683       c_2 = rans_const_c(2)
    643        c_3 = rans_const_c(3)   !> @todo clarify how to switch between different models
     684       c_3 = rans_const_c(3)   !> @todo Clarify how to switch between different models
    644685       c_4 = rans_const_c(4)
    645686
    646687       IF ( turbulent_inflow .OR. turbulent_outflow )  THEN
    647           message_string = 'turbulent inflow/outflow is not yet '//            &
    648                            'implemented for RANS mode'
     688          message_string = 'turbulent inflow/outflow is not yet '// 'implemented for RANS mode'
    649689          CALL message( 'tcm_check_parameters', 'PA0501', 1, 2, 0, 6, 0 )
    650690       ENDIF
    … …  
    653693!
    654694!--    LES mode
    655        c_0 = 0.1_wp    !according to Lilly (1967) and Deardorff (1980)
    656 
    657        dsig_e = 1.0_wp !assure to use K_m to calculate TKE instead
    658                        !of K_e which is used in RANS mode
     695       c_0 = 0.1_wp    !According to Lilly (1967) and Deardorff (1980)
     696
     697       dsig_e = 1.0_wp !Assure to use K_m to calculate TKE instead of K_e which is used in RANS mode
    659698
    660699    ENDIF
    … …  
    662701 END SUBROUTINE tcm_check_parameters
    663702
    664 !------------------------------------------------------------------------------!
     703!--------------------------------------------------------------------------------------------------!
    665704! Description:
    666705! ------------
    667706!> Check data output.
    668 !------------------------------------------------------------------------------!
     707!--------------------------------------------------------------------------------------------------!
    669708 SUBROUTINE tcm_check_data_output( var, unit )
    670709
    671710    IMPLICIT NONE
    672711
    673     CHARACTER (LEN=*) ::  unit     !< unit of output variable
    674     CHARACTER (LEN=*) ::  var      !< name of output variable
     712    CHARACTER(LEN=*) ::  unit  !< unit of output variable
     713    CHARACTER(LEN=*) ::  var   !< name of output variable
    675714
    676715
    … …  
    691730
    692731
    693 !------------------------------------------------------------------------------!
     732!--------------------------------------------------------------------------------------------------!
    694733! Description:
    695734! ------------
    696 !> Define appropriate grid for netcdf variables.
    697 !> It is called out from subroutine netcdf.
    698 !------------------------------------------------------------------------------!
     735!> Define appropriate grid for netcdf variables. It is called out from subroutine netcdf.
     736!--------------------------------------------------------------------------------------------------!
    699737 SUBROUTINE tcm_define_netcdf_grid( var, found, grid_x, grid_y, grid_z )
    700738
    701739    IMPLICIT NONE
    702740
    703     CHARACTER (LEN=*), INTENT(OUT) ::  grid_x   !< x grid of output variable
    704     CHARACTER (LEN=*), INTENT(OUT) ::  grid_y   !< y grid of output variable
    705     CHARACTER (LEN=*), INTENT(OUT) ::  grid_z   !< z grid of output variable
    706     CHARACTER (LEN=*), INTENT(IN)  ::  var      !< name of output variable
     741    CHARACTER(LEN=*), INTENT(OUT) ::  grid_x  !< x grid of output variable
     742    CHARACTER(LEN=*), INTENT(OUT) ::  grid_y  !< y grid of output variable
     743    CHARACTER(LEN=*), INTENT(OUT) ::  grid_z  !< z grid of output variable
     744    CHARACTER(LEN=*), INTENT(IN)  ::  var     !< name of output variable
    707745
    708746    LOGICAL, INTENT(OUT) ::  found   !< flag if output variable is found
    … …  
    740778
    741779
    742 !------------------------------------------------------------------------------!
     780!--------------------------------------------------------------------------------------------------!
    743781! Description:
    744782! ------------
    745783!> Average 3D data.
    746 !------------------------------------------------------------------------------!
     784!--------------------------------------------------------------------------------------------------!
    747785 SUBROUTINE tcm_3d_data_averaging( mode, variable )
    748786
    749787
    750     USE averaging,                                                             &
    751         ONLY:  diss_av, kh_av, km_av
    752 
    753     USE control_parameters,                                                    &
     788    USE averaging,                                                                                 &
     789        ONLY:  diss_av,                                                                            &
     790               kh_av,                                                                              &
     791               km_av
     792
     793    USE control_parameters,                                                                        &
    754794        ONLY:  average_count_3d
    755795
    756796    IMPLICIT NONE
    757797
    758     CHARACTER (LEN=*) ::  mode       !< flag defining mode 'allocate', 'sum' or 'average'
    759     CHARACTER (LEN=*) ::  variable   !< name of variable
    760 
    761     INTEGER(iwp) ::  i   !< loop index
    762     INTEGER(iwp) ::  j   !< loop index
    763     INTEGER(iwp) ::  k   !< loop index
     798    CHARACTER(LEN=*) ::  mode      !< flag defining mode 'allocate', 'sum' or 'average'
     799    CHARACTER(LEN=*) ::  variable  !< name of variable
     800
     801    INTEGER(iwp) ::  i  !< loop index
     802    INTEGER(iwp) ::  j  !< loop index
     803    INTEGER(iwp) ::  k  !< loop index
    764804
    765805    IF ( mode == 'allocate' )  THEN
    … …  
    841881                   DO  j = nysg, nyng
    842882                      DO  k = nzb, nzt+1
    843                          diss_av(k,j,i) = diss_av(k,j,i)                       &
    844                                         / REAL( average_count_3d, KIND=wp )
     883                         diss_av(k,j,i) = diss_av(k,j,i) / REAL( average_count_3d, KIND = wp )
    845884                      ENDDO
    846885                   ENDDO
    … …  
    853892                   DO  j = nysg, nyng
    854893                      DO  k = nzb, nzt+1
    855                          kh_av(k,j,i) = kh_av(k,j,i)                           &
    856                                         / REAL( average_count_3d, KIND=wp )
     894                         kh_av(k,j,i) = kh_av(k,j,i) / REAL( average_count_3d, KIND = wp )
    857895                      ENDDO
    858896                   ENDDO
    … …  
    865903                   DO  j = nysg, nyng
    866904                      DO  k = nzb, nzt+1
    867                          km_av(k,j,i) = km_av(k,j,i)                           &
    868                                         / REAL( average_count_3d, KIND=wp )
     905                         km_av(k,j,i) = km_av(k,j,i) / REAL( average_count_3d, KIND = wp )
    869906                      ENDDO
    870907                   ENDDO
    … …  
    879916
    880917
    881 !------------------------------------------------------------------------------!
     918!--------------------------------------------------------------------------------------------------!
    882919! Description:
    883920! ------------
    884921!> Define 2D output variables.
    885 !------------------------------------------------------------------------------!
    886  SUBROUTINE tcm_data_output_2d( av, variable, found, grid, mode, local_pf,     &
    887                                 nzb_do, nzt_do )
    888 
    889     USE averaging,                                                             &
    890         ONLY:  diss_av, kh_av, km_av
     922!--------------------------------------------------------------------------------------------------!
     923 SUBROUTINE tcm_data_output_2d( av, variable, found, grid, mode, local_pf, nzb_do, nzt_do )
     924
     925    USE averaging,                                                                                 &
     926        ONLY:  diss_av,                                                                            &
     927               kh_av,                                                                              &
     928               km_av
    891929
    892930    IMPLICIT NONE
    893931
    894     CHARACTER (LEN=*) ::  grid       !< name of vertical grid
    895     CHARACTER (LEN=*) ::  mode       !< either 'xy', 'xz' or 'yz'
    896     CHARACTER (LEN=*) ::  variable   !< name of variable
    897 
    898     INTEGER(iwp) ::  av        !< flag for (non-)average output
    899     INTEGER(iwp) ::  flag_nr   !< number of masking flag
    900     INTEGER(iwp) ::  i         !< loop index
    901     INTEGER(iwp) ::  j         !< loop index
    902     INTEGER(iwp) ::  k         !< loop index
    903     INTEGER(iwp) ::  nzb_do    !< vertical output index (bottom)
    904     INTEGER(iwp) ::  nzt_do    !< vertical output index (top)
     932    CHARACTER(LEN=*) ::  grid      !< name of vertical grid
     933    CHARACTER(LEN=*) ::  mode      !< either 'xy', 'xz' or 'yz'
     934    CHARACTER(LEN=*) ::  variable  !< name of variable
     935
     936    INTEGER(iwp) ::  av       !< flag for (non-)average output
     937    INTEGER(iwp) ::  flag_nr  !< number of masking flag
     938    INTEGER(iwp) ::  i        !< loop index
     939    INTEGER(iwp) ::  j        !< loop index
     940    INTEGER(iwp) ::  k        !< loop index
     941    INTEGER(iwp) ::  nzb_do   !< vertical output index (bottom)
     942    INTEGER(iwp) ::  nzt_do   !< vertical output index (top)
    905943
    906944    LOGICAL ::  found     !< flag if output variable is found
    … …  
    909947    REAL(wp) ::  fill_value = -9999.0_wp  !< value for the _FillValue attribute
    910948
    911     REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf !< local
    912        !< array to which output data is resorted to
     949    REAL(wp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf  !< local
     950     !< array to which output data is resorted to
    913951
    914952    REAL(wp), DIMENSION(:,:,:), POINTER ::  to_be_resorted  !< points to selected output variable
    … …  
    926964             to_be_resorted => diss
    927965          ELSE
    928              IF ( .NOT. ALLOCATED( diss_av ) ) THEN
     966             IF ( .NOT. ALLOCATED( diss_av ) )  THEN
    929967                ALLOCATE( diss_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    930968                diss_av = REAL( fill_value, KIND = wp )
    … …  
    938976             to_be_resorted => kh
    939977          ELSE
    940              IF ( .NOT. ALLOCATED( kh_av ) ) THEN
     978             IF ( .NOT. ALLOCATED( kh_av ) )  THEN
    941979                ALLOCATE( kh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    942980                kh_av = REAL( fill_value, KIND = wp )
    … …  
    950988             to_be_resorted => km
    951989          ELSE
    952              IF ( .NOT. ALLOCATED( km_av ) ) THEN
     990             IF ( .NOT. ALLOCATED( km_av ) )  THEN
    953991                ALLOCATE( km_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    954992                km_av = REAL( fill_value, KIND = wp )
    … …  
    9641002    END SELECT
    9651003
    966     IF ( found .AND. .NOT. resorted )  THEN
     1004    IF ( found  .AND.  .NOT. resorted )  THEN
    9671005       DO  i = nxl, nxr
    9681006          DO  j = nys, nyn
    9691007             DO  k = nzb_do, nzt_do
    970                 local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i),                &
    971                                   REAL( fill_value, KIND = wp ),               &
     1008                local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value, KIND = wp ),     &
    9721009                                  BTEST( wall_flags_total_0(k,j,i), flag_nr ) )
    9731010             ENDDO
    … …  
    9791016
    9801017
    981 !------------------------------------------------------------------------------!
     1018!--------------------------------------------------------------------------------------------------!
    9821019! Description:
    9831020! ------------
    9841021!> Define 3D output variables.
    985 !------------------------------------------------------------------------------!
     1022!--------------------------------------------------------------------------------------------------!
    9861023 SUBROUTINE tcm_data_output_3d( av, variable, found, local_pf, nzb_do, nzt_do )
    9871024
    9881025
    989     USE averaging,                                                             &
    990         ONLY:  diss_av, kh_av, km_av
     1026    USE averaging,                                                                                 &
     1027        ONLY:  diss_av,                                                                            &
     1028               kh_av,                                                                              &
     1029               km_av
    9911030
    9921031    IMPLICIT NONE
    9931032
    994     CHARACTER (LEN=*) ::  variable   !< name of variable
    995 
    996     INTEGER(iwp) ::  av        !< flag for (non-)average output
    997     INTEGER(iwp) ::  flag_nr   !< number of masking flag
    998     INTEGER(iwp) ::  i         !< loop index
    999     INTEGER(iwp) ::  j         !< loop index
    1000     INTEGER(iwp) ::  k         !< loop index
    1001     INTEGER(iwp) ::  nzb_do    !< lower limit of the data output (usually 0)
    1002     INTEGER(iwp) ::  nzt_do    !< vertical upper limit of the data output (usually nz_do3d)
     1033    CHARACTER(LEN=*) ::  variable  !< name of variable
     1034
     1035    INTEGER(iwp) ::  av       !< flag for (non-)average output
     1036    INTEGER(iwp) ::  flag_nr  !< number of masking flag
     1037    INTEGER(iwp) ::  i        !< loop index
     1038    INTEGER(iwp) ::  j        !< loop index
     1039    INTEGER(iwp) ::  k        !< loop index
     1040    INTEGER(iwp) ::  nzb_do   !< lower limit of the data output (usually 0)
     1041    INTEGER(iwp) ::  nzt_do   !< vertical upper limit of the data output (usually nz_do3d)
    10031042
    10041043    LOGICAL ::  found     !< flag if output variable is found
    … …  
    10071046    REAL(wp) ::  fill_value = -9999.0_wp  !< value for the _FillValue attribute
    10081047
    1009     REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf   !< local
     1048    REAL(sp), DIMENSION(nxl:nxr,nys:nyn,nzb_do:nzt_do) ::  local_pf  !< local
    10101049       !< array to which output data is resorted to
    10111050
    … …  
    10241063             to_be_resorted => diss
    10251064          ELSE
    1026              IF ( .NOT. ALLOCATED( diss_av ) ) THEN
     1065             IF ( .NOT. ALLOCATED( diss_av ) )  THEN
    10271066                ALLOCATE( diss_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    10281067                diss_av = REAL( fill_value, KIND = wp )
    … …  
    10351074             to_be_resorted => kh
    10361075          ELSE
    1037              IF ( .NOT. ALLOCATED( kh_av ) ) THEN
     1076             IF ( .NOT. ALLOCATED( kh_av ) )  THEN
    10381077                ALLOCATE( kh_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    10391078                kh_av = REAL( fill_value, KIND = wp )
    … …  
    10461085             to_be_resorted => km
    10471086          ELSE
    1048              IF ( .NOT. ALLOCATED( km_av ) ) THEN
     1087             IF ( .NOT. ALLOCATED( km_av ) )  THEN
    10491088                ALLOCATE( km_av(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    10501089                km_av = REAL( fill_value, KIND = wp )
    … …  
    10591098
    10601099
    1061     IF ( found .AND. .NOT. resorted )  THEN
     1100    IF ( found  .AND.  .NOT. resorted )  THEN
    10621101       DO  i = nxl, nxr
    10631102          DO  j = nys, nyn
    10641103             DO  k = nzb_do, nzt_do
    1065                 local_pf(i,j,k) = MERGE(                                 &
    1066                                    to_be_resorted(k,j,i),                &
    1067                                    REAL( fill_value, KIND = wp ),        &
    1068                                    BTEST( wall_flags_total_0(k,j,i), flag_nr ) )
     1104                local_pf(i,j,k) = MERGE( to_be_resorted(k,j,i), REAL( fill_value, KIND = wp ),     &
     1105                                  BTEST( wall_flags_total_0(k,j,i), flag_nr ) )
    10691106             ENDDO
    10701107          ENDDO
    … …  
    10761113
    10771114
    1078 !------------------------------------------------------------------------------!
     1115!--------------------------------------------------------------------------------------------------!
    10791116! Description:
    10801117! ------------
    10811118!> Allocate arrays and assign pointers.
    1082 !------------------------------------------------------------------------------!
     1119!--------------------------------------------------------------------------------------------------!
    10831120 SUBROUTINE tcm_init_arrays
    10841121
    1085     USE bulk_cloud_model_mod,                                                  &
     1122    USE bulk_cloud_model_mod,                                                                      &
    10861123        ONLY:  collision_turbulence
    10871124
    1088     USE pmc_interface,                                                         &
     1125    USE pmc_interface,                                                                             &
    10891126        ONLY:  nested_run
    10901127
    … …  
    11001137!
    11011138!-- Allocate arrays required for dissipation.
    1102 !-- Please note, if it is a nested run, arrays need to be allocated even if
    1103 !-- they do not necessarily need to be transferred, which is attributed to
    1104 !-- the design of the model coupler which allocates memory for each variable.
     1139!-- Please note, if it is a nested run, arrays need to be allocated even if they do not necessarily
     1140!-- need to be transferred, which is attributed to the design of the model coupler which allocates
     1141!-- memory for each variable.
    11051142    ALLOCATE( diss_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    11061143
    … …  
    11221159
    11231160
    1124 !------------------------------------------------------------------------------!
     1161!--------------------------------------------------------------------------------------------------!
    11251162! Description:
    11261163! ------------
    11271164!> Initialization of turbulence closure module.
    1128 !------------------------------------------------------------------------------!
     1165!--------------------------------------------------------------------------------------------------!
    11291166 SUBROUTINE tcm_init
    11301167
    1131     USE control_parameters,                                                    &
    1132         ONLY:  bc_dirichlet_l, complex_terrain, topography
    1133 
    1134     USE model_1d_mod,                                                          &
    1135         ONLY:  e1d, kh1d, km1d
     1168    USE control_parameters,                                                                        &
     1169        ONLY:  bc_dirichlet_l,                                                                     &
     1170               complex_terrain,                                                                    &
     1171               topography
     1172
     1173    USE model_1d_mod,                                                                              &
     1174        ONLY:  e1d,                                                                                &
     1175               kh1d,                                                                               &
     1176               km1d
    11361177
    11371178    IMPLICIT NONE
    11381179
    1139     INTEGER(iwp) :: i            !< loop index
    1140     INTEGER(iwp) :: j            !< loop index
    1141     INTEGER(iwp) :: k            !< loop index
    1142     INTEGER(iwp) :: nz_s_shift   !< lower shift index for scalars
    1143     INTEGER(iwp) :: nz_s_shift_l !< local lower shift index in case of turbulent inflow
     1180    INTEGER(iwp) :: i             !< loop index
     1181    INTEGER(iwp) :: j             !< loop index
     1182    INTEGER(iwp) :: k             !< loop index
     1183    INTEGER(iwp) :: nz_s_shift    !< lower shift index for scalars
     1184    INTEGER(iwp) :: nz_s_shift_l  !< local lower shift index in case of turbulent inflow
    11441185
    11451186!
    … …  
    11491190!
    11501191!-- Actions for initial runs
    1151     IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.            &
     1192    IF ( TRIM( initializing_actions ) /= 'read_restart_data'  .AND.                                &
    11521193         TRIM( initializing_actions ) /= 'cyclic_fill' )  THEN
    11531194
    11541195       IF ( INDEX( initializing_actions, 'set_1d-model_profiles' ) /= 0 )  THEN
    11551196
    1156           IF ( .NOT. rans_tke_e ) THEN
     1197          IF ( .NOT. rans_tke_e )  THEN
    11571198!
    11581199!--          Transfer initial profiles to the arrays of the 3D model
    … …  
    11731214          ELSE
    11741215!
    1175 !--          In case of TKE-e closure in RANS mode, do not use e, diss, and km
    1176 !--          profiles from 1D model. Instead, initialize with constant profiles
     1216!--          In case of TKE-e closure in RANS mode, do not use e, diss, and km profiles from 1D
     1217!--          model. Instead, initialize with constant profiles
    11771218             IF ( constant_diffusion )  THEN
    11781219                km = km_constant
    … …  
    11831224                   DO  j = nysg, nyng
    11841225                      DO  k = nzb+1, nzt
    1185                          km(k,j,i) = c_0 * SQRT( e_init ) * MIN( delta(k,j,i), &
    1186                                                                  ml_blackadar(k) )
     1226                         km(k,j,i) = c_0 * SQRT( e_init ) * MIN( delta(k,j,i), ml_blackadar(k) )
    11871227                      ENDDO
    11881228                   ENDDO
    … …  
    11941234             ELSE
    11951235                IF ( .NOT. ocean_mode )  THEN
    1196                    kh   = 0.01_wp   ! there must exist an initial diffusion, because
    1197                    km   = 0.01_wp   ! otherwise no TKE would be produced by the
    1198                                     ! production terms, as long as not yet
    1199                                     ! e = (u*/cm)**2 at k=nzb+1
     1236                   kh   = 0.01_wp   ! There must exist an initial diffusion, because otherwise no
     1237                   km   = 0.01_wp   ! TKE would be produced by the production terms, as long as not
     1238                                    ! yet e = (u*/cm)**2 at k=nzb+1
    12001239                ELSE
    12011240                   kh   = 0.00001_wp
    … …  
    12171256          ENDIF
    12181257
    1219        ELSEIF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0 .OR. &
     1258       ELSEIF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0 .OR.                   &
    12201259                INDEX( initializing_actions, 'inifor' ) /= 0 )  THEN
    12211260
    … …  
    12381277          ELSE
    12391278             IF ( .NOT. ocean_mode )  THEN
    1240                 kh   = 0.01_wp   ! there must exist an initial diffusion, because
    1241                 km   = 0.01_wp   ! otherwise no TKE would be produced by the
    1242                                  ! production terms, as long as not yet
     1279                kh   = 0.01_wp   ! There must exist an initial diffusion, because otherwise no TKE
     1280                km   = 0.01_wp   ! would be produced by the production terms, as long as not yet
    12431281                                 ! e = (u*/cm)**2 at k=nzb+1
    12441282             ELSE
    … …  
    12771315       ENDIF
    12781316
    1279     ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data'  .OR.         &
    1280              TRIM( initializing_actions ) == 'cyclic_fill' )                   &
    1281     THEN
    1282 
    1283 !
    1284 !--    In case of complex terrain and cyclic fill method as initialization,
    1285 !--    shift initial data in the vertical direction for each point in the
    1286 !--    x-y-plane depending on local surface height
    1287        IF ( complex_terrain  .AND.                                             &
    1288             TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
     1317    ELSEIF ( TRIM( initializing_actions ) == 'read_restart_data'  .OR.                             &
     1318             TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
     1319
     1320!
     1321!--    In case of complex terrain and cyclic fill method as initialization, shift initial data in
     1322!--    the vertical direction for each point in the x-y-plane depending on local surface height
     1323       IF ( complex_terrain  .AND.  TRIM( initializing_actions ) == 'cyclic_fill' )  THEN
    12891324          DO  i = nxlg, nxrg
    12901325             DO  j = nysg, nyng
    … …  
    13111346!
    13121347!--    Initialization of the turbulence recycling method
    1313        IF ( TRIM( initializing_actions ) == 'cyclic_fill'  .AND.               &
    1314             turbulent_inflow )  THEN
     1348       IF ( TRIM( initializing_actions ) == 'cyclic_fill'  .AND.  turbulent_inflow )  THEN
    13151349          mean_inflow_profiles(:,5) = hom_sum(:,8,0)   ! e
    13161350!
    1317 !--       In case of complex terrain, determine vertical displacement at inflow
    1318 !--       boundary and adjust mean inflow profiles
     1351!--       In case of complex terrain, determine vertical displacement at inflow boundary and adjust
     1352!--       mean inflow profiles
    13191353          IF ( complex_terrain )  THEN
    1320              IF ( nxlg <= 0 .AND. nxrg >= 0 .AND.  &
    1321                   nysg <= 0 .AND. nyng >= 0        )  THEN
     1354             IF ( nxlg <= 0  .AND.  nxrg >= 0  .AND.  nysg <= 0  .AND.  nyng >= 0 )  THEN
    13221355                nz_s_shift_l = topo_top_ind(0,0,0)
    13231356             ELSE
    … …  
    13251358             ENDIF
    13261359#if defined( __parallel )
    1327              CALL MPI_ALLREDUCE(nz_s_shift_l, nz_s_shift, 1, MPI_INTEGER,      &
    1328                                 MPI_MAX, comm2d, ierr)
     1360             CALL MPI_ALLREDUCE(nz_s_shift_l, nz_s_shift, 1, MPI_INTEGER, MPI_MAX, comm2d, ierr)
    13291361#else
    13301362             nz_s_shift = nz_s_shift_l
    13311363#endif
    1332              mean_inflow_profiles(nz_s_shift:nzt+1,5) =  &
    1333                 hom_sum(0:nzt+1-nz_s_shift,8,0)  ! e
     1364             mean_inflow_profiles(nz_s_shift:nzt+1,5) = hom_sum(0:nzt+1-nz_s_shift,8,0)  ! e
    13341365          ENDIF
    13351366!
    1336 !--       Use these mean profiles at the inflow (provided that Dirichlet
    1337 !--       conditions are used)
     1367!--       Use these mean profiles at the inflow (provided that Dirichlet conditions are used)
    13381368          IF ( bc_dirichlet_l )  THEN
    13391369             DO  j = nysg, nyng
    … …  
    13461376!
    13471377!--    Inside buildings set TKE back to zero
    1348        IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND.                &
    1349             topography /= 'flat' )  THEN
     1378       IF ( TRIM( initializing_actions ) == 'cyclic_fill' .AND. topography /= 'flat' )  THEN
    13501379!
    13511380!--       Inside buildings set TKE back to zero.
    1352 !--       Other scalars (km, kh,...) are ignored at present,
    1353 !--       maybe revise later.
     1381!--       Other scalars (km, kh,...) are ignored at present, maybe revise later.
    13541382          DO  i = nxlg, nxrg
    13551383             DO  j = nysg, nyng
    13561384                DO  k = nzb, nzt
    1357                    e(k,j,i)     = MERGE( e(k,j,i), 0.0_wp,                     &
    1358                                          BTEST( wall_flags_total_0(k,j,i), 0 ) )
     1385                   e(k,j,i) = MERGE( e(k,j,i), 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
    13591386                ENDDO
    13601387             ENDDO
    … …  
    13651392                DO  j = nysg, nyng
    13661393                   DO  k = nzb, nzt
    1367                       diss(k,j,i)    = MERGE( diss(k,j,i), 0.0_wp,             &
    1368                                          BTEST( wall_flags_total_0(k,j,i), 0 ) )
     1394                      diss(k,j,i) = MERGE( diss(k,j,i), 0.0_wp,                                    &
     1395                                    BTEST( wall_flags_total_0(k,j,i), 0 ) )
    13691396                   ENDDO
    13701397                ENDDO
    … …  
    13731400       ENDIF
    13741401!
    1375 !--    Initialize new time levels (only done in order to set boundary values
    1376 !--    including ghost points)
     1402!--    Initialize new time levels (only done in order to set boundary values including ghost points)
    13771403       e_p = e
    13781404!
    1379 !--    Allthough tendency arrays are set in prognostic_equations, they have
    1380 !--    to be predefined here because there they are used (but multiplied with 0)
    1381 !--    before they are set.
     1405!--    Although tendency arrays are set in prognostic_equations, they have to be predefined here
     1406!--    because there they are used (but multiplied with 0) before they are set.
    13821407       te_m = 0.0_wp
    13831408
    … …  
    13921417
    13931418
    1394 !------------------------------------------------------------------------------!
     1419!--------------------------------------------------------------------------------------------------!
    13951420! Description:
    13961421! ------------
    13971422!> Pre-computation of grid-dependent and near-wall mixing length.
    1398 !> @todo consider walls in horizontal direction at a distance further than a
    1399 !>       single grid point (RANS mode)
    1400 !------------------------------------------------------------------------------!
     1423!> @todo consider walls in horizontal direction at a distance further than a single grid point
     1424!> (RANS mode)
     1425!--------------------------------------------------------------------------------------------------!
    14011426 SUBROUTINE tcm_init_mixing_length
    14021427
    1403     USE arrays_3d,                                                             &
    1404         ONLY:  dzw, ug, vg, zu, zw
    1405 
    1406     USE control_parameters,                                                    &
    1407         ONLY:  f, message_string, wall_adjustment
    1408 
    1409     USE exchange_horiz_mod,                                                    &
    1410         ONLY:  exchange_horiz, exchange_horiz_int
    1411 
    1412     USE grid_variables,                                                        &
    1413         ONLY:  dx, dy
    1414 
    1415     USE indices,                                                               &
    1416         ONLY:  nbgp, nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb,  &
    1417                nzt, wall_flags_total_0
     1428    USE arrays_3d,                                                                                 &
     1429        ONLY:  dzw,                                                                                &
     1430               ug,                                                                                 &
     1431               vg,                                                                                 &
     1432               zu,                                                                                 &
     1433               zw
     1434
     1435    USE control_parameters,                                                                        &
     1436        ONLY:  f,                                                                                  &
     1437               message_string,                                                                     &
     1438               wall_adjustment
     1439
     1440    USE exchange_horiz_mod,                                                                        &
     1441        ONLY:  exchange_horiz,                                                                     &
     1442               exchange_horiz_int
     1443
     1444    USE grid_variables,                                                                            &
     1445        ONLY:  dx,                                                                                 &
     1446               dy
     1447
     1448    USE indices,                                                                                   &
     1449        ONLY:  nbgp,                                                                               &
     1450               nx,                                                                                 &
     1451               nxl,                                                                                &
     1452               nxlg,                                                                               &
     1453               nxr,                                                                                &
     1454               nxrg,                                                                               &
     1455               ny,                                                                                 &
     1456               nyn,                                                                                &
     1457               nyng,                                                                               &
     1458               nys,                                                                                &
     1459               nysg,                                                                               &
     1460               nzb,                                                                                &
     1461               nzt,                                                                                &
     1462               wall_flags_total_0
    14181463
    14191464    USE kinds
    … …  
    14221467    IMPLICIT NONE
    14231468
    1424     INTEGER(iwp) :: dist_dx     !< found distance devided by dx
    1425     INTEGER(iwp) :: i           !< index variable along x
    1426     INTEGER(iwp) :: ii          !< index variable along x
    1427     INTEGER(iwp) :: j           !< index variable along y
    1428     INTEGER(iwp) :: k           !< index variable along z
    1429     INTEGER(iwp) :: k_max_topo  !< index of maximum topography height
    1430     INTEGER(iwp) :: kk          !< index variable along z
    1431     INTEGER(iwp) :: rad_i       !< search radius in grid points along x
    1432     INTEGER(iwp) :: rad_i_l     !< possible search radius to the left
    1433     INTEGER(iwp) :: rad_i_r     !< possible search radius to the right
    1434     INTEGER(iwp) :: rad_j       !< search radius in grid points along y
    1435     INTEGER(iwp) :: rad_j_n     !< possible search radius to north
    1436     INTEGER(iwp) :: rad_j_s     !< possible search radius to south
    1437     INTEGER(iwp) :: rad_k       !< search radius in grid points along z
    1438     INTEGER(iwp) :: rad_k_b     !< search radius in grid points along negative z
    1439     INTEGER(iwp) :: rad_k_t     !< search radius in grid points along positive z
    1440 
    1441     INTEGER(KIND=1), DIMENSION(:,:), ALLOCATABLE :: vic_yz !< contains a quarter of a single yz-slice of vicinity
    1442 
    1443     INTEGER(KIND=1), DIMENSION(:,:,:), ALLOCATABLE :: vicinity !< contains topography information of the vicinity of (i/j/k)
    1444 
    1445     REAL(wp) :: distance_up            !< distance of grid box center to its boundary in upper direction
    1446     REAL(wp) :: distance_down          !< distance of grid box center to its boundary in lower direction
    1447     REAL(wp) :: distance_ns            !< distance of grid box center to its boundary in y direction
    1448     REAL(wp) :: distance_lr            !< distance of grid box center to its boundary in x direction
    1449     REAL(wp) :: distance_edge_yz_down  !< distance of grid box center to its boundary along yz diagonal (down)
    1450     REAL(wp) :: distance_edge_yz_up    !< distance of grid box center to its boundary along yz diagonal (up)
    1451     REAL(wp) :: distance_edge_xz_down  !< distance of grid box center to its boundary along xz diagonal
    1452     REAL(wp) :: distance_edge_xz_up    !< distance of grid box center to its boundary along xz diagonal (up)
    1453     REAL(wp) :: distance_edge_xy       !< distance of grid box center to its boundary along xy diagonal
    1454     REAL(wp) :: distance_corners_down  !< distance of grid box center to its upper corners
    1455     REAL(wp) :: distance_corners_up    !< distance of grid box center to its lower corners
    1456     REAL(wp) :: radius                 !< search radius in meter
     1469    INTEGER(iwp) ::  dist_dx     !< found distance devided by dx
     1470    INTEGER(iwp) ::  i           !< index variable along x
     1471    INTEGER(iwp) ::  ii          !< index variable along x
     1472    INTEGER(iwp) ::  j           !< index variable along y
     1473    INTEGER(iwp) ::  k           !< index variable along z
     1474    INTEGER(iwp) ::  k_max_topo  !< index of maximum topography height
     1475    INTEGER(iwp) ::  kk          !< index variable along z
     1476    INTEGER(iwp) ::  rad_i       !< search radius in grid points along x
     1477    INTEGER(iwp) ::  rad_i_l     !< possible search radius to the left
     1478    INTEGER(iwp) ::  rad_i_r     !< possible search radius to the right
     1479    INTEGER(iwp) ::  rad_j       !< search radius in grid points along y
     1480    INTEGER(iwp) ::  rad_j_n     !< possible search radius to north
     1481    INTEGER(iwp) ::  rad_j_s     !< possible search radius to south
     1482    INTEGER(iwp) ::  rad_k       !< search radius in grid points along z
     1483    INTEGER(iwp) ::  rad_k_b     !< search radius in grid points along negative z
     1484    INTEGER(iwp) ::  rad_k_t     !< search radius in grid points along positive z
     1485
     1486    INTEGER(KIND=1), DIMENSION(:,:), ALLOCATABLE ::  vic_yz !< contains a quarter of a single yz-slice of vicinity
     1487
     1488    INTEGER(KIND=1), DIMENSION(:,:,:), ALLOCATABLE ::  vicinity !< contains topography information of the vicinity of (i/j/k)
     1489
     1490    REAL(wp) ::  distance_up            !< distance of grid box center to its boundary in upper direction
     1491    REAL(wp) ::  distance_down          !< distance of grid box center to its boundary in lower direction
     1492    REAL(wp) ::  distance_ns            !< distance of grid box center to its boundary in y direction
     1493    REAL(wp) ::  distance_lr            !< distance of grid box center to its boundary in x direction
     1494    REAL(wp) ::  distance_edge_yz_down  !< distance of grid box center to its boundary along yz diagonal (down)
     1495    REAL(wp) ::  distance_edge_yz_up    !< distance of grid box center to its boundary along yz diagonal (up)
     1496    REAL(wp) ::  distance_edge_xz_down  !< distance of grid box center to its boundary along xz diagonal
     1497    REAL(wp) ::  distance_edge_xz_up    !< distance of grid box center to its boundary along xz diagonal (up)
     1498    REAL(wp) ::  distance_edge_xy       !< distance of grid box center to its boundary along xy diagonal
     1499    REAL(wp) ::  distance_corners_down  !< distance of grid box center to its upper corners
     1500    REAL(wp) ::  distance_corners_up    !< distance of grid box center to its lower corners
     1501    REAL(wp) ::  radius                 !< search radius in meter
    14571502
    14581503    REAL(wp), DIMENSION(nzb:nzt+1) ::  gridsize_geometric_mean  !< geometric mean of grid sizes dx, dy, dz
    … …  
    14761521       gridsize_geometric_mean(nzt+1) = gridsize_geometric_mean(nzt)
    14771522
    1478        IF ( ANY( gridsize_geometric_mean > 1.5_wp * dx * wall_adjustment_factor ) .OR. &
     1523       IF ( ANY( gridsize_geometric_mean > 1.5_wp * dx * wall_adjustment_factor )  .OR.            &
    14791524            ANY( gridsize_geometric_mean > 1.5_wp * dy * wall_adjustment_factor ) )  THEN
    1480             WRITE( message_string, * ) 'grid anisotropy exceeds threshold',    &
    1481                                       ' &starting from height level k = ', k, &
    1482                                       '.'
     1525            WRITE( message_string, * ) 'grid anisotropy exceeds threshold',                        &
     1526                                       ' &starting from height level k = ', k, '.'
    14831527           CALL message( 'init_grid', 'PA0202', 0, 1, 0, 6, 0 )
    14841528       ENDIF
    … …  
    14911535
    14921536!
    1493 !--    If Dai et al. (2020) closure is used, the distance to the wall (distance to nearest upward facing surface)
    1494 !--    must be stored
     1537!--    If Dai et al. (2020) closure is used, the distance to the wall (distance to nearest upward
     1538!--    facing surface) must be stored
    14951539       IF ( les_dai )  THEN
    14961540          DO  i = nxl, nxr
    … …  
    15051549       IF ( wall_adjustment )  THEN
    15061550!
    1507 !--       In case of topography, adjust mixing length if there is any wall at
    1508 !--       the surrounding grid boxes:
     1551!--       In case of topography, adjust mixing length if there is any wall at the surrounding grid
     1552!--       boxes:
    15091553          !> @todo check if this is correct also for the ocean case
    15101554          DO  i = nxl, nxr
    … …  
    15151559                   IF ( BTEST( wall_flags_total_0(k,j,i), 0 ) )  THEN
    15161560!
    1517 !--                   First, check if grid points directly next to current grid point
    1518 !--                   are surface grid points
     1561!--                   First, check if grid points directly next to current grid point are surface
     1562!--                   grid points
    15191563!--                   Check along...
    15201564!--                   ...vertical direction, down
    … …  
    15331577!
    15341578!--                   ...y-direction (north/south)
    1535                       IF ( .NOT. BTEST( wall_flags_total_0(k,j-1,i), 0 )  .OR. &
     1579                      IF ( .NOT. BTEST( wall_flags_total_0(k,j-1,i), 0 )  .OR.                     &
    15361580                           .NOT. BTEST( wall_flags_total_0(k,j+1,i), 0 ) )  THEN
    15371581                         distance_ns = 0.5_wp * dy
    … …  
    15411585!
    15421586!--                   ...x-direction (left/right)
    1543                       IF ( .NOT. BTEST( wall_flags_total_0(k,j,i-1), 0 )  .OR. &
     1587                      IF ( .NOT. BTEST( wall_flags_total_0(k,j,i-1), 0 )  .OR.                     &
    15441588                           .NOT. BTEST( wall_flags_total_0(k,j,i+1), 0 ) )  THEN
    15451589                         distance_lr = 0.5_wp * dx
    … …  
    15501594!--                   Now, check the edges along...
    15511595!--                   ...yz-direction (vertical edges, down)
    1552                       IF ( .NOT. BTEST( wall_flags_total_0(k-1,j-1,i), 0 )  .OR.  &
     1596                      IF ( .NOT. BTEST( wall_flags_total_0(k-1,j-1,i), 0 )  .OR.                   &
    15531597                           .NOT. BTEST( wall_flags_total_0(k-1,j+1,i), 0 )  )  THEN
    15541598                         distance_edge_yz_down = SQRT( 0.25_wp * dy**2 + ( zu(k) - zw(k-1) )**2 )
    … …  
    15581602!
    15591603!--                   ...yz-direction (vertical edges, up)
    1560                       IF ( .NOT. BTEST( wall_flags_total_0(k+1,j-1,i), 0 )  .OR.  &
     1604                      IF ( .NOT. BTEST( wall_flags_total_0(k+1,j-1,i), 0 )  .OR.                   &
    15611605                           .NOT. BTEST( wall_flags_total_0(k+1,j+1,i), 0 )  )  THEN
    15621606                         distance_edge_yz_up = SQRT( 0.25_wp * dy**2 + ( zw(k) - zu(k) )**2 )
    … …  
    15661610!
    15671611!--                   ...xz-direction (vertical edges, down)
    1568                       IF ( .NOT. BTEST( wall_flags_total_0(k-1,j,i-1), 0 )  .OR.  &
     1612                      IF ( .NOT. BTEST( wall_flags_total_0(k-1,j,i-1), 0 )  .OR.                   &
    15691613                           .NOT. BTEST( wall_flags_total_0(k-1,j,i+1), 0 )  )  THEN
    15701614                         distance_edge_xz_down = SQRT( 0.25_wp * dx**2 + ( zu(k) - zw(k-1) )**2 )
    … …  
    15741618!
    15751619!--                   ...xz-direction (vertical edges, up)
    1576                       IF ( .NOT. BTEST( wall_flags_total_0(k+1,j,i-1), 0 )  .OR.  &
     1620                      IF ( .NOT. BTEST( wall_flags_total_0(k+1,j,i-1), 0 )  .OR.                   &
    15771621                           .NOT. BTEST( wall_flags_total_0(k+1,j,i+1), 0 )  )  THEN
    15781622                         distance_edge_xz_up = SQRT( 0.25_wp * dx**2 + ( zw(k) - zu(k) )**2 )
    … …  
    15821626!
    15831627!--                   ...xy-direction (horizontal edges)
    1584                       IF ( .NOT. BTEST( wall_flags_total_0(k,j-1,i-1), 0 )  .OR. &
    1585                            .NOT. BTEST( wall_flags_total_0(k,j+1,i-1), 0 )  .OR. &
    1586                            .NOT. BTEST( wall_flags_total_0(k,j-1,i+1), 0 )  .OR. &
     1628                      IF ( .NOT. BTEST( wall_flags_total_0(k,j-1,i-1), 0 )  .OR.                   &
     1629                           .NOT. BTEST( wall_flags_total_0(k,j+1,i-1), 0 )  .OR.                   &
     1630                           .NOT. BTEST( wall_flags_total_0(k,j-1,i+1), 0 )  .OR.                   &
    15871631                           .NOT. BTEST( wall_flags_total_0(k,j+1,i+1), 0 ) )  THEN
    15881632                         distance_edge_xy = SQRT( 0.25_wp * ( dx**2 + dy**2 ) )
    … …  
    15931637!--                   Now, check the corners...
    15941638!--                   ...lower four corners
    1595                       IF ( .NOT. BTEST( wall_flags_total_0(k-1,j-1,i-1), 0 )  .OR. &
    1596                            .NOT. BTEST( wall_flags_total_0(k-1,j+1,i-1), 0 )  .OR. &
    1597                            .NOT. BTEST( wall_flags_total_0(k-1,j-1,i+1), 0 )  .OR. &
     1639                      IF ( .NOT. BTEST( wall_flags_total_0(k-1,j-1,i-1), 0 )  .OR.                 &
     1640                           .NOT. BTEST( wall_flags_total_0(k-1,j+1,i-1), 0 )  .OR.                 &
     1641                           .NOT. BTEST( wall_flags_total_0(k-1,j-1,i+1), 0 )  .OR.                 &
    15981642                           .NOT. BTEST( wall_flags_total_0(k-1,j+1,i+1), 0 ) )  THEN
    1599                          distance_corners_down = SQRT( 0.25_wp * ( dx**2 + dy**2 ) &
     1643                         distance_corners_down = SQRT( 0.25_wp * ( dx**2 + dy**2 )                 &
    16001644                                                       + ( zu(k) - zw(k-1) )**2 )
    16011645                      ELSE
    … …  
    16041648!
    16051649!--                   ...upper four corners
    1606                       IF ( .NOT. BTEST( wall_flags_total_0(k+1,j-1,i-1), 0 )  .OR. &
    1607                            .NOT. BTEST( wall_flags_total_0(k+1,j+1,i-1), 0 )  .OR. &
    1608                            .NOT. BTEST( wall_flags_total_0(k+1,j-1,i+1), 0 )  .OR. &
     1650                      IF ( .NOT. BTEST( wall_flags_total_0(k+1,j-1,i-1), 0 )  .OR.                 &
     1651                           .NOT. BTEST( wall_flags_total_0(k+1,j+1,i-1), 0 )  .OR.                 &
     1652                           .NOT. BTEST( wall_flags_total_0(k+1,j-1,i+1), 0 )  .OR.                 &
    16091653                           .NOT. BTEST( wall_flags_total_0(k+1,j+1,i+1), 0 ) )  THEN
    1610                          distance_corners_up = SQRT( 0.25_wp * ( dx**2 + dy**2 )   &
     1654                         distance_corners_up = SQRT( 0.25_wp * ( dx**2 + dy**2 )                   &
    16111655                                                     + ( zw(k) - zu(k) )**2 )
    16121656                      ELSE
    … …  
    16151659
    16161660!
    1617 !--                   Calculate the minimum distance from the wall and store it
    1618 !--                   temporarily in the array delta
    1619                       delta(k,j,i) = MIN(                                      &
    1620                          distance_up, distance_down, distance_ns, distance_lr, &
    1621                          distance_edge_yz_down, distance_edge_yz_up,           &
    1622                          distance_edge_xz_down, distance_edge_xz_up,           &
    1623                          distance_edge_xy,                                     &
    1624                          distance_corners_down, distance_corners_up )
    1625 
    1626 !
    1627 !--                   If Dai et al. (2020) closure is used, the distance to the wall
    1628 !--                   must be permanently stored
     1661!--                   Calculate the minimum distance from the wall and store it temporarily in the
     1662!--                   array delta
     1663                      delta(k,j,i) = MIN( distance_up, distance_down, distance_ns, distance_lr,    &
     1664                                          distance_edge_yz_down, distance_edge_yz_up,              &
     1665                                          distance_edge_xz_down, distance_edge_xz_up,              &
     1666                                          distance_edge_xy, distance_corners_down,                 &
     1667                                          distance_corners_up )
     1668
     1669!
     1670!--                   If Dai et al. (2020) closure is used, the distance to the wall must be
     1671!--                   permanently stored
    16291672                      IF ( les_dai  .AND.  delta(k,j,i) /= 9999999.9_wp )  THEN
    16301673                         distance_to_wall(k,j,i) = delta(k,j,i)
    … …  
    16351678                      delta(k,j,i) = wall_adjustment_factor * delta(k,j,i)
    16361679
    1637                   ENDIF  !if grid point belongs to atmosphere
    1638 
    1639 
    1640 
    1641 !
    1642 !--               The grid size (delta) is defined as the the minimum of the distance to
    1643 !--               the nearest wall * 1.8 and the geometric mean grid size.
     1680                  ENDIF  ! If grid point belongs to atmosphere
     1681
     1682
     1683
     1684!
     1685!--               The grid size (delta) is defined as the the minimum of the distance to the nearest
     1686!--               wall * 1.8 and the geometric mean grid size.
    16441687                  delta(k,j,i) = MIN( delta(k,j,i), gridsize_geometric_mean(k) )
    16451688
    … …  
    16591702!--    Calculate mixing length according to Blackadar (1962)
    16601703       IF ( f /= 0.0_wp )  THEN
    1661           length_scale_max = 2.7E-4_wp * SQRT( ug(nzt+1)**2 + vg(nzt+1)**2 )   &
    1662                            / ABS( f ) + 1.0E-10_wp
     1704          length_scale_max = 2.7E-4_wp * SQRT( ug(nzt+1)**2 + vg(nzt+1)**2 ) / ABS( f ) + 1.0E-10_wp
    16631705       ELSE
    16641706          length_scale_max = 30.0_wp
    … …  
    16771719          DO  j = nysg, nyng
    16781720             DO  k = nzb+1, nzt-1
    1679                 IF ( .NOT. BTEST( wall_flags_total_0(k,j,i), 0 ) .AND.  &
    1680                      k > k_max_topo )  &
     1721                IF ( .NOT. BTEST( wall_flags_total_0(k,j,i), 0 )  .AND.  k > k_max_topo )          &
    16811722                   k_max_topo = k
    16821723             ENDDO
    … …  
    16871728       delta(nzt+1,:,:) = ml_blackadar(nzt+1)
    16881729!
    1689 !--    Limit mixing length to either nearest wall or Blackadar mixing length.
    1690 !--    For that, analyze each grid point (i/j/k) ("analysed grid point") and
    1691 !--    search within its vicinity for the shortest distance to a wall by cal-
    1692 !--    culating the distance between the analysed grid point and the "viewed
    1693 !--    grid point" if it contains a wall (belongs to topography).
     1730!--    Limit mixing length to either nearest wall or Blackadar mixing length. For that, analyze each
     1731!--    grid point (i/j/k) ("analysed grid point") and search within its vicinity for the shortest
     1732!--    distance to a wall by calculating the distance between the analysed grid point and the
     1733!--    "viewed grid point" if it contains a wall (belongs to topography).
    16941734       DO  k = nzb+1, nzt
    16951735
    … …  
    17201760          rad_k = MAX( rad_k_b, rad_k_t )
    17211761!
    1722 !--       When analysed grid point lies above maximum topography, set search
    1723 !--       radius to 0 if the distance between the analysed grid point and max
    1724 !--       topography height is larger than the maximum search radius
     1762!--       When analysed grid point lies above maximum topography, set search radius to 0 if the
     1763!--       distance between the analysed grid point and max topography height is larger than the
     1764!--       maximum search radius
    17251765          IF ( zu(k-rad_k_b) > zu(k_max_topo) )  rad_k_b = 0
    17261766!
    … …  
    17291769
    17301770             !> @note shape of vicinity is larger in z direction
    1731              !>   Shape of vicinity is two grid points larger than actual search
    1732              !>   radius in vertical direction. The first and last grid point is
    1733              !>   always set to 1 to asure correct detection of topography. See
    1734              !>   function "shortest_distance" for details.
     1771             !>   Shape of vicinity is two grid points larger than actual search radius in vertical
     1772             !>   direction. The first and last grid point is always set to 1 to asure correct
     1773             !>   detection of topography. See function "shortest_distance" for details.
    17351774             !>   2018-03-16, gronemeier
    17361775             ALLOCATE( vicinity(-rad_k-1:rad_k+1,-rad_j:rad_j,-rad_i:rad_i) )
    … …  
    17481787                      vicinity(-rad_k:rad_k,:,:) = 0
    17491788!
    1750 !--                   Copy area surrounding analysed grid point into vicinity.
    1751 !--                   First, limit size of data copied to vicinity by the domain
    1752 !--                   border
     1789!--                   Copy area surrounding analysed grid point into vicinity. First, limit size of
     1790!--                   data copied to vicinity by the domain border
    17531791                      !> @note limit copied area to 1 grid point in hor. dir.
    1754                       !>   Ignore walls in horizontal direction which are
    1755                       !>   further away than a single grid point. This allows to
    1756                       !>   only search within local subdomain without the need
    1757                       !>   of global topography information.
    1758                       !>   The error made by this assumption are acceptable at
    1759                       !>   the moment.
     1792                      !>   Ignore walls in horizontal direction which are further away than a single
     1793                      !>   grid point. This allows to only search within local subdomain without the
     1794                      !>   need of global topography information. The error made by this assumption
     1795                      !>   are acceptable at the moment.
    17601796                      !>   2018-10-01, gronemeier
    17611797                      rad_i_l = MIN( 1, rad_i, i )
    … …  
    17651801                      rad_j_n = MIN( 1, rad_j, ny-j )
    17661802
    1767                       CALL copy_into_vicinity( k, j, i,           &
    1768                                                -rad_k_b, rad_k_t, &
    1769                                                -rad_j_s, rad_j_n, &
     1803                      CALL copy_into_vicinity( k, j, i, -rad_k_b, rad_k_t, -rad_j_s, rad_j_n,      &
    17701804                                               -rad_i_l, rad_i_r  )
    1771                       !> @note in case of cyclic boundaries, those parts of the
    1772                       !>   topography which lies beyond the domain borders but
    1773                       !>   still within the search radius still needs to be
    1774                       !>   copied into 'vicinity'. As the effective search
    1775                       !>   radius is limited to 1 at the moment, no further
    1776                       !>   copying is needed. Old implementation (prior to
    1777                       !>   2018-10-01) had this covered but used a global array.
    1778                       !>   2018-10-01, gronemeier
     1805                      !> @note in case of cyclic boundaries, those parts of the topography which
     1806                      !> lies beyond the domain borders but still within the search radius still
     1807                      !> needs to be copied into 'vicinity'. As the effective search radius is
     1808                      !> limited to 1 at the moment, no further copying is needed. Old
     1809                      !> implementation (prior to 2018-10-01) had this covered but used a global
     1810                      !> array.
     1811                      !> 2018-10-01, gronemeier
    17791812
    17801813!
    … …  
    17861819                         DO  ii = 0, dist_dx
    17871820!
    1788 !--                         Search along vertical direction only if below
    1789 !--                         maximum topography
    1790                             IF ( rad_k_t > 0 ) THEN
     1821!--                         Search along vertical direction only if below maximum topography
     1822                            IF ( rad_k_t > 0 )  THEN
    17911823!
    17921824!--                            Search for walls within octant (+++)
    17931825                               vic_yz = vicinity(0:rad_k+1,0:rad_j,ii)
    1794                                delta(k,j,i) = MIN( delta(k,j,i),             &
    1795                                        shortest_distance( vic_yz, .TRUE., ii ) )
     1826                               delta(k,j,i) = MIN( delta(k,j,i),                                   &
     1827                                                   shortest_distance( vic_yz, .TRUE., ii ) )
    17961828!
    17971829!--                            Search for walls within octant (+-+)
    1798 !--                            Switch order of array so that the analysed grid
    1799 !--                            point is always located at (0/0) (required by
    1800 !--                            shortest_distance").
     1830!--                            Switch order of array so that the analysed grid point is always
     1831!--                            located at (0/0) (required by shortest_distance").
    18011832                               vic_yz = vicinity(0:rad_k+1,0:-rad_j:-1,ii)
    1802                                delta(k,j,i) = MIN( delta(k,j,i),             &
    1803                                        shortest_distance( vic_yz, .TRUE., ii ) )
     1833                               delta(k,j,i) = MIN( delta(k,j,i),                                   &
     1834                                                   shortest_distance( vic_yz, .TRUE., ii ) )
    18041835
    18051836                            ENDIF
    … …  
    18071838!--                         Search for walls within octant (+--)
    18081839                            vic_yz = vicinity(0:-rad_k-1:-1,0:-rad_j:-1,ii)
    1809                             delta(k,j,i) = MIN( delta(k,j,i),                &
    1810                                       shortest_distance( vic_yz, .FALSE., ii ) )
     1840                            delta(k,j,i) = MIN( delta(k,j,i),                                      &
     1841                                                shortest_distance( vic_yz, .FALSE., ii ) )
    18111842!
    18121843!--                         Search for walls within octant (++-)
    18131844                            vic_yz = vicinity(0:-rad_k-1:-1,0:rad_j,ii)
    1814                             delta(k,j,i) = MIN( delta(k,j,i),                &
    1815                                       shortest_distance( vic_yz, .FALSE., ii ) )
     1845                            delta(k,j,i) = MIN( delta(k,j,i),                                      &
     1846                                                shortest_distance( vic_yz, .FALSE., ii ) )
    18161847!
    18171848!--                         Reduce search along x by already found distance
    … …  
    18231854                         DO  ii = 0, -dist_dx, -1
    18241855!
    1825 !--                         Search along vertical direction only if below
    1826 !--                         maximum topography
    1827                             IF ( rad_k_t > 0 ) THEN
     1856!--                         Search along vertical direction only if below maximum topography
     1857                            IF ( rad_k_t > 0 )  THEN
    18281858!
    18291859!--                            Search for walls within octant (-++)
    18301860                               vic_yz = vicinity(0:rad_k+1,0:rad_j,ii)
    1831                                delta(k,j,i) = MIN( delta(k,j,i),             &
    1832                                       shortest_distance( vic_yz, .TRUE., -ii ) )
     1861                               delta(k,j,i) = MIN( delta(k,j,i),                                   &
     1862                                                   shortest_distance( vic_yz, .TRUE., -ii ) )
    18331863!
    18341864!--                            Search for walls within octant (--+)
    1835 !--                            Switch order of array so that the analysed grid
    1836 !--                            point is always located at (0/0) (required by
    1837 !--                            shortest_distance").
     1865!--                            Switch order of array so that the analysed grid point is always
     1866!--                            located at (0/0) (required by shortest_distance").
    18381867                               vic_yz = vicinity(0:rad_k+1,0:-rad_j:-1,ii)
    1839                                delta(k,j,i) = MIN( delta(k,j,i),             &
    1840                                       shortest_distance( vic_yz, .TRUE., -ii ) )
     1868                               delta(k,j,i) = MIN( delta(k,j,i),                                   &
     1869                                                   shortest_distance( vic_yz, .TRUE., -ii ) )
    18411870
    18421871                            ENDIF
    … …  
    18441873!--                         Search for walls within octant (---)
    18451874                            vic_yz = vicinity(0:-rad_k-1:-1,0:-rad_j:-1,ii)
    1846                             delta(k,j,i) = MIN( delta(k,j,i),                &
    1847                                      shortest_distance( vic_yz, .FALSE., -ii ) )
     1875                            delta(k,j,i) = MIN( delta(k,j,i),                                      &
     1876                                                shortest_distance( vic_yz, .FALSE., -ii ) )
    18481877!
    18491878!--                         Search for walls within octant (-+-)
    18501879                            vic_yz = vicinity(0:-rad_k-1:-1,0:rad_j,ii)
    1851                             delta(k,j,i) = MIN( delta(k,j,i),                &
    1852                                      shortest_distance( vic_yz, .FALSE., -ii ) )
     1880                            delta(k,j,i) = MIN( delta(k,j,i),                                      &
     1881                                                shortest_distance( vic_yz, .FALSE., -ii ) )
    18531882!
    18541883!--                         Reduce search along x by already found distance
    … …  
    18711900             DEALLOCATE( vic_yz )
    18721901
    1873           ENDIF  !check vertical size of vicinity
     1902          ENDIF  !Check vertical size of vicinity
    18741903
    18751904       ENDDO  !k loop
    … …  
    18921921
    18931922    CONTAINS
    1894 !------------------------------------------------------------------------------!
     1923!--------------------------------------------------------------------------------------------------!
    18951924! Description:
    18961925! ------------
    1897 !> Calculate the shortest distance between position (i/j/k)=(0/0/0) and
    1898 !> (pos_i/jj/kk), where (jj/kk) is the position of the maximum of 'array'
    1899 !> closest to the origin (0/0) of 'array'.
    1900 !------------------------------------------------------------------------------!
     1926!> Calculate the shortest distance between position (i/j/k)=(0/0/0) and (pos_i/jj/kk), where (jj/kk)
     1927!> is the position of the maximum of 'array' closest to the origin (0/0) of 'array'.
     1928!--------------------------------------------------------------------------------------------------!
    19011929    REAL(wp) FUNCTION shortest_distance( array, orientation, pos_i )
    19021930
    19031931       IMPLICIT NONE
    19041932
    1905        LOGICAL, INTENT(IN) :: orientation    !< flag if array represents an array oriented upwards (true) or downwards (false)
    1906 
    1907        INTEGER(iwp), INTENT(IN) :: pos_i     !< x position of the yz-plane 'array'
    1908 
    1909        INTEGER(iwp) :: a                     !< loop index
    1910        INTEGER(iwp) :: b                     !< loop index
    1911        INTEGER(iwp) :: jj                    !< loop index
    1912 
    1913        INTEGER(KIND=1) :: maximum            !< maximum of array along z dimension
    1914 
    1915        INTEGER(iwp), DIMENSION(0:rad_j) :: loc_k !< location of closest wall along vertical dimension
    1916 
    1917        INTEGER(KIND=1), DIMENSION(0:rad_k+1,0:rad_j), INTENT(IN) :: array !< array containing a yz-plane at position pos_i
    1918 
    1919 !
    1920 !--    Get coordinate of first maximum along vertical dimension
    1921 !--    at each y position of array (similar to function maxloc but more stable).
     1933       INTEGER(iwp), INTENT(IN) ::  pos_i  !< x position of the yz-plane 'array'
     1934
     1935       INTEGER(iwp) ::  a   !< loop index
     1936       INTEGER(iwp) ::  b   !< loop index
     1937       INTEGER(iwp) ::  jj  !< loop index
     1938
     1939       INTEGER(KIND=1) ::  maximum  !< maximum of array along z dimension
     1940
     1941       INTEGER(iwp), DIMENSION(0:rad_j) ::  loc_k  !< location of closest wall along vertical dimension
     1942
     1943       INTEGER(KIND=1), DIMENSION(0:rad_k+1,0:rad_j), INTENT(IN) ::  array  !< array containing a yz-plane at position pos_i
     1944
     1945       LOGICAL, INTENT(IN) ::  orientation  !< flag if array represents an array oriented upwards (true) or downwards (false)
     1946
     1947!
     1948!--    Get coordinate of first maximum along vertical dimension at each y position of array
     1949!--    (similar to function maxloc but more stable).
    19221950       DO  a = 0, rad_j
    19231951          loc_k(a) = rad_k+1
    … …  
    19341962       shortest_distance = radius
    19351963!
    1936 !--    Calculate distance between position (0/0/0) and
    1937 !--    position (pos_i/jj/loc(jj)) and only save the shortest distance.
    1938        IF ( orientation ) THEN  !if array is oriented upwards
     1964!--    Calculate distance between position (0/0/0) and position (pos_i/jj/loc(jj)) and only save the
     1965!--    shortest distance.
     1966       IF ( orientation )  THEN  !if array is oriented upwards
    19391967          DO  jj = 0, rad_j
    1940              shortest_distance =                                               &
    1941                 MIN( shortest_distance,                                        &
    1942                      SQRT( MAX(REAL(pos_i, KIND=wp)*dx-0.5_wp*dx, 0.0_wp)**2   &
    1943                          + MAX(REAL(jj, KIND=wp)*dy-0.5_wp*dy, 0.0_wp)**2      &
    1944                          + MAX(zw(loc_k(jj)+k-1)-zu(k), 0.0_wp)**2             &
    1945                          )                                                     &
    1946                    )
     1968             shortest_distance =                                                                   &
     1969               MIN( shortest_distance,                                                             &
     1970                    SQRT( MAX( REAL( pos_i, KIND = wp ) * dx - 0.5_wp * dx, 0.0_wp)**2             &
     1971                        + MAX( REAL( jj, KIND = wp ) * dy - 0.5_wp * dy, 0.0_wp)**2                &
     1972                        + MAX( zw(loc_k(jj)+k-1) - zu(k), 0.0_wp)**2                               &
     1973                        )                                                                          &
     1974                  )
    19471975          ENDDO
    19481976       ELSE  !if array is oriented downwards
    19491977          !> @note MAX within zw required to circumvent error at domain border
    1950           !>   At the domain border, if non-cyclic boundary is present, the
    1951           !>   index for zw could be -1, which will be errorneous (zw(-1) does
    1952           !>   not exist). The MAX function limits the index to be at least 0.
     1978          !>   At the domain border, if non-cyclic boundary is present, the index for zw could be
     1979          !>   -1, which will be errorneous (zw(-1) does not exist). The MAX function limits the
     1980          !>   index to be at least 0.
    19531981          DO  jj = 0, rad_j
    1954              shortest_distance =                                               &
    1955                 MIN( shortest_distance,                                        &
    1956                      SQRT( MAX(REAL(pos_i, KIND=wp)*dx-0.5_wp*dx, 0.0_wp)**2   &
    1957                          + MAX(REAL(jj, KIND=wp)*dy-0.5_wp*dy, 0.0_wp)**2      &
    1958                          + MAX(zu(k)-zw(MAX(k-loc_k(jj),0_iwp)), 0.0_wp)**2    &
    1959                          )                                                     &
     1982             shortest_distance =                                                                   &
     1983                MIN( shortest_distance,                                                            &
     1984                     SQRT( MAX(REAL(pos_i, KIND=wp)*dx-0.5_wp*dx, 0.0_wp)**2                       &
     1985                         + MAX(REAL(jj, KIND=wp)*dy-0.5_wp*dy, 0.0_wp)**2                          &
     1986                         + MAX(zu(k)-zw(MAX(k-loc_k(jj),0_iwp)), 0.0_wp)**2                        &
     1987                         )                                                                         &
    19601988                   )
    19611989          ENDDO
    … …  
    19641992    END FUNCTION
    19651993
    1966 !------------------------------------------------------------------------------!
     1994!--------------------------------------------------------------------------------------------------!
    19671995! Description:
    19681996! ------------
    1969 !> Copy a subarray of size (kb:kt,js:jn,il:ir) centered around grid point
    1970 !> (kp,jp,ip) containing the first bit of wall_flags_total_0 into the array
    1971 !> 'vicinity'. Only copy first bit as this indicates the presence of topography.
    1972 !------------------------------------------------------------------------------!
     1997!> Copy a subarray of size (kb:kt,js:jn,il:ir) centered around grid point (kp,jp,ip) containing the
     1998!> first bit of wall_flags_total_0 into the array 'vicinity'. Only copy first bit as this indicates
     1999!> the presence of topography.
     2000!--------------------------------------------------------------------------------------------------!
    19732001    SUBROUTINE copy_into_vicinity( kp, jp, ip, kb, kt, js, jn, il, ir )
    19742002
    19752003       IMPLICIT NONE
    19762004
    1977        INTEGER(iwp), INTENT(IN) :: il !< left loop boundary
    1978        INTEGER(iwp), INTENT(IN) :: ip !< center position in x-direction
    1979        INTEGER(iwp), INTENT(IN) :: ir !< right loop boundary
    1980        INTEGER(iwp), INTENT(IN) :: jn !< northern loop boundary
    1981        INTEGER(iwp), INTENT(IN) :: jp !< center position in y-direction
    1982        INTEGER(iwp), INTENT(IN) :: js !< southern loop boundary
    1983        INTEGER(iwp), INTENT(IN) :: kb !< bottom loop boundary
    1984        INTEGER(iwp), INTENT(IN) :: kp !< center position in z-direction
    1985        INTEGER(iwp), INTENT(IN) :: kt !< top loop boundary
    1986 
    1987        INTEGER(iwp) :: i   !< loop index
    1988        INTEGER(iwp) :: j   !< loop index
    1989        INTEGER(iwp) :: k   !< loop index
     2005       INTEGER(iwp), INTENT(IN) ::  il !< left loop boundary
     2006       INTEGER(iwp), INTENT(IN) ::  ip !< center position in x-direction
     2007       INTEGER(iwp), INTENT(IN) ::  ir !< right loop boundary
     2008       INTEGER(iwp), INTENT(IN) ::  jn !< northern loop boundary
     2009       INTEGER(iwp), INTENT(IN) ::  jp !< center position in y-direction
     2010       INTEGER(iwp), INTENT(IN) ::  js !< southern loop boundary
     2011       INTEGER(iwp), INTENT(IN) ::  kb !< bottom loop boundary
     2012       INTEGER(iwp), INTENT(IN) ::  kp !< center position in z-direction
     2013       INTEGER(iwp), INTENT(IN) ::  kt !< top loop boundary
     2014
     2015       INTEGER(iwp) ::  i  !< loop index
     2016       INTEGER(iwp) ::  j  !< loop index
     2017       INTEGER(iwp) ::  k  !< loop index
    19902018
    19912019       DO  i = il, ir
    19922020          DO  j = js, jn
    19932021             DO  k = kb, kt
    1994                 vicinity(k,j,i) = MERGE( 0, 1,               &
    1995                        BTEST( wall_flags_total_0(kp+k,jp+j,ip+i), 0 ) )
     2022                vicinity(k,j,i) = MERGE( 0, 1, BTEST( wall_flags_total_0(kp+k,jp+j,ip+i), 0 ) )
    19962023             ENDDO
    19972024          ENDDO
    … …  
    20032030
    20042031
    2005 !------------------------------------------------------------------------------!
     2032!--------------------------------------------------------------------------------------------------!
    20062033! Description:
    20072034! ------------
    20082035!> Initialize virtual velocities used later in production_e.
    2009 !------------------------------------------------------------------------------!
     2036!--------------------------------------------------------------------------------------------------!
    20102037 SUBROUTINE production_e_init
    20112038
    2012     USE arrays_3d,                                                             &
    2013         ONLY:  drho_air_zw, zu
    2014 
    2015     USE control_parameters,                                                    &
     2039    USE arrays_3d,                                                                                 &
     2040        ONLY:  drho_air_zw,                                                                        &
     2041               zu
     2042
     2043    USE control_parameters,                                                                        &
    20162044        ONLY:  constant_flux_layer
    20172045
    2018     USE surface_layer_fluxes_mod,                                              &
     2046    USE surface_layer_fluxes_mod,                                                                  &
    20192047        ONLY:  phi_m
    20202048
    20212049    IMPLICIT NONE
    20222050
    2023     INTEGER(iwp) ::  i      !< grid index x-direction
    2024     INTEGER(iwp) ::  j      !< grid index y-direction
    2025     INTEGER(iwp) ::  k      !< grid index z-direction
    2026     INTEGER(iwp) ::  m      !< running index surface elements
    2027 
    2028     REAL(wp) ::  km_sfc     !< diffusion coefficient, used to compute virtual velocities
     2051    INTEGER(iwp) ::  i  !< grid index x-direction
     2052    INTEGER(iwp) ::  j  !< grid index y-direction
     2053    INTEGER(iwp) ::  k  !< grid index z-direction
     2054    INTEGER(iwp) ::  m  !< running index surface elements
     2055
     2056    REAL(wp) ::  km_sfc  !< diffusion coefficient, used to compute virtual velocities
    20292057
    20302058    IF ( constant_flux_layer )  THEN
    20312059!
    2032 !--    Calculate a virtual velocity at the surface in a way that the
    2033 !--    vertical velocity gradient at k = 1 (u(k+1)-u_0) matches the
    2034 !--    Prandtl law (-w'u'/km). This gradient is used in the TKE shear
    2035 !--    production term at k=1 (see production_e_ij).
    2036 !--    The velocity gradient has to be limited in case of too small km
    2037 !--    (otherwise the timestep may be significantly reduced by large
    2038 !--    surface winds).
    2039 !--    not available in case of non-cyclic boundary conditions.
     2060!--    Calculate a virtual velocity at the surface in a way that the vertical velocity gradient at
     2061!--    k = 1 (u(k+1)-u_0) matches the Prandtl law (-w'u'/km). This gradient is used in the TKE shear
     2062!--    production term at k=1 (see production_e_ij). The velocity gradient has to be limited in case
     2063!--    of too small km (otherwise the timestep may be significantly reduced by large surface winds).
     2064!--    Not available in case of non-cyclic boundary conditions.
    20402065!--    Default surfaces, upward-facing
    20412066       !$OMP PARALLEL DO PRIVATE(i,j,k,m)
    … …  
    20482073          k = surf_def_h(0)%k(m)
    20492074!
    2050 !--       Note, calculation of u_0 and v_0 is not fully accurate, as u/v
    2051 !--       and km are not on the same grid. Actually, a further
    2052 !--       interpolation of km onto the u/v-grid is necessary. However, the
     2075!--       Note, calculation of u_0 and v_0 is not fully accurate, as u/v and km are not on the same
     2076!--       grid. Actually, a further interpolation of km onto the u/v-grid is necessary. However, the
    20532077!--       effect of this error is negligible.
    2054           km_sfc = kappa * surf_def_h(0)%us(m) * surf_def_h(0)%z_mo(m) /       &
     2078          km_sfc = kappa * surf_def_h(0)%us(m) * surf_def_h(0)%z_mo(m) /                           &
    20552079                   phi_m( surf_def_h(0)%z_mo(m) / surf_def_h(0)%ol(m) )
    20562080
    2057           surf_def_h(0)%u_0(m) = u(k+1,j,i) + surf_def_h(0)%usws(m) *          &
    2058                                      drho_air_zw(k-1)               *          &
    2059                                      ( zu(k+1) - zu(k-1)    )       /          &
    2060                                      ( km_sfc  + 1.0E-20_wp )
    2061           surf_def_h(0)%v_0(m) = v(k+1,j,i) + surf_def_h(0)%vsws(m) *          &
    2062                                      drho_air_zw(k-1)               *          &
    2063                                      ( zu(k+1) - zu(k-1)    )       /          &
    2064                                      ( km_sfc  + 1.0E-20_wp )
    2065 
    2066           IF ( ABS( u(k+1,j,i) - surf_def_h(0)%u_0(m) )  >                     &
    2067                ABS( u(k+1,j,i) - u(k-1,j,i)           )                        &
    2068              )  surf_def_h(0)%u_0(m) = u(k-1,j,i)
    2069 
    2070           IF ( ABS( v(k+1,j,i) - surf_def_h(0)%v_0(m) )  >                     &
    2071                ABS( v(k+1,j,i) - v(k-1,j,i)           )                        &
    2072              )  surf_def_h(0)%v_0(m) = v(k-1,j,i)
     2081          surf_def_h(0)%u_0(m) = u(k+1,j,i) + surf_def_h(0)%usws(m) * drho_air_zw(k-1) *           &
     2082                                 ( zu(k+1) - zu(k-1) ) / ( km_sfc  + 1.0E-20_wp )
     2083          surf_def_h(0)%v_0(m) = v(k+1,j,i) + surf_def_h(0)%vsws(m) * drho_air_zw(k-1) *           &
     2084                                 ( zu(k+1) - zu(k-1) ) / ( km_sfc  + 1.0E-20_wp )
     2085
     2086          IF ( ABS( u(k+1,j,i) - surf_def_h(0)%u_0(m) ) > ABS( u(k+1,j,i) - u(k-1,j,i) ) )         &
     2087             surf_def_h(0)%u_0(m) = u(k-1,j,i)
     2088
     2089          IF ( ABS( v(k+1,j,i) - surf_def_h(0)%v_0(m) ) > ABS( v(k+1,j,i) - v(k-1,j,i) ) )         &
     2090             surf_def_h(0)%v_0(m) = v(k-1,j,i)
    20732091
    20742092       ENDDO
    … …  
    20842102          k = surf_def_h(1)%k(m)
    20852103!
    2086 !--       Note, calculation of u_0 and v_0 is not fully accurate, as u/v
    2087 !--       and km are not on the same grid. Actually, a further
    2088 !--       interpolation of km onto the u/v-grid is necessary. However, the
     2104!--       Note, calculation of u_0 and v_0 is not fully accurate, as u/v and km are not on the same
     2105!--       grid. Actually, a further interpolation of km onto the u/v-grid is necessary. However, the
    20892106!--       effect of this error is negligible.
    2090           surf_def_h(1)%u_0(m) = u(k-1,j,i) - surf_def_h(1)%usws(m) *          &
    2091                                      drho_air_zw(k-1) *                        &
    2092                                      ( zu(k+1)    - zu(k-1)    )  /            &
    2093                                      ( km(k,j,i)  + 1.0E-20_wp )
    2094           surf_def_h(1)%v_0(m) = v(k-1,j,i) - surf_def_h(1)%vsws(m) *          &
    2095                                      drho_air_zw(k-1) *                        &
    2096                                      ( zu(k+1)    - zu(k-1)    )  /            &
    2097                                      ( km(k,j,i)  + 1.0E-20_wp )
    2098 
    2099           IF ( ABS( surf_def_h(1)%u_0(m) - u(k-1,j,i) )  >                     &
    2100                ABS( u(k+1,j,i)           - u(k-1,j,i) )                        &
    2101              )  surf_def_h(1)%u_0(m) = u(k+1,j,i)
    2102 
    2103           IF ( ABS( surf_def_h(1)%v_0(m) - v(k-1,j,i) )  >                     &
    2104                ABS( v(k+1,j,i)           - v(k-1,j,i) )                        &
    2105              )  surf_def_h(1)%v_0(m) = v(k+1,j,i)
     2107          surf_def_h(1)%u_0(m) = u(k-1,j,i) - surf_def_h(1)%usws(m) * drho_air_zw(k-1) *           &
     2108                                     ( zu(k+1) - zu(k-1) ) / ( km(k,j,i) + 1.0E-20_wp )
     2109          surf_def_h(1)%v_0(m) = v(k-1,j,i) - surf_def_h(1)%vsws(m) * drho_air_zw(k-1) *           &
     2110                                     ( zu(k+1) - zu(k-1) ) / (km(k,j,i) + 1.0E-20_wp )
     2111
     2112          IF ( ABS( surf_def_h(1)%u_0(m) - u(k-1,j,i) ) > ABS( u(k+1,j,i) - u(k-1,j,i) ) )         &
     2113             surf_def_h(1)%u_0(m) = u(k+1,j,i)
     2114
     2115          IF ( ABS( surf_def_h(1)%v_0(m) - v(k-1,j,i) ) > ABS( v(k+1,j,i) - v(k-1,j,i) ) )         &
     2116             surf_def_h(1)%v_0(m) = v(k+1,j,i)
    21062117
    21072118       ENDDO
    … …  
    21172128          k = surf_lsm_h%k(m)
    21182129!
    2119 !--       Note, calculation of u_0 and v_0 is not fully accurate, as u/v
    2120 !--       and km are not on the same grid. Actually, a further
    2121 !--       interpolation of km onto the u/v-grid is necessary. However, the
     2130!--       Note, calculation of u_0 and v_0 is not fully accurate, as u/v and km are not on the same
     2131!--       grid. Actually, a further interpolation of km onto the u/v-grid is necessary. However, the
    21222132!--       effect of this error is negligible.
    2123           km_sfc = kappa * surf_lsm_h%us(m) * surf_lsm_h%z_mo(m) /             &
     2133          km_sfc = kappa * surf_lsm_h%us(m) * surf_lsm_h%z_mo(m) /                                 &
    21242134                   phi_m( surf_lsm_h%z_mo(m) / surf_lsm_h%ol(m) )
    21252135
    2126           surf_lsm_h%u_0(m) = u(k+1,j,i) + surf_lsm_h%usws(m)    *             &
    2127                                         drho_air_zw(k-1)         *             &
    2128                                         ( zu(k+1) - zu(k-1)    ) /             &
    2129                                         ( km_sfc  + 1.0E-20_wp )
    2130           surf_lsm_h%v_0(m) = v(k+1,j,i) + surf_lsm_h%vsws(m)    *             &
    2131                                         drho_air_zw(k-1)         *             &
    2132                                         ( zu(k+1) - zu(k-1)    ) /             &
    2133                                         ( km_sfc  + 1.0E-20_wp )
    2134 
    2135           IF ( ABS( u(k+1,j,i) - surf_lsm_h%u_0(m) )  >                        &
    2136                ABS( u(k+1,j,i) - u(k-1,j,i)   )                                &
    2137              )  surf_lsm_h%u_0(m) = u(k-1,j,i)
    2138 
    2139           IF ( ABS( v(k+1,j,i) - surf_lsm_h%v_0(m) )  >                        &
    2140                ABS( v(k+1,j,i) - v(k-1,j,i)   )                                &
    2141              )  surf_lsm_h%v_0(m) = v(k-1,j,i)
     2136          surf_lsm_h%u_0(m) = u(k+1,j,i) + surf_lsm_h%usws(m) * drho_air_zw(k-1) *                 &
     2137                              ( zu(k+1) - zu(k-1) ) / ( km_sfc  + 1.0E-20_wp )
     2138          surf_lsm_h%v_0(m) = v(k+1,j,i) + surf_lsm_h%vsws(m) * drho_air_zw(k-1) *                 &
     2139                                        ( zu(k+1) - zu(k-1)) / ( km_sfc  + 1.0E-20_wp )
     2140
     2141          IF ( ABS( u(k+1,j,i) - surf_lsm_h%u_0(m) ) > ABS( u(k+1,j,i) - u(k-1,j,i) ) )            &
     2142             surf_lsm_h%u_0(m) = u(k-1,j,i)
     2143
     2144          IF ( ABS( v(k+1,j,i) - surf_lsm_h%v_0(m) ) > ABS( v(k+1,j,i) - v(k-1,j,i) ) )            &
     2145             surf_lsm_h%v_0(m) = v(k-1,j,i)
    21422146
    21432147       ENDDO
    … …  
    21532157          k = surf_usm_h%k(m)
    21542158!
    2155 !--       Note, calculation of u_0 and v_0 is not fully accurate, as u/v
    2156 !--       and km are not on the same grid. Actually, a further
    2157 !--       interpolation of km onto the u/v-grid is necessary. However, the
     2159!--       Note, calculation of u_0 and v_0 is not fully accurate, as u/v and km are not on the same
     2160!--       grid. Actually, a further interpolation of km onto the u/v-grid is necessary. However, the
    21582161!--       effect of this error is negligible.
    2159           km_sfc = kappa * surf_usm_h%us(m) * surf_usm_h%z_mo(m) /             &
     2162          km_sfc = kappa * surf_usm_h%us(m) * surf_usm_h%z_mo(m) /                                 &
    21602163                   phi_m( surf_usm_h%z_mo(m) / surf_usm_h%ol(m) )
    21612164
    2162           surf_usm_h%u_0(m) = u(k+1,j,i) + surf_usm_h%usws(m)    *             &
    2163                                         drho_air_zw(k-1)         *             &
    2164                                         ( zu(k+1) - zu(k-1)    ) /             &
    2165                                         ( km_sfc  + 1.0E-20_wp )
    2166           surf_usm_h%v_0(m) = v(k+1,j,i) + surf_usm_h%vsws(m)    *             &
    2167                                         drho_air_zw(k-1)         *             &
    2168                                         ( zu(k+1) - zu(k-1)    ) /             &
    2169                                         ( km_sfc  + 1.0E-20_wp )
    2170 
    2171           IF ( ABS( u(k+1,j,i) - surf_usm_h%u_0(m) )  >                        &
    2172                ABS( u(k+1,j,i) - u(k-1,j,i)   )                                &
    2173              )  surf_usm_h%u_0(m) = u(k-1,j,i)
    2174 
    2175           IF ( ABS( v(k+1,j,i) - surf_usm_h%v_0(m) )  >                        &
    2176                ABS( v(k+1,j,i) - v(k-1,j,i)   )                                &
    2177              )  surf_usm_h%v_0(m) = v(k-1,j,i)
     2165          surf_usm_h%u_0(m) = u(k+1,j,i) + surf_usm_h%usws(m) * drho_air_zw(k-1) *                 &
     2166                                        ( zu(k+1) - zu(k-1) ) / ( km_sfc  + 1.0E-20_wp )
     2167          surf_usm_h%v_0(m) = v(k+1,j,i) + surf_usm_h%vsws(m) * drho_air_zw(k-1) *                 &
     2168                                        ( zu(k+1) - zu(k-1) ) / ( km_sfc  + 1.0E-20_wp )
     2169
     2170          IF ( ABS( u(k+1,j,i) - surf_usm_h%u_0(m) ) > ABS( u(k+1,j,i) - u(k-1,j,i) ) )            &
     2171             surf_usm_h%u_0(m) = u(k-1,j,i)
     2172
     2173          IF ( ABS( v(k+1,j,i) - surf_usm_h%v_0(m) ) > ABS( v(k+1,j,i) - v(k-1,j,i) ) )            &
     2174             surf_usm_h%v_0(m) = v(k-1,j,i)
    21782175
    21792176       ENDDO
    … …  
    21922189
    21932190
    2194     CHARACTER (LEN=*) ::  location !<
     2191    CHARACTER(LEN=*) ::  location !<
    21952192
    21962193!    INTEGER(iwp) ::  i !<
    … …  
    22002197!
    22012198!-- Here the module-specific actions follow
    2202 !-- No calls for single grid points are allowed at locations before and
    2203 !-- after the timestep, since these calls are not within an i,j-loop
     2199!-- No calls for single grid points are allowed at locations before and after the timestep, since
     2200!-- these calls are not within an i, j-loop
    22042201    SELECT CASE ( location )
    22052202
    … …  
    22582255
    22592256
    2260     CHARACTER (LEN=*) ::  location
    2261 
    2262     INTEGER(iwp) ::  i
    2263     INTEGER(iwp) ::  j
     2257    CHARACTER(LEN=*) ::  location  !<
     2258
     2259    INTEGER(iwp) ::  i  !<
     2260    INTEGER(iwp) ::  j  !<
    22642261
    22652262!
    … …  
    23012298
    23022299
    2303 !------------------------------------------------------------------------------!
     2300!--------------------------------------------------------------------------------------------------!
    23042301! Description:
    23052302! ------------
    2306 !> Prognostic equation for subgrid-scale TKE and TKE dissipation rate.
    2307 !> Vector-optimized version
    2308 !------------------------------------------------------------------------------!
     2303!> Prognostic equation for subgrid-scale TKE and TKE dissipation rate. Vector-optimized version
     2304!--------------------------------------------------------------------------------------------------!
    23092305 SUBROUTINE tcm_prognostic_equations
    23102306
    2311     USE control_parameters,                                                    &
    2312         ONLY:  scalar_advec, tsc
     2307    USE control_parameters,                                                                        &
     2308        ONLY:  scalar_advec,                                                                       &
     2309               tsc
    23132310
    23142311    IMPLICIT NONE
    23152312
    2316     INTEGER(iwp) ::  i       !< loop index
    2317     INTEGER(iwp) ::  j       !< loop index
    2318     INTEGER(iwp) ::  k       !< loop index
    2319 
    2320     REAL(wp)     ::  sbt     !< wheighting factor for sub-time step
    2321 
    2322 !
    2323 !-- If required, compute prognostic equation for turbulent kinetic
    2324 !-- energy (TKE)
     2313    INTEGER(iwp) ::  i  !< loop index
     2314    INTEGER(iwp) ::  j  !< loop index
     2315    INTEGER(iwp) ::  k  !< loop index
     2316
     2317    REAL(wp)     ::  sbt  !< wheighting factor for sub-time step
     2318
     2319!
     2320!-- If required, compute prognostic equation for turbulent kinetic energy (TKE)
    23252321    IF ( .NOT. constant_diffusion )  THEN
    23262322
    … …  
    23542350             IF ( timestep_scheme(1:5) == 'runge' )  THEN
    23552351                IF ( ws_scheme_sca )  THEN
    2356                    CALL advec_s_ws( advc_flags_s, e, 'e',                      &
    2357                                     bc_dirichlet_l  .OR.  bc_radiation_l,      &
    2358                                     bc_dirichlet_n  .OR.  bc_radiation_n,      &
    2359                                     bc_dirichlet_r  .OR.  bc_radiation_r,      &
    2360                                     bc_dirichlet_s  .OR.  bc_radiation_s )
     2352                   CALL advec_s_ws( advc_flags_s, e, 'e', bc_dirichlet_l  .OR.  bc_radiation_l,    &
     2353                                                          bc_dirichlet_n  .OR.  bc_radiation_n,    &
     2354                                                          bc_dirichlet_r  .OR.  bc_radiation_r,    &
     2355                                                          bc_dirichlet_s  .OR.  bc_radiation_s )
    23612356                ELSE
    23622357                   CALL advec_s_pw( e )
    … …  
    23892384!
    23902385!--    Prognostic equation for TKE.
    2391 !--    Eliminate negative TKE values, which can occur due to numerical
    2392 !--    reasons in the course of the integration. In such cases the old TKE
    2393 !--    value is reduced by 90%.
     2386!--    Eliminate negative TKE values, which can occur due to numerical reasons in the course of the
     2387!--    integration. In such cases the old TKE value is reduced by 90%.
    23942388       !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
    23952389       !$ACC PRESENT(e, tend, te_m, wall_flags_total_0) &
    … …  
    23982392       DO  i = nxl, nxr
    23992393          DO  j = nys, nyn
    2400              !following directive is required to vectorize on Intel19
     2394             !Following directive is required to vectorize on Intel19
    24012395             !DIR$ IVDEP
    24022396             DO  k = nzb+1, nzt
    2403                 e_p(k,j,i) = e(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +        &
    2404                                                  tsc(3) * te_m(k,j,i) )        &
    2405                                         )                                      &
    2406                                    * MERGE( 1.0_wp, 0.0_wp,                    &
    2407                                        BTEST( wall_flags_total_0(k,j,i), 0 )   &
    2408                                           )
     2397                e_p(k,j,i) = e(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + tsc(3) * te_m(k,j,i) ) )   &
     2398                             * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
    24092399                IF ( e_p(k,j,i) < 0.0_wp )  e_p(k,j,i) = 0.1_wp * e(k,j,i)
    24102400             ENDDO
    … …  
    24252415                ENDDO
    24262416             ENDDO
    2427           ELSEIF ( intermediate_timestep_count < &
    2428                    intermediate_timestep_count_max )  THEN
     2417          ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max )  THEN
    24292418             !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
    24302419             !$ACC PRESENT(tend, te_m)
    … …  
    24322421                DO  j = nys, nyn
    24332422                   DO  k = nzb+1, nzt
    2434                       te_m(k,j,i) =   -9.5625_wp * tend(k,j,i)                 &
    2435                                      + 5.3125_wp * te_m(k,j,i)
     2423                      te_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * te_m(k,j,i)
    24362424                   ENDDO
    24372425                ENDDO
    … …  
    24662454
    24672455!
    2468 !--    dissipation-tendency terms with no communication
     2456!--    Dissipation-tendency terms with no communication
    24692457       IF ( scalar_advec /= 'bc-scheme'  .OR.  use_upstream_for_tke )  THEN
    24702458          IF ( use_upstream_for_tke )  THEN
    … …  
    24752463             IF ( timestep_scheme(1:5) == 'runge' )  THEN
    24762464                IF ( ws_scheme_sca )  THEN
    2477                    CALL advec_s_ws( advc_flags_s, diss, 'diss',                &
    2478                                     bc_dirichlet_l  .OR.  bc_radiation_l,      &
    2479                                     bc_dirichlet_n  .OR.  bc_radiation_n,      &
    2480                                     bc_dirichlet_r  .OR.  bc_radiation_r,      &
     2465                   CALL advec_s_ws( advc_flags_s, diss, 'diss',                                    &
     2466                                    bc_dirichlet_l  .OR.  bc_radiation_l,                          &
     2467                                    bc_dirichlet_n  .OR.  bc_radiation_n,                          &
     2468                                    bc_dirichlet_r  .OR.  bc_radiation_r,                          &
    24812469                                    bc_dirichlet_s  .OR.  bc_radiation_s )
    24822470                ELSE
    … …  
    25022490!
    25032491!--    Prognostic equation for TKE dissipation.
    2504 !--    Eliminate negative dissipation values, which can occur due to numerical
    2505 !--    reasons in the course of the integration. In such cases the old
    2506 !--    dissipation value is reduced by 90%.
     2492!--    Eliminate negative dissipation values, which can occur due to numerical reasons in the course
     2493!--    of the integration. In such cases the old dissipation value is reduced by 90%.
    25072494       DO  i = nxl, nxr
    25082495          DO  j = nys, nyn
    25092496             DO  k = nzb+1, nzt
    2510                 diss_p(k,j,i) = diss(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) +  &
    2511                                                  tsc(3) * tdiss_m(k,j,i) )     &
    2512                                         )                                      &
    2513                                    * MERGE( 1.0_wp, 0.0_wp,                    &
    2514                                        BTEST( wall_flags_total_0(k,j,i), 0 )   &
    2515                                           )
    2516                 IF ( diss_p(k,j,i) < 0.0_wp )                                  &
    2517                    diss_p(k,j,i) = 0.1_wp * diss(k,j,i)
     2497                diss_p(k,j,i) = diss(k,j,i) + ( dt_3d * ( sbt * tend(k,j,i) + tsc(3)               &
     2498                                                          * tdiss_m(k,j,i) ) )                     &
     2499                                * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
     2500                IF ( diss_p(k,j,i) < 0.0_wp )  diss_p(k,j,i) = 0.1_wp * diss(k,j,i)
    25182501             ENDDO
    25192502          ENDDO
    … …  
    25312514                ENDDO
    25322515             ENDDO
    2533           ELSEIF ( intermediate_timestep_count < &
    2534                    intermediate_timestep_count_max )  THEN
     2516          ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max )  THEN
    25352517             DO  i = nxl, nxr
    25362518                DO  j = nys, nyn
    25372519                   DO  k = nzb+1, nzt
    2538                       tdiss_m(k,j,i) =   -9.5625_wp * tend(k,j,i)              &
    2539                                         + 5.3125_wp * tdiss_m(k,j,i)
     2520                      tdiss_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tdiss_m(k,j,i)
    25402521                   ENDDO
    25412522                ENDDO
    … …  
    25512532
    25522533
    2553 !------------------------------------------------------------------------------!
     2534!--------------------------------------------------------------------------------------------------!
    25542535! Description:
    25552536! ------------
    2556 !> Prognostic equation for subgrid-scale TKE and TKE dissipation rate.
    2557 !> Cache-optimized version
    2558 !------------------------------------------------------------------------------!
     2537!> Prognostic equation for subgrid-scale TKE and TKE dissipation rate. Cache-optimized version
     2538!--------------------------------------------------------------------------------------------------!
    25592539 SUBROUTINE tcm_prognostic_equations_ij( i, j, i_omp, tn )
    25602540
    2561     USE arrays_3d,                                                             &
    2562         ONLY:  diss_l_diss, diss_l_e, diss_s_diss, diss_s_e, flux_l_diss,      &
    2563                flux_l_e, flux_s_diss, flux_s_e
    2564 
    2565     USE control_parameters,                                                    &
     2541    USE arrays_3d,                                                                                 &
     2542        ONLY:  diss_l_diss,                                                                        &
     2543               diss_l_e,                                                                           &
     2544               diss_s_diss,                                                                        &
     2545               diss_s_e,                                                                           &
     2546               flux_l_diss,                                                                        &
     2547               flux_l_e,                                                                           &
     2548               flux_s_diss,                                                                        &
     2549               flux_s_e
     2550
     2551    USE control_parameters,                                                                        &
    25662552        ONLY:  tsc
    25672553
    25682554    IMPLICIT NONE
    25692555
    2570     INTEGER(iwp) ::  i       !< loop index x direction
    2571     INTEGER(iwp) ::  i_omp   !< first loop index of i-loop in prognostic_equations
    2572     INTEGER(iwp) ::  j       !< loop index y direction
    2573     INTEGER(iwp) ::  k       !< loop index z direction
    2574     INTEGER(iwp) ::  tn      !< task number of openmp task
    2575 
    2576 !
    2577 !-- If required, compute prognostic equation for turbulent kinetic
    2578 !-- energy (TKE)
     2556    INTEGER(iwp) ::  i      !< loop index x direction
     2557    INTEGER(iwp) ::  i_omp  !< first loop index of i-loop in prognostic_equations
     2558    INTEGER(iwp) ::  j      !< loop index y direction
     2559    INTEGER(iwp) ::  k      !< loop index z direction
     2560    INTEGER(iwp) ::  tn     !< task number of openmp task
     2561
     2562!
     2563!-- If required, compute prognostic equation for turbulent kinetic energy (TKE)
    25792564    IF ( .NOT. constant_diffusion )  THEN
    25802565
    … …  
    25822567!--    Tendency-terms for TKE
    25832568       tend(:,j,i) = 0.0_wp
    2584        IF ( timestep_scheme(1:5) == 'runge'  &
    2585            .AND.  .NOT. use_upstream_for_tke )  THEN
     2569       IF ( timestep_scheme(1:5) == 'runge'  .AND.  .NOT. use_upstream_for_tke )  THEN
    25862570           IF ( ws_scheme_sca )  THEN
    2587                CALL advec_s_ws( advc_flags_s,                                  &
    2588                                 i, j, e, 'e', flux_s_e, diss_s_e,              &
    2589                                 flux_l_e, diss_l_e , i_omp, tn,                &
    2590                                 bc_dirichlet_l  .OR.  bc_radiation_l,          &
    2591                                 bc_dirichlet_n  .OR.  bc_radiation_n,          &
    2592                                 bc_dirichlet_r  .OR.  bc_radiation_r,          &
     2571               CALL advec_s_ws( advc_flags_s, i, j, e, 'e', flux_s_e, diss_s_e, flux_l_e,          &
     2572                                diss_l_e , i_omp, tn,                                              &
     2573                                bc_dirichlet_l  .OR.  bc_radiation_l,                              &
     2574                                bc_dirichlet_n  .OR.  bc_radiation_n,                              &
     2575                                bc_dirichlet_r  .OR.  bc_radiation_r,                              &
    25932576                                bc_dirichlet_s  .OR.  bc_radiation_s )
    25942577           ELSE
    … …  
    26192602!
    26202603!--    Prognostic equation for TKE.
    2621 !--    Eliminate negative TKE values, which can occur due to numerical
    2622 !--    reasons in the course of the integration. In such cases the old
    2623 !--    TKE value is reduced by 90%.
     2604!--    Eliminate negative TKE values, which can occur due to numerical reasons in the course of the
     2605!--    integration. In such cases the old TKE value is reduced by 90%.
    26242606       DO  k = nzb+1, nzt
    2625           e_p(k,j,i) = e(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) +           &
    2626                                               tsc(3) * te_m(k,j,i) )           &
    2627                                   )                                            &
    2628                                  * MERGE( 1.0_wp, 0.0_wp,                      &
    2629                                     BTEST( wall_flags_total_0(k,j,i), 0 )      &
    2630                                         )
     2607          e_p(k,j,i) = e(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) + tsc(3) * te_m(k,j,i) ) )      &
     2608                       * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
    26312609          IF ( e_p(k,j,i) <= 0.0_wp )  e_p(k,j,i) = 0.1_wp * e(k,j,i)
    26322610       ENDDO
    … …  
    26392617                te_m(k,j,i) = tend(k,j,i)
    26402618             ENDDO
    2641           ELSEIF ( intermediate_timestep_count < &
    2642                    intermediate_timestep_count_max )  THEN
     2619          ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max )  THEN
    26432620             DO  k = nzb+1, nzt
    2644                 te_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +                     &
    2645                                  5.3125_wp * te_m(k,j,i)
     2621                te_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * te_m(k,j,i)
    26462622             ENDDO
    26472623          ENDIF
    … …  
    26562632!--    Tendency-terms for dissipation
    26572633       tend(:,j,i) = 0.0_wp
    2658        IF ( timestep_scheme(1:5) == 'runge'  &
    2659            .AND.  .NOT. use_upstream_for_tke )  THEN
     2634       IF ( timestep_scheme(1:5) == 'runge'  .AND.  .NOT. use_upstream_for_tke )  THEN
    26602635           IF ( ws_scheme_sca )  THEN
    2661                CALL advec_s_ws( advc_flags_s,                                  &
    2662                                 i, j, diss, 'diss', flux_s_diss, diss_s_diss,  &
    2663                                 flux_l_diss, diss_l_diss, i_omp, tn,           &
    2664                                 bc_dirichlet_l  .OR.  bc_radiation_l,          &
    2665                                 bc_dirichlet_n  .OR.  bc_radiation_n,          &
    2666                                 bc_dirichlet_r  .OR.  bc_radiation_r,          &
     2636               CALL advec_s_ws( advc_flags_s, i, j, diss, 'diss', flux_s_diss, diss_s_diss,        &
     2637                                flux_l_diss, diss_l_diss, i_omp, tn,                               &
     2638                                bc_dirichlet_l  .OR.  bc_radiation_l,                              &
     2639                                bc_dirichlet_n  .OR.  bc_radiation_n,                              &
     2640                                bc_dirichlet_r  .OR.  bc_radiation_r,                              &
    26672641                                bc_dirichlet_s  .OR.  bc_radiation_s )
    26682642           ELSE
    … …  
    26862660!
    26872661!--    Prognostic equation for TKE dissipation
    2688 !--    Eliminate negative dissipation values, which can occur due to
    2689 !--    numerical reasons in the course of the integration. In such cases
    2690 !--    the old dissipation value is reduced by 90%.
     2662!--    Eliminate negative dissipation values, which can occur due to numerical reasons in the course
     2663!--    of the integration. In such cases the old dissipation value is reduced by 90%.
    26912664       DO  k = nzb+1, nzt
    2692           diss_p(k,j,i) = diss(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) +     &
    2693                                                     tsc(3) * tdiss_m(k,j,i) )  &
    2694                                         )                                      &
    2695                                         * MERGE( 1.0_wp, 0.0_wp,               &
    2696                                           BTEST( wall_flags_total_0(k,j,i), 0 )&
    2697                                                )
     2665          diss_p(k,j,i) = diss(k,j,i) + ( dt_3d * ( tsc(2) * tend(k,j,i) + tsc(3)                  &
     2666                          * tdiss_m(k,j,i) ) )                                                     &
     2667                          * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
    26982668       ENDDO
    26992669
    … …  
    27052675                tdiss_m(k,j,i) = tend(k,j,i)
    27062676             ENDDO
    2707           ELSEIF ( intermediate_timestep_count < &
    2708                    intermediate_timestep_count_max )  THEN
     2677          ELSEIF ( intermediate_timestep_count < intermediate_timestep_count_max )  THEN
    27092678             DO  k = nzb+1, nzt
    2710                 tdiss_m(k,j,i) =   -9.5625_wp * tend(k,j,i) +                  &
    2711                                     5.3125_wp * tdiss_m(k,j,i)
     2679                tdiss_m(k,j,i) = -9.5625_wp * tend(k,j,i) + 5.3125_wp * tdiss_m(k,j,i)
    27122680             ENDDO
    27132681          ENDIF
    27142682       ENDIF
    27152683
    2716     ENDIF   ! dissipation equation
     2684    ENDIF   ! Dissipation equation
    27172685
    27182686 END SUBROUTINE tcm_prognostic_equations_ij
    27192687
    27202688
    2721 !------------------------------------------------------------------------------!
     2689!--------------------------------------------------------------------------------------------------!
    27222690! Description:
    27232691! ------------
    27242692!> Production terms (shear + buoyancy) of the TKE.
    27252693!> Vector-optimized version
    2726 !> @warning The case with constant_flux_layer = F and use_surface_fluxes = T is
    2727 !>          not considered well!
    2728 !------------------------------------------------------------------------------!
     2694!> @warning The case with constant_flux_layer = F and use_surface_fluxes = T is not considered well!
     2695!--------------------------------------------------------------------------------------------------!
    27292696 SUBROUTINE production_e( diss_production )
    27302697
    2731     USE arrays_3d,                                                             &
    2732         ONLY:  ddzw, dd2zu, drho_air_zw, q, ql, d_exner, exner
    2733 
    2734     USE control_parameters,                                                    &
    2735         ONLY:  cloud_droplets, constant_flux_layer, neutral,                   &
    2736                rho_reference, use_single_reference_value, use_surface_fluxes,  &
     2698    USE arrays_3d,                                                                                 &
     2699        ONLY:  ddzw,                                                                               &
     2700               dd2zu,                                                                              &
     2701               d_exner,                                                                            &
     2702               drho_air_zw,                                                                        &
     2703               exner,                                                                              &
     2704               q,                                                                                  &
     2705               ql
     2706
     2707
     2708
     2709    USE control_parameters,                                                                        &
     2710        ONLY:  cloud_droplets,                                                                     &
     2711               constant_flux_layer,                                                                &
     2712               neutral,                                                                            &
     2713               rho_reference,                                                                      &
     2714               use_single_reference_value,                                                         &
     2715               use_surface_fluxes,                                                                 &
    27372716               use_top_fluxes
    27382717
    2739     USE grid_variables,                                                        &
    2740         ONLY:  ddx, dx, ddy, dy
    2741 
    2742     USE bulk_cloud_model_mod,                                                  &
     2718    USE grid_variables,                                                                            &
     2719        ONLY:  ddx,                                                                                &
     2720               dx,                                                                                 &
     2721               ddy,                                                                                &
     2722               dy
     2723
     2724    USE bulk_cloud_model_mod,                                                                      &
    27432725        ONLY:  bulk_cloud_model
    27442726
    27452727    IMPLICIT NONE
    27462728
    2747     LOGICAL :: diss_production
    2748 
    2749     INTEGER(iwp) ::  i       !< running index x-direction
    2750     INTEGER(iwp) ::  j       !< running index y-direction
    2751     INTEGER(iwp) ::  k       !< running index z-direction
    2752     INTEGER(iwp) ::  l       !< running index for different surface type orientation
    2753     INTEGER(iwp) ::  m       !< running index surface elements
    2754     INTEGER(iwp) ::  surf_e  !< end index of surface elements at given i-j position
    2755     INTEGER(iwp) ::  surf_s  !< start index of surface elements at given i-j position
    2756     INTEGER(iwp) ::  flag_nr !< number of masking flag
     2729    LOGICAL ::  diss_production  !<
     2730
     2731    INTEGER(iwp) ::  flag_nr  !< number of masking flag
     2732    INTEGER(iwp) ::  i        !< running index x-direction
     2733    INTEGER(iwp) ::  j        !< running index y-direction
     2734    INTEGER(iwp) ::  k        !< running index z-direction
     2735    INTEGER(iwp) ::  l        !< running index for different surface type orientation
     2736    INTEGER(iwp) ::  m        !< running index surface elements
     2737    INTEGER(iwp) ::  surf_e   !< end index of surface elements at given i-j position
     2738    INTEGER(iwp) ::  surf_s   !< start index of surface elements at given i-j position
    27572739
    27582740    REAL(wp)     ::  def         !< ( du_i/dx_j + du_j/dx_i ) * du_i/dx_j
    … …  
    27612743    REAL(wp)     ::  k2          !< temporary factor
    27622744    REAL(wp)     ::  km_neutral  !< diffusion coefficient assuming neutral conditions - used to compute shear production at surfaces
     2745    REAL(wp)     ::  sign_dir    !< sign of wall-tke flux, depending on wall orientation
    27632746    REAL(wp)     ::  theta       !< virtual potential temperature
    27642747    REAL(wp)     ::  temp        !< theta * Exner-function
    2765     REAL(wp)     ::  sign_dir    !< sign of wall-tke flux, depending on wall orientation
    27662748    REAL(wp)     ::  usvs        !< momentum flux u"v"
    27672749    REAL(wp)     ::  vsus        !< momentum flux v"u"
    … …  
    27692751    REAL(wp)     ::  wsvs        !< momentum flux w"v"
    27702752
    2771     REAL(wp), DIMENSION(nzb+1:nzt) ::  dudx  !< Gradient of u-component in x-direction
    2772     REAL(wp), DIMENSION(nzb+1:nzt) ::  dudy  !< Gradient of u-component in y-direction
    2773     REAL(wp), DIMENSION(nzb+1:nzt) ::  dudz  !< Gradient of u-component in z-direction
    2774     REAL(wp), DIMENSION(nzb+1:nzt) ::  dvdx  !< Gradient of v-component in x-direction
    2775     REAL(wp), DIMENSION(nzb+1:nzt) ::  dvdy  !< Gradient of v-component in y-direction
    2776     REAL(wp), DIMENSION(nzb+1:nzt) ::  dvdz  !< Gradient of v-component in z-direction
    2777     REAL(wp), DIMENSION(nzb+1:nzt) ::  dwdx  !< Gradient of w-component in x-direction
    2778     REAL(wp), DIMENSION(nzb+1:nzt) ::  dwdy  !< Gradient of w-component in y-direction
    2779     REAL(wp), DIMENSION(nzb+1:nzt) ::  dwdz  !< Gradient of w-component in z-direction
     2753    REAL(wp), DIMENSION(nzb+1:nzt) ::  dudx      !< Gradient of u-component in x-direction
     2754    REAL(wp), DIMENSION(nzb+1:nzt) ::  dudy      !< Gradient of u-component in y-direction
     2755    REAL(wp), DIMENSION(nzb+1:nzt) ::  dudz      !< Gradient of u-component in z-direction
     2756    REAL(wp), DIMENSION(nzb+1:nzt) ::  dvdx      !< Gradient of v-component in x-direction
     2757    REAL(wp), DIMENSION(nzb+1:nzt) ::  dvdy      !< Gradient of v-component in y-direction
     2758    REAL(wp), DIMENSION(nzb+1:nzt) ::  dvdz      !< Gradient of v-component in z-direction
     2759    REAL(wp), DIMENSION(nzb+1:nzt) ::  dwdx      !< Gradient of w-component in x-direction
     2760    REAL(wp), DIMENSION(nzb+1:nzt) ::  dwdy      !< Gradient of w-component in y-direction
     2761    REAL(wp), DIMENSION(nzb+1:nzt) ::  dwdz      !< Gradient of w-component in z-direction
    27802762    REAL(wp), DIMENSION(nzb+1:nzt) ::  tmp_flux  !< temporary flux-array in z-direction
    27812763
    … …  
    27832765
    27842766!
    2785 !-- Calculate TKE production by shear. Calculate gradients at all grid
    2786 !-- points first, gradients at surface-bounded grid points will be
    2787 !-- overwritten further below.
     2767!-- Calculate TKE production by shear. Calculate gradients at all grid points first, gradients at
     2768!-- surface-bounded grid points will be overwritten further below.
    27882769    !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j, l) &
    27892770    !$ACC PRIVATE(surf_s, surf_e) &
    … …  
    28002781
    28012782             dudx(k) =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
    2802              dudy(k) = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) -                 &
    2803                                    u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
    2804              dudz(k) = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) -                 &
    2805                                    u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
    2806 
    2807              dvdx(k) = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) -                 &
    2808                                    v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
     2783             dudy(k) = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
     2784             dudz(k) = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
     2785
     2786             dvdx(k) = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
    28092787             dvdy(k) =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
    2810              dvdz(k) = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) -                 &
    2811                                      v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
    2812 
    2813              dwdx(k) = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) -                 &
    2814                                    w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
    2815              dwdy(k) = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) -                 &
    2816                                    w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
     2788             dvdz(k) = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
     2789
     2790             dwdx(k) = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
     2791             dwdy(k) = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
    28172792             dwdz(k) =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
    28182793
    … …  
    28322807!--          'bottom and wall: use u_0,v_0 and wall functions'
    28332808!
    2834 !--          Compute gradients at north- and south-facing surfaces.
    2835 !--          First, for default surfaces, then for urban surfaces.
    2836 !--          Note, so far no natural vertical surfaces implemented
     2809!--          Compute gradients at north- and south-facing surfaces. First, for default surfaces,
     2810!--          then for urban surfaces. Note, so far no natural vertical surfaces implemented
    28372811             DO  l = 0, 1
    28382812                surf_s = surf_def_v(l)%start_index(j,i)
    … …  
    28442818                   wsvs        = surf_def_v(l)%mom_flux_tke(1,m)
    28452819
    2846                    km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp         &
    2847                                       * 0.5_wp * dy
     2820                   km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp * 0.5_wp * dy
    28482821!
    28492822!--                -1.0 for right-facing wall, 1.0 for left-facing wall
    2850                    sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
    2851                               BTEST( wall_flags_total_0(k,j-1,i), flag_nr ) )
     2823                   sign_dir = MERGE( 1.0_wp, -1.0_wp, BTEST( wall_flags_total_0(k,j-1,i), flag_nr ) )
    28522824                   dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
    28532825                   dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
    … …  
    28632835                   wsvs        = surf_lsm_v(l)%mom_flux_tke(1,m)
    28642836
    2865                    km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp         &
    2866                                       * 0.5_wp * dy
     2837                   km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp * 0.5_wp * dy
    28672838!
    28682839!--                -1.0 for right-facing wall, 1.0 for left-facing wall
    2869                    sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
    2870                               BTEST( wall_flags_total_0(k,j-1,i), flag_nr ) )
     2840                   sign_dir = MERGE( 1.0_wp, -1.0_wp, BTEST( wall_flags_total_0(k,j-1,i), flag_nr ) )
    28712841                   dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
    28722842                   dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
    … …  
    28822852                   wsvs        = surf_usm_v(l)%mom_flux_tke(1,m)
    28832853
    2884                    km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp         &
    2885                                       * 0.5_wp * dy
     2854                   km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp * 0.5_wp * dy
    28862855!
    28872856!--                -1.0 for right-facing wall, 1.0 for left-facing wall
    2888                    sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
    2889                               BTEST( wall_flags_total_0(k,j-1,i), flag_nr ) )
     2857                   sign_dir = MERGE( 1.0_wp, -1.0_wp, BTEST( wall_flags_total_0(k,j-1,i), flag_nr ) )
    28902858                   dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
    28912859                   dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
    … …  
    29032871                   wsus  = surf_def_v(l)%mom_flux_tke(1,m)
    29042872
    2905                    km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp         &
    2906                                       * 0.5_wp * dx
     2873                   km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp * 0.5_wp * dx
    29072874!
    29082875!--                -1.0 for right-facing wall, 1.0 for left-facing wall
    2909                    sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
    2910                               BTEST( wall_flags_total_0(k,j,i-1), flag_nr ) )
     2876                   sign_dir = MERGE( 1.0_wp, -1.0_wp, BTEST( wall_flags_total_0(k,j,i-1), flag_nr ) )
    29112877                   dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
    29122878                   dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
    … …  
    29222888                   wsus  = surf_lsm_v(l)%mom_flux_tke(1,m)
    29232889
    2924                    km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp         &
    2925                                       * 0.5_wp * dx
     2890                   km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp * 0.5_wp * dx
    29262891!
    29272892!--                -1.0 for right-facing wall, 1.0 for left-facing wall
    2928                    sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
    2929                               BTEST( wall_flags_total_0(k,j,i-1), flag_nr ) )
     2893                   sign_dir = MERGE( 1.0_wp, -1.0_wp, BTEST( wall_flags_total_0(k,j,i-1), flag_nr ) )
    29302894                   dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
    29312895                   dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
    … …  
    29412905                   wsus  = surf_usm_v(l)%mom_flux_tke(1,m)
    29422906
    2943                    km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp         &
    2944                                       * 0.5_wp * dx
     2907                   km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp * 0.5_wp * dx
    29452908!
    29462909!--                -1.0 for right-facing wall, 1.0 for left-facing wall
    2947                    sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
    2948                               BTEST( wall_flags_total_0(k,j,i-1), flag_nr ) )
     2910                   sign_dir = MERGE( 1.0_wp, -1.0_wp, BTEST( wall_flags_total_0(k,j,i-1), flag_nr ) )
    29492911                   dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
    29502912                   dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
    … …  
    29592921                k = surf_def_h(0)%k(m)
    29602922!
    2961 !--             Please note, actually, an interpolation of u_0 and v_0
    2962 !--             onto the grid center would be required. However, this
    2963 !--             would require several data transfers between 2D-grid and
    2964 !--             wall type. The effect of this missing interpolation is
    2965 !--             negligible. (See also production_e_init).
     2923!--             Please note, actually, an interpolation of u_0 and v_0 onto the grid center would be
     2924!--             required. However, this would require several data transfers between 2D-grid and
     2925!--             wall type. The effect of this missing interpolation is negligible.
     2926!--            (See also production_e_init).
    29662927                dudz(k) = ( u(k+1,j,i) - surf_def_h(0)%u_0(m) ) * dd2zu(k)
    29672928                dvdz(k) = ( v(k+1,j,i) - surf_def_h(0)%v_0(m) ) * dd2zu(k)
    … …  
    29932954             ENDDO
    29942955!
    2995 !--          Compute gradients at downward-facing walls, only for
    2996 !--          non-natural default surfaces
     2956!--          Compute gradients at downward-facing walls, only for non-natural default surfaces
    29972957             surf_s = surf_def_h(1)%start_index(j,i)
    29982958             surf_e = surf_def_h(1)%end_index(j,i)
    … …  
    30132973          DO  k = nzb+1, nzt
    30142974
    3015              def = 2.0_wp * ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 ) +         &
    3016                               dudy(k)**2 + dvdx(k)**2 + dwdx(k)**2 +           &
    3017                               dwdy(k)**2 + dudz(k)**2 + dvdz(k)**2 +           &
    3018                    2.0_wp * ( dvdx(k)*dudy(k) + dwdx(k)*dudz(k) +              &
    3019                               dwdy(k)*dvdz(k) )
     2975             def = 2.0_wp * ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 ) + dudy(k)**2 + dvdx(k)**2 +   &
     2976                   dwdx(k)**2 + dwdy(k)**2 + dudz(k)**2 + dvdz(k)**2 +                             &
     2977                   2.0_wp * ( dvdx(k)*dudy(k) + dwdx(k)*dudz(k) + dwdy(k)*dvdz(k) )
    30202978
    30212979             IF ( def < 0.0_wp )  def = 0.0_wp
    … …  
    30312989
    30322990!--             RANS mode: Compute tendency for dissipation-rate-production from shear
    3033                 tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag *           &
    3034                               diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) * c_1
     2991                tend(k,j,i) = tend(k,j,i) + km(k,j,i) * def * flag * diss(k,j,i) /                 &
     2992                              ( e(k,j,i) + 1.0E-20_wp ) * c_1
    30352993
    30362994             ENDIF
    … …  
    30503008          IF ( ocean_mode )  THEN
    30513009!
    3052 !--          So far in the ocean no special treatment of density flux
    3053 !--          in the bottom and top surface layer
     3010!--          So far in the ocean no special treatment of density flux in the bottom and top surface
     3011!--          layer
    30543012             DO  i = nxl, nxr
    30553013                DO  j = nys, nyn
    … …  
    30593017                   ENDDO
    30603018!
    3061 !--                Treatment of near-surface grid points, at up- and down-
    3062 !--                ward facing surfaces
     3019!--                Treatment of near-surface grid points, at up- and down-ward facing surfaces
    30633020                   IF ( use_surface_fluxes )  THEN
    30643021                      DO  l = 0, 1
    … …  
    30853042!--                   Compute tendency for TKE-production from shear
    30863043                      DO  k = nzb+1, nzt
    3087                          flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_total_0(k,j,i),0) )
    3088                          tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
    3089                                        MERGE( rho_reference, prho(k,j,i),       &
     3044                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
     3045                         tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g /                    &
     3046                                       MERGE( rho_reference, prho(k,j,i),                          &
    30903047                                              use_single_reference_value ) )
    30913048                      ENDDO
    … …  
    30953052!--                   RANS mode: Compute tendency for dissipation-rate-production from shear
    30963053                      DO  k = nzb+1, nzt
    3097                          flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_total_0(k,j,i),0) )
    3098                          tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
    3099                                        MERGE( rho_reference, prho(k,j,i),       &
    3100                                               use_single_reference_value ) ) *  &
    3101                                        diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) *  &
    3102                                        c_3
     3054                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
     3055                         tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g /                    &
     3056                                       MERGE( rho_reference, prho(k,j,i),                          &
     3057                                              use_single_reference_value ) ) *                     &
     3058                                       diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) *  c_3
    31033059                      ENDDO
    31043060
    … …  
    31733129                     !$ACC LOOP PRIVATE(k, flag)
    31743130                      DO  k = nzb+1, nzt
    3175                          flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_total_0(k,j,i),0) )
    3176                          tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
    3177                                        MERGE( pt_reference, pt(k,j,i),          &
     3131                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
     3132                         tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g /                    &
     3133                                       MERGE( pt_reference, pt(k,j,i),                             &
    31783134                                              use_single_reference_value ) )
    31793135                      ENDDO
    … …  
    31833139!--                   RANS mode: Compute tendency for dissipation-rate-production from shear
    31843140                      DO  k = nzb+1, nzt
    3185                          flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_total_0(k,j,i),0) )
    3186                          tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
    3187                                        MERGE( pt_reference, pt(k,j,i),          &
    3188                                               use_single_reference_value ) ) *  &
    3189                                        diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) *  &
    3190                                        c_3
     3141                         flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
     3142                         tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g /                    &
     3143                                       MERGE( pt_reference, pt(k,j,i),                             &
     3144                                              use_single_reference_value ) ) *                     &
     3145                                       diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) * c_3
    31913146                      ENDDO
    31923147
    … …  
    31963151             ENDDO
    31973152
    3198           ENDIF ! from IF ( .NOT. ocean_mode )
    3199 
    3200        ELSE ! or IF ( humidity )  THEN
     3153          ENDIF ! From IF ( .NOT. ocean_mode )
     3154
     3155       ELSE ! Or IF ( humidity )  THEN
    32013156
    32023157          DO  i = nxl, nxr
    … …  
    32053160                DO  k = nzb+1, nzt
    32063161
    3207                    IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
     3162                   IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets )  THEN
    32083163                      k1 = 1.0_wp + 0.61_wp * q(k,j,i)
    32093164                      k2 = 0.61_wp * pt(k,j,i)
    3210                       tmp_flux(k) = -1.0_wp * kh(k,j,i) *                      &
    3211                                       ( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) +   &
    3212                                         k2 * ( q(k+1,j,i)  - q(k-1,j,i) )      &
    3213                                       ) * dd2zu(k)
     3165                      tmp_flux(k) = -1.0_wp * kh(k,j,i) * ( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) +   &
     3166                                                            k2 * ( q(k+1,j,i)  - q(k-1,j,i) ) )    &
     3167                                    * dd2zu(k)
    32143168                   ELSE IF ( bulk_cloud_model )  THEN
    32153169                      IF ( ql(k,j,i) == 0.0_wp )  THEN
    … …  
    32193173                         theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
    32203174                         temp  = theta * exner(k)
    3221                          k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *                  &
    3222                                        ( q(k,j,i) - ql(k,j,i) ) *              &
    3223                               ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /        &
    3224                               ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *         &
     3175                         k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * ( q(k,j,i) - ql(k,j,i) ) *           &
     3176                              ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /                            &
     3177                              ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *                             &
    32253178                              ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
    32263179                         k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
    32273180                      ENDIF
    3228                       tmp_flux(k) = -1.0_wp * kh(k,j,i) *                      &
    3229                                       ( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) +   &
    3230                                         k2 * ( q(k+1,j,i)  - q(k-1,j,i) )      &
    3231                                       ) * dd2zu(k)
     3181                      tmp_flux(k) = -1.0_wp * kh(k,j,i) * ( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) +   &
     3182                                                            k2 * ( q(k+1,j,i)  - q(k-1,j,i) ) )    &
     3183                                    * dd2zu(k)
    32323184                   ELSE IF ( cloud_droplets )  THEN
    32333185                      k1 = 1.0_wp + 0.61_wp * q(k,j,i) - ql(k,j,i)
    32343186                      k2 = 0.61_wp * pt(k,j,i)
    3235                       tmp_flux(k) = -1.0_wp * kh(k,j,i) * &
    3236                                       ( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) +   &
    3237                                         k2 * ( q(k+1,j,i)  - q(k-1,j,i) ) -    &
    3238                                         pt(k,j,i) * ( ql(k+1,j,i) -            &
    3239                                         ql(k-1,j,i) ) ) * dd2zu(k)
     3187                      tmp_flux(k) = -1.0_wp * kh(k,j,i) * ( k1 * ( pt(k+1,j,i) - pt(k-1,j,i) ) +   &
     3188                                                            k2 * ( q(k+1,j,i)  - q(k-1,j,i) ) -    &
     3189                                                            pt(k,j,i) * ( ql(k+1,j,i) -            &
     3190                                                            ql(k-1,j,i) ) ) * dd2zu(k)
    32403191                   ENDIF
    32413192
    … …  
    32513202                         k = surf_def_h(l)%k(m)
    32523203
    3253                          IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
     3204                         IF (  .NOT.  bulk_cloud_model  .AND.  .NOT.  cloud_droplets ) THEN
    32543205                            k1 = 1.0_wp + 0.61_wp * q(k,j,i)
    32553206                            k2 = 0.61_wp * pt(k,j,i)
    … …  
    32613212                               theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
    32623213                               temp  = theta * exner(k)
    3263                                k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *            &
    3264                                           ( q(k,j,i) - ql(k,j,i) ) *           &
    3265                                  ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /     &
    3266                                  ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *      &
    3267                                  ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
     3214                               k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * ( q(k,j,i) - ql(k,j,i) ) *     &
     3215                                    ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /                      &
     3216                                    ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *                       &
     3217                                    ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
    32683218                               k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
    32693219                            ENDIF
    … …  
    32733223                         ENDIF
    32743224
    3275                          tmp_flux(k) = ( k1 * surf_def_h(l)%shf(m) +           &
    3276                                          k2 * surf_def_h(l)%qsws(m)            &
    3277                                        ) * drho_air_zw(k-1)
     3225                         tmp_flux(k) = ( k1 * surf_def_h(l)%shf(m) + k2 * surf_def_h(l)%qsws(m) )  &
     3226                                         * drho_air_zw(k-1)
    32783227                      ENDDO
    32793228                   ENDDO
    … …  
    32853234                      k = surf_lsm_h%k(m)
    32863235
    3287                       IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
     3236                      IF (  .NOT.  bulk_cloud_model  .AND.  .NOT.  cloud_droplets ) THEN
    32883237                         k1 = 1.0_wp + 0.61_wp * q(k,j,i)
    32893238                         k2 = 0.61_wp * pt(k,j,i)
    … …  
    32953244                            theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
    32963245                            temp  = theta * exner(k)
    3297                             k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *               &
    3298                                           ( q(k,j,i) - ql(k,j,i) ) *           &
    3299                                  ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /     &
    3300                                  ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *      &
     3246                            k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * ( q(k,j,i) - ql(k,j,i) ) *        &
     3247                                 ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /                         &
     3248                                 ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *                          &
    33013249                                 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
    33023250                            k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
    … …  
    33073255                      ENDIF
    33083256
    3309                       tmp_flux(k) = ( k1 * surf_lsm_h%shf(m) +                 &
    3310                                       k2 * surf_lsm_h%qsws(m)                  &
    3311                                     ) * drho_air_zw(k-1)
     3257                      tmp_flux(k) = ( k1 * surf_lsm_h%shf(m) + k2 * surf_lsm_h%qsws(m) )           &
     3258                                    * drho_air_zw(k-1)
    33123259                   ENDDO
    33133260!
    … …  
    33183265                      k = surf_usm_h%k(m)
    33193266
    3320                       IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
     3267                      IF ( .NOT. bulk_cloud_model  .AND.  .NOT. cloud_droplets ) THEN
    33213268                         k1 = 1.0_wp + 0.61_wp * q(k,j,i)
    33223269                         k2 = 0.61_wp * pt(k,j,i)
    … …  
    33283275                            theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
    33293276                            temp  = theta * exner(k)
    3330                             k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *               &
    3331                                           ( q(k,j,i) - ql(k,j,i) ) *           &
    3332                                  ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /     &
    3333                                  ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *      &
     3277                            k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * ( q(k,j,i) - ql(k,j,i) ) *        &
     3278                                 ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /                         &
     3279                                 ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *                          &
    33343280                                 ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
    33353281                            k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
    … …  
    33403286                      ENDIF
    33413287
    3342                       tmp_flux(k) = ( k1 * surf_usm_h%shf(m) +                 &
    3343                                       k2 * surf_usm_h%qsws(m)                  &
    3344                                     ) * drho_air_zw(k-1)
     3288                      tmp_flux(k) = ( k1 * surf_usm_h%shf(m) + k2 * surf_usm_h%qsws(m) )           &
     3289                                    * drho_air_zw(k-1)
    33453290                   ENDDO
    33463291
    3347                 ENDIF ! from IF ( use_surface_fluxes )  THEN
     3292                ENDIF ! From IF ( use_surface_fluxes )  THEN
    33483293
    33493294                IF ( use_top_fluxes )  THEN
    … …  
    33543299                      k = surf_def_h(2)%k(m)
    33553300
    3356                       IF ( .NOT. bulk_cloud_model .AND. .NOT. cloud_droplets ) THEN
     3301                      IF ( .NOT. bulk_cloud_model  .AND.  .NOT. cloud_droplets ) THEN
    33573302                         k1 = 1.0_wp + 0.61_wp * q(k,j,i)
    33583303                         k2 = 0.61_wp * pt(k,j,i)
    … …  
    33643309                            theta = pt(k,j,i) + d_exner(k) * lv_d_cp * ql(k,j,i)
    33653310                            temp  = theta * exner(k)
    3366                             k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp *               &
    3367                                        ( q(k,j,i) - ql(k,j,i) ) *              &
    3368                               ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /        &
    3369                               ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *         &
     3311                            k1 = ( 1.0_wp - q(k,j,i) + 1.61_wp * ( q(k,j,i) - ql(k,j,i) ) *        &
     3312                              ( 1.0_wp + rd_d_rv * lv_d_rd / temp ) ) /                            &
     3313                              ( 1.0_wp + rd_d_rv * lv_d_rd * lv_d_cp *                             &
    33703314                              ( q(k,j,i) - ql(k,j,i) ) / ( temp * temp ) )
    33713315                            k2 = theta * ( lv_d_cp / temp * k1 - 1.0_wp )
    … …  
    33763320                      ENDIF
    33773321
    3378                       tmp_flux(k) = ( k1 * surf_def_h(2)%shf(m) +              &
    3379                                       k2 * surf_def_h(2)%qsws(m)               &
    3380                                     ) * drho_air_zw(k)
     3322                      tmp_flux(k) = ( k1 * surf_def_h(2)%shf(m) + k2 * surf_def_h(2)%qsws(m) )     &
     3323                                    * drho_air_zw(k)
    33813324
    33823325                   ENDDO
    33833326
    3384                 ENDIF ! from IF ( use_top_fluxes )  THEN
     3327                ENDIF ! From IF ( use_top_fluxes )  THEN
    33853328
    33863329                IF ( .NOT. diss_production )  THEN
    … …  
    33883331!--                Compute tendency for TKE-production from shear
    33893332                   DO  k = nzb+1, nzt
    3390                       flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_total_0(k,j,i),0) )
    3391                       tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
    3392                                     MERGE( vpt_reference, vpt(k,j,i),          &
    3393                                            use_single_reference_value ) )
     3333                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
     3334                      tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g /                       &
     3335                                    MERGE( vpt_reference, vpt(k,j,i), use_single_reference_value ) )
    33943336                   ENDDO
    33953337
    … …  
    33983340!--                RANS mode: Compute tendency for dissipation-rate-production from shear
    33993341                   DO  k = nzb+1, nzt
    3400                       flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_total_0(k,j,i),0) )
    3401                       tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g /   &
    3402                                     MERGE( vpt_reference, vpt(k,j,i),          &
    3403                                            use_single_reference_value ) ) *    &
    3404                                     diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) *    &
    3405                                     c_3
     3342                      flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
     3343                      tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) *                             &
     3344                                    ( g / MERGE( vpt_reference, vpt(k,j,i),                        &
     3345                                                  use_single_reference_value ) )                   &
     3346                                   * diss(k,j,i)/( e(k,j,i) + 1.0E-20_wp ) * c_3
    34063347                   ENDDO
    34073348
    … …  
    34183359
    34193360
    3420 !------------------------------------------------------------------------------!
     3361!--------------------------------------------------------------------------------------------------!
    34213362! Description:
    34223363! ------------
    34233364!> Production terms (shear + buoyancy) of the TKE.
    34243365!> Cache-optimized version
    3425 !> @warning The case with constant_flux_layer = F and use_surface_fluxes = T is
    3426 !>          not considered well!
    3427 !------------------------------------------------------------------------------!
     3366!> @warning The case with constant_flux_layer = F and use_surface_fluxes = T is not considered well!
     3367!--------------------------------------------------------------------------------------------------!
    34283368 SUBROUTINE production_e_ij( i, j, diss_production )
    34293369
    3430     USE arrays_3d,                                                             &
    3431         ONLY:  ddzw, dd2zu, drho_air_zw, q, ql, d_exner, exner
    3432 
    3433     USE control_parameters,                                                    &
    3434         ONLY:  cloud_droplets, constant_flux_layer, neutral,                   &
    3435                rho_reference, use_single_reference_value, use_surface_fluxes,  &
     3370    USE arrays_3d,                                                                                 &
     3371        ONLY:  ddzw,                                                                               &
     3372               dd2zu,                                                                              &
     3373               drho_air_zw,                                                                        &
     3374               d_exner,                                                                            &
     3375               exner,                                                                              &
     3376               q,                                                                                  &
     3377               ql
     3378
     3379
     3380
     3381    USE control_parameters,                                                                        &
     3382        ONLY:  cloud_droplets,                                                                     &
     3383               constant_flux_layer,                                                                &
     3384               neutral,                                                                            &
     3385               rho_reference,                                                                      &
     3386               use_single_reference_value,                                                         &
     3387               use_surface_fluxes,                                                                 &
    34363388               use_top_fluxes
    34373389
    3438     USE grid_variables,                                                        &
    3439         ONLY:  ddx, dx, ddy, dy
    3440 
    3441     USE bulk_cloud_model_mod,                                                  &
     3390    USE grid_variables,                                                                            &
     3391        ONLY:  ddx,                                                                                &
     3392               dx,                                                                                 &
     3393               ddy,                                                                                &
     3394               dy
     3395
     3396    USE bulk_cloud_model_mod,                                                                      &
    34423397        ONLY:  bulk_cloud_model
    34433398
    34443399    IMPLICIT NONE
    34453400
    3446     LOGICAL :: diss_production
    3447 
    3448     INTEGER(iwp) ::  i       !< running index x-direction
    3449     INTEGER(iwp) ::  j       !< running index y-direction
    3450     INTEGER(iwp) ::  k       !< running index z-direction
    3451     INTEGER(iwp) ::  l       !< running index for different surface type orientation
    3452     INTEGER(iwp) ::  m       !< running index surface elements
    3453     INTEGER(iwp) ::  surf_e  !< end index of surface elements at given i-j position
    3454     INTEGER(iwp) ::  surf_s  !< start index of surface elements at given i-j position
    3455     INTEGER(iwp) ::  flag_nr !< number of masking flag
     3401    LOGICAL ::  diss_production  !<
     3402
     3403    INTEGER(iwp) ::  flag_nr  !< number of masking flag
     3404    INTEGER(iwp) ::  i        !< running index x-direction
     3405    INTEGER(iwp) ::  j        !< running index y-direction
     3406    INTEGER(iwp) ::  k        !< running index z-direction
     3407    INTEGER(iwp) ::  l        !< running index for different surface type orientation
     3408    INTEGER(iwp) ::  m        !< running index surface elements
     3409    INTEGER(iwp) ::  surf_e   !< end index of surface elements at given i-j position
     3410    INTEGER(iwp) ::  surf_s   !< start index of surface elements at given i-j position
    34563411
    34573412    REAL(wp)     ::  def         !< ( du_i/dx_j + du_j/dx_i ) * du_i/dx_j
    … …  
    34603415    REAL(wp)     ::  k2          !< temporary factor
    34613416    REAL(wp)     ::  km_neutral  !< diffusion coefficient assuming neutral conditions - used to compute shear production at surfaces
     3417    REAL(wp)     ::  sign_dir    !< sign of wall-tke flux, depending on wall orientation
    34623418    REAL(wp)     ::  theta       !< virtual potential temperature
    34633419    REAL(wp)     ::  temp        !< theta * Exner-function
    3464     REAL(wp)     ::  sign_dir    !< sign of wall-tke flux, depending on wall orientation
    34653420    REAL(wp)     ::  usvs        !< momentum flux u"v"
    34663421    REAL(wp)     ::  vsus        !< momentum flux v"u"
    … …  
    34683423    REAL(wp)     ::  wsvs        !< momentum flux w"v"
    34693424
    3470     REAL(wp), DIMENSION(nzb+1:nzt) ::  dudx  !< Gradient of u-component in x-direction
    3471     REAL(wp), DIMENSION(nzb+1:nzt) ::  dudy  !< Gradient of u-component in y-direction
    3472     REAL(wp), DIMENSION(nzb+1:nzt) ::  dudz  !< Gradient of u-component in z-direction
    3473     REAL(wp), DIMENSION(nzb+1:nzt) ::  dvdx  !< Gradient of v-component in x-direction
    3474     REAL(wp), DIMENSION(nzb+1:nzt) ::  dvdy  !< Gradient of v-component in y-direction
    3475     REAL(wp), DIMENSION(nzb+1:nzt) ::  dvdz  !< Gradient of v-component in z-direction
    3476     REAL(wp), DIMENSION(nzb+1:nzt) ::  dwdx  !< Gradient of w-component in x-direction
    3477     REAL(wp), DIMENSION(nzb+1:nzt) ::  dwdy  !< Gradient of w-component in y-direction
    3478     REAL(wp), DIMENSION(nzb+1:nzt) ::  dwdz  !< Gradient of w-component in z-direction
     3425    REAL(wp), DIMENSION(nzb+1:nzt) ::  dudx      !< Gradient of u-component in x-direction
     3426    REAL(wp), DIMENSION(nzb+1:nzt) ::  dudy      !< Gradient of u-component in y-direction
     3427    REAL(wp), DIMENSION(nzb+1:nzt) ::  dudz      !< Gradient of u-component in z-direction
     3428    REAL(wp), DIMENSION(nzb+1:nzt) ::  dvdx      !< Gradient of v-component in x-direction
     3429    REAL(wp), DIMENSION(nzb+1:nzt) ::  dvdy      !< Gradient of v-component in y-direction
     3430    REAL(wp), DIMENSION(nzb+1:nzt) ::  dvdz      !< Gradient of v-component in z-direction
     3431    REAL(wp), DIMENSION(nzb+1:nzt) ::  dwdx      !< Gradient of w-component in x-direction
     3432    REAL(wp), DIMENSION(nzb+1:nzt) ::  dwdy      !< Gradient of w-component in y-direction
     3433    REAL(wp), DIMENSION(nzb+1:nzt) ::  dwdz      !< Gradient of w-component in z-direction
    34793434    REAL(wp), DIMENSION(nzb+1:nzt) ::  tmp_flux  !< temporary flux-array in z-direction
    34803435
    … …  
    34823437
    34833438!
    3484 !-- Calculate TKE production by shear. Calculate gradients at all grid
    3485 !-- points first, gradients at surface-bounded grid points will be
    3486 !-- overwritten further below.
     3439!-- Calculate TKE production by shear. Calculate gradients at all grid points first, gradients at
     3440!-- surface-bounded grid points will be overwritten further below.
    34873441    DO  k = nzb+1, nzt
    34883442
    34893443       dudx(k) =           ( u(k,j,i+1) - u(k,j,i)     ) * ddx
    3490        dudy(k) = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) -                 &
    3491                              u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
    3492        dudz(k) = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) -                 &
    3493                              u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
    3494 
    3495        dvdx(k) = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) -                 &
    3496                              v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
     3444       dudy(k) = 0.25_wp * ( u(k,j+1,i) + u(k,j+1,i+1) - u(k,j-1,i) - u(k,j-1,i+1) ) * ddy
     3445       dudz(k) = 0.5_wp  * ( u(k+1,j,i) + u(k+1,j,i+1) - u(k-1,j,i) - u(k-1,j,i+1) ) * dd2zu(k)
     3446
     3447       dvdx(k) = 0.25_wp * ( v(k,j,i+1) + v(k,j+1,i+1) - v(k,j,i-1) - v(k,j+1,i-1) ) * ddx
    34973448       dvdy(k) =           ( v(k,j+1,i) - v(k,j,i)     ) * ddy
    3498        dvdz(k) = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) -                 &
    3499                              v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
    3500 
    3501        dwdx(k) = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) -                 &
    3502                              w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
    3503        dwdy(k) = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) -                 &
    3504                              w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
     3449       dvdz(k) = 0.5_wp  * ( v(k+1,j,i) + v(k+1,j+1,i) - v(k-1,j,i) - v(k-1,j+1,i) ) * dd2zu(k)
     3450
     3451       dwdx(k) = 0.25_wp * ( w(k,j,i+1) + w(k-1,j,i+1) - w(k,j,i-1) - w(k-1,j,i-1) ) * ddx
     3452       dwdy(k) = 0.25_wp * ( w(k,j+1,i) + w(k-1,j+1,i) - w(k,j-1,i) - w(k-1,j-1,i) ) * ddy
    35053453       dwdz(k) =           ( w(k,j,i)   - w(k-1,j,i)   ) * ddzw(k)
    35063454
    … …  
    35173465!--    'bottom and wall: use u_0,v_0 and wall functions'
    35183466!
    3519 !--    Compute gradients at north- and south-facing surfaces.
    3520 !--    First, for default surfaces, then for urban surfaces.
    3521 !--    Note, so far no natural vertical surfaces implemented
     3467!--    Compute gradients at north- and south-facing surfaces. First, for default surfaces, then for
     3468!--    urban surfaces. Note, so far no natural vertical surfaces implemented
    35223469       DO  l = 0, 1
    35233470          surf_s = surf_def_v(l)%start_index(j,i)
    … …  
    35283475             wsvs        = surf_def_v(l)%mom_flux_tke(1,m)
    35293476
    3530              km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp         &
    3531                                 * 0.5_wp * dy
     3477             km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp * 0.5_wp * dy
    35323478!
    35333479!--          -1.0 for right-facing wall, 1.0 for left-facing wall
    3534              sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
    3535                         BTEST( wall_flags_total_0(k,j-1,i), flag_nr ) )
     3480             sign_dir = MERGE( 1.0_wp, -1.0_wp, BTEST( wall_flags_total_0(k,j-1,i), flag_nr ) )
    35363481             dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
    35373482             dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
    … …  
    35463491             wsvs        = surf_lsm_v(l)%mom_flux_tke(1,m)
    35473492
    3548              km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp         &
    3549                                 * 0.5_wp * dy
     3493             km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp * 0.5_wp * dy
    35503494!
    35513495!--          -1.0 for right-facing wall, 1.0 for left-facing wall
    3552              sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
    3553                         BTEST( wall_flags_total_0(k,j-1,i), flag_nr ) )
     3496             sign_dir = MERGE( 1.0_wp, -1.0_wp, BTEST( wall_flags_total_0(k,j-1,i), flag_nr ) )
    35543497             dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
    35553498             dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
    … …  
    35643507             wsvs        = surf_usm_v(l)%mom_flux_tke(1,m)
    35653508
    3566              km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp         &
    3567                                 * 0.5_wp * dy
     3509             km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp * 0.5_wp * dy
    35683510!
    35693511!--          -1.0 for right-facing wall, 1.0 for left-facing wall
    3570              sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
    3571                         BTEST( wall_flags_total_0(k,j-1,i), flag_nr ) )
     3512             sign_dir = MERGE( 1.0_wp, -1.0_wp, BTEST( wall_flags_total_0(k,j-1,i), flag_nr ) )
    35723513             dudy(k) = sign_dir * usvs / ( km_neutral + 1E-10_wp )
    35733514             dwdy(k) = sign_dir * wsvs / ( km_neutral + 1E-10_wp )
    … …  
    35843525             wsus  = surf_def_v(l)%mom_flux_tke(1,m)
    35853526
    3586              km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp         &
    3587                                 * 0.5_wp * dx
     3527             km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp * 0.5_wp * dx
    35883528!
    35893529!--          -1.0 for right-facing wall, 1.0 for left-facing wall
    3590              sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
    3591                         BTEST( wall_flags_total_0(k,j,i-1), flag_nr ) )
     3530             sign_dir = MERGE( 1.0_wp, -1.0_wp, BTEST( wall_flags_total_0(k,j,i-1), flag_nr ) )
    35923531             dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
    35933532             dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
    … …  
    36023541             wsus  = surf_lsm_v(l)%mom_flux_tke(1,m)
    36033542
    3604              km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp         &
    3605                                 * 0.5_wp * dx
     3543             km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp * 0.5_wp * dx
    36063544!
    36073545!--          -1.0 for right-facing wall, 1.0 for left-facing wall
    3608              sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
    3609                         BTEST( wall_flags_total_0(k,j,i-1), flag_nr ) )
     3546             sign_dir = MERGE( 1.0_wp, -1.0_wp, BTEST( wall_flags_total_0(k,j,i-1), flag_nr ) )
    36103547             dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
    36113548             dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
    … …  
    36203557             wsus  = surf_usm_v(l)%mom_flux_tke(1,m)
    36213558
    3622              km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp         &
    3623                                 * 0.5_wp * dx
     3559             km_neutral = kappa * ( vsus**2 + wsus**2 )**0.25_wp * 0.5_wp * dx
    36243560!
    36253561!--          -1.0 for right-facing wall, 1.0 for left-facing wall
    3626              sign_dir = MERGE( 1.0_wp, -1.0_wp,                          &
    3627                         BTEST( wall_flags_total_0(k,j,i-1), flag_nr ) )
     3562             sign_dir = MERGE( 1.0_wp, -1.0_wp, BTEST( wall_flags_total_0(k,j,i-1), flag_nr ) )
    36283563             dvdx(k) = sign_dir * vsus / ( km_neutral + 1E-10_wp )
    36293564             dwdx(k) = sign_dir * wsus / ( km_neutral + 1E-10_wp )
    … …  
    36373572          k = surf_def_h(0)%k(m)
    36383573!
    3639 !--       Please note, actually, an interpolation of u_0 and v_0
    3640 !--       onto the grid center would be required. However, this
    3641 !--       would require several data transfers between 2D-grid and
    3642 !--       wall type. The effect of this missing interpolation is
    3643 !--       negligible. (See also production_e_init).
     3574!--       Please note, actually, an interpolation of u_0 and v_0 onto the grid center would be
     3575!--       required. However, this would require several data transfers between 2D-grid and wall
     3576!--       type. The effect of this missing interpolation is negligible. (See also production_e_init).
    36443577          dudz(k) = ( u(k+1,j,i) - surf_def_h(0)%u_0(m) ) * dd2zu(k)
    36453578          dvdz(k) = ( v(k+1,j,i) - surf_def_h(0)%v_0(m) ) * dd2zu(k)
    … …  
    36693602       ENDDO
    36703603!
    3671 !--    Compute gradients at downward-facing walls, only for
    3672 !--    non-natural default surfaces
     3604!--    Compute gradients at downward-facing walls, only for non-natural default surfaces
    36733605       surf_s = surf_def_h(1)%start_index(j,i)
    36743606       surf_e = surf_def_h(1)%end_index(j,i)
    … …  
    36853617    DO  k = nzb+1, nzt
    36863618
    3687        def = 2.0_wp * ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 ) +         &
    3688                         dudy(k)**2 + dvdx(k)**2 + dwdx(k)**2 +           &
    3689                         dwdy(k)**2 + dudz(k)**2 + dvdz(k)**2 +           &
    3690              2.0_wp * ( dvdx(k)*dudy(k) + dwdx(k)*dudz(k) +              &
    3691                         dwdy(k)*dvdz(k) )
     3619       def = 2.0_wp * ( dudx(k)**2 + dvdy(k)**2 + dwdz(k)**2 ) +                                   &
     3620                        dudy(k)**2 + dvdx(k)**2 + dwdx(k)**2 +                                     &
     3621                        dwdy(k)**2 + dudz(k)**2 + dvdz(k)**2 +                                     &
     3622             2.0_wp * ( dvdx(k)*dudy(k) + dwdx(k)*dudz(k) + dwdy(k)*dvdz(k) )
    36923623
    36933624       IF ( def < 0.0_wp )  def = 0.0_wp
    … …  
    37033634
    37043635!--       RANS mode: Compute tendency for dissipation-rate-production from shear