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

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

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

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