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

Last change on this file since 1952 was 1933, checked in by hellstea, 8 years ago

last commit documented

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