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

Last change on this file since 1417 was 1410, checked in by suehring, 10 years ago

last commit documented

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