Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |