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

Last change on this file since 2556 was 2365, checked in by kanani, 7 years ago

Vertical nesting implemented (SadiqHuq?)

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