Home

Context Navigation

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

Last change on this file since 2805 was 2766, checked in by kanani, 7 years ago
Removal of chem directive, plus minor changes
Property svn:keywords set to `Id`
File size: 40.6 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 2766 2018-01-22 17:17:47Z suehring $
27	! Removed preprocessor directive __chem
28	!
29	! 2718 2018-01-02 08:49:38Z maronga
30	! Corrected "Former revisions" section
31	!
32	! 2696 2017-12-14 17:12:51Z kanani
33	! Change in file header (GPL part)
34	! Adjust boundary conditions for e and diss in case of TKE-e closure (TG)
35	! Implementation of chemistry module (FK)
36	!
37	! 2569 2017-10-20 11:54:42Z kanani
38	! Removed redundant code for ibc_s_b=1 and ibc_q_b=1
39	!
40	! 2365 2017-08-21 14:59:59Z kanani
41	! Vertical grid nesting implemented: exclude setting vertical velocity to zero
42	! on fine grid (SadiqHuq)
43	!
44	! 2320 2017-07-21 12:47:43Z suehring
45	! Remove unused control parameter large_scale_forcing from only-list
46	!
47	! 2292 2017-06-20 09:51:42Z schwenkel
48	! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
49	! includes two more prognostic equations for cloud drop concentration (nc)
50	! and cloud water content (qc).
51	!
52	! 2233 2017-05-30 18:08:54Z suehring
53	!
54	! 2232 2017-05-30 17:47:52Z suehring
55	! Set boundary conditions on topography top using flag method.
56	!
57	! 2118 2017-01-17 16:38:49Z raasch
58	! OpenACC directives removed
59	!
60	! 2000 2016-08-20 18:09:15Z knoop
61	! Forced header and separation lines into 80 columns
62	!
63	! 1992 2016-08-12 15:14:59Z suehring
64	! Adjustments for top boundary condition for passive scalar
65	!
66	! 1960 2016-07-12 16:34:24Z suehring
67	! Treat humidity and passive scalar separately
68	!
69	! 1823 2016-04-07 08:57:52Z hoffmann
70	! Initial version of purely vertical nesting introduced.
71	!
72	! 1822 2016-04-07 07:49:42Z hoffmann
73	! icloud_scheme removed. microphyisics_seifert added.
74	!
75	! 1764 2016-02-28 12:45:19Z raasch
76	! index bug for u_p at left outflow removed
77	!
78	! 1762 2016-02-25 12:31:13Z hellstea
79	! Introduction of nested domain feature
80	!
81	! 1742 2016-01-13 09:50:06Z raasch
82	! bugfix for outflow Neumann boundary conditions at bottom and top
83	!
84	! 1717 2015-11-11 15:09:47Z raasch
85	! Bugfix: index error in outflow conditions for left boundary
86	!
87	! 1682 2015-10-07 23:56:08Z knoop
88	! Code annotations made doxygen readable
89	!
90	! 1410 2014-05-23 12:16:18Z suehring
91	! Bugfix: set dirichlet boundary condition for passive_scalar at model domain
92	! top
93	!
94	! 1399 2014-05-07 11:16:25Z heinze
95	! Bugfix: set inflow boundary conditions also if no humidity or passive_scalar
96	! is used.
97	!
98	! 1398 2014-05-07 11:15:00Z heinze
99	! Dirichlet-condition at the top for u and v changed to u_init and v_init also
100	! for large_scale_forcing
101	!
102	! 1380 2014-04-28 12:40:45Z heinze
103	! Adjust Dirichlet-condition at the top for pt in case of nudging
104	!
105	! 1361 2014-04-16 15:17:48Z hoffmann
106	! Bottom and top boundary conditions of rain water content (qr) and
107	! rain drop concentration (nr) changed to Dirichlet
108	!
109	! 1353 2014-04-08 15:21:23Z heinze
110	! REAL constants provided with KIND-attribute
111	!
112	! 1320 2014-03-20 08:40:49Z raasch
113	! ONLY-attribute added to USE-statements,
114	! kind-parameters added to all INTEGER and REAL declaration statements,
115	! kinds are defined in new module kinds,
116	! revision history before 2012 removed,
117	! comment fields (!:) to be used for variable explanations added to
118	! all variable declaration statements
119	!
120	! 1257 2013-11-08 15:18:40Z raasch
121	! loop independent clauses added
122	!
123	! 1241 2013-10-30 11:36:58Z heinze
124	! Adjust ug and vg at each timestep in case of large_scale_forcing
125	!
126	! 1159 2013-05-21 11:58:22Z fricke
127	! Bugfix: Neumann boundary conditions for the velocity components at the
128	! outflow are in fact radiation boundary conditions using the maximum phase
129	! velocity that ensures numerical stability (CFL-condition).
130	! Hence, logical operator use_cmax is now used instead of bc_lr_dirneu/_neudir.
131	! Bugfix: In case of use_cmax at the outflow, u, v, w are replaced by
132	! u_p, v_p, w_p
133	!
134	! 1115 2013-03-26 18:16:16Z hoffmann
135	! boundary conditions of two-moment cloud scheme are restricted to Neumann-
136	! boundary-conditions
137	!
138	! 1113 2013-03-10 02:48:14Z raasch
139	! GPU-porting
140	! dummy argument "range" removed
141	! Bugfix: wrong index in loops of radiation boundary condition
142	!
143	! 1053 2012-11-13 17:11:03Z hoffmann
144	! boundary conditions for the two new prognostic equations (nr, qr) of the
145	! two-moment cloud scheme
146	!
147	! 1036 2012-10-22 13:43:42Z raasch
148	! code put under GPL (PALM 3.9)
149	!
150	! 996 2012-09-07 10:41:47Z raasch
151	! little reformatting
152	!
153	! 978 2012-08-09 08:28:32Z fricke
154	! Neumann boudnary conditions are added at the inflow boundary for the SGS-TKE.
155	! Outflow boundary conditions for the velocity components can be set to Neumann
156	! conditions or to radiation conditions with a horizontal averaged phase
157	! velocity.
158	!
159	! 875 2012-04-02 15:35:15Z gryschka
160	! Bugfix in case of dirichlet inflow bc at the right or north boundary
161	!
162	! Revision 1.1 1997/09/12 06:21:34 raasch
163	! Initial revision
164	!
165	!
166	! Description:
167	! ------------
168	!> Boundary conditions for the prognostic quantities.
169	!> One additional bottom boundary condition is applied for the TKE (=(u)*2)
170	!> in prandtl_fluxes. The cyclic lateral boundary conditions are implicitly
171	!> handled in routine exchange_horiz. Pressure boundary conditions are
172	!> explicitly set in routines pres, poisfft, poismg and sor.
173	!------------------------------------------------------------------------------!
174	SUBROUTINE boundary_conds
175
176
177	USE arrays_3d, &
178	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, &
179	diss_p, dzu, e_p, nc_p, nr_p, pt, pt_p, q, q_p, qc_p, qr_p, s, &
180	s_p, sa, sa_p, u, ug, u_init, u_m_l, u_m_n, u_m_r, u_m_s, u_p, &
181	v, vg, v_init, v_m_l, v_m_n, v_m_r, v_m_s, v_p, &
182	w, w_p, w_m_l, w_m_n, w_m_r, w_m_s, pt_init
183
184	USE chemistry_model_mod, &
185	ONLY: chem_boundary_conds
186
187	USE control_parameters, &
188	ONLY: air_chemistry, bc_pt_t_val, bc_q_t_val, bc_s_t_val, &
189	constant_diffusion, cloud_physics, coupling_mode, dt_3d, &
190	force_bound_l, force_bound_s, forcing, humidity, &
191	ibc_pt_b, ibc_pt_t, ibc_q_b, ibc_q_t, ibc_s_b, ibc_s_t, &
192	ibc_sa_t, ibc_uv_b, ibc_uv_t, inflow_l, inflow_n, inflow_r, &
193	inflow_s, intermediate_timestep_count, &
194	microphysics_morrison, microphysics_seifert, nest_domain, &
195	nest_bound_l, nest_bound_s, nudging, ocean, outflow_l, &
196	outflow_n, outflow_r, outflow_s, passive_scalar, rans_tke_e, &
197	tsc, use_cmax
198
199	USE grid_variables, &
200	ONLY: ddx, ddy, dx, dy
201
202	USE indices, &
203	ONLY: nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg, &
204	nzb, nzt, wall_flags_0
205
206	USE kinds
207
208	USE pegrid
209
210	USE pmc_interface, &
211	ONLY : nesting_mode
212
213	USE surface_mod, &
214	ONLY : bc_h
215
216	IMPLICIT NONE
217
218	INTEGER(iwp) :: i !< grid index x direction
219	INTEGER(iwp) :: j !< grid index y direction
220	INTEGER(iwp) :: k !< grid index z direction
221	INTEGER(iwp) :: kb !< variable to set respective boundary value, depends on facing.
222	INTEGER(iwp) :: l !< running index boundary type, for up- and downward-facing walls
223	INTEGER(iwp) :: m !< running index surface elements
224
225	REAL(wp) :: c_max !<
226	REAL(wp) :: denom !<
227
228
229	!
230	!-- Bottom boundary
231	IF ( ibc_uv_b == 1 ) THEN
232	u_p(nzb,:,:) = u_p(nzb+1,:,:)
233	v_p(nzb,:,:) = v_p(nzb+1,:,:)
234	ENDIF
235	!
236	!-- Set zero vertical velocity at topography top (l=0), or bottom (l=1) in case
237	!-- of downward-facing surfaces.
238	DO l = 0, 1
239	!
240	!-- Set kb, for upward-facing surfaces value at topography top (k-1) is set,
241	!-- for downward-facing surfaces at topography bottom (k+1).
242	kb = MERGE( -1, 1, l == 0 )
243	!$OMP PARALLEL DO PRIVATE( i, j, k )
244	DO m = 1, bc_h(l)%ns
245	i = bc_h(l)%i(m)
246	j = bc_h(l)%j(m)
247	k = bc_h(l)%k(m)
248	w_p(k+kb,j,i) = 0.0_wp
249	ENDDO
250	ENDDO
251
252	!
253	!-- Top boundary. A nested domain ( ibc_uv_t = 3 ) does not require settings.
254	IF ( ibc_uv_t == 0 ) THEN
255	u_p(nzt+1,:,:) = u_init(nzt+1)
256	v_p(nzt+1,:,:) = v_init(nzt+1)
257	ELSEIF ( ibc_uv_t == 1 ) THEN
258	u_p(nzt+1,:,:) = u_p(nzt,:,:)
259	v_p(nzt+1,:,:) = v_p(nzt,:,:)
260	ENDIF
261
262	!
263	!-- Vertical nesting: Vertical velocity not zero at the top of the fine grid
264	IF ( .NOT. nest_domain .AND. &
265	TRIM(coupling_mode) /= 'vnested_fine' ) THEN
266	w_p(nzt:nzt+1,:,:) = 0.0_wp !< nzt is not a prognostic level (but cf. pres)
267	ENDIF
268
269	!
270	!-- Temperature at bottom and top boundary.
271	!-- In case of coupled runs (ibc_pt_b = 2) the temperature is given by
272	!-- the sea surface temperature of the coupled ocean model.
273	!-- Dirichlet
274	IF ( ibc_pt_b == 0 ) THEN
275	DO l = 0, 1
276	!
277	!-- Set kb, for upward-facing surfaces value at topography top (k-1) is set,
278	!-- for downward-facing surfaces at topography bottom (k+1).
279	kb = MERGE( -1, 1, l == 0 )
280	!$OMP PARALLEL DO PRIVATE( i, j, k )
281	DO m = 1, bc_h(l)%ns
282	i = bc_h(l)%i(m)
283	j = bc_h(l)%j(m)
284	k = bc_h(l)%k(m)
285	pt_p(k+kb,j,i) = pt(k+kb,j,i)
286	ENDDO
287	ENDDO
288	!
289	!-- Neumann, zero-gradient
290	ELSEIF ( ibc_pt_b == 1 ) THEN
291	DO l = 0, 1
292	!
293	!-- Set kb, for upward-facing surfaces value at topography top (k-1) is set,
294	!-- for downward-facing surfaces at topography bottom (k+1).
295	kb = MERGE( -1, 1, l == 0 )
296	!$OMP PARALLEL DO PRIVATE( i, j, k )
297	DO m = 1, bc_h(l)%ns
298	i = bc_h(l)%i(m)
299	j = bc_h(l)%j(m)
300	k = bc_h(l)%k(m)
301	pt_p(k+kb,j,i) = pt_p(k,j,i)
302	ENDDO
303	ENDDO
304	ENDIF
305
306	!
307	!-- Temperature at top boundary
308	IF ( ibc_pt_t == 0 ) THEN
309	pt_p(nzt+1,:,:) = pt(nzt+1,:,:)
310	!
311	!-- In case of nudging adjust top boundary to pt which is
312	!-- read in from NUDGING-DATA
313	IF ( nudging ) THEN
314	pt_p(nzt+1,:,:) = pt_init(nzt+1)
315	ENDIF
316	ELSEIF ( ibc_pt_t == 1 ) THEN
317	pt_p(nzt+1,:,:) = pt_p(nzt,:,:)
318	ELSEIF ( ibc_pt_t == 2 ) THEN
319	pt_p(nzt+1,:,:) = pt_p(nzt,:,:) + bc_pt_t_val * dzu(nzt+1)
320	ENDIF
321
322	!
323	!-- Boundary conditions for TKE
324	!-- Generally Neumann conditions with de/dz=0 are assumed
325	IF ( .NOT. constant_diffusion ) THEN
326
327	IF ( .NOT. rans_tke_e ) THEN
328	DO l = 0, 1
329	!
330	!-- Set kb, for upward-facing surfaces value at topography top (k-1) is set,
331	!-- for downward-facing surfaces at topography bottom (k+1).
332	kb = MERGE( -1, 1, l == 0 )
333	!$OMP PARALLEL DO PRIVATE( i, j, k )
334	DO m = 1, bc_h(l)%ns
335	i = bc_h(l)%i(m)
336	j = bc_h(l)%j(m)
337	k = bc_h(l)%k(m)
338	e_p(k+kb,j,i) = e_p(k,j,i)
339	ENDDO
340	ENDDO
341	ENDIF
342
343	IF ( .NOT. nest_domain ) THEN
344	e_p(nzt+1,:,:) = e_p(nzt,:,:)
345	ENDIF
346	ENDIF
347
348	!
349	!-- Boundary conditions for TKE dissipation rate
350	IF ( rans_tke_e .AND. .NOT. nest_domain ) THEN
351	diss_p(nzt+1,:,:) = diss_p(nzt,:,:)
352	ENDIF
353
354	!
355	!-- Boundary conditions for salinity
356	IF ( ocean ) THEN
357	!
358	!-- Bottom boundary: Neumann condition because salinity flux is always
359	!-- given.
360	DO l = 0, 1
361	!
362	!-- Set kb, for upward-facing surfaces value at topography top (k-1) is set,
363	!-- for downward-facing surfaces at topography bottom (k+1).
364	kb = MERGE( -1, 1, l == 0 )
365	!$OMP PARALLEL DO PRIVATE( i, j, k )
366	DO m = 1, bc_h(l)%ns
367	i = bc_h(l)%i(m)
368	j = bc_h(l)%j(m)
369	k = bc_h(l)%k(m)
370	sa_p(k+kb,j,i) = sa_p(k,j,i)
371	ENDDO
372	ENDDO
373	!
374	!-- Top boundary: Dirichlet or Neumann
375	IF ( ibc_sa_t == 0 ) THEN
376	sa_p(nzt+1,:,:) = sa(nzt+1,:,:)
377	ELSEIF ( ibc_sa_t == 1 ) THEN
378	sa_p(nzt+1,:,:) = sa_p(nzt,:,:)
379	ENDIF
380
381	ENDIF
382
383	!
384	!-- Boundary conditions for total water content,
385	!-- bottom and top boundary (see also temperature)
386	IF ( humidity ) THEN
387	!
388	!-- Surface conditions for constant_humidity_flux
389	!-- Run loop over all non-natural and natural walls. Note, in wall-datatype
390	!-- the k coordinate belongs to the atmospheric grid point, therefore, set
391	!-- q_p at k-1
392	IF ( ibc_q_b == 0 ) THEN
393
394	DO l = 0, 1
395	!
396	!-- Set kb, for upward-facing surfaces value at topography top (k-1) is set,
397	!-- for downward-facing surfaces at topography bottom (k+1).
398	kb = MERGE( -1, 1, l == 0 )
399	!$OMP PARALLEL DO PRIVATE( i, j, k )
400	DO m = 1, bc_h(l)%ns
401	i = bc_h(l)%i(m)
402	j = bc_h(l)%j(m)
403	k = bc_h(l)%k(m)
404	q_p(k+kb,j,i) = q(k+kb,j,i)
405	ENDDO
406	ENDDO
407
408	ELSE
409
410	DO l = 0, 1
411	!
412	!-- Set kb, for upward-facing surfaces value at topography top (k-1) is set,
413	!-- for downward-facing surfaces at topography bottom (k+1).
414	kb = MERGE( -1, 1, l == 0 )
415	!$OMP PARALLEL DO PRIVATE( i, j, k )
416	DO m = 1, bc_h(l)%ns
417	i = bc_h(l)%i(m)
418	j = bc_h(l)%j(m)
419	k = bc_h(l)%k(m)
420	q_p(k+kb,j,i) = q_p(k,j,i)
421	ENDDO
422	ENDDO
423	ENDIF
424	!
425	!-- Top boundary
426	IF ( ibc_q_t == 0 ) THEN
427	q_p(nzt+1,:,:) = q(nzt+1,:,:)
428	ELSEIF ( ibc_q_t == 1 ) THEN
429	q_p(nzt+1,:,:) = q_p(nzt,:,:) + bc_q_t_val * dzu(nzt+1)
430	ENDIF
431
432	IF ( cloud_physics .AND. microphysics_morrison ) THEN
433	!
434	!-- Surface conditions cloud water (Dirichlet)
435	!-- Run loop over all non-natural and natural walls. Note, in wall-datatype
436	!-- the k coordinate belongs to the atmospheric grid point, therefore, set
437	!-- qr_p and nr_p at k-1
438	!$OMP PARALLEL DO PRIVATE( i, j, k )
439	DO m = 1, bc_h(0)%ns
440	i = bc_h(0)%i(m)
441	j = bc_h(0)%j(m)
442	k = bc_h(0)%k(m)
443	qc_p(k-1,j,i) = 0.0_wp
444	nc_p(k-1,j,i) = 0.0_wp
445	ENDDO
446	!
447	!-- Top boundary condition for cloud water (Dirichlet)
448	qc_p(nzt+1,:,:) = 0.0_wp
449	nc_p(nzt+1,:,:) = 0.0_wp
450
451	ENDIF
452
453	IF ( cloud_physics .AND. microphysics_seifert ) THEN
454	!
455	!-- Surface conditions rain water (Dirichlet)
456	!-- Run loop over all non-natural and natural walls. Note, in wall-datatype
457	!-- the k coordinate belongs to the atmospheric grid point, therefore, set
458	!-- qr_p and nr_p at k-1
459	!$OMP PARALLEL DO PRIVATE( i, j, k )
460	DO m = 1, bc_h(0)%ns
461	i = bc_h(0)%i(m)
462	j = bc_h(0)%j(m)
463	k = bc_h(0)%k(m)
464	qr_p(k-1,j,i) = 0.0_wp
465	nr_p(k-1,j,i) = 0.0_wp
466	ENDDO
467	!
468	!-- Top boundary condition for rain water (Dirichlet)
469	qr_p(nzt+1,:,:) = 0.0_wp
470	nr_p(nzt+1,:,:) = 0.0_wp
471
472	ENDIF
473	ENDIF
474	!
475	!-- Boundary conditions for scalar,
476	!-- bottom and top boundary (see also temperature)
477	IF ( passive_scalar ) THEN
478	!
479	!-- Surface conditions for constant_humidity_flux
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	!-- s_p at k-1
483	IF ( ibc_s_b == 0 ) THEN
484
485	DO l = 0, 1
486	!
487	!-- Set kb, for upward-facing surfaces value at topography top (k-1) is set,
488	!-- for downward-facing surfaces at topography bottom (k+1).
489	kb = MERGE( -1, 1, l == 0 )
490	!$OMP PARALLEL DO PRIVATE( i, j, k )
491	DO m = 1, bc_h(l)%ns
492	i = bc_h(l)%i(m)
493	j = bc_h(l)%j(m)
494	k = bc_h(l)%k(m)
495	s_p(k+kb,j,i) = s(k+kb,j,i)
496	ENDDO
497	ENDDO
498
499	ELSE
500
501	DO l = 0, 1
502	!
503	!-- Set kb, for upward-facing surfaces value at topography top (k-1) is set,
504	!-- for downward-facing surfaces at topography bottom (k+1).
505	kb = MERGE( -1, 1, l == 0 )
506	!$OMP PARALLEL DO PRIVATE( i, j, k )
507	DO m = 1, bc_h(l)%ns
508	i = bc_h(l)%i(m)
509	j = bc_h(l)%j(m)
510	k = bc_h(l)%k(m)
511	s_p(k+kb,j,i) = s_p(k,j,i)
512	ENDDO
513	ENDDO
514	ENDIF
515	!
516	!-- Top boundary condition
517	IF ( ibc_s_t == 0 ) THEN
518	s_p(nzt+1,:,:) = s(nzt+1,:,:)
519	ELSEIF ( ibc_s_t == 1 ) THEN
520	s_p(nzt+1,:,:) = s_p(nzt,:,:)
521	ELSEIF ( ibc_s_t == 2 ) THEN
522	s_p(nzt+1,:,:) = s_p(nzt,:,:) + bc_s_t_val * dzu(nzt+1)
523	ENDIF
524
525	ENDIF
526	!
527	!-- Top/bottom boundary conditions for chemical species
528	IF ( air_chemistry ) CALL chem_boundary_conds( 'set_bc_bottomtop' )
529	!
530	!-- In case of inflow or nest boundary at the south boundary the boundary for v
531	!-- is at nys and in case of inflow or nest boundary at the left boundary the
532	!-- boundary for u is at nxl. Since in prognostic_equations (cache optimized
533	!-- version) these levels are handled as a prognostic level, boundary values
534	!-- have to be restored here.
535	!-- For the SGS-TKE, Neumann boundary conditions are used at the inflow.
536	IF ( inflow_s ) THEN
537	v_p(:,nys,:) = v_p(:,nys-1,:)
538	IF ( .NOT. constant_diffusion ) e_p(:,nys-1,:) = e_p(:,nys,:)
539	ELSEIF ( inflow_n ) THEN
540	IF ( .NOT. constant_diffusion ) e_p(:,nyn+1,:) = e_p(:,nyn,:)
541	ELSEIF ( inflow_l ) THEN
542	u_p(:,:,nxl) = u_p(:,:,nxl-1)
543	IF ( .NOT. constant_diffusion ) e_p(:,:,nxl-1) = e_p(:,:,nxl)
544	ELSEIF ( inflow_r ) THEN
545	IF ( .NOT. constant_diffusion ) e_p(:,:,nxr+1) = e_p(:,:,nxr)
546	ENDIF
547
548	!
549	!-- The same restoration for u at i=nxl and v at j=nys as above must be made
550	!-- in case of nest boundaries. This must not be done in case of vertical nesting
551	!-- mode as in that case the lateral boundaries are actually cyclic.
552	IF ( nesting_mode /= 'vertical' .OR. forcing ) THEN
553	IF ( nest_bound_s .OR. force_bound_s ) THEN
554	v_p(:,nys,:) = v_p(:,nys-1,:)
555	ENDIF
556	IF ( nest_bound_l .OR. force_bound_l ) THEN
557	u_p(:,:,nxl) = u_p(:,:,nxl-1)
558	ENDIF
559	ENDIF
560
561	!
562	!-- Lateral boundary conditions for scalar quantities at the outflow
563	IF ( outflow_s ) THEN
564	pt_p(:,nys-1,:) = pt_p(:,nys,:)
565	IF ( .NOT. constant_diffusion ) e_p(:,nys-1,:) = e_p(:,nys,:)
566	IF ( rans_tke_e ) diss_p(:,nys-1,:) = diss_p(:,nys,:)
567	IF ( humidity ) THEN
568	q_p(:,nys-1,:) = q_p(:,nys,:)
569	IF ( cloud_physics .AND. microphysics_morrison ) THEN
570	qc_p(:,nys-1,:) = qc_p(:,nys,:)
571	nc_p(:,nys-1,:) = nc_p(:,nys,:)
572	ENDIF
573	IF ( cloud_physics .AND. microphysics_seifert ) THEN
574	qr_p(:,nys-1,:) = qr_p(:,nys,:)
575	nr_p(:,nys-1,:) = nr_p(:,nys,:)
576	ENDIF
577	ENDIF
578	IF ( passive_scalar ) s_p(:,nys-1,:) = s_p(:,nys,:)
579	ELSEIF ( outflow_n ) THEN
580	pt_p(:,nyn+1,:) = pt_p(:,nyn,:)
581	IF ( .NOT. constant_diffusion ) e_p(:,nyn+1,:) = e_p(:,nyn,:)
582	IF ( rans_tke_e ) diss_p(:,nyn+1,:) = diss_p(:,nyn,:)
583	IF ( humidity ) THEN
584	q_p(:,nyn+1,:) = q_p(:,nyn,:)
585	IF ( cloud_physics .AND. microphysics_morrison ) THEN
586	qc_p(:,nyn+1,:) = qc_p(:,nyn,:)
587	nc_p(:,nyn+1,:) = nc_p(:,nyn,:)
588	ENDIF
589	IF ( cloud_physics .AND. microphysics_seifert ) THEN
590	qr_p(:,nyn+1,:) = qr_p(:,nyn,:)
591	nr_p(:,nyn+1,:) = nr_p(:,nyn,:)
592	ENDIF
593	ENDIF
594	IF ( passive_scalar ) s_p(:,nyn+1,:) = s_p(:,nyn,:)
595	ELSEIF ( outflow_l ) THEN
596	pt_p(:,:,nxl-1) = pt_p(:,:,nxl)
597	IF ( .NOT. constant_diffusion ) e_p(:,:,nxl-1) = e_p(:,:,nxl)
598	IF ( rans_tke_e ) diss_p(:,:,nxl-1) = diss_p(:,:,nxl)
599	IF ( humidity ) THEN
600	q_p(:,:,nxl-1) = q_p(:,:,nxl)
601	IF ( cloud_physics .AND. microphysics_morrison ) THEN
602	qc_p(:,:,nxl-1) = qc_p(:,:,nxl)
603	nc_p(:,:,nxl-1) = nc_p(:,:,nxl)
604	ENDIF
605	IF ( cloud_physics .AND. microphysics_seifert ) THEN
606	qr_p(:,:,nxl-1) = qr_p(:,:,nxl)
607	nr_p(:,:,nxl-1) = nr_p(:,:,nxl)
608	ENDIF
609	ENDIF
610	IF ( passive_scalar ) s_p(:,:,nxl-1) = s_p(:,:,nxl)
611	ELSEIF ( outflow_r ) THEN
612	pt_p(:,:,nxr+1) = pt_p(:,:,nxr)
613	IF ( .NOT. constant_diffusion ) e_p(:,:,nxr+1) = e_p(:,:,nxr)
614	IF ( rans_tke_e ) diss_p(:,:,nxr+1) = diss_p(:,:,nxr)
615	IF ( humidity ) THEN
616	q_p(:,:,nxr+1) = q_p(:,:,nxr)
617	IF ( cloud_physics .AND. microphysics_morrison ) THEN
618	qc_p(:,:,nxr+1) = qc_p(:,:,nxr)
619	nc_p(:,:,nxr+1) = nc_p(:,:,nxr)
620	ENDIF
621	IF ( cloud_physics .AND. microphysics_seifert ) THEN
622	qr_p(:,:,nxr+1) = qr_p(:,:,nxr)
623	nr_p(:,:,nxr+1) = nr_p(:,:,nxr)
624	ENDIF
625	ENDIF
626	IF ( passive_scalar ) s_p(:,:,nxr+1) = s_p(:,:,nxr)
627	ENDIF
628
629	!
630	!-- Lateral boundary conditions for chemical species
631	IF ( air_chemistry ) CALL chem_boundary_conds( 'set_bc_lateral' )
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 |