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

Last change on this file since 3864 was 3864, checked in by monakurppa, 2 years ago

major changes in salsa: data input, format and performance

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