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

Last change on this file since 30 was 19, checked in by raasch, 18 years ago

preliminary version of modified boundary conditions at top

  • Property svn:keywords set to Id
File size: 9.7 KB
Line 
1 SUBROUTINE boundary_conds( range )
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6! Boundary conditions for e(nzt), pt(nzt), and q(nzt) removed because these
7! values are now calculated by the prognostic equation,
8! Dirichlet and zero gradient condition for pt established at top boundary
9!
10! Former revisions:
11! -----------------
12! $Id: boundary_conds.f90 19 2007-02-23 04:53:48Z raasch $
13! RCS Log replace by Id keyword, revision history cleaned up
14!
15! Revision 1.15  2006/02/23 09:54:55  raasch
16! Surface boundary conditions in case of topography: nzb replaced by
17! 2d-k-index-arrays (nzb_w_inner, etc.). Conditions for u and v remain
18! unchanged (still using nzb) because a non-flat topography must use a
19! Prandtl-layer, which don't requires explicit setting of the surface values.
20!
21! Revision 1.1  1997/09/12 06:21:34  raasch
22! Initial revision
23!
24!
25! Description:
26! ------------
27! Boundary conditions for the prognostic quantities (range='main').
28! In case of non-cyclic lateral boundaries the conditions for velocities at
29! the outflow are set after the pressure solver has been called (range=
30! 'outflow_uvw').
31! One additional bottom boundary condition is applied for the TKE (=(u*)**2)
32! in prandtl_fluxes. The cyclic lateral boundary conditions are implicitly
33! handled in routine exchange_horiz. Pressure boundary conditions are
34! explicitly set in routines pres, poisfft, poismg and sor.
35!------------------------------------------------------------------------------!
36
37    USE arrays_3d
38    USE control_parameters
39    USE grid_variables
40    USE indices
41    USE pegrid
42
43    IMPLICIT NONE
44
45    CHARACTER (LEN=*) ::  range
46
47    INTEGER ::  i, j, k
48
49
50    IF ( range == 'main')  THEN
51!
52!--    Bottom boundary
53       IF ( ibc_uv_b == 0 )  THEN
54          u(nzb,:,:) = -u(nzb+1,:,:)  ! satisfying the Dirichlet condition with
55          v(nzb,:,:) = -v(nzb+1,:,:)  ! an extra layer below the surface where
56       ELSE                           ! the u and v component change their sign
57          u(nzb,:,:) = u(nzb+1,:,:)
58          v(nzb,:,:) = v(nzb+1,:,:)
59       ENDIF
60       DO  i = nxl-1, nxr+1
61          DO  j = nys-1, nyn+1
62             w(nzb_w_inner(j,i),j,i) = 0.0
63          ENDDO
64       ENDDO
65
66!
67!--    Top boundary
68       IF ( ibc_uv_t == 0 )  THEN
69          u(nzt+1,:,:) = ug(nzt+1)
70          v(nzt+1,:,:) = vg(nzt+1)
71       ELSE
72          u(nzt+1,:,:) = u(nzt,:,:)
73          v(nzt+1,:,:) = v(nzt,:,:)
74       ENDIF
75       w(nzt:nzt+1,:,:) = 0.0  ! nzt is not a prognostic level (but cf. pres)
76
77!
78!--    Temperature at bottom boundary
79       IF ( ibc_pt_b == 0 )  THEN
80          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
81             DO  i = nxl-1, nxr+1
82                DO  j = nys-1, nyn+1
83                   pt(nzb_s_inner(j,i),j,i) = pt_m(nzb_s_inner(j,i),j,i)
84                ENDDO
85             ENDDO
86          ELSE
87             DO  i = nxl-1, nxr+1
88                DO  j = nys-1, nyn+1
89                   pt(nzb_s_inner(j,i),j,i) = pt_p(nzb_s_inner(j,i),j,i) 
90                                            ! pt_m is not used for Runge-Kutta
91                ENDDO
92             ENDDO
93          ENDIF
94       ELSE
95          DO  i = nxl-1, nxr+1
96             DO  j = nys-1, nyn+1
97                pt(nzb_s_inner(j,i),j,i) = pt(nzb_s_inner(j,i)+1,j,i)
98             ENDDO
99          ENDDO
100       ENDIF
101
102!
103!--    Temperature at top boundary
104       IF ( ibc_pt_t == 0 )  THEN
105          IF ( timestep_scheme(1:5) /= 'runge' )  THEN
106             pt(nzt+1,:,:) = pt_m(nzt+1,:,:)
107          ELSE
108             pt(nzt+1,:,:) = pt_p(nzt+1,:,:)  ! pt_m not used for Runge-Kutta
109          ENDIF
110       ELSEIF ( ibc_pt_t == 1 )  THEN
111          pt(nzt+1,:,:) = pt(nzt,:,:)
112       ELSEIF ( ibc_pt_t == 2 )  THEN
113          pt(nzt+1,:,:) = pt(nzt,:,:)   + bc_pt_t_val * dzu(nzt+1)
114       ENDIF
115
116!
117!--    Boundary conditions for TKE
118!--    Generally Neumann conditions with de/dz=0 are assumed
119       IF ( .NOT. constant_diffusion )  THEN
120          DO  i = nxl-1, nxr+1
121             DO  j = nys-1, nyn+1
122                e(nzb_s_inner(j,i),j,i) = e(nzb_s_inner(j,i)+1,j,i)
123             ENDDO
124          ENDDO
125          e(nzt+1,:,:) = e(nzt,:,:)
126       ENDIF
127
128!
129!--    Boundary conditions for total water content or scalar,
130!--    bottom and surface boundary (see also temperature)
131       IF ( moisture  .OR.  passive_scalar )  THEN
132!
133!--       Surface conditions for constant_moisture_flux
134          IF ( ibc_q_b == 0 ) THEN
135             IF ( timestep_scheme(1:5) /= 'runge' )  THEN
136                DO  i = nxl-1, nxr+1
137                   DO  j = nys-1, nyn+1
138                      q(nzb_s_inner(j,i),j,i) = q_m(nzb_s_inner(j,i),j,i)
139                   ENDDO
140                ENDDO
141             ELSE
142                DO  i = nxl-1, nxr+1
143                   DO  j = nys-1, nyn+1
144                      q(nzb_s_inner(j,i),j,i) = q_p(nzb_s_inner(j,i),j,i)
145                   ENDDO                      ! q_m is not used for Runge-Kutta
146                ENDDO
147             ENDIF
148          ELSE
149             DO  i = nxl-1, nxr+1
150                DO  j = nys-1, nyn+1
151                   q(nzb_s_inner(j,i),j,i) = q(nzb_s_inner(j,i)+1,j,i)
152                ENDDO
153             ENDDO
154          ENDIF
155!
156!--       Top boundary
157          q(nzt+1,:,:) = q(nzt,:,:)   + bc_q_t_val * dzu(nzt+1)
158       ENDIF
159
160!
161!--    Lateral boundary conditions at the inflow. Quasi Neumann conditions
162!--    are needed for the wall normal velocity in order to ensure zero
163!--    divergence. Dirichlet conditions are used for all other quantities.
164       IF ( inflow_s )  THEN
165          v(:,nys,:) = v(:,nys-1,:)
166       ELSEIF ( inflow_n ) THEN
167          v(:,nyn+vynp,:) = v(:,nyn+vynp+1,:)
168       ELSEIF ( inflow_l ) THEN
169          u(:,:,nxl) = u(:,:,nxl-1)
170       ELSEIF ( inflow_r ) THEN
171          u(:,:,nxr+uxrp) = u(:,:,nxr+uxrp+1)
172       ENDIF
173
174!
175!--    Lateral boundary conditions for scalar quantities at the outflow
176       IF ( outflow_s )  THEN
177          pt(:,nys-1,:)     = pt(:,nys,:)
178          IF ( .NOT. constant_diffusion     )  e(:,nys-1,:) = e(:,nys,:)
179          IF ( moisture .OR. passive_scalar )  q(:,nys-1,:) = q(:,nys,:)
180       ELSEIF ( outflow_n )  THEN
181          pt(:,nyn+1,:)     = pt(:,nyn,:)
182          IF ( .NOT. constant_diffusion     )  e(:,nyn+1,:) = e(:,nyn,:)
183          IF ( moisture .OR. passive_scalar )  q(:,nyn+1,:) = q(:,nyn,:)
184       ELSEIF ( outflow_l )  THEN
185          pt(:,:,nxl-1)     = pt(:,:,nxl)
186          IF ( .NOT. constant_diffusion     )  e(:,:,nxl-1) = e(:,:,nxl)
187          IF ( moisture .OR. passive_scalar )  q(:,:,nxl-1) = q(:,:,nxl)     
188       ELSEIF ( outflow_r )  THEN
189          pt(:,:,nxr+1)     = pt(:,:,nxr)
190          IF ( .NOT. constant_diffusion     )  e(:,:,nxr+1) = e(:,:,nxr)
191          IF ( moisture .OR. passive_scalar )  q(:,:,nxr+1) = q(:,:,nxr)
192       ENDIF
193
194    ENDIF
195
196    IF ( range == 'outflow_uvw' ) THEN
197!
198!--    Horizontal boundary conditions for the velocities at the outflow.
199!--    A Neumann condition is used for the wall normal velocity. The vertical
200!--    velocity is assumed as zero and a horizontal average along the wall is
201!--    used for the wall parallel horizontal velocity. The combination of all
202!--    three conditions ensures that the velocity field is free of divergence
203!--    at the outflow (uvmean_outflow_l is calculated in pres).
204       IF ( outflow_s ) THEN
205          v(:,nys-1,:) = v(:,nys,:)
206          w(:,nys-1,:) = 0.0
207!
208!--       Compute the mean horizontal wind parallel to and within the outflow
209!--       wall and use this as boundary condition for u
210#if defined( __parallel )
211          CALL MPI_ALLREDUCE( uvmean_outflow_l, uvmean_outflow, nzt-nzb+2, &
212                              MPI_REAL, MPI_SUM, comm1dx, ierr )
213          uvmean_outflow = uvmean_outflow / ( nx + 1.0 )
214#else
215          uvmean_outflow = uvmean_outflow_l / ( nx + 1.0 )
216#endif
217          DO  k = nzb, nzt+1
218             u(k,nys-1,:) = uvmean_outflow(k)
219          ENDDO
220       ENDIF
221
222       IF ( outflow_n ) THEN
223          v(:,nyn+vynp+1,:) = v(:,nyn+vynp,:)
224          w(:,nyn+1,:)      = 0.0
225!
226!--       Compute the mean horizontal wind parallel to and within the outflow
227!--       wall and use this as boundary condition for u
228#if defined( __parallel )
229          CALL MPI_ALLREDUCE( uvmean_outflow_l, uvmean_outflow, nzt-nzb+2, &
230                              MPI_REAL, MPI_SUM, comm1dx, ierr )
231          uvmean_outflow = uvmean_outflow / ( nx + 1.0 )
232#else
233          uvmean_outflow = uvmean_outflow_l / ( nx + 1.0 )
234#endif
235          DO  k = nzb, nzt+1
236             u(k,nyn+1,:) = uvmean_outflow(k)
237          ENDDO
238       ENDIF
239
240       IF ( outflow_l ) THEN
241          u(:,:,nxl-1) = u(:,:,nxl)
242          w(:,:,nxl-1) = 0.0
243!
244!--       Compute the mean horizontal wind parallel to and within the outflow
245!--       wall and use this as boundary condition for v
246#if defined( __parallel )
247          CALL MPI_ALLREDUCE( uvmean_outflow_l, uvmean_outflow, nzt-nzb+2, &
248                              MPI_REAL, MPI_SUM, comm1dy, ierr )
249          uvmean_outflow = uvmean_outflow / ( ny + 1.0 )
250#else
251          uvmean_outflow = uvmean_outflow_l / ( ny + 1.0 )
252#endif
253          DO  k = nzb, nzt+1
254             v(k,:,nxl-1) = uvmean_outflow(k)
255          ENDDO   
256
257       ENDIF
258
259       IF ( outflow_r ) THEN
260          u(:,:,nxr+uxrp+1) = u(:,:,nxr+uxrp)
261          w(:,:,nxr+1) = 0.0
262!
263!--       Compute the mean horizontal wind parallel to and within the outflow
264!--       wall and use this as boundary condition for v
265#if defined( __parallel )
266          CALL MPI_ALLREDUCE( uvmean_outflow_l, uvmean_outflow, nzt-nzb+2, &
267                              MPI_REAL, MPI_SUM, comm1dy, ierr )
268          uvmean_outflow = uvmean_outflow / ( ny + 1.0 )
269#else
270          uvmean_outflow = uvmean_outflow_l / ( ny + 1.0 )
271#endif
272          DO  k = nzb, nzt+1
273             v(k,:,nxr+1) = uvmean_outflow(k)
274          ENDDO   
275       ENDIF
276
277    ENDIF
278
279 
280 END SUBROUTINE boundary_conds
Note: See TracBrowser for help on using the repository browser.