- Timestamp:
- Jun 14, 2017 10:17:34 AM (7 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/pmc_interface_mod.f90 (modified) (21 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/pmc_interface_mod.f90
r2281 r2283 26 26 ! ----------------- 27 27 ! $Id$ 28 ! Bugfixes in inititalization of log-law correction concerning 29 ! new topography concept 30 ! 31 ! 2281 2017-06-13 11:34:50Z suehring 28 32 ! Bugfix, add pre-preprocessor directives to enable non-parrallel mode 29 33 ! … … 1320 1324 IMPLICIT NONE 1321 1325 1322 INTEGER(iwp) :: direction !: Wall normal index: 1=k, 2=j, 3=i. 1323 INTEGER(iwp) :: end_index !: End index of present surface data type 1324 INTEGER(iwp) :: i !: 1325 INTEGER(iwp) :: icorr !: 1326 INTEGER(iwp) :: inc !: Wall outward-normal index increment -1 1327 !: or 1, for direction=1, inc=1 always 1328 INTEGER(iwp) :: iw !: 1329 INTEGER(iwp) :: j !: 1330 INTEGER(iwp) :: jcorr !: 1331 INTEGER(iwp) :: jw !: 1332 INTEGER(iwp) :: k !: 1333 INTEGER(iwp) :: k_wall_u_ji !: 1334 INTEGER(iwp) :: k_wall_u_ji_p !: 1335 INTEGER(iwp) :: k_wall_u_ji_m !: 1336 INTEGER(iwp) :: k_wall_v_ji !: 1337 INTEGER(iwp) :: k_wall_v_ji_p !: 1338 INTEGER(iwp) :: k_wall_v_ji_m !: 1339 INTEGER(iwp) :: k_wall_w_ji !: 1340 INTEGER(iwp) :: k_wall_w_ji_p !: 1341 INTEGER(iwp) :: k_wall_w_ji_m !: 1342 INTEGER(iwp) :: kb !: 1343 INTEGER(iwp) :: kcorr !: 1344 INTEGER(iwp) :: lc !: 1345 INTEGER(iwp) :: m !: Running index for surface data type 1346 INTEGER(iwp) :: ni !: 1347 INTEGER(iwp) :: nj !: 1348 INTEGER(iwp) :: nk !: 1349 INTEGER(iwp) :: nzt_topo_max !: 1350 INTEGER(iwp) :: start_index !: Start index of present surface data type 1351 INTEGER(iwp) :: wall_index !: Index of the wall-node coordinate 1352 1353 REAL(wp) :: z0_topo !: roughness at vertical walls 1354 REAL(wp), ALLOCATABLE, DIMENSION(:) :: lcr !: 1326 INTEGER(iwp) :: direction !< Wall normal index: 1=k, 2=j, 3=i. 1327 INTEGER(iwp) :: i !< 1328 INTEGER(iwp) :: icorr !< 1329 INTEGER(iwp) :: inc !< Wall outward-normal index increment -1 1330 !< or 1, for direction=1, inc=1 always 1331 INTEGER(iwp) :: iw !< 1332 INTEGER(iwp) :: j !< 1333 INTEGER(iwp) :: jcorr !< 1334 INTEGER(iwp) :: jw !< 1335 INTEGER(iwp) :: k !< 1336 INTEGER(iwp) :: k_wall_u_ji !< topography top index on u-grid 1337 INTEGER(iwp) :: k_wall_u_ji_p !< topography top index on u-grid 1338 INTEGER(iwp) :: k_wall_u_ji_m !< topography top index on u-grid 1339 INTEGER(iwp) :: k_wall_v_ji !< topography top index on v-grid 1340 INTEGER(iwp) :: k_wall_v_ji_p !< topography top index on v-grid 1341 INTEGER(iwp) :: k_wall_v_ji_m !< topography top index on v-grid 1342 INTEGER(iwp) :: k_wall_w_ji !< topography top index on w-grid 1343 INTEGER(iwp) :: k_wall_w_ji_p !< topography top index on w-grid 1344 INTEGER(iwp) :: k_wall_w_ji_m !< topography top index on w-grid 1345 INTEGER(iwp) :: kb !< 1346 INTEGER(iwp) :: kcorr !< 1347 INTEGER(iwp) :: lc !< 1348 INTEGER(iwp) :: m !< Running index for surface data type 1349 INTEGER(iwp) :: ni !< 1350 INTEGER(iwp) :: nj !< 1351 INTEGER(iwp) :: nk !< 1352 INTEGER(iwp) :: nzt_topo_max !< 1353 INTEGER(iwp) :: wall_index !< Index of the wall-node coordinate 1354 1355 REAL(wp) :: z0_topo !< roughness at vertical walls 1356 REAL(wp), ALLOCATABLE, DIMENSION(:) :: lcr !< 1355 1357 1356 1358 ! … … 1452 1454 MAXLOC( & 1453 1455 MERGE( 1, 0, & 1454 BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 ) &1456 BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 ) & 1455 1457 ), DIM = 1 & 1456 1458 ) - 1 ) … … 1460 1462 MAXLOC( & 1461 1463 MERGE( 1, 0, & 1462 BTEST( wall_flags_0(nzb:nzb_max,j,i), 14 ) &1464 BTEST( wall_flags_0(nzb:nzb_max,j,i), 14 ) & 1463 1465 ), DIM = 1 & 1464 1466 ) - 1 ) … … 1468 1470 MAXLOC( & 1469 1471 MERGE( 1, 0, & 1470 BTEST( wall_flags_0(nzb:nzb_max,j,i), 16 ) &1472 BTEST( wall_flags_0(nzb:nzb_max,j,i), 16 ) & 1471 1473 ), DIM = 1 & 1472 1474 ) - 1 ) … … 1476 1478 MAXLOC( & 1477 1479 MERGE( 1, 0, & 1478 BTEST( wall_flags_0(nzb:nzb_max,j,i), 18 ) &1480 BTEST( wall_flags_0(nzb:nzb_max,j,i), 18 ) & 1479 1481 ), DIM = 1 & 1480 1482 ) - 1 ) … … 1494 1496 MAXLOC( & 1495 1497 MERGE( 1, 0, & 1496 BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 ) &1498 BTEST( wall_flags_0(nzb:nzb_max,j,i), 12 ) & 1497 1499 ), DIM = 1 & 1498 1500 ) - 1 ) … … 1502 1504 MAXLOC( & 1503 1505 MERGE( 1, 0, & 1504 BTEST( wall_flags_0(nzb:nzb_max,j,i), 14 ) &1506 BTEST( wall_flags_0(nzb:nzb_max,j,i), 14 ) & 1505 1507 ), DIM = 1 & 1506 1508 ) - 1 ) … … 1510 1512 MAXLOC( & 1511 1513 MERGE( 1, 0, & 1512 BTEST( wall_flags_0(nzb:nzb_max,j,i), 16 ) &1514 BTEST( wall_flags_0(nzb:nzb_max,j,i), 16 ) & 1513 1515 ), DIM = 1 & 1514 1516 ) - 1 ) … … 1518 1520 MAXLOC( & 1519 1521 MERGE( 1, 0, & 1520 BTEST( wall_flags_0(nzb:nzb_max,j,i), 18 ) &1522 BTEST( wall_flags_0(nzb:nzb_max,j,i), 18 ) & 1521 1523 ), DIM = 1 & 1522 1524 ) - 1 ) … … 1544 1546 ALLOCATE( lcr(0:ncorr-1) ) 1545 1547 lcr = 1.0_wp 1548 1549 z0_topo = roughness_length 1546 1550 1547 1551 ! … … 1572 1576 ! 1573 1577 !-- For loglaw correction the roughness z0 is required. z0, however, 1574 !-- is part of the surfacetypes now, so call subroutine according 1575 !-- to the present surface tpye. 1576 !-- Default surface type 1577 IF ( surf_def_h(0)%start_index(j,i) <= & 1578 surf_def_h(0)%end_index(j,i) ) THEN 1579 start_index = surf_def_h(0)%start_index(j,i) 1580 end_index = surf_def_h(0)%end_index(j,i) 1581 DO m = start_index, end_index 1582 k = surf_def_h(0)%k(m) 1583 wall_index = k - 1 1584 kb = k - 1 1585 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, & 1586 j, inc, wall_index, surf_def_h(0)%z0(m), & 1578 !-- is part of the surfacetypes now. Set default roughness instead. 1579 !-- Determine topography top index on u-grid 1580 kb = MAXLOC( MERGE( 1, 0, & 1581 BTEST( wall_flags_0(nzb:nzt_topo_nestbc_l,j,i), 14 ) & 1582 ), DIM = 1 & 1583 ) - 1 1584 k = kb + 1 1585 wall_index = kb 1586 1587 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, & 1588 j, inc, wall_index, z0_topo, & 1587 1589 kb, direction, ncorr ) 1588 ENDDO 1589 ! 1590 !-- Natural surface type 1591 ELSEIF ( surf_lsm_h%start_index(j,i) <= & 1592 surf_lsm_h%end_index(j,i) ) THEN 1593 start_index = surf_lsm_h%start_index(j,i) 1594 end_index = surf_lsm_h%end_index(j,i) 1595 DO m = start_index, end_index 1596 k = surf_lsm_h%k(m) 1597 wall_index = k - 1 1598 kb = k - 1 1599 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, & 1600 j, inc, wall_index, surf_lsm_h%z0(m), & 1601 kb, direction, ncorr ) 1602 ENDDO 1603 ! 1604 !-- Urban surface type 1605 ELSEIF ( surf_usm_h%start_index(j,i) <= & 1606 surf_usm_h%end_index(j,i) ) THEN 1607 start_index = surf_usm_h%start_index(j,i) 1608 end_index = surf_usm_h%end_index(j,i) 1609 DO m = start_index, end_index 1610 k = surf_usm_h%k(m) 1611 wall_index = k - 1 1612 kb = k - 1 1613 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, & 1614 j, inc, wall_index, surf_usm_h%z0(m), & 1615 kb, direction, ncorr ) 1616 ENDDO 1617 ENDIF 1590 1618 1591 logc_u_l(1,k,j) = lc 1619 1592 logc_ratio_u_l(1,0:ncorr-1,k,j) = lcr(0:ncorr-1) … … 1623 1596 i = -1 1624 1597 ! 1625 !-- For loglaw correction the roughness z0 is required. z0, however, 1626 !-- is part of the surfacetypes now, so call subroutine according 1627 !-- to the present surface tpye. 1628 !-- Default surface type 1629 IF ( surf_def_h(0)%start_index(j,i) <= & 1630 surf_def_h(0)%end_index(j,i) ) THEN 1631 start_index = surf_def_h(0)%start_index(j,i) 1632 end_index = surf_def_h(0)%end_index(j,i) 1633 DO m = start_index, end_index 1634 k = surf_def_h(0)%k(m) 1635 wall_index = k - 1 1636 kb = k - 1 1637 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, & 1638 j, inc, wall_index, surf_def_h(0)%z0(m), & 1598 !-- Determine topography top index on v-grid 1599 kb = MAXLOC( MERGE( 1, 0, & 1600 BTEST( wall_flags_0(nzb:nzt_topo_nestbc_l,j,i), 16 ) & 1601 ), DIM = 1 & 1602 ) - 1 1603 k = kb + 1 1604 wall_index = kb 1605 1606 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, & 1607 j, inc, wall_index, z0_topo, & 1639 1608 kb, direction, ncorr ) 1640 ENDDO 1641 ! 1642 !-- Natural surface type 1643 ELSEIF ( surf_lsm_h%start_index(j,i) <= & 1644 surf_lsm_h%end_index(j,i) ) THEN 1645 start_index = surf_lsm_h%start_index(j,i) 1646 end_index = surf_lsm_h%end_index(j,i) 1647 DO m = start_index, end_index 1648 k = surf_lsm_h%k(m) 1649 wall_index = k - 1 1650 kb = k - 1 1651 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, & 1652 j, inc, wall_index, surf_lsm_h%z0(m), & 1653 kb, direction, ncorr ) 1654 ENDDO 1655 ! 1656 !-- Urban surface type 1657 ELSEIF ( surf_usm_h%start_index(j,i) <= & 1658 surf_usm_h%end_index(j,i) ) THEN 1659 start_index = surf_usm_h%start_index(j,i) 1660 end_index = surf_usm_h%end_index(j,i) 1661 DO m = start_index, end_index 1662 k = surf_usm_h%k(m) 1663 wall_index = k - 1 1664 kb = k - 1 1665 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, & 1666 j, inc, wall_index, surf_usm_h%z0(m), & 1667 kb, direction, ncorr ) 1668 ENDDO 1669 ENDIF 1609 1670 1610 logc_v_l(1,k,j) = lc 1671 1611 logc_ratio_v_l(1,0:ncorr-1,k,j) = lcr(0:ncorr-1) … … 1703 1643 !-- is part of the surfacetypes now, so call subroutine according 1704 1644 !-- to the present surface tpye. 1705 !-- Default surface type 1706 IF ( surf_def_h(0)%start_index(j,i) <= & 1707 surf_def_h(0)%end_index(j,i) ) THEN 1708 start_index = surf_def_h(0)%start_index(j,i) 1709 end_index = surf_def_h(0)%end_index(j,i) 1710 DO m = start_index, end_index 1711 k = surf_def_h(0)%k(m) 1712 wall_index = k - 1 1713 kb = k - 1 1714 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, & 1715 j, inc, wall_index, surf_def_h(0)%z0(m), & 1716 kb, direction, ncorr ) 1717 ENDDO 1718 ! 1719 !-- Natural surface type 1720 ELSEIF ( surf_lsm_h%start_index(j,i) <= & 1721 surf_lsm_h%end_index(j,i) ) THEN 1722 start_index = surf_lsm_h%start_index(j,i) 1723 end_index = surf_lsm_h%end_index(j,i) 1724 DO m = start_index, end_index 1725 k = surf_lsm_h%k(m) 1726 wall_index = k - 1 1727 kb = k - 1 1728 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, & 1729 j, inc, wall_index, surf_lsm_h%z0(m), & 1730 kb, direction, ncorr ) 1731 ENDDO 1732 ! 1733 !-- Urban surface type 1734 ELSEIF ( surf_usm_h%start_index(j,i) <= & 1735 surf_usm_h%end_index(j,i) ) THEN 1736 start_index = surf_usm_h%start_index(j,i) 1737 end_index = surf_usm_h%end_index(j,i) 1738 DO m = start_index, end_index 1739 k = surf_usm_h%k(m) 1740 wall_index = k - 1 1741 kb = k - 1 1742 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, & 1743 j, inc, wall_index, surf_usm_h%z0(m), & 1744 kb, direction, ncorr ) 1745 ENDDO 1746 ENDIF 1645 !-- Determine topography top index on u-grid 1646 kb = MAXLOC( MERGE( 1, 0, & 1647 BTEST( wall_flags_0(nzb:nzt_topo_nestbc_l,j,i), 14 ) & 1648 ), DIM = 1 & 1649 ) - 1 1650 k = kb + 1 1651 wall_index = kb 1652 1653 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, & 1654 j, inc, wall_index, z0_topo, & 1655 kb, direction, ncorr ) 1747 1656 1748 1657 logc_u_r(1,k,j) = lc … … 1753 1662 i = nxr + 1 1754 1663 ! 1755 !-- For loglaw correction the roughness z0 is required. z0, however, 1756 !-- is part of the surfacetypes now, so call subroutine according 1757 !-- to the present surface tpye. 1758 !-- Default surface type 1759 IF ( surf_def_h(0)%start_index(j,i) <= & 1760 surf_def_h(0)%end_index(j,i) ) THEN 1761 start_index = surf_def_h(0)%start_index(j,i) 1762 end_index = surf_def_h(0)%end_index(j,i) 1763 DO m = start_index, end_index 1764 k = surf_def_h(0)%k(m) 1765 wall_index = k - 1 1766 kb = k - 1 1767 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, & 1768 j, inc, wall_index, surf_def_h(0)%z0(m), & 1769 kb, direction, ncorr ) 1770 ENDDO 1771 ! 1772 !-- Natural surface type 1773 ELSEIF ( surf_lsm_h%start_index(j,i) <= & 1774 surf_lsm_h%end_index(j,i) ) THEN 1775 start_index = surf_lsm_h%start_index(j,i) 1776 end_index = surf_lsm_h%end_index(j,i) 1777 DO m = start_index, end_index 1778 k = surf_lsm_h%k(m) 1779 wall_index = k - 1 1780 kb = k - 1 1781 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, & 1782 j, inc, wall_index, surf_lsm_h%z0(m), & 1783 kb, direction, ncorr ) 1784 ENDDO 1785 ! 1786 !-- Urban surface type 1787 ELSEIF ( surf_usm_h%start_index(j,i) <= & 1788 surf_usm_h%end_index(j,i) ) THEN 1789 start_index = surf_usm_h%start_index(j,i) 1790 end_index = surf_usm_h%end_index(j,i) 1791 DO m = start_index, end_index 1792 k = surf_usm_h%k(m) 1793 wall_index = k - 1 1794 kb = k - 1 1795 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, & 1796 j, inc, wall_index, surf_usm_h%z0(m), & 1797 kb, direction, ncorr ) 1798 ENDDO 1799 ENDIF 1664 !-- Determine topography top index on v-grid 1665 kb = MAXLOC( MERGE( 1, 0, & 1666 BTEST( wall_flags_0(nzb:nzt_topo_nestbc_l,j,i), 16 ) & 1667 ), DIM = 1 & 1668 ) - 1 1669 k = kb + 1 1670 wall_index = kb 1671 1672 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, & 1673 j, inc, wall_index, z0_topo, & 1674 kb, direction, ncorr ) 1675 1800 1676 logc_v_r(1,k,j) = lc 1801 1677 logc_ratio_v_r(1,0:ncorr-1,k,j) = lcr(0:ncorr-1) … … 1830 1706 j = -1 1831 1707 ! 1832 !-- For loglaw correction the roughness z0 is required. z0, however, 1833 !-- is part of the surfacetypes now, so call subroutine according 1834 !-- to the present surface tpye. 1835 !-- Default surface type 1836 IF ( surf_def_h(0)%start_index(j,i) <= & 1837 surf_def_h(0)%end_index(j,i) ) THEN 1838 start_index = surf_def_h(0)%start_index(j,i) 1839 end_index = surf_def_h(0)%end_index(j,i) 1840 DO m = start_index, end_index 1841 k = surf_def_h(0)%k(m) 1842 wall_index = k - 1 1843 kb = k - 1 1844 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, & 1845 j, inc, wall_index, surf_def_h(0)%z0(m), & 1846 kb, direction, ncorr ) 1847 ENDDO 1848 ! 1849 !-- Natural surface type 1850 ELSEIF ( surf_lsm_h%start_index(j,i) <= & 1851 surf_lsm_h%end_index(j,i) ) THEN 1852 start_index = surf_lsm_h%start_index(j,i) 1853 end_index = surf_lsm_h%end_index(j,i) 1854 DO m = start_index, end_index 1855 k = surf_lsm_h%k(m) 1856 wall_index = k - 1 1857 kb = k - 1 1858 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, & 1859 j, inc, wall_index, surf_lsm_h%z0(m), & 1860 kb, direction, ncorr ) 1861 ENDDO 1862 ! 1863 !-- Urban surface type 1864 ELSEIF ( surf_usm_h%start_index(j,i) <= & 1865 surf_usm_h%end_index(j,i) ) THEN 1866 start_index = surf_usm_h%start_index(j,i) 1867 end_index = surf_usm_h%end_index(j,i) 1868 DO m = start_index, end_index 1869 k = surf_usm_h%k(m) 1870 wall_index = k - 1 1871 kb = k - 1 1872 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, & 1873 j, inc, wall_index, surf_usm_h%z0(m), & 1874 kb, direction, ncorr ) 1875 ENDDO 1876 ENDIF 1708 !-- Determine topography top index on u-grid 1709 kb = MAXLOC( MERGE( 1, 0, & 1710 BTEST( wall_flags_0(nzb:nzt_topo_nestbc_l,j,i), 14 ) & 1711 ), DIM = 1 & 1712 ) - 1 1713 k = kb + 1 1714 wall_index = kb 1715 1716 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, & 1717 j, inc, wall_index, z0_topo, & 1718 kb, direction, ncorr ) 1719 1877 1720 logc_u_s(1,k,i) = lc 1878 1721 logc_ratio_u_s(1,0:ncorr-1,k,i) = lcr(0:ncorr-1) … … 1882 1725 j = 0 1883 1726 ! 1884 !-- For loglaw correction the roughness z0 is required. z0, however, 1885 !-- is part of the surfacetypes now, so call subroutine according 1886 !-- to the present surface tpye. 1887 !-- Default surface type 1888 IF ( surf_def_h(0)%start_index(j,i) <= & 1889 surf_def_h(0)%end_index(j,i) ) THEN 1890 start_index = surf_def_h(0)%start_index(j,i) 1891 end_index = surf_def_h(0)%end_index(j,i) 1892 DO m = start_index, end_index 1893 k = surf_def_h(0)%k(m) 1894 wall_index = k - 1 1895 kb = k - 1 1896 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, & 1897 j, inc, wall_index, surf_def_h(0)%z0(m), & 1898 kb, direction, ncorr ) 1899 ENDDO 1900 ! 1901 !-- Natural surface type 1902 ELSEIF ( surf_lsm_h%start_index(j,i) <= & 1903 surf_lsm_h%end_index(j,i) ) THEN 1904 start_index = surf_lsm_h%start_index(j,i) 1905 end_index = surf_lsm_h%end_index(j,i) 1906 DO m = start_index, end_index 1907 k = surf_lsm_h%k(m) 1908 wall_index = k - 1 1909 kb = k - 1 1910 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, & 1911 j, inc, wall_index, surf_lsm_h%z0(m), & 1912 kb, direction, ncorr ) 1913 ENDDO 1914 ! 1915 !-- Urban surface type 1916 ELSEIF ( surf_usm_h%start_index(j,i) <= & 1917 surf_usm_h%end_index(j,i) ) THEN 1918 start_index = surf_usm_h%start_index(j,i) 1919 end_index = surf_usm_h%end_index(j,i) 1920 DO m = start_index, end_index 1921 k = surf_usm_h%k(m) 1922 wall_index = k - 1 1923 kb = k - 1 1924 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, & 1925 j, inc, wall_index, surf_usm_h%z0(m), & 1926 kb, direction, ncorr ) 1927 ENDDO 1928 ENDIF 1727 !-- Determine topography top index on v-grid 1728 kb = MAXLOC( MERGE( 1, 0, & 1729 BTEST( wall_flags_0(nzb:nzt_topo_nestbc_l,j,i), 16 ) & 1730 ), DIM = 1 & 1731 ) - 1 1732 k = kb + 1 1733 wall_index = kb 1734 1735 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, & 1736 j, inc, wall_index, z0_topo, & 1737 kb, direction, ncorr ) 1738 1929 1739 logc_v_s(1,k,i) = lc 1930 1740 logc_ratio_v_s(1,0:ncorr-1,k,i) = lcr(0:ncorr-1) … … 1959 1769 j = nyn + 1 1960 1770 ! 1961 !-- For loglaw correction the roughness z0 is required. z0, however, 1962 !-- is part of the surfacetypes now, so call subroutine according 1963 !-- to the present surface tpye. 1964 !-- Default surface type 1965 IF ( surf_def_h(0)%start_index(j,i) <= & 1966 surf_def_h(0)%end_index(j,i) ) THEN 1967 start_index = surf_def_h(0)%start_index(j,i) 1968 end_index = surf_def_h(0)%end_index(j,i) 1969 DO m = start_index, end_index 1970 k = surf_def_h(0)%k(m) 1971 wall_index = k - 1 1972 kb = k - 1 1973 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, & 1974 j, inc, wall_index, surf_def_h(0)%z0(m), & 1975 kb, direction, ncorr ) 1976 ENDDO 1977 ! 1978 !-- Natural surface type 1979 ELSEIF ( surf_lsm_h%start_index(j,i) <= & 1980 surf_lsm_h%end_index(j,i) ) THEN 1981 start_index = surf_lsm_h%start_index(j,i) 1982 end_index = surf_lsm_h%end_index(j,i) 1983 DO m = start_index, end_index 1984 k = surf_lsm_h%k(m) 1985 wall_index = k - 1 1986 kb = k - 1 1987 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, & 1988 j, inc, wall_index, surf_lsm_h%z0(m), & 1989 kb, direction, ncorr ) 1990 ENDDO 1991 ! 1992 !-- Urban surface type 1993 ELSEIF ( surf_usm_h%start_index(j,i) <= & 1994 surf_usm_h%end_index(j,i) ) THEN 1995 start_index = surf_usm_h%start_index(j,i) 1996 end_index = surf_usm_h%end_index(j,i) 1997 DO m = start_index, end_index 1998 k = surf_usm_h%k(m) 1999 wall_index = k - 1 2000 kb = k - 1 2001 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, & 2002 j, inc, wall_index, surf_usm_h%z0(m), & 2003 kb, direction, ncorr ) 2004 ENDDO 2005 ENDIF 1771 !-- Determine topography top index on v-grid 1772 kb = MAXLOC( MERGE( 1, 0, & 1773 BTEST( wall_flags_0(nzb:nzt_topo_nestbc_l,j,i), 14 ) & 1774 ), DIM = 1 & 1775 ) - 1 1776 k = kb + 1 1777 wall_index = kb 1778 1779 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, & 1780 j, inc, wall_index, z0_topo, & 1781 kb, direction, ncorr ) 1782 2006 1783 logc_u_n(1,k,i) = lc 2007 1784 logc_ratio_u_n(1,0:ncorr-1,k,i) = lcr(0:ncorr-1) … … 2011 1788 j = nyn + 1 2012 1789 ! 2013 !-- For loglaw correction the roughness z0 is required. z0, however, 2014 !-- is part of the surfacetypes now, so call subroutine according 2015 !-- to the present surface tpye. 2016 !-- Default surface type 2017 IF ( surf_def_h(0)%start_index(j,i) <= & 2018 surf_def_h(0)%end_index(j,i) ) THEN 2019 start_index = surf_def_h(0)%start_index(j,i) 2020 end_index = surf_def_h(0)%end_index(j,i) 2021 DO m = start_index, end_index 2022 k = surf_def_h(0)%k(m) 2023 wall_index = k - 1 2024 kb = k - 1 2025 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, & 2026 j, inc, wall_index, surf_def_h(0)%z0(m), & 2027 kb, direction, ncorr ) 2028 ENDDO 2029 ! 2030 !-- Natural surface type 2031 ELSEIF ( surf_lsm_h%start_index(j,i) <= & 2032 surf_lsm_h%end_index(j,i) ) THEN 2033 start_index = surf_lsm_h%start_index(j,i) 2034 end_index = surf_lsm_h%end_index(j,i) 2035 DO m = start_index, end_index 2036 k = surf_lsm_h%k(m) 2037 wall_index = k - 1 2038 kb = k - 1 2039 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, & 2040 j, inc, wall_index, surf_lsm_h%z0(m), & 2041 kb, direction, ncorr ) 2042 ENDDO 2043 ! 2044 !-- Urban surface type 2045 ELSEIF ( surf_usm_h%start_index(j,i) <= & 2046 surf_usm_h%end_index(j,i) ) THEN 2047 start_index = surf_usm_h%start_index(j,i) 2048 end_index = surf_usm_h%end_index(j,i) 2049 DO m = start_index, end_index 2050 k = surf_usm_h%k(m) 2051 wall_index = k - 1 2052 kb = k - 1 2053 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, & 2054 j, inc, wall_index, surf_usm_h%z0(m), & 2055 kb, direction, ncorr ) 2056 ENDDO 2057 ENDIF 1790 !-- Determine topography top index on v-grid 1791 kb = MAXLOC( MERGE( 1, 0, & 1792 BTEST( wall_flags_0(nzb:nzt_topo_nestbc_l,j,i), 16 ) & 1793 ), DIM = 1 & 1794 ) - 1 1795 k = kb + 1 1796 wall_index = kb 1797 1798 CALL pmci_define_loglaw_correction_parameters( lc, lcr, k, & 1799 j, inc, wall_index, z0_topo, & 1800 kb, direction, ncorr ) 2058 1801 logc_v_n(1,k,i) = lc 2059 1802 logc_ratio_v_n(1,0:ncorr-1,k,i) = lcr(0:ncorr-1) … … 5671 5414 !-- For simplicity anterpolate within buildings and under elevated 5672 5415 !-- terrain too 5673 DO kk = kcb, kct -15416 DO kk = kcb, kct 5674 5417 ! 5675 5418 !-- ijfc and kfc are precomputed in pmci_init_anterp_tophat … … 5690 5433 5691 5434 ENDDO 5692 fc(kct,jj,ii) = fc(kct,jj,ii)5693 5435 5694 5436 ENDDO
Note: See TracChangeset
for help on using the changeset viewer.
