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

Last change on this file since 342 was 110, checked in by raasch, 17 years ago

New:
---
Allows runs for a coupled atmosphere-ocean LES,
coupling frequency is controlled by new d3par-parameter dt_coupling,
the coupling mode (atmosphere_to_ocean or ocean_to_atmosphere) for the
respective processes is read from environment variable coupling_mode,
which is set by the mpiexec-command,
communication between the two models is done using the intercommunicator
comm_inter,
local files opened by the ocean model get the additional suffic "_O".
Assume saturation at k=nzb_s_inner(j,i) for atmosphere coupled to ocean.

A momentum flux can be set as top boundary condition using the new
inipar parameter top_momentumflux_u|v.

Non-cyclic boundary conditions can be used along all horizontal directions.

Quantities w*p* and w"e can be output as vertical profiles.

Initial profiles are reset to constant profiles in case that initializing_actions /= 'set_constant_profiles'. (init_rankine)

Optionally calculate km and kh from initial TKE e_init.

Changed:


Remaining variables iran changed to iran_part (advec_particles, init_particles).

In case that the presure solver is not called for every Runge-Kutta substep
(call_psolver_at_all_substeps = .F.), it is called after the first substep
instead of the last. In that case, random perturbations are also added to the
velocity field after the first substep.

Initialization of km,kh = 0.00001 for ocean = .T. (for ocean = .F. it remains 0.01).

Allow data_output_pr= q, wq, w"q", w*q* for humidity = .T. (instead of cloud_physics = .T.).

Errors:


Bugs from code parts for non-cyclic boundary conditions are removed: loops for
u and v are starting from index nxlu, nysv, respectively. The radiation boundary
condition is used for every Runge-Kutta substep. Velocity phase speeds for
the radiation boundary conditions are calculated for the first Runge-Kutta
substep only and reused for the further substeps. New arrays c_u, c_v, and c_w
are defined for this purpose. Several index errors are removed from the
radiation boundary condition code parts. Upper bounds for calculating
u_0 and v_0 (in production_e) are nxr+1 and nyn+1 because otherwise these
values are not available in case of non-cyclic boundary conditions.

+dots_num_palm in module user, +module netcdf_control in user_init (both in user_interface)

Bugfix: wrong sign removed from the buoyancy production term in the case use_reference = .T. (production_e)

Bugfix: Error message concerning output of particle concentration (pc) modified (check_parameters).

Bugfix: Rayleigh damping for ocean fixed.

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