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

Last change on this file since 2065 was 2001, checked in by knoop, 8 years ago

last commit documented

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