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

Last change on this file since 1037 was 1037, checked in by raasch, 9 years ago

last commit documented

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