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

Last change on this file since 1321 was 1321, checked in by raasch, 10 years ago

last commit documented

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