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

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