Changeset 4473 for palm/trunk/SOURCE


Ignore:
Timestamp:
Mar 25, 2020 9:04:07 PM (5 years ago)
Author:
gronemeier
Message:

revised naming of mixing length and some further cleaning of mixing-length calculation

Location:
palm/trunk/SOURCE
Files:
3 edited

Legend:

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

    r4444 r4473  
    2525! -----------------
    2626! $Id$
     27! revised message if wall_adjustment is used
     28!
     29! 4444 2020-03-05 15:59:50Z raasch
    2730! bugfix: cpp-directives for serial mode added
    2831!
     
    15391542       IF ( e_init > 0.0_wp )  WRITE ( io, 455 )  e_init
    15401543       IF ( e_min > 0.0_wp )  WRITE ( io, 454 )  e_min
    1541        IF ( wall_adjustment )  WRITE ( io, 453 )  wall_adjustment_factor
     1544       IF ( wall_adjustment )  WRITE ( io, 453 )
    15421545    ENDIF
    15431546    IF ( rans_mode )  THEN
     
    18841887451 FORMAT ('    Diffusion coefficients are constant:'/ &
    18851888            '    Km = ',F6.2,' m**2/s   Kh = ',F6.2,' m**2/s   Pr = ',F5.2)
    1886 453 FORMAT ('    Mixing length is limited to',F5.2,' * z')
     1889453 FORMAT ('    Mixing length is limited close to surfaces')
    18871890454 FORMAT ('    TKE is not allowed to fall below ',E9.2,' (m/s)**2')
    18881891455 FORMAT ('    initial TKE is prescribed as ',E9.2,' (m/s)**2')
  • palm/trunk/SOURCE/modules.f90

    r4472 r4473  
    2525! -----------------
    2626! $Id$
     27! moved wall_adjustment_factor to turbulence_closure_mod
     28!
     29! 4472 2020-03-24 12:21:00Z Giersch
    2730! Additional switch added to activate calculations in flow_statistics for the
    2831! kolmogorov length scale
     
    961964    REAL(wp) ::  v_bulk = 0.0_wp                               !< namelist parameter
    962965    REAL(wp) ::  v_gtrans = 0.0_wp                             !< transformed wind component (galilei transformation)
    963     REAL(wp) ::  wall_adjustment_factor = 1.8_wp               !< adjustment factor for mixing length l
    964966    REAL(wp) ::  zeta_max = 20.0_wp                            !< namelist parameter
    965967    REAL(wp) ::  zeta_min = -20.0_wp                           !< namelist parameter
  • palm/trunk/SOURCE/turbulence_closure_mod.f90

    r4457 r4473  
    2525! -----------------
    2626! $Id$
     27! - rename l-grid to gridsize-geometric-mean
     28!          l-wall to ml-wall-adjusted
     29!          l-stable to ml-stratification
     30!          l-black to ml-blackadar
     31!          l-v to ml-local-profile
     32!          l-max to max-length-scale
     33!          l to ml
     34! - adjust some comments
     35! - corrected definition of wall-adjusted mixing length to include
     36!   gridsize-geometric-mean
     37! - moved definition of wall_adjustment_factor to this module
     38!
     39! 4457 2020-03-11 14:20:43Z raasch
    2740! use statement for exchange horiz added
    28 ! 
     41!
    2942! 4433 2020-02-28 22:14:43Z gronemeier
    3043! remove warning for newly implemented RANS mode
     
    97110! ------------
    98111!> This module contains the available turbulence closures for PALM.
    99 !>
    100112!>
    101113!> @todo test initialization for all possibilities
     
    181193
    182194
    183     REAL(wp) ::  c_0                !< constant used for diffusion coefficient and dissipation (dependent on mode RANS/LES)
    184     REAL(wp) ::  c_1                !< model constant for RANS mode
    185     REAL(wp) ::  c_2                !< model constant for RANS mode
    186     REAL(wp) ::  c_3                !< model constant for RANS mode
    187     REAL(wp) ::  c_4                !< model constant for RANS mode
    188     REAL(wp) ::  l_max              !< maximum length scale for Blackadar mixing length
    189     REAL(wp) ::  dsig_e = 1.0_wp    !< factor to calculate Ke from Km (1/sigma_e)
    190     REAL(wp) ::  dsig_diss = 1.0_wp !< factor to calculate K_diss from Km (1/sigma_diss)
    191 
    192     REAL(wp), DIMENSION(0:4) :: rans_const_c = &       !< model constants for RANS mode (namelist param)
     195    REAL(wp) ::  c_0                 !< constant used for diffusion coefficient and dissipation (dependent on mode RANS/LES)
     196    REAL(wp) ::  c_1                 !< model constant for RANS mode
     197    REAL(wp) ::  c_2                 !< model constant for RANS mode
     198    REAL(wp) ::  c_3                 !< model constant for RANS mode
     199    REAL(wp) ::  c_4                 !< model constant for RANS mode
     200    REAL(wp) ::  dsig_e = 1.0_wp     !< factor to calculate Ke from Km (1/sigma_e)
     201    REAL(wp) ::  dsig_diss = 1.0_wp  !< factor to calculate K_diss from Km (1/sigma_diss)
     202    REAL(wp) ::  length_scale_max    !< maximum length scale for Blackadar mixing length
     203    REAL(wp) ::  wall_adjustment_factor = 1.8_wp  !< adjustment factor for mixing length
     204
     205    REAL(wp), DIMENSION(0:4) :: rans_const_c = &        !< model constants for RANS mode (namelist param)
    193206       (/ 0.55_wp, 1.44_wp, 1.92_wp, 1.44_wp, 0.0_wp /) !> default values fit for standard-tke-e closure
    194207
    195     REAL(wp), DIMENSION(2) :: rans_const_sigma = &     !< model constants for RANS mode, sigma values (sigma_e, sigma_diss) (namelist param)
    196        (/ 1.0_wp, 1.30_wp /)
    197 
    198     REAL(wp), DIMENSION(:), ALLOCATABLE ::  l_black    !< mixing length according to Blackadar
    199 
    200     REAL(wp), DIMENSION(:), ALLOCATABLE ::  l_grid     !< geometric mean of grid sizes dx, dy, dz
    201 
    202     REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  l_wall !< near-wall mixing length
     208    REAL(wp), DIMENSION(2) :: rans_const_sigma = &      !< model constants for RANS mode, sigma values (namelist param)
     209       (/ 1.0_wp, 1.30_wp /)                            !> (sigma_e, sigma_diss)
     210
     211    REAL(wp), ALLOCATABLE, DIMENSION(:)     ::  ml_blackadar             !< mixing length according to Blackadar
     212    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::  ml_wall_adjusted         !< mixing length adjusted according to near-by walls
    203213
    204214!
     
    11471157                   DO  j = nysg, nyng
    11481158                      DO  k = nzb+1, nzt
    1149                          km(k,j,i) = c_0 * l_wall(k,j,i) * SQRT( e_init )
     1159                         km(k,j,i) = c_0 * SQRT( e_init ) * MIN( ml_wall_adjusted(k,j,i), &
     1160                                                                 ml_blackadar(k) )
    11501161                      ENDDO
    11511162                   ENDDO
     
    11911202                DO  j = nysg, nyng
    11921203                   DO  k = nzb+1, nzt
    1193                       km(k,j,i) = c_0 * l_wall(k,j,i) * SQRT( e_init )
     1204                      km(k,j,i) = c_0 * SQRT( e_init ) * ml_wall_adjusted(k,j,i)
    11941205                   ENDDO
    11951206                ENDDO
     
    13681379
    13691380    USE control_parameters,                                                    &
    1370         ONLY:  f, message_string, wall_adjustment, wall_adjustment_factor
     1381        ONLY:  f, message_string, wall_adjustment
    13711382
    13721383    USE exchange_horiz_mod,                                                    &
     
    14061417    INTEGER(KIND=1), DIMENSION(:,:,:), ALLOCATABLE :: vicinity !< contains topography information of the vicinity of (i/j/k)
    14071418
    1408     REAL(wp) :: radius          !< search radius in meter
    1409 
    1410     ALLOCATE( l_grid(1:nzt) )
    1411     ALLOCATE( l_wall(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     1419    REAL(wp) :: distance_up            !< distance of grid box center to its boundary in upper direction
     1420    REAL(wp) :: distance_down          !< distance of grid box center to its boundary in lower direction
     1421    REAL(wp) :: distance_ns            !< distance of grid box center to its boundary in y direction
     1422    REAL(wp) :: distance_lr            !< distance of grid box center to its boundary in x direction
     1423    REAL(wp) :: distance_edge_yz_down  !< distance of grid box center to its boundary along yz diagonal (down)
     1424    REAL(wp) :: distance_edge_yz_up    !< distance of grid box center to its boundary along yz diagonal (up)
     1425    REAL(wp) :: distance_edge_xz_down  !< distance of grid box center to its boundary along xz diagonal (down)
     1426    REAL(wp) :: distance_edge_xz_up    !< distance of grid box center to its boundary along xz diagonal (up)
     1427    REAL(wp) :: distance_edge_xy       !< distance of grid box center to its boundary along xy diagonal
     1428    REAL(wp) :: distance_corners_down  !< distance of grid box center to its upper corners
     1429    REAL(wp) :: distance_corners_up    !< distance of grid box center to its lower corners
     1430    REAL(wp) :: radius                 !< search radius in meter
     1431
     1432    REAL(wp), DIMENSION(nzb:nzt+1) ::  gridsize_geometric_mean  !< geometric mean of grid sizes dx, dy, dz
     1433
     1434    ! ALLOCATE( gridsize_geometric_mean(nzb:nzt+1) )
     1435    ALLOCATE( ml_wall_adjusted(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    14121436!
    14131437!-- Initialize the mixing length in case of an LES-simulation
    14141438    IF ( .NOT. rans_mode )  THEN
    14151439!
    1416 !--    Compute the grid-dependent mixing length.
    1417        DO  k = 1, nzt
    1418           l_grid(k)  = ( dx * dy * dzw(k) )**0.33333333333333_wp
     1440!--    Compute the geometric averaged grid size (used for limiting the mixing length)
     1441       DO  k = nzb+1, nzt
     1442          gridsize_geometric_mean(k)  = ( dx * dy * dzw(k) )**0.33333333333333_wp
    14191443       ENDDO
    1420 !
    1421 !--    Initialize near-wall mixing length l_wall only in the vertical direction
    1422 !--    for the moment, multiplication with wall_adjustment_factor further below
    1423        l_wall(nzb,:,:)   = l_grid(1)
    1424        DO  k = nzb+1, nzt
    1425           l_wall(k,:,:)  = l_grid(k)
     1444       gridsize_geometric_mean(nzb)   = gridsize_geometric_mean(1)
     1445       gridsize_geometric_mean(nzt+1) = gridsize_geometric_mean(nzt)
     1446
     1447       IF ( ANY( gridsize_geometric_mean > 1.5_wp * dx * wall_adjustment_factor ) .OR. &
     1448            ANY( gridsize_geometric_mean > 1.5_wp * dy * wall_adjustment_factor ) )  THEN
     1449            WRITE( message_string, * ) 'grid anisotropy exceeds threshold',    &
     1450                                      ' &starting from height level k = ', k, &
     1451                                      '.'
     1452           CALL message( 'init_grid', 'PA0202', 0, 1, 0, 6, 0 )
     1453       ENDIF
     1454
     1455!
     1456!--    Initialize near-wall mixing length ml_wall_adjusted
     1457       DO  k = nzb, nzt+1
     1458          ml_wall_adjusted(k,:,:) = gridsize_geometric_mean(k)
    14261459       ENDDO
    1427        l_wall(nzt+1,:,:) = l_grid(nzt)
    14281460
    14291461       IF ( wall_adjustment )  THEN
    1430 
    1431           DO  k = 1, nzt
    1432              IF ( l_grid(k) > 1.5_wp * dx * wall_adjustment_factor .OR.            &
    1433                   l_grid(k) > 1.5_wp * dy * wall_adjustment_factor )  THEN
    1434                 WRITE( message_string, * ) 'grid anisotropy exceeds ',             &
    1435                                            'threshold given by only local',        &
    1436                                            ' &horizontal reduction of near_wall ', &
    1437                                            'mixing length l_wall',                 &
    1438                                            ' &starting from height level k = ', k, &
    1439                                            '.'
    1440                 CALL message( 'init_grid', 'PA0202', 0, 1, 0, 6, 0 )
    1441                 EXIT
    1442              ENDIF
    1443           ENDDO
    1444 !
    1445 !--       In case of topography: limit near-wall mixing length l_wall further:
    1446 !--       Go through all points of the subdomain one by one and look for the closest
    1447 !--       surface.
    1448 !--       Is this correct in the ocean case?
     1462!
     1463!--       In case of topography, adjust mixing length if there is any wall at
     1464!--       the surrounding grid boxes:
     1465          !> @todo check if this is correct also for the ocean case
    14491466          DO  i = nxl, nxr
    14501467             DO  j = nys, nyn
     
    14541471                   IF ( BTEST( wall_flags_total_0(k,j,i), 0 ) )  THEN
    14551472!
    1456 !--                   Check for neighbouring grid-points.
    1457 !--                   Vertical distance, down
    1458                       IF ( .NOT. BTEST( wall_flags_total_0(k-1,j,i), 0 ) )             &
    1459                          l_wall(k,j,i) = MIN( l_grid(k), zu(k) - zw(k-1) )
    1460 !
    1461 !--                   Vertical distance, up
    1462                       IF ( .NOT. BTEST( wall_flags_total_0(k+1,j,i), 0 ) )             &
    1463                          l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), zw(k) - zu(k) )
    1464 !
    1465 !--                   y-distance
    1466                       IF ( .NOT. BTEST( wall_flags_total_0(k,j-1,i), 0 )  .OR.         &
    1467                            .NOT. BTEST( wall_flags_total_0(k,j+1,i), 0 ) )             &
    1468                          l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), 0.5_wp * dy )
    1469 !
    1470 !--                   x-distance
    1471                       IF ( .NOT. BTEST( wall_flags_total_0(k,j,i-1), 0 )  .OR.         &
    1472                            .NOT. BTEST( wall_flags_total_0(k,j,i+1), 0 ) )             &
    1473                          l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), 0.5_wp * dx )
    1474 !
    1475 !--                   yz-distance (vertical edges, down)
    1476                       IF ( .NOT. BTEST( wall_flags_total_0(k-1,j-1,i), 0 )  .OR.       &
    1477                            .NOT. BTEST( wall_flags_total_0(k-1,j+1,i), 0 )  )          &
    1478                          l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),                &
    1479                                               SQRT( 0.25_wp * dy**2 +                  &
    1480                                              ( zu(k) - zw(k-1) )**2 ) )
    1481 !
    1482 !--                   yz-distance (vertical edges, up)
    1483                       IF ( .NOT. BTEST( wall_flags_total_0(k+1,j-1,i), 0 )  .OR.       &
    1484                            .NOT. BTEST( wall_flags_total_0(k+1,j+1,i), 0 )  )          &
    1485                          l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),                &
    1486                                               SQRT( 0.25_wp * dy**2 +                  &
    1487                                              ( zw(k) - zu(k) )**2 ) )
    1488 !
    1489 !--                   xz-distance (vertical edges, down)
    1490                       IF ( .NOT. BTEST( wall_flags_total_0(k-1,j,i-1), 0 )  .OR.       &
    1491                            .NOT. BTEST( wall_flags_total_0(k-1,j,i+1), 0 )  )          &
    1492                          l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),                &
    1493                                               SQRT( 0.25_wp * dx**2 +                  &
    1494                                              ( zu(k) - zw(k-1) )**2 ) )
    1495 !
    1496 !--                   xz-distance (vertical edges, up)
    1497                       IF ( .NOT. BTEST( wall_flags_total_0(k+1,j,i-1), 0 )  .OR.       &
    1498                            .NOT. BTEST( wall_flags_total_0(k+1,j,i+1), 0 )  )          &
    1499                          l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),                &
    1500                                               SQRT( 0.25_wp * dx**2 +                  &
    1501                                              ( zw(k) - zu(k) )**2 ) )
    1502 !
    1503 !--                   xy-distance (horizontal edges)
    1504                       IF ( .NOT. BTEST( wall_flags_total_0(k,j-1,i-1), 0 )  .OR.        &
    1505                            .NOT. BTEST( wall_flags_total_0(k,j+1,i-1), 0 )  .OR.        &
    1506                            .NOT. BTEST( wall_flags_total_0(k,j-1,i+1), 0 )  .OR.        &
    1507                            .NOT. BTEST( wall_flags_total_0(k,j+1,i+1), 0 ) )            &
    1508                          l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),                 &
    1509                                               SQRT( 0.25_wp * ( dx**2 + dy**2 ) ) )
    1510 !
    1511 !--                   xyz distance (vertical and horizontal edges, down)
    1512                       IF ( .NOT. BTEST( wall_flags_total_0(k-1,j-1,i-1), 0 )  .OR.      &
    1513                            .NOT. BTEST( wall_flags_total_0(k-1,j+1,i-1), 0 )  .OR.      &
    1514                            .NOT. BTEST( wall_flags_total_0(k-1,j-1,i+1), 0 )  .OR.      &
    1515                            .NOT. BTEST( wall_flags_total_0(k-1,j+1,i+1), 0 ) )          &
    1516                          l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),                 &
    1517                                               SQRT( 0.25_wp * ( dx**2 + dy**2 )         &
    1518                                                     +  ( zu(k) - zw(k-1) )**2  ) )
    1519 !
    1520 !--                   xyz distance (vertical and horizontal edges, up)
    1521                       IF ( .NOT. BTEST( wall_flags_total_0(k+1,j-1,i-1), 0 )  .OR.      &
    1522                            .NOT. BTEST( wall_flags_total_0(k+1,j+1,i-1), 0 )  .OR.      &
    1523                            .NOT. BTEST( wall_flags_total_0(k+1,j-1,i+1), 0 )  .OR.      &
    1524                            .NOT. BTEST( wall_flags_total_0(k+1,j+1,i+1), 0 ) )          &
    1525                          l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),                 &
    1526                                               SQRT( 0.25_wp * ( dx**2 + dy**2 )         &
    1527                                                     +  ( zw(k) - zu(k) )**2  ) )
    1528 
    1529                    ENDIF
    1530 !
    1531 !--                Adjust mixing length by wall-adjustment factor and limit it by l_grid
    1532                    l_wall(k,j,i) = MIN( l_wall(k,j,i) * wall_adjustment_factor, l_grid(k) )
     1473!--                   First, check if grid points directly next to current grid point
     1474!--                   are surface grid points
     1475!--                   Check along...
     1476!--                   ...vertical direction, down
     1477                      IF ( .NOT. BTEST( wall_flags_total_0(k-1,j,i), 0 ) )  THEN
     1478                         distance_down =  zu(k) - zw(k-1)
     1479                      ELSE
     1480                         distance_down = 9999999.9_wp
     1481                      ENDIF
     1482!
     1483!--                   ...vertical direction, up
     1484                      IF ( .NOT. BTEST( wall_flags_total_0(k+1,j,i), 0 ) )  THEN
     1485                         distance_up = zw(k) - zu(k)
     1486                      ELSE
     1487                         distance_up = 9999999.9_wp
     1488                      ENDIF
     1489!
     1490!--                   ...y-direction (north/south)
     1491                      IF ( .NOT. BTEST( wall_flags_total_0(k,j-1,i), 0 )  .OR. &
     1492                           .NOT. BTEST( wall_flags_total_0(k,j+1,i), 0 ) )  THEN
     1493                         distance_ns = 0.5_wp * dy
     1494                      ELSE
     1495                         distance_ns = 9999999.9_wp
     1496                      ENDIF
     1497!
     1498!--                   ...x-direction (left/right)
     1499                      IF ( .NOT. BTEST( wall_flags_total_0(k,j,i-1), 0 )  .OR. &
     1500                           .NOT. BTEST( wall_flags_total_0(k,j,i+1), 0 ) )  THEN
     1501                         distance_lr = 0.5_wp * dx
     1502                      ELSE
     1503                         distance_lr = 9999999.9_wp
     1504                      ENDIF
     1505!
     1506!--                   Now, check the edges along...
     1507!--                   ...yz-direction (vertical edges, down)
     1508                      IF ( .NOT. BTEST( wall_flags_total_0(k-1,j-1,i), 0 )  .OR.  &
     1509                           .NOT. BTEST( wall_flags_total_0(k-1,j+1,i), 0 )  )  THEN
     1510                         distance_edge_yz_down = SQRT( 0.25_wp * dy**2 + ( zu(k) - zw(k-1) )**2 )
     1511                      ELSE
     1512                         distance_edge_yz_down = 9999999.9_wp
     1513                      ENDIF
     1514!
     1515!--                   ...yz-direction (vertical edges, up)
     1516                      IF ( .NOT. BTEST( wall_flags_total_0(k+1,j-1,i), 0 )  .OR.  &
     1517                           .NOT. BTEST( wall_flags_total_0(k+1,j+1,i), 0 )  )  THEN
     1518                         distance_edge_yz_up = SQRT( 0.25_wp * dy**2 + ( zw(k) - zu(k) )**2 )
     1519                      ELSE
     1520                         distance_edge_yz_up = 9999999.9_wp
     1521                      ENDIF
     1522!
     1523!--                   ...xz-direction (vertical edges, down)
     1524                      IF ( .NOT. BTEST( wall_flags_total_0(k-1,j,i-1), 0 )  .OR.  &
     1525                           .NOT. BTEST( wall_flags_total_0(k-1,j,i+1), 0 )  )  THEN
     1526                         distance_edge_xz_down = SQRT( 0.25_wp * dx**2 + ( zu(k) - zw(k-1) )**2 )
     1527                      ELSE
     1528                         distance_edge_xz_down = 9999999.9_wp
     1529                      ENDIF
     1530!
     1531!--                   ...xz-direction (vertical edges, up)
     1532                      IF ( .NOT. BTEST( wall_flags_total_0(k+1,j,i-1), 0 )  .OR.  &
     1533                           .NOT. BTEST( wall_flags_total_0(k+1,j,i+1), 0 )  )  THEN
     1534                         distance_edge_xz_up = SQRT( 0.25_wp * dx**2 + ( zw(k) - zu(k) )**2 )
     1535                      ELSE
     1536                         distance_edge_xz_up = 9999999.9_wp
     1537                      ENDIF
     1538!
     1539!--                   ...xy-direction (horizontal edges)
     1540                      IF ( .NOT. BTEST( wall_flags_total_0(k,j-1,i-1), 0 )  .OR. &
     1541                           .NOT. BTEST( wall_flags_total_0(k,j+1,i-1), 0 )  .OR. &
     1542                           .NOT. BTEST( wall_flags_total_0(k,j-1,i+1), 0 )  .OR. &
     1543                           .NOT. BTEST( wall_flags_total_0(k,j+1,i+1), 0 ) )  THEN
     1544                         distance_edge_xy = SQRT( 0.25_wp * ( dx**2 + dy**2 ) )
     1545                      ELSE
     1546                         distance_edge_xy = 9999999.9_wp
     1547                      ENDIF
     1548!
     1549!--                   Now, check the corners...
     1550!--                   ...lower four corners
     1551                      IF ( .NOT. BTEST( wall_flags_total_0(k-1,j-1,i-1), 0 )  .OR. &
     1552                           .NOT. BTEST( wall_flags_total_0(k-1,j+1,i-1), 0 )  .OR. &
     1553                           .NOT. BTEST( wall_flags_total_0(k-1,j-1,i+1), 0 )  .OR. &
     1554                           .NOT. BTEST( wall_flags_total_0(k-1,j+1,i+1), 0 ) )  THEN
     1555                         distance_corners_down = SQRT( 0.25_wp * ( dx**2 + dy**2 ) &
     1556                                                       + ( zu(k) - zw(k-1) )**2 )
     1557                      ELSE
     1558                         distance_corners_down = 9999999.9_wp
     1559                      ENDIF
     1560!
     1561!--                   ...upper four corners
     1562                      IF ( .NOT. BTEST( wall_flags_total_0(k+1,j-1,i-1), 0 )  .OR. &
     1563                           .NOT. BTEST( wall_flags_total_0(k+1,j+1,i-1), 0 )  .OR. &
     1564                           .NOT. BTEST( wall_flags_total_0(k+1,j-1,i+1), 0 )  .OR. &
     1565                           .NOT. BTEST( wall_flags_total_0(k+1,j+1,i+1), 0 ) )  THEN
     1566                         distance_corners_up = SQRT( 0.25_wp * ( dx**2 + dy**2 )   &
     1567                                                     + ( zw(k) - zu(k) )**2 )
     1568                      ELSE
     1569                         distance_corners_up = 9999999.9_wp
     1570                      ENDIF
     1571!
     1572!--                   Adjust mixing length by wall-adjustment factor
     1573                      ml_wall_adjusted(k,j,i) = wall_adjustment_factor * MIN(  &
     1574                         distance_up, distance_down, distance_ns, distance_lr, &
     1575                         distance_edge_yz_down, distance_edge_yz_up,           &
     1576                         distance_edge_xz_down, distance_edge_xz_up,           &
     1577                         distance_edge_xy,                                     &
     1578                         distance_corners_down, distance_corners_up )
     1579
     1580                  ENDIF  !if grid point belongs to atmosphere
     1581
     1582                  ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(k,j,i), &
     1583                                                 gridsize_geometric_mean(k) )
    15331584
    15341585                ENDDO  !k loop
     
    15381589       ENDIF  !if wall_adjustment
    15391590
    1540     ELSE
     1591    ELSE  !--> RANS mode
    15411592!
    15421593!--    Initialize the mixing length in case of a RANS simulation
    1543        ALLOCATE( l_black(nzb:nzt+1) )
     1594       ALLOCATE( ml_blackadar(nzb:nzt+1) )
    15441595
    15451596!
    15461597!--    Calculate mixing length according to Blackadar (1962)
    15471598       IF ( f /= 0.0_wp )  THEN
    1548           l_max = 2.7E-4_wp * SQRT( ug(nzt+1)**2 + vg(nzt+1)**2 ) /            &
    1549                   ABS( f ) + 1.0E-10_wp
     1599          length_scale_max = 2.7E-4_wp * SQRT( ug(nzt+1)**2 + vg(nzt+1)**2 )   &
     1600                           / ABS( f ) + 1.0E-10_wp
    15501601       ELSE
    1551           l_max = 30.0_wp
     1602          length_scale_max = 30.0_wp
    15521603       ENDIF
    15531604
    15541605       DO  k = nzb, nzt
    1555           l_black(k) = kappa * zu(k) / ( 1.0_wp + kappa * zu(k) / l_max )
     1606          ml_blackadar(k) = kappa * zu(k) / ( 1.0_wp + kappa * zu(k) / length_scale_max )
    15561607       ENDDO
    15571608
    1558        l_black(nzt+1) = l_black(nzt)
     1609       ml_blackadar(nzt+1) = ml_blackadar(nzt)
    15591610
    15601611!
     
    15711622       ENDDO
    15721623
    1573        l_wall(nzb,:,:) = l_black(nzb)
    1574        l_wall(nzt+1,:,:) = l_black(nzt+1)
     1624       ml_wall_adjusted(nzb,:,:) = ml_blackadar(nzb)
     1625       ml_wall_adjusted(nzt+1,:,:) = ml_blackadar(nzt+1)
    15751626!
    15761627!--    Limit mixing length to either nearest wall or Blackadar mixing length.
     
    15811632       DO  k = nzb+1, nzt
    15821633
    1583           radius = l_black(k)  ! radius within walls are searched
    1584 !
    1585 !--       Set l_wall to its default maximum value (l_back)
    1586           l_wall(k,:,:) = radius
     1634          radius = ml_blackadar(k)  ! radius within walls are searched
     1635!
     1636!--       Set ml_wall_adjusted to its default maximum value (ml_blackadar)
     1637          ml_wall_adjusted(k,:,:) = radius
    15871638
    15881639!
     
    16791730!--                            Search for walls within octant (+++)
    16801731                               vic_yz = vicinity(0:rad_k+1,0:rad_j,ii)
    1681                                l_wall(k,j,i) = MIN( l_wall(k,j,i),             &
     1732                               ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(k,j,i),             &
    16821733                                       shortest_distance( vic_yz, .TRUE., ii ) )
    16831734!
     
    16871738!--                            shortest_distance").
    16881739                               vic_yz = vicinity(0:rad_k+1,0:-rad_j:-1,ii)
    1689                                l_wall(k,j,i) = MIN( l_wall(k,j,i),             &
     1740                               ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(k,j,i),             &
    16901741                                       shortest_distance( vic_yz, .TRUE., ii ) )
    16911742
     
    16941745!--                         Search for walls within octant (+--)
    16951746                            vic_yz = vicinity(0:-rad_k-1:-1,0:-rad_j:-1,ii)
    1696                             l_wall(k,j,i) = MIN( l_wall(k,j,i),                &
     1747                            ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(k,j,i),                &
    16971748                                      shortest_distance( vic_yz, .FALSE., ii ) )
    16981749!
    16991750!--                         Search for walls within octant (++-)
    17001751                            vic_yz = vicinity(0:-rad_k-1:-1,0:rad_j,ii)
    1701                             l_wall(k,j,i) = MIN( l_wall(k,j,i),                &
     1752                            ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(k,j,i),                &
    17021753                                      shortest_distance( vic_yz, .FALSE., ii ) )
    17031754!
    17041755!--                         Reduce search along x by already found distance
    1705                             dist_dx = CEILING( l_wall(k,j,i) / dx )
     1756                            dist_dx = CEILING( ml_wall_adjusted(k,j,i) / dx )
    17061757
    17071758                         ENDDO
     
    17161767!--                            Search for walls within octant (-++)
    17171768                               vic_yz = vicinity(0:rad_k+1,0:rad_j,ii)
    1718                                l_wall(k,j,i) = MIN( l_wall(k,j,i),             &
     1769                               ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(k,j,i),             &
    17191770                                      shortest_distance( vic_yz, .TRUE., -ii ) )
    17201771!
     
    17241775!--                            shortest_distance").
    17251776                               vic_yz = vicinity(0:rad_k+1,0:-rad_j:-1,ii)
    1726                                l_wall(k,j,i) = MIN( l_wall(k,j,i),             &
     1777                               ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(k,j,i),             &
    17271778                                      shortest_distance( vic_yz, .TRUE., -ii ) )
    17281779
     
    17311782!--                         Search for walls within octant (---)
    17321783                            vic_yz = vicinity(0:-rad_k-1:-1,0:-rad_j:-1,ii)
    1733                             l_wall(k,j,i) = MIN( l_wall(k,j,i),                &
     1784                            ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(k,j,i),                &
    17341785                                     shortest_distance( vic_yz, .FALSE., -ii ) )
    17351786!
    17361787!--                         Search for walls within octant (-+-)
    17371788                            vic_yz = vicinity(0:-rad_k-1:-1,0:rad_j,ii)
    1738                             l_wall(k,j,i) = MIN( l_wall(k,j,i),                &
     1789                            ml_wall_adjusted(k,j,i) = MIN( ml_wall_adjusted(k,j,i),                &
    17391790                                     shortest_distance( vic_yz, .FALSE., -ii ) )
    17401791!
    17411792!--                         Reduce search along x by already found distance
    1742                             dist_dx = CEILING( l_wall(k,j,i) / dx )
     1793                            dist_dx = CEILING( ml_wall_adjusted(k,j,i) / dx )
    17431794
    17441795                         ENDDO
     
    17481799                   ELSE  !Check if (i,j,k) belongs to atmosphere
    17491800
    1750                       l_wall(k,j,i) = l_black(k)
     1801                      ml_wall_adjusted(k,j,i) = ml_blackadar(k)
    17511802
    17521803                   ENDIF
     
    17621813       ENDDO  !k loop
    17631814
    1764        !$ACC ENTER DATA COPYIN(l_black(nzb:nzt+1))
     1815       !$ACC ENTER DATA COPYIN(ml_blackadar(nzb:nzt+1))
    17651816
    17661817    ENDIF  !LES or RANS mode
    17671818
    17681819!
    1769 !-- Set lateral boundary conditions for l_wall
    1770     CALL exchange_horiz( l_wall, nbgp )
    1771 
    1772     !$ACC ENTER DATA COPYIN(l_grid(nzb:nzt+1)) &
    1773     !$ACC COPYIN(l_wall(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
     1820!-- Set lateral boundary conditions for ml_wall_adjusted
     1821    CALL exchange_horiz( ml_wall_adjusted, nbgp )
     1822
     1823    !$ACC ENTER DATA COPYIN(ml_wall_adjusted(nzb:nzt+1,nysg:nyng,nxlg:nxrg))
    17741824
    17751825    CONTAINS
     
    22662316       IF ( plant_canopy )  CALL pcm_tendency( 6 )
    22672317
    2268 !       CALL user_actions( 'e-tendency' ) ToDo: find general solution for circular dependency between modules
     2318       !> @todo find general solution for circular dependency between modules
     2319       ! CALL user_actions( 'e-tendency' )
    22692320
    22702321!
     
    23772428!
    23782429!--    Additional sink term for flows through plant canopies
    2379 !        IF ( plant_canopy )  CALL pcm_tendency( ? )         !> @todo not yet implemented
    2380 
    2381 !       CALL user_actions( 'e-tendency' ) ToDo: find general solution for circular dependency between modules
     2430       ! IF ( plant_canopy )  CALL pcm_tendency( ? )     !> @todo not yet implemented
     2431
     2432       ! CALL user_actions( 'e-tendency' ) !> @todo find general solution for circular dependency between modules
    23822433
    23832434!
     
    24962547       IF ( plant_canopy )  CALL pcm_tendency( i, j, 6 )
    24972548
    2498 !       CALL user_actions( i, j, 'e-tendency' ) ToDo: find general solution for circular dependency between modules
     2549       ! CALL user_actions( i, j, 'e-tendency' ) !> @todo: find general solution for circular dependency between modules
    24992550
    25002551!
     
    25632614!        IF ( plant_canopy )  CALL pcm_tendency( i, j, ? )     !> @todo not yet implemented
    25642615
    2565 !       CALL user_actions( i, j, 'diss-tendency' ) ToDo: find general solution for circular dependency between modules
     2616!       CALL user_actions( i, j, 'diss-tendency' ) !> @todo: find general solution for circular dependency between modules
    25662617
    25672618!
     
    27262777
    27272778                   km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp         &
    2728                                    * 0.5_wp * dy
     2779                                      * 0.5_wp * dy
    27292780!
    27302781!--                -1.0 for right-facing wall, 1.0 for left-facing wall
     
    27452796
    27462797                   km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp         &
    2747                                    * 0.5_wp * dy
     2798                                      * 0.5_wp * dy
    27482799!
    27492800!--                -1.0 for right-facing wall, 1.0 for left-facing wall
     
    27642815
    27652816                   km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp         &
    2766                                    * 0.5_wp * dy
     2817                                      * 0.5_wp * dy
    27672818!
    27682819!--                -1.0 for right-facing wall, 1.0 for left-facing wall
     
    34103461
    34113462             km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp         &
    3412                              * 0.5_wp * dy
     3463                                * 0.5_wp * dy
    34133464!
    34143465!--          -1.0 for right-facing wall, 1.0 for left-facing wall
     
    34283479
    34293480             km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp         &
    3430                              * 0.5_wp * dy
     3481                                * 0.5_wp * dy
    34313482!
    34323483!--          -1.0 for right-facing wall, 1.0 for left-facing wall
     
    34463497
    34473498             km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp         &
    3448                              * 0.5_wp * dy
     3499                                * 0.5_wp * dy
    34493500!
    34503501!--          -1.0 for right-facing wall, 1.0 for left-facing wall
     
    39113962                flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_total_0(k,j,i),0) )
    39123963                tend(k,j,i) = tend(k,j,i) + flag * tmp_flux(k) * ( g / &
    3913                               MERGE( vpt_reference, vpt(k,j,i),          &
     3964                              MERGE( vpt_reference, vpt(k,j,i),        &
    39143965                                     use_single_reference_value ) )
    39153966             ENDDO
     
    39714022    REAL(wp)     ::  duv2_dz2       !< squared vertical gradient of wind vector
    39724023    REAL(wp)     ::  dvar_dz        !< vertical gradient of var
    3973     REAL(wp)     ::  l              !< mixing length
     4024    REAL(wp)     ::  ml             !< mixing length
    39744025    REAL(wp)     ::  var_reference  !< reference temperature
    39754026
    3976     REAL(wp), DIMENSION(nzb+1:nzt) ::  l_stable  !< mixing length according to stratification
    3977     REAL(wp), DIMENSION(nzb+1:nzt) ::  rif       !< Richardson flux number
     4027    REAL(wp), DIMENSION(nzb+1:nzt) ::  ml_stratification  !< mixing length according to stratification
     4028    REAL(wp), DIMENSION(nzb+1:nzt) ::  rif                !< Richardson flux number
    39784029
    39794030    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) ::  var  !< temperature
     
    39844035
    39854036       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
    3986        !$ACC PRIVATE(l, l_stable, dvar_dz) &
     4037       !$ACC PRIVATE(ml, ml_stratification, dvar_dz) &
    39874038       !$ACC PRESENT(diss, e, var, wall_flags_total_0) &
    3988        !$ACC PRESENT(dd2zu, l_grid, l_wall)
     4039       !$ACC PRESENT(dd2zu, ml_wall_adjusted)
    39894040       DO  i = nxl, nxr
    39904041          DO  j = nys, nyn
     
    39954046                IF ( dvar_dz > 0.0_wp ) THEN
    39964047                   IF ( use_single_reference_value )  THEN
    3997                      l_stable(k) = 0.76_wp * SQRT( e(k,j,i) )                  &
     4048                     ml_stratification(k) = 0.76_wp * SQRT( e(k,j,i) )         &
    39984049                                 / SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
    39994050                   ELSE
    4000                      l_stable(k) = 0.76_wp * SQRT( e(k,j,i) )               &
     4051                     ml_stratification(k) = 0.76_wp * SQRT( e(k,j,i) )         &
    40014052                                 / SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
    40024053                   ENDIF
    40034054                ELSE
    4004                    l_stable(k) = l_grid(k)
     4055                   ml_stratification(k) = ml_wall_adjusted(k,j,i)
    40054056                ENDIF
    40064057
     
    40114062             DO  k = nzb+1, nzt
    40124063
    4013                 l  = MIN( l_wall(k,j,i), l_stable(k) )
    4014 
    4015                 diss(k,j,i) = ( 0.19_wp + 0.74_wp * l / l_wall(k,j,i) )                &
    4016                               * e(k,j,i) * SQRT( e(k,j,i) ) / l                        &
     4064                ml = MIN( ml_wall_adjusted(k,j,i), ml_stratification(k) )
     4065
     4066                diss(k,j,i) = ( 0.19_wp + 0.74_wp * ml / ml_wall_adjusted(k,j,i) )           &
     4067                              * e(k,j,i) * SQRT( e(k,j,i) ) / ml                             &
    40174068                              * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
    40184069
     
    40254076
    40264077      !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
    4027       !$ACC PRIVATE(l_stable, duv2_dz2, rif, dvar_dz) &
     4078      !$ACC PRIVATE(ml_stratification, duv2_dz2, rif, dvar_dz) &
    40284079      !$ACC PRESENT(diss, e, u, v, var, wall_flags_total_0) &
    4029       !$ACC PRESENT(dd2zu, l_black, l_wall)
     4080      !$ACC PRESENT(dd2zu, ml_blackadar, ml_wall_adjusted)
    40304081       DO  i = nxl, nxr
    40314082          DO  j = nys, nyn
     
    40624113             DO  k = nzb+1, nzt
    40634114                IF ( rif(k) >= 0.0_wp )  THEN
    4064                    l_stable(k) = MIN( l_black(k) / ( 1.0_wp + 5.0_wp * rif(k) ), l_wall(k,j,i) )
     4115                   ml_stratification(k) = ml_blackadar(k) / ( 1.0_wp + 5.0_wp * rif(k) )
    40654116                ELSE
    4066                    l_stable(k) = l_wall(k,j,i) * SQRT( 1.0_wp - 16.0_wp * rif(k) )
     4117                   ml_stratification(k) = ml_blackadar(k) * SQRT( 1.0_wp - 16.0_wp * rif(k) )
    40674118                ENDIF
    40684119             ENDDO
     
    40724123             DO  k = nzb+1, nzt
    40734124
    4074                 diss(k,j,i) = c_0**3 * e(k,j,i) * SQRT( e(k,j,i) ) / l_stable(k)     &
     4125                diss(k,j,i) = c_0**3 * e(k,j,i) * SQRT( e(k,j,i) )                         &
     4126                            / MIN( ml_stratification(k), ml_wall_adjusted(k,j,i) )         &
    40754127                            * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
    40764128
     
    41804232    REAL(wp)     ::  duv2_dz2       !< squared vertical gradient of wind vector
    41814233    REAL(wp)     ::  dvar_dz        !< vertical gradient of var
    4182     REAL(wp)     ::  l              !< mixing length
     4234    REAL(wp)     ::  ml             !< mixing length
    41834235    REAL(wp)     ::  var_reference  !< reference temperature
    41844236
    4185     REAL(wp), DIMENSION(nzb+1:nzt) ::  l_stable  !< mixing length according to stratification
    4186     REAL(wp), DIMENSION(nzb+1:nzt) ::  rif       !< Richardson flux number
     4237    REAL(wp), DIMENSION(nzb+1:nzt) ::  ml_stratification  !< mixing length according to stratification
     4238    REAL(wp), DIMENSION(nzb+1:nzt) ::  rif                !< Richardson flux number
    41874239
    41884240    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) ::  var  !< temperature
     
    41974249          IF ( dvar_dz > 0.0_wp ) THEN
    41984250             IF ( use_single_reference_value )  THEN
    4199                 l_stable(k) = 0.76_wp * SQRT( e(k,j,i) )                  &
     4251                ml_stratification(k) = 0.76_wp * SQRT( e(k,j,i) )              &
    42004252                            / SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
    42014253             ELSE
    4202                 l_stable(k) = 0.76_wp * SQRT( e(k,j,i) )               &
     4254                ml_stratification(k) = 0.76_wp * SQRT( e(k,j,i) )              &
    42034255                            / SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
    42044256             ENDIF
    42054257          ELSE
    4206              l_stable(k) = l_grid(k)
     4258             ml_stratification(k) = ml_wall_adjusted(k,j,i)
    42074259          ENDIF
    42084260       ENDDO
     
    42104262       !DIR$ IVDEP
    42114263       DO  k = nzb+1, nzt
    4212           l  = MIN( l_wall(k,j,i), l_stable(k) )
    4213 
    4214           diss(k,j,i) = ( 0.19_wp + 0.74_wp * l / l_wall(k,j,i) )              &
    4215                         * e(k,j,i) * SQRT( e(k,j,i) ) / l                      &
     4264          ml  = MIN( ml_wall_adjusted(k,j,i), ml_stratification(k) )
     4265
     4266          diss(k,j,i) = ( 0.19_wp + 0.74_wp * ml / ml_wall_adjusted(k,j,i) )           &
     4267                        * e(k,j,i) * SQRT( e(k,j,i) ) / ml                             &
    42164268                        * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
    42174269       ENDDO
     
    42484300       DO  k = nzb+1, nzt
    42494301          IF ( rif(k) >= 0.0_wp )  THEN
    4250              l_stable(k) = MIN( l_black(k) / ( 1.0_wp + 5.0_wp * rif(k) ), l_wall(k,j,i) )
     4302             ml_stratification(k) = ml_blackadar(k) / ( 1.0_wp + 5.0_wp * rif(k) )
    42514303          ELSE
    4252              l_stable(k) = l_wall(k,j,i) * SQRT( 1.0_wp - 16.0_wp * rif(k) )
     4304             ml_stratification(k) = ml_blackadar(k) * SQRT( 1.0_wp - 16.0_wp * rif(k) )
    42534305          ENDIF
    42544306
     
    42574309       !DIR$ IVDEP
    42584310       DO  k = nzb+1, nzt
    4259           diss(k,j,i) = c_0**3 * e(k,j,i) * SQRT( e(k,j,i) ) / l_stable(k)     &
     4311          diss(k,j,i) = c_0**3 * e(k,j,i) * SQRT( e(k,j,i) )                         &
     4312                      / MIN( ml_stratification(k), ml_wall_adjusted(k,j,i) )         &
    42604313                      * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
    42614314       ENDDO
     
    46524705    REAL(wp)     ::  duv2_dz2            !< squared vertical gradient of wind vector
    46534706    REAL(wp)     ::  dvar_dz             !< vertical gradient of var
    4654     REAL(wp)     ::  l                   !< mixing length (single height)
     4707    REAL(wp)     ::  ml                  !< mixing length (single height)
    46554708    REAL(wp)     ::  var_reference       !< reference temperature
    46564709
    4657     !DIR$ ATTRIBUTES ALIGN:64:: l_v, l_stable, rif
    4658     REAL(wp), DIMENSION(nzb+1:nzt) ::  l_v       !< mixing length (all heights)
    4659     REAL(wp), DIMENSION(nzb+1:nzt) ::  l_stable  !< mixing length according to stratification
    4660     REAL(wp), DIMENSION(nzb+1:nzt) ::  rif       !< Richardson flux number
     4710    !DIR$ ATTRIBUTES ALIGN:64:: ml_local_profile, ml_stratification, rif
     4711    REAL(wp), DIMENSION(nzb+1:nzt) ::  ml_local_profile   !< mixing length (all heights)
     4712    REAL(wp), DIMENSION(nzb+1:nzt) ::  ml_stratification  !< mixing length according to stratification
     4713    REAL(wp), DIMENSION(nzb+1:nzt) ::  rif                !< Richardson flux number
    46614714
    46624715    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) ::  var  !< temperature
     
    46744727!
    46754728!-- Compute the turbulent diffusion coefficient for momentum
    4676     !$OMP PARALLEL PRIVATE (i,j,k,l,sr,tn)
     4729    !$OMP PARALLEL PRIVATE (i,j,k,ml,sr,tn)
    46774730!$  tn = omp_get_thread_num()
    46784731
     
    46804733       !$OMP DO
    46814734       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
    4682        !$ACC PRIVATE(dvar_dz, l, l_stable, l_v) &
    4683        !$ACC PRESENT(wall_flags_total_0, var, dd2zu, e, l_wall, l_grid, rmask) &
    4684        !$ACC PRESENT(kh, km, sums_l_l)
     4735       !$ACC PRIVATE(dvar_dz, ml, ml_stratification, ml_local_profile) &
     4736       !$ACC PRESENT(wall_flags_total_0, var, dd2zu, e, ml_wall_adjusted) &
     4737       !$ACC PRESENT(kh, km, sums_l_l, rmask)
    46854738       DO  i = nxlg, nxrg
    46864739          DO  j = nysg, nyng
     
    46964749                IF ( dvar_dz > 0.0_wp ) THEN
    46974750                   IF ( use_single_reference_value )  THEN
    4698                       l_stable(k) = 0.76_wp * SQRT( e(k,j,i) )                  &
     4751                      ml_stratification(k) = 0.76_wp * SQRT( e(k,j,i) )         &
    46994752                                  / SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
    47004753                   ELSE
    4701                       l_stable(k) = 0.76_wp * SQRT( e(k,j,i) )               &
     4754                      ml_stratification(k) = 0.76_wp * SQRT( e(k,j,i) )         &
    47024755                                  / SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
    47034756                   ENDIF
    47044757                ELSE
    4705                    l_stable(k) = l_grid(k)
     4758                   ml_stratification(k) = ml_wall_adjusted(k,j,i)
    47064759                ENDIF
    47074760
     
    47124765             DO  k = nzb+1, nzt
    47134766
    4714                 l_v(k) = MIN( l_wall(k,j,i), l_stable(k) )                      &
    4715                        * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
    4716                 l = l_v(k)
     4767                ml = MIN( ml_wall_adjusted(k,j,i), ml_stratification(k) )         &
     4768                   * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
     4769                ml_local_profile(k) = ml
    47174770!
    47184771!--             Compute diffusion coefficients for momentum and heat
    4719                 km(k,j,i) = c_0 * l * SQRT( e(k,j,i) )
    4720                 kh(k,j,i) = ( 1.0_wp + 2.0_wp * l / l_wall(k,j,i) ) * km(k,j,i)
     4772                km(k,j,i) = c_0 * ml * SQRT( e(k,j,i) )
     4773                kh(k,j,i) = ( 1.0_wp + 2.0_wp * ml / ml_wall_adjusted(k,j,i) ) * km(k,j,i)
    47214774
    47224775             ENDDO
     
    47264779             DO  sr = 0, statistic_regions
    47274780                DO  k = nzb+1, nzt
    4728                    sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l_v(k) * rmask(j,i,sr)
     4781                   sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + ml_local_profile(k) * rmask(j,i,sr)
    47294782                ENDDO
    47304783             ENDDO
     
    47374790       !$OMP DO
    47384791       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
    4739        !$ACC PRIVATE(dvar_dz, duv2_dz2, l_stable, l_v, rif) &
    4740        !$ACC PRESENT(wall_flags_total_0, var, dd2zu, e, u, v, l_wall, l_black, rmask) &
    4741        !$ACC PRESENT(kh, km, sums_l_l)
     4792       !$ACC PRIVATE(dvar_dz, duv2_dz2, ml_stratification, ml_local_profile, rif) &
     4793       !$ACC PRESENT(wall_flags_total_0, var, dd2zu, e, u, v, ml_wall_adjusted, ml_blackadar) &
     4794       !$ACC PRESENT(kh, km, sums_l_l, rmask)
    47424795       DO  i = nxlg, nxrg
    47434796          DO  j = nysg, nyng
     
    47754828             DO  k = nzb+1, nzt
    47764829                IF ( rif(k) >= 0.0_wp )  THEN
    4777                    l_stable(k)  = MIN( l_black(k) / ( 1.0_wp + 5.0_wp * rif(k) ), l_wall(k,j,i) )
     4830                   ml_stratification(k) = ml_blackadar(k) / ( 1.0_wp + 5.0_wp * rif(k) )
    47784831                ELSE
    4779                    l_stable(k)  = l_wall(k,j,i)
     4832                   ml_stratification(k) = ml_blackadar(k)
    47804833                ENDIF
    47814834
     
    47864839             !DIR$ IVDEP
    47874840             DO  k = nzb+1, nzt
    4788                 l_v(k)    = l_stable(k) * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
    4789                 km(k,j,i) = c_0 * l_v(k) * SQRT( e(k,j,i) )
     4841                ml_local_profile(k) = MIN( ml_stratification(k), ml_wall_adjusted(k,j,i) ) &
     4842                                    * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
     4843                km(k,j,i) = c_0 * ml_local_profile(k) * SQRT( e(k,j,i) )
    47904844                kh(k,j,i) = km(k,j,i) / prandtl_number
    47914845             ENDDO
     
    47954849             DO  sr = 0, statistic_regions
    47964850                DO  k = nzb+1, nzt
    4797                    sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l_v(k) * rmask(j,i,sr)
     4851                   sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + ml_local_profile(k) * rmask(j,i,sr)
    47984852                ENDDO
    47994853             ENDDO
     
    48064860       !$OMP DO
    48074861       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
    4808        !$ACC PRIVATE(l_v) &
     4862       !$ACC PRIVATE(ml_local_profile) &
    48094863       !$ACC PRESENT(wall_flags_total_0, e, diss, rmask) &
    48104864       !$ACC PRESENT(kh, km, sums_l_l)
     
    48174871             DO  k = nzb+1, nzt
    48184872
    4819                 l_v(k) = c_0**3 * e(k,j,i) * SQRT(e(k,j,i)) / ( diss(k,j,i) + 1.0E-30_wp ) &
    4820                        * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
    4821 
    4822                 km(k,j,i) = c_0 * SQRT( e(k,j,i) ) * l_v(k)
     4873                ml_local_profile(k) = c_0**3 * e(k,j,i) * SQRT(e(k,j,i))                           &
     4874                                    / ( diss(k,j,i) + 1.0E-30_wp )                                 &
     4875                                    * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) )
     4876
     4877                km(k,j,i) = c_0 * SQRT( e(k,j,i) ) * ml_local_profile(k)
    48234878                kh(k,j,i) = km(k,j,i) / prandtl_number
    48244879
     
    48294884             DO  sr = 0, statistic_regions
    48304885                DO  k = nzb+1, nzt
    4831                    sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l_v(k) * rmask(j,i,sr)
     4886                   sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + ml_local_profile(k) * rmask(j,i,sr)
    48324887                ENDDO
    48334888             ENDDO
Note: See TracChangeset for help on using the changeset viewer.