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

Last change on this file since 1744 was 1744, checked in by raasch, 6 years ago

two last commits documented

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