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

Last change on this file since 104 was 102, checked in by raasch, 17 years ago

preliminary version for coupled runs

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