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

Last change on this file since 2236 was 2233, checked in by suehring, 7 years ago

last commit documented

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