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

Last change on this file since 3173 was 3129, checked in by gronemeier, 6 years ago

merge with branch rans: update of rans mode and data output

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