Ignore:
Timestamp:
Mar 7, 2007 12:33:47 PM (15 years ago)
Author:
raasch
Message:

preliminary version, several changes to be explained later

File:
1 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/wall_fluxes.f90

    r52 r53  
    2626    IMPLICIT NONE
    2727
    28     INTEGER ::  i, j, k, nzb_w, nzt_w
    29     REAL    ::  a, b, c1, c2, h1, h2, delta_p
     28    INTEGER ::  i, j, k, nzb_w, nzt_w, wall_index
     29    REAL    ::  a, b, c1, c2, h1, h2, zp
    3030
    3131    REAL ::  pts, pt_i, rifs, u_i, v_i, us_wall, vel_total, ws, wspts
     
    3434
    3535
    36     delta_p   = 0.5 * ( (a+c1) * dy + (b+c2) * dx )
    37     wall_flux = 0.0
     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 )
    3839
    3940!
     
    4142!-- the relevant variable is defined
    4243    DO  k = nzb_w, nzt_w
     44
    4345!
    4446!--    (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
     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
    5765!
    5866!--    (2) Compute wall-parallel absolute velocity vel_total
    5967       vel_total = SQRT( ws**2 + ( a+c1 ) * u_i**2 + ( b+c2 ) * v_i**2 )
     68
    6069!
    6170!--    (3) Compute wall friction velocity us_wall
    6271       IF ( rifs >= 0.0 )  THEN
     72
    6373!
    6474!--       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  &
     75          us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +               &
     76                                          5.0 * rifs * ( zp - z0(j,i) ) / zp  &
    6877                                        )
    6978       ELSE
     79
    7080!
    7181!--       Unstable stratification
    7282          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) ) )
     83          h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) )
     84
    7485!
    7586!--       If a borderline case occurs, the formula for stable stratification
     
    7788!--       argument of the logarithm.
    7889          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                                                    )
     90             us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +             &
     91                                            5.0 * rifs * ( zp - z0(j,i) ) / zp &
     92                                           )
    8393          ELSE
    84              us_wall = kappa * vel_total / (                           &
     94             us_wall = kappa * vel_total / (                                   &
    8595                            LOG( (1.0+h2) / (1.0-h2) * (1.0-h1) / (1.0+h1) ) + &
    8696                            2.0 * ( ATAN( h2 ) - ATAN( h1 ) )                  &
    87                                                    )
     97                                           )
    8898          ENDIF
    89        ENDIF
    90 !
    91 !--    (4) Compute delta_p/L (corresponds to neutral Richardson flux number
     99
     100       ENDIF
     101
     102!
     103!--    (4) Compute zp/L (corresponds to neutral Richardson flux number
    92104!--        rifs)
    93        rifs = -1.0 * delta_p * kappa * g * wspts / &
    94                        ( pt_i * ( us_wall**3 + 1E-30 ) )
     105       rifs = -1.0 * zp * kappa * g * wspts / ( pt_i * ( us_wall**3 + 1E-30 ) )
     106
    95107!
    96108!--    Limit the value range of the Richardson numbers.
     
    101113       IF ( rifs < rif_min )  rifs = rif_min
    102114       IF ( rifs > rif_max )  rifs = rif_max
     115
    103116!
    104117!--    (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
    105118       IF ( rifs >= 0.0 )  THEN
     119
    106120!
    107121!--       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
     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                         )
     127       ELSE
     128
    114129!
    115130!--       Unstable stratification
    116131          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) ) )
     132          h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) )
     133
    118134!
    119135!--       If a borderline case occurs, the formula for stable stratification
     
    121137!--       argument of the logarithm.
    122138          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                                                                             )
     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                            )
    128144          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                                                                             )
     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                            )
    134150          ENDIF
     151
    135152       ENDIF
    136153       wall_flux(k) = -wall_flux(k) * ABS( wall_flux(k) )
     
    138155!
    139156!--    store rifs for next time step
    140        rif_wall(k,j,i,NINT(a+2*b+3*c1+4*c2)) = rifs
     157       rif_wall(k,j,i,wall_index) = rifs
    141158
    142159    ENDDO
     160
    143161 END SUBROUTINE wall_fluxes
     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
     192    REAL    ::  a, b, c1, c2, h1, h2, zp
     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')
     226       IF ( rifs >= 0.0 )  THEN
     227
     228!
     229!--       Stable stratification (and neutral)
     230          wall_flux(k) = kappa *                                          &
     231                         ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / &
     232                         (  LOG( zp / z0(j,i) ) +                         &
     233                            5.0 * rifs * ( zp - z0(j,i) ) / zp            &
     234                         )
     235       ELSE
     236
     237!
     238!--       Unstable stratification
     239          h1 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs ) )
     240          h2 = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifs / zp * z0(j,i) ) )
     241
     242!
     243!--       If a borderline case occurs, the formula for stable stratification
     244!--       must be used anyway, or else a zero division would occur in the
     245!--       argument of the logarithm.
     246          IF ( h1 == 1.0  .OR.  h2 == 1.0 )  THEN
     247             wall_flux(k) = kappa *                                          &
     248                            ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / &
     249                            ( LOG( zp / z0(j,i) ) +                          &
     250                              5.0 * rifs * ( zp - z0(j,i) ) / zp             &
     251                            )
     252          ELSE
     253             wall_flux(k) = kappa *                                            &
     254                            ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) /   &
     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
Note: See TracChangeset for help on using the changeset viewer.