Home

Context Navigation

boundary_conds.f90 @ 91

Last change on this file since 91 was 77, checked in by raasch, 18 years ago

New:
---

particle reflection from vertical walls implemented, particle SGS model adjusted to walls

Wall functions for vertical walls now include diabatic conditions. New subroutines wall_fluxes, wall_fluxes_e. New 4D-array rif_wall.

new d3par-parameter netcdf_64bit_3d to switch on 64bit offset only for 3D files

new d3par-parameter dt_max to define the maximum value for the allowed timestep

new inipar-parameter loop_optimization to control the loop optimization method

new inipar-parameter pt_refrence. If given, this value is used as the reference that in buoyancy terms (otherwise, the instantaneous horizontally averaged temperature is used).

new user interface user_advec_particles

new initializing action "by_user" calls user_init_3d_model and allows the initial setting of all 3d arrays

topography height informations are stored on arrays zu_s_inner and zw_w_inner and output to the 2d/3d NetCDF files

samples added to the user interface which show how to add user-define time series quantities.

calculation/output of precipitation amount, precipitation rate and z0 (by setting "pra*", "prr*", "z0*" with data_output). The time interval on which the precipitation amount is defined is set by new d3par-parameter precipitation_amount_interval

unit 9 opened for debug output (file DEBUG_<pe#>)

Makefile, advec_particles, average_3d_data, buoyancy, calc_precipitation, check_open, check_parameters, data_output_2d, diffusion_e, diffusion_u, diffusion_v, diffusion_w, diffusivities, header, impact_of_latent_heat, init_particles, init_3d_model, modules, netcdf, parin, production_e, read_var_list, read_3d_binary, sum_up_3d_data, user_interface, write_var_list, write_3d_binary

New: wall_fluxes

Changed:

General revision of non-cyclic horizontal boundary conditions: radiation boundary conditions are now used instead of Neumann conditions at the outflow (calculation needs velocity values for t-dt, which are stored on new arrays u_m_l, u_m_r, etc.), calculation of mean outflow is not needed any more, volume flow control is added for the outflow boundary (currently only for the north boundary!!), additional gridpoints along x and y (uxrp, vynp) are not needed any more, routine "boundary_conds" now operates on timelevel t+dt and is not split in two parts (main, uvw_outflow) any more, Neumann boundary conditions at inflow/outflow in case of non-cyclic boundary conditions for all 2d-arrays that are handled by exchange_horiz_2d

The FFT-method for solving the Poisson-equation is now working with Neumann boundary conditions both at the bottom and the top. This requires adjustments of the tridiagonal coefficients and subtracting the horizontally averaged mean from the vertical velocity field.

+age_m in particle_type

Particles-package is now part of the default code ("-p particles" is not needed any more).

Move call of user_actions( 'after_integration' ) below increment of times
and counters. user_actions is now called for each statistic region and has as an argument the number of the respective region (sr)

d3par-parameter data_output_ts removed. Timeseries output for "profil" removed. Timeseries are now switched on by dt_dots. Timeseries data is collected in flow_statistics.

Initial velocities at nzb+1 are regarded for volume flow control in case they have been set zero before (to avoid small timesteps); see new internal parameters u/v_nzb_p1_for_vfc.

q is not allowed to become negative (prognostic_equations).

poisfft_init is only called if fft-solver is switched on (init_pegrid).

d3par-parameter moisture renamed to humidity.

Subversion global revision number is read from mrun and added to the run description header and to the run control (_rc) file.

vtk directives removed from main program.

The uitility routine interpret_config reads PALM environment variables from NAMELIST instead using the system call GETENV.

advec_u_pw, advec_u_up, advec_v_pw, advec_v_up, asselin_filter, check_parameters, coriolis, data_output_dvrp, data_output_ptseries, data_output_ts, data_output_2d, data_output_3d, diffusion_u, diffusion_v, exchange_horiz, exchange_horiz_2d, flow_statistics, header, init_grid, init_particles, init_pegrid, init_rankine, init_pt_anomaly, init_1d_model, init_3d_model, modules, palm, package_parin, parin, poisfft, poismg, prandtl_fluxes, pres, production_e, prognostic_equations, read_var_list, read_3d_binary, sor, swap_timelevel, time_integration, write_var_list, write_3d_binary

Errors:

Bugfix: preset of tendencies te_em, te_um, te_vm in init_1d_model

Bugfix in sample for reading user defined data from restart file (user_init)

Bugfix in setting diffusivities for cases with the outflow damping layer extending over more than one subdomain (init_3d_model)

Check for possible negative humidities in the initial humidity profile.

in Makefile, default suffixes removed from the suffix list to avoid calling of m2c in
# case of .mod files

Makefile
check_parameters, init_1d_model, init_3d_model, user_interface

Property svn:keywords set to Id

File size: 17.1 KB

Line
1	SUBROUTINE boundary_conds( range )
2
3	!------------------------------------------------------------------------------!
4	! Actual revisions:
5	! -----------------
6	!
7	!
8	! Former revisions:
9	! -----------------
10	! $Id: boundary_conds.f90 77 2007-03-29 04:26:56Z raasch $
11	!
12	! 75 2007-03-22 09:54:05Z raasch
13	! The "main" part sets conditions for time level t+dt instead of level t,
14	! outflow boundary conditions changed from Neumann to radiation condition,
15	! uxrp, vynp eliminated, moisture renamed humidity
16	!
17	! 19 2007-02-23 04:53:48Z raasch
18	! Boundary conditions for e(nzt), pt(nzt), and q(nzt) removed because these
19	! gridpoints are now calculated by the prognostic equation,
20	! Dirichlet and zero gradient condition for pt established at top boundary
21	!
22	! RCS Log replace by Id keyword, revision history cleaned up
23	!
24	! Revision 1.15 2006/02/23 09:54:55 raasch
25	! Surface boundary conditions in case of topography: nzb replaced by
26	! 2d-k-index-arrays (nzb_w_inner, etc.). Conditions for u and v remain
27	! unchanged (still using nzb) because a non-flat topography must use a
28	! Prandtl-layer, which don't requires explicit setting of the surface values.
29	!
30	! Revision 1.1 1997/09/12 06:21:34 raasch
31	! Initial revision
32	!
33	!
34	! Description:
35	! ------------
36	! Boundary conditions for the prognostic quantities (range='main').
37	! In case of non-cyclic lateral boundaries the conditions for velocities at
38	! the outflow are set after the pressure solver has been called (range=
39	! 'outflow_uvw').
40	! One additional bottom boundary condition is applied for the TKE (=(u)*2)
41	! in prandtl_fluxes. The cyclic lateral boundary conditions are implicitly
42	! handled in routine exchange_horiz. Pressure boundary conditions are
43	! explicitly set in routines pres, poisfft, poismg and sor.
44	!------------------------------------------------------------------------------!
45
46	USE arrays_3d
47	USE control_parameters
48	USE grid_variables
49	USE indices
50	USE pegrid
51
52	IMPLICIT NONE
53
54	CHARACTER (LEN=*) :: range
55
56	INTEGER :: i, j, k
57
58	REAL :: c_max, c_u, c_v, c_w, denom
59
60
61	IF ( range == 'main') THEN
62	!
63	!-- Bottom boundary
64	IF ( ibc_uv_b == 0 ) THEN
65	!
66	!-- Satisfying the Dirichlet condition with an extra layer below the
67	!-- surface where the u and v component change their sign
68	u_p(nzb,:,:) = -u_p(nzb+1,:,:)
69	v_p(nzb,:,:) = -v_p(nzb+1,:,:)
70	ELSE
71	u_p(nzb,:,:) = u_p(nzb+1,:,:)
72	v_p(nzb,:,:) = v_p(nzb+1,:,:)
73	ENDIF
74	DO i = nxl-1, nxr+1
75	DO j = nys-1, nyn+1
76	w_p(nzb_w_inner(j,i),j,i) = 0.0
77	ENDDO
78	ENDDO
79
80	!
81	!-- Top boundary
82	IF ( ibc_uv_t == 0 ) THEN
83	u_p(nzt+1,:,:) = ug(nzt+1)
84	v_p(nzt+1,:,:) = vg(nzt+1)
85	ELSE
86	u_p(nzt+1,:,:) = u_p(nzt,:,:)
87	v_p(nzt+1,:,:) = v_p(nzt,:,:)
88	ENDIF
89	w_p(nzt:nzt+1,:,:) = 0.0 ! nzt is not a prognostic level (but cf. pres)
90
91	!
92	!-- Temperature at bottom boundary
93	IF ( ibc_pt_b == 0 ) THEN
94	DO i = nxl-1, nxr+1
95	DO j = nys-1, nyn+1
96	pt_p(nzb_s_inner(j,i),j,i) = pt(nzb_s_inner(j,i),j,i)
97	ENDDO
98	ENDDO
99	ELSE
100	DO i = nxl-1, nxr+1
101	DO j = nys-1, nyn+1
102	pt_p(nzb_s_inner(j,i),j,i) = pt_p(nzb_s_inner(j,i)+1,j,i)
103	ENDDO
104	ENDDO
105	ENDIF
106
107	!
108	!-- Temperature at top boundary
109	IF ( ibc_pt_t == 0 ) THEN
110	pt_p(nzt+1,:,:) = pt(nzt+1,:,:)
111	ELSEIF ( ibc_pt_t == 1 ) THEN
112	pt_p(nzt+1,:,:) = pt_p(nzt,:,:)
113	ELSEIF ( ibc_pt_t == 2 ) THEN
114	pt_p(nzt+1,:,:) = pt_p(nzt,:,:) + bc_pt_t_val * dzu(nzt+1)
115	ENDIF
116
117	!
118	!-- Boundary conditions for TKE
119	!-- Generally Neumann conditions with de/dz=0 are assumed
120	IF ( .NOT. constant_diffusion ) THEN
121	DO i = nxl-1, nxr+1
122	DO j = nys-1, nyn+1
123	e_p(nzb_s_inner(j,i),j,i) = e_p(nzb_s_inner(j,i)+1,j,i)
124	ENDDO
125	ENDDO
126	e_p(nzt+1,:,:) = e_p(nzt,:,:)
127	ENDIF
128
129	!
130	!-- Boundary conditions for total water content or scalar,
131	!-- bottom and surface boundary (see also temperature)
132	IF ( humidity .OR. passive_scalar ) THEN
133	!
134	!-- Surface conditions for constant_humidity_flux
135	IF ( ibc_q_b == 0 ) THEN
136	DO i = nxl-1, nxr+1
137	DO j = nys-1, nyn+1
138	q_p(nzb_s_inner(j,i),j,i) = q(nzb_s_inner(j,i),j,i)
139	ENDDO
140	ENDDO
141	ELSE
142	DO i = nxl-1, nxr+1
143	DO j = nys-1, nyn+1
144	q_p(nzb_s_inner(j,i),j,i) = q_p(nzb_s_inner(j,i)+1,j,i)
145	ENDDO
146	ENDDO
147	ENDIF
148	!
149	!-- Top boundary
150	q_p(nzt+1,:,:) = q_p(nzt,:,:) + bc_q_t_val * dzu(nzt+1)
151	ENDIF
152
153	!
154	!-- Lateral boundary conditions at the inflow. Quasi Neumann conditions
155	!-- are needed for the wall normal velocity in order to ensure zero
156	!-- divergence. Dirichlet conditions are used for all other quantities.
157	IF ( inflow_s ) THEN
158	v_p(:,nys,:) = v_p(:,nys-1,:)
159	ELSEIF ( inflow_n ) THEN
160	v_p(:,nyn,:) = v_p(:,nyn+1,:)
161	ELSEIF ( inflow_l ) THEN
162	u_p(:,:,nxl) = u_p(:,:,nxl-1)
163	ELSEIF ( inflow_r ) THEN
164	u_p(:,:,nxr) = u_p(:,:,nxr+1)
165	ENDIF
166
167	!
168	!-- Lateral boundary conditions for scalar quantities at the outflow
169	IF ( outflow_s ) THEN
170	pt_p(:,nys-1,:) = pt_p(:,nys,:)
171	IF ( .NOT. constant_diffusion ) e_p(:,nys-1,:) = e_p(:,nys,:)
172	IF ( humidity .OR. passive_scalar ) q_p(:,nys-1,:) = q_p(:,nys,:)
173	ELSEIF ( outflow_n ) THEN
174	pt_p(:,nyn+1,:) = pt_p(:,nyn,:)
175	IF ( .NOT. constant_diffusion ) e_p(:,nyn+1,:) = e_p(:,nyn,:)
176	IF ( humidity .OR. passive_scalar ) q_p(:,nyn+1,:) = q_p(:,nyn,:)
177	ELSEIF ( outflow_l ) THEN
178	pt_p(:,:,nxl-1) = pt_p(:,:,nxl)
179	IF ( .NOT. constant_diffusion ) e_p(:,:,nxl-1) = e_p(:,:,nxl)
180	IF ( humidity .OR. passive_scalar ) q_p(:,:,nxl-1) = q_p(:,:,nxl)
181	ELSEIF ( outflow_r ) THEN
182	pt_p(:,:,nxr+1) = pt_p(:,:,nxr)
183	IF ( .NOT. constant_diffusion ) e_p(:,:,nxr+1) = e_p(:,:,nxr)
184	IF ( humidity .OR. passive_scalar ) q_p(:,:,nxr+1) = q_p(:,:,nxr)
185	ENDIF
186
187	ENDIF
188
189	!
190	!-- Radiation boundary condition for the velocities at the respective outflow
191	IF ( outflow_s .AND. &
192	intermediate_timestep_count == intermediate_timestep_count_max ) &
193	THEN
194
195	c_max = dy / dt_3d
196
197	DO i = nxl-1, nxr+1
198	DO k = nzb+1, nzt+1
199
200	!
201	!-- First calculate the phase speeds for u,v, and w
202	denom = u_m_s(k,0,i) - u_m_s(k,1,i)
203
204	IF ( denom /= 0.0 ) THEN
205	c_u = -c_max * ( u(k,0,i) - u_m_s(k,0,i) ) / denom
206	IF ( c_u < 0.0 ) THEN
207	c_u = 0.0
208	ELSEIF ( c_u > c_max ) THEN
209	c_u = c_max
210	ENDIF
211	ELSE
212	c_u = c_max
213	ENDIF
214	denom = v_m_s(k,0,i) - v_m_s(k,1,i)
215
216	IF ( denom /= 0.0 ) THEN
217	c_v = -c_max * ( v(k,0,i) - v_m_s(k,0,i) ) / denom
218	IF ( c_v < 0.0 ) THEN
219	c_v = 0.0
220	ELSEIF ( c_v > c_max ) THEN
221	c_v = c_max
222	ENDIF
223	ELSE
224	c_v = c_max
225	ENDIF
226
227	denom = w_m_s(k,0,i) - w_m_s(k,1,i)
228
229	IF ( denom /= 0.0 ) THEN
230	c_w = -c_max * ( w(k,0,i) - w_m_s(k,0,i) ) / denom
231	IF ( c_w < 0.0 ) THEN
232	c_w = 0.0
233	ELSEIF ( c_w > c_max ) THEN
234	c_w = c_max
235	ENDIF
236	ELSE
237	c_w = c_max
238	ENDIF
239
240	!
241	!-- Calculate the new velocities
242	u_p(k,-1,i) = u(k,-1,i) + dt_3d * c_u * &
243	( u(k,-1,i) - u(k,0,i) ) * ddy
244
245	v_p(k,-1,i) = v(k,-1,i) + dt_3d * c_v * &
246	( v(k,-1,i) - v_m_s(k,0,i) ) * ddy
247
248	w_p(k,-1,i) = w(k,-1,i) + dt_3d * c_w * &
249	( w(k,-1,i) - w(k,0,i) ) * ddy
250
251	!
252	!-- Save old timelevels for the next timestep
253	u_m_s(k,:,i) = u(k,-1:1,i)
254	v_m_s(k,:,i) = v(k,-1:1,i)
255	w_m_s(k,:,i) = w(k,-1:1,i)
256
257	ENDDO
258	ENDDO
259
260	!
261	!-- Bottom boundary at the outflow
262	IF ( ibc_uv_b == 0 ) THEN
263	u_p(nzb,-1,:) = -u_p(nzb+1,-1,:)
264	v_p(nzb,-1,:) = -v_p(nzb+1,-1,:)
265	ELSE
266	u_p(nzb,-1,:) = u_p(nzb+1,-1,:)
267	v_p(nzb,-1,:) = v_p(nzb+1,-1,:)
268	ENDIF
269	w_p(nzb,ny+1,:) = 0.0
270
271	!
272	!-- Top boundary at the outflow
273	IF ( ibc_uv_t == 0 ) THEN
274	u_p(nzt+1,-1,:) = ug(nzt+1)
275	v_p(nzt+1,-1,:) = vg(nzt+1)
276	ELSE
277	u_p(nzt+1,-1,:) = u(nzt,-1,:)
278	v_p(nzt+1,-1,:) = v(nzt,-1,:)
279	ENDIF
280	w_p(nzt:nzt+1,-1,:) = 0.0
281
282	ENDIF
283
284	IF ( outflow_n .AND. &
285	intermediate_timestep_count == intermediate_timestep_count_max ) &
286	THEN
287
288	c_max = dy / dt_3d
289
290	DO i = nxl-1, nxr+1
291	DO k = nzb+1, nzt+1
292
293	!
294	!-- First calculate the phase speeds for u,v, and w
295	denom = u_m_n(k,ny,i) - u_m_n(k,ny-1,i)
296
297	IF ( denom /= 0.0 ) THEN
298	c_u = -c_max * ( u(k,ny,i) - u_m_n(k,ny,i) ) / denom
299	IF ( c_u < 0.0 ) THEN
300	c_u = 0.0
301	ELSEIF ( c_u > c_max ) THEN
302	c_u = c_max
303	ENDIF
304	ELSE
305	c_u = c_max
306	ENDIF
307
308	denom = v_m_n(k,ny,i) - v_m_n(k,ny-1,i)
309
310	IF ( denom /= 0.0 ) THEN
311	c_v = -c_max * ( v(k,ny,i) - v_m_n(k,ny,i) ) / denom
312	IF ( c_v < 0.0 ) THEN
313	c_v = 0.0
314	ELSEIF ( c_v > c_max ) THEN
315	c_v = c_max
316	ENDIF
317	ELSE
318	c_v = c_max
319	ENDIF
320
321	denom = w_m_n(k,ny,i) - w_m_n(k,ny-1,i)
322
323	IF ( denom /= 0.0 ) THEN
324	c_w = -c_max * ( w(k,ny,i) - w_m_n(k,ny,i) ) / denom
325	IF ( c_w < 0.0 ) THEN
326	c_w = 0.0
327	ELSEIF ( c_w > c_max ) THEN
328	c_w = c_max
329	ENDIF
330	ELSE
331	c_w = c_max
332	ENDIF
333
334	!
335	!-- Calculate the new velocities
336	u_p(k,ny+1,i) = u(k,ny+1,i) - dt_3d * c_u * &
337	( u(k,ny+1,i) - u(k,ny,i) ) * ddy
338
339	v_p(k,ny+1,i) = v(k,ny+1,i) - dt_3d * c_v * &
340	( v(k,ny+1,i) - v(k,ny,i) ) * ddy
341
342	w_p(k,ny+1,i) = w(k,ny+1,i) - dt_3d * c_w * &
343	( w(k,ny+1,i) - w(k,ny,i) ) * ddy
344
345	!
346	!-- Swap timelevels for the next timestep
347	u_m_n(k,:,i) = u(k,ny-1:ny+1,i)
348	v_m_n(k,:,i) = v(k,ny-1:ny+1,i)
349	w_m_n(k,:,i) = w(k,ny-1:ny+1,i)
350
351	ENDDO
352	ENDDO
353
354	!
355	!-- Bottom boundary at the outflow
356	IF ( ibc_uv_b == 0 ) THEN
357	u_p(nzb,ny+1,:) = -u_p(nzb+1,ny+1,:)
358	v_p(nzb,ny+1,:) = -v_p(nzb+1,ny+1,:)
359	ELSE
360	u_p(nzb,ny+1,:) = u_p(nzb+1,ny+1,:)
361	v_p(nzb,ny+1,:) = v_p(nzb+1,ny+1,:)
362	ENDIF
363	w_p(nzb,ny+1,:) = 0.0
364
365	!
366	!-- Top boundary at the outflow
367	IF ( ibc_uv_t == 0 ) THEN
368	u_p(nzt+1,ny+1,:) = ug(nzt+1)
369	v_p(nzt+1,ny+1,:) = vg(nzt+1)
370	ELSE
371	u_p(nzt+1,ny+1,:) = u_p(nzt,nyn+1,:)
372	v_p(nzt+1,ny+1,:) = v_p(nzt,nyn+1,:)
373	ENDIF
374	w_p(nzt:nzt+1,ny+1,:) = 0.0
375
376	ENDIF
377
378	IF ( outflow_l .AND. &
379	intermediate_timestep_count == intermediate_timestep_count_max ) &
380	THEN
381
382	c_max = dx / dt_3d
383
384	DO j = nys-1, nyn+1
385	DO k = nzb+1, nzt+1
386
387	!
388	!-- First calculate the phase speeds for u,v, and w
389	denom = u_m_l(k,j,0) - u_m_l(k,j,1)
390
391	IF ( denom /= 0.0 ) THEN
392	c_u = -c_max * ( u(k,j,0) - u_m_r(k,j,0) ) / denom
393	IF ( c_u > 0.0 ) THEN
394	c_u = 0.0
395	ELSEIF ( c_u < -c_max ) THEN
396	c_u = -c_max
397	ENDIF
398	ELSE
399	c_u = -c_max
400	ENDIF
401
402	denom = v_m_l(k,j,0) - v_m_l(k,j,1)
403
404	IF ( denom /= 0.0 ) THEN
405	c_v = -c_max * ( v(k,j,0) - v_m_l(k,j,0) ) / denom
406	IF ( c_v < 0.0 ) THEN
407	c_v = 0.0
408	ELSEIF ( c_v > c_max ) THEN
409	c_v = c_max
410	ENDIF
411	ELSE
412	c_v = c_max
413	ENDIF
414
415	denom = w_m_l(k,j,0) - w_m_l(k,j,1)
416
417	IF ( denom /= 0.0 ) THEN
418	c_w = -c_max * ( w(k,j,0) - w_m_l(k,j,0) ) / denom
419	IF ( c_w < 0.0 ) THEN
420	c_w = 0.0
421	ELSEIF ( c_w > c_max ) THEN
422	c_w = c_max
423	ENDIF
424	ELSE
425	c_w = c_max
426	ENDIF
427
428	!
429	!-- Calculate the new velocities
430	u_p(k,j,-1) = u(k,j,-1) + dt_3d * c_u * &
431	( u(k,j,-1) - u(k,j,0) ) * ddx
432
433	v_p(k,j,-1) = v(k,j,-1) + dt_3d * c_v * &
434	( v(k,j,-1) - v(k,j,0) ) * ddx
435
436	w_p(k,j,-1) = w(k,j,-1) + dt_3d * c_w * &
437	( w(k,j,-1) - w(k,j,0) ) * ddx
438
439	!
440	!-- Swap timelevels for the next timestep
441	u_m_l(k,j,:) = u(k,j,-1:1)
442	v_m_l(k,j,:) = v(k,j,-1:1)
443	w_m_l(k,j,:) = w(k,j,-1:1)
444
445	ENDDO
446	ENDDO
447
448	!
449	!-- Bottom boundary at the outflow
450	IF ( ibc_uv_b == 0 ) THEN
451	u_p(nzb,:,-1) = -u_p(nzb+1,:,-1)
452	v_p(nzb,:,-1) = -v_p(nzb+1,:,-1)
453	ELSE
454	u_p(nzb,:,-1) = u_p(nzb+1,:,-1)
455	v_p(nzb,:,-1) = v_p(nzb+1,:,-1)
456	ENDIF
457	w_p(nzb,:,-1) = 0.0
458
459	!
460	!-- Top boundary at the outflow
461	IF ( ibc_uv_t == 0 ) THEN
462	u_p(nzt+1,:,-1) = ug(nzt+1)
463	v_p(nzt+1,:,-1) = vg(nzt+1)
464	ELSE
465	u_p(nzt+1,:,-1) = u_p(nzt,:,-1)
466	v_p(nzt+1,:,-1) = v_p(nzt,:,-1)
467	ENDIF
468	w_p(nzt:nzt+1,:,-1) = 0.0
469
470	ENDIF
471
472	IF ( outflow_r .AND. &
473	intermediate_timestep_count == intermediate_timestep_count_max ) &
474	THEN
475
476	c_max = dx / dt_3d
477
478	DO j = nys-1, nyn+1
479	DO k = nzb+1, nzt+1
480
481	!
482	!-- First calculate the phase speeds for u,v, and w
483	denom = u_m_r(k,j,nx) - u_m_r(k,j,nx-1)
484
485	IF ( denom /= 0.0 ) THEN
486	c_u = -c_max * ( u(k,j,nx) - u_m_r(k,j,nx) ) / denom
487	IF ( c_u < 0.0 ) THEN
488	c_u = 0.0
489	ELSEIF ( c_u > c_max ) THEN
490	c_u = c_max
491	ENDIF
492	ELSE
493	c_u = c_max
494	ENDIF
495
496	denom = v_m_r(k,j,nx) - v_m_r(k,j,nx-1)
497
498	IF ( denom /= 0.0 ) THEN
499	c_v = -c_max * ( v(k,j,nx) - v_m_r(k,j,nx) ) / denom
500	IF ( c_v < 0.0 ) THEN
501	c_v = 0.0
502	ELSEIF ( c_v > c_max ) THEN
503	c_v = c_max
504	ENDIF
505	ELSE
506	c_v = c_max
507	ENDIF
508
509	denom = w_m_r(k,j,nx) - w_m_r(k,j,nx-1)
510
511	IF ( denom /= 0.0 ) THEN
512	c_w = -c_max * ( w(k,j,nx) - w_m_r(k,j,nx) ) / denom
513	IF ( c_w < 0.0 ) THEN
514	c_w = 0.0
515	ELSEIF ( c_w > c_max ) THEN
516	c_w = c_max
517	ENDIF
518	ELSE
519	c_w = c_max
520	ENDIF
521
522	!
523	!-- Calculate the new velocities
524	u_p(k,j,nx+1) = u(k,j,nx+1) - dt_3d * c_u * &
525	( u(k,j,nx+1) - u(k,j,nx) ) * ddx
526
527	v_p(k,j,nx+1) = v(k,j,nx+1) - dt_3d * c_v * &
528	( v(k,j,nx+1) - v(k,j,nx) ) * ddx
529
530	w_p(k,j,nx+1) = w(k,j,nx+1) - dt_3d * c_w * &
531	( w(k,j,nx+1) - w(k,j,nx) ) * ddx
532
533	!
534	!-- Swap timelevels for the next timestep
535	u_m_r(k,j,:) = u(k,j,nx-1:nx+1)
536	v_m_r(k,j,:) = v(k,j,nx-1:nx+1)
537	w_m_r(k,j,:) = w(k,j,nx-1:nx+1)
538
539	ENDDO
540	ENDDO
541
542	!
543	!-- Bottom boundary at the outflow
544	IF ( ibc_uv_b == 0 ) THEN
545	u_p(nzb,:,nx+1) = -u_p(nzb+1,:,nx+1)
546	v_p(nzb,:,nx+1) = -v_p(nzb+1,:,nx+1)
547	ELSE
548	u_p(nzb,:,nx+1) = u_p(nzb+1,:,nx+1)
549	v_p(nzb,:,nx+1) = v_p(nzb+1,:,nx+1)
550	ENDIF
551	w_p(nzb,:,nx+1) = 0.0
552
553	!
554	!-- Top boundary at the outflow
555	IF ( ibc_uv_t == 0 ) THEN
556	u_p(nzt+1,:,nx+1) = ug(nzt+1)
557	v_p(nzt+1,:,nx+1) = vg(nzt+1)
558	ELSE
559	u_p(nzt+1,:,nx+1) = u_p(nzt,:,nx+1)
560	v_p(nzt+1,:,nx+1) = v_p(nzt,:,nx+1)
561	ENDIF
562	w(nzt:nzt+1,:,nx+1) = 0.0
563
564	ENDIF
565
566
567	END SUBROUTINE boundary_conds

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

Context Navigation

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

Download in other formats: