Changeset 4510 for palm/trunk/SOURCE/urban_surface_mod.f90
- Timestamp:
- Apr 29, 2020 2:19:18 PM (4 years ago)
- File:
-
- 1 edited
-
palm/trunk/SOURCE/urban_surface_mod.f90 (modified) (205 diffs)
Legend:
- Unmodified
- Added
- Removed
-
palm/trunk/SOURCE/urban_surface_mod.f90
r4509 r4510 18 18 ! Copyright 1997-2020 Leibniz Universitaet Hannover 19 19 !--------------------------------------------------------------------------------------------------! 20 ! 20 21 ! 21 22 ! Current revisions: … … 26 27 ! ----------------- 27 28 ! $Id$ 28 ! file re-formatted to follow the PALM coding standard 29 ! Further re-formatting to follow the PALM coding standard 30 ! 31 ! 4509 2020-04-26 15:57:55Z raasch 32 ! File re-formatted to follow the PALM coding standard 29 33 ! 30 34 ! 4500 2020-04-17 10:12:45Z suehring 31 35 ! Allocate array for wall heat flux, which is further used to aggregate tile 32 36 ! fractions in the surface output 33 ! 37 ! 34 38 ! 4495 2020-04-13 20:11:20Z raasch 35 ! restart data handling with MPI-IO added39 ! Restart data handling with MPI-IO added 36 40 ! 37 41 ! 4493 2020-04-10 09:49:43Z pavelkrc 38 42 ! J.Resler, 2020/03/19 39 ! - remove reading of deprecated input parameters c_surface and lambda_surf40 ! - and calculate them from parameters of the outer wall/roof layer43 ! - Remove reading of deprecated input parameters c_surface and lambda_surf 44 ! - And calculate them from parameters of the outer wall/roof layer 41 45 ! 42 46 ! 4481 2020-03-31 18:55:54Z maronga 43 ! use statement for exchange horiz added47 ! Use statement for exchange horiz added 44 48 ! 45 49 ! 4442 2020-03-04 19:21:13Z suehring … … 80 84 ! 81 85 ! 4227 2019-09-10 18:04:34Z gronemeier 82 ! implement new palm_date_time_mod86 ! Implement new palm_date_time_mod 83 87 ! 84 88 ! 4214 2019-09-02 15:57:02Z suehring … … 152 156 ! 153 157 ! 3832 2019-03-28 13:16:58Z raasch 154 ! instrumented with openmp directives158 ! Instrumented with openmp directives 155 159 ! 156 160 ! 3824 2019-03-27 15:56:16Z pavelkrc … … 159 163 ! 160 164 ! 3814 2019-03-26 08:40:31Z pavelkrc 161 ! unused subroutine commented out165 ! Unused subroutine commented out 162 166 ! 163 167 ! 3769 2019-02-28 10:16:49Z moh.hefny 164 ! removed unused variables168 ! Removed unused variables 165 169 ! 166 170 ! 3767 2019-02-27 08:18:02Z raasch 167 ! unused variables removed from rrd-subroutines parameter list171 ! Unused variables removed from rrd-subroutines parameter list 168 172 ! 169 173 ! 3748 2019-02-18 10:38:31Z suehring … … 172 176 ! 3745 2019-02-15 18:57:56Z suehring 173 177 ! - Remove internal flag indoor_model (is a global control parameter) 174 ! - add waste heat from buildings to the kinmatic heat flux175 ! - consider waste heat in restart data176 ! - remove unused USE statements178 ! - Add waste heat from buildings to the kinmatic heat flux 179 ! - Consider waste heat in restart data 180 ! - Remove unused USE statements 177 181 ! 178 182 ! 3744 2019-02-15 18:38:58Z suehring 179 ! fixed surface heat capacity in the building parameters convert the file back to unix format183 ! Fixed surface heat capacity in the building parameters convert the file back to unix format 180 184 ! 181 185 ! 3730 2019-02-11 11:26:47Z moh.hefny … … 186 190 ! 187 191 ! 3705 2019-01-29 19:56:39Z suehring 188 ! make nzb_wall public, required for virtual-measurements192 ! Make nzb_wall public, required for virtual-measurements 189 193 ! 190 194 ! 3704 2019-01-29 19:51:41Z suehring … … 208 212 !> Module for Urban Surface Model (USM) 209 213 !> The module includes: 210 !> 1. radiation model with direct/diffuse radiation, shading, reflections and integration with214 !> 1. Radiation model with direct/diffuse radiation, shading, reflections and integration with 211 215 !> plant canopy 212 !> 2. wall and wall surface model213 !> 3. surface layer energy balance214 !> 4. anthropogenic heat (only from transportation so far)215 !> 5. necessary auxiliary subroutines (reading inputs, writing outputs, restart simulations, ...)216 !> 2. Wall and wall surface model 217 !> 3. Surface layer energy balance 218 !> 4. Anthropogenic heat (only from transportation so far) 219 !> 5. Necessary auxiliary subroutines (reading inputs, writing outputs, restart simulations, ...) 216 220 !> It also makes use of standard radiation and integrates it into urban surface model. 217 221 !> … … 223 227 !> fraq(0,m) + fraq(1,m) = 0?! 224 228 !> @todo Use unit 90 for OPEN/CLOSE of input files (FK) 225 !> @todo remove reading of old csv inputs229 !> @todo Remove reading of old csv inputs 226 230 !--------------------------------------------------------------------------------------------------! 227 231 MODULE urban_surface_mod … … 418 422 /), (/ 4, 7 /) ) 419 423 ! 420 !-- value 9999999.9_wp -> generic available or user-defined value must be set otherwise421 !-- -> no generic variable and user setting is optional424 !-- Value 9999999.9_wp -> Generic available or user-defined value must be set otherwise 425 !-- -> No generic variable and user setting is optional 422 426 REAL(wp) :: alpha_vangenuchten = 9999999.9_wp !< NAMELIST alpha_vg 423 427 REAL(wp) :: field_capacity = 9999999.9_wp !< NAMELIST fc … … 430 434 431 435 ! 432 !-- configuration parameters (they can be setup in PALM config)436 !-- Configuration parameters (they can be setup in PALM config) 433 437 LOGICAL :: force_radiation_call_l = .FALSE. !< flag parameter for unscheduled radiation model calls 434 438 LOGICAL :: read_wall_temp_3d = .FALSE. !< … … 645 649 TYPE(surf_type_usm), TARGET :: tm_liq_usm_h_m !< liquid water reservoir tendency (m), horizontal surface elements 646 650 ! 647 !-- anthropogenic heat sources651 !-- Anthropogenic heat sources 648 652 INTEGER(iwp) :: naheatlayers = 1 !< number of layers of anthropogenic heat 649 653 … … 653 657 654 658 ! 655 !-- wall surface model656 !-- wall surface model constants659 !-- Wall surface model 660 !-- Wall surface model constants 657 661 INTEGER(iwp), PARAMETER :: nzb_wall = 0 !< inner side of the wall model (to be switched) 658 662 INTEGER(iwp), PARAMETER :: nzt_wall = 3 !< outer side of the wall model (to be switched) … … 680 684 681 685 ! 682 !-- surface and material model variables for walls, ground, roofs686 !-- Surface and material model variables for walls, ground, roofs 683 687 REAL(wp), DIMENSION(:), ALLOCATABLE :: zwn !< normalized wall layer depths (m) 684 688 REAL(wp), DIMENSION(:), ALLOCATABLE :: zwn_green !< normalized green layer depths (m) … … 715 719 ! 716 720 !-- Energy balance variables 717 !-- parameters of the land, roof and wall surfaces721 !-- Parameters of the land, roof and wall surfaces 718 722 REAL(wp), DIMENSION(:,:), POINTER :: fc_h !< 719 723 REAL(wp), DIMENSION(:,:), POINTER :: rootfr_h !< … … 763 767 ! 764 768 !-- Surface and material parameter classes (surface_type) 765 !-- albedo, emissivity, lambda_surf, roughness, thickness, volumetric heat capacity, thermal conductivity769 !-- Albedo, emissivity, lambda_surf, roughness, thickness, volumetric heat capacity, thermal conductivity 766 770 CHARACTER(12), DIMENSION(:), ALLOCATABLE :: surface_type_names !< names of wall types (used only for reports) 767 771 … … 790 794 791 795 ! 792 !-- interfaces of subroutines accessed from outside of this module796 !-- Interfaces of subroutines accessed from outside of this module 793 797 INTERFACE usm_3d_data_averaging 794 798 MODULE PROCEDURE usm_3d_data_averaging … … 923 927 ! 924 928 !-- Wall surface model 925 !-- allocate arrays for wall surface model and define pointers926 !-- allocate array of wall types and wall parameters929 !-- Allocate arrays for wall surface model and define pointers 930 !-- Allocate array of wall types and wall parameters 927 931 ALLOCATE ( surf_usm_h%surface_types(1:surf_usm_h%ns) ) 928 932 ALLOCATE ( surf_usm_h%building_type(1:surf_usm_h%ns) ) … … 984 988 985 989 ! 986 !-- wall and roof surface parameters. First for horizontal surfaces990 !-- Wall and roof surface parameters. First for horizontal surfaces 987 991 ALLOCATE ( surf_usm_h%isroof_surf(1:surf_usm_h%ns) ) 988 992 ALLOCATE ( surf_usm_h%lambda_surf(1:surf_usm_h%ns) ) … … 1020 1024 1021 1025 ! 1022 !-- allocate wall and roof material parameters. First for horizontal surfaces1026 !-- Allocate wall and roof material parameters. First for horizontal surfaces 1023 1027 ALLOCATE ( surf_usm_h%thickness_wall(1:surf_usm_h%ns) ) 1024 1028 ALLOCATE ( surf_usm_h%thickness_window(1:surf_usm_h%ns) ) … … 1055 1059 1056 1060 ! 1057 !-- allocate green wall and roof vegetation and soil parameters. First horizontal surfaces1061 !-- Allocate green wall and roof vegetation and soil parameters. First horizontal surfaces 1058 1062 ALLOCATE ( surf_usm_h%g_d(1:surf_usm_h%ns) ) 1059 1063 ALLOCATE ( surf_usm_h%c_liq(1:surf_usm_h%ns) ) … … 1077 1081 1078 1082 ! 1079 !-- allocate wall and roof layers sizes. For horizontal surfaces.1083 !-- Allocate wall and roof layers sizes. For horizontal surfaces. 1080 1084 ALLOCATE ( zwn(nzb_wall:nzt_wall) ) 1081 1085 ALLOCATE ( surf_usm_h%dz_wall(nzb_wall:nzt_wall+1,1:surf_usm_h%ns) ) … … 1326 1330 1327 1331 1328 IF ( variable(1:4) == 'usm_' ) THEN ! is such a check really required?1329 1330 ! 1331 !-- find the real name of the variable1332 IF ( variable(1:4) == 'usm_' ) THEN ! Is such a check really required? 1333 1334 ! 1335 !-- Find the real name of the variable 1332 1336 ids = -1 1333 1337 l = -1 … … 1343 1347 ENDIF 1344 1348 ENDDO 1345 l = idsint - 2 ! horizontal direction index - terrible hack !1349 l = idsint - 2 ! Horizontal direction index - terrible hack ! 1346 1350 IF ( l < 0 .OR. l > 3 ) THEN 1347 1351 l = -1 … … 1352 1356 IF ( var(1:11) == 'usm_t_wall_' .AND. len( TRIM( var ) ) >= 12 ) THEN 1353 1357 ! 1354 !-- wall layers1358 !-- Wall layers 1355 1359 READ( var(12:12), '(I1)', iostat=istat ) iwl 1356 1360 IF ( istat == 0 .AND. iwl >= nzb_wall .AND. iwl <= nzt_wall ) THEN … … 1358 1362 ELSE 1359 1363 ! 1360 !-- wrong wall layer index1364 !-- Wrong wall layer index 1361 1365 RETURN 1362 1366 ENDIF … … 1364 1368 IF ( var(1:13) == 'usm_t_window_' .AND. len( TRIM(var) ) >= 14 ) THEN 1365 1369 ! 1366 !-- wall layers1370 !-- Wall layers 1367 1371 READ( var(14:14), '(I1)', iostat=istat ) iwl 1368 1372 IF ( istat == 0 .AND. iwl >= nzb_wall .AND. iwl <= nzt_wall ) THEN … … 1370 1374 ELSE 1371 1375 ! 1372 !-- wrong window layer index1376 !-- Wrong window layer index 1373 1377 RETURN 1374 1378 ENDIF … … 1376 1380 IF ( var(1:12) == 'usm_t_green_' .AND. len( TRIM( var ) ) >= 13 ) THEN 1377 1381 ! 1378 !-- wall layers1382 !-- Wall layers 1379 1383 READ( var(13:13), '(I1)', iostat=istat ) iwl 1380 1384 IF ( istat == 0 .AND. iwl >= nzb_wall .AND. iwl <= nzt_wall ) THEN … … 1382 1386 ELSE 1383 1387 ! 1384 !-- wrong green layer index1388 !-- Wrong green layer index 1385 1389 RETURN 1386 1390 ENDIF … … 1388 1392 IF ( var(1:8) == 'usm_swc_' .AND. len( TRIM( var ) ) >= 9 ) THEN 1389 1393 ! 1390 !-- swc layers1394 !-- Swc layers 1391 1395 READ( var(9:9), '(I1)', iostat=istat ) iwl 1392 1396 IF ( istat == 0 .AND. iwl >= nzb_wall .AND. iwl <= nzt_wall ) THEN … … 1394 1398 ELSE 1395 1399 ! 1396 !-- wrong swc layer index1400 !-- Wrong swc layer index 1397 1401 RETURN 1398 1402 ENDIF … … 1405 1409 CASE ( 'usm_wshf' ) 1406 1410 ! 1407 !-- array of sensible heat flux from surfaces1408 !-- land surfaces1411 !-- Array of sensible heat flux from surfaces 1412 !-- Land surfaces 1409 1413 IF ( l == -1 ) THEN 1410 1414 IF ( .NOT. ALLOCATED( surf_usm_h%wshf_eb_av ) ) THEN … … 1421 1425 CASE ( 'usm_qsws' ) 1422 1426 ! 1423 !-- array of latent heat flux from surfaces1424 !-- land surfaces1427 !-- Array of latent heat flux from surfaces 1428 !-- Land surfaces 1425 1429 IF ( l == -1 .AND. .NOT. ALLOCATED( surf_usm_h%qsws_av ) ) THEN 1426 1430 ALLOCATE ( surf_usm_h%qsws_av(1:surf_usm_h%ns) ) … … 1435 1439 CASE ( 'usm_qsws_veg' ) 1436 1440 ! 1437 !-- array of latent heat flux from vegetation surfaces1438 !-- land surfaces1441 !-- Array of latent heat flux from vegetation surfaces 1442 !-- Land surfaces 1439 1443 IF ( l == -1 .AND. .NOT. ALLOCATED( surf_usm_h%qsws_veg_av ) ) THEN 1440 1444 ALLOCATE ( surf_usm_h%qsws_veg_av(1:surf_usm_h%ns) ) … … 1449 1453 CASE ( 'usm_qsws_liq' ) 1450 1454 ! 1451 !-- array of latent heat flux from surfaces with liquid1452 !-- land surfaces1455 !-- Array of latent heat flux from surfaces with liquid 1456 !-- Land surfaces 1453 1457 IF ( l == -1 .AND. .NOT. ALLOCATED( surf_usm_h%qsws_liq_av ) ) THEN 1454 1458 ALLOCATE ( surf_usm_h%qsws_liq_av(1:surf_usm_h%ns) ) … … 1466 1470 CASE ( 'usm_wghf' ) 1467 1471 ! 1468 !-- array of heat flux from ground (wall, roof, land)1472 !-- Array of heat flux from ground (wall, roof, land) 1469 1473 IF ( l == -1 ) THEN 1470 1474 IF ( .NOT. ALLOCATED( surf_usm_h%wghf_eb_av ) ) THEN … … 1481 1485 CASE ( 'usm_wghf_window' ) 1482 1486 ! 1483 !-- array of heat flux from window ground (wall, roof, land)1487 !-- Array of heat flux from window ground (wall, roof, land) 1484 1488 IF ( l == -1 ) THEN 1485 1489 IF ( .NOT. ALLOCATED( surf_usm_h%wghf_eb_window_av ) ) THEN … … 1496 1500 CASE ( 'usm_wghf_green' ) 1497 1501 ! 1498 !-- array of heat flux from green ground (wall, roof, land)1502 !-- Array of heat flux from green ground (wall, roof, land) 1499 1503 IF ( l == -1 ) THEN 1500 1504 IF ( .NOT. ALLOCATED( surf_usm_h%wghf_eb_green_av ) ) THEN … … 1511 1515 CASE ( 'usm_iwghf' ) 1512 1516 ! 1513 !-- array of heat flux from indoor ground (wall, roof, land)1517 !-- Array of heat flux from indoor ground (wall, roof, land) 1514 1518 IF ( l == -1 ) THEN 1515 1519 IF ( .NOT. ALLOCATED( surf_usm_h%iwghf_eb_av ) ) THEN … … 1526 1530 CASE ( 'usm_iwghf_window' ) 1527 1531 ! 1528 !-- array of heat flux from indoor window ground (wall, roof, land)1532 !-- Array of heat flux from indoor window ground (wall, roof, land) 1529 1533 IF ( l == -1 ) THEN 1530 1534 IF ( .NOT. ALLOCATED( surf_usm_h%iwghf_eb_window_av ) ) THEN … … 1541 1545 CASE ( 'usm_t_surf_wall' ) 1542 1546 ! 1543 !-- surface temperature for surfaces1547 !-- Surface temperature for surfaces 1544 1548 IF ( l == -1 ) THEN 1545 1549 IF ( .NOT. ALLOCATED( surf_usm_h%t_surf_wall_av ) ) THEN … … 1556 1560 CASE ( 'usm_t_surf_window' ) 1557 1561 ! 1558 !-- surface temperature for window surfaces1562 !-- Surface temperature for window surfaces 1559 1563 IF ( l == -1 ) THEN 1560 1564 IF ( .NOT. ALLOCATED( surf_usm_h%t_surf_window_av ) ) THEN … … 1571 1575 CASE ( 'usm_t_surf_green' ) 1572 1576 ! 1573 !-- surface temperature for green surfaces1577 !-- Surface temperature for green surfaces 1574 1578 IF ( l == -1 ) THEN 1575 1579 IF ( .NOT. ALLOCATED( surf_usm_h%t_surf_green_av ) ) THEN … … 1586 1590 CASE ( 'usm_theta_10cm' ) 1587 1591 ! 1588 !-- near surface (10cm) temperature for whole surfaces1592 !-- Near surface (10cm) temperature for whole surfaces 1589 1593 IF ( l == -1 ) THEN 1590 1594 IF ( .NOT. ALLOCATED( surf_usm_h%pt_10cm_av ) ) THEN … … 1601 1605 CASE ( 'usm_t_wall' ) 1602 1606 ! 1603 !-- wall temperature for iwl layer of walls and land1607 !-- Wall temperature for iwl layer of walls and land 1604 1608 IF ( l == -1 ) THEN 1605 1609 IF ( .NOT. ALLOCATED( surf_usm_h%t_wall_av ) ) THEN … … 1616 1620 CASE ( 'usm_t_window' ) 1617 1621 ! 1618 !-- window temperature for iwl layer of walls and land1622 !-- Window temperature for iwl layer of walls and land 1619 1623 IF ( l == -1 ) THEN 1620 1624 IF ( .NOT. ALLOCATED( surf_usm_h%t_window_av ) ) THEN … … 1631 1635 CASE ( 'usm_t_green' ) 1632 1636 ! 1633 !-- green temperature for iwl layer of walls and land1637 !-- Green temperature for iwl layer of walls and land 1634 1638 IF ( l == -1 ) THEN 1635 1639 IF ( .NOT. ALLOCATED( surf_usm_h%t_green_av ) ) THEN … … 1645 1649 CASE ( 'usm_swc' ) 1646 1650 ! 1647 !-- soil water content for iwl layer of walls and land1651 !-- Soil water content for iwl layer of walls and land 1648 1652 IF ( l == -1 .AND. .NOT. ALLOCATED( surf_usm_h%swc_av ) ) THEN 1649 1653 ALLOCATE ( surf_usm_h%swc_av(nzb_wall:nzt_wall,1:surf_usm_h%ns) ) … … 1667 1671 CASE ( 'usm_wshf' ) 1668 1672 ! 1669 !-- array of sensible heat flux from surfaces (land, roof, wall)1673 !-- Array of sensible heat flux from surfaces (land, roof, wall) 1670 1674 IF ( l == -1 ) THEN 1671 1675 DO m = 1, surf_usm_h%ns … … 1681 1685 CASE ( 'usm_qsws' ) 1682 1686 ! 1683 !-- array of latent heat flux from surfaces (land, roof, wall)1687 !-- Array of latent heat flux from surfaces (land, roof, wall) 1684 1688 IF ( l == -1 ) THEN 1685 1689 DO m = 1, surf_usm_h%ns … … 1695 1699 CASE ( 'usm_qsws_veg' ) 1696 1700 ! 1697 !-- array of latent heat flux from vegetation surfaces (land, roof, wall)1701 !-- Array of latent heat flux from vegetation surfaces (land, roof, wall) 1698 1702 IF ( l == -1 ) THEN 1699 1703 DO m = 1, surf_usm_h%ns … … 1709 1713 CASE ( 'usm_qsws_liq' ) 1710 1714 ! 1711 !-- array of latent heat flux from surfaces with liquid (land, roof, wall)1715 !-- Array of latent heat flux from surfaces with liquid (land, roof, wall) 1712 1716 IF ( l == -1 ) THEN 1713 1717 DO m = 1, surf_usm_h%ns … … 1724 1728 CASE ( 'usm_wghf' ) 1725 1729 ! 1726 !-- array of heat flux from ground (wall, roof, land)1730 !-- Array of heat flux from ground (wall, roof, land) 1727 1731 IF ( l == -1 ) THEN 1728 1732 DO m = 1, surf_usm_h%ns … … 1739 1743 CASE ( 'usm_wghf_window' ) 1740 1744 ! 1741 !-- array of heat flux from window ground (wall, roof, land)1745 !-- Array of heat flux from window ground (wall, roof, land) 1742 1746 IF ( l == -1 ) THEN 1743 1747 DO m = 1, surf_usm_h%ns … … 1754 1758 CASE ( 'usm_wghf_green' ) 1755 1759 ! 1756 !-- array of heat flux from green ground (wall, roof, land)1760 !-- Array of heat flux from green ground (wall, roof, land) 1757 1761 IF ( l == -1 ) THEN 1758 1762 DO m = 1, surf_usm_h%ns … … 1769 1773 CASE ( 'usm_iwghf' ) 1770 1774 ! 1771 !-- array of heat flux from indoor ground (wall, roof, land)1775 !-- Array of heat flux from indoor ground (wall, roof, land) 1772 1776 IF ( l == -1 ) THEN 1773 1777 DO m = 1, surf_usm_h%ns … … 1783 1787 CASE ( 'usm_iwghf_window' ) 1784 1788 ! 1785 !-- array of heat flux from indoor window ground (wall, roof, land)1789 !-- Array of heat flux from indoor window ground (wall, roof, land) 1786 1790 IF ( l == -1 ) THEN 1787 1791 DO m = 1, surf_usm_h%ns … … 1798 1802 CASE ( 'usm_t_surf_wall' ) 1799 1803 ! 1800 !-- surface temperature for surfaces1804 !-- Surface temperature for surfaces 1801 1805 IF ( l == -1 ) THEN 1802 1806 DO m = 1, surf_usm_h%ns … … 1812 1816 CASE ( 'usm_t_surf_window' ) 1813 1817 ! 1814 !-- surface temperature for window surfaces1818 !-- Surface temperature for window surfaces 1815 1819 IF ( l == -1 ) THEN 1816 1820 DO m = 1, surf_usm_h%ns … … 1827 1831 CASE ( 'usm_t_surf_green' ) 1828 1832 ! 1829 !-- surface temperature for green surfaces1833 !-- Surface temperature for green surfaces 1830 1834 IF ( l == -1 ) THEN 1831 1835 DO m = 1, surf_usm_h%ns … … 1842 1846 CASE ( 'usm_theta_10cm' ) 1843 1847 ! 1844 !-- near surface temperature for whole surfaces1848 !-- Near surface temperature for whole surfaces 1845 1849 IF ( l == -1 ) THEN 1846 1850 DO m = 1, surf_usm_h%ns … … 1857 1861 CASE ( 'usm_t_wall' ) 1858 1862 ! 1859 !-- wall temperature for iwl layer of walls and land1863 !-- Wall temperature for iwl layer of walls and land 1860 1864 IF ( l == -1 ) THEN 1861 1865 DO m = 1, surf_usm_h%ns … … 1872 1876 CASE ( 'usm_t_window' ) 1873 1877 ! 1874 !-- window temperature for iwl layer of walls and land1878 !-- Window temperature for iwl layer of walls and land 1875 1879 IF ( l == -1 ) THEN 1876 1880 DO m = 1, surf_usm_h%ns … … 1887 1891 CASE ( 'usm_t_green' ) 1888 1892 ! 1889 !-- green temperature for iwl layer of walls and land1893 !-- Green temperature for iwl layer of walls and land 1890 1894 IF ( l == -1 ) THEN 1891 1895 DO m = 1, surf_usm_h%ns … … 1901 1905 CASE ( 'usm_swc' ) 1902 1906 ! 1903 !-- soil water content for iwl layer of walls and land1907 !-- Soil water content for iwl layer of walls and land 1904 1908 IF ( l == -1 ) THEN 1905 1909 DO m = 1, surf_usm_h%ns … … 1920 1924 CASE ( 'usm_wshf' ) 1921 1925 ! 1922 !-- array of sensible heat flux from surfaces (land, roof, wall)1926 !-- Array of sensible heat flux from surfaces (land, roof, wall) 1923 1927 IF ( l == -1 ) THEN 1924 1928 DO m = 1, surf_usm_h%ns … … 1935 1939 CASE ( 'usm_qsws' ) 1936 1940 ! 1937 !-- array of latent heat flux from surfaces (land, roof, wall)1941 !-- Array of latent heat flux from surfaces (land, roof, wall) 1938 1942 IF ( l == -1 ) THEN 1939 1943 DO m = 1, surf_usm_h%ns … … 1950 1954 CASE ( 'usm_qsws_veg' ) 1951 1955 ! 1952 !-- array of latent heat flux from vegetation surfaces (land, roof, wall)1956 !-- Array of latent heat flux from vegetation surfaces (land, roof, wall) 1953 1957 IF ( l == -1 ) THEN 1954 1958 DO m = 1, surf_usm_h%ns … … 1965 1969 CASE ( 'usm_qsws_liq' ) 1966 1970 ! 1967 !-- array of latent heat flux from surfaces with liquid (land, roof, wall)1971 !-- Array of latent heat flux from surfaces with liquid (land, roof, wall) 1968 1972 IF ( l == -1 ) THEN 1969 1973 DO m = 1, surf_usm_h%ns … … 1980 1984 CASE ( 'usm_wghf' ) 1981 1985 ! 1982 !-- array of heat flux from ground (wall, roof, land)1986 !-- Array of heat flux from ground (wall, roof, land) 1983 1987 IF ( l == -1 ) THEN 1984 1988 DO m = 1, surf_usm_h%ns … … 1995 1999 CASE ( 'usm_wghf_window' ) 1996 2000 ! 1997 !-- array of heat flux from window ground (wall, roof, land)2001 !-- Array of heat flux from window ground (wall, roof, land) 1998 2002 IF ( l == -1 ) THEN 1999 2003 DO m = 1, surf_usm_h%ns … … 2010 2014 CASE ( 'usm_wghf_green' ) 2011 2015 ! 2012 !-- array of heat flux from green ground (wall, roof, land)2016 !-- Array of heat flux from green ground (wall, roof, land) 2013 2017 IF ( l == -1 ) THEN 2014 2018 DO m = 1, surf_usm_h%ns … … 2025 2029 CASE ( 'usm_iwghf' ) 2026 2030 ! 2027 !-- array of heat flux from indoor ground (wall, roof, land)2031 !-- Array of heat flux from indoor ground (wall, roof, land) 2028 2032 IF ( l == -1 ) THEN 2029 2033 DO m = 1, surf_usm_h%ns … … 2040 2044 CASE ( 'usm_iwghf_window' ) 2041 2045 ! 2042 !-- array of heat flux from indoor window ground (wall, roof, land)2046 !-- Array of heat flux from indoor window ground (wall, roof, land) 2043 2047 IF ( l == -1 ) THEN 2044 2048 DO m = 1, surf_usm_h%ns … … 2055 2059 CASE ( 'usm_t_surf_wall' ) 2056 2060 ! 2057 !-- surface temperature for surfaces2061 !-- Surface temperature for surfaces 2058 2062 IF ( l == -1 ) THEN 2059 2063 DO m = 1, surf_usm_h%ns … … 2070 2074 CASE ( 'usm_t_surf_window' ) 2071 2075 ! 2072 !-- surface temperature for window surfaces2076 !-- Surface temperature for window surfaces 2073 2077 IF ( l == -1 ) THEN 2074 2078 DO m = 1, surf_usm_h%ns … … 2085 2089 CASE ( 'usm_t_surf_green' ) 2086 2090 ! 2087 !-- surface temperature for green surfaces2091 !-- Surface temperature for green surfaces 2088 2092 IF ( l == -1 ) THEN 2089 2093 DO m = 1, surf_usm_h%ns … … 2100 2104 CASE ( 'usm_theta_10cm' ) 2101 2105 ! 2102 !-- near surface temperature for whole surfaces2106 !-- Near surface temperature for whole surfaces 2103 2107 IF ( l == -1 ) THEN 2104 2108 DO m = 1, surf_usm_h%ns … … 2116 2120 CASE ( 'usm_t_wall' ) 2117 2121 ! 2118 !-- wall temperature for iwl layer of walls and land2122 !-- Wall temperature for iwl layer of walls and land 2119 2123 IF ( l == -1 ) THEN 2120 2124 DO m = 1, surf_usm_h%ns … … 2131 2135 CASE ( 'usm_t_window' ) 2132 2136 ! 2133 !-- window temperature for iwl layer of walls and land2137 !-- Window temperature for iwl layer of walls and land 2134 2138 IF ( l == -1 ) THEN 2135 2139 DO m = 1, surf_usm_h%ns … … 2146 2150 CASE ( 'usm_t_green' ) 2147 2151 ! 2148 !-- green temperature for iwl layer of walls and land2152 !-- Green temperature for iwl layer of walls and land 2149 2153 IF ( l == -1 ) THEN 2150 2154 DO m = 1, surf_usm_h%ns … … 2161 2165 CASE ( 'usm_swc' ) 2162 2166 ! 2163 !-- soil water content for iwl layer of walls and land2167 !-- Soil water content for iwl layer of walls and land 2164 2168 IF ( l == -1 ) THEN 2165 2169 DO m = 1, surf_usm_h%ns … … 2281 2285 2282 2286 ! 2283 !-- check if variable exists2284 !-- directional variables2287 !-- Check if variable exists 2288 !-- Directional variables 2285 2289 DO i = 1, nl1 2286 2290 DO j = 1, nd … … 2294 2298 IF ( lfound ) GOTO 10 2295 2299 ! 2296 !-- directional layer variables2300 !-- Directional layer variables 2297 2301 DO i = 1, nl2 2298 2302 DO j = 1, nd … … 2388 2392 ENDIF 2389 2393 ! 2390 !-- naheatlayers2394 !-- Naheatlayers 2391 2395 IF ( naheatlayers > nzt ) THEN 2392 2396 message_string = 'number of anthropogenic heat layers "naheatlayers" can not be larger ' // & … … 2473 2477 IF ( var(1:11) == 'usm_t_wall_' .AND. len( TRIM( var ) ) >= 12 ) THEN 2474 2478 ! 2475 !-- wall layers2479 !-- Wall layers 2476 2480 READ( var(12:12), '(I1)', iostat = istat ) iwl 2477 2481 IF ( istat == 0 .AND. iwl >= nzb_wall .AND. iwl <= nzt_wall ) THEN … … 2481 2485 IF ( var(1:13) == 'usm_t_window_' .AND. len( TRIM( var ) ) >= 14 ) THEN 2482 2486 ! 2483 !-- window layers2487 !-- Window layers 2484 2488 READ( var(14:14), '(I1)', iostat = istat ) iwl 2485 2489 IF ( istat == 0 .AND. iwl >= nzb_wall .AND. iwl <= nzt_wall ) THEN … … 2489 2493 IF ( var(1:12) == 'usm_t_green_' .AND. len( TRIM( var ) ) >= 13 ) THEN 2490 2494 ! 2491 !-- green layers2495 !-- Green layers 2492 2496 READ( var(13:13), '(I1)', iostat = istat ) iwl 2493 2497 IF ( istat == 0 .AND. iwl >= nzb_wall .AND. iwl <= nzt_wall ) THEN … … 2497 2501 IF ( var(1:8) == 'usm_swc_' .AND. len( TRIM( var ) ) >= 9 ) THEN 2498 2502 ! 2499 !-- green layers soil water content2503 !-- Green layers soil water content 2500 2504 READ( var(9:9), '(I1)', iostat = istat ) iwl 2501 2505 IF ( istat == 0 .AND. iwl >= nzb_wall .AND. iwl <= nzt_wall ) THEN … … 2508 2512 CASE ( 'usm_surfz' ) 2509 2513 ! 2510 !-- array of surface height (z)2514 !-- Array of surface height (z) 2511 2515 IF ( idsint == iup_u ) THEN 2512 2516 DO m = 1, surf_usm_h%ns … … 2528 2532 CASE ( 'usm_surfcat' ) 2529 2533 ! 2530 !-- surface category2534 !-- Surface category 2531 2535 IF ( idsint == iup_u ) THEN 2532 2536 DO m = 1, surf_usm_h%ns … … 2548 2552 CASE ( 'usm_surfwintrans' ) 2549 2553 ! 2550 !-- transmissivity window tiles2554 !-- Transmissivity window tiles 2551 2555 IF ( idsint == iup_u ) THEN 2552 2556 DO m = 1, surf_usm_h%ns … … 2568 2572 CASE ( 'usm_wshf' ) 2569 2573 ! 2570 !-- array of sensible heat flux from surfaces2574 !-- Array of sensible heat flux from surfaces 2571 2575 IF ( av == 0 ) THEN 2572 2576 IF ( idsint == iup_u ) THEN … … 2608 2612 CASE ( 'usm_qsws' ) 2609 2613 ! 2610 !-- array of latent heat flux from surfaces2614 !-- Array of latent heat flux from surfaces 2611 2615 IF ( av == 0 ) THEN 2612 2616 IF ( idsint == iup_u ) THEN … … 2647 2651 CASE ( 'usm_qsws_veg' ) 2648 2652 ! 2649 !-- array of latent heat flux from vegetation surfaces2653 !-- Array of latent heat flux from vegetation surfaces 2650 2654 IF ( av == 0 ) THEN 2651 2655 IF ( idsint == iup_u ) THEN … … 2686 2690 CASE ( 'usm_qsws_liq' ) 2687 2691 ! 2688 !-- array of latent heat flux from surfaces with liquid2692 !-- Array of latent heat flux from surfaces with liquid 2689 2693 IF ( av == 0 ) THEN 2690 2694 IF ( idsint == iup_u ) THEN … … 2725 2729 CASE ( 'usm_wghf' ) 2726 2730 ! 2727 !-- array of heat flux from ground (land, wall, roof)2731 !-- Array of heat flux from ground (land, wall, roof) 2728 2732 IF ( av == 0 ) THEN 2729 2733 IF ( idsint == iup_u ) THEN … … 2764 2768 CASE ( 'usm_wghf_window' ) 2765 2769 ! 2766 !-- array of heat flux from window ground (land, wall, roof)2770 !-- Array of heat flux from window ground (land, wall, roof) 2767 2771 IF ( av == 0 ) THEN 2768 2772 IF ( idsint == iup_u ) THEN … … 2803 2807 CASE ( 'usm_wghf_green' ) 2804 2808 ! 2805 !-- array of heat flux from green ground (land, wall, roof)2809 !-- Array of heat flux from green ground (land, wall, roof) 2806 2810 IF ( av == 0 ) THEN 2807 2811 IF ( idsint == iup_u ) THEN … … 2842 2846 CASE ( 'usm_iwghf' ) 2843 2847 ! 2844 !-- array of heat flux from indoor ground (land, wall, roof)2848 !-- Array of heat flux from indoor ground (land, wall, roof) 2845 2849 IF ( av == 0 ) THEN 2846 2850 IF ( idsint == iup_u ) THEN … … 2881 2885 CASE ( 'usm_iwghf_window' ) 2882 2886 ! 2883 !-- array of heat flux from indoor window ground (land, wall, roof)2887 !-- Array of heat flux from indoor window ground (land, wall, roof) 2884 2888 IF ( av == 0 ) THEN 2885 2889 IF ( idsint == iup_u ) THEN … … 2920 2924 CASE ( 'usm_t_surf_wall' ) 2921 2925 ! 2922 !-- surface temperature for surfaces2926 !-- Surface temperature for surfaces 2923 2927 IF ( av == 0 ) THEN 2924 2928 IF ( idsint == iup_u ) THEN … … 2959 2963 CASE ( 'usm_t_surf_window' ) 2960 2964 ! 2961 !-- surface temperature for window surfaces2965 !-- Surface temperature for window surfaces 2962 2966 IF ( av == 0 ) THEN 2963 2967 IF ( idsint == iup_u ) THEN … … 3001 3005 CASE ( 'usm_t_surf_green' ) 3002 3006 ! 3003 !-- surface temperature for green surfaces3007 !-- Surface temperature for green surfaces 3004 3008 IF ( av == 0 ) THEN 3005 3009 IF ( idsint == iup_u ) THEN … … 3043 3047 CASE ( 'usm_theta_10cm' ) 3044 3048 ! 3045 !-- near surface temperature for whole surfaces3049 !-- Near surface temperature for whole surfaces 3046 3050 IF ( av == 0 ) THEN 3047 3051 IF ( idsint == iup_u ) THEN … … 3085 3089 CASE ( 'usm_t_wall' ) 3086 3090 ! 3087 !-- wall temperature for iwl layer of walls and land3091 !-- Wall temperature for iwl layer of walls and land 3088 3092 IF ( av == 0 ) THEN 3089 3093 IF ( idsint == iup_u ) THEN … … 3124 3128 CASE ( 'usm_t_window' ) 3125 3129 ! 3126 !-- window temperature for iwl layer of walls and land3130 !-- Window temperature for iwl layer of walls and land 3127 3131 IF ( av == 0 ) THEN 3128 3132 IF ( idsint == iup_u ) THEN … … 3163 3167 CASE ( 'usm_t_green' ) 3164 3168 ! 3165 !-- green temperature for iwl layer of walls and land3169 !-- Green temperature for iwl layer of walls and land 3166 3170 IF ( av == 0 ) THEN 3167 3171 IF ( idsint == iup_u ) THEN … … 3202 3206 CASE ( 'usm_swc' ) 3203 3207 ! 3204 !-- soil water content for iwl layer of walls and land3208 !-- Soil water content for iwl layer of walls and land 3205 3209 IF ( av == 0 ) THEN 3206 3210 IF ( idsint == iup_u ) THEN … … 3346 3350 IF (surf_usm_h%green_type_roof(m) == 2.0_wp ) THEN 3347 3351 ! 3348 !-- extensive green roof3349 !-- set ratio of substrate layer thickness, soil-type and LAI3352 !-- Extensive green roof 3353 !-- Set ratio of substrate layer thickness, soil-type and LAI 3350 3354 soil_type = 3 3351 3355 surf_usm_h%lai(m) = 2.0_wp … … 3357 3361 ELSE 3358 3362 ! 3359 !-- intensiv green roof3360 !-- set ratio of substrate layer thickness, soil-type and LAI3363 !-- Intensiv green roof 3364 !-- Set ratio of substrate layer thickness, soil-type and LAI 3361 3365 soil_type = 6 3362 3366 surf_usm_h%lai(m) = 4.0_wp … … 3560 3564 CALL cpu_log( log_point_s(78), 'usm_init', 'start' ) 3561 3565 ! 3562 !-- surface forcing has to be disabled for LSF in case of enabled urban surface module3566 !-- Surface forcing has to be disabled for LSF in case of enabled urban surface module 3563 3567 IF ( large_scale_forcing ) THEN 3564 3568 lsf_surf = .FALSE. … … 3694 3698 surf_usm_h%target_temp_winter(m) = building_pars(ind_indoor_target_temp_winter,building_type) 3695 3699 ! 3696 !-- emissivity of wall-, green- and window fraction3700 !-- Emissivity of wall-, green- and window fraction 3697 3701 surf_usm_h%emissivity(m,ind_veg_wall) = building_pars(ind_emis_wall_r,building_type) 3698 3702 surf_usm_h%emissivity(m,ind_pav_green) = building_pars(ind_emis_green_r,building_type) … … 3705 3709 surf_usm_h%z0q(m) = building_pars(ind_z0qh,building_type) 3706 3710 ! 3707 !-- albedo type for wall fraction, green fraction, window fraction3711 !-- Albedo type for wall fraction, green fraction, window fraction 3708 3712 surf_usm_h%albedo_type(m,ind_veg_wall) = INT( building_pars(ind_alb_wall_r,building_type) ) 3709 3713 surf_usm_h%albedo_type(m,ind_pav_green) = INT( building_pars(ind_alb_green_r,building_type) ) … … 3732 3736 DO m = 1, surf_usm_v(l)%ns 3733 3737 3734 surf_usm_v(l)%surface_types(m) = wall_category !< default category for root surface3738 surf_usm_v(l)%surface_types(m) = wall_category !< Default category for root surface 3735 3739 ! 3736 3740 !-- In order to distinguish between ground floor level and above-ground-floor level surfaces, … … 3850 3854 surf_usm_v(l)%target_temp_winter(m) = building_pars(ind_indoor_target_temp_winter,building_type) 3851 3855 ! 3852 !-- emissivity of wall-, green- and window fraction3856 !-- Emissivity of wall-, green- and window fraction 3853 3857 surf_usm_v(l)%emissivity(m,ind_veg_wall) = building_pars(ind_emis_wall,building_type) 3854 3858 surf_usm_v(l)%emissivity(m,ind_pav_green) = building_pars(ind_emis_green,building_type) … … 3943 3947 surf_usm_h%target_temp_winter(m) = building_pars(ind_indoor_target_temp_winter,st) 3944 3948 ! 3945 !-- emissivity of wall-, green- and window fraction3949 !-- Emissivity of wall-, green- and window fraction 3946 3950 surf_usm_h%emissivity(m,ind_veg_wall) = building_pars(ind_emis_wall_r,st) 3947 3951 surf_usm_h%emissivity(m,ind_pav_green) = building_pars(ind_emis_green_r,st) … … 3954 3958 surf_usm_h%z0q(m) = building_pars(ind_z0qh,st) 3955 3959 ! 3956 !-- albedo type for wall fraction, green fraction, window fraction3960 !-- Albedo type for wall fraction, green fraction, window fraction 3957 3961 surf_usm_h%albedo_type(m,ind_veg_wall) = INT( building_pars(ind_alb_wall_r,st) ) 3958 3962 surf_usm_h%albedo_type(m,ind_pav_green) = INT( building_pars(ind_alb_green_r,st) ) … … 4106 4110 surf_usm_v(l)%target_temp_winter(m) = building_pars(ind_indoor_target_temp_winter,st) 4107 4111 ! 4108 !-- emissivity of wall-, green- and window fraction4112 !-- Emissivity of wall-, green- and window fraction 4109 4113 surf_usm_v(l)%emissivity(m,ind_veg_wall) = building_pars(ind_emis_wall,st) 4110 4114 surf_usm_v(l)%emissivity(m,ind_pav_green) = building_pars(ind_emis_green,st) … … 4659 4663 ENDIF 4660 4664 4661 EXIT ! surface was found and processed4665 EXIT ! Surface was found and processed 4662 4666 ENDIF 4663 4667 ENDDO … … 4691 4695 surf_usm_v(l)%frac(m,ind_pav_green) = & 4692 4696 building_surface_pars_f%pars(ind_s_green_frac_r,is) 4693 !TODO clarify: why should _w and _r be on the same surface?4697 !TODO Clarify: why should _w and _r be on the same surface? 4694 4698 4695 4699 IF ( building_surface_pars_f%pars(ind_s_win_frac,is) /= & … … 4802 4806 ENDIF 4803 4807 4804 EXIT ! surface was found and processed4808 EXIT ! Surface was found and processed 4805 4809 ENDIF 4806 4810 ENDDO … … 4883 4887 CALL usm_init_material_model() 4884 4888 4885 !-- init skin layer properties (can be done after initialization of wall layers)4889 !-- Init skin layer properties (can be done after initialization of wall layers) 4886 4890 4887 4891 DO m = 1, surf_usm_h%ns … … 4924 4928 4925 4929 ! 4926 !-- init anthropogenic sources of heat4930 !-- Init anthropogenic sources of heat 4927 4931 IF ( usm_anthropogenic_heat ) THEN 4928 4932 ! 4929 !-- init anthropogenic sources of heat (from transportation for now)4933 !-- Init anthropogenic sources of heat (from transportation for now) 4930 4934 CALL usm_read_anthropogenic_heat() 4931 4935 ENDIF … … 5031 5035 ENDIF 5032 5036 ! 5033 !-- initial values for t_wall5034 !-- outer value is set to surface temperature, inner value is set to wall_inner_temperature5037 !-- Initial values for t_wall 5038 !-- Outer value is set to surface temperature, inner value is set to wall_inner_temperature 5035 5039 !-- and profile is logaritmic (linear in nz). 5036 5040 !-- Horizontal surfaces … … 5090 5094 5091 5095 ! 5092 !-- initialize prognostic values for the first timestep5096 !-- Initialize prognostic values for the first timestep 5093 5097 t_surf_wall_h_p = t_surf_wall_h 5094 5098 t_surf_wall_v_p = t_surf_wall_v … … 5182 5186 k = surf_usm_h%k(m) 5183 5187 ! 5184 !-- prognostic equation for ground/roof temperature t_wall_h5188 !-- Prognostic equation for ground/roof temperature t_wall_h 5185 5189 wtend(:) = 0.0_wp 5186 5190 wtend(nzb_wall) = ( 1.0_wp / surf_usm_h%rho_c_wall(nzb_wall,m) ) & … … 5206 5210 ) * surf_usm_h%ddz_wall_stag(nzb_wall,m) 5207 5211 ! 5208 !-- if indoor model is used inner wall layer is calculated by using iwghf (indoor wall ground heat flux)5212 !-- If indoor model is used inner wall layer is calculated by using iwghf (indoor wall ground heat flux) 5209 5213 IF ( indoor_model ) THEN 5210 5214 DO kw = nzb_wall+1, nzt_wall-1 … … 5245 5249 5246 5250 ! 5247 !-- during spinup the tempeature inside window layers is not calculated to make larger timesteps possible5251 !-- During spinup the tempeature inside window layers is not calculated to make larger timesteps possible 5248 5252 IF ( .NOT. during_spinup ) THEN 5249 5253 win_absorp = -log( surf_usm_h%transmissivity(m) ) / surf_usm_h%zw_window(nzt_wall,m) 5250 5254 ! 5251 !-- prognostic equation for ground/roof window temperature t_window_h takes absorption of5255 !-- Prognostic equation for ground/roof window temperature t_window_h takes absorption of 5252 5256 !-- shortwave radiation into account 5253 5257 wintend(:) = 0.0_wp … … 5312 5316 5313 5317 ! 5314 !-- calculate t_wall tendencies for the next Runge-Kutta step5318 !-- Calculate t_wall tendencies for the next Runge-Kutta step 5315 5319 IF ( timestep_scheme(1:5) == 'runge' ) THEN 5316 5320 IF ( intermediate_timestep_count == 1 ) THEN … … 5328 5332 IF ( .NOT. during_spinup ) THEN 5329 5333 ! 5330 !-- calculate t_window tendencies for the next Runge-Kutta step5334 !-- Calculate t_window tendencies for the next Runge-Kutta step 5331 5335 IF ( timestep_scheme(1:5) == 'runge' ) THEN 5332 5336 IF ( intermediate_timestep_count == 1 ) THEN … … 5356 5360 k = surf_usm_v(l)%k(m) 5357 5361 ! 5358 !-- prognostic equation for wall temperature t_wall_v5362 !-- Prognostic equation for wall temperature t_wall_v 5359 5363 wtend(:) = 0.0_wp 5360 5364 … … 5425 5429 surf_usm_v(l)%zw_window(nzt_wall,m) 5426 5430 ! 5427 !-- prognostic equation for window temperature t_window_v5431 !-- Prognostic equation for window temperature t_window_v 5428 5432 wintend(:) = 0.0_wp 5429 5433 wintend(nzb_wall) = ( 1.0_wp / surf_usm_v(l)%rho_c_window(nzb_wall,m) ) & … … 5492 5496 5493 5497 ! 5494 !-- calculate t_wall tendencies for the next Runge-Kutta step5498 !-- Calculate t_wall tendencies for the next Runge-Kutta step 5495 5499 IF ( timestep_scheme(1:5) == 'runge' ) THEN 5496 5500 IF ( intermediate_timestep_count == 1 ) THEN … … 5509 5513 IF ( .NOT. during_spinup ) THEN 5510 5514 ! 5511 !-- calculate t_window tendencies for the next Runge-Kutta step5515 !-- Calculate t_window tendencies for the next Runge-Kutta step 5512 5516 IF ( timestep_scheme(1:5) == 'runge' ) THEN 5513 5517 IF ( intermediate_timestep_count == 1 ) THEN … … 5837 5841 IF (surf_usm_v(l)%frac(m,ind_pav_green) > 0.0_wp) THEN 5838 5842 ! 5839 !-- no substrate layer for green walls / only groundbase green walls (ivy i.e.) -> green layers get5843 !-- No substrate layer for green walls / only groundbase green walls (ivy i.e.) -> Green layers get 5840 5844 !-- same temperature as first wall layer, therefore no temperature calculations for vertical green 5841 5845 !-- substrate layers now … … 5850 5854 ! t_green_v(l)%t(nzt_wall+1,m) = t_wall_v(l)%t(nzb_wall,m) 5851 5855 ! ! 5852 ! !-- prognostic equation for green temperature t_green_v5856 ! !-- Prognostic equation for green temperature t_green_v 5853 5857 ! gtend(:) = 0.0_wp 5854 5858 ! gtend(nzb_wall) = (1.0_wp / surf_usm_v(l)%rho_c_green(nzb_wall,m)) * & … … 5878 5882 ! 5879 5883 ! ! 5880 ! !-- calculate t_green tendencies for the next Runge-Kutta step5884 ! !-- Calculate t_green tendencies for the next Runge-Kutta step 5881 5885 ! IF ( timestep_scheme(1:5) == 'runge' ) THEN 5882 5886 ! IF ( intermediate_timestep_count == 1 ) THEN … … 6037 6041 IF ( ii == io_group ) THEN 6038 6042 6039 !-- open anthropogenic heat file6043 !-- Open anthropogenic heat file 6040 6044 OPEN( 151, file = 'ANTHROPOGENIC_HEAT' // TRIM( coupling_char ), action = 'read', & 6041 6045 status = 'old', form = 'formatted', err = 11 ) … … 6046 6050 IF ( i >= nxl .AND. i <= nxr .AND. j >= nys .AND. j <= nyn ) THEN 6047 6051 IF ( k <= naheatlayers .AND. k > topo_top_ind(j,i,0) ) THEN 6048 !-- write heat into the array6052 !-- Write heat into the array 6049 6053 aheat(k,j,i) = heat 6050 6054 ENDIF … … 6704 6708 6705 6709 ! 6706 !-- read categories of walls and their parameters6710 !-- Read categories of walls and their parameters 6707 6711 DO ii = 0, io_blocks-1 6708 6712 IF ( ii == io_group ) THEN 6709 6713 ! 6710 !-- open urban surface file6714 !-- Open urban surface file 6711 6715 OPEN( 151, file = 'SURFACE_PARAMETERS' // coupling_char, action = 'read', & 6712 6716 status = 'old', form = 'formatted', err = 15 ) 6713 6717 ! 6714 !-- first test and get n_surface_types6718 !-- First test and get n_surface_types 6715 6719 k = 0 6716 6720 l = 0 … … 6727 6731 ALLOCATE( surface_params(n_surface_params, n_surface_types) ) 6728 6732 ! 6729 !-- real reading6733 !-- Real reading 6730 6734 rewind( 151 ) 6731 6735 k = 0 … … 6749 6753 6750 6754 ! 6751 !-- read types of surfaces6755 !-- Read types of surfaces 6752 6756 usm_par = 0 6753 6757 DO ii = 0, io_blocks-1 … … 6755 6759 6756 6760 ! 6757 !-- open csv urban surface file6761 !-- Open csv urban surface file 6758 6762 OPEN( 151, file = 'URBAN_SURFACE' // TRIM( coupling_char ), action = 'read', & 6759 6763 status = 'old', form = 'formatted', err = 23 ) … … 6778 6782 IF ( i >= nxlg .AND. i <= nxrg .AND. j >= nysg .AND. j <= nyng ) THEN 6779 6783 ! 6780 !-- write integer variables into array6784 !-- Write integer variables into array 6781 6785 usm_par(:,j,i) = (/1, nz, roof, dirwe, dirsn, category, & 6782 6786 weheight1, wecat1, weheight2, wecat2, weheight3, wecat3, & 6783 6787 snheight1, sncat1, snheight2, sncat2, snheight3, sncat3 /) 6784 6788 ! 6785 !-- write real values into array6789 !-- Write real values into array 6786 6790 usm_val(:,j,i) = (/ albedo, thick, & 6787 6791 wealbedo1, wethick1, wealbedo2, wethick2, & … … 6807 6811 6808 6812 ! 6809 !-- check completeness and formal correctness of the data6813 !-- Check completeness and formal correctness of the data 6810 6814 DO i = nxlg, nxrg 6811 6815 DO j = nysg, nyng … … 6829 6833 ) ) THEN 6830 6834 ! 6831 !-- incorrect input data6835 !-- Incorrect input data 6832 6836 WRITE( message_string, '(A,2I5)' ) & 6833 6837 'missing or incorrect data in file URBAN_SURFACE' // TRIM( coupling_char ) // & … … 6852 6856 IF ( zu(kw) >= roof_height_limit ) THEN 6853 6857 surf_usm_h%isroof_surf(m) = .TRUE. 6854 surf_usm_h%surface_types(m) = roof_category !< default category for root surface6858 surf_usm_h%surface_types(m) = roof_category !< Default category for root surface 6855 6859 ELSE 6856 6860 surf_usm_h%isroof_surf(m) = .FALSE. 6857 surf_usm_h%surface_types(m) = land_category !< default category for land surface6861 surf_usm_h%surface_types(m) = land_category !< Default category for land surface 6858 6862 ENDIF 6859 6863 … … 6890 6894 IF ( ip == -99999 ) THEN 6891 6895 ! 6892 !-- land/roof category not found6896 !-- Land/roof category not found 6893 6897 WRITE(9, '(A, I5, A, 3I5)' ) 'land/roof category ', it, ' not found for i, j, k = ', & 6894 6898 iw, jw, kw … … 6907 6911 IF ( ip == -99999 ) THEN 6908 6912 ! 6909 !-- default land/roof category not found6913 !-- Default land/roof category not found 6910 6914 WRITE( 9, '(A, I5, A, 3I5)' ) 'Default land/roof category ', category, ' not found!' 6911 6915 FLUSH( 9 ) … … 6927 6931 ENDIF 6928 6932 ! 6929 !-- emissivity of the wall6933 !-- Emissivity of the wall 6930 6934 surf_usm_h%emissivity(m,:) = surface_params(iemiss, ip) 6931 6935 ! 6932 !-- heat conductivity λS between air and wall ( W mâ2 Kâ1 )6936 !-- Heat conductivity λS between air and wall ( W mâ2 Kâ1 ) 6933 6937 surf_usm_h%lambda_surf(m) = surface_params(ilambdas,ip) 6934 6938 surf_usm_h%lambda_surf_window(m) = surface_params(ilambdas,ip) 6935 6939 surf_usm_h%lambda_surf_green(m) = surface_params(ilambdas,ip) 6936 6940 ! 6937 !-- roughness length for momentum, heat and humidity6941 !-- Roughness length for momentum, heat and humidity 6938 6942 surf_usm_h%z0(m) = surface_params(irough,ip) 6939 6943 surf_usm_h%z0h(m) = surface_params(iroughh,ip) … … 6945 6949 surf_usm_h%c_surface_green(m) = surface_params(icsurf,ip) 6946 6950 ! 6947 !-- wall material parameters:6948 !-- thickness of the wall (m) missing values are replaced by default value for category6951 !-- Wall material parameters: 6952 !-- Thickness of the wall (m) missing values are replaced by default value for category 6949 6953 IF ( surf_usm_h%thickness_wall(m) <= 0.001_wp ) THEN 6950 6954 surf_usm_h%thickness_wall(m) = surface_params(ithick,ip) … … 6957 6961 ENDIF 6958 6962 ! 6959 !-- volumetric heat capacity rho*C of the wall ( J mâ3 Kâ1 )6963 !-- Volumetric heat capacity rho*C of the wall ( J mâ3 Kâ1 ) 6960 6964 surf_usm_h%rho_c_wall(:,m) = surface_params(irhoC,ip) 6961 6965 surf_usm_h%rho_c_window(:,m) = surface_params(irhoC,ip) 6962 6966 surf_usm_h%rho_c_green(:,m) = surface_params(irhoC,ip) 6963 6967 ! 6964 !-- thermal conductivity λH of the wall (W mâ1 Kâ1 )6968 !-- Thermal conductivity λH of the wall (W mâ1 Kâ1 ) 6965 6969 surf_usm_h%lambda_h(:,m) = surface_params(ilambdah,ip) 6966 6970 surf_usm_h%lambda_h_window(:,m) = surface_params(ilambdah,ip) … … 6981 6985 kw = surf_usm_v(l)%k(m) 6982 6986 6983 IF ( l == 3 ) THEN ! westward facing6987 IF ( l == 3 ) THEN ! Westward facing 6984 6988 iw = i 6985 6989 jw = j … … 7005 7009 IF ( iw < 0 .OR. jw < 0 ) THEN 7006 7010 ! 7007 !-- wall on west or south border of the domain - assign default category7011 !-- Wall on west or south border of the domain - assign default category 7008 7012 IF ( kw <= roof_height_limit ) THEN 7009 surf_usm_v(l)%surface_types(m) = wall_category !< default category for wall surface in wall zone7013 surf_usm_v(l)%surface_types(m) = wall_category !< Default category for wall surface in wall zone 7010 7014 ELSE 7011 surf_usm_v(l)%surface_types(m) = roof_category !< default category for wall surface in roof zone7015 surf_usm_v(l)%surface_types(m) = roof_category !< Default category for wall surface in roof zone 7012 7016 ENDIF 7013 7017 surf_usm_v(l)%albedo(m,:) = -1.0_wp … … 7018 7022 ELSE IF ( kw <= usm_par(ii,jw,iw) ) THEN 7019 7023 ! 7020 !-- pedestrian zone7024 !-- Pedestrian zone 7021 7025 IF ( usm_par(ii+1,jw,iw) == 0 ) THEN 7022 surf_usm_v(l)%surface_types(m) = pedestrian_category !< default category for wall surface in7023 !< pedestrian zone7026 surf_usm_v(l)%surface_types(m) = pedestrian_category !< Default category for wall surface in 7027 !< Pedestrian zone 7024 7028 surf_usm_v(l)%albedo(m,:) = -1.0_wp 7025 7029 surf_usm_v(l)%thickness_wall(m) = -1.0_wp … … 7037 7041 ELSE IF ( kw <= usm_par(ii+2,jw,iw) ) THEN 7038 7042 ! 7039 !-- wall zone7043 !-- Wall zone 7040 7044 IF ( usm_par(ii+3,jw,iw) == 0 ) THEN 7041 7045 surf_usm_v(l)%surface_types(m) = wall_category !< default category for wall surface … … 7055 7059 ELSE IF ( kw <= usm_par(ii+4,jw,iw) ) THEN 7056 7060 ! 7057 !-- roof zone7061 !-- Roof zone 7058 7062 IF ( usm_par(ii+5,jw,iw) == 0 ) THEN 7059 surf_usm_v(l)%surface_types(m) = roof_category !< default category for roof surface7063 surf_usm_v(l)%surface_types(m) = roof_category !< Default category for roof surface 7060 7064 surf_usm_v(l)%albedo(m,:) = -1.0_wp 7061 7065 surf_usm_v(l)%thickness_wall(m) = -1.0_wp … … 7081 7085 FLUSH( 9 ) 7082 7086 ! 7083 !-- supply the default category7087 !-- Supply the default category 7084 7088 IF ( kw <= roof_height_limit ) THEN 7085 surf_usm_v(l)%surface_types(m) = wall_category !< default category for wall surface in wall zone7089 surf_usm_v(l)%surface_types(m) = wall_category !< Default category for wall surface in wall zone 7086 7090 ELSE 7087 surf_usm_v(l)%surface_types(m) = roof_category !< default category for wall surface in roof zone7091 surf_usm_v(l)%surface_types(m) = roof_category !< Default category for wall surface in roof zone 7088 7092 ENDIF 7089 7093 surf_usm_v(l)%albedo(m,:) = -1.0_wp … … 7105 7109 IF ( ip == -99999 ) THEN 7106 7110 ! 7107 !-- wall category not found7111 !-- Wall category not found 7108 7112 WRITE( 9, '(A,I7,A,3I5)' ) 'wall category ', it, ' not found for i,j,k = ', iw, jw, kw 7109 7113 FLUSH(9) … … 7117 7121 IF ( ip == -99999 ) THEN 7118 7122 ! 7119 !-- default wall category not found7123 !-- Default wall category not found 7120 7124 WRITE ( 9, '(A,I5,A,3I5)' ) 'Default wall category', category, ' not found!' 7121 7125 FLUSH( 9 ) … … 7136 7140 ENDIF 7137 7141 ! 7138 !-- emissivity of the wall7142 !-- Emissivity of the wall 7139 7143 surf_usm_v(l)%emissivity(:,m) = surface_params(iemiss,ip) 7140 7144 ! 7141 !-- heat conductivity lambda S between air and wall ( W m-2 K-1 )7145 !-- Heat conductivity lambda S between air and wall ( W m-2 K-1 ) 7142 7146 surf_usm_v(l)%lambda_surf(m) = surface_params(ilambdas,ip) 7143 7147 surf_usm_v(l)%lambda_surf_window(m) = surface_params(ilambdas,ip) 7144 7148 surf_usm_v(l)%lambda_surf_green(m) = surface_params(ilambdas,ip) 7145 7149 ! 7146 !-- roughness length7150 !-- Roughness length 7147 7151 surf_usm_v(l)%z0(m) = surface_params(irough,ip) 7148 7152 surf_usm_v(l)%z0h(m) = surface_params(iroughh,ip) … … 7154 7158 surf_usm_v(l)%c_surface_green(m) = surface_params(icsurf,ip) 7155 7159 ! 7156 !-- wall material parameters:7157 !-- thickness of the wall (m)7158 !-- missing values are replaced by default value for category7160 !-- Wall material parameters: 7161 !-- Thickness of the wall (m) 7162 !-- Missing values are replaced by default value for category 7159 7163 IF ( surf_usm_v(l)%thickness_wall(m) <= 0.001_wp ) THEN 7160 7164 surf_usm_v(l)%thickness_wall(m) = surface_params(ithick,ip) … … 7167 7171 ENDIF 7168 7172 ! 7169 !-- volumetric heat capacity rho*C of the wall ( J m-3 K-1 )7173 !-- Volumetric heat capacity rho*C of the wall ( J m-3 K-1 ) 7170 7174 surf_usm_v(l)%rho_c_wall(:,m) = surface_params(irhoC,ip) 7171 7175 surf_usm_v(l)%rho_c_window(:,m) = surface_params(irhoC,ip) 7172 7176 surf_usm_v(l)%rho_c_green(:,m) = surface_params(irhoC,ip) 7173 7177 ! 7174 !-- thermal conductivity lambda H of the wall (W m-1 K-1 )7178 !-- Thermal conductivity lambda H of the wall (W m-1 K-1 ) 7175 7179 surf_usm_v(l)%lambda_h(:,m) = surface_params(ilambdah,ip) 7176 7180 surf_usm_v(l)%lambda_h_window(:,m) = surface_params(ilambdah,ip) … … 7189 7193 ENDDO 7190 7194 ! 7191 !-- apply for all particular surface grids. First for horizontal surfaces7195 !-- Apply for all particular surface grids. First for horizontal surfaces 7192 7196 DO m = 1, surf_usm_h%ns 7193 7197 surf_usm_h%zw(:,m) = zwn(:) * surf_usm_h%thickness_wall(m) … … 7264 7268 IF ( ii == io_group ) THEN 7265 7269 ! 7266 !-- open wall temperature file7270 !-- Open wall temperature file 7267 7271 OPEN( 152, file = 'WALL_TEMPERATURE' // coupling_char, action = 'read', & 7268 7272 status = 'old', form = 'formatted', err = 15 ) … … 7271 7275 iline = 1 7272 7276 DO 7273 rtwall = -9999.0_wp !< for incomplete lines7277 rtwall = -9999.0_wp !< For incomplete lines 7274 7278 READ( 152, *, err = 13, end = 14 ) i, j, k, d, rtsurf, rtwall 7275 7279 7276 IF ( nxl <= i .AND. i <= nxr .AND. nys <= j .AND. j <= nyn) THEN !< local processor7280 IF ( nxl <= i .AND. i <= nxr .AND. nys <= j .AND. j <= nyn) THEN !< Local processor 7277 7281 !-- identify surface id 7278 7282 isurfl = find_surface( i, j, k, d ) … … 7284 7288 ENDIF 7285 7289 ! 7286 !-- assign temperatures7290 !-- Assign temperatures 7287 7291 IF ( d == 0 ) THEN 7288 7292 t_surf_wall_h(isurfl) = rtsurf … … 7446 7450 ENDIF 7447 7451 ! 7448 !-- calculate rho * c_p coefficient at surface layer7452 !-- Calculate rho * c_p coefficient at surface layer 7449 7453 rho_cp = c_p * hyp(k) / ( r_d * surf_usm_h%pt1(m) * exner(k) ) 7450 7454 … … 7486 7490 7487 7491 ! 7488 !-- factor for shf_eb7492 !-- Factor for shf_eb 7489 7493 f_shf = rho_cp / surf_usm_h%r_a(m) 7490 7494 f_shf_window = rho_cp / surf_usm_h%r_a_window(m) … … 7497 7501 !-- ECMWF documentation 7498 7502 7499 !-- f1: correction for incoming shortwave radiation (stomata close at night)7503 !-- f1: Correction for incoming shortwave radiation (stomata close at night) 7500 7504 f1 = MIN( 1.0_wp, ( 0.004_wp * surf_usm_h%rad_sw_in(m) + 0.05_wp ) / & 7501 7505 (0.81_wp * ( 0.004_wp * surf_usm_h%rad_sw_in(m) + 1.0_wp ) ) ) 7502 7506 ! 7503 !-- f2: correction for soil moisture availability to plants (the integrated soil moisture must7507 !-- f2: Correction for soil moisture availability to plants (the integrated soil moisture must 7504 7508 !-- thus be considered here) f2 = 0 for very dry soils 7505 7509 m_total = 0.0_wp … … 7521 7525 / ( t_surf_green_h(m) - 35.86_wp ) ) 7522 7526 ! 7523 !-- f3: correction for vapour pressure deficit7527 !-- f3: Correction for vapour pressure deficit 7524 7528 IF ( surf_usm_h%g_d(m) /= 0.0_wp ) THEN 7525 7529 ! … … 7574 7578 ENDIF 7575 7579 ! 7576 !-- add LW up so that it can be removed in prognostic equation7580 !-- Add LW up so that it can be removed in prognostic equation 7577 7581 surf_usm_h%rad_net_l(m) = surf_usm_h%rad_sw_in(m) - surf_usm_h%rad_sw_out(m) + & 7578 7582 surf_usm_h%rad_lw_in(m) - surf_usm_h%rad_lw_out(m) 7579 7583 ! 7580 !-- numerator of the prognostic equation7584 !-- Numerator of the prognostic equation 7581 7585 !-- Todo: Adjust to tile approach. So far, emissivity for wall (element 0) is used 7582 7586 coef_1 = surf_usm_h%rad_net_l(m) + ( 3.0_wp + 1.0_wp ) & … … 7603 7607 ENDIF 7604 7608 ! 7605 !-- denominator of the prognostic equation7609 !-- Denominator of the prognostic equation 7606 7610 coef_2 = 4.0_wp * surf_usm_h%emissivity(m,ind_veg_wall) * sigma_sb * t_surf_wall_h(m)**3 & 7607 7611 + lambda_surface + f_shf / exner(k) … … 7619 7623 ENDIF 7620 7624 ! 7621 !-- implicit solution when the surface layer has no heat capacity, otherwise use RK3 scheme.7625 !-- Implicit solution when the surface layer has no heat capacity, otherwise use RK3 scheme. 7622 7626 t_surf_wall_h_p(m) = ( coef_1 * dt_3d * tsc(2) + surf_usm_h%c_surface(m) & 7623 7627 * t_surf_wall_h(m) ) & … … 7632 7636 / ( surf_usm_h%c_surface_green(m) + coef_green_2 * dt_3d * tsc(2) ) 7633 7637 ! 7634 !-- add RK3 term7638 !-- Add RK3 term 7635 7639 t_surf_wall_h_p(m) = t_surf_wall_h_p(m) + dt_3d * tsc(3) * surf_usm_h%tt_surface_wall_m(m) 7636 7640 … … 7650 7654 IF ( humidity ) surf_usm_h%vpt_surface(m) = surf_usm_h%pt_surface(m) 7651 7655 ! 7652 !-- calculate true tendency7656 !-- Calculate true tendency 7653 7657 stend_wall = ( t_surf_wall_h_p(m) - t_surf_wall_h(m) - dt_3d * tsc(3) * & 7654 7658 surf_usm_h%tt_surface_wall_m(m) ) / ( dt_3d * tsc(2) ) … … 7658 7662 surf_usm_h%tt_surface_green_m(m) ) / ( dt_3d * tsc(2) ) 7659 7663 ! 7660 !-- calculate t_surf tendencies for the next Runge-Kutta step7664 !-- Calculate t_surf tendencies for the next Runge-Kutta step 7661 7665 IF ( timestep_scheme(1:5) == 'runge' ) THEN 7662 7666 IF ( intermediate_timestep_count == 1 ) THEN … … 7674 7678 ENDIF 7675 7679 ! 7676 !-- in case of fast changes in the skin temperature, it is required to update the radiative7680 !-- In case of fast changes in the skin temperature, it is required to update the radiative 7677 7681 !-- fluxes in order to keep the solution stable 7678 7682 IF ( ( ( ABS( t_surf_wall_h_p(m) - t_surf_wall_h(m) ) > 1.0_wp ) .OR. & … … 7683 7687 ENDIF 7684 7688 ! 7685 !-- calculate fluxes7686 !-- rad_net_l is never used!7689 !-- Calculate fluxes 7690 !-- Rad_net_l is never used! 7687 7691 surf_usm_h%rad_net_l(m) = surf_usm_h%rad_net_l(m) + surf_usm_h%frac(m,ind_veg_wall) & 7688 7692 * sigma_sb * surf_usm_h%emissivity(m,ind_veg_wall) & … … 7702 7706 7703 7707 ! 7704 !-- ground/wall/roof surface heat flux7708 !-- Ground/wall/roof surface heat flux 7705 7709 surf_usm_h%wshf_eb(m) = - f_shf * ( surf_usm_h%pt1(m) - t_surf_wall_h_p(m) / exner(k) ) & 7706 7710 * surf_usm_h%frac(m,ind_veg_wall) - f_shf_window & … … 7710 7714 * surf_usm_h%frac(m,ind_pav_green) 7711 7715 ! 7712 !-- store kinematic surface heat fluxes for utilization in other processes diffusion_s,7716 !-- Store kinematic surface heat fluxes for utilization in other processes diffusion_s, 7713 7717 !-- surface_layer_fluxes,... 7714 7718 surf_usm_h%shf(m) = surf_usm_h%wshf_eb(m) / c_p … … 7862 7866 ENDIF 7863 7867 ! 7864 !-- calculate rho * c_p coefficient at wall layer7868 !-- Calculate rho * c_p coefficient at wall layer 7865 7869 rho_cp = c_p * hyp(k) / ( r_d * surf_usm_v(l)%pt1(m) * exner(k) ) 7866 7870 … … 7874 7878 !-- Calculation of r_a for vertical surfaces 7875 7879 !-- 7876 !-- heat transfer coefficient for forced convection along vertical walls follows formulation7880 !-- Heat transfer coefficient for forced convection along vertical walls follows formulation 7877 7881 !-- in TUF3d model (Krayenhoff & Voogt, 2006) 7878 7882 !-- … … 7880 7884 !-- httc = rw * (11.8 + 4.2 * Ueff) - 4.0 7881 7885 !-- 7882 !-- rw: wall patch roughness relative to 1.0 for concrete7883 !-- Ueff: effective wind speed7886 !-- rw: Wall patch roughness relative to 1.0 for concrete 7887 !-- Ueff: Effective wind speed 7884 7888 !-- - 4.0 is a reduction of Rowley et al (1930) formulation based on 7885 7889 !-- Cole and Sturrock (1977) 7886 7890 !-- 7887 7891 !-- Ucan: Canyon wind speed 7888 !-- wstar: convective velocity7889 !-- Qs: surface heat flux7890 !-- zH: height of the convective layer7892 !-- wstar: Convective velocity 7893 !-- Qs: Surface heat flux 7894 !-- zH: Height of the convective layer 7891 7895 !-- wstar = (g/Tcan*Qs*zH)**(1./3.) 7892 7896 !-- Effective velocity components must always be defined at scalar grid point. The wall … … 7917 7921 / (0.81_wp * (0.004_wp * surf_usm_v(l)%rad_sw_in(m) + 1.0_wp) ) ) 7918 7922 ! 7919 !-- f2: correction for soil moisture availability to plants (the integrated soil moisture7923 !-- f2: Correction for soil moisture availability to plants (the integrated soil moisture 7920 7924 !-- must thus be considered here) f2 = 0 for very dry soils 7921 7925 f2=1.0_wp … … 7926 7930 / ( t_surf_green_v_p(l)%t(m) - 35.86_wp ) ) 7927 7931 ! 7928 !-- f3: correction for vapour pressure deficit7932 !-- f3: Correction for vapour pressure deficit 7929 7933 IF ( surf_usm_v(l)%g_d(m) /= 0.0_wp ) THEN 7930 7934 ! … … 7972 7976 7973 7977 ! 7974 !-- add LW up so that it can be removed in prognostic equation7978 !-- Add LW up so that it can be removed in prognostic equation 7975 7979 surf_usm_v(l)%rad_net_l(m) = surf_usm_v(l)%rad_sw_in(m) - surf_usm_v(l)%rad_sw_out(m) & 7976 7980 + surf_usm_v(l)%rad_lw_in(m) - surf_usm_v(l)%rad_lw_out(m) 7977 7981 ! 7978 !-- numerator of the prognostic equation7979 coef_1 = surf_usm_v(l)%rad_net_l(m) + & ! coef +1 corresponds to -lwout7982 !-- Numerator of the prognostic equation 7983 coef_1 = surf_usm_v(l)%rad_net_l(m) + & ! Coef +1 corresponds to -lwout 7980 7984 ! included in calculation of radnet_l 7981 7985 ( 3.0_wp + 1.0_wp ) * surf_usm_v(l)%emissivity(m,ind_veg_wall) * & … … 7984 7988 lambda_surface * t_wall_v(l)%t(nzb_wall,m) 7985 7989 IF ( ( .NOT. during_spinup ) .AND. ( surf_usm_v(l)%frac(m,ind_wat_win) > 0.0_wp ) ) THEN 7986 coef_window_1 = surf_usm_v(l)%rad_net_l(m) + & ! coef +1 corresponds to -lwout7990 coef_window_1 = surf_usm_v(l)%rad_net_l(m) + & ! Coef +1 corresponds to -lwout 7987 7991 ! included in calculation of radnet_l 7988 7992 ( 3.0_wp + 1.0_wp ) * surf_usm_v(l)%emissivity(m,ind_wat_win) * & … … 7992 7996 ENDIF 7993 7997 IF ( ( humidity ) .AND. ( surf_usm_v(l)%frac(m,ind_pav_green) > 0.0_wp ) ) THEN 7994 coef_green_1 = surf_usm_v(l)%rad_net_l(m) + & ! coef +1 corresponds to -lwout7998 coef_green_1 = surf_usm_v(l)%rad_net_l(m) + & ! Coef +1 corresponds to -lwout 7995 7999 ! included in calculation of radnet_l 7996 8000 ( 3.0_wp + 1.0_wp ) * surf_usm_v(l)%emissivity(m,ind_pav_green) * sigma_sb * & … … 8000 8004 lambda_surface_green * t_wall_v(l)%t(nzb_wall,m) 8001 8005 ELSE 8002 coef_green_1 = surf_usm_v(l)%rad_net_l(m) + & ! coef +1 corresponds to -lwout included8006 coef_green_1 = surf_usm_v(l)%rad_net_l(m) + & ! Coef +1 corresponds to -lwout included 8003 8007 ! in calculation of radnet_l 8004 8008 ( 3.0_wp + 1.0_wp ) * surf_usm_v(l)%emissivity(m,ind_pav_green) * sigma_sb * & … … 8009 8013 8010 8014 ! 8011 !-- denominator of the prognostic equation8015 !-- Denominator of the prognostic equation 8012 8016 coef_2 = 4.0_wp * surf_usm_v(l)%emissivity(m,ind_veg_wall) * sigma_sb & 8013 8017 * t_surf_wall_v(l)%t(m)**3 + lambda_surface + f_shf / exner(k) … … 8025 8029 ENDIF 8026 8030 ! 8027 !-- implicit solution when the surface layer has no heat capacity, otherwise use RK3 scheme.8031 !-- Implicit solution when the surface layer has no heat capacity, otherwise use RK3 scheme. 8028 8032 t_surf_wall_v_p(l)%t(m) = ( coef_1 * dt_3d * tsc(2) + surf_usm_v(l)%c_surface(m) & 8029 8033 * t_surf_wall_v(l)%t(m) ) / ( surf_usm_v(l)%c_surface(m) & … … 8042 8046 + coef_green_2 * dt_3d * tsc(2) ) 8043 8047 ! 8044 !-- add RK3 term8048 !-- Add RK3 term 8045 8049 t_surf_wall_v_p(l)%t(m) = t_surf_wall_v_p(l)%t(m) + dt_3d * tsc(3) * & 8046 8050 surf_usm_v(l)%tt_surface_wall_m(m) … … 8062 8066 IF ( humidity ) surf_usm_v(l)%vpt_surface(m) = surf_usm_v(l)%pt_surface(m) 8063 8067 ! 8064 !-- calculate true tendency8068 !-- Calculate true tendency 8065 8069 stend_wall = ( t_surf_wall_v_p(l)%t(m) - t_surf_wall_v(l)%t(m) - dt_3d * tsc(3) * & 8066 8070 surf_usm_v(l)%tt_surface_wall_m(m) ) / ( dt_3d * tsc(2) ) … … 8071 8075 8072 8076 ! 8073 !-- calculate t_surf_* tendencies for the next Runge-Kutta step8077 !-- Calculate t_surf_* tendencies for the next Runge-Kutta step 8074 8078 IF ( timestep_scheme(1:5) == 'runge' ) THEN 8075 8079 IF ( intermediate_timestep_count == 1 ) THEN … … 8088 8092 8089 8093 ! 8090 !-- in case of fast changes in the skin temperature, it is required to update the radiative8094 !-- In case of fast changes in the skin temperature, it is required to update the radiative 8091 8095 !-- fluxes in order to keep the solution stable 8092 8096 … … 8099 8103 8100 8104 ! 8101 !-- calculate fluxes8102 !-- prognostic rad_net_l is used just for output!8105 !-- Calculate fluxes 8106 !-- Prognostic rad_net_l is used just for output! 8103 8107 surf_usm_v(l)%rad_net_l(m) = surf_usm_v(l)%frac(m,ind_veg_wall) * & 8104 8108 ( surf_usm_v(l)%rad_net_l(m) + 3.0_wp * sigma_sb * & … … 8124 8128 8125 8129 ! 8126 !-- ground/wall/roof surface heat flux8130 !-- Ground/wall/roof surface heat flux 8127 8131 surf_usm_v(l)%wshf_eb(m) = - f_shf * ( surf_usm_v(l)%pt1(m) - t_surf_wall_v_p(l)%t(m) & 8128 8132 / exner(k) ) * surf_usm_v(l)%frac(m,ind_veg_wall) & … … 8134 8138 8135 8139 ! 8136 !-- store kinematic surface heat fluxes for utilization in other processes diffusion_s,8140 !-- Store kinematic surface heat fluxes for utilization in other processes diffusion_s, 8137 8141 !-- surface_layer_fluxes,... 8138 8142 surf_usm_v(l)%shf(m) = surf_usm_v(l)%wshf_eb(m) / c_p … … 8221 8225 !-- TO_DO: activate, if testcase is available 8222 8226 !-- !$OMP PARALLEL DO PRIVATE (i, j, k, acoef, rho_cp) 8223 !-- it may also improve performance to move topo_top_ind before the k-loop8227 !-- It may also improve performance to move topo_top_ind before the k-loop 8224 8228 DO i = nxl, nxr 8225 8229 DO j = nys, nyn … … 8227 8231 IF ( k > topo_top_ind(j,i,0) ) THEN 8228 8232 ! 8229 !-- increase of pt in box i,j,k in time dt_3d given to anthropogenic heat8233 !-- Increase of pt in box i,j,k in time dt_3d given to anthropogenic heat 8230 8234 !-- aheat*acoef (W*m-2) 8231 8235 !-- linear interpolation of coeficient … … 8236 8240 IF ( aheat(k,j,i) > 0.0_wp ) THEN 8237 8241 ! 8238 !-- calculate rho * c_p coefficient at layer k8242 !-- Calculate rho * c_p coefficient at layer k 8239 8243 rho_cp = c_p * hyp(k) / ( r_d * pt(k+1,j,i) * exner(k) ) 8240 8244 pt(k,j,i) = pt(k,j,i) + aheat(k,j,i) * acoef * dt_3d / (exner(k) * rho_cp & … … 8303 8307 !! ( surf_usm_h%ts(m) * surf_usm_h%us(m) + 1.0E-10_wp ) 8304 8308 !! 8305 !! !- make sure that the resistance does not drop to zero8309 !! !- Make sure that the resistance does not drop to zero 8306 8310 !! IF ( ABS(surf_usm_h%r_a_green(m)) < 1.0E-10_wp ) surf_usm_h%r_a_green(m) = 1.0E-10_wp 8307 8311 !
Note: See TracChangeset
for help on using the changeset viewer.
