Changeset 4170 for palm/trunk
- Timestamp:
- Aug 19, 2019 5:12:31 PM (5 years ago)
- Location:
- palm/trunk/SOURCE
- Files:
-
- 2 edited
-
time_integration.f90 (modified) (5 diffs)
-
turbulence_closure_mod.f90 (modified) (28 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/time_integration.f90
r4144 r4170 25 25 ! ----------------- 26 26 ! $Id$ 27 ! copy diss, diss_p, tdiss_m to GPU 28 ! 29 ! 4144 2019-08-06 09:11:47Z raasch 27 30 ! relational operators .EQ., .NE., etc. replaced by ==, /=, etc. 28 31 ! … … 718 721 ONLY: d, dd2zu, ddzu, ddzw, drho_air, drho_air_zw, dzw, heatflux_output_conversion, kh, & 719 722 km, momentumflux_output_conversion, p, ptdf_x, ptdf_y, rdf, rdf_sc, rho_air, & 720 rho_air_zw, t e_m, tpt_m, tu_m, tv_m, tw_m, ug, u_init, u_stokes_zu, vg, v_init,&721 v_ stokes_zu, zu723 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 722 725 723 726 USE control_parameters, & … … 761 764 !$ACC DATA & 762 765 !$ACC COPY(d(nzb+1:nzt,nys:nyn,nxl:nxr)) & 766 !$ACC COPY(diss(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) & 763 767 !$ACC COPY(e(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) & 764 768 !$ACC COPY(u(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) & … … 771 775 772 776 !$ACC DATA & 777 !$ACC COPY(diss_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) & 773 778 !$ACC COPY(e_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) & 774 779 !$ACC COPY(u_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) & … … 777 782 !$ACC COPY(pt_p(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) & 778 783 !$ACC COPY(tend(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) & 784 !$ACC COPY(tdiss_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) & 779 785 !$ACC COPY(te_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) & 780 786 !$ACC COPY(tu_m(nzb:nzt+1,nysg:nyng,nxlg:nxrg)) & -
palm/trunk/SOURCE/turbulence_closure_mod.f90
r4168 r4170 25 25 ! ----------------- 26 26 ! $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 27 34 ! Replace function get_topography_top_index by topo_top_ind 28 35 ! … … 205 212 !> 206 213 !> @todo test initialization for all possibilities 207 !> 214 !> @todo add OpenMP directives whereever possible 208 215 !> @todo Check for random disturbances 209 216 !> @note <Enter notes on the module> … … 657 664 ENDDO 658 665 659 660 666 diss_p(nzt+1,:,:) = diss_p(nzt,:,:) 661 667 662 668 ENDIF 669 663 670 END SUBROUTINE tcm_boundary_conds 664 671 … … 1160 1167 ONLY: collision_turbulence 1161 1168 1162 USE particle_attributes, &1163 ONLY: use_sgs_for_particles, wang_kernel1164 1165 1169 USE pmc_interface, & 1166 1170 ONLY: nested_run … … 1180 1184 !-- they do not necessarily need to be transferred, which is attributed to 1181 1185 !-- 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) ) 1191 1191 ENDIF 1192 1192 … … 1195 1195 e => e_1; e_p => e_2; te_m => e_3 1196 1196 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 1201 1199 diss_p => diss_2; tdiss_m => diss_3 1202 ENDIF1203 1200 ENDIF 1204 1201 … … 1252 1249 e = 0.0_wp 1253 1250 ENDIF 1251 1252 diss = 0.0_wp 1254 1253 1255 1254 ELSE … … 1298 1297 ENDIF 1299 1298 1300 ELSEIF ( INDEX( initializing_actions, 'set_constant_profiles') /= 0 .OR. &1299 ELSEIF ( INDEX( initializing_actions, 'set_constant_profiles' ) /= 0 .OR. & 1301 1300 INDEX( initializing_actions, 'inifor' ) /= 0 ) THEN 1302 1301 … … 1340 1339 diss(nzb,:,:) = diss(nzb+1,:,:) 1341 1340 diss(nzt+1,:,:) = diss(nzt,:,:) 1341 ELSE 1342 diss = 0.0_wp 1342 1343 ENDIF 1343 1344 … … 1383 1384 ENDDO 1384 1385 ENDDO 1386 ELSE 1387 diss = 0.0_wp 1385 1388 ENDIF 1386 1389 ENDIF … … 1469 1472 1470 1473 1474 !------------------------------------------------------------------------------! 1471 1475 ! Description: 1472 ! ------------ -----------------------------------------------------------------!1476 ! ------------ 1473 1477 !> Pre-computation of grid-dependent and near-wall mixing length. 1474 1478 !> @todo consider walls in horizontal direction at a distance further than a … … 1481 1485 1482 1486 USE control_parameters, & 1483 ONLY: f, message_string, wall_adjustment _factor1487 ONLY: f, message_string, wall_adjustment, wall_adjustment_factor 1484 1488 1485 1489 USE grid_variables, & … … 1495 1499 IMPLICIT NONE 1496 1500 1497 INTEGER(iwp) :: dist_dx 1498 INTEGER(iwp) :: i 1499 INTEGER(iwp) :: ii 1500 INTEGER(iwp) :: j 1501 INTEGER(iwp) :: k 1502 INTEGER(iwp) :: k_max_topo = 0!< index of maximum topography height1503 INTEGER(iwp) :: kk 1504 INTEGER(iwp) :: rad_i 1505 INTEGER(iwp) :: rad_i_l 1506 INTEGER(iwp) :: rad_i_r 1507 INTEGER(iwp) :: rad_j 1508 INTEGER(iwp) :: rad_j_n 1509 INTEGER(iwp) :: rad_j_s 1510 INTEGER(iwp) :: rad_k 1511 INTEGER(iwp) :: rad_k_b 1512 INTEGER(iwp) :: rad_k_t 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 1513 1517 1514 1518 INTEGER(KIND=1), DIMENSION(:,:), ALLOCATABLE :: vic_yz !< contains a quarter of a single yz-slice of vicinity … … 1516 1520 INTEGER(KIND=1), DIMENSION(:,:,:), ALLOCATABLE :: vicinity !< contains topography information of the vicinity of (i/j/k) 1517 1521 1518 REAL(wp) :: radius 1522 REAL(wp) :: radius !< search radius in meter 1519 1523 1520 1524 ALLOCATE( l_grid(1:nzt) ) … … 1537 1541 l_wall(nzt+1,:,:) = l_grid(nzt) 1538 1542 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 1639 1557 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 1641 1653 1642 1654 ELSE 1643 1655 ! 1644 !-- Initialize the mixing length in case of a RANS simulation1656 !-- Initialize the mixing length in case of a RANS simulation 1645 1657 ALLOCATE( l_black(nzb:nzt+1) ) 1646 1658 … … 1662 1674 ! 1663 1675 !-- Get height level of highest topography within local subdomain 1676 k_max_topo = 0 1664 1677 DO i = nxlg, nxrg 1665 1678 DO j = nysg, nyng … … 4043 4056 SUBROUTINE diffusion_e( var, var_reference ) 4044 4057 4045 #if defined( _OPENACC )4046 4058 USE arrays_3d, & 4047 ONLY: dd zu, dd2zu, ddzw, drho_air, rho_air_zw4059 ONLY: dd2zu, ddzu, ddzw, drho_air, rho_air_zw 4048 4060 4049 4061 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 4056 4064 USE grid_variables, & 4057 4065 ONLY: ddx2, ddy2 … … 4070 4078 INTEGER(iwp) :: m !< running index surface elements 4071 4079 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 4073 4082 REAL(wp) :: l !< mixing length 4074 REAL(wp) :: ll !< adjusted l4075 4083 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 4097 4195 !$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) 4104 4198 DO i = nxl, nxr 4105 4199 DO j = nys, nyn 4106 4200 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) 4218 4221 4219 4222 ENDDO … … 4258 4261 4259 4262 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 4261 4267 4262 4268 USE grid_variables, & … … 4278 4284 INTEGER(iwp) :: surf_s !< Start index of surface elements at (j,i)-gridpoint 4279 4285 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 4281 4288 REAL(wp) :: l !< mixing length 4282 REAL(wp) :: ll !< adjusted l4283 4289 REAL(wp) :: var_reference !< reference temperature 4284 4290 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 4291 4376 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 4324 4378 tend(k,j,i) = tend(k,j,i) + ( & 4325 4379 ( & … … 4337 4391 * rho_air_zw(k-1) & 4338 4392 ) * 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) 4341 4397 4342 4398 ENDDO 4343 4399 4344 4400 ! 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. 4346 4408 IF ( .NOT. rans_tke_e .AND. ( use_sgs_for_particles .OR. wang_kernel & 4347 4409 .OR. collision_turbulence ) ) THEN 4348 DO k = nzb+1, nzt4349 diss(k,j,i) = dissipation(k) * MERGE( 1.0_wp, 0.0_wp, &4350 BTEST( wall_flags_0(k,j,i), 0 ) )4351 ENDDO4352 !4353 !-- Neumann boundary condition for dissipation diss(nzb,:,:) = diss(nzb+1,:,:)4354 !-- For each surface type determine start and end index (in case of elevated4355 !-- topography several up/downward facing surfaces may exist.4356 !-- Note, bc cannot be set in tcm_boundary conditions as the dissipation4357 !-- in LES mode is only a diagnostic quantity.4358 4410 surf_s = bc_h(0)%start_index(j,i) 4359 4411 surf_e = bc_h(0)%end_index(j,i) … … 4490 4542 ! Description: 4491 4543 ! ------------ 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: dd2zu4498 4499 USE control_parameters, &4500 ONLY: atmos_ocean_sign, use_single_reference_value, &4501 wall_adjustment, wall_adjustment_factor4502 4503 IMPLICIT NONE4504 4505 INTEGER(iwp) :: i !< loop index4506 INTEGER(iwp) :: j !< loop index4507 INTEGER(iwp) :: k !< loop index4508 4509 REAL(wp) :: dvar_dz !< vertical gradient of var4510 REAL(wp) :: l !< mixing length4511 REAL(wp) :: l_stable !< mixing length according to stratification4512 REAL(wp) :: ll !< adjusted l_grid4513 REAL(wp) :: var_reference !< var at reference height4514 4515 REAL(wp), DIMENSION(:,:,:), POINTER :: var !< temperature4516 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 ) THEN4520 IF ( use_single_reference_value ) THEN4521 l_stable = 0.76_wp * SQRT( e(k,j,i) ) &4522 / SQRT( g / var_reference * dvar_dz ) + 1E-5_wp4523 ELSE4524 l_stable = 0.76_wp * SQRT( e(k,j,i) ) &4525 / SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp4526 ENDIF4527 ELSE4528 l_stable = l_grid(k)4529 ENDIF4530 !4531 !-- Adjustment of the mixing length4532 IF ( wall_adjustment ) THEN4533 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 ELSE4536 l = MIN( l_grid(k), l_stable )4537 ll = l_grid(k)4538 ENDIF4539 4540 END SUBROUTINE mixing_length_les4541 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: dd2zu4552 4553 USE control_parameters, &4554 ONLY: atmos_ocean_sign, use_single_reference_value4555 4556 IMPLICIT NONE4557 4558 INTEGER(iwp) :: i !< loop index4559 INTEGER(iwp) :: j !< loop index4560 INTEGER(iwp) :: k !< loop index4561 4562 REAL(wp) :: duv2_dz2 !< squared vertical gradient of wind vector4563 REAL(wp) :: dvar_dz !< vertical gradient of var4564 REAL(wp) :: l !< mixing length4565 REAL(wp) :: l_diss !< mixing length for dissipation4566 REAL(wp) :: rif !< Richardson flux number4567 REAL(wp) :: var_reference !< var at reference height4568 4569 REAL(wp), DIMENSION(:,:,:), POINTER :: var !< temperature4570 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_wp4577 4578 IF ( use_single_reference_value ) THEN4579 rif = g / var_reference * dvar_dz / duv2_dz24580 ELSE4581 rif = g / var(k,j,i) * dvar_dz / duv2_dz24582 ENDIF4583 4584 rif = MAX( rif, -5.0_wp )4585 rif = MIN( rif, 1.0_wp )4586 4587 !4588 !-- Calculate diabatic mixing length using Dyer-profile functions4589 IF ( rif >= 0.0_wp ) THEN4590 l = MIN( l_black(k) / ( 1.0_wp + 5.0_wp * rif ), l_wall(k,j,i) )4591 l_diss = l4592 ELSE4593 !4594 !-- In case of unstable stratification, use mixing length of neutral case4595 !-- for l, but consider profile functions for l_diss4596 l = l_wall(k,j,i)4597 l_diss = l * SQRT( 1.0_wp - 16.0_wp * rif )4598 ENDIF4599 4600 END SUBROUTINE mixing_length_rans4601 4602 4603 !------------------------------------------------------------------------------!4604 ! Description:4605 ! ------------4606 4544 !> Computation of the turbulent diffusion coefficients for momentum and heat. 4607 4545 !> @bug unstable stratification is not properly considered for kh in rans mode. … … 4624 4562 REAL(wp) :: var_reference !< reference temperature 4625 4563 4626 REAL(wp), DIMENSION( :,:,:), POINTER:: var !< temperature4564 REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg), INTENT(IN) :: var !< temperature 4627 4565 4628 4566 … … 4797 4735 SUBROUTINE tcm_diffusivities_default( var, var_reference ) 4798 4736 4799 #if defined( _OPENACC )4800 4737 USE arrays_3d, & 4801 4738 ONLY: dd2zu 4802 4739 4803 4740 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 4807 4742 4808 4743 USE statistics, & … … 4818 4753 INTEGER(iwp) :: tn !< thread number 4819 4754 4820 REAL(wp) :: flag !< topography flag4821 REAL(wp) :: l !< mixing length4822 REAL(wp) :: l l !< adjusted mixing length4755 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) 4823 4758 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 4833 4766 4834 4767 ! … … 4844 4777 ! 4845 4778 !-- 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) 4847 4780 !$ tn = omp_get_thread_num() 4848 4781 4849 4782 IF ( les_dynamic .OR. les_mw ) THEN 4850 4783 !$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) & 4853 4786 !$ACC PRESENT(wall_flags_0, var, dd2zu, e, l_wall, l_grid, rmask) & 4854 4787 !$ACC PRESENT(kh, km, sums_l_l) 4855 4788 DO i = nxlg, nxrg 4856 4789 DO j = nysg, nyng 4790 !$ACC LOOP PRIVATE(k) 4857 4791 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). 4868 4798 dvar_dz = atmos_ocean_sign * ( var(k+1,j,i) - var(k-1,j,i) ) * dd2zu(k) 4869 4799 IF ( dvar_dz > 0.0_wp ) THEN 4870 4800 IF ( use_single_reference_value ) THEN 4871 l_stable = 0.76_wp * SQRT( e(k,j,i) )&4872 4801 l_stable(k) = 0.76_wp * SQRT( e(k,j,i) ) & 4802 / SQRT( g / var_reference * dvar_dz ) + 1E-5_wp 4873 4803 ELSE 4874 l_stable = 0.76_wp * SQRT( e(k,j,i) )&4875 4804 l_stable(k) = 0.76_wp * SQRT( e(k,j,i) ) & 4805 / SQRT( g / var(k,j,i) * dvar_dz ) + 1E-5_wp 4876 4806 ENDIF 4877 4807 ELSE 4878 l_stable = l_grid(k)4808 l_stable(k) = l_grid(k) 4879 4809 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) 4893 4820 ! 4894 4821 !-- 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) 4902 4832 ENDDO 4903 4904 ENDDO 4833 ENDDO 4834 4905 4835 ENDDO 4906 4836 ENDDO … … 4909 4839 4910 4840 !$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) 4911 4845 DO i = nxlg, nxrg 4912 4846 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) 4913 4878 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) 4928 4901 ENDDO 4929 4930 ENDDO 4902 ENDDO 4903 4931 4904 ENDDO 4932 4905 ENDDO … … 4935 4908 4936 4909 !$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) 4937 4914 DO i = nxlg, nxrg 4938 4915 DO j = nysg, nyng 4916 ! 4917 !-- Compute diffusion coefficients for momentum and heat 4918 !$ACC LOOP PRIVATE(k) 4919 !DIR$ IVDEP 4939 4920 DO k = nzb+1, nzt 4940 4921 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) 4952 4935 ENDDO 4953 4954 ENDDO 4936 ENDDO 4937 4955 4938 ENDDO 4956 4939 ENDDO
Note: See TracChangeset
for help on using the changeset viewer.
