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

Last change on this file since 3182 was 3182, checked in by suehring, 3 years ago

New Inifor features: grid stretching, improved command-interface, support start dates in different formats in both YYYYMMDD and YYYYMMDDHH, Ability to manually control input file prefixes (--radiation-prefix, --soil-preifx, --flow-prefix, --soilmoisture-prefix) for compatiblity with DWD forcast naming scheme; GNU-style short and long option; Prepared output of large-scale forcing profiles (no computation yet); Added preprocessor flag netcdf4 to switch output format between netCDF 3 and 4; Updated netCDF variable names and attributes to comply with PIDS v1.9; Inifor bugfixes: Improved compatibility with older Intel Intel compilers by avoiding implicit array allocation; Added origin_lon/_lat values and correct reference time in dynamic driver global attributes; corresponding PALM changes: adjustments to revised Inifor; variables names in dynamic driver adjusted; enable geostrophic forcing also in offline nested mode; variable names in LES-LES and COSMO offline nesting changed; lateral boundary flags for nesting, in- and outflow conditions renamed

  • Property svn:keywords set to Id
File size: 46.2 KB
Line 
1!> @file boundary_conds.f90
2!------------------------------------------------------------------------------!
3! This file is part of the PALM model system.
4!
5! PALM is free software: you can redistribute it and/or modify it under the
6! terms of the GNU General Public License as published by the Free Software
7! Foundation, either version 3 of the License, or (at your option) any later
8! version.
9!
10! PALM is distributed in the hope that it will be useful, but WITHOUT ANY
11! WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
12! A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
13!
14! You should have received a copy of the GNU General Public License along with
15! PALM. If not, see <http://www.gnu.org/licenses/>.
16!
17! Copyright 1997-2018 Leibniz Universitaet Hannover
18!------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22! Rename some variables concerning LES-LES as well as offline nesting
23!
24! Former revisions:
25! -----------------
26! $Id: boundary_conds.f90 3182 2018-07-27 13:36:03Z 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,     &
189               qr_p, s,                                                        & 
190               s_p, sa, sa_p, u, ug, u_init, u_m_l, u_m_n, u_m_r, u_m_s, u_p,  &
191               v, vg, v_init, v_m_l, v_m_n, v_m_r, v_m_s, v_p,                 &
192               w, w_p, w_m_l, w_m_n, w_m_r, w_m_s, pt_init, ddzu
193
194    USE chemistry_model_mod,                                                   &
195        ONLY:  chem_boundary_conds
196             
197    USE control_parameters,                                                    &
198        ONLY:  air_chemistry, bc_dirichlet_l, bc_dirichlet_n, bc_dirichlet_r,  &
199               bc_dirichlet_s, bc_radiation_l, bc_radiation_n, bc_radiation_r, &
200               bc_radiation_s, bc_pt_t_val, bc_q_t_val, bc_s_t_val,            &
201               child_domain, constant_diffusion, cloud_physics, coupling_mode, &
202               dt_3d, humidity, ibc_pt_b, ibc_pt_t, ibc_q_b, ibc_q_t, ibc_s_b, &
203               ibc_s_t,ibc_sa_t, ibc_uv_b, ibc_uv_t,                           &
204               intermediate_timestep_count, kappa,                             &
205               microphysics_morrison, microphysics_seifert,                    &
206               nesting_offline, nudging,                                       &
207               ocean, passive_scalar, rans_mode, rans_tke_e, tsc, use_cmax
208
209    USE grid_variables,                                                        &
210        ONLY:  ddx, ddy, dx, dy
211
212    USE indices,                                                               &
213        ONLY:  nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg,             &
214               nzb, nzt, wall_flags_0
215
216    USE kinds
217
218    USE pegrid
219
220    USE pmc_interface,                                                         &
221        ONLY : nesting_mode, rans_mode_parent
222
223    USE surface_mod,                                                           &
224        ONLY :  bc_h, surf_def_h, surf_def_v, surf_lsm_h, surf_lsm_v,          &
225                surf_usm_h, surf_usm_v
226
227    USE turbulence_closure_mod,                                                &
228        ONLY:  c_0
229
230    IMPLICIT NONE
231
232    INTEGER(iwp) ::  i  !< grid index x direction
233    INTEGER(iwp) ::  j  !< grid index y direction
234    INTEGER(iwp) ::  k  !< grid index z direction
235    INTEGER(iwp) ::  kb !< variable to set respective boundary value, depends on facing.
236    INTEGER(iwp) ::  l  !< running index boundary type, for up- and downward-facing walls
237    INTEGER(iwp) ::  m  !< running index surface elements
238
239    REAL(wp)    ::  c_max !<
240    REAL(wp)    ::  denom !<
241
242
243!
244!-- Bottom boundary
245    IF ( ibc_uv_b == 1 )  THEN
246       u_p(nzb,:,:) = u_p(nzb+1,:,:)
247       v_p(nzb,:,:) = v_p(nzb+1,:,:)
248    ENDIF
249!
250!-- Set zero vertical velocity at topography top (l=0), or bottom (l=1) in case
251!-- of downward-facing surfaces.
252    DO  l = 0, 1
253!
254!--    Set kb, for upward-facing surfaces value at topography top (k-1) is set,
255!--    for downward-facing surfaces at topography bottom (k+1).
256       kb = MERGE( -1, 1, l == 0 )
257       !$OMP PARALLEL DO PRIVATE( i, j, k )
258       DO  m = 1, bc_h(l)%ns
259          i = bc_h(l)%i(m)           
260          j = bc_h(l)%j(m)
261          k = bc_h(l)%k(m)
262          w_p(k+kb,j,i) = 0.0_wp
263       ENDDO
264    ENDDO
265
266!
267!-- Top boundary. A nested domain ( ibc_uv_t = 3 ) does not require settings.
268    IF ( ibc_uv_t == 0 )  THEN
269        u_p(nzt+1,:,:) = u_init(nzt+1)
270        v_p(nzt+1,:,:) = v_init(nzt+1)
271    ELSEIF ( ibc_uv_t == 1 )  THEN
272        u_p(nzt+1,:,:) = u_p(nzt,:,:)
273        v_p(nzt+1,:,:) = v_p(nzt,:,:)
274    ENDIF
275
276!
277!-- Vertical nesting: Vertical velocity not zero at the top of the fine grid
278    IF (  .NOT.  child_domain  .AND.                                           &
279                 TRIM(coupling_mode) /= 'vnested_fine' )  THEN
280       w_p(nzt:nzt+1,:,:) = 0.0_wp  !< nzt is not a prognostic level (but cf. pres)
281    ENDIF
282
283!
284!-- Temperature at bottom and top boundary.
285!-- In case of coupled runs (ibc_pt_b = 2) the temperature is given by
286!-- the sea surface temperature of the coupled ocean model.
287!-- Dirichlet
288    IF ( ibc_pt_b == 0 )  THEN
289       DO  l = 0, 1
290!
291!--       Set kb, for upward-facing surfaces value at topography top (k-1) is set,
292!--       for downward-facing surfaces at topography bottom (k+1).
293          kb = MERGE( -1, 1, l == 0 )
294          !$OMP PARALLEL DO PRIVATE( i, j, k )
295          DO  m = 1, bc_h(l)%ns
296             i = bc_h(l)%i(m)           
297             j = bc_h(l)%j(m)
298             k = bc_h(l)%k(m)
299             pt_p(k+kb,j,i) = pt(k+kb,j,i)
300          ENDDO
301       ENDDO
302!
303!-- Neumann, zero-gradient
304    ELSEIF ( ibc_pt_b == 1 )  THEN
305       DO  l = 0, 1
306!
307!--       Set kb, for upward-facing surfaces value at topography top (k-1) is set,
308!--       for downward-facing surfaces at topography bottom (k+1).
309          kb = MERGE( -1, 1, l == 0 )
310          !$OMP PARALLEL DO PRIVATE( i, j, k )
311          DO  m = 1, bc_h(l)%ns
312             i = bc_h(l)%i(m)           
313             j = bc_h(l)%j(m)
314             k = bc_h(l)%k(m)
315             pt_p(k+kb,j,i) = pt_p(k,j,i)
316          ENDDO
317       ENDDO
318    ENDIF
319
320!
321!-- Temperature at top boundary
322    IF ( ibc_pt_t == 0 )  THEN
323        pt_p(nzt+1,:,:) = pt(nzt+1,:,:)
324!
325!--     In case of nudging adjust top boundary to pt which is
326!--     read in from NUDGING-DATA
327        IF ( nudging )  THEN
328           pt_p(nzt+1,:,:) = pt_init(nzt+1)
329        ENDIF
330    ELSEIF ( ibc_pt_t == 1 )  THEN
331        pt_p(nzt+1,:,:) = pt_p(nzt,:,:)
332    ELSEIF ( ibc_pt_t == 2 )  THEN
333        pt_p(nzt+1,:,:) = pt_p(nzt,:,:) + bc_pt_t_val * dzu(nzt+1)
334    ENDIF
335
336!
337!-- Boundary conditions for TKE.
338!-- Generally Neumann conditions with de/dz=0 are assumed.
339    IF ( .NOT. constant_diffusion )  THEN
340
341       IF ( .NOT. rans_tke_e )  THEN
342          DO  l = 0, 1
343!
344!--         Set kb, for upward-facing surfaces value at topography top (k-1) is set,
345!--         for downward-facing surfaces at topography bottom (k+1).
346             kb = MERGE( -1, 1, l == 0 )
347             !$OMP PARALLEL DO PRIVATE( i, j, k )
348             DO  m = 1, bc_h(l)%ns
349                i = bc_h(l)%i(m)           
350                j = bc_h(l)%j(m)
351                k = bc_h(l)%k(m)
352                e_p(k+kb,j,i) = e_p(k,j,i)
353             ENDDO
354          ENDDO
355       ELSE
356!
357!--       Use wall function within constant-flux layer
358!--       Upward-facing surfaces
359!--       Default surfaces
360          DO  m = 1, surf_def_h(0)%ns
361             i = surf_def_h(0)%i(m)
362             j = surf_def_h(0)%j(m)
363             k = surf_def_h(0)%k(m)
364             e_p(k,j,i) = ( surf_def_h(0)%us(m) / c_0 )**2
365          ENDDO
366!
367!--       Natural surfaces
368          DO  m = 1, surf_lsm_h%ns
369             i = surf_lsm_h%i(m)
370             j = surf_lsm_h%j(m)
371             k = surf_lsm_h%k(m)
372             e_p(k,j,i) = ( surf_lsm_h%us(m) / c_0 )**2
373          ENDDO
374!
375!--       Urban surfaces
376          DO  m = 1, surf_usm_h%ns
377             i = surf_usm_h%i(m)
378             j = surf_usm_h%j(m)
379             k = surf_usm_h%k(m)
380             e_p(k,j,i) = ( surf_usm_h%us(m) / c_0 )**2
381          ENDDO
382!
383!--       Vertical surfaces
384          DO  l = 0, 3
385!
386!--          Default surfaces
387             DO  m = 1, surf_def_v(l)%ns
388                i = surf_def_v(l)%i(m)
389                j = surf_def_v(l)%j(m)
390                k = surf_def_v(l)%k(m)
391                e_p(k,j,i) = ( surf_def_v(l)%us(m) / c_0 )**2
392             ENDDO
393!
394!--          Natural surfaces
395             DO  m = 1, surf_lsm_v(l)%ns
396                i = surf_lsm_v(l)%i(m)
397                j = surf_lsm_v(l)%j(m)
398                k = surf_lsm_v(l)%k(m)
399                e_p(k,j,i) = ( surf_lsm_v(l)%us(m) / c_0 )**2
400             ENDDO
401!
402!--          Urban surfaces
403             DO  m = 1, surf_usm_v(l)%ns
404                i = surf_usm_v(l)%i(m)
405                j = surf_usm_v(l)%j(m)
406                k = surf_usm_v(l)%k(m)
407                e_p(k,j,i) = ( surf_usm_v(l)%us(m) / c_0 )**2
408             ENDDO
409          ENDDO
410       ENDIF
411
412       IF ( .NOT. child_domain )  THEN
413          e_p(nzt+1,:,:) = e_p(nzt,:,:)
414!
415!--    Nesting case: if parent operates in RANS mode and child in LES mode,
416!--    no TKE is transfered. This case, set Neumann conditions at lateral and
417!--    top child boundaries.
418!--    If not ( both either in RANS or in LES mode ), TKE boundary condition
419!--    is treated in the nesting.
420       ELSE
421
422          IF ( rans_mode_parent  .AND.  .NOT. rans_mode )  THEN
423
424             e_p(nzt+1,:,:) = e_p(nzt,:,:)
425             IF ( bc_dirichlet_l )  e_p(:,:,nxl-1) = e_p(:,:,nxl)
426             IF ( bc_dirichlet_r )  e_p(:,:,nxr+1) = e_p(:,:,nxr)
427             IF ( bc_dirichlet_s )  e_p(:,nys-1,:) = e_p(:,nys,:)
428             IF ( bc_dirichlet_n )  e_p(:,nyn+1,:) = e_p(:,nyn,:)
429
430          ENDIF
431       ENDIF
432    ENDIF
433
434!
435!-- Boundary conditions for TKE dissipation rate.
436    IF ( rans_tke_e )  THEN
437!
438!--    Use wall function within constant-flux layer
439!--    Upward-facing surfaces
440!--    Default surfaces
441       DO  m = 1, surf_def_h(0)%ns
442          i = surf_def_h(0)%i(m)
443          j = surf_def_h(0)%j(m)
444          k = surf_def_h(0)%k(m)
445          diss_p(k,j,i) = surf_def_h(0)%us(m)**3          &
446                        / ( kappa * surf_def_h(0)%z_mo(m) )
447       ENDDO
448!
449!--    Natural surfaces
450       DO  m = 1, surf_lsm_h%ns
451          i = surf_lsm_h%i(m)
452          j = surf_lsm_h%j(m)
453          k = surf_lsm_h%k(m)
454          diss_p(k,j,i) = surf_lsm_h%us(m)**3          &
455                        / ( kappa * surf_lsm_h%z_mo(m) )
456       ENDDO
457!
458!--    Urban surfaces
459       DO  m = 1, surf_usm_h%ns
460          i = surf_usm_h%i(m)
461          j = surf_usm_h%j(m)
462          k = surf_usm_h%k(m)
463          diss_p(k,j,i) = surf_usm_h%us(m)**3          &
464                        / ( kappa * surf_usm_h%z_mo(m) )
465       ENDDO
466!
467!--    Vertical surfaces
468       DO  l = 0, 3
469!
470!--       Default surfaces
471          DO  m = 1, surf_def_v(l)%ns
472             i = surf_def_v(l)%i(m)
473             j = surf_def_v(l)%j(m)
474             k = surf_def_v(l)%k(m)
475             diss_p(k,j,i) = surf_def_v(l)%us(m)**3          &
476                           / ( kappa * surf_def_v(l)%z_mo(m) )
477          ENDDO
478!
479!--       Natural surfaces
480          DO  m = 1, surf_lsm_v(l)%ns
481             i = surf_lsm_v(l)%i(m)
482             j = surf_lsm_v(l)%j(m)
483             k = surf_lsm_v(l)%k(m)
484             diss_p(k,j,i) = surf_lsm_v(l)%us(m)**3          &
485                           / ( kappa * surf_lsm_v(l)%z_mo(m) )
486          ENDDO
487!
488!--       Urban surfaces
489          DO  m = 1, surf_usm_v(l)%ns
490             i = surf_usm_v(l)%i(m)
491             j = surf_usm_v(l)%j(m)
492             k = surf_usm_v(l)%k(m)
493             diss_p(k,j,i) = surf_usm_v(l)%us(m)**3          &
494                           / ( kappa * surf_usm_v(l)%z_mo(m) )
495          ENDDO
496       ENDDO
497!
498!--    Limit change of diss to be between -90% and +100%. Also, set an absolute
499!--    minimum value
500       DO  i = nxl, nxr
501          DO  j = nys, nyn
502             DO  k = nzb, nzt+1
503                diss_p(k,j,i) = MAX( MIN( diss_p(k,j,i),          &
504                                          2.0_wp * diss(k,j,i) ), &
505                                     0.1_wp * diss(k,j,i),        &
506                                     0.0001_wp )
507             ENDDO
508          ENDDO
509       ENDDO
510
511       IF ( .NOT. child_domain )  THEN
512          diss_p(nzt+1,:,:) = diss_p(nzt,:,:)
513       ENDIF
514    ENDIF
515
516!
517!-- Boundary conditions for salinity
518    IF ( ocean )  THEN
519!
520!--    Bottom boundary: Neumann condition because salinity flux is always
521!--    given.
522       DO  l = 0, 1
523!
524!--       Set kb, for upward-facing surfaces value at topography top (k-1) is set,
525!--       for downward-facing surfaces at topography bottom (k+1).
526          kb = MERGE( -1, 1, l == 0 )
527          !$OMP PARALLEL DO PRIVATE( i, j, k )
528          DO  m = 1, bc_h(l)%ns
529             i = bc_h(l)%i(m)           
530             j = bc_h(l)%j(m)
531             k = bc_h(l)%k(m)
532             sa_p(k+kb,j,i) = sa_p(k,j,i)
533          ENDDO
534       ENDDO
535!
536!--    Top boundary: Dirichlet or Neumann
537       IF ( ibc_sa_t == 0 )  THEN
538           sa_p(nzt+1,:,:) = sa(nzt+1,:,:)
539       ELSEIF ( ibc_sa_t == 1 )  THEN
540           sa_p(nzt+1,:,:) = sa_p(nzt,:,:)
541       ENDIF
542
543    ENDIF
544
545!
546!-- Boundary conditions for total water content,
547!-- bottom and top boundary (see also temperature)
548    IF ( humidity )  THEN
549!
550!--    Surface conditions for constant_humidity_flux
551!--    Run loop over all non-natural and natural walls. Note, in wall-datatype
552!--    the k coordinate belongs to the atmospheric grid point, therefore, set
553!--    q_p at k-1
554       IF ( ibc_q_b == 0 ) THEN
555
556          DO  l = 0, 1
557!
558!--          Set kb, for upward-facing surfaces value at topography top (k-1) is set,
559!--          for downward-facing surfaces at topography bottom (k+1).
560             kb = MERGE( -1, 1, l == 0 )
561             !$OMP PARALLEL DO PRIVATE( i, j, k )
562             DO  m = 1, bc_h(l)%ns
563                i = bc_h(l)%i(m)           
564                j = bc_h(l)%j(m)
565                k = bc_h(l)%k(m)
566                q_p(k+kb,j,i) = q(k+kb,j,i)
567             ENDDO
568          ENDDO
569         
570       ELSE
571         
572          DO  l = 0, 1
573!
574!--          Set kb, for upward-facing surfaces value at topography top (k-1) is set,
575!--          for downward-facing surfaces at topography bottom (k+1).
576             kb = MERGE( -1, 1, l == 0 )
577             !$OMP PARALLEL DO PRIVATE( i, j, k )
578             DO  m = 1, bc_h(l)%ns
579                i = bc_h(l)%i(m)           
580                j = bc_h(l)%j(m)
581                k = bc_h(l)%k(m)
582                q_p(k+kb,j,i) = q_p(k,j,i)
583             ENDDO
584          ENDDO
585       ENDIF
586!
587!--    Top boundary
588       IF ( ibc_q_t == 0 ) THEN
589          q_p(nzt+1,:,:) = q(nzt+1,:,:)
590       ELSEIF ( ibc_q_t == 1 ) THEN
591          q_p(nzt+1,:,:) = q_p(nzt,:,:) + bc_q_t_val * dzu(nzt+1)
592       ENDIF
593
594       IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
595!             
596!--       Surface conditions cloud water (Dirichlet)
597!--       Run loop over all non-natural and natural walls. Note, in wall-datatype
598!--       the k coordinate belongs to the atmospheric grid point, therefore, set
599!--       qr_p and nr_p at k-1
600          !$OMP PARALLEL DO PRIVATE( i, j, k )
601          DO  m = 1, bc_h(0)%ns
602             i = bc_h(0)%i(m)           
603             j = bc_h(0)%j(m)
604             k = bc_h(0)%k(m)
605             qc_p(k-1,j,i) = 0.0_wp
606             nc_p(k-1,j,i) = 0.0_wp
607          ENDDO
608!
609!--       Top boundary condition for cloud water (Dirichlet)
610          qc_p(nzt+1,:,:) = 0.0_wp
611          nc_p(nzt+1,:,:) = 0.0_wp
612           
613       ENDIF
614
615       IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
616!             
617!--       Surface conditions rain water (Dirichlet)
618!--       Run loop over all non-natural and natural walls. Note, in wall-datatype
619!--       the k coordinate belongs to the atmospheric grid point, therefore, set
620!--       qr_p and nr_p at k-1
621          !$OMP PARALLEL DO PRIVATE( i, j, k )
622          DO  m = 1, bc_h(0)%ns
623             i = bc_h(0)%i(m)           
624             j = bc_h(0)%j(m)
625             k = bc_h(0)%k(m)
626             qr_p(k-1,j,i) = 0.0_wp
627             nr_p(k-1,j,i) = 0.0_wp
628          ENDDO
629!
630!--       Top boundary condition for rain water (Dirichlet)
631          qr_p(nzt+1,:,:) = 0.0_wp
632          nr_p(nzt+1,:,:) = 0.0_wp
633           
634       ENDIF
635    ENDIF
636!
637!-- Boundary conditions for scalar,
638!-- bottom and top boundary (see also temperature)
639    IF ( passive_scalar )  THEN
640!
641!--    Surface conditions for constant_humidity_flux
642!--    Run loop over all non-natural and natural walls. Note, in wall-datatype
643!--    the k coordinate belongs to the atmospheric grid point, therefore, set
644!--    s_p at k-1
645       IF ( ibc_s_b == 0 ) THEN
646         
647          DO  l = 0, 1
648!
649!--          Set kb, for upward-facing surfaces value at topography top (k-1) is set,
650!--          for downward-facing surfaces at topography bottom (k+1).
651             kb = MERGE( -1, 1, l == 0 )
652             !$OMP PARALLEL DO PRIVATE( i, j, k )
653             DO  m = 1, bc_h(l)%ns
654                i = bc_h(l)%i(m)           
655                j = bc_h(l)%j(m)
656                k = bc_h(l)%k(m)
657                s_p(k+kb,j,i) = s(k+kb,j,i)
658             ENDDO
659          ENDDO
660         
661       ELSE
662         
663          DO  l = 0, 1
664!
665!--          Set kb, for upward-facing surfaces value at topography top (k-1) is set,
666!--          for downward-facing surfaces at topography bottom (k+1).
667             kb = MERGE( -1, 1, l == 0 )
668             !$OMP PARALLEL DO PRIVATE( i, j, k )
669             DO  m = 1, bc_h(l)%ns
670                i = bc_h(l)%i(m)           
671                j = bc_h(l)%j(m)
672                k = bc_h(l)%k(m)
673                s_p(k+kb,j,i) = s_p(k,j,i)
674             ENDDO
675          ENDDO
676       ENDIF
677!
678!--    Top boundary condition
679       IF ( ibc_s_t == 0 )  THEN
680          s_p(nzt+1,:,:) = s(nzt+1,:,:)
681       ELSEIF ( ibc_s_t == 1 )  THEN
682          s_p(nzt+1,:,:) = s_p(nzt,:,:)
683       ELSEIF ( ibc_s_t == 2 )  THEN
684          s_p(nzt+1,:,:) = s_p(nzt,:,:) + bc_s_t_val * dzu(nzt+1)
685       ENDIF
686
687    ENDIF   
688!
689!-- Top/bottom boundary conditions for chemical species
690    IF ( air_chemistry )  CALL chem_boundary_conds( 'set_bc_bottomtop' )
691!
692!-- In case of inflow or nest boundary at the south boundary the boundary for v
693!-- is at nys and in case of inflow or nest boundary at the left boundary the
694!-- boundary for u is at nxl. Since in prognostic_equations (cache optimized
695!-- version) these levels are handled as a prognostic level, boundary values
696!-- have to be restored here.
697!-- For the SGS-TKE, Neumann boundary conditions are used at the inflow.
698    IF ( bc_dirichlet_s )  THEN
699       v_p(:,nys,:) = v_p(:,nys-1,:)
700       IF ( .NOT. constant_diffusion ) e_p(:,nys-1,:) = e_p(:,nys,:)
701    ELSEIF ( bc_dirichlet_n )  THEN
702       IF ( .NOT. constant_diffusion ) e_p(:,nyn+1,:) = e_p(:,nyn,:)
703    ELSEIF ( bc_dirichlet_l ) THEN
704       u_p(:,:,nxl) = u_p(:,:,nxl-1)
705       IF ( .NOT. constant_diffusion ) e_p(:,:,nxl-1) = e_p(:,:,nxl)
706    ELSEIF ( bc_dirichlet_r )  THEN
707       IF ( .NOT. constant_diffusion ) e_p(:,:,nxr+1) = e_p(:,:,nxr)
708    ENDIF
709
710!
711!-- The same restoration for u at i=nxl and v at j=nys as above must be made
712!-- in case of nest boundaries. This must not be done in case of vertical nesting
713!-- mode as in that case the lateral boundaries are actually cyclic.
714!-- @todo: Is this really needed? Boundary values will be overwritten in
715!--        coupler or by Inifor data.
716    IF ( nesting_mode /= 'vertical'  .OR.  nesting_offline )  THEN
717       IF ( bc_dirichlet_s )  THEN
718          v_p(:,nys,:) = v_p(:,nys-1,:)
719       ENDIF
720       IF ( bc_dirichlet_l )  THEN
721          u_p(:,:,nxl) = u_p(:,:,nxl-1)
722       ENDIF
723    ENDIF
724
725!
726!-- Lateral boundary conditions for scalar quantities at the outflow
727    IF ( bc_radiation_s )  THEN
728       pt_p(:,nys-1,:)     = pt_p(:,nys,:)
729       IF ( .NOT. constant_diffusion )  e_p(:,nys-1,:) = e_p(:,nys,:)
730       IF ( rans_tke_e )  diss_p(:,nys-1,:) = diss_p(:,nys,:)
731       IF ( humidity )  THEN
732          q_p(:,nys-1,:) = q_p(:,nys,:)
733          IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
734             qc_p(:,nys-1,:) = qc_p(:,nys,:)
735             nc_p(:,nys-1,:) = nc_p(:,nys,:)
736          ENDIF
737          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
738             qr_p(:,nys-1,:) = qr_p(:,nys,:)
739             nr_p(:,nys-1,:) = nr_p(:,nys,:)
740          ENDIF
741       ENDIF
742       IF ( passive_scalar )  s_p(:,nys-1,:) = s_p(:,nys,:)
743    ELSEIF ( bc_radiation_n )  THEN
744       pt_p(:,nyn+1,:)     = pt_p(:,nyn,:)
745       IF ( .NOT. constant_diffusion )  e_p(:,nyn+1,:) = e_p(:,nyn,:)
746       IF ( rans_tke_e )  diss_p(:,nyn+1,:) = diss_p(:,nyn,:)
747       IF ( humidity )  THEN
748          q_p(:,nyn+1,:) = q_p(:,nyn,:)
749          IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
750             qc_p(:,nyn+1,:) = qc_p(:,nyn,:)
751             nc_p(:,nyn+1,:) = nc_p(:,nyn,:)
752          ENDIF
753          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
754             qr_p(:,nyn+1,:) = qr_p(:,nyn,:)
755             nr_p(:,nyn+1,:) = nr_p(:,nyn,:)
756          ENDIF
757       ENDIF
758       IF ( passive_scalar )  s_p(:,nyn+1,:) = s_p(:,nyn,:)
759    ELSEIF ( bc_radiation_l )  THEN
760       pt_p(:,:,nxl-1)     = pt_p(:,:,nxl)
761       IF ( .NOT. constant_diffusion )  e_p(:,:,nxl-1) = e_p(:,:,nxl)
762       IF ( rans_tke_e )  diss_p(:,:,nxl-1) = diss_p(:,:,nxl)
763       IF ( humidity )  THEN
764          q_p(:,:,nxl-1) = q_p(:,:,nxl)
765          IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
766             qc_p(:,:,nxl-1) = qc_p(:,:,nxl)
767             nc_p(:,:,nxl-1) = nc_p(:,:,nxl)
768          ENDIF
769          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
770             qr_p(:,:,nxl-1) = qr_p(:,:,nxl)
771             nr_p(:,:,nxl-1) = nr_p(:,:,nxl)
772          ENDIF
773       ENDIF
774       IF ( passive_scalar )  s_p(:,:,nxl-1) = s_p(:,:,nxl)
775    ELSEIF ( bc_radiation_r )  THEN
776       pt_p(:,:,nxr+1)     = pt_p(:,:,nxr)
777       IF ( .NOT. constant_diffusion )  e_p(:,:,nxr+1) = e_p(:,:,nxr)
778       IF ( rans_tke_e )  diss_p(:,:,nxr+1) = diss_p(:,:,nxr)
779       IF ( humidity )  THEN
780          q_p(:,:,nxr+1) = q_p(:,:,nxr)
781          IF ( cloud_physics  .AND.  microphysics_morrison )  THEN
782             qc_p(:,:,nxr+1) = qc_p(:,:,nxr)
783             nc_p(:,:,nxr+1) = nc_p(:,:,nxr)
784          ENDIF
785          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
786             qr_p(:,:,nxr+1) = qr_p(:,:,nxr)
787             nr_p(:,:,nxr+1) = nr_p(:,:,nxr)
788          ENDIF
789       ENDIF
790       IF ( passive_scalar )  s_p(:,:,nxr+1) = s_p(:,:,nxr)
791    ENDIF
792
793!
794!-- Lateral boundary conditions for chemical species
795    IF ( air_chemistry )  CALL chem_boundary_conds( 'set_bc_lateral' )   
796
797!
798!-- Radiation boundary conditions for the velocities at the respective outflow.
799!-- The phase velocity is either assumed to the maximum phase velocity that
800!-- ensures numerical stability (CFL-condition) or calculated after
801!-- Orlanski(1976) and averaged along the outflow boundary.
802    IF ( bc_radiation_s )  THEN
803
804       IF ( use_cmax )  THEN
805          u_p(:,-1,:) = u(:,0,:)
806          v_p(:,0,:)  = v(:,1,:)
807          w_p(:,-1,:) = w(:,0,:)         
808       ELSEIF ( .NOT. use_cmax )  THEN
809
810          c_max = dy / dt_3d
811
812          c_u_m_l = 0.0_wp 
813          c_v_m_l = 0.0_wp
814          c_w_m_l = 0.0_wp
815
816          c_u_m = 0.0_wp 
817          c_v_m = 0.0_wp
818          c_w_m = 0.0_wp
819
820!
821!--       Calculate the phase speeds for u, v, and w, first local and then
822!--       average along the outflow boundary.
823          DO  k = nzb+1, nzt+1
824             DO  i = nxl, nxr
825
826                denom = u_m_s(k,0,i) - u_m_s(k,1,i)
827
828                IF ( denom /= 0.0_wp )  THEN
829                   c_u(k,i) = -c_max * ( u(k,0,i) - u_m_s(k,0,i) ) / ( denom * tsc(2) )
830                   IF ( c_u(k,i) < 0.0_wp )  THEN
831                      c_u(k,i) = 0.0_wp
832                   ELSEIF ( c_u(k,i) > c_max )  THEN
833                      c_u(k,i) = c_max
834                   ENDIF
835                ELSE
836                   c_u(k,i) = c_max
837                ENDIF
838
839                denom = v_m_s(k,1,i) - v_m_s(k,2,i)
840
841                IF ( denom /= 0.0_wp )  THEN
842                   c_v(k,i) = -c_max * ( v(k,1,i) - v_m_s(k,1,i) ) / ( denom * tsc(2) )
843                   IF ( c_v(k,i) < 0.0_wp )  THEN
844                      c_v(k,i) = 0.0_wp
845                   ELSEIF ( c_v(k,i) > c_max )  THEN
846                      c_v(k,i) = c_max
847                   ENDIF
848                ELSE
849                   c_v(k,i) = c_max
850                ENDIF
851
852                denom = w_m_s(k,0,i) - w_m_s(k,1,i)
853
854                IF ( denom /= 0.0_wp )  THEN
855                   c_w(k,i) = -c_max * ( w(k,0,i) - w_m_s(k,0,i) ) / ( denom * tsc(2) )
856                   IF ( c_w(k,i) < 0.0_wp )  THEN
857                      c_w(k,i) = 0.0_wp
858                   ELSEIF ( c_w(k,i) > c_max )  THEN
859                      c_w(k,i) = c_max
860                   ENDIF
861                ELSE
862                   c_w(k,i) = c_max
863                ENDIF
864
865                c_u_m_l(k) = c_u_m_l(k) + c_u(k,i)
866                c_v_m_l(k) = c_v_m_l(k) + c_v(k,i)
867                c_w_m_l(k) = c_w_m_l(k) + c_w(k,i)
868
869             ENDDO
870          ENDDO
871
872#if defined( __parallel )   
873          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
874          CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
875                              MPI_SUM, comm1dx, ierr )   
876          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
877          CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
878                              MPI_SUM, comm1dx, ierr ) 
879          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
880          CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
881                              MPI_SUM, comm1dx, ierr ) 
882#else
883          c_u_m = c_u_m_l
884          c_v_m = c_v_m_l
885          c_w_m = c_w_m_l
886#endif
887
888          c_u_m = c_u_m / (nx+1)
889          c_v_m = c_v_m / (nx+1)
890          c_w_m = c_w_m / (nx+1)
891
892!
893!--       Save old timelevels for the next timestep
894          IF ( intermediate_timestep_count == 1 )  THEN
895             u_m_s(:,:,:) = u(:,0:1,:)
896             v_m_s(:,:,:) = v(:,1:2,:)
897             w_m_s(:,:,:) = w(:,0:1,:)
898          ENDIF
899
900!
901!--       Calculate the new velocities
902          DO  k = nzb+1, nzt+1
903             DO  i = nxlg, nxrg
904                u_p(k,-1,i) = u(k,-1,i) - dt_3d * tsc(2) * c_u_m(k) *          &
905                                       ( u(k,-1,i) - u(k,0,i) ) * ddy
906
907                v_p(k,0,i)  = v(k,0,i)  - dt_3d * tsc(2) * c_v_m(k) *          &
908                                       ( v(k,0,i) - v(k,1,i) ) * ddy
909
910                w_p(k,-1,i) = w(k,-1,i) - dt_3d * tsc(2) * c_w_m(k) *          &
911                                       ( w(k,-1,i) - w(k,0,i) ) * ddy
912             ENDDO
913          ENDDO
914
915!
916!--       Bottom boundary at the outflow
917          IF ( ibc_uv_b == 0 )  THEN
918             u_p(nzb,-1,:) = 0.0_wp 
919             v_p(nzb,0,:)  = 0.0_wp 
920          ELSE                   
921             u_p(nzb,-1,:) =  u_p(nzb+1,-1,:)
922             v_p(nzb,0,:)  =  v_p(nzb+1,0,:)
923          ENDIF
924          w_p(nzb,-1,:) = 0.0_wp
925
926!
927!--       Top boundary at the outflow
928          IF ( ibc_uv_t == 0 )  THEN
929             u_p(nzt+1,-1,:) = u_init(nzt+1)
930             v_p(nzt+1,0,:)  = v_init(nzt+1)
931          ELSE
932             u_p(nzt+1,-1,:) = u_p(nzt,-1,:)
933             v_p(nzt+1,0,:)  = v_p(nzt,0,:)
934          ENDIF
935          w_p(nzt:nzt+1,-1,:) = 0.0_wp
936
937       ENDIF
938
939    ENDIF
940
941    IF ( bc_radiation_n )  THEN
942
943       IF ( use_cmax )  THEN
944          u_p(:,ny+1,:) = u(:,ny,:)
945          v_p(:,ny+1,:) = v(:,ny,:)
946          w_p(:,ny+1,:) = w(:,ny,:)         
947       ELSEIF ( .NOT. use_cmax )  THEN
948
949          c_max = dy / dt_3d
950
951          c_u_m_l = 0.0_wp 
952          c_v_m_l = 0.0_wp
953          c_w_m_l = 0.0_wp
954
955          c_u_m = 0.0_wp 
956          c_v_m = 0.0_wp
957          c_w_m = 0.0_wp
958
959!
960!--       Calculate the phase speeds for u, v, and w, first local and then
961!--       average along the outflow boundary.
962          DO  k = nzb+1, nzt+1
963             DO  i = nxl, nxr
964
965                denom = u_m_n(k,ny,i) - u_m_n(k,ny-1,i)
966
967                IF ( denom /= 0.0_wp )  THEN
968                   c_u(k,i) = -c_max * ( u(k,ny,i) - u_m_n(k,ny,i) ) / ( denom * tsc(2) )
969                   IF ( c_u(k,i) < 0.0_wp )  THEN
970                      c_u(k,i) = 0.0_wp
971                   ELSEIF ( c_u(k,i) > c_max )  THEN
972                      c_u(k,i) = c_max
973                   ENDIF
974                ELSE
975                   c_u(k,i) = c_max
976                ENDIF
977
978                denom = v_m_n(k,ny,i) - v_m_n(k,ny-1,i)
979
980                IF ( denom /= 0.0_wp )  THEN
981                   c_v(k,i) = -c_max * ( v(k,ny,i) - v_m_n(k,ny,i) ) / ( denom * tsc(2) )
982                   IF ( c_v(k,i) < 0.0_wp )  THEN
983                      c_v(k,i) = 0.0_wp
984                   ELSEIF ( c_v(k,i) > c_max )  THEN
985                      c_v(k,i) = c_max
986                   ENDIF
987                ELSE
988                   c_v(k,i) = c_max
989                ENDIF
990
991                denom = w_m_n(k,ny,i) - w_m_n(k,ny-1,i)
992
993                IF ( denom /= 0.0_wp )  THEN
994                   c_w(k,i) = -c_max * ( w(k,ny,i) - w_m_n(k,ny,i) ) / ( denom * tsc(2) )
995                   IF ( c_w(k,i) < 0.0_wp )  THEN
996                      c_w(k,i) = 0.0_wp
997                   ELSEIF ( c_w(k,i) > c_max )  THEN
998                      c_w(k,i) = c_max
999                   ENDIF
1000                ELSE
1001                   c_w(k,i) = c_max
1002                ENDIF
1003
1004                c_u_m_l(k) = c_u_m_l(k) + c_u(k,i)
1005                c_v_m_l(k) = c_v_m_l(k) + c_v(k,i)
1006                c_w_m_l(k) = c_w_m_l(k) + c_w(k,i)
1007
1008             ENDDO
1009          ENDDO
1010
1011#if defined( __parallel )   
1012          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
1013          CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
1014                              MPI_SUM, comm1dx, ierr )   
1015          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
1016          CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
1017                              MPI_SUM, comm1dx, ierr ) 
1018          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
1019          CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
1020                              MPI_SUM, comm1dx, ierr ) 
1021#else
1022          c_u_m = c_u_m_l
1023          c_v_m = c_v_m_l
1024          c_w_m = c_w_m_l
1025#endif
1026
1027          c_u_m = c_u_m / (nx+1)
1028          c_v_m = c_v_m / (nx+1)
1029          c_w_m = c_w_m / (nx+1)
1030
1031!
1032!--       Save old timelevels for the next timestep
1033          IF ( intermediate_timestep_count == 1 )  THEN
1034                u_m_n(:,:,:) = u(:,ny-1:ny,:)
1035                v_m_n(:,:,:) = v(:,ny-1:ny,:)
1036                w_m_n(:,:,:) = w(:,ny-1:ny,:)
1037          ENDIF
1038
1039!
1040!--       Calculate the new velocities
1041          DO  k = nzb+1, nzt+1
1042             DO  i = nxlg, nxrg
1043                u_p(k,ny+1,i) = u(k,ny+1,i) - dt_3d * tsc(2) * c_u_m(k) *      &
1044                                       ( u(k,ny+1,i) - u(k,ny,i) ) * ddy
1045
1046                v_p(k,ny+1,i) = v(k,ny+1,i)  - dt_3d * tsc(2) * c_v_m(k) *     &
1047                                       ( v(k,ny+1,i) - v(k,ny,i) ) * ddy
1048
1049                w_p(k,ny+1,i) = w(k,ny+1,i) - dt_3d * tsc(2) * c_w_m(k) *      &
1050                                       ( w(k,ny+1,i) - w(k,ny,i) ) * ddy
1051             ENDDO
1052          ENDDO
1053
1054!
1055!--       Bottom boundary at the outflow
1056          IF ( ibc_uv_b == 0 )  THEN
1057             u_p(nzb,ny+1,:) = 0.0_wp
1058             v_p(nzb,ny+1,:) = 0.0_wp   
1059          ELSE                   
1060             u_p(nzb,ny+1,:) =  u_p(nzb+1,ny+1,:)
1061             v_p(nzb,ny+1,:) =  v_p(nzb+1,ny+1,:)
1062          ENDIF
1063          w_p(nzb,ny+1,:) = 0.0_wp
1064
1065!
1066!--       Top boundary at the outflow
1067          IF ( ibc_uv_t == 0 )  THEN
1068             u_p(nzt+1,ny+1,:) = u_init(nzt+1)
1069             v_p(nzt+1,ny+1,:) = v_init(nzt+1)
1070          ELSE
1071             u_p(nzt+1,ny+1,:) = u_p(nzt,nyn+1,:)
1072             v_p(nzt+1,ny+1,:) = v_p(nzt,nyn+1,:)
1073          ENDIF
1074          w_p(nzt:nzt+1,ny+1,:) = 0.0_wp
1075
1076       ENDIF
1077
1078    ENDIF
1079
1080    IF ( bc_radiation_l )  THEN
1081
1082       IF ( use_cmax )  THEN
1083          u_p(:,:,0)  = u(:,:,1)
1084          v_p(:,:,-1) = v(:,:,0)
1085          w_p(:,:,-1) = w(:,:,0)         
1086       ELSEIF ( .NOT. use_cmax )  THEN
1087
1088          c_max = dx / dt_3d
1089
1090          c_u_m_l = 0.0_wp 
1091          c_v_m_l = 0.0_wp
1092          c_w_m_l = 0.0_wp
1093
1094          c_u_m = 0.0_wp 
1095          c_v_m = 0.0_wp
1096          c_w_m = 0.0_wp
1097
1098!
1099!--       Calculate the phase speeds for u, v, and w, first local and then
1100!--       average along the outflow boundary.
1101          DO  k = nzb+1, nzt+1
1102             DO  j = nys, nyn
1103
1104                denom = u_m_l(k,j,1) - u_m_l(k,j,2)
1105
1106                IF ( denom /= 0.0_wp )  THEN
1107                   c_u(k,j) = -c_max * ( u(k,j,1) - u_m_l(k,j,1) ) / ( denom * tsc(2) )
1108                   IF ( c_u(k,j) < 0.0_wp )  THEN
1109                      c_u(k,j) = 0.0_wp
1110                   ELSEIF ( c_u(k,j) > c_max )  THEN
1111                      c_u(k,j) = c_max
1112                   ENDIF
1113                ELSE
1114                   c_u(k,j) = c_max
1115                ENDIF
1116
1117                denom = v_m_l(k,j,0) - v_m_l(k,j,1)
1118
1119                IF ( denom /= 0.0_wp )  THEN
1120                   c_v(k,j) = -c_max * ( v(k,j,0) - v_m_l(k,j,0) ) / ( denom * tsc(2) )
1121                   IF ( c_v(k,j) < 0.0_wp )  THEN
1122                      c_v(k,j) = 0.0_wp
1123                   ELSEIF ( c_v(k,j) > c_max )  THEN
1124                      c_v(k,j) = c_max
1125                   ENDIF
1126                ELSE
1127                   c_v(k,j) = c_max
1128                ENDIF
1129
1130                denom = w_m_l(k,j,0) - w_m_l(k,j,1)
1131
1132                IF ( denom /= 0.0_wp )  THEN
1133                   c_w(k,j) = -c_max * ( w(k,j,0) - w_m_l(k,j,0) ) / ( denom * tsc(2) )
1134                   IF ( c_w(k,j) < 0.0_wp )  THEN
1135                      c_w(k,j) = 0.0_wp
1136                   ELSEIF ( c_w(k,j) > c_max )  THEN
1137                      c_w(k,j) = c_max
1138                   ENDIF
1139                ELSE
1140                   c_w(k,j) = c_max
1141                ENDIF
1142
1143                c_u_m_l(k) = c_u_m_l(k) + c_u(k,j)
1144                c_v_m_l(k) = c_v_m_l(k) + c_v(k,j)
1145                c_w_m_l(k) = c_w_m_l(k) + c_w(k,j)
1146
1147             ENDDO
1148          ENDDO
1149
1150#if defined( __parallel )   
1151          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
1152          CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
1153                              MPI_SUM, comm1dy, ierr )   
1154          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
1155          CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
1156                              MPI_SUM, comm1dy, ierr ) 
1157          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
1158          CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
1159                              MPI_SUM, comm1dy, ierr ) 
1160#else
1161          c_u_m = c_u_m_l
1162          c_v_m = c_v_m_l
1163          c_w_m = c_w_m_l
1164#endif
1165
1166          c_u_m = c_u_m / (ny+1)
1167          c_v_m = c_v_m / (ny+1)
1168          c_w_m = c_w_m / (ny+1)
1169
1170!
1171!--       Save old timelevels for the next timestep
1172          IF ( intermediate_timestep_count == 1 )  THEN
1173                u_m_l(:,:,:) = u(:,:,1:2)
1174                v_m_l(:,:,:) = v(:,:,0:1)
1175                w_m_l(:,:,:) = w(:,:,0:1)
1176          ENDIF
1177
1178!
1179!--       Calculate the new velocities
1180          DO  k = nzb+1, nzt+1
1181             DO  j = nysg, nyng
1182                u_p(k,j,0) = u(k,j,0) - dt_3d * tsc(2) * c_u_m(k) *            &
1183                                       ( u(k,j,0) - u(k,j,1) ) * ddx
1184
1185                v_p(k,j,-1) = v(k,j,-1) - dt_3d * tsc(2) * c_v_m(k) *          &
1186                                       ( v(k,j,-1) - v(k,j,0) ) * ddx
1187
1188                w_p(k,j,-1) = w(k,j,-1) - dt_3d * tsc(2) * c_w_m(k) *          &
1189                                       ( w(k,j,-1) - w(k,j,0) ) * ddx
1190             ENDDO
1191          ENDDO
1192
1193!
1194!--       Bottom boundary at the outflow
1195          IF ( ibc_uv_b == 0 )  THEN
1196             u_p(nzb,:,0)  = 0.0_wp 
1197             v_p(nzb,:,-1) = 0.0_wp
1198          ELSE                   
1199             u_p(nzb,:,0)  =  u_p(nzb+1,:,0)
1200             v_p(nzb,:,-1) =  v_p(nzb+1,:,-1)
1201          ENDIF
1202          w_p(nzb,:,-1) = 0.0_wp
1203
1204!
1205!--       Top boundary at the outflow
1206          IF ( ibc_uv_t == 0 )  THEN
1207             u_p(nzt+1,:,0)  = u_init(nzt+1)
1208             v_p(nzt+1,:,-1) = v_init(nzt+1)
1209          ELSE
1210             u_p(nzt+1,:,0)  = u_p(nzt,:,0)
1211             v_p(nzt+1,:,-1) = v_p(nzt,:,-1)
1212          ENDIF
1213          w_p(nzt:nzt+1,:,-1) = 0.0_wp
1214
1215       ENDIF
1216
1217    ENDIF
1218
1219    IF ( bc_radiation_r )  THEN
1220
1221       IF ( use_cmax )  THEN
1222          u_p(:,:,nx+1) = u(:,:,nx)
1223          v_p(:,:,nx+1) = v(:,:,nx)
1224          w_p(:,:,nx+1) = w(:,:,nx)         
1225       ELSEIF ( .NOT. use_cmax )  THEN
1226
1227          c_max = dx / dt_3d
1228
1229          c_u_m_l = 0.0_wp 
1230          c_v_m_l = 0.0_wp
1231          c_w_m_l = 0.0_wp
1232
1233          c_u_m = 0.0_wp 
1234          c_v_m = 0.0_wp
1235          c_w_m = 0.0_wp
1236
1237!
1238!--       Calculate the phase speeds for u, v, and w, first local and then
1239!--       average along the outflow boundary.
1240          DO  k = nzb+1, nzt+1
1241             DO  j = nys, nyn
1242
1243                denom = u_m_r(k,j,nx) - u_m_r(k,j,nx-1)
1244
1245                IF ( denom /= 0.0_wp )  THEN
1246                   c_u(k,j) = -c_max * ( u(k,j,nx) - u_m_r(k,j,nx) ) / ( denom * tsc(2) )
1247                   IF ( c_u(k,j) < 0.0_wp )  THEN
1248                      c_u(k,j) = 0.0_wp
1249                   ELSEIF ( c_u(k,j) > c_max )  THEN
1250                      c_u(k,j) = c_max
1251                   ENDIF
1252                ELSE
1253                   c_u(k,j) = c_max
1254                ENDIF
1255
1256                denom = v_m_r(k,j,nx) - v_m_r(k,j,nx-1)
1257
1258                IF ( denom /= 0.0_wp )  THEN
1259                   c_v(k,j) = -c_max * ( v(k,j,nx) - v_m_r(k,j,nx) ) / ( denom * tsc(2) )
1260                   IF ( c_v(k,j) < 0.0_wp )  THEN
1261                      c_v(k,j) = 0.0_wp
1262                   ELSEIF ( c_v(k,j) > c_max )  THEN
1263                      c_v(k,j) = c_max
1264                   ENDIF
1265                ELSE
1266                   c_v(k,j) = c_max
1267                ENDIF
1268
1269                denom = w_m_r(k,j,nx) - w_m_r(k,j,nx-1)
1270
1271                IF ( denom /= 0.0_wp )  THEN
1272                   c_w(k,j) = -c_max * ( w(k,j,nx) - w_m_r(k,j,nx) ) / ( denom * tsc(2) )
1273                   IF ( c_w(k,j) < 0.0_wp )  THEN
1274                      c_w(k,j) = 0.0_wp
1275                   ELSEIF ( c_w(k,j) > c_max )  THEN
1276                      c_w(k,j) = c_max
1277                   ENDIF
1278                ELSE
1279                   c_w(k,j) = c_max
1280                ENDIF
1281
1282                c_u_m_l(k) = c_u_m_l(k) + c_u(k,j)
1283                c_v_m_l(k) = c_v_m_l(k) + c_v(k,j)
1284                c_w_m_l(k) = c_w_m_l(k) + c_w(k,j)
1285
1286             ENDDO
1287          ENDDO
1288
1289#if defined( __parallel )   
1290          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
1291          CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
1292                              MPI_SUM, comm1dy, ierr )   
1293          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
1294          CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
1295                              MPI_SUM, comm1dy, ierr ) 
1296          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
1297          CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
1298                              MPI_SUM, comm1dy, ierr ) 
1299#else
1300          c_u_m = c_u_m_l
1301          c_v_m = c_v_m_l
1302          c_w_m = c_w_m_l
1303#endif
1304
1305          c_u_m = c_u_m / (ny+1)
1306          c_v_m = c_v_m / (ny+1)
1307          c_w_m = c_w_m / (ny+1)
1308
1309!
1310!--       Save old timelevels for the next timestep
1311          IF ( intermediate_timestep_count == 1 )  THEN
1312                u_m_r(:,:,:) = u(:,:,nx-1:nx)
1313                v_m_r(:,:,:) = v(:,:,nx-1:nx)
1314                w_m_r(:,:,:) = w(:,:,nx-1:nx)
1315          ENDIF
1316
1317!
1318!--       Calculate the new velocities
1319          DO  k = nzb+1, nzt+1
1320             DO  j = nysg, nyng
1321                u_p(k,j,nx+1) = u(k,j,nx+1) - dt_3d * tsc(2) * c_u_m(k) *      &
1322                                       ( u(k,j,nx+1) - u(k,j,nx) ) * ddx
1323
1324                v_p(k,j,nx+1) = v(k,j,nx+1) - dt_3d * tsc(2) * c_v_m(k) *      &
1325                                       ( v(k,j,nx+1) - v(k,j,nx) ) * ddx
1326
1327                w_p(k,j,nx+1) = w(k,j,nx+1) - dt_3d * tsc(2) * c_w_m(k) *      &
1328                                       ( w(k,j,nx+1) - w(k,j,nx) ) * ddx
1329             ENDDO
1330          ENDDO
1331
1332!
1333!--       Bottom boundary at the outflow
1334          IF ( ibc_uv_b == 0 )  THEN
1335             u_p(nzb,:,nx+1) = 0.0_wp
1336             v_p(nzb,:,nx+1) = 0.0_wp 
1337          ELSE                   
1338             u_p(nzb,:,nx+1) =  u_p(nzb+1,:,nx+1)
1339             v_p(nzb,:,nx+1) =  v_p(nzb+1,:,nx+1)
1340          ENDIF
1341          w_p(nzb,:,nx+1) = 0.0_wp
1342
1343!
1344!--       Top boundary at the outflow
1345          IF ( ibc_uv_t == 0 )  THEN
1346             u_p(nzt+1,:,nx+1) = u_init(nzt+1)
1347             v_p(nzt+1,:,nx+1) = v_init(nzt+1)
1348          ELSE
1349             u_p(nzt+1,:,nx+1) = u_p(nzt,:,nx+1)
1350             v_p(nzt+1,:,nx+1) = v_p(nzt,:,nx+1)
1351          ENDIF
1352          w_p(nzt:nzt+1,:,nx+1) = 0.0_wp
1353
1354       ENDIF
1355
1356    ENDIF
1357
1358 END SUBROUTINE boundary_conds
Note: See TracBrowser for help on using the repository browser.