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

Last change on this file since 3272 was 3241, checked in by raasch, 6 years ago

various changes to avoid compiler warnings (mainly removal of unused variables)

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