Home

Context Navigation

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

Last change on this file since 73 was 73, checked in by raasch, 17 years ago
preliminary changes for radiation conditions
Property svn:keywords set to `Id`
File size: 13.8 KB

Line
1	SUBROUTINE boundary_conds( range )
2
3	!------------------------------------------------------------------------------!
4	! Actual revisions:
5	! -----------------
6	! The "main" part sets conditions for time level t+dt insteat of level t,
7	! outflow boundary conditions changed from Neumann to radiation condition
8	!
9	! Former revisions:
10	! -----------------
11	! $Id: boundary_conds.f90 73 2007-03-20 08:33:14Z raasch $
12	!
13	! 19 2007-02-23 04:53:48Z raasch
14	! Boundary conditions for e(nzt), pt(nzt), and q(nzt) removed because these
15	! gridpoints are now calculated by the prognostic equation,
16	! Dirichlet and zero gradient condition for pt established at top boundary
17	!
18	! RCS Log replace by Id keyword, revision history cleaned up
19	!
20	! Revision 1.15 2006/02/23 09:54:55 raasch
21	! Surface boundary conditions in case of topography: nzb replaced by
22	! 2d-k-index-arrays (nzb_w_inner, etc.). Conditions for u and v remain
23	! unchanged (still using nzb) because a non-flat topography must use a
24	! Prandtl-layer, which don't requires explicit setting of the surface values.
25	!
26	! Revision 1.1 1997/09/12 06:21:34 raasch
27	! Initial revision
28	!
29	!
30	! Description:
31	! ------------
32	! Boundary conditions for the prognostic quantities (range='main').
33	! In case of non-cyclic lateral boundaries the conditions for velocities at
34	! the outflow are set after the pressure solver has been called (range=
35	! 'outflow_uvw').
36	! One additional bottom boundary condition is applied for the TKE (=(u)*2)
37	! in prandtl_fluxes. The cyclic lateral boundary conditions are implicitly
38	! handled in routine exchange_horiz. Pressure boundary conditions are
39	! explicitly set in routines pres, poisfft, poismg and sor.
40	!------------------------------------------------------------------------------!
41
42	USE arrays_3d
43	USE control_parameters
44	USE grid_variables
45	USE indices
46	USE pegrid
47
48	IMPLICIT NONE
49
50	CHARACTER (LEN=*) :: range
51
52	INTEGER :: i, j, k
53
54	REAL :: c_max, c_u, c_v, c_w, denom
55
56
57	IF ( range == 'main') THEN
58	!
59	!-- Bottom boundary
60	IF ( ibc_uv_b == 0 ) THEN
61	!
62	!-- Satisfying the Dirichlet condition with an extra layer below the
63	!-- surface where the u and v component change their sign
64	u_p(nzb,:,:) = -u_p(nzb+1,:,:)
65	v_p(nzb,:,:) = -v_p(nzb+1,:,:)
66	ELSE
67	u_p(nzb,:,:) = u_p(nzb+1,:,:)
68	v_p(nzb,:,:) = v_p(nzb+1,:,:)
69	ENDIF
70	DO i = nxl-1, nxr+1
71	DO j = nys-1, nyn+1
72	w_p(nzb_w_inner(j,i),j,i) = 0.0
73	ENDDO
74	ENDDO
75
76	!
77	!-- Top boundary
78	IF ( ibc_uv_t == 0 ) THEN
79	u_p(nzt+1,:,:) = ug(nzt+1)
80	v_p(nzt+1,:,:) = vg(nzt+1)
81	ELSE
82	u_p(nzt+1,:,:) = u_p(nzt,:,:)
83	v_p(nzt+1,:,:) = v_p(nzt,:,:)
84	ENDIF
85	w_p(nzt:nzt+1,:,:) = 0.0 ! nzt is not a prognostic level (but cf. pres)
86
87	!
88	!-- Temperature at bottom boundary
89	IF ( ibc_pt_b == 0 ) THEN
90	DO i = nxl-1, nxr+1
91	DO j = nys-1, nyn+1
92	pt_p(nzb_s_inner(j,i),j,i) = pt(nzb_s_inner(j,i),j,i)
93	ENDDO
94	ENDDO
95	ELSE
96	DO i = nxl-1, nxr+1
97	DO j = nys-1, nyn+1
98	pt_p(nzb_s_inner(j,i),j,i) = pt_p(nzb_s_inner(j,i)+1,j,i)
99	ENDDO
100	ENDDO
101	ENDIF
102
103	!
104	!-- Temperature at top boundary
105	IF ( ibc_pt_t == 0 ) THEN
106	pt_p(nzt+1,:,:) = pt(nzt+1,:,:)
107	ELSEIF ( ibc_pt_t == 1 ) THEN
108	pt_p(nzt+1,:,:) = pt_p(nzt,:,:)
109	ELSEIF ( ibc_pt_t == 2 ) THEN
110	pt_p(nzt+1,:,:) = pt_p(nzt,:,:) + bc_pt_t_val * dzu(nzt+1)
111	ENDIF
112
113	!
114	!-- Boundary conditions for TKE
115	!-- Generally Neumann conditions with de/dz=0 are assumed
116	IF ( .NOT. constant_diffusion ) THEN
117	DO i = nxl-1, nxr+1
118	DO j = nys-1, nyn+1
119	e_p(nzb_s_inner(j,i),j,i) = e_p(nzb_s_inner(j,i)+1,j,i)
120	ENDDO
121	ENDDO
122	e_p(nzt+1,:,:) = e_p(nzt,:,:)
123	ENDIF
124
125	!
126	!-- Boundary conditions for total water content or scalar,
127	!-- bottom and surface boundary (see also temperature)
128	IF ( moisture .OR. passive_scalar ) THEN
129	!
130	!-- Surface conditions for constant_moisture_flux
131	IF ( ibc_q_b == 0 ) THEN
132	DO i = nxl-1, nxr+1
133	DO j = nys-1, nyn+1
134	q_p(nzb_s_inner(j,i),j,i) = q(nzb_s_inner(j,i),j,i)
135	ENDDO
136	ENDDO
137	ELSE
138	DO i = nxl-1, nxr+1
139	DO j = nys-1, nyn+1
140	q_p(nzb_s_inner(j,i),j,i) = q_p(nzb_s_inner(j,i)+1,j,i)
141	ENDDO
142	ENDDO
143	ENDIF
144	!
145	!-- Top boundary
146	q_p(nzt+1,:,:) = q_p(nzt,:,:) + bc_q_t_val * dzu(nzt+1)
147	ENDIF
148
149	!
150	!-- Lateral boundary conditions at the inflow. Quasi Neumann conditions
151	!-- are needed for the wall normal velocity in order to ensure zero
152	!-- divergence. Dirichlet conditions are used for all other quantities.
153	IF ( inflow_s ) THEN
154	v_p(:,nys,:) = v_p(:,nys-1,:)
155	ELSEIF ( inflow_n ) THEN
156	v_p(:,nyn+vynp,:) = v_p(:,nyn+vynp+1,:)
157	ELSEIF ( inflow_l ) THEN
158	u_p(:,:,nxl) = u_p(:,:,nxl-1)
159	ELSEIF ( inflow_r ) THEN
160	u_p(:,:,nxr+uxrp) = u_p(:,:,nxr+uxrp+1)
161	ENDIF
162
163	!
164	!-- Lateral boundary conditions for scalar quantities at the outflow
165	IF ( outflow_s ) THEN
166	pt_p(:,nys-1,:) = pt_p(:,nys,:)
167	IF ( .NOT. constant_diffusion ) e_p(:,nys-1,:) = e_p(:,nys,:)
168	IF ( moisture .OR. passive_scalar ) q_p(:,nys-1,:) = q_p(:,nys,:)
169	ELSEIF ( outflow_n ) THEN
170	pt_p(:,nyn+1,:) = pt_p(:,nyn,:)
171	IF ( .NOT. constant_diffusion ) e_p(:,nyn+1,:) = e_p(:,nyn,:)
172	IF ( moisture .OR. passive_scalar ) q_p(:,nyn+1,:) = q_p(:,nyn,:)
173	ELSEIF ( outflow_l ) THEN
174	pt_p(:,:,nxl-1) = pt_p(:,:,nxl)
175	IF ( .NOT. constant_diffusion ) e_p(:,:,nxl-1) = e_p(:,:,nxl)
176	IF ( moisture .OR. passive_scalar ) q_p(:,:,nxl-1) = q_p(:,:,nxl)
177	ELSEIF ( outflow_r ) THEN
178	pt_p(:,:,nxr+1) = pt_p(:,:,nxr)
179	IF ( .NOT. constant_diffusion ) e_p(:,:,nxr+1) = e_p(:,:,nxr)
180	IF ( moisture .OR. passive_scalar ) q_p(:,:,nxr+1) = q_p(:,:,nxr)
181	ENDIF
182
183	ENDIF
184
185	IF ( range == 'outflow_uvw' ) THEN
186	!
187	!-- Radiation boundary condition for the velocities at the respective outflow
188	IF ( outflow_s ) THEN
189	! v(:,nys-1,:) = v(:,nys,:)
190	! w(:,nys-1,:) = 0.0
191	!!
192	!!-- Compute the mean horizontal wind parallel to and within the outflow
193	!!-- wall and use this as boundary condition for u
194	!#if defined( __parallel )
195	! CALL MPI_ALLREDUCE( uvmean_outflow_l, uvmean_outflow, nzt-nzb+2, &
196	! MPI_REAL, MPI_SUM, comm1dx, ierr )
197	! uvmean_outflow = uvmean_outflow / ( nx + 1.0 )
198	!#else
199	! uvmean_outflow = uvmean_outflow_l / ( nx + 1.0 )
200	!#endif
201	! DO k = nzb, nzt+1
202	! u(k,nys-1,:) = uvmean_outflow(k)
203	! ENDDO
204	ENDIF
205
206	IF ( outflow_n .AND. &
207	intermediate_timestep_count == intermediate_timestep_count_max ) &
208	THEN
209
210	c_max = dy / dt_3d
211
212	DO i = nxl-1, nxr+1
213	DO k = nzb+1, nzt+1
214
215	!
216	!-- First calculate the phase speeds for u,v, and w
217	denom = u_m_n(k,ny,i,-2) - u_m_n(k,ny-1,i,-2)
218
219	IF ( denom /= 0.0 ) THEN
220	c_u = -c_max * ( u_m_n(k,ny,i,-1)-u_m_n(k,ny,i,-2) ) / denom
221	IF ( c_u < 0.0 ) THEN
222	c_u = 0.0
223	ELSEIF ( c_u > c_max ) THEN
224	c_u = c_max
225	ENDIF
226	ELSE
227	c_u = c_max
228	ENDIF
229
230	denom = v_m_n(k,ny,i,-2) - v_m_n(k,ny-1,i,-2)
231
232	IF ( denom /= 0.0 ) THEN
233	c_v = -c_max * ( v_m_n(k,ny,i,-1)-v_m_n(k,ny,i,-2) ) / denom
234	IF ( c_v < 0.0 ) THEN
235	c_v = 0.0
236	ELSEIF ( c_v > c_max ) THEN
237	c_v = c_max
238	ENDIF
239	ELSE
240	c_v = c_max
241	ENDIF
242
243	denom = w_m_n(k,ny,i,-2) - w_m_n(k,ny-1,i,-2)
244
245	IF ( denom /= 0.0 ) THEN
246	c_w = -c_max * ( w_m_n(k,ny,i,-1)-w_m_n(k,ny,i,-2) ) / denom
247	IF ( c_w < 0.0 ) THEN
248	c_w = 0.0
249	ELSEIF ( c_w > c_max ) THEN
250	c_w = c_max
251	ENDIF
252	ELSE
253	c_w = c_max
254	ENDIF
255
256	!
257	!-- Calculate the new velocities
258	u(k,ny+1,i) = u_m_n(k,ny+1,i,-1) - dt_3d * c_u * &
259	( u_m_n(k,ny+1,i,-1) - u_m_n(k,ny,i,-1) ) * ddy
260
261	v(k,ny+1,i) = v_m_n(k,ny+1,i,-1) - dt_3d * c_v * &
262	( v_m_n(k,ny+1,i,-1) - v_m_n(k,ny,i,-1) ) * ddy
263
264	w(k,ny+1,i) = w_m_n(k,ny+1,i,-1) - dt_3d * c_w * &
265	( w_m_n(k,ny+1,i,-1) - w_m_n(k,ny,i,-1) ) * ddy
266
267	!
268	!-- Swap timelevels for the next timestep
269	u_m_n(k,:,i,-2) = u_m_n(k,:,i,-1)
270	u_m_n(k,:,i,-1) = u(k,ny-1:ny+1,i)
271	v_m_n(k,:,i,-2) = v_m_n(k,:,i,-1)
272	v_m_n(k,:,i,-1) = v(k,ny-1:ny+1,i)
273	w_m_n(k,:,i,-2) = w_m_n(k,:,i,-1)
274	w_m_n(k,:,i,-1) = w(k,ny-1:ny+1,i)
275
276	ENDDO
277	ENDDO
278
279	!
280	!-- Bottom boundary at the outflow
281	IF ( ibc_uv_b == 0 ) THEN
282	u(nzb,ny+1,:) = -u(nzb+1,ny+1,:)
283	v(nzb,ny+1,:) = -v(nzb+1,ny+1,:)
284	ELSE
285	u(nzb,ny+1,:) = u(nzb+1,ny+1,:)
286	v(nzb,ny+1,:) = v(nzb+1,ny+1,:)
287	ENDIF
288	w(nzb,ny+1,:) = 0.0
289
290	!
291	!-- Top boundary at the outflow
292	IF ( ibc_uv_t == 0 ) THEN
293	u(nzt+1,ny+1,:) = ug(nzt+1)
294	v(nzt+1,ny+1,:) = vg(nzt+1)
295	ELSE
296	u(nzt+1,ny+1,:) = u(nzt,nyn+1,:)
297	v(nzt+1,ny+1,:) = v(nzt,nyn+1,:)
298	ENDIF
299	w(nzt:nzt+1,ny+1,:) = 0.0
300
301	ENDIF
302
303	IF ( outflow_l ) THEN
304	! u(:,:,nxl-1) = u(:,:,nxl)
305	! w(:,:,nxl-1) = 0.0
306	!
307	!-- Compute the mean horizontal wind parallel to and within the outflow
308	!-- wall and use this as boundary condition for v
309	!#if defined( __parallel )
310	! CALL MPI_ALLREDUCE( uvmean_outflow_l, uvmean_outflow, nzt-nzb+2, &
311	! MPI_REAL, MPI_SUM, comm1dy, ierr )
312	! uvmean_outflow = uvmean_outflow / ( ny + 1.0 )
313	!#else
314	! uvmean_outflow = uvmean_outflow_l / ( ny + 1.0 )
315	!#endif
316	! DO k = nzb, nzt+1
317	! v(k,:,nxl-1) = uvmean_outflow(k)
318	! ENDDO
319	!
320	ENDIF
321
322	IF ( outflow_r .AND. &
323	intermediate_timestep_count == intermediate_timestep_count_max ) &
324	THEN
325
326	c_max = dx / dt_3d
327
328	DO j = nys-1, nyn+1
329	DO k = nzb+1, nzt+1
330
331	!
332	!-- First calculate the phase speeds for u,v, and w
333	denom = u_m_r(k,j,nx,-2) - u_m_r(k,j,nx-1,-2)
334
335	IF ( denom /= 0.0 ) THEN
336	c_u = -c_max * ( u_m_r(k,j,nx,-1)-u_m_r(k,j,nx,-2) ) / denom
337	IF ( c_u < 0.0 ) THEN
338	c_u = 0.0
339	ELSEIF ( c_u > c_max ) THEN
340	c_u = c_max
341	ENDIF
342	ELSE
343	c_u = c_max
344	ENDIF
345
346	denom = v_m_r(k,j,nx,-2) - v_m_r(k,j,nx-1,-2)
347
348	IF ( denom /= 0.0 ) THEN
349	c_v = -c_max * ( v_m_r(k,j,nx,-1)-v_m_r(k,j,nx,-2) ) / denom
350	IF ( c_v < 0.0 ) THEN
351	c_v = 0.0
352	ELSEIF ( c_v > c_max ) THEN
353	c_v = c_max
354	ENDIF
355	ELSE
356	c_v = c_max
357	ENDIF
358
359	denom = w_m_r(k,j,nx,-2) - w_m_r(k,j,nx-1,-2)
360
361	IF ( denom /= 0.0 ) THEN
362	c_w = -c_max * ( w_m_r(k,j,nx,-1)-w_m_n(k,j,nx,-2) ) / denom
363	IF ( c_w < 0.0 ) THEN
364	c_w = 0.0
365	ELSEIF ( c_w > c_max ) THEN
366	c_w = c_max
367	ENDIF
368	ELSE
369	c_w = c_max
370	ENDIF
371
372	!
373	!-- Calculate the new velocities
374	u(k,j,nx+1) = u_m_r(k,j,nx+1,-1) - dt_3d * c_u * &
375	( u_m_r(k,j,nx+1,-1) - u_m_r(k,j,nx,-1) ) * ddx
376
377	v(k,j,nx+1) = v_m_r(k,j,nx+1,-1) - dt_3d * c_v * &
378	( v_m_r(k,j,nx+1,-1) - v_m_r(k,j,nx,-1) ) * ddx
379
380	w(k,j,nx+1) = w_m_r(k,j,nx+1,-1) - dt_3d * c_w * &
381	( w_m_r(k,j,nx+1,-1) - w_m_r(k,j,nx,-1) ) * ddx
382
383	!
384	!-- Swap timelevels for the next timestep
385	u_m_r(k,j,:,-2) = u_m_r(k,j,:,-1)
386	u_m_r(k,j,:,-1) = u(k,j,nx-1:nx+1)
387	v_m_r(k,j,:,-2) = v_m_r(k,j,:,-1)
388	v_m_r(k,j,:,-1) = v(k,j,nx-1:nx+1)
389	w_m_r(k,j,:,-2) = w_m_r(k,j,:,-1)
390	w_m_r(k,j,:,-1) = w(k,j,nx-1:nx+1)
391
392	ENDDO
393	ENDDO
394
395	!
396	!-- Bottom boundary at the outflow
397	IF ( ibc_uv_b == 0 ) THEN
398	u(nzb,ny+1,:) = -u(nzb+1,ny+1,:)
399	v(nzb,ny+1,:) = -v(nzb+1,ny+1,:)
400	ELSE
401	u(nzb,ny+1,:) = u(nzb+1,ny+1,:)
402	v(nzb,ny+1,:) = v(nzb+1,ny+1,:)
403	ENDIF
404	w(nzb,ny+1,:) = 0.0
405
406	!
407	!-- Top boundary at the outflow
408	IF ( ibc_uv_t == 0 ) THEN
409	u(nzt+1,ny+1,:) = ug(nzt+1)
410	v(nzt+1,ny+1,:) = vg(nzt+1)
411	ELSE
412	u(nzt+1,ny+1,:) = u(nzt,nyn+1,:)
413	v(nzt+1,ny+1,:) = v(nzt,nyn+1,:)
414	ENDIF
415	w(nzt:nzt+1,ny+1,:) = 0.0
416
417	ENDIF
418
419	ENDIF
420
421
422	END SUBROUTINE boundary_conds

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |