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

Last change on this file since 97 was 95, checked in by raasch, 17 years ago

further preliminary uncomplete changes for ocean version

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