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

Last change on this file since 2709 was 2696, checked in by kanani, 7 years ago

Merge of branch palm4u into trunk

  • Property svn:keywords set to Id
File size: 40.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-2017 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22!
23!
24! Former revisions:
25! -----------------
26! $Id: boundary_conds.f90 2696 2017-12-14 17:12:51Z suehring $
27! Adjust boundary conditions for e and diss in case of TKE-e closure (TG)
28! Implementation of chemistry module (FK)
29!
30! 2569 2017-10-20 11:54:42Z kanani
31! Removed redundant code for ibc_s_b=1 and ibc_q_b=1
32!
33! 2365 2017-08-21 14:59:59Z kanani
34! Vertical grid nesting implemented: exclude setting vertical velocity to zero
35! on fine grid (SadiqHuq)
36!
37! 2320 2017-07-21 12:47:43Z suehring
38! Remove unused control parameter large_scale_forcing from only-list
39!
40! 2292 2017-06-20 09:51:42Z schwenkel
41! Implementation of new microphysic scheme: cloud_scheme = 'morrison'
42! includes two more prognostic equations for cloud drop concentration (nc) 
43! and cloud water content (qc).
44!
45! 2233 2017-05-30 18:08:54Z suehring
46!
47! 2232 2017-05-30 17:47:52Z suehring
48! Set boundary conditions on topography top using flag method.
49!
50! 2118 2017-01-17 16:38:49Z raasch
51! OpenACC directives removed
52!
53! 2000 2016-08-20 18:09:15Z knoop
54! Forced header and separation lines into 80 columns
55!
56! 1992 2016-08-12 15:14:59Z suehring
57! Adjustments for top boundary condition for passive scalar
58!
59! 1960 2016-07-12 16:34:24Z suehring
60! Treat humidity and passive scalar separately
61!
62! 1823 2016-04-07 08:57:52Z hoffmann
63! Initial version of purely vertical nesting introduced.
64!
65! 1822 2016-04-07 07:49:42Z hoffmann
66! icloud_scheme removed. microphyisics_seifert added.
67!
68! 1764 2016-02-28 12:45:19Z raasch
69! index bug for u_p at left outflow removed
70!
71! 1762 2016-02-25 12:31:13Z hellstea
72! Introduction of nested domain feature
73!
74! 1742 2016-01-13 09:50:06Z raasch
75! bugfix for outflow Neumann boundary conditions at bottom and top
76!
77! 1717 2015-11-11 15:09:47Z raasch
78! Bugfix: index error in outflow conditions for left boundary
79!
80! 1682 2015-10-07 23:56:08Z knoop
81! Code annotations made doxygen readable
82!
83! 1410 2014-05-23 12:16:18Z suehring
84! Bugfix: set dirichlet boundary condition for passive_scalar at model domain
85! top
86!
87! 1399 2014-05-07 11:16:25Z heinze
88! Bugfix: set inflow boundary conditions also if no humidity or passive_scalar
89! is used.
90!
91! 1398 2014-05-07 11:15:00Z heinze
92! Dirichlet-condition at the top for u and v changed to u_init and v_init also
93! for large_scale_forcing
94!
95! 1380 2014-04-28 12:40:45Z heinze
96! Adjust Dirichlet-condition at the top for pt in case of nudging
97!
98! 1361 2014-04-16 15:17:48Z hoffmann
99! Bottom and top boundary conditions of rain water content (qr) and
100! rain drop concentration (nr) changed to Dirichlet
101!
102! 1353 2014-04-08 15:21:23Z heinze
103! REAL constants provided with KIND-attribute
104
105! 1320 2014-03-20 08:40:49Z raasch
106! ONLY-attribute added to USE-statements,
107! kind-parameters added to all INTEGER and REAL declaration statements,
108! kinds are defined in new module kinds,
109! revision history before 2012 removed,
110! comment fields (!:) to be used for variable explanations added to
111! all variable declaration statements
112!
113! 1257 2013-11-08 15:18:40Z raasch
114! loop independent clauses added
115!
116! 1241 2013-10-30 11:36:58Z heinze
117! Adjust ug and vg at each timestep in case of large_scale_forcing
118!
119! 1159 2013-05-21 11:58:22Z fricke
120! Bugfix: Neumann boundary conditions for the velocity components at the
121! outflow are in fact radiation boundary conditions using the maximum phase
122! velocity that ensures numerical stability (CFL-condition).
123! Hence, logical operator use_cmax is now used instead of bc_lr_dirneu/_neudir.
124! Bugfix: In case of use_cmax at the outflow, u, v, w are replaced by
125! u_p, v_p, w_p 
126!
127! 1115 2013-03-26 18:16:16Z hoffmann
128! boundary conditions of two-moment cloud scheme are restricted to Neumann-
129! boundary-conditions
130!
131! 1113 2013-03-10 02:48:14Z raasch
132! GPU-porting
133! dummy argument "range" removed
134! Bugfix: wrong index in loops of radiation boundary condition
135!
136! 1053 2012-11-13 17:11:03Z hoffmann
137! boundary conditions for the two new prognostic equations (nr, qr) of the
138! two-moment cloud scheme
139!
140! 1036 2012-10-22 13:43:42Z raasch
141! code put under GPL (PALM 3.9)
142!
143! 996 2012-09-07 10:41:47Z raasch
144! little reformatting
145!
146! 978 2012-08-09 08:28:32Z fricke
147! Neumann boudnary conditions are added at the inflow boundary for the SGS-TKE.
148! Outflow boundary conditions for the velocity components can be set to Neumann
149! conditions or to radiation conditions with a horizontal averaged phase
150! velocity.
151!
152! 875 2012-04-02 15:35:15Z gryschka
153! Bugfix in case of dirichlet inflow bc at the right or north boundary
154!
155! Revision 1.1  1997/09/12 06:21:34  raasch
156! Initial revision
157!
158!
159! Description:
160! ------------
161!> Boundary conditions for the prognostic quantities.
162!> One additional bottom boundary condition is applied for the TKE (=(u*)**2)
163!> in prandtl_fluxes. The cyclic lateral boundary conditions are implicitly
164!> handled in routine exchange_horiz. Pressure boundary conditions are
165!> explicitly set in routines pres, poisfft, poismg and sor.
166!------------------------------------------------------------------------------!
167 SUBROUTINE boundary_conds
168 
169
170    USE arrays_3d,                                                             &
171        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,  &
172               diss_p, dzu, e_p, nc_p, nr_p, pt, pt_p, q, q_p, qc_p, qr_p, s,  & 
173               s_p, sa, sa_p, u, ug, u_init, u_m_l, u_m_n, u_m_r, u_m_s, u_p,  &
174               v, vg, v_init, v_m_l, v_m_n, v_m_r, v_m_s, v_p,                 &
175               w, w_p, w_m_l, w_m_n, w_m_r, w_m_s, pt_init
176
177#if defined( __chem )
178    USE chemistry_model_mod,                                                   &
179        ONLY:  chem_boundary_conds 
180#endif
181             
182    USE control_parameters,                                                    &
183        ONLY:  air_chemistry, bc_pt_t_val, bc_q_t_val, bc_s_t_val,             &
184               constant_diffusion, cloud_physics, coupling_mode, dt_3d,        &
185               force_bound_l, force_bound_s, forcing, humidity,                &
186               ibc_pt_b, ibc_pt_t, ibc_q_b, ibc_q_t, ibc_s_b, ibc_s_t,         &
187               ibc_sa_t, ibc_uv_b, ibc_uv_t, inflow_l, inflow_n, inflow_r,     &
188               inflow_s, intermediate_timestep_count,                          &
189               microphysics_morrison, microphysics_seifert, nest_domain,       &
190               nest_bound_l, nest_bound_s, nudging, ocean, outflow_l,          &
191               outflow_n, outflow_r, outflow_s, passive_scalar, rans_tke_e,    &
192               tsc, use_cmax
193
194    USE grid_variables,                                                        &
195        ONLY:  ddx, ddy, dx, dy
196
197    USE indices,                                                               &
198        ONLY:  nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg,             &
199               nzb, nzt, wall_flags_0
200
201    USE kinds
202
203    USE pegrid
204
205    USE pmc_interface,                                                         &
206        ONLY : nesting_mode
207
208    USE surface_mod,                                                           &
209        ONLY :  bc_h
210
211    IMPLICIT NONE
212
213    INTEGER(iwp) ::  i  !< grid index x direction
214    INTEGER(iwp) ::  j  !< grid index y direction
215    INTEGER(iwp) ::  k  !< grid index z direction
216    INTEGER(iwp) ::  kb !< variable to set respective boundary value, depends on facing.
217    INTEGER(iwp) ::  l  !< running index boundary type, for up- and downward-facing walls
218    INTEGER(iwp) ::  m  !< running index surface elements
219
220    REAL(wp)    ::  c_max !<
221    REAL(wp)    ::  denom !<
222
223
224!
225!-- Bottom boundary
226    IF ( ibc_uv_b == 1 )  THEN
227       u_p(nzb,:,:) = u_p(nzb+1,:,:)
228       v_p(nzb,:,:) = v_p(nzb+1,:,:)
229    ENDIF
230!
231!-- Set zero vertical velocity at topography top (l=0), or bottom (l=1) in case
232!-- of downward-facing surfaces.
233    DO  l = 0, 1
234!
235!--    Set kb, for upward-facing surfaces value at topography top (k-1) is set,
236!--    for downward-facing surfaces at topography bottom (k+1).
237       kb = MERGE( -1, 1, l == 0 )
238       !$OMP PARALLEL DO PRIVATE( i, j, k )
239       DO  m = 1, bc_h(l)%ns
240          i = bc_h(l)%i(m)           
241          j = bc_h(l)%j(m)
242          k = bc_h(l)%k(m)
243          w_p(k+kb,j,i) = 0.0_wp
244       ENDDO
245    ENDDO
246
247!
248!-- Top boundary. A nested domain ( ibc_uv_t = 3 ) does not require settings.
249    IF ( ibc_uv_t == 0 )  THEN
250        u_p(nzt+1,:,:) = u_init(nzt+1)
251        v_p(nzt+1,:,:) = v_init(nzt+1)
252    ELSEIF ( ibc_uv_t == 1 )  THEN
253        u_p(nzt+1,:,:) = u_p(nzt,:,:)
254        v_p(nzt+1,:,:) = v_p(nzt,:,:)
255    ENDIF
256
257!
258!-- Vertical nesting: Vertical velocity not zero at the top of the fine grid
259    IF (  .NOT.  nest_domain  .AND.                                            &
260                 TRIM(coupling_mode) /= 'vnested_fine' )  THEN
261       w_p(nzt:nzt+1,:,:) = 0.0_wp  !< nzt is not a prognostic level (but cf. pres)
262    ENDIF
263
264!
265!-- Temperature at bottom and top boundary.
266!-- In case of coupled runs (ibc_pt_b = 2) the temperature is given by
267!-- the sea surface temperature of the coupled ocean model.
268!-- Dirichlet
269    IF ( ibc_pt_b == 0 )  THEN
270       DO  l = 0, 1
271!
272!--       Set kb, for upward-facing surfaces value at topography top (k-1) is set,
273!--       for downward-facing surfaces at topography bottom (k+1).
274          kb = MERGE( -1, 1, l == 0 )
275          !$OMP PARALLEL DO PRIVATE( i, j, k )
276          DO  m = 1, bc_h(l)%ns
277             i = bc_h(l)%i(m)           
278             j = bc_h(l)%j(m)
279             k = bc_h(l)%k(m)
280             pt_p(k+kb,j,i) = pt(k+kb,j,i)
281          ENDDO
282       ENDDO
283!
284!-- Neumann, zero-gradient
285    ELSEIF ( ibc_pt_b == 1 )  THEN
286       DO  l = 0, 1
287!
288!--       Set kb, for upward-facing surfaces value at topography top (k-1) is set,
289!--       for downward-facing surfaces at topography bottom (k+1).
290          kb = MERGE( -1, 1, l == 0 )
291          !$OMP PARALLEL DO PRIVATE( i, j, k )
292          DO  m = 1, bc_h(l)%ns
293             i = bc_h(l)%i(m)           
294             j = bc_h(l)%j(m)
295             k = bc_h(l)%k(m)
296             pt_p(k+kb,j,i) = pt_p(k,j,i)
297          ENDDO
298       ENDDO
299    ENDIF
300
301!
302!-- Temperature at top boundary
303    IF ( ibc_pt_t == 0 )  THEN
304        pt_p(nzt+1,:,:) = pt(nzt+1,:,:)
305!
306!--     In case of nudging adjust top boundary to pt which is
307!--     read in from NUDGING-DATA
308        IF ( nudging )  THEN
309           pt_p(nzt+1,:,:) = pt_init(nzt+1)
310        ENDIF
311    ELSEIF ( ibc_pt_t == 1 )  THEN
312        pt_p(nzt+1,:,:) = pt_p(nzt,:,:)
313    ELSEIF ( ibc_pt_t == 2 )  THEN
314        pt_p(nzt+1,:,:) = pt_p(nzt,:,:) + bc_pt_t_val * dzu(nzt+1)
315    ENDIF
316
317!
318!-- Boundary conditions for TKE
319!-- Generally Neumann conditions with de/dz=0 are assumed
320    IF ( .NOT. constant_diffusion )  THEN
321
322       IF ( .NOT. rans_tke_e )  THEN
323          DO  l = 0, 1
324!
325!--         Set kb, for upward-facing surfaces value at topography top (k-1) is set,
326!--         for downward-facing surfaces at topography bottom (k+1).
327             kb = MERGE( -1, 1, l == 0 )
328             !$OMP PARALLEL DO PRIVATE( i, j, k )
329             DO  m = 1, bc_h(l)%ns
330                i = bc_h(l)%i(m)           
331                j = bc_h(l)%j(m)
332                k = bc_h(l)%k(m)
333                e_p(k+kb,j,i) = e_p(k,j,i)
334             ENDDO
335          ENDDO
336       ENDIF
337
338       IF ( .NOT. nest_domain )  THEN
339          e_p(nzt+1,:,:) = e_p(nzt,:,:)
340       ENDIF
341    ENDIF
342
343!
344!-- Boundary conditions for TKE dissipation rate
345    IF ( rans_tke_e .AND. .NOT. nest_domain )  THEN
346       diss_p(nzt+1,:,:) = diss_p(nzt,:,:)
347    ENDIF
348
349!
350!-- Boundary conditions for salinity
351    IF ( ocean )  THEN
352!
353!--    Bottom boundary: Neumann condition because salinity flux is always
354!--    given.
355       DO  l = 0, 1
356!
357!--       Set kb, for upward-facing surfaces value at topography top (k-1) is set,
358!--       for downward-facing surfaces at topography bottom (k+1).
359          kb = MERGE( -1, 1, l == 0 )
360          !$OMP PARALLEL DO PRIVATE( i, j, k )
361          DO  m = 1, bc_h(l)%ns
362             i = bc_h(l)%i(m)           
363             j = bc_h(l)%j(m)
364             k = bc_h(l)%k(m)
365             sa_p(k+kb,j,i) = sa_p(k,j,i)
366          ENDDO
367       ENDDO
368!
369!--    Top boundary: Dirichlet or Neumann
370       IF ( ibc_sa_t == 0 )  THEN
371           sa_p(nzt+1,:,:) = sa(nzt+1,:,:)
372       ELSEIF ( ibc_sa_t == 1 )  THEN
373           sa_p(nzt+1,:,:) = sa_p(nzt,:,:)
374       ENDIF
375
376    ENDIF
377
378!
379!-- Boundary conditions for total water content,
380!-- bottom and top boundary (see also temperature)
381    IF ( humidity )  THEN
382!
383!--    Surface conditions for constant_humidity_flux
384!--    Run loop over all non-natural and natural walls. Note, in wall-datatype
385!--    the k coordinate belongs to the atmospheric grid point, therefore, set
386!--    q_p at k-1
387       IF ( ibc_q_b == 0 ) THEN
388
389          DO  l = 0, 1
390!
391!--          Set kb, for upward-facing surfaces value at topography top (k-1) is set,
392!--          for downward-facing surfaces at topography bottom (k+1).
393             kb = MERGE( -1, 1, l == 0 )
394             !$OMP PARALLEL DO PRIVATE( i, j, k )
395             DO  m = 1, bc_h(l)%ns
396                i = bc_h(l)%i(m)           
397                j = bc_h(l)%j(m)
398                k = bc_h(l)%k(m)
399                q_p(k+kb,j,i) = q(k+kb,j,i)
400             ENDDO
401          ENDDO
402         
403       ELSE
404         
405          DO  l = 0, 1
406!
407!--          Set kb, for upward-facing surfaces value at topography top (k-1) is set,
408!--          for downward-facing surfaces at topography bottom (k+1).
409             kb = MERGE( -1, 1, l == 0 )
410             !$OMP PARALLEL DO PRIVATE( i, j, k )
411             DO  m = 1, bc_h(l)%ns
412                i = bc_h(l)%i(m)           
413                j = bc_h(l)%j(m)
414                k = bc_h(l)%k(m)
415                q_p(k+kb,j,i) = q_p(k,j,i)
416             ENDDO
417          ENDDO
418       ENDIF
419!
420!--    Top boundary
421       IF ( ibc_q_t == 0 ) THEN
422          q_p(nzt+1,:,:) = q(nzt+1,:,:)
423       ELSEIF ( ibc_q_t == 1 ) THEN
424          q_p(nzt+1,:,:) = q_p(nzt,:,:) + bc_q_t_val * dzu(nzt+1)
425       ENDIF
426
427       IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
428!             
429!--       Surface conditions cloud water (Dirichlet)
430!--       Run loop over all non-natural and natural walls. Note, in wall-datatype
431!--       the k coordinate belongs to the atmospheric grid point, therefore, set
432!--       qr_p and nr_p at k-1
433          !$OMP PARALLEL DO PRIVATE( i, j, k )
434          DO  m = 1, bc_h(0)%ns
435             i = bc_h(0)%i(m)           
436             j = bc_h(0)%j(m)
437             k = bc_h(0)%k(m)
438             qc_p(k-1,j,i) = 0.0_wp
439             nc_p(k-1,j,i) = 0.0_wp
440          ENDDO
441!
442!--       Top boundary condition for cloud water (Dirichlet)
443          qc_p(nzt+1,:,:) = 0.0_wp
444          nc_p(nzt+1,:,:) = 0.0_wp
445           
446       ENDIF
447
448       IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
449!             
450!--       Surface conditions rain water (Dirichlet)
451!--       Run loop over all non-natural and natural walls. Note, in wall-datatype
452!--       the k coordinate belongs to the atmospheric grid point, therefore, set
453!--       qr_p and nr_p at k-1
454          !$OMP PARALLEL DO PRIVATE( i, j, k )
455          DO  m = 1, bc_h(0)%ns
456             i = bc_h(0)%i(m)           
457             j = bc_h(0)%j(m)
458             k = bc_h(0)%k(m)
459             qr_p(k-1,j,i) = 0.0_wp
460             nr_p(k-1,j,i) = 0.0_wp
461          ENDDO
462!
463!--       Top boundary condition for rain water (Dirichlet)
464          qr_p(nzt+1,:,:) = 0.0_wp
465          nr_p(nzt+1,:,:) = 0.0_wp
466           
467       ENDIF
468    ENDIF
469!
470!-- Boundary conditions for scalar,
471!-- bottom and top boundary (see also temperature)
472    IF ( passive_scalar )  THEN
473!
474!--    Surface conditions for constant_humidity_flux
475!--    Run loop over all non-natural and natural walls. Note, in wall-datatype
476!--    the k coordinate belongs to the atmospheric grid point, therefore, set
477!--    s_p at k-1
478       IF ( ibc_s_b == 0 ) THEN
479         
480          DO  l = 0, 1
481!
482!--          Set kb, for upward-facing surfaces value at topography top (k-1) is set,
483!--          for downward-facing surfaces at topography bottom (k+1).
484             kb = MERGE( -1, 1, l == 0 )
485             !$OMP PARALLEL DO PRIVATE( i, j, k )
486             DO  m = 1, bc_h(l)%ns
487                i = bc_h(l)%i(m)           
488                j = bc_h(l)%j(m)
489                k = bc_h(l)%k(m)
490                s_p(k+kb,j,i) = s(k+kb,j,i)
491             ENDDO
492          ENDDO
493         
494       ELSE
495         
496          DO  l = 0, 1
497!
498!--          Set kb, for upward-facing surfaces value at topography top (k-1) is set,
499!--          for downward-facing surfaces at topography bottom (k+1).
500             kb = MERGE( -1, 1, l == 0 )
501             !$OMP PARALLEL DO PRIVATE( i, j, k )
502             DO  m = 1, bc_h(l)%ns
503                i = bc_h(l)%i(m)           
504                j = bc_h(l)%j(m)
505                k = bc_h(l)%k(m)
506                s_p(k+kb,j,i) = s_p(k,j,i)
507             ENDDO
508          ENDDO
509       ENDIF
510!
511!--    Top boundary condition
512       IF ( ibc_s_t == 0 )  THEN
513          s_p(nzt+1,:,:) = s(nzt+1,:,:)
514       ELSEIF ( ibc_s_t == 1 )  THEN
515          s_p(nzt+1,:,:) = s_p(nzt,:,:)
516       ELSEIF ( ibc_s_t == 2 )  THEN
517          s_p(nzt+1,:,:) = s_p(nzt,:,:) + bc_s_t_val * dzu(nzt+1)
518       ENDIF
519
520    ENDIF   
521!
522!-- Top/bottom boundary conditions for chemical species
523#if defined( __chem )
524    IF ( air_chemistry )  CALL chem_boundary_conds( 'set_bc_bottomtop' )
525#endif
526!
527!-- In case of inflow or nest boundary at the south boundary the boundary for v
528!-- is at nys and in case of inflow or nest boundary at the left boundary the
529!-- boundary for u is at nxl. Since in prognostic_equations (cache optimized
530!-- version) these levels are handled as a prognostic level, boundary values
531!-- have to be restored here.
532!-- For the SGS-TKE, Neumann boundary conditions are used at the inflow.
533    IF ( inflow_s )  THEN
534       v_p(:,nys,:) = v_p(:,nys-1,:)
535       IF ( .NOT. constant_diffusion ) e_p(:,nys-1,:) = e_p(:,nys,:)
536    ELSEIF ( inflow_n )  THEN
537       IF ( .NOT. constant_diffusion ) e_p(:,nyn+1,:) = e_p(:,nyn,:)
538    ELSEIF ( inflow_l ) THEN
539       u_p(:,:,nxl) = u_p(:,:,nxl-1)
540       IF ( .NOT. constant_diffusion ) e_p(:,:,nxl-1) = e_p(:,:,nxl)
541    ELSEIF ( inflow_r )  THEN
542       IF ( .NOT. constant_diffusion ) e_p(:,:,nxr+1) = e_p(:,:,nxr)
543    ENDIF
544
545!
546!-- The same restoration for u at i=nxl and v at j=nys as above must be made
547!-- in case of nest boundaries. This must not be done in case of vertical nesting
548!-- mode as in that case the lateral boundaries are actually cyclic.
549    IF ( nesting_mode /= 'vertical'  .OR.  forcing )  THEN
550       IF ( nest_bound_s  .OR.  force_bound_s )  THEN
551          v_p(:,nys,:) = v_p(:,nys-1,:)
552       ENDIF
553       IF ( nest_bound_l  .OR.  force_bound_l )  THEN
554          u_p(:,:,nxl) = u_p(:,:,nxl-1)
555       ENDIF
556    ENDIF
557
558!
559!-- Lateral boundary conditions for scalar quantities at the outflow
560    IF ( outflow_s )  THEN
561       pt_p(:,nys-1,:)     = pt_p(:,nys,:)
562       IF ( .NOT. constant_diffusion )  e_p(:,nys-1,:) = e_p(:,nys,:)
563       IF ( rans_tke_e )  diss_p(:,nys-1,:) = diss_p(:,nys,:)
564       IF ( humidity )  THEN
565          q_p(:,nys-1,:) = q_p(:,nys,:)
566          IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
567             qc_p(:,nys-1,:) = qc_p(:,nys,:)
568             nc_p(:,nys-1,:) = nc_p(:,nys,:)
569          ENDIF
570          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
571             qr_p(:,nys-1,:) = qr_p(:,nys,:)
572             nr_p(:,nys-1,:) = nr_p(:,nys,:)
573          ENDIF
574       ENDIF
575       IF ( passive_scalar )  s_p(:,nys-1,:) = s_p(:,nys,:)
576    ELSEIF ( outflow_n )  THEN
577       pt_p(:,nyn+1,:)     = pt_p(:,nyn,:)
578       IF ( .NOT. constant_diffusion )  e_p(:,nyn+1,:) = e_p(:,nyn,:)
579       IF ( rans_tke_e )  diss_p(:,nyn+1,:) = diss_p(:,nyn,:)
580       IF ( humidity )  THEN
581          q_p(:,nyn+1,:) = q_p(:,nyn,:)
582          IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
583             qc_p(:,nyn+1,:) = qc_p(:,nyn,:)
584             nc_p(:,nyn+1,:) = nc_p(:,nyn,:)
585          ENDIF
586          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
587             qr_p(:,nyn+1,:) = qr_p(:,nyn,:)
588             nr_p(:,nyn+1,:) = nr_p(:,nyn,:)
589          ENDIF
590       ENDIF
591       IF ( passive_scalar )  s_p(:,nyn+1,:) = s_p(:,nyn,:)
592    ELSEIF ( outflow_l )  THEN
593       pt_p(:,:,nxl-1)     = pt_p(:,:,nxl)
594       IF ( .NOT. constant_diffusion )  e_p(:,:,nxl-1) = e_p(:,:,nxl)
595       IF ( rans_tke_e )  diss_p(:,:,nxl-1) = diss_p(:,:,nxl)
596       IF ( humidity )  THEN
597          q_p(:,:,nxl-1) = q_p(:,:,nxl)
598          IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
599             qc_p(:,:,nxl-1) = qc_p(:,:,nxl)
600             nc_p(:,:,nxl-1) = nc_p(:,:,nxl)
601          ENDIF
602          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
603             qr_p(:,:,nxl-1) = qr_p(:,:,nxl)
604             nr_p(:,:,nxl-1) = nr_p(:,:,nxl)
605          ENDIF
606       ENDIF
607       IF ( passive_scalar )  s_p(:,:,nxl-1) = s_p(:,:,nxl)
608    ELSEIF ( outflow_r )  THEN
609       pt_p(:,:,nxr+1)     = pt_p(:,:,nxr)
610       IF ( .NOT. constant_diffusion )  e_p(:,:,nxr+1) = e_p(:,:,nxr)
611       IF ( rans_tke_e )  diss_p(:,:,nxr+1) = diss_p(:,:,nxr)
612       IF ( humidity )  THEN
613          q_p(:,:,nxr+1) = q_p(:,:,nxr)
614          IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
615             qc_p(:,:,nxr+1) = qc_p(:,:,nxr)
616             nc_p(:,:,nxr+1) = nc_p(:,:,nxr)
617          ENDIF
618          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
619             qr_p(:,:,nxr+1) = qr_p(:,:,nxr)
620             nr_p(:,:,nxr+1) = nr_p(:,:,nxr)
621          ENDIF
622       ENDIF
623       IF ( passive_scalar )  s_p(:,:,nxr+1) = s_p(:,:,nxr)
624    ENDIF
625
626!
627!-- Lateral boundary conditions for chemical species
628#if defined( __chem )
629    IF ( air_chemistry )  CALL chem_boundary_conds( 'set_bc_lateral' )   
630#endif
631
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.