1 | SUBROUTINE wall_fluxes( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 ) |
---|
2 | !------------------------------------------------------------------------------! |
---|
3 | ! Actual revisions: |
---|
4 | ! ----------------- |
---|
5 | ! |
---|
6 | ! |
---|
7 | ! Former revisions: |
---|
8 | ! ----------------- |
---|
9 | ! $Id$ |
---|
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 | !------------------------------------------------------------------------------! |
---|
18 | |
---|
19 | USE arrays_3d |
---|
20 | USE control_parameters |
---|
21 | USE grid_variables |
---|
22 | USE indices |
---|
23 | USE statistics |
---|
24 | USE user |
---|
25 | |
---|
26 | IMPLICIT NONE |
---|
27 | |
---|
28 | INTEGER :: i, j, k, nzb_w, nzt_w |
---|
29 | REAL :: a, b, c1, c2, h1, h2, delta_p |
---|
30 | |
---|
31 | REAL :: pts, pt_i, rifs, u_i, v_i, us_wall, vel_total, ws, wspts |
---|
32 | |
---|
33 | REAL, DIMENSION(nzb:nzt+1) :: wall_flux |
---|
34 | |
---|
35 | |
---|
36 | delta_p = 0.5 * ( (a+c1) * dy + (b+c2) * dx ) |
---|
37 | wall_flux = 0.0 |
---|
38 | |
---|
39 | ! |
---|
40 | !-- All subsequent variables are computed for the respective location where |
---|
41 | !-- the relevant variable is defined |
---|
42 | DO k = nzb_w, nzt_w |
---|
43 | ! |
---|
44 | !-- (1) Compute rifs, u_i, v_i, ws, pt' and w'pt' |
---|
45 | rifs = rif_wall(k,j,i,NINT(a+2*b+3*c1+4*c2)) |
---|
46 | u_i = a * u(k,j,i) + & |
---|
47 | c1 * 0.25 * ( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) ) |
---|
48 | v_i = b * v(k,j,i) + & |
---|
49 | c2 * 0.25 * ( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) ) |
---|
50 | ws = ( c1 + c2 ) * w(k,j,i) + & |
---|
51 | a * 0.25 * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) + & |
---|
52 | b * 0.25 * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) |
---|
53 | pt_i = 0.5 * ( pt(k,j,i) + & |
---|
54 | a * pt(k,j,i-1) + b * pt(k,j-1,i) + ( c1 + c2 ) * pt(k+1,j,i) ) |
---|
55 | pts = pt_i - hom(k,1,4,0) |
---|
56 | wspts = ws * pts |
---|
57 | ! |
---|
58 | !-- (2) Compute wall-parallel absolute velocity vel_total |
---|
59 | vel_total = SQRT( ws**2 + ( a+c1 ) * u_i**2 + ( b+c2 ) * v_i**2 ) |
---|
60 | ! |
---|
61 | !-- (3) Compute wall friction velocity us_wall |
---|
62 | IF ( rifs >= 0.0 ) THEN |
---|
63 | ! |
---|
64 | !-- Stable stratification (and neutral) |
---|
65 | us_wall = kappa * vel_total / ( & |
---|
66 | LOG( delta_p / z0(j,i) ) + & |
---|
67 | 5.0 * rifs * ( delta_p - z0(j,i) ) / delta_p & |
---|
68 | ) |
---|
69 | ELSE |
---|
70 | ! |
---|
71 | !-- Unstable stratification |
---|
72 | h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) ) |
---|
73 | h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / delta_p * z0(j,i) ) ) |
---|
74 | ! |
---|
75 | !-- If a borderline case occurs, the formula for stable stratification |
---|
76 | !-- must be used anyway, or else a zero division would occur in the |
---|
77 | !-- argument of the logarithm. |
---|
78 | IF ( h1 == 1.0 .OR. h2 == 1.0 ) THEN |
---|
79 | us_wall = kappa * vel_total / ( & |
---|
80 | LOG( delta_p / z0(j,i) ) + & |
---|
81 | 5.0 * rifs * ( delta_p - z0(j,i) ) / delta_p & |
---|
82 | ) |
---|
83 | ELSE |
---|
84 | us_wall = kappa * vel_total / ( & |
---|
85 | LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) + & |
---|
86 | 2.0 * ( ATAN( h2 ) - ATAN( h1 ) ) & |
---|
87 | ) |
---|
88 | ENDIF |
---|
89 | ENDIF |
---|
90 | ! |
---|
91 | !-- (4) Compute delta_p/L (corresponds to neutral Richardson flux number |
---|
92 | !-- rifs) |
---|
93 | rifs = -1.0 * delta_p * kappa * g * wspts / & |
---|
94 | ( pt_i * ( us_wall**3 + 1E-30 ) ) |
---|
95 | ! |
---|
96 | !-- Limit the value range of the Richardson numbers. |
---|
97 | !-- This is necessary for very small velocities (u,w --> 0), because |
---|
98 | !-- the absolute value of rif can then become very large, which in |
---|
99 | !-- consequence would result in very large shear stresses and very |
---|
100 | !-- small momentum fluxes (both are generally unrealistic). |
---|
101 | IF ( rifs < rif_min ) rifs = rif_min |
---|
102 | IF ( rifs > rif_max ) rifs = rif_max |
---|
103 | ! |
---|
104 | !-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u') |
---|
105 | IF ( rifs >= 0.0 ) THEN |
---|
106 | ! |
---|
107 | !-- Stable stratification (and neutral) |
---|
108 | wall_flux(k) = kappa * & |
---|
109 | ( a * u(k,j,i) + b * v(k,j,i) + (c1 + c2 ) * w(k,j,i) ) / ( & |
---|
110 | LOG( delta_p / z0(j,i) ) + & |
---|
111 | 5.0 * rifs * ( delta_p - z0(j,i) ) / delta_p & |
---|
112 | ) |
---|
113 | ELSE |
---|
114 | ! |
---|
115 | !-- Unstable stratification |
---|
116 | h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) ) |
---|
117 | h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / delta_p * z0(j,i) ) ) |
---|
118 | ! |
---|
119 | !-- If a borderline case occurs, the formula for stable stratification |
---|
120 | !-- must be used anyway, or else a zero division would occur in the |
---|
121 | !-- argument of the logarithm. |
---|
122 | IF ( h1 == 1.0 .OR. h2 == 1.0 ) THEN |
---|
123 | wall_flux(k) = kappa * & |
---|
124 | ( a * u(k,j,i) + b * v(k,j,i) + (c1 + c2 ) * w(k,j,i) ) / ( & |
---|
125 | LOG( delta_p / z0(j,i) ) + & |
---|
126 | 5.0 * rifs * ( delta_p - z0(j,i) ) / delta_p & |
---|
127 | ) |
---|
128 | ELSE |
---|
129 | wall_flux(k) = kappa * & |
---|
130 | ( a * u(k,j,i) + b * v(k,j,i) + (c1 + c2 ) * w(k,j,i) ) / ( & |
---|
131 | LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) + & |
---|
132 | 2.0 * ( ATAN( h2 ) - ATAN( h1 ) ) & |
---|
133 | ) |
---|
134 | ENDIF |
---|
135 | ENDIF |
---|
136 | wall_flux(k) = -wall_flux(k) * ABS( wall_flux(k) ) |
---|
137 | |
---|
138 | ! |
---|
139 | !-- store rifs for next time step |
---|
140 | rif_wall(k,j,i,NINT(a+2*b+3*c1+4*c2)) = rifs |
---|
141 | |
---|
142 | ENDDO |
---|
143 | END SUBROUTINE wall_fluxes |
---|