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

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

preliminary update of further changes, running

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