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

Last change on this file since 73 was 73, checked in by raasch, 15 years ago

preliminary changes for radiation conditions

  • Property svn:keywords set to Id
File size: 13.8 KB
Line 
1 SUBROUTINE boundary_conds( range )
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6! The "main" part sets conditions for time level t+dt insteat of level t,
7! outflow boundary conditions changed from Neumann to radiation condition
8!
9! Former revisions:
10! -----------------
11! $Id: boundary_conds.f90 73 2007-03-20 08:33:14Z raasch $
12!
13! 19 2007-02-23 04:53:48Z raasch
14! Boundary conditions for e(nzt), pt(nzt), and q(nzt) removed because these
15! gridpoints are now calculated by the prognostic equation,
16! Dirichlet and zero gradient condition for pt established at top boundary
17!
18! RCS Log replace by Id keyword, revision history cleaned up
19!
20! Revision 1.15  2006/02/23 09:54:55  raasch
21! Surface boundary conditions in case of topography: nzb replaced by
22! 2d-k-index-arrays (nzb_w_inner, etc.). Conditions for u and v remain
23! unchanged (still using nzb) because a non-flat topography must use a
24! Prandtl-layer, which don't requires explicit setting of the surface values.
25!
26! Revision 1.1  1997/09/12 06:21:34  raasch
27! Initial revision
28!
29!
30! Description:
31! ------------
32! Boundary conditions for the prognostic quantities (range='main').
33! In case of non-cyclic lateral boundaries the conditions for velocities at
34! the outflow are set after the pressure solver has been called (range=
35! 'outflow_uvw').
36! One additional bottom boundary condition is applied for the TKE (=(u*)**2)
37! in prandtl_fluxes. The cyclic lateral boundary conditions are implicitly
38! handled in routine exchange_horiz. Pressure boundary conditions are
39! explicitly set in routines pres, poisfft, poismg and sor.
40!------------------------------------------------------------------------------!
41
42    USE arrays_3d
43    USE control_parameters
44    USE grid_variables
45    USE indices
46    USE pegrid
47
48    IMPLICIT NONE
49
50    CHARACTER (LEN=*) ::  range
51
52    INTEGER ::  i, j, k
53
54    REAL    ::  c_max, c_u, c_v, c_w, denom
55
56
57    IF ( range == 'main')  THEN
58!
59!--    Bottom boundary
60       IF ( ibc_uv_b == 0 )  THEN
61!
62!--       Satisfying the Dirichlet condition with an extra layer below the
63!--       surface where the u and v component change their sign
64          u_p(nzb,:,:) = -u_p(nzb+1,:,:)
65          v_p(nzb,:,:) = -v_p(nzb+1,:,:)
66       ELSE
67          u_p(nzb,:,:) = u_p(nzb+1,:,:)
68          v_p(nzb,:,:) = v_p(nzb+1,:,:)
69       ENDIF
70       DO  i = nxl-1, nxr+1
71          DO  j = nys-1, nyn+1
72             w_p(nzb_w_inner(j,i),j,i) = 0.0
73          ENDDO
74       ENDDO
75
76!
77!--    Top boundary
78       IF ( ibc_uv_t == 0 )  THEN
79          u_p(nzt+1,:,:) = ug(nzt+1)
80          v_p(nzt+1,:,:) = vg(nzt+1)
81       ELSE
82          u_p(nzt+1,:,:) = u_p(nzt,:,:)
83          v_p(nzt+1,:,:) = v_p(nzt,:,:)
84       ENDIF
85       w_p(nzt:nzt+1,:,:) = 0.0  ! nzt is not a prognostic level (but cf. pres)
86
87!
88!--    Temperature at bottom boundary
89       IF ( ibc_pt_b == 0 )  THEN
90          DO  i = nxl-1, nxr+1
91             DO  j = nys-1, nyn+1
92                pt_p(nzb_s_inner(j,i),j,i) = pt(nzb_s_inner(j,i),j,i)
93             ENDDO
94          ENDDO
95       ELSE
96          DO  i = nxl-1, nxr+1
97             DO  j = nys-1, nyn+1
98                pt_p(nzb_s_inner(j,i),j,i) = pt_p(nzb_s_inner(j,i)+1,j,i)
99             ENDDO
100          ENDDO
101       ENDIF
102
103!
104!--    Temperature at top boundary
105       IF ( ibc_pt_t == 0 )  THEN
106          pt_p(nzt+1,:,:) = pt(nzt+1,:,:)
107       ELSEIF ( ibc_pt_t == 1 )  THEN
108          pt_p(nzt+1,:,:) = pt_p(nzt,:,:)
109       ELSEIF ( ibc_pt_t == 2 )  THEN
110          pt_p(nzt+1,:,:) = pt_p(nzt,:,:)   + bc_pt_t_val * dzu(nzt+1)
111       ENDIF
112
113!
114!--    Boundary conditions for TKE
115!--    Generally Neumann conditions with de/dz=0 are assumed
116       IF ( .NOT. constant_diffusion )  THEN
117          DO  i = nxl-1, nxr+1
118             DO  j = nys-1, nyn+1
119                e_p(nzb_s_inner(j,i),j,i) = e_p(nzb_s_inner(j,i)+1,j,i)
120             ENDDO
121          ENDDO
122          e_p(nzt+1,:,:) = e_p(nzt,:,:)
123       ENDIF
124
125!
126!--    Boundary conditions for total water content or scalar,
127!--    bottom and surface boundary (see also temperature)
128       IF ( moisture  .OR.  passive_scalar )  THEN
129!
130!--       Surface conditions for constant_moisture_flux
131          IF ( ibc_q_b == 0 ) THEN
132             DO  i = nxl-1, nxr+1
133                DO  j = nys-1, nyn+1
134                   q_p(nzb_s_inner(j,i),j,i) = q(nzb_s_inner(j,i),j,i)
135                ENDDO
136             ENDDO
137          ELSE
138             DO  i = nxl-1, nxr+1
139                DO  j = nys-1, nyn+1
140                   q_p(nzb_s_inner(j,i),j,i) = q_p(nzb_s_inner(j,i)+1,j,i)
141                ENDDO
142             ENDDO
143          ENDIF
144!
145!--       Top boundary
146          q_p(nzt+1,:,:) = q_p(nzt,:,:)   + bc_q_t_val * dzu(nzt+1)
147       ENDIF
148
149!
150!--    Lateral boundary conditions at the inflow. Quasi Neumann conditions
151!--    are needed for the wall normal velocity in order to ensure zero
152!--    divergence. Dirichlet conditions are used for all other quantities.
153       IF ( inflow_s )  THEN
154          v_p(:,nys,:) = v_p(:,nys-1,:)
155       ELSEIF ( inflow_n ) THEN
156          v_p(:,nyn+vynp,:) = v_p(:,nyn+vynp+1,:)
157       ELSEIF ( inflow_l ) THEN
158          u_p(:,:,nxl) = u_p(:,:,nxl-1)
159       ELSEIF ( inflow_r ) THEN
160          u_p(:,:,nxr+uxrp) = u_p(:,:,nxr+uxrp+1)
161       ENDIF
162
163!
164!--    Lateral boundary conditions for scalar quantities at the outflow
165       IF ( outflow_s )  THEN
166          pt_p(:,nys-1,:)     = pt_p(:,nys,:)
167          IF ( .NOT. constant_diffusion     )  e_p(:,nys-1,:) = e_p(:,nys,:)
168          IF ( moisture .OR. passive_scalar )  q_p(:,nys-1,:) = q_p(:,nys,:)
169       ELSEIF ( outflow_n )  THEN
170          pt_p(:,nyn+1,:)     = pt_p(:,nyn,:)
171          IF ( .NOT. constant_diffusion     )  e_p(:,nyn+1,:) = e_p(:,nyn,:)
172          IF ( moisture .OR. passive_scalar )  q_p(:,nyn+1,:) = q_p(:,nyn,:)
173       ELSEIF ( outflow_l )  THEN
174          pt_p(:,:,nxl-1)     = pt_p(:,:,nxl)
175          IF ( .NOT. constant_diffusion     )  e_p(:,:,nxl-1) = e_p(:,:,nxl)
176          IF ( moisture .OR. passive_scalar )  q_p(:,:,nxl-1) = q_p(:,:,nxl)     
177       ELSEIF ( outflow_r )  THEN
178          pt_p(:,:,nxr+1)     = pt_p(:,:,nxr)
179          IF ( .NOT. constant_diffusion     )  e_p(:,:,nxr+1) = e_p(:,:,nxr)
180          IF ( moisture .OR. passive_scalar )  q_p(:,:,nxr+1) = q_p(:,:,nxr)
181       ENDIF
182
183    ENDIF
184
185    IF ( range == 'outflow_uvw' ) THEN
186!
187!--    Radiation boundary condition for the velocities at the respective outflow
188       IF ( outflow_s ) THEN
189!          v(:,nys-1,:) = v(:,nys,:)
190!          w(:,nys-1,:) = 0.0
191!!
192!!--       Compute the mean horizontal wind parallel to and within the outflow
193!!--       wall and use this as boundary condition for u
194!#if defined( __parallel )
195!          CALL MPI_ALLREDUCE( uvmean_outflow_l, uvmean_outflow, nzt-nzb+2, &
196!                              MPI_REAL, MPI_SUM, comm1dx, ierr )
197!          uvmean_outflow = uvmean_outflow / ( nx + 1.0 )
198!#else
199!          uvmean_outflow = uvmean_outflow_l / ( nx + 1.0 )
200!#endif
201!          DO  k = nzb, nzt+1
202!             u(k,nys-1,:) = uvmean_outflow(k)
203!          ENDDO
204       ENDIF
205
206       IF ( outflow_n  .AND.                                                 &
207            intermediate_timestep_count == intermediate_timestep_count_max ) &
208       THEN
209
210          c_max = dy / dt_3d
211
212          DO i = nxl-1, nxr+1
213             DO k = nzb+1, nzt+1
214
215!
216!--             First calculate the phase speeds for u,v, and w
217                denom = u_m_n(k,ny,i,-2) - u_m_n(k,ny-1,i,-2)
218
219                IF ( denom /= 0.0 )  THEN
220                   c_u = -c_max * ( u_m_n(k,ny,i,-1)-u_m_n(k,ny,i,-2) ) / denom
221                   IF ( c_u < 0.0 )  THEN
222                      c_u = 0.0
223                   ELSEIF ( c_u > c_max )  THEN
224                      c_u = c_max
225                   ENDIF
226                ELSE
227                   c_u = c_max
228                ENDIF
229
230                denom = v_m_n(k,ny,i,-2) - v_m_n(k,ny-1,i,-2)
231
232                IF ( denom /= 0.0 )  THEN
233                   c_v = -c_max * ( v_m_n(k,ny,i,-1)-v_m_n(k,ny,i,-2) ) / denom
234                   IF ( c_v < 0.0 )  THEN
235                      c_v = 0.0
236                   ELSEIF ( c_v > c_max )  THEN
237                      c_v = c_max
238                   ENDIF
239                ELSE
240                   c_v = c_max
241                ENDIF
242
243                denom = w_m_n(k,ny,i,-2) - w_m_n(k,ny-1,i,-2)
244
245                IF ( denom /= 0.0 )  THEN
246                   c_w = -c_max * ( w_m_n(k,ny,i,-1)-w_m_n(k,ny,i,-2) ) / denom
247                   IF ( c_w < 0.0 )  THEN
248                      c_w = 0.0
249                   ELSEIF ( c_w > c_max )  THEN
250                      c_w = c_max
251                   ENDIF
252                ELSE
253                   c_w = c_max
254                ENDIF
255
256!
257!--             Calculate the new velocities
258                u(k,ny+1,i) = u_m_n(k,ny+1,i,-1) - dt_3d * c_u * &
259                              ( u_m_n(k,ny+1,i,-1) - u_m_n(k,ny,i,-1) ) * ddy
260
261                v(k,ny+1,i) = v_m_n(k,ny+1,i,-1) - dt_3d * c_v * &
262                              ( v_m_n(k,ny+1,i,-1) - v_m_n(k,ny,i,-1) ) * ddy
263
264                w(k,ny+1,i) = w_m_n(k,ny+1,i,-1) - dt_3d * c_w * &
265                              ( w_m_n(k,ny+1,i,-1) - w_m_n(k,ny,i,-1) ) * ddy
266
267!
268!--             Swap timelevels for the next timestep
269                u_m_n(k,:,i,-2) = u_m_n(k,:,i,-1)
270                u_m_n(k,:,i,-1) = u(k,ny-1:ny+1,i)
271                v_m_n(k,:,i,-2) = v_m_n(k,:,i,-1)
272                v_m_n(k,:,i,-1) = v(k,ny-1:ny+1,i)
273                w_m_n(k,:,i,-2) = w_m_n(k,:,i,-1)
274                w_m_n(k,:,i,-1) = w(k,ny-1:ny+1,i)
275
276             ENDDO
277          ENDDO
278
279!
280!--       Bottom boundary at the outflow
281          IF ( ibc_uv_b == 0 )  THEN
282             u(nzb,ny+1,:) = -u(nzb+1,ny+1,:) 
283             v(nzb,ny+1,:) = -v(nzb+1,ny+1,:) 
284          ELSE                   
285             u(nzb,ny+1,:) =  u(nzb+1,ny+1,:)
286             v(nzb,ny+1,:) =  v(nzb+1,ny+1,:)
287          ENDIF
288          w(nzb,ny+1,:) = 0.0
289
290!
291!--       Top boundary at the outflow
292          IF ( ibc_uv_t == 0 )  THEN
293             u(nzt+1,ny+1,:) = ug(nzt+1)
294             v(nzt+1,ny+1,:) = vg(nzt+1)
295          ELSE
296             u(nzt+1,ny+1,:) = u(nzt,nyn+1,:)
297             v(nzt+1,ny+1,:) = v(nzt,nyn+1,:)
298          ENDIF
299          w(nzt:nzt+1,ny+1,:) = 0.0
300
301       ENDIF
302
303       IF ( outflow_l ) THEN
304!          u(:,:,nxl-1) = u(:,:,nxl)
305!          w(:,:,nxl-1) = 0.0
306!
307!--       Compute the mean horizontal wind parallel to and within the outflow
308!--       wall and use this as boundary condition for v
309!#if defined( __parallel )
310!          CALL MPI_ALLREDUCE( uvmean_outflow_l, uvmean_outflow, nzt-nzb+2, &
311!                              MPI_REAL, MPI_SUM, comm1dy, ierr )
312!          uvmean_outflow = uvmean_outflow / ( ny + 1.0 )
313!#else
314!          uvmean_outflow = uvmean_outflow_l / ( ny + 1.0 )
315!#endif
316!          DO  k = nzb, nzt+1
317!             v(k,:,nxl-1) = uvmean_outflow(k)
318!          ENDDO   
319!
320       ENDIF
321
322       IF ( outflow_r  .AND.                                                 &
323            intermediate_timestep_count == intermediate_timestep_count_max ) &
324       THEN
325
326          c_max = dx / dt_3d
327
328          DO j = nys-1, nyn+1
329             DO k = nzb+1, nzt+1
330
331!
332!--             First calculate the phase speeds for u,v, and w
333                denom = u_m_r(k,j,nx,-2) - u_m_r(k,j,nx-1,-2)
334
335                IF ( denom /= 0.0 )  THEN
336                   c_u = -c_max * ( u_m_r(k,j,nx,-1)-u_m_r(k,j,nx,-2) ) / denom
337                   IF ( c_u < 0.0 )  THEN
338                      c_u = 0.0
339                   ELSEIF ( c_u > c_max )  THEN
340                      c_u = c_max
341                   ENDIF
342                ELSE
343                   c_u = c_max
344                ENDIF
345
346                denom = v_m_r(k,j,nx,-2) - v_m_r(k,j,nx-1,-2)
347
348                IF ( denom /= 0.0 )  THEN
349                   c_v = -c_max * ( v_m_r(k,j,nx,-1)-v_m_r(k,j,nx,-2) ) / denom
350                   IF ( c_v < 0.0 )  THEN
351                      c_v = 0.0
352                   ELSEIF ( c_v > c_max )  THEN
353                      c_v = c_max
354                   ENDIF
355                ELSE
356                   c_v = c_max
357                ENDIF
358
359                denom = w_m_r(k,j,nx,-2) - w_m_r(k,j,nx-1,-2)
360
361                IF ( denom /= 0.0 )  THEN
362                   c_w = -c_max * ( w_m_r(k,j,nx,-1)-w_m_n(k,j,nx,-2) ) / denom
363                   IF ( c_w < 0.0 )  THEN
364                      c_w = 0.0
365                   ELSEIF ( c_w > c_max )  THEN
366                      c_w = c_max
367                   ENDIF
368                ELSE
369                   c_w = c_max
370                ENDIF
371
372!
373!--             Calculate the new velocities
374                u(k,j,nx+1) = u_m_r(k,j,nx+1,-1) - dt_3d * c_u * &
375                              ( u_m_r(k,j,nx+1,-1) - u_m_r(k,j,nx,-1) ) * ddx
376
377                v(k,j,nx+1) = v_m_r(k,j,nx+1,-1) - dt_3d * c_v * &
378                              ( v_m_r(k,j,nx+1,-1) - v_m_r(k,j,nx,-1) ) * ddx
379
380                w(k,j,nx+1) = w_m_r(k,j,nx+1,-1) - dt_3d * c_w * &
381                              ( w_m_r(k,j,nx+1,-1) - w_m_r(k,j,nx,-1) ) * ddx
382
383!
384!--             Swap timelevels for the next timestep
385                u_m_r(k,j,:,-2) = u_m_r(k,j,:,-1)
386                u_m_r(k,j,:,-1) = u(k,j,nx-1:nx+1)
387                v_m_r(k,j,:,-2) = v_m_r(k,j,:,-1)
388                v_m_r(k,j,:,-1) = v(k,j,nx-1:nx+1)
389                w_m_r(k,j,:,-2) = w_m_r(k,j,:,-1)
390                w_m_r(k,j,:,-1) = w(k,j,nx-1:nx+1)
391
392             ENDDO
393          ENDDO
394
395!
396!--       Bottom boundary at the outflow
397          IF ( ibc_uv_b == 0 )  THEN
398             u(nzb,ny+1,:) = -u(nzb+1,ny+1,:) 
399             v(nzb,ny+1,:) = -v(nzb+1,ny+1,:) 
400          ELSE                   
401             u(nzb,ny+1,:) =  u(nzb+1,ny+1,:)
402             v(nzb,ny+1,:) =  v(nzb+1,ny+1,:)
403          ENDIF
404          w(nzb,ny+1,:) = 0.0
405
406!
407!--       Top boundary at the outflow
408          IF ( ibc_uv_t == 0 )  THEN
409             u(nzt+1,ny+1,:) = ug(nzt+1)
410             v(nzt+1,ny+1,:) = vg(nzt+1)
411          ELSE
412             u(nzt+1,ny+1,:) = u(nzt,nyn+1,:)
413             v(nzt+1,ny+1,:) = v(nzt,nyn+1,:)
414          ENDIF
415          w(nzt:nzt+1,ny+1,:) = 0.0
416
417       ENDIF
418
419    ENDIF
420
421 
422 END SUBROUTINE boundary_conds
Note: See TracBrowser for help on using the repository browser.