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

Last change on this file since 3456 was 3347, checked in by suehring, 6 years ago

Offline nesting revised and separated from large_scale_forcing_mod; Time-dependent synthetic turbulence generator; bugfixes in USM/LSM radiation model and init_pegrid

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