Changeset 187 for palm/trunk/SOURCE


Ignore:
Timestamp:
Aug 6, 2008 4:25:09 PM (16 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)
Location:
palm/trunk/SOURCE
Files:
4 edited

Legend:

Unmodified
Added
Removed
  • palm/trunk/SOURCE/CURRENT_MODIFICATIONS

    r184 r187  
    4242Changed:
    4343-------
     44Modification of the integrated version of the profile function for momentum
     45for unstable stratification; more consistent flux definitions. (wall_fluxes,
     46production_e)
     47
    4448Strict grid matching along z is not needed for mg-solver. (check_parameters)
    4549
     
    7983(init_1d_model)
    8084
    81 advec_s_ups, advec_u_ups, advec_v_ups, advec_w_ups, calc_spectra, check_open, check_parameters, cpu_statistics, init_1d_model, init_3d_model, modules, palm, parin, poisfft, read_var_list, read_3d_binary, transpose, write_var_list, write_3d_binary
     85advec_s_ups, advec_u_ups, advec_v_ups, advec_w_ups, calc_spectra, check_open, check_parameters, cpu_statistics, init_1d_model, init_3d_model, modules, palm, parin, poisfft, production_e, read_var_list, read_3d_binary, transpose, wall_fluxes, write_var_list, write_3d_binary
    8286
    8387
    8488Errors:
    8589------
     90Bugfix: change definition of us_wall from 1D to 2D. Tests showed that this decreases
     91u* by some 10% and increases TKE and momentum fluxes by some 10% because friction
     92was underestimated and momentum fluxes were wrongly calculated due to the bug.
     93(prandtl_fluxes, wall_fluxes)
     94
    8695Bugfix: calculation of horizontal fluxes at vertical walls (diffusion_s)
    8796
     
    106115
    107116Introduce prefix_chr to ensure unique dvrp_file path.
     117
    108118small bugfixes for user_interface sample code (comments):
    109119- initialize ustvst with 0.0 as it is now computed only until nxr and nyn
    110120- two ALLOCATE statements moved from user_read_restart_data back to user_init
    111121- remove 'READ (13) u2_av' statement in user_read_restart_data
     122
    112123Bugfix: remove IF statement in plant_canopy_model_ij (plant_canopy_model)
     124
    113125Bugfix: divide sums(k,8) (e) and sums(k,34) (e*) by ngp_2dh_s_inner(k,sr)
    114126(like other scalars) (flow_statistics)
     127
    115128Bugfix: dopr_time_count was written on the binary file, which caused that
    116129NetCDF files newly created by restart files (no append of existing files!)
    117130contained uneccessary time levels. (read_3d_binary, write_3d_binary)
     131
    118132Bugfix: extra '*' removed in user_statistics sample code (user_interface)
     133
    119134Bugfix: a stop command was missing in some cases of the parallel branch (local_stop)
     135
    120136Bugfix in volume flow control for non-cyclic boundary conditions (pres)
     137
    121138Bugfix: misplaced #endif directives (combine_plot_fields)
    122139
     140check_parameters, diffusion_s, flow_statistics, init_dvrp, init_3d_model, local_stop, plant_canopy_model, poismg, prandtl_fluxes, pres, read_3d_binary, user_interface, wall_fluxes, write_3d_binary
    123141
    124 check_parameters, diffusion_s, flow_statistics, init_dvrp, init_3d_model, local_stop, plant_canopy_model, poismg, pres, read_3d_binary, user_interface, write_3d_binary
    125 
  • palm/trunk/SOURCE/prandtl_fluxes.f90

    r110 r187  
    44! Actual revisions:
    55! -----------------
    6 !
     6! Bugfix: change definition of us_wall from 1D to 2D:
     7!         Modification of the evaluation of the vertical turbulent momentum
     8!         fluxes u'w' and v'w'; the first usws that is computed corresponds
     9!         to -u'w'/u* and not as priorily assumed to (-u'w')**0.5, the first
     10!         vsws that is computed corresponds to -v'w'/u* and not as priorily
     11!         assumed to (-v'w')**0.5. Therefore, the intermediate result for
     12!         usws has to be multiplied by -u* instead by itself in order to
     13!         get u'w'. Accordingly, the intermediate result for vsws has to be
     14!         multiplied by -u* instead by itself in order to get v'w'. As u*
     15!         is calculated for the position of a scalar an additional
     16!         interpolation of u* to the position of u and v, respectively,
     17!         is necessary. As u* is not determined for the ghost points on
     18!         each PE, an additional exchange of information from neighbouring
     19!         PEs is necessary.! 
     20! Change: Modification of the integrated version of the profile function for
     21!         momentum for unstable stratification
    722!
    823! Former revisions:
     
    8196!--             Unstable stratification
    8297                a = SQRT( 1.0 - 16.0 * rif(j,i) )
    83                 b = SQRT( 1.0 - 16.0 * rif(j,i) / z_p * z0(j,i) )
    84 !
    85 !--             If a borderline case occurs, the formula for stable
    86 !--             stratification must be used anyway, or else a zero division
    87 !--             would occur in the argument of the logarithm
    88                 IF ( a == 1.0  .OR.  b == 1.0 )  THEN
    89                    ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) / (          &
    90                                      LOG( z_p / z0(j,i) ) +                   &
    91                                      5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
    92                                                                    )
    93                 ELSE
    94                    ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) / (          &
    95                                  LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) &
    96                                                                    )
    97                 ENDIF
     98                b = SQRT( 1.0 - 16.0 * rif(j,i) * z0(j,i) / z_p )
     99
     100                ts(j,i) = kappa * ( pt(k+1,j,i) - pt(k,j,i) ) /  (         &
     101                          LOG( z_p / z0(j,i) ) -                           &
     102                          2.0 * LOG( ( 1.0 + a ) / ( 1.0 + b ) ) )
    98103             ENDIF
    99104
     
    168173!
    169174!--          Unstable stratification
    170              a = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
    171              b = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rif(j,i) / z_p * z0(j,i) ) )
    172 !
    173 !--          If a borderline case occurs, the formula for stable stratification
    174 !--          must be used anyway, or else a zero division would occur in the
    175 !--          argument of the logarithm.
    176              IF ( a == 1.0  .OR.  b == 1.0 )  THEN
    177                 us(j,i) = kappa * uv_total / (                                &
    178                                      LOG( z_p / z0(j,i) ) +                   &
    179                                      5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
    180                                              )
    181              ELSE
    182                 us(j,i) = kappa * uv_total / (                               &
    183                               LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) + &
    184                               2.0 * ( ATAN( b ) - ATAN( a ) )                &
    185                                              )
    186              ENDIF
     175             a = SQRT( SQRT( 1.0 - 16.0 * rif(j,i) ) )
     176             b = SQRT( SQRT( 1.0 - 16.0 * rif(j,i) / z_p * z0(j,i) ) )
     177
     178             us(j,i) = kappa * uv_total / (                                  &
     179                       LOG( z_p / z0(j,i) ) -                                &
     180                       LOG( ( 1.0 + a )**2 * ( 1.0 + a**2 ) / (              &
     181                            ( 1.0 + b )**2 * ( 1.0 + b**2 )   ) ) +          &
     182                            2.0 * ( ATAN( a ) - ATAN( b ) )                  &
     183                                           )
    187184          ENDIF
    188185       ENDDO
    189186    ENDDO
    190187
     188!
     189!-- Values of us at ghost point locations are needed for the evaluation of usws
     190!-- and vsws.
     191    CALL exchange_horiz_2d( us )
    191192!
    192193!-- Compute u'w' for the total model domain.
     
    212213!
    213214!--          Unstable stratification
    214              a = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifm ) )
    215              b = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifm / z_p * z0(j,i) ) )
    216 !
    217 !--          If a borderline case occurs, the formula for stable stratification
    218 !--          must be used anyway, or else a zero division would occur in the
    219 !--          argument of the logarithm.
    220              IF ( a == 1.0  .OR.  B == 1.0 )  THEN
    221                 usws(j,i) = kappa * u(k+1,j,i) / (                           &
    222                                         LOG( z_p / z0(j,i) ) +               &
    223                                         5.0 * rifm * ( z_p - z0(j,i) ) / z_p &
     215             a = SQRT( SQRT( 1.0 - 16.0 * rifm ) )
     216             b = SQRT( SQRT( 1.0 - 16.0 * rifm / z_p * z0(j,i) ) )
     217
     218             usws(j,i) = kappa * u(k+1,j,i) / (                           &
     219                         LOG( z_p / z0(j,i) ) -                           &
     220                         LOG( (1.0 + a )**2 * ( 1.0 + a**2 ) / (          &
     221                              (1.0 + b )**2 * ( 1.0 + b**2 )   ) ) +      &
     222                              2.0 * ( ATAN( a ) - ATAN( b ) )             &
    224223                                                 )
    225              ELSE
    226                 usws(j,i) = kappa * u(k+1,j,i) / (                           &
    227                               LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) + &
    228                               2.0 * ( ATAN( b ) - ATAN( a ) )                &
    229                                                  )
    230              ENDIF
    231224          ENDIF
    232           usws(j,i) = -usws(j,i) * ABS( usws(j,i) )
     225          usws(j,i) = -usws(j,i) * 0.5 * ( us(j,i-1) + us(j,i) )
    233226       ENDDO
    234227    ENDDO
     
    257250!
    258251!--          Unstable stratification
    259              a = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifm ) )
    260              b = 1.0 / SQRT( SQRT( 1.0 - 16.0 * rifm / z_p * z0(j,i) ) )
    261 !
    262 !--          If a borderline case occurs, the formula for stable stratification
    263 !--          must be used anyway, or else a zero division would occur in the
    264 !--          argument of the logarithm.
    265              IF ( a == 1.0  .OR.  b == 1.0 )  THEN
    266                 vsws(j,i) = kappa * v(k+1,j,i) / (                           &
    267                                         LOG( z_p / z0(j,i) ) +               &
    268                                         5.0 * rifm * ( z_p - z0(j,i) ) / z_p &
     252             a = SQRT( SQRT( 1.0 - 16.0 * rifm ) )
     253             b = SQRT( SQRT( 1.0 - 16.0 * rifm / z_p * z0(j,i) ) )
     254
     255             vsws(j,i) = kappa * v(k+1,j,i) / (                           &
     256                         LOG( z_p / z0(j,i) ) -                           &
     257                         LOG( (1.0 + a )**2 * ( 1.0 + a**2 ) / (          &
     258                              (1.0 + b )**2 * ( 1.0 + b**2 )   ) ) +      &
     259                              2.0 * ( ATAN( a ) - ATAN( b ) )             &
    269260                                                 )
    270              ELSE
    271                 vsws(j,i) = kappa * v(k+1,j,i) / (                           &
    272                               LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) + &
    273                               2.0 * ( ATAN( b ) - ATAN( a ) )                &
    274                                                  )
    275              ENDIF
    276261          ENDIF
    277           vsws(j,i) = -vsws(j,i) * ABS( vsws(j,i) )
     262          vsws(j,i) = -vsws(j,i) * 0.5 * ( us(j-1,i) + us(j,i) )
    278263       ENDDO
    279264    ENDDO
     
    317302!
    318303!--                Unstable stratification
    319                    a = SQRT( 1.0 - 16.0 * rif(j,i) )
    320                    b = SQRT( 1.0 - 16.0 * rif(j,i) / z_p * z0(j,i) )
    321 !
    322 !--                If a borderline case occurs, the formula for stable
    323 !--                stratification must be used anyway, or else a zero division
    324 !--                would occur in the argument of the logarithm.
    325                    IF ( a == 1.0  .OR.  b == 1.0 )  THEN
    326                       qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / (         &
    327                                      LOG( z_p / z0(j,i) ) +                   &
    328                                      5.0 * rif(j,i) * ( z_p - z0(j,i) ) / z_p &
    329                                                                     )
    330                    ELSE
    331                       qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) / (          &
    332                                   LOG( (1.0+b) / (1.0-b) * (1.0-a) / (1.0+a) ) &
    333                                                                     )
    334                    ENDIF
     304                   a = SQRT( 1.0 - 16.0 * rif(j,i) )
     305                   b = SQRT( 1.0 - 16.0 * rif(j,i) * z0(j,i) / z_p )
     306 
     307                   qs(j,i) = kappa * ( q(k+1,j,i) - q(k,j,i) ) /   (   &
     308                             LOG( z_p / z0(j,i) ) -                    &
     309                              2.0 * LOG( (1.0 + a ) / ( 1.0 + b ) ) )
    335310                ENDIF
    336311
     
    341316
    342317!
    343 !-- Exchange the boundaries for u* and the momentum fluxes (fluxes only for
    344 !-- completeness's sake).
    345     CALL exchange_horiz_2d( us )
     318!-- Exchange the boundaries for the momentum fluxes (only for sake of
     319!-- completeness)
    346320    CALL exchange_horiz_2d( usws )
    347321    CALL exchange_horiz_2d( vsws )
  • palm/trunk/SOURCE/production_e.f90

    r139 r187  
    44! Actual revisions:
    55! -----------------
    6 !
     6! Change: add 'minus' sign to fluxes obtained from subroutine wall_fluxes_e for
     7!         consistency with subroutine wall_fluxes
    78!
    89! Former revisions:
     
    164165                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
    165166                                          usvs, 1.0, 0.0, 0.0, 0.0 )
    166                       dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i)
    167 !                      dudy = wall_e_y(j,i) * usvs(k,j,i) / km(k,j,i)
     167                      dudy = - wall_e_y(j,i) * usvs(k) / km(k,j,i)
     168!                      dudy = - wall_e_y(j,i) * usvs(k,j,i) / km(k,j,i)
    168169                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
    169170                                          wsvs, 0.0, 0.0, 1.0, 0.0 )
    170                       dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i)
    171 !                      dwdy = wall_e_y(j,i) * wsvs(k,j,i) / km(k,j,i)
     171                      dwdy = - wall_e_y(j,i) * wsvs(k) / km(k,j,i)
     172!                      dwdy = - wall_e_y(j,i) * wsvs(k,j,i) / km(k,j,i)
    172173                   ELSE
    173174                      dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
     
    180181                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
    181182                                          vsus, 0.0, 1.0, 0.0, 0.0 )
    182                       dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i)
    183 !                      dvdx = wall_e_x(j,i) * vsus(k,j,i) / km(k,j,i)
     183                      dvdx = - wall_e_x(j,i) * vsus(k) / km(k,j,i)
     184!                      dvdx = - wall_e_x(j,i) * vsus(k,j,i) / km(k,j,i)
    184185                      CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
    185186                                          wsus, 0.0, 0.0, 0.0, 1.0 )
    186                       dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i)
    187 !                      dwdx = wall_e_x(j,i) * wsus(k,j,i) / km(k,j,i)
     187                      dwdx = - wall_e_x(j,i) * wsus(k) / km(k,j,i)
     188!                      dwdx = - wall_e_x(j,i) * wsus(k,j,i) / km(k,j,i)
    188189                   ELSE
    189190                      dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
     
    220221
    221222                      IF ( wall_e_y(j,i) /= 0.0 )  THEN
    222                          dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i)
    223 !                         dudy = wall_e_y(j,i) * usvs(k,j,i) / km(k,j,i)
    224                          dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i)
    225 !                         dwdy = wall_e_y(j,i) * wsvs(k,j,i) / km(k,j,i)
     223                         dudy = - wall_e_y(j,i) * usvs(k) / km(k,j,i)
     224!                         dudy = - wall_e_y(j,i) * usvs(k,j,i) / km(k,j,i)
     225                         dwdy = - wall_e_y(j,i) * wsvs(k) / km(k,j,i)
     226!                         dwdy = - wall_e_y(j,i) * wsvs(k,j,i) / km(k,j,i)
    226227                      ELSE
    227228                         dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
     
    232233
    233234                      IF ( wall_e_x(j,i) /= 0.0 )  THEN
    234                          dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i)
    235 !                         dvdx = wall_e_x(j,i) * vsus(k,j,i) / km(k,j,i)
    236                          dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i)
    237 !                         dwdx = wall_e_x(j,i) * wsus(k,j,i) / km(k,j,i)
     235                         dvdx = - wall_e_x(j,i) * vsus(k) / km(k,j,i)
     236!                         dvdx = - wall_e_x(j,i) * vsus(k,j,i) / km(k,j,i)
     237                         dwdx = - wall_e_x(j,i) * wsus(k) / km(k,j,i)
     238!                         dwdx = - wall_e_x(j,i) * wsus(k,j,i) / km(k,j,i)
    238239                      ELSE
    239240                         dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
     
    629630                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
    630631                                    usvs, 1.0, 0.0, 0.0, 0.0 )
    631                 dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i)
     632                dudy = - wall_e_y(j,i) * usvs(k) / km(k,j,i)
    632633                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
    633634                                    wsvs, 0.0, 0.0, 1.0, 0.0 )
    634                 dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i)
     635                dwdy = - wall_e_y(j,i) * wsvs(k) / km(k,j,i)
    635636             ELSE
    636637                dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
     
    643644                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
    644645                                    vsus, 0.0, 1.0, 0.0, 0.0 )
    645                 dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i)
     646                dvdx = - wall_e_x(j,i) * vsus(k) / km(k,j,i)
    646647                CALL wall_fluxes_e( i, j, k, nzb_diff_s_outer(j,i)-2, &
    647648                                    wsus, 0.0, 0.0, 0.0, 1.0 )
    648                 dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i)
     649                dwdx = - wall_e_x(j,i) * wsus(k) / km(k,j,i)
    649650             ELSE
    650651                dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
     
    679680
    680681                IF ( wall_e_y(j,i) /= 0.0 )  THEN
    681                    dudy = wall_e_y(j,i) * usvs(k) / km(k,j,i)
    682                    dwdy = wall_e_y(j,i) * wsvs(k) / km(k,j,i)
     682                   dudy = - wall_e_y(j,i) * usvs(k) / km(k,j,i)
     683                   dwdy = - wall_e_y(j,i) * wsvs(k) / km(k,j,i)
    683684                ELSE
    684685                   dudy = 0.25 * ( u(k,j+1,i) + u(k,j+1,i+1) - &
     
    689690
    690691                IF ( wall_e_x(j,i) /= 0.0 )  THEN
    691                    dvdx = wall_e_x(j,i) * vsus(k) / km(k,j,i)
    692                    dwdx = wall_e_x(j,i) * wsus(k) / km(k,j,i)
     692                   dvdx = - wall_e_x(j,i) * vsus(k) / km(k,j,i)
     693                   dwdx = - wall_e_x(j,i) * wsus(k) / km(k,j,i)
    693694                ELSE
    694695                   dvdx = 0.25 * ( v(k,j,i+1) + v(k,j+1,i+1) - &
  • 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.