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

Last change on this file since 2128 was 2119, checked in by raasch, 8 years ago

last commit documented

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