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

Last change on this file since 3189 was 3183, checked in by suehring, 6 years ago

last commit documented

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