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

Last change on this file since 1267 was 1258, checked in by raasch, 11 years ago

last commit documented

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