Home

Context Navigation

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

Last change on this file since 2000 was 2000, checked in by knoop, 8 years ago
Forced header and separation lines into 80 columns
Property svn:keywords set to `Id`
File size: 33.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |