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

Last change on this file since 1116 was 1116, checked in by hoffmann, 12 years ago

last commit documented

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