Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |