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

Last change on this file since 1972 was 1961, checked in by suehring, 8 years ago

last commit documented

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