source: palm/tags/release-3.2/SOURCE/boundary_conds.f90 @ 235

Last change on this file since 235 was 77, checked in by raasch, 17 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.