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

Last change on this file since 1625 was 1463, checked in by suehring, 10 years ago

last commit documented

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