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

Last change on this file since 1734 was 1718, checked in by raasch, 9 years ago

last commit documented

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