phys/module_sf_urban.F

   1 MODULE module_sf_urban
   2 #if defined(mpas)
   3 use mpas_atmphys_utilities, only: physics_message,physics_error_fatal
   4 #define FATAL_ERROR(M) call physics_error_fatal( M )
   5 #define WRITE_MESSAGE(M) call physics_message( M )
   6 #else
   7 use module_wrf_error
   8 #define FATAL_ERROR(M) call wrf_error_fatal( M )
   9 #define WRITE_MESSAGE(M) call wrf_message( M )
  10 #endif
  11
  12 !===============================================================================
  13 !     Single-Layer Urban Canopy Model for WRF Noah-LSM
  14 !     Original Version: 2002/11/06 by Hiroyuki Kusaka
  15 !     Last Update:      2006/08/24 by Fei Chen and Mukul Tewari (NCAR/RAL)
  16 !===============================================================================
  17
  18    CHARACTER(LEN=4)                :: LU_DATA_TYPE
  19
  20    INTEGER                         :: ICATE
  21
  22    REAL, ALLOCATABLE, DIMENSION(:) :: ZR_TBL
  23    REAL, ALLOCATABLE, DIMENSION(:) :: Z0C_TBL
  24    REAL, ALLOCATABLE, DIMENSION(:) :: Z0HC_TBL
  25    REAL, ALLOCATABLE, DIMENSION(:) :: ZDC_TBL
  26    REAL, ALLOCATABLE, DIMENSION(:) :: SVF_TBL
  27    REAL, ALLOCATABLE, DIMENSION(:) :: R_TBL
  28    REAL, ALLOCATABLE, DIMENSION(:) :: RW_TBL
  29    REAL, ALLOCATABLE, DIMENSION(:) :: HGT_TBL
  30    REAL, ALLOCATABLE, DIMENSION(:) :: AH_TBL
  31    REAL, ALLOCATABLE, DIMENSION(:) :: ALH_TBL
  32    REAL, ALLOCATABLE, DIMENSION(:) :: BETR_TBL
  33    REAL, ALLOCATABLE, DIMENSION(:) :: BETB_TBL
  34    REAL, ALLOCATABLE, DIMENSION(:) :: BETG_TBL
  35    REAL, ALLOCATABLE, DIMENSION(:) :: FRC_URB_TBL
  36
  37    REAL, ALLOCATABLE, DIMENSION(:) :: COP_TBL
  38    REAL, ALLOCATABLE, DIMENSION(:) :: BLDAC_FRC_TBL
  39    REAL, ALLOCATABLE, DIMENSION(:) :: COOLED_FRC_TBL
  40    REAL, ALLOCATABLE, DIMENSION(:) :: PWIN_TBL
  41    REAL, ALLOCATABLE, DIMENSION(:) :: BETA_TBL
  42    INTEGER, ALLOCATABLE, DIMENSION(:) :: SW_COND_TBL
  43    REAL, ALLOCATABLE, DIMENSION(:) :: TIME_ON_TBL
  44    REAL, ALLOCATABLE, DIMENSION(:) :: TIME_OFF_TBL
  45    REAL, ALLOCATABLE, DIMENSION(:) :: TARGTEMP_TBL
  46    REAL, ALLOCATABLE, DIMENSION(:) :: GAPTEMP_TBL
  47    REAL, ALLOCATABLE, DIMENSION(:) :: TARGHUM_TBL
  48    REAL, ALLOCATABLE, DIMENSION(:) :: GAPHUM_TBL
  49    REAL, ALLOCATABLE, DIMENSION(:) :: PERFLO_TBL
  50    REAL, ALLOCATABLE, DIMENSION(:) :: PV_FRAC_ROOF_TBL !GRZ
  51    REAL, ALLOCATABLE, DIMENSION(:) :: GR_FRAC_ROOF_TBL !GRZ
  52    INTEGER :: GR_FLAG_TBL !GRZ
  53    INTEGER :: GR_TYPE_TBL !GRZ
  54    REAL, DIMENSION(1:24)           :: IRHO_TBL
  55    REAL, ALLOCATABLE, DIMENSION(:) :: HSESF_TBL
  56
  57    REAL, ALLOCATABLE, DIMENSION(:) :: CAPR_TBL, CAPB_TBL, CAPG_TBL
  58    REAL, ALLOCATABLE, DIMENSION(:) :: AKSR_TBL, AKSB_TBL, AKSG_TBL
  59    REAL, ALLOCATABLE, DIMENSION(:) :: ALBR_TBL, ALBB_TBL, ALBG_TBL
  60    REAL, ALLOCATABLE, DIMENSION(:) :: EPSR_TBL, EPSB_TBL, EPSG_TBL
  61    REAL, ALLOCATABLE, DIMENSION(:) :: Z0R_TBL,  Z0B_TBL,  Z0G_TBL
  62    REAL, ALLOCATABLE, DIMENSION(:) :: SIGMA_ZED_TBL
  63    REAL, ALLOCATABLE, DIMENSION(:) :: Z0HB_TBL, Z0HG_TBL
  64    REAL, ALLOCATABLE, DIMENSION(:) :: TRLEND_TBL, TBLEND_TBL, TGLEND_TBL
  65    REAL, ALLOCATABLE, DIMENSION(:) :: AKANDA_URBAN_TBL
  66 !for BEP
  67
  68    ! MAXDIRS :: The maximum number of street directions we're allowed to define
  69    INTEGER, PARAMETER :: MAXDIRS = 3
  70    ! MAXHGTS :: The maximum number of building height bins we're allowed to define
  71    INTEGER, PARAMETER :: MAXHGTS = 50
  72
  73    INTEGER, ALLOCATABLE, DIMENSION(:)   :: NUMDIR_TBL
  74    REAL,    ALLOCATABLE, DIMENSION(:,:) :: STREET_DIRECTION_TBL
  75    REAL,    ALLOCATABLE, DIMENSION(:,:) :: STREET_WIDTH_TBL
  76    REAL,    ALLOCATABLE, DIMENSION(:,:) :: BUILDING_WIDTH_TBL
  77    INTEGER, ALLOCATABLE, DIMENSION(:)   :: NUMHGT_TBL
  78    REAL,    ALLOCATABLE, DIMENSION(:,:) :: HEIGHT_BIN_TBL
  79    REAL,    ALLOCATABLE, DIMENSION(:,:) :: HPERCENT_BIN_TBL
  80 !end BEP
  81    INTEGER                         :: BOUNDR_DATA,BOUNDB_DATA,BOUNDG_DATA
  82    INTEGER                         :: CH_SCHEME_DATA, TS_SCHEME_DATA
  83    INTEGER                         :: ahoption        ! Miao, 2007/01/17, cal. ah
  84    REAL, DIMENSION(1:24)           :: ahdiuprf        ! ah diurnal profile, tloc: 1-24
  85    REAL, DIMENSION(1:24)           :: hsequip_tbl
  86
  87 !===Yang, 2014/10/08, urban hydrological processes for single layer UCM===
  88    INTEGER                         :: IMP_SCHEME, IRI_SCHEME
  89    INTEGER                         :: alhoption       ! anthropogenic latent heat option
  90    INTEGER                         :: groption        ! anthropogenic latent heat option
  91    REAL                            :: fgr             ! green roof fraction
  92    REAL                            :: oasis           ! urban oasis parameter
  93    REAL, DIMENSION(1:4)            :: DZGR            ! Layer depth of green roof
  94    REAL, DIMENSION(1:4)            :: alhseason       ! seasonal variation of alh
  95    REAL, DIMENSION(1:48)           :: alhdiuprf       ! alh diurnal profile, tloc2: 1-48
  96    REAL, DIMENSION(1:3)            :: porimp          ! porosity of pavement over impervious surface
  97    REAL, DIMENSION(1:3)            :: dengimp         ! maximum water-holding depth of pavement
  98
  99 !===end hydrological processes===
 100
 101    INTEGER                         :: allocate_status
 102
 103 !   INTEGER                         :: num_roof_layers
 104 !   INTEGER                         :: num_wall_layers
 105 !   INTEGER                         :: num_road_layers
 106
 107    CHARACTER (LEN=256) , PRIVATE   :: mesg
 108
 109    CONTAINS
 110
 111 !===============================================================================
 112 !
 113 ! Author:
 114 !      Hiroyuki KUSAKA, PhD
 115 !      University of Tsukuba, JAPAN
 116 !      (CRIEPI, NCAR/MMM visiting scientist, 2002-2004)
 117 !      kusaka@ccs.tsukuba.ac.jp
 118 !
 119 ! Co-Researchers:
 120 !     Fei CHEN, PhD
 121 !      NCAR/RAP feichen@ucar.edu
 122 !     Mukul TEWARI, PhD
 123 !      NCAR/RAP mukul@ucar.edu
 124 !
 125 ! Purpose:
 126 !     Calculate surface temeprature, fluxes, canopy air temperature, and canopy wind
 127 !
 128 ! Subroutines:
 129 !     module_sf_urban
 130 !       |- urban
 131 !            |- read_param
 132 !            |- mos or jurges
 133 !            |- multi_layer or force_restore
 134 !       |- urban_param_init <-- URBPARM.TBL
 135 !       |- urban_var_init
 136 !
 137 ! Input Data from WRF [MKS unit]:
 138 !
 139 !     UTYPE  [-]     : Urban type. 1=Commercial/Industrial; 2=High-intensity residential;
 140 !                    : 3=low-intensity residential
 141 !     TA     [K]     : Potential temperature at 1st wrf level (absolute temp)
 142 !     QA     [kg/kg] : Mixing ratio at 1st atmospheric level
 143 !     UA     [m/s]   : Wind speed at 1st atmospheric level
 144 !     SSG    [W/m/m] : Short wave downward radiation at a flat surface
 145 !                      Note this is the total of direct and diffusive solar
 146 !                      downward radiation. If without two components, the
 147 !                      single solar downward can be used instead.
 148 !                      SSG = SSGD + SSGQ
 149 !     LSOLAR [-]     : Indicating the input type of solar downward radiation
 150 !                      True: both direct and diffusive solar radiation
 151 !                      are available
 152 !                      False: only total downward ridiation is available.
 153 !     SSGD   [W/m/m] : Direct solar radiation at a flat surface
 154 !                      if SSGD is not available, one can assume a ratio SRATIO
 155 !                      (e.g., 0.7), so that SSGD = SRATIO*SSG
 156 !     SSGQ   [W/m/m] : Diffuse solar radiation at a flat surface
 157 !                      If SSGQ is not available, SSGQ = SSG - SSGD
 158 !     LLG    [W/m/m] : Long wave downward radiation at a flat surface
 159 !     RAIN   [mm/h]  : Precipitation
 160 !     RHOO   [kg/m/m/m] : Air density
 161 !     ZA     [m]     : First atmospheric level
 162 !                      as a lowest boundary condition
 163 !     DECLIN [rad]   : solar declination
 164 !     COSZ           : = sin(fai)*sin(del)+cos(fai)*cos(del)*cos(omg)
 165 !     OMG    [rad]   : solar hour angle
 166 !     XLAT   [deg]   : latitude
 167 !     DELT   [sec]   : Time step
 168 !     ZNT    [m]     : Roughnes length
 169 !
 170 ! Output Data to WRF [MKS unit]:
 171 !
 172 !     TS  [K]            : Surface potential temperature (absolute temp)
 173 !     QS  [-]            : Surface humidity
 174 !
 175 !     SH  [W/m/m/]       : Sensible heat flux, = FLXTH*RHOO*CPP
 176 !     LH  [W/m/m]        : Latent heat flux, = FLXHUM*RHOO*ELL
 177 !     LH_INEMATIC [kg/m/m/sec]: Moisture Kinematic flux, = FLXHUM*RHOO
 178 !     SW  [W/m/m]        : Upward shortwave radiation flux,
 179 !                          = SSG-SNET*697.7*60. (697.7*60.=100.*100.*4.186)
 180 !     ALB [-]            : Time-varying albedo
 181 !     LW  [W/m/m]        : Upward longwave radiation flux,
 182 !                          = LNET*697.7*60.-LLG
 183 !     G   [W/m/m]        : Heat Flux into the Ground
 184 !     RN  [W/m/m]        : Net radiation
 185 !
 186 !     PSIM  [-]          : Diagnostic similarity stability function for momentum
 187 !     PSIH  [-]          : Diagnostic similarity stability function for heat
 188 !
 189 !     TC  [K]            : Diagnostic canopy air temperature
 190 !     QC  [-]            : Diagnostic canopy humidity
 191 !
 192 !     TH2 [K]            : Diagnostic potential temperature at 2 m
 193 !     Q2  [-]            : Diagnostic humidity at 2 m
 194 !     U10 [m/s]          : Diagnostic u wind component at 10 m
 195 !     V10 [m/s]          : Diagnostic v wind component at 10 m
 196 !
 197 !     CHS, CHS2 [m/s]    : CH*U at ZA, CH*U at 2 m (not used)
 198 !
 199 ! Important parameters:
 200 !
 201 ! Morphology of the urban canyon:
 202 ! These parameters assigned in the URBPARM.TBL
 203 !
 204 !    ZR  [m]             : roof level (building height)
 205 !    SIGMA_ZED [m]       : Standard Deviation of roof height
 206 !    ROOF_WIDTH [m]      : roof (i.e., building) width
 207 !    ROAD_WIDTH [m]      : road width
 208 !
 209 ! Parameters derived from the morphological terms above.
 210 ! These parameters are computed in the code.
 211 !
 212 !    HGT [-]             : normalized building height
 213 !    SVF [-]             : sky view factor
 214 !    R   [-]             : Normalized roof width (a.k.a. "building coverage ratio")
 215 !    RW  [-]             : = 1 - R
 216 !    Z0C [m]             : Roughness length above canyon for momentum (1/10 of ZR)
 217 !    Z0HC [m]            : Roughness length above canyon for heat (1/10 of Z0C)
 218 !    ZDC [m]             : Zero plane displacement height (1/5 of ZR)
 219 !
 220 ! Following parameter are assigned in run/URBPARM.TBL
 221 !
 222 !  AH  [ W m{-2} ]        : anthropogenic heat  ( W m{-2} in the table, converted internally to cal cm{-2} )
 223 !  ALH [ W m{-2} ]        : anthropogenic latent heat  ( W m{-2} in the table, converted internally to cal cm{-2} )
 224 !  CAPR[ J m{-3} K{-1} ]  : heat capacity of roof ( units converted in code to [ cal cm{-3} deg{-1} ] )
 225 !  CAPB[ J m{-3} K{-1} ]  : heat capacity of building wall ( units converted in code to [ cal cm{-3} deg{-1} ] )
 226 !  CAPG[ J m{-3} K{-1} ]  : heat capacity of road ( units converted in code to [ cal cm{-3} deg{-1} ] )
 227 !  AKSR [ J m{-1} s{-1} K{-1} ] : thermal conductivity of roof ( units converted in code to [ cal cm{-1} s{-1} deg{-1} ] )
 228 !  AKSB [ J m{-1} s{-1} K{-1} ] : thermal conductivity of building wall ( units converted in code to [ cal cm{-1} s{-1} deg{-1} ] )
 229 !  AKSG [ J m{-1} s{-1} K{-1} ] : thermal conductivity of road ( units converted in code to [ cal cm{-1} s{-1} deg{-1} ] )
 230 !  ALBR [-]               : surface albedo of roof
 231 !  ALBB [-]               : surface albedo of building wall
 232 !  ALBG [-]               : surface albedo of road
 233 !  EPSR [-]               : surface emissivity of roof
 234 !  EPSB [-]               : surface emissivity of building wall
 235 !  EPSG [-]               : surface emissivity of road
 236 !  Z0B [m]                : roughness length for momentum of building wall (only for CH_SCHEME = 1)
 237 !  Z0G [m]                : roughness length for momentum of road (only for CH_SCHEME = 1)
 238 !  Z0HB [m]               : roughness length for heat of building wall (only for CH_SCHEME = 1)
 239 !  Z0HG [m]               : roughness length for heat of road
 240 !  num_roof_layers        : number of layers within roof
 241 !  num_wall_layers        : number of layers within building walls
 242 !  num_road_layers        : number of layers within below road surface
 243 !   NOTE: for now, these layers are defined as same as the number of soil layers in namelist.input
 244 !  DZR [cm]               : thickness of each roof layer
 245 !  DZB [cm]               : thickness of each building wall layer
 246 !  DZG [cm]               : thickness of each ground layer
 247 !  BOUNDR [integer 1 or 2] : Boundary Condition for Roof Layer Temp [1: Zero-Flux, 2: T = Constant]
 248 !  BOUNDB [integer 1 or 2] : Boundary Condition for Building Wall Layer Temp [1: Zero-Flux, 2: T = Constant]
 249 !  BOUNDG [integer 1 or 2] : Boundary Condition for Road Layer Temp [1: Zero-Flux, 2: T = Constant]
 250 !  TRLEND [K]             : lower boundary condition of roof temperature
 251 !  TBLEND [K]             : lower boundary condition of building temperature
 252 !  TGLEND [K]             : lower boundary condition of ground temperature
 253 !  CH_SCHEME [integer 1 or 2] : Sfc exchange scheme used for building wall and road
 254 !                         [1: M-O Similarity Theory,  2: Empirical Form (recommend)]
 255 !  TS_SCHEME [integer 1 or 2] : Scheme for computing surface temperature (for roof, wall, and road)
 256 !                         [1: 4-layer model,  2: Force-Restore method]
 257 !  IMP_SCHEME[integer 1 or 2] : Evaporation scheme for impervious surfaces (roof, wall, and road)
 258 !                         [1: Hypothesized evaporation during large rainfall events
 259 !                         [2: Water-holding scheme over impervious surface
 260 !  IRI_SCHEME[integer 0 or 1] : Scheme for urban irrigation
 261 !                         [0: No irrigation,  1: Summertime (May-Sep) irrigation everyday at 9pm]
 262 !for BEP
 263 !  numdir [ - ]             : Number of street directions defined for a particular urban category
 264 !  street_direction [ deg ] : Direction of streets for a particular urban category and a particular street direction
 265 !  street_width [ m ]       : Width of street for a particular urban category and a particular street direction
 266 !  building_width [ m ]     : Width of buildings for a particular urban category and a particular street direction
 267 !  numhgt [ - ]             : Number of building height levels defined for a particular urban category
 268 !  height_bin [ m ]         : Building height bins defined for a particular urban category.
 269 !  hpercent_bin [ % ]       : Percentage of a particular urban category populated by buildings of particular height bins
 270 !end BEP
 271 ! Moved from URBPARM.TBL
 272 !
 273 !  BETR [-]            : minimum moisture availability of roof
 274 !  BETB [-]            : minimum moisture availability of building wall
 275 !  BETG [-]            : minimum moisture availability of road
 276 !  Z0R [m]                : roughness length for momentum of roof
 277 !  Z0HB [m]               : roughness length for heat of building wall (only for CH_SCHEME = 1)
 278 !  Z0HG [m]               : roughness length for heat of road
 279 !  num_roof_layers        : number of layers within roof
 280 !  num_wall_layers        : number of layers within building walls
 281 !  num_road_layers        : number of layers within below road surface
 282 !   NOTE: for now, these layers are defined as same as the number of soil layers in namelist.input
 283 !
 284 ! References:
 285 !     Kusaka and Kimura (2004) J.Appl.Meteor., vol.43, p1899-1910
 286 !     Kusaka and Kimura (2004) J.Meteor.Soc.Japan, vol.82, p45-65
 287 !     Kusaka et al. (2001) Bound.-Layer Meteor., vol.101, p329-358
 288 !
 289 ! History:
 290 !     2014/10,         modified by Jiachuan Yang (ASU)
 291 !     2006/06          modified by H. Kusaka (Univ. Tsukuba), M. Tewari
 292 !     2005/10/26,      modified by Fei Chen, Mukul Tewari
 293 !     2003/07/21 WRF , modified  by H. Kusaka of CRIEPI (NCAR/MMM)
 294 !     2001/08/26 PhD , modified  by H. Kusaka of CRIEPI (Univ.Tsukuba)
 295 !     1999/08/25 LCM , developed by H. Kusaka of CRIEPI (Univ.Tsukuba)
 296 !
 297 !===============================================================================
 298 !
 299 !  subroutine urban:
 300 !
 301 !===============================================================================
 302
 303    SUBROUTINE urban(LSOLAR,                                           & ! L
 304                     num_roof_layers,num_wall_layers,num_road_layers,  & ! I
 305                     DZR,DZB,DZG,                                      & ! I
 306                     UTYPE,TA,QA,UA,U1,V1,SSG,SSGD,SSGQ,LLG,RAIN,RHOO, & ! I
 307                     ZA,DECLIN,COSZ,OMG,XLAT,DELT,ZNT,                 & ! I
 308                     CHS, CHS2,                                        & ! I
 309                     TR, TB, TG, TC, QC, UC,                           & ! H
 310                     TRL,TBL,TGL,                                      & ! H
 311                     XXXR, XXXB, XXXG, XXXC,                           & ! H
 312                     TS,QS,SH,LH,LH_KINEMATIC,                         & ! O
 313                     SW,ALB,LW,G,RN,PSIM,PSIH,                         & ! O
 314                     GZ1OZ0,                                           & ! O
 315                     CMR_URB,CHR_URB,CMC_URB,CHC_URB,                  & ! I/O
 316                     U10,V10,TH2,Q2,UST,mh_urb,stdh_urb,lf_urb,        & ! O
 317                     lp_urb,hgt_urb,frc_urb,lb_urb,zo_check,           & ! O
 318                     CMCR,TGR,TGRL,SMR,CMGR_URB,CHGR_URB,jmonth,       & ! H
 319                     DRELR,DRELB,DRELG,FLXHUMR,FLXHUMB,FLXHUMG)
 320
 321    IMPLICIT NONE
 322
 323    REAL, PARAMETER    :: CP=0.24          ! heat capacity of dry air  [cgs unit]
 324    REAL, PARAMETER    :: EL=583.          ! latent heat of vaporation [cgs unit]
 325    REAL, PARAMETER    :: SIG=8.17E-11     ! stefun bolzman constant   [cgs unit]
 326    REAL, PARAMETER    :: SIG_SI=5.67E-8   !                           [MKS unit]
 327    REAL, PARAMETER    :: AK=0.4           ! kalman const.                [-]
 328    REAL, PARAMETER    :: PI=3.14159       ! pi                           [-]
 329    REAL, PARAMETER    :: TETENA=7.5       ! const. of Tetens Equation    [-]
 330    REAL, PARAMETER    :: TETENB=237.3     ! const. of Tetens Equation    [-]
 331    REAL, PARAMETER    :: SRATIO=0.75      ! ratio between direct/total solar [-]
 332
 333    REAL, PARAMETER    :: CPP=1004.5       ! heat capacity of dry air    [J/K/kg]
 334    REAL, PARAMETER    :: ELL=2.442E+06    ! latent heat of vaporization [J/kg]
 335    REAL, PARAMETER    :: XKA=2.4E-5
 336
 337 !-------------------------------------------------------------------------------
 338 ! C: configuration variables
 339 !-------------------------------------------------------------------------------
 340
 341    LOGICAL, INTENT(IN) :: LSOLAR  ! logical    [true=both, false=SSG only]
 342
 343 !  The following variables are also model configuration variables, but are
 344 !  defined in the URBAN.TBL and in the contains statement in the top of
 345 !  the module_urban_init, so we should not declare them here.
 346
 347   INTEGER, INTENT(IN) :: num_roof_layers
 348   INTEGER, INTENT(IN) :: num_wall_layers
 349   INTEGER, INTENT(IN) :: num_road_layers
 350
 351
 352   REAL, INTENT(IN), DIMENSION(1:num_roof_layers) :: DZR ! grid interval of roof layers [cm]
 353   REAL, INTENT(IN), DIMENSION(1:num_wall_layers) :: DZB ! grid interval of wall layers [cm]
 354   REAL, INTENT(IN), DIMENSION(1:num_road_layers) :: DZG ! grid interval of road layers [cm]
 355
 356 !-------------------------------------------------------------------------------
 357 ! I: input variables from LSM to Urban
 358 !-------------------------------------------------------------------------------
 359
 360    INTEGER, INTENT(IN) :: UTYPE ! urban type [1=Commercial/Industrial, 2=High-intensity residential,
 361                                 ! 3=low-intensity residential]
 362    INTEGER, INTENT(IN) :: jmonth! current month
 363    REAL, INTENT(IN)    :: TA   ! potential temp at 1st atmospheric level [K]
 364    REAL, INTENT(IN)    :: QA   ! mixing ratio at 1st atmospheric level  [kg/kg]
 365    REAL, INTENT(IN)    :: UA   ! wind speed at 1st atmospheric level    [m/s]
 366    REAL, INTENT(IN)    :: U1   ! u at 1st atmospheric level             [m/s]
 367    REAL, INTENT(IN)    :: V1   ! v at 1st atmospheric level             [m/s]
 368    REAL, INTENT(IN)    :: SSG  ! downward total short wave radiation    [W/m/m]
 369    REAL, INTENT(IN)    :: LLG  ! downward long wave radiation           [W/m/m]
 370    REAL, INTENT(IN)    :: RAIN ! precipitation                          [mm/h]
 371    REAL, INTENT(IN)    :: RHOO ! air density                            [kg/m^3]
 372    REAL, INTENT(IN)    :: ZA   ! first atmospheric level                [m]
 373    REAL, INTENT(IN)    :: DECLIN ! solar declination                    [rad]
 374    REAL, INTENT(IN)    :: COSZ ! sin(fai)*sin(del)+cos(fai)*cos(del)*cos(omg)
 375    REAL, INTENT(IN)    :: OMG  ! solar hour angle                       [rad]
 376
 377    REAL, INTENT(IN)    :: XLAT ! latitude                               [deg]
 378    REAL, INTENT(IN)    :: DELT ! time step                              [s]
 379    REAL, INTENT(IN)    :: CHS,CHS2 ! CH*U at za and 2 m             [m/s]
 380
 381    REAL, INTENT(INOUT) :: SSGD ! downward direct short wave radiation   [W/m/m]
 382    REAL, INTENT(INOUT) :: SSGQ ! downward diffuse short wave radiation  [W/m/m]
 383    REAL, INTENT(INOUT) :: CMR_URB
 384    REAL, INTENT(INOUT) :: CHR_URB
 385    REAL, INTENT(INOUT) :: CMC_URB
 386    REAL, INTENT(INOUT) :: CHC_URB
 387    REAL, INTENT(INOUT) :: ZNT  ! roughness length                       [m]    ! modified by danli
 388 !-------------------------------------------------------------------------------
 389 ! I: NUDAPT Input Parameters
 390 !-------------------------------------------------------------------------------
 391    REAL, INTENT(INOUT) :: mh_urb   ! mean building height [m]
 392    REAL, INTENT(INOUT) :: stdh_urb ! standard deviation of building height [m]
 393    REAL, INTENT(INOUT) :: hgt_urb  ! area weighted mean building height [m]
 394    REAL, INTENT(INOUT) :: lp_urb   ! plan area fraction [-]
 395    REAL, INTENT(INOUT) :: frc_urb  ! urban fraction [-]
 396    REAL, INTENT(INOUT) :: lb_urb   ! building surface to plan area ratio [-]
 397    REAL, INTENT(INOUT), DIMENSION(4) :: lf_urb   ! frontal area index [-]
 398    REAL, INTENT(INOUT) :: zo_check  ! check for printing ZOC
 399
 400 !-------------------------------------------------------------------------------
 401 ! O: output variables from Urban to LSM
 402 !-------------------------------------------------------------------------------
 403
 404    REAL, INTENT(OUT) :: TS     ! surface potential temperature    [K]
 405    REAL, INTENT(OUT) :: QS     ! surface humidity                 [K]
 406    REAL, INTENT(OUT) :: SH     ! sensible heat flux               [W/m/m]
 407    REAL, INTENT(OUT) :: LH     ! latent heat flux                 [W/m/m]
 408    REAL, INTENT(OUT) :: LH_KINEMATIC ! latent heat, kinetic     [kg/m/m/s]
 409    REAL, INTENT(OUT) :: SW     ! upward short wave radiation flux [W/m/m]
 410    REAL, INTENT(OUT) :: ALB    ! time-varying albedo            [fraction]
 411    REAL, INTENT(OUT) :: LW     ! upward long wave radiation flux  [W/m/m]
 412    REAL, INTENT(OUT) :: G      ! heat flux into the ground        [W/m/m]
 413    REAL, INTENT(OUT) :: RN     ! net radition                     [W/m/m]
 414    REAL, INTENT(OUT) :: PSIM   ! similality stability shear function for momentum
 415    REAL, INTENT(OUT) :: PSIH   ! similality stability shear function for heat
 416    REAL, INTENT(OUT) :: GZ1OZ0
 417    REAL, INTENT(OUT) :: U10    ! u at 10m                         [m/s]
 418    REAL, INTENT(OUT) :: V10    ! u at 10m                         [m/s]
 419    REAL, INTENT(OUT) :: TH2    ! potential temperature at 2 m     [K]
 420    REAL, INTENT(OUT) :: Q2     ! humidity at 2 m                  [-]
 421 !m   REAL, INTENT(OUT) :: CHS,CHS2 ! CH*U at za and 2 m             [m/s]
 422    REAL, INTENT(OUT) :: UST    ! friction velocity                [m/s]
 423
 424
 425 !-------------------------------------------------------------------------------
 426 ! H: Historical (state) variables of Urban : LSM <--> Urban
 427 !-------------------------------------------------------------------------------
 428
 429 ! TR: roof temperature              [K];      TRP: at previous time step [K]
 430 ! TB: building wall temperature     [K];      TBP: at previous time step [K]
 431 ! TG: road temperature              [K];      TGP: at previous time step [K]
 432 ! TC: urban-canopy air temperature  [K];      TCP: at previous time step [K]
 433 !                                                  (absolute temperature)
 434 ! QC: urban-canopy air mixing ratio [kg/kg];  QCP: at previous time step [kg/kg]
 435 !
 436 ! XXXR: Monin-Obkhov length for roof          [dimensionless]
 437 ! XXXB: Monin-Obkhov length for building wall [dimensionless]
 438 ! XXXG: Monin-Obkhov length for road          [dimensionless]
 439 ! XXXC: Monin-Obkhov length for urban-canopy  [dimensionless]
 440 !
 441 ! TRL, TBL, TGL: layer temperature [K] (absolute temperature)
 442
 443    REAL, INTENT(INOUT):: TR, TB, TG, TC, QC, UC
 444    REAL, INTENT(INOUT):: XXXR, XXXB, XXXG, XXXC
 445
 446    REAL, DIMENSION(1:num_roof_layers), INTENT(INOUT) :: TRL
 447    REAL, DIMENSION(1:num_wall_layers), INTENT(INOUT) :: TBL
 448    REAL, DIMENSION(1:num_road_layers), INTENT(INOUT) :: TGL
 449
 450 !===Yang,2014/10/08, urban hydrological variables for single layer UCM===
 451 ! FLXHUMR: evaporation over roof  [m/s];     FLXHUMRP: at previous time step [m/s]
 452 ! FLXHUMB: evaporation over wall  [m/s];     FLXHUMBP: at previous time step [m/s]
 453 ! FLXHUMG: evaporation over road  [m/s];     FLXHUMGP: at previous time step [m/s]
 454
 455 ! DRELR: water retention depth on roof [m];  DRELRP: at previous time stp [m]
 456 ! DRELB: water retention depth on wall [m];  DRELBP: at previous time stp [m]
 457 ! DRELG: water retention depth on road [m];  DRELGP: at previous time stp [m]
 458
 459 ! TGR: green roof surface temperature  [K];      TGRP: at previous time step [K]
 460 ! CMCR: Canopy intercepted water on green roof;  CMCRP: at previous time step
 461 ! SMR: soil moisture at each layer on roof [-];  SMRP: at previous time step
 462 ! TGRL:layer temperature on green roof [K]
 463
 464    REAL, INTENT(INOUT):: FLXHUMR,FLXHUMB,FLXHUMG,DRELR,DRELB,DRELG
 465    REAL, INTENT(INOUT):: TGR,CMCR,CHGR_URB,CMGR_URB
 466    REAL, DIMENSION(1:num_roof_layers), INTENT(INOUT) :: SMR
 467    REAL, DIMENSION(1:num_roof_layers), INTENT(INOUT) :: TGRL
 468
 469 !-------------------------------------------------------------------------------
 470 ! L:  Local variables from read_param
 471 !-------------------------------------------------------------------------------
 472
 473    REAL :: ZR, Z0C, Z0HC, ZDC, SVF, R, RW, HGT, AH, ALH
 474    REAL :: SIGMA_ZED
 475    REAL :: CAPR, CAPB, CAPG, AKSR, AKSB, AKSG, ALBR, ALBB, ALBG
 476    REAL :: EPSR, EPSB, EPSG, Z0R,  Z0B,  Z0G,  Z0HB, Z0HG
 477    REAL :: TRLEND,TBLEND,TGLEND
 478    REAL :: T1VR, T1VC,TH2V
 479    REAL :: RLMO_URB
 480    REAL :: AKANDA_URBAN
 481
 482    REAL :: TH2X                                                !m
 483
 484    INTEGER :: BOUNDR, BOUNDB, BOUNDG
 485    INTEGER :: CH_SCHEME, TS_SCHEME
 486
 487    LOGICAL :: SHADOW  ! [true=consider svf and shadow effects, false=consider svf effect only]
 488
 489 !for BEP
 490    INTEGER                        :: NUMDIR
 491    REAL,    DIMENSION ( MAXDIRS ) :: STREET_DIRECTION
 492    REAL,    DIMENSION ( MAXDIRS ) :: STREET_WIDTH
 493    REAL,    DIMENSION ( MAXDIRS ) :: BUILDING_WIDTH
 494    INTEGER                        :: NUMHGT
 495    REAL,    DIMENSION ( MAXHGTS ) :: HEIGHT_BIN
 496    REAL,    DIMENSION ( MAXHGTS ) :: HPERCENT_BIN
 497 !end BEP
 498 !-------------------------------------------------------------------------------
 499 ! L: Local variables
 500 !-------------------------------------------------------------------------------
 501
 502    REAL :: BETR, BETB, BETG
 503    REAL :: SX, SD, SQ, RX
 504    REAL :: UR, ZC, XLB, BB
 505    REAL :: Z, RIBB, RIBG, RIBC, BHR, BHB, BHG, BHC
 506    REAL :: TSC, LNET, SNET, FLXUV, THG, FLXTH, FLXHUM, FLXG
 507    REAL :: W, VFGS, VFGW, VFWG, VFWS, VFWW
 508    REAL :: HOUI1, HOUI2, HOUI3, HOUI4, HOUI5, HOUI6, HOUI7, HOUI8
 509    REAL :: SLX, SLX1, SLX2, SLX3, SLX4, SLX5, SLX6, SLX7, SLX8
 510    REAL :: FLXTHR, FLXTHB, FLXTHG
 511    REAL :: SR, SB, SG, RR, RB, RG
 512    REAL :: SR1, SR2, SB1, SB2, SG1, SG2, RR1, RR2, RB1, RB2, RG1, RG2
 513    REAL :: HR, HB, HG, ELER, ELEB, ELEG, G0R, G0B, G0G
 514    REAL :: ALPHAC, ALPHAR, ALPHAB, ALPHAG
 515    REAL :: CHC, CHR, CHB, CHG, CDC, CDR, CDB, CDG, CDGR
 516    REAL :: C1R, C1B, C1G, TE, TC1, TC2, QC1, QC2, QS0R, QS0B, QS0G,RHO,ES
 517
 518    REAL :: DESDT
 519    REAL :: F
 520    REAL :: DQS0RDTR
 521    REAL :: DRRDTR, DHRDTR, DELERDTR, DG0RDTR
 522    REAL :: DTR, DFDT
 523    REAL :: FX, FY, GF, GX, GY
 524    REAL :: DTCDTB, DTCDTG
 525    REAL :: DQCDTB, DQCDTG
 526    REAL :: DRBDTB1,  DRBDTG1,  DRBDTB2,  DRBDTG2
 527    REAL :: DRGDTB1,  DRGDTG1,  DRGDTB2,  DRGDTG2
 528    REAL :: DRBDTB,   DRBDTG,   DRGDTB,   DRGDTG
 529    REAL :: DHBDTB,   DHBDTG,   DHGDTB,   DHGDTG
 530    REAL :: DELEBDTB, DELEBDTG, DELEGDTG, DELEGDTB
 531    REAL :: DG0BDTB,  DG0BDTG,  DG0GDTG,  DG0GDTB
 532    REAL :: DQS0BDTB, DQS0GDTG
 533    REAL :: DTB, DTG, DTC
 534
 535    REAL :: THEATAZ    ! Solar Zenith Angle [rad]
 536    REAL :: THEATAS    ! = PI/2. - THETAZ
 537    REAL :: FAI        ! Latitude [rad]
 538    REAL :: CNT,SNT
 539    REAL :: PS         ! Surface Pressure [hPa]
 540    REAL :: TAV        ! Vertial Temperature [K]
 541
 542    REAL :: XXX, X, Z0, Z0H, CD, CH
 543    REAL :: XXX2, PSIM2, PSIH2, XXX10, PSIM10, PSIH10
 544    REAL :: PSIX, PSIT, PSIX2, PSIT2, PSIX10, PSIT10
 545
 546    REAL :: TRP, TBP, TGP, TCP, QCP, TST, QST
 547
 548    REAL :: WDR,HGT2,BW,DHGT
 549    REAL, parameter :: VonK = 0.4
 550    REAL :: lambda_f,alpha_macd,beta_macd,lambda_fr
 551
 552    INTEGER :: iteration, K, NUDAPT
 553    INTEGER :: tloc, tloc2, Kalh
 554
 555 !===Yang,2014/10/08, urban hydrological variables for single layer UCM===
 556    REAL :: FLXHUMRP, FLXHUMBP, FLXHUMGP
 557    REAL :: DRELRP, DRELBP, DRELGP
 558    REAL :: TGRP, CMCRP
 559    REAL, DIMENSION(1:num_roof_layers) :: ZSOILR, ETR, SMRP
 560 !===Define parameters for green roof===
 561    INTEGER :: KZ
 562    REAL :: RUNOFF1, RUNOFF2, RUNOFF3
 563    REAL :: SGR, SGR1, T1VGR, CHGR, ALPHAGR
 564    REAL :: FLXTHGR, FLXHUMGR, HGR, ELEGR, G0GR
 565    REAL :: QS0GR, EPGR, EDIR, ETTR, FV, DTGR, DRIP
 566 !   REAL :: DQS0GRDTGR, ETR, ECR,RAIN1, RAINDR, DEW, ETAR, BETGR
 567    REAL :: DQS0GRDTGR, ECR,RAIN1, RAINDR, DEW, ETAR, BETGR
 568 !   REAL :: DF1, RGR, RGRR, RCH, RR1, RR2, YY, ZZ1, SSOILR
 569    REAL :: DF1, RGR, RGRR, RCH, YY, ZZ1, SSOILR
 570    REAL :: DRRDTGR, DHRDTGR, DELERDTGR, DG0RDTGR, DFDVT
 571    real,parameter  :: SHDFAC   = 0.80   ! Vegetated area fraction of green roof vegetation
 572    real,parameter  :: ALBV     = 0.20   ! green roof albedo
 573    real,parameter  :: EPSV     = 0.93   ! green roof emissivity
 574    real,parameter  :: LAI      = 1.50   ! leaf area index on green roof
 575    real,parameter  :: CMCMAX   = 0.5E-3 ! Maximum canopy interception capacity
 576    real,parameter  :: SMCREF   = 0.329  ! Reference soil moisture
 577    real,parameter  :: SMCDRY   = 0.066  ! Residual  soil moisture
 578    real,parameter  :: SMCWLT   = 0.084  ! Wilting point
 579    real,parameter  :: SMCMAX   = 0.439  ! Saturated soil moisture
 580    real,parameter  :: RSMAX    = 5000   ! Maximum stomatal resistance
 581    real,parameter  :: RSMIN    = 100    ! Minimum stomatal resistance
 582    real,parameter  :: RGL      = 100    ! Radiation limit where photosynthesis begins
 583    real,parameter  :: CFACTR   = 0.5    ! Parameter used in the canopy inteception calculation
 584    real,parameter  :: DWSAT    = 0.143E-4 ! Saturated soil conductivity
 585    real,parameter  :: DKSAT    = 3.38E-6  ! Saturated soil diffusivity
 586    real,parameter  :: BEXP     = 5.25   ! B parameter in soil hydraulic calculation
 587    real,parameter  :: FXEXP    = 2.0    ! Parameter for computing direct soil evaporation
 588    real,parameter  :: ZBOT     = -2.0
 589    real,parameter  :: QUARTZ   = 0.40
 590    real,parameter  :: CSOIL    = 2.0E+6
 591    real,parameter  :: HS       = 36
 592    integer,parameter ::  NROOT = 2      ! Root depth layer of green roof
 593    integer,parameter ::  NGR   = 4      ! Layer of green roof
 594    integer,parameter ::  IMPR  = 1
 595    integer,parameter ::  IMPB  = 2
 596    integer,parameter ::  IMPG  = 3
 597
 598 !-------------------------------------------------------------------------------
 599 ! Set parameters
 600 !-------------------------------------------------------------------------------
 601
 602 ! Miao, 2007/01/17, cal. ah
 603    if(ahoption==1) then
 604      tloc=mod(int(OMG/PI*180./15.+12.+0.5 ),24)
 605      if(tloc.lt.0) tloc=tloc+24
 606      if(tloc==0) tloc=24
 607    endif
 608 ! Yang, 2014/10/08, cal. alh
 609    if(alhoption==1) then
 610      tloc2=mod(int((OMG/PI*180./15.+12.)*2.+0.5 ),48)
 611      if(tloc2.lt.0) tloc2=tloc2+48
 612      if(tloc2==0) tloc2=48
 613    endif
 614
 615    CALL read_param(UTYPE,ZR,SIGMA_ZED,Z0C,Z0HC,ZDC,SVF,R,RW,HGT,  &
 616                    AH,CAPR,CAPB,CAPG,AKSR,AKSB,AKSG,ALBR,ALBB,    &
 617                    ALBG,EPSR,EPSB,EPSG,Z0R,Z0B,Z0G,Z0HB,Z0HG,     &
 618                    BETR,BETB,BETG,TRLEND,TBLEND,TGLEND,           &
 619 !for BEP
 620                    NUMDIR, STREET_DIRECTION, STREET_WIDTH,        &
 621                    BUILDING_WIDTH, NUMHGT, HEIGHT_BIN,            &
 622                    HPERCENT_BIN,                                  &
 623 !end BEP
 624                    BOUNDR,BOUNDB,BOUNDG,CH_SCHEME,TS_SCHEME,      &
 625                    AKANDA_URBAN,ALH)
 626
 627 ! Glotfelty, 2012/07/05, NUDAPT Modification
 628
 629  if(mh_urb.gt.0.0)THEN
 630   !write(mesg,*) 'Mean Height NUDAPT',mh_urb
 631   !WRITE_MESSAGE(mesg)
 632   !write(mesg,*) 'Mean Height Table',ZR
 633   !WRITE_MESSAGE(mesg)
 634   if(zo_check.eq.1)THEN
 635    write(mesg,*) 'Mean Height NUDAPT',mh_urb
 636    WRITE_MESSAGE(mesg)
 637    write(mesg,*) 'Mean Height Table',ZR
 638    WRITE_MESSAGE(mesg)
 639    write(mesg,*) 'Roughness Length Table',Z0C
 640    WRITE_MESSAGE(mesg)
 641    write(mesg,*) 'Roof Roughness Length Table',Z0R
 642    WRITE_MESSAGE(mesg)
 643    write(mesg,*) 'Sky View Factor Table',SVF
 644    WRITE_MESSAGE(mesg)
 645    write(mesg,*) 'Normalized Height Table',HGT
 646    WRITE_MESSAGE(mesg)
 647    write(mesg,*) 'Plan Area Fraction', lp_urb
 648    WRITE_MESSAGE(mesg)
 649    write(mesg,*) 'Plan Area Fraction table', R
 650    WRITE_MESSAGE(mesg)
 651   end if
 652   !write(mesg,*) 'Area Weighted Mean Height',hgt_urb
 653   !WRITE_MESSAGE(mesg)
 654   !write(mesg,*) 'Plan Area Fraction', lp_urb
 655   !WRITE_MESSAGE(mesg)
 656   !write(mesg,*) 'STD Height', stdh_urb
 657   !WRITE_MESSAGE(mesg)
 658   !write(mesg,*) 'Frontal Area Index',lf_urb
 659   !WRITE_MESSAGE(mesg)
 660   !write(mesg,*) 'Urban Fraction',frc_urb
 661   !WRITE_MESSAGE(mesg)
 662   !write(mesg,*) 'Building Surf Ratio',lb_urb
 663   !WRITE_MESSAGE(mesg)
 664
 665  !Calculate Building Width and Street Width Based on BEP formulation
 666  if(lb_urb.gt.lp_urb)THEN
 667   BW=2.*hgt_urb*lp_urb/(lb_urb-lp_urb)
 668   SW=2.*hgt_urb*lp_urb*((frc_urb/lp_urb)-1.)/(lb_urb-lp_urb)
 669   !write(mesg,*) 'Building Width',BW
 670   !WRITE_MESSAGE(mesg)
 671   !write(mesg,*) 'Street Width',SW
 672   !WRITE_MESSAGE(mesg)
 673  elseif (SW.lt.0.0.or.BW.lt.0.0)then
 674   BW=BUILDING_WIDTH(1)
 675   SW=STREET_WIDTH(1)
 676  else
 677     BW=BUILDING_WIDTH(1)
 678   SW=STREET_WIDTH(1)
 679  end if
 680
 681  !Assign NUDAPT Parameters
 682    ZR = mh_urb
 683     R = lp_urb
 684     RW = 1.0-R
 685    HGT = mh_urb/(BW+SW)
 686    SIGMA_ZED = stdh_urb
 687
 688  !Calculate Wind Direction and Assign Appropriae lf_urb
 689    WDR = (180.0/PI)*ATAN2(U10,V10)
 690
 691    IF(WDR.ge.0.0.and.WDR.lt.22.5)THEN
 692      lambda_f = lf_urb(1)
 693    ELSEIF(WDR.ge.-22.5.and.WDR.lt.0.0)THEN
 694      lambda_f = lf_urb(1)
 695    ELSEIF(WDR.gt.157.5.and.WDR.le.180.0)THEN
 696      lambda_f = lf_urb(1)
 697    ELSEIF(WDR.lt.-157.5)THEN
 698      lambda_f = lf_urb(1)
 699    ELSEIF(WDR.gt.22.5.and.WDR.le.67.5)THEN
 700      lambda_f = lf_urb(2)
 701    ELSEIF(WDR.ge.-157.5.and.WDR.lt.-112.5)THEN
 702      lambda_f = lf_urb(2)
 703    ELSEIF(WDR.gt.67.5.and.WDR.le.112.5)THEN
 704      lambda_f = lf_urb(3)
 705    ELSEIF(WDR.ge.-112.5.and.WDR.lt.-67.5)THEN
 706      lambda_f = lf_urb(3)
 707    ELSEIF(WDR.gt.112.5.and.WDR.le.157.5)THEN
 708      lambda_f = lf_urb(4)
 709    ELSEIF(WDR.ge.-67.5.and.WDR.lt.-22.5)THEN
 710      lambda_f = lf_urb(4)
 711    ELSE
 712      lambda_f = lf_urb(1)
 713    ENDIF
 714
 715   !Calculate the following urban canyon geometry parameters following Macdonald's (1998) formulations
 716       Cd         = 1.2
 717       alpha_macd = 4.43
 718       beta_macd  = 1.0
 719
 720
 721       ZDC = ZR * ( 1.0 + ( alpha_macd ** ( -R ) )  * ( R - 1.0 ) )
 722
 723       Z0C = ZR * ( 1.0 - ZDC/ZR ) * &
 724            exp (-(0.5 * beta_macd * Cd / (VonK**2) * ( 1.0-ZDC/ZR) * lambda_f )**(-0.5))
 725
 726       if(zo_check.eq.1)THEN
 727         write(mesg,*) 'Roughness Length NUDAPT',Z0C
 728         WRITE_MESSAGE(mesg)
 729       end if
 730
 731       lambda_fr  = stdh_urb/(SW + BW)
 732
 733        Z0R = ZR * ( 1.0 - ZDC/ZR) &
 734              * exp ( -(0.5 * beta_macd * Cd / (VonK**2) &
 735               * ( 1.0-ZDC/ZR) * lambda_fr )**(-0.5))
 736
 737
 738
 739        Z0HC = 0.1 * Z0C
 740
 741   ! Calculate Sky View Factor
 742
 743      DHGT=HGT/100.
 744      HGT2=0.
 745      VFWS=0.
 746      HGT2=HGT-DHGT/2.
 747       do NUDAPT=1,99
 748          HGT2=HGT2-DHGT
 749          VFWS=VFWS+0.25*(1.-HGT2/SQRT(HGT2**2.+RW**2.))
 750       end do
 751
 752      VFWS=VFWS/99.
 753      VFWS=VFWS*2.
 754
 755      VFGS=1.-2.*VFWS*HGT/RW
 756      SVF=VFGS
 757
 758     if(zo_check.eq.1)THEN
 759       write(mesg,*) 'Roof Roughness Length NUDAPT',Z0R
 760       WRITE_MESSAGE(mesg)
 761       write(mesg,*) 'Sky View Factor NUDAPT',SVF
 762       WRITE_MESSAGE(mesg)
 763       write(mesg,*) 'normalized Height NUDAPT', HGT
 764       WRITE_MESSAGE(mesg)
 765     end if
 766
 767
 768  endif
 769
 770  !End NUDAPT Modification
 771
 772
 773 ! Miao, 2007/01/17, cal. ah
 774    if(ahoption==1) AH=AH*ahdiuprf(tloc)
 775
 776 ! Yang, 2014/10/08, cal. alh
 777    Kalh=0
 778    if(alhoption==1) THEN
 779      if(jmonth==3 .or. jmonth==4 .or. jmonth==5) Kalh=1
 780      if(jmonth==6 .or. jmonth==7 .or. jmonth==8) Kalh=2
 781      if(jmonth==9 .or. jmonth==10.or. jmonth==11)Kalh=3
 782      if(jmonth==12.or. jmonth==1 .or. jmonth==2) Kalh=4
 783    endif
 784    if(alhoption==1) ALH = ALH*alhdiuprf(tloc2)*alhseason(Kalh)
 785
 786    IF( ZDC+Z0C+2. >= ZA) THEN
 787     FATAL_ERROR("ZDC + Z0C + 2m is larger than the 1st WRF level - Stop in subroutine urban - change ZDC and Z0C" )
 788    END IF
 789
 790    IF(.NOT.LSOLAR) THEN
 791      SSGD = SRATIO*SSG
 792      SSGQ = SSG - SSGD
 793    ENDIF
 794    SSGD = SRATIO*SSG   ! No radiation scheme has SSGD and SSGQ.
 795    SSGQ = SSG - SSGD
 796
 797    W=2.*1.*HGT
 798    VFGS=SVF
 799    VFGW=1.-SVF
 800    VFWG=(1.-SVF)*(1.-R)/W
 801    VFWS=VFWG
 802    VFWW=1.-2.*VFWG
 803
 804 !-------------------------------------------------------------------------------
 805 ! Convert unit from MKS to cgs
 806 ! Renew surface and layer temperatures
 807 !-------------------------------------------------------------------------------
 808
 809    SX=(SSGD+SSGQ)/697.7/60.  ! downward short wave radition [ly/min]
 810    SD=SSGD/697.7/60.         ! downward direct short wave radiation
 811    SQ=SSGQ/697.7/60.         ! downward diffuse short wave radiation
 812    RX=LLG/697.7/60.          ! downward long wave radiation
 813    RHO=RHOO*0.001            ! air density at first atmospheric level
 814
 815    TRP=TR
 816    TBP=TB
 817    TGP=TG
 818    TCP=TC
 819    QCP=QC
 820
 821 !===Yang,2014/10/08, urban hydrological variables for single layer UCM===
 822    FLXHUMRP = FLXHUMR
 823    FLXHUMBP = FLXHUMB
 824    FLXHUMGP = FLXHUMG
 825    DRELRP = DRELR
 826    DRELBP = DRELB
 827    DRELGP = DRELG
 828    TGRP   = TGR
 829    CMCRP  = CMCR
 830    SMRP   = SMR
 831
 832 !===Yang,2014/10/08, urban irrigation, May-Sep, 9-10pm
 833    IF(IRI_SCHEME==1) THEN
 834        IF (tloc==21 .or. tloc==22) THEN
 835          IF(jmonth==5 .or. jmonth==6 .or. jmonth ==7 .or. &
 836           jmonth==8 .or. jmonth==9 ) THEN
 837             DO KZ = 1,2
 838                SMRP(KZ)= SMCREF
 839              END DO
 840        ENDIF
 841      ENDIF
 842    ENDIF
 843
 844    TAV=TA*(1.+0.61*QA)
 845    PS=RHOO*287.*TAV/100. ![hPa]
 846
 847 !-------------------------------------------------------------------------------
 848 ! Canopy wind
 849 !-------------------------------------------------------------------------------
 850
 851    IF ( ZR + 2. < ZA ) THEN
 852      UR=UA*LOG((ZR-ZDC)/Z0C)/LOG((ZA-ZDC)/Z0C)
 853      ZC=0.7*ZR
 854      XLB=0.4*(ZR-ZDC)
 855      ! BB formulation from Inoue (1963)
 856      BB = 0.4 * ZR / ( XLB * alog( ( ZR - ZDC ) / Z0C ) )
 857      UC=UR*EXP(-BB*(1.-ZC/ZR))
 858    ELSE
 859      ! PRINT *, 'Warning ZR + 2m  is larger than the 1st WRF level'
 860      ZC=ZA/2.
 861      UC=UA/2.
 862    END IF
 863
 864 !-------------------------------------------------------------------------------
 865 ! Net Short Wave Radiation at roof, wall, and road
 866 !-------------------------------------------------------------------------------
 867
 868    SHADOW = .false.
 869 !   SHADOW = .true.
 870
 871    IF (SSG > 0.0) THEN
 872
 873      IF(.NOT.SHADOW) THEN              ! no shadow effects model
 874
 875        SR1=SX*(1.-ALBR)
 876        SGR1=SX*(1.-ALBV)
 877        SG1=SX*VFGS*(1.-ALBG)
 878        SB1=SX*VFWS*(1.-ALBB)
 879        SG2=SB1*ALBB/(1.-ALBB)*VFGW*(1.-ALBG)
 880        SB2=SG1*ALBG/(1.-ALBG)*VFWG*(1.-ALBB)
 881
 882      ELSE                              ! shadow effects model
 883
 884        FAI=XLAT*PI/180.
 885
 886        THEATAS=ABS(ASIN(COSZ))
 887        THEATAZ=ABS(ACOS(COSZ))
 888
 889        SNT=COS(DECLIN)*SIN(OMG)/COS(THEATAS)
 890        CNT=(COSZ*SIN(FAI)-SIN(DECLIN))/COS(THEATAS)/COS(FAI)
 891
 892        HOUI1=(SNT*COS(PI/8.)    -CNT*SIN(PI/8.))
 893        HOUI2=(SNT*COS(2.*PI/8.) -CNT*SIN(2.*PI/8.))
 894        HOUI3=(SNT*COS(3.*PI/8.) -CNT*SIN(3.*PI/8.))
 895        HOUI4=(SNT*COS(4.*PI/8.) -CNT*SIN(4.*PI/8.))
 896        HOUI5=(SNT*COS(5.*PI/8.) -CNT*SIN(5.*PI/8.))
 897        HOUI6=(SNT*COS(6.*PI/8.) -CNT*SIN(6.*PI/8.))
 898        HOUI7=(SNT*COS(7.*PI/8.) -CNT*SIN(7.*PI/8.))
 899        HOUI8=(SNT*COS(8.*PI/8.) -CNT*SIN(8.*PI/8.))
 900
 901        SLX1=HGT*ABS(TAN(THEATAZ))*ABS(HOUI1)
 902        SLX2=HGT*ABS(TAN(THEATAZ))*ABS(HOUI2)
 903        SLX3=HGT*ABS(TAN(THEATAZ))*ABS(HOUI3)
 904        SLX4=HGT*ABS(TAN(THEATAZ))*ABS(HOUI4)
 905        SLX5=HGT*ABS(TAN(THEATAZ))*ABS(HOUI5)
 906        SLX6=HGT*ABS(TAN(THEATAZ))*ABS(HOUI6)
 907        SLX7=HGT*ABS(TAN(THEATAZ))*ABS(HOUI7)
 908        SLX8=HGT*ABS(TAN(THEATAZ))*ABS(HOUI8)
 909
 910        IF(SLX1 > RW) SLX1=RW
 911        IF(SLX2 > RW) SLX2=RW
 912        IF(SLX3 > RW) SLX3=RW
 913        IF(SLX4 > RW) SLX4=RW
 914        IF(SLX5 > RW) SLX5=RW
 915        IF(SLX6 > RW) SLX6=RW
 916        IF(SLX7 > RW) SLX7=RW
 917        IF(SLX8 > RW) SLX8=RW
 918
 919        SLX=(SLX1+SLX2+SLX3+SLX4+SLX5+SLX6+SLX7+SLX8)/8.
 920
 921        SR1=SD*(1.-ALBR)+SQ*(1.-ALBR)
 922        SGR1=SD*(1.-ALBV)+SQ*(1.-ALBV)
 923        SG1=SD*(RW-SLX)/RW*(1.-ALBG)+SQ*VFGS*(1.-ALBG)
 924        SB1=SD*SLX/W*(1.-ALBB)+SQ*VFWS*(1.-ALBB)
 925        SG2=SB1*ALBB/(1.-ALBB)*VFGW*(1.-ALBG)
 926        SB2=SG1*ALBG/(1.-ALBG)*VFWG*(1.-ALBB)
 927
 928      END IF
 929
 930      SR=SR1
 931      SGR=SGR1
 932      SG=SG1+SG2
 933      SB=SB1+SB2
 934
 935      IF (GROPTION ==1) THEN
 936      SNET=R*FGR*SGR+R*(1.-FGR)*SR+W*SB+RW*SG
 937      ELSE
 938      SNET=R*SR+W*SB+RW*SG
 939      ENDIF
 940
 941    ELSE
 942
 943      SR=0.
 944      SG=0.
 945      SGR=0.
 946      SB=0.
 947      SNET=0.
 948
 949    END IF
 950
 951 !-------------------------------------------------------------------------------
 952 ! Roof
 953 !-------------------------------------------------------------------------------
 954
 955 !-------------------------------------------------------------------------------
 956 ! CHR, CDR, BETR
 957 !-------------------------------------------------------------------------------
 958
 959    ! Z=ZA-ZDC
 960    ! BHR=LOG(Z0R/Z0HR)/0.4
 961    ! RIBR=(9.8*2./(TA+TRP))*(TA-TRP)*(Z+Z0R)/(UA*UA)
 962    ! CALL mos(XXXR,ALPHAR,CDR,BHR,RIBR,Z,Z0R,UA,TA,TRP,RHO)
 963
 964    ! Alternative option for MOST using SFCDIF routine from Noah
 965    ! Virtual temperatures needed by SFCDIF
 966    T1VR = TRP* (1.0+ 0.61 * QA)
 967    TH2V = (TA + ( 0.0098 * ZA)) * (1.0+ 0.61 * QA)
 968
 969    ! note that CHR_URB contains UA (=CHR_MOS*UA)
 970    RLMO_URB=0.0
 971    CALL SFCDIF_URB (ZA,Z0R,T1VR,TH2V,UA,AKANDA_URBAN,CMR_URB,CHR_URB,RLMO_URB,CDR)
 972    ALPHAR =  RHO*CP*CHR_URB
 973    CHR=ALPHAR/RHO/CP/UA
 974
 975 ! Yang, 03/12/2014 -- LH for impervious roof surface
 976    RAIN1 = RAIN * 0.001 /3600          ! CONVERT FROM mm/hr to m/s
 977    IF (IMP_SCHEME==1) then
 978    IF (RAIN > 1.) BETR=0.7
 979    ENDIF
 980
 981    IF (IMP_SCHEME==2) then
 982    IF (FLXHUMRP <= 0.) FLXHUMRP = 0.
 983 !  Compute water retention depth from previous time step
 984    DrelR = DrelRP+(RAIN1-FLXHUMRP)*DELT/porimp(IMPR)
 985    IF (RAIN > 0. .AND. DrelR < DrelRP) DrelR = DrelRP
 986
 987    IF (DrelR <= 0.) then
 988    DrelR  = 0.0
 989    BETR   = 0.0
 990    ELSEIf (DrelR <= dengimp(IMPR)) then
 991    BETR = DrelR/dengimp(IMPR)*porimp(IMPR)
 992    ELSE
 993    DrelR = dengimp(IMPR)
 994    BETR  = porimp(IMPR)
 995    ENDIF
 996
 997    IF ( BETR < 1.E-5 ) BETR = 0.0
 998    ENDIF
 999
1000    IF (TS_SCHEME == 1) THEN
1001
1002 !-------------------------------------------------------------------------------
1003 ! TR  Solving Non-Linear Equation by Newton-Rapson
1004 ! TRL Solving Heat Equation by Tri Diagonal Matrix Algorithm
1005 !-------------------------------------------------------------------------------
1006 !      TSC=TRP-273.15
1007 !      ES=EXP(19.482-4303.4/(TSC+243.5))        ! WMO
1008 !      ES=6.11*10.**(TETENA*TSC/(TETENB+TSC))   ! Tetens
1009 !      DESDT=( 6.1078*(2500.-2.4*TSC)/ &        ! Tetens
1010 !            (0.46151*(TSC+273.15)**2.) )*10.**(7.5*TSC/(237.3+TSC))
1011 !      ES=6.11*EXP((2.5*10.**6./461.51)*(TRP-273.15)/(273.15*TRP) ) ! Clausius-Clapeyron
1012 !      DESDT=(2.5*10.**6./461.51)*ES/(TRP**2.)                      ! Clausius-Clapeyron
1013 !      QS0R=0.622*ES/(PS-0.378*ES)
1014 !      DQS0RDTR = DESDT*0.622*PS/((PS-0.378*ES)**2.)
1015 !      DQS0RDTR = 17.269*(273.15-35.86)/((TRP-35.86)**2.)*QS0R
1016
1017 !    TRP=350.
1018
1019      DO ITERATION=1,20
1020
1021        ES=6.11*EXP( (2.5*10.**6./461.51)*(TRP-273.15)/(273.15*TRP) )
1022        DESDT=(2.5*10.**6./461.51)*ES/(TRP**2.)
1023        QS0R=0.622*ES/(PS-0.378*ES)
1024        DQS0RDTR = DESDT*0.622*PS/((PS-0.378*ES)**2.)
1025
1026        RR=EPSR*(RX-SIG*(TRP**4.)/60.)
1027        HR=RHO*CP*CHR*UA*(TRP-TA)*100.
1028        ELER=RHO*EL*CHR*UA*BETR*(QS0R-QA)*100.
1029        G0R=AKSR*(TRP-TRL(1))/(DZR(1)/2.)
1030
1031        F = SR + RR - HR - ELER - G0R
1032
1033        DRRDTR = (-4.*EPSR*SIG*TRP**3.)/60.
1034        DHRDTR = RHO*CP*CHR*UA*100.
1035        DELERDTR = RHO*EL*CHR*UA*BETR*DQS0RDTR*100.
1036        DG0RDTR =  2.*AKSR/DZR(1)
1037
1038        DFDT = DRRDTR - DHRDTR - DELERDTR - DG0RDTR
1039        DTR = F/DFDT
1040
1041        TR = TRP - DTR
1042        TRP = TR
1043
1044        IF( ABS(F) < 0.000001 .AND. ABS(DTR) < 0.000001 ) EXIT
1045
1046      END DO
1047
1048 ! multi-layer heat equation model
1049
1050      CALL multi_layer(num_roof_layers,BOUNDR,G0R,CAPR,AKSR,TRL,DZR,DELT,TRLEND)
1051
1052    ELSE
1053
1054      ES=6.11*EXP( (2.5*10.**6./461.51)*(TRP-273.15)/(273.15*TRP) )
1055      QS0R=0.622*ES/(PS-0.378*ES)
1056
1057      RR=EPSR*(RX-SIG*(TRP**4.)/60.)
1058      HR=RHO*CP*CHR*UA*(TRP-TA)*100.
1059      ELER=RHO*EL*CHR*UA*BETR*(QS0R-QA)*100.
1060      G0R=SR+RR-HR-ELER
1061
1062      CALL force_restore(CAPR,AKSR,DELT,SR,RR,HR,ELER,TRLEND,TRP,TR)
1063
1064      TRP=TR
1065
1066    END IF
1067
1068    FLXTHR=HR/RHO/CP/100.
1069    FLXHUMR=ELER/RHO/EL/100.
1070
1071 !-------------------------------------------------------------------------------
1072 ! Green Roof
1073 ! Must use multiple layers scheme (TS_SCHEME=1)
1074 !-------------------------------------------------------------------------------
1075    IF (GROPTION == 1) THEN
1076    T1VGR = TGRP* (1.0+ 0.61 * QA)
1077    RLMO_URB=0.0
1078    CALL SFCDIF_URB (ZA,Z0R,T1VGR,TH2V,UA,AKANDA_URBAN,CMGR_URB,CHGR_URB,RLMO_URB,CDGR)
1079    ALPHAGR =  RHO*CP*CHGR_URB
1080    CHGR=ALPHAGR/RHO/CP/UA
1081    RUNOFF1 = 0.0
1082    RUNOFF2 = 0.0
1083    RUNOFF3 = 0.0
1084
1085    KZ = 1
1086    ZSOILR (KZ) = - DZGR (KZ)
1087    DO KZ = 2,NGR
1088       ZSOILR (KZ) = - DZGR(KZ) + ZSOILR (KZ -1)
1089    END DO
1090
1091    DO ITERATION=1,100
1092        KZ=1
1093        ES=6.11*EXP( (2.5*10.**6./461.51)*(TGRP-273.15)/(273.15*TGRP) )
1094        DESDT=(2.5*10.**6./461.51)*ES/(TGRP**2.)
1095        QS0GR=0.622*ES/(PS-0.378*ES)
1096        DQS0GRDTGR = DESDT*0.622*PS/((PS-0.378*ES)**2.)
1097        EPGR=RHOO*CHGR*UA*(QS0GR-QA)   ! Potential evaporation [kg/m2/s]
1098
1099        IF (EPGR > 0.0) THEN
1100        ! Direct evaporation from soil on green roof
1101        CALL DIREVAP (EDIR,EPGR,SMRP(KZ),SHDFAC,SMCMAX,SMCDRY,FXEXP)
1102        ! Evapotranspiration and canopy intercepted evaporation
1103        CALL TRANSP  (ETTR,ETR,ECR,SHDFAC,EPGR,CMCRP,CFACTR,CMCMAX,LAI,RSMIN,RSMAX,RGL,SX, &
1104                      TGRP,TA,QA,SMRP,SMCWLT,SMCREF,CPP,PS,CHGR,EPSV,DELT,NROOT,NGR,DZGR, &
1105                      ZSOILR,HS)
1106        ! Update moisture in soil layers
1107        CALL SMFLX   (SMRP,SMR,NGR,CMCRP,CMCR,DELT,RAIN,ZSOILR,SMCMAX,BEXP,SMCWLT,DKSAT,&
1108                      DWSAT,SHDFAC,CMCMAX,RUNOFF1,RUNOFF2,RUNOFF3,EDIR,ECR,ETR,DRIP)
1109        else
1110        DEW  = - EPGR
1111        RAINDR = RAIN  + DEW * 3600.
1112        EDIR=0.0
1113        ECR =0.0
1114        ETTR=0.0
1115        CALL SMFLX (SMRP,SMR,NGR,CMCRP,CMCR,DELT,RAINDR,ZSOILR,SMCMAX,BEXP,SMCWLT,DKSAT,&
1116                    DWSAT,SHDFAC,CMCMAX,RUNOFF1,RUNOFF2,RUNOFF3,EDIR,ECR,ETR,DRIP)
1117        END IF
1118 ! ----------------------------------------------------------------------
1119 ! CONVERT MODELED EVAPOTRANSPIRATION FROM  M S-1  TO  KG M-2 S-1.
1120 ! ----------------------------------------------------------------------
1121        EDIR = EDIR  * 1000.0
1122        ETTR = ETTR  * 1000.0
1123        ECR  = ECR   * 1000.0
1124        ETAR = EDIR + ETTR + ECR
1125        IF (ETAR < 1.E-20) ETAR = 0.0
1126
1127        IF ( EPGR <= 0.0 ) THEN
1128          BETGR = 0.0
1129        ELSE
1130          BETGR = ETAR / EPGR
1131        END IF
1132        ELEGR= ETAR* RHO * EL /RHOO * 100
1133
1134        CALL TDFCND (DF1,SMR(KZ), QUARTZ, SMCMAX )
1135        DF1 = DF1 * EXP(-2.0 * SHDFAC)
1136        RGR = EPSV*(RX-SIG*(TGRP**4.)/60.)
1137        RGRR= (SGR+RGR) * 697.7 * 60.
1138        RCH = RHOO*CPP*CHGR
1139        RR1 = EPSV*(TA**4) * 6.48E-8 / (PS* CHGR) + 1.0
1140        IF (RAIN >  0.0) then
1141        RR2 = RR1 + RAIN / 3600 * 4.218E+3 / RCH
1142        else
1143        RR2 = RR1
1144        end if
1145        YY  = TA + (RGRR / RCH - BETGR * EPGR * ELL/ RCH) / RR2
1146        ZZ1 = DF1 / (-0.5 * ZSOILR (KZ) * RCH * RR2 ) + 1.0
1147
1148
1149        HGR=RHO*CP*CHGR*UA*(TGRP-TA)*100.
1150        RUNOFF3 = RUNOFF3/ DELT
1151        RUNOFF2 = RUNOFF2+ RUNOFF3
1152        G0GR    = DF1*(TGRP-TGRL(1))/(DZGR(1)/2.)/697.7/60
1153
1154        FV = SGR + RGR - HGR - ELEGR - G0GR
1155        DRRDTGR   = (-4.*EPSV*SIG*TGRP**3.)/60.
1156        DHRDTGR   = RHO*CP*CHGR*UA*100.
1157        DELERDTGR = RHO*EL*CHGR*UA*BETGR*DQS0GRDTGR*100.
1158        DG0RDTGR  = 2.*DF1/ DZGR(KZ) * ( 1.0 / 4.1868 ) * 1.E-4
1159        DFDVT = DRRDTGR - DHRDTGR - DELERDTGR - DG0RDTGR
1160        DTGR = FV/DFDVT/ 6
1161        TGR  = TGRP - DTGR
1162        TGRP = TGR
1163
1164        IF( ABS(FV) < 0.0001 .AND. ABS(DTGR) < 0.001 ) then
1165        EXIT
1166        ENDIF
1167      END DO
1168        ! Update temperature in soil layer
1169        CALL SHFLX (SSOILR,TGRL,SMR,SMCMAX,NGR,TGRP,DELT,YY,ZZ1,ZSOILR,       &
1170                   TRLEND,ZBOT,SMCWLT,DF1,QUARTZ,CSOIL,CAPR)
1171    FLXTHGR=HGR/RHO/CP/100.
1172    FLXHUMGR=ELEGR/RHO/EL/100.
1173 ELSE
1174    FLXTHGR=0.
1175    FLXHUMGR=0.
1176 ENDIF
1177
1178 !-------------------------------------------------------------------------------
1179 !  Wall and Road
1180 !-------------------------------------------------------------------------------
1181
1182 !-------------------------------------------------------------------------------
1183 !  CHC, CHB, CDB, BETB, CHG, CDG, BETG
1184 !-------------------------------------------------------------------------------
1185
1186    ! Z=ZA-ZDC
1187    ! BHC=LOG(Z0C/Z0HC)/0.4
1188    ! RIBC=(9.8*2./(TA+TCP))*(TA-TCP)*(Z+Z0C)/(UA*UA)
1189    !
1190    ! CALL mos(XXXC,ALPHAC,CDC,BHC,RIBC,Z,Z0C,UA,TA,TCP,RHO)
1191    ! Virtual temperatures needed by SFCDIF routine from Noah
1192
1193    T1VC = TCP* (1.0+ 0.61 * QA)
1194    RLMO_URB=0.0
1195    CALL SFCDIF_URB(ZA,Z0C,T1VC,TH2V,UA,AKANDA_URBAN,CMC_URB,CHC_URB,RLMO_URB,CDC)
1196    ALPHAC =  RHO*CP*CHC_URB
1197
1198    IF (CH_SCHEME == 1) THEN
1199
1200      Z=ZDC
1201      BHB=LOG(Z0B/Z0HB)/0.4
1202      BHG=LOG(Z0G/Z0HG)/0.4
1203      RIBB=(9.8*2./(TCP+TBP))*(TCP-TBP)*(Z+Z0B)/(UC*UC)
1204      RIBG=(9.8*2./(TCP+TGP))*(TCP-TGP)*(Z+Z0G)/(UC*UC)
1205
1206      CALL mos(XXXB,ALPHAB,CDB,BHB,RIBB,Z,Z0B,UC,TCP,TBP,RHO)
1207      CALL mos(XXXG,ALPHAG,CDG,BHG,RIBG,Z,Z0G,UC,TCP,TGP,RHO)
1208
1209    ELSE
1210
1211      ALPHAB=RHO*CP*(6.15+4.18*UC)/1200.
1212      IF(UC > 5.) ALPHAB=RHO*CP*(7.51*UC**0.78)/1200.
1213      ALPHAG=RHO*CP*(6.15+4.18*UC)/1200.
1214      IF(UC > 5.) ALPHAG=RHO*CP*(7.51*UC**0.78)/1200.
1215
1216    END IF
1217
1218    CHC=ALPHAC/RHO/CP/UA
1219    CHB=ALPHAB/RHO/CP/UC
1220    CHG=ALPHAG/RHO/CP/UC
1221
1222 !Yang 10/10/2013 -- LH from impervious wall and ground
1223  IF (IMP_SCHEME==1) then
1224    BETB=0.0
1225    IF(RAIN > 1.) BETG=0.7
1226  ENDIF
1227
1228  IF (IMP_SCHEME==2) then
1229    IF (FLXHUMBP <= 0.) FLXHUMBP = 0.
1230    IF (FLXHUMGP <= 0.) FLXHUMGP = 0.
1231 ! Compute water retention from previous time step for wall and ground
1232   DrelB = DrelBP+(RAIN1-FLXHUMBP)*DELT/porimp(IMPB)
1233     IF (RAIN > 0. .AND. DrelB < DrelBP) DrelB = DrelBP
1234   DrelG = DrelGP+(RAIN1-FLXHUMGP)*DELT/porimp(IMPG)
1235     IF (RAIN > 0. .AND. DrelG < DrelGP) DrelG = DrelGP
1236
1237   IF (DrelB <= 0.) then
1238   DrelB   = 0.0
1239   BETB    = 0.0
1240   ELSEIf (DrelB <= dengimp(IMPB)) then
1241   BETB  = DrelB/dengimp(IMPB)*porimp(IMPB)
1242   ELSE
1243   DrelB = dengimp(IMPB)
1244   BETB = porimp(IMPB)
1245   ENDIF
1246
1247   IF (DrelG <= 0.) then
1248   DrelG   = 0.0
1249   BETG    = 0.0
1250   ELSEIf (DrelG <= dengimp(IMPG)) then
1251   BETG  = DrelG/dengimp(IMPG)*porimp(IMPG)
1252   ELSE
1253   DrelG = dengimp(IMPG)
1254   BETG  = porimp(IMPG)
1255   ENDIF
1256
1257   if ( BETG < 1.E-5 ) BETG = 0.0
1258   if ( BETB < 1.E-5 ) BETB = 0.0
1259
1260 ENDIF
1261
1262    IF (TS_SCHEME == 1) THEN
1263
1264 !-------------------------------------------------------------------------------
1265 !  TB, TG  Solving Non-Linear Simultaneous Equation by Newton-Rapson
1266 !  TBL,TGL Solving Heat Equation by Tri Diagonal Matrix Algorithm
1267 !-------------------------------------------------------------------------------
1268
1269 !    TBP=350.
1270 !    TGP=350.
1271
1272      DO ITERATION=1,20
1273
1274        ES=6.11*EXP( (2.5*10.**6./461.51)*(TBP-273.15)/(273.15*TBP) )
1275        DESDT=(2.5*10.**6./461.51)*ES/(TBP**2.)
1276        QS0B=0.622*ES/(PS-0.378*ES)
1277        DQS0BDTB=DESDT*0.622*PS/((PS-0.378*ES)**2.)
1278
1279        ES=6.11*EXP( (2.5*10.**6./461.51)*(TGP-273.15)/(273.15*TGP) )
1280        DESDT=(2.5*10.**6./461.51)*ES/(TGP**2.)
1281        QS0G=0.622*ES/(PS-0.378*ES)
1282        DQS0GDTG=DESDT*0.22*PS/((PS-0.378*ES)**2.)
1283
1284        RG1=EPSG*( RX*VFGS          &
1285        +EPSB*VFGW*SIG*TBP**4./60.  &
1286        -SIG*TGP**4./60. )
1287
1288        RB1=EPSB*( RX*VFWS         &
1289        +EPSG*VFWG*SIG*TGP**4./60. &
1290        +EPSB*VFWW*SIG*TBP**4./60. &
1291        -SIG*TBP**4./60. )
1292
1293        RG2=EPSG*( (1.-EPSB)*(1.-SVF)*VFWS*RX                  &
1294        +(1.-EPSB)*(1.-SVF)*VFWG*EPSG*SIG*TGP**4./60.          &
1295        +EPSB*(1.-EPSB)*(1.-SVF)*(1.-2.*VFWS)*SIG*TBP**4./60. )
1296
1297        RB2=EPSB*( (1.-EPSG)*VFWG*VFGS*RX                          &
1298        +(1.-EPSG)*EPSB*VFGW*VFWG*SIG*(TBP**4.)/60.                &
1299        +(1.-EPSB)*VFWS*(1.-2.*VFWS)*RX                            &
1300        +(1.-EPSB)*VFWG*(1.-2.*VFWS)*EPSG*SIG*EPSG*TGP**4./60.     &
1301        +EPSB*(1.-EPSB)*(1.-2.*VFWS)*(1.-2.*VFWS)*SIG*TBP**4./60. )
1302
1303        RG=RG1+RG2
1304        RB=RB1+RB2
1305
1306        DRBDTB1=EPSB*(4.*EPSB*SIG*TB**3.*VFWW-4.*SIG*TB**3.)/60.
1307        DRBDTG1=EPSB*(4.*EPSG*SIG*TG**3.*VFWG)/60.
1308        DRBDTB2=EPSB*(4.*(1.-EPSG)*EPSB*SIG*TB**3.*VFGW*VFWG &
1309                +4.*EPSB*(1.-EPSB)*SIG*TB**3.*VFWW*VFWW)/60.
1310        DRBDTG2=EPSB*(4.*(1.-EPSB)*EPSG*SIG*TG**3.*VFWG*VFWW)/60.
1311
1312        DRGDTB1=EPSG*(4.*EPSB*SIG*TB**3.*VFGW)/60.
1313        DRGDTG1=EPSG*(-4.*SIG*TG**3.)/60.
1314        DRGDTB2=EPSG*(4.*EPSB*(1.-EPSB)*SIG*TB**3.*VFWW*VFGW)/60.
1315        DRGDTG2=EPSG*(4.*(1.-EPSB)*EPSG*SIG*TG**3.*VFWG*VFGW)/60.
1316
1317        DRBDTB=DRBDTB1+DRBDTB2
1318        DRBDTG=DRBDTG1+DRBDTG2
1319        DRGDTB=DRGDTB1+DRGDTB2
1320        DRGDTG=DRGDTG1+DRGDTG2
1321
1322        HB=RHO*CP*CHB*UC*(TBP-TCP)*100.
1323        HG=RHO*CP*CHG*UC*(TGP-TCP)*100.
1324
1325        DTCDTB=W*ALPHAB/(RW*ALPHAC+RW*ALPHAG+W*ALPHAB)
1326        DTCDTG=RW*ALPHAG/(RW*ALPHAC+RW*ALPHAG+W*ALPHAB)
1327
1328        DHBDTB=RHO*CP*CHB*UC*(1.-DTCDTB)*100.
1329        DHBDTG=RHO*CP*CHB*UC*(0.-DTCDTG)*100.
1330        DHGDTG=RHO*CP*CHG*UC*(1.-DTCDTG)*100.
1331        DHGDTB=RHO*CP*CHG*UC*(0.-DTCDTB)*100.
1332
1333        ELEB=RHO*EL*CHB*UC*BETB*(QS0B-QCP)*100.
1334        ELEG=RHO*EL*CHG*UC*BETG*(QS0G-QCP)*100.
1335
1336        DQCDTB=W*ALPHAB*BETB*DQS0BDTB/(RW*ALPHAC+RW*ALPHAG*BETG+W*ALPHAB*BETB)
1337        DQCDTG=RW*ALPHAG*BETG*DQS0GDTG/(RW*ALPHAC+RW*ALPHAG*BETG+W*ALPHAB*BETB)
1338
1339        DELEBDTB=RHO*EL*CHB*UC*BETB*(DQS0BDTB-DQCDTB)*100.
1340        DELEBDTG=RHO*EL*CHB*UC*BETB*(0.-DQCDTG)*100.
1341        DELEGDTG=RHO*EL*CHG*UC*BETG*(DQS0GDTG-DQCDTG)*100.
1342        DELEGDTB=RHO*EL*CHG*UC*BETG*(0.-DQCDTB)*100.
1343
1344        G0B=AKSB*(TBP-TBL(1))/(DZB(1)/2.)
1345        G0G=AKSG*(TGP-TGL(1))/(DZG(1)/2.)
1346
1347        DG0BDTB=2.*AKSB/DZB(1)
1348        DG0BDTG=0.
1349        DG0GDTG=2.*AKSG/DZG(1)
1350        DG0GDTB=0.
1351
1352        F = SB + RB - HB - ELEB - G0B
1353        FX = DRBDTB - DHBDTB - DELEBDTB - DG0BDTB
1354        FY = DRBDTG - DHBDTG - DELEBDTG - DG0BDTG
1355
1356        GF = SG + RG - HG - ELEG - G0G
1357        GX = DRGDTB - DHGDTB - DELEGDTB - DG0GDTB
1358        GY = DRGDTG - DHGDTG - DELEGDTG - DG0GDTG
1359
1360        DTB =  (GF*FY-F*GY)/(FX*GY-GX*FY)
1361        DTG = -(GF+GX*DTB)/GY
1362
1363        TB = TBP + DTB
1364        TG = TGP + DTG
1365
1366        TBP = TB
1367        TGP = TG
1368
1369        TC1=RW*ALPHAC+RW*ALPHAG+W*ALPHAB
1370        TC2=RW*ALPHAC*TA+RW*ALPHAG*TGP+W*ALPHAB*TBP
1371        TC=TC2/TC1
1372
1373        QC1=RW*ALPHAC+RW*ALPHAG*BETG+W*ALPHAB*BETB
1374        QC2=RW*ALPHAC*QA+RW*ALPHAG*BETG*QS0G+W*ALPHAB*BETB*QS0B
1375        QC=QC2/QC1
1376
1377        DTC=TCP - TC
1378        TCP=TC
1379        QCP=QC
1380
1381        IF( ABS(F) < 0.000001 .AND. ABS(DTB) < 0.000001 &
1382         .AND. ABS(GF) < 0.000001 .AND. ABS(DTG) < 0.000001 &
1383         .AND. ABS(DTC) < 0.000001) EXIT
1384
1385      END DO
1386
1387      CALL multi_layer(num_wall_layers,BOUNDB,G0B,CAPB,AKSB,TBL,DZB,DELT,TBLEND)
1388
1389      CALL multi_layer(num_road_layers,BOUNDG,G0G,CAPG,AKSG,TGL,DZG,DELT,TGLEND)
1390
1391    ELSE
1392
1393 !-------------------------------------------------------------------------------
1394 !  TB, TG  by Force-Restore Method
1395 !-------------------------------------------------------------------------------
1396
1397        ES=6.11*EXP((2.5*10.**6./461.51)*(TBP-273.15)/(273.15*TBP) )
1398        QS0B=0.622*ES/(PS-0.378*ES)
1399
1400        ES=6.11*EXP((2.5*10.**6./461.51)*(TGP-273.15)/(273.15*TGP) )
1401        QS0G=0.622*ES/(PS-0.378*ES)
1402
1403        RG1=EPSG*( RX*VFGS             &
1404        +EPSB*VFGW*SIG*TBP**4./60.     &
1405        -SIG*TGP**4./60. )
1406
1407        RB1=EPSB*( RX*VFWS &
1408        +EPSG*VFWG*SIG*TGP**4./60. &
1409        +EPSB*VFWW*SIG*TBP**4./60. &
1410        -SIG*TBP**4./60. )
1411
1412        RG2=EPSG*( (1.-EPSB)*(1.-SVF)*VFWS*RX                   &
1413        +(1.-EPSB)*(1.-SVF)*VFWG*EPSG*SIG*TGP**4./60.           &
1414        +EPSB*(1.-EPSB)*(1.-SVF)*(1.-2.*VFWS)*SIG*TBP**4./60. )
1415
1416        RB2=EPSB*( (1.-EPSG)*VFWG*VFGS*RX                          &
1417        +(1.-EPSG)*EPSB*VFGW*VFWG*SIG*(TBP**4.)/60.                &
1418        +(1.-EPSB)*VFWS*(1.-2.*VFWS)*RX                            &
1419        +(1.-EPSB)*VFWG*(1.-2.*VFWS)*EPSG*SIG*EPSG*TGP**4./60.     &
1420        +EPSB*(1.-EPSB)*(1.-2.*VFWS)*(1.-2.*VFWS)*SIG*TBP**4./60. )
1421
1422        RG=RG1+RG2
1423        RB=RB1+RB2
1424
1425        HB=RHO*CP*CHB*UC*(TBP-TCP)*100.
1426        ELEB=RHO*EL*CHB*UC*BETB*(QS0B-QCP)*100.
1427        G0B=SB+RB-HB-ELEB
1428
1429        HG=RHO*CP*CHG*UC*(TGP-TCP)*100.
1430        ELEG=RHO*EL*CHG*UC*BETG*(QS0G-QCP)*100.
1431        G0G=SG+RG-HG-ELEG
1432
1433        CALL force_restore(CAPB,AKSB,DELT,SB,RB,HB,ELEB,TBLEND,TBP,TB)
1434        CALL force_restore(CAPG,AKSG,DELT,SG,RG,HG,ELEG,TGLEND,TGP,TG)
1435
1436        TBP=TB
1437        TGP=TG
1438
1439        TC1=RW*ALPHAC+RW*ALPHAG+W*ALPHAB
1440        TC2=RW*ALPHAC*TA+RW*ALPHAG*TGP+W*ALPHAB*TBP
1441        TC=TC2/TC1
1442
1443        QC1=RW*ALPHAC+RW*ALPHAG*BETG+W*ALPHAB*BETB
1444        QC2=RW*ALPHAC*QA+RW*ALPHAG*BETG*QS0G+W*ALPHAB*BETB*QS0B
1445        QC=QC2/QC1
1446
1447        TCP=TC
1448        QCP=QC
1449
1450    END IF
1451
1452
1453    FLXTHB=HB/RHO/CP/100.
1454    FLXHUMB=ELEB/RHO/EL/100.
1455    FLXTHG=HG/RHO/CP/100.
1456    FLXHUMG=ELEG/RHO/EL/100.
1457
1458 !-------------------------------------------------------------------------------
1459 ! Total Fluxes from Urban Canopy
1460 !-------------------------------------------------------------------------------
1461 !===Yang, 2014/10/08, cal. ah. alh. green roof===
1462  if(groption==1) then
1463    if(ahoption==1) then
1464      FLXTH  = ((1.-FGR)*R*FLXTHR + FGR*R*FLXTHGR + W*FLXTHB + RW*FLXTHG)+ AH/RHOO/CPP
1465    else
1466      FLXTH  = ((1.-FGR)*R*FLXTHR + FGR*R*FLXTHGR + W*FLXTHB + RW*FLXTHG)
1467    endif
1468    if(alhoption==1) then
1469      FLXHUM  = ((1.-FGR)*R*FLXHUMR + FGR*R*FLXHUMGR + W*FLXHUMB + RW*FLXHUMG)+ ALH/RHOO/ELL
1470    else
1471      FLXHUM  = ((1.-FGR)*R*FLXHUMR + FGR*R*FLXHUMGR + W*FLXHUMB + RW*FLXHUMG)
1472    endif
1473    FLXUV  = ((1.-FGR)*R*CDR + FGR*R*CDGR + RW*CDC )*UA*UA
1474    FLXG =   ((1.-FGR)*R*G0R + FGR*R*G0GR+ W*G0B + RW*G0G)
1475    LNET =   (1.-FGR) * R * RR + FGR *R* RGR + W * RB +  RW * RG
1476  else
1477    if(ahoption==1) then
1478      FLXTH  = ( R*FLXTHR  + W*FLXTHB  + RW*FLXTHG ) + AH/RHOO/CPP
1479    else
1480      FLXTH  = ( R*FLXTHR  + W*FLXTHB  + RW*FLXTHG )
1481    endif
1482    if(alhoption==1) then
1483      FLXHUM = ( R*FLXHUMR + W*FLXHUMB + RW*FLXHUMG )+ ALH/RHOO/ELL
1484    else
1485      FLXHUM = ( R*FLXHUMR + W*FLXHUMB + RW*FLXHUMG )
1486    endif
1487    FLXUV  = ( R*CDR + RW*CDC )*UA*UA
1488    FLXG =   ( R*G0R + W*G0B + RW*G0G )
1489    LNET =     R*RR + W*RB + RW*RG
1490  endif
1491
1492 !----------------------------------------------------------------------------
1493 ! Convert Unit: FLUXES and u* T* q*  --> WRF
1494 !----------------------------------------------------------------------------
1495
1496    SH    = FLXTH  * RHOO * CPP    ! Sensible heat flux          [W/m/m]
1497    LH    = FLXHUM * RHOO * ELL    ! Latent heat flux            [W/m/m]
1498    LH_KINEMATIC = FLXHUM * RHOO   ! Latent heat, Kinematic      [kg/m/m/s]
1499    LW    = LLG - (LNET*697.7*60.) ! Upward longwave radiation   [W/m/m]
1500    SW    = SSG - (SNET*697.7*60.) ! Upward shortwave radiation  [W/m/m]
1501    ALB   = 0.
1502    IF( ABS(SSG) > 0.0001) ALB = SW/SSG ! Effective albedo [-]
1503    G = -FLXG*697.7*60.            ! [W/m/m]
1504    RN = (SNET+LNET)*697.7*60.     ! Net radiation [W/m/m]
1505
1506    UST = SQRT(FLXUV)              ! u* [m/s]
1507    TST = -FLXTH/UST               ! T* [K]
1508    QST = -FLXHUM/UST              ! q* [-]
1509
1510 !------------------------------------------------------
1511 !  diagnostic GRID AVERAGED  PSIM  PSIH  TS QS --> WRF
1512 !------------------------------------------------------
1513
1514    Z0 = Z0C
1515    Z0H = Z0HC
1516    Z = ZA - ZDC
1517    ZNT = Z0   ! add by Dan Li
1518
1519    XXX = 0.4*9.81*Z*TST/TA/UST/UST
1520
1521    IF ( XXX >= 1. ) XXX = 1.
1522    IF ( XXX <= -5. ) XXX = -5.
1523
1524    IF ( XXX > 0 ) THEN
1525      PSIM = -5. * XXX
1526      PSIH = -5. * XXX
1527    ELSE
1528      X = (1.-16.*XXX)**0.25
1529      PSIM = 2.*ALOG((1.+X)/2.) + ALOG((1.+X*X)/2.) - 2.*ATAN(X) + PI/2.
1530      PSIH = 2.*ALOG((1.+X*X)/2.)
1531    END IF
1532
1533    GZ1OZ0 = ALOG(Z/Z0)
1534    CD = 0.4**2./(ALOG(Z/Z0)-PSIM)**2.
1535 !
1536 !m   CH = 0.4**2./(ALOG(Z/Z0)-PSIM)/(ALOG(Z/Z0H)-PSIH)
1537 !m   CHS = 0.4*UST/(ALOG(Z/Z0H)-PSIH)
1538 !m   TS = TA + FLXTH/CH/UA    ! surface potential temp (flux temp)
1539 !m   QS = QA + FLXHUM/CH/UA   ! surface humidity
1540 !
1541    TS = TA + FLXTH/CHS    ! surface potential temp (flux temp)
1542    QS = QA + FLXHUM/CHS   ! surface humidity
1543
1544 !-------------------------------------------------------
1545 !  diagnostic  GRID AVERAGED  U10  V10  TH2  Q2 --> WRF
1546 !-------------------------------------------------------
1547
1548    XXX2 = (2./Z)*XXX
1549    IF ( XXX2 >= 1. ) XXX2 = 1.
1550    IF ( XXX2 <= -5. ) XXX2 = -5.
1551
1552    IF ( XXX2 > 0 ) THEN
1553       PSIM2 = -5. * XXX2
1554       PSIH2 = -5. * XXX2
1555    ELSE
1556       X = (1.-16.*XXX2)**0.25
1557       PSIM2 = 2.*ALOG((1.+X)/2.) + ALOG((1.+X*X)/2.) - 2.*ATAN(X) + 2.*ATAN(1.)
1558       PSIH2 = 2.*ALOG((1.+X*X)/2.)
1559    END IF
1560 !
1561 !m   CHS2 = 0.4*UST/(ALOG(2./Z0H)-PSIH2)
1562 !
1563
1564    XXX10 = (10./Z)*XXX
1565    IF ( XXX10 >= 1. ) XXX10 = 1.
1566    IF ( XXX10 <= -5. ) XXX10 = -5.
1567
1568    IF ( XXX10 > 0 ) THEN
1569       PSIM10 = -5. * XXX10
1570       PSIH10 = -5. * XXX10
1571    ELSE
1572       X = (1.-16.*XXX10)**0.25
1573       PSIM10 = 2.*ALOG((1.+X)/2.) + ALOG((1.+X*X)/2.) - 2.*ATAN(X) + 2.*ATAN(1.)
1574       PSIH10 = 2.*ALOG((1.+X*X)/2.)
1575    END IF
1576
1577    PSIX = ALOG(Z/Z0) - PSIM
1578    PSIT = ALOG(Z/Z0H) - PSIH
1579
1580    PSIX2 = ALOG(2./Z0) - PSIM2
1581    PSIT2 = ALOG(2./Z0H) - PSIH2
1582
1583    PSIX10 = ALOG(10./Z0) - PSIM10
1584    PSIT10 = ALOG(10./Z0H) - PSIH10
1585
1586    U10 = U1 * (PSIX10/PSIX)       ! u at 10 m [m/s]
1587    V10 = V1 * (PSIX10/PSIX)       ! v at 10 m [m/s]
1588
1589 !   TH2 = TS + (TA-TS)*(PSIT2/PSIT)   ! potential temp at 2 m [K]
1590 !  TH2 = TS + (TA-TS)*(PSIT2/PSIT)   ! Fei: this seems to be temp (not potential) at 2 m [K]
1591 !Fei: consistant with M-O theory
1592    TH2 = TS + (TA-TS) *(CHS/CHS2)
1593
1594    Q2 = QS + (QA-QS)*(PSIT2/PSIT)    ! humidity at 2 m       [-]
1595
1596 !   TS = (LW/SIG_SI/0.88)**0.25       ! Radiative temperature [K]
1597
1598    END SUBROUTINE urban
1599 !===============================================================================
1600 !
1601 !  mos
1602 !
1603 !===============================================================================
1604    SUBROUTINE mos(XXX,ALPHA,CD,B1,RIB,Z,Z0,UA,TA,TSF,RHO)
1605
1606 !  XXX:   z/L (requires iteration by Newton-Rapson method)
1607 !  B1:    Stanton number
1608 !  PSIM:  = PSIX of LSM
1609 !  PSIH:  = PSIT of LSM
1610
1611    IMPLICIT NONE
1612
1613    REAL, PARAMETER     :: CP=0.24
1614    REAL, INTENT(IN)    :: B1, Z, Z0, UA, TA, TSF, RHO
1615    REAL, INTENT(OUT)   :: ALPHA, CD
1616    REAL, INTENT(INOUT) :: XXX, RIB
1617    REAL                :: XXX0, X, X0, FAIH, DPSIM, DPSIH
1618    REAL                :: F, DF, XXXP, US, TS, AL, XKB, DD, PSIM, PSIH
1619    INTEGER             :: NEWT
1620    INTEGER, PARAMETER  :: NEWT_END=10
1621
1622    IF(RIB <= -15.) RIB=-15.
1623
1624    IF(RIB < 0.) THEN
1625
1626       DO NEWT=1,NEWT_END
1627
1628         IF(XXX >= 0.) XXX=-1.E-3
1629
1630         XXX0=XXX*Z0/(Z+Z0)
1631
1632         X=(1.-16.*XXX)**0.25
1633         X0=(1.-16.*XXX0)**0.25
1634
1635         PSIM=ALOG((Z+Z0)/Z0) &
1636             -ALOG((X+1.)**2.*(X**2.+1.)) &
1637             +2.*ATAN(X) &
1638             +ALOG((X+1.)**2.*(X0**2.+1.)) &
1639             -2.*ATAN(X0)
1640         FAIH=1./SQRT(1.-16.*XXX)
1641         PSIH=ALOG((Z+Z0)/Z0)+0.4*B1 &
1642             -2.*ALOG(SQRT(1.-16.*XXX)+1.) &
1643             +2.*ALOG(SQRT(1.-16.*XXX0)+1.)
1644
1645         DPSIM=(1.-16.*XXX)**(-0.25)/XXX &
1646              -(1.-16.*XXX0)**(-0.25)/XXX
1647         DPSIH=1./SQRT(1.-16.*XXX)/XXX &
1648              -1./SQRT(1.-16.*XXX0)/XXX
1649
1650         F=RIB*PSIM**2./PSIH-XXX
1651
1652         DF=RIB*(2.*DPSIM*PSIM*PSIH-DPSIH*PSIM**2.) &
1653           /PSIH**2.-1.
1654
1655         XXXP=XXX
1656         XXX=XXXP-F/DF
1657         IF(XXX <= -10.) XXX=-10.
1658
1659       END DO
1660
1661    ELSE IF(RIB >= 0.142857) THEN
1662
1663       XXX=0.714
1664       PSIM=ALOG((Z+Z0)/Z0)+7.*XXX
1665       PSIH=PSIM+0.4*B1
1666
1667    ELSE
1668
1669       AL=ALOG((Z+Z0)/Z0)
1670       XKB=0.4*B1
1671       DD=-4.*RIB*7.*XKB*AL+(AL+XKB)**2.
1672       IF(DD <= 0.) DD=0.
1673       XXX=(AL+XKB-2.*RIB*7.*AL-SQRT(DD))/(2.*(RIB*7.**2-7.))
1674       PSIM=ALOG((Z+Z0)/Z0)+7.*MIN(XXX,0.714)
1675       PSIH=PSIM+0.4*B1
1676
1677    END IF
1678
1679    US=0.4*UA/PSIM             ! u*
1680    IF(US <= 0.01) US=0.01
1681    TS=0.4*(TA-TSF)/PSIH       ! T*
1682
1683    CD=US*US/UA**2.            ! CD
1684    ALPHA=RHO*CP*0.4*US/PSIH   ! RHO*CP*CH*U
1685
1686    RETURN
1687    END SUBROUTINE mos
1688 !===============================================================================
1689 !
1690 !  louis79
1691 !
1692 !===============================================================================
1693    SUBROUTINE louis79(ALPHA,CD,RIB,Z,Z0,UA,RHO)
1694
1695    IMPLICIT NONE
1696
1697    REAL, PARAMETER     :: CP=0.24
1698    REAL, INTENT(IN)    :: Z, Z0, UA, RHO
1699    REAL, INTENT(OUT)   :: ALPHA, CD
1700    REAL, INTENT(INOUT) :: RIB
1701    REAL                :: A2, XX, CH, CMB, CHB
1702
1703    A2=(0.4/ALOG(Z/Z0))**2.
1704
1705    IF(RIB <= -15.) RIB=-15.
1706
1707    IF(RIB >= 0.0) THEN
1708       IF(RIB >= 0.142857) THEN
1709          XX=0.714
1710       ELSE
1711          XX=RIB*LOG(Z/Z0)/(1.-7.*RIB)
1712       END IF
1713       CH=0.16/0.74/(LOG(Z/Z0)+7.*MIN(XX,0.714))**2.
1714       CD=0.16/(LOG(Z/Z0)+7.*MIN(XX,0.714))**2.
1715    ELSE
1716       CMB=7.4*A2*9.4*SQRT(Z/Z0)
1717       CHB=5.3*A2*9.4*SQRT(Z/Z0)
1718       CH=A2/0.74*(1.-9.4*RIB/(1.+CHB*SQRT(-RIB)))
1719       CD=A2*(1.-9.4*RIB/(1.+CHB*SQRT(-RIB)))
1720    END IF
1721
1722    ALPHA=RHO*CP*CH*UA
1723
1724    RETURN
1725    END SUBROUTINE louis79
1726 !===============================================================================
1727 !
1728 !  louis82
1729 !
1730 !===============================================================================
1731    SUBROUTINE louis82(ALPHA,CD,RIB,Z,Z0,UA,RHO)
1732
1733    IMPLICIT NONE
1734
1735    REAL, PARAMETER     :: CP=0.24
1736    REAL, INTENT(IN)    :: Z, Z0, UA, RHO
1737    REAL, INTENT(OUT)   :: ALPHA, CD
1738    REAL, INTENT(INOUT) :: RIB
1739    REAL                :: A2, FM, FH, CH, CHH
1740
1741    A2=(0.4/ALOG(Z/Z0))**2.
1742
1743    IF(RIB <= -15.) RIB=-15.
1744
1745    IF(RIB >= 0.0) THEN
1746       FM=1./((1.+(2.*5.*RIB)/SQRT(1.+5.*RIB)))
1747       FH=1./(1.+(3.*5.*RIB)*SQRT(1.+5.*RIB))
1748       CH=A2*FH
1749       CD=A2*FM
1750    ELSE
1751       CHH=5.*3.*5.*A2*SQRT(Z/Z0)
1752       FM=1.-(2.*5.*RIB)/(1.+3.*5.*5.*A2*SQRT(Z/Z0+1.)*(-RIB))
1753       FH=1.-(3.*5.*RIB)/(1.+CHH*SQRT(-RIB))
1754       CH=A2*FH
1755       CD=A2*FM
1756    END IF
1757
1758    ALPHA=RHO*CP*CH*UA
1759
1760    RETURN
1761    END SUBROUTINE louis82
1762 !===============================================================================
1763 !
1764 !  multi_layer
1765 !
1766 !===============================================================================
1767    SUBROUTINE multi_layer(KM,BOUND,G0,CAP,AKS,TSL,DZ,DELT,TSLEND)
1768
1769    IMPLICIT NONE
1770
1771    REAL, INTENT(IN)                   :: G0
1772
1773    REAL, INTENT(IN)                   :: CAP
1774
1775    REAL, INTENT(IN)                   :: AKS
1776
1777    REAL, INTENT(IN)                   :: DELT      ! Time step [ s ]
1778
1779    REAL, INTENT(IN)                   :: TSLEND
1780
1781    INTEGER, INTENT(IN)                :: KM
1782
1783    INTEGER, INTENT(IN)                :: BOUND
1784
1785    REAL, DIMENSION(KM), INTENT(IN)    :: DZ
1786
1787    REAL, DIMENSION(KM), INTENT(INOUT) :: TSL
1788
1789    REAL, DIMENSION(KM)                :: A, B, C, D, X, P, Q
1790
1791    REAL                               :: DZEND
1792
1793    INTEGER                            :: K
1794
1795    DZEND=DZ(KM)
1796
1797    A(1) = 0.0
1798
1799    B(1) = CAP*DZ(1)/DELT &
1800           +2.*AKS/(DZ(1)+DZ(2))
1801    C(1) = -2.*AKS/(DZ(1)+DZ(2))
1802    D(1) = CAP*DZ(1)/DELT*TSL(1) + G0
1803
1804    DO K=2,KM-1
1805       A(K) = -2.*AKS/(DZ(K-1)+DZ(K))
1806       B(K) = CAP*DZ(K)/DELT + 2.*AKS/(DZ(K-1)+DZ(K)) + 2.*AKS/(DZ(K)+DZ(K+1))
1807       C(K) = -2.*AKS/(DZ(K)+DZ(K+1))
1808       D(K) = CAP*DZ(K)/DELT*TSL(K)
1809    END DO
1810
1811    IF(BOUND == 1) THEN                  ! Flux=0
1812       A(KM) = -2.*AKS/(DZ(KM-1)+DZ(KM))
1813       B(KM) = CAP*DZ(KM)/DELT + 2.*AKS/(DZ(KM-1)+DZ(KM))
1814       C(KM) = 0.0
1815       D(KM) = CAP*DZ(KM)/DELT*TSL(KM)
1816    ELSE                                 ! T=constant
1817       A(KM) = -2.*AKS/(DZ(KM-1)+DZ(KM))
1818       B(KM) = CAP*DZ(KM)/DELT + 2.*AKS/(DZ(KM-1)+DZ(KM)) + 2.*AKS/(DZ(KM)+DZEND)
1819       C(KM) = 0.0
1820       D(KM) = CAP*DZ(KM)/DELT*TSL(KM) + 2.*AKS*TSLEND/(DZ(KM)+DZEND)
1821    END IF
1822
1823    P(1) = -C(1)/B(1)
1824    Q(1) =  D(1)/B(1)
1825
1826    DO K=2,KM
1827       P(K) = -C(K)/(A(K)*P(K-1)+B(K))
1828       Q(K) = (-A(K)*Q(K-1)+D(K))/(A(K)*P(K-1)+B(K))
1829    END DO
1830
1831    X(KM) = Q(KM)
1832
1833    DO K=KM-1,1,-1
1834       X(K) = P(K)*X(K+1)+Q(K)
1835    END DO
1836
1837    DO K=1,KM
1838       TSL(K) = X(K)
1839    END DO
1840
1841    RETURN
1842    END SUBROUTINE multi_layer
1843 !===============================================================================
1844 !
1845 !  subroutine read_param
1846 !
1847 !===============================================================================
1848    SUBROUTINE read_param(UTYPE,                                        & ! in
1849                          ZR,SIGMA_ZED,Z0C,Z0HC,ZDC,SVF,R,RW,HGT,AH,    & ! out
1850                          CAPR,CAPB,CAPG,AKSR,AKSB,AKSG,ALBR,ALBB,ALBG, & ! out
1851                          EPSR,EPSB,EPSG,Z0R,Z0B,Z0G,Z0HB,Z0HG,         & ! out
1852                          BETR,BETB,BETG,TRLEND,TBLEND,TGLEND,          & ! out
1853 !for BEP
1854                          NUMDIR, STREET_DIRECTION, STREET_WIDTH,       & ! out
1855                          BUILDING_WIDTH, NUMHGT, HEIGHT_BIN,           & ! out
1856                          HPERCENT_BIN,                                 & ! out
1857 !end BEP
1858                          BOUNDR,BOUNDB,BOUNDG,CH_SCHEME,TS_SCHEME,     & ! out
1859                          AKANDA_URBAN,ALH)         ! out
1860
1861    INTEGER, INTENT(IN)  :: UTYPE
1862
1863    REAL, INTENT(OUT)    :: ZR,Z0C,Z0HC,ZDC,SVF,R,RW,HGT,AH,ALH,          &
1864                            CAPR,CAPB,CAPG,AKSR,AKSB,AKSG,ALBR,ALBB,ALBG, &
1865                            SIGMA_ZED,                                    &
1866                            EPSR,EPSB,EPSG,Z0R,Z0B,Z0G,Z0HB,Z0HG,         &
1867                            BETR,BETB,BETG,TRLEND,TBLEND,TGLEND
1868    REAL, INTENT(OUT)    :: AKANDA_URBAN
1869 !for BEP
1870    INTEGER,                     INTENT(OUT) :: NUMDIR
1871    REAL,    DIMENSION(MAXDIRS), INTENT(OUT) :: STREET_DIRECTION
1872    REAL,    DIMENSION(MAXDIRS), INTENT(OUT) :: STREET_WIDTH
1873    REAL,    DIMENSION(MAXDIRS), INTENT(OUT) :: BUILDING_WIDTH
1874    INTEGER,                     INTENT(OUT) :: NUMHGT
1875    REAL,    DIMENSION(MAXHGTS), INTENT(OUT) :: HEIGHT_BIN
1876    REAL,    DIMENSION(MAXHGTS), INTENT(OUT) :: HPERCENT_BIN
1877
1878 !end BEP
1879
1880    INTEGER, INTENT(OUT) :: BOUNDR,BOUNDB,BOUNDG,CH_SCHEME,TS_SCHEME
1881
1882    ZR =     ZR_TBL(UTYPE)
1883    SIGMA_ZED = SIGMA_ZED_TBL(UTYPE)
1884    Z0C=     Z0C_TBL(UTYPE)
1885    Z0HC=    Z0HC_TBL(UTYPE)
1886    ZDC=     ZDC_TBL(UTYPE)
1887    SVF=     SVF_TBL(UTYPE)
1888    R=       R_TBL(UTYPE)
1889    RW=      RW_TBL(UTYPE)
1890    HGT=     HGT_TBL(UTYPE)
1891    AH=      AH_TBL(UTYPE)
1892    ALH=     ALH_TBL(UTYPE)
1893    BETR=    BETR_TBL(UTYPE)
1894    BETB=    BETB_TBL(UTYPE)
1895    BETG=    BETG_TBL(UTYPE)
1896
1897 !m   FRC_URB= FRC_URB_TBL(UTYPE)
1898
1899    CAPR=    CAPR_TBL(UTYPE)
1900    CAPB=    CAPB_TBL(UTYPE)
1901    CAPG=    CAPG_TBL(UTYPE)
1902    AKSR=    AKSR_TBL(UTYPE)
1903    AKSB=    AKSB_TBL(UTYPE)
1904    AKSG=    AKSG_TBL(UTYPE)
1905    ALBR=    ALBR_TBL(UTYPE)
1906    ALBB=    ALBB_TBL(UTYPE)
1907    ALBG=    ALBG_TBL(UTYPE)
1908    EPSR=    EPSR_TBL(UTYPE)
1909    EPSB=    EPSB_TBL(UTYPE)
1910    EPSG=    EPSG_TBL(UTYPE)
1911    Z0R=     Z0R_TBL(UTYPE)
1912    Z0B=     Z0B_TBL(UTYPE)
1913    Z0G=     Z0G_TBL(UTYPE)
1914    Z0HB=    Z0HB_TBL(UTYPE)
1915    Z0HG=    Z0HG_TBL(UTYPE)
1916    TRLEND=  TRLEND_TBL(UTYPE)
1917    TBLEND=  TBLEND_TBL(UTYPE)
1918    TGLEND=  TGLEND_TBL(UTYPE)
1919    BOUNDR=  BOUNDR_DATA
1920    BOUNDB=  BOUNDB_DATA
1921    BOUNDG=  BOUNDG_DATA
1922    CH_SCHEME = CH_SCHEME_DATA
1923    TS_SCHEME = TS_SCHEME_DATA
1924    AKANDA_URBAN = AKANDA_URBAN_TBL(UTYPE)
1925
1926 !for BEP
1927
1928    STREET_DIRECTION = -1.E36
1929    STREET_WIDTH     = -1.E36
1930    BUILDING_WIDTH   = -1.E36
1931    HEIGHT_BIN       = -1.E36
1932    HPERCENT_BIN     = -1.E36
1933
1934    NUMDIR                     = NUMDIR_TBL ( UTYPE )
1935    STREET_DIRECTION(1:NUMDIR) = STREET_DIRECTION_TBL( 1:NUMDIR, UTYPE )
1936    STREET_WIDTH    (1:NUMDIR) = STREET_WIDTH_TBL    ( 1:NUMDIR, UTYPE )
1937    BUILDING_WIDTH  (1:NUMDIR) = BUILDING_WIDTH_TBL  ( 1:NUMDIR, UTYPE )
1938    NUMHGT                     = NUMHGT_TBL ( UTYPE )
1939    HEIGHT_BIN      (1:NUMHGT) = HEIGHT_BIN_TBL   ( 1:NUMHGT , UTYPE )
1940    HPERCENT_BIN    (1:NUMHGT) = HPERCENT_BIN_TBL ( 1:NUMHGT , UTYPE )
1941
1942 !end BEP
1943    END SUBROUTINE read_param
1944 !===============================================================================
1945 !
1946 ! subroutine urban_param_init: Read parameters from URBPARM.TBL
1947 !
1948 !===============================================================================
1949    SUBROUTINE urban_param_init(DZR,DZB,DZG,num_soil_layers, &
1950                                sf_urban_physics,use_wudapt_lcz)
1951 !                              num_roof_layers,num_wall_layers,num_road_layers)
1952
1953    IMPLICIT NONE
1954
1955    INTEGER, INTENT(IN) :: num_soil_layers
1956
1957 !   REAL, DIMENSION(1:num_roof_layers), INTENT(INOUT) :: DZR
1958 !   REAL, DIMENSION(1:num_wall_layers), INTENT(INOUT) :: DZB
1959 !   REAL, DIMENSION(1:num_road_layers), INTENT(INOUT) :: DZG
1960    REAL, DIMENSION(1:num_soil_layers), INTENT(INOUT) :: DZR
1961    REAL, DIMENSION(1:num_soil_layers), INTENT(INOUT) :: DZB
1962    REAL, DIMENSION(1:num_soil_layers), INTENT(INOUT) :: DZG
1963    INTEGER,                            INTENT(IN)    :: SF_URBAN_PHYSICS
1964    INTEGER,                            INTENT(IN)    :: USE_WUDAPT_LCZ !AndreaLCZ
1965
1966    INTEGER :: LC, K
1967    INTEGER :: IOSTATUS, ALLOCATE_STATUS
1968    INTEGER :: num_roof_layers
1969    INTEGER :: num_wall_layers
1970    INTEGER :: num_road_layers
1971    INTEGER :: dummy
1972    REAL    :: DHGT, HGT, VFWS, VFGS
1973
1974    REAL, allocatable, dimension(:) :: ROOF_WIDTH
1975    REAL, allocatable, dimension(:) :: ROAD_WIDTH
1976
1977    character(len=512) :: string
1978    character(len=128) :: name
1979    integer :: indx
1980
1981    real, parameter :: VonK = 0.4
1982    real :: lambda_p
1983    real :: lambda_f
1984    real :: Cd
1985    real :: alpha_macd
1986    real :: beta_macd
1987    real :: lambda_fr
1988
1989
1990 !for BEP
1991    real :: dummy_hgt
1992    real :: dummy_pct
1993    real :: pctsum
1994 !end BEP
1995    num_roof_layers = num_soil_layers
1996    num_wall_layers = num_soil_layers
1997    num_road_layers = num_soil_layers
1998
1999
2000    ICATE=0
2001
2002
2003   if(USE_WUDAPT_LCZ.eq.0)then   !AndreaLCZ
2004    OPEN (UNIT=11,                &
2005          FILE='URBPARM.TBL', &
2006          ACCESS='SEQUENTIAL',    &
2007          STATUS='OLD',           &
2008          ACTION='READ',          &
2009          POSITION='REWIND',      &
2010          IOSTAT=IOSTATUS)
2011
2012          IF (IOSTATUS > 0) THEN
2013          FATAL_ERROR('Error opening URBPARM.TBL. Make sure URBPARM.TBL (found in run/) is linked to the running directory.')
2014          ENDIF
2015
2016  else
2017    OPEN (UNIT=11,                &
2018          FILE='URBPARM_LCZ.TBL', &
2019          ACCESS='SEQUENTIAL',    &
2020          STATUS='OLD',           &
2021          ACTION='READ',          &
2022          POSITION='REWIND',      &
2023          IOSTAT=IOSTATUS)
2024
2025          IF (IOSTATUS > 0) THEN
2026          FATAL_ERROR('Error opening URBPARM_LCZ.TBL. Make sure URBPARM_LCZ.TBL (found in run/) is linked to the running directory.')
2027          ENDIF
2028  endif
2029
2030
2031    READLOOP : do
2032       read(11,'(A512)', iostat=iostatus) string
2033       if (iostatus /= 0) exit READLOOP
2034       if (string(1:1) == "#") cycle READLOOP
2035       if (trim(string) == "") cycle READLOOP
2036       indx = index(string,":")
2037       if (indx <= 0) cycle READLOOP
2038       name = trim(adjustl(string(1:indx-1)))
2039
2040       ! Here are the variables we expect to be defined in the URBPARM.TBL:
2041       if (name == "Number of urban categories") then
2042          read(string(indx+1:),*) icate
2043          IF (.not. ALLOCATED(ZR_TBL)) then
2044             ALLOCATE( ZR_TBL(ICATE), stat=allocate_status )
2045             if(allocate_status /= 0) FATAL_ERROR('Error allocating ZR_TBL in urban_param_init')
2046             ALLOCATE( SIGMA_ZED_TBL(ICATE), stat=allocate_status )
2047             if(allocate_status /= 0) FATAL_ERROR('Error allocating SIGMA_ZED_TBL in urban_param_init')
2048             ALLOCATE( Z0C_TBL(ICATE), stat=allocate_status )
2049             if(allocate_status /= 0) FATAL_ERROR('Error allocating Z0C_TBL in urban_param_init')
2050             ALLOCATE( Z0HC_TBL(ICATE), stat=allocate_status )
2051             if(allocate_status /= 0) FATAL_ERROR('Error allocating Z0HC_TBL in urban_param_init')
2052             ALLOCATE( ZDC_TBL(ICATE), stat=allocate_status )
2053             if(allocate_status /= 0) FATAL_ERROR('Error allocating ZDC_TBL in urban_param_init')
2054             ALLOCATE( SVF_TBL(ICATE), stat=allocate_status )
2055             if(allocate_status /= 0) FATAL_ERROR('Error allocating SVF_TBL in urban_param_init')
2056             ALLOCATE( R_TBL(ICATE), stat=allocate_status )
2057             if(allocate_status /= 0) FATAL_ERROR('Error allocating R_TBL in urban_param_init')
2058             ALLOCATE( RW_TBL(ICATE), stat=allocate_status )
2059             if(allocate_status /= 0) FATAL_ERROR('Error allocating RW_TBL in urban_param_init')
2060             ALLOCATE( HGT_TBL(ICATE), stat=allocate_status )
2061             if(allocate_status /= 0) FATAL_ERROR('Error allocating HGT_TBL in urban_param_init')
2062             ALLOCATE( AH_TBL(ICATE), stat=allocate_status )
2063             if(allocate_status /= 0) FATAL_ERROR('Error allocating AH_TBL in urban_param_init')
2064             ALLOCATE( ALH_TBL(ICATE), stat=allocate_status )
2065             if(allocate_status /= 0) FATAL_ERROR('Error allocating ALH_TBL in urban_param_init')
2066             ALLOCATE( BETR_TBL(ICATE), stat=allocate_status )
2067             if(allocate_status /= 0) FATAL_ERROR('Error allocating BETR_TBL in urban_param_init')
2068             ALLOCATE( BETB_TBL(ICATE), stat=allocate_status )
2069             if(allocate_status /= 0) FATAL_ERROR('Error allocating BETB_TBL in urban_param_init')
2070             ALLOCATE( BETG_TBL(ICATE), stat=allocate_status )
2071             if(allocate_status /= 0) FATAL_ERROR('Error allocating BETG_TBL in urban_param_init')
2072             ALLOCATE( CAPR_TBL(ICATE), stat=allocate_status )
2073             if(allocate_status /= 0) FATAL_ERROR('Error allocating CAPR_TBL in urban_param_init')
2074             ALLOCATE( CAPB_TBL(ICATE), stat=allocate_status )
2075             if(allocate_status /= 0) FATAL_ERROR('Error allocating CAPB_TBL in urban_param_init')
2076             ALLOCATE( CAPG_TBL(ICATE), stat=allocate_status )
2077             if(allocate_status /= 0) FATAL_ERROR('Error allocating CAPG_TBL in urban_param_init')
2078             ALLOCATE( AKSR_TBL(ICATE), stat=allocate_status )
2079             if(allocate_status /= 0) FATAL_ERROR('Error allocating AKSR_TBL in urban_param_init')
2080             ALLOCATE( AKSB_TBL(ICATE), stat=allocate_status )
2081             if(allocate_status /= 0) FATAL_ERROR('Error allocating AKSB_TBL in urban_param_init')
2082             ALLOCATE( AKSG_TBL(ICATE), stat=allocate_status )
2083             if(allocate_status /= 0) FATAL_ERROR('Error allocating AKSG_TBL in urban_param_init')
2084             ALLOCATE( ALBR_TBL(ICATE), stat=allocate_status )
2085             if(allocate_status /= 0) FATAL_ERROR('Error allocating ALBR_TBL in urban_param_init')
2086             ALLOCATE( ALBB_TBL(ICATE), stat=allocate_status )
2087             if(allocate_status /= 0) FATAL_ERROR('Error allocating ALBB_TBL in urban_param_init')
2088             ALLOCATE( ALBG_TBL(ICATE), stat=allocate_status )
2089             if(allocate_status /= 0) FATAL_ERROR('Error allocating ALBG_TBL in urban_param_init')
2090             ALLOCATE( EPSR_TBL(ICATE), stat=allocate_status )
2091             if(allocate_status /= 0) FATAL_ERROR('Error allocating EPSR_TBL in urban_param_init')
2092             ALLOCATE( EPSB_TBL(ICATE), stat=allocate_status )
2093             if(allocate_status /= 0) FATAL_ERROR('Error allocating EPSB_TBL in urban_param_init')
2094             ALLOCATE( EPSG_TBL(ICATE), stat=allocate_status )
2095             if(allocate_status /= 0) FATAL_ERROR('Error allocating EPSG_TBL in urban_param_init')
2096             ALLOCATE( Z0R_TBL(ICATE), stat=allocate_status )
2097             if(allocate_status /= 0) FATAL_ERROR('Error allocating Z0R_TBL in urban_param_init')
2098             ALLOCATE( Z0B_TBL(ICATE), stat=allocate_status )
2099             if(allocate_status /= 0) FATAL_ERROR('Error allocating Z0B_TBL in urban_param_init')
2100             ALLOCATE( Z0G_TBL(ICATE), stat=allocate_status )
2101             if(allocate_status /= 0) FATAL_ERROR('Error allocating Z0G_TBL in urban_param_init')
2102             ALLOCATE( AKANDA_URBAN_TBL(ICATE), stat=allocate_status )
2103             if(allocate_status /= 0) FATAL_ERROR('Error allocating AKANDA_URBAN_TBL in urban_param_init')
2104             ALLOCATE( Z0HB_TBL(ICATE), stat=allocate_status )
2105             if(allocate_status /= 0) FATAL_ERROR('Error allocating Z0HB_TBL in urban_param_init')
2106             ALLOCATE( Z0HG_TBL(ICATE), stat=allocate_status )
2107             if(allocate_status /= 0) FATAL_ERROR('Error allocating Z0HG_TBL in urban_param_init')
2108             ALLOCATE( TRLEND_TBL(ICATE), stat=allocate_status )
2109             if(allocate_status /= 0) FATAL_ERROR('Error allocating TRLEND_TBL in urban_param_init')
2110             ALLOCATE( TBLEND_TBL(ICATE), stat=allocate_status )
2111             if(allocate_status /= 0) FATAL_ERROR('Error allocating TBLEND_TBL in urban_param_init')
2112             ALLOCATE( TGLEND_TBL(ICATE), stat=allocate_status )
2113             if(allocate_status /= 0) FATAL_ERROR('Error allocating TGLEND_TBL in urban_param_init')
2114             ALLOCATE( FRC_URB_TBL(ICATE), stat=allocate_status )
2115             if(allocate_status /= 0) FATAL_ERROR('Error allocating FRC_URB_TBL in urban_param_init')
2116             ! ALLOCATE( ROOF_WIDTH(ICATE), stat=allocate_status )
2117             ! if(allocate_status /= 0) FATAL_ERROR('Error allocating ROOF_WIDTH in urban_param_init')
2118             ! ALLOCATE( ROAD_WIDTH(ICATE), stat=allocate_status )
2119             ! if(allocate_status /= 0) FATAL_ERROR('Error allocating ROAD_WIDTH in urban_param_init')
2120             !for BEP
2121             ALLOCATE( NUMDIR_TBL(ICATE), stat=allocate_status )
2122             if(allocate_status /= 0) FATAL_ERROR('Error allocating NUMDIR_TBL in urban_param_init')
2123             ALLOCATE( STREET_DIRECTION_TBL(MAXDIRS , ICATE), stat=allocate_status )
2124             if(allocate_status /= 0) FATAL_ERROR('Error allocating STREET_DIRECTION_TBL in urban_param_init')
2125             ALLOCATE( STREET_WIDTH_TBL(MAXDIRS , ICATE), stat=allocate_status )
2126             if(allocate_status /= 0) FATAL_ERROR('Error allocating STREET_WIDTH_TBL in urban_param_init')
2127             ALLOCATE( BUILDING_WIDTH_TBL(MAXDIRS , ICATE), stat=allocate_status )
2128             if(allocate_status /= 0) FATAL_ERROR('Error allocating BUILDING_WIDTH_TBL in urban_param_init')
2129             ALLOCATE( NUMHGT_TBL(ICATE), stat=allocate_status )
2130             if(allocate_status /= 0) FATAL_ERROR('Error allocating NUMHGT_TBL in urban_param_init')
2131             ALLOCATE( HEIGHT_BIN_TBL(MAXHGTS , ICATE), stat=allocate_status )
2132             if(allocate_status /= 0) FATAL_ERROR('Error allocating HEIGHT_BIN_TBL in urban_param_init')
2133             ALLOCATE( HPERCENT_BIN_TBL(MAXHGTS , ICATE), stat=allocate_status )
2134             if(allocate_status /= 0) FATAL_ERROR('Error allocating HPERCENT_BIN_TBL in urban_param_init')
2135             ALLOCATE( COP_TBL(ICATE), stat=allocate_status )
2136             if(allocate_status /= 0) FATAL_ERROR('Error allocating COP_TBL in urban_param_init')
2137             ALLOCATE( BLDAC_FRC_TBL(ICATE), stat=allocate_status )
2138             if(allocate_status /= 0) FATAL_ERROR('Error allocating BLDAC_FRC_TBL in urban_param_init')
2139             ALLOCATE( COOLED_FRC_TBL(ICATE), stat=allocate_status )
2140             if(allocate_status /= 0) FATAL_ERROR('Error allocating COOLED_FRC_TBL in urban_param_init')
2141             ALLOCATE( PWIN_TBL(ICATE), stat=allocate_status )
2142             if(allocate_status /= 0) FATAL_ERROR('Error allocating PWIN_TBL in urban_param_init')
2143             ALLOCATE( BETA_TBL(ICATE), stat=allocate_status )
2144             if(allocate_status /= 0) FATAL_ERROR('Error allocating BETA_TBL in urban_param_init')
2145             ALLOCATE( SW_COND_TBL(ICATE), stat=allocate_status )
2146             if(allocate_status /= 0) FATAL_ERROR('Error allocating SW_COND_TBL in urban_param_init')
2147             ALLOCATE( TIME_ON_TBL(ICATE), stat=allocate_status )
2148             if(allocate_status /= 0) FATAL_ERROR('Error allocating TIME_ON_TBL in urban_param_init')
2149             ALLOCATE( TIME_OFF_TBL(ICATE), stat=allocate_status )
2150             if(allocate_status /= 0) FATAL_ERROR('Error allocating TIME_OFF_TBL in urban_param_init')
2151             ALLOCATE( TARGTEMP_TBL(ICATE), stat=allocate_status )
2152             if(allocate_status /= 0) FATAL_ERROR('Error allocating TARGTEMP_TBL in urban_param_init')
2153             ALLOCATE( GAPTEMP_TBL(ICATE), stat=allocate_status )
2154             if(allocate_status /= 0) FATAL_ERROR('Error allocating GAPTEMP_TBL in urban_param_init')
2155             ALLOCATE( TARGHUM_TBL(ICATE), stat=allocate_status )
2156             if(allocate_status /= 0) FATAL_ERROR('Error allocating TARGHUM_TBL in urban_param_init')
2157             ALLOCATE( GAPHUM_TBL(ICATE), stat=allocate_status )
2158             if(allocate_status /= 0) FATAL_ERROR('Error allocating GAPHUM_TBL in urban_param_init')
2159             ALLOCATE( PERFLO_TBL(ICATE), stat=allocate_status )
2160             if(allocate_status /= 0) FATAL_ERROR('Error allocating PERFLO_TBL in urban_param_init')
2161             ALLOCATE( HSESF_TBL(ICATE), stat=allocate_status )
2162             if(allocate_status /= 0) FATAL_ERROR('Error allocating HSESF_TBL in urban_param_init')
2163             ALLOCATE( PV_FRAC_ROOF_TBL(ICATE), stat=allocate_status )
2164             if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating PV_FRAC_ROOF_TBL in urban_param_init')
2165             ALLOCATE( GR_FRAC_ROOF_TBL(ICATE), stat=allocate_status )
2166             if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating GR_FRAC_ROOF_TBL in urban_param_init')
2167
2168          endif
2169          numdir_tbl = 0
2170          street_direction_tbl = -1.E36
2171          street_width_tbl = 0
2172          building_width_tbl = 0
2173          numhgt_tbl = 0
2174          height_bin_tbl = -1.E36
2175          hpercent_bin_tbl = -1.E36
2176 !end BEP
2177
2178       else if (name == "ZR") then
2179          read(string(indx+1:),*) zr_tbl(1:icate)
2180       else if (name == "SIGMA_ZED") then
2181          read(string(indx+1:),*) sigma_zed_tbl(1:icate)
2182       else if (name == "ROOF_WIDTH") then
2183          ALLOCATE( ROOF_WIDTH(ICATE), stat=allocate_status )
2184          if(allocate_status /= 0) FATAL_ERROR('Error allocating ROOF_WIDTH in urban_param_init')
2185
2186          read(string(indx+1:),*) roof_width(1:icate)
2187       else if (name == "ROAD_WIDTH") then
2188          ALLOCATE( ROAD_WIDTH(ICATE), stat=allocate_status )
2189          if(allocate_status /= 0) FATAL_ERROR('Error allocating ROAD_WIDTH in urban_param_init')
2190          read(string(indx+1:),*) road_width(1:icate)
2191       else if (name == "AH") then
2192          read(string(indx+1:),*) ah_tbl(1:icate)
2193       else if (name == "ALH") then
2194          read(string(indx+1:),*) alh_tbl(1:icate)
2195       else if (name == "FRC_URB") then
2196          read(string(indx+1:),*) frc_urb_tbl(1:icate)
2197       else if (name == "CAPR") then
2198          read(string(indx+1:),*) capr_tbl(1:icate)
2199          ! Convert CAPR_TBL from J m{-3} K{-1} to cal cm{-3} deg{-1}
2200          capr_tbl = capr_tbl * ( 1.0 / 4.1868 ) * 1.E-6
2201       else if (name == "CAPB") then
2202          read(string(indx+1:),*) capb_tbl(1:icate)
2203          ! Convert CABR_TBL from J m{-3} K{-1} to cal cm{-3} deg{-1}
2204          capb_tbl = capb_tbl * ( 1.0 / 4.1868 ) * 1.E-6
2205       else if (name == "CAPG") then
2206          read(string(indx+1:),*) capg_tbl(1:icate)
2207          ! Convert CABG_TBL from J m{-3} K{-1} to cal cm{-3} deg{-1}
2208          capg_tbl = capg_tbl * ( 1.0 / 4.1868 ) * 1.E-6
2209       else if (name == "AKSR") then
2210          read(string(indx+1:),*) aksr_tbl(1:icate)
2211          ! Convert AKSR_TBL from J m{-1} s{-1} K{-1} to cal cm{-1} s{-1} deg{-1}
2212          AKSR_TBL = AKSR_TBL * ( 1.0 / 4.1868 ) * 1.E-2
2213       else if (name == "AKSB") then
2214          read(string(indx+1:),*) aksb_tbl(1:icate)
2215          ! Convert AKSB_TBL from J m{-1} s{-1} K{-1} to cal cm{-1} s{-1} deg{-1}
2216          AKSB_TBL = AKSB_TBL * ( 1.0 / 4.1868 ) * 1.E-2
2217       else if (name == "AKSG") then
2218          read(string(indx+1:),*) aksg_tbl(1:icate)
2219          ! Convert AKSG_TBL from J m{-1} s{-1} K{-1} to cal cm{-1} s{-1} deg{-1}
2220          AKSG_TBL = AKSG_TBL * ( 1.0 / 4.1868 ) * 1.E-2
2221       else if (name == "ALBR") then
2222          read(string(indx+1:),*) albr_tbl(1:icate)
2223       else if (name == "ALBB") then
2224          read(string(indx+1:),*) albb_tbl(1:icate)
2225       else if (name == "ALBG") then
2226          read(string(indx+1:),*) albg_tbl(1:icate)
2227       else if (name == "EPSR") then
2228          read(string(indx+1:),*) epsr_tbl(1:icate)
2229       else if (name == "EPSB") then
2230          read(string(indx+1:),*) epsb_tbl(1:icate)
2231       else if (name == "EPSG") then
2232          read(string(indx+1:),*) epsg_tbl(1:icate)
2233       else if (name == "AKANDA_URBAN") then
2234          read(string(indx+1:),*) akanda_urban_tbl(1:icate)
2235       else if (name == "Z0B") then
2236          read(string(indx+1:),*) z0b_tbl(1:icate)
2237       else if (name == "Z0G") then
2238          read(string(indx+1:),*) z0g_tbl(1:icate)
2239       else if (name == "DDZR") then
2240          read(string(indx+1:),*) dzr(1:num_roof_layers)
2241          ! Convert thicknesses from m to cm
2242          dzr = dzr * 100.0
2243       else if (name == "DDZB") then
2244          read(string(indx+1:),*) dzb(1:num_wall_layers)
2245          ! Convert thicknesses from m to cm
2246          dzb = dzb * 100.0
2247       else if (name == "DDZG") then
2248          read(string(indx+1:),*) dzg(1:num_road_layers)
2249          ! Convert thicknesses from m to cm
2250          dzg = dzg * 100.0
2251       else if (name == "BOUNDR") then
2252          read(string(indx+1:),*) boundr_data
2253       else if (name == "BOUNDB") then
2254          read(string(indx+1:),*) boundb_data
2255       else if (name == "BOUNDG") then
2256          read(string(indx+1:),*) boundg_data
2257       else if (name == "TRLEND") then
2258          read(string(indx+1:),*) trlend_tbl(1:icate)
2259       else if (name == "TBLEND") then
2260          read(string(indx+1:),*) tblend_tbl(1:icate)
2261       else if (name == "TGLEND") then
2262          read(string(indx+1:),*) tglend_tbl(1:icate)
2263       else if (name == "CH_SCHEME") then
2264          read(string(indx+1:),*) ch_scheme_data
2265       else if (name == "TS_SCHEME") then
2266          read(string(indx+1:),*) ts_scheme_data
2267       else if (name == "AHOPTION") then
2268          read(string(indx+1:),*) ahoption
2269       else if (name == "AHDIUPRF") then
2270          read(string(indx+1:),*) ahdiuprf(1:24)
2271       else if (name == "ALHOPTION") then
2272          read(string(indx+1:),*) alhoption
2273       else if (name == "ALHSEASON") then
2274          read(string(indx+1:),*) alhseason(1:4)
2275       else if (name == "ALHDIUPRF") then
2276          read(string(indx+1:),*) alhdiuprf(1:48)
2277       else if (name == "PORIMP") then
2278          read(string(indx+1:),*) porimp(1:3)
2279       else if (name == "DENGIMP") then
2280          read(string(indx+1:),*) dengimp(1:3)
2281       else if (name == "IMP_SCHEME") then
2282          read(string(indx+1:),*) imp_scheme
2283       else if (name == "IRI_SCHEME") then
2284          read(string(indx+1:),*) iri_scheme
2285       else if (name == "OASIS") then
2286          read(string(indx+1:),*) oasis
2287       else if (name == "GROPTION") then
2288          read(string(indx+1:),*) groption
2289       else if (name == "FGR") then
2290          read(string(indx+1:),*) fgr
2291       else if (name == "DZGR") then
2292          read(string(indx+1:),*) dzgr(1:4)
2293 !for BEP
2294       else if (name == "STREET PARAMETERS") then
2295
2296          STREETLOOP : do
2297             read(11,'(A512)', iostat=iostatus) string
2298             if (string(1:1) == "#") cycle STREETLOOP
2299             if (trim(string) == "") cycle STREETLOOP
2300             if (string == "END STREET PARAMETERS") exit STREETLOOP
2301             read(string, *) k ! , dirst, ws, bs
2302             numdir_tbl(k) = numdir_tbl(k) + 1
2303             read(string, *) k, street_direction_tbl(numdir_tbl(k),k), &
2304                                street_width_tbl(numdir_tbl(k),k), &
2305                                building_width_tbl(numdir_tbl(k),k)
2306          enddo STREETLOOP
2307
2308       else if (name == "BUILDING HEIGHTS") then
2309
2310          read(string(indx+1:),*) k
2311          HEIGHTLOOP : do
2312             read(11,'(A512)', iostat=iostatus) string
2313             if (string(1:1) == "#") cycle HEIGHTLOOP
2314             if (trim(string) == "") cycle HEIGHTLOOP
2315             if (string == "END BUILDING HEIGHTS") exit HEIGHTLOOP
2316             read(string,*) dummy_hgt, dummy_pct
2317             numhgt_tbl(k) = numhgt_tbl(k) + 1
2318             height_bin_tbl(numhgt_tbl(k), k) = dummy_hgt
2319             hpercent_bin_tbl(numhgt_tbl(k),k) = dummy_pct
2320
2321          enddo HEIGHTLOOP
2322          pctsum = sum ( hpercent_bin_tbl(:,k) , mask=(hpercent_bin_tbl(:,k)>-1.E25 ) )
2323          if ( pctsum /= 100.) then
2324             write (*,'(//,"Building height percentages for category ", I2, " must sum to 100.0")') k
2325             write (*,'("Currently, they sum to ", F6.2,/)') pctsum
2326             FATAL_ERROR('pctsum is not equal to 100.')
2327          endif
2328       else if ( name == "Z0R") then
2329          read(string(indx+1:),*) Z0R_tbl(1:icate)
2330       else if ( name == "COP") then
2331          read(string(indx+1:),*) cop_tbl(1:icate)
2332       else if ( name == "BLDAC_FRC") then
2333          read(string(indx+1:),*) bldac_frc_tbl(1:icate)
2334       else if ( name == "COOLED_FRC") then
2335          read(string(indx+1:),*) cooled_frc_tbl(1:icate)
2336       else if ( name == "PWIN") then
2337          read(string(indx+1:),*) pwin_tbl(1:icate)
2338       else if ( name == "BETA") then
2339          read(string(indx+1:),*) beta_tbl(1:icate)
2340       else if ( name == "SW_COND") then
2341          read(string(indx+1:),*) sw_cond_tbl(1:icate)
2342       else if ( name == "TIME_ON") then
2343          read(string(indx+1:),*) time_on_tbl(1:icate)
2344       else if ( name == "TIME_OFF") then
2345          read(string(indx+1:),*) time_off_tbl(1:icate)
2346       else if ( name == "TARGTEMP") then
2347          read(string(indx+1:),*) targtemp_tbl(1:icate)
2348       else if ( name == "GAPTEMP") then
2349          read(string(indx+1:),*) gaptemp_tbl(1:icate)
2350       else if ( name == "TARGHUM") then
2351          read(string(indx+1:),*) targhum_tbl(1:icate)
2352       else if ( name == "GAPHUM") then
2353          read(string(indx+1:),*) gaphum_tbl(1:icate)
2354       else if ( name == "PERFLO") then
2355          read(string(indx+1:),*) perflo_tbl(1:icate)
2356       else if (name == "HSEQUIP") then
2357          read(string(indx+1:),*) hsequip_tbl(1:24)
2358       else if (name == "HSEQUIP_SCALE_FACTOR") then
2359          read(string(indx+1:),*) hsesf_tbl(1:icate)
2360       else if (name == "IRHO") then
2361          read(string(indx+1:),*) IRHO_TBL(1:24)
2362        else if ( name == "PV_FRAC_ROOF") then
2363          read(string(indx+1:),*) pv_frac_roof_tbl(1:icate)
2364      else if ( name == "GR_FRAC_ROOF") then
2365          read(string(indx+1:),*) gr_frac_roof_tbl(1:icate)
2366       else if (name == "GR_FLAG") then
2367          read(string(indx+1:),*) gr_flag_tbl
2368       else if (name == "GR_TYPE") then
2369          read(string(indx+1:),*) gr_type_tbl
2370
2371 !end BEP
2372       else
2373          FATAL_ERROR('URBPARM.TBL: Unrecognized NAME = "'//trim(name)//'" in Subr URBAN_PARAM_INIT')
2374       endif
2375    enddo READLOOP
2376
2377    CLOSE(11)
2378
2379    ! Assign a few table values that do not need to come from URBPARM.TBL
2380
2381    Z0HB_TBL = 0.1 * Z0B_TBL
2382    Z0HG_TBL = 0.1 * Z0G_TBL
2383
2384    DO LC = 1, ICATE
2385
2386       ! HGT:  Normalized height
2387       HGT_TBL(LC) = ZR_TBL(LC) / ( ROAD_WIDTH(LC) + ROOF_WIDTH(LC) )
2388
2389       ! R:  Normalized Roof Width (a.k.a. "building coverage ratio")
2390       R_TBL(LC)  = ROOF_WIDTH(LC) / ( ROAD_WIDTH(LC) + ROOF_WIDTH(LC) )
2391
2392       RW_TBL(LC) = 1.0 - R_TBL(LC)
2393       BETR_TBL(LC) = 0.0
2394       BETB_TBL(LC) = 0.0
2395       BETG_TBL(LC) = 0.0
2396
2397       ! The following urban canyon geometry parameters are following Macdonald's (1998) formulations
2398
2399       ! Lambda_P :: Plan areal fraction, which corresponds to R for a 2-d canyon.
2400       ! Lambda_F :: Frontal area index, which corresponds to HGT for a 2-d canyon
2401       ! Cd       :: Drag coefficient ( 1.2 from Grimmond and Oke, 1998 )
2402       ! Alpha_macd :: Emperical coefficient ( 4.43 from Macdonald et al., 1998 )
2403       ! Beta_macd  :: Correction factor for the drag coefficient ( 1.0 from Macdonald et al., 1998 )
2404
2405       Lambda_P = R_TBL(LC)
2406       Lambda_F = HGT_TBL(LC)
2407       Cd         = 1.2
2408       alpha_macd = 4.43
2409       beta_macd  = 1.0
2410
2411
2412       ZDC_TBL(LC) = ZR_TBL(LC) * ( 1.0 + ( alpha_macd ** ( -Lambda_P ) )  * ( Lambda_P - 1.0 ) )
2413
2414       Z0C_TBL(LC) = ZR_TBL(LC) * ( 1.0 - ZDC_TBL(LC)/ZR_TBL(LC) ) * &
2415            exp (-(0.5 * beta_macd * Cd / (VonK**2) * ( 1.0-ZDC_TBL(LC)/ZR_TBL(LC) ) * Lambda_F )**(-0.5))
2416
2417       IF (SF_URBAN_PHYSICS == 1) THEN
2418          ! Include roof height variability in Macdonald
2419          ! to parameterize Z0R as a function of ZR_SD (Standard Deviation)
2420          Lambda_FR  = SIGMA_ZED_TBL(LC) / ( ROAD_WIDTH(LC) + ROOF_WIDTH(LC) )
2421          Z0R_TBL(LC) = ZR_TBL(LC) * ( 1.0 - ZDC_TBL(LC)/ZR_TBL(LC) ) &
2422               * exp ( -(0.5 * beta_macd * Cd / (VonK**2) &
2423               * ( 1.0-ZDC_TBL(LC)/ZR_TBL(LC) ) * Lambda_FR )**(-0.5))
2424       ENDIF
2425
2426       !
2427       ! Z0HC still one-tenth of Z0C, as before ?
2428       !
2429
2430       Z0HC_TBL(LC) = 0.1 * Z0C_TBL(LC)
2431
2432       !
2433       ! Calculate Sky View Factor:
2434       !
2435       DHGT=HGT_TBL(LC)/100.
2436       HGT=0.
2437       VFWS=0.
2438       HGT=HGT_TBL(LC)-DHGT/2.
2439       do k=1,99
2440          HGT=HGT-DHGT
2441          VFWS=VFWS+0.25*(1.-HGT/SQRT(HGT**2.+RW_TBL(LC)**2.))
2442       end do
2443
2444      VFWS=VFWS/99.
2445      VFWS=VFWS*2.
2446
2447      VFGS=1.-2.*VFWS*HGT_TBL(LC)/RW_TBL(LC)
2448      SVF_TBL(LC)=VFGS
2449    END DO
2450
2451    deallocate(roof_width)
2452    deallocate(road_width)
2453
2454    END SUBROUTINE urban_param_init
2455 !===========================================================================
2456 !
2457 ! subroutine urban_var_init: initialization of urban state variables
2458 !
2459 !===========================================================================
2460    SUBROUTINE urban_var_init(ISURBAN, TSURFACE0_URB,TLAYER0_URB,TDEEP0_URB,IVGTYP,  & ! in
2461                              ims,ime,jms,jme,kms,kme,num_soil_layers,         & ! in
2462                              LCZ_1,LCZ_2,LCZ_3,LCZ_4,LCZ_5,                &
2463                              LCZ_6,LCZ_7,LCZ_8,LCZ_9,LCZ_10,LCZ_11,        &
2464                              restart,sf_urban_physics,                     & !in
2465                              XXXR_URB2D,XXXB_URB2D,XXXG_URB2D,XXXC_URB2D,  & ! inout
2466                              TR_URB2D,TB_URB2D,TG_URB2D,TC_URB2D,QC_URB2D, & ! inout
2467                              TRL_URB3D,TBL_URB3D,TGL_URB3D,                & ! inout
2468                              SH_URB2D,LH_URB2D,G_URB2D,RN_URB2D,           & ! inout
2469                              TS_URB2D,                                     & ! inout
2470                              num_urban_ndm,                                & ! in
2471                              urban_map_zrd,                                & ! in
2472                              urban_map_zwd,                                & ! in
2473                              urban_map_gd,                                 & ! in
2474                              urban_map_zd,                                 & ! in
2475                              urban_map_zdf,                                & ! in
2476                              urban_map_bd,                                 & ! in
2477                              urban_map_wd,                                 & ! in
2478                              urban_map_gbd,                                & ! in
2479                              urban_map_fbd,                                & ! in
2480                              urban_map_zgrd,                                & ! in
2481                              num_urban_hi,                                 & ! in
2482                              TRB_URB4D,TW1_URB4D,TW2_URB4D,TGB_URB4D,      & ! inout
2483                              TLEV_URB3D,QLEV_URB3D,                        & ! inout
2484                              TW1LEV_URB3D,TW2LEV_URB3D,                    & ! inout
2485                              TGLEV_URB3D,TFLEV_URB3D,                      & ! inout
2486                              SF_AC_URB3D,LF_AC_URB3D,CM_AC_URB3D,          & ! inout
2487                              SFVENT_URB3D,LFVENT_URB3D,                    & ! inout
2488                              SFWIN1_URB3D,SFWIN2_URB3D,                    & ! inout
2489                              SFW1_URB3D,SFW2_URB3D,SFR_URB3D,SFG_URB3D,    & ! inout
2490                              EP_PV_URB3D,T_PV_URB3D,                       & !GRZ
2491                              TRV_URB4D,QR_URB4D,QGR_URB3D,TGR_URB3D,       & !GRZ
2492                              DRAIN_URB4D,DRAINGR_URB3D,SFRV_URB3D,         & !GRZ
2493                              LFRV_URB3D,DGR_URB3D,DG_URB3D,LFR_URB3D,LFG_URB3D,&!GRZ
2494                              SMOIS_URB,                                    &
2495                              LP_URB2D,HI_URB2D,LB_URB2D,                   & ! inout
2496                              HGT_URB2D,MH_URB2D,STDH_URB2D,                & ! inout
2497                              LF_URB2D,                                     & ! inout
2498                              CMCR_URB2D,TGR_URB2D,TGRL_URB3D,SMR_URB3D,    & ! inout
2499                              DRELR_URB2D,DRELB_URB2D,DRELG_URB2D,          & ! inout
2500                              FLXHUMR_URB2D, FLXHUMB_URB2D, FLXHUMG_URB2D,  & ! inout
2501                              A_U_BEP,A_V_BEP,A_T_BEP,A_Q_BEP,              & ! inout multi-layer urban
2502                              A_E_BEP,B_U_BEP,B_V_BEP,                      & ! inout multi-layer urban
2503                              B_T_BEP,B_Q_BEP,B_E_BEP,DLG_BEP,              & ! inout multi-layer urban
2504                              DL_U_BEP,SF_BEP,VL_BEP,                       & ! inout multi-layer urban
2505                              FRC_URB2D, UTYPE_URB2D,USE_WUDAPT_LCZ)                         ! inout
2506    IMPLICIT NONE
2507
2508    INTEGER, INTENT(IN) :: ISURBAN, sf_urban_physics,use_wudapt_lcz
2509    INTEGER, INTENT(IN) :: LCZ_1,LCZ_2,LCZ_3,LCZ_4,LCZ_5,LCZ_6,LCZ_7,LCZ_8,LCZ_9,LCZ_10,LCZ_11
2510    INTEGER, INTENT(IN) :: ims,ime,jms,jme,kms,kme,num_soil_layers
2511    INTEGER, INTENT(IN) :: num_urban_ndm
2512    INTEGER, INTENT(IN) :: urban_map_zrd
2513    INTEGER, INTENT(IN) :: urban_map_zwd
2514    INTEGER, INTENT(IN) :: urban_map_gd
2515    INTEGER, INTENT(IN) :: urban_map_zd
2516    INTEGER, INTENT(IN) :: urban_map_zdf
2517    INTEGER, INTENT(IN) :: urban_map_bd
2518    INTEGER, INTENT(IN) :: urban_map_wd
2519    INTEGER, INTENT(IN) :: urban_map_gbd
2520    INTEGER, INTENT(IN) :: urban_map_fbd
2521    INTEGER, INTENT(IN) :: urban_map_zgrd
2522    INTEGER, INTENT(IN) :: num_urban_hi                                !multi-layer urban
2523 !   INTEGER, INTENT(IN) :: num_roof_layers, num_wall_layers, num_road_layers
2524
2525    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN)                    :: TSURFACE0_URB
2526    REAL, DIMENSION( ims:ime, 1:num_soil_layers, jms:jme ), INTENT(IN) :: TLAYER0_URB
2527    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN)                    :: TDEEP0_URB
2528    INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(IN)                 :: IVGTYP
2529    LOGICAL , INTENT(IN) :: restart
2530
2531    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TR_URB2D
2532    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TB_URB2D
2533    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TG_URB2D
2534    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TC_URB2D
2535    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: QC_URB2D
2536    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXR_URB2D
2537    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXB_URB2D
2538    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXG_URB2D
2539    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXC_URB2D
2540
2541    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: DRELR_URB2D
2542    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: DRELB_URB2D
2543    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: DRELG_URB2D
2544    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: FLXHUMR_URB2D
2545    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: FLXHUMB_URB2D
2546    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: FLXHUMG_URB2D
2547    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CMCR_URB2D
2548    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TGR_URB2D
2549
2550 !   REAL, DIMENSION(ims:ime, 1:num_roof_layers, jms:jme), INTENT(INOUT) :: TRL_URB3D
2551 !   REAL, DIMENSION(ims:ime, 1:num_wall_layers, jms:jme), INTENT(INOUT) :: TBL_URB3D
2552 !   REAL, DIMENSION(ims:ime, 1:num_road_layers, jms:jme), INTENT(INOUT) :: TGL_URB3D
2553    REAL, DIMENSION(ims:ime, 1:num_soil_layers, jms:jme), INTENT(INOUT) :: TRL_URB3D
2554    REAL, DIMENSION(ims:ime, 1:num_soil_layers, jms:jme), INTENT(INOUT) :: TBL_URB3D
2555    REAL, DIMENSION(ims:ime, 1:num_soil_layers, jms:jme), INTENT(INOUT) :: TGL_URB3D
2556    REAL, DIMENSION(ims:ime, 1:num_soil_layers, jms:jme), INTENT(INOUT) :: TGRL_URB3D
2557    REAL, DIMENSION(ims:ime, 1:num_soil_layers, jms:jme), INTENT(INOUT) :: SMR_URB3D
2558
2559    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: SH_URB2D
2560    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: LH_URB2D
2561    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: G_URB2D
2562    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: RN_URB2D
2563    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TS_URB2D
2564
2565 ! multi-layer UCM variables
2566    REAL, DIMENSION(ims:ime, 1:urban_map_zrd, jms:jme), INTENT(INOUT) :: TRB_URB4D
2567    REAL, DIMENSION(ims:ime, 1:urban_map_zwd, jms:jme), INTENT(INOUT) :: TW1_URB4D
2568    REAL, DIMENSION(ims:ime, 1:urban_map_zwd, jms:jme), INTENT(INOUT) :: TW2_URB4D
2569    REAL, DIMENSION(ims:ime, 1:urban_map_gd , jms:jme), INTENT(INOUT) :: TGB_URB4D
2570    REAL, DIMENSION(ims:ime, 1:urban_map_bd , jms:jme), INTENT(INOUT) :: TLEV_URB3D
2571    REAL, DIMENSION(ims:ime, 1:urban_map_bd , jms:jme), INTENT(INOUT) :: QLEV_URB3D
2572    REAL, DIMENSION(ims:ime, 1:urban_map_wd , jms:jme), INTENT(INOUT) :: TW1LEV_URB3D
2573    REAL, DIMENSION(ims:ime, 1:urban_map_wd , jms:jme), INTENT(INOUT) :: TW2LEV_URB3D
2574    REAL, DIMENSION(ims:ime, 1:urban_map_gbd, jms:jme), INTENT(INOUT) :: TGLEV_URB3D
2575    REAL, DIMENSION(ims:ime, 1:urban_map_fbd, jms:jme), INTENT(INOUT) :: TFLEV_URB3D
2576    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: LF_AC_URB3D
2577    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: SF_AC_URB3D
2578    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CM_AC_URB3D
2579    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: SFVENT_URB3D
2580    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: LFVENT_URB3D
2581    REAL, DIMENSION( ims:ime, 1:urban_map_wd, jms:jme), INTENT(INOUT) :: SFWIN1_URB3D
2582    REAL, DIMENSION( ims:ime, 1:urban_map_wd, jms:jme), INTENT(INOUT) :: SFWIN2_URB3D
2583    REAL, DIMENSION(ims:ime, 1:urban_map_zd , jms:jme), INTENT(INOUT) :: SFW1_URB3D
2584    REAL, DIMENSION(ims:ime, 1:urban_map_zd , jms:jme), INTENT(INOUT) :: SFW2_URB3D
2585    REAL, DIMENSION(ims:ime, 1:urban_map_zdf, jms:jme), INTENT(INOUT) :: SFR_URB3D
2586    REAL, DIMENSION(ims:ime, 1:num_urban_ndm, jms:jme), INTENT(INOUT) :: SFG_URB3D
2587    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: EP_PV_URB3D!GRZ
2588    REAL, DIMENSION( ims:ime, 1:urban_map_zdf,jms:jme ), INTENT(INOUT) :: T_PV_URB3D!GRZ
2589    REAL, DIMENSION( ims:ime, 1:urban_map_zgrd, jms:jme),INTENT(INOUT) :: TRV_URB4D ! GRZ
2590    REAL, DIMENSION( ims:ime, 1:urban_map_zgrd, jms:jme),INTENT(INOUT) :: QR_URB4D ! GRZ
2591    REAL, DIMENSION( ims:ime,jms:jme), INTENT(INOUT) :: QGR_URB3D ! GRZ
2592    REAL, DIMENSION( ims:ime,jms:jme), INTENT(INOUT) :: TGR_URB3D ! GRZ
2593    REAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme),INTENT(INOUT) :: DRAIN_URB4D !GRZ
2594    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: DRAINGR_URB3D   !GRZ
2595    REAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme),INTENT(INOUT) :: SFRV_URB3D !GRZ
2596    REAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme),INTENT(INOUT) :: LFRV_URB3D ! GRZ
2597    REAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme ),INTENT(INOUT) :: DGR_URB3D !GRZ
2598    REAL, DIMENSION( ims:ime, 1:num_urban_ndm, jms:jme ),INTENT(INOUT) :: DG_URB3D !GRZ
2599    REAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme ),INTENT(INOUT) :: LFR_URB3D !GRZ
2600    REAL, DIMENSION( ims:ime, 1:num_urban_ndm, jms:jme ), INTENT(INOUT) :: LFG_URB3D !GRZ
2601    REAL, DIMENSION( ims:ime, 1:num_soil_layers, jms:jme ), INTENT(IN) ::SMOIS_URB
2602    REAL, DIMENSION( ims:ime,1:num_urban_hi , jms:jme), INTENT(INOUT) :: HI_URB2D
2603    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: LP_URB2D
2604    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: LB_URB2D
2605    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: HGT_URB2D
2606    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: MH_URB2D
2607    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: STDH_URB2D
2608    REAL, DIMENSION( ims:ime, 4,jms:jme ), INTENT(INOUT) :: LF_URB2D
2609    REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: A_U_BEP
2610    REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: A_V_BEP
2611    REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: A_T_BEP
2612    REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: A_Q_BEP
2613    REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: A_E_BEP
2614    REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: B_U_BEP
2615    REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: B_V_BEP
2616    REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: B_T_BEP
2617    REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: B_Q_BEP
2618    REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: B_E_BEP
2619    REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: VL_BEP
2620    REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: DLG_BEP
2621    REAL, DIMENSION(ims:ime, kms:kme,jms:jme),INTENT(INOUT) :: SF_BEP
2622    REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: DL_U_BEP
2623 !
2624    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: FRC_URB2D
2625    INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: UTYPE_URB2D
2626    INTEGER                                            :: UTYPE_URB
2627 !FS
2628    INTEGER :: SWITCH_URB
2629
2630    INTEGER :: I,J,K,CHECK
2631
2632    CHECK = 0
2633
2634    DO I=ims,ime
2635     DO J=jms,jme
2636
2637 !      XXXR_URB2D(I,J)=0.
2638 !      XXXB_URB2D(I,J)=0.
2639 !      XXXG_URB2D(I,J)=0.
2640 !      XXXC_URB2D(I,J)=0.
2641
2642       SH_URB2D(I,J)=0.
2643       LH_URB2D(I,J)=0.
2644       G_URB2D(I,J)=0.
2645       RN_URB2D(I,J)=0.
2646
2647 !m
2648 !FS   FRC_URB2D(I,J)=0.
2649       UTYPE_URB2D(I,J)=0
2650       SWITCH_URB=1
2651
2652             IF( IVGTYP(I,J) == ISURBAN)  THEN
2653                IF(use_wudapt_lcz==0) THEN
2654                   UTYPE_URB2D(I,J) = 2  ! for default. high-intensity
2655                ELSE
2656                   UTYPE_URB2D(I,J) = 5  ! for default. high-intensity
2657                ENDIF
2658             ELSE IF( IVGTYP(I,J) == LCZ_1) THEN
2659                UTYPE_URB2D(I,J) = 1
2660             ELSE IF( IVGTYP(I,J) == LCZ_2) THEN
2661                UTYPE_URB2D(I,J) = 2
2662             ELSE IF( IVGTYP(I,J) == LCZ_3) THEN
2663                UTYPE_URB2D(I,J) = 3
2664             ELSE IF( IVGTYP(I,J) == LCZ_4) THEN
2665                UTYPE_URB2D(I,J) = 4
2666             ELSE IF( IVGTYP(I,J) == LCZ_5) THEN
2667                UTYPE_URB2D(I,J) = 5
2668             ELSE IF( IVGTYP(I,J) == LCZ_6) THEN
2669                UTYPE_URB2D(I,J) = 6
2670             ELSE IF( IVGTYP(I,J) == LCZ_7) THEN
2671                UTYPE_URB2D(I,J) = 7
2672             ELSE IF( IVGTYP(I,J) == LCZ_8) THEN
2673                UTYPE_URB2D(I,J) = 8
2674             ELSE IF( IVGTYP(I,J) == LCZ_9) THEN
2675                UTYPE_URB2D(I,J) = 9
2676             ELSE IF( IVGTYP(I,J) == LCZ_10) THEN
2677                UTYPE_URB2D(I,J) = 10
2678             ELSE IF( IVGTYP(I,J) == LCZ_11) THEN
2679                UTYPE_URB2D(I,J) = 11
2680             ELSE
2681                SWITCH_URB=0
2682             ENDIF
2683
2684             IF (SWITCH_URB == 1) THEN
2685                UTYPE_URB = UTYPE_URB2D(I,J)  ! for default. high-intensity
2686                IF (HGT_URB2D(I,J)>0.) THEN
2687                   CONTINUE
2688                ELSE
2689                   WRITE(mesg,*) 'USING DEFAULT URBAN MORPHOLOGY'
2690                   WRITE_MESSAGE(mesg)
2691                   LP_URB2D(I,J)=0.
2692                   LB_URB2D(I,J)=0.
2693                   HGT_URB2D(I,J)=0.
2694                   IF ( sf_urban_physics == 1 ) THEN
2695                      MH_URB2D(I,J)=0.
2696                      STDH_URB2D(I,J)=0.
2697                      DO K=1,4
2698                        LF_URB2D(I,K,J)=0.
2699                      ENDDO
2700                   ELSE IF ( ( sf_urban_physics == 2 ) .or. ( sf_urban_physics == 3 ) ) THEN
2701                      DO K=1,num_urban_hi
2702                         HI_URB2D(I,K,J)=0.
2703                      ENDDO
2704                   ENDIF
2705                ENDIF
2706                 IF (FRC_URB2D(I,J)>0.and.FRC_URB2D(I,J)<=1.) THEN
2707                   CONTINUE
2708                 ELSE
2709                   WRITE(mesg,*) 'WARNING, FRC_URB2D = 0 BUT IVGTYP IS URBAN'
2710                   WRITE_MESSAGE(mesg)
2711                   WRITE(mesg,*) 'WARNING, THE URBAN FRACTION WILL BE READ FROM URBPARM.TBL'
2712                   WRITE_MESSAGE(mesg)
2713                   FRC_URB2D(I,J) = FRC_URB_TBL(UTYPE_URB)
2714                 ENDIF
2715             ELSE
2716                 FRC_URB2D(I,J)=0.
2717                 LP_URB2D(I,J)=0.
2718                 LB_URB2D(I,J)=0.
2719                 HGT_URB2D(I,J)=0.
2720                 IF ( sf_urban_physics == 1 ) THEN
2721                    MH_URB2D(I,J)=0.
2722                    STDH_URB2D(I,J)=0.
2723                    DO K=1,4
2724                       LF_URB2D(I,K,J)=0.
2725                    ENDDO
2726                 ELSE IF ( ( sf_urban_physics == 2 ) .or. ( sf_urban_physics == 3 ) ) THEN
2727                    DO K=1,num_urban_hi
2728                       HI_URB2D(I,K,J)=0.
2729                    ENDDO
2730                 ENDIF
2731             ENDIF
2732
2733
2734       QC_URB2D(I,J)=0.01
2735
2736        IF (.not.restart) THEN
2737
2738       XXXR_URB2D(I,J)=0.
2739       XXXB_URB2D(I,J)=0.
2740       XXXG_URB2D(I,J)=0.
2741       XXXC_URB2D(I,J)=0.
2742
2743       IF ( sf_urban_physics == 1 ) THEN
2744       DRELR_URB2D(I,J) = 0.
2745       DRELB_URB2D(I,J) = 0.
2746       DRELG_URB2D(I,J) = 0.
2747       FLXHUMR_URB2D(I,J) = 0.
2748       FLXHUMB_URB2D(I,J) = 0.
2749       FLXHUMG_URB2D(I,J) = 0.
2750       CMCR_URB2D(I,J)  = 0.
2751       TGR_URB2D(I,J)=TSURFACE0_URB(I,J)+0.
2752       ENDIF
2753
2754       TC_URB2D(I,J)=TSURFACE0_URB(I,J)+0.
2755       TR_URB2D(I,J)=TSURFACE0_URB(I,J)+0.
2756       TB_URB2D(I,J)=TSURFACE0_URB(I,J)+0.
2757       TG_URB2D(I,J)=TSURFACE0_URB(I,J)+0.
2758 !
2759       TS_URB2D(I,J)=TSURFACE0_URB(I,J)+0.
2760
2761 !     DO K=1,num_roof_layers
2762 !     DO K=1,num_soil_layers
2763 !         TRL_URB3D(I,1,J)=TLAYER0_URB(I,1,J)+0.
2764 !         TRL_URB3D(I,2,J)=TLAYER0_URB(I,2,J)+0.
2765 !         TRL_URB3D(I,3,J)=TLAYER0_URB(I,3,J)+0.
2766 !         TRL_URB3D(I,4,J)=TLAYER0_URB(I,4,J)+0.
2767
2768           TRL_URB3D(I,1,J)=TLAYER0_URB(I,1,J)+0.
2769           TRL_URB3D(I,2,J)=0.5*(TLAYER0_URB(I,1,J)+TLAYER0_URB(I,2,J))
2770           TRL_URB3D(I,3,J)=TLAYER0_URB(I,2,J)+0.
2771           TRL_URB3D(I,4,J)=TLAYER0_URB(I,2,J)+(TLAYER0_URB(I,3,J)-TLAYER0_URB(I,2,J))*0.29
2772
2773       IF ( sf_urban_physics == 1 ) THEN
2774           TGRL_URB3D(I,1,J)=TLAYER0_URB(I,1,J)+0.
2775           TGRL_URB3D(I,2,J)=0.5*(TLAYER0_URB(I,1,J)+TLAYER0_URB(I,2,J))
2776           TGRL_URB3D(I,3,J)=TLAYER0_URB(I,2,J)+0.
2777           TGRL_URB3D(I,4,J)=TLAYER0_URB(I,2,J)+(TLAYER0_URB(I,3,J)-TLAYER0_URB(I,2,J))*0.29
2778
2779           SMR_URB3D(I,1,J)=0.2
2780           SMR_URB3D(I,2,J)=0.2
2781           SMR_URB3D(I,3,J)=0.2
2782           SMR_URB3D(I,4,J)=0.
2783       ENDIF
2784 !     END DO
2785
2786 !      DO K=1,num_wall_layers
2787 !      DO K=1,num_soil_layers
2788 !m        TBL_URB3D(I,1,J)=TLAYER0_URB(I,1,J)+0.
2789 !m        TBL_URB3D(I,2,J)=TLAYER0_URB(I,2,J)+0.
2790 !m        TBL_URB3D(I,3,J)=TLAYER0_URB(I,3,J)+0.
2791 !m        TBL_URB3D(I,4,J)=TLAYER0_URB(I,4,J)+0.
2792
2793         TBL_URB3D(I,1,J)=TLAYER0_URB(I,1,J)+0.
2794         TBL_URB3D(I,2,J)=0.5*(TLAYER0_URB(I,1,J)+TLAYER0_URB(I,2,J))
2795         TBL_URB3D(I,3,J)=TLAYER0_URB(I,2,J)+0.
2796         TBL_URB3D(I,4,J)=TLAYER0_URB(I,2,J)+(TLAYER0_URB(I,3,J)-TLAYER0_URB(I,2,J))*0.29
2797
2798 !      END DO
2799
2800 !      DO K=1,num_road_layers
2801       DO K=1,num_soil_layers
2802         TGL_URB3D(I,K,J)=TLAYER0_URB(I,K,J)+0.
2803       END DO
2804
2805 ! multi-layer urban
2806        IF((SF_URBAN_PHYSICS.eq.2).OR.(SF_URBAN_PHYSICS.eq.3)) THEN
2807         IF (UTYPE_URB2D(I,J) > 0) THEN
2808            TRB_URB4D(I,:,J)=TBLEND_TBL(UTYPE_URB2D(I,J))
2809            TW1_URB4D(I,:,J)=TBLEND_TBL(UTYPE_URB2D(I,J))
2810            TW2_URB4D(I,:,J)=TBLEND_TBL(UTYPE_URB2D(I,J))
2811         ELSE
2812            TRB_URB4D(I,:,J)=tlayer0_urb(I,1,J)
2813            TW1_URB4D(I,:,J)=tlayer0_urb(I,1,J)
2814            TW2_URB4D(I,:,J)=tlayer0_urb(I,1,J)
2815         ENDIF
2816         TGB_URB4D(I,:,J)=tlayer0_urb(I,1,J)
2817         SFW1_URB3D(I,:,J)=0.
2818         SFW2_URB3D(I,:,J)=0.
2819         SFR_URB3D(I,:,J)=0.
2820         SFG_URB3D(I,:,J)=0.
2821
2822        ENDIF
2823
2824       if (SF_URBAN_PHYSICS.EQ.3) then
2825          LF_AC_URB3D(I,J)=0.
2826          SF_AC_URB3D(I,J)=0.
2827          CM_AC_URB3D(I,J)=0.
2828          SFVENT_URB3D(I,J)=0.
2829          LFVENT_URB3D(I,J)=0.
2830          EP_PV_URB3D(I,J)=0.
2831          T_PV_URB3D(I,:,J)=tlayer0_urb(I,1,J)
2832          TGR_URB3D(I,J)=tlayer0_urb(I,1,J)
2833          QR_URB4D(I,:,J)=SMOIS_URB(I,1,J)
2834          DRAIN_URB4D(I,:,J)=0. !GRZ
2835          SFRV_URB3D(I,:,J)=0.                !GRZ
2836          LFRV_URB3D(I,:,J)=0.                !GRZ
2837          DGR_URB3D(I,:,J)=0.  !GRZ
2838          DG_URB3D(I,:,J)=0.
2839          LFR_URB3D(I,:,J)=0.
2840          LFG_URB3D(I,:,J)=0.
2841          QGR_URB3D(I,J)=SMOIS_URB(I,1,J)  !GRZ
2842          TGR_URB3D(I,J)=0.
2843          DRAINGR_URB3D(I,J)=0. !GRZ
2844
2845        IF (UTYPE_URB2D(I,J) > 0) THEN
2846         TRV_URB4D(I,:,J)=TBLEND_TBL(UTYPE_URB2D(I,J))
2847        ELSE
2848        TRV_URB4D(I,:,J)=tlayer0_urb(I,1,J) !GRZ
2849       ENDIF
2850
2851             TLEV_URB3D(I,:,J)=tlayer0_urb(I,1,J)
2852             TW1LEV_URB3D(I,:,J)=tlayer0_urb(I,1,J)
2853             TW2LEV_URB3D(I,:,J)=tlayer0_urb(I,1,J)
2854             TGLEV_URB3D(I,:,J)=tlayer0_urb(I,1,J)
2855             TFLEV_URB3D(I,:,J)=tlayer0_urb(I,1,J)
2856             QLEV_URB3D(I,:,J)=0.01
2857             SFWIN1_URB3D(I,:,J)=0.
2858             SFWIN2_URB3D(I,:,J)=0.
2859 !rm         LF_AC_URB3D(I,J)=0.
2860 !rm         SF_AC_URB3D(I,J)=0.
2861 !rm         CM_AC_URB3D(I,J)=0.
2862 !rm         SFVENT_URB3D(I,J)=0.
2863 !rm         LFVENT_URB3D(I,J)=0.
2864
2865       endif
2866
2867 !      IF( sf_urban_physics .EQ. 2 )THEN
2868       IF((SF_URBAN_PHYSICS.eq.2).OR.(SF_URBAN_PHYSICS.eq.3)) THEN
2869       DO K= KMS,KME
2870           SF_BEP(I,K,J)=1.
2871           VL_BEP(I,K,J)=1.
2872           A_U_BEP(I,K,J)=0.
2873           A_V_BEP(I,K,J)=0.
2874           A_T_BEP(I,K,J)=0.
2875           A_E_BEP(I,K,J)=0.
2876           A_Q_BEP(I,K,J)=0.
2877           B_U_BEP(I,K,J)=0.
2878           B_V_BEP(I,K,J)=0.
2879           B_T_BEP(I,K,J)=0.
2880           B_E_BEP(I,K,J)=0.
2881           B_Q_BEP(I,K,J)=0.
2882           DLG_BEP(I,K,J)=0.
2883           DL_U_BEP(I,K,J)=0.
2884       END DO
2885       ENDIF       !sf_urban_physics=2
2886      ENDIF        !restart
2887
2888
2889       IF (CHECK.EQ.0)THEN
2890       IF(IVGTYP(I,J).EQ.1)THEN
2891         write(mesg,*) 'Sample of Urban settings'
2892         call wrf_message(mesg)
2893         write(mesg,*) 'TSURFACE0_URB',TSURFACE0_URB(I,J)
2894         WRITE_MESSAGE(mesg)
2895         write(mesg,*) 'TDEEP0_URB', TDEEP0_URB(I,J)
2896         WRITE_MESSAGE(mesg)
2897         write(mesg,*) 'IVGTYP',IVGTYP(I,J)
2898         WRITE_MESSAGE(mesg)
2899         write(mesg,*) 'TR_URB2D',TR_URB2D(I,J)
2900         WRITE_MESSAGE(mesg)
2901         write(mesg,*) 'TB_URB2D',TB_URB2D(I,J)
2902         WRITE_MESSAGE(mesg)
2903         write(mesg,*) 'TG_URB2D',TG_URB2D(I,J)
2904         WRITE_MESSAGE(mesg)
2905         write(mesg,*) 'TC_URB2D',TC_URB2D(I,J)
2906         WRITE_MESSAGE(mesg)
2907         write(mesg,*) 'QC_URB2D',QC_URB2D(I,J)
2908         WRITE_MESSAGE(mesg)
2909         write(mesg,*) 'XXXR_URB2D',XXXR_URB2D(I,J)
2910         WRITE_MESSAGE(mesg)
2911         write(mesg,*) 'SH_URB2D',SH_URB2D(I,J)
2912         WRITE_MESSAGE(mesg)
2913         write(mesg,*) 'LH_URB2D',LH_URB2D(I,J)
2914         WRITE_MESSAGE(mesg)
2915         write(mesg,*) 'G_URB2D',G_URB2D(I,J)
2916         WRITE_MESSAGE(mesg)
2917         write(mesg,*) 'RN_URB2D',RN_URB2D(I,J)
2918         WRITE_MESSAGE(mesg)
2919         write(mesg,*) 'TS_URB2D',TS_URB2D(I,J)
2920         WRITE_MESSAGE(mesg)
2921         write(mesg,*) 'LF_AC_URB3D', LF_AC_URB3D(I,J)
2922         WRITE_MESSAGE(mesg)
2923         write(mesg,*) 'SF_AC_URB3D', SF_AC_URB3D(I,J)
2924         WRITE_MESSAGE(mesg)
2925         write(mesg,*) 'CM_AC_URB3D', CM_AC_URB3D(I,J)
2926         WRITE_MESSAGE(mesg)
2927         write(mesg,*) 'SFVENT_URB3D', SFVENT_URB3D(I,J)
2928         WRITE_MESSAGE(mesg)
2929         write(mesg,*) 'LFVENT_URB3D', LFVENT_URB3D(I,J)
2930         WRITE_MESSAGE(mesg)
2931         write(mesg,*) 'FRC_URB2D', FRC_URB2D(I,J)
2932         WRITE_MESSAGE(mesg)
2933         write(mesg,*) 'UTYPE_URB2D', UTYPE_URB2D(I,J)
2934         WRITE_MESSAGE(mesg)
2935         write(mesg,*) 'I',I,'J',J
2936         WRITE_MESSAGE(mesg)
2937         write(mesg,*) 'num_urban_hi', num_urban_hi
2938         WRITE_MESSAGE(mesg)
2939         CHECK = 1
2940        END IF
2941        END IF
2942
2943     END DO
2944    END DO
2945    RETURN
2946    END SUBROUTINE urban_var_init
2947 !===========================================================================
2948 !
2949 !  force_restore
2950 !
2951 !===========================================================================
2952    SUBROUTINE force_restore(CAP,AKS,DELT,S,R,H,LE,TSLEND,TSP,TS)
2953
2954      REAL, INTENT(IN)  :: CAP,AKS,DELT,S,R,H,LE,TSLEND,TSP
2955      REAL, INTENT(OUT) :: TS
2956      REAL              :: C1,C2
2957
2958      C2=24.*3600./2./3.14159
2959      C1=SQRT(0.5*C2*CAP*AKS)
2960
2961      TS = TSP + DELT*( (S+R-H-LE)/C1 -(TSP-TSLEND)/C2 )
2962
2963    END SUBROUTINE force_restore
2964 !===========================================================================
2965 !
2966 !  bisection (not used)
2967 !
2968 !==============================================================================
2969    SUBROUTINE bisection(TSP,PS,S,EPS,RX,SIG,RHO,CP,CH,UA,QA,TA,EL,BET,AKS,TSL,DZ,TS)
2970
2971      REAL, INTENT(IN) :: TSP,PS,S,EPS,RX,SIG,RHO,CP,CH,UA,QA,TA,EL,BET,AKS,TSL,DZ
2972      REAL, INTENT(OUT) :: TS
2973      REAL :: ES,QS0,R,H,ELE,G0,F1,F
2974
2975      TS1 = TSP - 5.
2976      TS2 = TSP + 5.
2977
2978      DO ITERATION = 1,22
2979
2980        ES=6.11*EXP( (2.5*10.**6./461.51)*(TS1-273.15)/(273.15*TS1) )
2981        QS0=0.622*ES/(PS-0.378*ES)
2982        R=EPS*(RX-SIG*(TS1**4.)/60.)
2983        H=RHO*CP*CH*UA*(TS1-TA)*100.
2984        ELE=RHO*EL*CH*UA*BET*(QS0-QA)*100.
2985        G0=AKS*(TS1-TSL)/(DZ/2.)
2986        F1= S + R - H - ELE - G0
2987
2988        TS=0.5*(TS1+TS2)
2989
2990        ES=6.11*EXP( (2.5*10.**6./461.51)*(TS-273.15)/(273.15*TS) )
2991        QS0=0.622*ES/(PS-0.378*ES)
2992        R=EPS*(RX-SIG*(TS**4.)/60.)
2993        H=RHO*CP*CH*UA*(TS-TA)*100.
2994        ELE=RHO*EL*CH*UA*BET*(QS0-QA)*100.
2995        G0=AKS*(TS-TSL)/(DZ/2.)
2996        F = S + R - H - ELE - G0
2997
2998        IF (F1*F > 0.0) THEN
2999          TS1=TS
3000         ELSE
3001          TS2=TS
3002        END IF
3003
3004      END DO
3005
3006      RETURN
3007 END SUBROUTINE bisection
3008 !===========================================================================
3009
3010 SUBROUTINE SFCDIF_URB (ZLM,Z0,THZ0,THLM,SFCSPD,AKANDA,AKMS,AKHS,RLMO,CD)
3011
3012 ! ----------------------------------------------------------------------
3013 ! SUBROUTINE SFCDIF_URB (Urban version of SFCDIF_off)
3014 ! ----------------------------------------------------------------------
3015 ! CALCULATE SURFACE LAYER EXCHANGE COEFFICIENTS VIA ITERATIVE PROCESS.
3016 ! SEE CHEN ET AL (1997, BLM)
3017 ! ----------------------------------------------------------------------
3018
3019       IMPLICIT NONE
3020       REAL     WWST, WWST2, G, VKRM, EXCM, BETA, BTG, ELFC, WOLD, WNEW
3021       REAL     PIHF, EPSU2, EPSUST, EPSIT, EPSA, ZTMIN, ZTMAX, HPBL,     &
3022      & SQVISC
3023       REAL     RIC, RRIC, FHNEU, RFC,RLMO_THR, RFAC, ZZ, PSLMU, PSLMS, PSLHU,     &
3024      & PSLHS
3025       REAL     XX, PSPMU, YY, PSPMS, PSPHU, PSPHS, ZLM, Z0, THZ0, THLM
3026       REAL     SFCSPD, AKANDA, AKMS, AKHS, ZU, ZT, RDZ, CXCH
3027       REAL     DTHV, DU2, BTGH, WSTAR2, USTAR, ZSLU, ZSLT, RLOGU, RLOGT
3028       REAL     RLMO, ZETALT, ZETALU, ZETAU, ZETAT, XLU4, XLT4, XU4, XT4
3029 !CC   ......REAL ZTFC
3030
3031       REAL     XLU, XLT, XU, XT, PSMZ, SIMM, PSHZ, SIMH, USTARK, RLMN,  &
3032      &         RLMA
3033
3034       INTEGER  ITRMX, ILECH, ITR
3035       REAL,    INTENT(OUT) :: CD
3036       PARAMETER                                                         &
3037      &        (WWST = 1.2,WWST2 = WWST * WWST,G = 9.8,VKRM = 0.40,      &
3038      &         EXCM = 0.001                                             &
3039      &        ,BETA = 1./270.,BTG = BETA * G,ELFC = VKRM * BTG          &
3040      &                  ,WOLD =.15,WNEW = 1. - WOLD,ITRMX = 05,         &
3041      &                   PIHF = 3.14159265/2.)
3042       PARAMETER                                                         &
3043      &         (EPSU2 = 1.E-4,EPSUST = 0.07,EPSIT = 1.E-4,EPSA = 1.E-8  &
3044      &         ,ZTMIN = -5.,ZTMAX = 1.,HPBL = 1000.0                    &
3045      &          ,SQVISC = 258.2)
3046       PARAMETER                                                         &
3047      &       (RIC = 0.183,RRIC = 1.0/ RIC,FHNEU = 0.8,RFC = 0.191       &
3048      &        ,RLMO_THR = 0.001,RFAC = RIC / (FHNEU * RFC * RFC))
3049
3050 ! ----------------------------------------------------------------------
3051 ! NOTE: THE TWO CODE BLOCKS BELOW DEFINE FUNCTIONS
3052 ! ----------------------------------------------------------------------
3053 ! LECH'S SURFACE FUNCTIONS
3054 ! ----------------------------------------------------------------------
3055       PSLMU (ZZ)= -0.96* log (1.0-4.5* ZZ)
3056       PSLMS (ZZ)= ZZ * RRIC -2.076* (1. -1./ (ZZ +1.))
3057       PSLHU (ZZ)= -0.96* log (1.0-4.5* ZZ)
3058
3059 ! ----------------------------------------------------------------------
3060 ! PAULSON'S SURFACE FUNCTIONS
3061 ! ----------------------------------------------------------------------
3062       PSLHS (ZZ)= ZZ * RFAC -2.076* (1. -1./ (ZZ +1.))
3063       PSPMU (XX)= -2.* log ( (XX +1.)*0.5) - log ( (XX * XX +1.)*0.5)   &
3064      &        +2.* ATAN (XX)                                            &
3065      &- PIHF
3066       PSPMS (YY)= 5.* YY
3067       PSPHU (XX)= -2.* log ( (XX * XX +1.)*0.5)
3068
3069 ! ----------------------------------------------------------------------
3070 ! THIS ROUTINE SFCDIF CAN HANDLE BOTH OVER OPEN WATER (SEA, OCEAN) AND
3071 ! OVER SOLID SURFACE (LAND, SEA-ICE).
3072 ! ----------------------------------------------------------------------
3073       PSPHS (YY)= 5.* YY
3074
3075 ! ----------------------------------------------------------------------
3076 !     ZTFC: RATIO OF ZOH/ZOM  LESS OR EQUAL THAN 1
3077 !     C......ZTFC=0.1
3078 !     CZIL: CONSTANT C IN Zilitinkevich, S. S.1995,:NOTE ABOUT ZT
3079 ! ----------------------------------------------------------------------
3080       ILECH = 0
3081
3082 ! ----------------------------------------------------------------------
3083 !      ZILFC = - CZIL * VKRM * SQVISC
3084 !     C.......ZT=Z0*ZTFC
3085       ZU = Z0
3086       RDZ = 1./ ZLM
3087       CXCH = EXCM * RDZ
3088       DTHV = THLM - THZ0
3089
3090 ! ----------------------------------------------------------------------
3091 ! BELJARS CORRECTION OF USTAR
3092 ! ----------------------------------------------------------------------
3093       DU2 = MAX (SFCSPD * SFCSPD,EPSU2)
3094 !cc   If statements to avoid TANGENT LINEAR problems near zero
3095       BTGH = BTG * HPBL
3096       IF (BTGH * AKHS * DTHV .ne. 0.0) THEN
3097          WSTAR2 = WWST2* ABS (BTGH * AKHS * DTHV)** (2./3.)
3098       ELSE
3099          WSTAR2 = 0.0
3100       END IF
3101
3102 ! ----------------------------------------------------------------------
3103 ! ZILITINKEVITCH APPROACH FOR ZT
3104 ! ----------------------------------------------------------------------
3105       USTAR = MAX (SQRT (AKMS * SQRT (DU2+ WSTAR2)),EPSUST)
3106
3107 ! ----------------------------------------------------------------------
3108 ! KCL/TL Try Kanda approach instead (Kanda et al. 2007, JAMC)
3109 !      ZT = EXP (ZILFC * SQRT (USTAR * Z0))* Z0
3110       ZT = EXP (2.0-AKANDA*(SQVISC**2 * USTAR * Z0)**0.25)* Z0
3111
3112       ZSLU = ZLM + ZU
3113
3114       ZSLT = ZLM + ZT
3115       RLOGU = log (ZSLU / ZU)
3116
3117       RLOGT = log (ZSLT / ZT)
3118
3119       RLMO = ELFC * AKHS * DTHV / USTAR **3
3120 ! ----------------------------------------------------------------------
3121 ! 1./MONIN-OBUKKHOV LENGTH-SCALE
3122 ! ----------------------------------------------------------------------
3123       DO ITR = 1,ITRMX
3124          ZETALT = MAX (ZSLT * RLMO,ZTMIN)
3125          RLMO = ZETALT / ZSLT
3126          ZETALU = ZSLU * RLMO
3127          ZETAU = ZU * RLMO
3128
3129          ZETAT = ZT * RLMO
3130          IF (ILECH .eq. 0) THEN
3131             IF (RLMO .lt. 0.0)THEN
3132                XLU4 = 1. -16.* ZETALU
3133                XLT4 = 1. -16.* ZETALT
3134                XU4 = 1. -16.* ZETAU
3135
3136                XT4 = 1. -16.* ZETAT
3137                XLU = SQRT (SQRT (XLU4))
3138                XLT = SQRT (SQRT (XLT4))
3139                XU = SQRT (SQRT (XU4))
3140
3141                XT = SQRT (SQRT (XT4))
3142
3143                PSMZ = PSPMU (XU)
3144                SIMM = PSPMU (XLU) - PSMZ + RLOGU
3145                PSHZ = PSPHU (XT)
3146                SIMH = PSPHU (XLT) - PSHZ + RLOGT
3147             ELSE
3148                ZETALU = MIN (ZETALU,ZTMAX)
3149                ZETALT = MIN (ZETALT,ZTMAX)
3150                PSMZ = PSPMS (ZETAU)
3151                SIMM = PSPMS (ZETALU) - PSMZ + RLOGU
3152                PSHZ = PSPHS (ZETAT)
3153                SIMH = PSPHS (ZETALT) - PSHZ + RLOGT
3154             END IF
3155 ! ----------------------------------------------------------------------
3156 ! LECH'S FUNCTIONS
3157 ! ----------------------------------------------------------------------
3158          ELSE
3159             IF (RLMO .lt. 0.)THEN
3160                PSMZ = PSLMU (ZETAU)
3161                SIMM = PSLMU (ZETALU) - PSMZ + RLOGU
3162                PSHZ = PSLHU (ZETAT)
3163                SIMH = PSLHU (ZETALT) - PSHZ + RLOGT
3164             ELSE
3165                ZETALU = MIN (ZETALU,ZTMAX)
3166                ZETALT = MIN (ZETALT,ZTMAX)
3167                PSMZ = PSLMS (ZETAU)
3168                SIMM = PSLMS (ZETALU) - PSMZ + RLOGU
3169                PSHZ = PSLHS (ZETAT)
3170                SIMH = PSLHS (ZETALT) - PSHZ + RLOGT
3171             END IF
3172 ! ----------------------------------------------------------------------
3173 ! BELJAARS CORRECTION FOR USTAR
3174 ! ----------------------------------------------------------------------
3175          END IF
3176             USTAR = MAX (SQRT (AKMS * SQRT (DU2+ WSTAR2)),EPSUST)
3177             !KCL/TL
3178             !ZT = EXP (ZILFC * SQRT (USTAR * Z0))* Z0
3179             ZT = EXP (2.0-AKANDA*(SQVISC**2 * USTAR * Z0)**0.25)* Z0
3180             ZSLT = ZLM + ZT
3181             RLOGT = log (ZSLT / ZT)
3182             USTARK = USTAR * VKRM
3183             AKMS = MAX (USTARK / SIMM,CXCH)
3184             AKHS = MAX (USTARK / SIMH,CXCH)
3185 !
3186          IF (BTGH * AKHS * DTHV .ne. 0.0) THEN
3187             WSTAR2 = WWST2* ABS (BTGH * AKHS * DTHV)** (2./3.)
3188          ELSE
3189             WSTAR2 = 0.0
3190          END IF
3191 !-----------------------------------------------------------------------
3192          RLMN = ELFC * AKHS * DTHV / USTAR **3
3193 !-----------------------------------------------------------------------
3194 !     IF(ABS((RLMN-RLMO)/RLMA).LT.EPSIT)    GO TO 110
3195 !-----------------------------------------------------------------------
3196          RLMA = RLMO * WOLD+ RLMN * WNEW
3197 !-----------------------------------------------------------------------
3198          RLMO = RLMA
3199
3200       END DO
3201
3202       CD = USTAR*USTAR/SFCSPD**2
3203 ! ----------------------------------------------------------------------
3204   END SUBROUTINE SFCDIF_URB
3205 ! ----------------------------------------------------------------------
3206 !===========================================================================
3207 !  DIREVAP
3208 !  CALCULATE DIRECT SOIL EVAPORATION
3209 !===========================================================================
3210    SUBROUTINE DIREVAP (EDIR,ETP,SMC,SHDFAC,SMCMAX,SMCDRY,FXEXP)
3211
3212      REAL, INTENT(IN)  :: ETP,SMC,SHDFAC,SMCMAX,SMCDRY,FXEXP
3213      REAL, INTENT(OUT) :: EDIR
3214      REAL              :: FX, SRATIO
3215
3216 ! ----------------------------------------------------------------------
3217 ! FX > 1 REPRESENTS DEMAND CONTROL
3218 ! FX < 1 REPRESENTS FLUX CONTROL
3219 ! ----------------------------------------------------------------------
3220       SRATIO = (SMC - SMCDRY) / (SMCMAX - SMCDRY)
3221       IF (SRATIO > 0.) THEN
3222         FX = SRATIO**FXEXP
3223         FX = MAX ( MIN ( FX, 1. ) ,0. )
3224       ELSE
3225         FX = 0.
3226       ENDIF
3227       EDIR = FX * ( 1.0- SHDFAC ) * ETP * 0.001
3228
3229  END SUBROUTINE DIREVAP
3230 !===========================================================================
3231 !  TRANSP
3232 !  CALCULATE EVAPOTRANSPIRATION FOR VEGETATIO SURFACE
3233 !===========================================================================
3234
3235      SUBROUTINE TRANSP (ETT,ET,EC,SHDFAC,ETP1,CMC,CFACTR,CMCMAX,LAI,RSMIN,RSMAX,RGL,SX, &
3236                         TS,TA,QA,SMC,SMCWLT,SMCREF,CPP,PS,CH,EPSV,DELT, NROOT,NSOIL,    &
3237                         DZVR, ZSOIL, HS)
3238      INTEGER, INTENT(IN)   :: NROOT, NSOIL
3239      REAL, INTENT(IN)  :: SHDFAC,ETP1,CMC,CFACTR,CMCMAX,LAI,RSMIN,RSMAX,RGL,SX,TA
3240      REAL, INTENT(IN)  :: TS,QA, SMCWLT, SMCREF, CPP, PS,CH, EPSV, DELT, HS
3241      REAL, DIMENSION(1:NSOIL), INTENT(IN)   :: ZSOIL, DZVR, SMC
3242      REAL, DIMENSION(1:NSOIL), INTENT(INOUT):: ET
3243      REAL, INTENT(OUT) :: EC, ETT
3244      REAL              :: RC, RCS, RCT, RCQ, RCSOIL, FF, WS, SLV, DESDT
3245      REAL              :: SIGMA, PC, CMC2MS, SGX, DENOM, RTX, ETT1
3246      INTEGER           :: K
3247      REAL, DIMENSION(1:NROOT) ::  PART, GX
3248
3249      SLV    = 2.501E+6
3250      SIGMA  = 5.67E-8
3251      ETT    = 0.0
3252      DO K = 1, NSOIL
3253        ET(K) = 0.
3254      END DO
3255
3256 ! ----------------------------------------------------------------------
3257 ! INITIALIZE CANOPY RESISTANCE MULTIPLIER TERMS.
3258 ! ----------------------------------------------------------------------
3259       RCS = 0.0
3260       RCT = 0.0
3261       RCQ = 0.0
3262       RCSOIL = 0.0
3263
3264 ! ----------------------------------------------------------------------
3265 ! CONTRIBUTION DUE TO INCOMING SOLAR RADIATION
3266 ! ----------------------------------------------------------------------
3267       FF  = 0.55*2.0* SX*697.7 * 60/ (RGL * LAI)
3268       RCS = (FF + RSMIN / RSMAX) / (1.0+ FF)
3269       RCS = MAX (RCS,0.0001)
3270 ! ----------------------------------------------------------------------
3271 ! CONTRIBUTION DUE TO AIR TEMPERATURE AT FIRST MODEL LEVEL ABOVE GROUND
3272 ! RCT EXPRESSION FROM NOILHAN AND PLANTON (1989, MWR).
3273 ! ----------------------------------------------------------------------
3274       RCT = 1.0- 0.0016* ( (298 - TA)**2.0)
3275       RCT = MAX (RCT,0.0001)
3276 ! ----------------------------------------------------------------------
3277 ! CONTRIBUTION DUE TO VAPOR PRESSURE DEFICIT AT FIRST MODEL LEVEL.
3278 ! RCQ EXPRESSION FROM SSIB  (Niyogi and Raman, 1997)
3279 ! ----------------------------------------------------------------------
3280       EA = 6.11*EXP((2.5*10.**6./461.51)*(TA-273.15)/(273.15*TA) )
3281       WS = 0.622*EA/1013
3282       RCQ = 1.0/ (1.0+ HS * (WS - QA))
3283       RCQ = MAX (RCQ,0.01)
3284 ! ----------------------------------------------------------------------
3285 ! CONTRIBUTION DUE TO SOIL MOISTURE AVAILABILITY.
3286 ! DETERMINE CONTRIBUTION FROM EACH SOIL LAYER, THEN ADD THEM UP.
3287 ! ----------------------------------------------------------------------
3288       DO K = 1, NROOT
3289        GX(K) = (SMC(K) - SMCWLT) / (SMCREF - SMCWLT)
3290          IF (GX(K) >  1.) GX(K) = 1.
3291          IF (GX(K) <  0.) GX(K) = 0.
3292          PART (K) = ( -DZVR (K)/ ZSOIL (3)) * GX(K)
3293       END DO
3294
3295       SGX =0.0
3296       DO K = 1, NROOT
3297         SGX    = SGX    + GX (K)
3298         RCSOIL = RCSOIL + PART (K)
3299       END DO
3300       SGX =SGX / NROOT
3301
3302       RCSOIL = MAX (RCSOIL,0.0001)
3303
3304       RC = RSMIN / (LAI * RCS * RCT * RCQ * RCSOIL)
3305       DESDT = 0.622*SLV*EA/461.51/TA/TA/1013
3306       DELTA = (SLV / CPP)* DESDT
3307       RR = (4.* EPSV *SIGMA * 287.04 / CPP)* (TA **4.)/ (TS * CH) + 1.0
3308       PC = (RR + DELTA)/ (RR * (1. + RC * CH) + DELTA)
3309
3310       IF (CMC .ne. 0.0) THEN
3311          ETT1 = SHDFAC * PC * ETP1 * (1.0- (CMC / CMCMAX) ** CFACTR) * 0.001
3312       ELSE
3313          ETT1 = SHDFAC * PC * ETP1 * 0.001
3314       ENDIF
3315
3316       DENOM = 0.
3317       DO K = 1, NROOT
3318          RTX= (-DZVR (K)/ ZSOIL (3)) + GX(K) - SGX
3319          GX (K) = GX (K) * MAX ( RTX, 0. )
3320          DENOM  = DENOM + GX (K)
3321       END DO
3322       IF (DENOM .le. 0.0) DENOM =1.
3323
3324       DO K = 1, NROOT
3325          ET(K) = ETT1 * GX (K) / DENOM
3326          ETT   = ETT + ET (K)
3327       END DO
3328
3329
3330       IF (CMC > 0.0) THEN
3331       EC = SHDFAC * ( ( CMC / CMCMAX ) ** CFACTR ) * ETP1 * 0.001
3332       ELSE
3333       EC = 0.0
3334       END IF
3335       CMC2MS = CMC / DELT
3336       EC   = MIN ( CMC2MS, EC )
3337
3338   END SUBROUTINE TRANSP
3339 ! ----------------------------------------------------------------------
3340 ! SUBROUTINE SMFLX
3341 ! ----------------------------------------------------------------------
3342
3343        SUBROUTINE SMFLX (SMCP,SMC,NSOIL,CMCP,CMC,DT,PRCP1,ZSOIL,       &
3344      &                   SMCMAX,BEXP,SMCWLT,DKSAT,DWSAT,               &
3345      &                   SHDFAC,CMCMAX,RUNOFF1,RUNOFF2,RUNOFF3,        &
3346                          EDIR,EC,ET,DRIP)
3347
3348 ! CALCULATE SOIL MOISTURE FLUX.  THE SOIL MOISTURE CONTENT IS UPDATED WITH
3349 ! PROGNOSTIC EQNS. THE CANOPY MOISTURE CONTENT (CMC) IS ALSO UPDATED.
3350 ! ----------------------------------------------------------------------
3351       IMPLICIT NONE
3352
3353       INTEGER, INTENT(IN)   :: NSOIL
3354       INTEGER               :: I,K
3355       REAL, INTENT(IN)      :: BEXP, CMCMAX, DKSAT,DWSAT, DT, EC, EDIR,  &
3356                                PRCP1, SHDFAC, SMCMAX, SMCWLT
3357       REAL, INTENT(OUT)     :: DRIP, RUNOFF1, RUNOFF2, RUNOFF3
3358       REAL, INTENT(IN)      :: CMCP
3359       REAL, INTENT(OUT)     :: CMC
3360       REAL, DIMENSION(1:NSOIL), INTENT(IN)   :: ZSOIL, ET
3361       REAL, DIMENSION(1:NSOIL), INTENT(IN)   :: SMCP
3362       REAL, DIMENSION(1:NSOIL), INTENT(OUT)  :: SMC
3363       REAL, DIMENSION(1:NSOIL)               :: AI, BI, CI, STCF,RHSTS, RHSTT
3364       REAL                  :: EXCESS,PCPDRP,RHSCT,TRHSCT
3365
3366
3367 ! ----------------------------------------------------------------------
3368 ! ADD PRECIPITATION TO EXISTING CMC.IF RESULTING AMT EXCEEDS MAX CAPACITY,
3369 ! IT BECOMES DRIP AND WILL FALL TO THE GRND.
3370 ! ----------------------------------------------------------------------
3371       RHSCT = SHDFAC * PRCP1 * 0.001 /3600. - EC
3372       DRIP = 0.
3373       TRHSCT = DT * RHSCT
3374       EXCESS = CMCP + TRHSCT
3375
3376 ! ----------------------------------------------------------------------
3377 ! PCPDRP IS THE COMBINED PRCP1 AND DRIP (FROM CMCP) THAT GOES INTO THE
3378 ! SOIL
3379 ! ----------------------------------------------------------------------
3380       IF (EXCESS > CMCMAX) DRIP = EXCESS - CMCMAX
3381       PCPDRP = (1. - SHDFAC) * PRCP1 * 0.001 /3600. + DRIP / DT
3382
3383 ! ----------------------------------------------------------------------
3384 ! CALL SUBROUTINES SRT AND SSTEP TO SOLVE THE SOIL MOISTURE
3385 ! TENDENCY EQUATIONS.
3386 ! ----------------------------------------------------------------------
3387       CALL SRT (RHSTT,EDIR,ET,SMCP,NSOIL,PCPDRP,ZSOIL,DWSAT,DKSAT,    &
3388                    SMCMAX,BEXP,RUNOFF1,RUNOFF2,DT,SMCWLT,AI,BI,CI)
3389
3390       CALL SSTEP (SMCP,SMC,CMCP,CMC,RHSTT,RHSCT,DT,NSOIL,SMCMAX,     &
3391                         CMCMAX,RUNOFF3,ZSOIL,AI,BI,CI)
3392 ! ----------------------------------------------------------------------
3393   END SUBROUTINE SMFLX
3394 ! ----------------------------------------------------------------------
3395
3396       SUBROUTINE SRT (RHSTT,EDIR,ET,SMCP,NSOIL,PCPDRP,ZSOIL,DWSAT,    &
3397                       DKSAT,SMCMAX,BEXP,RUNOFF1,            &
3398                       RUNOFF2,DT,SMCWLT,AI,BI,CI)
3399
3400 ! ----------------------------------------------------------------------
3401 ! CALCULATE THE RIGHT HAND SIDE OF THE TIME TENDENCY TERM OF THE SOIL
3402 ! WATER DIFFUSION EQUATION.  ALSO TO COMPUTE ( PREPARE ) THE MATRIX
3403 ! COEFFICIENTS FOR THE TRI-DIAGONAL MATRIX OF THE IMPLICIT TIME SCHEME.
3404 ! ----------------------------------------------------------------------
3405       IMPLICIT NONE
3406       INTEGER, INTENT(IN)       :: NSOIL
3407       INTEGER                   :: K, KS
3408
3409       REAL, INTENT(IN)          :: BEXP, DKSAT, DT, DWSAT, EDIR,  &
3410                                    PCPDRP, SMCMAX, SMCWLT
3411       REAL, INTENT(OUT)         :: RUNOFF1, RUNOFF2
3412       REAL, DIMENSION(1:NSOIL), INTENT(IN)   :: SMCP, ZSOIL, ET
3413       REAL, DIMENSION(1:NSOIL), INTENT(OUT)  :: RHSTT
3414       REAL, DIMENSION(1:NSOIL), INTENT(OUT)  :: AI, BI, CI
3415       REAL, DIMENSION(1:NSOIL)  :: DDMAX
3416       REAL                      :: DD, DDT, DDZ, DDZ2, DENOM,       &
3417                                    DENOM2, DSMDZ, DSMDZ2, DT1,      &
3418                                    INFMAX,MXSMC,MXSMC2,NUMER,PDDUM, &
3419                                    PX,SMCAV, SSTT, PAR,         &
3420                                    VAL, WCND, WCND2, WDF, WDF2,KDT
3421
3422 ! ----------------------------------------------------------------------
3423 ! DETERMINE RAINFALL INFILTRATION RATE AND RUNOFF.  INCLUDE THE
3424 ! INFILTRATION FORMULE FROM SCHAAKE AND KOREN MODEL.
3425 ! MODIFIED BY Q DUAN
3426 ! ----------------------------------------------------------------------
3427
3428          PDDUM = PCPDRP
3429          RUNOFF1 = 0.0
3430          PAR = 2.0E-6
3431
3432          IF (PCPDRP /=  0.0) THEN
3433          SMCAV = SMCMAX - SMCWLT
3434          DDMAX (1) = - ZSOIL (1)* SMCAV
3435          DDMAX (1) = DDMAX (1)* (1.0- (SMCP (1) - SMCWLT)/ SMCAV)
3436          DDMAX (2) = (ZSOIL (1) - ZSOIL (2))* SMCAV
3437          DDMAX (2) = DDMAX (2)* (1.0- (SMCP (2) - SMCWLT)/ SMCAV)
3438          DDMAX (3) = (ZSOIL (2) - ZSOIL (3))* SMCAV
3439          DDMAX (3) = DDMAX (3)* (1.0- (SMCP (3) - SMCWLT)/ SMCAV)
3440
3441          DD = DDMAX(1)+DDMAX(2)+DDMAX(3)
3442          DT1 = DT/86400
3443          KDT = 3.0 * DKSAT / PAR
3444          VAL = (1. - EXP ( - KDT * DT1))
3445          DDT = DD * VAL
3446          PX = PCPDRP * DT
3447          IF (PX <  0.0) PX = 0.0
3448
3449          INFMAX = (PX * (DDT / (PX + DDT)))/ DT
3450          MXSMC = SMCP (1)
3451          CALL WDFCND (WDF,WCND,MXSMC,SMCMAX,BEXP,DKSAT,DWSAT)
3452          INFMAX = MAX (INFMAX,WCND)
3453          INFMAX = MIN (INFMAX,PX/DT)
3454
3455
3456          IF (PCPDRP >  INFMAX) THEN
3457           RUNOFF1  = PCPDRP - INFMAX
3458           PDDUM = INFMAX
3459          END IF
3460         END IF
3461 ! ----------------------------------------------------------------------
3462 ! TOP LAYER
3463 ! ----------------------------------------------------------------------
3464       CALL WDFCND (WDF,WCND,SMCP(1),SMCMAX,BEXP,DKSAT,DWSAT)
3465       DDZ = 1. / ( - .5 * ZSOIL (2) )
3466       AI (1) = 0.0
3467       BI (1) = WDF * DDZ / ( - ZSOIL (1) )
3468       CI (1) = - BI (1)
3469       DSMDZ = (SMCP (1) - SMCP (2) )/( - 0.5 * ZSOIL(2))
3470       RHSTT (1) = (WDF * DSMDZ + WCND- PDDUM + EDIR + ET(1))/ ZSOIL (1)
3471       SSTT = WDF * DSMDZ + WCND+ EDIR + ET(1)
3472
3473 ! ----------------------------------------------------------------------
3474 ! LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABV PROCESS
3475 ! ----------------------------------------------------------------------
3476       DDZ2 = 0.0
3477       DO K = 2,NSOIL-1
3478          DENOM2 = (ZSOIL (K -1) - ZSOIL (K))
3479          IF (K /= NSOIL-1) THEN
3480             MXSMC2 = SMCP (K)
3481             CALL WDFCND (WDF2,WCND2,MXSMC2,SMCMAX,BEXP,DKSAT,DWSAT)
3482             DENOM = (ZSOIL (K -1) - ZSOIL (K +1))
3483             DSMDZ2 = (SMCP (K) - SMCP (K +1)) / (DENOM * 0.5)
3484             DDZ2 = 2.0 / DENOM
3485             CI (K) = - WDF2 * DDZ2 / DENOM2
3486          ELSE
3487           CALL WDFCND (WDF2,WCND2,SMCP(NSOIL-1),SMCMAX,BEXP,DKSAT,DWSAT)
3488             DSMDZ2 = 0.0
3489             CI (K) = 0.0
3490          END IF
3491          NUMER = (WDF2 * DSMDZ2) - (WDF * DSMDZ)         &
3492                  - WCND+ ET(K)
3493          RHSTT (K) = NUMER / ( - DENOM2)
3494          AI (K) = - WDF * DDZ / DENOM2
3495          BI (K) = - ( AI (K) + CI (K) )
3496          IF (K .eq. NSOIL-1) THEN
3497             RUNOFF2 = 0.0
3498          END IF
3499          IF (K .ne. NSOIL-1) THEN
3500             WDF = WDF2
3501             WCND = WCND2
3502             DSMDZ = DSMDZ2
3503             DDZ = DDZ2
3504          END IF
3505       END DO
3506 ! ----------------------------------------------------------------------
3507   END SUBROUTINE SRT
3508 ! ----------------------------------------------------------------------
3509
3510       SUBROUTINE SSTEP (SMCP,SMC,CMCP,CMC,RHSTT,RHSCT,DT,     &
3511                         NSOIL,SMCMAX,CMCMAX,RUNOFF3,ZSOIL,        &
3512                         AI,BI,CI)
3513
3514 ! ----------------------------------------------------------------------
3515 ! SUBROUTINE SSTEP
3516 ! ----------------------------------------------------------------------
3517 ! CALCULATE/UPDATE SOIL MOISTURE CONTENT VALUES AND CANOPY MOISTURE
3518 ! CONTENT VALUES.
3519 ! ----------------------------------------------------------------------
3520       IMPLICIT NONE
3521       INTEGER, INTENT(IN)       :: NSOIL
3522       INTEGER                   :: I, K, KK11
3523
3524       REAL, INTENT(IN)          :: CMCMAX, DT, SMCMAX
3525       REAL, INTENT(OUT)         :: RUNOFF3
3526       REAL, INTENT(IN)          :: CMCP
3527       REAL, INTENT(OUT)         :: CMC
3528       REAL, DIMENSION(1:NSOIL), INTENT(IN)     :: SMCP, ZSOIL
3529       REAL, DIMENSION(1:NSOIL), INTENT(OUT)    :: SMC
3530       REAL, DIMENSION(1:NSOIL), INTENT(INOUT)  :: RHSTT
3531       REAL, DIMENSION(1:NSOIL), INTENT(INOUT)  :: AI, BI, CI
3532       REAL, DIMENSION(1:NSOIL)  :: RHSTTin, SMCOUT,SMCIN
3533       REAL, DIMENSION(1:NSOIL)  :: CIin
3534       REAL                      :: DDZ, RHSCT, WPLUS, STOT
3535
3536 ! ----------------------------------------------------------------------
3537 ! CREATE 'AMOUNT' VALUES OF VARIABLES TO BE INPUT TO THE
3538 ! TRI-DIAGONAL MATRIX ROUTINE.
3539 ! ----------------------------------------------------------------------
3540       DO K = 1,NSOIL-1
3541          RHSTT (K) = RHSTT (K) * DT
3542          AI (K) = AI (K) * DT
3543          BI (K) = 1. + BI (K) * DT
3544          CI (K) = CI (K) * DT
3545       END DO
3546 ! ----------------------------------------------------------------------
3547 ! COPY VALUES FOR INPUT VARIABLES BEFORE CALL TO ROSR12
3548 ! ----------------------------------------------------------------------
3549       DO K = 1,NSOIL-1
3550          RHSTTin (K) = RHSTT (K)
3551       END DO
3552       DO K = 1,NSOIL-1
3553          CIin (K) = CI (K)
3554       END DO
3555 ! ----------------------------------------------------------------------
3556 ! CALL ROSR12 TO SOLVE THE TRI-DIAGONAL MATRIX
3557 ! ----------------------------------------------------------------------
3558       CALL ROSR12 (CI,AI,BI,CIin,RHSTTin,RHSTT,NSOIL-1)
3559 ! ----------------------------------------------------------------------
3560 ! SUM THE PREVIOUS SMC VALUE AND THE MATRIX SOLUTION TO GET A
3561 ! NEW VALUE.  MIN ALLOWABLE VALUE OF SMC WILL BE 0.02.
3562 ! RUNOFF3: RUNOFF WITHIN SOIL LAYERS
3563 ! ----------------------------------------------------------------------
3564       WPLUS = 0.0
3565       RUNOFF3 = 0.
3566
3567       DDZ = - ZSOIL (1)
3568       DO K = 1,NSOIL-1
3569          IF (K /= 1) DDZ = ZSOIL (K - 1) - ZSOIL (K)
3570          SMCOUT (K) = SMCP (K) + CI (K) + WPLUS / DDZ
3571          STOT = SMCOUT (K)
3572          IF (STOT > SMCMAX) THEN
3573             IF (K .eq. 1) THEN
3574                DDZ = - ZSOIL (1)
3575             ELSE
3576                KK11 = K - 1
3577                DDZ = - ZSOIL (K) + ZSOIL (KK11)
3578             END IF
3579             WPLUS = (STOT - SMCMAX) * DDZ
3580          ELSE
3581             WPLUS = 0.
3582          END IF
3583          SMC (K) = MAX ( MIN (STOT,SMCMAX),0.066 )
3584       END DO
3585
3586 ! ----------------------------------------------------------------------
3587 ! UPDATE CANOPY WATER CONTENT/INTERCEPTION (CMC).  CONVERT RHSCT TO
3588 ! AN 'AMOUNT' VALUE AND ADD TO PREVIOUS CMC VALUE TO GET NEW CMC.
3589 ! ----------------------------------------------------------------------
3590       RUNOFF3 = WPLUS
3591       CMC = CMCP + DT * RHSCT
3592       IF (CMC < 1.E-20) CMC = 0.0
3593       CMC = MIN (CMC,CMCMAX)
3594
3595 ! ----------------------------------------------------------------------
3596   END SUBROUTINE SSTEP
3597 ! ----------------------------------------------------------------------
3598
3599       SUBROUTINE WDFCND (WDF,WCND,SMC,SMCMAX,BEXP,DKSAT,DWSAT)
3600
3601 ! ----------------------------------------------------------------------
3602 ! SUBROUTINE WDFCND
3603 ! ----------------------------------------------------------------------
3604 ! CALCULATE SOIL WATER DIFFUSIVITY AND SOIL HYDRAULIC CONDUCTIVITY.
3605 ! ----------------------------------------------------------------------
3606       IMPLICIT NONE
3607       REAL     BEXP
3608       REAL     DKSAT
3609       REAL     DWSAT
3610       REAL     EXPON
3611       REAL     FACTR1
3612       REAL     FACTR2
3613       REAL     SMC
3614       REAL     SMCMAX
3615       REAL     WCND
3616
3617 ! ----------------------------------------------------------------------
3618 !     CALC THE RATIO OF THE ACTUAL TO THE MAX PSBL SOIL H2O CONTENT
3619 ! ----------------------------------------------------------------------
3620       REAL     WDF
3621       FACTR1 = 0.05 / SMCMAX
3622
3623 ! ----------------------------------------------------------------------
3624 ! PREP AN EXPNTL COEF AND CALC THE SOIL WATER DIFFUSIVITY AND CONDUCTIVITY
3625 ! ----------------------------------------------------------------------
3626       FACTR2 = SMC / SMCMAX
3627       FACTR1 = MIN(FACTR1,FACTR2)
3628       EXPON  = BEXP + 2.0
3629       WDF    = DWSAT * FACTR2 ** EXPON
3630       EXPON  = (2.0 * BEXP) + 3.0
3631       WCND   = DKSAT * FACTR2 ** EXPON
3632
3633 ! ----------------------------------------------------------------------
3634   END SUBROUTINE WDFCND
3635 ! ----------------------------------------------------------------------
3636 ! SUBROUTINE ROSR12
3637 ! ----------------------------------------------------------------------
3638 ! INVERT (SOLVE) THE TRI-DIAGONAL MATRIX PROBLEM SHOWN BELOW:
3639 ! ###                                            ### ###  ###   ###  ###
3640 ! #B(1), C(1),  0  ,  0  ,  0  ,   . . .  ,    0   # #      #   #      #
3641 ! #A(2), B(2), C(2),  0  ,  0  ,   . . .  ,    0   # #      #   #      #
3642 ! # 0  , A(3), B(3), C(3),  0  ,   . . .  ,    0   # #      #   # D(3) #
3643 ! # 0  ,  0  , A(4), B(4), C(4),   . . .  ,    0   # # P(4) #   # D(4) #
3644 ! # 0  ,  0  ,  0  , A(5), B(5),   . . .  ,    0   # # P(5) #   # D(5) #
3645 ! # .                                          .   # #  .   # = #   .  #
3646 ! # .                                          .   # #  .   #   #   .  #
3647 ! # .                                          .   # #  .   #   #   .  #
3648 ! # 0  , . . . , 0 , A(M-2), B(M-2), C(M-2),   0   # #P(M-2)#   #D(M-2)#
3649 ! # 0  , . . . , 0 ,   0   , A(M-1), B(M-1), C(M-1)# #P(M-1)#   #D(M-1)#
3650 ! # 0  , . . . , 0 ,   0   ,   0   ,  A(M) ,  B(M) # # P(M) #   # D(M) #
3651 ! ###                                            ### ###  ###   ###  ###
3652 ! ----------------------------------------------------------------------
3653
3654   SUBROUTINE ROSR12 (P,A,B,C,D,DELTA,NSOIL)
3655       IMPLICIT NONE
3656
3657       INTEGER, INTENT(IN)   :: NSOIL
3658       INTEGER               :: K, KK
3659
3660       REAL, DIMENSION(1:NSOIL), INTENT(IN):: A, B, D
3661       REAL, DIMENSION(1:NSOIL),INTENT(INOUT):: C,P,DELTA
3662
3663 ! ----------------------------------------------------------------------
3664 ! INITIALIZE EQN COEF C FOR THE LOWEST SOIL LAYER
3665 ! ----------------------------------------------------------------------
3666       C (NSOIL) = 0.0
3667       P (1) = - C (1) / B (1)
3668       DELTA (1) = D (1) / B (1)
3669       DO K = 2,NSOIL
3670          P (K) = - C (K) * ( 1.0 / (B (K) + A (K) * P (K -1)) )
3671          DELTA (K) = (D (K) - A (K)* DELTA (K -1))* (1.0/ (B (K) + A (K)&
3672                     * P (K -1)))
3673       END DO
3674 ! ----------------------------------------------------------------------
3675 ! SET P TO DELTA FOR LOWEST SOIL LAYER
3676 ! ----------------------------------------------------------------------
3677       P (NSOIL) = DELTA (NSOIL)
3678
3679 ! ----------------------------------------------------------------------
3680 ! ADJUST P FOR SOIL LAYERS 2 THRU NSOIL
3681 ! ----------------------------------------------------------------------
3682       DO K = 2,NSOIL
3683          KK = NSOIL - K + 1
3684          P (KK) = P (KK) * P (KK +1) + DELTA (KK)
3685       END DO
3686 ! ----------------------------------------------------------------------
3687   END SUBROUTINE ROSR12
3688 !----------------------------------------------------------------------
3689
3690       SUBROUTINE SHFLX (SSOIL,STC,SMC,SMCMAX,NSOIL,T1,DT,YY,ZZ1,ZSOIL, &
3691                          TBOT,ZBOT,SMCWLT,DF1,QUARTZ,CSOIL,CAPR)
3692
3693 ! ----------------------------------------------------------------------
3694 ! SUBROUTINE SHFLX
3695 ! ----------------------------------------------------------------------
3696 ! UPDATE THE TEMPERATURE STATE OF THE SOIL COLUMN BASED ON THE THERMAL
3697 ! DIFFUSION EQUATION AND UPDATE THE FROZEN SOIL MOISTURE CONTENT BASED
3698 ! ON THE TEMPERATURE.
3699 ! ----------------------------------------------------------------------
3700       IMPLICIT NONE
3701
3702       INTEGER, INTENT(IN)   :: NSOIL
3703       INTEGER               :: I
3704
3705       REAL, INTENT(IN)      :: DF1,DT,SMCMAX, SMCWLT, TBOT,YY, ZBOT,ZZ1, QUARTZ
3706       REAL, INTENT(IN)      :: CSOIL, CAPR
3707       REAL, INTENT(INOUT)   :: T1
3708       REAL, INTENT(OUT)     :: SSOIL
3709       REAL, DIMENSION(1:NSOIL), INTENT(IN)    :: SMC,ZSOIL
3710       REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: STC
3711       REAL, DIMENSION(1:NSOIL)             :: AI, BI, CI, STCF,RHSTS
3712
3713 ! ----------------------------------------------------------------------
3714 ! HRT ROUTINE CALCS THE RIGHT HAND SIDE OF THE SOIL TEMP DIF EQN
3715 ! ----------------------------------------------------------------------
3716
3717       ! Land case
3718
3719       CALL HRT (RHSTS,STC,SMC,SMCMAX,NSOIL,ZSOIL,YY,ZZ1,TBOT,        &
3720                 ZBOT,DT,DF1,AI,BI,CI,QUARTZ,CSOIL,CAPR)
3721
3722       CALL HSTEP (STCF,STC,RHSTS,DT,NSOIL,AI,BI,CI)
3723
3724       DO I = 1,NSOIL
3725          STC (I) = STCF (I)
3726       ENDDO
3727
3728 ! ----------------------------------------------------------------------
3729 ! CALCULATE SURFACE SOIL HEAT FLUX
3730 ! ----------------------------------------------------------------------
3731       T1 = (YY + (ZZ1- 1.0) * STC (1)) / ZZ1
3732       SSOIL = DF1 * (STC (1) - T1) / (0.5 * ZSOIL (1))
3733
3734 ! ----------------------------------------------------------------------
3735   END SUBROUTINE SHFLX
3736 ! ----------------------------------------------------------------------
3737 ! SUBROUTINE HRT
3738 ! ----------------------------------------------------------------------
3739 ! CALCULATE THE RIGHT HAND SIDE OF THE TIME TENDENCY TERM OF THE SOIL
3740 ! THERMAL DIFFUSION EQUATION.  ALSO TO COMPUTE ( PREPARE ) THE MATRIX
3741 ! COEFFICIENTS FOR THE TRI-DIAGONAL MATRIX OF THE IMPLICIT TIME SCHEME.
3742 ! ----------------------------------------------------------------------
3743
3744       SUBROUTINE HRT (RHSTS,STC,SMC,SMCMAX,NSOIL,ZSOIL,YY,ZZ1,          &
3745                        TBOT,ZBOT,DT,DF1,AI,BI,CI,QUARTZ,CSOIL,CAPR)
3746
3747       IMPLICIT NONE
3748       LOGICAL              :: ITAVG
3749       INTEGER, INTENT(IN)  :: NSOIL
3750       INTEGER              :: I, K
3751
3752       REAL, INTENT(IN)     :: DF1, DT,SMCMAX ,TBOT,YY,ZZ1, ZBOT, QUARTZ, CSOIL, CAPR
3753       REAL, DIMENSION(1:NSOIL), INTENT(IN)   :: SMC,STC,ZSOIL
3754       REAL, DIMENSION(1:NSOIL), INTENT(OUT)  :: RHSTS
3755       REAL, DIMENSION(1:NSOIL), INTENT(OUT)  :: AI, BI,CI
3756       REAL                 :: DDZ, DDZ2, DENOM, DF1K, DTSDZ,DF1N,      &
3757                               DTSDZ2,HCPCT,QTOT,SSOIL,SICE,TAVG,TBK,     &
3758                               TBK1,TSNSR,TSURF
3759       REAL, PARAMETER      :: CAIR = 1004.0, CH2O = 4.2E6
3760
3761
3762 ! ----------------------------------------------------------------------
3763 ! INITIALIZE LOGICAL FOR SOIL LAYER TEMPERATURE AVERAGING.
3764 ! ----------------------------------------------------------------------
3765        ITAVG = .TRUE.
3766
3767 ! ----------------------------------------------------------------------
3768 ! TOP SOIL LAYER
3769 ! ----------------------------------------------------------------------
3770       HCPCT = SMC (1)* CH2O + (1.0- SMCMAX)* CSOIL + (SMCMAX - SMC (1))&
3771        * CAIR
3772       DDZ = 1.0 / ( -0.5 * ZSOIL (2) )
3773       AI (1) = 0.0
3774       CI (1) = (DF1 * DDZ) / (ZSOIL (1) * HCPCT)
3775
3776 ! ----------------------------------------------------------------------
3777 ! CALCULATE THE VERTICAL SOIL TEMP GRADIENT BTWN THE 1ST AND 2ND SOIL
3778 ! LAYERS.  THEN CALCULATE THE SUBSURFACE HEAT FLUX.
3779 ! ----------------------------------------------------------------------
3780       BI (1) = - CI (1) + DF1 / (0.5 * ZSOIL (1) * ZSOIL (1)* HCPCT *    &
3781        ZZ1)
3782       DTSDZ = (STC (1) - STC (2)) / (-0.5 * ZSOIL (2))
3783       SSOIL = DF1 * (STC (1) - YY) / (0.5 * ZSOIL (1) * ZZ1)
3784       DENOM = (ZSOIL (1) * HCPCT)
3785
3786 ! ----------------------------------------------------------------------
3787 ! NEXT CAPTURE THE VERTICAL DIFFERENCE OF THE HEAT FLUX AT TOP AND
3788 ! BOTTOM OF FIRST SOIL LAYER FOR USE IN HEAT FLUX CONSTRAINT
3789 ! ----------------------------------------------------------------------
3790       RHSTS (1) = (DF1 * DTSDZ - SSOIL) / DENOM
3791       QTOT = -1.0* RHSTS (1)* DENOM
3792       IF (ITAVG) THEN
3793          TSURF = (YY + (ZZ1-1) * STC (1)) / ZZ1
3794          CALL TBND (STC (1),STC (2),ZSOIL,ZBOT,1,NSOIL,TBK)
3795      ENDIF
3796       DDZ2 = 0.0
3797       DF1N = DF1
3798
3799 ! ----------------------------------------------------------------------
3800 ! LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABOVE PROCESS
3801 ! (EXCEPT SUBSFC OR "GROUND" HEAT FLUX NOT REPEATED IN LOWER LAYERS)
3802 ! ----------------------------------------------------------------------
3803       DO K = 2,NSOIL
3804 ! ----------------------------------------------------------------------
3805 ! THIS SECTION FOR LAYER 2 OR GREATER, BUT NOT LAST LAYER.
3806 ! ----------------------------------------------------------------------
3807          IF (K < NSOIL-1 ) THEN
3808          HCPCT = SMC (K)* CH2O + (1.0- SMCMAX)* CSOIL + (SMCMAX - SMC (  &
3809                 K))* CAIR
3810             CALL TDFCND  (DF1K, SMC(K), QUARTZ, SMCMAX)
3811             DENOM = 0.5 * ( ZSOIL (K -1) - ZSOIL (K +1) )
3812             DTSDZ2 = (STC (K) - STC (K +1) ) / DENOM
3813             DDZ2 = 2. / (ZSOIL (K -1) - ZSOIL (K +1))
3814
3815 ! ----------------------------------------------------------------------
3816 ! IF TEMPERATURE AVERAGING INVOKED (ITAVG=TRUE; ELSE SKIP):  CALCULATE
3817 ! TEMP AT BOTTOM OF LAYER.
3818 ! ----------------------------------------------------------------------
3819             CI (K) = - DF1K * DDZ2 / ( (ZSOIL (K -1) - ZSOIL (K)) *      &
3820        HCPCT)
3821             IF (ITAVG) THEN
3822                CALL TBND (STC (K),STC (K +1),ZSOIL,ZBOT,K,NSOIL,TBK1)
3823             END IF
3824
3825          ELSEIF (K == NSOIL-1) THEN
3826
3827          HCPCT = SMC (K)* CH2O + (1.0- SMCMAX)* CSOIL + (SMCMAX- SMC (  &
3828                 K))* CAIR
3829             CALL TDFCND  (DF1K, SMC(K), QUARTZ, SMCMAX)
3830             DENOM = 0.5 * ( ZSOIL (K -1) - ZSOIL (K +1) )
3831             DTSDZ2 = (STC (K) - STC (K +1) ) / DENOM
3832             DDZ2 = 2. / (ZSOIL (K -1) - ZSOIL (K +1))
3833 !-----------------------------------------------------------------------
3834 ! IF TEMPERATURE AVERAGING INVOKED (ITAVG=TRUE; ELSE SKIP):  CALCULATE
3835 ! TEMP AT BOTTOM OF LAST LAYER.
3836 ! ----------------------------------------------------------------------
3837             CI (K) = - DF1K * DDZ2 / ( (ZSOIL (K -1) - ZSOIL (K)) *      &
3838        HCPCT)
3839             IF (ITAVG) THEN
3840                CALL TBND (STC (K),TBOT,ZSOIL,ZBOT,K,NSOIL,TBK1)
3841             END IF
3842          ELSE
3843 ! ----------------------------------------------------------------------
3844 ! SPECIAL CASE OF BOTTOM LAYER (CONCRETE ROOF)
3845 ! ----------------------------------------------------------------------
3846             HCPCT = CAPR * 4.1868 * 1.E6
3847             DF1K  = 3.24
3848 ! ----------------------------------------------------------------------
3849 ! CALC THE VERTICAL TEMP GRADIENT THRU BOTTOM LAYER.
3850 ! ----------------------------------------------------------------------
3851             DENOM = .5 * (ZSOIL (K -1) + ZSOIL (K)) - ZBOT
3852             DTSDZ2 = (STC (K) - TBOT) / DENOM
3853 ! ----------------------------------------------------------------------
3854 ! IF TEMPERATURE AVERAGING INVOKED (ITAVG=TRUE; ELSE SKIP):  CALCULATE
3855 ! TEMP AT BOTTOM OF LAST LAYER.
3856 ! ----------------------------------------------------------------------
3857             CI (K) = 0.
3858             IF (ITAVG) THEN
3859                CALL TBND (STC (K),TBOT,ZSOIL,ZBOT,K,NSOIL,TBK1)
3860             END IF
3861 ! ----------------------------------------------------------------------
3862 ! THIS ENDS SPECIAL LOOP FOR BOTTOM LAYER.
3863          END IF
3864 ! ----------------------------------------------------------------------
3865 ! CALCULATE RHSTS FOR THIS LAYER AFTER CALC'NG A PARTIAL PRODUCT.
3866 ! ----------------------------------------------------------------------
3867          DENOM = ( ZSOIL (K) - ZSOIL (K -1) ) * HCPCT
3868          RHSTS (K) = ( DF1K * DTSDZ2- DF1N * DTSDZ ) / DENOM
3869          QTOT = -1.0* DENOM * RHSTS (K)
3870
3871 ! ----------------------------------------------------------------------
3872 ! CALC MATRIX COEFS, AI, AND BI FOR THIS LAYER.
3873 ! ----------------------------------------------------------------------
3874          AI (K) = - DF1N * DDZ / ( (ZSOIL (K -1) - ZSOIL (K)) * HCPCT)
3875
3876 ! ----------------------------------------------------------------------
3877 ! RESET VALUES OF DF1, DTSDZ, DDZ, AND TBK FOR LOOP TO NEXT SOIL LAYER.
3878 ! ----------------------------------------------------------------------
3879          BI (K) = - (AI (K) + CI (K))
3880          TBK = TBK1
3881          DF1N = DF1K
3882          DTSDZ = DTSDZ2
3883          DDZ = DDZ2
3884       END DO
3885 ! ----------------------------------------------------------------------
3886   END SUBROUTINE HRT
3887 ! ----------------------------------------------------------------------
3888
3889       SUBROUTINE HSTEP (STCOUT,STCIN,RHSTS,DT,NSOIL,AI,BI,CI)
3890 !      CALCULATE/UPDATE THE SOIL TEMPERATURE FIELD.
3891 ! ----------------------------------------------------------------------
3892       IMPLICIT NONE
3893       INTEGER, INTENT(IN)  :: NSOIL
3894       INTEGER              :: K
3895
3896       REAL, DIMENSION(1:NSOIL), INTENT(IN):: STCIN
3897       REAL, DIMENSION(1:NSOIL), INTENT(OUT):: STCOUT
3898       REAL, DIMENSION(1:NSOIL), INTENT(INOUT):: RHSTS
3899       REAL, DIMENSION(1:NSOIL), INTENT(INOUT):: AI,BI,CI
3900       REAL, DIMENSION(1:NSOIL) :: RHSTSin
3901       REAL, DIMENSION(1:NSOIL) :: CIin
3902       REAL                 :: DT
3903
3904 ! ----------------------------------------------------------------------
3905 ! CREATE FINITE DIFFERENCE VALUES FOR USE IN ROSR12 ROUTINE
3906 ! ----------------------------------------------------------------------
3907       DO K = 1,NSOIL
3908          RHSTS (K) = RHSTS (K) * DT
3909          AI (K) = AI (K) * DT
3910          BI (K) = 1. + BI (K) * DT
3911          CI (K) = CI (K) * DT
3912       END DO
3913 ! ----------------------------------------------------------------------
3914 ! COPY VALUES FOR INPUT VARIABLES BEFORE CALL TO ROSR12
3915 ! ----------------------------------------------------------------------
3916       DO K = 1,NSOIL
3917          RHSTSin (K) = RHSTS (K)
3918       END DO
3919       DO K = 1,NSOIL
3920          CIin (K) = CI (K)
3921       END DO
3922 ! ----------------------------------------------------------------------
3923 ! SOLVE THE TRI-DIAGONAL MATRIX EQUATION
3924 ! ----------------------------------------------------------------------
3925       CALL ROSR12 (CI,AI,BI,CIin,RHSTSin,RHSTS,NSOIL)
3926 ! ----------------------------------------------------------------------
3927 ! CALC/UPDATE THE SOIL TEMPS USING MATRIX SOLUTION
3928 ! ----------------------------------------------------------------------
3929       DO K = 1,NSOIL
3930          STCOUT (K) = STCIN (K) + CI (K)
3931       END DO
3932 ! ----------------------------------------------------------------------
3933   END SUBROUTINE HSTEP
3934 ! ----------------------------------------------------------------------
3935 ! ----------------------------------------------------------------------
3936
3937       SUBROUTINE TBND (TU,TB,ZSOIL,ZBOT,K,NSOIL,TBND1)
3938
3939 ! ----------------------------------------------------------------------
3940 ! SUBROUTINE TBND
3941 ! ----------------------------------------------------------------------
3942 ! CALCULATE TEMPERATURE ON THE BOUNDARY OF THE LAYER BY INTERPOLATION OF
3943 ! THE MIDDLE LAYER TEMPERATURES
3944 ! ----------------------------------------------------------------------
3945       IMPLICIT NONE
3946       INTEGER, INTENT(IN)       :: NSOIL
3947       INTEGER                   :: K
3948       REAL, INTENT(IN)          :: TB, TU, ZBOT
3949       REAL, INTENT(OUT)         :: TBND1
3950       REAL, DIMENSION(1:NSOIL), INTENT(IN)   :: ZSOIL
3951       REAL                      :: ZB, ZUP
3952
3953 ! ----------------------------------------------------------------------
3954 ! USE SURFACE TEMPERATURE ON THE TOP OF THE FIRST LAYER
3955 ! ----------------------------------------------------------------------
3956      IF (K == 1) THEN
3957          ZUP = 0.
3958       ELSE
3959          ZUP = ZSOIL (K -1)
3960       END IF
3961 ! ----------------------------------------------------------------------
3962 ! USE DEPTH OF THE CONSTANT BOTTOM TEMPERATURE WHEN INTERPOLATE
3963 ! TEMPERATURE INTO THE LAST LAYER BOUNDARY
3964 ! ----------------------------------------------------------------------
3965       IF (K ==  NSOIL) THEN
3966          ZB = 2.* ZBOT - ZSOIL (K)
3967       ELSE
3968          ZB = ZSOIL (K +1)
3969       END IF
3970 ! ----------------------------------------------------------------------
3971 ! LINEAR INTERPOLATION BETWEEN THE AVERAGE LAYER TEMPERATURES
3972 ! ----------------------------------------------------------------------
3973
3974       TBND1 = TU + (TB - TU)* (ZUP - ZSOIL (K))/ (ZUP - ZB)
3975 ! ----------------------------------------------------------------------
3976   END SUBROUTINE TBND
3977 ! ----------------------------------------------------------------------
3978       SUBROUTINE TDFCND (DF, SMC, QZ, SMCMAX)
3979 ! ----------------------------------------------------------------------
3980 ! CALCULATE THERMAL CONDUCTIVITY OF THE SOIL
3981 ! ----------------------------------------------------------------------
3982 ! PETERS-LIDARD APPROACH (PETERS-LIDARD et al., 1998)
3983 ! ----------------------------------------------------------------------
3984       IMPLICIT NONE
3985       REAL, INTENT(IN)          :: QZ,  SMC, SMCMAX
3986       REAL, INTENT(OUT)         :: DF
3987       REAL                      :: AKE, GAMMD, THKDRY, THKO,         &
3988                                    THKQTZ,THKSAT,THKS,THKW,SATRATIO
3989
3990 ! ----------------------------------------------------------------------
3991 ! IF THE SOIL HAS ANY MOISTURE CONTENT COMPUTE A PARTIAL SUM/PRODUCT
3992 ! OTHERWISE USE A CONSTANT VALUE WHICH WORKS WELL WITH MOST SOILS
3993 ! ----------------------------------------------------------------------
3994 !  THKW ......WATER THERMAL CONDUCTIVITY
3995 !  THKQTZ ....THERMAL CONDUCTIVITY FOR QUARTZ
3996 !  THKO ......THERMAL CONDUCTIVITY FOR OTHER SOIL COMPONENTS
3997 !  THKS ......THERMAL CONDUCTIVITY FOR THE SOLIDS COMBINED(QUARTZ+OTHER)
3998 !  SMCMAX ....POROSITY (= SMCMAX)
3999 !  QZ .........QUARTZ CONTENT (SOIL TYPE DEPENDENT)
4000 ! ----------------------------------------------------------------------
4001 ! USE AS IN PETERS-LIDARD, 1998 (MODIF. FROM JOHANSEN, 1975).
4002
4003 !                                  PABLO GRUNMANN, 08/17/98
4004 ! REFS.:
4005 !      FAROUKI, O.T.,1986: THERMAL PROPERTIES OF SOILS. SERIES ON ROCK
4006 !              AND SOIL MECHANICS, VOL. 11, TRANS TECH, 136 PP.
4007 !      JOHANSEN, O., 1975: THERMAL CONDUCTIVITY OF SOILS. PH.D. THESIS,
4008 !              UNIVERSITY OF TRONDHEIM,
4009 !      PETERS-LIDARD, C. D., ET AL., 1998: THE EFFECT OF SOIL THERMAL
4010 !              CONDUCTIVITY PARAMETERIZATION ON SURFACE ENERGY FLUXES
4011 !              AND TEMPERATURES. JOURNAL OF THE ATMOSPHERIC SCIENCES,
4012 !              VOL. 55, PP. 1209-1224.
4013 ! ----------------------------------------------------------------------
4014 ! NEEDS PARAMETERS
4015 ! POROSITY(SOIL TYPE):
4016 !      POROS = SMCMAX
4017 ! SATURATION RATIO:
4018 ! PARAMETERS  W/(M.K)
4019       SATRATIO = SMC / SMCMAX
4020 ! WATER CONDUCTIVITY:
4021       THKW = 0.57
4022 ! THERMAL CONDUCTIVITY OF "OTHER" SOIL COMPONENTS
4023 !      IF (QZ .LE. 0.2) THKO = 3.0
4024       THKO = 2.0
4025 ! QUARTZ' CONDUCTIVITY
4026       THKQTZ = 7.7
4027 ! SOLIDS' CONDUCTIVITY
4028       THKS = (THKQTZ ** QZ)* (THKO ** (1. - QZ))
4029
4030 ! SATURATED THERMAL CONDUCTIVITY
4031       THKSAT = THKS ** (1. - SMCMAX)* THKW ** (SMCMAX)
4032
4033 ! DRY DENSITY IN KG/M3
4034       GAMMD = (1. - SMCMAX)*2700.
4035
4036 ! DRY THERMAL CONDUCTIVITY IN W.M-1.K-1
4037       THKDRY = (0.135* GAMMD+ 64.7)/ (2700. - 0.947* GAMMD)
4038
4039 ! KERSTEN NUMBER (USING "FINE" FORMULA, VALID FOR SOILS CONTAINING AT
4040 ! LEAST 5% OF PARTICLES WITH DIAMETER LESS THAN 2.E-6 METERS.)
4041 ! (FOR "COARSE" FORMULA, SEE PETERS-LIDARD ET AL., 1998).
4042
4043          IF ( SATRATIO >  0.1 ) THEN
4044
4045             AKE = LOG10 (SATRATIO) + 1.0
4046
4047 ! USE K = KDRY
4048          ELSE
4049
4050             AKE = 0.0
4051          END IF
4052 !  THERMAL CONDUCTIVITY
4053
4054       DF = AKE * (THKSAT - THKDRY) + THKDRY
4055 ! ----------------------------------------------------------------------
4056   END SUBROUTINE TDFCND
4057 ! ----------------------------------------------------------------------
4058 !===========================================================================
4059 END MODULE module_sf_urban