source: palm/tags/release-3.1b/SOURCE/boundary_conds.f90 @ 256

Last change on this file since 256 was 4, checked in by raasch, 17 years ago

Id keyword set as property for all *.f90 files

  • Property svn:keywords set to Id
File size: 9.3 KB
Line 
1 SUBROUTINE boundary_conds( range )
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6!
7!
8! Former revisions:
9! -----------------
10! $Id: boundary_conds.f90 4 2007-02-13 11:33:16Z letzel $
11! RCS Log replace by Id keyword, revision history cleaned up
12!
13! Revision 1.15  2006/02/23 09:54:55  raasch
14! Surface boundary conditions in case of topography: nzb replaced by
15! 2d-k-index-arrays (nzb_w_inner, etc.). Conditions for u and v remain
16! unchanged (still using nzb) because a non-flat topography must use a
17! Prandtl-layer, which don't requires explicit setting of the surface values.
18!
19! Revision 1.1  1997/09/12 06:21:34  raasch
20! Initial revision
21!
22!
23! Description:
24! ------------
25! Boundary conditions for the prognostic quantities (range='main').
26! In case of non-cyclic lateral boundaries the conditions for velocities at
27! the outflow are set after the pressure solver has been called (range=
28! 'outflow_uvw').
29! One additional bottom boundary condition is applied for the TKE (=(u*)**2)
30! in prandtl_fluxes. The cyclic lateral boundary conditions are implicitly
31! handled in routine exchange_horiz. Pressure boundary conditions are
32! explicitly set in routines pres, poisfft, poismg and sor.
33!------------------------------------------------------------------------------!
34
35    USE arrays_3d
36    USE control_parameters
37    USE grid_variables
38    USE indices
39    USE pegrid
40
41    IMPLICIT NONE
42
43    CHARACTER (LEN=*) ::  range
44
45    INTEGER ::  i, j, k
46
47
48    IF ( range == 'main')  THEN
49!
50!--    Bottom boundary
51       IF ( ibc_uv_b == 0 )  THEN
52          u(nzb,:,:) = -u(nzb+1,:,:)  ! satisfying the Dirichlet condition with
53          v(nzb,:,:) = -v(nzb+1,:,:)  ! an extra layer below the surface where
54       ELSE                           ! the u and v component change their sign
55          u(nzb,:,:) = u(nzb+1,:,:)
56          v(nzb,:,:) = v(nzb+1,:,:)
57       ENDIF
58       DO  i = nxl-1, nxr+1
59          DO  j = nys-1, nyn+1
60             w(nzb_w_inner(j,i),j,i) = 0.0
61          ENDDO
62       ENDDO
63
64!
65!--    Top boundary
66       IF ( ibc_uv_t == 0 )  THEN
67          u(nzt+1,:,:) = ug(nzt+1)
68          v(nzt+1,:,:) = vg(nzt+1)
69       ELSE
70          u(nzt+1,:,:) = u(nzt,:,:)
71          v(nzt+1,:,:) = v(nzt,:,:)
72       ENDIF
73       w(nzt:nzt+1,:,:) = 0.0  ! nzt is not a prognostic level (but cf. pres)
74
75!
76!--    Temperature at bottom boundary
77       IF ( ibc_pt_b == 0 )  THEN
78          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
79             DO  i = nxl-1, nxr+1
80                DO  j = nys-1, nyn+1
81                   pt(nzb_s_inner(j,i),j,i) = pt_m(nzb_s_inner(j,i),j,i)
82                ENDDO
83             ENDDO
84          ELSE
85             DO  i = nxl-1, nxr+1
86                DO  j = nys-1, nyn+1
87                   pt(nzb_s_inner(j,i),j,i) = pt_p(nzb_s_inner(j,i),j,i) 
88                                            ! pt_m is not used for Runge-Kutta
89                ENDDO
90             ENDDO
91          ENDIF
92       ELSE
93          DO  i = nxl-1, nxr+1
94             DO  j = nys-1, nyn+1
95                pt(nzb_s_inner(j,i),j,i) = pt(nzb_s_inner(j,i)+1,j,i)
96             ENDDO
97          ENDDO
98       ENDIF
99
100!
101!--    Temperature at top boundary
102       IF ( ibc_pt_t == 1 )  THEN
103          pt(nzt,:,:)   = pt(nzt-1,:,:) + bc_pt_t_val * dzu(nzt)
104          pt(nzt+1,:,:) = pt(nzt,:,:)   + bc_pt_t_val * dzu(nzt+1)
105       ENDIF
106
107!
108!--    Boundary conditions for TKE
109!--    Generally Neumann conditions with de/dz=0 are assumed
110       IF ( .NOT. constant_diffusion )  THEN
111          DO  i = nxl-1, nxr+1
112             DO  j = nys-1, nyn+1
113                e(nzb_s_inner(j,i),j,i) = e(nzb_s_inner(j,i)+1,j,i)
114             ENDDO
115          ENDDO
116          e(nzt,:,:)   = e(nzt-1,:,:)
117          e(nzt+1,:,:) = e(nzt,:,:)
118       ENDIF
119
120!
121!--    Boundary conditions for total water content or scalar,
122!--    bottom and surface boundary (see also temperature)
123       IF ( moisture  .OR.  passive_scalar )  THEN
124!
125!--       Surface conditions for constant_moisture_flux
126          IF ( ibc_q_b == 0 ) THEN
127             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
128                DO  i = nxl-1, nxr+1
129                   DO  j = nys-1, nyn+1
130                      q(nzb_s_inner(j,i),j,i) = q_m(nzb_s_inner(j,i),j,i)
131                   ENDDO
132                ENDDO
133             ELSE
134                DO  i = nxl-1, nxr+1
135                   DO  j = nys-1, nyn+1
136                      q(nzb_s_inner(j,i),j,i) = q_p(nzb_s_inner(j,i),j,i)
137                   ENDDO                      ! q_m is not used for Runge-Kutta
138                ENDDO
139             ENDIF
140          ELSE
141             DO  i = nxl-1, nxr+1
142                DO  j = nys-1, nyn+1
143                   q(nzb_s_inner(j,i),j,i) = q(nzb_s_inner(j,i)+1,j,i)
144                ENDDO
145             ENDDO
146          ENDIF
147!
148!--       Top boundary
149          q(nzt,:,:)   = q(nzt-1,:,:) + bc_q_t_val * dzu(nzt)
150          q(nzt+1,:,:) = q(nzt,:,:)   + bc_q_t_val * dzu(nzt+1)
151       ENDIF
152
153!
154!--    Lateral boundary conditions at the inflow. Quasi Neumann conditions
155!--    are needed for the wall normal velocity in order to ensure zero
156!--    divergence. Dirichlet conditions are used for all other quantities.
157       IF ( inflow_s )  THEN
158          v(:,nys,:) = v(:,nys-1,:)
159       ELSEIF ( inflow_n ) THEN
160          v(:,nyn+vynp,:) = v(:,nyn+vynp+1,:)
161       ELSEIF ( inflow_l ) THEN
162          u(:,:,nxl) = u(:,:,nxl-1)
163       ELSEIF ( inflow_r ) THEN
164          u(:,:,nxr+uxrp) = u(:,:,nxr+uxrp+1)
165       ENDIF
166
167!
168!--    Lateral boundary conditions for scalar quantities at the outflow
169       IF ( outflow_s )  THEN
170          pt(:,nys-1,:)     = pt(:,nys,:)
171          IF ( .NOT. constant_diffusion     )  e(:,nys-1,:) = e(:,nys,:)
172          IF ( moisture .OR. passive_scalar )  q(:,nys-1,:) = q(:,nys,:)
173       ELSEIF ( outflow_n )  THEN
174          pt(:,nyn+1,:)     = pt(:,nyn,:)
175          IF ( .NOT. constant_diffusion     )  e(:,nyn+1,:) = e(:,nyn,:)
176          IF ( moisture .OR. passive_scalar )  q(:,nyn+1,:) = q(:,nyn,:)
177       ELSEIF ( outflow_l )  THEN
178          pt(:,:,nxl-1)     = pt(:,:,nxl)
179          IF ( .NOT. constant_diffusion     )  e(:,:,nxl-1) = e(:,:,nxl)
180          IF ( moisture .OR. passive_scalar )  q(:,:,nxl-1) = q(:,:,nxl)     
181       ELSEIF ( outflow_r )  THEN
182          pt(:,:,nxr+1)     = pt(:,:,nxr)
183          IF ( .NOT. constant_diffusion     )  e(:,:,nxr+1) = e(:,:,nxr)
184          IF ( moisture .OR. passive_scalar )  q(:,:,nxr+1) = q(:,:,nxr)
185       ENDIF
186
187    ENDIF
188
189    IF ( range == 'outflow_uvw' ) THEN
190!
191!--    Horizontal boundary conditions for the velocities at the outflow.
192!--    A Neumann condition is used for the wall normal velocity. The vertical
193!--    velocity is assumed as zero and a horizontal average along the wall is
194!--    used for the wall parallel horizontal velocity. The combination of all
195!--    three conditions ensures that the velocity field is free of divergence
196!--    at the outflow (uvmean_outflow_l is calculated in pres).
197       IF ( outflow_s ) THEN
198          v(:,nys-1,:) = v(:,nys,:)
199          w(:,nys-1,:) = 0.0
200!
201!--       Compute the mean horizontal wind parallel to and within the outflow
202!--       wall and use this as boundary condition for u
203#if defined( __parallel )
204          CALL MPI_ALLREDUCE( uvmean_outflow_l, uvmean_outflow, nzt-nzb+2, &
205                              MPI_REAL, MPI_SUM, comm1dx, ierr )
206          uvmean_outflow = uvmean_outflow / ( nx + 1.0 )
207#else
208          uvmean_outflow = uvmean_outflow_l / ( nx + 1.0 )
209#endif
210          DO  k = nzb, nzt+1
211             u(k,nys-1,:) = uvmean_outflow(k)
212          ENDDO
213       ENDIF
214
215       IF ( outflow_n ) THEN
216          v(:,nyn+vynp+1,:) = v(:,nyn+vynp,:)
217          w(:,nyn+1,:)      = 0.0
218!
219!--       Compute the mean horizontal wind parallel to and within the outflow
220!--       wall and use this as boundary condition for u
221#if defined( __parallel )
222          CALL MPI_ALLREDUCE( uvmean_outflow_l, uvmean_outflow, nzt-nzb+2, &
223                              MPI_REAL, MPI_SUM, comm1dx, ierr )
224          uvmean_outflow = uvmean_outflow / ( nx + 1.0 )
225#else
226          uvmean_outflow = uvmean_outflow_l / ( nx + 1.0 )
227#endif
228          DO  k = nzb, nzt+1
229             u(k,nyn+1,:) = uvmean_outflow(k)
230          ENDDO
231       ENDIF
232
233       IF ( outflow_l ) THEN
234          u(:,:,nxl-1) = u(:,:,nxl)
235          w(:,:,nxl-1) = 0.0
236!
237!--       Compute the mean horizontal wind parallel to and within the outflow
238!--       wall and use this as boundary condition for v
239#if defined( __parallel )
240          CALL MPI_ALLREDUCE( uvmean_outflow_l, uvmean_outflow, nzt-nzb+2, &
241                              MPI_REAL, MPI_SUM, comm1dy, ierr )
242          uvmean_outflow = uvmean_outflow / ( ny + 1.0 )
243#else
244          uvmean_outflow = uvmean_outflow_l / ( ny + 1.0 )
245#endif
246          DO  k = nzb, nzt+1
247             v(k,:,nxl-1) = uvmean_outflow(k)
248          ENDDO   
249
250       ENDIF
251
252       IF ( outflow_r ) THEN
253          u(:,:,nxr+uxrp+1) = u(:,:,nxr+uxrp)
254          w(:,:,nxr+1) = 0.0
255!
256!--       Compute the mean horizontal wind parallel to and within the outflow
257!--       wall and use this as boundary condition for v
258#if defined( __parallel )
259          CALL MPI_ALLREDUCE( uvmean_outflow_l, uvmean_outflow, nzt-nzb+2, &
260                              MPI_REAL, MPI_SUM, comm1dy, ierr )
261          uvmean_outflow = uvmean_outflow / ( ny + 1.0 )
262#else
263          uvmean_outflow = uvmean_outflow_l / ( ny + 1.0 )
264#endif
265          DO  k = nzb, nzt+1
266             v(k,:,nxr+1) = uvmean_outflow(k)
267          ENDDO   
268       ENDIF
269
270    ENDIF
271
272 
273 END SUBROUTINE boundary_conds
Note: See TracBrowser for help on using the repository browser.