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

Last change on this file since 56 was 56, checked in by raasch, 15 years ago

further checkin of preliminary changes

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