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(TOP, T, B, 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, 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, PAC_WIN ) 
!   ##########################################################################
!
!!****  *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 (TOP%CZ0H) 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_teb_option_n, ONLY : teb_options_t
USE modd_teb_n, ONLY : teb_t
USE modd_bem_n, ONLY : bem_t
!
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
!
TYPE(teb_options_t), INTENT(INOUT) :: TOP
TYPE(teb_t), INTENT(INOUT) :: T
TYPE(bem_t), INTENT(INOUT) :: B
!
 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)    :: 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
!
REAL, DIMENSION(:), INTENT(OUT)   :: PAC_WIN        ! aerodynamical conductance for window
!
!*      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(:) = t%XZ0_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(:) = (t%XWS_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(:) = (t%XWS_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) = t%XBLD(jj) * pts_roof(jj) + (1.-t%XBLD(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) =  t%XBLD(jj) * zq_roof(jj) + (1.-t%XBLD(jj)) * pq_canyon(jj)
!
ENDDO
!
!-------------------------------------------------------------------------------
!
!*      5.     Momentum drag coefficient
!              -------------------------
!
IF (.NOT. top%LCANOPY) THEN
  CALL urban_exch_coef(top%CZ0H, zz0_o_z0h, zts_town, zq_town, pexns, pexna, pta, pqa,     &
                       pzref+ t%XBLD_HEIGHT/3., puref+t%XBLD_HEIGHT/3., pvmod, t%XZ0_TOWN, &
                       pri, pcd, pcdn, zac, zra, zch                              )
ENDIF
!
!-------------------------------------------------------------------------------
!
!*      6.     Drag coefficient for heat fluxes between roofs and atmosphere
!              -------------------------------------------------------------
!
IF (top%CCH_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, t%XROUGH_ROOF) / prhoa / xcpd
ELSE
   CALL urban_exch_coef(top%CZ0H, 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)     = t%XWS_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. top%LCANOPY) THEN
  CALL urban_exch_coef('MASC95', 1., pt_canyon, pq_canyon, pexns, pexna, pta, pqa,    &
                        pzref+t%XBLD_HEIGHT-pz_lowcan, puref+t%XBLD_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 (top%CCH_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), t%XROUGH_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(top%CZ0H, zz0_o_z0h, pts_road, pq_lowcan, pexns, pexna,  &
                       pt_lowcan, pq_lowcan, pz_lowcan, pz_lowcan,              &
                       pu_lowcan+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) * t%XWALL_O_GRND(jj)

    IF (t%XROAD(jj) .GT. 0.) THEN
      zq0(jj)   = zq0(jj) &
            + (pts_road(jj) - pt_lowcan(jj)) * pac_road(jj) * t%XROAD  (jj)/(t%XROAD(jj)+t%XGARDEN(jj)) 
    ENDIF
    IF (t%XGARDEN(jj) .GT. 0.) THEN
      zq0(jj)   = zq0(jj) &
            + (pts_garden(jj) - pt_lowcan(jj)) * pac_garden(jj) * t%XGARDEN(jj)/(t%XROAD(jj)+t%XGARDEN(jj))
    ENDIF    
    !
    IF (zq0(jj) >= 0.) THEN
      zw_star(jj) = ( (xg * pexna(jj) / pta(jj)) * zq0(jj) * t%XBLD_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)     = t%XWS_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(b%XT_WIN1(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 (top%LCANOPY) 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