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

Last change on this file since 2320 was 2320, checked in by suehring, 4 years ago

large-scale forcing and nudging modularized

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