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

Last change on this file since 4125 was 4102, checked in by suehring, 5 years ago

Bugfix, set Neumann boundary conditions for the subgrid TKE at vertical walls instead of implicit Dirichlet conditions that always act as a sink term for the subgrid TKE. Therefore, add new data structure for vertical surfaces and revise the setting of the boundary grid point index space. Moreover, slightly revise setting of boundary conditions at upward- and downward facing surfaces. Finally, setting of boundary conditions for subgrid TKE and dissipation (in RANS mode) is now modularized. Update test case results.

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