Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |