urban_drag.F90

Go to the documentation of this file.

!SFX_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
!SFX_LIC This is part of the SURFEX software governed by the CeCILL-C licence
!SFX_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
!SFX_LIC for details. version 1.
!     #########
    SUBROUTINE urban_drag(HZ0H, HIMPLICIT_WIND, PTSTEP,                     &
                          pt_canyon, pq_canyon, pu_canyon,                  &
                          pt_lowcan, pq_lowcan, pu_lowcan, pz_lowcan,       &
                          pts_roof, pts_road, pts_wall, pts_garden,         &
                          pdelt_snow_roof, pdelt_snow_road,                 &
                          pexns, pexna, pta, pqa, pps, prhoa,               &
                          pzref, puref, pvmod,                              &
                          pz0_town,                                         &
                          pbld, pgarden, proad,                             &
                          pbld_height, pcan_hw_ratio,                       &
                          pwall_o_grnd,                                     &
                          pws_roof, pws_road,                               &
                          pws_roof_max, pws_road_max,                       &
                          ppew_a_coef, ppew_b_coef,                         &
                          ppew_a_coef_lowcan, ppew_b_coef_lowcan,           &
                          pqsat_roof, pqsat_road, pdelt_roof, pdelt_road,   &
                          pcd, pcdn, pac_roof, pac_roof_wat,                &
                          pac_wall, pac_road, pac_road_wat, pac_top,        &
                          pac_garden, pri,                                  &
                          puw_road, puw_roof, pduwdu_road, pduwdu_roof,     &
                          pustar_town, ocanopy, pts_win, pac_win, hch_bem,  &
                          prough_roof, prough_wall                          ) 
!   ##########################################################################
!
!!****  *URBAN_DRAG*  
!!
!!    PURPOSE
!!    -------
!
!     Computes the surface drag over artificial surfaces as towns, 
!     taking into account the canyon like geometry of urbanized areas.
!         
!     
!!**  METHOD
!!    ------
!
!
!
!
!!    EXTERNAL
!!    --------
!!
!!
!!    IMPLICIT ARGUMENTS
!!    ------------------
!!
!!      
!!    REFERENCE
!!    ---------
!!
!!      
!!    AUTHOR
!!    ------
!!
!!      V. Masson           * Meteo-France *
!!
!!    MODIFICATIONS
!!    -------------
!!      Original    20/01/98 
!!          01/00 (V. Masson)  separation of skimming, wake and isolated flows
!!          09/00 (V. Masson)  use of Z0 for roads
!!          12/02 (A. Lemonsu) convective speed w* in canyon
!             04 (A. Lemonsu) z0h=z0m for resistance canyon-atmosphere
!          03/08 (S. Leroyer) debug PU_CAN (1. * H/3)
!          12/08 (S. Leroyer) option (HZ0H) for z0h applied on roof, road and town
!!         09/12 B. Decharme new wind implicitation
!          11/11 (G. Pigeon) apply only urban_exch_coef when necessary if
!                            canopy/no canopy
!          09/12 (G. Pigeon) add new formulation for outdoor conv. coef for
!                            wall/roof/window
!!
!-------------------------------------------------------------------------------
!
!*       0.     DECLARATIONS
!               ------------
!
USE modd_surf_par, ONLY : xundef
USE modd_csts,ONLY : xlvtt, xpi, xcpd, xg, xkarman
!
!USE MODE_SBLS
USE mode_thermos
USE modi_urban_exch_coef
USE mode_conv_doe
!
USE yomhook   ,ONLY : lhook,   dr_hook
USE parkind1  ,ONLY : jprb
!
IMPLICIT NONE
!
!*      0.1    declarations of arguments
!
!
 CHARACTER(LEN=6),   INTENT(IN)    :: hz0h           ! TEB option for z0h roof & road
!
 CHARACTER(LEN=*),     INTENT(IN)  :: himplicit_wind   ! wind implicitation option
!                                                     ! 'OLD' = direct
!                                                     ! 'NEW' = Taylor serie, order 1
!
REAL,               INTENT(IN)    :: ptstep         ! time-step
REAL, DIMENSION(:), INTENT(IN)    :: pt_canyon      ! canyon air temperature
REAL, DIMENSION(:), INTENT(IN)    :: pq_canyon      ! canyon air specific humidity.
REAL, DIMENSION(:), INTENT(IN)    :: pu_canyon      ! hor. wind in canyon
REAL, DIMENSION(:), INTENT(IN)    :: pu_lowcan     ! wind near the road
REAL, DIMENSION(:), INTENT(IN)    :: pt_lowcan     ! temp. near the road
REAL, DIMENSION(:), INTENT(IN)    :: pq_lowcan     ! hum. near the road
REAL, DIMENSION(:), INTENT(IN)    :: pz_lowcan     ! height of atm. var. near the road
REAL, DIMENSION(:), INTENT(IN)    :: pts_roof       ! surface temperature
REAL, DIMENSION(:), INTENT(IN)    :: pts_road       ! surface temperature
REAL, DIMENSION(:), INTENT(IN)    :: pts_wall       ! surface temperature
REAL, DIMENSION(:), INTENT(IN)    :: pts_garden     ! surface temperature
REAL, DIMENSION(:), INTENT(IN)    :: pdelt_snow_roof! fraction of snow on roof
REAL, DIMENSION(:), INTENT(IN)    :: pdelt_snow_road! fraction of snow on road
REAL, DIMENSION(:), INTENT(IN)    :: pexns          ! surface exner function
REAL, DIMENSION(:), INTENT(IN)    :: pta            ! temperature at the lowest level
REAL, DIMENSION(:), INTENT(IN)    :: pqa            ! specific humidity
                                                    ! at the lowest level
REAL, DIMENSION(:), INTENT(IN)    :: pvmod          ! module of the horizontal wind
REAL, DIMENSION(:), INTENT(IN)    :: pps            ! pressure at the surface
REAL, DIMENSION(:), INTENT(IN)    :: pexna          ! exner function
                                                    ! at the lowest level
REAL, DIMENSION(:), INTENT(IN)    :: prhoa          ! air density
REAL, DIMENSION(:), INTENT(IN)    :: pzref          ! reference height of the first
                                                    ! atmospheric level (temperature)
REAL, DIMENSION(:), INTENT(IN)    :: puref          ! reference height of the first
                                                    ! atmospheric level (wind)
REAL, DIMENSION(:), INTENT(IN)    :: pz0_town       ! roughness length for momentum
REAL, DIMENSION(:), INTENT(IN)    :: pbld           ! fraction of buildings
REAL, DIMENSION(:), INTENT(IN)    :: proad          ! fraction of roads
REAL, DIMENSION(:), INTENT(IN)    :: pgarden        ! fraction of GARDEN areas
REAL, DIMENSION(:), INTENT(IN)    :: pbld_height    ! h
REAL, DIMENSION(:), INTENT(IN)    :: pcan_hw_ratio  ! h/W
REAL, DIMENSION(:), INTENT(IN)    :: pwall_o_grnd   ! wall surf. / (road+GARDEN area) surf.
!
REAL, DIMENSION(:), INTENT(IN)    :: pws_roof       ! roof water content (kg/m2)
REAL, DIMENSION(:), INTENT(IN)    :: pws_road       ! road water content (kg/m2)
REAL, DIMENSION(:), INTENT(IN)    :: pws_roof_max   ! maximum deepness of roof
REAL, DIMENSION(:), INTENT(IN)    :: pws_road_max   ! and water reservoirs (kg/m2)
!
!
REAL, DIMENSION(:), INTENT(IN)    :: ppew_a_coef    ! implicit coefficients (m2s/kg)
REAL, DIMENSION(:), INTENT(IN)    :: ppew_b_coef    ! for wind coupling     (m/s)
REAL, DIMENSION(:), INTENT(IN)    :: ppew_a_coef_lowcan ! implicit coefficients for wind coupling (m2s/kg)
REAL, DIMENSION(:), INTENT(IN)    :: ppew_b_coef_lowcan ! between low canyon wind and road (m/s)
!
REAL, DIMENSION(:), INTENT(OUT)   :: pqsat_roof     ! qsat(Ts)
REAL, DIMENSION(:), INTENT(OUT)   :: pqsat_road     ! qsat(Ts)
REAL, DIMENSION(:), INTENT(OUT)   :: pdelt_roof     ! water fraction on
REAL, DIMENSION(:), INTENT(OUT)   :: pdelt_road     ! snow-free surfaces
REAL, DIMENSION(:), INTENT(OUT)   :: pcd            ! drag coefficient
REAL, DIMENSION(:), INTENT(OUT)   :: pcdn           ! neutral drag coefficient
REAL, DIMENSION(:), INTENT(OUT)   :: pac_roof       ! aerodynamical conductance
REAL, DIMENSION(:), INTENT(OUT)   :: pac_roof_wat   ! aerodynamical conductance (for water)
REAL, DIMENSION(:), INTENT(OUT)   :: pac_wall       ! aerodynamical conductance
!                                                   ! between canyon air and
!                                                   ! walls 
REAL, DIMENSION(:), INTENT(OUT)   :: pac_road       ! aerodynamical conductance
!                                                   ! between canyon air and
!                                                   ! roads
REAL, DIMENSION(:), INTENT(OUT)   :: pac_road_wat   ! aerodynamical conductance
!                                                   ! between canyon air and
!                                                   ! road (for water)
REAL, DIMENSION(:), INTENT(OUT)   :: pac_top        ! aerodynamical conductance
!                                                   ! between canyon top and atm.
REAL, DIMENSION(:), INTENT(IN)    :: pac_garden     ! aerodynamical conductance
!                                                   ! between canyon air and GARDEN areas
REAL, DIMENSION(:), INTENT(OUT)   :: pri            ! Town Richardson number
!
REAL, DIMENSION(:), INTENT(OUT)   :: puw_road       ! Momentum flux for roads
REAL, DIMENSION(:), INTENT(OUT)   :: puw_roof       ! Momentum flux for roofs
REAL, DIMENSION(:), INTENT(OUT)   :: pduwdu_road    ! 
REAL, DIMENSION(:), INTENT(OUT)   :: pduwdu_roof    ! 
REAL, DIMENSION(:), INTENT(OUT)   :: pustar_town    ! Fraction velocity for town
!
LOGICAL, INTENT(IN)               :: ocanopy        ! is canopy active
REAL, DIMENSION(:), INTENT(OUT)   :: pac_win        ! aerodynamical conductance for window
REAL, DIMENSION(:), INTENT(IN)    :: pts_win        ! window outdoor surface temp
 CHARACTER(LEN=5), INTENT(IN)      :: hch_bem        ! BEM option for convective heat transfer coef.
REAL, DIMENSION(:), INTENT(IN)    :: prough_roof    ! roughness coef for the roof
REAL, DIMENSION(:), INTENT(IN)    :: prough_wall    ! roughness coef for the wall
!
!*      0.2    declarations of local variables
!
!
REAL, DIMENSION(SIZE(PTA)) :: zts_town     ! town averaged temp.
REAL, DIMENSION(SIZE(PTA)) :: zq_town      ! town averaged hum.
REAL, DIMENSION(SIZE(PTA)) :: zavdelt_roof ! averaged water frac.
REAL, DIMENSION(SIZE(PTA)) :: zq_roof      ! roof spec. hum.
REAL, DIMENSION(SIZE(PTA)) :: zz0_roof     ! roof roughness length
REAL, DIMENSION(SIZE(PTA)) :: zz0_road     ! road roughness length
REAL, DIMENSION(SIZE(PTA)) :: zw_can       ! ver. wind in canyon
REAL, DIMENSION(SIZE(PTA)) :: zri          ! Richardson number
REAL, DIMENSION(SIZE(PTA)) :: zle_max      ! maximum latent heat flux available
REAL, DIMENSION(SIZE(PTA)) :: zle          ! actual latent heat flux
REAL, DIMENSION(SIZE(PTA)) :: zra_roof     ! aerodynamical resistance
REAL, DIMENSION(SIZE(PTA)) :: zch_roof     ! drag coefficient for heat
REAL, DIMENSION(SIZE(PTA)) :: zra_top      ! aerodynamical resistance
REAL, DIMENSION(SIZE(PTA)) :: zch_top      ! drag coefficient for heat
REAL, DIMENSION(SIZE(PTA)) :: zra_road     ! aerodynamical resistance
REAL, DIMENSION(SIZE(PTA)) :: zch_road     ! drag coeifficient for heat
REAL, DIMENSION(SIZE(PTA)) :: zcd_road     ! road  surf. exchange coefficient
REAL, DIMENSION(SIZE(PTA)) :: zac          ! town aerodynamical conductance (not used)
REAL, DIMENSION(SIZE(PTA)) :: zra          ! town aerodynamical resistance  (not used)
REAL, DIMENSION(SIZE(PTA)) :: zch          ! town drag coefficient for heat (not used)
REAL, DIMENSION(SIZE(PTA)) :: zcd          ! any surf. exchange coefficient (not used)
REAL, DIMENSION(SIZE(PTA)) :: zcdn         ! any surf. neutral exch. coef.  (not used)
!
REAL, DIMENSION(SIZE(PTA)) :: zu_star, zw_star !! 
REAL, DIMENSION(SIZE(PTA)) :: zq0              !! 
!
REAL, DIMENSION(SIZE(PTA)) :: zustar2      ! square of friction velocity (m2/s2
REAL, DIMENSION(SIZE(PTA)) :: zvmod        ! module of the horizontal wind at t+1
!
! for calculation of momentum fluxes
REAL, DIMENSION(SIZE(PTA)) :: zlmo         ! Monin-Obukhov length
REAL, DIMENSION(SIZE(PTA)) :: zustar_road  ! friction velocity for roads
REAL, DIMENSION(SIZE(PTA)) :: zustar_roof  ! friction velocity for roofs
REAL, DIMENSION(SIZE(PTA)) :: zustar_town  !
!
REAL, DIMENSION(SIZE(PTA)) :: zz0_top      ! roughness length for zac_top calculation
REAL, DIMENSION(SIZE(PTA)) :: zchtcn_win   ! natural convective heat transfer coef. for window [W/(m2.K)]
REAL, DIMENSION(SIZE(PTA)) :: zchtcn_roof  ! natural convective heat transfer coef. for roof [W/(m2.K)]
REAL, DIMENSION(SIZE(PTA)) :: zchtcs_roof  ! forced convective heat transfer coef. for smooth roof [W/(m2.K)]
REAL, DIMENSION(SIZE(PTA)) :: zchtcn_wall  ! natural convective heat transfer coef. for wall [W/(m2.K)]
REAL, DIMENSION(SIZE(PTA)) :: zchtcs_wall  ! forced natural convective heat transfer coef. for smooth wall [W/(m2.K)]
!
INTEGER                   ::  jloop, jj            !! 
!
REAL :: zz0_o_z0h = 200.  ! z0/z0h ratio used in Mascart (1995) formulation.
!                         ! It is set to the maximum value acceptable by
!                         ! formulation. Observed values are often larger in cities.
REAL(KIND=JPRB) :: zhook_handle
!-------------------------------------------------------------------------------
!
!
IF (lhook) CALL dr_hook('URBAN_DRAG',0,zhook_handle)
!
zz0_roof(:)    = 0.15                      ! z0 for roofs
zz0_road(:)    = min(0.05,0.1*pz_lowcan(:))! z0 for roads
!
zz0_top(:) = pz0_town(:)
!
pcd(:) = xundef
pcdn(:) = xundef
pac_top(:) = xundef
pri(:) = xundef
!
puw_road(:) = xundef
puw_roof(:) = xundef
pduwdu_road(:) = xundef
pduwdu_roof(:) = xundef
pustar_town(:) = xundef
!
!-------------------------------------------------------------------------------
!
!*      1.     roof and road saturation specific humidity
!              ------------------------------------------
!
pqsat_roof(:) = qsat(pts_roof(:), pps(:))
!
pqsat_road(:) = qsat(pts_road(:), pps(:))
!
!-------------------------------------------------------------------------------
!
!*      2.     fraction of water on roofs
!              --------------------------
!
pdelt_roof=1.
!
!*      2.1    general case
!              ------------
!
WHERE (pqsat_roof(:) >= pqa(:) )
  pdelt_roof(:) = (pws_roof(:)/pws_roof_max)**(2./3.)
END WHERE
!
!*      2.2    dew deposition on roofs (PDELT_ROOF=1)
!              -----------------------
!
!-------------------------------------------------------------------------------
!
!*      3.     fraction of water on roads
!              --------------------------
!
pdelt_road=1.
!
!*      3.1    general case
!              ------------
!
WHERE (pqsat_road(:) >= pq_canyon(:) )
  pdelt_road(:) = (pws_road(:)/pws_road_max)**(2./3.)
END WHERE
!
!*      3.2    dew deposition on roads (PDELT_ROAD=1)
!              -----------------------
!
!-------------------------------------------------------------------------------
!
!*      4.     Drag coefficient for momentum between roof level and atmosphere
!              ---------------------------------------------------------------
!
!
DO jj=1,SIZE(pta)
!
!*      4.1    Averaged temperature at roof level
!              ----------------------------------
!
  zts_town(jj) = pbld(jj) * pts_roof(jj) + (1.-pbld(jj)) * pt_canyon(jj)
!
!*      4.2    Averaged water fraction on roofs
!              -------------------------------
!
  zavdelt_roof(jj) = pdelt_roof(jj) * pdelt_snow_roof(jj)
!
!*      4.3    Roof specific humidity
!              ----------------------
!
  zq_roof(jj) = pqsat_roof(jj) * zavdelt_roof(jj)
!
!*      4.4    Averaged Saturation specific humidity
!              -------------------------------------
!
  zq_town(jj) =       pbld(jj) * zq_roof(jj) + (1.-pbld(jj)) * pq_canyon(jj)
!
ENDDO
!
!-------------------------------------------------------------------------------
!
!*      5.     Momentum drag coefficient
!              -------------------------
!
IF (.NOT. ocanopy) THEN
  CALL urban_exch_coef(hz0h, zz0_o_z0h, zts_town, zq_town, pexns, pexna, pta, pqa, &
                        pzref+ pbld_height/3.,puref+pbld_height/3.,pvmod,pz0_town, &
                        pri, pcd, pcdn, zac, zra, zch                              )
ENDIF
!
!-------------------------------------------------------------------------------
!
!*      6.     Drag coefficient for heat fluxes between roofs and atmosphere
!              -------------------------------------------------------------
!
IF (hch_bem == "DOE-2") THEN
   zchtcn_roof = chtc_up_doe(pts_roof, pta)
   zchtcs_roof = chtc_smooth_wind_doe(zchtcn_roof, pvmod)
   pac_roof = chtc_rough_doe(zchtcn_roof, zchtcs_roof, prough_roof) / prhoa / xcpd
ELSE
   CALL urban_exch_coef(hz0h, zz0_o_z0h, pts_roof, zq_roof, pexns, pexna, pta, pqa, &
                        pzref, puref, pvmod, zz0_roof,                              &
                        zri, zcd, zcdn, pac_roof, zra_roof, zch_roof               )
ENDIF
!
!
DO jj=1,SIZE(pta)
  zle_max(jj)     = pws_roof(jj) / ptstep * xlvtt
  zle(jj)     =(pqsat_roof(jj) - pqa(jj))                     &
                   * pac_roof(jj) * pdelt_roof(jj) * xlvtt * prhoa(jj)
!
  pac_roof_wat(jj) = pac_roof(jj)
!
  IF (pdelt_roof(jj)==0.) pac_roof_wat(jj)=0.
!
  IF (zle(jj)>0.) pac_roof_wat(jj) = pac_roof(jj) * min( 1. , zle_max(jj)/zle(jj) )
!
ENDDO
!-------------------------------------------------------------------------------
!
!*      7.     Drag coefficient for heat fluxes between canyon and atmosphere
!              --------------------------------------------------------------
!
!* Because air/air exchanges are considered, roughness length for heat is set
!  equal to roughness length for momentum.
!
IF (.NOT. ocanopy) THEN
  CALL urban_exch_coef('MASC95', 1., pt_canyon, pq_canyon, pexns, pexna, pta, pqa, &
                        pzref+pbld_height-pz_lowcan,puref+pbld_height-pz_lowcan,   &
                        pvmod, zz0_top,                                            &
                        zri, zcd, zcdn, pac_top, zra_top, zch_top                  )
ENDIF
!
!-------------------------------------------------------------------------------
!
!*      8.     Drag coefficient for heat fluxes between walls, road and canyon
!              ---------------------------------------------------------------
!
!*      8.1    aerodynamical conductance for walls
!              -----------------------------------
!
IF (hch_bem == "DOE-2") THEN
  DO jj=1,SIZE(pta)
    zchtcn_wall(jj) = chtc_vert_doe(pts_wall(jj), pt_canyon(jj))
    zchtcs_wall(jj) = 0.5 * (chtc_smooth_lee_doe(zchtcn_wall(jj), pu_canyon(jj)) + &
                             chtc_smooth_wind_doe(zchtcn_wall(jj), pu_canyon(jj)) )
                      
    pac_wall(jj) = chtc_rough_doe(zchtcn_wall(jj), zchtcs_wall(jj), prough_wall(jj)) / xcpd / prhoa(jj)
  END DO
ELSE
  pac_wall(:) = ( 11.8 + 4.2 * pu_canyon(:) ) / xcpd / prhoa(:)
END IF
!
!*      8.2    aerodynamical conductance for roads
!              -----------------------------------
!
zw_star(:) = 0.
zq0(:)     = 0.
!
!
DO jloop=1,3
 !
  zw_can(:)   = zw_star(:)
  !
  !
  CALL urban_exch_coef(hz0h, zz0_o_z0h, pts_road, pq_lowcan, pexns, pexna,  &
                        pt_lowcan, pq_lowcan,                               &
                        pz_lowcan, pz_lowcan, pu_canyon+zw_can, zz0_road,   &
                        zri, zcd_road, zcdn, pac_road, zra_road, zch_road   )
  !
  DO jj=1,SIZE(pta)

    zq0(jj)     = (pts_wall(jj) - pt_canyon(jj)) * pac_wall(jj) * pwall_o_grnd(jj)

    IF (proad(jj)+pgarden(jj) .GT. 0.) THEN
      zq0(jj)   = zq0(jj) &
                +(pts_road(jj) - pt_lowcan(jj)) * pac_road(jj) * proad(jj)/(proad(jj)+pgarden(jj)) &
                +(pts_garden(jj) - pt_lowcan(jj)) * pac_garden(jj) * pgarden(jj)/(proad(jj)+pgarden(jj))
    ENDIF
    !
    IF (zq0(jj) >= 0.) THEN
      zw_star(jj) = ( (xg * pexna(jj) / pta(jj)) * zq0(jj) * pbld_height(jj)) ** (1/3.)
    ELSE
      zw_star(jj) = 0.
    ENDIF
!
  ENDDO
!
END DO
!
!
!*      8.4    aerodynamical conductance for water limited by available water
!              --------------------------------------------------------------
!
DO jj=1,SIZE(pta)
  !
  zle_max(jj)     = pws_road(jj) / ptstep * xlvtt
  zle(jj)     = ( pqsat_road(jj) - pq_lowcan(jj) )                   &
                   *   pac_road(jj) * pdelt_road(jj) * xlvtt * prhoa(jj)
  !
  pac_road_wat(jj) = pac_road(jj)
  !
  IF (pdelt_road(jj)==0.) pac_road_wat(jj) = 0.
  !
  IF (zle(jj)>0.) pac_road_wat(jj) = pac_road(jj) * min( 1. , zle_max(jj)/zle(jj) )
  !
  !
  !*      8.5    aerodynamical conductance for window
  !              ------------------------------------
  !
  zchtcn_win(jj) = chtc_vert_doe(pts_win(jj), pt_canyon(jj))
  !
  pac_win(jj) = 0.5 * (chtc_smooth_lee_doe(zchtcn_win(jj), pu_canyon(jj)) + &
                   chtc_smooth_wind_doe(zchtcn_win(jj), pu_canyon(jj)) ) &
                   / prhoa(jj) / xcpd
  !
  !-------------------------------------------------------------------------------
  !
  !*      9.     Momentum fluxes
  !              ---------------
  !
  !*      9.1    For roads
  !              ---------
  !
  !* road friction
  !
  IF (ocanopy) THEN
    !
    zustar2(jj)=xundef
    !
    IF(himplicit_wind=='OLD')THEN
      !   old implicitation
      zustar2(jj) = (zcd_road(jj)*pu_lowcan(jj)*ppew_b_coef_lowcan(jj))/              &
                    (1.0-prhoa(jj)*zcd_road(jj)*pu_lowcan(jj)*ppew_a_coef_lowcan(jj))
    ELSE
      !   new implicitation
      zustar2(jj) = (zcd_road(jj)*pu_lowcan(jj)*(2.*ppew_b_coef_lowcan(jj)-pu_lowcan(jj)))/  &
                    (1.0-2.0*prhoa(jj)*zcd_road(jj)*pu_lowcan(jj)*ppew_a_coef_lowcan(jj))
      !                   
      zvmod(jj) = prhoa(jj)*ppew_a_coef_lowcan(jj)*zustar2(jj) + ppew_b_coef_lowcan(jj)
      zvmod(jj) = max(zvmod(jj),0.)
      !
      IF(ppew_a_coef_lowcan(jj)/= 0.)THEN
        zustar2(jj) = max( ( zvmod(jj) - ppew_b_coef_lowcan(jj) ) / (prhoa(jj)*ppew_a_coef_lowcan(jj)), 0.)
      ENDIF
      !              
    ENDIF
    !
    puw_road(jj) = - zustar2(jj)
    !
    pduwdu_road(jj) = 0. ! implicitation already taken into account in PUW_ROAD
    !
    !*      9.2    For roofs
    !              ---------
    !
    !* roof friction
    !* neutral case, as guess
    !
    !
    zustar_roof(jj) = pvmod(jj) * xkarman / log(pzref(jj)/zz0_roof(jj))
    !
    puw_roof(jj)    = - zustar_roof(jj)**2
    pduwdu_roof(jj) = 0.
    IF (pvmod(jj)/=0.) pduwdu_roof(jj) = 2. * puw_roof(jj) / pvmod(jj)
    !
  ELSE
    !
    !*      9.3    For town
    !              --------
    !
    zustar2(jj)=xundef
    !  
    IF(himplicit_wind=='OLD')THEN
      !   old implicitation
      zustar2(jj) = (pcd(jj)*pvmod(jj)*ppew_b_coef(jj))/            &
                    (1.0-prhoa(jj)*pcd(jj)*pvmod(jj)*ppew_a_coef(jj))
    ELSE
      !   new implicitation
      zustar2(jj) = (pcd(jj)*pvmod(jj)*(2.*ppew_b_coef(jj)-pvmod(jj)))/ &
                    (1.0-2.0*prhoa(jj)*pcd(jj)*pvmod(jj)*ppew_a_coef(jj)) 
      !                   
      zvmod(jj) = prhoa(jj)*ppew_a_coef(jj)*zustar2(jj) + ppew_b_coef(jj)
      zvmod(jj) = max(zvmod(jj),0.)
      !
      IF(ppew_a_coef(jj)/= 0.)THEN
        zustar2(jj) = max( ( zvmod(jj) - ppew_b_coef(jj) ) / (prhoa(jj)*ppew_a_coef(jj)), 0.)
      ENDIF
      !                        
    ENDIF
    !
    pustar_town(jj) = sqrt(zustar2(jj))
    !
  ENDIF
  !
ENDDO
!
IF (lhook) CALL dr_hook('URBAN_DRAG',1,zhook_handle)
!-------------------------------------------------------------------------------
!
END SUBROUTINE urban_drag