Ignore:
Timestamp:
Aug 6, 2008 4:25:09 PM (13 years ago)
Author:
letzel
Message:
  • new: descriptions of plant canopy model and prandtl layer (trunk/DOC/misc)
  • changed: more consistent flux definitions; modification of the integrated version of the profile function for momentum for unstable stratification (wall_fluxes, production_e)
  • bugfix: change definition of us_wall from 1D to 2D (prandtl_fluxes, wall_fluxes)
File:
1 edited

Legend:

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

    r75 r187  
    33! Actual revisions:
    44! -----------------
    5 !
     5! Bugfix: change definition of us_wall from 1D to 2D:
     6!         Modification of the evaluation of the vertical turbulent momentum
     7!         fluxes u'w' and v'w'; the first usws that is computed corresponds
     8!         to -u'w'/u* and not as priorily assumed to (-u'w')**0.5, the first
     9!         vsws that is computed corresponds to -v'w'/u* and not as priorily
     10!         assumed to (-v'w')**0.5. Therefore, the intermediate result for
     11!         usws has to be multiplied by -u* instead by itself in order to
     12!         get u'w'. Accordingly, the intermediate result for vsws has to be
     13!         multiplied by -u* instead by itself in order to get v'w'.
     14!         This requires the calculation of us_wall (and vel_total, u_i, v_i, ws)
     15!         also in wall_fluxes_e.
     16! Bugfix: storage of rifs to rifs_wall in wall_fluxes_e removed
     17! Change: add 'minus' sign to fluxes produced by subroutine wall_fluxes_e for
     18!         consistency with subroutine wall_fluxes
     19! Change: Modification of the integrated version of the profile function for
     20!         momentum for unstable stratification
    621!
    722! Former revisions:
     
    6984!
    7085!--             All subsequent variables are computed for the respective
    71 !--             location where the relevant variable is defined
     86!--             location where the respective flux is defined.
    7287                DO  k = nzb_uvw_inner(j,i)+1, nzb_uvw_outer(j,i)
    7388
     
    109124!
    110125!--                   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 
     126                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
     127                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
     128
     129                      us_wall = kappa * vel_total / (                          &
     130                           LOG( zp / z0(j,i) ) -                               &
     131                           LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
     132                                ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
     133                                2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
     134                                                    )
    129135                   ENDIF
    130136
     
    160166!
    161167!--                   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 
     168                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
     169                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
     170
     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                           LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
     175                                ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
     176                                2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
     177                                                                            )
    183178                   ENDIF
    184                    wall_flux(k,j,i) = -wall_flux(k,j,i) * ABS(wall_flux(k,j,i))
     179                   wall_flux(k,j,i) = -wall_flux(k,j,i) * us_wall
    185180
    186181!
     
    226221!
    227222!--    All subsequent variables are computed for the respective location where
    228 !--    the relevant variable is defined
     223!--    the respective flux is defined.
    229224       DO  k = nzb_w, nzt_w
    230225
     
    266261!
    267262!--          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 
     263             h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
     264             h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
     265
     266             us_wall = kappa * vel_total / (                          &
     267                  LOG( zp / z0(j,i) ) -                               &
     268                  LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
     269                       ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
     270                       2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
     271                                           )
    286272          ENDIF
    287273
     
    315301!
    316302!--          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 
     303             h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
     304             h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
     305
     306             wall_flux(k) = kappa *                               &
     307                  ( a*u(k,j,i) + b*v(k,j,i) + (c1+c2)*w(k,j,i) ) / (  &
     308                  LOG( zp / z0(j,i) ) -                               &
     309                  LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
     310                       ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
     311                       2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
     312                                                                   )
    338313          ENDIF
    339           wall_flux(k) = -wall_flux(k) * ABS( wall_flux(k) )
     314          wall_flux(k) = -wall_flux(k) * us_wall
    340315
    341316!
     
    371346
    372347       INTEGER ::  i, j, k, kk, wall_index
    373        REAL    ::  a, b, c1, c2, h1, h2, vel_zp, zp
     348       REAL    ::  a, b, c1, c2, h1, h2, u_i, v_i, us_wall, vel_total, vel_zp, &
     349                   ws, zp
    374350
    375351       REAL ::  rifs
     
    388364             IF ( wall(j,i) /= 0.0 )  THEN
    389365!
    390 !--             All subsequent variables are computed for the respective
    391 !--             location where the relevant variable is defined
     366!--             All subsequent variables are computed for scalar locations.
    392367                DO  k = nzb_diff_s_inner(j,i)-1, nzb_diff_s_outer(j,i)-2
    393 
    394 !
    395 !--                (1) Compute rifs
     368!
     369!--                (1) Compute rifs, u_i, v_i, and ws
    396370                   IF ( k == nzb_diff_s_inner(j,i)-1 )  THEN
    397371                      kk = nzb_diff_s_inner(j,i)-1
     
    404378                                 )
    405379
    406 !
    407 !--                Skip (2) to (4) of wall_fluxes, because here rifs is
    408 !--                already available from (1)
    409 
     380                   u_i   = 0.5 * ( u(k,j,i) + u(k,j,i+1) )
     381                   v_i   = 0.5 * ( v(k,j,i) + v(k,j+1,i) )
     382                   ws    = 0.5 * ( w(k,j,i) + w(k-1,j,i) )
     383!
     384!--                (2) Compute wall-parallel absolute velocity vel_total and
     385!--                interpolate appropriate velocity component vel_zp.
     386                   vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 )
     387                   vel_zp = 0.5 * ( a * u_i + b * v_i + (c1+c2) * ws )
     388!
     389!--                (3) Compute wall friction velocity us_wall
     390                   IF ( rifs >= 0.0 )  THEN
     391
     392!
     393!--                   Stable stratification (and neutral)
     394                      us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +    &
     395                                            5.0 * rifs * ( zp - z0(j,i) ) / zp &
     396                                                    )
     397                   ELSE
     398
     399!
     400!--                   Unstable stratification
     401                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
     402                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
     403
     404                      us_wall = kappa * vel_total / (                          &
     405                           LOG( zp / z0(j,i) ) -                               &
     406                           LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
     407                                ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
     408                                2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
     409                                                    )
     410                   ENDIF
     411
     412!
     413!--                Skip step (4) of wall_fluxes, because here rifs is already
     414!--                available from (1)
    410415!
    411416!--                (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                                   )
    416417
    417418                   IF ( rifs >= 0.0 )  THEN
     
    425426!
    426427!--                   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 
     428                      h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
     429                      h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
     430
     431                      wall_flux(k,j,i) = kappa * vel_zp / (                    &
     432                           LOG( zp / z0(j,i) ) -                               &
     433                           LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
     434                                ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
     435                                2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
     436                                                          )
    446437                   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
     438                   wall_flux(k,j,i) = - wall_flux(k,j,i) * us_wall
    452439
    453440                ENDDO
     
    476463
    477464       INTEGER ::  i, j, k, kk, nzb_w, nzt_w, wall_index
    478        REAL    ::  a, b, c1, c2, h1, h2, vel_zp, zp
     465       REAL    ::  a, b, c1, c2, h1, h2, u_i, v_i, us_wall, vel_total, vel_zp, &
     466                   ws, zp
    479467
    480468       REAL ::  rifs
     
    488476
    489477!
    490 !--    All subsequent variables are computed for the respective location where
    491 !--    the relevant variable is defined
     478!--    All subsequent variables are computed for scalar locations.
    492479       DO  k = nzb_w, nzt_w
    493480
    494481!
    495 !--       (1) Compute rifs
     482!--       (1) Compute rifs, u_i, v_i, and ws
    496483          IF ( k == nzb_w )  THEN
    497484             kk = nzb_w
     
    504491                        )
    505492
    506 !
    507 !--       Skip (2) to (4) of wall_fluxes, because here rifs is already available
    508 !--       from (1)
    509 
     493          u_i   = 0.5 * ( u(k,j,i) + u(k,j,i+1) )
     494          v_i   = 0.5 * ( v(k,j,i) + v(k,j+1,i) )
     495          ws    = 0.5 * ( w(k,j,i) + w(k-1,j,i) )
     496!
     497!--       (2) Compute wall-parallel absolute velocity vel_total and
     498!--       interpolate appropriate velocity component vel_zp.
     499          vel_total = SQRT( ws**2 + (a+c1) * u_i**2 + (b+c2) * v_i**2 )
     500          vel_zp = 0.5 * ( a * u_i + b * v_i + (c1+c2) * ws )
     501!
     502!--       (3) Compute wall friction velocity us_wall
     503          IF ( rifs >= 0.0 )  THEN
     504
     505!
     506!--          Stable stratification (and neutral)
     507             us_wall = kappa * vel_total / ( LOG( zp / z0(j,i) ) +             &
     508                                            5.0 * rifs * ( zp - z0(j,i) ) / zp &
     509                                           )
     510          ELSE
     511
     512!
     513!--          Unstable stratification
     514             h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
     515             h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
     516
     517             us_wall = kappa * vel_total / (                          &
     518                  LOG( zp / z0(j,i) ) -                               &
     519                  LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
     520                       ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
     521                       2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
     522                                           )
     523          ENDIF
     524
     525!
     526!--       Skip step (4) of wall_fluxes, because here rifs is already
     527!--       available from (1)
    510528!
    511529!--       (5) Compute wall_flux (u'v', v'u', w'v', or w'u')
     530!--       First interpolate the velocity (this is different from
     531!--       subroutine wall_fluxes because fluxes in subroutine
     532!--       wall_fluxes_e are defined at scalar locations).
    512533          vel_zp = 0.5 * (       a * ( u(k,j,i) + u(k,j,i+1) ) +  &
    513534                                 b * ( v(k,j,i) + v(k,j+1,i) ) +  &
     
    525546!
    526547!--          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 
     548             h1 = SQRT( SQRT( 1.0 - 16.0 * rifs ) )
     549             h2 = SQRT( SQRT( 1.0 - 16.0 * rifs * z0(j,i) / zp ) )
     550
     551             wall_flux(k) = kappa * vel_zp / (                        &
     552                  LOG( zp / z0(j,i) ) -                               &
     553                  LOG( ( 1.0 + h1 )**2 * ( 1.0 + h1**2 ) / (          &
     554                       ( 1.0 + h2 )**2 * ( 1.0 + h2**2 )   ) ) +      &
     555                       2.0 * ( ATAN( h1 ) - ATAN( h2 ) )              &
     556                                                 )
    546557          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
     558          wall_flux(k) = - wall_flux(k) * us_wall
    552559
    553560       ENDDO
Note: See TracChangeset for help on using the changeset viewer.