Home

Context Navigation

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

Last change on this file since 2703 was 2696, checked in by kanani, 7 years ago
Merge of branch palm4u into trunk
Property svn:keywords set to `Id`
File size: 40.5 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |