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

Last change on this file since 2724 was 2718, checked in by maronga, 7 years ago

deleting of deprecated files; headers updated where needed

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