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

Last change on this file since 1361 was 1361, checked in by hoffmann, 10 years ago

improved version of two-moment cloud physics

  • Property svn:keywords set to Id
File size: 30.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-2014 Leibniz Universitaet Hannover
18!--------------------------------------------------------------------------------!
19!
20! Current revisions:
21! -----------------
22! Bottom and top boundary conditions of rain water content (qr) and
23! rain drop concentration (nr) changed to Dirichlet
24!
25! Former revisions:
26! -----------------
27! $Id: boundary_conds.f90 1361 2014-04-16 15:17:48Z hoffmann $
28!
29! 1353 2014-04-08 15:21:23Z heinze
30! REAL constants provided with KIND-attribute
31
32! 1320 2014-03-20 08:40:49Z raasch
33! ONLY-attribute added to USE-statements,
34! kind-parameters added to all INTEGER and REAL declaration statements,
35! kinds are defined in new module kinds,
36! revision history before 2012 removed,
37! comment fields (!:) to be used for variable explanations added to
38! all variable declaration statements
39!
40! 1257 2013-11-08 15:18:40Z raasch
41! loop independent clauses added
42!
43! 1241 2013-10-30 11:36:58Z heinze
44! Adjust ug and vg at each timestep in case of large_scale_forcing
45!
46! 1159 2013-05-21 11:58:22Z fricke
47! Bugfix: Neumann boundary conditions for the velocity components at the
48! outflow are in fact radiation boundary conditions using the maximum phase
49! velocity that ensures numerical stability (CFL-condition).
50! Hence, logical operator use_cmax is now used instead of bc_lr_dirneu/_neudir.
51! Bugfix: In case of use_cmax at the outflow, u, v, w are replaced by
52! u_p, v_p, w_p 
53!
54! 1115 2013-03-26 18:16:16Z hoffmann
55! boundary conditions of two-moment cloud scheme are restricted to Neumann-
56! boundary-conditions
57!
58! 1113 2013-03-10 02:48:14Z raasch
59! GPU-porting
60! dummy argument "range" removed
61! Bugfix: wrong index in loops of radiation boundary condition
62!
63! 1053 2012-11-13 17:11:03Z hoffmann
64! boundary conditions for the two new prognostic equations (nr, qr) of the
65! two-moment cloud scheme
66!
67! 1036 2012-10-22 13:43:42Z raasch
68! code put under GPL (PALM 3.9)
69!
70! 996 2012-09-07 10:41:47Z raasch
71! little reformatting
72!
73! 978 2012-08-09 08:28:32Z fricke
74! Neumann boudnary conditions are added at the inflow boundary for the SGS-TKE.
75! Outflow boundary conditions for the velocity components can be set to Neumann
76! conditions or to radiation conditions with a horizontal averaged phase
77! velocity.
78!
79! 875 2012-04-02 15:35:15Z gryschka
80! Bugfix in case of dirichlet inflow bc at the right or north boundary
81!
82! Revision 1.1  1997/09/12 06:21:34  raasch
83! Initial revision
84!
85!
86! Description:
87! ------------
88! Boundary conditions for the prognostic quantities.
89! One additional bottom boundary condition is applied for the TKE (=(u*)**2)
90! in prandtl_fluxes. The cyclic lateral boundary conditions are implicitly
91! handled in routine exchange_horiz. Pressure boundary conditions are
92! explicitly set in routines pres, poisfft, poismg and sor.
93!------------------------------------------------------------------------------!
94
95    USE arrays_3d,                                                             &
96        ONLY:  c_u, c_u_m, c_u_m_l, c_v, c_v_m, c_v_m_l, c_w, c_w_m, c_w_m_l,  &
97               dzu, e_p, nr_p, pt, pt_p, q, q_p, qr_p, sa, sa_p,               &
98               u, ug, u_init, u_m_l, u_m_n, u_m_r, u_m_s, u_p,                 &
99               v, vg, v_init, v_m_l, v_m_n, v_m_r, v_m_s, v_p,                 &
100               w, w_p, w_m_l, w_m_n, w_m_r, w_m_s
101
102    USE control_parameters,                                                    &
103        ONLY:  bc_pt_t_val, bc_q_t_val, constant_diffusion,                    &
104               cloud_physics, dt_3d, humidity,                                 &
105               ibc_pt_b, ibc_pt_t, ibc_q_b, ibc_sa_t, ibc_uv_b, ibc_uv_t,      &
106               icloud_scheme, inflow_l, inflow_n, inflow_r, inflow_s,          &
107               intermediate_timestep_count, large_scale_forcing, ocean,        &
108               outflow_l, outflow_n, outflow_r, outflow_s, passive_scalar,     &
109               precipitation, tsc, use_cmax
110
111    USE grid_variables,                                                        &
112        ONLY:  ddx, ddy, dx, dy
113
114    USE indices,                                                               &
115        ONLY:  nx, nxl, nxlg, nxr, nxrg, ny, nyn, nyng, nys, nysg,             &
116               nzb, nzb_s_inner, nzb_w_inner, nzt
117
118    USE kinds
119
120    USE pegrid
121
122
123    IMPLICIT NONE
124
125    INTEGER(iwp) ::  i !:
126    INTEGER(iwp) ::  j !:
127    INTEGER(iwp) ::  k !:
128
129    REAL(wp)    ::  c_max !:
130    REAL(wp)    ::  denom !:
131
132
133!
134!-- Bottom boundary
135    IF ( ibc_uv_b == 1 )  THEN
136       !$acc kernels present( u_p, v_p )
137       u_p(nzb,:,:) = u_p(nzb+1,:,:)
138       v_p(nzb,:,:) = v_p(nzb+1,:,:)
139       !$acc end kernels
140    ENDIF
141
142    !$acc kernels present( nzb_w_inner, w_p )
143    DO  i = nxlg, nxrg
144       DO  j = nysg, nyng
145          w_p(nzb_w_inner(j,i),j,i) = 0.0_wp
146       ENDDO
147    ENDDO
148    !$acc end kernels
149
150!
151!-- Top boundary
152    IF ( ibc_uv_t == 0 )  THEN
153       !$acc kernels present( u_init, u_p, v_init, v_p )
154        u_p(nzt+1,:,:) = u_init(nzt+1)
155        v_p(nzt+1,:,:) = v_init(nzt+1)
156        IF ( large_scale_forcing) THEN
157           u_p(nzt+1,:,:) = ug(nzt+1)
158           v_p(nzt+1,:,:) = vg(nzt+1)
159        END IF
160       !$acc end kernels
161    ELSE
162       !$acc kernels present( u_p, v_p )
163        u_p(nzt+1,:,:) = u_p(nzt,:,:)
164        v_p(nzt+1,:,:) = v_p(nzt,:,:)
165       !$acc end kernels
166    ENDIF
167    !$acc kernels present( w_p )
168    w_p(nzt:nzt+1,:,:) = 0.0_wp  ! nzt is not a prognostic level (but cf. pres)
169    !$acc end kernels
170
171!
172!-- Temperature at bottom boundary.
173!-- In case of coupled runs (ibc_pt_b = 2) the temperature is given by
174!-- the sea surface temperature of the coupled ocean model.
175    IF ( ibc_pt_b == 0 )  THEN
176       !$acc kernels present( nzb_s_inner, pt, pt_p )
177       !$acc loop independent
178       DO  i = nxlg, nxrg
179          !$acc loop independent
180          DO  j = nysg, nyng
181             pt_p(nzb_s_inner(j,i),j,i) = pt(nzb_s_inner(j,i),j,i)
182          ENDDO
183       ENDDO
184       !$acc end kernels
185    ELSEIF ( ibc_pt_b == 1 )  THEN
186       !$acc kernels present( nzb_s_inner, pt_p )
187       !$acc loop independent
188       DO  i = nxlg, nxrg
189          !$acc loop independent
190          DO  j = nysg, nyng
191             pt_p(nzb_s_inner(j,i),j,i) = pt_p(nzb_s_inner(j,i)+1,j,i)
192          ENDDO
193       ENDDO
194      !$acc end kernels
195    ENDIF
196
197!
198!-- Temperature at top boundary
199    IF ( ibc_pt_t == 0 )  THEN
200       !$acc kernels present( pt, pt_p )
201        pt_p(nzt+1,:,:) = pt(nzt+1,:,:)
202       !$acc end kernels
203    ELSEIF ( ibc_pt_t == 1 )  THEN
204       !$acc kernels present( pt_p )
205        pt_p(nzt+1,:,:) = pt_p(nzt,:,:)
206       !$acc end kernels
207    ELSEIF ( ibc_pt_t == 2 )  THEN
208       !$acc kernels present( dzu, pt_p )
209        pt_p(nzt+1,:,:) = pt_p(nzt,:,:)   + bc_pt_t_val * dzu(nzt+1)
210       !$acc end kernels
211    ENDIF
212
213!
214!-- Boundary conditions for TKE
215!-- Generally Neumann conditions with de/dz=0 are assumed
216    IF ( .NOT. constant_diffusion )  THEN
217       !$acc kernels present( e_p, nzb_s_inner )
218       !$acc loop independent
219       DO  i = nxlg, nxrg
220          !$acc loop independent
221          DO  j = nysg, nyng
222             e_p(nzb_s_inner(j,i),j,i) = e_p(nzb_s_inner(j,i)+1,j,i)
223          ENDDO
224       ENDDO
225       e_p(nzt+1,:,:) = e_p(nzt,:,:)
226       !$acc end kernels
227    ENDIF
228
229!
230!-- Boundary conditions for salinity
231    IF ( ocean )  THEN
232!
233!--    Bottom boundary: Neumann condition because salinity flux is always
234!--    given
235       DO  i = nxlg, nxrg
236          DO  j = nysg, nyng
237             sa_p(nzb_s_inner(j,i),j,i) = sa_p(nzb_s_inner(j,i)+1,j,i)
238          ENDDO
239       ENDDO
240
241!
242!--    Top boundary: Dirichlet or Neumann
243       IF ( ibc_sa_t == 0 )  THEN
244           sa_p(nzt+1,:,:) = sa(nzt+1,:,:)
245       ELSEIF ( ibc_sa_t == 1 )  THEN
246           sa_p(nzt+1,:,:) = sa_p(nzt,:,:)
247       ENDIF
248
249    ENDIF
250
251!
252!-- Boundary conditions for total water content or scalar,
253!-- bottom and top boundary (see also temperature)
254    IF ( humidity  .OR.  passive_scalar )  THEN
255!
256!--    Surface conditions for constant_humidity_flux
257       IF ( ibc_q_b == 0 ) THEN
258          DO  i = nxlg, nxrg
259             DO  j = nysg, nyng
260                q_p(nzb_s_inner(j,i),j,i) = q(nzb_s_inner(j,i),j,i)
261             ENDDO
262          ENDDO
263       ELSE
264          DO  i = nxlg, nxrg
265             DO  j = nysg, nyng
266                q_p(nzb_s_inner(j,i),j,i) = q_p(nzb_s_inner(j,i)+1,j,i)
267             ENDDO
268          ENDDO
269       ENDIF
270!
271!--    Top boundary
272       q_p(nzt+1,:,:) = q_p(nzt,:,:)   + bc_q_t_val * dzu(nzt+1)
273
274       IF ( cloud_physics  .AND.  icloud_scheme == 0  .AND.  precipitation )  THEN
275!             
276!--       Surface conditions rain water (Dirichlet)
277          DO  i = nxlg, nxrg
278             DO  j = nysg, nyng
279                qr_p(nzb_s_inner(j,i),j,i) = 0.0_wp
280                nr_p(nzb_s_inner(j,i),j,i) = 0.0_wp
281             ENDDO
282          ENDDO
283!
284!--       Top boundary condition for rain water (Dirichlet)
285          qr_p(nzt+1,:,:) = 0.0_wp
286          nr_p(nzt+1,:,:) = 0.0_wp
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_wp 
372          c_v_m_l = 0.0_wp
373          c_w_m_l = 0.0_wp
374
375          c_u_m = 0.0_wp 
376          c_v_m = 0.0_wp
377          c_w_m = 0.0_wp
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_wp )  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_wp )  THEN
390                      c_u(k,i) = 0.0_wp
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_wp )  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_wp )  THEN
403                      c_v(k,i) = 0.0_wp
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_wp )  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_wp )  THEN
416                      c_w(k,i) = 0.0_wp
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_wp 
478             v_p(nzb,0,:)  = 0.0_wp 
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_wp
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_wp
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_wp 
511          c_v_m_l = 0.0_wp
512          c_w_m_l = 0.0_wp
513
514          c_u_m = 0.0_wp 
515          c_v_m = 0.0_wp
516          c_w_m = 0.0_wp
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_wp )  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_wp )  THEN
529                      c_u(k,i) = 0.0_wp
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_wp )  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_wp )  THEN
542                      c_v(k,i) = 0.0_wp
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_wp )  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_wp )  THEN
555                      c_w(k,i) = 0.0_wp
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_wp
617             v_p(nzb,ny+1,:) = 0.0_wp   
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_wp
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_wp
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_wp 
650          c_v_m_l = 0.0_wp
651          c_w_m_l = 0.0_wp
652
653          c_u_m = 0.0_wp 
654          c_v_m = 0.0_wp
655          c_w_m = 0.0_wp
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_wp )  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_wp )  THEN
668                      c_u(k,j) = 0.0_wp
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_wp )  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_wp )  THEN
681                      c_v(k,j) = 0.0_wp
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_wp )  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_wp )  THEN
694                      c_w(k,j) = 0.0_wp
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_wp 
756             v_p(nzb,:,-1) = 0.0_wp
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_wp
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_wp
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_wp 
789          c_v_m_l = 0.0_wp
790          c_w_m_l = 0.0_wp
791
792          c_u_m = 0.0_wp 
793          c_v_m = 0.0_wp
794          c_w_m = 0.0_wp
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_wp )  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_wp )  THEN
807                      c_u(k,j) = 0.0_wp
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_wp )  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_wp )  THEN
820                      c_v(k,j) = 0.0_wp
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_wp )  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_wp )  THEN
833                      c_w(k,j) = 0.0_wp
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_wp
895             v_p(nzb,:,nx+1) = 0.0_wp 
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_wp
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_wp
912
913       ENDIF
914
915    ENDIF
916
917 END SUBROUTINE boundary_conds
Note: See TracBrowser for help on using the repository browser.