Changeset 4170 for palm/trunk


Ignore:
Timestamp:
Aug 19, 2019 5:12:31 PM (5 years ago)
Author:
gronemeier
Message:

changes in turbulence_closure_mod:

  • add performance optimizations according to K. Ketelsen to diffusion_e and tcm_diffusivities_default
  • bugfix in calculating l_wall for vertical walls
  • bugfix in using l_wall in initialization (consider wall_adjustment_factor)
  • always initialize diss and save the dissipation to that array

related changes in time_integration:

  • copy diss, diss_p, tdiss_m to GPU
Location:
palm/trunk/SOURCE
Files:
2 edited

Legend:

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

    r4144 r4170  
    2525! -----------------
    2626! $Id$
     27! copy diss, diss_p, tdiss_m to GPU
     28!
     29! 4144 2019-08-06 09:11:47Z raasch
    2730! relational operators .EQ., .NE., etc. replaced by ==, /=, etc.
    2831!
     
    718721        ONLY:  d, dd2zu, ddzu, ddzw, drho_air, drho_air_zw, dzw, heatflux_output_conversion, kh,   &
    719722               km, momentumflux_output_conversion, p, ptdf_x, ptdf_y, rdf, rdf_sc, rho_air,        &
    720                rho_air_zw, te_m, tpt_m, tu_m, tv_m, tw_m, ug, u_init, u_stokes_zu, vg, v_init,     &
    721                v_stokes_zu, zu
     723               rho_air_zw, tdiss_m, te_m, tpt_m, tu_m, tv_m, tw_m, ug, u_init, u_stokes_zu, vg,    &
     724               v_init, v_stokes_zu, zu
    722725
    723726    USE control_parameters,                                                                        &
     
    761764!$ACC DATA &
    762765!$ACC COPY(d(nzb+1:nzt,nys:nyn,nxl:nxr)) &
     766!$ACC COPY(diss(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
    763767!$ACC COPY(e(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
    764768!$ACC COPY(u(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
     
    771775
    772776!$ACC DATA &
     777!$ACC COPY(diss_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
    773778!$ACC COPY(e_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
    774779!$ACC COPY(u_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
     
    777782!$ACC COPY(pt_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
    778783!$ACC COPY(tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
     784!$ACC COPY(tdiss_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
    779785!$ACC COPY(te_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
    780786!$ACC COPY(tu_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) &
  • palm/trunk/SOURCE/turbulence_closure_mod.f90

    r4168 r4170  
    2525! -----------------
    2626! $Id$
     27! - add performance optimizations according to K. Ketelsen
     28!   to diffusion_e and tcm_diffusivities_default
     29! - bugfix in calculating l_wall for vertical walls
     30! - bugfix in using l_wall in initialization (consider wall_adjustment_factor)
     31! - always initialize diss and save the dissipation to that array
     32!
     33! 4168 2019-08-16 13:50:17Z suehring
    2734! Replace function get_topography_top_index by topo_top_ind
    2835!
     
    205212!>
    206213!> @todo test initialization for all possibilities
    207 !>       add OpenMP directives whereever possible
     214!> @todo add OpenMP directives whereever possible
    208215!> @todo Check for random disturbances
    209216!> @note <Enter notes on the module>
     
    657664       ENDDO
    658665
    659        
    660666       diss_p(nzt+1,:,:) = diss_p(nzt,:,:)
    661      
     667
    662668    ENDIF
     669
    663670 END SUBROUTINE tcm_boundary_conds
    664671 
     
    11601167        ONLY:  collision_turbulence
    11611168
    1162     USE particle_attributes,                                                   &
    1163         ONLY:  use_sgs_for_particles, wang_kernel
    1164 
    11651169    USE pmc_interface,                                                         &
    11661170        ONLY:  nested_run
     
    11801184!-- they do not necessarily need to be transferred, which is attributed to
    11811185!-- the design of the model coupler which allocates memory for each variable.
    1182     IF ( rans_mode  .OR.  use_sgs_for_particles  .OR.  wang_kernel  .OR.       &
    1183          collision_turbulence  .OR.  nested_run )  THEN
    1184 
    1185        ALLOCATE( diss_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    1186        IF ( rans_tke_e  .OR.  nested_run )  THEN
    1187           ALLOCATE( diss_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    1188           ALLOCATE( diss_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    1189        ENDIF
    1190 
     1186    ALLOCATE( diss_1(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     1187
     1188    IF ( rans_tke_e  .OR.  nested_run )  THEN
     1189       ALLOCATE( diss_2(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
     1190       ALLOCATE( diss_3(nzb:nzt+1,nysg:nyng,nxlg:nxrg) )
    11911191    ENDIF
    11921192
     
    11951195    e  => e_1;   e_p  => e_2;   te_m  => e_3
    11961196
    1197     IF ( rans_mode  .OR.  use_sgs_for_particles  .OR.     &
    1198          wang_kernel  .OR.  collision_turbulence  .OR.  nested_run )  THEN
    1199        diss => diss_1
    1200        IF ( rans_tke_e  .OR.  nested_run )  THEN
     1197    diss => diss_1
     1198    IF ( rans_tke_e  .OR.  nested_run )  THEN
    12011199       diss_p => diss_2; tdiss_m => diss_3
    1202        ENDIF
    12031200    ENDIF
    12041201
     
    12521249                e = 0.0_wp
    12531250             ENDIF
     1251
     1252             diss = 0.0_wp
    12541253
    12551254          ELSE
     
    12981297          ENDIF
    12991298
    1300        ELSEIF ( INDEX(initializing_actions, 'set_constant_profiles') /= 0 .OR. &
     1299       ELSEIF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0 .OR. &
    13011300                INDEX( initializing_actions, 'inifor' ) /= 0 )  THEN
    13021301
     
    13401339             diss(nzb,:,:) = diss(nzb+1,:,:)
    13411340             diss(nzt+1,:,:) = diss(nzt,:,:)
     1341          ELSE
     1342             diss = 0.0_wp
    13421343          ENDIF
    13431344
     
    13831384                ENDDO
    13841385             ENDDO
     1386          ELSE
     1387             diss = 0.0_wp
    13851388          ENDIF
    13861389       ENDIF
     
    14691472
    14701473
     1474!------------------------------------------------------------------------------!
    14711475! Description:
    1472 ! -----------------------------------------------------------------------------!
     1476! ------------
    14731477!> Pre-computation of grid-dependent and near-wall mixing length.
    14741478!> @todo consider walls in horizontal direction at a distance further than a
     
    14811485
    14821486    USE control_parameters,                                                    &
    1483         ONLY:  f, message_string, wall_adjustment_factor
     1487        ONLY:  f, message_string, wall_adjustment, wall_adjustment_factor
    14841488
    14851489    USE grid_variables,                                                        &
     
    14951499    IMPLICIT NONE
    14961500
    1497     INTEGER(iwp) :: dist_dx        !< found distance devided by dx
    1498     INTEGER(iwp) :: i              !< index variable along x
    1499     INTEGER(iwp) :: ii             !< index variable along x
    1500     INTEGER(iwp) :: j              !< index variable along y
    1501     INTEGER(iwp) :: k              !< index variable along z
    1502     INTEGER(iwp) :: k_max_topo = 0 !< index of maximum topography height
    1503     INTEGER(iwp) :: kk             !< index variable along z
    1504     INTEGER(iwp) :: rad_i          !< search radius in grid points along x
    1505     INTEGER(iwp) :: rad_i_l        !< possible search radius to the left
    1506     INTEGER(iwp) :: rad_i_r        !< possible search radius to the right
    1507     INTEGER(iwp) :: rad_j          !< search radius in grid points along y
    1508     INTEGER(iwp) :: rad_j_n        !< possible search radius to north
    1509     INTEGER(iwp) :: rad_j_s        !< possible search radius to south
    1510     INTEGER(iwp) :: rad_k          !< search radius in grid points along z
    1511     INTEGER(iwp) :: rad_k_b        !< search radius in grid points along negative z
    1512     INTEGER(iwp) :: rad_k_t        !< search radius in grid points along positive z
     1501    INTEGER(iwp) :: dist_dx     !< found distance devided by dx
     1502    INTEGER(iwp) :: i           !< index variable along x
     1503    INTEGER(iwp) :: ii          !< index variable along x
     1504    INTEGER(iwp) :: j           !< index variable along y
     1505    INTEGER(iwp) :: k           !< index variable along z
     1506    INTEGER(iwp) :: k_max_topo !< index of maximum topography height
     1507    INTEGER(iwp) :: kk          !< index variable along z
     1508    INTEGER(iwp) :: rad_i       !< search radius in grid points along x
     1509    INTEGER(iwp) :: rad_i_l     !< possible search radius to the left
     1510    INTEGER(iwp) :: rad_i_r     !< possible search radius to the right
     1511    INTEGER(iwp) :: rad_j       !< search radius in grid points along y
     1512    INTEGER(iwp) :: rad_j_n     !< possible search radius to north
     1513    INTEGER(iwp) :: rad_j_s     !< possible search radius to south
     1514    INTEGER(iwp) :: rad_k       !< search radius in grid points along z
     1515    INTEGER(iwp) :: rad_k_b     !< search radius in grid points along negative z
     1516    INTEGER(iwp) :: rad_k_t     !< search radius in grid points along positive z
    15131517
    15141518    INTEGER(KIND=1), DIMENSION(:,:), ALLOCATABLE :: vic_yz !< contains a quarter of a single yz-slice of vicinity
     
    15161520    INTEGER(KIND=1), DIMENSION(:,:,:), ALLOCATABLE :: vicinity !< contains topography information of the vicinity of (i/j/k)
    15171521
    1518     REAL(wp) :: radius           !< search radius in meter
     1522    REAL(wp) :: radius          !< search radius in meter
    15191523
    15201524    ALLOCATE( l_grid(1:nzt) )
     
    15371541       l_wall(nzt+1,:,:) = l_grid(nzt)
    15381542
    1539        DO  k = 1, nzt
    1540           IF ( l_grid(k) > 1.5_wp * dx * wall_adjustment_factor .OR.            &
    1541                l_grid(k) > 1.5_wp * dy * wall_adjustment_factor )  THEN
    1542              WRITE( message_string, * ) 'grid anisotropy exceeds ',             &
    1543                                         'threshold given by only local',        &
    1544                                         ' &horizontal reduction of near_wall ', &
    1545                                         'mixing length l_wall',                 &
    1546                                         ' &starting from height level k = ', k, &
    1547                                         '.'
    1548              CALL message( 'init_grid', 'PA0202', 0, 1, 0, 6, 0 )
    1549              EXIT
    1550           ENDIF
    1551        ENDDO
    1552 !
    1553 !--    In case of topography: limit near-wall mixing length l_wall further:
    1554 !--    Go through all points of the subdomain one by one and look for the closest
    1555 !--    surface.
    1556 !--    Is this correct in the ocean case?
    1557        DO  i = nxl, nxr
    1558           DO  j = nys, nyn
    1559              DO  k = nzb+1, nzt
    1560 !
    1561 !--             Check if current gridpoint belongs to the atmosphere
    1562                 IF ( BTEST( wall_flags_0(k,j,i), 0 ) )  THEN
    1563 !
    1564 !--                Check for neighbouring grid-points.
    1565 !--                Vertical distance, down
    1566                    IF ( .NOT. BTEST( wall_flags_0(k-1,j,i), 0 ) )              &
    1567                       l_wall(k,j,i) = MIN( l_grid(k), zu(k) - zw(k-1) )
    1568 !
    1569 !--                Vertical distance, up
    1570                    IF ( .NOT. BTEST( wall_flags_0(k+1,j,i), 0 ) )              &
    1571                       l_wall(k,j,i) = MIN( l_grid(k), zw(k) - zu(k) )
    1572 !
    1573 !--                y-distance
    1574                    IF ( .NOT. BTEST( wall_flags_0(k,j-1,i), 0 )  .OR.          &
    1575                         .NOT. BTEST( wall_flags_0(k,j+1,i), 0 ) )              &
    1576                       l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), 0.5_wp * dy )
    1577 !
    1578 !--                x-distance
    1579                    IF ( .NOT. BTEST( wall_flags_0(k,j,i-1), 0 )  .OR.          &
    1580                         .NOT. BTEST( wall_flags_0(k,j,i+1), 0 ) )              &
    1581                       l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), 0.5_wp * dx )
    1582 !
    1583 !--                 yz-distance (vertical edges, down)
    1584                     IF ( .NOT. BTEST( wall_flags_0(k-1,j-1,i), 0 )  .OR.       &
    1585                          .NOT. BTEST( wall_flags_0(k-1,j+1,i), 0 )  )          &
    1586                       l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),           &
    1587                                            SQRT( 0.25_wp * dy**2 +             &
    1588                                           ( zu(k) - zw(k-1) )**2 ) )
    1589 !
    1590 !--                  yz-distance (vertical edges, up)
    1591                     IF ( .NOT. BTEST( wall_flags_0(k+1,j-1,i), 0 )  .OR.       &
    1592                          .NOT. BTEST( wall_flags_0(k+1,j+1,i), 0 )  )          &
    1593                       l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),           &
    1594                                            SQRT( 0.25_wp * dy**2 +             &
    1595                                           ( zw(k) - zu(k) )**2 ) )
    1596 !
    1597 !--                 xz-distance (vertical edges, down)
    1598                     IF ( .NOT. BTEST( wall_flags_0(k-1,j,i-1), 0 )  .OR.       &
    1599                          .NOT. BTEST( wall_flags_0(k-1,j,i+1), 0 )  )          &
    1600                       l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),           &
    1601                                            SQRT( 0.25_wp * dx**2 +             &
    1602                                           ( zu(k) - zw(k-1) )**2 ) )
    1603 !
    1604 !--                 xz-distance (vertical edges, up)
    1605                     IF ( .NOT. BTEST( wall_flags_0(k+1,j,i-1), 0 )  .OR.       &
    1606                          .NOT. BTEST( wall_flags_0(k+1,j,i+1), 0 )  )          &
    1607                      l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),            &
    1608                                            SQRT( 0.25_wp * dx**2 +             &
    1609                                           ( zw(k) - zu(k) )**2 ) )
    1610 !
    1611 !--                xy-distance (horizontal edges)
    1612                    IF ( .NOT. BTEST( wall_flags_0(k,j-1,i-1), 0 )  .OR.        &
    1613                         .NOT. BTEST( wall_flags_0(k,j+1,i-1), 0 )  .OR.        &
    1614                         .NOT. BTEST( wall_flags_0(k,j-1,i+1), 0 )  .OR.        &
    1615                         .NOT. BTEST( wall_flags_0(k,j+1,i+1), 0 ) )            &
    1616                       l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),           &
    1617                                            SQRT( 0.25_wp * ( dx**2 + dy**2 ) ) )
    1618 !
    1619 !--                xyz distance (vertical and horizontal edges, down)
    1620                    IF ( .NOT. BTEST( wall_flags_0(k-1,j-1,i-1), 0 )  .OR.      &
    1621                         .NOT. BTEST( wall_flags_0(k-1,j+1,i-1), 0 )  .OR.      &
    1622                         .NOT. BTEST( wall_flags_0(k-1,j-1,i+1), 0 )  .OR.      &
    1623                         .NOT. BTEST( wall_flags_0(k-1,j+1,i+1), 0 ) )          &
    1624                       l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),           &
    1625                                            SQRT( 0.25_wp * ( dx**2 + dy**2 )   &
    1626                                                  +  ( zu(k) - zw(k-1) )**2  ) )
    1627 !
    1628 !--                xyz distance (vertical and horizontal edges, up)
    1629                    IF ( .NOT. BTEST( wall_flags_0(k+1,j-1,i-1), 0 )  .OR.      &
    1630                         .NOT. BTEST( wall_flags_0(k+1,j+1,i-1), 0 )  .OR.      &
    1631                         .NOT. BTEST( wall_flags_0(k+1,j-1,i+1), 0 )  .OR.      &
    1632                         .NOT. BTEST( wall_flags_0(k+1,j+1,i+1), 0 ) )          &
    1633                       l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),           &
    1634                                            SQRT( 0.25_wp * ( dx**2 + dy**2 )   &
    1635                                                  +  ( zw(k) - zu(k) )**2  ) )
    1636 
    1637                 ENDIF
    1638              ENDDO
     1543       IF ( wall_adjustment )  THEN
     1544
     1545          DO  k = 1, nzt
     1546             IF ( l_grid(k) > 1.5_wp * dx * wall_adjustment_factor .OR.            &
     1547                  l_grid(k) > 1.5_wp * dy * wall_adjustment_factor )  THEN
     1548                WRITE( message_string, * ) 'grid anisotropy exceeds ',             &
     1549                                           'threshold given by only local',        &
     1550                                           ' &horizontal reduction of near_wall ', &
     1551                                           'mixing length l_wall',                 &
     1552                                           ' &starting from height level k = ', k, &
     1553                                           '.'
     1554                CALL message( 'init_grid', 'PA0202', 0, 1, 0, 6, 0 )
     1555                EXIT
     1556             ENDIF
    16391557          ENDDO
    1640        ENDDO
     1558!
     1559!--       In case of topography: limit near-wall mixing length l_wall further:
     1560!--       Go through all points of the subdomain one by one and look for the closest
     1561!--       surface.
     1562!--       Is this correct in the ocean case?
     1563          DO  i = nxl, nxr
     1564             DO  j = nys, nyn
     1565                DO  k = nzb+1, nzt
     1566!
     1567!--                Check if current gridpoint belongs to the atmosphere
     1568                   IF ( BTEST( wall_flags_0(k,j,i), 0 ) )  THEN
     1569!
     1570!--                   Check for neighbouring grid-points.
     1571!--                   Vertical distance, down
     1572                      IF ( .NOT. BTEST( wall_flags_0(k-1,j,i), 0 ) )             &
     1573                         l_wall(k,j,i) = MIN( l_grid(k), zu(k) - zw(k-1) )
     1574!
     1575!--                   Vertical distance, up
     1576                      IF ( .NOT. BTEST( wall_flags_0(k+1,j,i), 0 ) )             &
     1577                         l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), zw(k) - zu(k) )
     1578!
     1579!--                   y-distance
     1580                      IF ( .NOT. BTEST( wall_flags_0(k,j-1,i), 0 )  .OR.         &
     1581                           .NOT. BTEST( wall_flags_0(k,j+1,i), 0 ) )             &
     1582                         l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), 0.5_wp * dy )
     1583!
     1584!--                   x-distance
     1585                      IF ( .NOT. BTEST( wall_flags_0(k,j,i-1), 0 )  .OR.         &
     1586                           .NOT. BTEST( wall_flags_0(k,j,i+1), 0 ) )             &
     1587                         l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k), 0.5_wp * dx )
     1588!
     1589!--                   yz-distance (vertical edges, down)
     1590                      IF ( .NOT. BTEST( wall_flags_0(k-1,j-1,i), 0 )  .OR.       &
     1591                           .NOT. BTEST( wall_flags_0(k-1,j+1,i), 0 )  )          &
     1592                         l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),          &
     1593                                              SQRT( 0.25_wp * dy**2 +            &
     1594                                             ( zu(k) - zw(k-1) )**2 ) )
     1595!
     1596!--                   yz-distance (vertical edges, up)
     1597                      IF ( .NOT. BTEST( wall_flags_0(k+1,j-1,i), 0 )  .OR.       &
     1598                           .NOT. BTEST( wall_flags_0(k+1,j+1,i), 0 )  )          &
     1599                         l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),          &
     1600                                              SQRT( 0.25_wp * dy**2 +            &
     1601                                             ( zw(k) - zu(k) )**2 ) )
     1602!
     1603!--                   xz-distance (vertical edges, down)
     1604                      IF ( .NOT. BTEST( wall_flags_0(k-1,j,i-1), 0 )  .OR.       &
     1605                           .NOT. BTEST( wall_flags_0(k-1,j,i+1), 0 )  )          &
     1606                         l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),          &
     1607                                              SQRT( 0.25_wp * dx**2 +            &
     1608                                             ( zu(k) - zw(k-1) )**2 ) )
     1609!
     1610!--                   xz-distance (vertical edges, up)
     1611                      IF ( .NOT. BTEST( wall_flags_0(k+1,j,i-1), 0 )  .OR.       &
     1612                           .NOT. BTEST( wall_flags_0(k+1,j,i+1), 0 )  )          &
     1613                         l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),          &
     1614                                              SQRT( 0.25_wp * dx**2 +            &
     1615                                             ( zw(k) - zu(k) )**2 ) )
     1616!
     1617!--                   xy-distance (horizontal edges)
     1618                      IF ( .NOT. BTEST( wall_flags_0(k,j-1,i-1), 0 )  .OR.        &
     1619                           .NOT. BTEST( wall_flags_0(k,j+1,i-1), 0 )  .OR.        &
     1620                           .NOT. BTEST( wall_flags_0(k,j-1,i+1), 0 )  .OR.        &
     1621                           .NOT. BTEST( wall_flags_0(k,j+1,i+1), 0 ) )            &
     1622                         l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),           &
     1623                                              SQRT( 0.25_wp * ( dx**2 + dy**2 ) ) )
     1624!
     1625!--                   xyz distance (vertical and horizontal edges, down)
     1626                      IF ( .NOT. BTEST( wall_flags_0(k-1,j-1,i-1), 0 )  .OR.      &
     1627                           .NOT. BTEST( wall_flags_0(k-1,j+1,i-1), 0 )  .OR.      &
     1628                           .NOT. BTEST( wall_flags_0(k-1,j-1,i+1), 0 )  .OR.      &
     1629                           .NOT. BTEST( wall_flags_0(k-1,j+1,i+1), 0 ) )          &
     1630                         l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),           &
     1631                                              SQRT( 0.25_wp * ( dx**2 + dy**2 )   &
     1632                                                    +  ( zu(k) - zw(k-1) )**2  ) )
     1633!
     1634!--                   xyz distance (vertical and horizontal edges, up)
     1635                      IF ( .NOT. BTEST( wall_flags_0(k+1,j-1,i-1), 0 )  .OR.      &
     1636                           .NOT. BTEST( wall_flags_0(k+1,j+1,i-1), 0 )  .OR.      &
     1637                           .NOT. BTEST( wall_flags_0(k+1,j-1,i+1), 0 )  .OR.      &
     1638                           .NOT. BTEST( wall_flags_0(k+1,j+1,i+1), 0 ) )          &
     1639                         l_wall(k,j,i) = MIN( l_wall(k,j,i), l_grid(k),           &
     1640                                              SQRT( 0.25_wp * ( dx**2 + dy**2 )   &
     1641                                                    +  ( zw(k) - zu(k) )**2  ) )
     1642
     1643                   ENDIF
     1644!
     1645!--                Adjust mixing length by wall-adjustment factor and limit it by l_grid
     1646                   l_wall(k,j,i) = MIN( l_wall(k,j,i) * wall_adjustment_factor, l_grid(k) )
     1647
     1648                ENDDO  !k loop
     1649             ENDDO  !j loop
     1650          ENDDO  !i loop
     1651
     1652       ENDIF  !if wall_adjustment
    16411653
    16421654    ELSE
    16431655!
    1644 !-- Initialize the mixing length in case of a RANS simulation
     1656!--    Initialize the mixing length in case of a RANS simulation
    16451657       ALLOCATE( l_black(nzb:nzt+1) )
    16461658
     
    16621674!
    16631675!--    Get height level of highest topography within local subdomain
     1676       k_max_topo = 0
    16641677       DO  i = nxlg, nxrg
    16651678          DO  j = nysg, nyng
     
    40434056 SUBROUTINE diffusion_e( var, var_reference )
    40444057
    4045 #if defined( _OPENACC )
    40464058    USE arrays_3d,                                                             &
    4047         ONLY:  ddzu, dd2zu, ddzw, drho_air, rho_air_zw
     4059        ONLY:  dd2zu, ddzu, ddzw, drho_air, rho_air_zw
    40484060
    40494061    USE control_parameters,                                                    &
    4050         ONLY:  atmos_ocean_sign, use_single_reference_value,                   &
    4051                wall_adjustment, wall_adjustment_factor
    4052 #else
    4053     USE arrays_3d,                                                             &
    4054         ONLY:  ddzu, ddzw, drho_air, rho_air_zw
    4055 #endif
     4062        ONLY:  atmos_ocean_sign, use_single_reference_value
     4063
    40564064    USE grid_variables,                                                        &
    40574065        ONLY:  ddx2, ddy2
     
    40704078    INTEGER(iwp) ::  m              !< running index surface elements
    40714079
    4072     REAL(wp)     ::  flag           !< flag to mask topography
     4080    REAL(wp)     ::  duv2_dz2       !< squared vertical gradient of wind vector
     4081    REAL(wp)     ::  dvar_dz        !< vertical gradient of var
    40734082    REAL(wp)     ::  l              !< mixing length
    4074     REAL(wp)     ::  ll             !< adjusted l
    40754083    REAL(wp)     ::  var_reference  !< reference temperature
    4076 #if defined( _OPENACC )
    4077 !
    4078 !-- From mixing_length_les:
    4079     REAL(wp)     :: l_stable        !< mixing length according to stratification
    4080 !
    4081 !-- From mixing_length_rans:
    4082     REAL(wp)     :: duv2_dz2        !< squared vertical gradient of wind vector
    4083     REAL(wp)     :: l_diss          !< mixing length for dissipation
    4084     REAL(wp)     :: rif             !< Richardson flux number
    4085 !
    4086 !-- From both:
    4087     REAL(wp)     :: dvar_dz         !< vertical gradient of var
    4088 #endif
    4089 
    4090     REAL(wp)     ::  dissipation  !< TKE dissipation
    4091 
    4092     REAL(wp), DIMENSION(:,:,:), POINTER ::  var  !< temperature
    4093 
    4094 
    4095 !
    4096 !-- Calculate the tendency terms
     4084
     4085    REAL(wp), DIMENSION(nzb+1:nzt) ::  l_stable  !< mixing length according to stratification
     4086    REAL(wp), DIMENSION(nzb+1:nzt) ::  rif       !< Richardson flux number
     4087
     4088    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) ::  var  !< temperature
     4089
     4090!
     4091!-- Calculate the dissipation
     4092    IF ( les_dynamic .OR. les_mw )  THEN
     4093
     4094       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
     4095       !$ACC PRIVATE(l, l_stable, dvar_dz) &
     4096       !$ACC PRESENT(diss, e, var, wall_flags_0) &
     4097       !$ACC PRESENT(dd2zu, l_grid, l_wall)
     4098       DO  i = nxl, nxr
     4099          DO  j = nys, nyn
     4100             !$ACC LOOP PRIVATE(k)
     4101             DO  k = nzb+1, nzt
     4102
     4103                dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
     4104                IF ( dvar_dz > 0.0_wp ) THEN
     4105                   IF ( use_single_reference_value )  THEN
     4106                     l_stable(k) = 0.76_wp * SQRT( e(k,j,i) )                  &
     4107                                 / SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
     4108                   ELSE
     4109                     l_stable(k) = 0.76_wp * SQRT( e(k,j,i) )               &
     4110                                 / SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
     4111                   ENDIF
     4112                ELSE
     4113                   l_stable(k) = l_grid(k)
     4114                ENDIF
     4115
     4116             ENDDO
     4117
     4118             !$ACC LOOP PRIVATE(k)
     4119             !DIR$ IVDEP
     4120             DO  k = nzb+1, nzt
     4121
     4122                l  = MIN( l_wall(k,j,i), l_stable(k) )
     4123
     4124                diss(k,j,i) = ( 0.19_wp + 0.74_wp * l / l_wall(k,j,i) )                &
     4125                              * e(k,j,i) * SQRT( e(k,j,i) ) / l                        &
     4126                              * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
     4127
     4128             ENDDO
     4129
     4130          ENDDO
     4131       ENDDO
     4132
     4133    ELSEIF ( rans_tke_l )  THEN
     4134
     4135      !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
     4136      !$ACC PRIVATE(l_stable, duv2_dz2, rif, dvar_dz) &
     4137      !$ACC PRESENT(diss, e, u, v, var, wall_flags_0) &
     4138      !$ACC PRESENT(dd2zu, l_black, l_wall)
     4139       DO  i = nxl, nxr
     4140          DO  j = nys, nyn
     4141!
     4142!--          Calculate Richardson-flux number
     4143             IF ( use_single_reference_value )  THEN
     4144                !$ACC LOOP PRIVATE(k)
     4145                DO  k = nzb+1, nzt
     4146
     4147                   dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
     4148
     4149                   duv2_dz2 =   ( ( u(k+1,j,i) - u(k-1,j,i) ) * dd2zu(k) )**2  &
     4150                              + ( ( v(k+1,j,i) - v(k-1,j,i) ) * dd2zu(k) )**2  &
     4151                              + 1E-30_wp
     4152
     4153                   rif(k) = MIN( MAX( g / var_reference * dvar_dz / duv2_dz2, -5.0_wp ),  1.0_wp )
     4154                ENDDO
     4155             ELSE
     4156                !$ACC LOOP PRIVATE(k)
     4157                DO  k = nzb+1, nzt
     4158
     4159                   dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
     4160
     4161                   duv2_dz2 =   ( ( u(k+1,j,i) - u(k-1,j,i) ) * dd2zu(k) )**2  &
     4162                              + ( ( v(k+1,j,i) - v(k-1,j,i) ) * dd2zu(k) )**2  &
     4163                              + 1E-30_wp
     4164
     4165                   rif(k) = MIN( MAX( g / var(k,j,i) * dvar_dz / duv2_dz2, -5.0_wp ),  1.0_wp )
     4166                ENDDO
     4167             ENDIF
     4168!
     4169!--          Calculate diabatic mixing length using Dyer-profile functions
     4170             !$ACC LOOP PRIVATE(k)
     4171             DO  k = nzb+1, nzt
     4172                IF ( rif(k) >= 0.0_wp )  THEN
     4173                   l_stable(k) = MIN( l_black(k) / ( 1.0_wp + 5.0_wp * rif(k) ), l_wall(k,j,i) )
     4174                ELSE
     4175                   l_stable(k) = l_wall(k,j,i) * SQRT( 1.0_wp - 16.0_wp * rif(k) )
     4176                ENDIF
     4177             ENDDO
     4178
     4179             !$ACC LOOP PRIVATE(k)
     4180             !DIR$ IVDEP
     4181             DO  k = nzb+1, nzt
     4182
     4183                diss(k,j,i) = c_0**3 * e(k,j,i) * SQRT( e(k,j,i) ) / l_stable(k)     &
     4184                            * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
     4185
     4186             ENDDO
     4187
     4188          ENDDO
     4189       ENDDO
     4190
     4191!-- Note, in case of rans_tke_e, the dissipation is already calculated
     4192!-- in prognostic_equations
     4193    ENDIF
     4194
    40974195    !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(i, j, k) &
    4098     !$ACC PRIVATE(flag, l, ll, dissipation) &
    4099     !$ACC PRIVATE(l_stable, duv2_dz2, l_diss, rif, dvar_dz) &
    4100     !$ACC PRESENT(e, u, v, km, var, wall_flags_0) &
    4101     !$ACC PRESENT(ddzu, dd2zu, ddzw, rho_air_zw, drho_air) &
    4102     !$ACC PRESENT(l_black, l_grid, l_wall) &
    4103     !$ACC PRESENT(diss, tend)
     4196    !$ACC PRESENT(diss, e, km, tend, wall_flags_0) &
     4197    !$ACC PRESENT(ddzu, ddzw, rho_air_zw, drho_air)
    41044198    DO  i = nxl, nxr
    41054199       DO  j = nys, nyn
    41064200          DO  k = nzb+1, nzt
    4107 !
    4108 !--          Predetermine flag to mask topography
    4109              flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
    4110 
    4111 !
    4112 !--          Calculate dissipation
    4113              IF ( les_dynamic .OR. les_mw )  THEN
    4114 
    4115 #if defined( _OPENACC )
    4116 !
    4117 !--             Cannot call subroutine mixing_length_les because after adding all required
    4118 !--             OpenACC directives the execution crashed reliably when inside the called
    4119 !--             subroutine. I'm not sure why that is...
    4120                 dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
    4121                 IF ( dvar_dz > 0.0_wp ) THEN
    4122                    IF ( use_single_reference_value )  THEN
    4123                       l_stable = 0.76_wp * SQRT( e(k,j,i) )                                &
    4124                                          / SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
    4125                    ELSE
    4126                       l_stable = 0.76_wp * SQRT( e(k,j,i) )                                &
    4127                                          / SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
    4128                    ENDIF
    4129                 ELSE
    4130                    l_stable = l_grid(k)
    4131                 ENDIF
    4132 !
    4133 !--             Adjustment of the mixing length
    4134                 IF ( wall_adjustment )  THEN
    4135                    l  = MIN( wall_adjustment_factor * l_wall(k,j,i), l_grid(k), l_stable )
    4136                    ll = MIN( wall_adjustment_factor * l_wall(k,j,i), l_grid(k) )
    4137                 ELSE
    4138                    l  = MIN( l_grid(k), l_stable )
    4139                    ll = l_grid(k)
    4140                 ENDIF
    4141 #else
    4142                 CALL mixing_length_les( i, j, k, l, ll, var, var_reference )
    4143 #endif
    4144 
    4145                 dissipation = ( 0.19_wp + 0.74_wp * l / ll )                   &
    4146                               * e(k,j,i) * SQRT( e(k,j,i) ) / l
    4147 
    4148              ELSEIF ( rans_tke_l )  THEN
    4149 
    4150 #if defined( _OPENACC )
    4151 !
    4152 !--             Same reason as above...
    4153                 dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
    4154 
    4155                 duv2_dz2 =   ( ( u(k+1,j,i) - u(k-1,j,i) ) * dd2zu(k) )**2                 &
    4156                            + ( ( v(k+1,j,i) - v(k-1,j,i) ) * dd2zu(k) )**2                 &
    4157                            + 1E-30_wp
    4158 
    4159                 IF ( use_single_reference_value )  THEN
    4160                    rif = g / var_reference * dvar_dz / duv2_dz2
    4161                 ELSE
    4162                    rif = g / var(k,j,i) * dvar_dz / duv2_dz2
    4163                 ENDIF
    4164 
    4165                 rif = MAX( rif, -5.0_wp )
    4166                 rif = MIN( rif,  1.0_wp )
    4167 
    4168 !
    4169 !--             Calculate diabatic mixing length using Dyer-profile functions
    4170                 IF ( rif >= 0.0_wp )  THEN
    4171                    l      = MIN( l_black(k) / ( 1.0_wp + 5.0_wp * rif ), l_wall(k,j,i) )
    4172                    l_diss = l
    4173                 ELSE
    4174 !
    4175 !--                In case of unstable stratification, use mixing length of neutral case
    4176 !--                for l, but consider profile functions for l_diss
    4177                    l      = l_wall(k,j,i)
    4178                    l_diss = l * SQRT( 1.0_wp - 16.0_wp * rif )
    4179                 ENDIF
    4180 #else
    4181                 CALL mixing_length_rans( i, j, k, l, ll, var, var_reference )
    4182 #endif
    4183 
    4184                 dissipation = c_0**3 * e(k,j,i) * SQRT( e(k,j,i) ) / ll
    4185 
    4186                 diss(k,j,i) = dissipation * flag
    4187 
    4188              ELSEIF ( rans_tke_e )  THEN
    4189 
    4190                 dissipation = diss(k,j,i)
    4191 
    4192              ENDIF
    4193 
    4194              tend(k,j,i) = tend(k,j,i) + (                                     &
    4195                                            (                                   &
    4196                        ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )     &
    4197                      - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )     &
    4198                                            ) * ddx2  * flag                    &
    4199                                          + (                                   &
    4200                        ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )     &
    4201                      - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )     &
    4202                                            ) * ddy2  * flag                    &
    4203                                          + (                                   &
    4204             ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1)    &
    4205                                                           * rho_air_zw(k)      &
    4206           - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)      &
    4207                                                           * rho_air_zw(k-1)    &
    4208                                            ) * ddzw(k) * drho_air(k)           &
    4209                                          ) * flag * dsig_e                     &
    4210                           - dissipation * flag
    4211 !
    4212 !--          Store dissipation if needed for calculating the sgs particle
    4213 !--          velocities
    4214              IF ( .NOT. rans_tke_e .AND. ( use_sgs_for_particles  .OR.         &
    4215                   wang_kernel  .OR.  collision_turbulence  ) )  THEN
    4216                 diss(k,j,i) = dissipation * flag
    4217              ENDIF
     4201
     4202             tend(k,j,i) = tend(k,j,i) + (                                            &
     4203                                           (                                          &
     4204                       ( km(k,j,i)+km(k,j,i+1) ) * ( e(k,j,i+1)-e(k,j,i) )            &
     4205                     - ( km(k,j,i)+km(k,j,i-1) ) * ( e(k,j,i)-e(k,j,i-1) )            &
     4206                                           ) * ddx2                                   &
     4207                                         + (                                          &
     4208                       ( km(k,j,i)+km(k,j+1,i) ) * ( e(k,j+1,i)-e(k,j,i) )            &
     4209                     - ( km(k,j,i)+km(k,j-1,i) ) * ( e(k,j,i)-e(k,j-1,i) )            &
     4210                                           ) * ddy2                                   &
     4211                                         + (                                          &
     4212            ( km(k,j,i)+km(k+1,j,i) ) * ( e(k+1,j,i)-e(k,j,i) ) * ddzu(k+1)           &
     4213                                                          * rho_air_zw(k)             &
     4214          - ( km(k,j,i)+km(k-1,j,i) ) * ( e(k,j,i)-e(k-1,j,i) ) * ddzu(k)             &
     4215                                                          * rho_air_zw(k-1)           &
     4216                                           ) * ddzw(k) * drho_air(k)                  &
     4217                                         ) * dsig_e                                   &
     4218                                           * MERGE( 1.0_wp, 0.0_wp,                   &
     4219                                                    BTEST( wall_flags_0(k,j,i), 0 ) ) &
     4220                          - diss(k,j,i)
    42184221
    42194222          ENDDO
     
    42584261
    42594262    USE arrays_3d,                                                             &
    4260         ONLY:  ddzu, ddzw, drho_air, rho_air_zw
     4263        ONLY:  dd2zu, ddzu, ddzw, drho_air, rho_air_zw
     4264
     4265    USE control_parameters,                                                    &
     4266        ONLY:  atmos_ocean_sign, use_single_reference_value
    42614267
    42624268    USE grid_variables,                                                        &
     
    42784284    INTEGER(iwp) ::  surf_s         !< Start index of surface elements at (j,i)-gridpoint
    42794285
    4280     REAL(wp)     ::  flag           !< flag to mask topography
     4286    REAL(wp)     ::  duv2_dz2       !< squared vertical gradient of wind vector
     4287    REAL(wp)     ::  dvar_dz        !< vertical gradient of var
    42814288    REAL(wp)     ::  l              !< mixing length
    4282     REAL(wp)     ::  ll             !< adjusted l
    42834289    REAL(wp)     ::  var_reference  !< reference temperature
    42844290
    4285     REAL(wp), DIMENSION(nzb+1:nzt) ::  dissipation  !< dissipation of TKE
    4286 
    4287     REAL(wp), DIMENSION(:,:,:), POINTER ::  var     !< temperature
    4288 
    4289 !
    4290 !-- Calculate the mixing length (for dissipation)
     4291    REAL(wp), DIMENSION(nzb+1:nzt) ::  l_stable  !< mixing length according to stratification
     4292    REAL(wp), DIMENSION(nzb+1:nzt) ::  rif       !< Richardson flux number
     4293
     4294    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) ::  var  !< temperature
     4295
     4296!
     4297!-- Calculate the mixing length and dissipation...
     4298!-- ...in case of LES
     4299    IF ( les_dynamic .OR. les_mw )  THEN
     4300
     4301       DO  k = nzb+1, nzt
     4302          dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
     4303          IF ( dvar_dz > 0.0_wp ) THEN
     4304             IF ( use_single_reference_value )  THEN
     4305                l_stable(k) = 0.76_wp * SQRT( e(k,j,i) )                  &
     4306                            / SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
     4307             ELSE
     4308                l_stable(k) = 0.76_wp * SQRT( e(k,j,i) )               &
     4309                            / SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
     4310             ENDIF
     4311          ELSE
     4312             l_stable(k) = l_grid(k)
     4313          ENDIF
     4314       ENDDO
     4315
     4316       !DIR$ IVDEP
     4317       DO  k = nzb+1, nzt
     4318          l  = MIN( l_wall(k,j,i), l_stable(k) )
     4319
     4320          diss(k,j,i) = ( 0.19_wp + 0.74_wp * l / l_wall(k,j,i) )              &
     4321                        * e(k,j,i) * SQRT( e(k,j,i) ) / l                      &
     4322                        * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
     4323       ENDDO
     4324
     4325!
     4326!-- ...in case of RANS
     4327    ELSEIF ( rans_tke_l )  THEN
     4328
     4329!
     4330!--    Calculate Richardson-flux number
     4331       IF ( use_single_reference_value )  THEN
     4332          DO  k = nzb+1, nzt
     4333             dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
     4334
     4335             duv2_dz2 =   ( ( u(k+1,j,i) - u(k-1,j,i) ) * dd2zu(k) )**2           &
     4336                        + ( ( v(k+1,j,i) - v(k-1,j,i) ) * dd2zu(k) )**2           &
     4337                        + 1E-30_wp
     4338
     4339             rif(k) = MIN( MAX( g / var_reference * dvar_dz / duv2_dz2, -5.0_wp ),  1.0_wp )
     4340          ENDDO
     4341       ELSE
     4342          DO  k = nzb+1, nzt
     4343             dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
     4344
     4345             duv2_dz2 =   ( ( u(k+1,j,i) - u(k-1,j,i) ) * dd2zu(k) )**2           &
     4346                        + ( ( v(k+1,j,i) - v(k-1,j,i) ) * dd2zu(k) )**2           &
     4347                        + 1E-30_wp
     4348
     4349             rif(k) = MIN( MAX( g / var(k,j,i) * dvar_dz / duv2_dz2, -5.0_wp ), 1.0_wp )
     4350          ENDDO
     4351       ENDIF
     4352!
     4353!--    Calculate diabatic mixing length using Dyer-profile functions
     4354       DO  k = nzb+1, nzt
     4355          IF ( rif(k) >= 0.0_wp )  THEN
     4356             l_stable(k) = MIN( l_black(k) / ( 1.0_wp + 5.0_wp * rif(k) ), l_wall(k,j,i) )
     4357          ELSE
     4358             l_stable(k) = l_wall(k,j,i) * SQRT( 1.0_wp - 16.0_wp * rif(k) )
     4359          ENDIF
     4360
     4361       ENDDO
     4362
     4363       !DIR$ IVDEP
     4364       DO  k = nzb+1, nzt
     4365          diss(k,j,i) = c_0**3 * e(k,j,i) * SQRT( e(k,j,i) ) / l_stable(k)     &
     4366                      * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
     4367       ENDDO
     4368
     4369!-- Note, in case of rans_tke_e, the dissipation is already calculated
     4370!-- in prognostic_equations
     4371    ENDIF
     4372
     4373!
     4374!-- Calculate the tendency term
     4375    !DIR$ IVDEP
    42914376    DO  k = nzb+1, nzt
    4292 !
    4293 !--    Predetermine flag to mask topography
    4294        flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
    4295 
    4296 !
    4297 !--    Calculate dissipation...
    4298 !--    ...in case of LES
    4299        IF ( les_dynamic .OR. les_mw )  THEN
    4300 
    4301           CALL mixing_length_les( i, j, k, l, ll, var, var_reference )
    4302 
    4303           dissipation(k) = ( 0.19_wp + 0.74_wp * l / ll )                      &
    4304                            * e(k,j,i) * SQRT( e(k,j,i) ) / l
    4305 
    4306 !
    4307 !--    ...in case of RANS
    4308        ELSEIF ( rans_tke_l )  THEN
    4309 
    4310           CALL mixing_length_rans( i, j, k, l, ll, var, var_reference  )
    4311 
    4312           dissipation(k) = c_0**3 * e(k,j,i) * SQRT( e(k,j,i) ) / ll
    4313 
    4314           diss(k,j,i) = dissipation(k) * flag
    4315 
    4316        ELSEIF ( rans_tke_e )  THEN
    4317 
    4318           dissipation(k) = diss(k,j,i)
    4319 
    4320        ENDIF
    4321 
    4322 !
    4323 !--    Calculate the tendency term
     4377
    43244378       tend(k,j,i) = tend(k,j,i) + (                                           &
    43254379                                      (                                        &
     
    43374391                                                         * rho_air_zw(k-1)     &
    43384392                                      ) * ddzw(k) * drho_air(k)                &
    4339                                    ) * flag * dsig_e                           &
    4340                                  - dissipation(k) * flag
     4393                                   ) * dsig_e                                  &
     4394                                     * MERGE( 1.0_wp, 0.0_wp,                  &
     4395                                              BTEST( wall_flags_0(k,j,i), 0 ) )&
     4396                                 - diss(k,j,i)
    43414397
    43424398    ENDDO
    43434399
    43444400!
    4345 !-- Store dissipation if needed for calculating the sgs particle velocities
     4401!-- Set boundary conditions of dissipation if needed for calculating the sgs
     4402!-- particle velocities.
     4403!-- Neumann boundary condition for dissipation diss(nzb,:,:) = diss(nzb+1,:,:)
     4404!-- For each surface type determine start and end index (in case of elevated
     4405!-- topography several up/downward facing surfaces may exist.
     4406!-- Note, bc cannot be set in tcm_boundary conditions as the dissipation
     4407!-- in LES mode is only a diagnostic quantity.
    43464408    IF ( .NOT. rans_tke_e .AND.  ( use_sgs_for_particles  .OR.  wang_kernel    &
    43474409          .OR.  collision_turbulence ) )  THEN
    4348        DO  k = nzb+1, nzt
    4349           diss(k,j,i) = dissipation(k) * MERGE( 1.0_wp, 0.0_wp,                &
    4350                                                BTEST( wall_flags_0(k,j,i), 0 ) )
    4351        ENDDO
    4352 !
    4353 !--    Neumann boundary condition for dissipation diss(nzb,:,:) = diss(nzb+1,:,:)
    4354 !--    For each surface type determine start and end index (in case of elevated
    4355 !--    topography several up/downward facing surfaces may exist.
    4356 !--    Note, bc cannot be set in tcm_boundary conditions as the dissipation
    4357 !--    in LES mode is only a diagnostic quantity.
    43584410       surf_s = bc_h(0)%start_index(j,i)
    43594411       surf_e = bc_h(0)%end_index(j,i)
     
    44904542! Description:
    44914543! ------------
    4492 !> Calculate mixing length for LES mode.
    4493 !------------------------------------------------------------------------------!
    4494  SUBROUTINE mixing_length_les( i, j, k, l, ll, var, var_reference )
    4495 
    4496     USE arrays_3d,                                                             &
    4497         ONLY:  dd2zu
    4498 
    4499     USE control_parameters,                                                    &
    4500         ONLY:  atmos_ocean_sign, use_single_reference_value,                   &
    4501                wall_adjustment, wall_adjustment_factor
    4502 
    4503     IMPLICIT NONE
    4504 
    4505     INTEGER(iwp) :: i   !< loop index
    4506     INTEGER(iwp) :: j   !< loop index
    4507     INTEGER(iwp) :: k   !< loop index
    4508 
    4509     REAL(wp)     :: dvar_dz         !< vertical gradient of var
    4510     REAL(wp)     :: l               !< mixing length
    4511     REAL(wp)     :: l_stable        !< mixing length according to stratification
    4512     REAL(wp)     :: ll              !< adjusted l_grid
    4513     REAL(wp)     :: var_reference   !< var at reference height
    4514 
    4515     REAL(wp), DIMENSION(:,:,:), POINTER ::  var  !< temperature
    4516 
    4517 
    4518     dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
    4519     IF ( dvar_dz > 0.0_wp ) THEN
    4520        IF ( use_single_reference_value )  THEN
    4521           l_stable = 0.76_wp * SQRT( e(k,j,i) )                                &
    4522                              / SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
    4523        ELSE
    4524           l_stable = 0.76_wp * SQRT( e(k,j,i) )                                &
    4525                              / SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
    4526        ENDIF
    4527     ELSE
    4528        l_stable = l_grid(k)
    4529     ENDIF
    4530 !
    4531 !-- Adjustment of the mixing length
    4532     IF ( wall_adjustment )  THEN
    4533        l  = MIN( wall_adjustment_factor * l_wall(k,j,i), l_grid(k), l_stable )
    4534        ll = MIN( wall_adjustment_factor * l_wall(k,j,i), l_grid(k) )
    4535     ELSE
    4536        l  = MIN( l_grid(k), l_stable )
    4537        ll = l_grid(k)
    4538     ENDIF
    4539 
    4540  END SUBROUTINE mixing_length_les
    4541 
    4542 
    4543 !------------------------------------------------------------------------------!
    4544 ! Description:
    4545 ! ------------
    4546 !> Calculate mixing length for RANS mode.
    4547 !------------------------------------------------------------------------------!
    4548  SUBROUTINE mixing_length_rans( i, j, k, l, l_diss, var, var_reference  )
    4549 
    4550     USE arrays_3d,                                                             &
    4551         ONLY:  dd2zu
    4552 
    4553     USE control_parameters,                                                    &
    4554         ONLY:  atmos_ocean_sign, use_single_reference_value
    4555 
    4556     IMPLICIT NONE
    4557 
    4558     INTEGER(iwp) :: i   !< loop index
    4559     INTEGER(iwp) :: j   !< loop index
    4560     INTEGER(iwp) :: k   !< loop index
    4561 
    4562     REAL(wp)     :: duv2_dz2        !< squared vertical gradient of wind vector
    4563     REAL(wp)     :: dvar_dz         !< vertical gradient of var
    4564     REAL(wp)     :: l               !< mixing length
    4565     REAL(wp)     :: l_diss          !< mixing length for dissipation
    4566     REAL(wp)     :: rif             !< Richardson flux number
    4567     REAL(wp)     :: var_reference   !< var at reference height
    4568 
    4569     REAL(wp), DIMENSION(:,:,:), POINTER ::  var  !< temperature
    4570 
    4571 
    4572     dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
    4573 
    4574     duv2_dz2 =   ( ( u(k+1,j,i) - u(k-1,j,i) ) * dd2zu(k) )**2                 &
    4575                + ( ( v(k+1,j,i) - v(k-1,j,i) ) * dd2zu(k) )**2                 &
    4576                + 1E-30_wp
    4577 
    4578     IF ( use_single_reference_value )  THEN
    4579        rif = g / var_reference * dvar_dz / duv2_dz2
    4580     ELSE
    4581        rif = g / var(k,j,i) * dvar_dz / duv2_dz2
    4582     ENDIF
    4583 
    4584     rif = MAX( rif, -5.0_wp )
    4585     rif = MIN( rif,  1.0_wp )
    4586 
    4587 !
    4588 !-- Calculate diabatic mixing length using Dyer-profile functions
    4589     IF ( rif >= 0.0_wp )  THEN
    4590        l      = MIN( l_black(k) / ( 1.0_wp + 5.0_wp * rif ), l_wall(k,j,i) )
    4591        l_diss = l
    4592     ELSE
    4593 !
    4594 !--    In case of unstable stratification, use mixing length of neutral case
    4595 !--    for l, but consider profile functions for l_diss
    4596        l      = l_wall(k,j,i)
    4597        l_diss = l * SQRT( 1.0_wp - 16.0_wp * rif )
    4598     ENDIF
    4599 
    4600  END SUBROUTINE mixing_length_rans
    4601 
    4602 
    4603 !------------------------------------------------------------------------------!
    4604 ! Description:
    4605 ! ------------
    46064544!> Computation of the turbulent diffusion coefficients for momentum and heat.
    46074545!> @bug unstable stratification is not properly considered for kh in rans mode.
     
    46244562    REAL(wp) ::  var_reference  !< reference temperature
    46254563
    4626     REAL(wp), DIMENSION(:,:,:), POINTER ::  var  !< temperature
     4564    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) ::  var  !< temperature
    46274565
    46284566
     
    47974735 SUBROUTINE tcm_diffusivities_default( var, var_reference )
    47984736
    4799 #if defined( _OPENACC )
    48004737    USE arrays_3d,                                                             &
    48014738        ONLY:  dd2zu
    48024739
    48034740    USE control_parameters,                                                    &
    4804         ONLY:  atmos_ocean_sign, use_single_reference_value,                   &
    4805                wall_adjustment, wall_adjustment_factor
    4806 #endif
     4741        ONLY:  atmos_ocean_sign, use_single_reference_value
    48074742
    48084743    USE statistics,                                                            &
     
    48184753    INTEGER(iwp) ::  tn                  !< thread number
    48194754
    4820     REAL(wp)     ::  flag                !< topography flag
    4821     REAL(wp)     ::  l                   !< mixing length
    4822     REAL(wp)     ::  ll                  !< adjusted mixing length
     4755    REAL(wp)     ::  duv2_dz2            !< squared vertical gradient of wind vector
     4756    REAL(wp)     ::  dvar_dz             !< vertical gradient of var
     4757    REAL(wp)     ::  l                   !< mixing length (single height)
    48234758    REAL(wp)     ::  var_reference       !< reference temperature
    4824 #if defined( _OPENACC )
    4825 !
    4826 !-- From mixing_length_les:
    4827     REAL(wp)     :: l_stable        !< mixing length according to stratification
    4828     REAL(wp)     :: dvar_dz         !< vertical gradient of var
    4829 #endif
    4830 
    4831     REAL(wp), DIMENSION(:,:,:), POINTER ::  var  !< temperature
    4832 
     4759
     4760    !DIR$ ATTRIBUTES ALIGN:64:: l_v, l_stable, rif
     4761    REAL(wp), DIMENSION(nzb+1:nzt) ::  l_v       !< mixing length (all heights)
     4762    REAL(wp), DIMENSION(nzb+1:nzt) ::  l_stable  !< mixing length according to stratification
     4763    REAL(wp), DIMENSION(nzb+1:nzt) ::  rif       !< Richardson flux number
     4764
     4765    REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) ::  var  !< temperature
    48334766
    48344767!
     
    48444777!
    48454778!-- Compute the turbulent diffusion coefficient for momentum
    4846     !$OMP PARALLEL PRIVATE (i,j,k,l,ll,sr,tn,flag)
     4779    !$OMP PARALLEL PRIVATE (i,j,k,l,sr,tn)
    48474780!$  tn = omp_get_thread_num()
    48484781
    48494782    IF ( les_dynamic .OR. les_mw )  THEN
    48504783       !$OMP DO
    4851        !$ACC PARALLEL LOOP COLLAPSE(3) PRIVATE(sr) &
    4852        !$ACC PRIVATE(flag, dvar_dz, l_stable, l, ll) &
     4784       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
     4785       !$ACC PRIVATE(dvar_dz, l, l_stable, l_v) &
    48534786       !$ACC PRESENT(wall_flags_0, var, dd2zu, e, l_wall, l_grid, rmask) &
    48544787       !$ACC PRESENT(kh, km, sums_l_l)
    48554788       DO  i = nxlg, nxrg
    48564789          DO  j = nysg, nyng
     4790             !$ACC LOOP PRIVATE(k)
    48574791             DO  k = nzb+1, nzt
    4858 
    4859                 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
    4860 
    4861 !
    4862 !--             Determine the mixing length for LES closure
    4863 #if defined( _OPENACC )
    4864 !
    4865 !--             Cannot call subroutine mixing_length_les because after adding all required
    4866 !--             OpenACC directives the execution crashed reliably when inside the called
    4867 !--             subroutine. I'm not sure why that is...
     4792!
     4793!--             Determine the mixing length
     4794!--             @note The following code cannot be transferred to a subroutine
     4795!--             due to errors when using OpenACC directives. The execution
     4796!--             crashes reliably if a subroutine is called at this point (the
     4797!--             reasong for this behaviour is unknown, however).
    48684798                dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
    48694799                IF ( dvar_dz > 0.0_wp ) THEN
    48704800                   IF ( use_single_reference_value )  THEN
    4871                       l_stable = 0.76_wp * SQRT( e(k,j,i) )                                &
    4872                                          / SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
     4801                      l_stable(k) = 0.76_wp * SQRT( e(k,j,i) )                  &
     4802                                  / SQRT( g / var_reference * dvar_dz ) + 1E-5_wp
    48734803                   ELSE
    4874                       l_stable = 0.76_wp * SQRT( e(k,j,i) )                                &
    4875                                          / SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
     4804                      l_stable(k) = 0.76_wp * SQRT( e(k,j,i) )               &
     4805                                  / SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp
    48764806                   ENDIF
    48774807                ELSE
    4878                    l_stable = l_grid(k)
     4808                   l_stable(k) = l_grid(k)
    48794809                ENDIF
    4880 !
    4881 !--             Adjustment of the mixing length
    4882                 IF ( wall_adjustment )  THEN
    4883                    l  = MIN( wall_adjustment_factor * l_wall(k,j,i), l_grid(k), l_stable )
    4884                    ll = MIN( wall_adjustment_factor * l_wall(k,j,i), l_grid(k) )
    4885                 ELSE
    4886                    l  = MIN( l_grid(k), l_stable )
    4887                    ll = l_grid(k)
    4888                 ENDIF
    4889 #else
    4890                 CALL mixing_length_les( i, j, k, l, ll, var, var_reference )
    4891 #endif
    4892 
     4810
     4811             ENDDO
     4812
     4813             !$ACC LOOP PRIVATE(k)
     4814             !DIR$ IVDEP
     4815             DO  k = nzb+1, nzt
     4816
     4817                l_v(k) = MIN( l_wall(k,j,i), l_stable(k) )                      &
     4818                       * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
     4819                l = l_v(k)
    48934820!
    48944821!--             Compute diffusion coefficients for momentum and heat
    4895                 km(k,j,i) = c_0 * l * SQRT( e(k,j,i) ) * flag
    4896                 kh(k,j,i) = ( 1.0_wp + 2.0_wp * l / ll ) * km(k,j,i) * flag
    4897 !
    4898 !--             Summation for averaged profile (cf. flow_statistics)
    4899                 DO  sr = 0, statistic_regions
    4900                    sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l * rmask(j,i,sr)   &
    4901                                                              * flag
     4822                km(k,j,i) = c_0 * l * SQRT( e(k,j,i) )
     4823                kh(k,j,i) = ( 1.0_wp + 2.0_wp * l / l_wall(k,j,i) ) * km(k,j,i)
     4824
     4825             ENDDO
     4826!
     4827!--          Summation for averaged profile (cf. flow_statistics)
     4828             !$ACC LOOP PRIVATE(sr, k)
     4829             DO  sr = 0, statistic_regions
     4830                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)
    49024832                ENDDO
    4903 
    4904              ENDDO
     4833             ENDDO
     4834
    49054835          ENDDO
    49064836       ENDDO
     
    49094839
    49104840       !$OMP DO
     4841       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
     4842       !$ACC PRIVATE(dvar_dz, duv2_dz2, l_stable, l_v, rif) &
     4843       !$ACC PRESENT(wall_flags_0, var, dd2zu, e, u, v, l_wall, l_black, rmask) &
     4844       !$ACC PRESENT(kh, km, sums_l_l)
    49114845       DO  i = nxlg, nxrg
    49124846          DO  j = nysg, nyng
     4847!
     4848!--          Calculate Richardson-flux number
     4849             IF ( use_single_reference_value )  THEN
     4850                !$ACC LOOP PRIVATE(k)
     4851                !DIR$ IVDEP
     4852                DO  k = nzb+1, nzt
     4853                   dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
     4854
     4855                   duv2_dz2 =   ( ( u(k+1,j,i) - u(k-1,j,i) ) * dd2zu(k) )**2  &
     4856                              + ( ( v(k+1,j,i) - v(k-1,j,i) ) * dd2zu(k) )**2  &
     4857                              + 1E-30_wp
     4858
     4859                   rif(k) = MIN( MAX( g / var_reference * dvar_dz / duv2_dz2, -5.0_wp ),  1.0_wp )
     4860                ENDDO
     4861             ELSE
     4862                !$ACC LOOP PRIVATE(k)
     4863                !DIR$ IVDEP
     4864                DO  k = nzb+1, nzt
     4865                   dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k)
     4866
     4867                   duv2_dz2 =   ( ( u(k+1,j,i) - u(k-1,j,i) ) * dd2zu(k) )**2  &
     4868                              + ( ( v(k+1,j,i) - v(k-1,j,i) ) * dd2zu(k) )**2  &
     4869                              + 1E-30_wp
     4870
     4871                   rif(k) = MIN( MAX( g / var(k,j,i) * dvar_dz / duv2_dz2, -5.0_wp ),  1.0_wp )
     4872                ENDDO
     4873             ENDIF
     4874!
     4875!--          Calculate diabatic mixing length using Dyer-profile functions
     4876!--          In case of unstable stratification, use mixing length of neutral case
     4877             !$ACC LOOP PRIVATE(k)
    49134878             DO  k = nzb+1, nzt
    4914 
    4915                 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
    4916 !
    4917 !--             Mixing length for RANS mode with TKE-l closure
    4918                 CALL mixing_length_rans( i, j, k, l, ll, var, var_reference )
    4919 !
    4920 !--             Compute diffusion coefficients for momentum and heat
    4921                 km(k,j,i) = c_0 * l * SQRT( e(k,j,i) ) * flag
    4922                 kh(k,j,i) = km(k,j,i) / prandtl_number * flag
    4923 !
    4924 !--             Summation for averaged profile (cf. flow_statistics)
    4925                 DO  sr = 0, statistic_regions
    4926                    sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l * rmask(j,i,sr)   &
    4927                                                              * flag
     4879                IF ( rif(k) >= 0.0_wp )  THEN
     4880                   l_stable(k)  = MIN( l_black(k) / ( 1.0_wp + 5.0_wp * rif(k) ), l_wall(k,j,i) )
     4881                ELSE
     4882                   l_stable(k)  = l_wall(k,j,i)
     4883                ENDIF
     4884
     4885             ENDDO
     4886!
     4887!--          Compute diffusion coefficients for momentum and heat
     4888             !$ACC LOOP PRIVATE(k)
     4889             !DIR$ IVDEP
     4890             DO  k = nzb+1, nzt
     4891                l_v(k)    = l_stable(k) * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
     4892                km(k,j,i) = c_0 * l_v(k) * SQRT( e(k,j,i) )
     4893                kh(k,j,i) = km(k,j,i) / prandtl_number
     4894             ENDDO
     4895!
     4896!--          Summation for averaged profile (cf. flow_statistics)
     4897             !$ACC LOOP PRIVATE(sr, k)
     4898             DO  sr = 0, statistic_regions
     4899                DO  k = nzb+1, nzt
     4900                   sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l_v(k) * rmask(j,i,sr)
    49284901                ENDDO
    4929 
    4930              ENDDO
     4902             ENDDO
     4903
    49314904          ENDDO
    49324905       ENDDO
     
    49354908
    49364909       !$OMP DO
     4910       !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) &
     4911       !$ACC PRIVATE(l_v) &
     4912       !$ACC PRESENT(wall_flags_0, e, diss, rmask) &
     4913       !$ACC PRESENT(kh, km, sums_l_l)
    49374914       DO  i = nxlg, nxrg
    49384915          DO  j = nysg, nyng
     4916!
     4917!--          Compute diffusion coefficients for momentum and heat
     4918             !$ACC LOOP PRIVATE(k)
     4919             !DIR$ IVDEP
    49394920             DO  k = nzb+1, nzt
    49404921
    4941                 flag = MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
    4942 !
    4943 !--             Compute diffusion coefficients for momentum and heat
    4944                 km(k,j,i) = c_0**4 * e(k,j,i)**2 / ( diss(k,j,i) + 1.0E-30_wp ) * flag
    4945                 kh(k,j,i) = km(k,j,i) / prandtl_number * flag
    4946 !
    4947 !--             Summation for averaged profile of mixing length (cf. flow_statistics)
    4948                 DO  sr = 0, statistic_regions
    4949                    sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) +                     &
    4950                       c_0**3 * e(k,j,i) * SQRT(e(k,j,i)) /                     &
    4951                       ( diss(k,j,i) + 1.0E-30_wp ) * rmask(j,i,sr) * flag
     4922                l_v(k) = c_0**3 * e(k,j,i) * SQRT(e(k,j,i)) / ( diss(k,j,i) + 1.0E-30_wp ) &
     4923                       * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_0(k,j,i), 0 ) )
     4924
     4925                km(k,j,i) = c_0 * SQRT( e(k,j,i) ) * l_v(k)
     4926                kh(k,j,i) = km(k,j,i) / prandtl_number
     4927
     4928             ENDDO
     4929!
     4930!--          Summation for averaged profile of mixing length (cf. flow_statistics)
     4931             !$ACC LOOP PRIVATE(sr, k)
     4932             DO  sr = 0, statistic_regions
     4933                DO  k = nzb+1, nzt
     4934                   sums_l_l(k,sr,tn) = sums_l_l(k,sr,tn) + l_v(k) * rmask(j,i,sr)
    49524935                ENDDO
    4953 
    4954              ENDDO
     4936             ENDDO
     4937
    49554938          ENDDO
    49564939       ENDDO
Note: See TracChangeset for help on using the changeset viewer.