Home

Context Navigation

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

Last change on this file since 2487 was 2365, checked in by kanani, 7 years ago
Vertical nesting implemented (SadiqHuq?)
Property svn:keywords set to `Id`
File size: 39.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |