wrfv2_fire/phys/module_sf_urban.F

   1 MODULE module_sf_urban
   2
   3 !===============================================================================
   4 !     Single-Layer Urban Canopy Model for WRF Noah-LSM
   5 !     Original Version: 2002/11/06 by Hiroyuki Kusaka
   6 !     Last Update:      2006/08/24 by Fei Chen and Mukul Tewari (NCAR/RAL)
   7 !===============================================================================
   8
   9    CHARACTER(LEN=4)                :: LU_DATA_TYPE
  10
  11    INTEGER                         :: ICATE
  12
  13    REAL, ALLOCATABLE, DIMENSION(:) :: ZR_TBL
  14    REAL, ALLOCATABLE, DIMENSION(:) :: Z0C_TBL
  15    REAL, ALLOCATABLE, DIMENSION(:) :: Z0HC_TBL
  16    REAL, ALLOCATABLE, DIMENSION(:) :: ZDC_TBL
  17    REAL, ALLOCATABLE, DIMENSION(:) :: SVF_TBL
  18    REAL, ALLOCATABLE, DIMENSION(:) :: R_TBL
  19    REAL, ALLOCATABLE, DIMENSION(:) :: RW_TBL
  20    REAL, ALLOCATABLE, DIMENSION(:) :: HGT_TBL
  21    REAL, ALLOCATABLE, DIMENSION(:) :: CDS_TBL
  22    REAL, ALLOCATABLE, DIMENSION(:) :: AS_TBL
  23    REAL, ALLOCATABLE, DIMENSION(:) :: AH_TBL
  24    REAL, ALLOCATABLE, DIMENSION(:) :: BETR_TBL
  25    REAL, ALLOCATABLE, DIMENSION(:) :: BETB_TBL
  26    REAL, ALLOCATABLE, DIMENSION(:) :: BETG_TBL
  27    REAL, ALLOCATABLE, DIMENSION(:) :: FRC_URB_TBL
  28
  29    REAL                            :: CAPR_DATA, CAPB_DATA, CAPG_DATA
  30    REAL                            :: AKSR_DATA, AKSB_DATA, AKSG_DATA
  31    REAL                            :: ALBR_DATA, ALBB_DATA, ALBG_DATA
  32    REAL                            :: EPSR_DATA, EPSB_DATA, EPSG_DATA
  33    REAL                            :: Z0R_DATA,  Z0B_DATA,  Z0G_DATA
  34    REAL                            :: Z0HR_DATA, Z0HB_DATA, Z0HG_DATA
  35    REAL                            :: TRLEND_DATA, TBLEND_DATA, TGLEND_DATA
  36
  37    INTEGER                         :: BOUNDR_DATA,BOUNDB_DATA,BOUNDG_DATA
  38    INTEGER                         :: CH_SCHEME_DATA, TS_SCHEME_DATA
  39    INTEGER                         :: ahoption        ! Miao, 2007/01/17, cal. ah
  40    REAL, DIMENSION(1:24)           :: ahdiuprf        ! ah diurnal profile, tloc: 1-24
  41
  42    INTEGER                         :: allocate_status
  43
  44 !   INTEGER                         :: num_roof_layers
  45 !   INTEGER                         :: num_wall_layers
  46 !   INTEGER                         :: num_road_layers
  47
  48    CONTAINS
  49
  50 !===============================================================================
  51 !
  52 ! Author:
  53 !      Hiroyuki KUSAKA, PhD
  54 !      University of Tsukuba, JAPAN
  55 !      (CRIEPI, NCAR/MMM visiting scientist, 2002-2004)
  56 !      kusaka@ccs.tsukuba.ac.jp
  57 !
  58 ! Co-Researchers:
  59 !     Fei CHEN, PhD
  60 !      NCAR/RAP feichen@ucar.edu
  61 !     Mukul TEWARI, PhD
  62 !      NCAR/RAP mukul@ucar.edu
  63 !
  64 ! Purpose:
  65 !     Calculate surface temeprature, fluxes, canopy air temperature, and canopy wind
  66 !
  67 ! Subroutines:
  68 !     module_sf_urban
  69 !       |- urban
  70 !            |- read_param
  71 !            |- mos or jurges
  72 !            |- multi_layer or force_restore
  73 !       |- urban_param_init <-- urban_param.tbl
  74 !       |- urban_var_init
  75 !
  76 ! Input Data from WRF [MKS unit]:
  77 !
  78 !     UTYPE  [-]     : Urban type. 1=urban, 2=suburban, 3=rural
  79 !     TA     [K]     : Potential temperature at 1st wrf level (absolute temp)
  80 !     QA     [kg/kg] : Mixing ratio at 1st atmospheric level
  81 !     UA     [m/s]   : Wind speed at 1st atmospheric level
  82 !     SSG    [W/m/m] : Short wave downward radiation at a flat surface
  83 !                      Note this is the total of direct and diffusive solar
  84 !                      downward radiation. If without two components, the
  85 !                      single solar downward can be used instead.
  86 !                      SSG = SSGD + SSGQ
  87 !     LSOLAR [-]     : Indicating the input type of solar downward radiation
  88 !                      True: both direct and diffusive solar radiation
  89 !                      are available
  90 !                      False: only total downward ridiation is available.
  91 !     SSGD   [W/m/m] : Direct solar radiation at a flat surface
  92 !                      if SSGD is not available, one can assume a ratio SRATIO
  93 !                      (e.g., 0.7), so that SSGD = SRATIO*SSG
  94 !     SSGQ   [W/m/m] : Diffuse solar radiation at a flat surface
  95 !                      If SSGQ is not available, SSGQ = SSG - SSGD
  96 !     LLG    [W/m/m] : Long wave downward radiation at a flat surface
  97 !     RAIN   [mm/h]  : Precipitation
  98 !     RHOO   [kg/m/m/m] : Air density
  99 !     ZA     [m]     : First atmospheric level
 100 !                      as a lowest boundary condition
 101 !     DECLIN [rad]   : solar declination
 102 !     COSZ           : = sin(fai)*sin(del)+cos(fai)*cos(del)*cos(omg)
 103 !     OMG    [rad]   : solar hour angle
 104 !     XLAT   [deg]   : latitude
 105 !     DELT   [sec]   : Time step
 106 !     ZNT    [m]     : Roughnes length
 107 !
 108 ! Output Data to WRF [MKS unit]:
 109 !
 110 !     TS  [K]            : Surface potential temperature (absolute temp)
 111 !     QS  [-]            : Surface humidity
 112 !
 113 !     SH  [W/m/m/]       : Sensible heat flux, = FLXTH*RHOO*CPP
 114 !     LH  [W/m/m]        : Latent heat flux, = FLXHUM*RHOO*ELL
 115 !     LH_INEMATIC [kg/m/m/sec]: Moisture Kinematic flux, = FLXHUM*RHOO
 116 !     SW  [W/m/m]        : Upward shortwave radiation flux,
 117 !                          = SSG-SNET*697.7*60. (697.7*60.=100.*100.*4.186)
 118 !     ALB [-]            : Time-varying albedo
 119 !     LW  [W/m/m]        : Upward longwave radiation flux,
 120 !                          = LNET*697.7*60.-LLG
 121 !     G   [W/m/m]        : Heat Flux into the Ground
 122 !     RN  [W/m/m]        : Net radiation
 123 !
 124 !     PSIM  [-]          : Diagnostic similarity stability function for momentum
 125 !     PSIH  [-]          : Diagnostic similarity stability function for heat
 126 !
 127 !     TC  [K]            : Diagnostic canopy air temperature
 128 !     QC  [-]            : Diagnostic canopy humidity
 129 !
 130 !     TH2 [K]            : Diagnostic potential temperature at 2 m
 131 !     Q2  [-]            : Diagnostic humidity at 2 m
 132 !     U10 [m/s]          : Diagnostic u wind component at 10 m
 133 !     V10 [m/s]          : Diagnostic v wind component at 10 m
 134 !
 135 !     CHS, CHS2 [m/s]    : CH*U at ZA, CH*U at 2 m (not used)
 136 !
 137 ! Important parameters:
 138 ! Following parameter are assigned in run/urban_param.tbl
 139 !
 140 !  ZR  [m]             : roof level (building height)
 141 !  Z0C [m]             : Roughness length above canyon for momentum (1/10 of ZR)
 142 !  Z0HC [m]            : Roughness length above canyon for heat (1/10 of Z0C)
 143 !  ZDC [m]             : Zero plane displacement height (1/5 of ZR)
 144 !  SVF [-]             : sky view factor. Calculated again in urban_param_init
 145 !  R   [-]             : building coverage ratio
 146 !  RW  [-]             : = 1 - R
 147 !  HGT [-]             : normalized building height
 148 !  CDS [-]             : drag coefficient by buildings
 149 !  AS  [1/m]           : buildings volumetric parameter
 150 !  AH  [cal/cm/cm]     : anthropogenic heat
 151 !  BETR [-]            : minimum moisture availability of roof
 152 !  BETB [-]            : minimum moisture availability of building wall
 153 !  BETG [-]            : minimum moisture availability of road
 154 !  CAPR[cal/cm/cm/cm/degC]: heat capacity of roof
 155 !  CAPB[cal/cm/cm/cm/degC]: heat capacity of building wall
 156 !  CAPG[cal/cm/cm/cm/degC]: heat capacity of road
 157 !  AKSR [cal/cm/sec/degC] : thermal conductivity of roof
 158 !  AKSB [cal/cm/sec/degC] : thermal conductivity of building wall
 159 !  AKSG [cal/cm/sec/degC] : thermal conductivity of road
 160 !  ALBR [-]               : surface albedo of roof
 161 !  ALBB [-]               : surface albedo of building wall
 162 !  ALBG [-]               : surface albedo of road
 163 !  EPSR [-]               : surface emissivity of roof
 164 !  EPSB [-]               : surface emissivity of building wall
 165 !  EPSG [-]               : surface emissivity of road
 166 !  Z0R [m]                : roughness length for momentum of roof
 167 !  Z0B [m]                : roughness length for momentum of building wall (only for CH_SCHEME = 1)
 168 !  Z0G [m]                : roughness length for momentum of road (only for CH_SCHEME = 1)
 169 !  Z0HR [m]               : roughness length for heat of roof
 170 !  Z0HB [m]               : roughness length for heat of building wall (only for CH_SCHEME = 1)
 171 !  Z0HG [m]               : roughness length for heat of road
 172 !  num_roof_layers        : number of layers within roof
 173 !  num_wall_layers        : number of layers within building walls
 174 !  num_road_layers        : number of layers within below road surface
 175 !   NOTE: for now, these layers are defined as same as the number of soil layers in namelist.input
 176 !  DZR [cm]               : thickness of each roof layer
 177 !  DZB [cm]               : thickness of each building wall layer
 178 !  DZG [cm]               : thickness of each ground layer
 179 !  BOUNDR [integer 1 or 2] : Boundary Condition for Roof Layer Temp [1: Zero-Flux, 2: T = Constant]
 180 !  BOUNDB [integer 1 or 2] : Boundary Condition for Building Wall Layer Temp [1: Zero-Flux, 2: T = Constant]
 181 !  BOUNDG [integer 1 or 2] : Boundary Condition for Road Layer Temp [1: Zero-Flux, 2: T = Constant]
 182 !  TRLEND [K]             : lower boundary condition of roof temperature
 183 !  TBLEND [K]             : lower boundary condition of building temperature
 184 !  TGLEND [K]             : lower boundary condition of ground temperature
 185 !  CH_SCHEME [integer 1 or 2] : Sfc exchange scheme used for building wall and road
 186 !                         [1: M-O Similarity Theory,  2: Empirical Form (recommend)]
 187 !  TS_SCHEME [integer 1 or 2] : Scheme for computing surface temperature (for roof, wall, and road)
 188 !                         [1: 4-layer model,  2: Force-Restore method]
 189 !
 190 !
 191 ! References:
 192 !     Kusaka and Kimura (2004) J.Appl.Meteor., vol.43, p1899-1910
 193 !     Kusaka and Kimura (2004) J.Meteor.Soc.Japan, vol.82, p45-65
 194 !     Kusaka et al. (2001) Bound.-Layer Meteor., vol.101, p329-358
 195 !
 196 ! History:
 197 !     2006/06          modified by H. Kusaka (Univ. Tsukuba), M. Tewari
 198 !     2005/10/26,      modified by Fei Chen, Mukul Tewari
 199 !     2003/07/21 WRF , modified  by H. Kusaka of CRIEPI (NCAR/MMM)
 200 !     2001/08/26 PhD , modified  by H. Kusaka of CRIEPI (Univ.Tsukuba)
 201 !     1999/08/25 LCM , developed by H. Kusaka of CRIEPI (Univ.Tsukuba)
 202 !
 203 !===============================================================================
 204 !
 205 !  subroutine urban:
 206 !
 207 !===============================================================================
 208
 209    SUBROUTINE urban(LSOLAR,                                           & ! L
 210                     num_roof_layers,num_wall_layers,num_road_layers,  & ! I
 211                     DZR,DZB,DZG,                                      & ! I
 212                     UTYPE,TA,QA,UA,U1,V1,SSG,SSGD,SSGQ,LLG,RAIN,RHOO, & ! I
 213                     ZA,DECLIN,COSZ,OMG,XLAT,DELT,ZNT,                 & ! I
 214                     CHS, CHS2,                                        & ! I
 215                     TR, TB, TG, TC, QC, UC,                           & ! H
 216                     TRL,TBL,TGL,                                      & ! H
 217                     XXXR, XXXB, XXXG, XXXC,                           & ! H
 218                     TS,QS,SH,LH,LH_KINEMATIC,                         & ! O
 219                     SW,ALB,LW,G,RN,PSIM,PSIH,                         & ! O
 220                     GZ1OZ0,                                           & ! O
 221                     U10,V10,TH2,Q2,UST                                & ! O
 222                     )
 223
 224    IMPLICIT NONE
 225
 226    REAL, PARAMETER    :: CP=0.24          ! heat capacity of dry air  [cgs unit]
 227    REAL, PARAMETER    :: EL=583.          ! latent heat of vaporation [cgs unit]
 228    REAL, PARAMETER    :: SIG=8.17E-11     ! stefun bolzman constant   [cgs unit]
 229    REAL, PARAMETER    :: SIG_SI=5.67E-8   !                           [MKS unit]
 230    REAL, PARAMETER    :: AK=0.4           ! kalman const.                [-]
 231    REAL, PARAMETER    :: PI=3.14159       ! pi                           [-]
 232    REAL, PARAMETER    :: TETENA=7.5       ! const. of Tetens Equation    [-]
 233    REAL, PARAMETER    :: TETENB=237.3     ! const. of Tetens Equation    [-]
 234    REAL, PARAMETER    :: SRATIO=0.75      ! ratio between direct/total solar [-]
 235
 236    REAL, PARAMETER    :: CPP=1004.5       ! heat capacity of dry air    [J/K/kg]
 237    REAL, PARAMETER    :: ELL=2.442E+06    ! latent heat of vaporization [J/kg]
 238    REAL, PARAMETER    :: XKA=2.4E-5
 239
 240 !-------------------------------------------------------------------------------
 241 ! C: configuration variables
 242 !-------------------------------------------------------------------------------
 243
 244    LOGICAL, INTENT(IN) :: LSOLAR  ! logical    [true=both, false=SSG only]
 245
 246 !  The following variables are also model configuration variables, but are
 247 !  defined in the URBAN.TBL and in the contains statement in the top of
 248 !  the module_urban_init, so we should not declare them here.
 249
 250   INTEGER, INTENT(IN) :: num_roof_layers
 251   INTEGER, INTENT(IN) :: num_wall_layers
 252   INTEGER, INTENT(IN) :: num_road_layers
 253
 254
 255   REAL, INTENT(IN), DIMENSION(1:num_roof_layers) :: DZR ! grid interval of roof layers [cm]
 256   REAL, INTENT(IN), DIMENSION(1:num_wall_layers) :: DZB ! grid interval of wall layers [cm]
 257   REAL, INTENT(IN), DIMENSION(1:num_road_layers) :: DZG ! grid interval of road layers [cm]
 258
 259 !-------------------------------------------------------------------------------
 260 ! I: input variables from LSM to Urban
 261 !-------------------------------------------------------------------------------
 262
 263    INTEGER, INTENT(IN) :: UTYPE ! urban type [urban=1, suburban=2, rural=3]
 264
 265    REAL, INTENT(IN)    :: TA   ! potential temp at 1st atmospheric level [K]
 266    REAL, INTENT(IN)    :: QA   ! mixing ratio at 1st atmospheric level  [kg/kg]
 267    REAL, INTENT(IN)    :: UA   ! wind speed at 1st atmospheric level    [m/s]
 268    REAL, INTENT(IN)    :: U1   ! u at 1st atmospheric level             [m/s]
 269    REAL, INTENT(IN)    :: V1   ! v at 1st atmospheric level             [m/s]
 270    REAL, INTENT(IN)    :: SSG  ! downward total short wave radiation    [W/m/m]
 271    REAL, INTENT(IN)    :: LLG  ! downward long wave radiation           [W/m/m]
 272    REAL, INTENT(IN)    :: RAIN ! precipitation                          [mm/h]
 273    REAL, INTENT(IN)    :: RHOO ! air density                            [kg/m^3]
 274    REAL, INTENT(IN)    :: ZA   ! first atmospheric level                [m]
 275    REAL, INTENT(IN)    :: DECLIN ! solar declination                    [rad]
 276    REAL, INTENT(IN)    :: COSZ ! sin(fai)*sin(del)+cos(fai)*cos(del)*cos(omg)
 277    REAL, INTENT(IN)    :: OMG  ! solar hour angle                       [rad]
 278
 279    REAL, INTENT(IN)    :: XLAT ! latitude                               [deg]
 280    REAL, INTENT(IN)    :: DELT ! time step                              [s]
 281    REAL, INTENT(IN)    :: ZNT  ! roughness length                       [m]
 282    REAL, INTENT(IN)    :: CHS,CHS2 ! CH*U at za and 2 m             [m/s]
 283
 284    REAL, INTENT(INOUT) :: SSGD ! downward direct short wave radiation   [W/m/m]
 285    REAL, INTENT(INOUT) :: SSGQ ! downward diffuse short wave radiation  [W/m/m]
 286
 287 !-------------------------------------------------------------------------------
 288 ! O: output variables from Urban to LSM
 289 !-------------------------------------------------------------------------------
 290
 291    REAL, INTENT(OUT) :: TS     ! surface potential temperature    [K]
 292    REAL, INTENT(OUT) :: QS     ! surface humidity                 [K]
 293    REAL, INTENT(OUT) :: SH     ! sensible heat flux               [W/m/m]
 294    REAL, INTENT(OUT) :: LH     ! latent heat flux                 [W/m/m]
 295    REAL, INTENT(OUT) :: LH_KINEMATIC ! latent heat, kinetic     [kg/m/m/s]
 296    REAL, INTENT(OUT) :: SW     ! upward short wave radiation flux [W/m/m]
 297    REAL, INTENT(OUT) :: ALB    ! time-varying albedo            [fraction]
 298    REAL, INTENT(OUT) :: LW     ! upward long wave radiation flux  [W/m/m]
 299    REAL, INTENT(OUT) :: G      ! heat flux into the ground        [W/m/m]
 300    REAL, INTENT(OUT) :: RN     ! net radition                     [W/m/m]
 301    REAL, INTENT(OUT) :: PSIM   ! similality stability shear function for momentum
 302    REAL, INTENT(OUT) :: PSIH   ! similality stability shear function for heat
 303    REAL, INTENT(OUT) :: GZ1OZ0
 304    REAL, INTENT(OUT) :: U10    ! u at 10m                         [m/s]
 305    REAL, INTENT(OUT) :: V10    ! u at 10m                         [m/s]
 306    REAL, INTENT(OUT) :: TH2    ! potential temperature at 2 m     [K]
 307    REAL, INTENT(OUT) :: Q2     ! humidity at 2 m                  [-]
 308 !m   REAL, INTENT(OUT) :: CHS,CHS2 ! CH*U at za and 2 m             [m/s]
 309    REAL, INTENT(OUT) :: UST    ! friction velocity                [m/s]
 310
 311
 312 !-------------------------------------------------------------------------------
 313 ! H: Historical (state) variables of Urban : LSM <--> Urban
 314 !-------------------------------------------------------------------------------
 315
 316 ! TR: roof temperature              [K];      TRP: at previous time step [K]
 317 ! TB: building wall temperature     [K];      TBP: at previous time step [K]
 318 ! TG: road temperature              [K];      TGP: at previous time step [K]
 319 ! TC: urban-canopy air temperature  [K];      TCP: at previous time step [K]
 320 !                                                  (absolute temperature)
 321 ! QC: urban-canopy air mixing ratio [kg/kg];  QCP: at previous time step [kg/kg]
 322 !
 323 ! XXXR: Monin-Obkhov length for roof          [dimensionless]
 324 ! XXXB: Monin-Obkhov length for building wall [dimensionless]
 325 ! XXXG: Monin-Obkhov length for road          [dimensionless]
 326 ! XXXC: Monin-Obkhov length for urban-canopy  [dimensionless]
 327 !
 328 ! TRL, TBL, TGL: layer temperature [K] (absolute temperature)
 329
 330    REAL, INTENT(INOUT):: TR, TB, TG, TC, QC, UC
 331    REAL, INTENT(INOUT):: XXXR, XXXB, XXXG, XXXC
 332
 333    REAL, DIMENSION(1:num_roof_layers), INTENT(INOUT) :: TRL
 334    REAL, DIMENSION(1:num_wall_layers), INTENT(INOUT) :: TBL
 335    REAL, DIMENSION(1:num_road_layers), INTENT(INOUT) :: TGL
 336
 337 !-------------------------------------------------------------------------------
 338 ! L:  Local variables from read_param
 339 !-------------------------------------------------------------------------------
 340
 341    REAL :: ZR, Z0C, Z0HC, ZDC, SVF, R, RW, HGT, CDS, AS, AH
 342    REAL :: CAPR, CAPB, CAPG, AKSR, AKSB, AKSG, ALBR, ALBB, ALBG
 343    REAL :: EPSR, EPSB, EPSG, Z0R,  Z0B,  Z0G,  Z0HR, Z0HB, Z0HG
 344    REAL :: TRLEND,TBLEND,TGLEND
 345
 346    REAL :: TH2X                                                !m
 347
 348    INTEGER :: BOUNDR, BOUNDB, BOUNDG
 349    INTEGER :: CH_SCHEME, TS_SCHEME
 350
 351    LOGICAL :: SHADOW  ! [true=consider svf and shadow effects, false=consider svf effect only]
 352
 353 !-------------------------------------------------------------------------------
 354 ! L: Local variables
 355 !-------------------------------------------------------------------------------
 356
 357    REAL :: BETR, BETB, BETG
 358    REAL :: SX, SD, SQ, RX
 359    REAL :: UR, ZC, XLB, BB
 360    REAL :: Z, RIBR, RIBB, RIBG, RIBC, BHR, BHB, BHG, BHC
 361    REAL :: TSC, LNET, SNET, FLXUV, THG, FLXTH, FLXHUM, FLXG
 362    REAL :: W, VFGS, VFGW, VFWG, VFWS, VFWW
 363    REAL :: HOUI1, HOUI2, HOUI3, HOUI4, HOUI5, HOUI6, HOUI7, HOUI8
 364    REAL :: SLX, SLX1, SLX2, SLX3, SLX4, SLX5, SLX6, SLX7, SLX8
 365    REAL :: FLXTHR, FLXTHB, FLXTHG, FLXHUMR, FLXHUMB, FLXHUMG
 366    REAL :: SR, SB, SG, RR, RB, RG
 367    REAL :: SR1, SR2, SB1, SB2, SG1, SG2, RR1, RR2, RB1, RB2, RG1, RG2
 368    REAL :: HR, HB, HG, ELER, ELEB, ELEG, G0R, G0B, G0G
 369    REAL :: ALPHAC, ALPHAR, ALPHAB, ALPHAG
 370    REAL :: CHC, CHR, CHB, CHG, CDC, CDR, CDB, CDG
 371    REAL :: C1R, C1B, C1G, TE, TC1, TC2, QC1, QC2, QS0R, QS0B, QS0G,RHO,ES
 372
 373    REAL :: DESDT
 374    REAL :: F
 375    REAL :: DQS0RDTR
 376    REAL :: DRRDTR, DHRDTR, DELERDTR, DG0RDTR
 377    REAL :: DTR, DFDT
 378    REAL :: FX, FY, GF, GX, GY
 379    REAL :: DTCDTB, DTCDTG
 380    REAL :: DQCDTB, DQCDTG
 381    REAL :: DRBDTB1,  DRBDTG1,  DRBDTB2,  DRBDTG2
 382    REAL :: DRGDTB1,  DRGDTG1,  DRGDTB2,  DRGDTG2
 383    REAL :: DRBDTB,   DRBDTG,   DRGDTB,   DRGDTG
 384    REAL :: DHBDTB,   DHBDTG,   DHGDTB,   DHGDTG
 385    REAL :: DELEBDTB, DELEBDTG, DELEGDTG, DELEGDTB
 386    REAL :: DG0BDTB,  DG0BDTG,  DG0GDTG,  DG0GDTB
 387    REAL :: DQS0BDTB, DQS0GDTG
 388    REAL :: DTB, DTG, DTC
 389
 390    REAL :: THEATAZ    ! Solar Zenith Angle [rad]
 391    REAL :: THEATAS    ! = PI/2. - THETAZ
 392    REAL :: FAI        ! Latitude [rad]
 393    REAL :: CNT,SNT
 394    REAL :: PS         ! Surface Pressure [hPa]
 395    REAL :: TAV        ! Vertial Temperature [K]
 396
 397    REAL :: XXX, X, Z0, Z0H, CD, CH
 398    REAL :: XXX2, PSIM2, PSIH2, XXX10, PSIM10, PSIH10
 399    REAL :: PSIX, PSIT, PSIX2, PSIT2, PSIX10, PSIT10
 400
 401    REAL :: TRP, TBP, TGP, TCP, QCP, TST, QST
 402
 403    INTEGER :: iteration, K
 404    INTEGER :: tloc
 405
 406 !-------------------------------------------------------------------------------
 407 ! Set parameters
 408 !-------------------------------------------------------------------------------
 409
 410 ! Miao, 2007/01/17, cal. ah
 411    if(ahoption==1) then
 412      tloc=mod(int(OMG/PI*180./15.+12.+0.5 ),24)
 413      if(tloc==0) tloc=24
 414    endif
 415
 416    CALL read_param(UTYPE,ZR,Z0C,Z0HC,ZDC,SVF,R,RW,HGT,CDS,AS,AH,  &
 417                    CAPR,CAPB,CAPG,AKSR,AKSB,AKSG,ALBR,ALBB,ALBG,  &
 418                    EPSR,EPSB,EPSG,Z0R,Z0B,Z0G,Z0HR,Z0HB,Z0HG,     &
 419                    BETR,BETB,BETG,TRLEND,TBLEND,TGLEND,           &
 420                    BOUNDR,BOUNDB,BOUNDG,CH_SCHEME,TS_SCHEME)
 421
 422 ! Miao, 2007/01/17, cal. ah
 423    if(ahoption==1) AH=AH*ahdiuprf(tloc)
 424
 425    IF( ZDC+Z0C+2. >= ZA) THEN
 426     PRINT *, 'ZDC + Z0C + 2m is larger than the 1st WRF level'
 427     PRINT *, 'Stop in the subroutine urban - change ZDC and Z0C'
 428     STOP
 429    END IF
 430
 431    IF(.NOT.LSOLAR) THEN
 432      SSGD = SRATIO*SSG
 433      SSGQ = SSG - SSGD
 434    ENDIF
 435    SSGD = SRATIO*SSG   ! No radiation scheme has SSGD and SSGQ.
 436    SSGQ = SSG - SSGD
 437
 438    W=2.*1.*HGT
 439    VFGS=SVF
 440    VFGW=1.-SVF
 441    VFWG=(1.-SVF)*(1.-R)/W
 442    VFWS=VFWG
 443    VFWW=1.-2.*VFWG
 444
 445 !-------------------------------------------------------------------------------
 446 ! Convert unit from MKS to cgs
 447 ! Renew surface and layer temperatures
 448 !-------------------------------------------------------------------------------
 449
 450    SX=(SSGD+SSGQ)/697.7/60.  ! downward short wave radition [ly/min]
 451    SD=SSGD/697.7/60.         ! downward direct short wave radiation
 452    SQ=SSGQ/697.7/60.         ! downward diffiusion short wave radiation
 453    RX=LLG/697.7/60.          ! downward long wave radiation
 454    RHO=RHOO*0.001            ! air density at first atmospheric level
 455
 456    TRP=TR
 457    TBP=TB
 458    TGP=TG
 459    TCP=TC
 460    QCP=QC
 461
 462    TAV=TA*(1.+0.61*QA)
 463    PS=RHOO*287.*TAV/100. ![hPa]
 464
 465 !-------------------------------------------------------------------------------
 466 ! Canopy wind
 467 !-------------------------------------------------------------------------------
 468
 469    IF ( ZR + 2. < ZA ) THEN
 470      UR=UA*LOG((ZR-ZDC)/Z0C)/LOG((ZA-ZDC)/Z0C)
 471      ZC=0.7*ZR
 472 !     ZC=0.5*ZR
 473      XLB=0.4*(ZR-ZDC)
 474      BB=ZR*(CDS*AS/(2.*XLB**2.))**(1./3.)
 475      UC=UR*EXP(-BB*(1.-ZC/ZR))
 476    ELSE
 477      print *,'ZR=',ZR, 'ZA=',ZA
 478      PRINT *, 'Warning ZR + 2m  is larger than the 1st WRF level'
 479      ZC=ZA/2.
 480      UC=UA/2.
 481    END IF
 482
 483 !-------------------------------------------------------------------------------
 484 ! Net Short Wave Radiation at roof, wall, and road
 485 !-------------------------------------------------------------------------------
 486
 487    SHADOW = .false.
 488 !   SHADOW = .true.
 489
 490    IF (SSG > 0.0) THEN
 491
 492      IF(.NOT.SHADOW) THEN              ! no shadow effects model
 493
 494        SR1=SX*(1.-ALBR)
 495        SG1=SX*VFGS*(1.-ALBG)
 496        SB1=SX*VFWS*(1.-ALBB)
 497        SG2=SB1*ALBB/(1.-ALBB)*VFGW*(1.-ALBG)
 498        SB2=SG1*ALBG/(1.-ALBG)*VFWG*(1.-ALBB)
 499
 500      ELSE                              ! shadow effects model
 501
 502        FAI=XLAT*PI/180.
 503
 504        THEATAS=ABS(ASIN(COSZ))
 505        THEATAZ=ABS(ACOS(COSZ))
 506
 507        SNT=COS(DECLIN)*SIN(OMG)/COS(THEATAS)
 508        CNT=(COSZ*SIN(FAI)-SIN(DECLIN))/COS(THEATAS)/COS(FAI)
 509
 510        HOUI1=(SNT*COS(PI/8.)    -CNT*SIN(PI/8.))
 511        HOUI2=(SNT*COS(2.*PI/8.) -CNT*SIN(2.*PI/8.))
 512        HOUI3=(SNT*COS(3.*PI/8.) -CNT*SIN(3.*PI/8.))
 513        HOUI4=(SNT*COS(4.*PI/8.) -CNT*SIN(4.*PI/8.))
 514        HOUI5=(SNT*COS(5.*PI/8.) -CNT*SIN(5.*PI/8.))
 515        HOUI6=(SNT*COS(6.*PI/8.) -CNT*SIN(6.*PI/8.))
 516        HOUI7=(SNT*COS(7.*PI/8.) -CNT*SIN(7.*PI/8.))
 517        HOUI8=(SNT*COS(8.*PI/8.) -CNT*SIN(8.*PI/8.))
 518
 519        SLX1=HGT*ABS(TAN(THEATAZ))*ABS(HOUI1)
 520        SLX2=HGT*ABS(TAN(THEATAZ))*ABS(HOUI2)
 521        SLX3=HGT*ABS(TAN(THEATAZ))*ABS(HOUI3)
 522        SLX4=HGT*ABS(TAN(THEATAZ))*ABS(HOUI4)
 523        SLX5=HGT*ABS(TAN(THEATAZ))*ABS(HOUI5)
 524        SLX6=HGT*ABS(TAN(THEATAZ))*ABS(HOUI6)
 525        SLX7=HGT*ABS(TAN(THEATAZ))*ABS(HOUI7)
 526        SLX8=HGT*ABS(TAN(THEATAZ))*ABS(HOUI8)
 527
 528        IF(SLX1 > RW) SLX1=RW
 529        IF(SLX2 > RW) SLX2=RW
 530        IF(SLX3 > RW) SLX3=RW
 531        IF(SLX4 > RW) SLX4=RW
 532        IF(SLX5 > RW) SLX5=RW
 533        IF(SLX6 > RW) SLX6=RW
 534        IF(SLX7 > RW) SLX7=RW
 535        IF(SLX8 > RW) SLX8=RW
 536
 537        SLX=(SLX1+SLX2+SLX3+SLX4+SLX5+SLX6+SLX7+SLX8)/8.
 538
 539      END IF
 540
 541      SR=SR1
 542      SG=SG1+SG2
 543      SB=SB1+SB2
 544
 545      SNET=R*SR+W*SB+RW*SG
 546
 547    ELSE
 548
 549      SR=0.
 550      SG=0.
 551      SB=0.
 552      SNET=0.
 553
 554    END IF
 555
 556 !-------------------------------------------------------------------------------
 557 ! Roof
 558 !-------------------------------------------------------------------------------
 559
 560 !-------------------------------------------------------------------------------
 561 ! CHR, CDR, BETR
 562 !-------------------------------------------------------------------------------
 563
 564    Z=ZA-ZDC
 565    BHR=LOG(Z0R/Z0HR)/0.4
 566    RIBR=(9.8*2./(TA+TRP))*(TA-TRP)*(Z+Z0R)/(UA*UA)
 567
 568    CALL mos(XXXR,ALPHAR,CDR,BHR,RIBR,Z,Z0R,UA,TA,TRP,RHO)
 569
 570    CHR=ALPHAR/RHO/CP/UA
 571
 572    IF(RAIN > 1.) BETR=0.7
 573
 574    IF (TS_SCHEME == 1) THEN
 575
 576 !-------------------------------------------------------------------------------
 577 ! TR  Solving Non-Linear Equation by Newton-Rapson
 578 ! TRL Solving Heat Equation by Tri Diagonal Matrix Algorithm
 579 !-------------------------------------------------------------------------------
 580 !      TSC=TRP-273.15
 581 !      ES=EXP(19.482-4303.4/(TSC+243.5))        ! WMO
 582 !      ES=6.11*10.**(TETENA*TSC/(TETENB+TSC))   ! Tetens
 583 !      DESDT=( 6.1078*(2500.-2.4*TSC)/ &        ! Tetens
 584 !            (0.46151*(TSC+273.15)**2.) )*10.**(7.5*TSC/(237.3+TSC))
 585 !      ES=6.11*EXP((2.5*10.**6./461.51)*(TRP-273.15)/(273.15*TRP) ) ! Clausius-Clapeyron
 586 !      DESDT=(2.5*10.**6./461.51)*ES/(TRP**2.)                      ! Clausius-Clapeyron
 587 !      QS0R=0.622*ES/(PS-0.378*ES)
 588 !      DQS0RDTR = DESDT*0.622*PS/((PS-0.378*ES)**2.)
 589 !      DQS0RDTR = 17.269*(273.15-35.86)/((TRP-35.86)**2.)*QS0R
 590
 591 !    TRP=350.
 592
 593      DO ITERATION=1,20
 594
 595        ES=6.11*EXP( (2.5*10.**6./461.51)*(TRP-273.15)/(273.15*TRP) )
 596        DESDT=(2.5*10.**6./461.51)*ES/(TRP**2.)
 597        QS0R=0.622*ES/(PS-0.378*ES)
 598        DQS0RDTR = DESDT*0.622*PS/((PS-0.378*ES)**2.)
 599
 600        RR=EPSR*(RX-SIG*(TRP**4.)/60.)
 601        HR=RHO*CP*CHR*UA*(TRP-TA)*100.
 602        ELER=RHO*EL*CHR*UA*BETR*(QS0R-QA)*100.
 603        G0R=AKSR*(TRP-TRL(1))/(DZR(1)/2.)
 604
 605        F = SR + RR - HR - ELER - G0R
 606
 607        DRRDTR = (-4.*EPSR*SIG*TRP**3.)/60.
 608        DHRDTR = RHO*CP*CHR*UA*100.
 609        DELERDTR = RHO*EL*CHR*UA*BETR*DQS0RDTR*100.
 610        DG0RDTR =  2.*AKSR/DZR(1)
 611
 612        DFDT = DRRDTR - DHRDTR - DELERDTR - DG0RDTR
 613        DTR = F/DFDT
 614
 615        TR = TRP - DTR
 616        TRP = TR
 617
 618        IF( ABS(F) < 0.000001 .AND. ABS(DTR) < 0.000001 ) EXIT
 619
 620      END DO
 621
 622 ! multi-layer heat equation model
 623
 624      CALL multi_layer(num_roof_layers,BOUNDR,G0R,CAPR,AKSR,TRL,DZR,DELT,TRLEND)
 625
 626    ELSE
 627
 628      ES=6.11*EXP( (2.5*10.**6./461.51)*(TRP-273.15)/(273.15*TRP) )
 629      QS0R=0.622*ES/(PS-0.378*ES)
 630
 631      RR=EPSR*(RX-SIG*(TRP**4.)/60.)
 632      HR=RHO*CP*CHR*UA*(TRP-TA)*100.
 633      ELER=RHO*EL*CHR*UA*BETR*(QS0R-QA)*100.
 634      G0R=SR+RR-HR-ELER
 635
 636      CALL force_restore(CAPR,AKSR,DELT,SR,RR,HR,ELER,TRLEND,TRP,TR)
 637
 638      TRP=TR
 639
 640    END IF
 641
 642    FLXTHR=HR/RHO/CP/100.
 643    FLXHUMR=ELER/RHO/EL/100.
 644
 645 !-------------------------------------------------------------------------------
 646 !  Wall and Road
 647 !-------------------------------------------------------------------------------
 648
 649 !-------------------------------------------------------------------------------
 650 !  CHC, CHB, CDB, BETB, CHG, CDG, BETG
 651 !-------------------------------------------------------------------------------
 652
 653    Z=ZA-ZDC
 654    BHC=LOG(Z0C/Z0HC)/0.4
 655    RIBC=(9.8*2./(TA+TCP))*(TA-TCP)*(Z+Z0C)/(UA*UA)
 656
 657    CALL mos(XXXC,ALPHAC,CDC,BHC,RIBC,Z,Z0C,UA,TA,TCP,RHO)
 658
 659    IF (CH_SCHEME == 1) THEN
 660
 661      Z=ZDC
 662      BHB=LOG(Z0B/Z0HB)/0.4
 663      BHG=LOG(Z0G/Z0HG)/0.4
 664      RIBB=(9.8*2./(TCP+TBP))*(TCP-TBP)*(Z+Z0B)/(UC*UC)
 665      RIBG=(9.8*2./(TCP+TGP))*(TCP-TGP)*(Z+Z0G)/(UC*UC)
 666
 667      CALL mos(XXXB,ALPHAB,CDB,BHB,RIBB,Z,Z0B,UC,TCP,TBP,RHO)
 668      CALL mos(XXXG,ALPHAG,CDG,BHG,RIBG,Z,Z0G,UC,TCP,TGP,RHO)
 669
 670    ELSE
 671
 672      ALPHAB=RHO*CP*(6.15+4.18*UC)/1200.
 673      IF(UC > 5.) ALPHAB=RHO*CP*(7.51*UC**0.78)/1200.
 674      ALPHAG=RHO*CP*(6.15+4.18*UC)/1200.
 675      IF(UC > 5.) ALPHAG=RHO*CP*(7.51*UC**0.78)/1200.
 676
 677    END IF
 678
 679    CHC=ALPHAC/RHO/CP/UA
 680    CHB=ALPHAB/RHO/CP/UC
 681    CHG=ALPHAG/RHO/CP/UC
 682
 683    BETB=0.0
 684    IF(RAIN > 1.) BETG=0.7
 685
 686    IF (TS_SCHEME == 1) THEN
 687
 688 !-------------------------------------------------------------------------------
 689 !  TB, TG  Solving Non-Linear Simultaneous Equation by Newton-Rapson
 690 !  TBL,TGL Solving Heat Equation by Tri Diagonal Matrix Algorithm
 691 !-------------------------------------------------------------------------------
 692
 693 !    TBP=350.
 694 !    TGP=350.
 695
 696      DO ITERATION=1,20
 697
 698        ES=6.11*EXP( (2.5*10.**6./461.51)*(TBP-273.15)/(273.15*TBP) )
 699        DESDT=(2.5*10.**6./461.51)*ES/(TBP**2.)
 700        QS0B=0.622*ES/(PS-0.378*ES)
 701        DQS0BDTB=DESDT*0.622*PS/((PS-0.378*ES)**2.)
 702
 703        ES=6.11*EXP( (2.5*10.**6./461.51)*(TGP-273.15)/(273.15*TGP) )
 704        DESDT=(2.5*10.**6./461.51)*ES/(TGP**2.)
 705        QS0G=0.622*ES/(PS-0.378*ES)
 706        DQS0GDTG=DESDT*0.22*PS/((PS-0.378*ES)**2.)
 707
 708        RG1=EPSG*( RX*VFGS          &
 709        +EPSB*VFGW*SIG*TBP**4./60.  &
 710        -SIG*TGP**4./60. )
 711
 712        RB1=EPSB*( RX*VFWS         &
 713        +EPSG*VFWG*SIG*TGP**4./60. &
 714        +EPSB*VFWW*SIG*TBP**4./60. &
 715        -SIG*TBP**4./60. )
 716
 717        RG2=EPSG*( (1.-EPSB)*(1.-SVF)*VFWS*RX                  &
 718        +(1.-EPSB)*(1.-SVF)*VFWG*EPSG*SIG*TGP**4./60.          &
 719        +EPSB*(1.-EPSB)*(1.-SVF)*(1.-2.*VFWS)*SIG*TBP**4./60. )
 720
 721        RB2=EPSB*( (1.-EPSG)*VFWG*VFGS*RX                          &
 722        +(1.-EPSG)*EPSB*VFGW*VFWG*SIG*(TBP**4.)/60.                &
 723        +(1.-EPSB)*VFWS*(1.-2.*VFWS)*RX                            &
 724        +(1.-EPSB)*VFWG*(1.-2.*VFWS)*EPSG*SIG*EPSG*TGP**4./60.     &
 725        +EPSB*(1.-EPSB)*(1.-2.*VFWS)*(1.-2.*VFWS)*SIG*TBP**4./60. )
 726
 727        RG=RG1+RG2
 728        RB=RB1+RB2
 729
 730        DRBDTB1=EPSB*(4.*EPSB*SIG*TB**3.*VFWW-4.*SIG*TB**3.)/60.
 731        DRBDTG1=EPSB*(4.*EPSG*SIG*TG**3.*VFWG)/60.
 732        DRBDTB2=EPSB*(4.*(1.-EPSG)*EPSB*SIG*TB**3.*VFGW*VFWG &
 733                +4.*EPSB*(1.-EPSB)*SIG*TB**3.*VFWW*VFWW)/60.
 734        DRBDTG2=EPSB*(4.*(1.-EPSB)*EPSG*SIG*TG**3.*VFWG*VFWW)/60.
 735
 736        DRGDTB1=EPSG*(4.*EPSB*SIG*TB**3.*VFGW)/60.
 737        DRGDTG1=EPSG*(-4.*SIG*TG**3.)/60.
 738        DRGDTB2=EPSG*(4.*EPSB*(1.-EPSB)*SIG*TB**3.*VFWW*VFGW)/60.
 739        DRGDTG2=EPSG*(4.*(1.-EPSB)*EPSG*SIG*TG**3.*VFWG*VFGW)/60.
 740
 741        DRBDTB=DRBDTB1+DRBDTB2
 742        DRBDTG=DRBDTG1+DRBDTG2
 743        DRGDTB=DRGDTB1+DRGDTB2
 744        DRGDTG=DRGDTG1+DRGDTG2
 745
 746        HB=RHO*CP*CHB*UC*(TBP-TCP)*100.
 747        HG=RHO*CP*CHG*UC*(TGP-TCP)*100.
 748
 749        DTCDTB=W*ALPHAB/(RW*ALPHAC+RW*ALPHAG+W*ALPHAB)
 750        DTCDTG=RW*ALPHAG/(RW*ALPHAC+RW*ALPHAG+W*ALPHAB)
 751
 752        DHBDTB=RHO*CP*CHB*UC*(1.-DTCDTB)*100.
 753        DHBDTG=RHO*CP*CHB*UC*(0.-DTCDTG)*100.
 754        DHGDTG=RHO*CP*CHG*UC*(1.-DTCDTG)*100.
 755        DHGDTB=RHO*CP*CHG*UC*(0.-DTCDTB)*100.
 756
 757        ELEB=RHO*EL*CHB*UC*BETB*(QS0B-QCP)*100.
 758        ELEG=RHO*EL*CHG*UC*BETG*(QS0G-QCP)*100.
 759
 760        DQCDTB=W*ALPHAB*BETB*DQS0BDTB/(RW*ALPHAC+RW*ALPHAG*BETG+W*ALPHAB*BETB)
 761        DQCDTG=RW*ALPHAG*BETG*DQS0GDTG/(RW*ALPHAC+RW*ALPHAG*BETG+W*ALPHAB*BETB)
 762
 763        DELEBDTB=RHO*EL*CHB*UC*BETB*(DQS0BDTB-DQCDTB)*100.
 764        DELEBDTG=RHO*EL*CHB*UC*BETB*(0.-DQCDTG)*100.
 765        DELEGDTG=RHO*EL*CHG*UC*BETG*(DQS0GDTG-DQCDTG)*100.
 766        DELEGDTB=RHO*EL*CHG*UC*BETG*(0.-DQCDTB)*100.
 767
 768        G0B=AKSB*(TBP-TBL(1))/(DZB(1)/2.)
 769        G0G=AKSG*(TGP-TGL(1))/(DZG(1)/2.)
 770
 771        DG0BDTB=2.*AKSB/DZB(1)
 772        DG0BDTG=0.
 773        DG0GDTG=2.*AKSG/DZG(1)
 774        DG0GDTB=0.
 775
 776        F = SB + RB - HB - ELEB - G0B
 777        FX = DRBDTB - DHBDTB - DELEBDTB - DG0BDTB
 778        FY = DRBDTG - DHBDTG - DELEBDTG - DG0BDTG
 779
 780        GF = SG + RG - HG - ELEG - G0G
 781        GX = DRGDTB - DHGDTB - DELEGDTB - DG0GDTB
 782        GY = DRGDTG - DHGDTG - DELEGDTG - DG0GDTG
 783
 784        DTB =  (GF*FY-F*GY)/(FX*GY-GX*FY)
 785        DTG = -(GF+GX*DTB)/GY
 786
 787        TB = TBP + DTB
 788        TG = TGP + DTG
 789
 790        TBP = TB
 791        TGP = TG
 792
 793        TC1=RW*ALPHAC+RW*ALPHAG+W*ALPHAB
 794        TC2=RW*ALPHAC*TA+RW*ALPHAG*TGP+W*ALPHAB*TBP
 795        TC=TC2/TC1
 796
 797        QC1=RW*ALPHAC+RW*ALPHAG*BETG+W*ALPHAB*BETB
 798        QC2=RW*ALPHAC*QA+RW*ALPHAG*BETG*QS0G+W*ALPHAB*BETB*QS0B
 799        QC=QC2/QC1
 800
 801        DTC=TCP - TC
 802        TCP=TC
 803        QCP=QC
 804
 805        IF( ABS(F) < 0.000001 .AND. ABS(DTB) < 0.000001 &
 806         .AND. ABS(GF) < 0.000001 .AND. ABS(DTG) < 0.000001 &
 807         .AND. ABS(DTC) < 0.000001) EXIT
 808
 809      END DO
 810
 811      CALL multi_layer(num_wall_layers,BOUNDB,G0B,CAPB,AKSB,TBL,DZB,DELT,TBLEND)
 812
 813      CALL multi_layer(num_road_layers,BOUNDG,G0G,CAPG,AKSG,TGL,DZG,DELT,TGLEND)
 814
 815    ELSE
 816
 817 !-------------------------------------------------------------------------------
 818 !  TB, TG  by Force-Restore Method
 819 !-------------------------------------------------------------------------------
 820
 821        ES=6.11*EXP((2.5*10.**6./461.51)*(TBP-273.15)/(273.15*TBP) )
 822        QS0B=0.622*ES/(PS-0.378*ES)
 823
 824        ES=6.11*EXP((2.5*10.**6./461.51)*(TGP-273.15)/(273.15*TGP) )
 825        QS0G=0.622*ES/(PS-0.378*ES)
 826
 827        RG1=EPSG*( RX*VFGS             &
 828        +EPSB*VFGW*SIG*TBP**4./60.     &
 829        -SIG*TGP**4./60. )
 830
 831        RB1=EPSB*( RX*VFWS &
 832        +EPSG*VFWG*SIG*TGP**4./60. &
 833        +EPSB*VFWW*SIG*TBP**4./60. &
 834        -SIG*TBP**4./60. )
 835
 836        RG2=EPSG*( (1.-EPSB)*(1.-SVF)*VFWS*RX                   &
 837        +(1.-EPSB)*(1.-SVF)*VFWG*EPSG*SIG*TGP**4./60.           &
 838        +EPSB*(1.-EPSB)*(1.-SVF)*(1.-2.*VFWS)*SIG*TBP**4./60. )
 839
 840        RB2=EPSB*( (1.-EPSG)*VFWG*VFGS*RX                          &
 841        +(1.-EPSG)*EPSB*VFGW*VFWG*SIG*(TBP**4.)/60.                &
 842        +(1.-EPSB)*VFWS*(1.-2.*VFWS)*RX                            &
 843        +(1.-EPSB)*VFWG*(1.-2.*VFWS)*EPSG*SIG*EPSG*TGP**4./60.     &
 844        +EPSB*(1.-EPSB)*(1.-2.*VFWS)*(1.-2.*VFWS)*SIG*TBP**4./60. )
 845
 846        RG=RG1+RG2
 847        RB=RB1+RB2
 848
 849        HB=RHO*CP*CHB*UC*(TBP-TCP)*100.
 850        ELEB=RHO*EL*CHB*UC*BETB*(QS0B-QCP)*100.
 851        G0B=SB+RB-HB-ELEB
 852
 853        HG=RHO*CP*CHG*UC*(TGP-TCP)*100.
 854        ELEG=RHO*EL*CHG*UC*BETG*(QS0G-QCP)*100.
 855        G0G=SG+RG-HG-ELEG
 856
 857        CALL force_restore(CAPB,AKSB,DELT,SB,RB,HB,ELEB,TBLEND,TBP,TB)
 858        CALL force_restore(CAPG,AKSG,DELT,SG,RG,HG,ELEG,TGLEND,TGP,TG)
 859
 860        TBP=TB
 861        TGP=TG
 862
 863        TC1=RW*ALPHAC+RW*ALPHAG+W*ALPHAB
 864        TC2=RW*ALPHAC*TA+RW*ALPHAG*TGP+W*ALPHAB*TBP
 865        TC=TC2/TC1
 866
 867        QC1=RW*ALPHAC+RW*ALPHAG*BETG+W*ALPHAB*BETB
 868        QC2=RW*ALPHAC*QA+RW*ALPHAG*BETG*QS0G+W*ALPHAB*BETB*QS0B
 869        QC=QC2/QC1
 870
 871        TCP=TC
 872        QCP=QC
 873
 874    END IF
 875
 876    FLXTHB=HB/RHO/CP/100.
 877    FLXHUMB=ELEB/RHO/EL/100.
 878    FLXTHG=HG/RHO/CP/100.
 879    FLXHUMG=ELEG/RHO/EL/100.
 880
 881 !-------------------------------------------------------------------------------
 882 ! Total Fulxes from Urban Canopy
 883 !-------------------------------------------------------------------------------
 884
 885    FLXUV  = ( R*CDR + RW*CDC )*UA*UA
 886 ! Miao, 2007/01/17, cal. ah
 887    if(ahoption==1) then
 888      FLXTH  = ( R*FLXTHR  + W*FLXTHB  + RW*FLXTHG ) + AH/RHOO/CPP
 889    else
 890      FLXTH  = ( R*FLXTHR  + W*FLXTHB  + RW*FLXTHG )
 891    endif
 892    FLXHUM = ( R*FLXHUMR + W*FLXHUMB + RW*FLXHUMG )
 893    FLXG =   ( R*G0R + W*G0B + RW*G0G )
 894    LNET =     R*RR + W*RB + RW*RG
 895
 896 !----------------------------------------------------------------------------
 897 ! Convert Unit: FLUXES and u* T* q*  --> WRF
 898 !----------------------------------------------------------------------------
 899
 900    SH    = FLXTH  * RHOO * CPP    ! Sensible heat flux          [W/m/m]
 901    LH    = FLXHUM * RHOO * ELL    ! Latent heat flux            [W/m/m]
 902    LH_KINEMATIC = FLXHUM * RHOO   ! Latent heat, Kinematic      [kg/m/m/s]
 903    LW    = LLG - (LNET*697.7*60.) ! Upward longwave radiation   [W/m/m]
 904    SW    = SSG - (SNET*697.7*60.) ! Upward shortwave radiation  [W/m/m]
 905    ALB   = 0.
 906    IF( ABS(SSG) > 0.0001) ALB = SW/SSG ! Effective albedo [-]
 907    G = -FLXG*697.7*60.            ! [W/m/m]
 908    RN = (SNET+LNET)*697.7*60.     ! Net radiation [W/m/m]
 909
 910    UST = SQRT(FLXUV)              ! u* [m/s]
 911    TST = -FLXTH/UST               ! T* [K]
 912    QST = -FLXHUM/UST              ! q* [-]
 913
 914 !------------------------------------------------------
 915 !  diagnostic GRID AVERAGED  PSIM  PSIH  TS QS --> WRF
 916 !------------------------------------------------------
 917
 918    Z0 = Z0C
 919    Z0H = Z0HC
 920    Z = ZA - ZDC
 921
 922    XXX = 0.4*9.81*Z*TST/TA/UST/UST
 923
 924    IF ( XXX >= 1. ) XXX = 1.
 925    IF ( XXX <= -5. ) XXX = -5.
 926
 927    IF ( XXX > 0 ) THEN
 928      PSIM = -5. * XXX
 929      PSIH = -5. * XXX
 930    ELSE
 931      X = (1.-16.*XXX)**0.25
 932      PSIM = 2.*ALOG((1.+X)/2.) + ALOG((1.+X*X)/2.) - 2.*ATAN(X) + PI/2.
 933      PSIH = 2.*ALOG((1.+X*X)/2.)
 934    END IF
 935
 936    GZ1OZ0 = ALOG(Z/Z0)
 937    CD = 0.4**2./(ALOG(Z/Z0)-PSIM)**2.
 938 !
 939 !m   CH = 0.4**2./(ALOG(Z/Z0)-PSIM)/(ALOG(Z/Z0H)-PSIH)
 940 !m   CHS = 0.4*UST/(ALOG(Z/Z0H)-PSIH)
 941 !m   TS = TA + FLXTH/CH/UA    ! surface potential temp (flux temp)
 942 !m   QS = QA + FLXHUM/CH/UA   ! surface humidity
 943 !
 944    TS = TA + FLXTH/CHS    ! surface potential temp (flux temp)
 945    QS = QA + FLXHUM/CHS   ! surface humidity
 946
 947 !-------------------------------------------------------
 948 !  diagnostic  GRID AVERAGED  U10  V10  TH2  Q2 --> WRF
 949 !-------------------------------------------------------
 950
 951    XXX2 = (2./Z)*XXX
 952    IF ( XXX2 >= 1. ) XXX2 = 1.
 953    IF ( XXX2 <= -5. ) XXX2 = -5.
 954
 955    IF ( XXX2 > 0 ) THEN
 956       PSIM2 = -5. * XXX2
 957       PSIH2 = -5. * XXX2
 958    ELSE
 959       X = (1.-16.*XXX2)**0.25
 960       PSIM2 = 2.*ALOG((1.+X)/2.) + ALOG((1.+X*X)/2.) - 2.*ATAN(X) + 2.*ATAN(1.)
 961       PSIH2 = 2.*ALOG((1.+X*X)/2.)
 962    END IF
 963 !
 964 !m   CHS2 = 0.4*UST/(ALOG(2./Z0H)-PSIH2)
 965 !
 966
 967    XXX10 = (10./Z)*XXX
 968    IF ( XXX10 >= 1. ) XXX10 = 1.
 969    IF ( XXX10 <= -5. ) XXX10 = -5.
 970
 971    IF ( XXX10 > 0 ) THEN
 972       PSIM10 = -5. * XXX10
 973       PSIH10 = -5. * XXX10
 974    ELSE
 975       X = (1.-16.*XXX10)**0.25
 976       PSIM10 = 2.*ALOG((1.+X)/2.) + ALOG((1.+X*X)/2.) - 2.*ATAN(X) + 2.*ATAN(1.)
 977       PSIH10 = 2.*ALOG((1.+X*X)/2.)
 978    END IF
 979
 980    PSIX = ALOG(Z/Z0) - PSIM
 981    PSIT = ALOG(Z/Z0H) - PSIH
 982
 983    PSIX2 = ALOG(2./Z0) - PSIM2
 984    PSIT2 = ALOG(2./Z0H) - PSIH2
 985
 986    PSIX10 = ALOG(10./Z0) - PSIM10
 987    PSIT10 = ALOG(10./Z0H) - PSIH10
 988
 989    U10 = U1 * (PSIX10/PSIX)       ! u at 10 m [m/s]
 990    V10 = V1 * (PSIX10/PSIX)       ! v at 10 m [m/s]
 991
 992 !   TH2 = TS + (TA-TS)*(PSIT2/PSIT)   ! potential temp at 2 m [K]
 993 !  TH2 = TS + (TA-TS)*(PSIT2/PSIT)   ! Fei: this seems to be temp (not potential) at 2 m [K]
 994 !Fei: consistant with M-O theory
 995    TH2 = TS + (TA-TS) *(CHS/CHS2)
 996
 997    Q2 = QS + (QA-QS)*(PSIT2/PSIT)    ! humidity at 2 m       [-]
 998
 999 !   TS = (LW/SIG_SI/0.88)**0.25       ! Radiative temperature [K]
1000
1001    RETURN
1002
1003    END SUBROUTINE urban
1004 !===============================================================================
1005 !
1006 !  mos
1007 !
1008 !===============================================================================
1009    SUBROUTINE mos(XXX,ALPHA,CD,B1,RIB,Z,Z0,UA,TA,TSF,RHO)
1010
1011 !  XXX:   z/L (requires iteration by Newton-Rapson method)
1012 !  B1:    Stanton number
1013 !  PSIM:  = PSIX of LSM
1014 !  PSIH:  = PSIT of LSM
1015
1016    IMPLICIT NONE
1017
1018    REAL, PARAMETER     :: CP=0.24
1019    REAL, INTENT(IN)    :: B1, Z, Z0, UA, TA, TSF, RHO
1020    REAL, INTENT(OUT)   :: ALPHA, CD
1021    REAL, INTENT(INOUT) :: XXX, RIB
1022    REAL                :: XXX0, X, X0, FAIH, DPSIM, DPSIH
1023    REAL                :: F, DF, XXXP, US, TS, AL, XKB, DD, PSIM, PSIH
1024    INTEGER             :: NEWT
1025    INTEGER, PARAMETER  :: NEWT_END=10
1026
1027    IF(RIB <= -15.) RIB=-15.
1028
1029    IF(RIB < 0.) THEN
1030
1031       DO NEWT=1,NEWT_END
1032
1033         IF(XXX >= 0.) XXX=-1.E-3
1034
1035         XXX0=XXX*Z0/(Z+Z0)
1036
1037         X=(1.-16.*XXX)**0.25
1038         X0=(1.-16.*XXX0)**0.25
1039
1040         PSIM=ALOG((Z+Z0)/Z0) &
1041             -ALOG((X+1.)**2.*(X**2.+1.)) &
1042             +2.*ATAN(X) &
1043             +ALOG((X+1.)**2.*(X0**2.+1.)) &
1044             -2.*ATAN(X0)
1045         FAIH=1./SQRT(1.-16.*XXX)
1046         PSIH=ALOG((Z+Z0)/Z0)+0.4*B1 &
1047             -2.*ALOG(SQRT(1.-16.*XXX)+1.) &
1048             +2.*ALOG(SQRT(1.-16.*XXX0)+1.)
1049
1050         DPSIM=(1.-16.*XXX)**(-0.25)/XXX &
1051              -(1.-16.*XXX0)**(-0.25)/XXX
1052         DPSIH=1./SQRT(1.-16.*XXX)/XXX &
1053              -1./SQRT(1.-16.*XXX0)/XXX
1054
1055         F=RIB*PSIM**2./PSIH-XXX
1056
1057         DF=RIB*(2.*DPSIM*PSIM*PSIH-DPSIH*PSIM**2.) &
1058           /PSIH**2.-1.
1059
1060         XXXP=XXX
1061         XXX=XXXP-F/DF
1062         IF(XXX <= -10.) XXX=-10.
1063
1064       END DO
1065
1066    ELSE IF(RIB >= 0.142857) THEN
1067
1068       XXX=0.714
1069       PSIM=ALOG((Z+Z0)/Z0)+7.*XXX
1070       PSIH=PSIM+0.4*B1
1071
1072    ELSE
1073
1074       AL=ALOG((Z+Z0)/Z0)
1075       XKB=0.4*B1
1076       DD=-4.*RIB*7.*XKB*AL+(AL+XKB)**2.
1077       IF(DD <= 0.) DD=0.
1078       XXX=(AL+XKB-2.*RIB*7.*AL-SQRT(DD))/(2.*(RIB*7.**2-7.))
1079       PSIM=ALOG((Z+Z0)/Z0)+7.*MIN(XXX,0.714)
1080       PSIH=PSIM+0.4*B1
1081
1082    END IF
1083
1084    US=0.4*UA/PSIM             ! u*
1085    IF(US <= 0.01) US=0.01
1086    TS=0.4*(TA-TSF)/PSIH       ! T*
1087
1088    CD=US*US/UA**2.            ! CD
1089    ALPHA=RHO*CP*0.4*US/PSIH   ! RHO*CP*CH*U
1090
1091    RETURN
1092    END SUBROUTINE mos
1093 !===============================================================================
1094 !
1095 !  louis79
1096 !
1097 !===============================================================================
1098    SUBROUTINE louis79(ALPHA,CD,RIB,Z,Z0,UA,RHO)
1099
1100    IMPLICIT NONE
1101
1102    REAL, PARAMETER     :: CP=0.24
1103    REAL, INTENT(IN)    :: Z, Z0, UA, RHO
1104    REAL, INTENT(OUT)   :: ALPHA, CD
1105    REAL, INTENT(INOUT) :: RIB
1106    REAL                :: A2, XX, CH, CMB, CHB
1107
1108    A2=(0.4/ALOG(Z/Z0))**2.
1109
1110    IF(RIB <= -15.) RIB=-15.
1111
1112    IF(RIB >= 0.0) THEN
1113       IF(RIB >= 0.142857) THEN
1114          XX=0.714
1115       ELSE
1116          XX=RIB*LOG(Z/Z0)/(1.-7.*RIB)
1117       END IF
1118       CH=0.16/0.74/(LOG(Z/Z0)+7.*MIN(XX,0.714))**2.
1119       CD=0.16/(LOG(Z/Z0)+7.*MIN(XX,0.714))**2.
1120    ELSE
1121       CMB=7.4*A2*9.4*SQRT(Z/Z0)
1122       CHB=5.3*A2*9.4*SQRT(Z/Z0)
1123       CH=A2/0.74*(1.-9.4*RIB/(1.+CHB*SQRT(-RIB)))
1124       CD=A2*(1.-9.4*RIB/(1.+CHB*SQRT(-RIB)))
1125    END IF
1126
1127    ALPHA=RHO*CP*CH*UA
1128
1129    RETURN
1130    END SUBROUTINE louis79
1131 !===============================================================================
1132 !
1133 !  louis82
1134 !
1135 !===============================================================================
1136    SUBROUTINE louis82(ALPHA,CD,RIB,Z,Z0,UA,RHO)
1137
1138    IMPLICIT NONE
1139
1140    REAL, PARAMETER     :: CP=0.24
1141    REAL, INTENT(IN)    :: Z, Z0, UA, RHO
1142    REAL, INTENT(OUT)   :: ALPHA, CD
1143    REAL, INTENT(INOUT) :: RIB
1144    REAL                :: A2, FM, FH, CH, CHH
1145
1146    A2=(0.4/ALOG(Z/Z0))**2.
1147
1148    IF(RIB <= -15.) RIB=-15.
1149
1150    IF(RIB >= 0.0) THEN
1151       FM=1./((1.+(2.*5.*RIB)/SQRT(1.+5.*RIB)))
1152       FH=1./(1.+(3.*5.*RIB)*SQRT(1.+5.*RIB))
1153       CH=A2*FH
1154       CD=A2*FM
1155    ELSE
1156       CHH=5.*3.*5.*A2*SQRT(Z/Z0)
1157       FM=1.-(2.*5.*RIB)/(1.+3.*5.*5.*A2*SQRT(Z/Z0+1.)*(-RIB))
1158       FH=1.-(3.*5.*RIB)/(1.+CHH*SQRT(-RIB))
1159       CH=A2*FH
1160       CD=A2*FM
1161    END IF
1162
1163    ALPHA=RHO*CP*CH*UA
1164
1165    RETURN
1166    END SUBROUTINE louis82
1167 !===============================================================================
1168 !
1169 !  multi_layer
1170 !
1171 !===============================================================================
1172    SUBROUTINE multi_layer(KM,BOUND,G0,CAP,AKS,TSL,DZ,DELT,TSLEND)
1173
1174    IMPLICIT NONE
1175
1176    REAL, INTENT(IN)                   :: G0, CAP, AKS, DELT,TSLEND
1177
1178    INTEGER, INTENT(IN)                :: KM, BOUND
1179
1180    REAL, DIMENSION(KM), INTENT(IN)    :: DZ
1181
1182    REAL, DIMENSION(KM), INTENT(INOUT) :: TSL
1183
1184    REAL, DIMENSION(KM)                :: A, B, C, D, X, P, Q
1185
1186    REAL                               :: DZEND
1187
1188    INTEGER                            :: K
1189
1190    DZEND=DZ(KM)
1191
1192    A(1) = 0.0
1193
1194    B(1) = CAP*DZ(1)/DELT &
1195           +2.*AKS/(DZ(1)+DZ(2))
1196    C(1) = -2.*AKS/(DZ(1)+DZ(2))
1197    D(1) = CAP*DZ(1)/DELT*TSL(1) + G0
1198
1199    DO K=2,KM-1
1200       A(K) = -2.*AKS/(DZ(K-1)+DZ(K))
1201       B(K) = CAP*DZ(K)/DELT + 2.*AKS/(DZ(K-1)+DZ(K)) + 2.*AKS/(DZ(K)+DZ(K+1))
1202       C(K) = -2.*AKS/(DZ(K)+DZ(K+1))
1203       D(K) = CAP*DZ(K)/DELT*TSL(K)
1204    END DO
1205
1206    IF(BOUND == 1) THEN                  ! Flux=0
1207       A(KM) = -2.*AKS/(DZ(KM-1)+DZ(KM))
1208       B(KM) = CAP*DZ(KM)/DELT + 2.*AKS/(DZ(KM-1)+DZ(KM))
1209       C(KM) = 0.0
1210       D(KM) = CAP*DZ(KM)/DELT*TSL(KM)
1211    ELSE                                 ! T=constant
1212       A(KM) = -2.*AKS/(DZ(KM-1)+DZ(KM))
1213       B(KM) = CAP*DZ(KM)/DELT + 2.*AKS/(DZ(KM-1)+DZ(KM)) + 2.*AKS/(DZ(KM)+DZEND)
1214       C(KM) = 0.0
1215       D(KM) = CAP*DZ(KM)/DELT*TSL(KM) + 2.*AKS*TSLEND/(DZ(KM)+DZEND)
1216    END IF
1217
1218    P(1) = -C(1)/B(1)
1219    Q(1) =  D(1)/B(1)
1220
1221    DO K=2,KM
1222       P(K) = -C(K)/(A(K)*P(K-1)+B(K))
1223       Q(K) = (-A(K)*Q(K-1)+D(K))/(A(K)*P(K-1)+B(K))
1224    END DO
1225
1226    X(KM) = Q(KM)
1227
1228    DO K=KM-1,1,-1
1229       X(K) = P(K)*X(K+1)+Q(K)
1230    END DO
1231
1232    DO K=1,KM
1233       TSL(K) = X(K)
1234    END DO
1235
1236    RETURN
1237    END SUBROUTINE multi_layer
1238 !===============================================================================
1239 !
1240 !  subroutine read_param
1241 !
1242 !===============================================================================
1243    SUBROUTINE read_param(UTYPE,                                        & ! in
1244                          ZR,Z0C,Z0HC,ZDC,SVF,R,RW,HGT,CDS,AS,AH,       & ! out
1245                          CAPR,CAPB,CAPG,AKSR,AKSB,AKSG,ALBR,ALBB,ALBG, & ! out
1246                          EPSR,EPSB,EPSG,Z0R,Z0B,Z0G,Z0HR,Z0HB,Z0HG,    & ! out
1247                          BETR,BETB,BETG,TRLEND,TBLEND,TGLEND,          & ! out
1248                          BOUNDR,BOUNDB,BOUNDG,CH_SCHEME,TS_SCHEME)       ! out
1249
1250    INTEGER, INTENT(IN)  :: UTYPE
1251
1252    REAL, INTENT(OUT)    :: ZR,Z0C,Z0HC,ZDC,SVF,R,RW,HGT,CDS,AS,AH,       &
1253                            CAPR,CAPB,CAPG,AKSR,AKSB,AKSG,ALBR,ALBB,ALBG, &
1254                            EPSR,EPSB,EPSG,Z0R,Z0B,Z0G,Z0HR,Z0HB,Z0HG,    &
1255                            BETR,BETB,BETG,TRLEND,TBLEND,TGLEND
1256
1257    INTEGER, INTENT(OUT) :: BOUNDR,BOUNDB,BOUNDG,CH_SCHEME,TS_SCHEME
1258
1259    ZR =     ZR_TBL(UTYPE)
1260    Z0C=     Z0C_TBL(UTYPE)
1261    Z0HC=    Z0HC_TBL(UTYPE)
1262    ZDC=     ZDC_TBL(UTYPE)
1263    SVF=     SVF_TBL(UTYPE)
1264    R=       R_TBL(UTYPE)
1265    RW=      RW_TBL(UTYPE)
1266    HGT=     HGT_TBL(UTYPE)
1267    CDS=     CDS_TBL(UTYPE)
1268    AS=      AS_TBL(UTYPE)
1269    AH=      AH_TBL(UTYPE)
1270    BETR=    BETR_TBL(UTYPE)
1271    BETB=    BETB_TBL(UTYPE)
1272    BETG=    BETG_TBL(UTYPE)
1273
1274 !m   FRC_URB= FRC_URB_TBL(UTYPE)
1275
1276    CAPR=    CAPR_DATA
1277    CAPB=    CAPB_DATA
1278    CAPG=    CAPG_DATA
1279    AKSR=    AKSR_DATA
1280    AKSB=    AKSB_DATA
1281    AKSG=    AKSG_DATA
1282    ALBR=    ALBR_DATA
1283    ALBB=    ALBB_DATA
1284    ALBG=    ALBG_DATA
1285    EPSR=    EPSR_DATA
1286    EPSB=    EPSB_DATA
1287    EPSG=    EPSG_DATA
1288    Z0R=     Z0R_DATA
1289    Z0B=     Z0B_DATA
1290    Z0G=     Z0G_DATA
1291    Z0HR=    Z0HR_DATA
1292    Z0HB=    Z0HB_DATA
1293    Z0HG=    Z0HG_DATA
1294    TRLEND=  TRLEND_DATA
1295    TBLEND=  TBLEND_DATA
1296    TGLEND=  TGLEND_DATA
1297    BOUNDR=  BOUNDR_DATA
1298    BOUNDB=  BOUNDB_DATA
1299    BOUNDG=  BOUNDG_DATA
1300    CH_SCHEME = CH_SCHEME_DATA
1301    TS_SCHEME = TS_SCHEME_DATA
1302
1303    RETURN
1304    END SUBROUTINE read_param
1305 !===============================================================================
1306 !
1307 ! subroutine urban_param_init: Read parameters from urban_param.tbl
1308 !
1309 !===============================================================================
1310    SUBROUTINE urban_param_init(DZR,DZB,DZG,num_soil_layers &
1311                                )
1312 !                              num_roof_layers,num_wall_layers,num_road_layers)
1313
1314    IMPLICIT NONE
1315
1316    INTEGER, INTENT(IN) :: num_soil_layers
1317
1318 !   REAL, DIMENSION(1:num_roof_layers), INTENT(INOUT) :: DZR
1319 !   REAL, DIMENSION(1:num_wall_layers), INTENT(INOUT) :: DZB
1320 !   REAL, DIMENSION(1:num_road_layers), INTENT(INOUT) :: DZG
1321    REAL, DIMENSION(1:num_soil_layers), INTENT(INOUT) :: DZR
1322    REAL, DIMENSION(1:num_soil_layers), INTENT(INOUT) :: DZB
1323    REAL, DIMENSION(1:num_soil_layers), INTENT(INOUT) :: DZG
1324
1325    INTEGER :: INDEX, LC, K
1326    INTEGER :: IOSTATUS, ALLOCATE_STATUS
1327    INTEGER :: num_roof_layers
1328    INTEGER :: num_wall_layers
1329    INTEGER :: num_road_layers
1330    INTEGER :: dummy
1331    REAL    :: DHGT, HGT, VFWS, VFGS
1332
1333    OPEN (UNIT=11,                &
1334          FILE='urban_param.tbl', &
1335          ACCESS='SEQUENTIAL',    &
1336          STATUS='OLD',           &
1337          ACTION='READ',          &
1338          POSITION='REWIND',      &
1339          IOSTAT=IOSTATUS)
1340
1341    IF (IOSTATUS > 0) STOP 'ERROR OPEN urban_param.tbl'
1342
1343    READ(11,*)
1344    READ(11,'(A4)') LU_DATA_TYPE
1345
1346    READ(11,*) ICATE
1347    ALLOCATE( ZR_TBL(ICATE), stat=allocate_status )
1348    if(allocate_status == 0) THEN
1349    ALLOCATE( Z0C_TBL(ICATE), stat=allocate_status )
1350     if(allocate_status /= 0) stop 'error allocate Z0C_TBL in urban_param_init'
1351    IF( .NOT. ALLOCATED( Z0HC_TBL ) ) &
1352    ALLOCATE( Z0HC_TBL(ICATE), stat=allocate_status )
1353     if(allocate_status /= 0) stop 'error allocate Z0HC_TBL in urban_param_init'
1354    IF( .NOT. ALLOCATED( ZDC_TBL ) ) &
1355    ALLOCATE( ZDC_TBL(ICATE), stat=allocate_status )
1356     if(allocate_status /= 0) stop 'error allocate ZDC_TBL in urban_param_init'
1357    IF( .NOT. ALLOCATED( SVF_TBL ) ) &
1358    ALLOCATE( SVF_TBL(ICATE), stat=allocate_status )
1359     if(allocate_status /= 0) stop 'error allocate SVF_TBL in urban_param_init'
1360    IF( .NOT. ALLOCATED( R_TBL ) ) &
1361    ALLOCATE( R_TBL(ICATE), stat=allocate_status )
1362     if(allocate_status /= 0) stop 'error allocate R_TBL in urban_param_init'
1363    IF( .NOT. ALLOCATED( RW_TBL ) ) &
1364    ALLOCATE( RW_TBL(ICATE), stat=allocate_status )
1365     if(allocate_status /= 0) stop 'error allocate RW_TBL in urban_param_init'
1366    IF( .NOT. ALLOCATED( HGT_TBL ) ) &
1367    ALLOCATE( HGT_TBL(ICATE), stat=allocate_status )
1368     if(allocate_status /= 0) stop 'error allocate HGT_TBL in urban_param_init'
1369    IF( .NOT. ALLOCATED( CDS_TBL ) ) &
1370    ALLOCATE( CDS_TBL(ICATE), stat=allocate_status )
1371     if(allocate_status /= 0) stop 'error allocate CDS_TBL in urban_param_init'
1372    IF( .NOT. ALLOCATED( AS_TBL ) ) &
1373    ALLOCATE( AS_TBL(ICATE), stat=allocate_status )
1374     if(allocate_status /= 0) stop 'error allocate AS_TBL in urban_param_init'
1375    IF( .NOT. ALLOCATED( AH_TBL ) ) &
1376    ALLOCATE( AH_TBL(ICATE), stat=allocate_status )
1377     if(allocate_status /= 0) stop 'error allocate AH_TBL in urban_param_init'
1378    IF( .NOT. ALLOCATED( BETR_TBL ) ) &
1379    ALLOCATE( BETR_TBL(ICATE), stat=allocate_status )
1380     if(allocate_status /= 0) stop 'error allocate BETR_TBL in urban_param_init'
1381    IF( .NOT. ALLOCATED( BETB_TBL ) ) &
1382    ALLOCATE( BETB_TBL(ICATE), stat=allocate_status )
1383     if(allocate_status /= 0) stop 'error allocate BETB_TBL in urban_param_init'
1384    IF( .NOT. ALLOCATED( BETG_TBL ) ) &
1385    ALLOCATE( BETG_TBL(ICATE), stat=allocate_status )
1386     if(allocate_status /= 0) stop 'error allocate BETG_TBL in urban_param_init'
1387    ALLOCATE( FRC_URB_TBL(ICATE), stat=allocate_status )
1388     if(allocate_status /= 0) stop 'error allocate FRC_URB_TBL in urban_param_init'
1389
1390 ENDIF
1391
1392    DO LC = 1, ICATE
1393       READ(11,*) INDEX,        &
1394                  ZR_TBL(LC),   &
1395                  Z0C_TBL(LC),  &
1396                  Z0HC_TBL(LC), &
1397                  ZDC_TBL(LC),  &
1398                  SVF_TBL(LC),  &
1399                  R_TBL(LC),    &
1400                  RW_TBL(LC),   &
1401                  HGT_TBL(LC),  &
1402                  CDS_TBL(LC),  &
1403                  AS_TBL(LC),   &
1404                  AH_TBL(LC),   &
1405                  BETR_TBL(LC), &
1406                  BETB_TBL(LC), &
1407                  BETG_TBL(LC), &
1408                  FRC_URB_TBL(LC)
1409    END DO
1410
1411    READ(11,*)
1412    READ(11,*) CAPR_DATA
1413    READ(11,*)
1414    READ(11,*) CAPB_DATA
1415    READ(11,*)
1416    READ(11,*) CAPG_DATA
1417    READ(11,*)
1418    READ(11,*) AKSR_DATA
1419    READ(11,*)
1420    READ(11,*) AKSB_DATA
1421    READ(11,*)
1422    READ(11,*) AKSG_DATA
1423    READ(11,*)
1424    READ(11,*) ALBR_DATA
1425    READ(11,*)
1426    READ(11,*) ALBB_DATA
1427    READ(11,*)
1428    READ(11,*) ALBG_DATA
1429    READ(11,*)
1430    READ(11,*) EPSR_DATA
1431    READ(11,*)
1432    READ(11,*) EPSB_DATA
1433    READ(11,*)
1434    READ(11,*) EPSG_DATA
1435    READ(11,*)
1436    READ(11,*) Z0R_DATA
1437    READ(11,*)
1438    READ(11,*) Z0B_DATA
1439    READ(11,*)
1440    READ(11,*) Z0G_DATA
1441    READ(11,*)
1442    READ(11,*) Z0HR_DATA
1443    READ(11,*)
1444    READ(11,*) Z0HB_DATA
1445    READ(11,*)
1446    READ(11,*) Z0HG_DATA
1447    READ(11,*)
1448 !   READ(11,*) num_roof_layers
1449    READ(11,*) dummy
1450    READ(11,*)
1451 !   READ(11,*) num_wall_layers
1452    READ(11,*) dummy
1453    READ(11,*)
1454 !   READ(11,*) num_road_layers
1455    READ(11,*) dummy
1456
1457    num_roof_layers = num_soil_layers
1458    num_wall_layers = num_soil_layers
1459    num_road_layers = num_soil_layers
1460
1461    DO K=1,num_roof_layers
1462       READ(11,*)
1463       READ(11,*) DZR(K)
1464    END DO
1465
1466    DO K=1,num_wall_layers
1467       READ(11,*)
1468       READ(11,*) DZB(K)
1469    END DO
1470
1471    DO K=1,num_road_layers
1472       READ(11,*)
1473       READ(11,*) DZG(K)
1474    END DO
1475
1476    READ(11,*)
1477    READ(11,*) BOUNDR_DATA
1478    READ(11,*)
1479    READ(11,*) BOUNDB_DATA
1480    READ(11,*)
1481    READ(11,*) BOUNDG_DATA
1482    READ(11,*)
1483    READ(11,*) TRLEND_DATA
1484    READ(11,*)
1485    READ(11,*) TBLEND_DATA
1486    READ(11,*)
1487    READ(11,*) TGLEND_DATA
1488    READ(11,*)
1489    READ(11,*) CH_SCHEME_DATA
1490    READ(11,*)
1491    READ(11,*) TS_SCHEME_DATA
1492 ! Miao, 2007/01/17, cal. ah
1493    READ(11,*)
1494    READ(11,*) ahoption
1495    if(ahoption==1) then
1496      READ(11,*)
1497      READ(11,*) (ahdiuprf(k),k=1,24)
1498    endif
1499
1500    CLOSE(11)
1501
1502 ! Calculate Sky View Factor
1503
1504    DO LC = 1, ICATE
1505       DHGT=HGT_TBL(LC)/100.
1506       HGT=0.
1507       VFWS=0.
1508       HGT=HGT_TBL(LC)-DHGT/2.
1509       do k=1,99
1510          HGT=HGT-DHGT
1511          VFWS=VFWS+0.25*(1.-HGT/SQRT(HGT**2.+RW_TBL(LC)**2.))
1512       end do
1513
1514      VFWS=VFWS/99.
1515      VFWS=VFWS*2.
1516
1517      VFGS=1.-2.*VFWS*HGT_TBL(LC)/RW_TBL(LC)
1518      SVF_TBL(LC)=VFGS
1519    END DO
1520
1521    END SUBROUTINE urban_param_init
1522 !===========================================================================
1523 !
1524 ! subroutine urban_var_init: initialization of urban state variables
1525 !
1526 !===========================================================================
1527    SUBROUTINE urban_var_init(TSURFACE0_URB,TLAYER0_URB,TDEEP0_URB,IVGTYP,     & ! in
1528                              ims,ime,jms,jme,num_soil_layers,                 & ! in
1529 !                             num_roof_layers,num_wall_layers,num_road_layers, & ! in
1530                              restart,                                      & !in
1531                              XXXR_URB2D,XXXB_URB2D,XXXG_URB2D,XXXC_URB2D,  & ! inout
1532                              TR_URB2D,TB_URB2D,TG_URB2D,TC_URB2D,QC_URB2D, & ! inout
1533                              TRL_URB3D,TBL_URB3D,TGL_URB3D,                & ! inout
1534                              SH_URB2D,LH_URB2D,G_URB2D,RN_URB2D,           & ! inout
1535                              TS_URB2D, FRC_URB2D, UTYPE_URB2D)               ! inout
1536    IMPLICIT NONE
1537
1538    INTEGER, INTENT(IN) :: ims,ime,jms,jme,num_soil_layers
1539 !   INTEGER, INTENT(IN) :: num_roof_layers, num_wall_layers, num_road_layers
1540
1541    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN)                    :: TSURFACE0_URB
1542    REAL, DIMENSION( ims:ime, 1:num_soil_layers, jms:jme ), INTENT(IN) :: TLAYER0_URB
1543    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN)                    :: TDEEP0_URB
1544    INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(IN)                 :: IVGTYP
1545    LOGICAL , INTENT(IN) :: restart
1546
1547    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TR_URB2D
1548    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TB_URB2D
1549    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TG_URB2D
1550    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TC_URB2D
1551    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: QC_URB2D
1552    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXR_URB2D
1553    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXB_URB2D
1554    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXG_URB2D
1555    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXC_URB2D
1556
1557 !   REAL, DIMENSION(ims:ime, 1:num_roof_layers, jms:jme), INTENT(INOUT) :: TRL_URB3D
1558 !   REAL, DIMENSION(ims:ime, 1:num_wall_layers, jms:jme), INTENT(INOUT) :: TBL_URB3D
1559 !   REAL, DIMENSION(ims:ime, 1:num_road_layers, jms:jme), INTENT(INOUT) :: TGL_URB3D
1560    REAL, DIMENSION(ims:ime, 1:num_soil_layers, jms:jme), INTENT(INOUT) :: TRL_URB3D
1561    REAL, DIMENSION(ims:ime, 1:num_soil_layers, jms:jme), INTENT(INOUT) :: TBL_URB3D
1562    REAL, DIMENSION(ims:ime, 1:num_soil_layers, jms:jme), INTENT(INOUT) :: TGL_URB3D
1563
1564    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: SH_URB2D
1565    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: LH_URB2D
1566    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: G_URB2D
1567    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: RN_URB2D
1568    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TS_URB2D
1569 !
1570    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: FRC_URB2D
1571    INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: UTYPE_URB2D
1572    INTEGER                                            :: UTYPE_URB
1573
1574    INTEGER :: I,J,K
1575
1576    DO I=ims,ime
1577     DO J=jms,jme
1578
1579 !      XXXR_URB2D(I,J)=0.
1580 !      XXXB_URB2D(I,J)=0.
1581 !      XXXG_URB2D(I,J)=0.
1582 !      XXXC_URB2D(I,J)=0.
1583
1584       SH_URB2D(I,J)=0.
1585       LH_URB2D(I,J)=0.
1586       G_URB2D(I,J)=0.
1587       RN_URB2D(I,J)=0.
1588 !m
1589       FRC_URB2D(I,J)=0.
1590       UTYPE_URB2D(I,J)=0.
1591
1592             IF( IVGTYP(I,J) == 1)  THEN
1593                UTYPE_URB2D(I,J) = 2  ! for default. high-density
1594                UTYPE_URB = UTYPE_URB2D(I,J)  ! for default. high-density
1595                FRC_URB2D(I,J) = FRC_URB_TBL(UTYPE_URB)
1596             ENDIF
1597             IF( IVGTYP(I,J) == 31) THEN
1598                UTYPE_URB2D(I,J) = 3  ! low-density residential
1599                UTYPE_URB = UTYPE_URB2D(I,J)  ! low-density residential
1600                FRC_URB2D(I,J) = FRC_URB_TBL(UTYPE_URB)
1601             ENDIF
1602             IF( IVGTYP(I,J) == 32) THEN
1603                UTYPE_URB2D(I,J) = 2  ! high-density
1604                UTYPE_URB = UTYPE_URB2D(I,J)  ! high-density
1605                FRC_URB2D(I,J) = FRC_URB_TBL(UTYPE_URB)
1606             ENDIF
1607             IF( IVGTYP(I,J) == 33) THEN
1608                UTYPE_URB2D(I,J) = 1  !  Commercial/Industrial/Transportation
1609                UTYPE_URB = UTYPE_URB2D(I,J)  !  Commercial/Industrial/Transportation
1610                FRC_URB2D(I,J) = FRC_URB_TBL(UTYPE_URB)
1611             ENDIF
1612
1613
1614       QC_URB2D(I,J)=0.01
1615
1616        IF (.not.restart) THEN
1617
1618       XXXR_URB2D(I,J)=0.
1619       XXXB_URB2D(I,J)=0.
1620       XXXG_URB2D(I,J)=0.
1621       XXXC_URB2D(I,J)=0.
1622
1623
1624       TC_URB2D(I,J)=TSURFACE0_URB(I,J)+0.
1625       TR_URB2D(I,J)=TSURFACE0_URB(I,J)+0.
1626       TB_URB2D(I,J)=TSURFACE0_URB(I,J)+0.
1627       TG_URB2D(I,J)=TSURFACE0_URB(I,J)+0.
1628 !
1629       TS_URB2D(I,J)=TSURFACE0_URB(I,J)+0.
1630
1631 !      DO K=1,num_roof_layers
1632 !     DO K=1,num_soil_layers
1633 !         TRL_URB3D(I,1,J)=TLAYER0_URB(I,1,J)+0.
1634 !         TRL_URB3D(I,2,J)=TLAYER0_URB(I,2,J)+0.
1635 !         TRL_URB3D(I,3,J)=TLAYER0_URB(I,3,J)+0.
1636 !         TRL_URB3D(I,4,J)=TLAYER0_URB(I,4,J)+0.
1637
1638           TRL_URB3D(I,1,J)=TLAYER0_URB(I,1,J)+0.
1639           TRL_URB3D(I,2,J)=0.5*(TLAYER0_URB(I,1,J)+TLAYER0_URB(I,2,J))
1640           TRL_URB3D(I,3,J)=TLAYER0_URB(I,2,J)+0.
1641           TRL_URB3D(I,4,J)=TLAYER0_URB(I,2,J)+(TLAYER0_URB(I,3,J)-TLAYER0_URB(I,2,J))*0.29
1642 !     END DO
1643
1644 !      DO K=1,num_wall_layers
1645 !      DO K=1,num_soil_layers
1646 !m        TBL_URB3D(I,1,J)=TLAYER0_URB(I,1,J)+0.
1647 !m        TBL_URB3D(I,2,J)=TLAYER0_URB(I,2,J)+0.
1648 !m        TBL_URB3D(I,3,J)=TLAYER0_URB(I,3,J)+0.
1649 !m        TBL_URB3D(I,4,J)=TLAYER0_URB(I,4,J)+0.
1650
1651         TBL_URB3D(I,1,J)=TLAYER0_URB(I,1,J)+0.
1652         TBL_URB3D(I,2,J)=0.5*(TLAYER0_URB(I,1,J)+TLAYER0_URB(I,2,J))
1653         TBL_URB3D(I,3,J)=TLAYER0_URB(I,2,J)+0.
1654         TBL_URB3D(I,4,J)=TLAYER0_URB(I,2,J)+(TLAYER0_URB(I,3,J)-TLAYER0_URB(I,2,J))*0.29
1655 !      END DO
1656
1657 !      DO K=1,num_road_layers
1658       DO K=1,num_soil_layers
1659         TGL_URB3D(I,K,J)=TLAYER0_URB(I,K,J)+0.
1660       END DO
1661
1662      ENDIF        !restart
1663     END DO
1664    END DO
1665
1666    RETURN
1667    END SUBROUTINE urban_var_init
1668 !===========================================================================
1669 !
1670 !  force_restore
1671 !
1672 !===========================================================================
1673    SUBROUTINE force_restore(CAP,AKS,DELT,S,R,H,LE,TSLEND,TSP,TS)
1674
1675      REAL, INTENT(IN)  :: CAP,AKS,DELT,S,R,H,LE,TSLEND,TSP
1676      REAL, INTENT(OUT) :: TS
1677      REAL              :: C1,C2
1678
1679      C2=24.*3600./2./3.14159
1680      C1=SQRT(0.5*C2*CAP*AKS)
1681
1682      TS = TSP + DELT*( (S+R-H-LE)/C1 -(TSP-TSLEND)/C2 )
1683
1684    END SUBROUTINE force_restore
1685 !===========================================================================
1686 !
1687 !  bisection (not used)
1688 !
1689 !==============================================================================
1690    SUBROUTINE bisection(TSP,PS,S,EPS,RX,SIG,RHO,CP,CH,UA,QA,TA,EL,BET,AKS,TSL,DZ,TS)
1691
1692      REAL, INTENT(IN) :: TSP,PS,S,EPS,RX,SIG,RHO,CP,CH,UA,QA,TA,EL,BET,AKS,TSL,DZ
1693      REAL, INTENT(OUT) :: TS
1694      REAL :: ES,QS0,R,H,ELE,G0,F1,F
1695
1696      TS1 = TSP - 5.
1697      TS2 = TSP + 5.
1698
1699      DO ITERATION = 1,22
1700
1701        ES=6.11*EXP( (2.5*10.**6./461.51)*(TS1-273.15)/(273.15*TS1) )
1702        QS0=0.622*ES/(PS-0.378*ES)
1703        R=EPS*(RX-SIG*(TS1**4.)/60.)
1704        H=RHO*CP*CH*UA*(TS1-TA)*100.
1705        ELE=RHO*EL*CH*UA*BET*(QS0-QA)*100.
1706        G0=AKS*(TS1-TSL)/(DZ/2.)
1707        F1= S + R - H - ELE - G0
1708
1709        TS=0.5*(TS1+TS2)
1710
1711        ES=6.11*EXP( (2.5*10.**6./461.51)*(TS-273.15)/(273.15*TS) )
1712        QS0=0.622*ES/(PS-0.378*ES)
1713        R=EPS*(RX-SIG*(TS**4.)/60.)
1714        H=RHO*CP*CH*UA*(TS-TA)*100.
1715        ELE=RHO*EL*CH*UA*BET*(QS0-QA)*100.
1716        G0=AKS*(TS-TSL)/(DZ/2.)
1717        F = S + R - H - ELE - G0
1718
1719        IF (F1*F > 0.0) THEN
1720          TS1=TS
1721         ELSE
1722          TS2=TS
1723        END IF
1724
1725      END DO
1726
1727      RETURN
1728 END SUBROUTINE bisection
1729 !===========================================================================
1730 END MODULE module_sf_urban