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

Last change on this file since 3540 was 3467, checked in by suehring, 6 years ago

Branch salsa @3446 re-integrated into trunk

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