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

Last change on this file since 3630 was 3589, checked in by suehring, 6 years ago

Remove erroneous UTF encoding; last commit documented

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