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

Last change on this file since 2292 was 2292, checked in by schwenkel, 4 years ago

implementation of new bulk microphysics scheme

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