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

Last change on this file since 3274 was 3274, checked in by knoop, 3 years ago

Modularization of all bulk cloud physics code components

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