Home

Context Navigation

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

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

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

Download in other formats:

| Impressum | ©Leibniz Universität Hannover |