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

 SUBROUTINE boundary_conds( range )

!------------------------------------------------------------------------------!
! Actual revisions:
! -----------------
! 
!
! Former revisions:
! -----------------
! $Log: boundary_conds.f90,v $
! Revision 1.15  2006/02/23 09:54:55  raasch
! Surface boundary conditions in case of topography: nzb replaced by
! 2d-k-index-arrays (nzb_w_inner, etc.). Conditions for u and v remain
! unchanged (still using nzb) because a non-flat topography must use a
! Prandtl-layer, which don't requires explicit setting of the surface values.
!
! Revision 1.14  2005/06/29 09:49:40  steinfeld
! Baroclinicity considered in the Dirichlet boundary condition for u and v at
! the top boundary
!
! Revision 1.13  2005/04/23 08:35:30  raasch
! Error removed in Dirichlet bottom boundary conditions for pt and q in case
! of Runge-Kutta schemes. So far, new surface values have been taken from
! timelevel t-dt, which does not exist in case of Runge-Kutta schemes.
!
! Revision 1.12  2005/03/26 14:58:24  raasch
! Non-cyclic boundary conditions included, argument range added
!
! Revision 1.11  2001/03/30 06:54:18  raasch
! Timelevel t+dt replaced by timelevel t,
! Translation of remaining German identifiers (variables, subroutines, etc.)
!
! Revision 1.10  2001/01/29 12:19:27  raasch
! Passive scalar is considered
!
! Revision 1.9  2001/01/22 05:25:56  raasch
! Module test_variables removed
!
! Revision 1.8  2000/07/04 14:07:59  raasch
! Missing diriclet boundary conditions for temperature and total water
! content added
!
! Revision 1.7  2000/04/13 13:56:05  schroeter
! Boundaray conditions for total water content
!
! Revision 1.6  2000/01/20  10:44:48  10:44:48  letzel (Marcus Letzel)
! All comments translated into English
! 
! Revision 1.5  2000/01/10 09:28:53  raasch
! Variablenuebergabe jetzt per Modul
! Randbedingungen fuer w im Rahmen der pointer-Einfuehrung
!
! Revision 1.4  1998/07/06 12:07:05  raasch
! + USE test_variables
!
! Revision 1.3  1998/03/10 07:19:37  raasch
! Beschreibung ergaenzt um Randbedingung fuer TKE
!
! Revision 1.2  1997/09/19 07:39:03  raasch
! Randbedingungen fuer TKE
!
! Revision 1.1  1997/09/12 06:21:34  raasch
! Initial revision
!
!
! Description:
! ------------
! Boundary conditions for the prognostic quantities (range='main').
! In case of non-cyclic lateral boundaries the conditions for velocities at
! the outflow are set after the pressure solver has been called (range=
! 'outflow_uvw').
! 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.
!------------------------------------------------------------------------------!

    USE arrays_3d
    USE control_parameters
    USE grid_variables
    USE indices
    USE pegrid

    IMPLICIT NONE

    CHARACTER (LEN=*) ::  range

    INTEGER ::  i, j, k


    IF ( range == 'main')  THEN
!
!--    Bottom boundary
       IF ( ibc_uv_b == 0 )  THEN
          u(nzb,:,:) = -u(nzb+1,:,:)  ! satisfying the Dirichlet condition with
          v(nzb,:,:) = -v(nzb+1,:,:)  ! an extra layer below the surface where
       ELSE                           ! the u and v component change their sign
          u(nzb,:,:) = u(nzb+1,:,:)
          v(nzb,:,:) = v(nzb+1,:,:)
       ENDIF
       DO  i = nxl-1, nxr+1
          DO  j = nys-1, nyn+1
             w(nzb_w_inner(j,i),j,i) = 0.0
          ENDDO
       ENDDO

!
!--    Top boundary
       IF ( ibc_uv_t == 0 )  THEN
          u(nzt+1,:,:) = ug(nzt+1)
          v(nzt+1,:,:) = vg(nzt+1)
       ELSE
          u(nzt+1,:,:) = u(nzt,:,:)
          v(nzt+1,:,:) = v(nzt,:,:)
       ENDIF
       w(nzt:nzt+1,:,:) = 0.0  ! nzt is not a prognostic level (but cf. pres)

!
!--    Temperature at bottom boundary
       IF ( ibc_pt_b == 0 )  THEN
          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
             DO  i = nxl-1, nxr+1
                DO  j = nys-1, nyn+1
                   pt(nzb_s_inner(j,i),j,i) = pt_m(nzb_s_inner(j,i),j,i)
                ENDDO
             ENDDO
          ELSE
             DO  i = nxl-1, nxr+1
                DO  j = nys-1, nyn+1
                   pt(nzb_s_inner(j,i),j,i) = pt_p(nzb_s_inner(j,i),j,i)  
                                            ! pt_m is not used for Runge-Kutta
                ENDDO
             ENDDO
          ENDIF
       ELSE
          DO  i = nxl-1, nxr+1
             DO  j = nys-1, nyn+1
                pt(nzb_s_inner(j,i),j,i) = pt(nzb_s_inner(j,i)+1,j,i)
             ENDDO
          ENDDO
       ENDIF

!
!--    Temperature at top boundary
       IF ( ibc_pt_t == 1 )  THEN
          pt(nzt,:,:)   = pt(nzt-1,:,:) + bc_pt_t_val * dzu(nzt)
          pt(nzt+1,:,:) = pt(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
          DO  i = nxl-1, nxr+1
             DO  j = nys-1, nyn+1
                e(nzb_s_inner(j,i),j,i) = e(nzb_s_inner(j,i)+1,j,i)
             ENDDO
          ENDDO
          e(nzt,:,:)   = e(nzt-1,:,:)
          e(nzt+1,:,:) = e(nzt,:,:)
       ENDIF

!
!--    Boundary conditions for total water content or scalar,
!--    bottom and surface boundary (see also temperature)
       IF ( moisture  .OR.  passive_scalar )  THEN
!
!--       Surface conditions for constant_moisture_flux
          IF ( ibc_q_b == 0 ) THEN
             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
                DO  i = nxl-1, nxr+1
                   DO  j = nys-1, nyn+1
                      q(nzb_s_inner(j,i),j,i) = q_m(nzb_s_inner(j,i),j,i)
                   ENDDO
                ENDDO
             ELSE
                DO  i = nxl-1, nxr+1
                   DO  j = nys-1, nyn+1
                      q(nzb_s_inner(j,i),j,i) = q_p(nzb_s_inner(j,i),j,i)
                   ENDDO                      ! q_m is not used for Runge-Kutta
                ENDDO
             ENDIF
          ELSE
             DO  i = nxl-1, nxr+1
                DO  j = nys-1, nyn+1
                   q(nzb_s_inner(j,i),j,i) = q(nzb_s_inner(j,i)+1,j,i)
                ENDDO
             ENDDO
          ENDIF
!
!--       Top boundary
          q(nzt,:,:)   = q(nzt-1,:,:) + bc_q_t_val * dzu(nzt)
          q(nzt+1,:,:) = q(nzt,:,:)   + bc_q_t_val * dzu(nzt+1)
       ENDIF

!
!--    Lateral boundary conditions at the inflow. Quasi Neumann conditions
!--    are needed for the wall normal velocity in order to ensure zero
!--    divergence. Dirichlet conditions are used for all other quantities.
       IF ( inflow_s )  THEN
          v(:,nys,:) = v(:,nys-1,:)
       ELSEIF ( inflow_n ) THEN
          v(:,nyn+vynp,:) = v(:,nyn+vynp+1,:)
       ELSEIF ( inflow_l ) THEN
          u(:,:,nxl) = u(:,:,nxl-1)
       ELSEIF ( inflow_r ) THEN
          u(:,:,nxr+uxrp) = u(:,:,nxr+uxrp+1)
       ENDIF

!
!--    Lateral boundary conditions for scalar quantities at the outflow
       IF ( outflow_s )  THEN
          pt(:,nys-1,:)     = pt(:,nys,:)
          IF ( .NOT. constant_diffusion     )  e(:,nys-1,:) = e(:,nys,:)
          IF ( moisture .OR. passive_scalar )  q(:,nys-1,:) = q(:,nys,:)
       ELSEIF ( outflow_n )  THEN
          pt(:,nyn+1,:)     = pt(:,nyn,:)
          IF ( .NOT. constant_diffusion     )  e(:,nyn+1,:) = e(:,nyn,:)
          IF ( moisture .OR. passive_scalar )  q(:,nyn+1,:) = q(:,nyn,:)
       ELSEIF ( outflow_l )  THEN
          pt(:,:,nxl-1)     = pt(:,:,nxl)
          IF ( .NOT. constant_diffusion     )  e(:,:,nxl-1) = e(:,:,nxl)
          IF ( moisture .OR. passive_scalar )  q(:,:,nxl-1) = q(:,:,nxl)      
       ELSEIF ( outflow_r )  THEN
          pt(:,:,nxr+1)     = pt(:,:,nxr)
          IF ( .NOT. constant_diffusion     )  e(:,:,nxr+1) = e(:,:,nxr)
          IF ( moisture .OR. passive_scalar )  q(:,:,nxr+1) = q(:,:,nxr)
       ENDIF

    ENDIF

    IF ( range == 'outflow_uvw' ) THEN
!
!--    Horizontal boundary conditions for the velocities at the outflow.
!--    A Neumann condition is used for the wall normal velocity. The vertical
!--    velocity is assumed as zero and a horizontal average along the wall is
!--    used for the wall parallel horizontal velocity. The combination of all
!--    three conditions ensures that the velocity field is free of divergence
!--    at the outflow (uvmean_outflow_l is calculated in pres).
       IF ( outflow_s ) THEN
          v(:,nys-1,:) = v(:,nys,:)
          w(:,nys-1,:) = 0.0
!
!--       Compute the mean horizontal wind parallel to and within the outflow
!--       wall and use this as boundary condition for u
#if defined( __parallel )
          CALL MPI_ALLREDUCE( uvmean_outflow_l, uvmean_outflow, nzt-nzb+2, &
                              MPI_REAL, MPI_SUM, comm1dx, ierr )
          uvmean_outflow = uvmean_outflow / ( nx + 1.0 )
#else
          uvmean_outflow = uvmean_outflow_l / ( nx + 1.0 )
#endif
          DO  k = nzb, nzt+1
             u(k,nys-1,:) = uvmean_outflow(k)
          ENDDO
       ENDIF

       IF ( outflow_n ) THEN
          v(:,nyn+vynp+1,:) = v(:,nyn+vynp,:)
          w(:,nyn+1,:)      = 0.0
!
!--       Compute the mean horizontal wind parallel to and within the outflow
!--       wall and use this as boundary condition for u
#if defined( __parallel )
          CALL MPI_ALLREDUCE( uvmean_outflow_l, uvmean_outflow, nzt-nzb+2, &
                              MPI_REAL, MPI_SUM, comm1dx, ierr )
          uvmean_outflow = uvmean_outflow / ( nx + 1.0 )
#else
          uvmean_outflow = uvmean_outflow_l / ( nx + 1.0 )
#endif
          DO  k = nzb, nzt+1
             u(k,nyn+1,:) = uvmean_outflow(k)
          ENDDO
       ENDIF

       IF ( outflow_l ) THEN
          u(:,:,nxl-1) = u(:,:,nxl)
          w(:,:,nxl-1) = 0.0
!
!--       Compute the mean horizontal wind parallel to and within the outflow
!--       wall and use this as boundary condition for v
#if defined( __parallel )
          CALL MPI_ALLREDUCE( uvmean_outflow_l, uvmean_outflow, nzt-nzb+2, &
                              MPI_REAL, MPI_SUM, comm1dy, ierr )
          uvmean_outflow = uvmean_outflow / ( ny + 1.0 )
#else
          uvmean_outflow = uvmean_outflow_l / ( ny + 1.0 )
#endif
          DO  k = nzb, nzt+1
             v(k,:,nxl-1) = uvmean_outflow(k)
          ENDDO    

       ENDIF

       IF ( outflow_r ) THEN
          u(:,:,nxr+uxrp+1) = u(:,:,nxr+uxrp)
          w(:,:,nxr+1) = 0.0
!
!--       Compute the mean horizontal wind parallel to and within the outflow
!--       wall and use this as boundary condition for v
#if defined( __parallel )
          CALL MPI_ALLREDUCE( uvmean_outflow_l, uvmean_outflow, nzt-nzb+2, &
                              MPI_REAL, MPI_SUM, comm1dy, ierr )
          uvmean_outflow = uvmean_outflow / ( ny + 1.0 )
#else
          uvmean_outflow = uvmean_outflow_l / ( ny + 1.0 )
#endif
          DO  k = nzb, nzt+1
             v(k,:,nxr+1) = uvmean_outflow(k)
          ENDDO    
       ENDIF

    ENDIF

 
 END SUBROUTINE boundary_conds