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

Last change on this file since 1381 was 1381, checked in by heinze, 8 years ago

last commit documented

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