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

!> @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-2019 Leibniz Universitaet Hannover
!------------------------------------------------------------------------------!
!
! Current revisions:
! -----------------
! 
! 
! Former revisions:
! -----------------
! $Id: boundary_conds.f90 4268 2019-10-17 11:29:38Z schwenkel $
! Removing bulk cloud variables to respective module
! 
! 4182 2019-08-22 15:20:23Z scharf
! Corrected "Former revisions" section
! 
! 4102 2019-07-17 16:00:03Z suehring
! - Revise setting for boundary conditions at horizontal walls, use the offset
!   index that belongs to the data structure instead of pre-calculating 
!   the offset index for each facing.
! - Set boundary conditions for bulk-cloud quantities also at downward facing
!   surfaces
! 
! 4087 2019-07-11 11:35:23Z gronemeier
! Add comment
! 
! 4086 2019-07-11 05:55:44Z gronemeier
! Bugfix: use constant-flux layer condition for e in all rans modes
! 
! 3879 2019-04-08 20:25:23Z knoop
! Bugfix, do not set boundary conditions for potential temperature in neutral
! runs.
! 
! 3655 2019-01-07 16:51:22Z knoop
! OpenACC port for SPEC
!
! 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,  &
               dzu, pt, pt_init, pt_p, q,                                      &
               q_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 chemistry_model_mod,                                                   &
        ONLY:  chem_boundary_conds

    USE control_parameters,                                                    &
        ONLY:  air_chemistry, bc_dirichlet_l,                                  &
               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, 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, neutral, nudging, ocean_mode, passive_scalar,  &
               tsc, salsa, 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
 
    USE salsa_mod,                                                             &
        ONLY:  salsa_boundary_conds

    USE surface_mod,                                                           &
        ONLY :  bc_h

    USE turbulence_closure_mod,                                                &
        ONLY:  tcm_boundary_conds

    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) ::  l  !< running index boundary type, for up- and downward-facing walls
    INTEGER(iwp) ::  m  !< running index surface elements

    REAL(wp)    ::  c_max !< maximum phase velocity allowed by CFL criterion, used for outflow boundary condition
    REAL(wp)    ::  denom !< horizontal gradient of velocity component normal to the outflow boundary 


!
!-- 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
       !$OMP PARALLEL DO PRIVATE( i, j, k )
       !$ACC PARALLEL LOOP PRIVATE(i, j, k) &
       !$ACC PRESENT(bc_h, w_p)
       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+bc_h(l)%koff,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
        !$ACC KERNELS PRESENT(u_p, v_p, u_init, v_init)
        u_p(nzt+1,:,:) = u_init(nzt+1)
        v_p(nzt+1,:,:) = v_init(nzt+1)
        !$ACC END KERNELS
    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
       !$ACC KERNELS PRESENT(w_p)
       w_p(nzt:nzt+1,:,:) = 0.0_wp  !< nzt is not a prognostic level (but cf. pres)
       !$ACC END KERNELS
    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 ( .NOT. neutral )  THEN
       IF ( ibc_pt_b == 0 )  THEN
          DO  l = 0, 1
             !$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+bc_h(l)%koff,j,i) = pt(k+bc_h(l)%koff,j,i)
             ENDDO
          ENDDO
!      
!--    Neumann, zero-gradient
       ELSEIF ( ibc_pt_b == 1 )  THEN
          DO  l = 0, 1
             !$OMP PARALLEL DO PRIVATE( i, j, k )
             !$ACC PARALLEL LOOP PRIVATE(i, j, k) &
             !$ACC PRESENT(bc_h, pt_p)
             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+bc_h(l)%koff,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
           !$ACC KERNELS PRESENT(pt_p, dzu)
           pt_p(nzt+1,:,:) = pt_p(nzt,:,:) + bc_pt_t_val * dzu(nzt+1)
           !$ACC END KERNELS
       ENDIF
    ENDIF
!
!-- Boundary conditions for salinity
    IF ( ocean_mode )  THEN
!
!--    Bottom boundary: Neumann condition because salinity flux is always
!--    given.
       DO  l = 0, 1
          !$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+bc_h(l)%koff,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
             !$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+bc_h(l)%koff,j,i) = q(k+bc_h(l)%koff,j,i)
             ENDDO
          ENDDO

       ELSE

          DO  l = 0, 1
             !$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+bc_h(l)%koff,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
    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
             !$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+bc_h(l)%koff,j,i) = s(k+bc_h(l)%koff,j,i)
             ENDDO
          ENDDO

       ELSE

          DO  l = 0, 1
             !$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+bc_h(l)%koff,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  
!
!-- Set boundary conditions for subgrid TKE and dissipation (RANS mode)
    CALL tcm_boundary_conds
!
!-- 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.
    IF ( bc_dirichlet_s )  THEN
       v_p(:,nys,:) = v_p(:,nys-1,:)
    ELSEIF ( bc_dirichlet_l ) THEN
       u_p(:,:,nxl) = u_p(:,:,nxl-1)
    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. 
!-- Lateral oundary conditions for TKE and dissipation are set 
!-- in tcm_boundary_conds.
    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.
!-- Lateral oundary conditions for TKE and dissipation are set 
!-- in tcm_boundary_conds.
    IF ( bc_radiation_s )  THEN
       pt_p(:,nys-1,:)     = pt_p(:,nys,:)
       IF ( humidity )  THEN
          q_p(:,nys-1,:) = q_p(:,nys,:)
       ENDIF
       IF ( passive_scalar )  s_p(:,nys-1,:) = s_p(:,nys,:)
    ELSEIF ( bc_radiation_n )  THEN
       pt_p(:,nyn+1,:)     = pt_p(:,nyn,:)
       IF ( humidity )  THEN
          q_p(:,nyn+1,:) = q_p(:,nyn,:)
       ENDIF
       IF ( passive_scalar )  s_p(:,nyn+1,:) = s_p(:,nyn,:)
    ELSEIF ( bc_radiation_l )  THEN
       pt_p(:,:,nxl-1)     = pt_p(:,:,nxl)
       IF ( humidity )  THEN
          q_p(:,:,nxl-1) = q_p(:,:,nxl)
       ENDIF
       IF ( passive_scalar )  s_p(:,:,nxl-1) = s_p(:,:,nxl)
    ELSEIF ( bc_radiation_r )  THEN
       pt_p(:,:,nxr+1)     = pt_p(:,:,nxr)
       IF ( humidity )  THEN
          q_p(:,:,nxr+1) = q_p(:,:,nxr)
       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