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

Last change on this file since 2118 was 2118, checked in by raasch, 5 years ago

all OpenACC directives and related parts removed from the code

  • 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! OpenACC directives removed
23!
24! Former revisions:
25! -----------------
26! $Id: boundary_conds.f90 2118 2017-01-17 16:38:49Z raasch $
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       u_p(nzb,:,:) = u_p(nzb+1,:,:)
191       v_p(nzb,:,:) = v_p(nzb+1,:,:)
192    ENDIF
193
194    DO  i = nxlg, nxrg
195       DO  j = nysg, nyng
196          w_p(nzb_w_inner(j,i),j,i) = 0.0_wp
197       ENDDO
198    ENDDO
199
200!
201!-- Top boundary. A nested domain ( ibc_uv_t = 3 ) does not require settings.
202    IF ( ibc_uv_t == 0 )  THEN
203        u_p(nzt+1,:,:) = u_init(nzt+1)
204        v_p(nzt+1,:,:) = v_init(nzt+1)
205    ELSEIF ( ibc_uv_t == 1 )  THEN
206        u_p(nzt+1,:,:) = u_p(nzt,:,:)
207        v_p(nzt+1,:,:) = v_p(nzt,:,:)
208    ENDIF
209
210    IF ( .NOT. nest_domain )  THEN
211       w_p(nzt:nzt+1,:,:) = 0.0_wp  ! nzt is not a prognostic level (but cf. pres)
212    ENDIF
213
214!
215!-- Temperature at bottom boundary.
216!-- In case of coupled runs (ibc_pt_b = 2) the temperature is given by
217!-- the sea surface temperature of the coupled ocean model.
218    IF ( ibc_pt_b == 0 )  THEN
219       DO  i = nxlg, nxrg
220          DO  j = nysg, nyng
221             pt_p(nzb_s_inner(j,i),j,i) = pt(nzb_s_inner(j,i),j,i)
222          ENDDO
223       ENDDO
224    ELSEIF ( ibc_pt_b == 1 )  THEN
225       DO  i = nxlg, nxrg
226          DO  j = nysg, nyng
227             pt_p(nzb_s_inner(j,i),j,i) = pt_p(nzb_s_inner(j,i)+1,j,i)
228          ENDDO
229       ENDDO
230    ENDIF
231
232!
233!-- Temperature at top boundary
234    IF ( ibc_pt_t == 0 )  THEN
235        pt_p(nzt+1,:,:) = pt(nzt+1,:,:)
236!
237!--     In case of nudging adjust top boundary to pt which is
238!--     read in from NUDGING-DATA
239        IF ( nudging )  THEN
240           pt_p(nzt+1,:,:) = pt_init(nzt+1)
241        ENDIF
242    ELSEIF ( ibc_pt_t == 1 )  THEN
243        pt_p(nzt+1,:,:) = pt_p(nzt,:,:)
244    ELSEIF ( ibc_pt_t == 2 )  THEN
245        pt_p(nzt+1,:,:) = pt_p(nzt,:,:) + bc_pt_t_val * dzu(nzt+1)
246    ENDIF
247
248!
249!-- Boundary conditions for TKE
250!-- Generally Neumann conditions with de/dz=0 are assumed
251    IF ( .NOT. constant_diffusion )  THEN
252       DO  i = nxlg, nxrg
253          DO  j = nysg, nyng
254             e_p(nzb_s_inner(j,i),j,i) = e_p(nzb_s_inner(j,i)+1,j,i)
255          ENDDO
256       ENDDO
257       IF ( .NOT. nest_domain )  THEN
258          e_p(nzt+1,:,:) = e_p(nzt,:,:)
259       ENDIF
260    ENDIF
261
262!
263!-- Boundary conditions for salinity
264    IF ( ocean )  THEN
265!
266!--    Bottom boundary: Neumann condition because salinity flux is always
267!--    given
268       DO  i = nxlg, nxrg
269          DO  j = nysg, nyng
270             sa_p(nzb_s_inner(j,i),j,i) = sa_p(nzb_s_inner(j,i)+1,j,i)
271          ENDDO
272       ENDDO
273
274!
275!--    Top boundary: Dirichlet or Neumann
276       IF ( ibc_sa_t == 0 )  THEN
277           sa_p(nzt+1,:,:) = sa(nzt+1,:,:)
278       ELSEIF ( ibc_sa_t == 1 )  THEN
279           sa_p(nzt+1,:,:) = sa_p(nzt,:,:)
280       ENDIF
281
282    ENDIF
283
284!
285!-- Boundary conditions for total water content,
286!-- bottom and top boundary (see also temperature)
287    IF ( humidity )  THEN
288!
289!--    Surface conditions for constant_humidity_flux
290       IF ( ibc_q_b == 0 ) THEN
291          DO  i = nxlg, nxrg
292             DO  j = nysg, nyng
293                q_p(nzb_s_inner(j,i),j,i) = q(nzb_s_inner(j,i),j,i)
294             ENDDO
295          ENDDO
296       ELSE
297          DO  i = nxlg, nxrg
298             DO  j = nysg, nyng
299                q_p(nzb_s_inner(j,i),j,i) = q_p(nzb_s_inner(j,i)+1,j,i)
300             ENDDO
301          ENDDO
302       ENDIF
303!
304!--    Top boundary
305       IF ( ibc_q_t == 0 ) THEN
306          q_p(nzt+1,:,:) = q(nzt+1,:,:)
307       ELSEIF ( ibc_q_t == 1 ) THEN
308          q_p(nzt+1,:,:) = q_p(nzt,:,:) + bc_q_t_val * dzu(nzt+1)
309       ENDIF
310
311       IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
312!             
313!--       Surface conditions rain water (Dirichlet)
314          DO  i = nxlg, nxrg
315             DO  j = nysg, nyng
316                qr_p(nzb_s_inner(j,i),j,i) = 0.0_wp
317                nr_p(nzb_s_inner(j,i),j,i) = 0.0_wp
318             ENDDO
319          ENDDO
320!
321!--       Top boundary condition for rain water (Dirichlet)
322          qr_p(nzt+1,:,:) = 0.0_wp
323          nr_p(nzt+1,:,:) = 0.0_wp
324           
325       ENDIF
326    ENDIF
327!
328!-- Boundary conditions for scalar,
329!-- bottom and top boundary (see also temperature)
330    IF ( passive_scalar )  THEN
331!
332!--    Surface conditions for constant_humidity_flux
333       IF ( ibc_s_b == 0 ) THEN
334          DO  i = nxlg, nxrg
335             DO  j = nysg, nyng
336                s_p(nzb_s_inner(j,i),j,i) = s(nzb_s_inner(j,i),j,i)
337             ENDDO
338          ENDDO
339       ELSE
340          DO  i = nxlg, nxrg
341             DO  j = nysg, nyng
342                s_p(nzb_s_inner(j,i),j,i) = s_p(nzb_s_inner(j,i)+1,j,i)
343             ENDDO
344          ENDDO
345       ENDIF
346!
347!--    Top boundary condition
348       IF ( ibc_s_t == 0 )  THEN
349          s_p(nzt+1,:,:) = s(nzt+1,:,:)
350       ELSEIF ( ibc_s_t == 1 )  THEN
351          s_p(nzt+1,:,:) = s_p(nzt,:,:)
352       ELSEIF ( ibc_s_t == 2 )  THEN
353          s_p(nzt+1,:,:) = s_p(nzt,:,:) + bc_s_t_val * dzu(nzt+1)
354       ENDIF
355
356    ENDIF   
357!
358!-- In case of inflow or nest boundary at the south boundary the boundary for v
359!-- is at nys and in case of inflow or nest boundary at the left boundary the
360!-- boundary for u is at nxl. Since in prognostic_equations (cache optimized
361!-- version) these levels are handled as a prognostic level, boundary values
362!-- have to be restored here.
363!-- For the SGS-TKE, Neumann boundary conditions are used at the inflow.
364    IF ( inflow_s )  THEN
365       v_p(:,nys,:) = v_p(:,nys-1,:)
366       IF ( .NOT. constant_diffusion ) e_p(:,nys-1,:) = e_p(:,nys,:)
367    ELSEIF ( inflow_n )  THEN
368       IF ( .NOT. constant_diffusion ) e_p(:,nyn+1,:) = e_p(:,nyn,:)
369    ELSEIF ( inflow_l ) THEN
370       u_p(:,:,nxl) = u_p(:,:,nxl-1)
371       IF ( .NOT. constant_diffusion ) e_p(:,:,nxl-1) = e_p(:,:,nxl)
372    ELSEIF ( inflow_r )  THEN
373       IF ( .NOT. constant_diffusion ) e_p(:,:,nxr+1) = e_p(:,:,nxr)
374    ENDIF
375
376!
377!-- The same restoration for u at i=nxl and v at j=nys as above must be made
378!-- in case of nest boundaries. This must not be done in case of vertical nesting
379!-- mode as in that case the lateral boundaries are actually cyclic.
380    IF ( nesting_mode /= 'vertical' )  THEN
381       IF ( nest_bound_s )  THEN
382          v_p(:,nys,:) = v_p(:,nys-1,:)
383       ENDIF
384       IF ( nest_bound_l )  THEN
385          u_p(:,:,nxl) = u_p(:,:,nxl-1)
386       ENDIF
387    ENDIF
388
389!
390!-- Lateral boundary conditions for scalar quantities at the outflow
391    IF ( outflow_s )  THEN
392       pt_p(:,nys-1,:)     = pt_p(:,nys,:)
393       IF ( .NOT. constant_diffusion     )  e_p(:,nys-1,:) = e_p(:,nys,:)
394       IF ( humidity )  THEN
395          q_p(:,nys-1,:) = q_p(:,nys,:)
396          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
397             qr_p(:,nys-1,:) = qr_p(:,nys,:)
398             nr_p(:,nys-1,:) = nr_p(:,nys,:)
399          ENDIF
400       ENDIF
401       IF ( passive_scalar )  s_p(:,nys-1,:) = s_p(:,nys,:)
402    ELSEIF ( outflow_n )  THEN
403       pt_p(:,nyn+1,:)     = pt_p(:,nyn,:)
404       IF ( .NOT. constant_diffusion     )  e_p(:,nyn+1,:) = e_p(:,nyn,:)
405       IF ( humidity )  THEN
406          q_p(:,nyn+1,:) = q_p(:,nyn,:)
407          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
408             qr_p(:,nyn+1,:) = qr_p(:,nyn,:)
409             nr_p(:,nyn+1,:) = nr_p(:,nyn,:)
410          ENDIF
411       ENDIF
412       IF ( passive_scalar )  s_p(:,nyn+1,:) = s_p(:,nyn,:)
413    ELSEIF ( outflow_l )  THEN
414       pt_p(:,:,nxl-1)     = pt_p(:,:,nxl)
415       IF ( .NOT. constant_diffusion     )  e_p(:,:,nxl-1) = e_p(:,:,nxl)
416       IF ( humidity )  THEN
417          q_p(:,:,nxl-1) = q_p(:,:,nxl)
418          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
419             qr_p(:,:,nxl-1) = qr_p(:,:,nxl)
420             nr_p(:,:,nxl-1) = nr_p(:,:,nxl)
421          ENDIF
422       ENDIF
423       IF ( passive_scalar )  s_p(:,:,nxl-1) = s_p(:,:,nxl)
424    ELSEIF ( outflow_r )  THEN
425       pt_p(:,:,nxr+1)     = pt_p(:,:,nxr)
426       IF ( .NOT. constant_diffusion     )  e_p(:,:,nxr+1) = e_p(:,:,nxr)
427       IF ( humidity )  THEN
428          q_p(:,:,nxr+1) = q_p(:,:,nxr)
429          IF ( cloud_physics  .AND.  microphysics_seifert )  THEN
430             qr_p(:,:,nxr+1) = qr_p(:,:,nxr)
431             nr_p(:,:,nxr+1) = nr_p(:,:,nxr)
432          ENDIF
433       ENDIF
434       IF ( passive_scalar )  s_p(:,:,nxr+1) = s_p(:,:,nxr)
435    ENDIF
436
437!
438!-- Radiation boundary conditions for the velocities at the respective outflow.
439!-- The phase velocity is either assumed to the maximum phase velocity that
440!-- ensures numerical stability (CFL-condition) or calculated after
441!-- Orlanski(1976) and averaged along the outflow boundary.
442    IF ( outflow_s )  THEN
443
444       IF ( use_cmax )  THEN
445          u_p(:,-1,:) = u(:,0,:)
446          v_p(:,0,:)  = v(:,1,:)
447          w_p(:,-1,:) = w(:,0,:)         
448       ELSEIF ( .NOT. use_cmax )  THEN
449
450          c_max = dy / dt_3d
451
452          c_u_m_l = 0.0_wp 
453          c_v_m_l = 0.0_wp
454          c_w_m_l = 0.0_wp
455
456          c_u_m = 0.0_wp 
457          c_v_m = 0.0_wp
458          c_w_m = 0.0_wp
459
460!
461!--       Calculate the phase speeds for u, v, and w, first local and then
462!--       average along the outflow boundary.
463          DO  k = nzb+1, nzt+1
464             DO  i = nxl, nxr
465
466                denom = u_m_s(k,0,i) - u_m_s(k,1,i)
467
468                IF ( denom /= 0.0_wp )  THEN
469                   c_u(k,i) = -c_max * ( u(k,0,i) - u_m_s(k,0,i) ) / ( denom * tsc(2) )
470                   IF ( c_u(k,i) < 0.0_wp )  THEN
471                      c_u(k,i) = 0.0_wp
472                   ELSEIF ( c_u(k,i) > c_max )  THEN
473                      c_u(k,i) = c_max
474                   ENDIF
475                ELSE
476                   c_u(k,i) = c_max
477                ENDIF
478
479                denom = v_m_s(k,1,i) - v_m_s(k,2,i)
480
481                IF ( denom /= 0.0_wp )  THEN
482                   c_v(k,i) = -c_max * ( v(k,1,i) - v_m_s(k,1,i) ) / ( denom * tsc(2) )
483                   IF ( c_v(k,i) < 0.0_wp )  THEN
484                      c_v(k,i) = 0.0_wp
485                   ELSEIF ( c_v(k,i) > c_max )  THEN
486                      c_v(k,i) = c_max
487                   ENDIF
488                ELSE
489                   c_v(k,i) = c_max
490                ENDIF
491
492                denom = w_m_s(k,0,i) - w_m_s(k,1,i)
493
494                IF ( denom /= 0.0_wp )  THEN
495                   c_w(k,i) = -c_max * ( w(k,0,i) - w_m_s(k,0,i) ) / ( denom * tsc(2) )
496                   IF ( c_w(k,i) < 0.0_wp )  THEN
497                      c_w(k,i) = 0.0_wp
498                   ELSEIF ( c_w(k,i) > c_max )  THEN
499                      c_w(k,i) = c_max
500                   ENDIF
501                ELSE
502                   c_w(k,i) = c_max
503                ENDIF
504
505                c_u_m_l(k) = c_u_m_l(k) + c_u(k,i)
506                c_v_m_l(k) = c_v_m_l(k) + c_v(k,i)
507                c_w_m_l(k) = c_w_m_l(k) + c_w(k,i)
508
509             ENDDO
510          ENDDO
511
512#if defined( __parallel )   
513          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
514          CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
515                              MPI_SUM, comm1dx, ierr )   
516          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
517          CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_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_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
521                              MPI_SUM, comm1dx, ierr ) 
522#else
523          c_u_m = c_u_m_l
524          c_v_m = c_v_m_l
525          c_w_m = c_w_m_l
526#endif
527
528          c_u_m = c_u_m / (nx+1)
529          c_v_m = c_v_m / (nx+1)
530          c_w_m = c_w_m / (nx+1)
531
532!
533!--       Save old timelevels for the next timestep
534          IF ( intermediate_timestep_count == 1 )  THEN
535             u_m_s(:,:,:) = u(:,0:1,:)
536             v_m_s(:,:,:) = v(:,1:2,:)
537             w_m_s(:,:,:) = w(:,0:1,:)
538          ENDIF
539
540!
541!--       Calculate the new velocities
542          DO  k = nzb+1, nzt+1
543             DO  i = nxlg, nxrg
544                u_p(k,-1,i) = u(k,-1,i) - dt_3d * tsc(2) * c_u_m(k) *          &
545                                       ( u(k,-1,i) - u(k,0,i) ) * ddy
546
547                v_p(k,0,i)  = v(k,0,i)  - dt_3d * tsc(2) * c_v_m(k) *          &
548                                       ( v(k,0,i) - v(k,1,i) ) * ddy
549
550                w_p(k,-1,i) = w(k,-1,i) - dt_3d * tsc(2) * c_w_m(k) *          &
551                                       ( w(k,-1,i) - w(k,0,i) ) * ddy
552             ENDDO
553          ENDDO
554
555!
556!--       Bottom boundary at the outflow
557          IF ( ibc_uv_b == 0 )  THEN
558             u_p(nzb,-1,:) = 0.0_wp 
559             v_p(nzb,0,:)  = 0.0_wp 
560          ELSE                   
561             u_p(nzb,-1,:) =  u_p(nzb+1,-1,:)
562             v_p(nzb,0,:)  =  v_p(nzb+1,0,:)
563          ENDIF
564          w_p(nzb,-1,:) = 0.0_wp
565
566!
567!--       Top boundary at the outflow
568          IF ( ibc_uv_t == 0 )  THEN
569             u_p(nzt+1,-1,:) = u_init(nzt+1)
570             v_p(nzt+1,0,:)  = v_init(nzt+1)
571          ELSE
572             u_p(nzt+1,-1,:) = u_p(nzt,-1,:)
573             v_p(nzt+1,0,:)  = v_p(nzt,0,:)
574          ENDIF
575          w_p(nzt:nzt+1,-1,:) = 0.0_wp
576
577       ENDIF
578
579    ENDIF
580
581    IF ( outflow_n )  THEN
582
583       IF ( use_cmax )  THEN
584          u_p(:,ny+1,:) = u(:,ny,:)
585          v_p(:,ny+1,:) = v(:,ny,:)
586          w_p(:,ny+1,:) = w(:,ny,:)         
587       ELSEIF ( .NOT. use_cmax )  THEN
588
589          c_max = dy / dt_3d
590
591          c_u_m_l = 0.0_wp 
592          c_v_m_l = 0.0_wp
593          c_w_m_l = 0.0_wp
594
595          c_u_m = 0.0_wp 
596          c_v_m = 0.0_wp
597          c_w_m = 0.0_wp
598
599!
600!--       Calculate the phase speeds for u, v, and w, first local and then
601!--       average along the outflow boundary.
602          DO  k = nzb+1, nzt+1
603             DO  i = nxl, nxr
604
605                denom = u_m_n(k,ny,i) - u_m_n(k,ny-1,i)
606
607                IF ( denom /= 0.0_wp )  THEN
608                   c_u(k,i) = -c_max * ( u(k,ny,i) - u_m_n(k,ny,i) ) / ( denom * tsc(2) )
609                   IF ( c_u(k,i) < 0.0_wp )  THEN
610                      c_u(k,i) = 0.0_wp
611                   ELSEIF ( c_u(k,i) > c_max )  THEN
612                      c_u(k,i) = c_max
613                   ENDIF
614                ELSE
615                   c_u(k,i) = c_max
616                ENDIF
617
618                denom = v_m_n(k,ny,i) - v_m_n(k,ny-1,i)
619
620                IF ( denom /= 0.0_wp )  THEN
621                   c_v(k,i) = -c_max * ( v(k,ny,i) - v_m_n(k,ny,i) ) / ( denom * tsc(2) )
622                   IF ( c_v(k,i) < 0.0_wp )  THEN
623                      c_v(k,i) = 0.0_wp
624                   ELSEIF ( c_v(k,i) > c_max )  THEN
625                      c_v(k,i) = c_max
626                   ENDIF
627                ELSE
628                   c_v(k,i) = c_max
629                ENDIF
630
631                denom = w_m_n(k,ny,i) - w_m_n(k,ny-1,i)
632
633                IF ( denom /= 0.0_wp )  THEN
634                   c_w(k,i) = -c_max * ( w(k,ny,i) - w_m_n(k,ny,i) ) / ( denom * tsc(2) )
635                   IF ( c_w(k,i) < 0.0_wp )  THEN
636                      c_w(k,i) = 0.0_wp
637                   ELSEIF ( c_w(k,i) > c_max )  THEN
638                      c_w(k,i) = c_max
639                   ENDIF
640                ELSE
641                   c_w(k,i) = c_max
642                ENDIF
643
644                c_u_m_l(k) = c_u_m_l(k) + c_u(k,i)
645                c_v_m_l(k) = c_v_m_l(k) + c_v(k,i)
646                c_w_m_l(k) = c_w_m_l(k) + c_w(k,i)
647
648             ENDDO
649          ENDDO
650
651#if defined( __parallel )   
652          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
653          CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
654                              MPI_SUM, comm1dx, ierr )   
655          IF ( collective_wait )  CALL MPI_BARRIER( comm1dx, ierr )
656          CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_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_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
660                              MPI_SUM, comm1dx, ierr ) 
661#else
662          c_u_m = c_u_m_l
663          c_v_m = c_v_m_l
664          c_w_m = c_w_m_l
665#endif
666
667          c_u_m = c_u_m / (nx+1)
668          c_v_m = c_v_m / (nx+1)
669          c_w_m = c_w_m / (nx+1)
670
671!
672!--       Save old timelevels for the next timestep
673          IF ( intermediate_timestep_count == 1 )  THEN
674                u_m_n(:,:,:) = u(:,ny-1:ny,:)
675                v_m_n(:,:,:) = v(:,ny-1:ny,:)
676                w_m_n(:,:,:) = w(:,ny-1:ny,:)
677          ENDIF
678
679!
680!--       Calculate the new velocities
681          DO  k = nzb+1, nzt+1
682             DO  i = nxlg, nxrg
683                u_p(k,ny+1,i) = u(k,ny+1,i) - dt_3d * tsc(2) * c_u_m(k) *      &
684                                       ( u(k,ny+1,i) - u(k,ny,i) ) * ddy
685
686                v_p(k,ny+1,i) = v(k,ny+1,i)  - dt_3d * tsc(2) * c_v_m(k) *     &
687                                       ( v(k,ny+1,i) - v(k,ny,i) ) * ddy
688
689                w_p(k,ny+1,i) = w(k,ny+1,i) - dt_3d * tsc(2) * c_w_m(k) *      &
690                                       ( w(k,ny+1,i) - w(k,ny,i) ) * ddy
691             ENDDO
692          ENDDO
693
694!
695!--       Bottom boundary at the outflow
696          IF ( ibc_uv_b == 0 )  THEN
697             u_p(nzb,ny+1,:) = 0.0_wp
698             v_p(nzb,ny+1,:) = 0.0_wp   
699          ELSE                   
700             u_p(nzb,ny+1,:) =  u_p(nzb+1,ny+1,:)
701             v_p(nzb,ny+1,:) =  v_p(nzb+1,ny+1,:)
702          ENDIF
703          w_p(nzb,ny+1,:) = 0.0_wp
704
705!
706!--       Top boundary at the outflow
707          IF ( ibc_uv_t == 0 )  THEN
708             u_p(nzt+1,ny+1,:) = u_init(nzt+1)
709             v_p(nzt+1,ny+1,:) = v_init(nzt+1)
710          ELSE
711             u_p(nzt+1,ny+1,:) = u_p(nzt,nyn+1,:)
712             v_p(nzt+1,ny+1,:) = v_p(nzt,nyn+1,:)
713          ENDIF
714          w_p(nzt:nzt+1,ny+1,:) = 0.0_wp
715
716       ENDIF
717
718    ENDIF
719
720    IF ( outflow_l )  THEN
721
722       IF ( use_cmax )  THEN
723          u_p(:,:,0)  = u(:,:,1)
724          v_p(:,:,-1) = v(:,:,0)
725          w_p(:,:,-1) = w(:,:,0)         
726       ELSEIF ( .NOT. use_cmax )  THEN
727
728          c_max = dx / dt_3d
729
730          c_u_m_l = 0.0_wp 
731          c_v_m_l = 0.0_wp
732          c_w_m_l = 0.0_wp
733
734          c_u_m = 0.0_wp 
735          c_v_m = 0.0_wp
736          c_w_m = 0.0_wp
737
738!
739!--       Calculate the phase speeds for u, v, and w, first local and then
740!--       average along the outflow boundary.
741          DO  k = nzb+1, nzt+1
742             DO  j = nys, nyn
743
744                denom = u_m_l(k,j,1) - u_m_l(k,j,2)
745
746                IF ( denom /= 0.0_wp )  THEN
747                   c_u(k,j) = -c_max * ( u(k,j,1) - u_m_l(k,j,1) ) / ( denom * tsc(2) )
748                   IF ( c_u(k,j) < 0.0_wp )  THEN
749                      c_u(k,j) = 0.0_wp
750                   ELSEIF ( c_u(k,j) > c_max )  THEN
751                      c_u(k,j) = c_max
752                   ENDIF
753                ELSE
754                   c_u(k,j) = c_max
755                ENDIF
756
757                denom = v_m_l(k,j,0) - v_m_l(k,j,1)
758
759                IF ( denom /= 0.0_wp )  THEN
760                   c_v(k,j) = -c_max * ( v(k,j,0) - v_m_l(k,j,0) ) / ( denom * tsc(2) )
761                   IF ( c_v(k,j) < 0.0_wp )  THEN
762                      c_v(k,j) = 0.0_wp
763                   ELSEIF ( c_v(k,j) > c_max )  THEN
764                      c_v(k,j) = c_max
765                   ENDIF
766                ELSE
767                   c_v(k,j) = c_max
768                ENDIF
769
770                denom = w_m_l(k,j,0) - w_m_l(k,j,1)
771
772                IF ( denom /= 0.0_wp )  THEN
773                   c_w(k,j) = -c_max * ( w(k,j,0) - w_m_l(k,j,0) ) / ( denom * tsc(2) )
774                   IF ( c_w(k,j) < 0.0_wp )  THEN
775                      c_w(k,j) = 0.0_wp
776                   ELSEIF ( c_w(k,j) > c_max )  THEN
777                      c_w(k,j) = c_max
778                   ENDIF
779                ELSE
780                   c_w(k,j) = c_max
781                ENDIF
782
783                c_u_m_l(k) = c_u_m_l(k) + c_u(k,j)
784                c_v_m_l(k) = c_v_m_l(k) + c_v(k,j)
785                c_w_m_l(k) = c_w_m_l(k) + c_w(k,j)
786
787             ENDDO
788          ENDDO
789
790#if defined( __parallel )   
791          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
792          CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
793                              MPI_SUM, comm1dy, ierr )   
794          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
795          CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_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_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
799                              MPI_SUM, comm1dy, ierr ) 
800#else
801          c_u_m = c_u_m_l
802          c_v_m = c_v_m_l
803          c_w_m = c_w_m_l
804#endif
805
806          c_u_m = c_u_m / (ny+1)
807          c_v_m = c_v_m / (ny+1)
808          c_w_m = c_w_m / (ny+1)
809
810!
811!--       Save old timelevels for the next timestep
812          IF ( intermediate_timestep_count == 1 )  THEN
813                u_m_l(:,:,:) = u(:,:,1:2)
814                v_m_l(:,:,:) = v(:,:,0:1)
815                w_m_l(:,:,:) = w(:,:,0:1)
816          ENDIF
817
818!
819!--       Calculate the new velocities
820          DO  k = nzb+1, nzt+1
821             DO  j = nysg, nyng
822                u_p(k,j,0) = u(k,j,0) - dt_3d * tsc(2) * c_u_m(k) *            &
823                                       ( u(k,j,0) - u(k,j,1) ) * ddx
824
825                v_p(k,j,-1) = v(k,j,-1) - dt_3d * tsc(2) * c_v_m(k) *          &
826                                       ( v(k,j,-1) - v(k,j,0) ) * ddx
827
828                w_p(k,j,-1) = w(k,j,-1) - dt_3d * tsc(2) * c_w_m(k) *          &
829                                       ( w(k,j,-1) - w(k,j,0) ) * ddx
830             ENDDO
831          ENDDO
832
833!
834!--       Bottom boundary at the outflow
835          IF ( ibc_uv_b == 0 )  THEN
836             u_p(nzb,:,0)  = 0.0_wp 
837             v_p(nzb,:,-1) = 0.0_wp
838          ELSE                   
839             u_p(nzb,:,0)  =  u_p(nzb+1,:,0)
840             v_p(nzb,:,-1) =  v_p(nzb+1,:,-1)
841          ENDIF
842          w_p(nzb,:,-1) = 0.0_wp
843
844!
845!--       Top boundary at the outflow
846          IF ( ibc_uv_t == 0 )  THEN
847             u_p(nzt+1,:,0)  = u_init(nzt+1)
848             v_p(nzt+1,:,-1) = v_init(nzt+1)
849          ELSE
850             u_p(nzt+1,:,0)  = u_p(nzt,:,0)
851             v_p(nzt+1,:,-1) = v_p(nzt,:,-1)
852          ENDIF
853          w_p(nzt:nzt+1,:,-1) = 0.0_wp
854
855       ENDIF
856
857    ENDIF
858
859    IF ( outflow_r )  THEN
860
861       IF ( use_cmax )  THEN
862          u_p(:,:,nx+1) = u(:,:,nx)
863          v_p(:,:,nx+1) = v(:,:,nx)
864          w_p(:,:,nx+1) = w(:,:,nx)         
865       ELSEIF ( .NOT. use_cmax )  THEN
866
867          c_max = dx / dt_3d
868
869          c_u_m_l = 0.0_wp 
870          c_v_m_l = 0.0_wp
871          c_w_m_l = 0.0_wp
872
873          c_u_m = 0.0_wp 
874          c_v_m = 0.0_wp
875          c_w_m = 0.0_wp
876
877!
878!--       Calculate the phase speeds for u, v, and w, first local and then
879!--       average along the outflow boundary.
880          DO  k = nzb+1, nzt+1
881             DO  j = nys, nyn
882
883                denom = u_m_r(k,j,nx) - u_m_r(k,j,nx-1)
884
885                IF ( denom /= 0.0_wp )  THEN
886                   c_u(k,j) = -c_max * ( u(k,j,nx) - u_m_r(k,j,nx) ) / ( denom * tsc(2) )
887                   IF ( c_u(k,j) < 0.0_wp )  THEN
888                      c_u(k,j) = 0.0_wp
889                   ELSEIF ( c_u(k,j) > c_max )  THEN
890                      c_u(k,j) = c_max
891                   ENDIF
892                ELSE
893                   c_u(k,j) = c_max
894                ENDIF
895
896                denom = v_m_r(k,j,nx) - v_m_r(k,j,nx-1)
897
898                IF ( denom /= 0.0_wp )  THEN
899                   c_v(k,j) = -c_max * ( v(k,j,nx) - v_m_r(k,j,nx) ) / ( denom * tsc(2) )
900                   IF ( c_v(k,j) < 0.0_wp )  THEN
901                      c_v(k,j) = 0.0_wp
902                   ELSEIF ( c_v(k,j) > c_max )  THEN
903                      c_v(k,j) = c_max
904                   ENDIF
905                ELSE
906                   c_v(k,j) = c_max
907                ENDIF
908
909                denom = w_m_r(k,j,nx) - w_m_r(k,j,nx-1)
910
911                IF ( denom /= 0.0_wp )  THEN
912                   c_w(k,j) = -c_max * ( w(k,j,nx) - w_m_r(k,j,nx) ) / ( denom * tsc(2) )
913                   IF ( c_w(k,j) < 0.0_wp )  THEN
914                      c_w(k,j) = 0.0_wp
915                   ELSEIF ( c_w(k,j) > c_max )  THEN
916                      c_w(k,j) = c_max
917                   ENDIF
918                ELSE
919                   c_w(k,j) = c_max
920                ENDIF
921
922                c_u_m_l(k) = c_u_m_l(k) + c_u(k,j)
923                c_v_m_l(k) = c_v_m_l(k) + c_v(k,j)
924                c_w_m_l(k) = c_w_m_l(k) + c_w(k,j)
925
926             ENDDO
927          ENDDO
928
929#if defined( __parallel )   
930          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
931          CALL MPI_ALLREDUCE( c_u_m_l(nzb+1), c_u_m(nzb+1), nzt-nzb, MPI_REAL, &
932                              MPI_SUM, comm1dy, ierr )   
933          IF ( collective_wait )  CALL MPI_BARRIER( comm1dy, ierr )
934          CALL MPI_ALLREDUCE( c_v_m_l(nzb+1), c_v_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_w_m_l(nzb+1), c_w_m(nzb+1), nzt-nzb, MPI_REAL, &
938                              MPI_SUM, comm1dy, ierr ) 
939#else
940          c_u_m = c_u_m_l
941          c_v_m = c_v_m_l
942          c_w_m = c_w_m_l
943#endif
944
945          c_u_m = c_u_m / (ny+1)
946          c_v_m = c_v_m / (ny+1)
947          c_w_m = c_w_m / (ny+1)
948
949!
950!--       Save old timelevels for the next timestep
951          IF ( intermediate_timestep_count == 1 )  THEN
952                u_m_r(:,:,:) = u(:,:,nx-1:nx)
953                v_m_r(:,:,:) = v(:,:,nx-1:nx)
954                w_m_r(:,:,:) = w(:,:,nx-1:nx)
955          ENDIF
956
957!
958!--       Calculate the new velocities
959          DO  k = nzb+1, nzt+1
960             DO  j = nysg, nyng
961                u_p(k,j,nx+1) = u(k,j,nx+1) - dt_3d * tsc(2) * c_u_m(k) *      &
962                                       ( u(k,j,nx+1) - u(k,j,nx) ) * ddx
963
964                v_p(k,j,nx+1) = v(k,j,nx+1) - dt_3d * tsc(2) * c_v_m(k) *      &
965                                       ( v(k,j,nx+1) - v(k,j,nx) ) * ddx
966
967                w_p(k,j,nx+1) = w(k,j,nx+1) - dt_3d * tsc(2) * c_w_m(k) *      &
968                                       ( w(k,j,nx+1) - w(k,j,nx) ) * ddx
969             ENDDO
970          ENDDO
971
972!
973!--       Bottom boundary at the outflow
974          IF ( ibc_uv_b == 0 )  THEN
975             u_p(nzb,:,nx+1) = 0.0_wp
976             v_p(nzb,:,nx+1) = 0.0_wp 
977          ELSE                   
978             u_p(nzb,:,nx+1) =  u_p(nzb+1,:,nx+1)
979             v_p(nzb,:,nx+1) =  v_p(nzb+1,:,nx+1)
980          ENDIF
981          w_p(nzb,:,nx+1) = 0.0_wp
982
983!
984!--       Top boundary at the outflow
985          IF ( ibc_uv_t == 0 )  THEN
986             u_p(nzt+1,:,nx+1) = u_init(nzt+1)
987             v_p(nzt+1,:,nx+1) = v_init(nzt+1)
988          ELSE
989             u_p(nzt+1,:,nx+1) = u_p(nzt,:,nx+1)
990             v_p(nzt+1,:,nx+1) = v_p(nzt,:,nx+1)
991          ENDIF
992          w_p(nzt:nzt+1,:,nx+1) = 0.0_wp
993
994       ENDIF
995
996    ENDIF
997
998 END SUBROUTINE boundary_conds
Note: See TracBrowser for help on using the repository browser.