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

Last change on this file since 3634 was 3634, checked in by knoop, 3 years ago

OpenACC port for SPEC

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