SUBROUTINE boundary_conds( range ) !------------------------------------------------------------------------------! ! Actual revisions: ! ----------------- ! ! ! Former revisions: ! ----------------- ! $Id: boundary_conds.f90 39 2007-03-01 12:46:59Z kanani $ ! ! 19 2007-02-23 04:53:48Z raasch ! Boundary conditions for e(nzt), pt(nzt), and q(nzt) removed because these ! gridpoints are now calculated by the prognostic equation, ! Dirichlet and zero gradient condition for pt established at top boundary ! ! RCS Log replace by Id keyword, revision history cleaned up ! ! 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.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 == 0 ) THEN IF ( timestep_scheme(1:5) /= 'runge' ) THEN pt(nzt+1,:,:) = pt_m(nzt+1,:,:) ELSE pt(nzt+1,:,:) = pt_p(nzt+1,:,:) ! pt_m not used for Runge-Kutta ENDIF ELSEIF ( ibc_pt_t == 1 ) THEN pt(nzt+1,:,:) = pt(nzt,:,:) ELSEIF ( ibc_pt_t == 2 ) THEN 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+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+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