Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |