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

Last change on this file since 1207 was 1160, checked in by fricke, 11 years ago

last commit documented

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