Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |