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

Last change on this file since 75 was 75, checked in by raasch, 14 years ago

preliminary update for changes concerning non-cyclic boundary conditions

  • Property svn:keywords set to Id
File size: 17.1 KB
Line 
1 SUBROUTINE boundary_conds( range )
2
3!------------------------------------------------------------------------------!
4! Actual revisions:
5! -----------------
6! The "main" part sets conditions for time level t+dt instead of level t,
7! outflow boundary conditions changed from Neumann to radiation condition,
8! uxrp, vynp eliminated, moisture renamed humidity
9!
10! Former revisions:
11! -----------------
12! $Id: boundary_conds.f90 75 2007-03-22 09:54:05Z raasch $
13!
14! 19 2007-02-23 04:53:48Z raasch
15! Boundary conditions for e(nzt), pt(nzt), and q(nzt) removed because these
16! gridpoints are now calculated by the prognostic equation,
17! Dirichlet and zero gradient condition for pt established at top boundary
18!
19! RCS Log replace by Id keyword, revision history cleaned up
20!
21! Revision 1.15  2006/02/23 09:54:55  raasch
22! Surface boundary conditions in case of topography: nzb replaced by
23! 2d-k-index-arrays (nzb_w_inner, etc.). Conditions for u and v remain
24! unchanged (still using nzb) because a non-flat topography must use a
25! Prandtl-layer, which don't requires explicit setting of the surface values.
26!
27! Revision 1.1  1997/09/12 06:21:34  raasch
28! Initial revision
29!
30!
31! Description:
32! ------------
33! Boundary conditions for the prognostic quantities (range='main').
34! In case of non-cyclic lateral boundaries the conditions for velocities at
35! the outflow are set after the pressure solver has been called (range=
36! 'outflow_uvw').
37! One additional bottom boundary condition is applied for the TKE (=(u*)**2)
38! in prandtl_fluxes. The cyclic lateral boundary conditions are implicitly
39! handled in routine exchange_horiz. Pressure boundary conditions are
40! explicitly set in routines pres, poisfft, poismg and sor.
41!------------------------------------------------------------------------------!
42
43    USE arrays_3d
44    USE control_parameters
45    USE grid_variables
46    USE indices
47    USE pegrid
48
49    IMPLICIT NONE
50
51    CHARACTER (LEN=*) ::  range
52
53    INTEGER ::  i, j, k
54
55    REAL    ::  c_max, c_u, c_v, c_w, denom
56
57
58    IF ( range == 'main')  THEN
59!
60!--    Bottom boundary
61       IF ( ibc_uv_b == 0 )  THEN
62!
63!--       Satisfying the Dirichlet condition with an extra layer below the
64!--       surface where the u and v component change their sign
65          u_p(nzb,:,:) = -u_p(nzb+1,:,:)
66          v_p(nzb,:,:) = -v_p(nzb+1,:,:)
67       ELSE
68          u_p(nzb,:,:) = u_p(nzb+1,:,:)
69          v_p(nzb,:,:) = v_p(nzb+1,:,:)
70       ENDIF
71       DO  i = nxl-1, nxr+1
72          DO  j = nys-1, nyn+1
73             w_p(nzb_w_inner(j,i),j,i) = 0.0
74          ENDDO
75       ENDDO
76
77!
78!--    Top boundary
79       IF ( ibc_uv_t == 0 )  THEN
80          u_p(nzt+1,:,:) = ug(nzt+1)
81          v_p(nzt+1,:,:) = vg(nzt+1)
82       ELSE
83          u_p(nzt+1,:,:) = u_p(nzt,:,:)
84          v_p(nzt+1,:,:) = v_p(nzt,:,:)
85       ENDIF
86       w_p(nzt:nzt+1,:,:) = 0.0  ! nzt is not a prognostic level (but cf. pres)
87
88!
89!--    Temperature at bottom boundary
90       IF ( ibc_pt_b == 0 )  THEN
91          DO  i = nxl-1, nxr+1
92             DO  j = nys-1, nyn+1
93                pt_p(nzb_s_inner(j,i),j,i) = pt(nzb_s_inner(j,i),j,i)
94             ENDDO
95          ENDDO
96       ELSE
97          DO  i = nxl-1, nxr+1
98             DO  j = nys-1, nyn+1
99                pt_p(nzb_s_inner(j,i),j,i) = pt_p(nzb_s_inner(j,i)+1,j,i)
100             ENDDO
101          ENDDO
102       ENDIF
103
104!
105!--    Temperature at top boundary
106       IF ( ibc_pt_t == 0 )  THEN
107          pt_p(nzt+1,:,:) = pt(nzt+1,:,:)
108       ELSEIF ( ibc_pt_t == 1 )  THEN
109          pt_p(nzt+1,:,:) = pt_p(nzt,:,:)
110       ELSEIF ( ibc_pt_t == 2 )  THEN
111          pt_p(nzt+1,:,:) = pt_p(nzt,:,:)   + bc_pt_t_val * dzu(nzt+1)
112       ENDIF
113
114!
115!--    Boundary conditions for TKE
116!--    Generally Neumann conditions with de/dz=0 are assumed
117       IF ( .NOT. constant_diffusion )  THEN
118          DO  i = nxl-1, nxr+1
119             DO  j = nys-1, nyn+1
120                e_p(nzb_s_inner(j,i),j,i) = e_p(nzb_s_inner(j,i)+1,j,i)
121             ENDDO
122          ENDDO
123          e_p(nzt+1,:,:) = e_p(nzt,:,:)
124       ENDIF
125
126!
127!--    Boundary conditions for total water content or scalar,
128!--    bottom and surface boundary (see also temperature)
129       IF ( humidity  .OR.  passive_scalar )  THEN
130!
131!--       Surface conditions for constant_humidity_flux
132          IF ( ibc_q_b == 0 ) THEN
133             DO  i = nxl-1, nxr+1
134                DO  j = nys-1, nyn+1
135                   q_p(nzb_s_inner(j,i),j,i) = q(nzb_s_inner(j,i),j,i)
136                ENDDO
137             ENDDO
138          ELSE
139             DO  i = nxl-1, nxr+1
140                DO  j = nys-1, nyn+1
141                   q_p(nzb_s_inner(j,i),j,i) = q_p(nzb_s_inner(j,i)+1,j,i)
142                ENDDO
143             ENDDO
144          ENDIF
145!
146!--       Top boundary
147          q_p(nzt+1,:,:) = q_p(nzt,:,:)   + bc_q_t_val * dzu(nzt+1)
148       ENDIF
149
150!
151!--    Lateral boundary conditions at the inflow. Quasi Neumann conditions
152!--    are needed for the wall normal velocity in order to ensure zero
153!--    divergence. Dirichlet conditions are used for all other quantities.
154       IF ( inflow_s )  THEN
155          v_p(:,nys,:) = v_p(:,nys-1,:)
156       ELSEIF ( inflow_n ) THEN
157          v_p(:,nyn,:) = v_p(:,nyn+1,:)
158       ELSEIF ( inflow_l ) THEN
159          u_p(:,:,nxl) = u_p(:,:,nxl-1)
160       ELSEIF ( inflow_r ) THEN
161          u_p(:,:,nxr) = u_p(:,:,nxr+1)
162       ENDIF
163
164!
165!--    Lateral boundary conditions for scalar quantities at the outflow
166       IF ( outflow_s )  THEN
167          pt_p(:,nys-1,:)     = pt_p(:,nys,:)
168          IF ( .NOT. constant_diffusion     )  e_p(:,nys-1,:) = e_p(:,nys,:)
169          IF ( humidity .OR. passive_scalar )  q_p(:,nys-1,:) = q_p(:,nys,:)
170       ELSEIF ( outflow_n )  THEN
171          pt_p(:,nyn+1,:)     = pt_p(:,nyn,:)
172          IF ( .NOT. constant_diffusion     )  e_p(:,nyn+1,:) = e_p(:,nyn,:)
173          IF ( humidity .OR. passive_scalar )  q_p(:,nyn+1,:) = q_p(:,nyn,:)
174       ELSEIF ( outflow_l )  THEN
175          pt_p(:,:,nxl-1)     = pt_p(:,:,nxl)
176          IF ( .NOT. constant_diffusion     )  e_p(:,:,nxl-1) = e_p(:,:,nxl)
177          IF ( humidity .OR. passive_scalar )  q_p(:,:,nxl-1) = q_p(:,:,nxl)
178       ELSEIF ( outflow_r )  THEN
179          pt_p(:,:,nxr+1)     = pt_p(:,:,nxr)
180          IF ( .NOT. constant_diffusion     )  e_p(:,:,nxr+1) = e_p(:,:,nxr)
181          IF ( humidity .OR. passive_scalar )  q_p(:,:,nxr+1) = q_p(:,:,nxr)
182       ENDIF
183
184    ENDIF
185
186!
187!-- Radiation boundary condition for the velocities at the respective outflow
188    IF ( outflow_s  .AND.                                                 &
189         intermediate_timestep_count == intermediate_timestep_count_max ) &
190    THEN
191
192       c_max = dy / dt_3d
193
194       DO i = nxl-1, nxr+1
195          DO k = nzb+1, nzt+1
196
197!
198!--          First calculate the phase speeds for u,v, and w
199             denom = u_m_s(k,0,i) - u_m_s(k,1,i)
200
201             IF ( denom /= 0.0 )  THEN
202                c_u = -c_max * ( u(k,0,i) - u_m_s(k,0,i) ) / denom
203                IF ( c_u < 0.0 )  THEN
204                   c_u = 0.0
205                ELSEIF ( c_u > c_max )  THEN
206                   c_u = c_max
207                ENDIF
208             ELSE
209                c_u = c_max
210             ENDIF
211             denom = v_m_s(k,0,i) - v_m_s(k,1,i)
212
213             IF ( denom /= 0.0 )  THEN
214                c_v = -c_max * ( v(k,0,i) - v_m_s(k,0,i) ) / denom
215                IF ( c_v < 0.0 )  THEN
216                   c_v = 0.0
217                ELSEIF ( c_v > c_max )  THEN
218                   c_v = c_max
219                ENDIF
220             ELSE
221                c_v = c_max
222             ENDIF
223
224             denom = w_m_s(k,0,i) - w_m_s(k,1,i)
225
226             IF ( denom /= 0.0 )  THEN
227                c_w = -c_max * ( w(k,0,i) - w_m_s(k,0,i) ) / denom
228                IF ( c_w < 0.0 )  THEN
229                   c_w = 0.0
230                ELSEIF ( c_w > c_max )  THEN
231                   c_w = c_max
232                ENDIF
233             ELSE
234                c_w = c_max
235             ENDIF
236
237!
238!--          Calculate the new velocities
239             u_p(k,-1,i) = u(k,-1,i) + dt_3d * c_u * &
240                                       ( u(k,-1,i) - u(k,0,i) ) * ddy
241
242             v_p(k,-1,i) = v(k,-1,i) + dt_3d * c_v * &
243                                       ( v(k,-1,i) - v_m_s(k,0,i) ) * ddy
244
245             w_p(k,-1,i) = w(k,-1,i) + dt_3d * c_w * &
246                                       ( w(k,-1,i) - w(k,0,i) ) * ddy
247
248!
249!--          Save old timelevels for the next timestep
250             u_m_s(k,:,i) = u(k,-1:1,i)
251             v_m_s(k,:,i) = v(k,-1:1,i)
252             w_m_s(k,:,i) = w(k,-1:1,i)
253
254          ENDDO
255       ENDDO
256
257!
258!--    Bottom boundary at the outflow
259       IF ( ibc_uv_b == 0 )  THEN
260          u_p(nzb,-1,:) = -u_p(nzb+1,-1,:) 
261          v_p(nzb,-1,:) = -v_p(nzb+1,-1,:) 
262       ELSE                   
263          u_p(nzb,-1,:) =  u_p(nzb+1,-1,:)
264          v_p(nzb,-1,:) =  v_p(nzb+1,-1,:)
265       ENDIF
266       w_p(nzb,ny+1,:) = 0.0
267
268!
269!--    Top boundary at the outflow
270       IF ( ibc_uv_t == 0 )  THEN
271          u_p(nzt+1,-1,:) = ug(nzt+1)
272          v_p(nzt+1,-1,:) = vg(nzt+1)
273       ELSE
274          u_p(nzt+1,-1,:) = u(nzt,-1,:)
275          v_p(nzt+1,-1,:) = v(nzt,-1,:)
276       ENDIF
277       w_p(nzt:nzt+1,-1,:) = 0.0
278
279    ENDIF
280
281    IF ( outflow_n  .AND.                                                 &
282         intermediate_timestep_count == intermediate_timestep_count_max ) &
283    THEN
284
285       c_max = dy / dt_3d
286
287       DO i = nxl-1, nxr+1
288          DO k = nzb+1, nzt+1
289
290!
291!--          First calculate the phase speeds for u,v, and w
292             denom = u_m_n(k,ny,i) - u_m_n(k,ny-1,i)
293
294             IF ( denom /= 0.0 )  THEN
295                c_u = -c_max * ( u(k,ny,i) - u_m_n(k,ny,i) ) / denom
296                IF ( c_u < 0.0 )  THEN
297                   c_u = 0.0
298                ELSEIF ( c_u > c_max )  THEN
299                   c_u = c_max
300                ENDIF
301             ELSE
302                c_u = c_max
303             ENDIF
304
305             denom = v_m_n(k,ny,i) - v_m_n(k,ny-1,i)
306
307             IF ( denom /= 0.0 )  THEN
308                c_v = -c_max * ( v(k,ny,i) - v_m_n(k,ny,i) ) / denom
309                IF ( c_v < 0.0 )  THEN
310                   c_v = 0.0
311                ELSEIF ( c_v > c_max )  THEN
312                   c_v = c_max
313                ENDIF
314             ELSE
315                c_v = c_max
316             ENDIF
317
318             denom = w_m_n(k,ny,i) - w_m_n(k,ny-1,i)
319
320             IF ( denom /= 0.0 )  THEN
321                c_w = -c_max * ( w(k,ny,i) - w_m_n(k,ny,i) ) / denom
322                IF ( c_w < 0.0 )  THEN
323                   c_w = 0.0
324                ELSEIF ( c_w > c_max )  THEN
325                   c_w = c_max
326                ENDIF
327             ELSE
328                c_w = c_max
329             ENDIF
330
331!
332!--          Calculate the new velocities
333             u_p(k,ny+1,i) = u(k,ny+1,i) - dt_3d * c_u * &
334                                           ( u(k,ny+1,i) - u(k,ny,i) ) * ddy
335
336             v_p(k,ny+1,i) = v(k,ny+1,i) - dt_3d * c_v * &
337                                           ( v(k,ny+1,i) - v(k,ny,i) ) * ddy
338
339             w_p(k,ny+1,i) = w(k,ny+1,i) - dt_3d * c_w * &
340                                           ( w(k,ny+1,i) - w(k,ny,i) ) * ddy
341
342!
343!--          Swap timelevels for the next timestep
344             u_m_n(k,:,i) = u(k,ny-1:ny+1,i)
345             v_m_n(k,:,i) = v(k,ny-1:ny+1,i)
346             w_m_n(k,:,i) = w(k,ny-1:ny+1,i)
347
348          ENDDO
349       ENDDO
350
351!
352!--    Bottom boundary at the outflow
353       IF ( ibc_uv_b == 0 )  THEN
354          u_p(nzb,ny+1,:) = -u_p(nzb+1,ny+1,:) 
355          v_p(nzb,ny+1,:) = -v_p(nzb+1,ny+1,:) 
356       ELSE                   
357          u_p(nzb,ny+1,:) =  u_p(nzb+1,ny+1,:)
358          v_p(nzb,ny+1,:) =  v_p(nzb+1,ny+1,:)
359       ENDIF
360       w_p(nzb,ny+1,:) = 0.0
361
362!
363!--    Top boundary at the outflow
364       IF ( ibc_uv_t == 0 )  THEN
365          u_p(nzt+1,ny+1,:) = ug(nzt+1)
366          v_p(nzt+1,ny+1,:) = vg(nzt+1)
367       ELSE
368          u_p(nzt+1,ny+1,:) = u_p(nzt,nyn+1,:)
369          v_p(nzt+1,ny+1,:) = v_p(nzt,nyn+1,:)
370       ENDIF
371       w_p(nzt:nzt+1,ny+1,:) = 0.0
372
373    ENDIF
374
375    IF ( outflow_l  .AND.                                                 &
376         intermediate_timestep_count == intermediate_timestep_count_max ) &
377    THEN
378
379       c_max = dx / dt_3d
380
381       DO j = nys-1, nyn+1
382          DO k = nzb+1, nzt+1
383
384!
385!--          First calculate the phase speeds for u,v, and w
386             denom = u_m_l(k,j,0) - u_m_l(k,j,1)
387
388             IF ( denom /= 0.0 )  THEN
389                c_u = -c_max * ( u(k,j,0) - u_m_r(k,j,0) ) / denom
390                IF ( c_u > 0.0 )  THEN
391                   c_u = 0.0
392                ELSEIF ( c_u < -c_max )  THEN
393                   c_u = -c_max
394                ENDIF
395             ELSE
396                c_u = -c_max
397             ENDIF
398
399             denom = v_m_l(k,j,0) - v_m_l(k,j,1)
400
401             IF ( denom /= 0.0 )  THEN
402                c_v = -c_max * ( v(k,j,0) - v_m_l(k,j,0) ) / denom
403                IF ( c_v < 0.0 )  THEN
404                   c_v = 0.0
405                ELSEIF ( c_v > c_max )  THEN
406                   c_v = c_max
407                ENDIF
408             ELSE
409                c_v = c_max
410             ENDIF
411
412             denom = w_m_l(k,j,0) - w_m_l(k,j,1)
413
414             IF ( denom /= 0.0 )  THEN
415                c_w = -c_max * ( w(k,j,0) - w_m_l(k,j,0) ) / denom
416                IF ( c_w < 0.0 )  THEN
417                   c_w = 0.0
418                ELSEIF ( c_w > c_max )  THEN
419                   c_w = c_max
420                ENDIF
421             ELSE
422                c_w = c_max
423             ENDIF
424
425!
426!--          Calculate the new velocities
427             u_p(k,j,-1) = u(k,j,-1) + dt_3d * c_u * &
428                                       ( u(k,j,-1) - u(k,j,0) ) * ddx
429
430             v_p(k,j,-1) = v(k,j,-1) + dt_3d * c_v * &
431                                       ( v(k,j,-1) - v(k,j,0) ) * ddx
432
433             w_p(k,j,-1) = w(k,j,-1) + dt_3d * c_w * &
434                                       ( w(k,j,-1) - w(k,j,0) ) * ddx
435
436!
437!--          Swap timelevels for the next timestep
438             u_m_l(k,j,:) = u(k,j,-1:1)
439             v_m_l(k,j,:) = v(k,j,-1:1)
440             w_m_l(k,j,:) = w(k,j,-1:1)
441
442          ENDDO
443       ENDDO
444
445!
446!--    Bottom boundary at the outflow
447       IF ( ibc_uv_b == 0 )  THEN
448          u_p(nzb,:,-1) = -u_p(nzb+1,:,-1) 
449          v_p(nzb,:,-1) = -v_p(nzb+1,:,-1) 
450       ELSE                   
451          u_p(nzb,:,-1) =  u_p(nzb+1,:,-1)
452          v_p(nzb,:,-1) =  v_p(nzb+1,:,-1)
453       ENDIF
454       w_p(nzb,:,-1) = 0.0
455
456!
457!--    Top boundary at the outflow
458       IF ( ibc_uv_t == 0 )  THEN
459          u_p(nzt+1,:,-1) = ug(nzt+1)
460          v_p(nzt+1,:,-1) = vg(nzt+1)
461       ELSE
462          u_p(nzt+1,:,-1) = u_p(nzt,:,-1)
463          v_p(nzt+1,:,-1) = v_p(nzt,:,-1)
464       ENDIF
465       w_p(nzt:nzt+1,:,-1) = 0.0
466
467    ENDIF
468
469    IF ( outflow_r  .AND.                                                 &
470         intermediate_timestep_count == intermediate_timestep_count_max ) &
471    THEN
472
473       c_max = dx / dt_3d
474
475       DO j = nys-1, nyn+1
476          DO k = nzb+1, nzt+1
477
478!
479!--          First calculate the phase speeds for u,v, and w
480             denom = u_m_r(k,j,nx) - u_m_r(k,j,nx-1)
481
482             IF ( denom /= 0.0 )  THEN
483                c_u = -c_max * ( u(k,j,nx) - u_m_r(k,j,nx) ) / denom
484                IF ( c_u < 0.0 )  THEN
485                   c_u = 0.0
486                ELSEIF ( c_u > c_max )  THEN
487                   c_u = c_max
488                ENDIF
489             ELSE
490                c_u = c_max
491             ENDIF
492
493             denom = v_m_r(k,j,nx) - v_m_r(k,j,nx-1)
494
495             IF ( denom /= 0.0 )  THEN
496                c_v = -c_max * ( v(k,j,nx) - v_m_r(k,j,nx) ) / denom
497                IF ( c_v < 0.0 )  THEN
498                   c_v = 0.0
499                ELSEIF ( c_v > c_max )  THEN
500                   c_v = c_max
501                ENDIF
502             ELSE
503                c_v = c_max
504             ENDIF
505
506             denom = w_m_r(k,j,nx) - w_m_r(k,j,nx-1)
507
508             IF ( denom /= 0.0 )  THEN
509                c_w = -c_max * ( w(k,j,nx) - w_m_r(k,j,nx) ) / denom
510                IF ( c_w < 0.0 )  THEN
511                   c_w = 0.0
512                ELSEIF ( c_w > c_max )  THEN
513                   c_w = c_max
514                ENDIF
515             ELSE
516                c_w = c_max
517             ENDIF
518
519!
520!--          Calculate the new velocities
521             u_p(k,j,nx+1) = u(k,j,nx+1) - dt_3d * c_u * &
522                                           ( u(k,j,nx+1) - u(k,j,nx) ) * ddx
523
524             v_p(k,j,nx+1) = v(k,j,nx+1) - dt_3d * c_v * &
525                                           ( v(k,j,nx+1) - v(k,j,nx) ) * ddx
526
527             w_p(k,j,nx+1) = w(k,j,nx+1) - dt_3d * c_w * &
528                                           ( w(k,j,nx+1) - w(k,j,nx) ) * ddx
529
530!
531!--          Swap timelevels for the next timestep
532             u_m_r(k,j,:) = u(k,j,nx-1:nx+1)
533             v_m_r(k,j,:) = v(k,j,nx-1:nx+1)
534             w_m_r(k,j,:) = w(k,j,nx-1:nx+1)
535
536          ENDDO
537       ENDDO
538
539!
540!--    Bottom boundary at the outflow
541       IF ( ibc_uv_b == 0 )  THEN
542          u_p(nzb,:,nx+1) = -u_p(nzb+1,:,nx+1) 
543          v_p(nzb,:,nx+1) = -v_p(nzb+1,:,nx+1) 
544       ELSE                   
545          u_p(nzb,:,nx+1) =  u_p(nzb+1,:,nx+1)
546          v_p(nzb,:,nx+1) =  v_p(nzb+1,:,nx+1)
547       ENDIF
548       w_p(nzb,:,nx+1) = 0.0
549
550!
551!--    Top boundary at the outflow
552       IF ( ibc_uv_t == 0 )  THEN
553          u_p(nzt+1,:,nx+1) = ug(nzt+1)
554          v_p(nzt+1,:,nx+1) = vg(nzt+1)
555       ELSE
556          u_p(nzt+1,:,nx+1) = u_p(nzt,:,nx+1)
557          v_p(nzt+1,:,nx+1) = v_p(nzt,:,nx+1)
558       ENDIF
559       w(nzt:nzt+1,:,nx+1) = 0.0
560
561    ENDIF
562
563 
564 END SUBROUTINE boundary_conds
Note: See TracBrowser for help on using the repository browser.