Skip to content
Snippets Groups Projects
coupling_tebn.F90 58.1 KiB
Newer Older
  • Learn to ignore specific revisions
  • !SURFEX_LIC Copyright 1994-2014 Meteo-France 
    !SURFEX_LIC This is part of the SURFEX software governed by the CeCILL-C  licence
    !SURFEX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt
    !SURFEX_LIC for details. version 1.
    
    5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292 293 294 295 296 297 298 299 300 301 302 303 304 305 306 307 308 309 310 311 312 313 314 315 316 317 318 319 320 321 322 323 324 325 326 327 328 329 330 331 332 333 334 335 336 337 338 339 340 341 342 343 344 345 346 347 348 349 350 351 352 353 354 355 356 357 358 359 360 361 362 363 364 365 366 367 368 369 370 371 372 373 374 375 376 377 378 379 380 381 382 383 384 385 386 387 388 389 390 391 392 393 394 395 396 397 398 399 400 401 402 403 404 405 406 407 408 409 410 411 412 413 414 415 416 417 418 419 420 421 422 423 424 425 426 427 428 429 430 431 432 433 434 435 436 437 438 439 440 441 442 443 444 445 446 447 448 449 450 451 452 453 454 455 456 457 458 459 460 461 462 463 464 465 466 467 468 469 470 471 472 473 474 475 476 477 478 479 480 481 482 483 484 485 486 487 488 489 490 491 492 493 494 495 496 497 498 499 500 501 502 503 504 505 506 507 508 509 510 511 512 513 514 515 516 517 518 519 520 521 522 523 524 525 526 527 528 529 530 531 532 533 534 535 536 537 538 539 540 541 542 543 544 545 546 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 563 564 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594 595 596 597 598 599 600 601 602 603 604 605 606 607 608 609 610 611 612 613 614 615 616 617 618 619 620 621 622 623 624 625 626 627 628 629 630 631 632 633 634 635 636 637 638 639 640 641 642 643 644 645 646 647 648 649 650 651 652 653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705 706 707 708 709 710 711 712 713 714 715 716 717 718 719 720 721 722 723 724 725 726 727 728 729 730 731 732 733 734 735 736 737 738 739 740 741 742 743 744 745 746 747 748 749 750 751 752 753 754 755 756 757 758 759 760 761 762 763 764 765 766 767 768 769 770 771 772 773 774 775 776 777 778 779 780 781 782 783 784 785 786 787 788 789 790 791 792 793 794 795 796 797 798 799 800 801 802 803 804 805 806 807 808 809 810 811 812 813 814 815 816 817 818 819 820 821 822 823 824 825 826 827 828 829 830 831 832 833 834 835 836 837 838 839 840 841 842 843 844 845 846 847 848 849 850 851 852 853 854 855 856 857 858 859 860 861 862 863 864 865 866 867 868 869 870 871 872 873 874 875 876 877 878 879 880 881 882 883 884 885 886 887 888 889 890 891 892 893 894 895 896 897 898 899 900 901 902 903 904 905 906 907 908 909 910 911 912 913 914 915 916 917 918 919 920 921 922 923 924 925 926 927 928 929 930 931 932 933 934 935 936 937 938 939 940 941 942 943 944 945 946 947 948 949 950 951 952 953 954 955 956 957 958 959 960 961 962 963 964 965 966 967 968 969 970 971 972 973 974 975 976 977 978 979 980 981 982 983 984 985 986 987 988 989 990 991 992 993 994 995 996 997 998 999 1000
    !     ###############################################################################
    SUBROUTINE COUPLING_TEB_n(HPROGRAM, HCOUPLING,                                             &
                   PTSTEP, KYEAR, KMONTH, KDAY, PTIME, KI, KSV, KSW, PTSUN, PZENITH, PAZIM,    &
                   PZREF, PUREF, PZS, PU, PV, PQA, PTA, PRHOA, PSV, PCO2, HSV,                 &
                   PRAIN, PSNOW, PLW, PDIR_SW, PSCA_SW, PSW_BANDS, PPS, PPA,                   &
                   PSFTQ, PSFTH, PSFTS, PSFCO2, PSFU, PSFV,                                    &
                   PTRAD, PDIR_ALB, PSCA_ALB, PEMIS,                                           &
                   PPEW_A_COEF, PPEW_B_COEF,                                                   &
                   PPET_A_COEF, PPEQ_A_COEF, PPET_B_COEF, PPEQ_B_COEF,                         &
                   HTEST                                                                       )
    !     ###############################################################################
    !
    !!****  *COUPLING_TEB_n * - Driver for TEB 
    !!
    !!    PURPOSE
    !!    -------
    !
    !!**  METHOD
    !!    ------
    !!
    !!    REFERENCE
    !!    ---------
    !!      
    !!
    !!    AUTHOR
    !!    ------
    !!     V. Masson 
    !!
    !!    MODIFICATIONS
    !!    -------------
    !!      Original    01/2004
    !!                  10/2005 (G.Pigeon) transfer of domestic heating
    !!      S. Riette   06/2009 Initialisation of XT, XQ, XU and XTKE on canopy levels
    !!      S. Riette   01/2010 Use of interpol_sbl to compute 10m wind diagnostic
    !!      G. Pigeon   09/2012 CCH_BEM, ROUGH_WALL, ROUGH_ROOF for building conv. coef
    !!      G. Pigeon   10/2012 XF_WIN_WIN as arg. of TEB_GARDEN
    !!      B. Decharme 09/2012 New wind implicitation
    !!      J. Escobar  09/2012 KI not allowed without-interface , replace by KI
    !!---------------------------------------------------------------
    !
    !
    USE MODD_CSTS,         ONLY : XRD, XCPD, XP00, XLVTT, XPI, XKARMAN, XG
    USE MODD_SURF_PAR,     ONLY : XUNDEF
    !
    USE MODD_SURF_ATM,     ONLY : CIMPLICIT_WIND
    !
    USE MODD_TEB_n,        ONLY : LGARDEN, LGREENROOF,                                     &
                                  CBEM, TTIME,LCANOPY,CZ0H,CROAD_DIR,CWALL_OPT,            &
                                  XT_CANYON, XQ_CANYON,                                    &
                                  XT_ROOF, XT_ROAD, XT_WALL_A, XT_WALL_B,                  &
                                  XWS_ROOF, XWS_ROAD,                                      &
                                  TSNOW_ROOF, TSNOW_ROAD,                                  &
                                  XH_TRAFFIC, XLE_TRAFFIC, XH_INDUSTRY, XLE_INDUSTRY,      &
                                  XZ0_TOWN, XBLD, XGARDEN, XROAD_DIR, XROAD, XGREENROOF,   &
                                  XBLD_HEIGHT, XWALL_O_HOR, XCAN_HW_RATIO,                 &
                                  XROAD_O_GRND, XGARDEN_O_GRND, XWALL_O_GRND,              &
                                  XALB_ROOF, XEMIS_ROOF, XHC_ROOF,XTC_ROOF, XD_ROOF,       &
                                  XALB_ROAD, XEMIS_ROAD, XHC_ROAD,XTC_ROAD, XD_ROAD,       &
                                  XALB_WALL, XEMIS_WALL, XHC_WALL,XTC_WALL, XD_WALL,       &
                                  XSVF_ROAD, XSVF_WALL,                                    &
                                  XSVF_GARDEN, XWALL_O_BLD,                                &
                                  XQSAT_ROOF, XQSAT_ROAD, XDELT_ROOF, XDELT_ROAD,          &
                                  NTEB_PATCH, XTEB_PATCH, CCH_BEM, XROUGH_ROOF, XROUGH_WALL                       
    !
    USE MODD_BEM_n,        ONLY : XHC_FLOOR, XTC_FLOOR, XD_FLOOR, XTCOOL_TARGET,           &
                                  XTHEAT_TARGET, XF_WASTE_CAN, XEFF_HEAT, XTI_BLD,         &
                                  XT_FLOOR, XT_MASS, XQIN, XQIN_FRAD, XSHGC, XSHGC_SH,     &
                                  XU_WIN, XGR, XINF, CCOOL_COIL, CHEAT_COIL,               &
                                  XF_WATER_COND, XAUX_MAX, XQIN_FLAT,                      &
                                  XHR_TARGET, XT_WIN2, XQI_BLD, XV_VENT, XCAP_SYS_HEAT,    &
                                  XCAP_SYS_RAT, XT_ADP, XM_SYS_RAT, XCOP_RAT, XT_WIN1,     &
                                  XALB_WIN, XABS_WIN, XT_SIZE_MAX, XT_SIZE_MIN, XUGG_WIN,  &
                                  LSHADE, CNATVENT, LSHAD_DAY, LNATVENT_NIGHT,             &
                                  XN_FLOOR, XGLAZ_O_BLD, XMASS_O_BLD, XFLOOR_HW_RATIO,     &
                                  XF_FLOOR_MASS, XF_FLOOR_WALL, XF_FLOOR_WIN,              &
                                  XF_FLOOR_ROOF, XF_WALL_FLOOR, XF_WALL_MASS,              &
                                  XF_WALL_WIN, XF_WIN_FLOOR, XF_WIN_MASS, XF_WIN_WALL,     &
                                  XF_MASS_FLOOR, XF_MASS_WALL, XF_MASS_WIN, &
                                  XTRAN_WIN, XF_WIN_WIN
                                   
    USE MODD_CH_TEB_n,     ONLY : CSV, CCH_DRY_DEP, XDEP, NBEQ, NSV_CHSBEG, NSV_CHSEND,    &
                                  NSV_DSTBEG, NSV_DSTEND, NAEREQ, NDSTEQ, NSLTEQ,          &
                                  NSV_AERBEG, NSV_AEREND, NSV_SLTBEG, NSV_SLTEND
    USE MODD_TEB_CANOPY_n, ONLY : XZ, XU, NLVL, XTKE, XT, XQ,                              &
                                  XLMO, XLM, XLEPS,XZF, XDZ, XDZF, XP
    USE MODD_DIAG_TEB_n,   ONLY : N2M, XZON10M, XMER10M
    USE MODD_DIAG_UTCI_TEB_n, ONLY : LUTCI, XUTCI_IN, XUTCI_OUTSUN,          &
                                     XUTCI_OUTSHADE, XTRAD_SUN, XTRAD_SHADE
    USE MODD_DST_n,        ONLY : XEMISRADIUS_DST, XEMISSIG_DST
    USE MODD_SLT_n,        ONLY : XEMISRADIUS_SLT, XEMISSIG_SLT
    USE MODD_DST_SURF
    USE MODD_SLT_SURF
    !
    !
    USE MODE_DSLT_SURF
    USE MODE_THERMOS
    USE MODE_SBLS
    !
    USE MODI_GOTO_TEB
    USE MODI_AVERAGE_RAD
    USE MODI_SM10
    USE MODI_ADD_FORECAST_TO_DATE_SURF
    USE MODI_DIAG_INLINE_TEB_n
    USE MODI_DIAG_MISC_TEB_n
    USE MODI_CH_AER_DEP
    USE MODI_CH_DEP_TOWN
    USE MODI_DSLT_DEP
    USE MODI_TEB_GARDEN
    USE MODI_TEB_CANOPY
    ! 
    USE YOMHOOK   ,ONLY : LHOOK,   DR_HOOK
    USE PARKIND1  ,ONLY : JPRB
    !
    USE MODI_ABOR1_SFX
    USE MODI_CANOPY_EVOL
    USE MODI_CANOPY_GRID_UPDATE
    USE MODI_UTCI_TEB
    USE MODI_CIRCUMSOLAR_RAD
    !
    IMPLICIT NONE
    !
    !*      0.1    declarations of arguments
    !
     CHARACTER(LEN=6),    INTENT(IN)  :: HPROGRAM  ! program calling surf. schemes
     CHARACTER(LEN=1),    INTENT(IN)  :: HCOUPLING ! type of coupling
                                                  ! 'E' : explicit
                                                  ! 'I' : implicit
    INTEGER,             INTENT(IN)  :: KYEAR     ! current year (UTC)
    INTEGER,             INTENT(IN)  :: KMONTH    ! current month (UTC)
    INTEGER,             INTENT(IN)  :: KDAY      ! current day (UTC)
    REAL,                INTENT(IN)  :: PTIME     ! current time since midnight (UTC, s)
    INTEGER,             INTENT(IN)  :: KI        ! number of points
    INTEGER,             INTENT(IN)  :: KSV       ! number of scalars
    INTEGER,             INTENT(IN)  :: KSW       ! number of short-wave spectral bands
    REAL, DIMENSION(KI), INTENT(IN)  :: PTSUN     ! solar time                    (s from midnight)
    REAL,                INTENT(IN)  :: PTSTEP    ! atmospheric time-step                 (s)
    REAL, DIMENSION(KI), INTENT(IN)  :: PZREF     ! height of T,q forcing                 (m)
    REAL, DIMENSION(KI), INTENT(IN)  :: PUREF     ! height of wind forcing                (m)
    !
    REAL, DIMENSION(KI), INTENT(IN)  :: PTA       ! air temperature forcing               (K)
    REAL, DIMENSION(KI), INTENT(IN)  :: PQA       ! air humidity forcing                  (kg/m3)
    REAL, DIMENSION(KI), INTENT(IN)  :: PRHOA     ! air density                           (kg/m3)
    REAL, DIMENSION(KI,KSV),INTENT(IN) :: PSV     ! scalar variables
    !                                             ! chemistry:   first char. in HSV: '#'  (molecule/m3)
    !                                             !
     CHARACTER(LEN=6), DIMENSION(KSV),INTENT(IN):: HSV  ! name of all scalar variables
    REAL, DIMENSION(KI), INTENT(IN)  :: PU        ! zonal wind                            (m/s)
    REAL, DIMENSION(KI), INTENT(IN)  :: PV        ! meridian wind                         (m/s)
    REAL, DIMENSION(KI,KSW),INTENT(IN) :: PDIR_SW ! direct  solar radiation (on horizontal surf.)
    !                                             !                                       (W/m2)
    REAL, DIMENSION(KI,KSW),INTENT(IN) :: PSCA_SW ! diffuse solar radiation (on horizontal surf.)
    !                                             !                                       (W/m2)
    REAL, DIMENSION(KSW),INTENT(IN)  :: PSW_BANDS ! mean wavelength of each shortwave band (m)
    REAL, DIMENSION(KI), INTENT(IN)  :: PZENITH   ! zenithal angle       (radian from the vertical)
    REAL, DIMENSION(KI), INTENT(IN)  :: PAZIM     ! azimuthal angle      (radian from North, clockwise)
    REAL, DIMENSION(KI), INTENT(IN)  :: PLW       ! longwave radiation (on horizontal surf.)
    !                                             !                                       (W/m2)
    REAL, DIMENSION(KI), INTENT(IN)  :: PPS       ! pressure at atmospheric model surface (Pa)
    REAL, DIMENSION(KI), INTENT(IN)  :: PPA       ! pressure at forcing level             (Pa)
    REAL, DIMENSION(KI), INTENT(IN)  :: PZS       ! atmospheric model orography           (m)
    REAL, DIMENSION(KI), INTENT(IN)  :: PCO2      ! CO2 concentration in the air          (kg/m3)
    REAL, DIMENSION(KI), INTENT(IN)  :: PSNOW     ! snow precipitation                    (kg/m2/s)
    REAL, DIMENSION(KI), INTENT(IN)  :: PRAIN     ! liquid precipitation                  (kg/m2/s)
    !
    !
    REAL, DIMENSION(KI), INTENT(OUT) :: PSFTH     ! flux of heat                          (W/m2)
    REAL, DIMENSION(KI), INTENT(OUT) :: PSFTQ     ! flux of water vapor                   (kg/m2/s)
    REAL, DIMENSION(KI), INTENT(OUT) :: PSFU      ! zonal momentum flux                   (Pa)
    REAL, DIMENSION(KI), INTENT(OUT) :: PSFV      ! meridian momentum flux                (Pa)
    REAL, DIMENSION(KI), INTENT(OUT) :: PSFCO2    ! flux of CO2                           (kg/m2/s)
    REAL, DIMENSION(KI,KSV),INTENT(OUT):: PSFTS   ! flux of scalar var.                   (kg/m2/s)
    !
    REAL, DIMENSION(KI), INTENT(OUT) :: PTRAD     ! radiative temperature                 (K)
    REAL, DIMENSION(KI,KSW),INTENT(OUT):: PDIR_ALB! direct albedo for each spectral band  (-)
    REAL, DIMENSION(KI,KSW),INTENT(OUT):: PSCA_ALB! diffuse albedo for each spectral band (-)
    REAL, DIMENSION(KI), INTENT(OUT) :: PEMIS     ! emissivity                            (-)
    !
    REAL, DIMENSION(KI), INTENT(IN) :: PPEW_A_COEF! implicit coefficients
    REAL, DIMENSION(KI), INTENT(IN) :: PPEW_B_COEF! needed if HCOUPLING='I'
    REAL, DIMENSION(KI), INTENT(IN) :: PPET_A_COEF
    REAL, DIMENSION(KI), INTENT(IN) :: PPEQ_A_COEF
    REAL, DIMENSION(KI), INTENT(IN) :: PPET_B_COEF
    REAL, DIMENSION(KI), INTENT(IN) :: PPEQ_B_COEF
     CHARACTER(LEN=2),    INTENT(IN) :: HTEST ! must be equal to 'OK'
    !
    !
    !*      0.2    declarations of local variables
    !
    INTEGER                     :: JSWB        ! loop counter on shortwave spectral bands
    !         
    REAL, DIMENSION(KI)  :: ZQA         ! specific humidity                 (kg/kg)
    REAL, DIMENSION(KI)  :: ZEXNA       ! Exner function at forcing level
    REAL, DIMENSION(KI)  :: ZEXNS       ! Exner function at surface level
    REAL, DIMENSION(KI)  :: ZWIND       ! wind
    !
    ! Ouput Diagnostics:
    !
    REAL, DIMENSION(KI)  :: ZU_CANYON   ! wind in canyon
    REAL, DIMENSION(KI)  :: ZT_CANYON   ! temperature in canyon
    REAL, DIMENSION(KI)  :: ZQ_CANYON   ! specific humidity in canyon
    REAL, DIMENSION(KI)  :: ZT_CAN      ! temperature in canyon       (evolving in TEB)
    REAL, DIMENSION(KI)  :: ZQ_CAN      ! specific humidity in canyon (evolving in TEB)
    !
    REAL, DIMENSION(KI)  :: ZRN_ROOF    ! net radiation on roof
    REAL, DIMENSION(KI)  :: ZH_ROOF     ! sensible heat flux on roof
    REAL, DIMENSION(KI)  :: ZLE_ROOF    ! latent heat flux on roof
    REAL, DIMENSION(KI)  :: ZLEW_ROOF   ! latent heat flux on snowfree roof
    REAL, DIMENSION(KI)  :: ZGFLUX_ROOF ! storage flux in roof
    REAL, DIMENSION(KI)  :: ZRUNOFF_ROOF! water runoff from roof
    REAL, DIMENSION(KI)  :: ZRN_ROAD    ! net radiation on road
    REAL, DIMENSION(KI)  :: ZH_ROAD     ! sensible heat flux on road
    REAL, DIMENSION(KI)  :: ZLE_ROAD    ! latent heat flux on road
    REAL, DIMENSION(KI)  :: ZLEW_ROAD   ! latent heat flux on snowfree road
    REAL, DIMENSION(KI)  :: ZGFLUX_ROAD ! storage flux in road
    REAL, DIMENSION(KI)  :: ZRUNOFF_ROAD! water runoff from road
    REAL, DIMENSION(KI)  :: ZRN_WALL_A  ! net radiation on walls
    REAL, DIMENSION(KI)  :: ZH_WALL_A   ! sensible heat flux on walls
    REAL, DIMENSION(KI)  :: ZLE_WALL_A  ! latent heat flux on walls
    REAL, DIMENSION(KI)  :: ZGFLUX_WALL_A!storage flux in walls
    REAL, DIMENSION(KI)  :: ZRN_WALL_B  ! net radiation on walls
    REAL, DIMENSION(KI)  :: ZH_WALL_B   ! sensible heat flux on walls
    REAL, DIMENSION(KI)  :: ZLE_WALL_B  ! latent heat flux on walls
    REAL, DIMENSION(KI)  :: ZGFLUX_WALL_B!storage flux in walls
    REAL, DIMENSION(KI)  :: ZRN_GARDEN  ! net radiation on green areas
    REAL, DIMENSION(KI)  :: ZH_GARDEN   ! sensible heat flux on green areas
    REAL, DIMENSION(KI)  :: ZLE_GARDEN  ! latent heat flux on green areas
    REAL, DIMENSION(KI)  :: ZGFLUX_GARDEN!storage flux in green areas
    REAL, DIMENSION(KI)  :: ZRN_GREENROOF! net radiation on green roofs
    REAL, DIMENSION(KI)  :: ZH_GREENROOF ! sensible heat flux on green roofs
    REAL, DIMENSION(KI)  :: ZLE_GREENROOF! latent heat flux on green roofs
    REAL, DIMENSION(KI)  :: ZGFLUX_GREENROOF    ! storage flux in green roofs
    REAL, DIMENSION(KI)  :: ZG_GREENROOF_ROOF   ! heat flux between base of greenroof
    REAL, DIMENSION(KI)  :: ZRUNOFF_GREENROOF   ! water runoff from green roof
    REAL, DIMENSION(KI)  :: ZDRAIN_GREENROOF    ! water drainage from green roof
    REAL, DIMENSION(KI)  :: ZRN_STRLROOF        ! net radiation on structural roof
    REAL, DIMENSION(KI)  :: ZH_STRLROOF         ! sensible heat flux on structural roof
    REAL, DIMENSION(KI)  :: ZLE_STRLROOF        ! latent heat flux on structural roof
    REAL, DIMENSION(KI)  :: ZGFLUX_STRLROOF     ! storage flux in structural roof
    REAL, DIMENSION(KI)  :: ZRN_BLT     ! net radiation on built surf 
    REAL, DIMENSION(KI)  :: ZH_BLT      ! sensible heat flux on built surf 
    REAL, DIMENSION(KI)  :: ZLE_BLT     ! latent heat flux on built surf 
    REAL, DIMENSION(KI)  :: ZGFLUX_BLT  ! storage flux in built surf 
    REAL, DIMENSION(KI)  :: ZRN_GRND    ! net radiation on ground built surf
    REAL, DIMENSION(KI)  :: ZH_GRND     ! sensible heat flux on ground built surf
    REAL, DIMENSION(KI)  :: ZLE_GRND    ! latent heat flux on ground built surf
    REAL, DIMENSION(KI)  :: ZGFLUX_GRND ! storage flux in ground built surf
    REAL, DIMENSION(KI)  :: ZRNSNOW_ROOF  ! net radiation over snow
    REAL, DIMENSION(KI)  :: ZHSNOW_ROOF   ! sensible heat flux over snow
    REAL, DIMENSION(KI)  :: ZLESNOW_ROOF  ! latent heat flux over snow
    REAL, DIMENSION(KI)  :: ZGSNOW_ROOF   ! flux under the snow
    REAL, DIMENSION(KI)  :: ZMELT_ROOF    ! snow melt
    REAL, DIMENSION(KI)  :: ZRNSNOW_ROAD  ! net radiation over snow
    REAL, DIMENSION(KI)  :: ZHSNOW_ROAD   ! sensible heat flux over snow
    REAL, DIMENSION(KI)  :: ZLESNOW_ROAD  ! latent heat flux over snow
    REAL, DIMENSION(KI)  :: ZGSNOW_ROAD   ! flux under the snow
    REAL, DIMENSION(KI)  :: ZMELT_ROAD    ! snow melt
    !
    REAL, DIMENSION(KI)  :: ZTRAD         ! radiative temperature for current patch
    REAL, DIMENSION(KI)  :: ZEMIS         ! emissivity for current patch
    REAL, DIMENSION(KI,NTEB_PATCH) :: ZTRAD_PATCH ! radiative temperature for each patch
    REAL, DIMENSION(KI,NTEB_PATCH) :: ZEMIS_PATCH ! emissivity for each patch
    REAL, DIMENSION(KI,KSW,NTEB_PATCH) :: ZDIR_ALB_PATCH ! direct albedo per wavelength and patch
    REAL, DIMENSION(KI,KSW,NTEB_PATCH) :: ZSCA_ALB_PATCH ! diffuse albedo per wavelength and patch
    !
    REAL, DIMENSION(KI)  :: ZRN           ! net radiation over town
    REAL, DIMENSION(KI)  :: ZH            ! sensible heat flux over town
    REAL, DIMENSION(KI)  :: ZLE           ! latent heat flux over town
    REAL, DIMENSION(KI)  :: ZGFLUX        ! flux through the ground
    REAL, DIMENSION(KI)  :: ZSFCO2        ! CO2 flux over town
    REAL, DIMENSION(KI)  :: ZQF_BLD       ! domestic heating
    REAL, DIMENSION(KI)  :: ZFLX_BLD      ! flux from bld
    REAL, DIMENSION(KI)  :: ZDQS_TOWN     ! storage inside town materials
    REAL, DIMENSION(KI)  :: ZQF_TOWN      ! total anthropogenic heat
    REAL, DIMENSION(KI)  :: ZEVAP         ! evaporation (km/m2/s)
    REAL, DIMENSION(KI)  :: ZRUNOFF       ! runoff over the ground
    REAL, DIMENSION(KI)  :: ZCD           ! drag coefficient
    REAL, DIMENSION(KI)  :: ZCDN          ! neutral drag coefficient
    REAL, DIMENSION(KI)  :: ZCH           ! heat drag
    REAL, DIMENSION(KI)  :: ZRI           ! Richardson number
    REAL, DIMENSION(KI)  :: ZUW_GRND      ! momentum flux for ground built surf
    REAL, DIMENSION(KI)  :: ZUW_ROOF      ! momentum flux for roofs
    REAL, DIMENSION(KI)  :: ZDUWDU_GRND   !
    REAL, DIMENSION(KI)  :: ZDUWDU_ROOF   !
    REAL, DIMENSION(KI)  :: ZUSTAR        ! friction velocity
    REAL, DIMENSION(KI)  :: ZSFU          ! momentum flux for patch (U direction)
    REAL, DIMENSION(KI)  :: ZSFV          ! momentum flux for patch (V direction)
    REAL, DIMENSION(KI)  :: ZAVG_DIR_ALB  ! direct albedo of town
    REAL, DIMENSION(KI)  :: ZAVG_SCA_ALB  ! diffuse albedo of town
    REAL, DIMENSION(KI)  :: ZAVG_T_CANYON ! temperature in canyon for town 
    REAL, DIMENSION(KI)  :: ZAVG_Q_CANYON ! specific humidity in canyon for town
    !
    REAL, DIMENSION(KI)  :: ZAVG_CD       ! aggregated drag coefficient
    REAL, DIMENSION(KI)  :: ZAVG_CDN      ! aggregated neutral drag coefficient
    REAL, DIMENSION(KI)  :: ZAVG_RI       ! aggregated Richardson number
    REAL, DIMENSION(KI)  :: ZAVG_CH       ! aggregated Heat transfer coefficient
    !
    REAL, DIMENSION(KI)  :: ZDIR_ALB      ! direct albedo of town
    REAL, DIMENSION(KI)  :: ZSCA_ALB      ! diffuse albedo of town
    !
    REAL, DIMENSION(KI)  :: ZH_TRAFFIC    ! anthropogenic sensible
    !                                            ! heat fluxes due to traffic
    REAL, DIMENSION(KI)  :: ZLE_TRAFFIC   ! anthropogenic latent
    !                                            ! heat fluxes due to traffic
    REAL, DIMENSION(KI)  :: ZRESA_TOWN    ! aerodynamical resistance
    REAL, DIMENSION(KI)  :: ZAC_ROAD      ! road aerodynamical conductance
    REAL, DIMENSION(KI)  :: ZAC_GARDEN    ! green area aerodynamical conductance
    REAL, DIMENSION(KI)  :: ZAC_GRND      ! ground built surf aerodynamical conductance
    REAL, DIMENSION(KI)  :: ZAC_GREENROOF ! green roof aerodynamical conductance
    REAL, DIMENSION(KI)  :: ZAC_ROAD_WAT  ! road water aerodynamical conductance
    REAL, DIMENSION(KI)  :: ZAC_GARDEN_WAT! green area water aerodynamical conductance
    REAL, DIMENSION(KI)  :: ZAC_GRND_WAT  ! ground built surf water aerodynamical conductance
    REAL, DIMENSION(KI)  :: ZAC_GREENROOF_WAT! green roof water aerodynamical conductance
    REAL, DIMENSION(KI,1):: ZESNOW_GARDEN    ! green area snow emissivity
    !
    REAL                        :: ZBEGIN_TRAFFIC_TIME ! start traffic time (solar time, s)
    REAL                        :: ZEND_TRAFFIC_TIME   ! end traffic time   (solar time, s)
    REAL, DIMENSION(KI)  :: ZDIR_SW       ! total direct SW
    REAL, DIMENSION(KI)  :: ZSCA_SW       ! total diffuse SW
    REAL, DIMENSION(KI)  :: ZPEW_A_COEF   ! implicit coefficients
    REAL, DIMENSION(KI)  :: ZPEW_B_COEF   ! needed if HCOUPLING='I'
    
    !***** CANOPY  *****
    REAL, DIMENSION(KI)        :: ZSFLUX_U  ! Surface flux u'w' (m2/s2)
    REAL, DIMENSION(KI)        :: ZSFLUX_T  ! Surface flux w'T' (mK/s)
    REAL, DIMENSION(KI)        :: ZSFLUX_Q  ! Surface flux w'q' (kgm2/s)
    REAL, DIMENSION(KI,NLVL)   :: ZFORC_U   ! tendency due to drag force for wind
    REAL, DIMENSION(KI,NLVL)   :: ZDFORC_UDU! formal derivative of
    !                                              ! tendency due to drag force for wind
    REAL, DIMENSION(KI,NLVL)   :: ZFORC_E   ! tendency due to drag force for TKE
    REAL, DIMENSION(KI,NLVL)   :: ZDFORC_EDE! formal derivative of
    !                                              ! tendency due to drag force for TKE
    REAL, DIMENSION(KI,NLVL)   :: ZFORC_T   ! tendency due to drag force for Temp
    REAL, DIMENSION(KI,NLVL)   :: ZDFORC_TDT! formal derivative of
    !                                              ! tendency due to drag force for Temp
    REAL, DIMENSION(KI,NLVL)   :: ZFORC_Q   ! tendency due to drag force for hum
    REAL, DIMENSION(KI,NLVL)   :: ZDFORC_QDQ! formal derivative of
    !                                              ! tendency due to drag force for hum.
    
    REAL, DIMENSION(KI)        :: ZAVG_UW_GRND
    REAL, DIMENSION(KI)        :: ZAVG_DUWDU_GRND
    REAL, DIMENSION(KI)        :: ZAVG_UW_ROOF
    REAL, DIMENSION(KI)        :: ZAVG_DUWDU_ROOF
    REAL, DIMENSION(KI)        :: ZAVG_H_GRND
    REAL, DIMENSION(KI)        :: ZAVG_H_WALL
    REAL, DIMENSION(KI)        :: ZAVG_H_ROOF
    REAL, DIMENSION(KI)        :: ZAVG_E_GRND
    REAL, DIMENSION(KI)        :: ZAVG_E_ROOF
    REAL, DIMENSION(KI)        :: ZAVG_AC_GRND
    REAL, DIMENSION(KI)        :: ZAVG_AC_GRND_WAT
    REAL, DIMENSION(KI)        :: ZAVG_Z0_TOWN
    REAL, DIMENSION(KI)        :: ZAVG_RESA_TOWN
    REAL, DIMENSION(KI)        :: ZAVG_USTAR        ! town avegared Ustar
    REAL, DIMENSION(KI)        :: ZAVG_BLD          ! town averaged building fraction
    REAL, DIMENSION(KI)        :: ZAVG_BLD_HEIGHT   ! town averaged building height
    REAL, DIMENSION(KI)        :: ZAVG_WALL_O_HOR   ! town averaged Wall/hor ratio
    REAL, DIMENSION(KI)        :: ZAVG_CAN_HW_RATIO ! town averaged road aspect ratio
    REAL, DIMENSION(KI)        :: ZAVG_H
    REAL, DIMENSION(KI)        :: ZAVG_LE
    REAL, DIMENSION(KI)        :: ZAVG_RN
    REAL, DIMENSION(KI)        :: ZAVG_GFLUX
    REAL, DIMENSION(KI)        :: ZAVG_REF_SW_GRND
    REAL, DIMENSION(KI)        :: ZAVG_REF_SW_FAC
    REAL, DIMENSION(KI)        :: ZAVG_SCA_SW
    REAL, DIMENSION(KI)        :: ZAVG_DIR_SW 
    REAL, DIMENSION(KI)        :: ZAVG_EMIT_LW_FAC
    REAL, DIMENSION(KI)        :: ZAVG_EMIT_LW_GRND
    REAL, DIMENSION(KI)        :: ZAVG_T_RAD_IND
    REAL, DIMENSION(KI)        :: ZT_LOWCAN  ! temperature at lowest canyon level (K)
    REAL, DIMENSION(KI)        :: ZQ_LOWCAN  ! humidity    at lowest canyon level (kg/kg)
    REAL, DIMENSION(KI)        :: ZU_LOWCAN  ! wind        at lowest canyon level (m/s)
    REAL, DIMENSION(KI)        :: ZZ_LOWCAN  ! height      of lowest canyon level (m)
    REAL, DIMENSION(KI)        :: ZPEW_A_COEF_LOWCAN   ! implicit coefficients for wind coupling
    REAL, DIMENSION(KI)        :: ZPEW_B_COEF_LOWCAN   ! between first canopy level and road
    REAL, DIMENSION(KI)        :: ZTA        ! temperature at canyon level just above roof (K)
    REAL, DIMENSION(KI)        :: ZPA        ! pressure    at canyon level just above roof (K)
    REAL, DIMENSION(KI)        :: ZUA        ! wind        at canyon level just above roof (m/s)
    REAL, DIMENSION(KI)        :: ZUREF      ! height      of canyon level just above roof (m)
    REAL, DIMENSION(KI)        :: ZZREF      ! height      of canyon level just above roof (m)
    REAL, DIMENSION(KI)        :: ZLAMBDA_F  ! frontal density (-)
    REAL, DIMENSION(KI)        :: ZLMO       ! Monin-Obukhov length at canopy height (m)
    REAL, DIMENSION(KI,NLVL)   :: ZL         ! Mixing length generic profile at mid levels
    !
    ! absorbed solar and infra-red radiation by road, wall and roof
    !                                                      
    REAL, DIMENSION(KI) :: ZABS_SW_ROAD
    REAL, DIMENSION(KI) :: ZABS_SW_WALL_A
    REAL, DIMENSION(KI) :: ZABS_SW_WALL_B
    REAL, DIMENSION(KI) :: ZABS_SW_ROOF
    REAL, DIMENSION(KI) :: ZABS_SW_GARDEN
    REAL, DIMENSION(KI) :: ZABS_SW_GREENROOF
    REAL, DIMENSION(KI) :: ZABS_SW_SNOW_ROAD
    REAL, DIMENSION(KI) :: ZABS_SW_SNOW_ROOF
    REAL, DIMENSION(KI) :: ZABS_LW_SNOW_ROAD
    REAL, DIMENSION(KI) :: ZABS_LW_SNOW_ROOF
    REAL, DIMENSION(KI) :: ZABS_LW_ROAD
    REAL, DIMENSION(KI) :: ZABS_LW_WALL_A
    REAL, DIMENSION(KI) :: ZABS_LW_WALL_B
    REAL, DIMENSION(KI) :: ZABS_LW_ROOF
    REAL, DIMENSION(KI) :: ZABS_LW_GARDEN 
    REAL, DIMENSION(KI) :: ZABS_LW_GREENROOF
    !
    REAL, DIMENSION(KI)        :: ZU_UTCI ! wind speed for the UTCI calculation (m/s) 
    
    REAL, DIMENSION(KI)        :: ZALFAU   ! V+(1) = alfa u'w'(1) + beta
    REAL, DIMENSION(KI)        :: ZBETAU   ! V+(1) = alfa u'w'(1) + beta
    REAL, DIMENSION(KI)        :: ZALFAT   ! Th+(1) = alfa w'th'(1) + beta
    REAL, DIMENSION(KI)        :: ZBETAT   ! Th+(1) = alfa w'th'(1) + beta
    REAL, DIMENSION(KI)        :: ZALFAQ   ! Q+(1) = alfa w'q'(1) + beta
    REAL, DIMENSION(KI)        :: ZBETAQ   ! Q+(1) = alfa w'q'(1) + beta
    !***** CANOPY  *****
    REAL, DIMENSION(KI)        :: ZWAKE      ! reduction of average wind speed
    !                                              ! in canyon due to direction average.
    ! new local variables after BEM
    !
    REAL, DIMENSION(KI) :: ZCAP_SYS
    REAL, DIMENSION(KI) :: ZM_SYS
    REAL, DIMENSION(KI) :: ZCOP
    REAL, DIMENSION(KI) :: ZQ_SYS
    REAL, DIMENSION(KI) :: ZT_SYS
    REAL, DIMENSION(KI) :: ZTR_SW_WIN
    REAL, DIMENSION(KI) :: ZFAN_POWER
    REAL, DIMENSION(KI) :: ZABS_SW_WIN
    REAL, DIMENSION(KI) :: ZABS_LW_WIN
    REAL, DIMENSION(KI) :: ZH_BLD_COOL
    REAL, DIMENSION(KI) :: ZT_BLD_COOL
    REAL, DIMENSION(KI) :: ZH_BLD_HEAT
    REAL, DIMENSION(KI) :: ZLE_BLD_COOL
    REAL, DIMENSION(KI) :: ZLE_BLD_HEAT  
    REAL, DIMENSION(KI) :: ZH_WASTE
    REAL, DIMENSION(KI) :: ZLE_WASTE
    REAL, DIMENSION(KI) :: ZHVAC_COOL  
    REAL, DIMENSION(KI) :: ZHVAC_HEAT
    
    !new local variables for UTCI calculation
    REAL, DIMENSION(KI) :: ZEMIT_LW_GRND
    REAL, DIMENSION(KI) :: ZEMIT_LW_FAC
    REAL, DIMENSION(KI) :: ZT_RAD_IND   ! Indoor mean radiant temperature [K]
    REAL, DIMENSION(KI) :: ZREF_SW_GRND ! total solar rad reflected from ground
    REAL, DIMENSION(KI) :: ZREF_SW_FAC  ! total solar rad reflected from facade
    REAL, DIMENSION(KI) :: ZHU_BLD
    REAL, DIMENSION(KI) :: ZAVG_TI_BLD
    REAL, DIMENSION(KI) :: ZAVG_QI_BLD
    REAL, DIMENSION(KI) :: ZF1_o_B
    REAL, DIMENSION(KI,SIZE(PDIR_SW,2))  :: ZDIR_SWB ! total direct SW per band
    REAL, DIMENSION(KI,SIZE(PSCA_SW,2))  :: ZSCA_SWB ! total diffuse SW per band
    !
    REAL, DIMENSION(KI)        :: ZCOEF
    !
    REAL                       :: ZCONVERTFACM0_SLT, ZCONVERTFACM0_DST
    REAL                       :: ZCONVERTFACM3_SLT, ZCONVERTFACM3_DST
    REAL                       :: ZCONVERTFACM6_SLT, ZCONVERTFACM6_DST
    !
    INTEGER                           :: JI
    INTEGER                           :: JLAYER
    INTEGER                           :: JJ
    !
    ! number of TEB patches
    !
    INTEGER                    :: JTEB_PATCH ! loop counter
    !
    REAL(KIND=JPRB) :: ZHOOK_HANDLE
    !
    !-------------------------------------------------------------------------------------
    ! Preliminaries:
    !-------------------------------------------------------------------------------------
    IF (LHOOK) CALL DR_HOOK('COUPLING_TEB_N',0,ZHOOK_HANDLE)
    IF (HTEST/='OK') THEN
      CALL ABOR1_SFX('COUPLING_TEBN: FATAL ERROR DURING ARGUMENT TRANSFER')
    END IF
    
    !-------------------------------------------------------------------------------------
    !
    ! scalar fluxes
    !
    PSFTS(:,:) = 0.
    !
    ! broadband radiative fluxes
    !
    ZDIR_SW(:) = 0.
    ZSCA_SW(:) = 0.
    DO JSWB=1,KSW
      !add directionnal contrib from scattered radiation
      CALL CIRCUMSOLAR_RAD(PDIR_SW(:,JSWB), PSCA_SW(:,JSWB), PZENITH, ZF1_o_B)
      ZDIR_SWB(:,JSWB) = PDIR_SW(:,JSWB) + PSCA_SW(:,JSWB) * ZF1_o_B
      ZSCA_SWB(:,JSWB) = PSCA_SW(:,JSWB) * (1. - ZF1_o_B)
      !add directionnal contrib from scattered radiation
      DO JJ=1,SIZE(PDIR_SW,1)
        ZDIR_SW(JJ) = ZDIR_SW(JJ) + ZDIR_SWB(JJ,JSWB)
        ZSCA_SW(JJ) = ZSCA_SW(JJ) + ZSCA_SWB(JJ,JSWB)
      ENDDO
    END DO
    !
    DO JJ=1,KI
    ! specific humidity (conversion from kg/m3 to kg/kg)
    !
      ZQA(JJ) = PQA(JJ) / PRHOA(JJ)
    !
    ! wind
    !
      ZWIND(JJ) = SQRT(PU(JJ)**2+PV(JJ)**2)
    !
    ENDDO
    ! method of wind coupling
    !
    IF (HCOUPLING=='I') THEN
      ZPEW_A_COEF = PPEW_A_COEF
      ZPEW_B_COEF = PPEW_B_COEF
    ELSE
      ZPEW_A_COEF =  0.
      ZPEW_B_COEF =  ZWIND
    END IF
    !
    !
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    ! Time evolution
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    !
    TTIME%TIME = TTIME%TIME + PTSTEP
     CALL ADD_FORECAST_TO_DATE_SURF(TTIME%TDATE%YEAR,TTIME%TDATE%MONTH,TTIME%TDATE%DAY,TTIME%TIME)
    !
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    !  Anthropogenic fluxes (except building heating)
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    !
    ZBEGIN_TRAFFIC_TIME = 21600.
    ZEND_TRAFFIC_TIME   = 64800.
    !
    WHERE(       PTSUN>ZBEGIN_TRAFFIC_TIME   &
          .AND.  PTSUN<ZEND_TRAFFIC_TIME     )
      ZH_TRAFFIC  (:) = XH_TRAFFIC   (:)
      ZLE_TRAFFIC (:) = XLE_TRAFFIC  (:)
    ELSEWHERE
      ZH_TRAFFIC  (:) = 0.
      ZLE_TRAFFIC (:) = 0.   
    END WHERE
    !
    !--------------------------------------------------------------------------------------
    !  Canyon forcing for TEB
    !--------------------------------------------------------------------------------------
    !-------------------------------------------------------------------------------------
    ! Town averaged quantities to force canopy atmospheric layers
    !-------------------------------------------------------------------------------------
    
    DO JTEB_PATCH=1,NTEB_PATCH
      CALL GOTO_TEB(JTEB_PATCH)
      CALL ADD_PATCH_CONTRIB(JTEB_PATCH,ZAVG_BLD,         XBLD         )
      CALL ADD_PATCH_CONTRIB(JTEB_PATCH,ZAVG_BLD_HEIGHT,  XBLD_HEIGHT  )
      CALL ADD_PATCH_CONTRIB(JTEB_PATCH,ZAVG_WALL_O_HOR,  XWALL_O_HOR  )
      CALL ADD_PATCH_CONTRIB(JTEB_PATCH,ZAVG_CAN_HW_RATIO,XCAN_HW_RATIO)
      CALL ADD_PATCH_CONTRIB(JTEB_PATCH,ZAVG_Z0_TOWN,     XZ0_TOWN     )
    END DO
    !
    IF (LCANOPY) THEN
    !-------------------------------------------------------------------------------------
    ! Updates canopy vertical grid as a function of forcing height
    !-------------------------------------------------------------------------------------
    !
    !* determines where is the forcing level and modifies the upper levels of the canopy grid
    !
      CALL CANOPY_GRID_UPDATE(KI,NLVL,ZAVG_BLD_HEIGHT,ZAVG_BLD_HEIGHT+PUREF,XZ,XZF,XDZ,XDZF)
    !
    !* Initialisations of T, Q, TKE and wind at first time step
    !
    
      IF(ANY(XT(:,:) == XUNDEF)) THEN
        DO JLAYER=1,NLVL
          XT(:,JLAYER) = PTA(:)
          XQ(:,JLAYER) = PQA(:)
          XU(:,JLAYER) = 2./XPI * ZWIND(:)                                  &
                  * LOG( (          2.* XBLD_HEIGHT(:)/3.) / XZ0_TOWN(:))   &
                  / LOG( (PUREF(:)+ 2.* XBLD_HEIGHT(:)/3.) / XZ0_TOWN(:))
        END  DO
        XTKE(:,:) = 1.
      ENDIF
    !
    !* default forcing above roof: forcing level
    ZUREF(:)     = PUREF(:)
    ZZREF(:)     = PZREF(:)
    ZUA(:)       = XU(:,NLVL)
    ZTA(:)       = XT(:,NLVL)
    ZQA(:)       = XQ(:,NLVL)/PRHOA(:)
    ZPA(:)       = XP(:,NLVL)
    !* for the time being, only one value is kept for wall in-canyon forcing, in the middle of the canyon
    ZU_CANYON(:) = ZUA(:)
    ZT_CANYON(:) = ZTA(:)
    ZQ_CANYON(:) = ZQA(:)
      DO JLAYER=1,NLVL-1
        DO JI=1,KI
          !* finds middle canyon layer
          IF (XZ(JI,JLAYER)<ZAVG_BLD_HEIGHT(JI)/2. .AND. XZ(JI,JLAYER+1)>=ZAVG_BLD_HEIGHT(JI)/2.) THEN
            ZCOEF(JI) = (ZAVG_BLD_HEIGHT(JI)/2.-XZ(JI,JLAYER))/(XZ(JI,JLAYER+1)-XZ(JI,JLAYER))
            ZU_CANYON(JI) = XU(JI,JLAYER) + ZCOEF(JI) * (XU(JI,JLAYER+1)-XU(JI,JLAYER))
            ZT_CANYON(JI) = XT(JI,JLAYER) + ZCOEF(JI) * (XT(JI,JLAYER+1)-XT(JI,JLAYER))
            ZQ_CANYON(JI) =(XQ(JI,JLAYER) + ZCOEF(JI) * (XQ(JI,JLAYER+1)-XQ(JI,JLAYER)))/PRHOA(JI)
          END IF
          !* finds layer just above roof (at least 1m above roof)
          IF (XZ(JI,JLAYER)<ZAVG_BLD_HEIGHT(JI)+1. .AND. XZ(JI,JLAYER+1)>=ZAVG_BLD_HEIGHT(JI)+1.) THEN
            ZUREF(JI) = XZ(JI,JLAYER+1) - ZAVG_BLD_HEIGHT(JI)
            ZZREF(JI) = XZ(JI,JLAYER+1) - ZAVG_BLD_HEIGHT(JI)
            ZTA  (JI) = XT(JI,JLAYER+1)
            ZQA  (JI) = XQ(JI,JLAYER+1)/PRHOA(JI)
            !ZUA  (JI) = XU(JI,JLAYER+1)
            ZUA  (JI) = MAX(XU(JI,JLAYER+1) - 2.*SQRT(XTKE(JI,JLAYER+1)) , XU(JI,JLAYER+1)/3.)
            ZPA  (JI) = XP(JI,JLAYER+1)
            ZLMO (JI) = XLMO(JI,JLAYER+1)
          END IF
        END DO
      END DO
      ZU_CANYON= MAX(ZU_CANYON,0.2)
      ZU_LOWCAN=XU(:,1)
      ZT_LOWCAN=XT(:,1)
      ZQ_LOWCAN=XQ(:,1) / PRHOA(:)
      ZZ_LOWCAN=XZ(:,1)
      WHERE(ZPA==XUNDEF) ZPA = PPA   ! security for first time step
    !
    !-------------------------------------------------------------------------------------
    ! determine the vertical profile for mixing and dissipative lengths (at full levels)
    !-------------------------------------------------------------------------------------
    !
    ! frontal density
      ZLAMBDA_F(:) = ZAVG_CAN_HW_RATIO*ZAVG_BLD / (0.5*XPI)
    !
      CALL SM10(XZ,ZAVG_BLD_HEIGHT,ZLAMBDA_F,ZL)
    !
    !-------------------------------------------------------------------------------------
    ! computes coefficients for implicitation
    !-------------------------------------------------------------------------------------
    !
      ZAVG_UW_GRND(:)      = 0.
      ZAVG_DUWDU_GRND(:)   = 0.
      ZAVG_UW_ROOF(:)      = 0.
      ZAVG_DUWDU_ROOF(:)   = 0.
      ZAVG_H_GRND(:)       = 0.
      ZAVG_H_WALL(:)       = 0.
      ZAVG_H_ROOF(:)       = 0.
      ZAVG_E_GRND(:)       = 0.
      ZAVG_E_ROOF(:)       = 0.
      ZAVG_AC_GRND(:)      = 0.
      ZAVG_AC_GRND_WAT(:)  = 0.
      ZSFLUX_U(:)          = 0.
      ZSFLUX_T(:)          = 0.
      ZSFLUX_Q(:)          = 0.
    !
      DO JLAYER=1,NLVL-1
          !* Monin-Obuhkov theory not used inside the urban canopy
          ! => neutral mixing  if layer is below : (roof level +1 meter)
          WHERE (XZ(:,JLAYER)<=ZAVG_BLD_HEIGHT(:)+1.) XLMO(:,JLAYER) = XUNDEF
      ENDDO
    !
    !
    !* computes tendencies on wind and Tke due to canopy
     CALL TEB_CANOPY(KI,NLVL,XZ,XZF,XDZ,XDZF,ZAVG_BLD,ZAVG_BLD_HEIGHT,ZAVG_WALL_O_HOR,     &
                    PPA,PRHOA,XU,                                                         &
                    ZAVG_DUWDU_GRND, ZAVG_UW_ROOF, ZAVG_DUWDU_ROOF,                       &
                    ZAVG_H_WALL,ZAVG_H_ROOF,ZAVG_E_ROOF,ZAVG_AC_GRND,ZAVG_AC_GRND_WAT,    &
                    ZFORC_U,ZDFORC_UDU,ZFORC_E,ZDFORC_EDE,ZFORC_T,ZDFORC_TDT,ZFORC_Q,ZDFORC_QDQ)
    !
    !* computes coefficients for implicitation
      CALL CANOPY_EVOL(KI,NLVL,PTSTEP,1,                         &
                         ZL,ZWIND,PTA,PQA,PPA,PRHOA,             &
                         ZSFLUX_U,ZSFLUX_T,ZSFLUX_Q,             &
                         ZFORC_U,ZDFORC_UDU,ZFORC_E,ZDFORC_EDE,  &
                         ZFORC_T,ZDFORC_TDT,ZFORC_Q,ZDFORC_QDQ,  &
                         XZ,XZF,XDZ,XDZF,XU,XTKE,XT,XQ,XLMO,     &
                         XLM,XLEPS,XP,ZAVG_USTAR,                &
                         ZALFAU,ZBETAU,ZALFAT,ZBETAT,ZALFAQ,ZBETAQ)
    !
      ZPEW_A_COEF_LOWCAN = - ZALFAU / PRHOA
      ZPEW_B_COEF_LOWCAN = ZBETAU  
    !
    !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
    ELSE              ! no canopy case
    !- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - 
      DO JI=1,KI
    !* skimming flow for h/w>1 (maximum effect of direction on wind in the canyon);
    !* isolated flow for h/w<0.5 (wind is the same in large streets for all dir.)
    !* wake flow between.
    !
        ZWAKE(JI)= 1. + (2./XPI-1.) * 2. * (ZAVG_CAN_HW_RATIO(JI)-0.5)
        ZWAKE(JI)= MAX(MIN(ZWAKE(JI),1.),2./XPI)
    !
    !* Estimation of canyon wind speed from wind just above roof level
    !  (at 1.33h). Wind at 1.33h is estimated using the log law.
    !
       IF (ZAVG_BLD_HEIGHT(JI) .GT. 0.) THEN
        ZU_CANYON(JI) = ZWAKE(JI) * EXP(-ZAVG_CAN_HW_RATIO(JI)/4.) * ZWIND(JI)     &
                  * LOG( (           2.* ZAVG_BLD_HEIGHT(JI)/3.) / ZAVG_Z0_TOWN(JI))   &
                  / LOG( (PUREF(JI)+ 2.* ZAVG_BLD_HEIGHT(JI)/3.) / ZAVG_Z0_TOWN(JI))
        ZZ_LOWCAN(JI) = ZAVG_BLD_HEIGHT(JI) / 2.
       ELSE
        ZU_CANYON(JI) = ZWIND(JI)
        ZZ_LOWCAN(JI) = PZREF(JI)
       ENDIF
     END DO
    !
    !* Without SBL scheme, canyon air is assumed at mid height
      ZU_LOWCAN=ZU_CANYON
      ZT_LOWCAN=XT_CANYON
      ZQ_LOWCAN=XQ_CANYON
      ZT_CANYON=XT_CANYON
      ZQ_CANYON=XQ_CANYON
      ZUREF    =PUREF
      ZZREF    =PZREF
      ZTA      =PTA
      ZUA      =ZWIND
      ZPA      =PPA
      ZPEW_A_COEF_LOWCAN =  0.
      ZPEW_B_COEF_LOWCAN =  ZU_CANYON
    END IF
    !
    ! Exner functions
    !
    ZEXNS     (:) = (PPS(:)/XP00)**(XRD/XCPD)
    ZEXNA     (:) = (ZPA(:)/XP00)**(XRD/XCPD)
    
    !--------------------------------------------------------------------------------------
    ! Over Urban surfaces/towns:
    !--------------------------------------------------------------------------------------
    !
    DO JTEB_PATCH=1,NTEB_PATCH
     CALL GOTO_TEB(JTEB_PATCH)
    !
    ZT_CAN=ZT_CANYON
    ZQ_CAN=ZQ_CANYON
    !
    IF (LCANOPY) THEN
      XT_CANYON(:) = ZT_CANYON(:)
      XQ_CANYON(:) = ZQ_CANYON(:)
    END IF
    !
    ZLESNOW_ROOF(:) = 0.
    ZLESNOW_ROAD(:) = 0.
    ZG_GREENROOF_ROOF(:) = 0.
    !
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    !  Reinitialize shading of windows when changing day
    ! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
    !
    IF (CBEM=='BEM') &
    WHERE (PTSUN .LT. PTSTEP + 1E-3) LSHAD_DAY(:) = .FALSE.
    !
    !
     CALL TEB_GARDEN      (LGARDEN, LGREENROOF, CZ0H, CIMPLICIT_WIND, CROAD_DIR, CWALL_OPT,&
                          TTIME, PTSUN, ZT_CAN, ZQ_CAN, ZU_CANYON,                         &
                          ZT_LOWCAN, ZQ_LOWCAN, ZU_LOWCAN, ZZ_LOWCAN,                      &
                          XTI_BLD,                                                         &
                          XT_ROOF, XT_ROAD, XT_WALL_A, XT_WALL_B, XWS_ROOF,XWS_ROAD,       &
                          TSNOW_ROOF%SCHEME,                                               &
                          TSNOW_ROOF%WSNOW(:,:,1), TSNOW_ROOF%T(:,:,1),                    &
                          TSNOW_ROOF%RHO(:,:,1), TSNOW_ROOF%ALB(:,1),                      &
                          TSNOW_ROOF%TS(:,1), TSNOW_ROOF%EMIS(:,1),                        &
                          TSNOW_ROAD%SCHEME,                                               &
                          TSNOW_ROAD%WSNOW(:,:,1), TSNOW_ROAD%T(:,:,1),                    &
                          TSNOW_ROAD%RHO(:,:,1), TSNOW_ROAD%ALB(:,1),                      &
                          TSNOW_ROAD%TS(:,1), TSNOW_ROAD%EMIS(:,1),                        &
                          ZPEW_A_COEF, ZPEW_B_COEF,                                        &
                          ZPEW_A_COEF_LOWCAN, ZPEW_B_COEF_LOWCAN,                          &
                          PPS, ZPA, ZEXNS, ZEXNA, ZTA, ZQA, PRHOA, PCO2,                   &
                          PLW, ZDIR_SWB, ZSCA_SWB, PSW_BANDS, KSW, PZENITH, PAZIM,         &
                          PRAIN, PSNOW, ZZREF, ZUREF, ZUA,                                 &
                          ZH_TRAFFIC, ZLE_TRAFFIC, XH_INDUSTRY, XLE_INDUSTRY,              &
                          PTSTEP,                                                          &
                          XZ0_TOWN,                                                        &
                          XBLD, XGARDEN, XROAD_DIR, XROAD, XGREENROOF,                     &
                          XBLD_HEIGHT, XWALL_O_HOR, XCAN_HW_RATIO,                         &
                          XROAD_O_GRND, XGARDEN_O_GRND, XWALL_O_GRND,                      &
                          XALB_ROOF, XEMIS_ROOF,                                           &
                          XHC_ROOF,XTC_ROOF,XD_ROOF,                                       &
                          XALB_ROAD, XEMIS_ROAD, XSVF_ROAD,                                &
                          XHC_ROAD,XTC_ROAD,XD_ROAD,                                       &
                          XALB_WALL, XEMIS_WALL, XSVF_WALL,                                &
                          XSVF_GARDEN,                                                     &
                          XHC_WALL,XTC_WALL,XD_WALL,                                       &
                          ZRN_ROOF, ZH_ROOF, ZLE_ROOF, ZLEW_ROOF, ZGFLUX_ROOF,             &
                          ZRUNOFF_ROOF,                                                    &
                          ZRN_ROAD, ZH_ROAD, ZLE_ROAD, ZLEW_ROAD, ZGFLUX_ROAD,             &
                          ZRUNOFF_ROAD,                                                    &
                          ZRN_WALL_A, ZH_WALL_A, ZLE_WALL_A, ZGFLUX_WALL_A,                &
                          ZRN_WALL_B, ZH_WALL_B, ZLE_WALL_B, ZGFLUX_WALL_B,                &
                          ZRN_GARDEN,ZH_GARDEN,ZLE_GARDEN,ZGFLUX_GARDEN,                   &
                          ZRN_GREENROOF,ZH_GREENROOF,ZLE_GREENROOF,ZGFLUX_GREENROOF,       &
                          ZRN_STRLROOF,ZH_STRLROOF,ZLE_STRLROOF,ZGFLUX_STRLROOF,           &
                          ZRN_BLT,ZH_BLT,ZLE_BLT,ZGFLUX_BLT,                               &
                          ZRNSNOW_ROOF, ZHSNOW_ROOF, ZLESNOW_ROOF, ZGSNOW_ROOF,            &
                          ZMELT_ROOF,                                                      &
                          ZRNSNOW_ROAD, ZHSNOW_ROAD, ZLESNOW_ROAD, ZGSNOW_ROAD,            &
                          ZMELT_ROAD,                                                      &
                          ZRN_GRND, ZH_GRND, ZLE_GRND, ZGFLUX_GRND,                        &
                          ZRN, ZH, ZLE, ZGFLUX, ZEVAP, ZRUNOFF, ZSFCO2,                    &
                          ZUW_GRND, ZUW_ROOF, ZDUWDU_GRND, ZDUWDU_ROOF,                    &
                          ZUSTAR, ZCD, ZCDN, ZCH, ZRI,                                     &
                          ZTRAD, ZEMIS, ZDIR_ALB, ZSCA_ALB, ZRESA_TOWN, ZDQS_TOWN,         &
                          ZQF_TOWN, ZQF_BLD,                                               &
                          ZFLX_BLD, ZAC_ROAD, ZAC_GARDEN, ZAC_GREENROOF,                   &
                          ZAC_ROAD_WAT, ZAC_GARDEN_WAT, ZAC_GREENROOF_WAT,                 &
                          ZABS_SW_ROOF,ZABS_LW_ROOF,                                       &
                          ZABS_SW_SNOW_ROOF,ZABS_LW_SNOW_ROOF,                             &
                          ZABS_SW_ROAD,ZABS_LW_ROAD,                                       &
                          ZABS_SW_SNOW_ROAD,ZABS_LW_SNOW_ROAD,                             &
                          ZABS_SW_WALL_A, ZABS_LW_WALL_A,                                  &
                          ZABS_SW_WALL_B, ZABS_LW_WALL_B,                                  &
                          ZABS_SW_GARDEN,ZABS_LW_GARDEN,                                   &
                          ZABS_SW_GREENROOF,ZABS_LW_GREENROOF, ZG_GREENROOF_ROOF,          &
                          ZRUNOFF_GREENROOF, ZDRAIN_GREENROOF,                             &
                          CCOOL_COIL, XF_WATER_COND, CHEAT_COIL, CNATVENT,                 &
                          KDAY, XAUX_MAX, XT_FLOOR, XT_MASS, ZH_BLD_COOL,                  &
                          ZT_BLD_COOL, ZH_BLD_HEAT, ZLE_BLD_COOL, ZLE_BLD_HEAT, ZH_WASTE,  &
                          ZLE_WASTE, XF_WASTE_CAN, ZHVAC_COOL, ZHVAC_HEAT, XQIN, XQIN_FRAD,&
                          XQIN_FLAT, XGR, XEFF_HEAT, XINF, XTCOOL_TARGET,                  &
                          XTHEAT_TARGET, XHR_TARGET, XT_WIN2, XQI_BLD, XV_VENT,            &
                          XCAP_SYS_HEAT, XCAP_SYS_RAT, XT_ADP, XM_SYS_RAT, XCOP_RAT,       &
                          ZCAP_SYS, ZM_SYS, ZCOP, ZQ_SYS, ZT_SYS, ZTR_SW_WIN, ZFAN_POWER,  &
                          XHC_FLOOR, XTC_FLOOR, XD_FLOOR, XT_WIN1, ZABS_SW_WIN,            &
                          ZABS_LW_WIN, XSHGC, XSHGC_SH, XUGG_WIN, XALB_WIN, XABS_WIN,      &
                          ZEMIT_LW_FAC, ZEMIT_LW_GRND, ZT_RAD_IND, ZREF_SW_GRND,           &
                          ZREF_SW_FAC, ZHU_BLD, PTIME, LSHADE, LSHAD_DAY, LNATVENT_NIGHT,  &
                          CBEM, XN_FLOOR, XWALL_O_BLD, XGLAZ_O_BLD, XMASS_O_BLD,           &
                          XFLOOR_HW_RATIO,                                                 &
                          XF_FLOOR_MASS, XF_FLOOR_WALL, XF_FLOOR_WIN,                      &
                          XF_FLOOR_ROOF, XF_WALL_FLOOR, XF_WALL_MASS,                      &
                          XF_WALL_WIN, XF_WIN_FLOOR, XF_WIN_MASS, XF_WIN_WALL,             &
                          XF_MASS_FLOOR, XF_MASS_WALL, XF_MASS_WIN, LCANOPY, XTRAN_WIN,    &
                          CCH_BEM, XROUGH_ROOF, XROUGH_WALL, XF_WIN_WIN                    )
    
    
    !
    IF (.NOT. LCANOPY) THEN
      CALL ADD_PATCH_CONTRIB(JTEB_PATCH,ZAVG_T_CANYON,ZT_CAN)
      CALL ADD_PATCH_CONTRIB(JTEB_PATCH,ZAVG_Q_CANYON,ZQ_CAN)
    !
    ! Momentum fluxes
    !
      ZSFU = 0.
      ZSFV = 0.
      DO JJ=1,SIZE(PU)
        IF (ZWIND(JJ)>0.) THEN
          ZCOEF(JJ) = - PRHOA(JJ) * ZUSTAR(JJ)**2 / ZWIND(JJ)
          ZSFU(JJ) = ZCOEF(JJ) * PU(JJ)
          ZSFV(JJ) = ZCOEF(JJ) * PV(JJ)
        ENDIF
      ENDDO
      CALL ADD_PATCH_CONTRIB(JTEB_PATCH,PSFU,ZSFU)
      CALL ADD_PATCH_CONTRIB(JTEB_PATCH,PSFV,ZSFV)
    !
    ENDIF
    !
    !-------------------------------------------------------------------------------------
    ! Outputs:
    !-------------------------------------------------------------------------------------
    !
    ! Grid box average fluxes/properties: Arguments and standard diagnostics
    !
     CALL ADD_PATCH_CONTRIB(JTEB_PATCH,PSFTH,ZH)
     CALL ADD_PATCH_CONTRIB(JTEB_PATCH,PSFTQ,ZEVAP)
     CALL ADD_PATCH_CONTRIB(JTEB_PATCH,PSFCO2,ZSFCO2)
    !
    !
    ! Albedo for each wavelength and patch
    !
    DO JSWB=1,SIZE(PSW_BANDS)
      DO JJ=1,SIZE(ZDIR_ALB)
        ZDIR_ALB_PATCH(JJ,JSWB,JTEB_PATCH) = ZDIR_ALB(JJ)
        ZSCA_ALB_PATCH(JJ,JSWB,JTEB_PATCH) = ZSCA_ALB(JJ)
      ENDDO
    END DO
    !
    ! emissivity and radiative temperature
    !
    ZEMIS_PATCH(:,JTEB_PATCH) = ZEMIS
    ZTRAD_PATCH(:,JTEB_PATCH) = ZTRAD
    !
    ! computes some aggregated diagnostics
    !
     IF (.NOT. LCANOPY) THEN
       CALL ADD_PATCH_CONTRIB(JTEB_PATCH,ZAVG_CD ,ZCD )
       CALL ADD_PATCH_CONTRIB(JTEB_PATCH,ZAVG_CDN,ZCDN)
       CALL ADD_PATCH_CONTRIB(JTEB_PATCH,ZAVG_RI ,ZRI )
     ENDIF
     CALL ADD_PATCH_CONTRIB(JTEB_PATCH,ZAVG_CH ,ZCH )
     CALL ADD_PATCH_CONTRIB(JTEB_PATCH,ZAVG_RN ,ZRN )
     CALL ADD_PATCH_CONTRIB(JTEB_PATCH,ZAVG_H  ,ZH  )
     CALL ADD_PATCH_CONTRIB(JTEB_PATCH,ZAVG_LE ,ZLE )
     CALL ADD_PATCH_CONTRIB(JTEB_PATCH,ZAVG_GFLUX ,ZGFLUX )
    !
    !* warning: aerodynamical resistance does not yet take into account gardens
     CALL ADD_PATCH_CONTRIB(JTEB_PATCH,ZAVG_RESA_TOWN,1./ZRESA_TOWN)
    IF (JTEB_PATCH==NTEB_PATCH) ZAVG_RESA_TOWN = 1./ZAVG_RESA_TOWN
    !
    !-------------------------------------------------------------------------------------
    ! Diagnostics on each patch
    !-------------------------------------------------------------------------------------
    !
     CALL DIAG_MISC_TEB_n(PTSTEP, ZDQS_TOWN, ZQF_BLD, ZQF_TOWN, ZFLX_BLD,           &
                         ZRN_ROAD, ZH_ROAD, ZLE_ROAD, ZGFLUX_ROAD,                 &
                         ZRN_WALL_A, ZH_WALL_A, ZGFLUX_WALL_A,                     &
                         ZRN_WALL_B, ZH_WALL_B, ZGFLUX_WALL_B,                     &
                         ZRN_ROOF, ZH_ROOF, ZLE_ROOF, ZGFLUX_ROOF, ZRUNOFF,        &
                         ZRN_STRLROOF, ZH_STRLROOF, ZLE_STRLROOF, ZGFLUX_STRLROOF, &
                         ZRN_GREENROOF, ZH_GREENROOF,                              &
                         ZLE_GREENROOF, ZGFLUX_GREENROOF, ZG_GREENROOF_ROOF,       &
                         ZRUNOFF_GREENROOF, ZDRAIN_GREENROOF,                      &
                         ZRN_GARDEN,ZH_GARDEN,ZLE_GARDEN,ZGFLUX_GARDEN,            &
                         ZRN_BLT,ZH_BLT,ZLE_BLT,ZGFLUX_BLT,                        &
                         ZABS_SW_ROOF,ZABS_LW_ROOF,                                &
                         ZABS_SW_SNOW_ROOF,ZABS_LW_SNOW_ROOF,                      &
                         ZABS_SW_ROAD,ZABS_LW_ROAD,                                &
                         ZABS_SW_SNOW_ROAD,ZABS_LW_SNOW_ROAD,                      &
                         ZABS_SW_WALL_A, ZABS_LW_WALL_A, ZABS_SW_WALL_B,           &
                         ZABS_LW_WALL_B,                                           &
                         ZABS_SW_GARDEN,ZABS_LW_GARDEN,                            &
                         ZABS_SW_GREENROOF,ZABS_LW_GREENROOF,                      &
                         ZH_BLD_COOL, ZT_BLD_COOL,                                 &     
                         ZH_BLD_HEAT, ZLE_BLD_COOL, ZLE_BLD_HEAT,                  &
                         ZH_WASTE, ZLE_WASTE, ZHVAC_COOL,                          &
                         ZHVAC_HEAT, ZCAP_SYS, ZM_SYS, ZCOP,                       &
                         ZQ_SYS, ZT_SYS, ZTR_SW_WIN, ZFAN_POWER,                   &
                         ZABS_SW_WIN, ZABS_LW_WIN                                  )
    !
    !
    !-------------------------------------------------------------------------------------
    ! Computes averaged parameters necessary for UTCI
    !-------------------------------------------------------------------------------------
    !
    IF (N2M >0 .AND. LUTCI) THEN
      CALL ADD_PATCH_CONTRIB(JTEB_PATCH,ZAVG_REF_SW_GRND ,ZREF_SW_GRND )
      CALL ADD_PATCH_CONTRIB(JTEB_PATCH,ZAVG_REF_SW_FAC  ,ZREF_SW_FAC  )
      CALL ADD_PATCH_CONTRIB(JTEB_PATCH,ZAVG_SCA_SW      ,ZSCA_SW      )
      CALL ADD_PATCH_CONTRIB(JTEB_PATCH,ZAVG_DIR_SW      ,ZDIR_SW      )
      CALL ADD_PATCH_CONTRIB(JTEB_PATCH,ZAVG_EMIT_LW_FAC ,ZEMIT_LW_FAC )
      CALL ADD_PATCH_CONTRIB(JTEB_PATCH,ZAVG_EMIT_LW_GRND,ZEMIT_LW_GRND)
      CALL ADD_PATCH_CONTRIB(JTEB_PATCH,ZAVG_T_RAD_IND   ,ZT_RAD_IND   )
      CALL ADD_PATCH_CONTRIB(JTEB_PATCH,ZAVG_TI_BLD      ,XTI_BLD      )
      CALL ADD_PATCH_CONTRIB(JTEB_PATCH,ZAVG_QI_BLD      ,XQI_BLD      )
    END IF
    !
    !-------------------------------------------------------------------------------------
    ! Use of the canopy version of TEB
    !-------------------------------------------------------------------------------------
    !
    IF (LCANOPY) THEN
    !-------------------------------------------------------------------------------------
    ! Town averaged quantities to force canopy atmospheric layers
    !-------------------------------------------------------------------------------------
    
     CALL ADD_PATCH_CONTRIB(JTEB_PATCH,ZAVG_DUWDU_GRND ,ZDUWDU_GRND )
     CALL ADD_PATCH_CONTRIB(JTEB_PATCH,ZAVG_UW_ROOF ,ZUW_ROOF)
     CALL ADD_PATCH_CONTRIB(JTEB_PATCH,ZAVG_DUWDU_ROOF ,ZDUWDU_ROOF)
     CALL ADD_PATCH_CONTRIB(JTEB_PATCH,ZAVG_H_WALL ,0.5*(ZH_WALL_A+ZH_WALL_B))
     CALL ADD_PATCH_CONTRIB(JTEB_PATCH,ZAVG_H_ROOF ,(ZH_ROOF+XH_INDUSTRY))
     CALL ADD_PATCH_CONTRIB(JTEB_PATCH,ZAVG_E_ROOF ,(ZLE_ROOF+XLE_INDUSTRY)/XLVTT)
    !
    !-------------------------------------------------------------------------------------
    ! Computes the impact of canopy and surfaces on air
    !-------------------------------------------------------------------------------------
    !
    ZAC_GRND    (:) = (XROAD(:)*ZAC_ROAD    (:) + XGARDEN(:)*ZAC_GARDEN    (:)) / (XROAD(:)+XGARDEN(:))
    ZAC_GRND_WAT(:) = (XROAD(:)*ZAC_ROAD_WAT(:) + XGARDEN(:)*ZAC_GARDEN_WAT(:)) / (XROAD(:)+XGARDEN(:))
    !
     CALL ADD_PATCH_CONTRIB(JTEB_PATCH,ZAVG_AC_GRND     ,ZAC_GRND    )
     CALL ADD_PATCH_CONTRIB(JTEB_PATCH,ZAVG_AC_GRND_WAT ,ZAC_GRND_WAT)
     CALL ADD_PATCH_CONTRIB(JTEB_PATCH,ZSFLUX_U ,ZUW_GRND * (1.-XBLD))
     CALL ADD_PATCH_CONTRIB(JTEB_PATCH,ZSFLUX_T ,ZH_GRND  * (1.-XBLD)/XCPD/PRHOA)
     CALL ADD_PATCH_CONTRIB(JTEB_PATCH,ZSFLUX_Q ,ZLE_GRND * (1.-XBLD)/XLVTT)
    !
    
    END IF
    !
    !-------------------------------------------------------------------------------------
    ! end of loop on TEB patches
    END DO
    !-------------------------------------------------------------------------------------
    !
    !-------------------------------------------------------------------------------------
    !* Evolution of canopy air if canopy option is active
    !-------------------------------------------------------------------------------------
    !
    IF (LCANOPY) THEN
    !
    !-------------------------------------------------------------------------------------
    !* Impact of TEB fluxes on the air
    !-------------------------------------------------------------------------------------
    !
     CALL TEB_CANOPY(KI,NLVL,XZ,XZF,XDZ,XDZF,ZAVG_BLD,ZAVG_BLD_HEIGHT,ZAVG_WALL_O_HOR,     &
                    PPA,PRHOA,XU,                                                         &
                    ZAVG_DUWDU_GRND, ZAVG_UW_ROOF, ZAVG_DUWDU_ROOF,                       &
                    ZAVG_H_WALL,ZAVG_H_ROOF,ZAVG_E_ROOF,ZAVG_AC_GRND,ZAVG_AC_GRND_WAT,    &
                    ZFORC_U,ZDFORC_UDU,ZFORC_E,ZDFORC_EDE,ZFORC_T,ZDFORC_TDT,ZFORC_Q,ZDFORC_QDQ)
    !
    !-------------------------------------------------------------------------------------
    !* Evolution of canopy air due to these impacts
    !-------------------------------------------------------------------------------------
    !
     CALL CANOPY_EVOL(KI,NLVL,PTSTEP,2,                                            &
                     ZL,ZWIND,PTA,PQA,PPA,PRHOA,                                  &
                     ZSFLUX_U,ZSFLUX_T,ZSFLUX_Q,                                  &
                     ZFORC_U,ZDFORC_UDU,ZFORC_E,ZDFORC_EDE,                       &