source: palm/trunk/SOURCE/wall_fluxes.f90 @ 329

Last change on this file since 329 was 198, checked in by raasch, 16 years ago

file headers updated for the next release 3.5

  • Property svn:keywords set to Id
File size: 20.8 KB
Line 
1 MODULE wall_fluxes_mod
2!------------------------------------------------------------------------------!
3! Actual revisions:
4! -----------------
5!
6!
7! Former revisions:
8! -----------------
9! $Id: wall_fluxes.f90 198 2008-09-17 08:55:28Z heinze $
10!
11! 187 2008-08-06 16:25:09Z letzel
12! Bugfix: Modification of the evaluation of the vertical turbulent momentum
13! fluxes u'w' and v'w (see prandtl_fluxes), this requires the calculation of
14! us_wall (and vel_total, u_i, v_i, ws) also in wall_fluxes_e.
15! Bugfix: change definition of us_wall from 1D to 2D
16! Bugfix: storage of rifs to rifs_wall in wall_fluxes_e removed
17! Change: add 'minus' sign to fluxes produced by subroutine wall_fluxes_e for
18! consistency with subroutine wall_fluxes
19! Change: Modification of the integrated version of the profile function for
20! momentum for unstable stratification
21!
22! Initial version (2007/03/07)
23!
24! Description:
25! ------------
26! Calculates momentum fluxes at vertical walls assuming Monin-Obukhov
27! similarity.
28! Indices: usvs a=1, vsus b=1, wsvs c1=1, wsus c2=1 (other=0).
29! The all-gridpoint version of wall_fluxes_e is not used so far, because
30! it gives slightly different results from the ij-version for some unknown
31! reason.
32!------------------------------------------------------------------------------!
33    PRIVATE
34    PUBLIC wall_fluxes, wall_fluxes_e
35   
36    INTERFACE wall_fluxes
37       MODULE PROCEDURE wall_fluxes
38       MODULE PROCEDURE wall_fluxes_ij
39    END INTERFACE wall_fluxes
40   
41    INTERFACE wall_fluxes_e
42       MODULE PROCEDURE wall_fluxes_e
43       MODULE PROCEDURE wall_fluxes_e_ij
44    END INTERFACE wall_fluxes_e
45 
46 CONTAINS
47
48!------------------------------------------------------------------------------!
49! Call for all grid points
50!------------------------------------------------------------------------------!
51    SUBROUTINE wall_fluxes( wall_flux, a, b, c1, c2, nzb_uvw_inner, &
52                            nzb_uvw_outer, wall )
53
54       USE arrays_3d
55       USE control_parameters
56       USE grid_variables
57       USE indices
58       USE statistics
59
60       IMPLICIT NONE
61
62       INTEGER ::  i, j, k, wall_index
63
64       INTEGER, DIMENSION(nys-1:nyn+1,nxl-1:nxr+1) ::  nzb_uvw_inner, &
65                                                       nzb_uvw_outer
66       REAL ::  a, b, c1, c2, h1, h2, zp
67       REAL ::  pts, pt_i, rifs, u_i, v_i, us_wall, vel_total, ws, wspts
68
69       REAL, DIMENSION(nys-1:nyn+1,nxl-1:nxr+1)   ::  wall
70       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wall_flux
71
72
73       zp         = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
74       wall_flux  = 0.0
75       wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )
76
77       DO  i = nxl, nxr
78          DO  j = nys, nyn
79
80             IF ( wall(j,i) /= 0.0 )  THEN
81!
82!--             All subsequent variables are computed for the respective
83!--             location where the respective flux is defined.
84                DO  k = nzb_uvw_inner(j,i)+1, nzb_uvw_outer(j,i)
85
86!
87!--                (1) Compute rifs, u_i, v_i, ws, pt' and w'pt'
88                   rifs  = rif_wall(k,j,i,wall_index)
89
90                   u_i   = a * u(k,j,i) + c1 * 0.25 * &
91                           ( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )
92
93                   v_i   = b * v(k,j,i) + c2 * 0.25 * &
94                           ( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) )
95
96                   ws    = ( c1 + c2 ) * w(k,j,i) + 0.25 * (                   &
97                     a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) &
98                   + b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) &
99                                                           )
100                   pt_i  = 0.5 * ( pt(k,j,i) + a *  pt(k,j,i-1) + &
101                                   b * pt(k,j-1,i) + ( c1 + c2 ) * pt(k+1,j,i) )
102
103                   pts   = pt_i - hom(k,1,4,0)
104                   wspts = ws * pts
105
106!
107!--                (2) Compute wall-parallel absolute velocity vel_total
108                   vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 )
109
110!
111!--                (3) Compute wall friction velocity us_wall
112                   IF ( rifs >= 0.0 )  THEN
113
114!
115!--                   Stable stratification (and neutral)
116                      us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +    &
117                                            5.0 * rifs * ( zp - z0(j,i) ) / zp &
118                                                    )
119                   ELSE
120
121!
122!--                   Unstable stratification
123                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
124                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
125
126                      us_wall = kappa * vel_total / (                          &
127                           LOG( zp / z0(j,i) ) -                               &
128                           LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
129                                ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
130                                2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
131                                                    )
132                   ENDIF
133
134!
135!--                (4) Compute zp/L (corresponds to neutral Richardson flux
136!--                    number rifs)
137                   rifs = -1.0 * zp * kappa * g * wspts / ( pt_i * &
138                                                        ( us_wall**3 + 1E-30 ) )
139
140!
141!--                Limit the value range of the Richardson numbers.
142!--                This is necessary for very small velocities (u,w --> 0),
143!--                because the absolute value of rif can then become very
144!--                large, which in consequence would result in very large
145!--                shear stresses and very small momentum fluxes (both are
146!--                generally unrealistic).
147                   IF ( rifs < rif_min )  rifs = rif_min
148                   IF ( rifs > rif_max )  rifs = rif_max
149
150!
151!--                (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
152                   IF ( rifs >= 0.0 )  THEN
153
154!
155!--                   Stable stratification (and neutral)
156                      wall_flux(k,j,i) = kappa *                               &
157                              ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / &
158                              (  LOG( zp / z0(j,i) ) +                         &
159                                 5.0 * rifs * ( zp - z0(j,i) ) / zp            &
160                              )
161                   ELSE
162
163!
164!--                   Unstable stratification
165                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
166                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
167
168                      wall_flux(k,j,i) = kappa *                               &
169                           ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / (  &
170                           LOG( zp / z0(j,i) ) -                               &
171                           LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
172                                ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
173                                2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
174                                                                            )
175                   ENDIF
176                   wall_flux(k,j,i) = -wall_flux(k,j,i) * us_wall
177
178!
179!--                store rifs for next time step
180                   rif_wall(k,j,i,wall_index) = rifs
181
182                ENDDO
183
184             ENDIF
185
186          ENDDO
187       ENDDO
188
189    END SUBROUTINE wall_fluxes
190
191
192
193!------------------------------------------------------------------------------!
194! Call for all grid point i,j
195!------------------------------------------------------------------------------!
196    SUBROUTINE wall_fluxes_ij( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 )
197
198       USE arrays_3d
199       USE control_parameters
200       USE grid_variables
201       USE indices
202       USE statistics
203
204       IMPLICIT NONE
205
206       INTEGER ::  i, j, k, nzb_w, nzt_w, wall_index
207       REAL    ::  a, b, c1, c2, h1, h2, zp
208
209       REAL ::  pts, pt_i, rifs, u_i, v_i, us_wall, vel_total, ws, wspts
210
211       REAL, DIMENSION(nzb:nzt+1) ::  wall_flux
212
213
214       zp         = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
215       wall_flux  = 0.0
216       wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )
217
218!
219!--    All subsequent variables are computed for the respective location where
220!--    the respective flux is defined.
221       DO  k = nzb_w, nzt_w
222
223!
224!--       (1) Compute rifs, u_i, v_i, ws, pt' and w'pt'
225          rifs  = rif_wall(k,j,i,wall_index)
226
227          u_i   = a * u(k,j,i) + c1 * 0.25 * &
228                  ( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) )
229
230          v_i   = b * v(k,j,i) + c2 * 0.25 * &
231                  ( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) )
232
233          ws    = ( c1 + c2 ) * w(k,j,i) + 0.25 * (                            &
234                     a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) &
235                   + b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) &
236                                                  )
237          pt_i  = 0.5 * ( pt(k,j,i) + a *  pt(k,j,i-1) + b * pt(k,j-1,i)  &
238                          + ( c1 + c2 ) * pt(k+1,j,i) )
239
240          pts   = pt_i - hom(k,1,4,0)
241          wspts = ws * pts
242
243!
244!--       (2) Compute wall-parallel absolute velocity vel_total
245          vel_total = SQRT( ws**2 + ( a+c1 ) * u_i**2 + ( b+c2 ) * v_i**2 )
246
247!
248!--       (3) Compute wall friction velocity us_wall
249          IF ( rifs >= 0.0 )  THEN
250
251!
252!--          Stable stratification (and neutral)
253             us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +             &
254                                            5.0 * rifs * ( zp - z0(j,i) ) / zp &
255                                           )
256          ELSE
257
258!
259!--          Unstable stratification
260             h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
261             h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
262
263             us_wall = kappa * vel_total / (                          &
264                  LOG( zp / z0(j,i) ) -                               &
265                  LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
266                       ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
267                       2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
268                                           )
269          ENDIF
270
271!
272!--       (4) Compute zp/L (corresponds to neutral Richardson flux number
273!--           rifs)
274          rifs = -1.0 * zp * kappa * g * wspts / ( pt_i * (us_wall**3 + 1E-30) )
275
276!
277!--       Limit the value range of the Richardson numbers.
278!--       This is necessary for very small velocities (u,w --> 0), because
279!--       the absolute value of rif can then become very large, which in
280!--       consequence would result in very large shear stresses and very
281!--       small momentum fluxes (both are generally unrealistic).
282          IF ( rifs < rif_min )  rifs = rif_min
283          IF ( rifs > rif_max )  rifs = rif_max
284
285!
286!--       (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
287          IF ( rifs >= 0.0 )  THEN
288
289!
290!--          Stable stratification (and neutral)
291             wall_flux(k) = kappa *                                          &
292                            ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / &
293                            (  LOG( zp / z0(j,i) ) +                         &
294                               5.0 * rifs * ( zp - z0(j,i) ) / zp            &
295                            )
296          ELSE
297
298!
299!--          Unstable stratification
300             h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
301             h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
302
303             wall_flux(k) = kappa *                               &
304                  ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / (  &
305                  LOG( zp / z0(j,i) ) -                               &
306                  LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
307                       ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
308                       2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
309                                                                   )
310          ENDIF
311          wall_flux(k) = -wall_flux(k) * us_wall
312
313!
314!--       store rifs for next time step
315          rif_wall(k,j,i,wall_index) = rifs
316
317       ENDDO
318
319    END SUBROUTINE wall_fluxes_ij
320
321
322
323!------------------------------------------------------------------------------!
324! Call for all grid points
325!------------------------------------------------------------------------------!
326    SUBROUTINE wall_fluxes_e( wall_flux, a, b, c1, c2, wall )
327
328!------------------------------------------------------------------------------!
329! Description:
330! ------------
331! Calculates momentum fluxes at vertical walls for routine production_e
332! assuming Monin-Obukhov similarity.
333! Indices: usvs a=1, vsus b=1, wsvs c1=1, wsus c2=1 (other=0).
334!------------------------------------------------------------------------------!
335
336       USE arrays_3d
337       USE control_parameters
338       USE grid_variables
339       USE indices
340       USE statistics
341
342       IMPLICIT NONE
343
344       INTEGER ::  i, j, k, kk, wall_index
345       REAL    ::  a, b, c1, c2, h1, h2, u_i, v_i, us_wall, vel_total, vel_zp, &
346                   ws, zp
347
348       REAL ::  rifs
349
350       REAL, DIMENSION(nys-1:nyn+1,nxl-1:nxr+1)   ::  wall
351       REAL, DIMENSION(nzb:nzt+1,nys:nyn,nxl:nxr) ::  wall_flux
352
353
354       zp         = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
355       wall_flux  = 0.0
356       wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )
357
358       DO  i = nxl, nxr
359          DO  j = nys, nyn
360
361             IF ( wall(j,i) /= 0.0 )  THEN
362!
363!--             All subsequent variables are computed for scalar locations.
364                DO  k = nzb_diff_s_inner(j,i)-1, nzb_diff_s_outer(j,i)-2
365!
366!--                (1) Compute rifs, u_i, v_i, and ws
367                   IF ( k == nzb_diff_s_inner(j,i)-1 )  THEN
368                      kk = nzb_diff_s_inner(j,i)-1
369                   ELSE
370                      kk = k-1
371                   ENDIF
372                   rifs  = 0.5 * ( rif_wall(k,j,i,wall_index) +                &
373                          a * rif_wall(k,j,i+1,1) +  b * rif_wall(k,j+1,i,2) + &
374                          c1 * rif_wall(kk,j,i,3) + c2 * rif_wall(kk,j,i,4)    &
375                                 )
376
377                   u_i   = 0.5 * ( u(k,j,i) + u(k,j,i+1) )
378                   v_i   = 0.5 * ( v(k,j,i) + v(k,j+1,i) )
379                   ws    = 0.5 * ( w(k,j,i) + w(k-1,j,i) )
380!
381!--                (2) Compute wall-parallel absolute velocity vel_total and
382!--                interpolate appropriate velocity component vel_zp.
383                   vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 )
384                   vel_zp = 0.5 * ( a * u_i + b * v_i + (c1+c2) * ws )
385!
386!--                (3) Compute wall friction velocity us_wall
387                   IF ( rifs >= 0.0 )  THEN
388
389!
390!--                   Stable stratification (and neutral)
391                      us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +    &
392                                            5.0 * rifs * ( zp - z0(j,i) ) / zp &
393                                                    )
394                   ELSE
395
396!
397!--                   Unstable stratification
398                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
399                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
400
401                      us_wall = kappa * vel_total / (                          &
402                           LOG( zp / z0(j,i) ) -                               &
403                           LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
404                                ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
405                                2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
406                                                    )
407                   ENDIF
408
409!
410!--                Skip step (4) of wall_fluxes, because here rifs is already
411!--                available from (1)
412!
413!--                (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
414
415                   IF ( rifs >= 0.0 )  THEN
416
417!
418!--                   Stable stratification (and neutral)
419                      wall_flux(k,j,i) = kappa *  vel_zp / &
420                          ( LOG( zp/z0(j,i) ) + 5.0*rifs * ( zp-z0(j,i) ) / zp )
421                   ELSE
422
423!
424!--                   Unstable stratification
425                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
426                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
427
428                      wall_flux(k,j,i) = kappa * vel_zp / (                    &
429                           LOG( zp / z0(j,i) ) -                               &
430                           LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
431                                ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
432                                2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
433                                                          )
434                   ENDIF
435                   wall_flux(k,j,i) = - wall_flux(k,j,i) * us_wall
436
437                ENDDO
438
439             ENDIF
440
441          ENDDO
442       ENDDO
443
444    END SUBROUTINE wall_fluxes_e
445
446
447
448!------------------------------------------------------------------------------!
449! Call for grid point i,j
450!------------------------------------------------------------------------------!
451    SUBROUTINE wall_fluxes_e_ij( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 )
452
453       USE arrays_3d
454       USE control_parameters
455       USE grid_variables
456       USE indices
457       USE statistics
458
459       IMPLICIT NONE
460
461       INTEGER ::  i, j, k, kk, nzb_w, nzt_w, wall_index
462       REAL    ::  a, b, c1, c2, h1, h2, u_i, v_i, us_wall, vel_total, vel_zp, &
463                   ws, zp
464
465       REAL ::  rifs
466
467       REAL, DIMENSION(nzb:nzt+1) ::  wall_flux
468
469
470       zp         = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
471       wall_flux  = 0.0
472       wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 )
473
474!
475!--    All subsequent variables are computed for scalar locations.
476       DO  k = nzb_w, nzt_w
477
478!
479!--       (1) Compute rifs, u_i, v_i, and ws
480          IF ( k == nzb_w )  THEN
481             kk = nzb_w
482          ELSE
483             kk = k-1
484          ENDIF
485          rifs  = 0.5 * ( rif_wall(k,j,i,wall_index) +                         &
486                          a * rif_wall(k,j,i+1,1) +  b * rif_wall(k,j+1,i,2) + &
487                          c1 * rif_wall(kk,j,i,3) + c2 * rif_wall(kk,j,i,4)    &
488                        )
489
490          u_i   = 0.5 * ( u(k,j,i) + u(k,j,i+1) )
491          v_i   = 0.5 * ( v(k,j,i) + v(k,j+1,i) )
492          ws    = 0.5 * ( w(k,j,i) + w(k-1,j,i) )
493!
494!--       (2) Compute wall-parallel absolute velocity vel_total and
495!--       interpolate appropriate velocity component vel_zp.
496          vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 )
497          vel_zp = 0.5 * ( a * u_i + b * v_i + (c1+c2) * ws )
498!
499!--       (3) Compute wall friction velocity us_wall
500          IF ( rifs >= 0.0 )  THEN
501
502!
503!--          Stable stratification (and neutral)
504             us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +             &
505                                            5.0 * rifs * ( zp - z0(j,i) ) / zp &
506                                           )
507          ELSE
508
509!
510!--          Unstable stratification
511             h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
512             h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
513
514             us_wall = kappa * vel_total / (                          &
515                  LOG( zp / z0(j,i) ) -                               &
516                  LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
517                       ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
518                       2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
519                                           )
520          ENDIF
521
522!
523!--       Skip step (4) of wall_fluxes, because here rifs is already
524!--       available from (1)
525!
526!--       (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
527!--       First interpolate the velocity (this is different from
528!--       subroutine wall_fluxes because fluxes in subroutine
529!--       wall_fluxes_e are defined at scalar locations).
530          vel_zp = 0.5 * (       a * ( u(k,j,i) + u(k,j,i+1) ) +  &
531                                 b * ( v(k,j,i) + v(k,j+1,i) ) +  &
532                           (c1+c2) * ( w(k,j,i) + w(k-1,j,i) )    &
533                         )
534
535          IF ( rifs >= 0.0 )  THEN
536
537!
538!--          Stable stratification (and neutral)
539             wall_flux(k) = kappa *  vel_zp / &
540                          ( LOG( zp/z0(j,i) ) + 5.0*rifs * ( zp-z0(j,i) ) / zp )
541          ELSE
542
543!
544!--          Unstable stratification
545             h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
546             h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
547
548             wall_flux(k) = kappa * vel_zp / (                        &
549                  LOG( zp / z0(j,i) ) -                               &
550                  LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
551                       ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
552                       2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
553                                                 )
554          ENDIF
555          wall_flux(k) = - wall_flux(k) * us_wall
556
557       ENDDO
558
559    END SUBROUTINE wall_fluxes_e_ij
560
561 END MODULE wall_fluxes_mod
Note: See TracBrowser for help on using the repository browser.