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

Last change on this file since 3587 was 3582, checked in by suehring, 5 years ago

Merge branch salsa with trunk

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