Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |