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

Last change on this file since 4087 was 4087, checked in by gronemeier, 2 years ago

Add comment (boundary_conds.f90)

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