[52] | 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: wall_fluxes.f90 55 2007-03-08 04:12:41Z 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 | !------------------------------------------------------------------------------! |
---|
| 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 | |
---|
[53] | 28 | INTEGER :: i, j, k, nzb_w, nzt_w, wall_index |
---|
| 29 | REAL :: a, b, c1, c2, h1, h2, zp |
---|
[52] | 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 | |
---|
[53] | 36 | zp = 0.5 * ( (a+c1) * dy + (b+c2) * dx ) |
---|
| 37 | wall_flux = 0.0 |
---|
| 38 | wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 ) |
---|
[52] | 39 | |
---|
| 40 | ! |
---|
| 41 | !-- All subsequent variables are computed for the respective location where |
---|
| 42 | !-- the relevant variable is defined |
---|
| 43 | DO k = nzb_w, nzt_w |
---|
[53] | 44 | |
---|
[52] | 45 | ! |
---|
| 46 | !-- (1) Compute rifs, u_i, v_i, ws, pt' and w'pt' |
---|
[53] | 47 | rifs = rif_wall(k,j,i,wall_index) |
---|
| 48 | |
---|
| 49 | u_i = a * u(k,j,i) + c1 * 0.25 * & |
---|
| 50 | ( u(k+1,j,i+1) + u(k+1,j,i) + u(k,j,i+1) + u(k,j,i) ) |
---|
| 51 | |
---|
| 52 | v_i = b * v(k,j,i) + c2 * 0.25 * & |
---|
| 53 | ( v(k+1,j+1,i) + v(k+1,j,i) + v(k,j+1,i) + v(k,j,i) ) |
---|
| 54 | |
---|
| 55 | ws = ( c1 + c2 ) * w(k,j,i) + 0.25 * ( & |
---|
| 56 | a * ( w(k-1,j,i-1) + w(k-1,j,i) + w(k,j,i-1) + w(k,j,i) ) & |
---|
| 57 | + b * ( w(k-1,j-1,i) + w(k-1,j,i) + w(k,j-1,i) + w(k,j,i) ) & |
---|
| 58 | ) |
---|
| 59 | pt_i = 0.5 * ( pt(k,j,i) + a * pt(k,j,i-1) + b * pt(k,j-1,i) & |
---|
| 60 | + ( c1 + c2 ) * pt(k+1,j,i) ) |
---|
| 61 | |
---|
| 62 | pts = pt_i - hom(k,1,4,0) |
---|
| 63 | wspts = ws * pts |
---|
| 64 | |
---|
[52] | 65 | ! |
---|
| 66 | !-- (2) Compute wall-parallel absolute velocity vel_total |
---|
| 67 | vel_total = SQRT( ws**2 + ( a+c1 ) * u_i**2 + ( b+c2 ) * v_i**2 ) |
---|
[53] | 68 | |
---|
[52] | 69 | ! |
---|
| 70 | !-- (3) Compute wall friction velocity us_wall |
---|
| 71 | IF ( rifs >= 0.0 ) THEN |
---|
[53] | 72 | |
---|
[52] | 73 | ! |
---|
| 74 | !-- Stable stratification (and neutral) |
---|
[53] | 75 | us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + & |
---|
| 76 | 5.0 * rifs * ( zp - z0(j,i) ) / zp & |
---|
[52] | 77 | ) |
---|
| 78 | ELSE |
---|
[53] | 79 | |
---|
[52] | 80 | ! |
---|
| 81 | !-- Unstable stratification |
---|
| 82 | h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) ) |
---|
[53] | 83 | h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) ) |
---|
| 84 | |
---|
[52] | 85 | ! |
---|
| 86 | !-- If a borderline case occurs, the formula for stable stratification |
---|
| 87 | !-- must be used anyway, or else a zero division would occur in the |
---|
| 88 | !-- argument of the logarithm. |
---|
| 89 | IF ( h1 == 1.0 .OR. h2 == 1.0 ) THEN |
---|
[53] | 90 | us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) + & |
---|
| 91 | 5.0 * rifs * ( zp - z0(j,i) ) / zp & |
---|
| 92 | ) |
---|
[52] | 93 | ELSE |
---|
[53] | 94 | us_wall = kappa * vel_total / ( & |
---|
[52] | 95 | LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) + & |
---|
| 96 | 2.0 * ( ATAN( h2 ) - ATAN( h1 ) ) & |
---|
[53] | 97 | ) |
---|
[52] | 98 | ENDIF |
---|
[53] | 99 | |
---|
[52] | 100 | ENDIF |
---|
[53] | 101 | |
---|
[52] | 102 | ! |
---|
[53] | 103 | !-- (4) Compute zp/L (corresponds to neutral Richardson flux number |
---|
[52] | 104 | !-- rifs) |
---|
[53] | 105 | rifs = -1.0 * zp * kappa * g * wspts / ( pt_i * ( us_wall**3 + 1E-30 ) ) |
---|
| 106 | |
---|
[52] | 107 | ! |
---|
| 108 | !-- Limit the value range of the Richardson numbers. |
---|
| 109 | !-- This is necessary for very small velocities (u,w --> 0), because |
---|
| 110 | !-- the absolute value of rif can then become very large, which in |
---|
| 111 | !-- consequence would result in very large shear stresses and very |
---|
| 112 | !-- small momentum fluxes (both are generally unrealistic). |
---|
| 113 | IF ( rifs < rif_min ) rifs = rif_min |
---|
| 114 | IF ( rifs > rif_max ) rifs = rif_max |
---|
[53] | 115 | |
---|
[52] | 116 | ! |
---|
| 117 | !-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u') |
---|
| 118 | IF ( rifs >= 0.0 ) THEN |
---|
[53] | 119 | |
---|
[52] | 120 | ! |
---|
| 121 | !-- Stable stratification (and neutral) |
---|
[53] | 122 | wall_flux(k) = kappa * & |
---|
| 123 | ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / & |
---|
| 124 | ( LOG( zp / z0(j,i) ) + & |
---|
| 125 | 5.0 * rifs * ( zp - z0(j,i) ) / zp & |
---|
| 126 | ) |
---|
[52] | 127 | ELSE |
---|
[53] | 128 | |
---|
[52] | 129 | ! |
---|
| 130 | !-- Unstable stratification |
---|
| 131 | h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) ) |
---|
[53] | 132 | h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) ) |
---|
| 133 | |
---|
[52] | 134 | ! |
---|
| 135 | !-- If a borderline case occurs, the formula for stable stratification |
---|
| 136 | !-- must be used anyway, or else a zero division would occur in the |
---|
| 137 | !-- argument of the logarithm. |
---|
| 138 | IF ( h1 == 1.0 .OR. h2 == 1.0 ) THEN |
---|
[53] | 139 | wall_flux(k) = kappa * & |
---|
| 140 | ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / & |
---|
| 141 | ( LOG( zp / z0(j,i) ) + & |
---|
| 142 | 5.0 * rifs * ( zp - z0(j,i) ) / zp & |
---|
| 143 | ) |
---|
[52] | 144 | ELSE |
---|
[53] | 145 | wall_flux(k) = kappa * & |
---|
| 146 | ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / & |
---|
| 147 | ( LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) & |
---|
| 148 | + 2.0 * ( ATAN( h2 ) - ATAN( h1 ) ) & |
---|
| 149 | ) |
---|
[52] | 150 | ENDIF |
---|
[53] | 151 | |
---|
[52] | 152 | ENDIF |
---|
| 153 | wall_flux(k) = -wall_flux(k) * ABS( wall_flux(k) ) |
---|
| 154 | |
---|
| 155 | ! |
---|
| 156 | !-- store rifs for next time step |
---|
[53] | 157 | rif_wall(k,j,i,wall_index) = rifs |
---|
[52] | 158 | |
---|
| 159 | ENDDO |
---|
[53] | 160 | |
---|
[52] | 161 | END SUBROUTINE wall_fluxes |
---|
[53] | 162 | |
---|
| 163 | |
---|
| 164 | |
---|
| 165 | SUBROUTINE wall_fluxes_e( i, j, nzb_w, nzt_w, wall_flux, a, b, c1, c2 ) |
---|
| 166 | !------------------------------------------------------------------------------! |
---|
| 167 | ! Actual revisions: |
---|
| 168 | ! ----------------- |
---|
| 169 | ! |
---|
| 170 | ! |
---|
| 171 | ! Former revisions: |
---|
| 172 | ! ----------------- |
---|
| 173 | ! Initial version (2007/03/07) |
---|
| 174 | ! |
---|
| 175 | ! Description: |
---|
| 176 | ! ------------ |
---|
| 177 | ! Calculates momentum fluxes at vertical walls for routine production_e |
---|
| 178 | ! assuming Monin-Obukhov similarity. |
---|
| 179 | ! Indices: usvs a=1, vsus b=1, wsvs c1=1, wsus c2=1 (other=0). |
---|
| 180 | !------------------------------------------------------------------------------! |
---|
| 181 | |
---|
| 182 | USE arrays_3d |
---|
| 183 | USE control_parameters |
---|
| 184 | USE grid_variables |
---|
| 185 | USE indices |
---|
| 186 | USE statistics |
---|
| 187 | USE user |
---|
| 188 | |
---|
| 189 | IMPLICIT NONE |
---|
| 190 | |
---|
| 191 | INTEGER :: i, j, k, kk, nzb_w, nzt_w, wall_index |
---|
[55] | 192 | REAL :: a, b, c1, c2, h1, h2, vel_zp, zp |
---|
[53] | 193 | |
---|
| 194 | REAL :: rifs |
---|
| 195 | |
---|
| 196 | REAL, DIMENSION(nzb:nzt+1) :: wall_flux |
---|
| 197 | |
---|
| 198 | |
---|
| 199 | zp = 0.5 * ( (a+c1) * dy + (b+c2) * dx ) |
---|
| 200 | wall_flux = 0.0 |
---|
| 201 | wall_index = NINT( a+ 2*b + 3*c1 + 4*c2 ) |
---|
| 202 | |
---|
| 203 | ! |
---|
| 204 | !-- All subsequent variables are computed for the respective location where |
---|
| 205 | !-- the relevant variable is defined |
---|
| 206 | DO k = nzb_w, nzt_w |
---|
| 207 | |
---|
| 208 | ! |
---|
| 209 | !-- (1) Compute rifs |
---|
| 210 | IF ( k == nzb_w ) THEN |
---|
| 211 | kk = nzb_w |
---|
| 212 | ELSE |
---|
| 213 | kk = k-1 |
---|
| 214 | ENDIF |
---|
| 215 | rifs = 0.5 * ( rif_wall(k,j,i,wall_index) + & |
---|
| 216 | a * rif_wall(k,j,i+1,1) + b * rif_wall(k,j+1,i,2) + & |
---|
| 217 | c1 * rif_wall(kk,j,i,3) + c2 * rif_wall(kk,j,i,4) & |
---|
| 218 | ) |
---|
| 219 | |
---|
| 220 | ! |
---|
| 221 | !-- Skip (2) to (4) of wall_fluxes, because here rifs is already available |
---|
| 222 | !-- from (1) |
---|
| 223 | |
---|
| 224 | ! |
---|
| 225 | !-- (5) Compute wall_flux (u'v', v'u', w'v', or w'u') |
---|
[55] | 226 | vel_zp = 0.5 * ( a * ( u(k,j,i) + u(k,j,i+1) ) + & |
---|
| 227 | b * ( v(k,j,i) + v(k,j+1,i) ) + & |
---|
| 228 | (c1+c2) * ( w(k,j,i) + w(k-1,j,i) ) & |
---|
| 229 | ) |
---|
| 230 | |
---|
[53] | 231 | IF ( rifs >= 0.0 ) THEN |
---|
| 232 | |
---|
| 233 | ! |
---|
| 234 | !-- Stable stratification (and neutral) |
---|
[55] | 235 | wall_flux(k) = kappa * vel_zp / & |
---|
| 236 | ( LOG( zp/z0(j,i) ) + 5.0*rifs * ( zp-z0(j,i) ) / zp ) |
---|
[53] | 237 | ELSE |
---|
| 238 | |
---|
| 239 | ! |
---|
| 240 | !-- Unstable stratification |
---|
| 241 | h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) ) |
---|
| 242 | h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) ) |
---|
| 243 | |
---|
| 244 | ! |
---|
| 245 | !-- If a borderline case occurs, the formula for stable stratification |
---|
| 246 | !-- must be used anyway, or else a zero division would occur in the |
---|
| 247 | !-- argument of the logarithm. |
---|
| 248 | IF ( h1 == 1.0 .OR. h2 == 1.0 ) THEN |
---|
[55] | 249 | wall_flux(k) = kappa * vel_zp / & |
---|
[53] | 250 | ( LOG( zp / z0(j,i) ) + & |
---|
| 251 | 5.0 * rifs * ( zp - z0(j,i) ) / zp & |
---|
| 252 | ) |
---|
| 253 | ELSE |
---|
[55] | 254 | wall_flux(k) = kappa * vel_zp / & |
---|
[53] | 255 | ( LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) & |
---|
| 256 | + 2.0 * ( ATAN( h2 ) - ATAN( h1 ) ) & |
---|
| 257 | ) |
---|
| 258 | ENDIF |
---|
| 259 | |
---|
| 260 | ENDIF |
---|
| 261 | wall_flux(k) = wall_flux(k) * ABS( wall_flux(k) ) |
---|
| 262 | |
---|
| 263 | ! |
---|
| 264 | !-- store rifs for next time step |
---|
| 265 | rif_wall(k,j,i,wall_index) = rifs |
---|
| 266 | |
---|
| 267 | ENDDO |
---|
| 268 | |
---|
| 269 | END SUBROUTINE wall_fluxes_e |
---|