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

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

last commit documented

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