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

Last change on this file since 1320 was 1320, checked in by raasch, 7 years ago

ONLY-attribute added to USE-statements,
kind-parameters added to all INTEGER and REAL declaration statements,
kinds are defined in new module kinds,
old module precision_kind is removed,
revision history before 2012 removed,
comment fields (!:) to be used for variable explanations added to all variable declaration statements

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