Changeset 3551 for palm/trunk/SOURCE
- Timestamp:
- Nov 21, 2018 5:05:46 PM (6 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/advec_ws.f90 (modified) (93 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/advec_ws.f90
r3547 r3551 25 25 ! ----------------- 26 26 ! $Id$ 27 ! - Computation of vertical fluxes separated from computation of horizontal 28 ! fluxes. Loops are splitted in order to avoid indirect indexing and allow 29 ! for better vectorization. 30 ! - Accelerate code by remove type-conversions of ibits 31 ! - replace pointer definition in scalar routine by simple explicit definition 32 ! of the passed array 33 ! 34 ! 3547 2018-11-21 13:21:24Z suehring 27 35 ! variables documented 28 36 ! … … 193 201 ! 194 202 ! 888 2012-04-20 15:03:46Z suehring 195 ! Number of IBITS() calls with identical arguments is reduced.203 ! Number of REAL( IBITS() calls with identical arguments is reduced. 196 204 ! 197 205 ! 862 2012-03-26 14:21:38Z suehring … … 1200 1208 1201 1209 INTEGER(iwp) :: i !< grid index along x-direction 1202 INTEGER(iwp) :: ibit0 !< flag indicating 1st-order scheme along x-direction1203 INTEGER(iwp) :: ibit1 !< flag indicating 3rd-order scheme along x-direction1204 INTEGER(iwp) :: ibit2 !< flag indicating 5th-order scheme along x-direction1205 INTEGER(iwp) :: ibit3 !< flag indicating 1st-order scheme along y-direction1206 INTEGER(iwp) :: ibit4 !< flag indicating 3rd-order scheme along y-direction1207 INTEGER(iwp) :: ibit5 !< flag indicating 5th-order scheme along y-direction1208 INTEGER(iwp) :: ibit6 !< flag indicating 1st-order scheme along z-direction1209 INTEGER(iwp) :: ibit7 !< flag indicating 3rd-order scheme along z-direction1210 INTEGER(iwp) :: ibit8 !< flag indicating 5th-order scheme along z-direction1211 1210 INTEGER(iwp) :: i_omp !< leftmost index on subdomain, or in case of OpenMP, on thread 1212 1211 INTEGER(iwp) :: j !< grid index along y-direction … … 1217 1216 INTEGER(iwp) :: tn !< number of OpenMP thread 1218 1217 1218 REAL(wp) :: ibit0 !< flag indicating 1st-order scheme along x-direction 1219 REAL(wp) :: ibit1 !< flag indicating 3rd-order scheme along x-direction 1220 REAL(wp) :: ibit2 !< flag indicating 5th-order scheme along x-direction 1221 REAL(wp) :: ibit3 !< flag indicating 1st-order scheme along y-direction 1222 REAL(wp) :: ibit4 !< flag indicating 3rd-order scheme along y-direction 1223 REAL(wp) :: ibit5 !< flag indicating 5th-order scheme along y-direction 1224 REAL(wp) :: ibit6 !< flag indicating 1st-order scheme along z-direction 1225 REAL(wp) :: ibit7 !< flag indicating 3rd-order scheme along z-direction 1226 REAL(wp) :: ibit8 !< flag indicating 5th-order scheme along z-direction 1219 1227 REAL(wp) :: diss_d !< artificial dissipation term at grid box bottom 1220 1228 REAL(wp) :: div !< diverence on scalar grid … … 1222 1230 REAL(wp) :: u_comp !< advection velocity along x-direction 1223 1231 REAL(wp) :: v_comp !< advection velocity along y-direction 1224 1225 #if defined( __nopointer ) 1226 REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: sk !< advected scalar 1227 #else 1228 REAL(wp), DIMENSION(:,:,:), POINTER :: sk !< advected scalar 1229 #endif 1232 ! 1233 !-- sk is an array from parameter list. It should not be a pointer, because 1234 !-- in that case the compiler can not assume a stride 1 and cannot perform 1235 !-- a strided one vector load. Adding the CONTIGUOUS keyword makes things 1236 !-- even worse, because the compiler cannot assume strided one in the 1237 !-- caller side. 1238 REAL(wp), INTENT(IN),DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: sk !< advected scalar 1239 1230 1240 REAL(wp), DIMENSION(nzb:nzt+1) :: diss_n !< discretized artificial dissipation at northward-side of the grid box 1231 1241 REAL(wp), DIMENSION(nzb:nzt+1) :: diss_r !< discretized artificial dissipation at rightward-side of the grid box … … 1249 1259 DO k = nzb+1, nzb_max 1250 1260 1251 ibit5 = IBITS(advc_flags_1(k,j-1,i),5,1)1252 ibit4 = IBITS(advc_flags_1(k,j-1,i),4,1)1253 ibit3 = IBITS(advc_flags_1(k,j-1,i),3,1)1261 ibit5 = REAL( IBITS(advc_flags_1(k,j-1,i),5,1), KIND = wp ) 1262 ibit4 = REAL( IBITS(advc_flags_1(k,j-1,i),4,1), KIND = wp ) 1263 ibit3 = REAL( IBITS(advc_flags_1(k,j-1,i),3,1), KIND = wp ) 1254 1264 1255 1265 v_comp = v(k,j,i) - v_gtrans + v_stokes_zu(k) … … 1310 1320 DO k = nzb+1, nzb_max 1311 1321 1312 ibit2 = IBITS(advc_flags_1(k,j,i-1),2,1)1313 ibit1 = IBITS(advc_flags_1(k,j,i-1),1,1)1314 ibit0 = IBITS(advc_flags_1(k,j,i-1),0,1)1322 ibit2 = REAL( IBITS(advc_flags_1(k,j,i-1),2,1), KIND = wp ) 1323 ibit1 = REAL( IBITS(advc_flags_1(k,j,i-1),1,1), KIND = wp ) 1324 ibit0 = REAL( IBITS(advc_flags_1(k,j,i-1),0,1), KIND = wp ) 1315 1325 1316 1326 u_comp = u(k,j,i) - u_gtrans + u_stokes_zu(k) … … 1365 1375 1366 1376 ENDIF 1367 1368 flux_t(0) = 0.0_wp1369 diss_t(0) = 0.0_wp1370 1377 ! 1371 1378 !-- Now compute the fluxes and tendency terms for the horizontal and … … 1377 1384 !-- flux at the end. 1378 1385 1379 ibit2 = IBITS(advc_flags_1(k,j,i),2,1)1380 ibit1 = IBITS(advc_flags_1(k,j,i),1,1)1381 ibit0 = IBITS(advc_flags_1(k,j,i),0,1)1386 ibit2 = REAL( IBITS(advc_flags_1(k,j,i),2,1), KIND = wp ) 1387 ibit1 = REAL( IBITS(advc_flags_1(k,j,i),1,1), KIND = wp ) 1388 ibit0 = REAL( IBITS(advc_flags_1(k,j,i),0,1), KIND = wp ) 1382 1389 1383 1390 u_comp = u(k,j,i+1) - u_gtrans + u_stokes_zu(k) … … 1412 1419 ) 1413 1420 1414 ibit5 = IBITS(advc_flags_1(k,j,i),5,1)1415 ibit4 = IBITS(advc_flags_1(k,j,i),4,1)1416 ibit3 = IBITS(advc_flags_1(k,j,i),3,1)1421 ibit5 = REAL( IBITS(advc_flags_1(k,j,i),5,1), KIND = wp ) 1422 ibit4 = REAL( IBITS(advc_flags_1(k,j,i),4,1), KIND = wp ) 1423 ibit3 = REAL( IBITS(advc_flags_1(k,j,i),3,1), KIND = wp ) 1417 1424 1418 1425 v_comp = v(k,j+1,i) - v_gtrans + v_stokes_zu(k) … … 1446 1453 ( sk(k,j+3,i) - sk(k,j-2,i) ) & 1447 1454 ) 1448 !1449 !-- k index has to be modified near bottom and top, else array1450 !-- subscripts will be exceeded.1451 ibit8 = IBITS(advc_flags_1(k,j,i),8,1)1452 ibit7 = IBITS(advc_flags_1(k,j,i),7,1)1453 ibit6 = IBITS(advc_flags_1(k,j,i),6,1)1454 1455 k_ppp = k + 3 * ibit81456 k_pp = k + 2 * ( 1 - ibit6 )1457 k_mm = k - 2 * ibit81458 1459 1460 flux_t(k) = w(k,j,i) * rho_air_zw(k) * ( &1461 ( 37.0_wp * ibit8 * adv_sca_5 &1462 + 7.0_wp * ibit7 * adv_sca_3 &1463 + ibit6 * adv_sca_1 &1464 ) * &1465 ( sk(k+1,j,i) + sk(k,j,i) ) &1466 - ( 8.0_wp * ibit8 * adv_sca_5 &1467 + ibit7 * adv_sca_3 &1468 ) * &1469 ( sk(k_pp,j,i) + sk(k-1,j,i) ) &1470 + ( ibit8 * adv_sca_5 &1471 ) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) &1472 )1473 1474 diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * ( &1475 ( 10.0_wp * ibit8 * adv_sca_5 &1476 + 3.0_wp * ibit7 * adv_sca_3 &1477 + ibit6 * adv_sca_1 &1478 ) * &1479 ( sk(k+1,j,i) - sk(k,j,i) ) &1480 - ( 5.0_wp * ibit8 * adv_sca_5 &1481 + ibit7 * adv_sca_3 &1482 ) * &1483 ( sk(k_pp,j,i) - sk(k-1,j,i) ) &1484 + ( ibit8 * adv_sca_5 &1485 ) * &1486 ( sk(k_ppp,j,i) - sk(k_mm,j,i) ) &1487 )1488 ENDDO1489 1490 DO k = nzb+1, nzb_max1491 1492 flux_d = flux_t(k-1)1493 diss_d = diss_t(k-1)1494 !1495 !-- Calculate the divergence of the velocity field. A respective1496 !-- correction is needed to overcome numerical instabilities caused1497 !-- by a not sufficient reduction of divergences near topography.1498 div = ( u(k,j,i+1) * ( ibit0 + ibit1 + ibit2 ) &1499 - u(k,j,i) * ( IBITS(advc_flags_1(k,j,i-1),0,1) &1500 + IBITS(advc_flags_1(k,j,i-1),1,1) &1501 + IBITS(advc_flags_1(k,j,i-1),2,1) &1502 ) &1503 ) * ddx &1504 + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 ) &1505 - v(k,j,i) * ( IBITS(advc_flags_1(k,j-1,i),3,1) &1506 + IBITS(advc_flags_1(k,j-1,i),4,1) &1507 + IBITS(advc_flags_1(k,j-1,i),5,1) &1508 ) &1509 ) * ddy &1510 + ( w(k,j,i) * rho_air_zw(k) * &1511 ( ibit6 + ibit7 + ibit8 ) &1512 - w(k-1,j,i) * rho_air_zw(k-1) * &1513 ( IBITS(advc_flags_1(k-1,j,i),6,1) &1514 + IBITS(advc_flags_1(k-1,j,i),7,1) &1515 + IBITS(advc_flags_1(k-1,j,i),8,1) &1516 ) &1517 ) * drho_air(k) * ddzw(k)1518 1519 1520 tend(k,j,i) = tend(k,j,i) - ( &1521 ( flux_r(k) + diss_r(k) - swap_flux_x_local(k,j,tn) - &1522 swap_diss_x_local(k,j,tn) ) * ddx &1523 + ( flux_n(k) + diss_n(k) - swap_flux_y_local(k,tn) - &1524 swap_diss_y_local(k,tn) ) * ddy &1525 + ( ( flux_t(k) + diss_t(k) ) - &1526 ( flux_d + diss_d ) &1527 ) * drho_air(k) * ddzw(k) &1528 ) + sk(k,j,i) * div1529 1530 swap_flux_y_local(k,tn) = flux_n(k)1531 swap_diss_y_local(k,tn) = diss_n(k)1532 swap_flux_x_local(k,j,tn) = flux_r(k)1533 swap_diss_x_local(k,j,tn) = diss_r(k)1534 1535 1455 ENDDO 1536 1456 ! … … 1559 1479 - 5.0_wp * ( sk(k,j+2,i) - sk(k,j-1,i) ) & 1560 1480 + ( sk(k,j+3,i) - sk(k,j-2,i) ) ) * adv_sca_5 1481 1482 ENDDO 1483 ! 1484 !-- Now, compute vertical fluxes. Split loop into a part treating the 1485 !-- lowest 2 grid points with indirect indexing, a main loop without 1486 !-- indirect indexing, and a loop for the uppermost 2 grip points with 1487 !-- indirect indexing. This allows better vectorization for the main loop. 1488 !-- First, compute the flux at model surface, which need has to be 1489 !-- calculated explicetely for the tendency at 1490 !-- the first w-level. For topography wall this is done implicitely by 1491 !-- advc_flags_1. 1492 flux_t(nzb) = 0.0_wp 1493 diss_t(nzb) = 0.0_wp 1494 DO k = nzb+1, nzb+2 1495 ibit8 = REAL( IBITS(advc_flags_1(k,j,i),8,1), KIND = wp ) 1496 ibit7 = REAL( IBITS(advc_flags_1(k,j,i),7,1), KIND = wp ) 1497 ibit6 = REAL( IBITS(advc_flags_1(k,j,i),6,1), KIND = wp ) 1561 1498 ! 1562 1499 !-- k index has to be modified near bottom and top, else array 1563 1500 !-- subscripts will be exceeded. 1564 ibit8 = IBITS(advc_flags_1(k,j,i),8,1)1565 ibit7 = IBITS(advc_flags_1(k,j,i),7,1)1566 ibit6 = IBITS(advc_flags_1(k,j,i),6,1)1567 1568 1501 k_ppp = k + 3 * ibit8 1569 1502 k_pp = k + 2 * ( 1 - ibit6 ) 1570 1503 k_mm = k - 2 * ibit8 1571 1504 1505 1572 1506 flux_t(k) = w(k,j,i) * rho_air_zw(k) * ( & 1573 ( 37.0_wp * ibit8 * adv_sca_5&1574 + 7.0_wp * ibit7 * adv_sca_3&1575 + ibit6 * adv_sca_1&1576 ) *&1577 ( sk(k+1,j,i) + sk(k,j,i) )&1578 - ( 8.0_wp * ibit8 * adv_sca_5&1507 ( 37.0_wp * ibit8 * adv_sca_5 & 1508 + 7.0_wp * ibit7 * adv_sca_3 & 1509 + ibit6 * adv_sca_1 & 1510 ) * & 1511 ( sk(k+1,j,i) + sk(k,j,i) ) & 1512 - ( 8.0_wp * ibit8 * adv_sca_5 & 1579 1513 + ibit7 * adv_sca_3 & 1580 ) *&1581 ( sk(k_pp,j,i) + sk(k-1,j,i) )&1582 + ( ibit8 * adv_sca_5&1583 ) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) )&1514 ) * & 1515 ( sk(k_pp,j,i) + sk(k-1,j,i) ) & 1516 + ( ibit8 * adv_sca_5 & 1517 ) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) & 1584 1518 ) 1585 1519 1586 1520 diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * ( & 1587 ( 10.0_wp * ibit8 * adv_sca_5&1588 + 3.0_wp * ibit7 * adv_sca_3&1589 + ibit6 * adv_sca_1&1590 ) *&1521 ( 10.0_wp * ibit8 * adv_sca_5 & 1522 + 3.0_wp * ibit7 * adv_sca_3 & 1523 + ibit6 * adv_sca_1 & 1524 ) * & 1591 1525 ( sk(k+1,j,i) - sk(k,j,i) ) & 1592 - ( 5.0_wp * ibit8 * adv_sca_5&1593 + ibit7 * adv_sca_3&1594 ) *&1526 - ( 5.0_wp * ibit8 * adv_sca_5 & 1527 + ibit7 * adv_sca_3 & 1528 ) * & 1595 1529 ( sk(k_pp,j,i) - sk(k-1,j,i) ) & 1596 + ( ibit8 * adv_sca_5&1597 ) *&1530 + ( ibit8 * adv_sca_5 & 1531 ) * & 1598 1532 ( sk(k_ppp,j,i) - sk(k_mm,j,i) ) & 1599 1533 ) 1600 1601 ENDDO 1602 1603 DO k = nzb_max+1, nzt 1534 ENDDO 1535 1536 DO k = nzb+3, nzt-2 1537 ibit8 = REAL( IBITS(advc_flags_1(k,j,i),8,1), KIND = wp ) 1538 ibit7 = REAL( IBITS(advc_flags_1(k,j,i),7,1), KIND = wp ) 1539 ibit6 = REAL( IBITS(advc_flags_1(k,j,i),6,1), KIND = wp ) 1540 1541 flux_t(k) = w(k,j,i) * rho_air_zw(k) * ( & 1542 ( 37.0_wp * ibit8 * adv_sca_5 & 1543 + 7.0_wp * ibit7 * adv_sca_3 & 1544 + ibit6 * adv_sca_1 & 1545 ) * & 1546 ( sk(k+1,j,i) + sk(k,j,i) ) & 1547 - ( 8.0_wp * ibit8 * adv_sca_5 & 1548 + ibit7 * adv_sca_3 & 1549 ) * & 1550 ( sk(k+2,j,i) + sk(k-1,j,i) ) & 1551 + ( ibit8 * adv_sca_5 & 1552 ) * ( sk(k+3,j,i)+ sk(k-2,j,i) ) & 1553 ) 1554 1555 diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * ( & 1556 ( 10.0_wp * ibit8 * adv_sca_5 & 1557 + 3.0_wp * ibit7 * adv_sca_3 & 1558 + ibit6 * adv_sca_1 & 1559 ) * & 1560 ( sk(k+1,j,i) - sk(k,j,i) ) & 1561 - ( 5.0_wp * ibit8 * adv_sca_5 & 1562 + ibit7 * adv_sca_3 & 1563 ) * & 1564 ( sk(k+2,j,i) - sk(k-1,j,i) ) & 1565 + ( ibit8 * adv_sca_5 & 1566 ) * & 1567 ( sk(k+3,j,i) - sk(k-2,j,i) ) & 1568 ) 1569 ENDDO 1570 1571 DO k = nzt-1, nzt 1572 ibit8 = REAL( IBITS(advc_flags_1(k,j,i),8,1), KIND = wp ) 1573 ibit7 = REAL( IBITS(advc_flags_1(k,j,i),7,1), KIND = wp ) 1574 ibit6 = REAL( IBITS(advc_flags_1(k,j,i),6,1), KIND = wp ) 1575 ! 1576 !-- k index has to be modified near bottom and top, else array 1577 !-- subscripts will be exceeded. 1578 k_ppp = k + 3 * ibit8 1579 k_pp = k + 2 * ( 1 - ibit6 ) 1580 k_mm = k - 2 * ibit8 1581 1582 1583 flux_t(k) = w(k,j,i) * rho_air_zw(k) * ( & 1584 ( 37.0_wp * ibit8 * adv_sca_5 & 1585 + 7.0_wp * ibit7 * adv_sca_3 & 1586 + ibit6 * adv_sca_1 & 1587 ) * & 1588 ( sk(k+1,j,i) + sk(k,j,i) ) & 1589 - ( 8.0_wp * ibit8 * adv_sca_5 & 1590 + ibit7 * adv_sca_3 & 1591 ) * & 1592 ( sk(k_pp,j,i) + sk(k-1,j,i) ) & 1593 + ( ibit8 * adv_sca_5 & 1594 ) * ( sk(k_ppp,j,i)+ sk(k_mm,j,i) ) & 1595 ) 1596 1597 diss_t(k) = -ABS( w(k,j,i) ) * rho_air_zw(k) * ( & 1598 ( 10.0_wp * ibit8 * adv_sca_5 & 1599 + 3.0_wp * ibit7 * adv_sca_3 & 1600 + ibit6 * adv_sca_1 & 1601 ) * & 1602 ( sk(k+1,j,i) - sk(k,j,i) ) & 1603 - ( 5.0_wp * ibit8 * adv_sca_5 & 1604 + ibit7 * adv_sca_3 & 1605 ) * & 1606 ( sk(k_pp,j,i) - sk(k-1,j,i) ) & 1607 + ( ibit8 * adv_sca_5 & 1608 ) * & 1609 ( sk(k_ppp,j,i) - sk(k_mm,j,i) ) & 1610 ) 1611 ENDDO 1612 1613 DO k = nzb+1, nzt 1604 1614 1605 1615 flux_d = flux_t(k-1) … … 1610 1620 !-- correction is needed to overcome numerical instabilities introduced 1611 1621 !-- by a not sufficient reduction of divergences near topography. 1612 div = ( u(k,j,i+1) - u(k,j,i) ) * ddx & 1613 + ( v(k,j+1,i) - v(k,j,i) ) * ddy & 1614 + ( w(k,j,i) * rho_air_zw(k) - & 1615 w(k-1,j,i) * rho_air_zw(k-1) & 1622 div = ( u(k,j,i+1) * ( ibit0 + ibit1 + ibit2 ) & 1623 - u(k,j,i) * ( & 1624 REAL( IBITS(advc_flags_1(k,j,i-1),0,1), KIND = wp ) & 1625 + REAL( IBITS(advc_flags_1(k,j,i-1),1,1), KIND = wp ) & 1626 + REAL( IBITS(advc_flags_1(k,j,i-1),2,1), KIND = wp ) & 1627 ) & 1628 ) * ddx & 1629 + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 ) & 1630 - v(k,j,i) * ( & 1631 REAL( IBITS(advc_flags_1(k,j-1,i),3,1), KIND = wp ) & 1632 + REAL( IBITS(advc_flags_1(k,j-1,i),4,1), KIND = wp ) & 1633 + REAL( IBITS(advc_flags_1(k,j-1,i),5,1), KIND = wp ) & 1634 ) & 1635 ) * ddy & 1636 + ( w(k,j,i) * rho_air_zw(k) * & 1637 ( ibit6 + ibit7 + ibit8 ) & 1638 - w(k-1,j,i) * rho_air_zw(k-1) * & 1639 ( & 1640 REAL( IBITS(advc_flags_1(k-1,j,i),6,1), KIND = wp ) & 1641 + REAL( IBITS(advc_flags_1(k-1,j,i),7,1), KIND = wp ) & 1642 + REAL( IBITS(advc_flags_1(k-1,j,i),8,1), KIND = wp ) & 1643 ) & 1616 1644 ) * drho_air(k) * ddzw(k) 1617 1645 … … 1777 1805 1778 1806 INTEGER(iwp) :: i !< grid index along x-direction 1779 INTEGER(iwp) :: ibit9 !< flag indicating 1st-order scheme along x-direction1780 INTEGER(iwp) :: ibit10 !< flag indicating 3rd-order scheme along x-direction1781 INTEGER(iwp) :: ibit11 !< flag indicating 5th-order scheme along x-direction1782 INTEGER(iwp) :: ibit12 !< flag indicating 1st-order scheme along y-direction1783 INTEGER(iwp) :: ibit13 !< flag indicating 3rd-order scheme along y-direction1784 INTEGER(iwp) :: ibit14 !< flag indicating 5th-order scheme along y-direction1785 INTEGER(iwp) :: ibit15 !< flag indicating 1st-order scheme along z-direction1786 INTEGER(iwp) :: ibit16 !< flag indicating 3rd-order scheme along z-direction1787 INTEGER(iwp) :: ibit17 !< flag indicating 5th-order scheme along z-direction1788 1807 INTEGER(iwp) :: i_omp !< leftmost index on subdomain, or in case of OpenMP, on thread 1789 1808 INTEGER(iwp) :: j !< grid index along y-direction … … 1794 1813 INTEGER(iwp) :: tn !< number of OpenMP thread 1795 1814 1815 REAL(wp) :: ibit9 !< flag indicating 1st-order scheme along x-direction 1816 REAL(wp) :: ibit10 !< flag indicating 3rd-order scheme along x-direction 1817 REAL(wp) :: ibit11 !< flag indicating 5th-order scheme along x-direction 1818 REAL(wp) :: ibit12 !< flag indicating 1st-order scheme along y-direction 1819 REAL(wp) :: ibit13 !< flag indicating 3rd-order scheme along y-direction 1820 REAL(wp) :: ibit14 !< flag indicating 5th-order scheme along y-direction 1821 REAL(wp) :: ibit15 !< flag indicating 1st-order scheme along z-direction 1822 REAL(wp) :: ibit16 !< flag indicating 3rd-order scheme along z-direction 1823 REAL(wp) :: ibit17 !< flag indicating 5th-order scheme along z-direction 1796 1824 REAL(wp) :: diss_d !< artificial dissipation term at grid box bottom 1797 1825 REAL(wp) :: div !< diverence on u-grid … … 1819 1847 DO k = nzb+1, nzb_max 1820 1848 1821 ibit14 = IBITS(advc_flags_1(k,j-1,i),14,1)1822 ibit13 = IBITS(advc_flags_1(k,j-1,i),13,1)1823 ibit12 = IBITS(advc_flags_1(k,j-1,i),12,1)1849 ibit14 = REAL( IBITS(advc_flags_1(k,j-1,i),14,1), KIND = wp ) 1850 ibit13 = REAL( IBITS(advc_flags_1(k,j-1,i),13,1), KIND = wp ) 1851 ibit12 = REAL( IBITS(advc_flags_1(k,j-1,i),12,1), KIND = wp ) 1824 1852 1825 1853 v_comp(k) = v(k,j,i) + v(k,j,i-1) - gv … … 1877 1905 DO k = nzb+1, nzb_max 1878 1906 1879 ibit11 = IBITS(advc_flags_1(k,j,i-1),11,1)1880 ibit10 = IBITS(advc_flags_1(k,j,i-1),10,1)1881 ibit9 = IBITS(advc_flags_1(k,j,i-1),9,1)1907 ibit11 = REAL( IBITS(advc_flags_1(k,j,i-1),11,1), KIND = wp ) 1908 ibit10 = REAL( IBITS(advc_flags_1(k,j,i-1),10,1), KIND = wp ) 1909 ibit9 = REAL( IBITS(advc_flags_1(k,j,i-1),9,1), KIND = wp ) 1882 1910 1883 1911 u_comp_l = u(k,j,i) + u(k,j,i-1) - gu … … 1929 1957 1930 1958 ENDIF 1931 1932 flux_t(0) = 0.0_wp1933 diss_t(0) = 0.0_wp1934 w_comp(0) = 0.0_wp1935 1959 ! 1936 1960 !-- Now compute the fluxes tendency terms for the horizontal and … … 1938 1962 DO k = nzb+1, nzb_max 1939 1963 1940 ibit11 = IBITS(advc_flags_1(k,j,i),11,1)1941 ibit10 = IBITS(advc_flags_1(k,j,i),10,1)1942 ibit9 = IBITS(advc_flags_1(k,j,i),9,1)1964 ibit11 = REAL( IBITS(advc_flags_1(k,j,i),11,1), KIND = wp ) 1965 ibit10 = REAL( IBITS(advc_flags_1(k,j,i),10,1), KIND = wp ) 1966 ibit9 = REAL( IBITS(advc_flags_1(k,j,i),9,1), KIND = wp ) 1943 1967 1944 1968 u_comp(k) = u(k,j,i+1) + u(k,j,i) … … 1973 1997 ) 1974 1998 1975 ibit14 = IBITS(advc_flags_1(k,j,i),14,1)1976 ibit13 = IBITS(advc_flags_1(k,j,i),13,1)1977 ibit12 = IBITS(advc_flags_1(k,j,i),12,1)1999 ibit14 = REAL( IBITS(advc_flags_1(k,j,i),14,1), KIND = wp ) 2000 ibit13 = REAL( IBITS(advc_flags_1(k,j,i),13,1), KIND = wp ) 2001 ibit12 = REAL( IBITS(advc_flags_1(k,j,i),12,1), KIND = wp ) 1978 2002 1979 2003 v_comp(k) = v(k,j+1,i) + v(k,j+1,i-1) - gv … … 2007 2031 ( u(k,j+3,i) - u(k,j-2,i) ) & 2008 2032 ) 2033 ENDDO 2034 2035 DO k = nzb_max+1, nzt 2036 2037 u_comp(k) = u(k,j,i+1) + u(k,j,i) 2038 flux_r(k) = ( u_comp(k) - gu ) * ( & 2039 37.0_wp * ( u(k,j,i+1) + u(k,j,i) ) & 2040 - 8.0_wp * ( u(k,j,i+2) + u(k,j,i-1) ) & 2041 + ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_5 2042 diss_r(k) = - ABS( u_comp(k) - gu ) * ( & 2043 10.0_wp * ( u(k,j,i+1) - u(k,j,i) ) & 2044 - 5.0_wp * ( u(k,j,i+2) - u(k,j,i-1) ) & 2045 + ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_5 2046 2047 v_comp(k) = v(k,j+1,i) + v(k,j+1,i-1) - gv 2048 flux_n(k) = v_comp(k) * ( & 2049 37.0_wp * ( u(k,j+1,i) + u(k,j,i) ) & 2050 - 8.0_wp * ( u(k,j+2,i) + u(k,j-1,i) ) & 2051 + ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_5 2052 diss_n(k) = - ABS( v_comp(k) ) * ( & 2053 10.0_wp * ( u(k,j+1,i) - u(k,j,i) ) & 2054 - 5.0_wp * ( u(k,j+2,i) - u(k,j-1,i) ) & 2055 + ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_5 2056 2057 ENDDO 2058 ! 2059 !-- Now, compute vertical fluxes. Split loop into a part treating the 2060 !-- lowest 2 grid points with indirect indexing, a main loop without 2061 !-- indirect indexing, and a loop for the uppermost 2 grip points with 2062 !-- indirect indexing. This allows better vectorization for the main loop. 2063 !-- First, compute the flux at model surface, which need has to be 2064 !-- calculated explicetely for the tendency at 2065 !-- the first w-level. For topography wall this is done implicitely by 2066 !-- advc_flags_1. 2067 flux_t(nzb) = 0.0_wp 2068 diss_t(nzb) = 0.0_wp 2069 w_comp(nzb) = 0.0_wp 2070 DO k = nzb+1, nzb+2 2009 2071 ! 2010 2072 !-- k index has to be modified near bottom and top, else array 2011 2073 !-- subscripts will be exceeded. 2012 ibit17 = IBITS(advc_flags_1(k,j,i),17,1)2013 ibit16 = IBITS(advc_flags_1(k,j,i),16,1)2014 ibit15 = IBITS(advc_flags_1(k,j,i),15,1)2074 ibit17 = REAL( IBITS(advc_flags_1(k,j,i),17,1), KIND = wp ) 2075 ibit16 = REAL( IBITS(advc_flags_1(k,j,i),16,1), KIND = wp ) 2076 ibit15 = REAL( IBITS(advc_flags_1(k,j,i),15,1), KIND = wp ) 2015 2077 2016 2078 k_ppp = k + 3 * ibit17 … … 2049 2111 ) 2050 2112 ENDDO 2051 2052 DO k = nzb+1, nzb_max 2113 2114 DO k = nzb+3, nzt-2 2115 2116 ibit17 = REAL( IBITS(advc_flags_1(k,j,i),17,1), KIND = wp ) 2117 ibit16 = REAL( IBITS(advc_flags_1(k,j,i),16,1), KIND = wp ) 2118 ibit15 = REAL( IBITS(advc_flags_1(k,j,i),15,1), KIND = wp ) 2119 2120 w_comp(k) = w(k,j,i) + w(k,j,i-1) 2121 flux_t(k) = w_comp(k) * rho_air_zw(k) * ( & 2122 ( 37.0_wp * ibit17 * adv_mom_5 & 2123 + 7.0_wp * ibit16 * adv_mom_3 & 2124 + ibit15 * adv_mom_1 & 2125 ) * & 2126 ( u(k+1,j,i) + u(k,j,i) ) & 2127 - ( 8.0_wp * ibit17 * adv_mom_5 & 2128 + ibit16 * adv_mom_3 & 2129 ) * & 2130 ( u(k+2,j,i) + u(k-1,j,i) ) & 2131 + ( ibit17 * adv_mom_5 & 2132 ) * & 2133 ( u(k+3,j,i) + u(k-2,j,i) ) & 2134 ) 2135 2136 diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * ( & 2137 ( 10.0_wp * ibit17 * adv_mom_5 & 2138 + 3.0_wp * ibit16 * adv_mom_3 & 2139 + ibit15 * adv_mom_1 & 2140 ) * & 2141 ( u(k+1,j,i) - u(k,j,i) ) & 2142 - ( 5.0_wp * ibit17 * adv_mom_5 & 2143 + ibit16 * adv_mom_3 & 2144 ) * & 2145 ( u(k+2,j,i) - u(k-1,j,i) ) & 2146 + ( ibit17 * adv_mom_5 & 2147 ) * & 2148 ( u(k+3,j,i) - u(k-2,j,i) ) & 2149 ) 2150 ENDDO 2151 2152 DO k = nzt-1, nzt 2153 ! 2154 !-- k index has to be modified near bottom and top, else array 2155 !-- subscripts will be exceeded. 2156 ibit17 = REAL( IBITS(advc_flags_1(k,j,i),17,1), KIND = wp ) 2157 ibit16 = REAL( IBITS(advc_flags_1(k,j,i),16,1), KIND = wp ) 2158 ibit15 = REAL( IBITS(advc_flags_1(k,j,i),15,1), KIND = wp ) 2159 2160 k_ppp = k + 3 * ibit17 2161 k_pp = k + 2 * ( 1 - ibit15 ) 2162 k_mm = k - 2 * ibit17 2163 2164 w_comp(k) = w(k,j,i) + w(k,j,i-1) 2165 flux_t(k) = w_comp(k) * rho_air_zw(k) * ( & 2166 ( 37.0_wp * ibit17 * adv_mom_5 & 2167 + 7.0_wp * ibit16 * adv_mom_3 & 2168 + ibit15 * adv_mom_1 & 2169 ) * & 2170 ( u(k+1,j,i) + u(k,j,i) ) & 2171 - ( 8.0_wp * ibit17 * adv_mom_5 & 2172 + ibit16 * adv_mom_3 & 2173 ) * & 2174 ( u(k_pp,j,i) + u(k-1,j,i) ) & 2175 + ( ibit17 * adv_mom_5 & 2176 ) * & 2177 ( u(k_ppp,j,i) + u(k_mm,j,i) ) & 2178 ) 2179 2180 diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * ( & 2181 ( 10.0_wp * ibit17 * adv_mom_5 & 2182 + 3.0_wp * ibit16 * adv_mom_3 & 2183 + ibit15 * adv_mom_1 & 2184 ) * & 2185 ( u(k+1,j,i) - u(k,j,i) ) & 2186 - ( 5.0_wp * ibit17 * adv_mom_5 & 2187 + ibit16 * adv_mom_3 & 2188 ) * & 2189 ( u(k_pp,j,i) - u(k-1,j,i) ) & 2190 + ( ibit17 * adv_mom_5 & 2191 ) * & 2192 ( u(k_ppp,j,i) - u(k_mm,j,i) ) & 2193 ) 2194 ENDDO 2195 2196 DO k = nzb+1, nzt 2053 2197 2054 2198 flux_d = flux_t(k-1) … … 2057 2201 !-- Calculate the divergence of the velocity field. A respective 2058 2202 !-- correction is needed to overcome numerical instabilities introduced 2059 !-- by a not sufficient reduction of divergences near topography. 2203 !-- by a not sufficient reduction of divergences near topography. 2060 2204 div = ( ( u_comp(k) * ( ibit9 + ibit10 + ibit11 ) & 2061 2205 - ( u(k,j,i) + u(k,j,i-1) ) & 2062 * ( IBITS(advc_flags_1(k,j,i-1),9,1) & 2063 + IBITS(advc_flags_1(k,j,i-1),10,1) & 2064 + IBITS(advc_flags_1(k,j,i-1),11,1) & 2206 * ( & 2207 REAL( IBITS(advc_flags_1(k,j,i-1),9,1), KIND = wp ) & 2208 + REAL( IBITS(advc_flags_1(k,j,i-1),10,1), KIND = wp ) & 2209 + REAL( IBITS(advc_flags_1(k,j,i-1),11,1), KIND = wp ) & 2065 2210 ) & 2066 2211 ) * ddx & 2067 2212 + ( ( v_comp(k) + gv ) * ( ibit12 + ibit13 + ibit14 ) & 2068 2213 - ( v(k,j,i) + v(k,j,i-1 ) ) & 2069 * ( IBITS(advc_flags_1(k,j-1,i),12,1) & 2070 + IBITS(advc_flags_1(k,j-1,i),13,1) & 2071 + IBITS(advc_flags_1(k,j-1,i),14,1) & 2214 * ( & 2215 REAL( IBITS(advc_flags_1(k,j-1,i),12,1), KIND = wp ) & 2216 + REAL( IBITS(advc_flags_1(k,j-1,i),13,1), KIND = wp ) & 2217 + REAL( IBITS(advc_flags_1(k,j-1,i),14,1), KIND = wp ) & 2072 2218 ) & 2073 2219 ) * ddy & 2074 2220 + ( w_comp(k) * rho_air_zw(k) * ( ibit15 + ibit16 + ibit17 )& 2075 2221 - w_comp(k-1) * rho_air_zw(k-1) & 2076 * ( IBITS(advc_flags_1(k-1,j,i),15,1) & 2077 + IBITS(advc_flags_1(k-1,j,i),16,1) & 2078 + IBITS(advc_flags_1(k-1,j,i),17,1) & 2222 * ( & 2223 REAL( IBITS(advc_flags_1(k-1,j,i),15,1), KIND = wp ) & 2224 + REAL( IBITS(advc_flags_1(k-1,j,i),16,1), KIND = wp ) & 2225 + REAL( IBITS(advc_flags_1(k-1,j,i),17,1), KIND = wp ) & 2079 2226 ) & 2080 2227 ) * drho_air(k) * ddzw(k) & 2081 2228 ) * 0.5_wp 2082 2229 2083 2084 tend(k,j,i) = tend(k,j,i) - ( & 2230 tend(k,j,i) = tend(k,j,i) - ( & 2085 2231 ( flux_r(k) + diss_r(k) & 2086 2232 - flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx & … … 2088 2234 - flux_s_u(k,tn) - diss_s_u(k,tn) ) * ddy & 2089 2235 + ( ( flux_t(k) + diss_t(k) ) & 2090 - ( flux_d + diss_d )&2236 - ( flux_d + diss_d ) & 2091 2237 ) * drho_air(k) * ddzw(k) & 2092 2238 ) + div * u(k,j,i) … … 2111 2257 sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) & 2112 2258 + ( flux_t(k) & 2113 * ( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) &2114 / ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) ) ) &2115 + diss_t(k) &2116 * ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) &2117 / ( ABS( w_comp(k) ) + 1.0E-20_wp ) &2118 ) * weight_substep(intermediate_timestep_count)2119 ENDDO2120 2121 DO k = nzb_max+1, nzt2122 2123 u_comp(k) = u(k,j,i+1) + u(k,j,i)2124 flux_r(k) = ( u_comp(k) - gu ) * ( &2125 37.0_wp * ( u(k,j,i+1) + u(k,j,i) ) &2126 - 8.0_wp * ( u(k,j,i+2) + u(k,j,i-1) ) &2127 + ( u(k,j,i+3) + u(k,j,i-2) ) ) * adv_mom_52128 diss_r(k) = - ABS( u_comp(k) - gu ) * ( &2129 10.0_wp * ( u(k,j,i+1) - u(k,j,i) ) &2130 - 5.0_wp * ( u(k,j,i+2) - u(k,j,i-1) ) &2131 + ( u(k,j,i+3) - u(k,j,i-2) ) ) * adv_mom_52132 2133 v_comp(k) = v(k,j+1,i) + v(k,j+1,i-1) - gv2134 flux_n(k) = v_comp(k) * ( &2135 37.0_wp * ( u(k,j+1,i) + u(k,j,i) ) &2136 - 8.0_wp * ( u(k,j+2,i) + u(k,j-1,i) ) &2137 + ( u(k,j+3,i) + u(k,j-2,i) ) ) * adv_mom_52138 diss_n(k) = - ABS( v_comp(k) ) * ( &2139 10.0_wp * ( u(k,j+1,i) - u(k,j,i) ) &2140 - 5.0_wp * ( u(k,j+2,i) - u(k,j-1,i) ) &2141 + ( u(k,j+3,i) - u(k,j-2,i) ) ) * adv_mom_52142 !2143 !-- k index has to be modified near bottom and top, else array2144 !-- subscripts will be exceeded.2145 ibit17 = IBITS(advc_flags_1(k,j,i),17,1)2146 ibit16 = IBITS(advc_flags_1(k,j,i),16,1)2147 ibit15 = IBITS(advc_flags_1(k,j,i),15,1)2148 2149 k_ppp = k + 3 * ibit172150 k_pp = k + 2 * ( 1 - ibit15 )2151 k_mm = k - 2 * ibit172152 2153 w_comp(k) = w(k,j,i) + w(k,j,i-1)2154 flux_t(k) = w_comp(k) * rho_air_zw(k) * ( &2155 ( 37.0_wp * ibit17 * adv_mom_5 &2156 + 7.0_wp * ibit16 * adv_mom_3 &2157 + ibit15 * adv_mom_1 &2158 ) * &2159 ( u(k+1,j,i) + u(k,j,i) ) &2160 - ( 8.0_wp * ibit17 * adv_mom_5 &2161 + ibit16 * adv_mom_3 &2162 ) * &2163 ( u(k_pp,j,i) + u(k-1,j,i) ) &2164 + ( ibit17 * adv_mom_5 &2165 ) * &2166 ( u(k_ppp,j,i) + u(k_mm,j,i) ) &2167 )2168 2169 diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * ( &2170 ( 10.0_wp * ibit17 * adv_mom_5 &2171 + 3.0_wp * ibit16 * adv_mom_3 &2172 + ibit15 * adv_mom_1 &2173 ) * &2174 ( u(k+1,j,i) - u(k,j,i) ) &2175 - ( 5.0_wp * ibit17 * adv_mom_5 &2176 + ibit16 * adv_mom_3 &2177 ) * &2178 ( u(k_pp,j,i) - u(k-1,j,i) ) &2179 + ( ibit17 * adv_mom_5 &2180 ) * &2181 ( u(k_ppp,j,i) - u(k_mm,j,i) ) &2182 )2183 2184 ENDDO2185 2186 DO k = nzb_max+1, nzt2187 2188 flux_d = flux_t(k-1)2189 diss_d = diss_t(k-1)2190 !2191 !-- Calculate the divergence of the velocity field. A respective2192 !-- correction is needed to overcome numerical instabilities introduced2193 !-- by a not sufficient reduction of divergences near topography.2194 div = ( ( u_comp(k) - ( u(k,j,i) + u(k,j,i-1) ) ) * ddx &2195 + ( v_comp(k) + gv - ( v(k,j,i) + v(k,j,i-1 ) ) ) * ddy &2196 + ( w_comp(k) * rho_air_zw(k) - &2197 w_comp(k-1) * rho_air_zw(k-1) &2198 ) * drho_air(k) * ddzw(k) &2199 ) * 0.5_wp2200 2201 tend(k,j,i) = tend(k,j,i) - ( &2202 ( flux_r(k) + diss_r(k) &2203 - flux_l_u(k,j,tn) - diss_l_u(k,j,tn) ) * ddx &2204 + ( flux_n(k) + diss_n(k) &2205 - flux_s_u(k,tn) - diss_s_u(k,tn) ) * ddy &2206 + ( ( flux_t(k) + diss_t(k) ) &2207 - ( flux_d + diss_d ) &2208 ) * drho_air(k) * ddzw(k) &2209 ) + div * u(k,j,i)2210 2211 flux_l_u(k,j,tn) = flux_r(k)2212 diss_l_u(k,j,tn) = diss_r(k)2213 flux_s_u(k,tn) = flux_n(k)2214 diss_s_u(k,tn) = diss_n(k)2215 !2216 !-- Statistical Evaluation of u'u'. The factor has to be applied for2217 !-- right evaluation when gallilei_trans = .T. .2218 sums_us2_ws_l(k,tn) = sums_us2_ws_l(k,tn) &2219 + ( flux_r(k) &2220 * ( u_comp(k) - 2.0_wp * hom(k,1,1,0) ) &2221 / ( u_comp(k) - gu + SIGN( 1.0E-20_wp, u_comp(k) - gu ) ) &2222 + diss_r(k) &2223 * ABS( u_comp(k) - 2.0_wp * hom(k,1,1,0) ) &2224 / ( ABS( u_comp(k) - gu ) + 1.0E-20_wp ) &2225 ) * weight_substep(intermediate_timestep_count)2226 !2227 !-- Statistical Evaluation of w'u'.2228 sums_wsus_ws_l(k,tn) = sums_wsus_ws_l(k,tn) &2229 + ( flux_t(k) &2230 2259 * ( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) & 2231 2260 / ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) ) ) & … … 2272 2301 2273 2302 INTEGER(iwp) :: i !< grid index along x-direction 2274 INTEGER(iwp) :: ibit18 !< flag indicating 1st-order scheme along x-direction2275 INTEGER(iwp) :: ibit19 !< flag indicating 3rd-order scheme along x-direction2276 INTEGER(iwp) :: ibit20 !< flag indicating 5th-order scheme along x-direction2277 INTEGER(iwp) :: ibit21 !< flag indicating 1st-order scheme along y-direction2278 INTEGER(iwp) :: ibit22 !< flag indicating 3rd-order scheme along y-direction2279 INTEGER(iwp) :: ibit23 !< flag indicating 3rd-order scheme along y-direction2280 INTEGER(iwp) :: ibit24 !< flag indicating 1st-order scheme along z-direction2281 INTEGER(iwp) :: ibit25 !< flag indicating 3rd-order scheme along z-direction2282 INTEGER(iwp) :: ibit26 !< flag indicating 3rd-order scheme along z-direction2283 2303 INTEGER(iwp) :: i_omp !< leftmost index on subdomain, or in case of OpenMP, on thread 2284 2304 INTEGER(iwp) :: j !< grid index along y-direction … … 2289 2309 INTEGER(iwp) :: tn !< number of OpenMP thread 2290 2310 2291 REAL(wp) :: diss_d !< artificial dissipation term at grid box bottom 2292 REAL(wp) :: div !< divergence on v-grid 2293 REAL(wp) :: flux_d !< 6th-order flux at grid box bottom 2294 REAL(wp) :: gu !< Galilei-transformation velocity along x 2295 REAL(wp) :: gv !< Galilei-transformation velocity along y 2296 REAL(wp) :: v_comp_l !< advection velocity along y on leftmost grid point on subdomain 2311 REAL(wp) :: ibit18 !< flag indicating 1st-order scheme along x-direction 2312 REAL(wp) :: ibit19 !< flag indicating 3rd-order scheme along x-direction 2313 REAL(wp) :: ibit20 !< flag indicating 5th-order scheme along x-direction 2314 REAL(wp) :: ibit21 !< flag indicating 1st-order scheme along y-direction 2315 REAL(wp) :: ibit22 !< flag indicating 3rd-order scheme along y-direction 2316 REAL(wp) :: ibit23 !< flag indicating 3rd-order scheme along y-direction 2317 REAL(wp) :: ibit24 !< flag indicating 1st-order scheme along z-direction 2318 REAL(wp) :: ibit25 !< flag indicating 3rd-order scheme along z-direction 2319 REAL(wp) :: ibit26 !< flag indicating 3rd-order scheme along z-direction 2320 REAL(wp) :: diss_d !< artificial dissipation term at grid box bottom 2321 REAL(wp) :: div !< divergence on v-grid 2322 REAL(wp) :: flux_d !< 6th-order flux at grid box bottom 2323 REAL(wp) :: gu !< Galilei-transformation velocity along x 2324 REAL(wp) :: gv !< Galilei-transformation velocity along y 2325 REAL(wp) :: v_comp_l !< advection velocity along y on leftmost grid point on subdomain 2297 2326 2298 2327 REAL(wp), DIMENSION(nzb:nzt+1) :: diss_n !< discretized artificial dissipation at northward-side of the grid box … … 2315 2344 DO k = nzb+1, nzb_max 2316 2345 2317 ibit20 = IBITS(advc_flags_1(k,j,i-1),20,1)2318 ibit19 = IBITS(advc_flags_1(k,j,i-1),19,1)2319 ibit18 = IBITS(advc_flags_1(k,j,i-1),18,1)2346 ibit20 = REAL( IBITS(advc_flags_1(k,j,i-1),20,1), KIND = wp ) 2347 ibit19 = REAL( IBITS(advc_flags_1(k,j,i-1),19,1), KIND = wp ) 2348 ibit18 = REAL( IBITS(advc_flags_1(k,j,i-1),18,1), KIND = wp ) 2320 2349 2321 2350 u_comp(k) = u(k,j-1,i) + u(k,j,i) - gu … … 2373 2402 DO k = nzb+1, nzb_max 2374 2403 2375 ibit23 = IBITS(advc_flags_1(k,j-1,i),23,1)2376 ibit22 = IBITS(advc_flags_1(k,j-1,i),22,1)2377 ibit21 = IBITS(advc_flags_1(k,j-1,i),21,1)2404 ibit23 = REAL( IBITS(advc_flags_1(k,j-1,i),23,1), KIND = wp ) 2405 ibit22 = REAL( IBITS(advc_flags_1(k,j-1,i),22,1), KIND = wp ) 2406 ibit21 = REAL( IBITS(advc_flags_1(k,j-1,i),21,1), KIND = wp ) 2378 2407 2379 2408 v_comp_l = v(k,j,i) + v(k,j-1,i) - gv … … 2425 2454 2426 2455 ENDIF 2427 2428 flux_t(0) = 0.0_wp2429 diss_t(0) = 0.0_wp2430 w_comp(0) = 0.0_wp2431 2456 ! 2432 2457 !-- Now compute the fluxes and tendency terms for the horizontal and … … 2434 2459 DO k = nzb+1, nzb_max 2435 2460 2436 ibit20 = IBITS(advc_flags_1(k,j,i),20,1)2437 ibit19 = IBITS(advc_flags_1(k,j,i),19,1)2438 ibit18 = IBITS(advc_flags_1(k,j,i),18,1)2461 ibit20 = REAL( IBITS(advc_flags_1(k,j,i),20,1), KIND = wp ) 2462 ibit19 = REAL( IBITS(advc_flags_1(k,j,i),19,1), KIND = wp ) 2463 ibit18 = REAL( IBITS(advc_flags_1(k,j,i),18,1), KIND = wp ) 2439 2464 2440 2465 u_comp(k) = u(k,j-1,i+1) + u(k,j,i+1) - gu … … 2469 2494 ) 2470 2495 2471 ibit23 = IBITS(advc_flags_1(k,j,i),23,1)2472 ibit22 = IBITS(advc_flags_1(k,j,i),22,1)2473 ibit21 = IBITS(advc_flags_1(k,j,i),21,1)2496 ibit23 = REAL( IBITS(advc_flags_1(k,j,i),23,1), KIND = wp ) 2497 ibit22 = REAL( IBITS(advc_flags_1(k,j,i),22,1), KIND = wp ) 2498 ibit21 = REAL( IBITS(advc_flags_1(k,j,i),21,1), KIND = wp ) 2474 2499 2475 2500 … … 2504 2529 ( v(k,j+3,i) - v(k,j-2,i) ) & 2505 2530 ) 2531 ENDDO 2532 2533 DO k = nzb_max+1, nzt 2534 2535 u_comp(k) = u(k,j-1,i+1) + u(k,j,i+1) - gu 2536 flux_r(k) = u_comp(k) * ( & 2537 37.0_wp * ( v(k,j,i+1) + v(k,j,i) ) & 2538 - 8.0_wp * ( v(k,j,i+2) + v(k,j,i-1) ) & 2539 + ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_5 2540 2541 diss_r(k) = - ABS( u_comp(k) ) * ( & 2542 10.0_wp * ( v(k,j,i+1) - v(k,j,i) ) & 2543 - 5.0_wp * ( v(k,j,i+2) - v(k,j,i-1) ) & 2544 + ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_5 2545 2546 2547 v_comp(k) = v(k,j+1,i) + v(k,j,i) 2548 flux_n(k) = ( v_comp(k) - gv ) * ( & 2549 37.0_wp * ( v(k,j+1,i) + v(k,j,i) ) & 2550 - 8.0_wp * ( v(k,j+2,i) + v(k,j-1,i) ) & 2551 + ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_5 2552 2553 diss_n(k) = - ABS( v_comp(k) - gv ) * ( & 2554 10.0_wp * ( v(k,j+1,i) - v(k,j,i) ) & 2555 - 5.0_wp * ( v(k,j+2,i) - v(k,j-1,i) ) & 2556 + ( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_5 2557 ENDDO 2558 ! 2559 !-- Now, compute vertical fluxes. Split loop into a part treating the 2560 !-- lowest 2 grid points with indirect indexing, a main loop without 2561 !-- indirect indexing, and a loop for the uppermost 2 grip points with 2562 !-- indirect indexing. This allows better vectorization for the main loop. 2563 !-- First, compute the flux at model surface, which need has to be 2564 !-- calculated explicetely for the tendency at 2565 !-- the first w-level. For topography wall this is done implicitely by 2566 !-- advc_flags_1. 2567 flux_t(nzb) = 0.0_wp 2568 diss_t(nzb) = 0.0_wp 2569 w_comp(nzb) = 0.0_wp 2570 2571 DO k = nzb+1, nzb+2 2506 2572 ! 2507 2573 !-- k index has to be modified near bottom and top, else array 2508 2574 !-- subscripts will be exceeded. 2509 ibit26 = IBITS(advc_flags_1(k,j,i),26,1)2510 ibit25 = IBITS(advc_flags_1(k,j,i),25,1)2511 ibit24 = IBITS(advc_flags_1(k,j,i),24,1)2575 ibit26 = REAL( IBITS(advc_flags_1(k,j,i),26,1), KIND = wp ) 2576 ibit25 = REAL( IBITS(advc_flags_1(k,j,i),25,1), KIND = wp ) 2577 ibit24 = REAL( IBITS(advc_flags_1(k,j,i),24,1), KIND = wp ) 2512 2578 2513 2579 k_ppp = k + 3 * ibit26 … … 2546 2612 ) 2547 2613 ENDDO 2548 2549 DO k = nzb+1, nzb_max 2614 2615 DO k = nzb+3, nzt-2 2616 2617 ibit26 = REAL( IBITS(advc_flags_1(k,j,i),26,1), KIND = wp ) 2618 ibit25 = REAL( IBITS(advc_flags_1(k,j,i),25,1), KIND = wp ) 2619 ibit24 = REAL( IBITS(advc_flags_1(k,j,i),24,1), KIND = wp ) 2620 2621 w_comp(k) = w(k,j-1,i) + w(k,j,i) 2622 flux_t(k) = w_comp(k) * rho_air_zw(k) * ( & 2623 ( 37.0_wp * ibit26 * adv_mom_5 & 2624 + 7.0_wp * ibit25 * adv_mom_3 & 2625 + ibit24 * adv_mom_1 & 2626 ) * & 2627 ( v(k+1,j,i) + v(k,j,i) ) & 2628 - ( 8.0_wp * ibit26 * adv_mom_5 & 2629 + ibit25 * adv_mom_3 & 2630 ) * & 2631 ( v(k+2,j,i) + v(k-1,j,i) ) & 2632 + ( ibit26 * adv_mom_5 & 2633 ) * & 2634 ( v(k+3,j,i) + v(k-2,j,i) ) & 2635 ) 2636 2637 diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * ( & 2638 ( 10.0_wp * ibit26 * adv_mom_5 & 2639 + 3.0_wp * ibit25 * adv_mom_3 & 2640 + ibit24 * adv_mom_1 & 2641 ) * & 2642 ( v(k+1,j,i) - v(k,j,i) ) & 2643 - ( 5.0_wp * ibit26 * adv_mom_5 & 2644 + ibit25 * adv_mom_3 & 2645 ) * & 2646 ( v(k+2,j,i) - v(k-1,j,i) ) & 2647 + ( ibit26 * adv_mom_5 & 2648 ) * & 2649 ( v(k+3,j,i) - v(k-2,j,i) ) & 2650 ) 2651 ENDDO 2652 2653 DO k = nzt-1, nzt 2654 ! 2655 !-- k index has to be modified near bottom and top, else array 2656 !-- subscripts will be exceeded. 2657 ibit26 = REAL( IBITS(advc_flags_1(k,j,i),26,1), KIND = wp ) 2658 ibit25 = REAL( IBITS(advc_flags_1(k,j,i),25,1), KIND = wp ) 2659 ibit24 = REAL( IBITS(advc_flags_1(k,j,i),24,1), KIND = wp ) 2660 2661 k_ppp = k + 3 * ibit26 2662 k_pp = k + 2 * ( 1 - ibit24 ) 2663 k_mm = k - 2 * ibit26 2664 2665 w_comp(k) = w(k,j-1,i) + w(k,j,i) 2666 flux_t(k) = w_comp(k) * rho_air_zw(k) * ( & 2667 ( 37.0_wp * ibit26 * adv_mom_5 & 2668 + 7.0_wp * ibit25 * adv_mom_3 & 2669 + ibit24 * adv_mom_1 & 2670 ) * & 2671 ( v(k+1,j,i) + v(k,j,i) ) & 2672 - ( 8.0_wp * ibit26 * adv_mom_5 & 2673 + ibit25 * adv_mom_3 & 2674 ) * & 2675 ( v(k_pp,j,i) + v(k-1,j,i) ) & 2676 + ( ibit26 * adv_mom_5 & 2677 ) * & 2678 ( v(k_ppp,j,i) + v(k_mm,j,i) ) & 2679 ) 2680 2681 diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * ( & 2682 ( 10.0_wp * ibit26 * adv_mom_5 & 2683 + 3.0_wp * ibit25 * adv_mom_3 & 2684 + ibit24 * adv_mom_1 & 2685 ) * & 2686 ( v(k+1,j,i) - v(k,j,i) ) & 2687 - ( 5.0_wp * ibit26 * adv_mom_5 & 2688 + ibit25 * adv_mom_3 & 2689 ) * & 2690 ( v(k_pp,j,i) - v(k-1,j,i) ) & 2691 + ( ibit26 * adv_mom_5 & 2692 ) * & 2693 ( v(k_ppp,j,i) - v(k_mm,j,i) ) & 2694 ) 2695 ENDDO 2696 2697 DO k = nzb+1, nzt 2550 2698 2551 2699 flux_d = flux_t(k-1) … … 2554 2702 !-- Calculate the divergence of the velocity field. A respective 2555 2703 !-- correction is needed to overcome numerical instabilities introduced 2556 !-- by a not sufficient reduction of divergences near topography. 2704 !-- by a not sufficient reduction of divergences near topography. 2557 2705 div = ( ( ( u_comp(k) + gu ) & 2558 2706 * ( ibit18 + ibit19 + ibit20 ) & 2559 2707 - ( u(k,j-1,i) + u(k,j,i) ) & 2560 * ( IBITS(advc_flags_1(k,j,i-1),18,1) & 2561 + IBITS(advc_flags_1(k,j,i-1),19,1) & 2562 + IBITS(advc_flags_1(k,j,i-1),20,1) & 2708 * ( & 2709 REAL( IBITS(advc_flags_1(k,j,i-1),18,1), KIND = wp ) & 2710 + REAL( IBITS(advc_flags_1(k,j,i-1),19,1), KIND = wp ) & 2711 + REAL( IBITS(advc_flags_1(k,j,i-1),20,1), KIND = wp ) & 2563 2712 ) & 2564 2713 ) * ddx & … … 2566 2715 * ( ibit21 + ibit22 + ibit23 ) & 2567 2716 - ( v(k,j,i) + v(k,j-1,i) ) & 2568 * ( IBITS(advc_flags_1(k,j-1,i),21,1) & 2569 + IBITS(advc_flags_1(k,j-1,i),22,1) & 2570 + IBITS(advc_flags_1(k,j-1,i),23,1) & 2717 * ( & 2718 REAL( IBITS(advc_flags_1(k,j-1,i),21,1), KIND = wp ) & 2719 + REAL( IBITS(advc_flags_1(k,j-1,i),22,1), KIND = wp ) & 2720 + REAL( IBITS(advc_flags_1(k,j-1,i),23,1), KIND = wp ) & 2571 2721 ) & 2572 2722 ) * ddy & 2573 2723 + ( w_comp(k) * rho_air_zw(k) * ( ibit24 + ibit25 + ibit26 )& 2574 2724 - w_comp(k-1) * rho_air_zw(k-1) & 2575 * ( IBITS(advc_flags_1(k-1,j,i),24,1) & 2576 + IBITS(advc_flags_1(k-1,j,i),25,1) & 2577 + IBITS(advc_flags_1(k-1,j,i),26,1) & 2725 * ( & 2726 REAL( IBITS(advc_flags_1(k-1,j,i),24,1), KIND = wp ) & 2727 + REAL( IBITS(advc_flags_1(k-1,j,i),25,1), KIND = wp ) & 2728 + REAL( IBITS(advc_flags_1(k-1,j,i),26,1), KIND = wp ) & 2578 2729 ) & 2579 2730 ) * drho_air(k) * ddzw(k) & 2580 2731 ) * 0.5_wp 2581 2582 2732 2583 2733 tend(k,j,i) = tend(k,j,i) - ( & … … 2595 2745 flux_s_v(k,tn) = flux_n(k) 2596 2746 diss_s_v(k,tn) = diss_n(k) 2597 2598 2747 ! 2599 2748 !-- Statistical Evaluation of v'v'. The factor has to be applied for … … 2619 2768 2620 2769 ENDDO 2621 2622 DO k = nzb_max+1, nzt2623 2624 u_comp(k) = u(k,j-1,i+1) + u(k,j,i+1) - gu2625 flux_r(k) = u_comp(k) * ( &2626 37.0_wp * ( v(k,j,i+1) + v(k,j,i) ) &2627 - 8.0_wp * ( v(k,j,i+2) + v(k,j,i-1) ) &2628 + ( v(k,j,i+3) + v(k,j,i-2) ) ) * adv_mom_52629 2630 diss_r(k) = - ABS( u_comp(k) ) * ( &2631 10.0_wp * ( v(k,j,i+1) - v(k,j,i) ) &2632 - 5.0_wp * ( v(k,j,i+2) - v(k,j,i-1) ) &2633 + ( v(k,j,i+3) - v(k,j,i-2) ) ) * adv_mom_52634 2635 2636 v_comp(k) = v(k,j+1,i) + v(k,j,i)2637 flux_n(k) = ( v_comp(k) - gv ) * ( &2638 37.0_wp * ( v(k,j+1,i) + v(k,j,i) ) &2639 - 8.0_wp * ( v(k,j+2,i) + v(k,j-1,i) ) &2640 + ( v(k,j+3,i) + v(k,j-2,i) ) ) * adv_mom_52641 2642 diss_n(k) = - ABS( v_comp(k) - gv ) * ( &2643 10.0_wp * ( v(k,j+1,i) - v(k,j,i) ) &2644 - 5.0_wp * ( v(k,j+2,i) - v(k,j-1,i) ) &2645 + ( v(k,j+3,i) - v(k,j-2,i) ) ) * adv_mom_52646 !2647 !-- k index has to be modified near bottom and top, else array2648 !-- subscripts will be exceeded.2649 ibit26 = IBITS(advc_flags_1(k,j,i),26,1)2650 ibit25 = IBITS(advc_flags_1(k,j,i),25,1)2651 ibit24 = IBITS(advc_flags_1(k,j,i),24,1)2652 2653 k_ppp = k + 3 * ibit262654 k_pp = k + 2 * ( 1 - ibit24 )2655 k_mm = k - 2 * ibit262656 2657 w_comp(k) = w(k,j-1,i) + w(k,j,i)2658 flux_t(k) = w_comp(k) * rho_air_zw(k) * ( &2659 ( 37.0_wp * ibit26 * adv_mom_5 &2660 + 7.0_wp * ibit25 * adv_mom_3 &2661 + ibit24 * adv_mom_1 &2662 ) * &2663 ( v(k+1,j,i) + v(k,j,i) ) &2664 - ( 8.0_wp * ibit26 * adv_mom_5 &2665 + ibit25 * adv_mom_3 &2666 ) * &2667 ( v(k_pp,j,i) + v(k-1,j,i) ) &2668 + ( ibit26 * adv_mom_5 &2669 ) * &2670 ( v(k_ppp,j,i) + v(k_mm,j,i) ) &2671 )2672 2673 diss_t(k) = - ABS( w_comp(k) ) * rho_air_zw(k) * ( &2674 ( 10.0_wp * ibit26 * adv_mom_5 &2675 + 3.0_wp * ibit25 * adv_mom_3 &2676 + ibit24 * adv_mom_1 &2677 ) * &2678 ( v(k+1,j,i) - v(k,j,i) ) &2679 - ( 5.0_wp * ibit26 * adv_mom_5 &2680 + ibit25 * adv_mom_3 &2681 ) * &2682 ( v(k_pp,j,i) - v(k-1,j,i) ) &2683 + ( ibit26 * adv_mom_5 &2684 ) * &2685 ( v(k_ppp,j,i) - v(k_mm,j,i) ) &2686 )2687 ENDDO2688 2689 DO k = nzb_max+1, nzt2690 2691 flux_d = flux_t(k-1)2692 diss_d = diss_t(k-1)2693 !2694 !-- Calculate the divergence of the velocity field. A respective2695 !-- correction is needed to overcome numerical instabilities introduced2696 !-- by a not sufficient reduction of divergences near topography.2697 div = ( ( u_comp(k) + gu - ( u(k,j-1,i) + u(k,j,i) ) ) * ddx &2698 + ( v_comp(k) - ( v(k,j,i) + v(k,j-1,i) ) ) * ddy &2699 + ( w_comp(k) * rho_air_zw(k) - &2700 w_comp(k-1) * rho_air_zw(k-1) &2701 ) * drho_air(k) * ddzw(k) &2702 ) * 0.5_wp2703 2704 tend(k,j,i) = tend(k,j,i) - ( &2705 ( flux_r(k) + diss_r(k) &2706 - flux_l_v(k,j,tn) - diss_l_v(k,j,tn) ) * ddx &2707 + ( flux_n(k) + diss_n(k) &2708 - flux_s_v(k,tn) - diss_s_v(k,tn) ) * ddy &2709 + ( ( flux_t(k) + diss_t(k) ) &2710 - ( flux_d + diss_d ) &2711 ) * drho_air(k) * ddzw(k) &2712 ) + v(k,j,i) * div2713 2714 flux_l_v(k,j,tn) = flux_r(k)2715 diss_l_v(k,j,tn) = diss_r(k)2716 flux_s_v(k,tn) = flux_n(k)2717 diss_s_v(k,tn) = diss_n(k)2718 !2719 !-- Statistical Evaluation of v'v'. The factor has to be applied for2720 !-- right evaluation when gallilei_trans = .T. .2721 sums_vs2_ws_l(k,tn) = sums_vs2_ws_l(k,tn) &2722 + ( flux_n(k) &2723 * ( v_comp(k) - 2.0_wp * hom(k,1,2,0) ) &2724 / ( v_comp(k) - gv + SIGN( 1.0E-20_wp, v_comp(k) - gv ) ) &2725 + diss_n(k) &2726 * ABS( v_comp(k) - 2.0_wp * hom(k,1,2,0) ) &2727 / ( ABS( v_comp(k) - gv ) + 1.0E-20_wp ) &2728 ) * weight_substep(intermediate_timestep_count)2729 !2730 !-- Statistical Evaluation of w'u'.2731 sums_wsvs_ws_l(k,tn) = sums_wsvs_ws_l(k,tn) &2732 + ( flux_t(k) &2733 * ( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) &2734 / ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) ) ) &2735 + diss_t(k) &2736 * ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) &2737 / ( ABS( w_comp(k) ) + 1.0E-20_wp ) &2738 ) * weight_substep(intermediate_timestep_count)2739 2740 ENDDO2741 2770 sums_vs2_ws_l(nzb,tn) = sums_vs2_ws_l(nzb+1,tn) 2742 2771 … … 2774 2803 2775 2804 INTEGER(iwp) :: i !< grid index along x-direction 2776 INTEGER(iwp) :: ibit27 !< flag indicating 1st-order scheme along x-direction2777 INTEGER(iwp) :: ibit28 !< flag indicating 3rd-order scheme along x-direction2778 INTEGER(iwp) :: ibit29 !< flag indicating 5th-order scheme along x-direction2779 INTEGER(iwp) :: ibit30 !< flag indicating 1st-order scheme along y-direction2780 INTEGER(iwp) :: ibit31 !< flag indicating 3rd-order scheme along y-direction2781 INTEGER(iwp) :: ibit32 !< flag indicating 5th-order scheme along y-direction2782 INTEGER(iwp) :: ibit33 !< flag indicating 1st-order scheme along z-direction2783 INTEGER(iwp) :: ibit34 !< flag indicating 3rd-order scheme along z-direction2784 INTEGER(iwp) :: ibit35 !< flag indicating 5th-order scheme along z-direction2785 2805 INTEGER(iwp) :: i_omp !< leftmost index on subdomain, or in case of OpenMP, on thread 2786 2806 INTEGER(iwp) :: j !< grid index along y-direction … … 2791 2811 INTEGER(iwp) :: tn !< number of OpenMP thread 2792 2812 2813 REAL(wp) :: ibit27 !< flag indicating 1st-order scheme along x-direction 2814 REAL(wp) :: ibit28 !< flag indicating 3rd-order scheme along x-direction 2815 REAL(wp) :: ibit29 !< flag indicating 5th-order scheme along x-direction 2816 REAL(wp) :: ibit30 !< flag indicating 1st-order scheme along y-direction 2817 REAL(wp) :: ibit31 !< flag indicating 3rd-order scheme along y-direction 2818 REAL(wp) :: ibit32 !< flag indicating 5th-order scheme along y-direction 2819 REAL(wp) :: ibit33 !< flag indicating 1st-order scheme along z-direction 2820 REAL(wp) :: ibit34 !< flag indicating 3rd-order scheme along z-direction 2821 REAL(wp) :: ibit35 !< flag indicating 5th-order scheme along z-direction 2793 2822 REAL(wp) :: diss_d !< discretized artificial dissipation at top of the grid box 2794 2823 REAL(wp) :: div !< divergence on w-grid … … 2815 2844 2816 2845 DO k = nzb+1, nzb_max 2817 ibit32 = IBITS(advc_flags_2(k,j-1,i),0,1)2818 ibit31 = IBITS(advc_flags_1(k,j-1,i),31,1)2819 ibit30 = IBITS(advc_flags_1(k,j-1,i),30,1)2846 ibit32 = REAL( IBITS(advc_flags_2(k,j-1,i),0,1), KIND = wp ) 2847 ibit31 = REAL( IBITS(advc_flags_1(k,j-1,i),31,1), KIND = wp ) 2848 ibit30 = REAL( IBITS(advc_flags_1(k,j-1,i),30,1), KIND = wp ) 2820 2849 2821 2850 v_comp(k) = v(k+1,j,i) + v(k,j,i) - gv … … 2859 2888 - 8.0_wp * ( w(k,j+1,i) +w(k,j-2,i) ) & 2860 2889 + ( w(k,j+2,i) + w(k,j-3,i) ) ) * adv_mom_5 2861 diss_s_w(k,tn) = - ABS( v_comp(k) ) * ( &2890 diss_s_w(k,tn) = - ABS( v_comp(k) ) * ( & 2862 2891 10.0_wp * ( w(k,j,i) - w(k,j-1,i) ) & 2863 2892 - 5.0_wp * ( w(k,j+1,i) - w(k,j-2,i) ) & … … 2873 2902 DO k = nzb+1, nzb_max 2874 2903 2875 ibit29 = IBITS(advc_flags_1(k,j,i-1),29,1)2876 ibit28 = IBITS(advc_flags_1(k,j,i-1),28,1)2877 ibit27 = IBITS(advc_flags_1(k,j,i-1),27,1)2904 ibit29 = REAL( IBITS(advc_flags_1(k,j,i-1),29,1), KIND = wp ) 2905 ibit28 = REAL( IBITS(advc_flags_1(k,j,i-1),28,1), KIND = wp ) 2906 ibit27 = REAL( IBITS(advc_flags_1(k,j,i-1),27,1), KIND = wp ) 2878 2907 2879 2908 u_comp(k) = u(k+1,j,i) + u(k,j,i) - gu … … 2926 2955 ENDIF 2927 2956 ! 2928 !-- The lower flux has to be calculated explicetely for the tendency at 2957 !-- Now compute the fluxes and tendency terms for the horizontal 2958 !-- and vertical parts. 2959 DO k = nzb+1, nzb_max 2960 2961 ibit29 = REAL( IBITS(advc_flags_1(k,j,i),29,1), KIND = wp ) 2962 ibit28 = REAL( IBITS(advc_flags_1(k,j,i),28,1), KIND = wp ) 2963 ibit27 = REAL( IBITS(advc_flags_1(k,j,i),27,1), KIND = wp ) 2964 2965 u_comp(k) = u(k+1,j,i+1) + u(k,j,i+1) - gu 2966 flux_r(k) = u_comp(k) * ( & 2967 ( 37.0_wp * ibit29 * adv_mom_5 & 2968 + 7.0_wp * ibit28 * adv_mom_3 & 2969 + ibit27 * adv_mom_1 & 2970 ) * & 2971 ( w(k,j,i+1) + w(k,j,i) ) & 2972 - ( 8.0_wp * ibit29 * adv_mom_5 & 2973 + ibit28 * adv_mom_3 & 2974 ) * & 2975 ( w(k,j,i+2) + w(k,j,i-1) ) & 2976 + ( ibit29 * adv_mom_5 & 2977 ) * & 2978 ( w(k,j,i+3) + w(k,j,i-2) ) & 2979 ) 2980 2981 diss_r(k) = - ABS( u_comp(k) ) * ( & 2982 ( 10.0_wp * ibit29 * adv_mom_5 & 2983 + 3.0_wp * ibit28 * adv_mom_3 & 2984 + ibit27 * adv_mom_1 & 2985 ) * & 2986 ( w(k,j,i+1) - w(k,j,i) ) & 2987 - ( 5.0_wp * ibit29 * adv_mom_5 & 2988 + ibit28 * adv_mom_3 & 2989 ) * & 2990 ( w(k,j,i+2) - w(k,j,i-1) ) & 2991 + ( ibit29 * adv_mom_5 & 2992 ) * & 2993 ( w(k,j,i+3) - w(k,j,i-2) ) & 2994 ) 2995 2996 ibit32 = REAL( IBITS(advc_flags_2(k,j,i),0,1), KIND = wp ) 2997 ibit31 = REAL( IBITS(advc_flags_1(k,j,i),31,1), KIND = wp ) 2998 ibit30 = REAL( IBITS(advc_flags_1(k,j,i),30,1), KIND = wp ) 2999 3000 v_comp(k) = v(k+1,j+1,i) + v(k,j+1,i) - gv 3001 flux_n(k) = v_comp(k) * ( & 3002 ( 37.0_wp * ibit32 * adv_mom_5 & 3003 + 7.0_wp * ibit31 * adv_mom_3 & 3004 + ibit30 * adv_mom_1 & 3005 ) * & 3006 ( w(k,j+1,i) + w(k,j,i) ) & 3007 - ( 8.0_wp * ibit32 * adv_mom_5 & 3008 + ibit31 * adv_mom_3 & 3009 ) * & 3010 ( w(k,j+2,i) + w(k,j-1,i) ) & 3011 + ( ibit32 * adv_mom_5 & 3012 ) * & 3013 ( w(k,j+3,i) + w(k,j-2,i) ) & 3014 ) 3015 3016 diss_n(k) = - ABS( v_comp(k) ) * ( & 3017 ( 10.0_wp * ibit32 * adv_mom_5 & 3018 + 3.0_wp * ibit31 * adv_mom_3 & 3019 + ibit30 * adv_mom_1 & 3020 ) * & 3021 ( w(k,j+1,i) - w(k,j,i) ) & 3022 - ( 5.0_wp * ibit32 * adv_mom_5 & 3023 + ibit31 * adv_mom_3 & 3024 ) * & 3025 ( w(k,j+2,i) - w(k,j-1,i) ) & 3026 + ( ibit32 * adv_mom_5 & 3027 ) * & 3028 ( w(k,j+3,i) - w(k,j-2,i) ) & 3029 ) 3030 ENDDO 3031 3032 DO k = nzb_max+1, nzt 3033 3034 u_comp(k) = u(k+1,j,i+1) + u(k,j,i+1) - gu 3035 flux_r(k) = u_comp(k) * ( & 3036 37.0_wp * ( w(k,j,i+1) + w(k,j,i) ) & 3037 - 8.0_wp * ( w(k,j,i+2) + w(k,j,i-1) ) & 3038 + ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_5 3039 3040 diss_r(k) = - ABS( u_comp(k) ) * ( & 3041 10.0_wp * ( w(k,j,i+1) - w(k,j,i) ) & 3042 - 5.0_wp * ( w(k,j,i+2) - w(k,j,i-1) ) & 3043 + ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_5 3044 3045 v_comp(k) = v(k+1,j+1,i) + v(k,j+1,i) - gv 3046 flux_n(k) = v_comp(k) * ( & 3047 37.0_wp * ( w(k,j+1,i) + w(k,j,i) ) & 3048 - 8.0_wp * ( w(k,j+2,i) + w(k,j-1,i) ) & 3049 + ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_5 3050 3051 diss_n(k) = - ABS( v_comp(k) ) * ( & 3052 10.0_wp * ( w(k,j+1,i) - w(k,j,i) ) & 3053 - 5.0_wp * ( w(k,j+2,i) - w(k,j-1,i) ) & 3054 + ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_5 3055 ENDDO 3056 3057 ! 3058 !-- Now, compute vertical fluxes. Split loop into a part treating the 3059 !-- lowest 2 grid points with indirect indexing, a main loop without 3060 !-- indirect indexing, and a loop for the uppermost 2 grip points with 3061 !-- indirect indexing. This allows better vectorization for the main loop. 3062 !-- First, compute the flux at model surface, which need has to be 3063 !-- calculated explicetely for the tendency at 2929 3064 !-- the first w-level. For topography wall this is done implicitely by 2930 3065 !-- advc_flags_1. … … 2933 3068 flux_t(0) = w_comp(k) * ( w(k,j,i) + w(k-1,j,i) ) * adv_mom_1 2934 3069 diss_t(0) = -ABS(w_comp(k)) * ( w(k,j,i) - w(k-1,j,i) ) * adv_mom_1 2935 ! 2936 !-- Now compute the fluxes and tendency terms for the horizontal 2937 !-- and vertical parts. 2938 DO k = nzb+1, nzb_max 2939 2940 ibit29 = IBITS(advc_flags_1(k,j,i),29,1) 2941 ibit28 = IBITS(advc_flags_1(k,j,i),28,1) 2942 ibit27 = IBITS(advc_flags_1(k,j,i),27,1) 2943 2944 u_comp(k) = u(k+1,j,i+1) + u(k,j,i+1) - gu 2945 flux_r(k) = u_comp(k) * ( & 2946 ( 37.0_wp * ibit29 * adv_mom_5 & 2947 + 7.0_wp * ibit28 * adv_mom_3 & 2948 + ibit27 * adv_mom_1 & 2949 ) * & 2950 ( w(k,j,i+1) + w(k,j,i) ) & 2951 - ( 8.0_wp * ibit29 * adv_mom_5 & 2952 + ibit28 * adv_mom_3 & 2953 ) * & 2954 ( w(k,j,i+2) + w(k,j,i-1) ) & 2955 + ( ibit29 * adv_mom_5 & 2956 ) * & 2957 ( w(k,j,i+3) + w(k,j,i-2) ) & 2958 ) 2959 2960 diss_r(k) = - ABS( u_comp(k) ) * ( & 2961 ( 10.0_wp * ibit29 * adv_mom_5 & 2962 + 3.0_wp * ibit28 * adv_mom_3 & 2963 + ibit27 * adv_mom_1 & 2964 ) * & 2965 ( w(k,j,i+1) - w(k,j,i) ) & 2966 - ( 5.0_wp * ibit29 * adv_mom_5 & 2967 + ibit28 * adv_mom_3 & 2968 ) * & 2969 ( w(k,j,i+2) - w(k,j,i-1) ) & 2970 + ( ibit29 * adv_mom_5 & 2971 ) * & 2972 ( w(k,j,i+3) - w(k,j,i-2) ) & 2973 ) 2974 2975 ibit32 = IBITS(advc_flags_2(k,j,i),0,1) 2976 ibit31 = IBITS(advc_flags_1(k,j,i),31,1) 2977 ibit30 = IBITS(advc_flags_1(k,j,i),30,1) 2978 2979 v_comp(k) = v(k+1,j+1,i) + v(k,j+1,i) - gv 2980 flux_n(k) = v_comp(k) * ( & 2981 ( 37.0_wp * ibit32 * adv_mom_5 & 2982 + 7.0_wp * ibit31 * adv_mom_3 & 2983 + ibit30 * adv_mom_1 & 2984 ) * & 2985 ( w(k,j+1,i) + w(k,j,i) ) & 2986 - ( 8.0_wp * ibit32 * adv_mom_5 & 2987 + ibit31 * adv_mom_3 & 2988 ) * & 2989 ( w(k,j+2,i) + w(k,j-1,i) ) & 2990 + ( ibit32 * adv_mom_5 & 2991 ) * & 2992 ( w(k,j+3,i) + w(k,j-2,i) ) & 2993 ) 2994 2995 diss_n(k) = - ABS( v_comp(k) ) * ( & 2996 ( 10.0_wp * ibit32 * adv_mom_5 & 2997 + 3.0_wp * ibit31 * adv_mom_3 & 2998 + ibit30 * adv_mom_1 & 2999 ) * & 3000 ( w(k,j+1,i) - w(k,j,i) ) & 3001 - ( 5.0_wp * ibit32 * adv_mom_5 & 3002 + ibit31 * adv_mom_3 & 3003 ) * & 3004 ( w(k,j+2,i) - w(k,j-1,i) ) & 3005 + ( ibit32 * adv_mom_5 & 3006 ) * & 3007 ( w(k,j+3,i) - w(k,j-2,i) ) & 3008 ) 3070 3071 DO k = nzb+1, nzb+2 3009 3072 ! 3010 3073 !-- k index has to be modified near bottom and top, else array 3011 3074 !-- subscripts will be exceeded. 3012 ibit35 = IBITS(advc_flags_2(k,j,i),3,1)3013 ibit34 = IBITS(advc_flags_2(k,j,i),2,1)3014 ibit33 = IBITS(advc_flags_2(k,j,i),1,1)3075 ibit35 = REAL( IBITS(advc_flags_2(k,j,i),3,1), KIND = wp ) 3076 ibit34 = REAL( IBITS(advc_flags_2(k,j,i),2,1), KIND = wp ) 3077 ibit33 = REAL( IBITS(advc_flags_2(k,j,i),1,1), KIND = wp ) 3015 3078 3016 3079 k_ppp = k + 3 * ibit35 … … 3049 3112 ) 3050 3113 ENDDO 3051 3052 DO k = nzb+1, nzb_max 3114 3115 DO k = nzb+3, nzt-2 3116 3117 ibit35 = REAL( IBITS(advc_flags_2(k,j,i),3,1), KIND = wp ) 3118 ibit34 = REAL( IBITS(advc_flags_2(k,j,i),2,1), KIND = wp ) 3119 ibit33 = REAL( IBITS(advc_flags_2(k,j,i),1,1), KIND = wp ) 3120 3121 w_comp(k) = w(k+1,j,i) + w(k,j,i) 3122 flux_t(k) = w_comp(k) * rho_air(k+1) * ( & 3123 ( 37.0_wp * ibit35 * adv_mom_5 & 3124 + 7.0_wp * ibit34 * adv_mom_3 & 3125 + ibit33 * adv_mom_1 & 3126 ) * & 3127 ( w(k+1,j,i) + w(k,j,i) ) & 3128 - ( 8.0_wp * ibit35 * adv_mom_5 & 3129 + ibit34 * adv_mom_3 & 3130 ) * & 3131 ( w(k+2,j,i) + w(k-1,j,i) ) & 3132 + ( ibit35 * adv_mom_5 & 3133 ) * & 3134 ( w(k+3,j,i) + w(k-2,j,i) ) & 3135 ) 3136 3137 diss_t(k) = - ABS( w_comp(k) ) * rho_air(k+1) * ( & 3138 ( 10.0_wp * ibit35 * adv_mom_5 & 3139 + 3.0_wp * ibit34 * adv_mom_3 & 3140 + ibit33 * adv_mom_1 & 3141 ) * & 3142 ( w(k+1,j,i) - w(k,j,i) ) & 3143 - ( 5.0_wp * ibit35 * adv_mom_5 & 3144 + ibit34 * adv_mom_3 & 3145 ) * & 3146 ( w(k+2,j,i) - w(k-1,j,i) ) & 3147 + ( ibit35 * adv_mom_5 & 3148 ) * & 3149 ( w(k+3,j,i) - w(k-2,j,i) ) & 3150 ) 3151 ENDDO 3152 3153 DO k = nzt-1, nzt 3154 ! 3155 !-- k index has to be modified near bottom and top, else array 3156 !-- subscripts will be exceeded. 3157 ibit35 = REAL( IBITS(advc_flags_2(k,j,i),3,1), KIND = wp ) 3158 ibit34 = REAL( IBITS(advc_flags_2(k,j,i),2,1), KIND = wp ) 3159 ibit33 = REAL( IBITS(advc_flags_2(k,j,i),1,1), KIND = wp ) 3160 3161 k_ppp = k + 3 * ibit35 3162 k_pp = k + 2 * ( 1 - ibit33 ) 3163 k_mm = k - 2 * ibit35 3164 3165 w_comp(k) = w(k+1,j,i) + w(k,j,i) 3166 flux_t(k) = w_comp(k) * rho_air(k+1) * ( & 3167 ( 37.0_wp * ibit35 * adv_mom_5 & 3168 + 7.0_wp * ibit34 * adv_mom_3 & 3169 + ibit33 * adv_mom_1 & 3170 ) * & 3171 ( w(k+1,j,i) + w(k,j,i) ) & 3172 - ( 8.0_wp * ibit35 * adv_mom_5 & 3173 + ibit34 * adv_mom_3 & 3174 ) * & 3175 ( w(k_pp,j,i) + w(k-1,j,i) ) & 3176 + ( ibit35 * adv_mom_5 & 3177 ) * & 3178 ( w(k_ppp,j,i) + w(k_mm,j,i) ) & 3179 ) 3180 3181 diss_t(k) = - ABS( w_comp(k) ) * rho_air(k+1) * ( & 3182 ( 10.0_wp * ibit35 * adv_mom_5 & 3183 + 3.0_wp * ibit34 * adv_mom_3 & 3184 + ibit33 * adv_mom_1 & 3185 ) * & 3186 ( w(k+1,j,i) - w(k,j,i) ) & 3187 - ( 5.0_wp * ibit35 * adv_mom_5 & 3188 + ibit34 * adv_mom_3 & 3189 ) * & 3190 ( w(k_pp,j,i) - w(k-1,j,i) ) & 3191 + ( ibit35 * adv_mom_5 & 3192 ) * & 3193 ( w(k_ppp,j,i) - w(k_mm,j,i) ) & 3194 ) 3195 ENDDO 3196 3197 DO k = nzb+1, nzt 3053 3198 3054 3199 flux_d = flux_t(k-1) … … 3057 3202 !-- Calculate the divergence of the velocity field. A respective 3058 3203 !-- correction is needed to overcome numerical instabilities introduced 3059 !-- by a not sufficient reduction of divergences near topography. 3204 !-- by a not sufficient reduction of divergences near topography. 3060 3205 div = ( ( ( u_comp(k) + gu ) * ( ibit27 + ibit28 + ibit29 ) & 3061 3206 - ( u(k+1,j,i) + u(k,j,i) ) & 3062 * ( IBITS(advc_flags_1(k,j,i-1),27,1) & 3063 + IBITS(advc_flags_1(k,j,i-1),28,1) & 3064 + IBITS(advc_flags_1(k,j,i-1),29,1) & 3207 * ( & 3208 REAL( IBITS(advc_flags_1(k,j,i-1),27,1), KIND = wp ) & 3209 + REAL( IBITS(advc_flags_1(k,j,i-1),28,1), KIND = wp ) & 3210 + REAL( IBITS(advc_flags_1(k,j,i-1),29,1), KIND = wp ) & 3065 3211 ) & 3066 3212 ) * ddx & 3067 3213 + ( ( v_comp(k) + gv ) * ( ibit30 + ibit31 + ibit32 ) & 3068 3214 - ( v(k+1,j,i) + v(k,j,i) ) & 3069 * ( IBITS(advc_flags_1(k,j-1,i),30,1) & 3070 + IBITS(advc_flags_1(k,j-1,i),31,1) & 3071 + IBITS(advc_flags_2(k,j-1,i),0,1) & 3215 * ( & 3216 REAL( IBITS(advc_flags_1(k,j-1,i),30,1), KIND = wp ) & 3217 + REAL( IBITS(advc_flags_1(k,j-1,i),31,1), KIND = wp ) & 3218 + REAL( IBITS(advc_flags_2(k,j-1,i),0,1), KIND = wp ) & 3072 3219 ) & 3073 3220 ) * ddy & … … 3075 3222 * ( ibit33 + ibit34 + ibit35 ) & 3076 3223 - ( w(k,j,i) + w(k-1,j,i) ) * rho_air(k) & 3077 * ( IBITS(advc_flags_2(k-1,j,i),1,1) & 3078 + IBITS(advc_flags_2(k-1,j,i),2,1) & 3079 + IBITS(advc_flags_2(k-1,j,i),3,1) & 3224 * ( & 3225 REAL( IBITS(advc_flags_2(k-1,j,i),1,1), KIND = wp ) & 3226 + REAL( IBITS(advc_flags_2(k-1,j,i),2,1), KIND = wp ) & 3227 + REAL( IBITS(advc_flags_2(k-1,j,i),3,1), KIND = wp ) & 3080 3228 ) & 3081 ) * drho_air_zw(k) * ddzu(k+1) &3082 ) * 0.5_wp3083 3084 3085 tend(k,j,i) = tend(k,j,i) - ( &3086 ( flux_r(k) + diss_r(k) &3087 - flux_l_w(k,j,tn) - diss_l_w(k,j,tn) ) * ddx &3088 + ( flux_n(k) + diss_n(k) &3089 - flux_s_w(k,tn) - diss_s_w(k,tn) ) * ddy &3090 + ( ( flux_t(k) + diss_t(k) ) &3091 - ( flux_d + diss_d ) &3092 ) * drho_air_zw(k) * ddzu(k+1) &3093 ) + div * w(k,j,i)3094 3095 flux_l_w(k,j,tn) = flux_r(k)3096 diss_l_w(k,j,tn) = diss_r(k)3097 flux_s_w(k,tn) = flux_n(k)3098 diss_s_w(k,tn) = diss_n(k)3099 !3100 !-- Statistical Evaluation of w'w'.3101 sums_ws2_ws_l(k,tn) = sums_ws2_ws_l(k,tn) &3102 + ( flux_t(k) &3103 * ( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) &3104 / ( w_comp(k) + SIGN( 1.0E-20_wp, w_comp(k) ) ) &3105 + diss_t(k) &3106 * ABS( w_comp(k) - 2.0_wp * hom(k,1,3,0) ) &3107 / ( ABS( w_comp(k) ) + 1.0E-20_wp ) &3108 ) * weight_substep(intermediate_timestep_count)3109 3110 ENDDO3111 3112 DO k = nzb_max+1, nzt3113 3114 u_comp(k) = u(k+1,j,i+1) + u(k,j,i+1) - gu3115 flux_r(k) = u_comp(k) * ( &3116 37.0_wp * ( w(k,j,i+1) + w(k,j,i) ) &3117 - 8.0_wp * ( w(k,j,i+2) + w(k,j,i-1) ) &3118 + ( w(k,j,i+3) + w(k,j,i-2) ) ) * adv_mom_53119 3120 diss_r(k) = - ABS( u_comp(k) ) * ( &3121 10.0_wp * ( w(k,j,i+1) - w(k,j,i) ) &3122 - 5.0_wp * ( w(k,j,i+2) - w(k,j,i-1) ) &3123 + ( w(k,j,i+3) - w(k,j,i-2) ) ) * adv_mom_53124 3125 v_comp(k) = v(k+1,j+1,i) + v(k,j+1,i) - gv3126 flux_n(k) = v_comp(k) * ( &3127 37.0_wp * ( w(k,j+1,i) + w(k,j,i) ) &3128 - 8.0_wp * ( w(k,j+2,i) + w(k,j-1,i) ) &3129 + ( w(k,j+3,i) + w(k,j-2,i) ) ) * adv_mom_53130 3131 diss_n(k) = - ABS( v_comp(k) ) * ( &3132 10.0_wp * ( w(k,j+1,i) - w(k,j,i) ) &3133 - 5.0_wp * ( w(k,j+2,i) - w(k,j-1,i) ) &3134 + ( w(k,j+3,i) - w(k,j-2,i) ) ) * adv_mom_53135 !3136 !-- k index has to be modified near bottom and top, else array3137 !-- subscripts will be exceeded.3138 ibit35 = IBITS(advc_flags_2(k,j,i),3,1)3139 ibit34 = IBITS(advc_flags_2(k,j,i),2,1)3140 ibit33 = IBITS(advc_flags_2(k,j,i),1,1)3141 3142 k_ppp = k + 3 * ibit353143 k_pp = k + 2 * ( 1 - ibit33 )3144 k_mm = k - 2 * ibit353145 3146 w_comp(k) = w(k+1,j,i) + w(k,j,i)3147 flux_t(k) = w_comp(k) * rho_air(k+1) * ( &3148 ( 37.0_wp * ibit35 * adv_mom_5 &3149 + 7.0_wp * ibit34 * adv_mom_3 &3150 + ibit33 * adv_mom_1 &3151 ) * &3152 ( w(k+1,j,i) + w(k,j,i) ) &3153 - ( 8.0_wp * ibit35 * adv_mom_5 &3154 + ibit34 * adv_mom_3 &3155 ) * &3156 ( w(k_pp,j,i) + w(k-1,j,i) ) &3157 + ( ibit35 * adv_mom_5 &3158 ) * &3159 ( w(k_ppp,j,i) + w(k_mm,j,i) ) &3160 )3161 3162 diss_t(k) = - ABS( w_comp(k) ) * rho_air(k+1) * ( &3163 ( 10.0_wp * ibit35 * adv_mom_5 &3164 + 3.0_wp * ibit34 * adv_mom_3 &3165 + ibit33 * adv_mom_1 &3166 ) * &3167 ( w(k+1,j,i) - w(k,j,i) ) &3168 - ( 5.0_wp * ibit35 * adv_mom_5 &3169 + ibit34 * adv_mom_3 &3170 ) * &3171 ( w(k_pp,j,i) - w(k-1,j,i) ) &3172 + ( ibit35 * adv_mom_5 &3173 ) * &3174 ( w(k_ppp,j,i) - w(k_mm,j,i) ) &3175 )3176 ENDDO3177 3178 DO k = nzb_max+1, nzt3179 3180 flux_d = flux_t(k-1)3181 diss_d = diss_t(k-1)3182 !3183 !-- Calculate the divergence of the velocity field. A respective3184 !-- correction is needed to overcome numerical instabilities introduced3185 !-- by a not sufficient reduction of divergences near topography.3186 div = ( ( u_comp(k) + gu - ( u(k+1,j,i) + u(k,j,i) ) ) * ddx &3187 + ( v_comp(k) + gv - ( v(k+1,j,i) + v(k,j,i) ) ) * ddy &3188 + ( w_comp(k) * rho_air(k+1) - &3189 ( w(k,j,i) + w(k-1,j,i) ) * rho_air(k) &3190 3229 ) * drho_air_zw(k) * ddzu(k+1) & 3191 3230 ) * 0.5_wp … … 3260 3299 3261 3300 INTEGER(iwp) :: i !< grid index along x-direction 3262 INTEGER(iwp) :: ibit0 !< flag indicating 1st-order scheme along x-direction3263 INTEGER(iwp) :: ibit1 !< flag indicating 3rd-order scheme along x-direction3264 INTEGER(iwp) :: ibit2 !< flag indicating 5th-order scheme along x-direction3265 INTEGER(iwp) :: ibit3 !< flag indicating 1st-order scheme along y-direction3266 INTEGER(iwp) :: ibit4 !< flag indicating 3rd-order scheme along y-direction3267 INTEGER(iwp) :: ibit5 !< flag indicating 5th-order scheme along y-direction3268 INTEGER(iwp) :: ibit6 !< flag indicating 1st-order scheme along z-direction3269 INTEGER(iwp) :: ibit7 !< flag indicating 3rd-order scheme along z-direction3270 INTEGER(iwp) :: ibit8 !< flag indicating 5th-order scheme along z-direction3271 3301 INTEGER(iwp) :: j !< grid index along y-direction 3272 3302 INTEGER(iwp) :: k !< grid index along z-direction … … 3276 3306 INTEGER(iwp) :: tn = 0 !< number of OpenMP thread 3277 3307 3278 #if defined( __nopointer ) 3279 REAL(wp), DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: sk !< advected scalar 3280 #else 3281 REAL(wp), DIMENSION(:,:,:), POINTER :: sk !< advected scalar 3282 #endif 3283 3308 ! 3309 !-- sk is an array from parameter list. It should not be a pointer, because 3310 !-- in that case the compiler can not assume a stride 1 and cannot perform 3311 !-- a strided one vector load. Adding the CONTIGUOUS keyword makes things 3312 !-- even worse, because the compiler cannot assume strided one in the 3313 !-- caller side. 3314 REAL(wp), INTENT(IN),DIMENSION(nzb:nzt+1,nysg:nyng,nxlg:nxrg) :: sk !< advected scalar 3315 3316 REAL(wp) :: ibit0 !< flag indicating 1st-order scheme along x-direction 3317 REAL(wp) :: ibit1 !< flag indicating 3rd-order scheme along x-direction 3318 REAL(wp) :: ibit2 !< flag indicating 5th-order scheme along x-direction 3319 REAL(wp) :: ibit3 !< flag indicating 1st-order scheme along y-direction 3320 REAL(wp) :: ibit4 !< flag indicating 3rd-order scheme along y-direction 3321 REAL(wp) :: ibit5 !< flag indicating 5th-order scheme along y-direction 3322 REAL(wp) :: ibit6 !< flag indicating 1st-order scheme along z-direction 3323 REAL(wp) :: ibit7 !< flag indicating 3rd-order scheme along z-direction 3324 REAL(wp) :: ibit8 !< flag indicating 5th-order scheme along z-direction 3284 3325 REAL(wp) :: diss_d !< artificial dissipation term at grid box bottom 3285 3326 REAL(wp) :: div !< diverence on scalar grid … … 3309 3350 DO k = nzb+1, nzb_max 3310 3351 3311 ibit2 = IBITS(advc_flags_1(k,j,i-1),2,1)3312 ibit1 = IBITS(advc_flags_1(k,j,i-1),1,1)3313 ibit0 = IBITS(advc_flags_1(k,j,i-1),0,1)3352 ibit2 = REAL( IBITS(advc_flags_1(k,j,i-1),2,1), KIND = wp ) 3353 ibit1 = REAL( IBITS(advc_flags_1(k,j,i-1),1,1), KIND = wp ) 3354 ibit0 = REAL( IBITS(advc_flags_1(k,j,i-1),0,1), KIND = wp ) 3314 3355 3315 3356 u_comp = u(k,j,i) - u_gtrans + u_stokes_zu(k) … … 3370 3411 DO k = nzb+1, nzb_max 3371 3412 3372 ibit5 = IBITS(advc_flags_1(k,j-1,i),5,1)3373 ibit4 = IBITS(advc_flags_1(k,j-1,i),4,1)3374 ibit3 = IBITS(advc_flags_1(k,j-1,i),3,1)3413 ibit5 = REAL( IBITS(advc_flags_1(k,j-1,i),5,1), KIND = wp ) 3414 ibit4 = REAL( IBITS(advc_flags_1(k,j-1,i),4,1), KIND = wp ) 3415 ibit3 = REAL( IBITS(advc_flags_1(k,j-1,i),3,1), KIND = wp ) 3375 3416 3376 3417 v_comp = v(k,j,i) - v_gtrans + v_stokes_zu(k) … … 3433 3474 DO k = nzb+1, nzb_max 3434 3475 3435 ibit2 = IBITS(advc_flags_1(k,j,i),2,1)3436 ibit1 = IBITS(advc_flags_1(k,j,i),1,1)3437 ibit0 = IBITS(advc_flags_1(k,j,i),0,1)3476 ibit2 = REAL( IBITS(advc_flags_1(k,j,i),2,1), KIND = wp ) 3477 ibit1 = REAL( IBITS(advc_flags_1(k,j,i),1,1), KIND = wp ) 3478 ibit0 = REAL( IBITS(advc_flags_1(k,j,i),0,1), KIND = wp ) 3438 3479 3439 3480 u_comp = u(k,j,i+1) - u_gtrans + u_stokes_zu(k) … … 3468 3509 ) 3469 3510 3470 ibit5 = IBITS(advc_flags_1(k,j,i),5,1)3471 ibit4 = IBITS(advc_flags_1(k,j,i),4,1)3472 ibit3 = IBITS(advc_flags_1(k,j,i),3,1)3511 ibit5 = REAL( IBITS(advc_flags_1(k,j,i),5,1), KIND = wp ) 3512 ibit4 = REAL( IBITS(advc_flags_1(k,j,i),4,1), KIND = wp ) 3513 ibit3 = REAL( IBITS(advc_flags_1(k,j,i),3,1), KIND = wp ) 3473 3514 3474 3515 v_comp = v(k,j+1,i) - v_gtrans + v_stokes_zu(k) … … 3505 3546 !-- k index has to be modified near bottom and top, else array 3506 3547 !-- subscripts will be exceeded. 3507 ibit8 = IBITS(advc_flags_1(k,j,i),8,1)3508 ibit7 = IBITS(advc_flags_1(k,j,i),7,1)3509 ibit6 = IBITS(advc_flags_1(k,j,i),6,1)3548 ibit8 = REAL( IBITS(advc_flags_1(k,j,i),8,1), KIND = wp ) 3549 ibit7 = REAL( IBITS(advc_flags_1(k,j,i),7,1), KIND = wp ) 3550 ibit6 = REAL( IBITS(advc_flags_1(k,j,i),6,1), KIND = wp ) 3510 3551 3511 3552 k_ppp = k + 3 * ibit8 … … 3547 3588 !-- by a not sufficient reduction of divergences near topography. 3548 3589 div = ( u(k,j,i+1) * ( ibit0 + ibit1 + ibit2 ) & 3549 - u(k,j,i) * ( IBITS(advc_flags_1(k,j,i-1),0,1) & 3550 + IBITS(advc_flags_1(k,j,i-1),1,1) & 3551 + IBITS(advc_flags_1(k,j,i-1),2,1) & 3590 - u(k,j,i) * ( & 3591 REAL( IBITS(advc_flags_1(k,j,i-1),0,1), KIND = wp ) & 3592 + REAL( IBITS(advc_flags_1(k,j,i-1),1,1), KIND = wp ) & 3593 + REAL( IBITS(advc_flags_1(k,j,i-1),2,1), KIND = wp ) & 3552 3594 ) & 3553 3595 ) * ddx & 3554 3596 + ( v(k,j+1,i) * ( ibit3 + ibit4 + ibit5 ) & 3555 - v(k,j,i) * ( IBITS(advc_flags_1(k,j-1,i),3,1) & 3556 + IBITS(advc_flags_1(k,j-1,i),4,1) & 3557 + IBITS(advc_flags_1(k,j-1,i),5,1) & 3597 - v(k,j,i) * ( & 3598 REAL( IBITS(advc_flags_1(k,j-1,i),3,1), KIND = wp ) & 3599 + REAL( IBITS(advc_flags_1(k,j-1,i),4,1), KIND = wp ) & 3600 + REAL( IBITS(advc_flags_1(k,j-1,i),5,1), KIND = wp ) & 3558 3601 ) & 3559 3602 ) * ddy & … … 3561 3604 ( ibit6 + ibit7 + ibit8 ) & 3562 3605 - w(k-1,j,i) * rho_air_zw(k-1) * & 3563 ( IBITS(advc_flags_1(k-1,j,i),6,1) & 3564 + IBITS(advc_flags_1(k-1,j,i),7,1) & 3565 + IBITS(advc_flags_1(k-1,j,i),8,1) & 3606 ( & 3607 REAL( IBITS(advc_flags_1(k-1,j,i),6,1), KIND = wp ) & 3608 + REAL( IBITS(advc_flags_1(k-1,j,i),7,1), KIND = wp ) & 3609 + REAL( IBITS(advc_flags_1(k-1,j,i),8,1), KIND = wp ) & 3566 3610 ) & 3567 3611 ) * drho_air(k) * ddzw(k) … … 3611 3655 !-- k index has to be modified near bottom and top, else array 3612 3656 !-- subscripts will be exceeded. 3613 ibit8 = IBITS(advc_flags_1(k,j,i),8,1)3614 ibit7 = IBITS(advc_flags_1(k,j,i),7,1)3615 ibit6 = IBITS(advc_flags_1(k,j,i),6,1)3657 ibit8 = REAL( IBITS(advc_flags_1(k,j,i),8,1), KIND = wp ) 3658 ibit7 = REAL( IBITS(advc_flags_1(k,j,i),7,1), KIND = wp ) 3659 ibit6 = REAL( IBITS(advc_flags_1(k,j,i),6,1), KIND = wp ) 3616 3660 3617 3661 k_ppp = k + 3 * ibit8 … … 3817 3861 3818 3862 INTEGER(iwp) :: i !< grid index along x-direction 3819 INTEGER(iwp) :: ibit9 !< flag indicating 1st-order scheme along x-direction3820 INTEGER(iwp) :: ibit10 !< flag indicating 3rd-order scheme along x-direction3821 INTEGER(iwp) :: ibit11 !< flag indicating 5th-order scheme along x-direction3822 INTEGER(iwp) :: ibit12 !< flag indicating 1st-order scheme along y-direction3823 INTEGER(iwp) :: ibit13 !< flag indicating 3rd-order scheme along y-direction3824 INTEGER(iwp) :: ibit14 !< flag indicating 5th-order scheme along y-direction3825 INTEGER(iwp) :: ibit15 !< flag indicating 1st-order scheme along z-direction3826 INTEGER(iwp) :: ibit16 !< flag indicating 3rd-order scheme along z-direction3827 INTEGER(iwp) :: ibit17 !< flag indicating 5th-order scheme along z-direction3828 3863 INTEGER(iwp) :: j !< grid index along y-direction 3829 3864 INTEGER(iwp) :: k !< grid index along z-direction … … 3833 3868 INTEGER(iwp) :: tn = 0 !< number of OpenMP thread 3834 3869 3870 REAL(wp) :: ibit9 !< flag indicating 1st-order scheme along x-direction 3871 REAL(wp) :: ibit10 !< flag indicating 3rd-order scheme along x-direction 3872 REAL(wp) :: ibit11 !< flag indicating 5th-order scheme along x-direction 3873 REAL(wp) :: ibit12 !< flag indicating 1st-order scheme along y-direction 3874 REAL(wp) :: ibit13 !< flag indicating 3rd-order scheme along y-direction 3875 REAL(wp) :: ibit14 !< flag indicating 5th-order scheme along y-direction 3876 REAL(wp) :: ibit15 !< flag indicating 1st-order scheme along z-direction 3877 REAL(wp) :: ibit16 !< flag indicating 3rd-order scheme along z-direction 3878 REAL(wp) :: ibit17 !< flag indicating 5th-order scheme along z-direction 3835 3879 REAL(wp) :: diss_d !< artificial dissipation term at grid box bottom 3836 3880 REAL(wp) :: div !< diverence on u-grid … … 3864 3908 DO k = nzb+1, nzb_max 3865 3909 3866 ibit11 = IBITS(advc_flags_1(k,j,i-1),11,1)3867 ibit10 = IBITS(advc_flags_1(k,j,i-1),10,1)3868 ibit9 = IBITS(advc_flags_1(k,j,i-1),9,1)3910 ibit11 = REAL( IBITS(advc_flags_1(k,j,i-1),11,1), KIND = wp ) 3911 ibit10 = REAL( IBITS(advc_flags_1(k,j,i-1),10,1), KIND = wp ) 3912 ibit9 = REAL( IBITS(advc_flags_1(k,j,i-1),9,1), KIND = wp ) 3869 3913 3870 3914 u_comp(k) = u(k,j,i) + u(k,j,i-1) - gu 3871 3915 swap_flux_x_local_u(k,j) = u_comp(k) * ( & 3872 ( 37.0_wp * ibit11 * adv_mom_5 3873 + 7.0_wp * ibit10 * adv_mom_3 3874 + ibit9 * adv_mom_1 3916 ( 37.0_wp * ibit11 * adv_mom_5 & 3917 + 7.0_wp * ibit10 * adv_mom_3 & 3918 + ibit9 * adv_mom_1 & 3875 3919 ) * & 3876 3920 ( u(k,j,i) + u(k,j,i-1) ) & 3877 - ( 8.0_wp * ibit11 * adv_mom_5 3878 + ibit10 * adv_mom_3 3921 - ( 8.0_wp * ibit11 * adv_mom_5 & 3922 + ibit10 * adv_mom_3 & 3879 3923 ) * & 3880 3924 ( u(k,j,i+1) + u(k,j,i-2) ) & 3881 + ( ibit11 * adv_mom_5 3925 + ( ibit11 * adv_mom_5 & 3882 3926 ) * & 3883 3927 ( u(k,j,i+2) + u(k,j,i-3) ) & … … 3885 3929 3886 3930 swap_diss_x_local_u(k,j) = - ABS( u_comp(k) ) * ( & 3887 ( 10.0_wp * ibit11 * adv_mom_5 3888 + 3.0_wp * ibit10 * adv_mom_3 3889 + ibit9 * adv_mom_1 3931 ( 10.0_wp * ibit11 * adv_mom_5 & 3932 + 3.0_wp * ibit10 * adv_mom_3 & 3933 + ibit9 * adv_mom_1 & 3890 3934 ) * & 3891 3935 ( u(k,j,i) - u(k,j,i-1) ) & 3892 - ( 5.0_wp * ibit11 * adv_mom_5 3893 + ibit10 * adv_mom_3 3936 - ( 5.0_wp * ibit11 * adv_mom_5 & 3937 + ibit10 * adv_mom_3 & 3894 3938 ) * & 3895 3939 ( u(k,j,i+1) - u(k,j,i-2) ) & 3896 + ( ibit11 * adv_mom_5 3940 + ( ibit11 * adv_mom_5 & 3897 3941 ) * & 3898 3942 ( u(k,j,i+2) - u(k,j,i-3) ) & … … 3905 3949 u_comp(k) = u(k,j,i) + u(k,j,i-1) - gu 3906 3950 swap_flux_x_local_u(k,j) = u_comp(k) * ( & 3907 37.0_wp * ( u(k,j,i) + u(k,j,i-1) ) 3908 - 8.0_wp * ( u(k,j,i+1) + u(k,j,i-2) ) 3951 37.0_wp * ( u(k,j,i) + u(k,j,i-1) ) & 3952 - 8.0_wp * ( u(k,j,i+1) + u(k,j,i-2) ) & 3909 3953 + ( u(k,j,i+2) + u(k,j,i-3) ) ) * adv_mom_5 3910 3954 swap_diss_x_local_u(k,j) = - ABS(u_comp(k)) * ( & 3911 10.0_wp * ( u(k,j,i) - u(k,j,i-1) ) 3912 - 5.0_wp * ( u(k,j,i+1) - u(k,j,i-2) ) 3955 10.0_wp * ( u(k,j,i) - u(k,j,i-1) ) & 3956 - 5.0_wp * ( u(k,j,i+1) - u(k,j,i-2) ) & 3913 3957 + ( u(k,j,i+2) - u(k,j,i-3) ) ) * adv_mom_5 3914 3958 … … 3922 3966 DO k = nzb+1, nzb_max 3923 3967 3924 ibit14 = IBITS(advc_flags_1(k,j-1,i),14,1)3925 ibit13 = IBITS(advc_flags_1(k,j-1,i),13,1)3926 ibit12 = IBITS(advc_flags_1(k,j-1,i),12,1)3968 ibit14 = REAL( IBITS(advc_flags_1(k,j-1,i),14,1), KIND = wp ) 3969 ibit13 = REAL( IBITS(advc_flags_1(k,j-1,i),13,1), KIND = wp ) 3970 ibit12 = REAL( IBITS(advc_flags_1(k,j-1,i),12,1), KIND = wp ) 3927 3971 3928 3972 v_comp = v(k,j,i) + v(k,j,i-1) - gv 3929 3973 swap_flux_y_local_u(k) = v_comp * ( & 3930 ( 37.0_wp * ibit14 * adv_mom_5 3931 + 7.0_wp * ibit13 * adv_mom_3 3932 + ibit12 * adv_mom_1 3974 ( 37.0_wp * ibit14 * adv_mom_5 & 3975 + 7.0_wp * ibit13 * adv_mom_3 & 3976 + ibit12 * adv_mom_1 & 3933 3977 ) * & 3934 3978 ( u(k,j,i) + u(k,j-1,i) ) & 3935 - ( 8.0_wp * ibit14 * adv_mom_5 3936 + ibit13 * adv_mom_3 3979 - ( 8.0_wp * ibit14 * adv_mom_5 & 3980 + ibit13 * adv_mom_3 & 3937 3981 ) * & 3938 3982 ( u(k,j+1,i) + u(k,j-2,i) ) & 3939 + ( ibit14 * adv_mom_5 3983 + ( ibit14 * adv_mom_5 & 3940 3984 ) * & 3941 3985 ( u(k,j+2,i) + u(k,j-3,i) ) & … … 3943 3987 3944 3988 swap_diss_y_local_u(k) = - ABS ( v_comp ) * ( & 3945 ( 10.0_wp * ibit14 * adv_mom_5 3946 + 3.0_wp * ibit13 * adv_mom_3 3947 + ibit12 * adv_mom_1 3989 ( 10.0_wp * ibit14 * adv_mom_5 & 3990 + 3.0_wp * ibit13 * adv_mom_3 & 3991 + ibit12 * adv_mom_1 & 3948 3992 ) * & 3949 3993 ( u(k,j,i) - u(k,j-1,i) ) & 3950 - ( 5.0_wp * ibit14 * adv_mom_5 3951 + ibit13 * adv_mom_3 3994 - ( 5.0_wp * ibit14 * adv_mom_5 & 3995 + ibit13 * adv_mom_3 & 3952 3996 ) * & 3953 3997 ( u(k,j+1,i) - u(k,j-2,i) ) & 3954 + ( ibit14 * adv_mom_5 3998 + ( ibit14 * adv_mom_5 & 3955 3999 ) * & 3956 4000 ( u(k,j+2,i) - u(k,j-3,i) ) & … … 3983 4027 DO k = nzb+1, nzb_max 3984 4028 3985 ibit11 = IBITS(advc_flags_1(k,j,i),11,1)3986 ibit10 = IBITS(advc_flags_1(k,j,i),10,1)3987 ibit9 = IBITS(advc_flags_1(k,j,i),9,1)4029 ibit11 = REAL( IBITS(advc_flags_1(k,j,i),11,1), KIND = wp ) 4030 ibit10 = REAL( IBITS(advc_flags_1(k,j,i),10,1), KIND = wp ) 4031 ibit9 = REAL( IBITS(advc_flags_1(k,j,i),9,1), KIND = wp ) 3988 4032 3989 4033 u_comp(k) = u(k,j,i+1) + u(k,j,i) 3990 4034 flux_r(k) = ( u_comp(k) - gu ) * ( & 3991 ( 37.0_wp * ibit11 * adv_mom_5 3992 + 7.0_wp * ibit10 * adv_mom_3 3993 + ibit9 * adv_mom_1 4035 ( 37.0_wp * ibit11 * adv_mom_5 & 4036 + 7.0_wp * ibit10 * adv_mom_3 & 4037 + ibit9 * adv_mom_1 & 3994 4038 ) * & 3995 4039 ( u(k,j,i+1) + u(k,j,i) ) & 3996 - ( 8.0_wp * ibit11 * adv_mom_5 3997 + ibit10 * adv_mom_3 4040 - ( 8.0_wp * ibit11 * adv_mom_5 & 4041 + ibit10 * adv_mom_3 & 3998 4042 ) * & 3999 4043 ( u(k,j,i+2) + u(k,j,i-1) ) & 4000 + ( ibit11 * adv_mom_5 4044 + ( ibit11 * adv_mom_5 & 4001 4045 ) * & 4002 4046 ( u(k,j,i+3) + u(k,j,i-2) ) & … … 4004 4048 4005 4049 diss_r(k) = - ABS( u_comp(k) - gu ) * ( & 4006 ( 10.0_wp * ibit11 * adv_mom_5 4007 + 3.0_wp * ibit10 * adv_mom_3 4008 + ibit9 * adv_mom_1 4050 ( 10.0_wp * ibit11 * adv_mom_5 & 4051 + 3.0_wp * ibit10 * adv_mom_3 & 4052 + ibit9 * adv_mom_1 & 4009 4053 ) * & 4010 4054 ( u(k,j,i+1) - u(k,j,i) ) & 4011 - ( 5.0_wp * ibit11 * adv_mom_5 4012 + ibit10 * adv_mom_3 4055 - ( 5.0_wp * ibit11 * adv_mom_5 & 4056 + ibit10 * adv_mom_3 & 4013 4057 ) * & 4014 4058 ( u(k,j,i+2) - u(k,j,i-1) ) & 4015 + ( ibit11 * adv_mom_5 4059 + ( ibit11 * adv_mom_5 & 4016 4060 ) * & 4017 4061 ( u(k,j,i+3) - u(k,j,i-2) ) & 4018 4062 ) 4019 4063 4020 ibit14 = IBITS(advc_flags_1(k,j,i),14,1)4021 ibit13 = IBITS(advc_flags_1(k,j,i),13,1)4022 ibit12 = IBITS(advc_flags_1(k,j,i),12,1)4064 ibit14 = REAL( IBITS(advc_flags_1(k,j,i),14,1), KIND = wp ) 4065 ibit13 = REAL( IBITS(advc_flags_1(k,j,i),13,1), KIND = wp ) 4066 ibit12 = REAL( IBITS(advc_flags_1(k,j,i),12,1), KIND = wp ) 4023 4067 4024 4068 v_comp = v(k,j+1,i) + v(k,j+1,i-1) - gv 4025 4069 flux_n(k) = v_comp * ( & 4026 ( 37.0_wp * ibit14 * adv_mom_5 4027 + 7.0_wp * ibit13 * adv_mom_3 4028 + ibit12 * adv_mom_1 4070 ( 37.0_wp * ibit14 * adv_mom_5 & 4071 + 7.0_wp * ibit13 * adv_mom_3 & 4072 + ibit12 * adv_mom_1 & 4029 4073 ) * & 4030 4074 ( u(k,j+1,i) + u(k,j,i) ) & 4031 - ( 8.0_wp * ibit14 * adv_mom_5 4032 + ibit13 * adv_mom_3 4075 - ( 8.0_wp * ibit14 * adv_mom_5 & 4076 + ibit13 * adv_mom_3 & 4033 4077 ) * & 4034 4078 ( u(k,j+2,i) + u(k,j-1,i) ) & 4035 + ( ibit14 * adv_mom_5 4079 + ( ibit14 * adv_mom_5 & 4036 4080 ) * & 4037 4081 ( u(k,j+3,i) + u(k,j-2,i) ) & … … 4039 4083 4040 4084 diss_n(k) = - ABS ( v_comp ) * ( & 4041 ( 10.0_wp * ibit14 * adv_mom_5 4042 + 3.0_wp * ibit13 * adv_mom_3 4043 + ibit12 * adv_mom_1 4085 ( 10.0_wp * ibit14 * adv_mom_5 & 4086 + 3.0_wp * ibit13 * adv_mom_3 & 4087 + ibit12 * adv_mom_1 & 4044 4088 ) * & 4045 4089 ( u(k,j+1,i) - u(k,j,i) ) & 4046 - ( 5.0_wp * ibit14 * adv_mom_5 4047 + ibit13 * adv_mom_3 4090 - ( 5.0_wp * ibit14 * adv_mom_5 & 4091 + ibit13 * adv_mom_3 & 4048 4092 ) * & 4049 4093 ( u(k,j+2,i) - u(k,j-1,i) ) & 4050 + ( ibit14 * adv_mom_5 4094 + ( ibit14 * adv_mom_5 & 4051 4095 ) * & 4052 4096 ( u(k,j+3,i) - u(k,j-2,i) ) & … … 4055 4099 !-- k index has to be modified near bottom and top, else array 4056 4100 !-- subscripts will be exceeded. 4057 ibit17 = IBITS(advc_flags_1(k,j,i),17,1)4058 ibit16 = IBITS(advc_flags_1(k,j,i),16,1)4059 ibit15 = IBITS(advc_flags_1(k,j,i),15,1)4101 ibit17 = REAL( IBITS(advc_flags_1(k,j,i),17,1), KIND = wp ) 4102 ibit16 = REAL( IBITS(advc_flags_1(k,j,i),16,1), KIND = wp ) 4103 ibit15 = REAL( IBITS(advc_flags_1(k,j,i),15,1), KIND = wp ) 4060 4104 4061 4105 k_ppp = k + 3 * ibit17 … … 4065 4109 w_comp = w(k,j,i) + w(k,j,i-1) 4066 4110 flux_t(k) = w_comp * rho_air_zw(k) * ( & 4067 ( 37.0_wp * ibit17 * adv_mom_5 4068 + 7.0_wp * ibit16 * adv_mom_3 4069 + ibit15 * adv_mom_1 4111 ( 37.0_wp * ibit17 * adv_mom_5 & 4112 + 7.0_wp * ibit16 * adv_mom_3 & 4113 + ibit15 * adv_mom_1 & 4070 4114 ) * & 4071 4115 ( u(k+1,j,i) + u(k,j,i) ) & 4072 - ( 8.0_wp * ibit17 * adv_mom_5 4073 + ibit16 * adv_mom_3 4116 - ( 8.0_wp * ibit17 * adv_mom_5 & 4117 + ibit16 * adv_mom_3 & 4074 4118 ) * & 4075 4119 ( u(k_pp,j,i) + u(k-1,j,i) ) & 4076 + ( ibit17 * adv_mom_5 4120 + ( ibit17 * adv_mom_5 & 4077 4121 ) * & 4078 4122 ( u(k_ppp,j,i) + u(k_mm,j,i) ) & … … 4080 4124 4081 4125 diss_t(k) = - ABS( w_comp ) * rho_air_zw(k) * ( & 4082 ( 10.0_wp * ibit17 * adv_mom_5 4083 + 3.0_wp * ibit16 * adv_mom_3 4084 + ibit15 * adv_mom_1 4126 ( 10.0_wp * ibit17 * adv_mom_5 & 4127 + 3.0_wp * ibit16 * adv_mom_3 & 4128 + ibit15 * adv_mom_1 & 4085 4129 ) * & 4086 4130 ( u(k+1,j,i) - u(k,j,i) ) & 4087 - ( 5.0_wp * ibit17 * adv_mom_5 4088 + ibit16 * adv_mom_3 4131 - ( 5.0_wp * ibit17 * adv_mom_5 & 4132 + ibit16 * adv_mom_3 & 4089 4133 ) * & 4090 4134 ( u(k_pp,j,i) - u(k-1,j,i) ) & 4091 + ( ibit17 * adv_mom_5 4135 + ( ibit17 * adv_mom_5 & 4092 4136 ) * & 4093 4137 ( u(k_ppp,j,i) - u(k_mm,j,i) ) & … … 4099 4143 div = ( ( u_comp(k) * ( ibit9 + ibit10 + ibit11 ) & 4100 4144 - ( u(k,j,i) + u(k,j,i-1) ) & 4101 * ( IBITS(advc_flags_1(k,j,i-1),9,1) & 4102 + IBITS(advc_flags_1(k,j,i-1),10,1) & 4103 + IBITS(advc_flags_1(k,j,i-1),11,1) & 4145 * ( & 4146 REAL( IBITS(advc_flags_1(k,j,i-1),9,1), KIND = wp ) & 4147 + REAL( IBITS(advc_flags_1(k,j,i-1),10,1), KIND = wp ) & 4148 + REAL( IBITS(advc_flags_1(k,j,i-1),11,1), KIND = wp ) & 4104 4149 ) & 4105 4150 ) * ddx & 4106 4151 + ( ( v_comp + gv ) * ( ibit12 + ibit13 + ibit14 ) & 4107 4152 - ( v(k,j,i) + v(k,j,i-1 ) ) & 4108 * ( IBITS(advc_flags_1(k,j-1,i),12,1) & 4109 + IBITS(advc_flags_1(k,j-1,i),13,1) & 4110 + IBITS(advc_flags_1(k,j-1,i),14,1) & 4153 * ( & 4154 REAL( IBITS(advc_flags_1(k,j-1,i),12,1), KIND = wp ) & 4155 + REAL( IBITS(advc_flags_1(k,j-1,i),13,1), KIND = wp ) & 4156 + REAL( IBITS(advc_flags_1(k,j-1,i),14,1), KIND = wp ) & 4111 4157 ) & 4112 4158 ) * ddy & 4113 4159 + ( w_comp * rho_air_zw(k) * ( ibit15 + ibit16 + ibit17 ) & 4114 4160 - ( w(k-1,j,i) + w(k-1,j,i-1) ) * rho_air_zw(k-1) & 4115 * ( IBITS(advc_flags_1(k-1,j,i),15,1) & 4116 + IBITS(advc_flags_1(k-1,j,i),16,1) & 4117 + IBITS(advc_flags_1(k-1,j,i),17,1) & 4161 * ( & 4162 REAL( IBITS(advc_flags_1(k-1,j,i),15,1), KIND = wp ) & 4163 + REAL( IBITS(advc_flags_1(k-1,j,i),16,1), KIND = wp ) & 4164 + REAL( IBITS(advc_flags_1(k-1,j,i),17,1), KIND = wp ) & 4118 4165 ) & 4119 4166 ) * drho_air(k) * ddzw(k) & … … 4186 4233 !-- k index has to be modified near bottom and top, else array 4187 4234 !-- subscripts will be exceeded. 4188 ibit17 = IBITS(advc_flags_1(k,j,i),17,1)4189 ibit16 = IBITS(advc_flags_1(k,j,i),16,1)4190 ibit15 = IBITS(advc_flags_1(k,j,i),15,1)4235 ibit17 = REAL( IBITS(advc_flags_1(k,j,i),17,1), KIND = wp ) 4236 ibit16 = REAL( IBITS(advc_flags_1(k,j,i),16,1), KIND = wp ) 4237 ibit15 = REAL( IBITS(advc_flags_1(k,j,i),15,1), KIND = wp ) 4191 4238 4192 4239 k_ppp = k + 3 * ibit17 … … 4309 4356 4310 4357 INTEGER(iwp) :: i !< grid index along x-direction 4311 INTEGER(iwp) :: ibit18 !< flag indicating 1st-order scheme along x-direction4312 INTEGER(iwp) :: ibit19 !< flag indicating 3rd-order scheme along x-direction4313 INTEGER(iwp) :: ibit20 !< flag indicating 5th-order scheme along x-direction4314 INTEGER(iwp) :: ibit21 !< flag indicating 1st-order scheme along y-direction4315 INTEGER(iwp) :: ibit22 !< flag indicating 3rd-order scheme along y-direction4316 INTEGER(iwp) :: ibit23 !< flag indicating 5th-order scheme along y-direction4317 INTEGER(iwp) :: ibit24 !< flag indicating 1st-order scheme along z-direction4318 INTEGER(iwp) :: ibit25 !< flag indicating 3rd-order scheme along z-direction4319 INTEGER(iwp) :: ibit26 !< flag indicating 5th-order scheme along z-direction4320 4358 INTEGER(iwp) :: j !< grid index along y-direction 4321 4359 INTEGER(iwp) :: k !< grid index along z-direction … … 4325 4363 INTEGER(iwp) :: tn = 0 !< number of OpenMP thread 4326 4364 4365 REAL(wp) :: ibit18 !< flag indicating 1st-order scheme along x-direction 4366 REAL(wp) :: ibit19 !< flag indicating 3rd-order scheme along x-direction 4367 REAL(wp) :: ibit20 !< flag indicating 5th-order scheme along x-direction 4368 REAL(wp) :: ibit21 !< flag indicating 1st-order scheme along y-direction 4369 REAL(wp) :: ibit22 !< flag indicating 3rd-order scheme along y-direction 4370 REAL(wp) :: ibit23 !< flag indicating 5th-order scheme along y-direction 4371 REAL(wp) :: ibit24 !< flag indicating 1st-order scheme along z-direction 4372 REAL(wp) :: ibit25 !< flag indicating 3rd-order scheme along z-direction 4373 REAL(wp) :: ibit26 !< flag indicating 5th-order scheme along z-direction 4327 4374 REAL(wp) :: diss_d !< artificial dissipation term at grid box bottom 4328 4375 REAL(wp) :: div !< diverence on v-grid … … 4355 4402 DO k = nzb+1, nzb_max 4356 4403 4357 ibit20 = IBITS(advc_flags_1(k,j,i-1),20,1)4358 ibit19 = IBITS(advc_flags_1(k,j,i-1),19,1)4359 ibit18 = IBITS(advc_flags_1(k,j,i-1),18,1)4404 ibit20 = REAL( IBITS(advc_flags_1(k,j,i-1),20,1), KIND = wp ) 4405 ibit19 = REAL( IBITS(advc_flags_1(k,j,i-1),19,1), KIND = wp ) 4406 ibit18 = REAL( IBITS(advc_flags_1(k,j,i-1),18,1), KIND = wp ) 4360 4407 4361 4408 u_comp = u(k,j-1,i) + u(k,j,i) - gu … … 4413 4460 DO k = nzb+1, nzb_max 4414 4461 4415 ibit23 = IBITS(advc_flags_1(k,j-1,i),23,1)4416 ibit22 = IBITS(advc_flags_1(k,j-1,i),22,1)4417 ibit21 = IBITS(advc_flags_1(k,j-1,i),21,1)4462 ibit23 = REAL( IBITS(advc_flags_1(k,j-1,i),23,1), KIND = wp ) 4463 ibit22 = REAL( IBITS(advc_flags_1(k,j-1,i),22,1), KIND = wp ) 4464 ibit21 = REAL( IBITS(advc_flags_1(k,j-1,i),21,1), KIND = wp ) 4418 4465 4419 4466 v_comp(k) = v(k,j,i) + v(k,j-1,i) - gv … … 4473 4520 DO k = nzb+1, nzb_max 4474 4521 4475 ibit20 = IBITS(advc_flags_1(k,j,i),20,1)4476 ibit19 = IBITS(advc_flags_1(k,j,i),19,1)4477 ibit18 = IBITS(advc_flags_1(k,j,i),18,1)4522 ibit20 = REAL( IBITS(advc_flags_1(k,j,i),20,1), KIND = wp ) 4523 ibit19 = REAL( IBITS(advc_flags_1(k,j,i),19,1), KIND = wp ) 4524 ibit18 = REAL( IBITS(advc_flags_1(k,j,i),18,1), KIND = wp ) 4478 4525 4479 4526 u_comp = u(k,j-1,i+1) + u(k,j,i+1) - gu … … 4508 4555 ) 4509 4556 4510 ibit23 = IBITS(advc_flags_1(k,j,i),23,1)4511 ibit22 = IBITS(advc_flags_1(k,j,i),22,1)4512 ibit21 = IBITS(advc_flags_1(k,j,i),21,1)4557 ibit23 = REAL( IBITS(advc_flags_1(k,j,i),23,1), KIND = wp ) 4558 ibit22 = REAL( IBITS(advc_flags_1(k,j,i),22,1), KIND = wp ) 4559 ibit21 = REAL( IBITS(advc_flags_1(k,j,i),21,1), KIND = wp ) 4513 4560 4514 4561 v_comp(k) = v(k,j+1,i) + v(k,j,i) … … 4545 4592 !-- k index has to be modified near bottom and top, else array 4546 4593 !-- subscripts will be exceeded. 4547 ibit26 = IBITS(advc_flags_1(k,j,i),26,1)4548 ibit25 = IBITS(advc_flags_1(k,j,i),25,1)4549 ibit24 = IBITS(advc_flags_1(k,j,i),24,1)4594 ibit26 = REAL( IBITS(advc_flags_1(k,j,i),26,1), KIND = wp ) 4595 ibit25 = REAL( IBITS(advc_flags_1(k,j,i),25,1), KIND = wp ) 4596 ibit24 = REAL( IBITS(advc_flags_1(k,j,i),24,1), KIND = wp ) 4550 4597 4551 4598 k_ppp = k + 3 * ibit26 … … 4590 4637 * ( ibit18 + ibit19 + ibit20 ) & 4591 4638 - ( u(k,j-1,i) + u(k,j,i) ) & 4592 * ( IBITS(advc_flags_1(k,j,i-1),18,1) & 4593 + IBITS(advc_flags_1(k,j,i-1),19,1) & 4594 + IBITS(advc_flags_1(k,j,i-1),20,1) & 4639 * ( & 4640 REAL( IBITS(advc_flags_1(k,j,i-1),18,1), KIND = wp ) & 4641 + REAL( IBITS(advc_flags_1(k,j,i-1),19,1), KIND = wp ) & 4642 + REAL( IBITS(advc_flags_1(k,j,i-1),20,1), KIND = wp ) & 4595 4643 ) & 4596 4644 ) * ddx & … … 4598 4646 * ( ibit21 + ibit22 + ibit23 ) & 4599 4647 - ( v(k,j,i) + v(k,j-1,i) ) & 4600 * ( IBITS(advc_flags_1(k,j-1,i),21,1) & 4601 + IBITS(advc_flags_1(k,j-1,i),22,1) & 4602 + IBITS(advc_flags_1(k,j-1,i),23,1) & 4648 * ( & 4649 REAL( IBITS(advc_flags_1(k,j-1,i),21,1), KIND = wp ) & 4650 + REAL( IBITS(advc_flags_1(k,j-1,i),22,1), KIND = wp ) & 4651 + REAL( IBITS(advc_flags_1(k,j-1,i),23,1), KIND = wp ) & 4603 4652 ) & 4604 4653 ) * ddy & … … 4606 4655 * ( ibit24 + ibit25 + ibit26 ) & 4607 4656 - ( w(k-1,j-1,i) + w(k-1,j,i) ) * rho_air_zw(k-1) & 4608 * ( IBITS(advc_flags_1(k-1,j,i),24,1) & 4609 + IBITS(advc_flags_1(k-1,j,i),25,1) & 4610 + IBITS(advc_flags_1(k-1,j,i),26,1) & 4657 * ( & 4658 REAL( IBITS(advc_flags_1(k-1,j,i),24,1), KIND = wp ) & 4659 + REAL( IBITS(advc_flags_1(k-1,j,i),25,1), KIND = wp ) & 4660 + REAL( IBITS(advc_flags_1(k-1,j,i),26,1), KIND = wp ) & 4611 4661 ) & 4612 4662 ) * drho_air(k) * ddzw(k) & … … 4684 4734 !-- k index has to be modified near bottom and top, else array 4685 4735 !-- subscripts will be exceeded. 4686 ibit26 = IBITS(advc_flags_1(k,j,i),26,1)4687 ibit25 = IBITS(advc_flags_1(k,j,i),25,1)4688 ibit24 = IBITS(advc_flags_1(k,j,i),24,1)4736 ibit26 = REAL( IBITS(advc_flags_1(k,j,i),26,1), KIND = wp ) 4737 ibit25 = REAL( IBITS(advc_flags_1(k,j,i),25,1), KIND = wp ) 4738 ibit24 = REAL( IBITS(advc_flags_1(k,j,i),24,1), KIND = wp ) 4689 4739 4690 4740 k_ppp = k + 3 * ibit26 … … 4811 4861 4812 4862 INTEGER(iwp) :: i !< grid index along x-direction 4813 INTEGER(iwp) :: ibit27 !< flag indicating 1st-order scheme along x-direction4814 INTEGER(iwp) :: ibit28 !< flag indicating 3rd-order scheme along x-direction4815 INTEGER(iwp) :: ibit29 !< flag indicating 5th-order scheme along x-direction4816 INTEGER(iwp) :: ibit30 !< flag indicating 1st-order scheme along y-direction4817 INTEGER(iwp) :: ibit31 !< flag indicating 3rd-order scheme along y-direction4818 INTEGER(iwp) :: ibit32 !< flag indicating 5th-order scheme along y-direction4819 INTEGER(iwp) :: ibit33 !< flag indicating 1st-order scheme along z-direction4820 INTEGER(iwp) :: ibit34 !< flag indicating 3rd-order scheme along z-direction4821 INTEGER(iwp) :: ibit35 !< flag indicating 5th-order scheme along z-direction4822 4863 INTEGER(iwp) :: j !< grid index along y-direction 4823 4864 INTEGER(iwp) :: k !< grid index along z-direction … … 4827 4868 INTEGER(iwp) :: tn = 0 !< number of OpenMP thread 4828 4869 4870 REAL(wp) :: ibit27 !< flag indicating 1st-order scheme along x-direction 4871 REAL(wp) :: ibit28 !< flag indicating 3rd-order scheme along x-direction 4872 REAL(wp) :: ibit29 !< flag indicating 5th-order scheme along x-direction 4873 REAL(wp) :: ibit30 !< flag indicating 1st-order scheme along y-direction 4874 REAL(wp) :: ibit31 !< flag indicating 3rd-order scheme along y-direction 4875 REAL(wp) :: ibit32 !< flag indicating 5th-order scheme along y-direction 4876 REAL(wp) :: ibit33 !< flag indicating 1st-order scheme along z-direction 4877 REAL(wp) :: ibit34 !< flag indicating 3rd-order scheme along z-direction 4878 REAL(wp) :: ibit35 !< flag indicating 5th-order scheme along z-direction 4829 4879 REAL(wp) :: diss_d !< artificial dissipation term at grid box bottom 4830 4880 REAL(wp) :: div !< divergence on w-grid … … 4857 4907 DO k = nzb+1, nzb_max 4858 4908 4859 ibit29 = IBITS(advc_flags_1(k,j,i-1),29,1)4860 ibit28 = IBITS(advc_flags_1(k,j,i-1),28,1)4861 ibit27 = IBITS(advc_flags_1(k,j,i-1),27,1)4909 ibit29 = REAL( IBITS(advc_flags_1(k,j,i-1),29,1), KIND = wp ) 4910 ibit28 = REAL( IBITS(advc_flags_1(k,j,i-1),28,1), KIND = wp ) 4911 ibit27 = REAL( IBITS(advc_flags_1(k,j,i-1),27,1), KIND = wp ) 4862 4912 4863 4913 u_comp = u(k+1,j,i) + u(k,j,i) - gu … … 4915 4965 DO k = nzb+1, nzb_max 4916 4966 4917 ibit32 = IBITS(advc_flags_2(k,j-1,i),0,1)4918 ibit31 = IBITS(advc_flags_1(k,j-1,i),31,1)4919 ibit30 = IBITS(advc_flags_1(k,j-1,i),30,1)4967 ibit32 = REAL( IBITS(advc_flags_2(k,j-1,i),0,1), KIND = wp ) 4968 ibit31 = REAL( IBITS(advc_flags_1(k,j-1,i),31,1), KIND = wp ) 4969 ibit30 = REAL( IBITS(advc_flags_1(k,j-1,i),30,1), KIND = wp ) 4920 4970 4921 4971 v_comp = v(k+1,j,i) + v(k,j,i) - gv … … 4981 5031 DO k = nzb+1, nzb_max 4982 5032 4983 ibit29 = IBITS(advc_flags_1(k,j,i),29,1)4984 ibit28 = IBITS(advc_flags_1(k,j,i),28,1)4985 ibit27 = IBITS(advc_flags_1(k,j,i),27,1)5033 ibit29 = REAL( IBITS(advc_flags_1(k,j,i),29,1), KIND = wp ) 5034 ibit28 = REAL( IBITS(advc_flags_1(k,j,i),28,1), KIND = wp ) 5035 ibit27 = REAL( IBITS(advc_flags_1(k,j,i),27,1), KIND = wp ) 4986 5036 4987 5037 u_comp = u(k+1,j,i+1) + u(k,j,i+1) - gu … … 5016 5066 ) 5017 5067 5018 ibit32 = IBITS(advc_flags_2(k,j,i),0,1)5019 ibit31 = IBITS(advc_flags_1(k,j,i),31,1)5020 ibit30 = IBITS(advc_flags_1(k,j,i),30,1)5068 ibit32 = REAL( IBITS(advc_flags_2(k,j,i),0,1), KIND = wp ) 5069 ibit31 = REAL( IBITS(advc_flags_1(k,j,i),31,1), KIND = wp ) 5070 ibit30 = REAL( IBITS(advc_flags_1(k,j,i),30,1), KIND = wp ) 5021 5071 5022 5072 v_comp = v(k+1,j+1,i) + v(k,j+1,i) - gv … … 5053 5103 !-- k index has to be modified near bottom and top, else array 5054 5104 !-- subscripts will be exceeded. 5055 ibit35 = IBITS(advc_flags_2(k,j,i),3,1)5056 ibit34 = IBITS(advc_flags_2(k,j,i),2,1)5057 ibit33 = IBITS(advc_flags_2(k,j,i),1,1)5105 ibit35 = REAL( IBITS(advc_flags_2(k,j,i),3,1), KIND = wp ) 5106 ibit34 = REAL( IBITS(advc_flags_2(k,j,i),2,1), KIND = wp ) 5107 ibit33 = REAL( IBITS(advc_flags_2(k,j,i),1,1), KIND = wp ) 5058 5108 5059 5109 k_ppp = k + 3 * ibit35 … … 5097 5147 div = ( ( ( u_comp + gu ) * ( ibit27 + ibit28 + ibit29 ) & 5098 5148 - ( u(k+1,j,i) + u(k,j,i) ) & 5099 * ( IBITS(advc_flags_1(k,j,i-1),27,1) & 5100 + IBITS(advc_flags_1(k,j,i-1),28,1) & 5101 + IBITS(advc_flags_1(k,j,i-1),29,1) & 5149 * ( & 5150 REAL( IBITS(advc_flags_1(k,j,i-1),27,1), KIND = wp ) & 5151 + REAL( IBITS(advc_flags_1(k,j,i-1),28,1), KIND = wp ) & 5152 + REAL( IBITS(advc_flags_1(k,j,i-1),29,1), KIND = wp ) & 5102 5153 ) & 5103 5154 ) * ddx & 5104 5155 + ( ( v_comp + gv ) * ( ibit30 + ibit31 + ibit32 ) & 5105 5156 - ( v(k+1,j,i) + v(k,j,i) ) & 5106 * ( IBITS(advc_flags_1(k,j-1,i),30,1) & 5107 + IBITS(advc_flags_1(k,j-1,i),31,1) & 5108 + IBITS(advc_flags_2(k,j-1,i),0,1) & 5157 * ( & 5158 REAL( IBITS(advc_flags_1(k,j-1,i),30,1), KIND = wp ) & 5159 + REAL( IBITS(advc_flags_1(k,j-1,i),31,1), KIND = wp ) & 5160 + REAL( IBITS(advc_flags_2(k,j-1,i),0,1), KIND = wp ) & 5109 5161 ) & 5110 5162 ) * ddy & 5111 5163 + ( w_comp * rho_air(k+1) * ( ibit33 + ibit34 + ibit35 ) & 5112 5164 - ( w(k,j,i) + w(k-1,j,i) ) * rho_air(k) & 5113 * ( IBITS(advc_flags_2(k-1,j,i),1,1) & 5114 + IBITS(advc_flags_2(k-1,j,i),2,1) & 5115 + IBITS(advc_flags_2(k-1,j,i),3,1) & 5165 * ( & 5166 REAL( IBITS(advc_flags_2(k-1,j,i),1,1), KIND = wp ) & 5167 + REAL( IBITS(advc_flags_2(k-1,j,i),2,1), KIND = wp ) & 5168 + REAL( IBITS(advc_flags_2(k-1,j,i),3,1), KIND = wp ) & 5116 5169 ) & 5117 5170 ) * drho_air_zw(k) * ddzu(k+1) & … … 5176 5229 !-- k index has to be modified near bottom and top, else array 5177 5230 !-- subscripts will be exceeded. 5178 ibit35 = IBITS(advc_flags_2(k,j,i),3,1)5179 ibit34 = IBITS(advc_flags_2(k,j,i),2,1)5180 ibit33 = IBITS(advc_flags_2(k,j,i),1,1)5231 ibit35 = REAL( IBITS(advc_flags_2(k,j,i),3,1), KIND = wp ) 5232 ibit34 = REAL( IBITS(advc_flags_2(k,j,i),2,1), KIND = wp ) 5233 ibit33 = REAL( IBITS(advc_flags_2(k,j,i),1,1), KIND = wp ) 5181 5234 5182 5235 k_ppp = k + 3 * ibit35
Note: See TracChangeset
for help on using the changeset viewer.
