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

Last change on this file since 1004 was 997, checked in by raasch, 12 years ago

last commit documented

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