Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |