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

Last change on this file since 2938 was 2938, checked in by suehring, 3 years ago

Nesting in RANS-LES and RANS-RANS mode enabled; synthetic turbulence generator at all lateral boundaries in nesting or non-cyclic forcing mode; revised Inifor initialization in nesting mode

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