Home

Context Navigation

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

Last change on this file since 4125 was 4102, checked in by suehring, 5 years ago
Bugfix, set Neumann boundary conditions for the subgrid TKE at vertical walls instead of implicit Dirichlet conditions that always act as a sink term for the subgrid TKE. Therefore, add new data structure for vertical surfaces and revise the setting of the boundary grid point index space. Moreover, slightly revise setting of boundary conditions at upward- and downward facing surfaces. Finally, setting of boundary conditions for subgrid TKE and dissipation (in RANS mode) is now modularized. Update test case results.
Property svn:keywords set to `Id`
File size: 40.6 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |