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

Last change on this file since 1140 was 1116, checked in by hoffmann, 11 years ago

last commit documented

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