Home

Context Navigation

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

Last change on this file since 2284 was 2233, checked in by suehring, 8 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 37.5 KB

Line
1	!> @file boundary_conds.f90
2	!------------------------------------------------------------------------------!
3	! This file is part of PALM.
4	!
5	! PALM is free software: you can redistribute it and/or modify it under the
6	! terms of the GNU General Public License as published by the Free Software
7	! Foundation, either version 3 of the License, or (at your option) any later
8	! 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-2017 Leibniz Universitaet Hannover
18	!------------------------------------------------------------------------------!
19	!
20	! Current revisions:
21	! -----------------
22	!
23	!
24	! Former revisions:
25	! -----------------
26	! $Id: boundary_conds.f90 2233 2017-05-30 18:08:54Z witha $
27	!
28	! 2232 2017-05-30 17:47:52Z suehring
29	! Set boundary conditions on topography top using flag method.
30	!
31	! 2118 2017-01-17 16:38:49Z raasch
32	! OpenACC directives removed
33	!
34	! 2000 2016-08-20 18:09:15Z knoop
35	! Forced header and separation lines into 80 columns
36	!
37	! 1992 2016-08-12 15:14:59Z suehring
38	! Adjustments for top boundary condition for passive scalar
39	!
40	! 1960 2016-07-12 16:34:24Z suehring
41	! Treat humidity and passive scalar separately
42	!
43	! 1823 2016-04-07 08:57:52Z hoffmann
44	! Initial version of purely vertical nesting introduced.
45	!
46	! 1822 2016-04-07 07:49:42Z hoffmann
47	! icloud_scheme removed. microphyisics_seifert added.
48	!
49	! 1764 2016-02-28 12:45:19Z raasch
50	! index bug for u_p at left outflow removed
51	!
52	! 1762 2016-02-25 12:31:13Z hellstea
53	! Introduction of nested domain feature
54	!
55	! 1742 2016-01-13 09:50:06Z raasch
56	! bugfix for outflow Neumann boundary conditions at bottom and top
57	!
58	! 1717 2015-11-11 15:09:47Z raasch
59	! Bugfix: index error in outflow conditions for left boundary
60	!
61	! 1682 2015-10-07 23:56:08Z knoop
62	! Code annotations made doxygen readable
63	!
64	! 1410 2014-05-23 12:16:18Z suehring
65	! Bugfix: set dirichlet boundary condition for passive_scalar at model domain
66	! top
67	!
68	! 1399 2014-05-07 11:16:25Z heinze
69	! Bugfix: set inflow boundary conditions also if no humidity or passive_scalar
70	! is used.
71	!
72	! 1398 2014-05-07 11:15:00Z heinze
73	! Dirichlet-condition at the top for u and v changed to u_init and v_init also
74	! for large_scale_forcing
75	!
76	! 1380 2014-04-28 12:40:45Z heinze
77	! Adjust Dirichlet-condition at the top for pt in case of nudging
78	!
79	! 1361 2014-04-16 15:17:48Z hoffmann
80	! Bottom and top boundary conditions of rain water content (qr) and
81	! rain drop concentration (nr) changed to Dirichlet
82	!
83	! 1353 2014-04-08 15:21:23Z heinze
84	! REAL constants provided with KIND-attribute
85	!
86	! 1320 2014-03-20 08:40:49Z raasch
87	! ONLY-attribute added to USE-statements,
88	! kind-parameters added to all INTEGER and REAL declaration statements,
89	! kinds are defined in new module kinds,
90	! revision history before 2012 removed,
91	! comment fields (!:) to be used for variable explanations added to
92	! all variable declaration statements
93	!
94	! 1257 2013-11-08 15:18:40Z raasch
95	! loop independent clauses added
96	!
97	! 1241 2013-10-30 11:36:58Z heinze
98	! Adjust ug and vg at each timestep in case of large_scale_forcing
99	!
100	! 1159 2013-05-21 11:58:22Z fricke
101	! Bugfix: Neumann boundary conditions for the velocity components at the
102	! outflow are in fact radiation boundary conditions using the maximum phase
103	! velocity that ensures numerical stability (CFL-condition).
104	! Hence, logical operator use_cmax is now used instead of bc_lr_dirneu/_neudir.
105	! Bugfix: In case of use_cmax at the outflow, u, v, w are replaced by
106	! u_p, v_p, w_p
107	!
108	! 1115 2013-03-26 18:16:16Z hoffmann
109	! boundary conditions of two-moment cloud scheme are restricted to Neumann-
110	! boundary-conditions
111	!
112	! 1113 2013-03-10 02:48:14Z raasch
113	! GPU-porting
114	! dummy argument "range" removed
115	! Bugfix: wrong index in loops of radiation boundary condition
116	!
117	! 1053 2012-11-13 17:11:03Z hoffmann
118	! boundary conditions for the two new prognostic equations (nr, qr) of the
119	! two-moment cloud scheme
120	!
121	! 1036 2012-10-22 13:43:42Z raasch
122	! code put under GPL (PALM 3.9)
123	!
124	! 996 2012-09-07 10:41:47Z raasch
125	! little reformatting
126	!
127	! 978 2012-08-09 08:28:32Z fricke
128	! Neumann boudnary conditions are added at the inflow boundary for the SGS-TKE.
129	! Outflow boundary conditions for the velocity components can be set to Neumann
130	! conditions or to radiation conditions with a horizontal averaged phase
131	! velocity.
132	!
133	! 875 2012-04-02 15:35:15Z gryschka
134	! Bugfix in case of dirichlet inflow bc at the right or north boundary
135	!
136	! Revision 1.1 1997/09/12 06:21:34 raasch
137	! Initial revision
138	!
139	!
140	! Description:
141	! ------------
142	!> Boundary conditions for the prognostic quantities.
143	!> One additional bottom boundary condition is applied for the TKE (=(u)*2)
144	!> in prandtl_fluxes. The cyclic lateral boundary conditions are implicitly
145	!> handled in routine exchange_horiz. Pressure boundary conditions are
146	!> explicitly set in routines pres, poisfft, poismg and sor.
147	!------------------------------------------------------------------------------!
148	SUBROUTINE boundary_conds
149
150
151	USE arrays_3d, &
152	ONLY: c_u, c_u_m, c_u_m_l, c_v, c_v_m, c_v_m_l, c_w, c_w_m, c_w_m_l, &
153	dzu, e_p, nr_p, pt, pt_p, q, q_p, qr_p, s, s_p, sa, sa_p, &
154	u, ug, u_init, u_m_l, u_m_n, u_m_r, u_m_s, u_p, &
155	v, vg, v_init, v_m_l, v_m_n, v_m_r, v_m_s, v_p, &
156	w, w_p, w_m_l, w_m_n, w_m_r, w_m_s, pt_init
157
158	USE control_parameters, &
159	ONLY: bc_pt_t_val, bc_q_t_val, bc_s_t_val, constant_diffusion, &
160	cloud_physics, dt_3d, humidity, &
161	ibc_pt_b, ibc_pt_t, ibc_q_b, ibc_q_t, ibc_s_b, ibc_s_t, &
162	ibc_sa_t, ibc_uv_b, ibc_uv_t, inflow_l, inflow_n, inflow_r, &
163	inflow_s, intermediate_timestep_count, large_scale_forcing, &
164	microphysics_seifert, nest_domain, nest_bound_l, nest_bound_s, &
165	nudging, ocean, outflow_l, outflow_n, outflow_r, outflow_s, &
166	passive_scalar, tsc, use_cmax
167
168	USE grid_variables, &
169	ONLY: ddx, ddy, dx, dy
170
171	USE indices, &
172	ONLY: nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, &
173	nzb, nzt, wall_flags_0
174
175	USE kinds
176
177	USE pegrid
178
179	USE pmc_interface, &
180	ONLY : nesting_mode
181
182	USE surface_mod, &
183	ONLY : bc_h
184
185	IMPLICIT NONE
186
187	INTEGER(iwp) :: i !< grid index x direction
188	INTEGER(iwp) :: j !< grid index y direction
189	INTEGER(iwp) :: k !< grid index z direction
190	INTEGER(iwp) :: kb !< variable to set respective boundary value, depends on facing.
191	INTEGER(iwp) :: l !< running index boundary type, for up- and downward-facing walls
192	INTEGER(iwp) :: m !< running index surface elements
193
194	REAL(wp) :: c_max !<
195	REAL(wp) :: denom !<
196
197
198	!
199	!-- Bottom boundary
200	IF ( ibc_uv_b == 1 ) THEN
201	u_p(nzb,:,:) = u_p(nzb+1,:,:)
202	v_p(nzb,:,:) = v_p(nzb+1,:,:)
203	ENDIF
204	!
205	!-- Set zero vertical velocity at topography top (l=0), or bottom (l=1) in case
206	!-- of downward-facing surfaces.
207	DO l = 0, 1
208	!
209	!-- Set kb, for upward-facing surfaces value at topography top (k-1) is set,
210	!-- for downward-facing surfaces at topography bottom (k+1).
211	kb = MERGE( -1, 1, l == 0 )
212	!$OMP PARALLEL DO PRIVATE( i, j, k )
213	DO m = 1, bc_h(l)%ns
214	i = bc_h(l)%i(m)
215	j = bc_h(l)%j(m)
216	k = bc_h(l)%k(m)
217	w_p(k+kb,j,i) = 0.0_wp
218	ENDDO
219	ENDDO
220
221	!
222	!-- Top boundary. A nested domain ( ibc_uv_t = 3 ) does not require settings.
223	IF ( ibc_uv_t == 0 ) THEN
224	u_p(nzt+1,:,:) = u_init(nzt+1)
225	v_p(nzt+1,:,:) = v_init(nzt+1)
226	ELSEIF ( ibc_uv_t == 1 ) THEN
227	u_p(nzt+1,:,:) = u_p(nzt,:,:)
228	v_p(nzt+1,:,:) = v_p(nzt,:,:)
229	ENDIF
230
231	IF ( .NOT. nest_domain ) THEN
232	w_p(nzt:nzt+1,:,:) = 0.0_wp ! nzt is not a prognostic level (but cf. pres)
233	ENDIF
234
235	!
236	!-- Temperature at bottom and top boundary.
237	!-- In case of coupled runs (ibc_pt_b = 2) the temperature is given by
238	!-- the sea surface temperature of the coupled ocean model.
239	!-- Dirichlet
240	IF ( ibc_pt_b == 0 ) THEN
241	DO l = 0, 1
242	!
243	!-- Set kb, for upward-facing surfaces value at topography top (k-1) is set,
244	!-- for downward-facing surfaces at topography bottom (k+1).
245	kb = MERGE( -1, 1, l == 0 )
246	!$OMP PARALLEL DO PRIVATE( i, j, k )
247	DO m = 1, bc_h(l)%ns
248	i = bc_h(l)%i(m)
249	j = bc_h(l)%j(m)
250	k = bc_h(l)%k(m)
251	pt_p(k+kb,j,i) = pt(k+kb,j,i)
252	ENDDO
253	ENDDO
254	!
255	!-- Neumann, zero-gradient
256	ELSEIF ( ibc_pt_b == 1 ) THEN
257	DO l = 0, 1
258	!
259	!-- Set kb, for upward-facing surfaces value at topography top (k-1) is set,
260	!-- for downward-facing surfaces at topography bottom (k+1).
261	kb = MERGE( -1, 1, l == 0 )
262	!$OMP PARALLEL DO PRIVATE( i, j, k )
263	DO m = 1, bc_h(l)%ns
264	i = bc_h(l)%i(m)
265	j = bc_h(l)%j(m)
266	k = bc_h(l)%k(m)
267	pt_p(k+kb,j,i) = pt_p(k,j,i)
268	ENDDO
269	ENDDO
270	ENDIF
271
272	!
273	!-- Temperature at top boundary
274	IF ( ibc_pt_t == 0 ) THEN
275	pt_p(nzt+1,:,:) = pt(nzt+1,:,:)
276	!
277	!-- In case of nudging adjust top boundary to pt which is
278	!-- read in from NUDGING-DATA
279	IF ( nudging ) THEN
280	pt_p(nzt+1,:,:) = pt_init(nzt+1)
281	ENDIF
282	ELSEIF ( ibc_pt_t == 1 ) THEN
283	pt_p(nzt+1,:,:) = pt_p(nzt,:,:)
284	ELSEIF ( ibc_pt_t == 2 ) THEN
285	pt_p(nzt+1,:,:) = pt_p(nzt,:,:) + bc_pt_t_val * dzu(nzt+1)
286	ENDIF
287
288	!
289	!-- Boundary conditions for TKE
290	!-- Generally Neumann conditions with de/dz=0 are assumed
291	IF ( .NOT. constant_diffusion ) THEN
292
293	DO l = 0, 1
294	!
295	!-- Set kb, for upward-facing surfaces value at topography top (k-1) is set,
296	!-- for downward-facing surfaces at topography bottom (k+1).
297	kb = MERGE( -1, 1, l == 0 )
298	!$OMP PARALLEL DO PRIVATE( i, j, k )
299	DO m = 1, bc_h(l)%ns
300	i = bc_h(l)%i(m)
301	j = bc_h(l)%j(m)
302	k = bc_h(l)%k(m)
303	e_p(k+kb,j,i) = e_p(k,j,i)
304	ENDDO
305	ENDDO
306
307	IF ( .NOT. nest_domain ) THEN
308	e_p(nzt+1,:,:) = e_p(nzt,:,:)
309	ENDIF
310	ENDIF
311
312	!
313	!-- Boundary conditions for salinity
314	IF ( ocean ) THEN
315	!
316	!-- Bottom boundary: Neumann condition because salinity flux is always
317	!-- given.
318	DO l = 0, 1
319	!
320	!-- Set kb, for upward-facing surfaces value at topography top (k-1) is set,
321	!-- for downward-facing surfaces at topography bottom (k+1).
322	kb = MERGE( -1, 1, l == 0 )
323	!$OMP PARALLEL DO PRIVATE( i, j, k )
324	DO m = 1, bc_h(l)%ns
325	i = bc_h(l)%i(m)
326	j = bc_h(l)%j(m)
327	k = bc_h(l)%k(m)
328	sa_p(k+kb,j,i) = sa_p(k,j,i)
329	ENDDO
330	ENDDO
331	!
332	!-- Top boundary: Dirichlet or Neumann
333	IF ( ibc_sa_t == 0 ) THEN
334	sa_p(nzt+1,:,:) = sa(nzt+1,:,:)
335	ELSEIF ( ibc_sa_t == 1 ) THEN
336	sa_p(nzt+1,:,:) = sa_p(nzt,:,:)
337	ENDIF
338
339	ENDIF
340
341	!
342	!-- Boundary conditions for total water content,
343	!-- bottom and top boundary (see also temperature)
344	IF ( humidity ) THEN
345	!
346	!-- Surface conditions for constant_humidity_flux
347	!-- Run loop over all non-natural and natural walls. Note, in wall-datatype
348	!-- the k coordinate belongs to the atmospheric grid point, therefore, set
349	!-- q_p at k-1
350	IF ( ibc_q_b == 0 ) THEN
351
352	DO l = 0, 1
353	!
354	!-- Set kb, for upward-facing surfaces value at topography top (k-1) is set,
355	!-- for downward-facing surfaces at topography bottom (k+1).
356	kb = MERGE( -1, 1, l == 0 )
357	!$OMP PARALLEL DO PRIVATE( i, j, k )
358	DO m = 1, bc_h(l)%ns
359	i = bc_h(l)%i(m)
360	j = bc_h(l)%j(m)
361	k = bc_h(l)%k(m)
362	q_p(k+kb,j,i) = q(k+kb,j,i)
363	ENDDO
364	ENDDO
365
366	ELSE
367	!$OMP PARALLEL DO PRIVATE( i, j, k )
368	DO m = 1, bc_h(0)%ns
369	i = bc_h(0)%i(m)
370	j = bc_h(0)%j(m)
371	k = bc_h(0)%k(m)
372	q_p(k-1,j,i) = q_p(k,j,i)
373	ENDDO
374
375	DO l = 0, 1
376	!
377	!-- Set kb, for upward-facing surfaces value at topography top (k-1) is set,
378	!-- for downward-facing surfaces at topography bottom (k+1).
379	kb = MERGE( -1, 1, l == 0 )
380	!$OMP PARALLEL DO PRIVATE( i, j, k )
381	DO m = 1, bc_h(l)%ns
382	i = bc_h(l)%i(m)
383	j = bc_h(l)%j(m)
384	k = bc_h(l)%k(m)
385	q_p(k+kb,j,i) = q_p(k,j,i)
386	ENDDO
387	ENDDO
388	ENDIF
389	!
390	!-- Top boundary
391	IF ( ibc_q_t == 0 ) THEN
392	q_p(nzt+1,:,:) = q(nzt+1,:,:)
393	ELSEIF ( ibc_q_t == 1 ) THEN
394	q_p(nzt+1,:,:) = q_p(nzt,:,:) + bc_q_t_val * dzu(nzt+1)
395	ENDIF
396
397	IF ( cloud_physics .AND. microphysics_seifert ) THEN
398	!
399	!-- Surface conditions rain water (Dirichlet)
400	!-- Run loop over all non-natural and natural walls. Note, in wall-datatype
401	!-- the k coordinate belongs to the atmospheric grid point, therefore, set
402	!-- qr_p and nr_p at k-1
403	!$OMP PARALLEL DO PRIVATE( i, j, k )
404	DO m = 1, bc_h(0)%ns
405	i = bc_h(0)%i(m)
406	j = bc_h(0)%j(m)
407	k = bc_h(0)%k(m)
408	qr_p(k-1,j,i) = 0.0_wp
409	nr_p(k-1,j,i) = 0.0_wp
410	ENDDO
411	!
412	!-- Top boundary condition for rain water (Dirichlet)
413	qr_p(nzt+1,:,:) = 0.0_wp
414	nr_p(nzt+1,:,:) = 0.0_wp
415
416	ENDIF
417	ENDIF
418	!
419	!-- Boundary conditions for scalar,
420	!-- bottom and top boundary (see also temperature)
421	IF ( passive_scalar ) THEN
422	!
423	!-- Surface conditions for constant_humidity_flux
424	!-- Run loop over all non-natural and natural walls. Note, in wall-datatype
425	!-- the k coordinate belongs to the atmospheric grid point, therefore, set
426	!-- s_p at k-1
427	IF ( ibc_s_b == 0 ) THEN
428
429	DO l = 0, 1
430	!
431	!-- Set kb, for upward-facing surfaces value at topography top (k-1) is set,
432	!-- for downward-facing surfaces at topography bottom (k+1).
433	kb = MERGE( -1, 1, l == 0 )
434	!$OMP PARALLEL DO PRIVATE( i, j, k )
435	DO m = 1, bc_h(l)%ns
436	i = bc_h(l)%i(m)
437	j = bc_h(l)%j(m)
438	k = bc_h(l)%k(m)
439	s_p(k+kb,j,i) = s(k+kb,j,i)
440	ENDDO
441	ENDDO
442
443	ELSE
444	!$OMP PARALLEL DO PRIVATE( i, j, k )
445	DO m = 1, bc_h(0)%ns
446	i = bc_h(0)%i(m)
447	j = bc_h(0)%j(m)
448	k = bc_h(0)%k(m)
449	s_p(k-1,j,i) = s_p(k,j,i)
450	ENDDO
451
452	DO l = 0, 1
453	!
454	!-- Set kb, for upward-facing surfaces value at topography top (k-1) is set,
455	!-- for downward-facing surfaces at topography bottom (k+1).
456	kb = MERGE( -1, 1, l == 0 )
457	!$OMP PARALLEL DO PRIVATE( i, j, k )
458	DO m = 1, bc_h(l)%ns
459	i = bc_h(l)%i(m)
460	j = bc_h(l)%j(m)
461	k = bc_h(l)%k(m)
462	s_p(k+kb,j,i) = s_p(k,j,i)
463	ENDDO
464	ENDDO
465	ENDIF
466	!
467	!-- Top boundary condition
468	IF ( ibc_s_t == 0 ) THEN
469	s_p(nzt+1,:,:) = s(nzt+1,:,:)
470	ELSEIF ( ibc_s_t == 1 ) THEN
471	s_p(nzt+1,:,:) = s_p(nzt,:,:)
472	ELSEIF ( ibc_s_t == 2 ) THEN
473	s_p(nzt+1,:,:) = s_p(nzt,:,:) + bc_s_t_val * dzu(nzt+1)
474	ENDIF
475
476	ENDIF
477	!
478	!-- In case of inflow or nest boundary at the south boundary the boundary for v
479	!-- is at nys and in case of inflow or nest boundary at the left boundary the
480	!-- boundary for u is at nxl. Since in prognostic_equations (cache optimized
481	!-- version) these levels are handled as a prognostic level, boundary values
482	!-- have to be restored here.
483	!-- For the SGS-TKE, Neumann boundary conditions are used at the inflow.
484	IF ( inflow_s ) THEN
485	v_p(:,nys,:) = v_p(:,nys-1,:)
486	IF ( .NOT. constant_diffusion ) e_p(:,nys-1,:) = e_p(:,nys,:)
487	ELSEIF ( inflow_n ) THEN
488	IF ( .NOT. constant_diffusion ) e_p(:,nyn+1,:) = e_p(:,nyn,:)
489	ELSEIF ( inflow_l ) THEN
490	u_p(:,:,nxl) = u_p(:,:,nxl-1)
491	IF ( .NOT. constant_diffusion ) e_p(:,:,nxl-1) = e_p(:,:,nxl)
492	ELSEIF ( inflow_r ) THEN
493	IF ( .NOT. constant_diffusion ) e_p(:,:,nxr+1) = e_p(:,:,nxr)
494	ENDIF
495
496	!
497	!-- The same restoration for u at i=nxl and v at j=nys as above must be made
498	!-- in case of nest boundaries. This must not be done in case of vertical nesting
499	!-- mode as in that case the lateral boundaries are actually cyclic.
500	IF ( nesting_mode /= 'vertical' ) THEN
501	IF ( nest_bound_s ) THEN
502	v_p(:,nys,:) = v_p(:,nys-1,:)
503	ENDIF
504	IF ( nest_bound_l ) THEN
505	u_p(:,:,nxl) = u_p(:,:,nxl-1)
506	ENDIF
507	ENDIF
508
509	!
510	!-- Lateral boundary conditions for scalar quantities at the outflow
511	IF ( outflow_s ) THEN
512	pt_p(:,nys-1,:) = pt_p(:,nys,:)
513	IF ( .NOT. constant_diffusion ) e_p(:,nys-1,:) = e_p(:,nys,:)
514	IF ( humidity ) THEN
515	q_p(:,nys-1,:) = q_p(:,nys,:)
516	IF ( cloud_physics .AND. microphysics_seifert ) THEN
517	qr_p(:,nys-1,:) = qr_p(:,nys,:)
518	nr_p(:,nys-1,:) = nr_p(:,nys,:)
519	ENDIF
520	ENDIF
521	IF ( passive_scalar ) s_p(:,nys-1,:) = s_p(:,nys,:)
522	ELSEIF ( outflow_n ) THEN
523	pt_p(:,nyn+1,:) = pt_p(:,nyn,:)
524	IF ( .NOT. constant_diffusion ) e_p(:,nyn+1,:) = e_p(:,nyn,:)
525	IF ( humidity ) THEN
526	q_p(:,nyn+1,:) = q_p(:,nyn,:)
527	IF ( cloud_physics .AND. microphysics_seifert ) THEN
528	qr_p(:,nyn+1,:) = qr_p(:,nyn,:)
529	nr_p(:,nyn+1,:) = nr_p(:,nyn,:)
530	ENDIF
531	ENDIF
532	IF ( passive_scalar ) s_p(:,nyn+1,:) = s_p(:,nyn,:)
533	ELSEIF ( outflow_l ) THEN
534	pt_p(:,:,nxl-1) = pt_p(:,:,nxl)
535	IF ( .NOT. constant_diffusion ) e_p(:,:,nxl-1) = e_p(:,:,nxl)
536	IF ( humidity ) THEN
537	q_p(:,:,nxl-1) = q_p(:,:,nxl)
538	IF ( cloud_physics .AND. microphysics_seifert ) THEN
539	qr_p(:,:,nxl-1) = qr_p(:,:,nxl)
540	nr_p(:,:,nxl-1) = nr_p(:,:,nxl)
541	ENDIF
542	ENDIF
543	IF ( passive_scalar ) s_p(:,:,nxl-1) = s_p(:,:,nxl)
544	ELSEIF ( outflow_r ) THEN
545	pt_p(:,:,nxr+1) = pt_p(:,:,nxr)
546	IF ( .NOT. constant_diffusion ) e_p(:,:,nxr+1) = e_p(:,:,nxr)
547	IF ( humidity ) THEN
548	q_p(:,:,nxr+1) = q_p(:,:,nxr)
549	IF ( cloud_physics .AND. microphysics_seifert ) THEN
550	qr_p(:,:,nxr+1) = qr_p(:,:,nxr)
551	nr_p(:,:,nxr+1) = nr_p(:,:,nxr)
552	ENDIF
553	ENDIF
554	IF ( passive_scalar ) s_p(:,:,nxr+1) = s_p(:,:,nxr)
555	ENDIF
556
557	!
558	!-- Radiation boundary conditions for the velocities at the respective outflow.
559	!-- The phase velocity is either assumed to the maximum phase velocity that
560	!-- ensures numerical stability (CFL-condition) or calculated after
561	!-- Orlanski(1976) and averaged along the outflow boundary.
562	IF ( outflow_s ) THEN
563
564	IF ( use_cmax ) THEN
565	u_p(:,-1,:) = u(:,0,:)
566	v_p(:,0,:) = v(:,1,:)
567	w_p(:,-1,:) = w(:,0,:)
568	ELSEIF ( .NOT. use_cmax ) THEN
569
570	c_max = dy / dt_3d
571
572	c_u_m_l = 0.0_wp
573	c_v_m_l = 0.0_wp
574	c_w_m_l = 0.0_wp
575
576	c_u_m = 0.0_wp
577	c_v_m = 0.0_wp
578	c_w_m = 0.0_wp
579
580	!
581	!-- Calculate the phase speeds for u, v, and w, first local and then
582	!-- average along the outflow boundary.
583	DO k = nzb+1, nzt+1
584	DO i = nxl, nxr
585
586	denom = u_m_s(k,0,i) - u_m_s(k,1,i)
587
588	IF ( denom /= 0.0_wp ) THEN
589	c_u(k,i) = -c_max * ( u(k,0,i) - u_m_s(k,0,i) ) / ( denom * tsc(2) )
590	IF ( c_u(k,i) < 0.0_wp ) THEN
591	c_u(k,i) = 0.0_wp
592	ELSEIF ( c_u(k,i) > c_max ) THEN
593	c_u(k,i) = c_max
594	ENDIF
595	ELSE
596	c_u(k,i) = c_max
597	ENDIF
598
599	denom = v_m_s(k,1,i) - v_m_s(k,2,i)
600
601	IF ( denom /= 0.0_wp ) THEN
602	c_v(k,i) = -c_max * ( v(k,1,i) - v_m_s(k,1,i) ) / ( denom * tsc(2) )
603	IF ( c_v(k,i) < 0.0_wp ) THEN
604	c_v(k,i) = 0.0_wp
605	ELSEIF ( c_v(k,i) > c_max ) THEN
606	c_v(k,i) = c_max
607	ENDIF
608	ELSE
609	c_v(k,i) = c_max
610	ENDIF
611
612	denom = w_m_s(k,0,i) - w_m_s(k,1,i)
613
614	IF ( denom /= 0.0_wp ) THEN
615	c_w(k,i) = -c_max * ( w(k,0,i) - w_m_s(k,0,i) ) / ( denom * tsc(2) )
616	IF ( c_w(k,i) < 0.0_wp ) THEN
617	c_w(k,i) = 0.0_wp
618	ELSEIF ( c_w(k,i) > c_max ) THEN
619	c_w(k,i) = c_max
620	ENDIF
621	ELSE
622	c_w(k,i) = c_max
623	ENDIF
624
625	c_u_m_l(k) = c_u_m_l(k) + c_u(k,i)
626	c_v_m_l(k) = c_v_m_l(k) + c_v(k,i)
627	c_w_m_l(k) = c_w_m_l(k) + c_w(k,i)
628
629	ENDDO
630	ENDDO
631
632	#if defined( __parallel )
633	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
634	CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
635	MPI_SUM, comm1dx, ierr )
636	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
637	CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
638	MPI_SUM, comm1dx, ierr )
639	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
640	CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
641	MPI_SUM, comm1dx, ierr )
642	#else
643	c_u_m = c_u_m_l
644	c_v_m = c_v_m_l
645	c_w_m = c_w_m_l
646	#endif
647
648	c_u_m = c_u_m / (nx+1)
649	c_v_m = c_v_m / (nx+1)
650	c_w_m = c_w_m / (nx+1)
651
652	!
653	!-- Save old timelevels for the next timestep
654	IF ( intermediate_timestep_count == 1 ) THEN
655	u_m_s(:,:,:) = u(:,0:1,:)
656	v_m_s(:,:,:) = v(:,1:2,:)
657	w_m_s(:,:,:) = w(:,0:1,:)
658	ENDIF
659
660	!
661	!-- Calculate the new velocities
662	DO k = nzb+1, nzt+1
663	DO i = nxlg, nxrg
664	u_p(k,-1,i) = u(k,-1,i) - dt_3d * tsc(2) * c_u_m(k) * &
665	( u(k,-1,i) - u(k,0,i) ) * ddy
666
667	v_p(k,0,i) = v(k,0,i) - dt_3d * tsc(2) * c_v_m(k) * &
668	( v(k,0,i) - v(k,1,i) ) * ddy
669
670	w_p(k,-1,i) = w(k,-1,i) - dt_3d * tsc(2) * c_w_m(k) * &
671	( w(k,-1,i) - w(k,0,i) ) * ddy
672	ENDDO
673	ENDDO
674
675	!
676	!-- Bottom boundary at the outflow
677	IF ( ibc_uv_b == 0 ) THEN
678	u_p(nzb,-1,:) = 0.0_wp
679	v_p(nzb,0,:) = 0.0_wp
680	ELSE
681	u_p(nzb,-1,:) = u_p(nzb+1,-1,:)
682	v_p(nzb,0,:) = v_p(nzb+1,0,:)
683	ENDIF
684	w_p(nzb,-1,:) = 0.0_wp
685
686	!
687	!-- Top boundary at the outflow
688	IF ( ibc_uv_t == 0 ) THEN
689	u_p(nzt+1,-1,:) = u_init(nzt+1)
690	v_p(nzt+1,0,:) = v_init(nzt+1)
691	ELSE
692	u_p(nzt+1,-1,:) = u_p(nzt,-1,:)
693	v_p(nzt+1,0,:) = v_p(nzt,0,:)
694	ENDIF
695	w_p(nzt:nzt+1,-1,:) = 0.0_wp
696
697	ENDIF
698
699	ENDIF
700
701	IF ( outflow_n ) THEN
702
703	IF ( use_cmax ) THEN
704	u_p(:,ny+1,:) = u(:,ny,:)
705	v_p(:,ny+1,:) = v(:,ny,:)
706	w_p(:,ny+1,:) = w(:,ny,:)
707	ELSEIF ( .NOT. use_cmax ) THEN
708
709	c_max = dy / dt_3d
710
711	c_u_m_l = 0.0_wp
712	c_v_m_l = 0.0_wp
713	c_w_m_l = 0.0_wp
714
715	c_u_m = 0.0_wp
716	c_v_m = 0.0_wp
717	c_w_m = 0.0_wp
718
719	!
720	!-- Calculate the phase speeds for u, v, and w, first local and then
721	!-- average along the outflow boundary.
722	DO k = nzb+1, nzt+1
723	DO i = nxl, nxr
724
725	denom = u_m_n(k,ny,i) - u_m_n(k,ny-1,i)
726
727	IF ( denom /= 0.0_wp ) THEN
728	c_u(k,i) = -c_max * ( u(k,ny,i) - u_m_n(k,ny,i) ) / ( denom * tsc(2) )
729	IF ( c_u(k,i) < 0.0_wp ) THEN
730	c_u(k,i) = 0.0_wp
731	ELSEIF ( c_u(k,i) > c_max ) THEN
732	c_u(k,i) = c_max
733	ENDIF
734	ELSE
735	c_u(k,i) = c_max
736	ENDIF
737
738	denom = v_m_n(k,ny,i) - v_m_n(k,ny-1,i)
739
740	IF ( denom /= 0.0_wp ) THEN
741	c_v(k,i) = -c_max * ( v(k,ny,i) - v_m_n(k,ny,i) ) / ( denom * tsc(2) )
742	IF ( c_v(k,i) < 0.0_wp ) THEN
743	c_v(k,i) = 0.0_wp
744	ELSEIF ( c_v(k,i) > c_max ) THEN
745	c_v(k,i) = c_max
746	ENDIF
747	ELSE
748	c_v(k,i) = c_max
749	ENDIF
750
751	denom = w_m_n(k,ny,i) - w_m_n(k,ny-1,i)
752
753	IF ( denom /= 0.0_wp ) THEN
754	c_w(k,i) = -c_max * ( w(k,ny,i) - w_m_n(k,ny,i) ) / ( denom * tsc(2) )
755	IF ( c_w(k,i) < 0.0_wp ) THEN
756	c_w(k,i) = 0.0_wp
757	ELSEIF ( c_w(k,i) > c_max ) THEN
758	c_w(k,i) = c_max
759	ENDIF
760	ELSE
761	c_w(k,i) = c_max
762	ENDIF
763
764	c_u_m_l(k) = c_u_m_l(k) + c_u(k,i)
765	c_v_m_l(k) = c_v_m_l(k) + c_v(k,i)
766	c_w_m_l(k) = c_w_m_l(k) + c_w(k,i)
767
768	ENDDO
769	ENDDO
770
771	#if defined( __parallel )
772	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
773	CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
774	MPI_SUM, comm1dx, ierr )
775	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
776	CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
777	MPI_SUM, comm1dx, ierr )
778	IF ( collective_wait ) CALL MPI_BARRIER( comm1dx, ierr )
779	CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
780	MPI_SUM, comm1dx, ierr )
781	#else
782	c_u_m = c_u_m_l
783	c_v_m = c_v_m_l
784	c_w_m = c_w_m_l
785	#endif
786
787	c_u_m = c_u_m / (nx+1)
788	c_v_m = c_v_m / (nx+1)
789	c_w_m = c_w_m / (nx+1)
790
791	!
792	!-- Save old timelevels for the next timestep
793	IF ( intermediate_timestep_count == 1 ) THEN
794	u_m_n(:,:,:) = u(:,ny-1:ny,:)
795	v_m_n(:,:,:) = v(:,ny-1:ny,:)
796	w_m_n(:,:,:) = w(:,ny-1:ny,:)
797	ENDIF
798
799	!
800	!-- Calculate the new velocities
801	DO k = nzb+1, nzt+1
802	DO i = nxlg, nxrg
803	u_p(k,ny+1,i) = u(k,ny+1,i) - dt_3d * tsc(2) * c_u_m(k) * &
804	( u(k,ny+1,i) - u(k,ny,i) ) * ddy
805
806	v_p(k,ny+1,i) = v(k,ny+1,i) - dt_3d * tsc(2) * c_v_m(k) * &
807	( v(k,ny+1,i) - v(k,ny,i) ) * ddy
808
809	w_p(k,ny+1,i) = w(k,ny+1,i) - dt_3d * tsc(2) * c_w_m(k) * &
810	( w(k,ny+1,i) - w(k,ny,i) ) * ddy
811	ENDDO
812	ENDDO
813
814	!
815	!-- Bottom boundary at the outflow
816	IF ( ibc_uv_b == 0 ) THEN
817	u_p(nzb,ny+1,:) = 0.0_wp
818	v_p(nzb,ny+1,:) = 0.0_wp
819	ELSE
820	u_p(nzb,ny+1,:) = u_p(nzb+1,ny+1,:)
821	v_p(nzb,ny+1,:) = v_p(nzb+1,ny+1,:)
822	ENDIF
823	w_p(nzb,ny+1,:) = 0.0_wp
824
825	!
826	!-- Top boundary at the outflow
827	IF ( ibc_uv_t == 0 ) THEN
828	u_p(nzt+1,ny+1,:) = u_init(nzt+1)
829	v_p(nzt+1,ny+1,:) = v_init(nzt+1)
830	ELSE
831	u_p(nzt+1,ny+1,:) = u_p(nzt,nyn+1,:)
832	v_p(nzt+1,ny+1,:) = v_p(nzt,nyn+1,:)
833	ENDIF
834	w_p(nzt:nzt+1,ny+1,:) = 0.0_wp
835
836	ENDIF
837
838	ENDIF
839
840	IF ( outflow_l ) THEN
841
842	IF ( use_cmax ) THEN
843	u_p(:,:,0) = u(:,:,1)
844	v_p(:,:,-1) = v(:,:,0)
845	w_p(:,:,-1) = w(:,:,0)
846	ELSEIF ( .NOT. use_cmax ) THEN
847
848	c_max = dx / dt_3d
849
850	c_u_m_l = 0.0_wp
851	c_v_m_l = 0.0_wp
852	c_w_m_l = 0.0_wp
853
854	c_u_m = 0.0_wp
855	c_v_m = 0.0_wp
856	c_w_m = 0.0_wp
857
858	!
859	!-- Calculate the phase speeds for u, v, and w, first local and then
860	!-- average along the outflow boundary.
861	DO k = nzb+1, nzt+1
862	DO j = nys, nyn
863
864	denom = u_m_l(k,j,1) - u_m_l(k,j,2)
865
866	IF ( denom /= 0.0_wp ) THEN
867	c_u(k,j) = -c_max * ( u(k,j,1) - u_m_l(k,j,1) ) / ( denom * tsc(2) )
868	IF ( c_u(k,j) < 0.0_wp ) THEN
869	c_u(k,j) = 0.0_wp
870	ELSEIF ( c_u(k,j) > c_max ) THEN
871	c_u(k,j) = c_max
872	ENDIF
873	ELSE
874	c_u(k,j) = c_max
875	ENDIF
876
877	denom = v_m_l(k,j,0) - v_m_l(k,j,1)
878
879	IF ( denom /= 0.0_wp ) THEN
880	c_v(k,j) = -c_max * ( v(k,j,0) - v_m_l(k,j,0) ) / ( denom * tsc(2) )
881	IF ( c_v(k,j) < 0.0_wp ) THEN
882	c_v(k,j) = 0.0_wp
883	ELSEIF ( c_v(k,j) > c_max ) THEN
884	c_v(k,j) = c_max
885	ENDIF
886	ELSE
887	c_v(k,j) = c_max
888	ENDIF
889
890	denom = w_m_l(k,j,0) - w_m_l(k,j,1)
891
892	IF ( denom /= 0.0_wp ) THEN
893	c_w(k,j) = -c_max * ( w(k,j,0) - w_m_l(k,j,0) ) / ( denom * tsc(2) )
894	IF ( c_w(k,j) < 0.0_wp ) THEN
895	c_w(k,j) = 0.0_wp
896	ELSEIF ( c_w(k,j) > c_max ) THEN
897	c_w(k,j) = c_max
898	ENDIF
899	ELSE
900	c_w(k,j) = c_max
901	ENDIF
902
903	c_u_m_l(k) = c_u_m_l(k) + c_u(k,j)
904	c_v_m_l(k) = c_v_m_l(k) + c_v(k,j)
905	c_w_m_l(k) = c_w_m_l(k) + c_w(k,j)
906
907	ENDDO
908	ENDDO
909
910	#if defined( __parallel )
911	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
912	CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
913	MPI_SUM, comm1dy, ierr )
914	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
915	CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
916	MPI_SUM, comm1dy, ierr )
917	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
918	CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
919	MPI_SUM, comm1dy, ierr )
920	#else
921	c_u_m = c_u_m_l
922	c_v_m = c_v_m_l
923	c_w_m = c_w_m_l
924	#endif
925
926	c_u_m = c_u_m / (ny+1)
927	c_v_m = c_v_m / (ny+1)
928	c_w_m = c_w_m / (ny+1)
929
930	!
931	!-- Save old timelevels for the next timestep
932	IF ( intermediate_timestep_count == 1 ) THEN
933	u_m_l(:,:,:) = u(:,:,1:2)
934	v_m_l(:,:,:) = v(:,:,0:1)
935	w_m_l(:,:,:) = w(:,:,0:1)
936	ENDIF
937
938	!
939	!-- Calculate the new velocities
940	DO k = nzb+1, nzt+1
941	DO j = nysg, nyng
942	u_p(k,j,0) = u(k,j,0) - dt_3d * tsc(2) * c_u_m(k) * &
943	( u(k,j,0) - u(k,j,1) ) * ddx
944
945	v_p(k,j,-1) = v(k,j,-1) - dt_3d * tsc(2) * c_v_m(k) * &
946	( v(k,j,-1) - v(k,j,0) ) * ddx
947
948	w_p(k,j,-1) = w(k,j,-1) - dt_3d * tsc(2) * c_w_m(k) * &
949	( w(k,j,-1) - w(k,j,0) ) * ddx
950	ENDDO
951	ENDDO
952
953	!
954	!-- Bottom boundary at the outflow
955	IF ( ibc_uv_b == 0 ) THEN
956	u_p(nzb,:,0) = 0.0_wp
957	v_p(nzb,:,-1) = 0.0_wp
958	ELSE
959	u_p(nzb,:,0) = u_p(nzb+1,:,0)
960	v_p(nzb,:,-1) = v_p(nzb+1,:,-1)
961	ENDIF
962	w_p(nzb,:,-1) = 0.0_wp
963
964	!
965	!-- Top boundary at the outflow
966	IF ( ibc_uv_t == 0 ) THEN
967	u_p(nzt+1,:,0) = u_init(nzt+1)
968	v_p(nzt+1,:,-1) = v_init(nzt+1)
969	ELSE
970	u_p(nzt+1,:,0) = u_p(nzt,:,0)
971	v_p(nzt+1,:,-1) = v_p(nzt,:,-1)
972	ENDIF
973	w_p(nzt:nzt+1,:,-1) = 0.0_wp
974
975	ENDIF
976
977	ENDIF
978
979	IF ( outflow_r ) THEN
980
981	IF ( use_cmax ) THEN
982	u_p(:,:,nx+1) = u(:,:,nx)
983	v_p(:,:,nx+1) = v(:,:,nx)
984	w_p(:,:,nx+1) = w(:,:,nx)
985	ELSEIF ( .NOT. use_cmax ) THEN
986
987	c_max = dx / dt_3d
988
989	c_u_m_l = 0.0_wp
990	c_v_m_l = 0.0_wp
991	c_w_m_l = 0.0_wp
992
993	c_u_m = 0.0_wp
994	c_v_m = 0.0_wp
995	c_w_m = 0.0_wp
996
997	!
998	!-- Calculate the phase speeds for u, v, and w, first local and then
999	!-- average along the outflow boundary.
1000	DO k = nzb+1, nzt+1
1001	DO j = nys, nyn
1002
1003	denom = u_m_r(k,j,nx) - u_m_r(k,j,nx-1)
1004
1005	IF ( denom /= 0.0_wp ) THEN
1006	c_u(k,j) = -c_max * ( u(k,j,nx) - u_m_r(k,j,nx) ) / ( denom * tsc(2) )
1007	IF ( c_u(k,j) < 0.0_wp ) THEN
1008	c_u(k,j) = 0.0_wp
1009	ELSEIF ( c_u(k,j) > c_max ) THEN
1010	c_u(k,j) = c_max
1011	ENDIF
1012	ELSE
1013	c_u(k,j) = c_max
1014	ENDIF
1015
1016	denom = v_m_r(k,j,nx) - v_m_r(k,j,nx-1)
1017
1018	IF ( denom /= 0.0_wp ) THEN
1019	c_v(k,j) = -c_max * ( v(k,j,nx) - v_m_r(k,j,nx) ) / ( denom * tsc(2) )
1020	IF ( c_v(k,j) < 0.0_wp ) THEN
1021	c_v(k,j) = 0.0_wp
1022	ELSEIF ( c_v(k,j) > c_max ) THEN
1023	c_v(k,j) = c_max
1024	ENDIF
1025	ELSE
1026	c_v(k,j) = c_max
1027	ENDIF
1028
1029	denom = w_m_r(k,j,nx) - w_m_r(k,j,nx-1)
1030
1031	IF ( denom /= 0.0_wp ) THEN
1032	c_w(k,j) = -c_max * ( w(k,j,nx) - w_m_r(k,j,nx) ) / ( denom * tsc(2) )
1033	IF ( c_w(k,j) < 0.0_wp ) THEN
1034	c_w(k,j) = 0.0_wp
1035	ELSEIF ( c_w(k,j) > c_max ) THEN
1036	c_w(k,j) = c_max
1037	ENDIF
1038	ELSE
1039	c_w(k,j) = c_max
1040	ENDIF
1041
1042	c_u_m_l(k) = c_u_m_l(k) + c_u(k,j)
1043	c_v_m_l(k) = c_v_m_l(k) + c_v(k,j)
1044	c_w_m_l(k) = c_w_m_l(k) + c_w(k,j)
1045
1046	ENDDO
1047	ENDDO
1048
1049	#if defined( __parallel )
1050	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
1051	CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
1052	MPI_SUM, comm1dy, ierr )
1053	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
1054	CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
1055	MPI_SUM, comm1dy, ierr )
1056	IF ( collective_wait ) CALL MPI_BARRIER( comm1dy, ierr )
1057	CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
1058	MPI_SUM, comm1dy, ierr )
1059	#else
1060	c_u_m = c_u_m_l
1061	c_v_m = c_v_m_l
1062	c_w_m = c_w_m_l
1063	#endif
1064
1065	c_u_m = c_u_m / (ny+1)
1066	c_v_m = c_v_m / (ny+1)
1067	c_w_m = c_w_m / (ny+1)
1068
1069	!
1070	!-- Save old timelevels for the next timestep
1071	IF ( intermediate_timestep_count == 1 ) THEN
1072	u_m_r(:,:,:) = u(:,:,nx-1:nx)
1073	v_m_r(:,:,:) = v(:,:,nx-1:nx)
1074	w_m_r(:,:,:) = w(:,:,nx-1:nx)
1075	ENDIF
1076
1077	!
1078	!-- Calculate the new velocities
1079	DO k = nzb+1, nzt+1
1080	DO j = nysg, nyng
1081	u_p(k,j,nx+1) = u(k,j,nx+1) - dt_3d * tsc(2) * c_u_m(k) * &
1082	( u(k,j,nx+1) - u(k,j,nx) ) * ddx
1083
1084	v_p(k,j,nx+1) = v(k,j,nx+1) - dt_3d * tsc(2) * c_v_m(k) * &
1085	( v(k,j,nx+1) - v(k,j,nx) ) * ddx
1086
1087	w_p(k,j,nx+1) = w(k,j,nx+1) - dt_3d * tsc(2) * c_w_m(k) * &
1088	( w(k,j,nx+1) - w(k,j,nx) ) * ddx
1089	ENDDO
1090	ENDDO
1091
1092	!
1093	!-- Bottom boundary at the outflow
1094	IF ( ibc_uv_b == 0 ) THEN
1095	u_p(nzb,:,nx+1) = 0.0_wp
1096	v_p(nzb,:,nx+1) = 0.0_wp
1097	ELSE
1098	u_p(nzb,:,nx+1) = u_p(nzb+1,:,nx+1)
1099	v_p(nzb,:,nx+1) = v_p(nzb+1,:,nx+1)
1100	ENDIF
1101	w_p(nzb,:,nx+1) = 0.0_wp
1102
1103	!
1104	!-- Top boundary at the outflow
1105	IF ( ibc_uv_t == 0 ) THEN
1106	u_p(nzt+1,:,nx+1) = u_init(nzt+1)
1107	v_p(nzt+1,:,nx+1) = v_init(nzt+1)
1108	ELSE
1109	u_p(nzt+1,:,nx+1) = u_p(nzt,:,nx+1)
1110	v_p(nzt+1,:,nx+1) = v_p(nzt,:,nx+1)
1111	ENDIF
1112	w_p(nzt:nzt+1,:,nx+1) = 0.0_wp
1113
1114	ENDIF
1115
1116	ENDIF
1117
1118	END SUBROUTINE boundary_conds

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |