Home

Context Navigation

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

Last change on this file since 75 was 75, checked in by raasch, 17 years ago
preliminary update for changes concerning non-cyclic boundary conditions
Property svn:keywords set to `Id`
File size: 17.1 KB

Line
1	SUBROUTINE boundary_conds( range )
2
3	!------------------------------------------------------------------------------!
4	! Actual revisions:
5	! -----------------
6	! The "main" part sets conditions for time level t+dt instead of level t,
7	! outflow boundary conditions changed from Neumann to radiation condition,
8	! uxrp, vynp eliminated, moisture renamed humidity
9	!
10	! Former revisions:
11	! -----------------
12	! $Id: boundary_conds.f90 75 2007-03-22 09:54:05Z raasch $
13	!
14	! 19 2007-02-23 04:53:48Z raasch
15	! Boundary conditions for e(nzt), pt(nzt), and q(nzt) removed because these
16	! gridpoints are now calculated by the prognostic equation,
17	! Dirichlet and zero gradient condition for pt established at top boundary
18	!
19	! RCS Log replace by Id keyword, revision history cleaned up
20	!
21	! Revision 1.15 2006/02/23 09:54:55 raasch
22	! Surface boundary conditions in case of topography: nzb replaced by
23	! 2d-k-index-arrays (nzb_w_inner, etc.). Conditions for u and v remain
24	! unchanged (still using nzb) because a non-flat topography must use a
25	! Prandtl-layer, which don't requires explicit setting of the surface values.
26	!
27	! Revision 1.1 1997/09/12 06:21:34 raasch
28	! Initial revision
29	!
30	!
31	! Description:
32	! ------------
33	! Boundary conditions for the prognostic quantities (range='main').
34	! In case of non-cyclic lateral boundaries the conditions for velocities at
35	! the outflow are set after the pressure solver has been called (range=
36	! 'outflow_uvw').
37	! One additional bottom boundary condition is applied for the TKE (=(u)*2)
38	! in prandtl_fluxes. The cyclic lateral boundary conditions are implicitly
39	! handled in routine exchange_horiz. Pressure boundary conditions are
40	! explicitly set in routines pres, poisfft, poismg and sor.
41	!------------------------------------------------------------------------------!
42
43	USE arrays_3d
44	USE control_parameters
45	USE grid_variables
46	USE indices
47	USE pegrid
48
49	IMPLICIT NONE
50
51	CHARACTER (LEN=*) :: range
52
53	INTEGER :: i, j, k
54
55	REAL :: c_max, c_u, c_v, c_w, denom
56
57
58	IF ( range == 'main') THEN
59	!
60	!-- Bottom boundary
61	IF ( ibc_uv_b == 0 ) THEN
62	!
63	!-- Satisfying the Dirichlet condition with an extra layer below the
64	!-- surface where the u and v component change their sign
65	u_p(nzb,:,:) = -u_p(nzb+1,:,:)
66	v_p(nzb,:,:) = -v_p(nzb+1,:,:)
67	ELSE
68	u_p(nzb,:,:) = u_p(nzb+1,:,:)
69	v_p(nzb,:,:) = v_p(nzb+1,:,:)
70	ENDIF
71	DO i = nxl-1, nxr+1
72	DO j = nys-1, nyn+1
73	w_p(nzb_w_inner(j,i),j,i) = 0.0
74	ENDDO
75	ENDDO
76
77	!
78	!-- Top boundary
79	IF ( ibc_uv_t == 0 ) THEN
80	u_p(nzt+1,:,:) = ug(nzt+1)
81	v_p(nzt+1,:,:) = vg(nzt+1)
82	ELSE
83	u_p(nzt+1,:,:) = u_p(nzt,:,:)
84	v_p(nzt+1,:,:) = v_p(nzt,:,:)
85	ENDIF
86	w_p(nzt:nzt+1,:,:) = 0.0 ! nzt is not a prognostic level (but cf. pres)
87
88	!
89	!-- Temperature at bottom boundary
90	IF ( ibc_pt_b == 0 ) THEN
91	DO i = nxl-1, nxr+1
92	DO j = nys-1, nyn+1
93	pt_p(nzb_s_inner(j,i),j,i) = pt(nzb_s_inner(j,i),j,i)
94	ENDDO
95	ENDDO
96	ELSE
97	DO i = nxl-1, nxr+1
98	DO j = nys-1, nyn+1
99	pt_p(nzb_s_inner(j,i),j,i) = pt_p(nzb_s_inner(j,i)+1,j,i)
100	ENDDO
101	ENDDO
102	ENDIF
103
104	!
105	!-- Temperature at top boundary
106	IF ( ibc_pt_t == 0 ) THEN
107	pt_p(nzt+1,:,:) = pt(nzt+1,:,:)
108	ELSEIF ( ibc_pt_t == 1 ) THEN
109	pt_p(nzt+1,:,:) = pt_p(nzt,:,:)
110	ELSEIF ( ibc_pt_t == 2 ) THEN
111	pt_p(nzt+1,:,:) = pt_p(nzt,:,:) + bc_pt_t_val * dzu(nzt+1)
112	ENDIF
113
114	!
115	!-- Boundary conditions for TKE
116	!-- Generally Neumann conditions with de/dz=0 are assumed
117	IF ( .NOT. constant_diffusion ) THEN
118	DO i = nxl-1, nxr+1
119	DO j = nys-1, nyn+1
120	e_p(nzb_s_inner(j,i),j,i) = e_p(nzb_s_inner(j,i)+1,j,i)
121	ENDDO
122	ENDDO
123	e_p(nzt+1,:,:) = e_p(nzt,:,:)
124	ENDIF
125
126	!
127	!-- Boundary conditions for total water content or scalar,
128	!-- bottom and surface boundary (see also temperature)
129	IF ( humidity .OR. passive_scalar ) THEN
130	!
131	!-- Surface conditions for constant_humidity_flux
132	IF ( ibc_q_b == 0 ) THEN
133	DO i = nxl-1, nxr+1
134	DO j = nys-1, nyn+1
135	q_p(nzb_s_inner(j,i),j,i) = q(nzb_s_inner(j,i),j,i)
136	ENDDO
137	ENDDO
138	ELSE
139	DO i = nxl-1, nxr+1
140	DO j = nys-1, nyn+1
141	q_p(nzb_s_inner(j,i),j,i) = q_p(nzb_s_inner(j,i)+1,j,i)
142	ENDDO
143	ENDDO
144	ENDIF
145	!
146	!-- Top boundary
147	q_p(nzt+1,:,:) = q_p(nzt,:,:) + bc_q_t_val * dzu(nzt+1)
148	ENDIF
149
150	!
151	!-- Lateral boundary conditions at the inflow. Quasi Neumann conditions
152	!-- are needed for the wall normal velocity in order to ensure zero
153	!-- divergence. Dirichlet conditions are used for all other quantities.
154	IF ( inflow_s ) THEN
155	v_p(:,nys,:) = v_p(:,nys-1,:)
156	ELSEIF ( inflow_n ) THEN
157	v_p(:,nyn,:) = v_p(:,nyn+1,:)
158	ELSEIF ( inflow_l ) THEN
159	u_p(:,:,nxl) = u_p(:,:,nxl-1)
160	ELSEIF ( inflow_r ) THEN
161	u_p(:,:,nxr) = u_p(:,:,nxr+1)
162	ENDIF
163
164	!
165	!-- Lateral boundary conditions for scalar quantities at the outflow
166	IF ( outflow_s ) THEN
167	pt_p(:,nys-1,:) = pt_p(:,nys,:)
168	IF ( .NOT. constant_diffusion ) e_p(:,nys-1,:) = e_p(:,nys,:)
169	IF ( humidity .OR. passive_scalar ) q_p(:,nys-1,:) = q_p(:,nys,:)
170	ELSEIF ( outflow_n ) THEN
171	pt_p(:,nyn+1,:) = pt_p(:,nyn,:)
172	IF ( .NOT. constant_diffusion ) e_p(:,nyn+1,:) = e_p(:,nyn,:)
173	IF ( humidity .OR. passive_scalar ) q_p(:,nyn+1,:) = q_p(:,nyn,:)
174	ELSEIF ( outflow_l ) THEN
175	pt_p(:,:,nxl-1) = pt_p(:,:,nxl)
176	IF ( .NOT. constant_diffusion ) e_p(:,:,nxl-1) = e_p(:,:,nxl)
177	IF ( humidity .OR. passive_scalar ) q_p(:,:,nxl-1) = q_p(:,:,nxl)
178	ELSEIF ( outflow_r ) THEN
179	pt_p(:,:,nxr+1) = pt_p(:,:,nxr)
180	IF ( .NOT. constant_diffusion ) e_p(:,:,nxr+1) = e_p(:,:,nxr)
181	IF ( humidity .OR. passive_scalar ) q_p(:,:,nxr+1) = q_p(:,:,nxr)
182	ENDIF
183
184	ENDIF
185
186	!
187	!-- Radiation boundary condition for the velocities at the respective outflow
188	IF ( outflow_s .AND. &
189	intermediate_timestep_count == intermediate_timestep_count_max ) &
190	THEN
191
192	c_max = dy / dt_3d
193
194	DO i = nxl-1, nxr+1
195	DO k = nzb+1, nzt+1
196
197	!
198	!-- First calculate the phase speeds for u,v, and w
199	denom = u_m_s(k,0,i) - u_m_s(k,1,i)
200
201	IF ( denom /= 0.0 ) THEN
202	c_u = -c_max * ( u(k,0,i) - u_m_s(k,0,i) ) / denom
203	IF ( c_u < 0.0 ) THEN
204	c_u = 0.0
205	ELSEIF ( c_u > c_max ) THEN
206	c_u = c_max
207	ENDIF
208	ELSE
209	c_u = c_max
210	ENDIF
211	denom = v_m_s(k,0,i) - v_m_s(k,1,i)
212
213	IF ( denom /= 0.0 ) THEN
214	c_v = -c_max * ( v(k,0,i) - v_m_s(k,0,i) ) / denom
215	IF ( c_v < 0.0 ) THEN
216	c_v = 0.0
217	ELSEIF ( c_v > c_max ) THEN
218	c_v = c_max
219	ENDIF
220	ELSE
221	c_v = c_max
222	ENDIF
223
224	denom = w_m_s(k,0,i) - w_m_s(k,1,i)
225
226	IF ( denom /= 0.0 ) THEN
227	c_w = -c_max * ( w(k,0,i) - w_m_s(k,0,i) ) / denom
228	IF ( c_w < 0.0 ) THEN
229	c_w = 0.0
230	ELSEIF ( c_w > c_max ) THEN
231	c_w = c_max
232	ENDIF
233	ELSE
234	c_w = c_max
235	ENDIF
236
237	!
238	!-- Calculate the new velocities
239	u_p(k,-1,i) = u(k,-1,i) + dt_3d * c_u * &
240	( u(k,-1,i) - u(k,0,i) ) * ddy
241
242	v_p(k,-1,i) = v(k,-1,i) + dt_3d * c_v * &
243	( v(k,-1,i) - v_m_s(k,0,i) ) * ddy
244
245	w_p(k,-1,i) = w(k,-1,i) + dt_3d * c_w * &
246	( w(k,-1,i) - w(k,0,i) ) * ddy
247
248	!
249	!-- Save old timelevels for the next timestep
250	u_m_s(k,:,i) = u(k,-1:1,i)
251	v_m_s(k,:,i) = v(k,-1:1,i)
252	w_m_s(k,:,i) = w(k,-1:1,i)
253
254	ENDDO
255	ENDDO
256
257	!
258	!-- Bottom boundary at the outflow
259	IF ( ibc_uv_b == 0 ) THEN
260	u_p(nzb,-1,:) = -u_p(nzb+1,-1,:)
261	v_p(nzb,-1,:) = -v_p(nzb+1,-1,:)
262	ELSE
263	u_p(nzb,-1,:) = u_p(nzb+1,-1,:)
264	v_p(nzb,-1,:) = v_p(nzb+1,-1,:)
265	ENDIF
266	w_p(nzb,ny+1,:) = 0.0
267
268	!
269	!-- Top boundary at the outflow
270	IF ( ibc_uv_t == 0 ) THEN
271	u_p(nzt+1,-1,:) = ug(nzt+1)
272	v_p(nzt+1,-1,:) = vg(nzt+1)
273	ELSE
274	u_p(nzt+1,-1,:) = u(nzt,-1,:)
275	v_p(nzt+1,-1,:) = v(nzt,-1,:)
276	ENDIF
277	w_p(nzt:nzt+1,-1,:) = 0.0
278
279	ENDIF
280
281	IF ( outflow_n .AND. &
282	intermediate_timestep_count == intermediate_timestep_count_max ) &
283	THEN
284
285	c_max = dy / dt_3d
286
287	DO i = nxl-1, nxr+1
288	DO k = nzb+1, nzt+1
289
290	!
291	!-- First calculate the phase speeds for u,v, and w
292	denom = u_m_n(k,ny,i) - u_m_n(k,ny-1,i)
293
294	IF ( denom /= 0.0 ) THEN
295	c_u = -c_max * ( u(k,ny,i) - u_m_n(k,ny,i) ) / denom
296	IF ( c_u < 0.0 ) THEN
297	c_u = 0.0
298	ELSEIF ( c_u > c_max ) THEN
299	c_u = c_max
300	ENDIF
301	ELSE
302	c_u = c_max
303	ENDIF
304
305	denom = v_m_n(k,ny,i) - v_m_n(k,ny-1,i)
306
307	IF ( denom /= 0.0 ) THEN
308	c_v = -c_max * ( v(k,ny,i) - v_m_n(k,ny,i) ) / denom
309	IF ( c_v < 0.0 ) THEN
310	c_v = 0.0
311	ELSEIF ( c_v > c_max ) THEN
312	c_v = c_max
313	ENDIF
314	ELSE
315	c_v = c_max
316	ENDIF
317
318	denom = w_m_n(k,ny,i) - w_m_n(k,ny-1,i)
319
320	IF ( denom /= 0.0 ) THEN
321	c_w = -c_max * ( w(k,ny,i) - w_m_n(k,ny,i) ) / denom
322	IF ( c_w < 0.0 ) THEN
323	c_w = 0.0
324	ELSEIF ( c_w > c_max ) THEN
325	c_w = c_max
326	ENDIF
327	ELSE
328	c_w = c_max
329	ENDIF
330
331	!
332	!-- Calculate the new velocities
333	u_p(k,ny+1,i) = u(k,ny+1,i) - dt_3d * c_u * &
334	( u(k,ny+1,i) - u(k,ny,i) ) * ddy
335
336	v_p(k,ny+1,i) = v(k,ny+1,i) - dt_3d * c_v * &
337	( v(k,ny+1,i) - v(k,ny,i) ) * ddy
338
339	w_p(k,ny+1,i) = w(k,ny+1,i) - dt_3d * c_w * &
340	( w(k,ny+1,i) - w(k,ny,i) ) * ddy
341
342	!
343	!-- Swap timelevels for the next timestep
344	u_m_n(k,:,i) = u(k,ny-1:ny+1,i)
345	v_m_n(k,:,i) = v(k,ny-1:ny+1,i)
346	w_m_n(k,:,i) = w(k,ny-1:ny+1,i)
347
348	ENDDO
349	ENDDO
350
351	!
352	!-- Bottom boundary at the outflow
353	IF ( ibc_uv_b == 0 ) THEN
354	u_p(nzb,ny+1,:) = -u_p(nzb+1,ny+1,:)
355	v_p(nzb,ny+1,:) = -v_p(nzb+1,ny+1,:)
356	ELSE
357	u_p(nzb,ny+1,:) = u_p(nzb+1,ny+1,:)
358	v_p(nzb,ny+1,:) = v_p(nzb+1,ny+1,:)
359	ENDIF
360	w_p(nzb,ny+1,:) = 0.0
361
362	!
363	!-- Top boundary at the outflow
364	IF ( ibc_uv_t == 0 ) THEN
365	u_p(nzt+1,ny+1,:) = ug(nzt+1)
366	v_p(nzt+1,ny+1,:) = vg(nzt+1)
367	ELSE
368	u_p(nzt+1,ny+1,:) = u_p(nzt,nyn+1,:)
369	v_p(nzt+1,ny+1,:) = v_p(nzt,nyn+1,:)
370	ENDIF
371	w_p(nzt:nzt+1,ny+1,:) = 0.0
372
373	ENDIF
374
375	IF ( outflow_l .AND. &
376	intermediate_timestep_count == intermediate_timestep_count_max ) &
377	THEN
378
379	c_max = dx / dt_3d
380
381	DO j = nys-1, nyn+1
382	DO k = nzb+1, nzt+1
383
384	!
385	!-- First calculate the phase speeds for u,v, and w
386	denom = u_m_l(k,j,0) - u_m_l(k,j,1)
387
388	IF ( denom /= 0.0 ) THEN
389	c_u = -c_max * ( u(k,j,0) - u_m_r(k,j,0) ) / denom
390	IF ( c_u > 0.0 ) THEN
391	c_u = 0.0
392	ELSEIF ( c_u < -c_max ) THEN
393	c_u = -c_max
394	ENDIF
395	ELSE
396	c_u = -c_max
397	ENDIF
398
399	denom = v_m_l(k,j,0) - v_m_l(k,j,1)
400
401	IF ( denom /= 0.0 ) THEN
402	c_v = -c_max * ( v(k,j,0) - v_m_l(k,j,0) ) / denom
403	IF ( c_v < 0.0 ) THEN
404	c_v = 0.0
405	ELSEIF ( c_v > c_max ) THEN
406	c_v = c_max
407	ENDIF
408	ELSE
409	c_v = c_max
410	ENDIF
411
412	denom = w_m_l(k,j,0) - w_m_l(k,j,1)
413
414	IF ( denom /= 0.0 ) THEN
415	c_w = -c_max * ( w(k,j,0) - w_m_l(k,j,0) ) / denom
416	IF ( c_w < 0.0 ) THEN
417	c_w = 0.0
418	ELSEIF ( c_w > c_max ) THEN
419	c_w = c_max
420	ENDIF
421	ELSE
422	c_w = c_max
423	ENDIF
424
425	!
426	!-- Calculate the new velocities
427	u_p(k,j,-1) = u(k,j,-1) + dt_3d * c_u * &
428	( u(k,j,-1) - u(k,j,0) ) * ddx
429
430	v_p(k,j,-1) = v(k,j,-1) + dt_3d * c_v * &
431	( v(k,j,-1) - v(k,j,0) ) * ddx
432
433	w_p(k,j,-1) = w(k,j,-1) + dt_3d * c_w * &
434	( w(k,j,-1) - w(k,j,0) ) * ddx
435
436	!
437	!-- Swap timelevels for the next timestep
438	u_m_l(k,j,:) = u(k,j,-1:1)
439	v_m_l(k,j,:) = v(k,j,-1:1)
440	w_m_l(k,j,:) = w(k,j,-1:1)
441
442	ENDDO
443	ENDDO
444
445	!
446	!-- Bottom boundary at the outflow
447	IF ( ibc_uv_b == 0 ) THEN
448	u_p(nzb,:,-1) = -u_p(nzb+1,:,-1)
449	v_p(nzb,:,-1) = -v_p(nzb+1,:,-1)
450	ELSE
451	u_p(nzb,:,-1) = u_p(nzb+1,:,-1)
452	v_p(nzb,:,-1) = v_p(nzb+1,:,-1)
453	ENDIF
454	w_p(nzb,:,-1) = 0.0
455
456	!
457	!-- Top boundary at the outflow
458	IF ( ibc_uv_t == 0 ) THEN
459	u_p(nzt+1,:,-1) = ug(nzt+1)
460	v_p(nzt+1,:,-1) = vg(nzt+1)
461	ELSE
462	u_p(nzt+1,:,-1) = u_p(nzt,:,-1)
463	v_p(nzt+1,:,-1) = v_p(nzt,:,-1)
464	ENDIF
465	w_p(nzt:nzt+1,:,-1) = 0.0
466
467	ENDIF
468
469	IF ( outflow_r .AND. &
470	intermediate_timestep_count == intermediate_timestep_count_max ) &
471	THEN
472
473	c_max = dx / dt_3d
474
475	DO j = nys-1, nyn+1
476	DO k = nzb+1, nzt+1
477
478	!
479	!-- First calculate the phase speeds for u,v, and w
480	denom = u_m_r(k,j,nx) - u_m_r(k,j,nx-1)
481
482	IF ( denom /= 0.0 ) THEN
483	c_u = -c_max * ( u(k,j,nx) - u_m_r(k,j,nx) ) / denom
484	IF ( c_u < 0.0 ) THEN
485	c_u = 0.0
486	ELSEIF ( c_u > c_max ) THEN
487	c_u = c_max
488	ENDIF
489	ELSE
490	c_u = c_max
491	ENDIF
492
493	denom = v_m_r(k,j,nx) - v_m_r(k,j,nx-1)
494
495	IF ( denom /= 0.0 ) THEN
496	c_v = -c_max * ( v(k,j,nx) - v_m_r(k,j,nx) ) / denom
497	IF ( c_v < 0.0 ) THEN
498	c_v = 0.0
499	ELSEIF ( c_v > c_max ) THEN
500	c_v = c_max
501	ENDIF
502	ELSE
503	c_v = c_max
504	ENDIF
505
506	denom = w_m_r(k,j,nx) - w_m_r(k,j,nx-1)
507
508	IF ( denom /= 0.0 ) THEN
509	c_w = -c_max * ( w(k,j,nx) - w_m_r(k,j,nx) ) / denom
510	IF ( c_w < 0.0 ) THEN
511	c_w = 0.0
512	ELSEIF ( c_w > c_max ) THEN
513	c_w = c_max
514	ENDIF
515	ELSE
516	c_w = c_max
517	ENDIF
518
519	!
520	!-- Calculate the new velocities
521	u_p(k,j,nx+1) = u(k,j,nx+1) - dt_3d * c_u * &
522	( u(k,j,nx+1) - u(k,j,nx) ) * ddx
523
524	v_p(k,j,nx+1) = v(k,j,nx+1) - dt_3d * c_v * &
525	( v(k,j,nx+1) - v(k,j,nx) ) * ddx
526
527	w_p(k,j,nx+1) = w(k,j,nx+1) - dt_3d * c_w * &
528	( w(k,j,nx+1) - w(k,j,nx) ) * ddx
529
530	!
531	!-- Swap timelevels for the next timestep
532	u_m_r(k,j,:) = u(k,j,nx-1:nx+1)
533	v_m_r(k,j,:) = v(k,j,nx-1:nx+1)
534	w_m_r(k,j,:) = w(k,j,nx-1:nx+1)
535
536	ENDDO
537	ENDDO
538
539	!
540	!-- Bottom boundary at the outflow
541	IF ( ibc_uv_b == 0 ) THEN
542	u_p(nzb,:,nx+1) = -u_p(nzb+1,:,nx+1)
543	v_p(nzb,:,nx+1) = -v_p(nzb+1,:,nx+1)
544	ELSE
545	u_p(nzb,:,nx+1) = u_p(nzb+1,:,nx+1)
546	v_p(nzb,:,nx+1) = v_p(nzb+1,:,nx+1)
547	ENDIF
548	w_p(nzb,:,nx+1) = 0.0
549
550	!
551	!-- Top boundary at the outflow
552	IF ( ibc_uv_t == 0 ) THEN
553	u_p(nzt+1,:,nx+1) = ug(nzt+1)
554	v_p(nzt+1,:,nx+1) = vg(nzt+1)
555	ELSE
556	u_p(nzt+1,:,nx+1) = u_p(nzt,:,nx+1)
557	v_p(nzt+1,:,nx+1) = v_p(nzt,:,nx+1)
558	ENDIF
559	w(nzt:nzt+1,:,nx+1) = 0.0
560
561	ENDIF
562
563
564	END SUBROUTINE boundary_conds

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |