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

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