Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |