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

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

variables documented

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