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

Last change on this file since 2053 was 2001, checked in by knoop, 8 years ago

last commit documented

  • Property svn:keywords set to Id
File size: 33.8 KB
RevLine 
[1682]1!> @file boundary_conds.f90
[2000]2!------------------------------------------------------------------------------!
[1036]3! This file is part of PALM.
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!
[1818]17! Copyright 1997-2016 Leibniz Universitaet Hannover
[2000]18!------------------------------------------------------------------------------!
[1036]19!
[484]20! Current revisions:
[1]21! -----------------
[1933]22!
[2001]23!
[1321]24! Former revisions:
25! -----------------
26! $Id: boundary_conds.f90 2001 2016-08-20 18:41:22Z gronemeier $
27!
[2001]28! 2000 2016-08-20 18:09:15Z knoop
29! Forced header and separation lines into 80 columns
30!
[1993]31! 1992 2016-08-12 15:14:59Z suehring
32! Adjustments for top boundary condition for passive scalar
33!
[1961]34! 1960 2016-07-12 16:34:24Z suehring
35! Treat humidity and passive scalar separately
36!
[1933]37! 1823 2016-04-07 08:57:52Z hoffmann
38! Initial version of purely vertical nesting introduced.
39!
[1823]40! 1822 2016-04-07 07:49:42Z hoffmann
41! icloud_scheme removed. microphyisics_seifert added.
42!
[1765]43! 1764 2016-02-28 12:45:19Z raasch
44! index bug for u_p at left outflow removed
45!
[1763]46! 1762 2016-02-25 12:31:13Z hellstea
47! Introduction of nested domain feature
48!
[1744]49! 1742 2016-01-13 09:50:06Z raasch
50! bugfix for outflow Neumann boundary conditions at bottom and top
51!
[1718]52! 1717 2015-11-11 15:09:47Z raasch
53! Bugfix: index error in outflow conditions for left boundary
54!
[1683]55! 1682 2015-10-07 23:56:08Z knoop
56! Code annotations made doxygen readable
57!
[1717]58! 1410 2014-05-23 12:16:18Z suehring
[1463]59! Bugfix: set dirichlet boundary condition for passive_scalar at model domain
60! top
61!
[1410]62! 1399 2014-05-07 11:16:25Z heinze
63! Bugfix: set inflow boundary conditions also if no humidity or passive_scalar
64! is used.
65!
[1399]66! 1398 2014-05-07 11:15:00Z heinze
67! Dirichlet-condition at the top for u and v changed to u_init and v_init also
68! for large_scale_forcing
69!
[1381]70! 1380 2014-04-28 12:40:45Z heinze
71! Adjust Dirichlet-condition at the top for pt in case of nudging
72!
[1362]73! 1361 2014-04-16 15:17:48Z hoffmann
74! Bottom and top boundary conditions of rain water content (qr) and
75! rain drop concentration (nr) changed to Dirichlet
76!
[1354]77! 1353 2014-04-08 15:21:23Z heinze
78! REAL constants provided with KIND-attribute
79
[1321]80! 1320 2014-03-20 08:40:49Z raasch
[1320]81! ONLY-attribute added to USE-statements,
82! kind-parameters added to all INTEGER and REAL declaration statements,
83! kinds are defined in new module kinds,
84! revision history before 2012 removed,
85! comment fields (!:) to be used for variable explanations added to
86! all variable declaration statements
[1160]87!
[1258]88! 1257 2013-11-08 15:18:40Z raasch
89! loop independent clauses added
90!
[1242]91! 1241 2013-10-30 11:36:58Z heinze
92! Adjust ug and vg at each timestep in case of large_scale_forcing
93!
[1160]94! 1159 2013-05-21 11:58:22Z fricke
[1159]95! Bugfix: Neumann boundary conditions for the velocity components at the
96! outflow are in fact radiation boundary conditions using the maximum phase
97! velocity that ensures numerical stability (CFL-condition).
98! Hence, logical operator use_cmax is now used instead of bc_lr_dirneu/_neudir.
99! Bugfix: In case of use_cmax at the outflow, u, v, w are replaced by
100! u_p, v_p, w_p 
[1116]101!
102! 1115 2013-03-26 18:16:16Z hoffmann
103! boundary conditions of two-moment cloud scheme are restricted to Neumann-
104! boundary-conditions
105!
[1114]106! 1113 2013-03-10 02:48:14Z raasch
107! GPU-porting
108! dummy argument "range" removed
109! Bugfix: wrong index in loops of radiation boundary condition
[1113]110!
[1054]111! 1053 2012-11-13 17:11:03Z hoffmann
112! boundary conditions for the two new prognostic equations (nr, qr) of the
113! two-moment cloud scheme
114!
[1037]115! 1036 2012-10-22 13:43:42Z raasch
116! code put under GPL (PALM 3.9)
117!
[997]118! 996 2012-09-07 10:41:47Z raasch
119! little reformatting
120!
[979]121! 978 2012-08-09 08:28:32Z fricke
122! Neumann boudnary conditions are added at the inflow boundary for the SGS-TKE.
123! Outflow boundary conditions for the velocity components can be set to Neumann
124! conditions or to radiation conditions with a horizontal averaged phase
125! velocity.
126!
[876]127! 875 2012-04-02 15:35:15Z gryschka
128! Bugfix in case of dirichlet inflow bc at the right or north boundary
129!
[1]130! Revision 1.1  1997/09/12 06:21:34  raasch
131! Initial revision
132!
133!
134! Description:
135! ------------
[1682]136!> Boundary conditions for the prognostic quantities.
137!> One additional bottom boundary condition is applied for the TKE (=(u*)**2)
138!> in prandtl_fluxes. The cyclic lateral boundary conditions are implicitly
139!> handled in routine exchange_horiz. Pressure boundary conditions are
140!> explicitly set in routines pres, poisfft, poismg and sor.
[1]141!------------------------------------------------------------------------------!
[1682]142 SUBROUTINE boundary_conds
143 
[1]144
[1320]145    USE arrays_3d,                                                             &
146        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,  &
[1960]147               dzu, e_p, nr_p, pt, pt_p, q, q_p, qr_p, s, s_p, sa, sa_p,       &
[1320]148               u, ug, u_init, u_m_l, u_m_n, u_m_r, u_m_s, u_p,                 &
149               v, vg, v_init, v_m_l, v_m_n, v_m_r, v_m_s, v_p,                 &
[1960]150               w, w_p, w_m_l, w_m_n, w_m_r, w_m_s, pt_init
[1320]151
152    USE control_parameters,                                                    &
[1960]153        ONLY:  bc_pt_t_val, bc_q_t_val, bc_s_t_val, constant_diffusion,        &
[1320]154               cloud_physics, dt_3d, humidity,                                 &
[1960]155               ibc_pt_b, ibc_pt_t, ibc_q_b, ibc_q_t, ibc_s_b, ibc_s_t,         &
156               ibc_sa_t, ibc_uv_b, ibc_uv_t, inflow_l, inflow_n, inflow_r,     &
157               inflow_s, intermediate_timestep_count, large_scale_forcing,     &
[1822]158               microphysics_seifert, nest_domain, nest_bound_l, nest_bound_s,  &
159               nudging, ocean, outflow_l, outflow_n, outflow_r, outflow_s,     &
160               passive_scalar, tsc, use_cmax
[1320]161
162    USE grid_variables,                                                        &
163        ONLY:  ddx, ddy, dx, dy
164
165    USE indices,                                                               &
166        ONLY:  nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg,             &
167               nzb, nzb_s_inner, nzb_w_inner, nzt
168
169    USE kinds
170
[1]171    USE pegrid
172
[1933]173    USE pmc_interface,                                                         &
174        ONLY : nesting_mode
[1320]175
[1933]176
[1]177    IMPLICIT NONE
178
[1682]179    INTEGER(iwp) ::  i !<
180    INTEGER(iwp) ::  j !<
181    INTEGER(iwp) ::  k !<
[1]182
[1682]183    REAL(wp)    ::  c_max !<
184    REAL(wp)    ::  denom !<
[1]185
[73]186
[1]187!
[1113]188!-- Bottom boundary
189    IF ( ibc_uv_b == 1 )  THEN
190       !$acc kernels present( u_p, v_p )
191       u_p(nzb,:,:) = u_p(nzb+1,:,:)
192       v_p(nzb,:,:) = v_p(nzb+1,:,:)
193       !$acc end kernels
194    ENDIF
195
196    !$acc kernels present( nzb_w_inner, w_p )
197    DO  i = nxlg, nxrg
198       DO  j = nysg, nyng
[1353]199          w_p(nzb_w_inner(j,i),j,i) = 0.0_wp
[1113]200       ENDDO
201    ENDDO
202    !$acc end kernels
203
204!
[1762]205!-- Top boundary. A nested domain ( ibc_uv_t = 3 ) does not require settings.
[1113]206    IF ( ibc_uv_t == 0 )  THEN
207       !$acc kernels present( u_init, u_p, v_init, v_p )
208        u_p(nzt+1,:,:) = u_init(nzt+1)
209        v_p(nzt+1,:,:) = v_init(nzt+1)
210       !$acc end kernels
[1762]211    ELSEIF ( ibc_uv_t == 1 )  THEN
[1113]212       !$acc kernels present( u_p, v_p )
213        u_p(nzt+1,:,:) = u_p(nzt,:,:)
214        v_p(nzt+1,:,:) = v_p(nzt,:,:)
215       !$acc end kernels
216    ENDIF
217
[1762]218    IF ( .NOT. nest_domain )  THEN
219       !$acc kernels present( w_p )
220       w_p(nzt:nzt+1,:,:) = 0.0_wp  ! nzt is not a prognostic level (but cf. pres)
221       !$acc end kernels
222    ENDIF
223
[1113]224!
225!-- Temperature at bottom boundary.
226!-- In case of coupled runs (ibc_pt_b = 2) the temperature is given by
227!-- the sea surface temperature of the coupled ocean model.
228    IF ( ibc_pt_b == 0 )  THEN
229       !$acc kernels present( nzb_s_inner, pt, pt_p )
[1257]230       !$acc loop independent
[667]231       DO  i = nxlg, nxrg
[1257]232          !$acc loop independent
[667]233          DO  j = nysg, nyng
[1113]234             pt_p(nzb_s_inner(j,i),j,i) = pt(nzb_s_inner(j,i),j,i)
[1]235          ENDDO
236       ENDDO
[1113]237       !$acc end kernels
238    ELSEIF ( ibc_pt_b == 1 )  THEN
239       !$acc kernels present( nzb_s_inner, pt_p )
[1257]240       !$acc loop independent
[1113]241       DO  i = nxlg, nxrg
[1257]242          !$acc loop independent
[1113]243          DO  j = nysg, nyng
244             pt_p(nzb_s_inner(j,i),j,i) = pt_p(nzb_s_inner(j,i)+1,j,i)
245          ENDDO
246       ENDDO
247      !$acc end kernels
248    ENDIF
[1]249
250!
[1113]251!-- Temperature at top boundary
252    IF ( ibc_pt_t == 0 )  THEN
253       !$acc kernels present( pt, pt_p )
254        pt_p(nzt+1,:,:) = pt(nzt+1,:,:)
[1380]255!
256!--     In case of nudging adjust top boundary to pt which is
257!--     read in from NUDGING-DATA
258        IF ( nudging )  THEN
259           pt_p(nzt+1,:,:) = pt_init(nzt+1)
260        ENDIF
[1113]261       !$acc end kernels
262    ELSEIF ( ibc_pt_t == 1 )  THEN
263       !$acc kernels present( pt_p )
264        pt_p(nzt+1,:,:) = pt_p(nzt,:,:)
265       !$acc end kernels
266    ELSEIF ( ibc_pt_t == 2 )  THEN
267       !$acc kernels present( dzu, pt_p )
[1992]268        pt_p(nzt+1,:,:) = pt_p(nzt,:,:) + bc_pt_t_val * dzu(nzt+1)
[1113]269       !$acc end kernels
270    ENDIF
[1]271
272!
[1113]273!-- Boundary conditions for TKE
274!-- Generally Neumann conditions with de/dz=0 are assumed
275    IF ( .NOT. constant_diffusion )  THEN
276       !$acc kernels present( e_p, nzb_s_inner )
[1257]277       !$acc loop independent
[1113]278       DO  i = nxlg, nxrg
[1257]279          !$acc loop independent
[1113]280          DO  j = nysg, nyng
281             e_p(nzb_s_inner(j,i),j,i) = e_p(nzb_s_inner(j,i)+1,j,i)
[73]282          ENDDO
[1113]283       ENDDO
[1762]284       IF ( .NOT. nest_domain )  THEN
285          e_p(nzt+1,:,:) = e_p(nzt,:,:)
286       ENDIF
[1113]287       !$acc end kernels
288    ENDIF
289
290!
291!-- Boundary conditions for salinity
292    IF ( ocean )  THEN
293!
294!--    Bottom boundary: Neumann condition because salinity flux is always
295!--    given
296       DO  i = nxlg, nxrg
297          DO  j = nysg, nyng
298             sa_p(nzb_s_inner(j,i),j,i) = sa_p(nzb_s_inner(j,i)+1,j,i)
[1]299          ENDDO
[1113]300       ENDDO
[1]301
302!
[1113]303!--    Top boundary: Dirichlet or Neumann
304       IF ( ibc_sa_t == 0 )  THEN
305           sa_p(nzt+1,:,:) = sa(nzt+1,:,:)
306       ELSEIF ( ibc_sa_t == 1 )  THEN
307           sa_p(nzt+1,:,:) = sa_p(nzt,:,:)
[1]308       ENDIF
309
[1113]310    ENDIF
311
[1]312!
[1960]313!-- Boundary conditions for total water content,
[1113]314!-- bottom and top boundary (see also temperature)
[1960]315    IF ( humidity )  THEN
[1113]316!
317!--    Surface conditions for constant_humidity_flux
318       IF ( ibc_q_b == 0 ) THEN
[667]319          DO  i = nxlg, nxrg
320             DO  j = nysg, nyng
[1113]321                q_p(nzb_s_inner(j,i),j,i) = q(nzb_s_inner(j,i),j,i)
[1]322             ENDDO
323          ENDDO
[1113]324       ELSE
[667]325          DO  i = nxlg, nxrg
326             DO  j = nysg, nyng
[1113]327                q_p(nzb_s_inner(j,i),j,i) = q_p(nzb_s_inner(j,i)+1,j,i)
[95]328             ENDDO
329          ENDDO
[1113]330       ENDIF
[95]331!
[1113]332!--    Top boundary
[1462]333       IF ( ibc_q_t == 0 ) THEN
334          q_p(nzt+1,:,:) = q(nzt+1,:,:)
335       ELSEIF ( ibc_q_t == 1 ) THEN
[1992]336          q_p(nzt+1,:,:) = q_p(nzt,:,:) + bc_q_t_val * dzu(nzt+1)
[1462]337       ENDIF
[95]338
[1822]339       IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
[1113]340!             
[1361]341!--       Surface conditions rain water (Dirichlet)
[1115]342          DO  i = nxlg, nxrg
343             DO  j = nysg, nyng
[1361]344                qr_p(nzb_s_inner(j,i),j,i) = 0.0_wp
345                nr_p(nzb_s_inner(j,i),j,i) = 0.0_wp
[73]346             ENDDO
[1115]347          ENDDO
[1]348!
[1361]349!--       Top boundary condition for rain water (Dirichlet)
350          qr_p(nzt+1,:,:) = 0.0_wp
351          nr_p(nzt+1,:,:) = 0.0_wp
[1115]352           
[1]353       ENDIF
[1409]354    ENDIF
[1]355!
[1960]356!-- Boundary conditions for scalar,
357!-- bottom and top boundary (see also temperature)
358    IF ( passive_scalar )  THEN
359!
360!--    Surface conditions for constant_humidity_flux
361       IF ( ibc_s_b == 0 ) THEN
362          DO  i = nxlg, nxrg
363             DO  j = nysg, nyng
364                s_p(nzb_s_inner(j,i),j,i) = s(nzb_s_inner(j,i),j,i)
365             ENDDO
366          ENDDO
367       ELSE
368          DO  i = nxlg, nxrg
369             DO  j = nysg, nyng
370                s_p(nzb_s_inner(j,i),j,i) = s_p(nzb_s_inner(j,i)+1,j,i)
371             ENDDO
372          ENDDO
373       ENDIF
374!
[1992]375!--    Top boundary condition
376       IF ( ibc_s_t == 0 )  THEN
[1960]377          s_p(nzt+1,:,:) = s(nzt+1,:,:)
[1992]378       ELSEIF ( ibc_s_t == 1 )  THEN
379          s_p(nzt+1,:,:) = s_p(nzt,:,:)
380       ELSEIF ( ibc_s_t == 2 )  THEN
381          s_p(nzt+1,:,:) = s_p(nzt,:,:) + bc_s_t_val * dzu(nzt+1)
[1960]382       ENDIF
383
384    ENDIF   
385!
[1762]386!-- In case of inflow or nest boundary at the south boundary the boundary for v
387!-- is at nys and in case of inflow or nest boundary at the left boundary the
388!-- boundary for u is at nxl. Since in prognostic_equations (cache optimized
389!-- version) these levels are handled as a prognostic level, boundary values
390!-- have to be restored here.
[1409]391!-- For the SGS-TKE, Neumann boundary conditions are used at the inflow.
392    IF ( inflow_s )  THEN
393       v_p(:,nys,:) = v_p(:,nys-1,:)
394       IF ( .NOT. constant_diffusion ) e_p(:,nys-1,:) = e_p(:,nys,:)
395    ELSEIF ( inflow_n )  THEN
396       IF ( .NOT. constant_diffusion ) e_p(:,nyn+1,:) = e_p(:,nyn,:)
397    ELSEIF ( inflow_l ) THEN
398       u_p(:,:,nxl) = u_p(:,:,nxl-1)
399       IF ( .NOT. constant_diffusion ) e_p(:,:,nxl-1) = e_p(:,:,nxl)
400    ELSEIF ( inflow_r )  THEN
401       IF ( .NOT. constant_diffusion ) e_p(:,:,nxr+1) = e_p(:,:,nxr)
402    ENDIF
[1]403
404!
[1762]405!-- The same restoration for u at i=nxl and v at j=nys as above must be made
[1933]406!-- in case of nest boundaries. This must not be done in case of vertical nesting
407!-- mode as in that case the lateral boundaries are actually cyclic.
408    IF ( nesting_mode /= 'vertical' )  THEN
409       IF ( nest_bound_s )  THEN
410          v_p(:,nys,:) = v_p(:,nys-1,:)
411       ENDIF
412       IF ( nest_bound_l )  THEN
413          u_p(:,:,nxl) = u_p(:,:,nxl-1)
414       ENDIF
[1762]415    ENDIF
416
417!
[1409]418!-- Lateral boundary conditions for scalar quantities at the outflow
419    IF ( outflow_s )  THEN
420       pt_p(:,nys-1,:)     = pt_p(:,nys,:)
421       IF ( .NOT. constant_diffusion     )  e_p(:,nys-1,:) = e_p(:,nys,:)
[1960]422       IF ( humidity )  THEN
[1409]423          q_p(:,nys-1,:) = q_p(:,nys,:)
[1822]424          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
[1409]425             qr_p(:,nys-1,:) = qr_p(:,nys,:)
426             nr_p(:,nys-1,:) = nr_p(:,nys,:)
[1053]427          ENDIF
[1409]428       ENDIF
[1960]429       IF ( passive_scalar )  s_p(:,nys-1,:) = s_p(:,nys,:)
[1409]430    ELSEIF ( outflow_n )  THEN
431       pt_p(:,nyn+1,:)     = pt_p(:,nyn,:)
432       IF ( .NOT. constant_diffusion     )  e_p(:,nyn+1,:) = e_p(:,nyn,:)
[1960]433       IF ( humidity )  THEN
[1409]434          q_p(:,nyn+1,:) = q_p(:,nyn,:)
[1822]435          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
[1409]436             qr_p(:,nyn+1,:) = qr_p(:,nyn,:)
437             nr_p(:,nyn+1,:) = nr_p(:,nyn,:)
[1053]438          ENDIF
[1409]439       ENDIF
[1960]440       IF ( passive_scalar )  s_p(:,nyn+1,:) = s_p(:,nyn,:)
[1409]441    ELSEIF ( outflow_l )  THEN
442       pt_p(:,:,nxl-1)     = pt_p(:,:,nxl)
443       IF ( .NOT. constant_diffusion     )  e_p(:,:,nxl-1) = e_p(:,:,nxl)
[1960]444       IF ( humidity )  THEN
[1409]445          q_p(:,:,nxl-1) = q_p(:,:,nxl)
[1822]446          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
[1409]447             qr_p(:,:,nxl-1) = qr_p(:,:,nxl)
448             nr_p(:,:,nxl-1) = nr_p(:,:,nxl)
[1053]449          ENDIF
[1409]450       ENDIF
[1960]451       IF ( passive_scalar )  s_p(:,:,nxl-1) = s_p(:,:,nxl)
[1409]452    ELSEIF ( outflow_r )  THEN
453       pt_p(:,:,nxr+1)     = pt_p(:,:,nxr)
454       IF ( .NOT. constant_diffusion     )  e_p(:,:,nxr+1) = e_p(:,:,nxr)
[1960]455       IF ( humidity )  THEN
[1409]456          q_p(:,:,nxr+1) = q_p(:,:,nxr)
[1822]457          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
[1409]458             qr_p(:,:,nxr+1) = qr_p(:,:,nxr)
459             nr_p(:,:,nxr+1) = nr_p(:,:,nxr)
[1053]460          ENDIF
[1]461       ENDIF
[1960]462       IF ( passive_scalar )  s_p(:,:,nxr+1) = s_p(:,:,nxr)
[1]463    ENDIF
464
465!
[1159]466!-- Radiation boundary conditions for the velocities at the respective outflow.
467!-- The phase velocity is either assumed to the maximum phase velocity that
468!-- ensures numerical stability (CFL-condition) or calculated after
469!-- Orlanski(1976) and averaged along the outflow boundary.
[106]470    IF ( outflow_s )  THEN
[75]471
[1159]472       IF ( use_cmax )  THEN
473          u_p(:,-1,:) = u(:,0,:)
474          v_p(:,0,:)  = v(:,1,:)
475          w_p(:,-1,:) = w(:,0,:)         
476       ELSEIF ( .NOT. use_cmax )  THEN
[75]477
[978]478          c_max = dy / dt_3d
[75]479
[1353]480          c_u_m_l = 0.0_wp 
481          c_v_m_l = 0.0_wp
482          c_w_m_l = 0.0_wp
[978]483
[1353]484          c_u_m = 0.0_wp 
485          c_v_m = 0.0_wp
486          c_w_m = 0.0_wp
[978]487
[75]488!
[996]489!--       Calculate the phase speeds for u, v, and w, first local and then
490!--       average along the outflow boundary.
491          DO  k = nzb+1, nzt+1
492             DO  i = nxl, nxr
[75]493
[106]494                denom = u_m_s(k,0,i) - u_m_s(k,1,i)
495
[1353]496                IF ( denom /= 0.0_wp )  THEN
[996]497                   c_u(k,i) = -c_max * ( u(k,0,i) - u_m_s(k,0,i) ) / ( denom * tsc(2) )
[1353]498                   IF ( c_u(k,i) < 0.0_wp )  THEN
499                      c_u(k,i) = 0.0_wp
[106]500                   ELSEIF ( c_u(k,i) > c_max )  THEN
501                      c_u(k,i) = c_max
502                   ENDIF
503                ELSE
504                   c_u(k,i) = c_max
[75]505                ENDIF
506
[106]507                denom = v_m_s(k,1,i) - v_m_s(k,2,i)
508
[1353]509                IF ( denom /= 0.0_wp )  THEN
[996]510                   c_v(k,i) = -c_max * ( v(k,1,i) - v_m_s(k,1,i) ) / ( denom * tsc(2) )
[1353]511                   IF ( c_v(k,i) < 0.0_wp )  THEN
512                      c_v(k,i) = 0.0_wp
[106]513                   ELSEIF ( c_v(k,i) > c_max )  THEN
514                      c_v(k,i) = c_max
515                   ENDIF
516                ELSE
517                   c_v(k,i) = c_max
[75]518                ENDIF
519
[106]520                denom = w_m_s(k,0,i) - w_m_s(k,1,i)
[75]521
[1353]522                IF ( denom /= 0.0_wp )  THEN
[996]523                   c_w(k,i) = -c_max * ( w(k,0,i) - w_m_s(k,0,i) ) / ( denom * tsc(2) )
[1353]524                   IF ( c_w(k,i) < 0.0_wp )  THEN
525                      c_w(k,i) = 0.0_wp
[106]526                   ELSEIF ( c_w(k,i) > c_max )  THEN
527                      c_w(k,i) = c_max
528                   ENDIF
529                ELSE
530                   c_w(k,i) = c_max
[75]531                ENDIF
[106]532
[978]533                c_u_m_l(k) = c_u_m_l(k) + c_u(k,i)
534                c_v_m_l(k) = c_v_m_l(k) + c_v(k,i)
535                c_w_m_l(k) = c_w_m_l(k) + c_w(k,i)
[106]536
[978]537             ENDDO
538          ENDDO
[75]539
[978]540#if defined( __parallel )   
541          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
542          CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
543                              MPI_SUM, comm1dx, ierr )   
544          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
545          CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
546                              MPI_SUM, comm1dx, ierr ) 
547          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
548          CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
549                              MPI_SUM, comm1dx, ierr ) 
550#else
551          c_u_m = c_u_m_l
552          c_v_m = c_v_m_l
553          c_w_m = c_w_m_l
554#endif
555
556          c_u_m = c_u_m / (nx+1)
557          c_v_m = c_v_m / (nx+1)
558          c_w_m = c_w_m / (nx+1)
559
[75]560!
[978]561!--       Save old timelevels for the next timestep
562          IF ( intermediate_timestep_count == 1 )  THEN
563             u_m_s(:,:,:) = u(:,0:1,:)
564             v_m_s(:,:,:) = v(:,1:2,:)
565             w_m_s(:,:,:) = w(:,0:1,:)
566          ENDIF
567
568!
569!--       Calculate the new velocities
[996]570          DO  k = nzb+1, nzt+1
571             DO  i = nxlg, nxrg
[978]572                u_p(k,-1,i) = u(k,-1,i) - dt_3d * tsc(2) * c_u_m(k) *          &
[75]573                                       ( u(k,-1,i) - u(k,0,i) ) * ddy
574
[978]575                v_p(k,0,i)  = v(k,0,i)  - dt_3d * tsc(2) * c_v_m(k) *          &
[106]576                                       ( v(k,0,i) - v(k,1,i) ) * ddy
[75]577
[978]578                w_p(k,-1,i) = w(k,-1,i) - dt_3d * tsc(2) * c_w_m(k) *          &
[75]579                                       ( w(k,-1,i) - w(k,0,i) ) * ddy
[978]580             ENDDO
[75]581          ENDDO
582
583!
[978]584!--       Bottom boundary at the outflow
585          IF ( ibc_uv_b == 0 )  THEN
[1353]586             u_p(nzb,-1,:) = 0.0_wp 
587             v_p(nzb,0,:)  = 0.0_wp 
[978]588          ELSE                   
589             u_p(nzb,-1,:) =  u_p(nzb+1,-1,:)
590             v_p(nzb,0,:)  =  v_p(nzb+1,0,:)
591          ENDIF
[1353]592          w_p(nzb,-1,:) = 0.0_wp
[73]593
[75]594!
[978]595!--       Top boundary at the outflow
596          IF ( ibc_uv_t == 0 )  THEN
597             u_p(nzt+1,-1,:) = u_init(nzt+1)
598             v_p(nzt+1,0,:)  = v_init(nzt+1)
599          ELSE
[1742]600             u_p(nzt+1,-1,:) = u_p(nzt,-1,:)
601             v_p(nzt+1,0,:)  = v_p(nzt,0,:)
[978]602          ENDIF
[1353]603          w_p(nzt:nzt+1,-1,:) = 0.0_wp
[978]604
[75]605       ENDIF
[73]606
[75]607    ENDIF
[73]608
[106]609    IF ( outflow_n )  THEN
[73]610
[1159]611       IF ( use_cmax )  THEN
612          u_p(:,ny+1,:) = u(:,ny,:)
613          v_p(:,ny+1,:) = v(:,ny,:)
614          w_p(:,ny+1,:) = w(:,ny,:)         
615       ELSEIF ( .NOT. use_cmax )  THEN
[75]616
[978]617          c_max = dy / dt_3d
[75]618
[1353]619          c_u_m_l = 0.0_wp 
620          c_v_m_l = 0.0_wp
621          c_w_m_l = 0.0_wp
[978]622
[1353]623          c_u_m = 0.0_wp 
624          c_v_m = 0.0_wp
625          c_w_m = 0.0_wp
[978]626
[1]627!
[996]628!--       Calculate the phase speeds for u, v, and w, first local and then
629!--       average along the outflow boundary.
630          DO  k = nzb+1, nzt+1
631             DO  i = nxl, nxr
[73]632
[106]633                denom = u_m_n(k,ny,i) - u_m_n(k,ny-1,i)
634
[1353]635                IF ( denom /= 0.0_wp )  THEN
[996]636                   c_u(k,i) = -c_max * ( u(k,ny,i) - u_m_n(k,ny,i) ) / ( denom * tsc(2) )
[1353]637                   IF ( c_u(k,i) < 0.0_wp )  THEN
638                      c_u(k,i) = 0.0_wp
[106]639                   ELSEIF ( c_u(k,i) > c_max )  THEN
640                      c_u(k,i) = c_max
641                   ENDIF
642                ELSE
643                   c_u(k,i) = c_max
[73]644                ENDIF
645
[106]646                denom = v_m_n(k,ny,i) - v_m_n(k,ny-1,i)
[73]647
[1353]648                IF ( denom /= 0.0_wp )  THEN
[996]649                   c_v(k,i) = -c_max * ( v(k,ny,i) - v_m_n(k,ny,i) ) / ( denom * tsc(2) )
[1353]650                   IF ( c_v(k,i) < 0.0_wp )  THEN
651                      c_v(k,i) = 0.0_wp
[106]652                   ELSEIF ( c_v(k,i) > c_max )  THEN
653                      c_v(k,i) = c_max
654                   ENDIF
655                ELSE
656                   c_v(k,i) = c_max
[73]657                ENDIF
658
[106]659                denom = w_m_n(k,ny,i) - w_m_n(k,ny-1,i)
[73]660
[1353]661                IF ( denom /= 0.0_wp )  THEN
[996]662                   c_w(k,i) = -c_max * ( w(k,ny,i) - w_m_n(k,ny,i) ) / ( denom * tsc(2) )
[1353]663                   IF ( c_w(k,i) < 0.0_wp )  THEN
664                      c_w(k,i) = 0.0_wp
[106]665                   ELSEIF ( c_w(k,i) > c_max )  THEN
666                      c_w(k,i) = c_max
667                   ENDIF
668                ELSE
669                   c_w(k,i) = c_max
[73]670                ENDIF
[106]671
[978]672                c_u_m_l(k) = c_u_m_l(k) + c_u(k,i)
673                c_v_m_l(k) = c_v_m_l(k) + c_v(k,i)
674                c_w_m_l(k) = c_w_m_l(k) + c_w(k,i)
[106]675
[978]676             ENDDO
677          ENDDO
[73]678
[978]679#if defined( __parallel )   
680          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
681          CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
682                              MPI_SUM, comm1dx, ierr )   
683          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
684          CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
685                              MPI_SUM, comm1dx, ierr ) 
686          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
687          CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
688                              MPI_SUM, comm1dx, ierr ) 
689#else
690          c_u_m = c_u_m_l
691          c_v_m = c_v_m_l
692          c_w_m = c_w_m_l
693#endif
694
695          c_u_m = c_u_m / (nx+1)
696          c_v_m = c_v_m / (nx+1)
697          c_w_m = c_w_m / (nx+1)
698
[73]699!
[978]700!--       Save old timelevels for the next timestep
701          IF ( intermediate_timestep_count == 1 )  THEN
702                u_m_n(:,:,:) = u(:,ny-1:ny,:)
703                v_m_n(:,:,:) = v(:,ny-1:ny,:)
704                w_m_n(:,:,:) = w(:,ny-1:ny,:)
705          ENDIF
[73]706
[978]707!
708!--       Calculate the new velocities
[996]709          DO  k = nzb+1, nzt+1
710             DO  i = nxlg, nxrg
[978]711                u_p(k,ny+1,i) = u(k,ny+1,i) - dt_3d * tsc(2) * c_u_m(k) *      &
712                                       ( u(k,ny+1,i) - u(k,ny,i) ) * ddy
[73]713
[978]714                v_p(k,ny+1,i) = v(k,ny+1,i)  - dt_3d * tsc(2) * c_v_m(k) *     &
715                                       ( v(k,ny+1,i) - v(k,ny,i) ) * ddy
[73]716
[978]717                w_p(k,ny+1,i) = w(k,ny+1,i) - dt_3d * tsc(2) * c_w_m(k) *      &
718                                       ( w(k,ny+1,i) - w(k,ny,i) ) * ddy
719             ENDDO
[1]720          ENDDO
721
722!
[978]723!--       Bottom boundary at the outflow
724          IF ( ibc_uv_b == 0 )  THEN
[1353]725             u_p(nzb,ny+1,:) = 0.0_wp
726             v_p(nzb,ny+1,:) = 0.0_wp   
[978]727          ELSE                   
728             u_p(nzb,ny+1,:) =  u_p(nzb+1,ny+1,:)
729             v_p(nzb,ny+1,:) =  v_p(nzb+1,ny+1,:)
730          ENDIF
[1353]731          w_p(nzb,ny+1,:) = 0.0_wp
[73]732
733!
[978]734!--       Top boundary at the outflow
735          IF ( ibc_uv_t == 0 )  THEN
736             u_p(nzt+1,ny+1,:) = u_init(nzt+1)
737             v_p(nzt+1,ny+1,:) = v_init(nzt+1)
738          ELSE
739             u_p(nzt+1,ny+1,:) = u_p(nzt,nyn+1,:)
740             v_p(nzt+1,ny+1,:) = v_p(nzt,nyn+1,:)
741          ENDIF
[1353]742          w_p(nzt:nzt+1,ny+1,:) = 0.0_wp
[978]743
[1]744       ENDIF
745
[75]746    ENDIF
747
[106]748    IF ( outflow_l )  THEN
[75]749
[1159]750       IF ( use_cmax )  THEN
[1717]751          u_p(:,:,0)  = u(:,:,1)
752          v_p(:,:,-1) = v(:,:,0)
[1159]753          w_p(:,:,-1) = w(:,:,0)         
754       ELSEIF ( .NOT. use_cmax )  THEN
[75]755
[978]756          c_max = dx / dt_3d
[75]757
[1353]758          c_u_m_l = 0.0_wp 
759          c_v_m_l = 0.0_wp
760          c_w_m_l = 0.0_wp
[978]761
[1353]762          c_u_m = 0.0_wp 
763          c_v_m = 0.0_wp
764          c_w_m = 0.0_wp
[978]765
[1]766!
[996]767!--       Calculate the phase speeds for u, v, and w, first local and then
768!--       average along the outflow boundary.
769          DO  k = nzb+1, nzt+1
770             DO  j = nys, nyn
[75]771
[106]772                denom = u_m_l(k,j,1) - u_m_l(k,j,2)
773
[1353]774                IF ( denom /= 0.0_wp )  THEN
[996]775                   c_u(k,j) = -c_max * ( u(k,j,1) - u_m_l(k,j,1) ) / ( denom * tsc(2) )
[1353]776                   IF ( c_u(k,j) < 0.0_wp )  THEN
777                      c_u(k,j) = 0.0_wp
[107]778                   ELSEIF ( c_u(k,j) > c_max )  THEN
779                      c_u(k,j) = c_max
[106]780                   ENDIF
781                ELSE
[107]782                   c_u(k,j) = c_max
[75]783                ENDIF
784
[106]785                denom = v_m_l(k,j,0) - v_m_l(k,j,1)
[75]786
[1353]787                IF ( denom /= 0.0_wp )  THEN
[996]788                   c_v(k,j) = -c_max * ( v(k,j,0) - v_m_l(k,j,0) ) / ( denom * tsc(2) )
[1353]789                   IF ( c_v(k,j) < 0.0_wp )  THEN
790                      c_v(k,j) = 0.0_wp
[106]791                   ELSEIF ( c_v(k,j) > c_max )  THEN
792                      c_v(k,j) = c_max
793                   ENDIF
794                ELSE
795                   c_v(k,j) = c_max
[75]796                ENDIF
797
[106]798                denom = w_m_l(k,j,0) - w_m_l(k,j,1)
[75]799
[1353]800                IF ( denom /= 0.0_wp )  THEN
[996]801                   c_w(k,j) = -c_max * ( w(k,j,0) - w_m_l(k,j,0) ) / ( denom * tsc(2) )
[1353]802                   IF ( c_w(k,j) < 0.0_wp )  THEN
803                      c_w(k,j) = 0.0_wp
[106]804                   ELSEIF ( c_w(k,j) > c_max )  THEN
805                      c_w(k,j) = c_max
806                   ENDIF
807                ELSE
808                   c_w(k,j) = c_max
[75]809                ENDIF
[106]810
[978]811                c_u_m_l(k) = c_u_m_l(k) + c_u(k,j)
812                c_v_m_l(k) = c_v_m_l(k) + c_v(k,j)
813                c_w_m_l(k) = c_w_m_l(k) + c_w(k,j)
[106]814
[978]815             ENDDO
816          ENDDO
[75]817
[978]818#if defined( __parallel )   
819          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
820          CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
821                              MPI_SUM, comm1dy, ierr )   
822          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
823          CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
824                              MPI_SUM, comm1dy, ierr ) 
825          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
826          CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
827                              MPI_SUM, comm1dy, ierr ) 
828#else
829          c_u_m = c_u_m_l
830          c_v_m = c_v_m_l
831          c_w_m = c_w_m_l
832#endif
833
834          c_u_m = c_u_m / (ny+1)
835          c_v_m = c_v_m / (ny+1)
836          c_w_m = c_w_m / (ny+1)
837
[73]838!
[978]839!--       Save old timelevels for the next timestep
840          IF ( intermediate_timestep_count == 1 )  THEN
841                u_m_l(:,:,:) = u(:,:,1:2)
842                v_m_l(:,:,:) = v(:,:,0:1)
843                w_m_l(:,:,:) = w(:,:,0:1)
844          ENDIF
845
846!
847!--       Calculate the new velocities
[996]848          DO  k = nzb+1, nzt+1
[1113]849             DO  j = nysg, nyng
[978]850                u_p(k,j,0) = u(k,j,0) - dt_3d * tsc(2) * c_u_m(k) *            &
[106]851                                       ( u(k,j,0) - u(k,j,1) ) * ddx
[75]852
[978]853                v_p(k,j,-1) = v(k,j,-1) - dt_3d * tsc(2) * c_v_m(k) *          &
[75]854                                       ( v(k,j,-1) - v(k,j,0) ) * ddx
855
[978]856                w_p(k,j,-1) = w(k,j,-1) - dt_3d * tsc(2) * c_w_m(k) *          &
[75]857                                       ( w(k,j,-1) - w(k,j,0) ) * ddx
[978]858             ENDDO
[75]859          ENDDO
860
861!
[978]862!--       Bottom boundary at the outflow
863          IF ( ibc_uv_b == 0 )  THEN
[1353]864             u_p(nzb,:,0)  = 0.0_wp 
865             v_p(nzb,:,-1) = 0.0_wp
[978]866          ELSE                   
867             u_p(nzb,:,0)  =  u_p(nzb+1,:,0)
868             v_p(nzb,:,-1) =  v_p(nzb+1,:,-1)
869          ENDIF
[1353]870          w_p(nzb,:,-1) = 0.0_wp
[1]871
[75]872!
[978]873!--       Top boundary at the outflow
874          IF ( ibc_uv_t == 0 )  THEN
[1764]875             u_p(nzt+1,:,0)  = u_init(nzt+1)
[978]876             v_p(nzt+1,:,-1) = v_init(nzt+1)
877          ELSE
[1764]878             u_p(nzt+1,:,0)  = u_p(nzt,:,0)
[978]879             v_p(nzt+1,:,-1) = v_p(nzt,:,-1)
880          ENDIF
[1353]881          w_p(nzt:nzt+1,:,-1) = 0.0_wp
[978]882
[75]883       ENDIF
[73]884
[75]885    ENDIF
[73]886
[106]887    IF ( outflow_r )  THEN
[73]888
[1159]889       IF ( use_cmax )  THEN
890          u_p(:,:,nx+1) = u(:,:,nx)
891          v_p(:,:,nx+1) = v(:,:,nx)
892          w_p(:,:,nx+1) = w(:,:,nx)         
893       ELSEIF ( .NOT. use_cmax )  THEN
[75]894
[978]895          c_max = dx / dt_3d
[75]896
[1353]897          c_u_m_l = 0.0_wp 
898          c_v_m_l = 0.0_wp
899          c_w_m_l = 0.0_wp
[978]900
[1353]901          c_u_m = 0.0_wp 
902          c_v_m = 0.0_wp
903          c_w_m = 0.0_wp
[978]904
[1]905!
[996]906!--       Calculate the phase speeds for u, v, and w, first local and then
907!--       average along the outflow boundary.
908          DO  k = nzb+1, nzt+1
909             DO  j = nys, nyn
[73]910
[106]911                denom = u_m_r(k,j,nx) - u_m_r(k,j,nx-1)
912
[1353]913                IF ( denom /= 0.0_wp )  THEN
[996]914                   c_u(k,j) = -c_max * ( u(k,j,nx) - u_m_r(k,j,nx) ) / ( denom * tsc(2) )
[1353]915                   IF ( c_u(k,j) < 0.0_wp )  THEN
916                      c_u(k,j) = 0.0_wp
[106]917                   ELSEIF ( c_u(k,j) > c_max )  THEN
918                      c_u(k,j) = c_max
919                   ENDIF
920                ELSE
921                   c_u(k,j) = c_max
[73]922                ENDIF
923
[106]924                denom = v_m_r(k,j,nx) - v_m_r(k,j,nx-1)
[73]925
[1353]926                IF ( denom /= 0.0_wp )  THEN
[996]927                   c_v(k,j) = -c_max * ( v(k,j,nx) - v_m_r(k,j,nx) ) / ( denom * tsc(2) )
[1353]928                   IF ( c_v(k,j) < 0.0_wp )  THEN
929                      c_v(k,j) = 0.0_wp
[106]930                   ELSEIF ( c_v(k,j) > c_max )  THEN
931                      c_v(k,j) = c_max
932                   ENDIF
933                ELSE
934                   c_v(k,j) = c_max
[73]935                ENDIF
936
[106]937                denom = w_m_r(k,j,nx) - w_m_r(k,j,nx-1)
[73]938
[1353]939                IF ( denom /= 0.0_wp )  THEN
[996]940                   c_w(k,j) = -c_max * ( w(k,j,nx) - w_m_r(k,j,nx) ) / ( denom * tsc(2) )
[1353]941                   IF ( c_w(k,j) < 0.0_wp )  THEN
942                      c_w(k,j) = 0.0_wp
[106]943                   ELSEIF ( c_w(k,j) > c_max )  THEN
944                      c_w(k,j) = c_max
945                   ENDIF
946                ELSE
947                   c_w(k,j) = c_max
[73]948                ENDIF
[106]949
[978]950                c_u_m_l(k) = c_u_m_l(k) + c_u(k,j)
951                c_v_m_l(k) = c_v_m_l(k) + c_v(k,j)
952                c_w_m_l(k) = c_w_m_l(k) + c_w(k,j)
[106]953
[978]954             ENDDO
955          ENDDO
[73]956
[978]957#if defined( __parallel )   
958          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
959          CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
960                              MPI_SUM, comm1dy, ierr )   
961          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
962          CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_m(nzb+1), nzt-nzb, MPI_REAL, &
963                              MPI_SUM, comm1dy, ierr ) 
964          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
965          CALL MPI_ALLREDUCE( c_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
966                              MPI_SUM, comm1dy, ierr ) 
967#else
968          c_u_m = c_u_m_l
969          c_v_m = c_v_m_l
970          c_w_m = c_w_m_l
971#endif
972
973          c_u_m = c_u_m / (ny+1)
974          c_v_m = c_v_m / (ny+1)
975          c_w_m = c_w_m / (ny+1)
976
[73]977!
[978]978!--       Save old timelevels for the next timestep
979          IF ( intermediate_timestep_count == 1 )  THEN
980                u_m_r(:,:,:) = u(:,:,nx-1:nx)
981                v_m_r(:,:,:) = v(:,:,nx-1:nx)
982                w_m_r(:,:,:) = w(:,:,nx-1:nx)
983          ENDIF
[73]984
[978]985!
986!--       Calculate the new velocities
[996]987          DO  k = nzb+1, nzt+1
[1113]988             DO  j = nysg, nyng
[978]989                u_p(k,j,nx+1) = u(k,j,nx+1) - dt_3d * tsc(2) * c_u_m(k) *      &
990                                       ( u(k,j,nx+1) - u(k,j,nx) ) * ddx
[73]991
[978]992                v_p(k,j,nx+1) = v(k,j,nx+1) - dt_3d * tsc(2) * c_v_m(k) *      &
993                                       ( v(k,j,nx+1) - v(k,j,nx) ) * ddx
[73]994
[978]995                w_p(k,j,nx+1) = w(k,j,nx+1) - dt_3d * tsc(2) * c_w_m(k) *      &
996                                       ( w(k,j,nx+1) - w(k,j,nx) ) * ddx
997             ENDDO
[73]998          ENDDO
999
1000!
[978]1001!--       Bottom boundary at the outflow
1002          IF ( ibc_uv_b == 0 )  THEN
[1353]1003             u_p(nzb,:,nx+1) = 0.0_wp
1004             v_p(nzb,:,nx+1) = 0.0_wp 
[978]1005          ELSE                   
1006             u_p(nzb,:,nx+1) =  u_p(nzb+1,:,nx+1)
1007             v_p(nzb,:,nx+1) =  v_p(nzb+1,:,nx+1)
1008          ENDIF
[1353]1009          w_p(nzb,:,nx+1) = 0.0_wp
[73]1010
1011!
[978]1012!--       Top boundary at the outflow
1013          IF ( ibc_uv_t == 0 )  THEN
1014             u_p(nzt+1,:,nx+1) = u_init(nzt+1)
1015             v_p(nzt+1,:,nx+1) = v_init(nzt+1)
1016          ELSE
1017             u_p(nzt+1,:,nx+1) = u_p(nzt,:,nx+1)
1018             v_p(nzt+1,:,nx+1) = v_p(nzt,:,nx+1)
1019          ENDIF
[1742]1020          w_p(nzt:nzt+1,:,nx+1) = 0.0_wp
[978]1021
[1]1022       ENDIF
1023
1024    ENDIF
1025
1026 END SUBROUTINE boundary_conds
Note: See TracBrowser for help on using the repository browser.