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

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

update of GPL copyright

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