Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |