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

Last change on this file since 2766 was 2766, checked in by kanani, 4 years ago

Removal of chem directive, plus minor changes

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