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

Last change on this file since 1734 was 1718, checked in by raasch, 9 years ago

last commit documented

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