Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |