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

Last change on this file since 108 was 107, checked in by raasch, 17 years ago

further bugfix for non-cyclic BCs

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