Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |