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

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

code has been put under the GNU General Public License (v3)

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