Home

Context Navigation

boundary_conds.f90 @ 1257

Last change on this file since 1257 was 1257, checked in by raasch, 10 years ago

New:
---

openACC porting of timestep calculation
(modules, timestep, time_integration)

Changed:

openACC loop directives and vector clauses removed (because they do not give any performance improvement with PGI
compiler versions > 13.6)
(advec_ws, buoyancy, coriolis, diffusion_e, diffusion_s, diffusion_u, diffusion_v, diffusion_w, diffusivities, exchange_horiz, fft_xy, pres, production_e, transpose, tridia_solver, wall_fluxes)

openACC loop independent clauses added
(boundary_conds, prandtl_fluxes, pres)

openACC declare create statements moved after FORTRAN declaration statement
(diffusion_u, diffusion_v, diffusion_w, fft_xy, poisfft, production_e, tridia_solver)

openACC end parallel replaced by end parallel loop
(flow_statistics, pres)

openACC "kernels do" replaced by "kernels loop"
(prandtl_fluxes)

output format for theta* changed to avoid output of *
(run_control)

Errors:

bugfix for calculation of advective timestep (old version may cause wrong timesteps in case of
vertixcally stretched grids)
Attention: standard run-control output has changed!
(timestep)

Property svn:keywords set to Id

File size: 30.2 KB

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

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

Context Navigation

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

Download in other formats: