Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |