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

Last change on this file since 1253 was 1242, checked in by heinze, 11 years ago

Last commmit documented

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