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

Last change on this file since 3294 was 3294, checked in by raasch, 3 years ago

modularization of the ocean code

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