Changeset 2283


Ignore:
Timestamp:
Jun 14, 2017 10:17:34 AM (7 years ago)
Author:
suehring
Message:

Bugfixes in inititalization of log-law correction concerning new topography concept

File:
1 edited

Legend:

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

    r2281 r2283  
    2626! -----------------
    2727! $Id$
     28! Bugfixes in inititalization of log-law correction concerning
     29! new topography concept
     30!
     31! 2281 2017-06-13 11:34:50Z suehring
    2832! Bugfix, add pre-preprocessor directives to enable non-parrallel mode
    2933!
     
    13201324       IMPLICIT NONE
    13211325
    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   !<
    13551357
    13561358!
     
    14521454                  MAXLOC(                                                      &
    14531455                          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 )    &
    14551457                               ), DIM = 1                                      &
    14561458                        ) - 1          )
     
    14601462                  MAXLOC(                                                      &
    14611463                          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 )    &
    14631465                               ), DIM = 1                                      &
    14641466                        ) - 1          )
     
    14681470                  MAXLOC(                                                      &
    14691471                          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 )    &
    14711473                               ), DIM = 1                                      &
    14721474                        ) - 1          )
     
    14761478                  MAXLOC(                                                      &
    14771479                          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 )    &
    14791481                               ), DIM = 1                                      &
    14801482                        ) - 1          )
     
    14941496                  MAXLOC(                                                      &
    14951497                          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 )    &
    14971499                               ), DIM = 1                                      &
    14981500                        ) - 1          )
     
    15021504                  MAXLOC(                                                      &
    15031505                          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 )    &
    15051507                               ), DIM = 1                                      &
    15061508                        ) - 1          )
     
    15101512                  MAXLOC(                                                      &
    15111513                          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 )    &
    15131515                               ), DIM = 1                                      &
    15141516                        ) - 1          )
     
    15181520                  MAXLOC(                                                      &
    15191521                          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 )    &
    15211523                               ), DIM = 1                                      &
    15221524                        ) - 1          )
     
    15441546       ALLOCATE( lcr(0:ncorr-1) )
    15451547       lcr = 1.0_wp
     1548
     1549       z0_topo = roughness_length
    15461550
    15471551!
     
    15721576!
    15731577!--          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,                       &
    15871589                            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
    16181591             logc_u_l(1,k,j) = lc
    16191592             logc_ratio_u_l(1,0:ncorr-1,k,j) = lcr(0:ncorr-1)
     
    16231596             i   = -1
    16241597!
    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,                       &
    16391608                            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
    16701610             logc_v_l(1,k,j) = lc
    16711611             logc_ratio_v_l(1,0:ncorr-1,k,j) = lcr(0:ncorr-1)
     
    17031643!--          is part of the surfacetypes now, so call subroutine according
    17041644!--          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 )
    17471656
    17481657             logc_u_r(1,k,j) = lc
     
    17531662             i   = nxr + 1
    17541663!
    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
    18001676             logc_v_r(1,k,j) = lc
    18011677             logc_ratio_v_r(1,0:ncorr-1,k,j) = lcr(0:ncorr-1)
     
    18301706             j   = -1
    18311707!
    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
    18771720             logc_u_s(1,k,i) = lc
    18781721             logc_ratio_u_s(1,0:ncorr-1,k,i) = lcr(0:ncorr-1)
     
    18821725             j   = 0
    18831726!
    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
    19291739             logc_v_s(1,k,i) = lc
    19301740             logc_ratio_v_s(1,0:ncorr-1,k,i) = lcr(0:ncorr-1)
     
    19591769             j   = nyn + 1
    19601770!
    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
    20061783             logc_u_n(1,k,i) = lc
    20071784             logc_ratio_u_n(1,0:ncorr-1,k,i) = lcr(0:ncorr-1)
     
    20111788             j   = nyn + 1
    20121789!
    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 )
    20581801             logc_v_n(1,k,i) = lc
    20591802             logc_ratio_v_n(1,0:ncorr-1,k,i) = lcr(0:ncorr-1)
     
    56715414!--          For simplicity anterpolate within buildings and under elevated
    56725415!--          terrain too
    5673              DO  kk = kcb, kct-1
     5416             DO  kk = kcb, kct
    56745417!
    56755418!--             ijfc and kfc are precomputed in pmci_init_anterp_tophat
     
    56905433
    56915434             ENDDO
    5692              fc(kct,jj,ii) = fc(kct,jj,ii)
    56935435
    56945436          ENDDO
Note: See TracChangeset for help on using the changeset viewer.