source: palm/trunk/SOURCE/boundary_conds.f90 @ 3467

!> @file boundary_conds.f90
!------------------------------------------------------------------------------!
! This file is part of the PALM model system.
!
! PALM is free software: you can redistribute it and/or modify it under the
! terms of the GNU General Public License as published by the Free Software
! Foundation, either version 3 of the License, or (at your option) any later
! version.
!
! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
!
! You should have received a copy of the GNU General Public License along with
! PALM. If not, see <http://www.gnu.org/licenses/>.
!
! Copyright 1997-2018 Leibniz Universitaet Hannover
!------------------------------------------------------------------------------!
!
! Current revisions:
! -----------------
! 
! 
! Former revisions:
! -----------------
! $Id: boundary_conds.f90 3467 2018-10-30 19:05:21Z suehring $
! Implementation of a new aerosol module salsa.
! 
! 3347 2018-10-15 14:21:08Z suehring
! Bugfix in setting boundary conditions in offline nesting
! 
! 3341 2018-10-15 10:31:27Z suehring
! changes concerning modularization of ocean option
! 
! 3274 2018-09-24 15:42:55Z knoop
! Modularization of all bulk cloud physics code components
! 
! 3241 2018-09-12 15:02:00Z raasch
! unused variables removed
! 
! 3183 2018-07-27 14:25:55Z suehring
! Rename some variables concerning LES-LES as well as offline nesting
! 
! 3182 2018-07-27 13:36:03Z suehring
! - Use wall function for e_p and diss_p in case of rans_tke_e
! - move limitation of diss_p from tcm_prognostic
! 
! 2938 2018-03-27 15:52:42Z suehring
! Set boundary condition for TKE and TKE dissipation rate in case of nesting
! and if parent model operates in RANS mode but child model in LES mode.
! mode
! 
! 2793 2018-02-07 10:54:33Z suehring
! Removed preprocessor directive __chem
! 
! 2718 2018-01-02 08:49:38Z maronga
! Corrected "Former revisions" section
! 
! 2696 2017-12-14 17:12:51Z kanani
! Change in file header (GPL part)
! Adjust boundary conditions for e and diss in case of TKE-e closure (TG)
! Implementation of chemistry module (FK)
! 
! 2569 2017-10-20 11:54:42Z kanani
! Removed redundant code for ibc_s_b=1 and ibc_q_b=1
! 
! 2365 2017-08-21 14:59:59Z kanani
! Vertical grid nesting implemented: exclude setting vertical velocity to zero 
! on fine grid (SadiqHuq)
! 
! 2320 2017-07-21 12:47:43Z suehring
! Remove unused control parameter large_scale_forcing from only-list
! 
! 2292 2017-06-20 09:51:42Z schwenkel
! Implementation of new microphysic scheme: cloud_scheme = 'morrison' 
! includes two more prognostic equations for cloud drop concentration (nc)  
! and cloud water content (qc). 
! 
! 2233 2017-05-30 18:08:54Z suehring
!
! 2232 2017-05-30 17:47:52Z suehring
! Set boundary conditions on topography top using flag method.
! 
! 2118 2017-01-17 16:38:49Z raasch
! OpenACC directives removed
! 
! 2000 2016-08-20 18:09:15Z knoop
! Forced header and separation lines into 80 columns
! 
! 1992 2016-08-12 15:14:59Z suehring
! Adjustments for top boundary condition for passive scalar
! 
! 1960 2016-07-12 16:34:24Z suehring
! Treat humidity and passive scalar separately
! 
! 1823 2016-04-07 08:57:52Z hoffmann
! Initial version of purely vertical nesting introduced.
!
! 1822 2016-04-07 07:49:42Z hoffmann
! icloud_scheme removed. microphyisics_seifert added.
!
! 1764 2016-02-28 12:45:19Z raasch
! index bug for u_p at left outflow removed
!
! 1762 2016-02-25 12:31:13Z hellstea
! Introduction of nested domain feature
!
! 1742 2016-01-13 09:50:06Z raasch
! bugfix for outflow Neumann boundary conditions at bottom and top
!
! 1717 2015-11-11 15:09:47Z raasch
! Bugfix: index error in outflow conditions for left boundary
!
! 1682 2015-10-07 23:56:08Z knoop
! Code annotations made doxygen readable
! 
! 1410 2014-05-23 12:16:18Z suehring
! Bugfix: set dirichlet boundary condition for passive_scalar at model domain 
! top 
!
! 1399 2014-05-07 11:16:25Z heinze
! Bugfix: set inflow boundary conditions also if no humidity or passive_scalar
! is used.
! 
! 1398 2014-05-07 11:15:00Z heinze
! Dirichlet-condition at the top for u and v changed to u_init and v_init also 
! for large_scale_forcing
! 
! 1380 2014-04-28 12:40:45Z heinze
! Adjust Dirichlet-condition at the top for pt in case of nudging
! 
! 1361 2014-04-16 15:17:48Z hoffmann
! Bottom and top boundary conditions of rain water content (qr) and 
! rain drop concentration (nr) changed to Dirichlet
! 
! 1353 2014-04-08 15:21:23Z heinze
! REAL constants provided with KIND-attribute
!  
! 1320 2014-03-20 08:40:49Z raasch
! ONLY-attribute added to USE-statements,
! kind-parameters added to all INTEGER and REAL declaration statements,
! kinds are defined in new module kinds,
! revision history before 2012 removed,
! comment fields (!:) to be used for variable explanations added to
! all variable declaration statements
!
! 1257 2013-11-08 15:18:40Z raasch
! loop independent clauses added
!
! 1241 2013-10-30 11:36:58Z heinze
! Adjust ug and vg at each timestep in case of large_scale_forcing
!
! 1159 2013-05-21 11:58:22Z fricke
! Bugfix: Neumann boundary conditions for the velocity components at the
! outflow are in fact radiation boundary conditions using the maximum phase
! velocity that ensures numerical stability (CFL-condition).
! Hence, logical operator use_cmax is now used instead of bc_lr_dirneu/_neudir.
! Bugfix: In case of use_cmax at the outflow, u, v, w are replaced by
! u_p, v_p, w_p  
!
! 1115 2013-03-26 18:16:16Z hoffmann
! boundary conditions of two-moment cloud scheme are restricted to Neumann-
! boundary-conditions
!
! 1113 2013-03-10 02:48:14Z raasch
! GPU-porting
! dummy argument "range" removed
! Bugfix: wrong index in loops of radiation boundary condition
!
! 1053 2012-11-13 17:11:03Z hoffmann
! boundary conditions for the two new prognostic equations (nr, qr) of the 
! two-moment cloud scheme
!
! 1036 2012-10-22 13:43:42Z raasch
! code put under GPL (PALM 3.9)
!
! 996 2012-09-07 10:41:47Z raasch
! little reformatting
!
! 978 2012-08-09 08:28:32Z fricke
! Neumann boudnary conditions are added at the inflow boundary for the SGS-TKE.
! Outflow boundary conditions for the velocity components can be set to Neumann
! conditions or to radiation conditions with a horizontal averaged phase
! velocity. 
!
! 875 2012-04-02 15:35:15Z gryschka
! Bugfix in case of dirichlet inflow bc at the right or north boundary
!
! Revision 1.1  1997/09/12 06:21:34  raasch
! Initial revision
!
!
! Description:
! ------------
!> Boundary conditions for the prognostic quantities.
!> One additional bottom boundary condition is applied for the TKE (=(u*)**2)
!> in prandtl_fluxes. The cyclic lateral boundary conditions are implicitly
!> handled in routine exchange_horiz. Pressure boundary conditions are
!> explicitly set in routines pres, poisfft, poismg and sor.
!------------------------------------------------------------------------------!
 SUBROUTINE boundary_conds
 

    USE arrays_3d,                                                             &
        ONLY:  c_u, c_u_m, c_u_m_l, c_v, c_v_m, c_v_m_l, c_w, c_w_m, c_w_m_l,  &
               diss, diss_p, dzu, e_p, nc_p, nr_p, pt, pt_init, pt_p, q,       &
               q_p, qc_p, qr_p, s, s_p, sa, sa_p, u, u_init, u_m_l, u_m_n,     &
               u_m_r, u_m_s, u_p, v, v_init, v_m_l, v_m_n, v_m_r, v_m_s, v_p,  &
               w, w_p, w_m_l, w_m_n, w_m_r, w_m_s

    USE basic_constants_and_equations_mod,                                     &
        ONLY:  kappa

    USE bulk_cloud_model_mod,                                                  &
        ONLY: bulk_cloud_model, microphysics_morrison, microphysics_seifert

    USE chemistry_model_mod,                                                   &
        ONLY:  chem_boundary_conds 
              
    USE control_parameters,                                                    &
        ONLY:  air_chemistry, bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r,  &
               bc_dirichlet_s, bc_radiation_l, bc_radiation_n, bc_radiation_r, &
               bc_radiation_s, bc_pt_t_val, bc_q_t_val, bc_s_t_val,            &
               child_domain, constant_diffusion, coupling_mode, &
               dt_3d, humidity, ibc_pt_b, ibc_pt_t, ibc_q_b, ibc_q_t, ibc_s_b, &
               ibc_s_t, ibc_uv_b, ibc_uv_t,                           &
               intermediate_timestep_count, nesting_offline, nudging,          &
               ocean_mode, passive_scalar, rans_mode, rans_tke_e, tsc, use_cmax

    USE grid_variables,                                                        &
        ONLY:  ddx, ddy, dx, dy

    USE indices,                                                               &
        ONLY:  nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, nzb, nzt

    USE kinds

    USE ocean_mod,                                                             &
        ONLY:  ibc_sa_t

    USE pegrid

    USE pmc_interface,                                                         &
        ONLY : nesting_mode, rans_mode_parent
 
    USE salsa_mod,                                                             &
        ONLY:  salsa, salsa_boundary_conds        

    USE surface_mod,                                                           &
        ONLY :  bc_h, surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v,          &
                surf_usm_h, surf_usm_v

    USE turbulence_closure_mod,                                                &
        ONLY:  c_0

    IMPLICIT NONE

    INTEGER(iwp) ::  i  !< grid index x direction
    INTEGER(iwp) ::  j  !< grid index y direction
    INTEGER(iwp) ::  k  !< grid index z direction
    INTEGER(iwp) ::  kb !< variable to set respective boundary value, depends on facing. 
    INTEGER(iwp) ::  l  !< running index boundary type, for up- and downward-facing walls
    INTEGER(iwp) ::  m  !< running index surface elements

    REAL(wp)    ::  c_max !<
    REAL(wp)    ::  denom !<


!
!-- Bottom boundary 
    IF ( ibc_uv_b == 1 )  THEN
       u_p(nzb,:,:) = u_p(nzb+1,:,:)
       v_p(nzb,:,:) = v_p(nzb+1,:,:)
    ENDIF
!
!-- Set zero vertical velocity at topography top (l=0), or bottom (l=1) in case
!-- of downward-facing surfaces. 
    DO  l = 0, 1
!
!--    Set kb, for upward-facing surfaces value at topography top (k-1) is set,
!--    for downward-facing surfaces at topography bottom (k+1). 
       kb = MERGE( -1, 1, l == 0 )
       !$OMP PARALLEL DO PRIVATE( i, j, k )
       DO  m = 1, bc_h(l)%ns
          i = bc_h(l)%i(m)            
          j = bc_h(l)%j(m)
          k = bc_h(l)%k(m)
          w_p(k+kb,j,i) = 0.0_wp
       ENDDO
    ENDDO

!
!-- Top boundary. A nested domain ( ibc_uv_t = 3 ) does not require settings.
    IF ( ibc_uv_t == 0 )  THEN
        u_p(nzt+1,:,:) = u_init(nzt+1)
        v_p(nzt+1,:,:) = v_init(nzt+1)
    ELSEIF ( ibc_uv_t == 1 )  THEN
        u_p(nzt+1,:,:) = u_p(nzt,:,:)
        v_p(nzt+1,:,:) = v_p(nzt,:,:)
    ENDIF

!
!-- Vertical nesting: Vertical velocity not zero at the top of the fine grid
    IF (  .NOT.  child_domain  .AND.  .NOT.  nesting_offline  .AND.            &
                 TRIM(coupling_mode) /= 'vnested_fine' )  THEN
       w_p(nzt:nzt+1,:,:) = 0.0_wp  !< nzt is not a prognostic level (but cf. pres)
    ENDIF

!
!-- Temperature at bottom and top boundary.
!-- In case of coupled runs (ibc_pt_b = 2) the temperature is given by
!-- the sea surface temperature of the coupled ocean model.
!-- Dirichlet
    IF ( ibc_pt_b == 0 )  THEN
       DO  l = 0, 1
!
!--       Set kb, for upward-facing surfaces value at topography top (k-1) is set,
!--       for downward-facing surfaces at topography bottom (k+1). 
          kb = MERGE( -1, 1, l == 0 )
          !$OMP PARALLEL DO PRIVATE( i, j, k )
          DO  m = 1, bc_h(l)%ns
             i = bc_h(l)%i(m)            
             j = bc_h(l)%j(m)
             k = bc_h(l)%k(m)
             pt_p(k+kb,j,i) = pt(k+kb,j,i)
          ENDDO
       ENDDO
!
!-- Neumann, zero-gradient
    ELSEIF ( ibc_pt_b == 1 )  THEN
       DO  l = 0, 1
!
!--       Set kb, for upward-facing surfaces value at topography top (k-1) is set,
!--       for downward-facing surfaces at topography bottom (k+1). 
          kb = MERGE( -1, 1, l == 0 )
          !$OMP PARALLEL DO PRIVATE( i, j, k )
          DO  m = 1, bc_h(l)%ns
             i = bc_h(l)%i(m)            
             j = bc_h(l)%j(m)
             k = bc_h(l)%k(m)
             pt_p(k+kb,j,i) = pt_p(k,j,i)
          ENDDO
       ENDDO
    ENDIF

!
!-- Temperature at top boundary
    IF ( ibc_pt_t == 0 )  THEN
        pt_p(nzt+1,:,:) = pt(nzt+1,:,:)
!
!--     In case of nudging adjust top boundary to pt which is
!--     read in from NUDGING-DATA
        IF ( nudging )  THEN
           pt_p(nzt+1,:,:) = pt_init(nzt+1)
        ENDIF
    ELSEIF ( ibc_pt_t == 1 )  THEN
        pt_p(nzt+1,:,:) = pt_p(nzt,:,:)
    ELSEIF ( ibc_pt_t == 2 )  THEN
        pt_p(nzt+1,:,:) = pt_p(nzt,:,:) + bc_pt_t_val * dzu(nzt+1)
    ENDIF

!
!-- Boundary conditions for TKE.
!-- Generally Neumann conditions with de/dz=0 are assumed.
    IF ( .NOT. constant_diffusion )  THEN

       IF ( .NOT. rans_tke_e )  THEN
          DO  l = 0, 1
!
!--         Set kb, for upward-facing surfaces value at topography top (k-1) is set,
!--         for downward-facing surfaces at topography bottom (k+1). 
             kb = MERGE( -1, 1, l == 0 )
             !$OMP PARALLEL DO PRIVATE( i, j, k )
             DO  m = 1, bc_h(l)%ns
                i = bc_h(l)%i(m)            
                j = bc_h(l)%j(m)
                k = bc_h(l)%k(m)
                e_p(k+kb,j,i) = e_p(k,j,i)
             ENDDO
          ENDDO
       ELSE
!
!--       Use wall function within constant-flux layer
!--       Upward-facing surfaces
!--       Default surfaces
          DO  m = 1, surf_def_h(0)%ns
             i = surf_def_h(0)%i(m)
             j = surf_def_h(0)%j(m)
             k = surf_def_h(0)%k(m)
             e_p(k,j,i) = ( surf_def_h(0)%us(m) / c_0 )**2
          ENDDO
!
!--       Natural surfaces
          DO  m = 1, surf_lsm_h%ns
             i = surf_lsm_h%i(m)
             j = surf_lsm_h%j(m)
             k = surf_lsm_h%k(m)
             e_p(k,j,i) = ( surf_lsm_h%us(m) / c_0 )**2
          ENDDO
!
!--       Urban surfaces
          DO  m = 1, surf_usm_h%ns
             i = surf_usm_h%i(m)
             j = surf_usm_h%j(m)
             k = surf_usm_h%k(m)
             e_p(k,j,i) = ( surf_usm_h%us(m) / c_0 )**2
          ENDDO
!
!--       Vertical surfaces
          DO  l = 0, 3
!
!--          Default surfaces
             DO  m = 1, surf_def_v(l)%ns
                i = surf_def_v(l)%i(m)
                j = surf_def_v(l)%j(m)
                k = surf_def_v(l)%k(m)
                e_p(k,j,i) = ( surf_def_v(l)%us(m) / c_0 )**2
             ENDDO
!
!--          Natural surfaces
             DO  m = 1, surf_lsm_v(l)%ns
                i = surf_lsm_v(l)%i(m)
                j = surf_lsm_v(l)%j(m)
                k = surf_lsm_v(l)%k(m)
                e_p(k,j,i) = ( surf_lsm_v(l)%us(m) / c_0 )**2
             ENDDO
!
!--          Urban surfaces
             DO  m = 1, surf_usm_v(l)%ns
                i = surf_usm_v(l)%i(m)
                j = surf_usm_v(l)%j(m)
                k = surf_usm_v(l)%k(m)
                e_p(k,j,i) = ( surf_usm_v(l)%us(m) / c_0 )**2
             ENDDO
          ENDDO
       ENDIF

       IF ( .NOT. child_domain )  THEN
          e_p(nzt+1,:,:) = e_p(nzt,:,:)
!
!--    Nesting case: if parent operates in RANS mode and child in LES mode,
!--    no TKE is transfered. This case, set Neumann conditions at lateral and 
!--    top child boundaries. 
!--    If not ( both either in RANS or in LES mode ), TKE boundary condition
!--    is treated in the nesting. 
       ELSE

          IF ( rans_mode_parent  .AND.  .NOT. rans_mode )  THEN

             e_p(nzt+1,:,:) = e_p(nzt,:,:)
             IF ( bc_dirichlet_l )  e_p(:,:,nxl-1) = e_p(:,:,nxl)
             IF ( bc_dirichlet_r )  e_p(:,:,nxr+1) = e_p(:,:,nxr)
             IF ( bc_dirichlet_s )  e_p(:,nys-1,:) = e_p(:,nys,:)
             IF ( bc_dirichlet_n )  e_p(:,nyn+1,:) = e_p(:,nyn,:)

          ENDIF
       ENDIF
    ENDIF

!
!-- Boundary conditions for TKE dissipation rate. 
    IF ( rans_tke_e )  THEN
!
!--    Use wall function within constant-flux layer
!--    Upward-facing surfaces
!--    Default surfaces
       DO  m = 1, surf_def_h(0)%ns
          i = surf_def_h(0)%i(m)
          j = surf_def_h(0)%j(m)
          k = surf_def_h(0)%k(m)
          diss_p(k,j,i) = surf_def_h(0)%us(m)**3          &
                        / ( kappa * surf_def_h(0)%z_mo(m) )
       ENDDO
!
!--    Natural surfaces
       DO  m = 1, surf_lsm_h%ns
          i = surf_lsm_h%i(m)
          j = surf_lsm_h%j(m)
          k = surf_lsm_h%k(m)
          diss_p(k,j,i) = surf_lsm_h%us(m)**3          &
                        / ( kappa * surf_lsm_h%z_mo(m) )
       ENDDO
!
!--    Urban surfaces
       DO  m = 1, surf_usm_h%ns
          i = surf_usm_h%i(m)
          j = surf_usm_h%j(m)
          k = surf_usm_h%k(m)
          diss_p(k,j,i) = surf_usm_h%us(m)**3          &
                        / ( kappa * surf_usm_h%z_mo(m) )
       ENDDO
!
!--    Vertical surfaces
       DO  l = 0, 3
!
!--       Default surfaces
          DO  m = 1, surf_def_v(l)%ns
             i = surf_def_v(l)%i(m)
             j = surf_def_v(l)%j(m)
             k = surf_def_v(l)%k(m)
             diss_p(k,j,i) = surf_def_v(l)%us(m)**3          &
                           / ( kappa * surf_def_v(l)%z_mo(m) )
          ENDDO
!
!--       Natural surfaces
          DO  m = 1, surf_lsm_v(l)%ns
             i = surf_lsm_v(l)%i(m)
             j = surf_lsm_v(l)%j(m)
             k = surf_lsm_v(l)%k(m)
             diss_p(k,j,i) = surf_lsm_v(l)%us(m)**3          &
                           / ( kappa * surf_lsm_v(l)%z_mo(m) )
          ENDDO
!
!--       Urban surfaces
          DO  m = 1, surf_usm_v(l)%ns
             i = surf_usm_v(l)%i(m)
             j = surf_usm_v(l)%j(m)
             k = surf_usm_v(l)%k(m)
             diss_p(k,j,i) = surf_usm_v(l)%us(m)**3          &
                           / ( kappa * surf_usm_v(l)%z_mo(m) )
          ENDDO
       ENDDO
!
!--    Limit change of diss to be between -90% and +100%. Also, set an absolute
!--    minimum value
       DO  i = nxl, nxr
          DO  j = nys, nyn
             DO  k = nzb, nzt+1
                diss_p(k,j,i) = MAX( MIN( diss_p(k,j,i),          &
                                          2.0_wp * diss(k,j,i) ), &
                                     0.1_wp * diss(k,j,i),        &
                                     0.0001_wp )
             ENDDO
          ENDDO
       ENDDO

       IF ( .NOT. child_domain )  THEN
          diss_p(nzt+1,:,:) = diss_p(nzt,:,:)
       ENDIF
    ENDIF

!
!-- Boundary conditions for salinity
    IF ( ocean_mode )  THEN
!
!--    Bottom boundary: Neumann condition because salinity flux is always
!--    given.
       DO  l = 0, 1
!
!--       Set kb, for upward-facing surfaces value at topography top (k-1) is set,
!--       for downward-facing surfaces at topography bottom (k+1). 
          kb = MERGE( -1, 1, l == 0 )
          !$OMP PARALLEL DO PRIVATE( i, j, k )
          DO  m = 1, bc_h(l)%ns
             i = bc_h(l)%i(m)            
             j = bc_h(l)%j(m)
             k = bc_h(l)%k(m)
             sa_p(k+kb,j,i) = sa_p(k,j,i)
          ENDDO
       ENDDO
!
!--    Top boundary: Dirichlet or Neumann
       IF ( ibc_sa_t == 0 )  THEN
           sa_p(nzt+1,:,:) = sa(nzt+1,:,:)
       ELSEIF ( ibc_sa_t == 1 )  THEN
           sa_p(nzt+1,:,:) = sa_p(nzt,:,:)
       ENDIF

    ENDIF

!
!-- Boundary conditions for total water content,
!-- bottom and top boundary (see also temperature)
    IF ( humidity )  THEN
!
!--    Surface conditions for constant_humidity_flux
!--    Run loop over all non-natural and natural walls. Note, in wall-datatype
!--    the k coordinate belongs to the atmospheric grid point, therefore, set
!--    q_p at k-1
       IF ( ibc_q_b == 0 ) THEN

          DO  l = 0, 1
!
!--          Set kb, for upward-facing surfaces value at topography top (k-1) is set,
!--          for downward-facing surfaces at topography bottom (k+1). 
             kb = MERGE( -1, 1, l == 0 )
             !$OMP PARALLEL DO PRIVATE( i, j, k )
             DO  m = 1, bc_h(l)%ns
                i = bc_h(l)%i(m)            
                j = bc_h(l)%j(m)
                k = bc_h(l)%k(m)
                q_p(k+kb,j,i) = q(k+kb,j,i)
             ENDDO
          ENDDO
          
       ELSE
          
          DO  l = 0, 1
!
!--          Set kb, for upward-facing surfaces value at topography top (k-1) is set,
!--          for downward-facing surfaces at topography bottom (k+1). 
             kb = MERGE( -1, 1, l == 0 )
             !$OMP PARALLEL DO PRIVATE( i, j, k )
             DO  m = 1, bc_h(l)%ns
                i = bc_h(l)%i(m)            
                j = bc_h(l)%j(m)
                k = bc_h(l)%k(m)
                q_p(k+kb,j,i) = q_p(k,j,i)
             ENDDO
          ENDDO
       ENDIF
!
!--    Top boundary
       IF ( ibc_q_t == 0 ) THEN
          q_p(nzt+1,:,:) = q(nzt+1,:,:)
       ELSEIF ( ibc_q_t == 1 ) THEN
          q_p(nzt+1,:,:) = q_p(nzt,:,:) + bc_q_t_val * dzu(nzt+1)
       ENDIF

       IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
!             
!--       Surface conditions cloud water (Dirichlet)
!--       Run loop over all non-natural and natural walls. Note, in wall-datatype
!--       the k coordinate belongs to the atmospheric grid point, therefore, set
!--       qr_p and nr_p at k-1
          !$OMP PARALLEL DO PRIVATE( i, j, k )
          DO  m = 1, bc_h(0)%ns
             i = bc_h(0)%i(m)            
             j = bc_h(0)%j(m)
             k = bc_h(0)%k(m)
             qc_p(k-1,j,i) = 0.0_wp
             nc_p(k-1,j,i) = 0.0_wp
          ENDDO
!
!--       Top boundary condition for cloud water (Dirichlet)
          qc_p(nzt+1,:,:) = 0.0_wp
          nc_p(nzt+1,:,:) = 0.0_wp
           
       ENDIF

       IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
!             
!--       Surface conditions rain water (Dirichlet)
!--       Run loop over all non-natural and natural walls. Note, in wall-datatype
!--       the k coordinate belongs to the atmospheric grid point, therefore, set
!--       qr_p and nr_p at k-1
          !$OMP PARALLEL DO PRIVATE( i, j, k )
          DO  m = 1, bc_h(0)%ns
             i = bc_h(0)%i(m)            
             j = bc_h(0)%j(m)
             k = bc_h(0)%k(m)
             qr_p(k-1,j,i) = 0.0_wp
             nr_p(k-1,j,i) = 0.0_wp
          ENDDO
!
!--       Top boundary condition for rain water (Dirichlet)
          qr_p(nzt+1,:,:) = 0.0_wp
          nr_p(nzt+1,:,:) = 0.0_wp
           
       ENDIF
    ENDIF
!
!-- Boundary conditions for scalar,
!-- bottom and top boundary (see also temperature)
    IF ( passive_scalar )  THEN
!
!--    Surface conditions for constant_humidity_flux
!--    Run loop over all non-natural and natural walls. Note, in wall-datatype
!--    the k coordinate belongs to the atmospheric grid point, therefore, set
!--    s_p at k-1
       IF ( ibc_s_b == 0 ) THEN
          
          DO  l = 0, 1
!
!--          Set kb, for upward-facing surfaces value at topography top (k-1) is set,
!--          for downward-facing surfaces at topography bottom (k+1). 
             kb = MERGE( -1, 1, l == 0 )
             !$OMP PARALLEL DO PRIVATE( i, j, k )
             DO  m = 1, bc_h(l)%ns
                i = bc_h(l)%i(m)            
                j = bc_h(l)%j(m)
                k = bc_h(l)%k(m)
                s_p(k+kb,j,i) = s(k+kb,j,i)
             ENDDO
          ENDDO
          
       ELSE
          
          DO  l = 0, 1
!
!--          Set kb, for upward-facing surfaces value at topography top (k-1) is set,
!--          for downward-facing surfaces at topography bottom (k+1). 
             kb = MERGE( -1, 1, l == 0 )
             !$OMP PARALLEL DO PRIVATE( i, j, k )
             DO  m = 1, bc_h(l)%ns
                i = bc_h(l)%i(m)            
                j = bc_h(l)%j(m)
                k = bc_h(l)%k(m)
                s_p(k+kb,j,i) = s_p(k,j,i)
             ENDDO
          ENDDO
       ENDIF
!
!--    Top boundary condition
       IF ( ibc_s_t == 0 )  THEN
          s_p(nzt+1,:,:) = s(nzt+1,:,:)
       ELSEIF ( ibc_s_t == 1 )  THEN
          s_p(nzt+1,:,:) = s_p(nzt,:,:)
       ELSEIF ( ibc_s_t == 2 )  THEN
          s_p(nzt+1,:,:) = s_p(nzt,:,:) + bc_s_t_val * dzu(nzt+1)
       ENDIF

    ENDIF    
!
!-- Top/bottom boundary conditions for chemical species
    IF ( air_chemistry )  CALL chem_boundary_conds( 'set_bc_bottomtop' )
!
!-- In case of inflow or nest boundary at the south boundary the boundary for v
!-- is at nys and in case of inflow or nest boundary at the left boundary the
!-- boundary for u is at nxl. Since in prognostic_equations (cache optimized
!-- version) these levels are handled as a prognostic level, boundary values
!-- have to be restored here.
!-- For the SGS-TKE, Neumann boundary conditions are used at the inflow.
    IF ( bc_dirichlet_s )  THEN
       v_p(:,nys,:) = v_p(:,nys-1,:)
       IF ( .NOT. constant_diffusion ) e_p(:,nys-1,:) = e_p(:,nys,:)
    ELSEIF ( bc_dirichlet_n )  THEN
       IF ( .NOT. constant_diffusion ) e_p(:,nyn+1,:) = e_p(:,nyn,:)
    ELSEIF ( bc_dirichlet_l ) THEN
       u_p(:,:,nxl) = u_p(:,:,nxl-1)
       IF ( .NOT. constant_diffusion ) e_p(:,:,nxl-1) = e_p(:,:,nxl)
    ELSEIF ( bc_dirichlet_r )  THEN
       IF ( .NOT. constant_diffusion ) e_p(:,:,nxr+1) = e_p(:,:,nxr)
    ENDIF

!
!-- The same restoration for u at i=nxl and v at j=nys as above must be made
!-- in case of nest boundaries. This must not be done in case of vertical nesting 
!-- mode as in that case the lateral boundaries are actually cyclic. 
    IF ( nesting_mode /= 'vertical'  .OR.  nesting_offline )  THEN
       IF ( bc_dirichlet_s )  THEN
          v_p(:,nys,:) = v_p(:,nys-1,:)
       ENDIF
       IF ( bc_dirichlet_l )  THEN
          u_p(:,:,nxl) = u_p(:,:,nxl-1)
       ENDIF
    ENDIF

!
!-- Lateral boundary conditions for scalar quantities at the outflow
    IF ( bc_radiation_s )  THEN
       pt_p(:,nys-1,:)     = pt_p(:,nys,:)
       IF ( .NOT. constant_diffusion )  e_p(:,nys-1,:) = e_p(:,nys,:)
       IF ( rans_tke_e )  diss_p(:,nys-1,:) = diss_p(:,nys,:)
       IF ( humidity )  THEN
          q_p(:,nys-1,:) = q_p(:,nys,:)
          IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
             qc_p(:,nys-1,:) = qc_p(:,nys,:)
             nc_p(:,nys-1,:) = nc_p(:,nys,:)
          ENDIF
          IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
             qr_p(:,nys-1,:) = qr_p(:,nys,:)
             nr_p(:,nys-1,:) = nr_p(:,nys,:)
          ENDIF
       ENDIF
       IF ( passive_scalar )  s_p(:,nys-1,:) = s_p(:,nys,:)
    ELSEIF ( bc_radiation_n )  THEN
       pt_p(:,nyn+1,:)     = pt_p(:,nyn,:)
       IF ( .NOT. constant_diffusion )  e_p(:,nyn+1,:) = e_p(:,nyn,:)
       IF ( rans_tke_e )  diss_p(:,nyn+1,:) = diss_p(:,nyn,:)
       IF ( humidity )  THEN
          q_p(:,nyn+1,:) = q_p(:,nyn,:)
          IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
             qc_p(:,nyn+1,:) = qc_p(:,nyn,:)
             nc_p(:,nyn+1,:) = nc_p(:,nyn,:)
          ENDIF
          IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
             qr_p(:,nyn+1,:) = qr_p(:,nyn,:)
             nr_p(:,nyn+1,:) = nr_p(:,nyn,:)
          ENDIF
       ENDIF
       IF ( passive_scalar )  s_p(:,nyn+1,:) = s_p(:,nyn,:)
    ELSEIF ( bc_radiation_l )  THEN
       pt_p(:,:,nxl-1)     = pt_p(:,:,nxl)
       IF ( .NOT. constant_diffusion )  e_p(:,:,nxl-1) = e_p(:,:,nxl)
       IF ( rans_tke_e )  diss_p(:,:,nxl-1) = diss_p(:,:,nxl)
       IF ( humidity )  THEN
          q_p(:,:,nxl-1) = q_p(:,:,nxl)
          IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
             qc_p(:,:,nxl-1) = qc_p(:,:,nxl)
             nc_p(:,:,nxl-1) = nc_p(:,:,nxl)
          ENDIF
          IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
             qr_p(:,:,nxl-1) = qr_p(:,:,nxl)
             nr_p(:,:,nxl-1) = nr_p(:,:,nxl)
          ENDIF
       ENDIF
       IF ( passive_scalar )  s_p(:,:,nxl-1) = s_p(:,:,nxl)
    ELSEIF ( bc_radiation_r )  THEN
       pt_p(:,:,nxr+1)     = pt_p(:,:,nxr)
       IF ( .NOT. constant_diffusion )  e_p(:,:,nxr+1) = e_p(:,:,nxr)
       IF ( rans_tke_e )  diss_p(:,:,nxr+1) = diss_p(:,:,nxr)
       IF ( humidity )  THEN
          q_p(:,:,nxr+1) = q_p(:,:,nxr)
          IF ( bulk_cloud_model  .AND.  microphysics_morrison )  THEN
             qc_p(:,:,nxr+1) = qc_p(:,:,nxr)
             nc_p(:,:,nxr+1) = nc_p(:,:,nxr)
          ENDIF
          IF ( bulk_cloud_model  .AND.  microphysics_seifert )  THEN
             qr_p(:,:,nxr+1) = qr_p(:,:,nxr)
             nr_p(:,:,nxr+1) = nr_p(:,:,nxr)
          ENDIF
       ENDIF
       IF ( passive_scalar )  s_p(:,:,nxr+1) = s_p(:,:,nxr)
    ENDIF

!
!-- Lateral boundary conditions for chemical species
    IF ( air_chemistry )  CALL chem_boundary_conds( 'set_bc_lateral' )    

!
!-- Radiation boundary conditions for the velocities at the respective outflow.
!-- The phase velocity is either assumed to the maximum phase velocity that
!-- ensures numerical stability (CFL-condition) or calculated after
!-- Orlanski(1976) and averaged along the outflow boundary.
    IF ( bc_radiation_s )  THEN

       IF ( use_cmax )  THEN
          u_p(:,-1,:) = u(:,0,:)
          v_p(:,0,:)  = v(:,1,:)
          w_p(:,-1,:) = w(:,0,:)          
       ELSEIF ( .NOT. use_cmax )  THEN

          c_max = dy / dt_3d

          c_u_m_l = 0.0_wp 
          c_v_m_l = 0.0_wp
          c_w_m_l = 0.0_wp

          c_u_m = 0.0_wp 
          c_v_m = 0.0_wp
          c_w_m = 0.0_wp

!
!--       Calculate the phase speeds for u, v, and w, first local and then
!--       average along the outflow boundary.
          DO  k = nzb+1, nzt+1
             DO  i = nxl, nxr

                denom = u_m_s(k,0,i) - u_m_s(k,1,i)

                IF ( denom /= 0.0_wp )  THEN
                   c_u(k,i) = -c_max * ( u(k,0,i) - u_m_s(k,0,i) ) / ( denom * tsc(2) )
                   IF ( c_u(k,i) < 0.0_wp )  THEN
                      c_u(k,i) = 0.0_wp
                   ELSEIF ( c_u(k,i) > c_max )  THEN
                      c_u(k,i) = c_max
                   ENDIF
                ELSE
                   c_u(k,i) = c_max
                ENDIF

                denom = v_m_s(k,1,i) - v_m_s(k,2,i)

                IF ( denom /= 0.0_wp )  THEN
                   c_v(k,i) = -c_max * ( v(k,1,i) - v_m_s(k,1,i) ) / ( denom * tsc(2) )
                   IF ( c_v(k,i) < 0.0_wp )  THEN
                      c_v(k,i) = 0.0_wp
                   ELSEIF ( c_v(k,i) > c_max )  THEN
                      c_v(k,i) = c_max
                   ENDIF
                ELSE
                   c_v(k,i) = c_max
                ENDIF

                denom = w_m_s(k,0,i) - w_m_s(k,1,i)

                IF ( denom /= 0.0_wp )  THEN
                   c_w(k,i) = -c_max * ( w(k,0,i) - w_m_s(k,0,i) ) / ( denom * tsc(2) )
                   IF ( c_w(k,i) < 0.0_wp )  THEN
                      c_w(k,i) = 0.0_wp
                   ELSEIF ( c_w(k,i) > c_max )  THEN
                      c_w(k,i) = c_max
                   ENDIF
                ELSE
                   c_w(k,i) = c_max
                ENDIF

                c_u_m_l(k) = c_u_m_l(k) + c_u(k,i)
                c_v_m_l(k) = c_v_m_l(k) + c_v(k,i)
                c_w_m_l(k) = c_w_m_l(k) + c_w(k,i)

             ENDDO
          ENDDO

#if defined( __parallel )   
          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
          CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
                              MPI_SUM, comm1dx, ierr )   
          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
          CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
                              MPI_SUM, comm1dx, ierr )  
          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
          CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
                              MPI_SUM, comm1dx, ierr )  
#else
          c_u_m = c_u_m_l
          c_v_m = c_v_m_l
          c_w_m = c_w_m_l
#endif

          c_u_m = c_u_m / (nx+1)
          c_v_m = c_v_m / (nx+1)
          c_w_m = c_w_m / (nx+1)

!
!--       Save old timelevels for the next timestep
          IF ( intermediate_timestep_count == 1 )  THEN
             u_m_s(:,:,:) = u(:,0:1,:)
             v_m_s(:,:,:) = v(:,1:2,:)
             w_m_s(:,:,:) = w(:,0:1,:)
          ENDIF

!
!--       Calculate the new velocities
          DO  k = nzb+1, nzt+1
             DO  i = nxlg, nxrg
                u_p(k,-1,i) = u(k,-1,i) - dt_3d * tsc(2) * c_u_m(k) *          &
                                       ( u(k,-1,i) - u(k,0,i) ) * ddy

                v_p(k,0,i)  = v(k,0,i)  - dt_3d * tsc(2) * c_v_m(k) *          &
                                       ( v(k,0,i) - v(k,1,i) ) * ddy

                w_p(k,-1,i) = w(k,-1,i) - dt_3d * tsc(2) * c_w_m(k) *          &
                                       ( w(k,-1,i) - w(k,0,i) ) * ddy
             ENDDO
          ENDDO

!
!--       Bottom boundary at the outflow
          IF ( ibc_uv_b == 0 )  THEN
             u_p(nzb,-1,:) = 0.0_wp 
             v_p(nzb,0,:)  = 0.0_wp  
          ELSE                    
             u_p(nzb,-1,:) =  u_p(nzb+1,-1,:)
             v_p(nzb,0,:)  =  v_p(nzb+1,0,:)
          ENDIF
          w_p(nzb,-1,:) = 0.0_wp

!
!--       Top boundary at the outflow
          IF ( ibc_uv_t == 0 )  THEN
             u_p(nzt+1,-1,:) = u_init(nzt+1)
             v_p(nzt+1,0,:)  = v_init(nzt+1)
          ELSE
             u_p(nzt+1,-1,:) = u_p(nzt,-1,:)
             v_p(nzt+1,0,:)  = v_p(nzt,0,:)
          ENDIF
          w_p(nzt:nzt+1,-1,:) = 0.0_wp

       ENDIF

    ENDIF

    IF ( bc_radiation_n )  THEN

       IF ( use_cmax )  THEN
          u_p(:,ny+1,:) = u(:,ny,:)
          v_p(:,ny+1,:) = v(:,ny,:)
          w_p(:,ny+1,:) = w(:,ny,:)          
       ELSEIF ( .NOT. use_cmax )  THEN

          c_max = dy / dt_3d

          c_u_m_l = 0.0_wp 
          c_v_m_l = 0.0_wp
          c_w_m_l = 0.0_wp

          c_u_m = 0.0_wp 
          c_v_m = 0.0_wp
          c_w_m = 0.0_wp

!
!--       Calculate the phase speeds for u, v, and w, first local and then
!--       average along the outflow boundary.
          DO  k = nzb+1, nzt+1
             DO  i = nxl, nxr

                denom = u_m_n(k,ny,i) - u_m_n(k,ny-1,i)

                IF ( denom /= 0.0_wp )  THEN
                   c_u(k,i) = -c_max * ( u(k,ny,i) - u_m_n(k,ny,i) ) / ( denom * tsc(2) )
                   IF ( c_u(k,i) < 0.0_wp )  THEN
                      c_u(k,i) = 0.0_wp
                   ELSEIF ( c_u(k,i) > c_max )  THEN
                      c_u(k,i) = c_max
                   ENDIF
                ELSE
                   c_u(k,i) = c_max
                ENDIF

                denom = v_m_n(k,ny,i) - v_m_n(k,ny-1,i)

                IF ( denom /= 0.0_wp )  THEN
                   c_v(k,i) = -c_max * ( v(k,ny,i) - v_m_n(k,ny,i) ) / ( denom * tsc(2) )
                   IF ( c_v(k,i) < 0.0_wp )  THEN
                      c_v(k,i) = 0.0_wp
                   ELSEIF ( c_v(k,i) > c_max )  THEN
                      c_v(k,i) = c_max
                   ENDIF
                ELSE
                   c_v(k,i) = c_max
                ENDIF

                denom = w_m_n(k,ny,i) - w_m_n(k,ny-1,i)

                IF ( denom /= 0.0_wp )  THEN
                   c_w(k,i) = -c_max * ( w(k,ny,i) - w_m_n(k,ny,i) ) / ( denom * tsc(2) )
                   IF ( c_w(k,i) < 0.0_wp )  THEN
                      c_w(k,i) = 0.0_wp
                   ELSEIF ( c_w(k,i) > c_max )  THEN
                      c_w(k,i) = c_max
                   ENDIF
                ELSE
                   c_w(k,i) = c_max
                ENDIF

                c_u_m_l(k) = c_u_m_l(k) + c_u(k,i)
                c_v_m_l(k) = c_v_m_l(k) + c_v(k,i)
                c_w_m_l(k) = c_w_m_l(k) + c_w(k,i)

             ENDDO
          ENDDO

#if defined( __parallel )   
          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
          CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
                              MPI_SUM, comm1dx, ierr )   
          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
          CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
                              MPI_SUM, comm1dx, ierr )  
          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
          CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
                              MPI_SUM, comm1dx, ierr )  
#else
          c_u_m = c_u_m_l
          c_v_m = c_v_m_l
          c_w_m = c_w_m_l
#endif

          c_u_m = c_u_m / (nx+1)
          c_v_m = c_v_m / (nx+1)
          c_w_m = c_w_m / (nx+1)

!
!--       Save old timelevels for the next timestep
          IF ( intermediate_timestep_count == 1 )  THEN
                u_m_n(:,:,:) = u(:,ny-1:ny,:)
                v_m_n(:,:,:) = v(:,ny-1:ny,:)
                w_m_n(:,:,:) = w(:,ny-1:ny,:)
          ENDIF

!
!--       Calculate the new velocities
          DO  k = nzb+1, nzt+1
             DO  i = nxlg, nxrg
                u_p(k,ny+1,i) = u(k,ny+1,i) - dt_3d * tsc(2) * c_u_m(k) *      &
                                       ( u(k,ny+1,i) - u(k,ny,i) ) * ddy

                v_p(k,ny+1,i) = v(k,ny+1,i)  - dt_3d * tsc(2) * c_v_m(k) *     &
                                       ( v(k,ny+1,i) - v(k,ny,i) ) * ddy

                w_p(k,ny+1,i) = w(k,ny+1,i) - dt_3d * tsc(2) * c_w_m(k) *      &
                                       ( w(k,ny+1,i) - w(k,ny,i) ) * ddy
             ENDDO
          ENDDO

!
!--       Bottom boundary at the outflow
          IF ( ibc_uv_b == 0 )  THEN
             u_p(nzb,ny+1,:) = 0.0_wp
             v_p(nzb,ny+1,:) = 0.0_wp   
          ELSE                    
             u_p(nzb,ny+1,:) =  u_p(nzb+1,ny+1,:)
             v_p(nzb,ny+1,:) =  v_p(nzb+1,ny+1,:)
          ENDIF
          w_p(nzb,ny+1,:) = 0.0_wp

!
!--       Top boundary at the outflow
          IF ( ibc_uv_t == 0 )  THEN
             u_p(nzt+1,ny+1,:) = u_init(nzt+1)
             v_p(nzt+1,ny+1,:) = v_init(nzt+1)
          ELSE
             u_p(nzt+1,ny+1,:) = u_p(nzt,nyn+1,:)
             v_p(nzt+1,ny+1,:) = v_p(nzt,nyn+1,:)
          ENDIF
          w_p(nzt:nzt+1,ny+1,:) = 0.0_wp

       ENDIF

    ENDIF

    IF ( bc_radiation_l )  THEN

       IF ( use_cmax )  THEN
          u_p(:,:,0)  = u(:,:,1)
          v_p(:,:,-1) = v(:,:,0)
          w_p(:,:,-1) = w(:,:,0)          
       ELSEIF ( .NOT. use_cmax )  THEN

          c_max = dx / dt_3d

          c_u_m_l = 0.0_wp 
          c_v_m_l = 0.0_wp
          c_w_m_l = 0.0_wp

          c_u_m = 0.0_wp 
          c_v_m = 0.0_wp
          c_w_m = 0.0_wp

!
!--       Calculate the phase speeds for u, v, and w, first local and then
!--       average along the outflow boundary.
          DO  k = nzb+1, nzt+1
             DO  j = nys, nyn

                denom = u_m_l(k,j,1) - u_m_l(k,j,2)

                IF ( denom /= 0.0_wp )  THEN
                   c_u(k,j) = -c_max * ( u(k,j,1) - u_m_l(k,j,1) ) / ( denom * tsc(2) )
                   IF ( c_u(k,j) < 0.0_wp )  THEN
                      c_u(k,j) = 0.0_wp
                   ELSEIF ( c_u(k,j) > c_max )  THEN
                      c_u(k,j) = c_max
                   ENDIF
                ELSE
                   c_u(k,j) = c_max
                ENDIF

                denom = v_m_l(k,j,0) - v_m_l(k,j,1)

                IF ( denom /= 0.0_wp )  THEN
                   c_v(k,j) = -c_max * ( v(k,j,0) - v_m_l(k,j,0) ) / ( denom * tsc(2) )
                   IF ( c_v(k,j) < 0.0_wp )  THEN
                      c_v(k,j) = 0.0_wp
                   ELSEIF ( c_v(k,j) > c_max )  THEN
                      c_v(k,j) = c_max
                   ENDIF
                ELSE
                   c_v(k,j) = c_max
                ENDIF

                denom = w_m_l(k,j,0) - w_m_l(k,j,1)

                IF ( denom /= 0.0_wp )  THEN
                   c_w(k,j) = -c_max * ( w(k,j,0) - w_m_l(k,j,0) ) / ( denom * tsc(2) )
                   IF ( c_w(k,j) < 0.0_wp )  THEN
                      c_w(k,j) = 0.0_wp
                   ELSEIF ( c_w(k,j) > c_max )  THEN
                      c_w(k,j) = c_max
                   ENDIF
                ELSE
                   c_w(k,j) = c_max
                ENDIF

                c_u_m_l(k) = c_u_m_l(k) + c_u(k,j)
                c_v_m_l(k) = c_v_m_l(k) + c_v(k,j)
                c_w_m_l(k) = c_w_m_l(k) + c_w(k,j)

             ENDDO
          ENDDO

#if defined( __parallel )   
          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
          CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
                              MPI_SUM, comm1dy, ierr )   
          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
          CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
                              MPI_SUM, comm1dy, ierr )  
          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
          CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
                              MPI_SUM, comm1dy, ierr )  
#else
          c_u_m = c_u_m_l
          c_v_m = c_v_m_l
          c_w_m = c_w_m_l
#endif

          c_u_m = c_u_m / (ny+1)
          c_v_m = c_v_m / (ny+1)
          c_w_m = c_w_m / (ny+1)

!
!--       Save old timelevels for the next timestep
          IF ( intermediate_timestep_count == 1 )  THEN
                u_m_l(:,:,:) = u(:,:,1:2)
                v_m_l(:,:,:) = v(:,:,0:1)
                w_m_l(:,:,:) = w(:,:,0:1)
          ENDIF

!
!--       Calculate the new velocities
          DO  k = nzb+1, nzt+1
             DO  j = nysg, nyng
                u_p(k,j,0) = u(k,j,0) - dt_3d * tsc(2) * c_u_m(k) *            &
                                       ( u(k,j,0) - u(k,j,1) ) * ddx

                v_p(k,j,-1) = v(k,j,-1) - dt_3d * tsc(2) * c_v_m(k) *          &
                                       ( v(k,j,-1) - v(k,j,0) ) * ddx

                w_p(k,j,-1) = w(k,j,-1) - dt_3d * tsc(2) * c_w_m(k) *          &
                                       ( w(k,j,-1) - w(k,j,0) ) * ddx
             ENDDO
          ENDDO

!
!--       Bottom boundary at the outflow
          IF ( ibc_uv_b == 0 )  THEN
             u_p(nzb,:,0)  = 0.0_wp 
             v_p(nzb,:,-1) = 0.0_wp
          ELSE                    
             u_p(nzb,:,0)  =  u_p(nzb+1,:,0)
             v_p(nzb,:,-1) =  v_p(nzb+1,:,-1)
          ENDIF
          w_p(nzb,:,-1) = 0.0_wp

!
!--       Top boundary at the outflow
          IF ( ibc_uv_t == 0 )  THEN
             u_p(nzt+1,:,0)  = u_init(nzt+1)
             v_p(nzt+1,:,-1) = v_init(nzt+1)
          ELSE
             u_p(nzt+1,:,0)  = u_p(nzt,:,0)
             v_p(nzt+1,:,-1) = v_p(nzt,:,-1)
          ENDIF
          w_p(nzt:nzt+1,:,-1) = 0.0_wp

       ENDIF

    ENDIF

    IF ( bc_radiation_r )  THEN

       IF ( use_cmax )  THEN
          u_p(:,:,nx+1) = u(:,:,nx)
          v_p(:,:,nx+1) = v(:,:,nx)
          w_p(:,:,nx+1) = w(:,:,nx)          
       ELSEIF ( .NOT. use_cmax )  THEN

          c_max = dx / dt_3d

          c_u_m_l = 0.0_wp 
          c_v_m_l = 0.0_wp
          c_w_m_l = 0.0_wp

          c_u_m = 0.0_wp 
          c_v_m = 0.0_wp
          c_w_m = 0.0_wp

!
!--       Calculate the phase speeds for u, v, and w, first local and then
!--       average along the outflow boundary.
          DO  k = nzb+1, nzt+1
             DO  j = nys, nyn

                denom = u_m_r(k,j,nx) - u_m_r(k,j,nx-1)

                IF ( denom /= 0.0_wp )  THEN
                   c_u(k,j) = -c_max * ( u(k,j,nx) - u_m_r(k,j,nx) ) / ( denom * tsc(2) )
                   IF ( c_u(k,j) < 0.0_wp )  THEN
                      c_u(k,j) = 0.0_wp
                   ELSEIF ( c_u(k,j) > c_max )  THEN
                      c_u(k,j) = c_max
                   ENDIF
                ELSE
                   c_u(k,j) = c_max
                ENDIF

                denom = v_m_r(k,j,nx) - v_m_r(k,j,nx-1)

                IF ( denom /= 0.0_wp )  THEN
                   c_v(k,j) = -c_max * ( v(k,j,nx) - v_m_r(k,j,nx) ) / ( denom * tsc(2) )
                   IF ( c_v(k,j) < 0.0_wp )  THEN
                      c_v(k,j) = 0.0_wp
                   ELSEIF ( c_v(k,j) > c_max )  THEN
                      c_v(k,j) = c_max
                   ENDIF
                ELSE
                   c_v(k,j) = c_max
                ENDIF

                denom = w_m_r(k,j,nx) - w_m_r(k,j,nx-1)

                IF ( denom /= 0.0_wp )  THEN
                   c_w(k,j) = -c_max * ( w(k,j,nx) - w_m_r(k,j,nx) ) / ( denom * tsc(2) )
                   IF ( c_w(k,j) < 0.0_wp )  THEN
                      c_w(k,j) = 0.0_wp
                   ELSEIF ( c_w(k,j) > c_max )  THEN
                      c_w(k,j) = c_max
                   ENDIF
                ELSE
                   c_w(k,j) = c_max
                ENDIF

                c_u_m_l(k) = c_u_m_l(k) + c_u(k,j)
                c_v_m_l(k) = c_v_m_l(k) + c_v(k,j)
                c_w_m_l(k) = c_w_m_l(k) + c_w(k,j)

             ENDDO
          ENDDO

#if defined( __parallel )   
          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
          CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
                              MPI_SUM, comm1dy, ierr )   
          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
          CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
                              MPI_SUM, comm1dy, ierr )  
          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
          CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
                              MPI_SUM, comm1dy, ierr )  
#else
          c_u_m = c_u_m_l
          c_v_m = c_v_m_l
          c_w_m = c_w_m_l
#endif

          c_u_m = c_u_m / (ny+1)
          c_v_m = c_v_m / (ny+1)
          c_w_m = c_w_m / (ny+1)

!
!--       Save old timelevels for the next timestep
          IF ( intermediate_timestep_count == 1 )  THEN
                u_m_r(:,:,:) = u(:,:,nx-1:nx)
                v_m_r(:,:,:) = v(:,:,nx-1:nx)
                w_m_r(:,:,:) = w(:,:,nx-1:nx)
          ENDIF

!
!--       Calculate the new velocities
          DO  k = nzb+1, nzt+1
             DO  j = nysg, nyng
                u_p(k,j,nx+1) = u(k,j,nx+1) - dt_3d * tsc(2) * c_u_m(k) *      &
                                       ( u(k,j,nx+1) - u(k,j,nx) ) * ddx

                v_p(k,j,nx+1) = v(k,j,nx+1) - dt_3d * tsc(2) * c_v_m(k) *      &
                                       ( v(k,j,nx+1) - v(k,j,nx) ) * ddx

                w_p(k,j,nx+1) = w(k,j,nx+1) - dt_3d * tsc(2) * c_w_m(k) *      &
                                       ( w(k,j,nx+1) - w(k,j,nx) ) * ddx
             ENDDO
          ENDDO

!
!--       Bottom boundary at the outflow
          IF ( ibc_uv_b == 0 )  THEN
             u_p(nzb,:,nx+1) = 0.0_wp
             v_p(nzb,:,nx+1) = 0.0_wp 
          ELSE                    
             u_p(nzb,:,nx+1) =  u_p(nzb+1,:,nx+1)
             v_p(nzb,:,nx+1) =  v_p(nzb+1,:,nx+1)
          ENDIF
          w_p(nzb,:,nx+1) = 0.0_wp

!
!--       Top boundary at the outflow
          IF ( ibc_uv_t == 0 )  THEN
             u_p(nzt+1,:,nx+1) = u_init(nzt+1)
             v_p(nzt+1,:,nx+1) = v_init(nzt+1)
          ELSE
             u_p(nzt+1,:,nx+1) = u_p(nzt,:,nx+1)
             v_p(nzt+1,:,nx+1) = v_p(nzt,:,nx+1)
          ENDIF
          w_p(nzt:nzt+1,:,nx+1) = 0.0_wp

       ENDIF

    ENDIF
    
    IF ( salsa )  THEN
       CALL salsa_boundary_conds
    ENDIF 

 END SUBROUTINE boundary_conds