Home

Context Navigation

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

Last change on this file since 1372 was 1362, checked in by hoffmann, 11 years ago
last commit documented
Property svn:keywords set to `Id`
File size: 30.7 KB

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |