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

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

last commit documented

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