Changeset 4473 for palm/trunk/SOURCE
- Timestamp:
- Mar 25, 2020 9:04:07 PM (5 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 3 edited
-
header.f90 (modified) (3 diffs)
-
modules.f90 (modified) (2 diffs)
-
turbulence_closure_mod.f90 (modified) (55 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/header.f90
r4444 r4473 25 25 ! ----------------- 26 26 ! $Id$ 27 ! revised message if wall_adjustment is used 28 ! 29 ! 4444 2020-03-05 15:59:50Z raasch 27 30 ! bugfix: cpp-directives for serial mode added 28 31 ! … … 1539 1542 IF ( e_init > 0.0_wp ) WRITE ( io, 455 ) e_init 1540 1543 IF ( e_min > 0.0_wp ) WRITE ( io, 454 ) e_min 1541 IF ( wall_adjustment ) WRITE ( io, 453 ) wall_adjustment_factor1544 IF ( wall_adjustment ) WRITE ( io, 453 ) 1542 1545 ENDIF 1543 1546 IF ( rans_mode ) THEN … … 1884 1887 451 FORMAT (' Diffusion coefficients are constant:'/ & 1885 1888 ' 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')1889 453 FORMAT (' Mixing length is limited close to surfaces') 1887 1890 454 FORMAT (' TKE is not allowed to fall below ',E9.2,' (m/s)**2') 1888 1891 455 FORMAT (' initial TKE is prescribed as ',E9.2,' (m/s)**2') -
palm/trunk/SOURCE/modules.f90
r4472 r4473 25 25 ! ----------------- 26 26 ! $Id$ 27 ! moved wall_adjustment_factor to turbulence_closure_mod 28 ! 29 ! 4472 2020-03-24 12:21:00Z Giersch 27 30 ! Additional switch added to activate calculations in flow_statistics for the 28 31 ! kolmogorov length scale … … 961 964 REAL(wp) :: v_bulk = 0.0_wp !< namelist parameter 962 965 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 l964 966 REAL(wp) :: zeta_max = 20.0_wp !< namelist parameter 965 967 REAL(wp) :: zeta_min = -20.0_wp !< namelist parameter -
palm/trunk/SOURCE/turbulence_closure_mod.f90
r4457 r4473 25 25 ! ----------------- 26 26 ! $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 27 40 ! use statement for exchange horiz added 28 ! 41 ! 29 42 ! 4433 2020-02-28 22:14:43Z gronemeier 30 43 ! remove warning for newly implemented RANS mode … … 97 110 ! ------------ 98 111 !> This module contains the available turbulence closures for PALM. 99 !>100 112 !> 101 113 !> @todo test initialization for all possibilities … … 181 193 182 194 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) 193 206 (/ 0.55_wp, 1.44_wp, 1.92_wp, 1.44_wp, 0.0_wp /) !> default values fit for standard-tke-e closure 194 207 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 203 213 204 214 ! … … 1147 1157 DO j = nysg, nyng 1148 1158 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) ) 1150 1161 ENDDO 1151 1162 ENDDO … … 1191 1202 DO j = nysg, nyng 1192 1203 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) 1194 1205 ENDDO 1195 1206 ENDDO … … 1368 1379 1369 1380 USE control_parameters, & 1370 ONLY: f, message_string, wall_adjustment , wall_adjustment_factor1381 ONLY: f, message_string, wall_adjustment 1371 1382 1372 1383 USE exchange_horiz_mod, & … … 1406 1417 INTEGER(KIND=1), DIMENSION(:,:,:), ALLOCATABLE :: vicinity !< contains topography information of the vicinity of (i/j/k) 1407 1418 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) ) 1412 1436 ! 1413 1437 !-- Initialize the mixing length in case of an LES-simulation 1414 1438 IF ( .NOT. rans_mode ) THEN 1415 1439 ! 1416 !-- Compute the g rid-dependent mixing length.1417 DO k = 1, nzt1418 l_grid(k) = ( dx * dy * dzw(k) )**0.33333333333333_wp1440 !-- 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 1419 1443 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) 1426 1459 ENDDO 1427 l_wall(nzt+1,:,:) = l_grid(nzt)1428 1460 1429 1461 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 1449 1466 DO i = nxl, nxr 1450 1467 DO j = nys, nyn … … 1454 1471 IF ( BTEST( wall_flags_total_0(k,j,i), 0 ) ) THEN 1455 1472 ! 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) ) 1533 1584 1534 1585 ENDDO !k loop … … 1538 1589 ENDIF !if wall_adjustment 1539 1590 1540 ELSE 1591 ELSE !--> RANS mode 1541 1592 ! 1542 1593 !-- Initialize the mixing length in case of a RANS simulation 1543 ALLOCATE( l_black(nzb:nzt+1) )1594 ALLOCATE( ml_blackadar(nzb:nzt+1) ) 1544 1595 1545 1596 ! 1546 1597 !-- Calculate mixing length according to Blackadar (1962) 1547 1598 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_wp1599 length_scale_max = 2.7E-4_wp * SQRT( ug(nzt+1)**2 + vg(nzt+1)**2 ) & 1600 / ABS( f ) + 1.0E-10_wp 1550 1601 ELSE 1551 l _max = 30.0_wp1602 length_scale_max = 30.0_wp 1552 1603 ENDIF 1553 1604 1554 1605 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 ) 1556 1607 ENDDO 1557 1608 1558 l_black(nzt+1) = l_black(nzt)1609 ml_blackadar(nzt+1) = ml_blackadar(nzt) 1559 1610 1560 1611 ! … … 1571 1622 ENDDO 1572 1623 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) 1575 1626 ! 1576 1627 !-- Limit mixing length to either nearest wall or Blackadar mixing length. … … 1581 1632 DO k = nzb+1, nzt 1582 1633 1583 radius = l_black(k) ! radius within walls are searched1584 ! 1585 !-- Set l_wall to its default maximum value (l_back)1586 l_wall(k,:,:) = radius1634 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 1587 1638 1588 1639 ! … … 1679 1730 !-- Search for walls within octant (+++) 1680 1731 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), & 1682 1733 shortest_distance( vic_yz, .TRUE., ii ) ) 1683 1734 ! … … 1687 1738 !-- shortest_distance"). 1688 1739 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), & 1690 1741 shortest_distance( vic_yz, .TRUE., ii ) ) 1691 1742 … … 1694 1745 !-- Search for walls within octant (+--) 1695 1746 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), & 1697 1748 shortest_distance( vic_yz, .FALSE., ii ) ) 1698 1749 ! 1699 1750 !-- Search for walls within octant (++-) 1700 1751 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), & 1702 1753 shortest_distance( vic_yz, .FALSE., ii ) ) 1703 1754 ! 1704 1755 !-- 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 ) 1706 1757 1707 1758 ENDDO … … 1716 1767 !-- Search for walls within octant (-++) 1717 1768 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), & 1719 1770 shortest_distance( vic_yz, .TRUE., -ii ) ) 1720 1771 ! … … 1724 1775 !-- shortest_distance"). 1725 1776 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), & 1727 1778 shortest_distance( vic_yz, .TRUE., -ii ) ) 1728 1779 … … 1731 1782 !-- Search for walls within octant (---) 1732 1783 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), & 1734 1785 shortest_distance( vic_yz, .FALSE., -ii ) ) 1735 1786 ! 1736 1787 !-- Search for walls within octant (-+-) 1737 1788 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), & 1739 1790 shortest_distance( vic_yz, .FALSE., -ii ) ) 1740 1791 ! 1741 1792 !-- 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 ) 1743 1794 1744 1795 ENDDO … … 1748 1799 ELSE !Check if (i,j,k) belongs to atmosphere 1749 1800 1750 l_wall(k,j,i) = l_black(k)1801 ml_wall_adjusted(k,j,i) = ml_blackadar(k) 1751 1802 1752 1803 ENDIF … … 1762 1813 ENDDO !k loop 1763 1814 1764 !$ACC ENTER DATA COPYIN( l_black(nzb:nzt+1))1815 !$ACC ENTER DATA COPYIN(ml_blackadar(nzb:nzt+1)) 1765 1816 1766 1817 ENDIF !LES or RANS mode 1767 1818 1768 1819 ! 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)) 1774 1824 1775 1825 CONTAINS … … 2266 2316 IF ( plant_canopy ) CALL pcm_tendency( 6 ) 2267 2317 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' ) 2269 2320 2270 2321 ! … … 2377 2428 ! 2378 2429 !-- Additional sink term for flows through plant canopies 2379 ! IF ( plant_canopy ) CALL pcm_tendency( ? )!> @todo not yet implemented2380 2381 ! CALL user_actions( 'e-tendency' ) ToDo:find general solution for circular dependency between modules2430 ! 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 2382 2433 2383 2434 ! … … 2496 2547 IF ( plant_canopy ) CALL pcm_tendency( i, j, 6 ) 2497 2548 2498 ! CALL user_actions( i, j, 'e-tendency' ) ToDo: find general solution for circular dependency between modules2549 ! CALL user_actions( i, j, 'e-tendency' ) !> @todo: find general solution for circular dependency between modules 2499 2550 2500 2551 ! … … 2563 2614 ! IF ( plant_canopy ) CALL pcm_tendency( i, j, ? ) !> @todo not yet implemented 2564 2615 2565 ! CALL user_actions( i, j, 'diss-tendency' ) ToDo: find general solution for circular dependency between modules2616 ! CALL user_actions( i, j, 'diss-tendency' ) !> @todo: find general solution for circular dependency between modules 2566 2617 2567 2618 ! … … 2726 2777 2727 2778 km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp & 2728 * 0.5_wp * dy2779 * 0.5_wp * dy 2729 2780 ! 2730 2781 !-- -1.0 for right-facing wall, 1.0 for left-facing wall … … 2745 2796 2746 2797 km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp & 2747 * 0.5_wp * dy2798 * 0.5_wp * dy 2748 2799 ! 2749 2800 !-- -1.0 for right-facing wall, 1.0 for left-facing wall … … 2764 2815 2765 2816 km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp & 2766 * 0.5_wp * dy2817 * 0.5_wp * dy 2767 2818 ! 2768 2819 !-- -1.0 for right-facing wall, 1.0 for left-facing wall … … 3410 3461 3411 3462 km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp & 3412 * 0.5_wp * dy3463 * 0.5_wp * dy 3413 3464 ! 3414 3465 !-- -1.0 for right-facing wall, 1.0 for left-facing wall … … 3428 3479 3429 3480 km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp & 3430 * 0.5_wp * dy3481 * 0.5_wp * dy 3431 3482 ! 3432 3483 !-- -1.0 for right-facing wall, 1.0 for left-facing wall … … 3446 3497 3447 3498 km_neutral = kappa * ( usvs**2 + wsvs**2 )**0.25_wp & 3448 * 0.5_wp * dy3499 * 0.5_wp * dy 3449 3500 ! 3450 3501 !-- -1.0 for right-facing wall, 1.0 for left-facing wall … … 3911 3962 flag = MERGE( 1.0_wp, 0.0_wp, BTEST(wall_flags_total_0(k,j,i),0) ) 3912 3963 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), & 3914 3965 use_single_reference_value ) ) 3915 3966 ENDDO … … 3971 4022 REAL(wp) :: duv2_dz2 !< squared vertical gradient of wind vector 3972 4023 REAL(wp) :: dvar_dz !< vertical gradient of var 3973 REAL(wp) :: l!< mixing length4024 REAL(wp) :: ml !< mixing length 3974 4025 REAL(wp) :: var_reference !< reference temperature 3975 4026 3976 REAL(wp), DIMENSION(nzb+1:nzt) :: l_stable!< mixing length according to stratification3977 REAL(wp), DIMENSION(nzb+1:nzt) :: rif !< Richardson flux number4027 REAL(wp), DIMENSION(nzb+1:nzt) :: ml_stratification !< mixing length according to stratification 4028 REAL(wp), DIMENSION(nzb+1:nzt) :: rif !< Richardson flux number 3978 4029 3979 4030 REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) :: var !< temperature … … 3984 4035 3985 4036 !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) & 3986 !$ACC PRIVATE( l, l_stable, dvar_dz) &4037 !$ACC PRIVATE(ml, ml_stratification, dvar_dz) & 3987 4038 !$ACC PRESENT(diss, e, var, wall_flags_total_0) & 3988 !$ACC PRESENT(dd2zu, l_grid, l_wall)4039 !$ACC PRESENT(dd2zu, ml_wall_adjusted) 3989 4040 DO i = nxl, nxr 3990 4041 DO j = nys, nyn … … 3995 4046 IF ( dvar_dz > 0.0_wp ) THEN 3996 4047 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) ) & 3998 4049 / SQRT( g / var_reference * dvar_dz ) + 1E-5_wp 3999 4050 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) ) & 4001 4052 / SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp 4002 4053 ENDIF 4003 4054 ELSE 4004 l_stable(k) = l_grid(k)4055 ml_stratification(k) = ml_wall_adjusted(k,j,i) 4005 4056 ENDIF 4006 4057 … … 4011 4062 DO k = nzb+1, nzt 4012 4063 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 & 4017 4068 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 4018 4069 … … 4025 4076 4026 4077 !$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) & 4028 4079 !$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) 4030 4081 DO i = nxl, nxr 4031 4082 DO j = nys, nyn … … 4062 4113 DO k = nzb+1, nzt 4063 4114 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) ) 4065 4116 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) ) 4067 4118 ENDIF 4068 4119 ENDDO … … 4072 4123 DO k = nzb+1, nzt 4073 4124 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) ) & 4075 4127 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 4076 4128 … … 4180 4232 REAL(wp) :: duv2_dz2 !< squared vertical gradient of wind vector 4181 4233 REAL(wp) :: dvar_dz !< vertical gradient of var 4182 REAL(wp) :: l!< mixing length4234 REAL(wp) :: ml !< mixing length 4183 4235 REAL(wp) :: var_reference !< reference temperature 4184 4236 4185 REAL(wp), DIMENSION(nzb+1:nzt) :: l_stable!< mixing length according to stratification4186 REAL(wp), DIMENSION(nzb+1:nzt) :: rif !< Richardson flux number4237 REAL(wp), DIMENSION(nzb+1:nzt) :: ml_stratification !< mixing length according to stratification 4238 REAL(wp), DIMENSION(nzb+1:nzt) :: rif !< Richardson flux number 4187 4239 4188 4240 REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) :: var !< temperature … … 4197 4249 IF ( dvar_dz > 0.0_wp ) THEN 4198 4250 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) ) & 4200 4252 / SQRT( g / var_reference * dvar_dz ) + 1E-5_wp 4201 4253 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) ) & 4203 4255 / SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp 4204 4256 ENDIF 4205 4257 ELSE 4206 l_stable(k) = l_grid(k)4258 ml_stratification(k) = ml_wall_adjusted(k,j,i) 4207 4259 ENDIF 4208 4260 ENDDO … … 4210 4262 !DIR$ IVDEP 4211 4263 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 & 4216 4268 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 4217 4269 ENDDO … … 4248 4300 DO k = nzb+1, nzt 4249 4301 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) ) 4251 4303 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) ) 4253 4305 ENDIF 4254 4306 … … 4257 4309 !DIR$ IVDEP 4258 4310 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) ) & 4260 4313 * MERGE( 1.0_wp, 0.0_wp, BTEST( wall_flags_total_0(k,j,i), 0 ) ) 4261 4314 ENDDO … … 4652 4705 REAL(wp) :: duv2_dz2 !< squared vertical gradient of wind vector 4653 4706 REAL(wp) :: dvar_dz !< vertical gradient of var 4654 REAL(wp) :: l!< mixing length (single height)4707 REAL(wp) :: ml !< mixing length (single height) 4655 4708 REAL(wp) :: var_reference !< reference temperature 4656 4709 4657 !DIR$ ATTRIBUTES ALIGN:64:: l_v, l_stable, rif4658 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 stratification4660 REAL(wp), DIMENSION(nzb+1:nzt) :: rif !< Richardson flux number4710 !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 4661 4714 4662 4715 REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) :: var !< temperature … … 4674 4727 ! 4675 4728 !-- 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) 4677 4730 !$ tn = omp_get_thread_num() 4678 4731 … … 4680 4733 !$OMP DO 4681 4734 !$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) 4685 4738 DO i = nxlg, nxrg 4686 4739 DO j = nysg, nyng … … 4696 4749 IF ( dvar_dz > 0.0_wp ) THEN 4697 4750 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) ) & 4699 4752 / SQRT( g / var_reference * dvar_dz ) + 1E-5_wp 4700 4753 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) ) & 4702 4755 / SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp 4703 4756 ENDIF 4704 4757 ELSE 4705 l_stable(k) = l_grid(k)4758 ml_stratification(k) = ml_wall_adjusted(k,j,i) 4706 4759 ENDIF 4707 4760 … … 4712 4765 DO k = nzb+1, nzt 4713 4766 4714 l_v(k) = MIN( l_wall(k,j,i), l_stable(k) )&4715 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 4717 4770 ! 4718 4771 !-- 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) 4721 4774 4722 4775 ENDDO … … 4726 4779 DO sr = 0, statistic_regions 4727 4780 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) 4729 4782 ENDDO 4730 4783 ENDDO … … 4737 4790 !$OMP DO 4738 4791 !$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) 4742 4795 DO i = nxlg, nxrg 4743 4796 DO j = nysg, nyng … … 4775 4828 DO k = nzb+1, nzt 4776 4829 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) ) 4778 4831 ELSE 4779 l_stable(k) = l_wall(k,j,i)4832 ml_stratification(k) = ml_blackadar(k) 4780 4833 ENDIF 4781 4834 … … 4786 4839 !DIR$ IVDEP 4787 4840 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) ) 4790 4844 kh(k,j,i) = km(k,j,i) / prandtl_number 4791 4845 ENDDO … … 4795 4849 DO sr = 0, statistic_regions 4796 4850 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) 4798 4852 ENDDO 4799 4853 ENDDO … … 4806 4860 !$OMP DO 4807 4861 !$ACC PARALLEL LOOP COLLAPSE(2) PRIVATE(i, j) & 4808 !$ACC PRIVATE( l_v) &4862 !$ACC PRIVATE(ml_local_profile) & 4809 4863 !$ACC PRESENT(wall_flags_total_0, e, diss, rmask) & 4810 4864 !$ACC PRESENT(kh, km, sums_l_l) … … 4817 4871 DO k = nzb+1, nzt 4818 4872 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) 4823 4878 kh(k,j,i) = km(k,j,i) / prandtl_number 4824 4879 … … 4829 4884 DO sr = 0, statistic_regions 4830 4885 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) 4832 4887 ENDDO 4833 4888 ENDDO
Note: See TracChangeset
for help on using the changeset viewer.
