Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |