source: palm/tags/release-3.8/SOURCE/boundary_conds.f90 @ 4343

Last change on this file since 4343 was 668, checked in by suehring, 11 years ago

last commit documented

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