Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |