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

Last change on this file since 768 was 768, checked in by raasch, 11 years ago

last commit documented

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