1 | SUBROUTINE boundary_conds( range ) |
---|
2 | |
---|
3 | !------------------------------------------------------------------------------! |
---|
4 | ! Actual revisions: |
---|
5 | ! ----------------- |
---|
6 | ! |
---|
7 | ! |
---|
8 | ! Former revisions: |
---|
9 | ! ----------------- |
---|
10 | ! $Log: boundary_conds.f90,v $ |
---|
11 | ! Revision 1.15 2006/02/23 09:54:55 raasch |
---|
12 | ! Surface boundary conditions in case of topography: nzb replaced by |
---|
13 | ! 2d-k-index-arrays (nzb_w_inner, etc.). Conditions for u and v remain |
---|
14 | ! unchanged (still using nzb) because a non-flat topography must use a |
---|
15 | ! Prandtl-layer, which don't requires explicit setting of the surface values. |
---|
16 | ! |
---|
17 | ! Revision 1.14 2005/06/29 09:49:40 steinfeld |
---|
18 | ! Baroclinicity considered in the Dirichlet boundary condition for u and v at |
---|
19 | ! the top boundary |
---|
20 | ! |
---|
21 | ! Revision 1.13 2005/04/23 08:35:30 raasch |
---|
22 | ! Error removed in Dirichlet bottom boundary conditions for pt and q in case |
---|
23 | ! of Runge-Kutta schemes. So far, new surface values have been taken from |
---|
24 | ! timelevel t-dt, which does not exist in case of Runge-Kutta schemes. |
---|
25 | ! |
---|
26 | ! Revision 1.12 2005/03/26 14:58:24 raasch |
---|
27 | ! Non-cyclic boundary conditions included, argument range added |
---|
28 | ! |
---|
29 | ! Revision 1.11 2001/03/30 06:54:18 raasch |
---|
30 | ! Timelevel t+dt replaced by timelevel t, |
---|
31 | ! Translation of remaining German identifiers (variables, subroutines, etc.) |
---|
32 | ! |
---|
33 | ! Revision 1.10 2001/01/29 12:19:27 raasch |
---|
34 | ! Passive scalar is considered |
---|
35 | ! |
---|
36 | ! Revision 1.9 2001/01/22 05:25:56 raasch |
---|
37 | ! Module test_variables removed |
---|
38 | ! |
---|
39 | ! Revision 1.8 2000/07/04 14:07:59 raasch |
---|
40 | ! Missing diriclet boundary conditions for temperature and total water |
---|
41 | ! content added |
---|
42 | ! |
---|
43 | ! Revision 1.7 2000/04/13 13:56:05 schroeter |
---|
44 | ! Boundaray conditions for total water content |
---|
45 | ! |
---|
46 | ! Revision 1.6 2000/01/20 10:44:48 10:44:48 letzel (Marcus Letzel) |
---|
47 | ! All comments translated into English |
---|
48 | ! |
---|
49 | ! Revision 1.5 2000/01/10 09:28:53 raasch |
---|
50 | ! Variablenuebergabe jetzt per Modul |
---|
51 | ! Randbedingungen fuer w im Rahmen der pointer-Einfuehrung |
---|
52 | ! |
---|
53 | ! Revision 1.4 1998/07/06 12:07:05 raasch |
---|
54 | ! + USE test_variables |
---|
55 | ! |
---|
56 | ! Revision 1.3 1998/03/10 07:19:37 raasch |
---|
57 | ! Beschreibung ergaenzt um Randbedingung fuer TKE |
---|
58 | ! |
---|
59 | ! Revision 1.2 1997/09/19 07:39:03 raasch |
---|
60 | ! Randbedingungen fuer TKE |
---|
61 | ! |
---|
62 | ! Revision 1.1 1997/09/12 06:21:34 raasch |
---|
63 | ! Initial revision |
---|
64 | ! |
---|
65 | ! |
---|
66 | ! Description: |
---|
67 | ! ------------ |
---|
68 | ! Boundary conditions for the prognostic quantities (range='main'). |
---|
69 | ! In case of non-cyclic lateral boundaries the conditions for velocities at |
---|
70 | ! the outflow are set after the pressure solver has been called (range= |
---|
71 | ! 'outflow_uvw'). |
---|
72 | ! One additional bottom boundary condition is applied for the TKE (=(u*)**2) |
---|
73 | ! in prandtl_fluxes. The cyclic lateral boundary conditions are implicitly |
---|
74 | ! handled in routine exchange_horiz. Pressure boundary conditions are |
---|
75 | ! explicitly set in routines pres, poisfft, poismg and sor. |
---|
76 | !------------------------------------------------------------------------------! |
---|
77 | |
---|
78 | USE arrays_3d |
---|
79 | USE control_parameters |
---|
80 | USE grid_variables |
---|
81 | USE indices |
---|
82 | USE pegrid |
---|
83 | |
---|
84 | IMPLICIT NONE |
---|
85 | |
---|
86 | CHARACTER (LEN=*) :: range |
---|
87 | |
---|
88 | INTEGER :: i, j, k |
---|
89 | |
---|
90 | |
---|
91 | IF ( range == 'main') THEN |
---|
92 | ! |
---|
93 | !-- Bottom boundary |
---|
94 | IF ( ibc_uv_b == 0 ) THEN |
---|
95 | u(nzb,:,:) = -u(nzb+1,:,:) ! satisfying the Dirichlet condition with |
---|
96 | v(nzb,:,:) = -v(nzb+1,:,:) ! an extra layer below the surface where |
---|
97 | ELSE ! the u and v component change their sign |
---|
98 | u(nzb,:,:) = u(nzb+1,:,:) |
---|
99 | v(nzb,:,:) = v(nzb+1,:,:) |
---|
100 | ENDIF |
---|
101 | DO i = nxl-1, nxr+1 |
---|
102 | DO j = nys-1, nyn+1 |
---|
103 | w(nzb_w_inner(j,i),j,i) = 0.0 |
---|
104 | ENDDO |
---|
105 | ENDDO |
---|
106 | |
---|
107 | ! |
---|
108 | !-- Top boundary |
---|
109 | IF ( ibc_uv_t == 0 ) THEN |
---|
110 | u(nzt+1,:,:) = ug(nzt+1) |
---|
111 | v(nzt+1,:,:) = vg(nzt+1) |
---|
112 | ELSE |
---|
113 | u(nzt+1,:,:) = u(nzt,:,:) |
---|
114 | v(nzt+1,:,:) = v(nzt,:,:) |
---|
115 | ENDIF |
---|
116 | w(nzt:nzt+1,:,:) = 0.0 ! nzt is not a prognostic level (but cf. pres) |
---|
117 | |
---|
118 | ! |
---|
119 | !-- Temperature at bottom boundary |
---|
120 | IF ( ibc_pt_b == 0 ) THEN |
---|
121 | IF ( timestep_scheme(1:5) /= 'runge' ) THEN |
---|
122 | DO i = nxl-1, nxr+1 |
---|
123 | DO j = nys-1, nyn+1 |
---|
124 | pt(nzb_s_inner(j,i),j,i) = pt_m(nzb_s_inner(j,i),j,i) |
---|
125 | ENDDO |
---|
126 | ENDDO |
---|
127 | ELSE |
---|
128 | DO i = nxl-1, nxr+1 |
---|
129 | DO j = nys-1, nyn+1 |
---|
130 | pt(nzb_s_inner(j,i),j,i) = pt_p(nzb_s_inner(j,i),j,i) |
---|
131 | ! pt_m is not used for Runge-Kutta |
---|
132 | ENDDO |
---|
133 | ENDDO |
---|
134 | ENDIF |
---|
135 | ELSE |
---|
136 | DO i = nxl-1, nxr+1 |
---|
137 | DO j = nys-1, nyn+1 |
---|
138 | pt(nzb_s_inner(j,i),j,i) = pt(nzb_s_inner(j,i)+1,j,i) |
---|
139 | ENDDO |
---|
140 | ENDDO |
---|
141 | ENDIF |
---|
142 | |
---|
143 | ! |
---|
144 | !-- Temperature at top boundary |
---|
145 | IF ( ibc_pt_t == 1 ) THEN |
---|
146 | pt(nzt,:,:) = pt(nzt-1,:,:) + bc_pt_t_val * dzu(nzt) |
---|
147 | pt(nzt+1,:,:) = pt(nzt,:,:) + bc_pt_t_val * dzu(nzt+1) |
---|
148 | ENDIF |
---|
149 | |
---|
150 | ! |
---|
151 | !-- Boundary conditions for TKE |
---|
152 | !-- Generally Neumann conditions with de/dz=0 are assumed |
---|
153 | IF ( .NOT. constant_diffusion ) THEN |
---|
154 | DO i = nxl-1, nxr+1 |
---|
155 | DO j = nys-1, nyn+1 |
---|
156 | e(nzb_s_inner(j,i),j,i) = e(nzb_s_inner(j,i)+1,j,i) |
---|
157 | ENDDO |
---|
158 | ENDDO |
---|
159 | e(nzt,:,:) = e(nzt-1,:,:) |
---|
160 | e(nzt+1,:,:) = e(nzt,:,:) |
---|
161 | ENDIF |
---|
162 | |
---|
163 | ! |
---|
164 | !-- Boundary conditions for total water content or scalar, |
---|
165 | !-- bottom and surface boundary (see also temperature) |
---|
166 | IF ( moisture .OR. passive_scalar ) THEN |
---|
167 | ! |
---|
168 | !-- Surface conditions for constant_moisture_flux |
---|
169 | IF ( ibc_q_b == 0 ) THEN |
---|
170 | IF ( timestep_scheme(1:5) /= 'runge' ) THEN |
---|
171 | DO i = nxl-1, nxr+1 |
---|
172 | DO j = nys-1, nyn+1 |
---|
173 | q(nzb_s_inner(j,i),j,i) = q_m(nzb_s_inner(j,i),j,i) |
---|
174 | ENDDO |
---|
175 | ENDDO |
---|
176 | ELSE |
---|
177 | DO i = nxl-1, nxr+1 |
---|
178 | DO j = nys-1, nyn+1 |
---|
179 | q(nzb_s_inner(j,i),j,i) = q_p(nzb_s_inner(j,i),j,i) |
---|
180 | ENDDO ! q_m is not used for Runge-Kutta |
---|
181 | ENDDO |
---|
182 | ENDIF |
---|
183 | ELSE |
---|
184 | DO i = nxl-1, nxr+1 |
---|
185 | DO j = nys-1, nyn+1 |
---|
186 | q(nzb_s_inner(j,i),j,i) = q(nzb_s_inner(j,i)+1,j,i) |
---|
187 | ENDDO |
---|
188 | ENDDO |
---|
189 | ENDIF |
---|
190 | ! |
---|
191 | !-- Top boundary |
---|
192 | q(nzt,:,:) = q(nzt-1,:,:) + bc_q_t_val * dzu(nzt) |
---|
193 | q(nzt+1,:,:) = q(nzt,:,:) + bc_q_t_val * dzu(nzt+1) |
---|
194 | ENDIF |
---|
195 | |
---|
196 | ! |
---|
197 | !-- Lateral boundary conditions at the inflow. Quasi Neumann conditions |
---|
198 | !-- are needed for the wall normal velocity in order to ensure zero |
---|
199 | !-- divergence. Dirichlet conditions are used for all other quantities. |
---|
200 | IF ( inflow_s ) THEN |
---|
201 | v(:,nys,:) = v(:,nys-1,:) |
---|
202 | ELSEIF ( inflow_n ) THEN |
---|
203 | v(:,nyn+vynp,:) = v(:,nyn+vynp+1,:) |
---|
204 | ELSEIF ( inflow_l ) THEN |
---|
205 | u(:,:,nxl) = u(:,:,nxl-1) |
---|
206 | ELSEIF ( inflow_r ) THEN |
---|
207 | u(:,:,nxr+uxrp) = u(:,:,nxr+uxrp+1) |
---|
208 | ENDIF |
---|
209 | |
---|
210 | ! |
---|
211 | !-- Lateral boundary conditions for scalar quantities at the outflow |
---|
212 | IF ( outflow_s ) THEN |
---|
213 | pt(:,nys-1,:) = pt(:,nys,:) |
---|
214 | IF ( .NOT. constant_diffusion ) e(:,nys-1,:) = e(:,nys,:) |
---|
215 | IF ( moisture .OR. passive_scalar ) q(:,nys-1,:) = q(:,nys,:) |
---|
216 | ELSEIF ( outflow_n ) THEN |
---|
217 | pt(:,nyn+1,:) = pt(:,nyn,:) |
---|
218 | IF ( .NOT. constant_diffusion ) e(:,nyn+1,:) = e(:,nyn,:) |
---|
219 | IF ( moisture .OR. passive_scalar ) q(:,nyn+1,:) = q(:,nyn,:) |
---|
220 | ELSEIF ( outflow_l ) THEN |
---|
221 | pt(:,:,nxl-1) = pt(:,:,nxl) |
---|
222 | IF ( .NOT. constant_diffusion ) e(:,:,nxl-1) = e(:,:,nxl) |
---|
223 | IF ( moisture .OR. passive_scalar ) q(:,:,nxl-1) = q(:,:,nxl) |
---|
224 | ELSEIF ( outflow_r ) THEN |
---|
225 | pt(:,:,nxr+1) = pt(:,:,nxr) |
---|
226 | IF ( .NOT. constant_diffusion ) e(:,:,nxr+1) = e(:,:,nxr) |
---|
227 | IF ( moisture .OR. passive_scalar ) q(:,:,nxr+1) = q(:,:,nxr) |
---|
228 | ENDIF |
---|
229 | |
---|
230 | ENDIF |
---|
231 | |
---|
232 | IF ( range == 'outflow_uvw' ) THEN |
---|
233 | ! |
---|
234 | !-- Horizontal boundary conditions for the velocities at the outflow. |
---|
235 | !-- A Neumann condition is used for the wall normal velocity. The vertical |
---|
236 | !-- velocity is assumed as zero and a horizontal average along the wall is |
---|
237 | !-- used for the wall parallel horizontal velocity. The combination of all |
---|
238 | !-- three conditions ensures that the velocity field is free of divergence |
---|
239 | !-- at the outflow (uvmean_outflow_l is calculated in pres). |
---|
240 | IF ( outflow_s ) THEN |
---|
241 | v(:,nys-1,:) = v(:,nys,:) |
---|
242 | w(:,nys-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 u |
---|
246 | #if defined( __parallel ) |
---|
247 | CALL MPI_ALLREDUCE( uvmean_outflow_l, uvmean_outflow, nzt-nzb+2, & |
---|
248 | MPI_REAL, MPI_SUM, comm1dx, ierr ) |
---|
249 | uvmean_outflow = uvmean_outflow / ( nx + 1.0 ) |
---|
250 | #else |
---|
251 | uvmean_outflow = uvmean_outflow_l / ( nx + 1.0 ) |
---|
252 | #endif |
---|
253 | DO k = nzb, nzt+1 |
---|
254 | u(k,nys-1,:) = uvmean_outflow(k) |
---|
255 | ENDDO |
---|
256 | ENDIF |
---|
257 | |
---|
258 | IF ( outflow_n ) THEN |
---|
259 | v(:,nyn+vynp+1,:) = v(:,nyn+vynp,:) |
---|
260 | w(:,nyn+1,:) = 0.0 |
---|
261 | ! |
---|
262 | !-- Compute the mean horizontal wind parallel to and within the outflow |
---|
263 | !-- wall and use this as boundary condition for u |
---|
264 | #if defined( __parallel ) |
---|
265 | CALL MPI_ALLREDUCE( uvmean_outflow_l, uvmean_outflow, nzt-nzb+2, & |
---|
266 | MPI_REAL, MPI_SUM, comm1dx, ierr ) |
---|
267 | uvmean_outflow = uvmean_outflow / ( nx + 1.0 ) |
---|
268 | #else |
---|
269 | uvmean_outflow = uvmean_outflow_l / ( nx + 1.0 ) |
---|
270 | #endif |
---|
271 | DO k = nzb, nzt+1 |
---|
272 | u(k,nyn+1,:) = uvmean_outflow(k) |
---|
273 | ENDDO |
---|
274 | ENDIF |
---|
275 | |
---|
276 | IF ( outflow_l ) THEN |
---|
277 | u(:,:,nxl-1) = u(:,:,nxl) |
---|
278 | w(:,:,nxl-1) = 0.0 |
---|
279 | ! |
---|
280 | !-- Compute the mean horizontal wind parallel to and within the outflow |
---|
281 | !-- wall and use this as boundary condition for v |
---|
282 | #if defined( __parallel ) |
---|
283 | CALL MPI_ALLREDUCE( uvmean_outflow_l, uvmean_outflow, nzt-nzb+2, & |
---|
284 | MPI_REAL, MPI_SUM, comm1dy, ierr ) |
---|
285 | uvmean_outflow = uvmean_outflow / ( ny + 1.0 ) |
---|
286 | #else |
---|
287 | uvmean_outflow = uvmean_outflow_l / ( ny + 1.0 ) |
---|
288 | #endif |
---|
289 | DO k = nzb, nzt+1 |
---|
290 | v(k,:,nxl-1) = uvmean_outflow(k) |
---|
291 | ENDDO |
---|
292 | |
---|
293 | ENDIF |
---|
294 | |
---|
295 | IF ( outflow_r ) THEN |
---|
296 | u(:,:,nxr+uxrp+1) = u(:,:,nxr+uxrp) |
---|
297 | w(:,:,nxr+1) = 0.0 |
---|
298 | ! |
---|
299 | !-- Compute the mean horizontal wind parallel to and within the outflow |
---|
300 | !-- wall and use this as boundary condition for v |
---|
301 | #if defined( __parallel ) |
---|
302 | CALL MPI_ALLREDUCE( uvmean_outflow_l, uvmean_outflow, nzt-nzb+2, & |
---|
303 | MPI_REAL, MPI_SUM, comm1dy, ierr ) |
---|
304 | uvmean_outflow = uvmean_outflow / ( ny + 1.0 ) |
---|
305 | #else |
---|
306 | uvmean_outflow = uvmean_outflow_l / ( ny + 1.0 ) |
---|
307 | #endif |
---|
308 | DO k = nzb, nzt+1 |
---|
309 | v(k,:,nxr+1) = uvmean_outflow(k) |
---|
310 | ENDDO |
---|
311 | ENDIF |
---|
312 | |
---|
313 | ENDIF |
---|
314 | |
---|
315 | |
---|
316 | END SUBROUTINE boundary_conds |
---|