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

Last change on this file since 51 was 39, checked in by raasch, 18 years ago

comments prepared for 3.1c

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