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

Last change on this file since 1713 was 1683, checked in by knoop, 9 years ago

last commit documented

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