Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |