Home

Context Navigation

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

Last change on this file since 1070 was 1054, checked in by hoffmann, 12 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 29.2 KB

Line
1	SUBROUTINE boundary_conds( range )
2
3	!--------------------------------------------------------------------------------!
4	! This file is part of PALM.
5	!
6	! PALM is free software: you can redistribute it and/or modify it under the terms
7	! of the GNU General Public License as published by the Free Software Foundation,
8	! either version 3 of the License, or (at your option) any later version.
9	!
10	! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11	! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12	! A PARTICULAR PURPOSE. See the GNU General Public License for more details.
13	!
14	! You should have received a copy of the GNU General Public License along with
15	! PALM. If not, see <http://www.gnu.org/licenses/>.
16	!
17	! Copyright 1997-2012 Leibniz University Hannover
18	!--------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	! Former revisions:
24	! -----------------
25	! $Id: boundary_conds.f90 1054 2012-11-13 17:30:09Z maronga $
26	!
27	! 1053 2012-11-13 17:11:03Z hoffmann
28	! boundary conditions for the two new prognostic equations (nr, qr) of the
29	! two-moment cloud scheme
30	!
31	! 1036 2012-10-22 13:43:42Z raasch
32	! code put under GPL (PALM 3.9)
33	!
34	! 996 2012-09-07 10:41:47Z raasch
35	! little reformatting
36	!
37	! 978 2012-08-09 08:28:32Z fricke
38	! Neumann boudnary conditions are added at the inflow boundary for the SGS-TKE.
39	! Outflow boundary conditions for the velocity components can be set to Neumann
40	! conditions or to radiation conditions with a horizontal averaged phase
41	! velocity.
42	!
43	! 875 2012-04-02 15:35:15Z gryschka
44	! Bugfix in case of dirichlet inflow bc at the right or north boundary
45	!
46	! 767 2011-10-14 06:39:12Z raasch
47	! ug,vg replaced by u_init,v_init as the Dirichlet top boundary condition
48	!
49	! 667 2010-12-23 12:06:00Z suehring/gryschka
50	! nxl-1, nxr+1, nys-1, nyn+1 replaced by nxlg, nxrg, nysg, nyng
51	! Removed mirror boundary conditions for u and v at the bottom in case of
52	! ibc_uv_b == 0. Instead, dirichelt boundary conditions (u=v=0) are set
53	! in init_3d_model
54	!
55	! 107 2007-08-17 13:54:45Z raasch
56	! Boundary conditions for temperature adjusted for coupled runs,
57	! bugfixes for the radiation boundary conditions at the outflow: radiation
58	! conditions are used for every substep, phase speeds are calculated for the
59	! first Runge-Kutta substep only and then reused, several index values changed
60	!
61	! 95 2007-06-02 16:48:38Z raasch
62	! Boundary conditions for salinity added
63	!
64	! 75 2007-03-22 09:54:05Z raasch
65	! The "main" part sets conditions for time level t+dt instead of level t,
66	! outflow boundary conditions changed from Neumann to radiation condition,
67	! uxrp, vynp eliminated, moisture renamed humidity
68	!
69	! 19 2007-02-23 04:53:48Z raasch
70	! Boundary conditions for e(nzt), pt(nzt), and q(nzt) removed because these
71	! gridpoints are now calculated by the prognostic equation,
72	! Dirichlet and zero gradient condition for pt established at top boundary
73	!
74	! RCS Log replace by Id keyword, revision history cleaned up
75	!
76	! Revision 1.15 2006/02/23 09:54:55 raasch
77	! Surface boundary conditions in case of topography: nzb replaced by
78	! 2d-k-index-arrays (nzb_w_inner, etc.). Conditions for u and v remain
79	! unchanged (still using nzb) because a non-flat topography must use a
80	! Prandtl-layer, which don't requires explicit setting of the surface values.
81	!
82	! Revision 1.1 1997/09/12 06:21:34 raasch
83	! Initial revision
84	!
85	!
86	! Description:
87	! ------------
88	! Boundary conditions for the prognostic quantities (range='main').
89	! In case of non-cyclic lateral boundaries the conditions for velocities at
90	! the outflow are set after the pressure solver has been called (range=
91	! 'outflow_uvw').
92	! One additional bottom boundary condition is applied for the TKE (=(u)*2)
93	! in prandtl_fluxes. The cyclic lateral boundary conditions are implicitly
94	! handled in routine exchange_horiz. Pressure boundary conditions are
95	! explicitly set in routines pres, poisfft, poismg and sor.
96	!------------------------------------------------------------------------------!
97
98	USE arrays_3d
99	USE control_parameters
100	USE grid_variables
101	USE indices
102	USE pegrid
103
104	IMPLICIT NONE
105
106	CHARACTER (LEN=*) :: range
107
108	INTEGER :: i, j, k
109
110	REAL :: c_max, denom
111
112
113	IF ( range == 'main') THEN
114	!
115	!-- Bottom boundary
116	IF ( ibc_uv_b == 1 ) THEN
117	u_p(nzb,:,:) = u_p(nzb+1,:,:)
118	v_p(nzb,:,:) = v_p(nzb+1,:,:)
119	ENDIF
120	DO i = nxlg, nxrg
121	DO j = nysg, nyng
122	w_p(nzb_w_inner(j,i),j,i) = 0.0
123	ENDDO
124	ENDDO
125
126	!
127	!-- Top boundary
128	IF ( ibc_uv_t == 0 ) THEN
129	u_p(nzt+1,:,:) = u_init(nzt+1)
130	v_p(nzt+1,:,:) = v_init(nzt+1)
131	ELSE
132	u_p(nzt+1,:,:) = u_p(nzt,:,:)
133	v_p(nzt+1,:,:) = v_p(nzt,:,:)
134	ENDIF
135	w_p(nzt:nzt+1,:,:) = 0.0 ! nzt is not a prognostic level (but cf. pres)
136
137	!
138	!-- Temperature at bottom boundary.
139	!-- In case of coupled runs (ibc_pt_b = 2) the temperature is given by
140	!-- the sea surface temperature of the coupled ocean model.
141	IF ( ibc_pt_b == 0 ) THEN
142	DO i = nxlg, nxrg
143	DO j = nysg, nyng
144	pt_p(nzb_s_inner(j,i),j,i) = pt(nzb_s_inner(j,i),j,i)
145	ENDDO
146	ENDDO
147	ELSEIF ( ibc_pt_b == 1 ) THEN
148	DO i = nxlg, nxrg
149	DO j = nysg, nyng
150	pt_p(nzb_s_inner(j,i),j,i) = pt_p(nzb_s_inner(j,i)+1,j,i)
151	ENDDO
152	ENDDO
153	ENDIF
154
155	!
156	!-- Temperature at top boundary
157	IF ( ibc_pt_t == 0 ) THEN
158	pt_p(nzt+1,:,:) = pt(nzt+1,:,:)
159	ELSEIF ( ibc_pt_t == 1 ) THEN
160	pt_p(nzt+1,:,:) = pt_p(nzt,:,:)
161	ELSEIF ( ibc_pt_t == 2 ) THEN
162	pt_p(nzt+1,:,:) = pt_p(nzt,:,:) + bc_pt_t_val * dzu(nzt+1)
163	ENDIF
164
165	!
166	!-- Boundary conditions for TKE
167	!-- Generally Neumann conditions with de/dz=0 are assumed
168	IF ( .NOT. constant_diffusion ) THEN
169	DO i = nxlg, nxrg
170	DO j = nysg, nyng
171	e_p(nzb_s_inner(j,i),j,i) = e_p(nzb_s_inner(j,i)+1,j,i)
172	ENDDO
173	ENDDO
174	e_p(nzt+1,:,:) = e_p(nzt,:,:)
175	ENDIF
176
177	!
178	!-- Boundary conditions for salinity
179	IF ( ocean ) THEN
180	!
181	!-- Bottom boundary: Neumann condition because salinity flux is always
182	!-- given
183	DO i = nxlg, nxrg
184	DO j = nysg, nyng
185	sa_p(nzb_s_inner(j,i),j,i) = sa_p(nzb_s_inner(j,i)+1,j,i)
186	ENDDO
187	ENDDO
188
189	!
190	!-- Top boundary: Dirichlet or Neumann
191	IF ( ibc_sa_t == 0 ) THEN
192	sa_p(nzt+1,:,:) = sa(nzt+1,:,:)
193	ELSEIF ( ibc_sa_t == 1 ) THEN
194	sa_p(nzt+1,:,:) = sa_p(nzt,:,:)
195	ENDIF
196
197	ENDIF
198
199	!
200	!-- Boundary conditions for total water content or scalar,
201	!-- bottom and top boundary (see also temperature)
202	IF ( humidity .OR. passive_scalar ) THEN
203	!
204	!-- Surface conditions for constant_humidity_flux
205	IF ( ibc_q_b == 0 ) THEN
206	DO i = nxlg, nxrg
207	DO j = nysg, nyng
208	q_p(nzb_s_inner(j,i),j,i) = q(nzb_s_inner(j,i),j,i)
209	ENDDO
210	ENDDO
211	ELSE
212	DO i = nxlg, nxrg
213	DO j = nysg, nyng
214	q_p(nzb_s_inner(j,i),j,i) = q_p(nzb_s_inner(j,i)+1,j,i)
215	ENDDO
216	ENDDO
217	ENDIF
218	!
219	!-- Top boundary
220	q_p(nzt+1,:,:) = q_p(nzt,:,:) + bc_q_t_val * dzu(nzt+1)
221
222	IF ( cloud_physics .AND. icloud_scheme == 0 ) THEN
223	!
224	!-- Surface conditions for constant_humidity_flux
225	IF ( ibc_qr_b == 0 ) THEN
226	DO i = nxlg, nxrg
227	DO j = nysg, nyng
228	qr_p(nzb_s_inner(j,i),j,i) = qr(nzb_s_inner(j,i),j,i)
229	ENDDO
230	ENDDO
231	ELSE
232	DO i = nxlg, nxrg
233	DO j = nysg, nyng
234	qr_p(nzb_s_inner(j,i),j,i) = qr_p(nzb_s_inner(j,i)+1,j,i)
235	ENDDO
236	ENDDO
237	ENDIF
238	!
239	!-- Top boundary
240	qr_p(nzt+1,:,:) = qr_p(nzt,:,:) + bc_qr_t_val * dzu(nzt+1)
241	!
242	!-- Surface conditions for constant_humidity_flux
243	IF ( ibc_nr_b == 0 ) THEN
244	DO i = nxlg, nxrg
245	DO j = nysg, nyng
246	nr_p(nzb_s_inner(j,i),j,i) = nr(nzb_s_inner(j,i),j,i)
247	ENDDO
248	ENDDO
249	ELSE
250	DO i = nxlg, nxrg
251	DO j = nysg, nyng
252	nr_p(nzb_s_inner(j,i),j,i) = nr_p(nzb_s_inner(j,i)+1,j,i)
253	ENDDO
254	ENDDO
255	ENDIF
256	!
257	!-- Top boundary
258	nr_p(nzt+1,:,:) = nr_p(nzt,:,:) + bc_nr_t_val * dzu(nzt+1)
259	ENDIF
260
261	ENDIF
262
263	!
264	!-- In case of inflow at the south boundary the boundary for v is at nys
265	!-- and in case of inflow at the left boundary the boundary for u is at nxl.
266	!-- Since in prognostic_equations (cache optimized version) these levels are
267	!-- handled as a prognostic level, boundary values have to be restored here.
268	!-- For the SGS-TKE, Neumann boundary conditions are used at the inflow.
269	IF ( inflow_s ) THEN
270	v_p(:,nys,:) = v_p(:,nys-1,:)
271	IF ( .NOT. constant_diffusion ) e_p(:,nys-1,:) = e_p(:,nys,:)
272	ELSEIF ( inflow_n ) THEN
273	IF ( .NOT. constant_diffusion ) e_p(:,nyn+1,:) = e_p(:,nyn,:)
274	ELSEIF ( inflow_l ) THEN
275	u_p(:,:,nxl) = u_p(:,:,nxl-1)
276	IF ( .NOT. constant_diffusion ) e_p(:,:,nxl-1) = e_p(:,:,nxl)
277	ELSEIF ( inflow_r ) THEN
278	IF ( .NOT. constant_diffusion ) e_p(:,:,nxr+1) = e_p(:,:,nxr)
279	ENDIF
280
281	!
282	!-- Lateral boundary conditions for scalar quantities at the outflow
283	IF ( outflow_s ) THEN
284	pt_p(:,nys-1,:) = pt_p(:,nys,:)
285	IF ( .NOT. constant_diffusion ) e_p(:,nys-1,:) = e_p(:,nys,:)
286	IF ( humidity .OR. passive_scalar ) THEN
287	q_p(:,nys-1,:) = q_p(:,nys,:)
288	IF ( cloud_physics .AND. icloud_scheme == 0 ) THEN
289	qr_p(:,nys-1,:) = qr_p(:,nys,:)
290	nr_p(:,nys-1,:) = nr_p(:,nys,:)
291	ENDIF
292	ENDIF
293	ELSEIF ( outflow_n ) THEN
294	pt_p(:,nyn+1,:) = pt_p(:,nyn,:)
295	IF ( .NOT. constant_diffusion ) e_p(:,nyn+1,:) = e_p(:,nyn,:)
296	IF ( humidity .OR. passive_scalar ) THEN
297	q_p(:,nyn+1,:) = q_p(:,nyn,:)
298	IF ( cloud_physics .AND. icloud_scheme == 0 ) THEN
299	qr_p(:,nyn+1,:) = qr_p(:,nyn,:)
300	nr_p(:,nyn+1,:) = nr_p(:,nyn,:)
301	ENDIF
302	ENDIF
303	ELSEIF ( outflow_l ) THEN
304	pt_p(:,:,nxl-1) = pt_p(:,:,nxl)
305	IF ( .NOT. constant_diffusion ) e_p(:,:,nxl-1) = e_p(:,:,nxl)
306	IF ( humidity .OR. passive_scalar ) THEN
307	q_p(:,:,nxl-1) = q_p(:,:,nxl)
308	IF ( cloud_physics .AND. icloud_scheme == 0 ) THEN
309	qr_p(:,:,nxl-1) = qr_p(:,:,nxl)
310	nr_p(:,:,nxl-1) = nr_p(:,:,nxl)
311	ENDIF
312	ENDIF
313	ELSEIF ( outflow_r ) THEN
314	pt_p(:,:,nxr+1) = pt_p(:,:,nxr)
315	IF ( .NOT. constant_diffusion ) e_p(:,:,nxr+1) = e_p(:,:,nxr)
316	IF ( humidity .OR. passive_scalar ) THEN
317	q_p(:,:,nxr+1) = q_p(:,:,nxr)
318	IF ( cloud_physics .AND. icloud_scheme == 0 ) THEN
319	qr_p(:,:,nxr+1) = qr_p(:,:,nxr)
320	nr_p(:,:,nxr+1) = nr_p(:,:,nxr)
321	ENDIF
322	ENDIF
323	ENDIF
324
325	ENDIF
326
327	!
328	!-- Neumann or Radiation boundary condition for the velocities at the
329	!-- respective outflow
330	IF ( outflow_s ) THEN
331
332	IF ( bc_ns_dirneu ) THEN
333	u(:,-1,:) = u(:,0,:)
334	v(:,0,:) = v(:,1,:)
335	w(:,-1,:) = w(:,0,:)
336	ELSEIF ( bc_ns_dirrad ) THEN
337
338	c_max = dy / dt_3d
339
340	c_u_m_l = 0.0
341	c_v_m_l = 0.0
342	c_w_m_l = 0.0
343
344	c_u_m = 0.0
345	c_v_m = 0.0
346	c_w_m = 0.0
347
348	!
349	!-- Calculate the phase speeds for u, v, and w, first local and then
350	!-- average along the outflow boundary.
351	DO k = nzb+1, nzt+1
352	DO i = nxl, nxr
353
354	denom = u_m_s(k,0,i) - u_m_s(k,1,i)
355
356	IF ( denom /= 0.0 ) THEN
357	c_u(k,i) = -c_max * ( u(k,0,i) - u_m_s(k,0,i) ) / ( denom * tsc(2) )
358	IF ( c_u(k,i) < 0.0 ) THEN
359	c_u(k,i) = 0.0
360	ELSEIF ( c_u(k,i) > c_max ) THEN
361	c_u(k,i) = c_max
362	ENDIF
363	ELSE
364	c_u(k,i) = c_max
365	ENDIF
366
367	denom = v_m_s(k,1,i) - v_m_s(k,2,i)
368
369	IF ( denom /= 0.0 ) THEN
370	c_v(k,i) = -c_max * ( v(k,1,i) - v_m_s(k,1,i) ) / ( denom * tsc(2) )
371	IF ( c_v(k,i) < 0.0 ) THEN
372	c_v(k,i) = 0.0
373	ELSEIF ( c_v(k,i) > c_max ) THEN
374	c_v(k,i) = c_max
375	ENDIF
376	ELSE
377	c_v(k,i) = c_max
378	ENDIF
379
380	denom = w_m_s(k,0,i) - w_m_s(k,1,i)
381
382	IF ( denom /= 0.0 ) THEN
383	c_w(k,i) = -c_max * ( w(k,0,i) - w_m_s(k,0,i) ) / ( denom * tsc(2) )
384	IF ( c_w(k,i) < 0.0 ) THEN
385	c_w(k,i) = 0.0
386	ELSEIF ( c_w(k,i) > c_max ) THEN
387	c_w(k,i) = c_max
388	ENDIF
389	ELSE
390	c_w(k,i) = c_max
391	ENDIF
392
393	c_u_m_l(k) = c_u_m_l(k) + c_u(k,i)
394	c_v_m_l(k) = c_v_m_l(k) + c_v(k,i)
395	c_w_m_l(k) = c_w_m_l(k) + c_w(k,i)
396
397	ENDDO
398	ENDDO
399
400	#if defined( __parallel )
401	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
402	CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
403	MPI_SUM, comm1dx, ierr )
404	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
405	CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
406	MPI_SUM, comm1dx, ierr )
407	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
408	CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
409	MPI_SUM, comm1dx, ierr )
410	#else
411	c_u_m = c_u_m_l
412	c_v_m = c_v_m_l
413	c_w_m = c_w_m_l
414	#endif
415
416	c_u_m = c_u_m / (nx+1)
417	c_v_m = c_v_m / (nx+1)
418	c_w_m = c_w_m / (nx+1)
419
420	!
421	!-- Save old timelevels for the next timestep
422	IF ( intermediate_timestep_count == 1 ) THEN
423	u_m_s(:,:,:) = u(:,0:1,:)
424	v_m_s(:,:,:) = v(:,1:2,:)
425	w_m_s(:,:,:) = w(:,0:1,:)
426	ENDIF
427
428	!
429	!-- Calculate the new velocities
430	DO k = nzb+1, nzt+1
431	DO i = nxlg, nxrg
432	u_p(k,-1,i) = u(k,-1,i) - dt_3d * tsc(2) * c_u_m(k) * &
433	( u(k,-1,i) - u(k,0,i) ) * ddy
434
435	v_p(k,0,i) = v(k,0,i) - dt_3d * tsc(2) * c_v_m(k) * &
436	( v(k,0,i) - v(k,1,i) ) * ddy
437
438	w_p(k,-1,i) = w(k,-1,i) - dt_3d * tsc(2) * c_w_m(k) * &
439	( w(k,-1,i) - w(k,0,i) ) * ddy
440	ENDDO
441	ENDDO
442
443	!
444	!-- Bottom boundary at the outflow
445	IF ( ibc_uv_b == 0 ) THEN
446	u_p(nzb,-1,:) = 0.0
447	v_p(nzb,0,:) = 0.0
448	ELSE
449	u_p(nzb,-1,:) = u_p(nzb+1,-1,:)
450	v_p(nzb,0,:) = v_p(nzb+1,0,:)
451	ENDIF
452	w_p(nzb,-1,:) = 0.0
453
454	!
455	!-- Top boundary at the outflow
456	IF ( ibc_uv_t == 0 ) THEN
457	u_p(nzt+1,-1,:) = u_init(nzt+1)
458	v_p(nzt+1,0,:) = v_init(nzt+1)
459	ELSE
460	u_p(nzt+1,-1,:) = u(nzt,-1,:)
461	v_p(nzt+1,0,:) = v(nzt,0,:)
462	ENDIF
463	w_p(nzt:nzt+1,-1,:) = 0.0
464
465	ENDIF
466
467	ENDIF
468
469	IF ( outflow_n ) THEN
470
471	IF ( bc_ns_neudir ) THEN
472	u(:,ny+1,:) = u(:,ny,:)
473	v(:,ny+1,:) = v(:,ny,:)
474	w(:,ny+1,:) = w(:,ny,:)
475	ELSEIF ( bc_ns_dirrad ) THEN
476
477	c_max = dy / dt_3d
478
479	c_u_m_l = 0.0
480	c_v_m_l = 0.0
481	c_w_m_l = 0.0
482
483	c_u_m = 0.0
484	c_v_m = 0.0
485	c_w_m = 0.0
486
487	!
488	!-- Calculate the phase speeds for u, v, and w, first local and then
489	!-- average along the outflow boundary.
490	DO k = nzb+1, nzt+1
491	DO i = nxl, nxr
492
493	denom = u_m_n(k,ny,i) - u_m_n(k,ny-1,i)
494
495	IF ( denom /= 0.0 ) THEN
496	c_u(k,i) = -c_max * ( u(k,ny,i) - u_m_n(k,ny,i) ) / ( denom * tsc(2) )
497	IF ( c_u(k,i) < 0.0 ) THEN
498	c_u(k,i) = 0.0
499	ELSEIF ( c_u(k,i) > c_max ) THEN
500	c_u(k,i) = c_max
501	ENDIF
502	ELSE
503	c_u(k,i) = c_max
504	ENDIF
505
506	denom = v_m_n(k,ny,i) - v_m_n(k,ny-1,i)
507
508	IF ( denom /= 0.0 ) THEN
509	c_v(k,i) = -c_max * ( v(k,ny,i) - v_m_n(k,ny,i) ) / ( denom * tsc(2) )
510	IF ( c_v(k,i) < 0.0 ) THEN
511	c_v(k,i) = 0.0
512	ELSEIF ( c_v(k,i) > c_max ) THEN
513	c_v(k,i) = c_max
514	ENDIF
515	ELSE
516	c_v(k,i) = c_max
517	ENDIF
518
519	denom = w_m_n(k,ny,i) - w_m_n(k,ny-1,i)
520
521	IF ( denom /= 0.0 ) THEN
522	c_w(k,i) = -c_max * ( w(k,ny,i) - w_m_n(k,ny,i) ) / ( denom * tsc(2) )
523	IF ( c_w(k,i) < 0.0 ) THEN
524	c_w(k,i) = 0.0
525	ELSEIF ( c_w(k,i) > c_max ) THEN
526	c_w(k,i) = c_max
527	ENDIF
528	ELSE
529	c_w(k,i) = c_max
530	ENDIF
531
532	c_u_m_l(k) = c_u_m_l(k) + c_u(k,i)
533	c_v_m_l(k) = c_v_m_l(k) + c_v(k,i)
534	c_w_m_l(k) = c_w_m_l(k) + c_w(k,i)
535
536	ENDDO
537	ENDDO
538
539	#if defined( __parallel )
540	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
541	CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
542	MPI_SUM, comm1dx, ierr )
543	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
544	CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
545	MPI_SUM, comm1dx, ierr )
546	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
547	CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
548	MPI_SUM, comm1dx, ierr )
549	#else
550	c_u_m = c_u_m_l
551	c_v_m = c_v_m_l
552	c_w_m = c_w_m_l
553	#endif
554
555	c_u_m = c_u_m / (nx+1)
556	c_v_m = c_v_m / (nx+1)
557	c_w_m = c_w_m / (nx+1)
558
559	!
560	!-- Save old timelevels for the next timestep
561	IF ( intermediate_timestep_count == 1 ) THEN
562	u_m_n(:,:,:) = u(:,ny-1:ny,:)
563	v_m_n(:,:,:) = v(:,ny-1:ny,:)
564	w_m_n(:,:,:) = w(:,ny-1:ny,:)
565	ENDIF
566
567	!
568	!-- Calculate the new velocities
569	DO k = nzb+1, nzt+1
570	DO i = nxlg, nxrg
571	u_p(k,ny+1,i) = u(k,ny+1,i) - dt_3d * tsc(2) * c_u_m(k) * &
572	( u(k,ny+1,i) - u(k,ny,i) ) * ddy
573
574	v_p(k,ny+1,i) = v(k,ny+1,i) - dt_3d * tsc(2) * c_v_m(k) * &
575	( v(k,ny+1,i) - v(k,ny,i) ) * ddy
576
577	w_p(k,ny+1,i) = w(k,ny+1,i) - dt_3d * tsc(2) * c_w_m(k) * &
578	( w(k,ny+1,i) - w(k,ny,i) ) * ddy
579	ENDDO
580	ENDDO
581
582	!
583	!-- Bottom boundary at the outflow
584	IF ( ibc_uv_b == 0 ) THEN
585	u_p(nzb,ny+1,:) = 0.0
586	v_p(nzb,ny+1,:) = 0.0
587	ELSE
588	u_p(nzb,ny+1,:) = u_p(nzb+1,ny+1,:)
589	v_p(nzb,ny+1,:) = v_p(nzb+1,ny+1,:)
590	ENDIF
591	w_p(nzb,ny+1,:) = 0.0
592
593	!
594	!-- Top boundary at the outflow
595	IF ( ibc_uv_t == 0 ) THEN
596	u_p(nzt+1,ny+1,:) = u_init(nzt+1)
597	v_p(nzt+1,ny+1,:) = v_init(nzt+1)
598	ELSE
599	u_p(nzt+1,ny+1,:) = u_p(nzt,nyn+1,:)
600	v_p(nzt+1,ny+1,:) = v_p(nzt,nyn+1,:)
601	ENDIF
602	w_p(nzt:nzt+1,ny+1,:) = 0.0
603
604	ENDIF
605
606	ENDIF
607
608	IF ( outflow_l ) THEN
609
610	IF ( bc_lr_neudir ) THEN
611	u(:,:,-1) = u(:,:,0)
612	v(:,:,0) = v(:,:,1)
613	w(:,:,-1) = w(:,:,0)
614	ELSEIF ( bc_ns_dirrad ) THEN
615
616	c_max = dx / dt_3d
617
618	c_u_m_l = 0.0
619	c_v_m_l = 0.0
620	c_w_m_l = 0.0
621
622	c_u_m = 0.0
623	c_v_m = 0.0
624	c_w_m = 0.0
625
626	!
627	!-- Calculate the phase speeds for u, v, and w, first local and then
628	!-- average along the outflow boundary.
629	DO k = nzb+1, nzt+1
630	DO j = nys, nyn
631
632	denom = u_m_l(k,j,1) - u_m_l(k,j,2)
633
634	IF ( denom /= 0.0 ) THEN
635	c_u(k,j) = -c_max * ( u(k,j,1) - u_m_l(k,j,1) ) / ( denom * tsc(2) )
636	IF ( c_u(k,j) < 0.0 ) THEN
637	c_u(k,j) = 0.0
638	ELSEIF ( c_u(k,j) > c_max ) THEN
639	c_u(k,j) = c_max
640	ENDIF
641	ELSE
642	c_u(k,j) = c_max
643	ENDIF
644
645	denom = v_m_l(k,j,0) - v_m_l(k,j,1)
646
647	IF ( denom /= 0.0 ) THEN
648	c_v(k,j) = -c_max * ( v(k,j,0) - v_m_l(k,j,0) ) / ( denom * tsc(2) )
649	IF ( c_v(k,j) < 0.0 ) THEN
650	c_v(k,j) = 0.0
651	ELSEIF ( c_v(k,j) > c_max ) THEN
652	c_v(k,j) = c_max
653	ENDIF
654	ELSE
655	c_v(k,j) = c_max
656	ENDIF
657
658	denom = w_m_l(k,j,0) - w_m_l(k,j,1)
659
660	IF ( denom /= 0.0 ) THEN
661	c_w(k,j) = -c_max * ( w(k,j,0) - w_m_l(k,j,0) ) / ( denom * tsc(2) )
662	IF ( c_w(k,j) < 0.0 ) THEN
663	c_w(k,j) = 0.0
664	ELSEIF ( c_w(k,j) > c_max ) THEN
665	c_w(k,j) = c_max
666	ENDIF
667	ELSE
668	c_w(k,j) = c_max
669	ENDIF
670
671	c_u_m_l(k) = c_u_m_l(k) + c_u(k,j)
672	c_v_m_l(k) = c_v_m_l(k) + c_v(k,j)
673	c_w_m_l(k) = c_w_m_l(k) + c_w(k,j)
674
675	ENDDO
676	ENDDO
677
678	#if defined( __parallel )
679	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
680	CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
681	MPI_SUM, comm1dy, ierr )
682	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
683	CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
684	MPI_SUM, comm1dy, ierr )
685	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
686	CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
687	MPI_SUM, comm1dy, ierr )
688	#else
689	c_u_m = c_u_m_l
690	c_v_m = c_v_m_l
691	c_w_m = c_w_m_l
692	#endif
693
694	c_u_m = c_u_m / (ny+1)
695	c_v_m = c_v_m / (ny+1)
696	c_w_m = c_w_m / (ny+1)
697
698	!
699	!-- Save old timelevels for the next timestep
700	IF ( intermediate_timestep_count == 1 ) THEN
701	u_m_l(:,:,:) = u(:,:,1:2)
702	v_m_l(:,:,:) = v(:,:,0:1)
703	w_m_l(:,:,:) = w(:,:,0:1)
704	ENDIF
705
706	!
707	!-- Calculate the new velocities
708	DO k = nzb+1, nzt+1
709	DO i = nxlg, nxrg
710	u_p(k,j,0) = u(k,j,0) - dt_3d * tsc(2) * c_u_m(k) * &
711	( u(k,j,0) - u(k,j,1) ) * ddx
712
713	v_p(k,j,-1) = v(k,j,-1) - dt_3d * tsc(2) * c_v_m(k) * &
714	( v(k,j,-1) - v(k,j,0) ) * ddx
715
716	w_p(k,j,-1) = w(k,j,-1) - dt_3d * tsc(2) * c_w_m(k) * &
717	( w(k,j,-1) - w(k,j,0) ) * ddx
718	ENDDO
719	ENDDO
720
721	!
722	!-- Bottom boundary at the outflow
723	IF ( ibc_uv_b == 0 ) THEN
724	u_p(nzb,:,0) = 0.0
725	v_p(nzb,:,-1) = 0.0
726	ELSE
727	u_p(nzb,:,0) = u_p(nzb+1,:,0)
728	v_p(nzb,:,-1) = v_p(nzb+1,:,-1)
729	ENDIF
730	w_p(nzb,:,-1) = 0.0
731
732	!
733	!-- Top boundary at the outflow
734	IF ( ibc_uv_t == 0 ) THEN
735	u_p(nzt+1,:,-1) = u_init(nzt+1)
736	v_p(nzt+1,:,-1) = v_init(nzt+1)
737	ELSE
738	u_p(nzt+1,:,-1) = u_p(nzt,:,-1)
739	v_p(nzt+1,:,-1) = v_p(nzt,:,-1)
740	ENDIF
741	w_p(nzt:nzt+1,:,-1) = 0.0
742
743	ENDIF
744
745	ENDIF
746
747	IF ( outflow_r ) THEN
748
749	IF ( bc_lr_dirneu ) THEN
750	u(:,:,nx+1) = u(:,:,nx)
751	v(:,:,nx+1) = v(:,:,nx)
752	w(:,:,nx+1) = w(:,:,nx)
753	ELSEIF ( bc_ns_dirrad ) THEN
754
755	c_max = dx / dt_3d
756
757	c_u_m_l = 0.0
758	c_v_m_l = 0.0
759	c_w_m_l = 0.0
760
761	c_u_m = 0.0
762	c_v_m = 0.0
763	c_w_m = 0.0
764
765	!
766	!-- Calculate the phase speeds for u, v, and w, first local and then
767	!-- average along the outflow boundary.
768	DO k = nzb+1, nzt+1
769	DO j = nys, nyn
770
771	denom = u_m_r(k,j,nx) - u_m_r(k,j,nx-1)
772
773	IF ( denom /= 0.0 ) THEN
774	c_u(k,j) = -c_max * ( u(k,j,nx) - u_m_r(k,j,nx) ) / ( denom * tsc(2) )
775	IF ( c_u(k,j) < 0.0 ) THEN
776	c_u(k,j) = 0.0
777	ELSEIF ( c_u(k,j) > c_max ) THEN
778	c_u(k,j) = c_max
779	ENDIF
780	ELSE
781	c_u(k,j) = c_max
782	ENDIF
783
784	denom = v_m_r(k,j,nx) - v_m_r(k,j,nx-1)
785
786	IF ( denom /= 0.0 ) THEN
787	c_v(k,j) = -c_max * ( v(k,j,nx) - v_m_r(k,j,nx) ) / ( denom * tsc(2) )
788	IF ( c_v(k,j) < 0.0 ) THEN
789	c_v(k,j) = 0.0
790	ELSEIF ( c_v(k,j) > c_max ) THEN
791	c_v(k,j) = c_max
792	ENDIF
793	ELSE
794	c_v(k,j) = c_max
795	ENDIF
796
797	denom = w_m_r(k,j,nx) - w_m_r(k,j,nx-1)
798
799	IF ( denom /= 0.0 ) THEN
800	c_w(k,j) = -c_max * ( w(k,j,nx) - w_m_r(k,j,nx) ) / ( denom * tsc(2) )
801	IF ( c_w(k,j) < 0.0 ) THEN
802	c_w(k,j) = 0.0
803	ELSEIF ( c_w(k,j) > c_max ) THEN
804	c_w(k,j) = c_max
805	ENDIF
806	ELSE
807	c_w(k,j) = c_max
808	ENDIF
809
810	c_u_m_l(k) = c_u_m_l(k) + c_u(k,j)
811	c_v_m_l(k) = c_v_m_l(k) + c_v(k,j)
812	c_w_m_l(k) = c_w_m_l(k) + c_w(k,j)
813
814	ENDDO
815	ENDDO
816
817	#if defined( __parallel )
818	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
819	CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
820	MPI_SUM, comm1dy, ierr )
821	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
822	CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
823	MPI_SUM, comm1dy, ierr )
824	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
825	CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
826	MPI_SUM, comm1dy, ierr )
827	#else
828	c_u_m = c_u_m_l
829	c_v_m = c_v_m_l
830	c_w_m = c_w_m_l
831	#endif
832
833	c_u_m = c_u_m / (ny+1)
834	c_v_m = c_v_m / (ny+1)
835	c_w_m = c_w_m / (ny+1)
836
837	!
838	!-- Save old timelevels for the next timestep
839	IF ( intermediate_timestep_count == 1 ) THEN
840	u_m_r(:,:,:) = u(:,:,nx-1:nx)
841	v_m_r(:,:,:) = v(:,:,nx-1:nx)
842	w_m_r(:,:,:) = w(:,:,nx-1:nx)
843	ENDIF
844
845	!
846	!-- Calculate the new velocities
847	DO k = nzb+1, nzt+1
848	DO i = nxlg, nxrg
849	u_p(k,j,nx+1) = u(k,j,nx+1) - dt_3d * tsc(2) * c_u_m(k) * &
850	( u(k,j,nx+1) - u(k,j,nx) ) * ddx
851
852	v_p(k,j,nx+1) = v(k,j,nx+1) - dt_3d * tsc(2) * c_v_m(k) * &
853	( v(k,j,nx+1) - v(k,j,nx) ) * ddx
854
855	w_p(k,j,nx+1) = w(k,j,nx+1) - dt_3d * tsc(2) * c_w_m(k) * &
856	( w(k,j,nx+1) - w(k,j,nx) ) * ddx
857	ENDDO
858	ENDDO
859
860	!
861	!-- Bottom boundary at the outflow
862	IF ( ibc_uv_b == 0 ) THEN
863	u_p(nzb,:,nx+1) = 0.0
864	v_p(nzb,:,nx+1) = 0.0
865	ELSE
866	u_p(nzb,:,nx+1) = u_p(nzb+1,:,nx+1)
867	v_p(nzb,:,nx+1) = v_p(nzb+1,:,nx+1)
868	ENDIF
869	w_p(nzb,:,nx+1) = 0.0
870
871	!
872	!-- Top boundary at the outflow
873	IF ( ibc_uv_t == 0 ) THEN
874	u_p(nzt+1,:,nx+1) = u_init(nzt+1)
875	v_p(nzt+1,:,nx+1) = v_init(nzt+1)
876	ELSE
877	u_p(nzt+1,:,nx+1) = u_p(nzt,:,nx+1)
878	v_p(nzt+1,:,nx+1) = v_p(nzt,:,nx+1)
879	ENDIF
880	w(nzt:nzt+1,:,nx+1) = 0.0
881
882	ENDIF
883
884	ENDIF
885
886	END SUBROUTINE boundary_conds

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |