PageRenderTime 73ms CodeModel.GetById 24ms RepoModel.GetById 0ms app.codeStats 1ms

/wrfv2_fire/phys/module_sf_urban.F

http://github.com/jbeezley/wrf-fire
FORTRAN Legacy | 2421 lines | 1481 code | 375 blank | 565 comment | 237 complexity | da74e6c6f59556f308ac988182e383d8 MD5 | raw file
Possible License(s): AGPL-1.0

Large files files are truncated, but you can click here to view the full file

    

  1. MODULE module_sf_urban
    2. 

  2. 

  3. !===============================================================================
    4. 

  4. !     Single-Layer Urban Canopy Model for WRF Noah-LSM
    5. 

  5. !     Original Version: 2002/11/06 by Hiroyuki Kusaka
    6. 

  6. !     Last Update:      2006/08/24 by Fei Chen and Mukul Tewari (NCAR/RAL)  
    7. 

  7. !===============================================================================
    8. 

  8. 

  9.    CHARACTER(LEN=4)                :: LU_DATA_TYPE
    10. 

  10. 

  11.    INTEGER                         :: ICATE
    12. 

  12. 

  13.    REAL, ALLOCATABLE, DIMENSION(:) :: ZR_TBL
    14. 

  14.    REAL, ALLOCATABLE, DIMENSION(:) :: Z0C_TBL
    15. 

  15.    REAL, ALLOCATABLE, DIMENSION(:) :: Z0HC_TBL
    16. 

  16.    REAL, ALLOCATABLE, DIMENSION(:) :: ZDC_TBL
    17. 

  17.    REAL, ALLOCATABLE, DIMENSION(:) :: SVF_TBL
    18. 

  18.    REAL, ALLOCATABLE, DIMENSION(:) :: R_TBL
    19. 

  19.    REAL, ALLOCATABLE, DIMENSION(:) :: RW_TBL
    20. 

  20.    REAL, ALLOCATABLE, DIMENSION(:) :: HGT_TBL
    21. 

  21.    REAL, ALLOCATABLE, DIMENSION(:) :: AH_TBL
    22. 

  22.    REAL, ALLOCATABLE, DIMENSION(:) :: BETR_TBL
    23. 

  23.    REAL, ALLOCATABLE, DIMENSION(:) :: BETB_TBL
    24. 

  24.    REAL, ALLOCATABLE, DIMENSION(:) :: BETG_TBL
    25. 

  25.    REAL, ALLOCATABLE, DIMENSION(:) :: FRC_URB_TBL
    26. 

  26. 

  27.    REAL, ALLOCATABLE, DIMENSION(:) :: COP_TBL
    28. 

  28.    REAL, ALLOCATABLE, DIMENSION(:) :: PWIN_TBL
    29. 

  29.    REAL, ALLOCATABLE, DIMENSION(:) :: BETA_TBL
    30. 

  30.    INTEGER, ALLOCATABLE, DIMENSION(:) :: SW_COND_TBL
    31. 

  31.    REAL, ALLOCATABLE, DIMENSION(:) :: TIME_ON_TBL
    32. 

  32.    REAL, ALLOCATABLE, DIMENSION(:) :: TIME_OFF_TBL
    33. 

  33.    REAL, ALLOCATABLE, DIMENSION(:) :: TARGTEMP_TBL
    34. 

  34.    REAL, ALLOCATABLE, DIMENSION(:) :: GAPTEMP_TBL
    35. 

  35.    REAL, ALLOCATABLE, DIMENSION(:) :: TARGHUM_TBL
    36. 

  36.    REAL, ALLOCATABLE, DIMENSION(:) :: GAPHUM_TBL
    37. 

  37.    REAL, ALLOCATABLE, DIMENSION(:) :: PERFLO_TBL
    38. 

  38.    REAL, ALLOCATABLE, DIMENSION(:) :: HSESF_TBL
    39. 

  39. 

  40.    REAL, ALLOCATABLE, DIMENSION(:) :: CAPR_TBL, CAPB_TBL, CAPG_TBL
    41. 

  41.    REAL, ALLOCATABLE, DIMENSION(:) :: AKSR_TBL, AKSB_TBL, AKSG_TBL
    42. 

  42.    REAL, ALLOCATABLE, DIMENSION(:) :: ALBR_TBL, ALBB_TBL, ALBG_TBL
    43. 

  43.    REAL, ALLOCATABLE, DIMENSION(:) :: EPSR_TBL, EPSB_TBL, EPSG_TBL
    44. 

  44.    REAL, ALLOCATABLE, DIMENSION(:) :: Z0R_TBL,  Z0B_TBL,  Z0G_TBL
    45. 

  45.    REAL, ALLOCATABLE, DIMENSION(:) :: SIGMA_ZED_TBL
    46. 

  46.    REAL, ALLOCATABLE, DIMENSION(:) :: Z0HB_TBL, Z0HG_TBL
    47. 

  47.    REAL, ALLOCATABLE, DIMENSION(:) :: TRLEND_TBL, TBLEND_TBL, TGLEND_TBL
    48. 

  48.    REAL, ALLOCATABLE, DIMENSION(:) :: AKANDA_URBAN_TBL
    49. 

  49. !for BEP
    50. 

  50. 

  51.    ! MAXDIRS :: The maximum number of street directions we're allowed to define
    52. 

  52.    INTEGER, PARAMETER :: MAXDIRS = 3
    53. 

  53.    ! MAXHGTS :: The maximum number of building height bins we're allowed to define
    54. 

  54.    INTEGER, PARAMETER :: MAXHGTS = 50
    55. 

  55. 

  56.    INTEGER, ALLOCATABLE, DIMENSION(:)   :: NUMDIR_TBL
    57. 

  57.    REAL,    ALLOCATABLE, DIMENSION(:,:) :: STREET_DIRECTION_TBL
    58. 

  58.    REAL,    ALLOCATABLE, DIMENSION(:,:) :: STREET_WIDTH_TBL
    59. 

  59.    REAL,    ALLOCATABLE, DIMENSION(:,:) :: BUILDING_WIDTH_TBL
    60. 

  60.    INTEGER, ALLOCATABLE, DIMENSION(:)   :: NUMHGT_TBL
    61. 

  61.    REAL,    ALLOCATABLE, DIMENSION(:,:) :: HEIGHT_BIN_TBL
    62. 

  62.    REAL,    ALLOCATABLE, DIMENSION(:,:) :: HPERCENT_BIN_TBL
    63. 

  63. !end BEP
    64. 

  64.    INTEGER                         :: BOUNDR_DATA,BOUNDB_DATA,BOUNDG_DATA
    65. 

  65.    INTEGER                         :: CH_SCHEME_DATA, TS_SCHEME_DATA
    66. 

  66.    INTEGER                         :: ahoption        ! Miao, 2007/01/17, cal. ah
    67. 

  67.    REAL, DIMENSION(1:24)           :: ahdiuprf        ! ah diurnal profile, tloc: 1-24
    68. 

  68.    REAL, DIMENSION(1:24)           :: hsequip_tbl
    69. 

  69. 

  70.    INTEGER                         :: allocate_status 
    71. 

  71. 

  72. !   INTEGER                         :: num_roof_layers
    73. 

  73. !   INTEGER                         :: num_wall_layers
    74. 

  74. !   INTEGER                         :: num_road_layers
    75. 

  75. 

  76.    CONTAINS
    77. 

  77. 

  78. !===============================================================================
    79. 

  79. !
    80. 

  80. ! Author:
    81. 

  81. !      Hiroyuki KUSAKA, PhD
    82. 

  82. !      University of Tsukuba, JAPAN
    83. 

  83. !      (CRIEPI, NCAR/MMM visiting scientist, 2002-2004)
    84. 

  84. !      kusaka@ccs.tsukuba.ac.jp
    85. 

  85. !
    86. 

  86. ! Co-Researchers:
    87. 

  87. !     Fei CHEN, PhD
    88. 

  88. !      NCAR/RAP feichen@ucar.edu
    89. 

  89. !     Mukul TEWARI, PhD
    90. 

  90. !      NCAR/RAP mukul@ucar.edu
    91. 

  91. !
    92. 

  92. ! Purpose:
    93. 

  93. !     Calculate surface temeprature, fluxes, canopy air temperature, and canopy wind
    94. 

  94. !
    95. 

  95. ! Subroutines:
    96. 

  96. !     module_sf_urban
    97. 

  97. !       |- urban
    98. 

  98. !            |- read_param
    99. 

  99. !            |- mos or jurges
    100. 

  100. !            |- multi_layer or force_restore
    101. 

  101. !       |- urban_param_init <-- URBPARM.TBL
    102. 

  102. !       |- urban_var_init
    103. 

  103. !
    104. 

  104. ! Input Data from WRF [MKS unit]:
    105. 

  105. !
    106. 

  106. !     UTYPE  [-]     : Urban type. 1=Commercial/Industrial; 2=High-intensity residential; 
    107. 

  107. !                    : 3=low-intensity residential 
    108. 

  108. !     TA     [K]     : Potential temperature at 1st wrf level (absolute temp)
    109. 

  109. !     QA     [kg/kg] : Mixing ratio at 1st atmospheric level
    110. 

  110. !     UA     [m/s]   : Wind speed at 1st atmospheric level
    111. 

  111. !     SSG    [W/m/m] : Short wave downward radiation at a flat surface
    112. 

  112. !                      Note this is the total of direct and diffusive solar
    113. 

  113. !                      downward radiation. If without two components, the
    114. 

  114. !                      single solar downward can be used instead.
    115. 

  115. !                      SSG = SSGD + SSGQ
    116. 

  116. !     LSOLAR [-]     : Indicating the input type of solar downward radiation
    117. 

  117. !                      True: both direct and diffusive solar radiation 
    118. 

  118. !                      are available
    119. 

  119. !                      False: only total downward ridiation is available.
    120. 

  120. !     SSGD   [W/m/m] : Direct solar radiation at a flat surface
    121. 

  121. !                      if SSGD is not available, one can assume a ratio SRATIO
    122. 

  122. !                      (e.g., 0.7), so that SSGD = SRATIO*SSG 
    123. 

  123. !     SSGQ   [W/m/m] : Diffuse solar radiation at a flat surface
    124. 

  124. !                      If SSGQ is not available, SSGQ = SSG - SSGD
    125. 

  125. !     LLG    [W/m/m] : Long wave downward radiation at a flat surface 
    126. 

  126. !     RAIN   [mm/h]  : Precipitation
    127. 

  127. !     RHOO   [kg/m/m/m] : Air density
    128. 

  128. !     ZA     [m]     : First atmospheric level
    129. 

  129. !                      as a lowest boundary condition
    130. 

  130. !     DECLIN [rad]   : solar declination
    131. 

  131. !     COSZ           : = sin(fai)*sin(del)+cos(fai)*cos(del)*cos(omg)
    132. 

  132. !     OMG    [rad]   : solar hour angle
    133. 

  133. !     XLAT   [deg]   : latitude
    134. 

  134. !     DELT   [sec]   : Time step
    135. 

  135. !     ZNT    [m]     : Roughnes length
    136. 

  136. !
    137. 

  137. ! Output Data to WRF [MKS unit]:
    138. 

  138. !
    139. 

  139. !     TS  [K]            : Surface potential temperature (absolute temp)
    140. 

  140. !     QS  [-]            : Surface humidity
    141. 

  141. !
    142. 

  142. !     SH  [W/m/m/]       : Sensible heat flux, = FLXTH*RHOO*CPP
    143. 

  143. !     LH  [W/m/m]        : Latent heat flux, = FLXHUM*RHOO*ELL
    144. 

  144. !     LH_INEMATIC [kg/m/m/sec]: Moisture Kinematic flux, = FLXHUM*RHOO
    145. 

  145. !     SW  [W/m/m]        : Upward shortwave radiation flux, 
    146. 

  146. !                          = SSG-SNET*697.7*60. (697.7*60.=100.*100.*4.186)
    147. 

  147. !     ALB [-]            : Time-varying albedo
    148. 

  148. !     LW  [W/m/m]        : Upward longwave radiation flux, 
    149. 

  149. !                          = LNET*697.7*60.-LLG
    150. 

  150. !     G   [W/m/m]        : Heat Flux into the Ground
    151. 

  151. !     RN  [W/m/m]        : Net radiation
    152. 

  152. !
    153. 

  153. !     PSIM  [-]          : Diagnostic similarity stability function for momentum
    154. 

  154. !     PSIH  [-]          : Diagnostic similarity stability function for heat
    155. 

  155. !
    156. 

  156. !     TC  [K]            : Diagnostic canopy air temperature
    157. 

  157. !     QC  [-]            : Diagnostic canopy humidity
    158. 

  158. !
    159. 

  159. !     TH2 [K]            : Diagnostic potential temperature at 2 m
    160. 

  160. !     Q2  [-]            : Diagnostic humidity at 2 m
    161. 

  161. !     U10 [m/s]          : Diagnostic u wind component at 10 m
    162. 

  162. !     V10 [m/s]          : Diagnostic v wind component at 10 m
    163. 

  163. !
    164. 

  164. !     CHS, CHS2 [m/s]    : CH*U at ZA, CH*U at 2 m (not used)
    165. 

  165. !
    166. 

  166. ! Important parameters:
    167. 

  167. !
    168. 

  168. ! Morphology of the urban canyon:
    169. 

  169. ! These parameters assigned in the URBPARM.TBL
    170. 

  170. !
    171. 

  171. !    ZR  [m]             : roof level (building height)
    172. 

  172. !    SIGMA_ZED [m]       : Standard Deviation of roof height
    173. 

  173. !    ROOF_WIDTH [m]      : roof (i.e., building) width
    174. 

  174. !    ROAD_WIDTH [m]      : road width
    175. 

  175. !
    176. 

  176. ! Parameters derived from the morphological terms above.
    177. 

  177. ! These parameters are computed in the code.
    178. 

  178. !
    179. 

  179. !    HGT [-]             : normalized building height
    180. 

  180. !    SVF [-]             : sky view factor
    181. 

  181. !    R   [-]             : Normalized roof width (a.k.a. "building coverage ratio")
    182. 

  182. !    RW  [-]             : = 1 - R
    183. 

  183. !    Z0C [m]             : Roughness length above canyon for momentum (1/10 of ZR)
    184. 

  184. !    Z0HC [m]            : Roughness length above canyon for heat (1/10 of Z0C)
    185. 

  185. !    ZDC [m]             : Zero plane displacement height (1/5 of ZR)
    186. 

  186. !
    187. 

  187. ! Following parameter are assigned in run/URBPARM.TBL
    188. 

  188. !
    189. 

  189. !  AH  [ W m{-2} ]        : anthropogenic heat  ( W m{-2} in the table, converted internally to cal cm{-2} )
    190. 

  190. !  CAPR[ J m{-3} K{-1} ]  : heat capacity of roof ( units converted in code to [ cal cm{-3} deg{-1} ] )
    191. 

  191. !  CAPB[ J m{-3} K{-1} ]  : heat capacity of building wall ( units converted in code to [ cal cm{-3} deg{-1} ] )
    192. 

  192. !  CAPG[ J m{-3} K{-1} ]  : heat capacity of road ( units converted in code to [ cal cm{-3} deg{-1} ] )
    193. 

  193. !  AKSR [ J m{-1} s{-1} K{-1} ] : thermal conductivity of roof ( units converted in code to [ cal cm{-1} s{-1} deg{-1} ] )
    194. 

  194. !  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} ] )
    195. 

  195. !  AKSG [ J m{-1} s{-1} K{-1} ] : thermal conductivity of road ( units converted in code to [ cal cm{-1} s{-1} deg{-1} ] )
    196. 

  196. !  ALBR [-]               : surface albedo of roof
    197. 

  197. !  ALBB [-]               : surface albedo of building wall
    198. 

  198. !  ALBG [-]               : surface albedo of road
    199. 

  199. !  EPSR [-]               : surface emissivity of roof
    200. 

  200. !  EPSB [-]               : surface emissivity of building wall
    201. 

  201. !  EPSG [-]               : surface emissivity of road
    202. 

  202. !  Z0B [m]                : roughness length for momentum of building wall (only for CH_SCHEME = 1)
    203. 

  203. !  Z0G [m]                : roughness length for momentum of road (only for CH_SCHEME = 1)
    204. 

  204. !  Z0HB [m]               : roughness length for heat of building wall (only for CH_SCHEME = 1)
    205. 

  205. !  Z0HG [m]               : roughness length for heat of road
    206. 

  206. !  num_roof_layers        : number of layers within roof
    207. 

  207. !  num_wall_layers        : number of layers within building walls
    208. 

  208. !  num_road_layers        : number of layers within below road surface
    209. 

  209. !   NOTE: for now, these layers are defined as same as the number of soil layers in namelist.input
    210. 

  210. !  DZR [cm]               : thickness of each roof layer
    211. 

  211. !  DZB [cm]               : thickness of each building wall layer
    212. 

  212. !  DZG [cm]               : thickness of each ground layer
    213. 

  213. !  BOUNDR [integer 1 or 2] : Boundary Condition for Roof Layer Temp [1: Zero-Flux, 2: T = Constant]
    214. 

  214. !  BOUNDB [integer 1 or 2] : Boundary Condition for Building Wall Layer Temp [1: Zero-Flux, 2: T = Constant]
    215. 

  215. !  BOUNDG [integer 1 or 2] : Boundary Condition for Road Layer Temp [1: Zero-Flux, 2: T = Constant]
    216. 

  216. !  TRLEND [K]             : lower boundary condition of roof temperature
    217. 

  217. !  TBLEND [K]             : lower boundary condition of building temperature
    218. 

  218. !  TGLEND [K]             : lower boundary condition of ground temperature
    219. 

  219. !  CH_SCHEME [integer 1 or 2] : Sfc exchange scheme used for building wall and road
    220. 

  220. !                         [1: M-O Similarity Theory,  2: Empirical Form (recommend)]
    221. 

  221. !  TS_SCHEME [integer 1 or 2] : Scheme for computing surface temperature (for roof, wall, and road)
    222. 

  222. !                         [1: 4-layer model,  2: Force-Restore method]
    223. 

  223. !
    224. 

  224. !for BEP
    225. 

  225. !  numdir [ - ]             : Number of street directions defined for a particular urban category
    226. 

  226. !  street_direction [ deg ] : Direction of streets for a particular urban category and a particular street direction
    227. 

  227. !  street_width [ m ]       : Width of street for a particular urban category and a particular street direction
    228. 

  228. !  building_width [ m ]     : Width of buildings for a particular urban category and a particular street direction
    229. 

  229. !  numhgt [ - ]             : Number of building height levels defined for a particular urban category
    230. 

  230. !  height_bin [ m ]         : Building height bins defined for a particular urban category.
    231. 

  231. !  hpercent_bin [ % ]       : Percentage of a particular urban category populated by buildings of particular height bins
    232. 

  232. !end BEP
    233. 

  233. ! Moved from URBPARM.TBL
    234. 

  234. !
    235. 

  235. !  BETR [-]            : minimum moisture availability of roof
    236. 

  236. !  BETB [-]            : minimum moisture availability of building wall
    237. 

  237. !  BETG [-]            : minimum moisture availability of road
    238. 

  238. !  Z0R [m]                : roughness length for momentum of roof
    239. 

  239. !  Z0HB [m]               : roughness length for heat of building wall (only for CH_SCHEME = 1)
    240. 

  240. !  Z0HG [m]               : roughness length for heat of road
    241. 

  241. !  num_roof_layers        : number of layers within roof
    242. 

  242. !  num_wall_layers        : number of layers within building walls
    243. 

  243. !  num_road_layers        : number of layers within below road surface
    244. 

  244. !   NOTE: for now, these layers are defined as same as the number of soil layers in namelist.input
    245. 

  245. !
    246. 

  246. ! References:
    247. 

  247. !     Kusaka and Kimura (2004) J.Appl.Meteor., vol.43, p1899-1910
    248. 

  248. !     Kusaka and Kimura (2004) J.Meteor.Soc.Japan, vol.82, p45-65
    249. 

  249. !     Kusaka et al. (2001) Bound.-Layer Meteor., vol.101, p329-358
    250. 

  250. !
    251. 

  251. ! History:
    252. 

  252. !     2006/06          modified by H. Kusaka (Univ. Tsukuba), M. Tewari 
    253. 

  253. !     2005/10/26,      modified by Fei Chen, Mukul Tewari
    254. 

  254. !     2003/07/21 WRF , modified  by H. Kusaka of CRIEPI (NCAR/MMM)
    255. 

  255. !     2001/08/26 PhD , modified  by H. Kusaka of CRIEPI (Univ.Tsukuba)
    256. 

  256. !     1999/08/25 LCM , developed by H. Kusaka of CRIEPI (Univ.Tsukuba)
    257. 

  257. !
    258. 

  258. !===============================================================================
    259. 

  259. !
    260. 

  260. !  subroutine urban:
    261. 

  261. !
    262. 

  262. !===============================================================================
    263. 

  263. 

  264.    SUBROUTINE urban(LSOLAR,                                           & ! L
    265. 

  265.                     num_roof_layers,num_wall_layers,num_road_layers,  & ! I
    266. 

  266.                     DZR,DZB,DZG,                                      & ! I
    267. 

  267.                     UTYPE,TA,QA,UA,U1,V1,SSG,SSGD,SSGQ,LLG,RAIN,RHOO, & ! I
    268. 

  268.                     ZA,DECLIN,COSZ,OMG,XLAT,DELT,ZNT,                 & ! I
    269. 

  269.                     CHS, CHS2,                                        & ! I
    270. 

  270.                     TR, TB, TG, TC, QC, UC,                           & ! H
    271. 

  271.                     TRL,TBL,TGL,                                      & ! H
    272. 

  272.                     XXXR, XXXB, XXXG, XXXC,                           & ! H
    273. 

  273.                     TS,QS,SH,LH,LH_KINEMATIC,                         & ! O
    274. 

  274.                     SW,ALB,LW,G,RN,PSIM,PSIH,                         & ! O
    275. 

  275.                     GZ1OZ0,                                           & ! O
    276. 

  276.                     CMR_URB,CHR_URB,CMC_URB,CHC_URB,                  & ! I/O
    277. 

  277.                     U10,V10,TH2,Q2,UST                                & ! O
    278. 

  278.                     )
    279. 

  279. 

  280.    IMPLICIT NONE
    281. 

  281. 

  282.    REAL, PARAMETER    :: CP=0.24          ! heat capacity of dry air  [cgs unit]
    283. 

  283.    REAL, PARAMETER    :: EL=583.          ! latent heat of vaporation [cgs unit]
    284. 

  284.    REAL, PARAMETER    :: SIG=8.17E-11     ! stefun bolzman constant   [cgs unit]
    285. 

  285.    REAL, PARAMETER    :: SIG_SI=5.67E-8   !                           [MKS unit]
    286. 

  286.    REAL, PARAMETER    :: AK=0.4           ! kalman const.                [-]
    287. 

  287.    REAL, PARAMETER    :: PI=3.14159       ! pi                           [-]
    288. 

  288.    REAL, PARAMETER    :: TETENA=7.5       ! const. of Tetens Equation    [-]
    289. 

  289.    REAL, PARAMETER    :: TETENB=237.3     ! const. of Tetens Equation    [-]
    290. 

  290.    REAL, PARAMETER    :: SRATIO=0.75      ! ratio between direct/total solar [-]
    291. 

  291. 

  292.    REAL, PARAMETER    :: CPP=1004.5       ! heat capacity of dry air    [J/K/kg]
    293. 

  293.    REAL, PARAMETER    :: ELL=2.442E+06    ! latent heat of vaporization [J/kg]
    294. 

  294.    REAL, PARAMETER    :: XKA=2.4E-5
    295. 

  295. 

  296. !-------------------------------------------------------------------------------
    297. 

  297. ! C: configuration variables
    298. 

  298. !-------------------------------------------------------------------------------
    299. 

  299. 

  300.    LOGICAL, INTENT(IN) :: LSOLAR  ! logical    [true=both, false=SSG only]
    301. 

  301. 

  302. !  The following variables are also model configuration variables, but are 
    303. 

  303. !  defined in the URBAN.TBL and in the contains statement in the top of 
    304. 

  304. !  the module_urban_init, so we should not declare them here.
    305. 

  305. 

  306.   INTEGER, INTENT(IN) :: num_roof_layers
    307. 

  307.   INTEGER, INTENT(IN) :: num_wall_layers
    308. 

  308.   INTEGER, INTENT(IN) :: num_road_layers
    309. 

  309. 

  310. 

  311.   REAL, INTENT(IN), DIMENSION(1:num_roof_layers) :: DZR ! grid interval of roof layers [cm]
    312. 

  312.   REAL, INTENT(IN), DIMENSION(1:num_wall_layers) :: DZB ! grid interval of wall layers [cm]
    313. 

  313.   REAL, INTENT(IN), DIMENSION(1:num_road_layers) :: DZG ! grid interval of road layers [cm]
    314. 

  314. 

  315. !-------------------------------------------------------------------------------
    316. 

  316. ! I: input variables from LSM to Urban
    317. 

  317. !-------------------------------------------------------------------------------
    318. 

  318. 

  319.    INTEGER, INTENT(IN) :: UTYPE ! urban type [1=Commercial/Industrial, 2=High-intensity residential, 
    320. 

  320.                                 ! 3=low-intensity residential]
    321. 

  321. 

  322.    REAL, INTENT(IN)    :: TA   ! potential temp at 1st atmospheric level [K]
    323. 

  323.    REAL, INTENT(IN)    :: QA   ! mixing ratio at 1st atmospheric level  [kg/kg]
    324. 

  324.    REAL, INTENT(IN)    :: UA   ! wind speed at 1st atmospheric level    [m/s]
    325. 

  325.    REAL, INTENT(IN)    :: U1   ! u at 1st atmospheric level             [m/s]
    326. 

  326.    REAL, INTENT(IN)    :: V1   ! v at 1st atmospheric level             [m/s]
    327. 

  327.    REAL, INTENT(IN)    :: SSG  ! downward total short wave radiation    [W/m/m]
    328. 

  328.    REAL, INTENT(IN)    :: LLG  ! downward long wave radiation           [W/m/m]
    329. 

  329.    REAL, INTENT(IN)    :: RAIN ! precipitation                          [mm/h]
    330. 

  330.    REAL, INTENT(IN)    :: RHOO ! air density                            [kg/m^3]
    331. 

  331.    REAL, INTENT(IN)    :: ZA   ! first atmospheric level                [m]
    332. 

  332.    REAL, INTENT(IN)    :: DECLIN ! solar declination                    [rad]
    333. 

  333.    REAL, INTENT(IN)    :: COSZ ! sin(fai)*sin(del)+cos(fai)*cos(del)*cos(omg)
    334. 

  334.    REAL, INTENT(IN)    :: OMG  ! solar hour angle                       [rad]
    335. 

  335. 

  336.    REAL, INTENT(IN)    :: XLAT ! latitude                               [deg]
    337. 

  337.    REAL, INTENT(IN)    :: DELT ! time step                              [s]
    338. 

  338.    REAL, INTENT(IN)    :: ZNT  ! roughness length                       [m]
    339. 

  339.    REAL, INTENT(IN)    :: CHS,CHS2 ! CH*U at za and 2 m             [m/s]
    340. 

  340. 

  341.    REAL, INTENT(INOUT) :: SSGD ! downward direct short wave radiation   [W/m/m]
    342. 

  342.    REAL, INTENT(INOUT) :: SSGQ ! downward diffuse short wave radiation  [W/m/m]
    343. 

  343.    REAL, INTENT(INOUT) :: CMR_URB
    344. 

  344.    REAL, INTENT(INOUT) :: CHR_URB
    345. 

  345.    REAL, INTENT(INOUT) :: CMC_URB
    346. 

  346.    REAL, INTENT(INOUT) :: CHC_URB
    347. 

  347. 

  348. !-------------------------------------------------------------------------------
    349. 

  349. ! O: output variables from Urban to LSM 
    350. 

  350. !-------------------------------------------------------------------------------
    351. 

  351. 

  352.    REAL, INTENT(OUT) :: TS     ! surface potential temperature    [K]
    353. 

  353.    REAL, INTENT(OUT) :: QS     ! surface humidity                 [K]
    354. 

  354.    REAL, INTENT(OUT) :: SH     ! sensible heat flux               [W/m/m]
    355. 

  355.    REAL, INTENT(OUT) :: LH     ! latent heat flux                 [W/m/m]
    356. 

  356.    REAL, INTENT(OUT) :: LH_KINEMATIC ! latent heat, kinetic     [kg/m/m/s]
    357. 

  357.    REAL, INTENT(OUT) :: SW     ! upward short wave radiation flux [W/m/m]
    358. 

  358.    REAL, INTENT(OUT) :: ALB    ! time-varying albedo            [fraction]
    359. 

  359.    REAL, INTENT(OUT) :: LW     ! upward long wave radiation flux  [W/m/m]
    360. 

  360.    REAL, INTENT(OUT) :: G      ! heat flux into the ground        [W/m/m]
    361. 

  361.    REAL, INTENT(OUT) :: RN     ! net radition                     [W/m/m]
    362. 

  362.    REAL, INTENT(OUT) :: PSIM   ! similality stability shear function for momentum
    363. 

  363.    REAL, INTENT(OUT) :: PSIH   ! similality stability shear function for heat
    364. 

  364.    REAL, INTENT(OUT) :: GZ1OZ0   
    365. 

  365.    REAL, INTENT(OUT) :: U10    ! u at 10m                         [m/s]
    366. 

  366.    REAL, INTENT(OUT) :: V10    ! u at 10m                         [m/s]
    367. 

  367.    REAL, INTENT(OUT) :: TH2    ! potential temperature at 2 m     [K]
    368. 

  368.    REAL, INTENT(OUT) :: Q2     ! humidity at 2 m                  [-]
    369. 

  369. !m   REAL, INTENT(OUT) :: CHS,CHS2 ! CH*U at za and 2 m             [m/s]
    370. 

  370.    REAL, INTENT(OUT) :: UST    ! friction velocity                [m/s]
    371. 

  371. 

  372. 

  373. !-------------------------------------------------------------------------------
    374. 

  374. ! H: Historical (state) variables of Urban : LSM <--> Urban
    375. 

  375. !-------------------------------------------------------------------------------
    376. 

  376. 

  377. ! TR: roof temperature              [K];      TRP: at previous time step [K]
    378. 

  378. ! TB: building wall temperature     [K];      TBP: at previous time step [K]
    379. 

  379. ! TG: road temperature              [K];      TGP: at previous time step [K]
    380. 

  380. ! TC: urban-canopy air temperature  [K];      TCP: at previous time step [K]
    381. 

  381. !                                                  (absolute temperature)
    382. 

  382. ! QC: urban-canopy air mixing ratio [kg/kg];  QCP: at previous time step [kg/kg]
    383. 

  383. !
    384. 

  384. ! XXXR: Monin-Obkhov length for roof          [dimensionless]
    385. 

  385. ! XXXB: Monin-Obkhov length for building wall [dimensionless]
    386. 

  386. ! XXXG: Monin-Obkhov length for road          [dimensionless]
    387. 

  387. ! XXXC: Monin-Obkhov length for urban-canopy  [dimensionless]
    388. 

  388. !
    389. 

  389. ! TRL, TBL, TGL: layer temperature [K] (absolute temperature)
    390. 

  390. 

  391.    REAL, INTENT(INOUT):: TR, TB, TG, TC, QC, UC
    392. 

  392.    REAL, INTENT(INOUT):: XXXR, XXXB, XXXG, XXXC
    393. 

  393. 

  394.    REAL, DIMENSION(1:num_roof_layers), INTENT(INOUT) :: TRL
    395. 

  395.    REAL, DIMENSION(1:num_wall_layers), INTENT(INOUT) :: TBL
    396. 

  396.    REAL, DIMENSION(1:num_road_layers), INTENT(INOUT) :: TGL
    397. 

  397. 

  398. !-------------------------------------------------------------------------------
    399. 

  399. ! L:  Local variables from read_param
    400. 

  400. !-------------------------------------------------------------------------------
    401. 

  401. 

  402.    REAL :: ZR, Z0C, Z0HC, ZDC, SVF, R, RW, HGT, AH
    403. 

  403.    REAL :: SIGMA_ZED
    404. 

  404.    REAL :: CAPR, CAPB, CAPG, AKSR, AKSB, AKSG, ALBR, ALBB, ALBG 
    405. 

  405.    REAL :: EPSR, EPSB, EPSG, Z0R,  Z0B,  Z0G,  Z0HB, Z0HG
    406. 

  406.    REAL :: TRLEND,TBLEND,TGLEND
    407. 

  407.    REAL :: T1VR, T1VC,TH2V
    408. 

  408.    REAL :: RLMO_URB
    409. 

  409.    REAL :: AKANDA_URBAN
    410. 

  410.    
    411. 

  411.    REAL :: TH2X                                                !m
    412. 

  412.    
    413. 

  413.    INTEGER :: BOUNDR, BOUNDB, BOUNDG
    414. 

  414.    INTEGER :: CH_SCHEME, TS_SCHEME
    415. 

  415. 

  416.    LOGICAL :: SHADOW  ! [true=consider svf and shadow effects, false=consider svf effect only]
    417. 

  417. 

  418. !for BEP
    419. 

  419.    INTEGER                        :: NUMDIR
    420. 

  420.    REAL,    DIMENSION ( MAXDIRS ) :: STREET_DIRECTION
    421. 

  421.    REAL,    DIMENSION ( MAXDIRS ) :: STREET_WIDTH
    422. 

  422.    REAL,    DIMENSION ( MAXDIRS ) :: BUILDING_WIDTH
    423. 

  423.    INTEGER                        :: NUMHGT
    424. 

  424.    REAL,    DIMENSION ( MAXHGTS ) :: HEIGHT_BIN
    425. 

  425.    REAL,    DIMENSION ( MAXHGTS ) :: HPERCENT_BIN
    426. 

  426. !end BEP
    427. 

  427. !-------------------------------------------------------------------------------
    428. 

  428. ! L: Local variables
    429. 

  429. !-------------------------------------------------------------------------------
    430. 

  430. 

  431.    REAL :: BETR, BETB, BETG
    432. 

  432.    REAL :: SX, SD, SQ, RX
    433. 

  433.    REAL :: UR, ZC, XLB, BB
    434. 

  434.    REAL :: Z, RIBB, RIBG, RIBC, BHR, BHB, BHG, BHC
    435. 

  435.    REAL :: TSC, LNET, SNET, FLXUV, THG, FLXTH, FLXHUM, FLXG
    436. 

  436.    REAL :: W, VFGS, VFGW, VFWG, VFWS, VFWW
    437. 

  437.    REAL :: HOUI1, HOUI2, HOUI3, HOUI4, HOUI5, HOUI6, HOUI7, HOUI8
    438. 

  438.    REAL :: SLX, SLX1, SLX2, SLX3, SLX4, SLX5, SLX6, SLX7, SLX8
    439. 

  439.    REAL :: FLXTHR, FLXTHB, FLXTHG, FLXHUMR, FLXHUMB, FLXHUMG
    440. 

  440.    REAL :: SR, SB, SG, RR, RB, RG
    441. 

  441.    REAL :: SR1, SR2, SB1, SB2, SG1, SG2, RR1, RR2, RB1, RB2, RG1, RG2
    442. 

  442.    REAL :: HR, HB, HG, ELER, ELEB, ELEG, G0R, G0B, G0G
    443. 

  443.    REAL :: ALPHAC, ALPHAR, ALPHAB, ALPHAG
    444. 

  444.    REAL :: CHC, CHR, CHB, CHG, CDC, CDR, CDB, CDG
    445. 

  445.    REAL :: C1R, C1B, C1G, TE, TC1, TC2, QC1, QC2, QS0R, QS0B, QS0G,RHO,ES
    446. 

  446. 

  447.    REAL :: DESDT
    448. 

  448.    REAL :: F
    449. 

  449.    REAL :: DQS0RDTR
    450. 

  450.    REAL :: DRRDTR, DHRDTR, DELERDTR, DG0RDTR
    451. 

  451.    REAL :: DTR, DFDT
    452. 

  452.    REAL :: FX, FY, GF, GX, GY
    453. 

  453.    REAL :: DTCDTB, DTCDTG
    454. 

  454.    REAL :: DQCDTB, DQCDTG
    455. 

  455.    REAL :: DRBDTB1,  DRBDTG1,  DRBDTB2,  DRBDTG2
    456. 

  456.    REAL :: DRGDTB1,  DRGDTG1,  DRGDTB2,  DRGDTG2
    457. 

  457.    REAL :: DRBDTB,   DRBDTG,   DRGDTB,   DRGDTG
    458. 

  458.    REAL :: DHBDTB,   DHBDTG,   DHGDTB,   DHGDTG
    459. 

  459.    REAL :: DELEBDTB, DELEBDTG, DELEGDTG, DELEGDTB
    460. 

  460.    REAL :: DG0BDTB,  DG0BDTG,  DG0GDTG,  DG0GDTB
    461. 

  461.    REAL :: DQS0BDTB, DQS0GDTG
    462. 

  462.    REAL :: DTB, DTG, DTC
    463. 

  463. 

  464.    REAL :: THEATAZ    ! Solar Zenith Angle [rad]
    465. 

  465.    REAL :: THEATAS    ! = PI/2. - THETAZ
    466. 

  466.    REAL :: FAI        ! Latitude [rad]
    467. 

  467.    REAL :: CNT,SNT
    468. 

  468.    REAL :: PS         ! Surface Pressure [hPa]
    469. 

  469.    REAL :: TAV        ! Vertial Temperature [K] 
    470. 

  470. 

  471.    REAL :: XXX, X, Z0, Z0H, CD, CH
    472. 

  472.    REAL :: XXX2, PSIM2, PSIH2, XXX10, PSIM10, PSIH10
    473. 

  473.    REAL :: PSIX, PSIT, PSIX2, PSIT2, PSIX10, PSIT10
    474. 

  474. 

  475.    REAL :: TRP, TBP, TGP, TCP, QCP, TST, QST
    476. 

  476. 

  477.    INTEGER :: iteration, K 
    478. 

  478.    INTEGER :: tloc
    479. 

  479. 

  480. !-------------------------------------------------------------------------------
    481. 

  481. ! Set parameters
    482. 

  482. !-------------------------------------------------------------------------------
    483. 

  483. 

  484. ! Miao, 2007/01/17, cal. ah
    485. 

  485.    if(ahoption==1) then
    486. 

  486.      tloc=mod(int(OMG/PI*180./15.+12.+0.5 ),24)
    487. 

  487.      if(tloc==0) tloc=24
    488. 

  488.    endif
    489. 

  489. 

  490.    CALL read_param(UTYPE,ZR,SIGMA_ZED,Z0C,Z0HC,ZDC,SVF,R,RW,HGT,  &
    491. 

  491.                    AH,CAPR,CAPB,CAPG,AKSR,AKSB,AKSG,ALBR,ALBB,    &
    492. 

  492.                    ALBG,EPSR,EPSB,EPSG,Z0R,Z0B,Z0G,Z0HB,Z0HG,     &
    493. 

  493.                    BETR,BETB,BETG,TRLEND,TBLEND,TGLEND,           &
    494. 

  494. !for BEP
    495. 

  495.                    NUMDIR, STREET_DIRECTION, STREET_WIDTH,        & 
    496. 

  496.                    BUILDING_WIDTH, NUMHGT, HEIGHT_BIN,            & 
    497. 

  497.                    HPERCENT_BIN,                                  & 
    498. 

  498. !end BEP
    499. 

  499.                    BOUNDR,BOUNDB,BOUNDG,CH_SCHEME,TS_SCHEME,      &
    500. 

  500.                    AKANDA_URBAN)
    501. 

  501. 

  502. ! Miao, 2007/01/17, cal. ah
    503. 

  503.    if(ahoption==1) AH=AH*ahdiuprf(tloc)
    504. 

  504. 

  505.    IF( ZDC+Z0C+2. >= ZA) THEN
    506. 

  506.     CALL wrf_error_fatal ("ZDC + Z0C + 2m is larger than the 1st WRF level "// &
    507. 

  507.                            "Stop in subroutine urban - change ZDC and Z0C" ) 
    508. 

  508.    END IF
    509. 

  509. 

  510.    IF(.NOT.LSOLAR) THEN
    511. 

  511.      SSGD = SRATIO*SSG
    512. 

  512.      SSGQ = SSG - SSGD
    513. 

  513.    ENDIF
    514. 

  514.    SSGD = SRATIO*SSG   ! No radiation scheme has SSGD and SSGQ.
    515. 

  515.    SSGQ = SSG - SSGD
    516. 

  516. 

  517.    W=2.*1.*HGT
    518. 

  518.    VFGS=SVF
    519. 

  519.    VFGW=1.-SVF
    520. 

  520.    VFWG=(1.-SVF)*(1.-R)/W
    521. 

  521.    VFWS=VFWG
    522. 

  522.    VFWW=1.-2.*VFWG
    523. 

  523. 

  524. !-------------------------------------------------------------------------------
    525. 

  525. ! Convert unit from MKS to cgs
    526. 

  526. ! Renew surface and layer temperatures
    527. 

  527. !-------------------------------------------------------------------------------
    528. 

  528. 

  529.    SX=(SSGD+SSGQ)/697.7/60.  ! downward short wave radition [ly/min]
    530. 

  530.    SD=SSGD/697.7/60.         ! downward direct short wave radiation
    531. 

  531.    SQ=SSGQ/697.7/60.         ! downward diffuse short wave radiation
    532. 

  532.    RX=LLG/697.7/60.          ! downward long wave radiation
    533. 

  533.    RHO=RHOO*0.001            ! air density at first atmospheric level
    534. 

  534. 

  535.    TRP=TR
    536. 

  536.    TBP=TB
    537. 

  537.    TGP=TG
    538. 

  538.    TCP=TC
    539. 

  539.    QCP=QC
    540. 

  540. 

  541.    TAV=TA*(1.+0.61*QA)
    542. 

  542.    PS=RHOO*287.*TAV/100. ![hPa]
    543. 

  543. 

  544. !-------------------------------------------------------------------------------
    545. 

  545. ! Canopy wind
    546. 

  546. !-------------------------------------------------------------------------------
    547. 

  547. 

  548.    IF ( ZR + 2. < ZA ) THEN
    549. 

  549.      UR=UA*LOG((ZR-ZDC)/Z0C)/LOG((ZA-ZDC)/Z0C)
    550. 

  550.      ZC=0.7*ZR
    551. 

  551.      XLB=0.4*(ZR-ZDC)
    552. 

  552.      ! BB formulation from Inoue (1963)
    553. 

  553.      BB = 0.4 * ZR / ( XLB * alog( ( ZR - ZDC ) / Z0C ) )
    554. 

  554.      UC=UR*EXP(-BB*(1.-ZC/ZR))
    555. 

  555.    ELSE
    556. 

  556.      ! PRINT *, 'Warning ZR + 2m  is larger than the 1st WRF level' 
    557. 

  557.      ZC=ZA/2.
    558. 

  558.      UC=UA/2.
    559. 

  559.    END IF
    560. 

  560. 

  561. !-------------------------------------------------------------------------------
    562. 

  562. ! Net Short Wave Radiation at roof, wall, and road
    563. 

  563. !-------------------------------------------------------------------------------
    564. 

  564. 

  565.    SHADOW = .false.
    566. 

  566. !   SHADOW = .true.
    567. 

  567. 

  568.    IF (SSG > 0.0) THEN
    569. 

  569. 

  570.      IF(.NOT.SHADOW) THEN              ! no shadow effects model
    571. 

  571. 

  572.        SR1=SX*(1.-ALBR)
    573. 

  573.        SG1=SX*VFGS*(1.-ALBG)
    574. 

  574.        SB1=SX*VFWS*(1.-ALBB)
    575. 

  575.        SG2=SB1*ALBB/(1.-ALBB)*VFGW*(1.-ALBG)
    576. 

  576.        SB2=SG1*ALBG/(1.-ALBG)*VFWG*(1.-ALBB)
    577. 

  577. 

  578.      ELSE                              ! shadow effects model
    579. 

  579. 

  580.        FAI=XLAT*PI/180.
    581. 

  581. 

  582.        THEATAS=ABS(ASIN(COSZ))
    583. 

  583.        THEATAZ=ABS(ACOS(COSZ))
    584. 

  584. 

  585.        SNT=COS(DECLIN)*SIN(OMG)/COS(THEATAS)
    586. 

  586.        CNT=(COSZ*SIN(FAI)-SIN(DECLIN))/COS(THEATAS)/COS(FAI)
    587. 

  587. 

  588.        HOUI1=(SNT*COS(PI/8.)    -CNT*SIN(PI/8.))
    589. 

  589.        HOUI2=(SNT*COS(2.*PI/8.) -CNT*SIN(2.*PI/8.))
    590. 

  590.        HOUI3=(SNT*COS(3.*PI/8.) -CNT*SIN(3.*PI/8.))
    591. 

  591.        HOUI4=(SNT*COS(4.*PI/8.) -CNT*SIN(4.*PI/8.))
    592. 

  592.        HOUI5=(SNT*COS(5.*PI/8.) -CNT*SIN(5.*PI/8.))
    593. 

  593.        HOUI6=(SNT*COS(6.*PI/8.) -CNT*SIN(6.*PI/8.))
    594. 

  594.        HOUI7=(SNT*COS(7.*PI/8.) -CNT*SIN(7.*PI/8.))
    595. 

  595.        HOUI8=(SNT*COS(8.*PI/8.) -CNT*SIN(8.*PI/8.))
    596. 

  596. 

  597.        SLX1=HGT*ABS(TAN(THEATAZ))*ABS(HOUI1)
    598. 

  598.        SLX2=HGT*ABS(TAN(THEATAZ))*ABS(HOUI2)
    599. 

  599.        SLX3=HGT*ABS(TAN(THEATAZ))*ABS(HOUI3)
    600. 

  600.        SLX4=HGT*ABS(TAN(THEATAZ))*ABS(HOUI4)
    601. 

  601.        SLX5=HGT*ABS(TAN(THEATAZ))*ABS(HOUI5)
    602. 

  602.        SLX6=HGT*ABS(TAN(THEATAZ))*ABS(HOUI6)
    603. 

  603.        SLX7=HGT*ABS(TAN(THEATAZ))*ABS(HOUI7)
    604. 

  604.        SLX8=HGT*ABS(TAN(THEATAZ))*ABS(HOUI8)
    605. 

  605. 

  606.        IF(SLX1 > RW) SLX1=RW
    607. 

  607.        IF(SLX2 > RW) SLX2=RW
    608. 

  608.        IF(SLX3 > RW) SLX3=RW
    609. 

  609.        IF(SLX4 > RW) SLX4=RW
    610. 

  610.        IF(SLX5 > RW) SLX5=RW
    611. 

  611.        IF(SLX6 > RW) SLX6=RW
    612. 

  612.        IF(SLX7 > RW) SLX7=RW
    613. 

  613.        IF(SLX8 > RW) SLX8=RW
    614. 

  614. 

  615.        SLX=(SLX1+SLX2+SLX3+SLX4+SLX5+SLX6+SLX7+SLX8)/8.
    616. 

  616. 

  617.        SR1=SD*(1.-ALBR)+SQ*(1.-ALBR)
    618. 

  618.        SG1=SD*(RW-SLX)/RW*(1.-ALBG)+SQ*VFGS*(1.-ALBG)
    619. 

  619.        SB1=SD*SLX/W*(1.-ALBB)+SQ*VFWS*(1.-ALBB)
    620. 

  620.        SG2=SB1*ALBB/(1.-ALBB)*VFGW*(1.-ALBG)
    621. 

  621.        SB2=SG1*ALBG/(1.-ALBG)*VFWG*(1.-ALBB)
    622. 

  622. 

  623.      END IF
    624. 

  624. 

  625.      SR=SR1
    626. 

  626.      SG=SG1+SG2
    627. 

  627.      SB=SB1+SB2
    628. 

  628. 

  629.      SNET=R*SR+W*SB+RW*SG
    630. 

  630. 

  631.    ELSE
    632. 

  632. 

  633.      SR=0.
    634. 

  634.      SG=0.
    635. 

  635.      SB=0.
    636. 

  636.      SNET=0.
    637. 

  637. 

  638.    END IF
    639. 

  639. 

  640. !-------------------------------------------------------------------------------
    641. 

  641. ! Roof
    642. 

  642. !-------------------------------------------------------------------------------
    643. 

  643. 

  644. !-------------------------------------------------------------------------------
    645. 

  645. ! CHR, CDR, BETR
    646. 

  646. !-------------------------------------------------------------------------------
    647. 

  647. 

  648.    ! Z=ZA-ZDC
    649. 

  649.    ! BHR=LOG(Z0R/Z0HR)/0.4
    650. 

  650.    ! RIBR=(9.8*2./(TA+TRP))*(TA-TRP)*(Z+Z0R)/(UA*UA)
    651. 

  651.    ! CALL mos(XXXR,ALPHAR,CDR,BHR,RIBR,Z,Z0R,UA,TA,TRP,RHO)
    652. 

  652. 

  653.    ! Alternative option for MOST using SFCDIF routine from Noah
    654. 

  654.    ! Virtual temperatures needed by SFCDIF
    655. 

  655.    T1VR = TRP* (1.0+ 0.61 * QA) 
    656. 

  656.    TH2V = (TA + ( 0.0098 * ZA)) * (1.0+ 0.61 * QA)
    657. 

  657.   
    658. 

  658.    ! note that CHR_URB contains UA (=CHR_MOS*UA)
    659. 

  659.    RLMO_URB=0.0
    660. 

  660.    CALL SFCDIF_URB (ZA,Z0R,T1VR,TH2V,UA,AKANDA_URBAN,CMR_URB,CHR_URB,RLMO_URB,CDR)
    661. 

  661.    ALPHAR =  RHO*CP*CHR_URB
    662. 

  662.    CHR=ALPHAR/RHO/CP/UA
    663. 

  663. 

  664.    IF(RAIN > 1.) BETR=0.7  
    665. 

  665. 

  666.    IF (TS_SCHEME == 1) THEN 
    667. 

  667. 

  668. !-------------------------------------------------------------------------------
    669. 

  669. ! TR  Solving Non-Linear Equation by Newton-Rapson
    670. 

  670. ! TRL Solving Heat Equation by Tri Diagonal Matrix Algorithm  
    671. 

  671. !-------------------------------------------------------------------------------
    672. 

  672. !      TSC=TRP-273.15                          
    673. 

  673. !      ES=EXP(19.482-4303.4/(TSC+243.5))        ! WMO
    674. 

  674. !      ES=6.11*10.**(TETENA*TSC/(TETENB+TSC))   ! Tetens
    675. 

  675. !      DESDT=( 6.1078*(2500.-2.4*TSC)/ &        ! Tetens
    676. 

  676. !            (0.46151*(TSC+273.15)**2.) )*10.**(7.5*TSC/(237.3+TSC)) 
    677. 

  677. !      ES=6.11*EXP((2.5*10.**6./461.51)*(TRP-273.15)/(273.15*TRP) ) ! Clausius-Clapeyron 
    678. 

  678. !      DESDT=(2.5*10.**6./461.51)*ES/(TRP**2.)                      ! Clausius-Clapeyron
    679. 

  679. !      QS0R=0.622*ES/(PS-0.378*ES) 
    680. 

  680. !      DQS0RDTR = DESDT*0.622*PS/((PS-0.378*ES)**2.) 
    681. 

  681. !      DQS0RDTR = 17.269*(273.15-35.86)/((TRP-35.86)**2.)*QS0R
    682. 

  682. 

  683. !    TRP=350.
    684. 

  684. 

  685.      DO ITERATION=1,20
    686. 

  686. 

  687.        ES=6.11*EXP( (2.5*10.**6./461.51)*(TRP-273.15)/(273.15*TRP) )
    688. 

  688.        DESDT=(2.5*10.**6./461.51)*ES/(TRP**2.)
    689. 

  689.        QS0R=0.622*ES/(PS-0.378*ES) 
    690. 

  690.        DQS0RDTR = DESDT*0.622*PS/((PS-0.378*ES)**2.) 
    691. 

  691. 

  692.        RR=EPSR*(RX-SIG*(TRP**4.)/60.)
    693. 

  693.        HR=RHO*CP*CHR*UA*(TRP-TA)*100.
    694. 

  694.        ELER=RHO*EL*CHR*UA*BETR*(QS0R-QA)*100.
    695. 

  695.        G0R=AKSR*(TRP-TRL(1))/(DZR(1)/2.)
    696. 

  696. 

  697.        F = SR + RR - HR - ELER - G0R
    698. 

  698. 

  699.        DRRDTR = (-4.*EPSR*SIG*TRP**3.)/60.
    700. 

  700.        DHRDTR = RHO*CP*CHR*UA*100.
    701. 

  701.        DELERDTR = RHO*EL*CHR*UA*BETR*DQS0RDTR*100.
    702. 

  702.        DG0RDTR =  2.*AKSR/DZR(1)
    703. 

  703. 

  704.        DFDT = DRRDTR - DHRDTR - DELERDTR - DG0RDTR 
    705. 

  705.        DTR = F/DFDT
    706. 

  706. 

  707.        TR = TRP - DTR
    708. 

  708.        TRP = TR
    709. 

  709. 

  710.        IF( ABS(F) < 0.000001 .AND. ABS(DTR) < 0.000001 ) EXIT
    711. 

  711. 

  712.      END DO
    713. 

  713. 

  714. ! multi-layer heat equation model
    715. 

  715. 

  716.      CALL multi_layer(num_roof_layers,BOUNDR,G0R,CAPR,AKSR,TRL,DZR,DELT,TRLEND)
    717. 

  717. 

  718.    ELSE
    719. 

  719. 

  720.      ES=6.11*EXP( (2.5*10.**6./461.51)*(TRP-273.15)/(273.15*TRP) )
    721. 

  721.      QS0R=0.622*ES/(PS-0.378*ES)        
    722. 

  722. 

  723.      RR=EPSR*(RX-SIG*(TRP**4.)/60.)
    724. 

  724.      HR=RHO*CP*CHR*UA*(TRP-TA)*100.
    725. 

  725.      ELER=RHO*EL*CHR*UA*BETR*(QS0R-QA)*100.
    726. 

  726.      G0R=SR+RR-HR-ELER
    727. 

  727. 

  728.      CALL force_restore(CAPR,AKSR,DELT,SR,RR,HR,ELER,TRLEND,TRP,TR)
    729. 

  729. 

  730.      TRP=TR
    731. 

  731. 

  732.    END IF
    733. 

  733. 

  734.    FLXTHR=HR/RHO/CP/100.
    735. 

  735.    FLXHUMR=ELER/RHO/EL/100.
    736. 

  736. 

  737. !-------------------------------------------------------------------------------
    738. 

  738. !  Wall and Road 
    739. 

  739. !-------------------------------------------------------------------------------
    740. 

  740. 

  741. !-------------------------------------------------------------------------------
    742. 

  742. !  CHC, CHB, CDB, BETB, CHG, CDG, BETG
    743. 

  743. !-------------------------------------------------------------------------------
    744. 

  744. 

  745.    ! Z=ZA-ZDC
    746. 

  746.    ! BHC=LOG(Z0C/Z0HC)/0.4
    747. 

  747.    ! RIBC=(9.8*2./(TA+TCP))*(TA-TCP)*(Z+Z0C)/(UA*UA)
    748. 

  748.    !
    749. 

  749.    ! CALL mos(XXXC,ALPHAC,CDC,BHC,RIBC,Z,Z0C,UA,TA,TCP,RHO)
    750. 

  750.    ! Virtual temperatures needed by SFCDIF routine from Noah
    751. 

  751. 

  752.    T1VC = TCP* (1.0+ 0.61 * QA) 
    753. 

  753.    RLMO_URB=0.0
    754. 

  754.    CALL SFCDIF_URB(ZA,Z0C,T1VC,TH2V,UA,AKANDA_URBAN,CMC_URB,CHC_URB,RLMO_URB,CDC)
    755. 

  755.    ALPHAC =  RHO*CP*CHC_URB
    756. 

  756. 

  757.    IF (CH_SCHEME == 1) THEN 
    758. 

  758. 

  759.      Z=ZDC
    760. 

  760.      BHB=LOG(Z0B/Z0HB)/0.4
    761. 

  761.      BHG=LOG(Z0G/Z0HG)/0.4
    762. 

  762.      RIBB=(9.8*2./(TCP+TBP))*(TCP-TBP)*(Z+Z0B)/(UC*UC)
    763. 

  763.      RIBG=(9.8*2./(TCP+TGP))*(TCP-TGP)*(Z+Z0G)/(UC*UC)
    764. 

  764. 

  765.      CALL mos(XXXB,ALPHAB,CDB,BHB,RIBB,Z,Z0B,UC,TCP,TBP,RHO)
    766. 

  766.      CALL mos(XXXG,ALPHAG,CDG,BHG,RIBG,Z,Z0G,UC,TCP,TGP,RHO)
    767. 

  767. 

  768.    ELSE
    769. 

  769. 

  770.      ALPHAB=RHO*CP*(6.15+4.18*UC)/1200.
    771. 

  771.      IF(UC > 5.) ALPHAB=RHO*CP*(7.51*UC**0.78)/1200.
    772. 

  772.      ALPHAG=RHO*CP*(6.15+4.18*UC)/1200.
    773. 

  773.      IF(UC > 5.) ALPHAG=RHO*CP*(7.51*UC**0.78)/1200.
    774. 

  774. 

  775.    END IF
    776. 

  776. 

  777.    CHC=ALPHAC/RHO/CP/UA
    778. 

  778.    CHB=ALPHAB/RHO/CP/UC
    779. 

  779.    CHG=ALPHAG/RHO/CP/UC
    780. 

  780. 

  781.    BETB=0.0
    782. 

  782.    IF(RAIN > 1.) BETG=0.7
    783. 

  783. 

  784.    IF (TS_SCHEME == 1) THEN
    785. 

  785. 

  786. !-------------------------------------------------------------------------------
    787. 

  787. !  TB, TG  Solving Non-Linear Simultaneous Equation by Newton-Rapson
    788. 

  788. !  TBL,TGL Solving Heat Equation by Tri Diagonal Matrix Algorithm  
    789. 

  789. !-------------------------------------------------------------------------------
    790. 

  790. 

  791. !    TBP=350.
    792. 

  792. !    TGP=350.
    793. 

  793. 

  794.      DO ITERATION=1,20
    795. 

  795. 

  796.        ES=6.11*EXP( (2.5*10.**6./461.51)*(TBP-273.15)/(273.15*TBP) ) 
    797. 

  797.        DESDT=(2.5*10.**6./461.51)*ES/(TBP**2.)
    798. 

  798.        QS0B=0.622*ES/(PS-0.378*ES) 
    799. 

  799.        DQS0BDTB=DESDT*0.622*PS/((PS-0.378*ES)**2.)
    800. 

  800. 

  801.        ES=6.11*EXP( (2.5*10.**6./461.51)*(TGP-273.15)/(273.15*TGP) ) 
    802. 

  802.        DESDT=(2.5*10.**6./461.51)*ES/(TGP**2.)
    803. 

  803.        QS0G=0.622*ES/(PS-0.378*ES)        
    804. 

  804.        DQS0GDTG=DESDT*0.22*PS/((PS-0.378*ES)**2.) 
    805. 

  805. 

  806.        RG1=EPSG*( RX*VFGS          &
    807. 

  807.        +EPSB*VFGW*SIG*TBP**4./60.  &
    808. 

  808.        -SIG*TGP**4./60. )
    809. 

  809. 

  810.        RB1=EPSB*( RX*VFWS         &
    811. 

  811.        +EPSG*VFWG*SIG*TGP**4./60. &
    812. 

  812.        +EPSB*VFWW*SIG*TBP**4./60. &
    813. 

  813.        -SIG*TBP**4./60. )
    814. 

  814. 

  815.        RG2=EPSG*( (1.-EPSB)*(1.-SVF)*VFWS*RX                  &
    816. 

  816.        +(1.-EPSB)*(1.-SVF)*VFWG*EPSG*SIG*TGP**4./60.          &
    817. 

  817.        +EPSB*(1.-EPSB)*(1.-SVF)*(1.-2.*VFWS)*SIG*TBP**4./60. )
    818. 

  818. 

  819.        RB2=EPSB*( (1.-EPSG)*VFWG*VFGS*RX                          &
    820. 

  820.        +(1.-EPSG)*EPSB*VFGW*VFWG*SIG*(TBP**4.)/60.                &
    821. 

  821.        +(1.-EPSB)*VFWS*(1.-2.*VFWS)*RX                            &
    822. 

  822.        +(1.-EPSB)*VFWG*(1.-2.*VFWS)*EPSG*SIG*EPSG*TGP**4./60.     &
    823. 

  823.        +EPSB*(1.-EPSB)*(1.-2.*VFWS)*(1.-2.*VFWS)*SIG*TBP**4./60. )
    824. 

  824. 

  825.        RG=RG1+RG2
    826. 

  826.        RB=RB1+RB2
    827. 

  827. 

  828.        DRBDTB1=EPSB*(4.*EPSB*SIG*TB**3.*VFWW-4.*SIG*TB**3.)/60.
    829. 

  829.        DRBDTG1=EPSB*(4.*EPSG*SIG*TG**3.*VFWG)/60.
    830. 

  830.        DRBDTB2=EPSB*(4.*(1.-EPSG)*EPSB*SIG*TB**3.*VFGW*VFWG &
    831. 

  831.                +4.*EPSB*(1.-EPSB)*SIG*TB**3.*VFWW*VFWW)/60.
    832. 

  832.        DRBDTG2=EPSB*(4.*(1.-EPSB)*EPSG*SIG*TG**3.*VFWG*VFWW)/60.
    833. 

  833. 

  834.        DRGDTB1=EPSG*(4.*EPSB*SIG*TB**3.*VFGW)/60.
    835. 

  835.        DRGDTG1=EPSG*(-4.*SIG*TG**3.)/60.
    836. 

  836.        DRGDTB2=EPSG*(4.*EPSB*(1.-EPSB)*SIG*TB**3.*VFWW*VFGW)/60.
    837. 

  837.        DRGDTG2=EPSG*(4.*(1.-EPSB)*EPSG*SIG*TG**3.*VFWG*VFGW)/60.
    838. 

  838. 

  839.        DRBDTB=DRBDTB1+DRBDTB2
    840. 

  840.        DRBDTG=DRBDTG1+DRBDTG2
    841. 

  841.        DRGDTB=DRGDTB1+DRGDTB2
    842. 

  842.        DRGDTG=DRGDTG1+DRGDTG2
    843. 

  843. 

  844.        HB=RHO*CP*CHB*UC*(TBP-TCP)*100.
    845. 

  845.        HG=RHO*CP*CHG*UC*(TGP-TCP)*100.
    846. 

  846. 

  847.        DTCDTB=W*ALPHAB/(RW*ALPHAC+RW*ALPHAG+W*ALPHAB)
    848. 

  848.        DTCDTG=RW*ALPHAG/(RW*ALPHAC+RW*ALPHAG+W*ALPHAB)
    849. 

  849. 

  850.        DHBDTB=RHO*CP*CHB*UC*(1.-DTCDTB)*100.
    851. 

  851.        DHBDTG=RHO*CP*CHB*UC*(0.-DTCDTG)*100.
    852. 

  852.        DHGDTG=RHO*CP*CHG*UC*(1.-DTCDTG)*100.
    853. 

  853.        DHGDTB=RHO*CP*CHG*UC*(0.-DTCDTB)*100.
    854. 

  854. 

  855.        ELEB=RHO*EL*CHB*UC*BETB*(QS0B-QCP)*100.
    856. 

  856.        ELEG=RHO*EL*CHG*UC*BETG*(QS0G-QCP)*100.
    857. 

  857. 

  858.        DQCDTB=W*ALPHAB*BETB*DQS0BDTB/(RW*ALPHAC+RW*ALPHAG*BETG+W*ALPHAB*BETB)
    859. 

  859.        DQCDTG=RW*ALPHAG*BETG*DQS0GDTG/(RW*ALPHAC+RW*ALPHAG*BETG+W*ALPHAB*BETB)
    860. 

  860. 

  861.        DELEBDTB=RHO*EL*CHB*UC*BETB*(DQS0BDTB-DQCDTB)*100.
    862. 

  862.        DELEBDTG=RHO*EL*CHB*UC*BETB*(0.-DQCDTG)*100.
    863. 

  863.        DELEGDTG=RHO*EL*CHG*UC*BETG*(DQS0GDTG-DQCDTG)*100.
    864. 

  864.        DELEGDTB=RHO*EL*CHG*UC*BETG*(0.-DQCDTB)*100.
    865. 

  865. 

  866.        G0B=AKSB*(TBP-TBL(1))/(DZB(1)/2.)
    867. 

  867.        G0G=AKSG*(TGP-TGL(1))/(DZG(1)/2.)
    868. 

  868. 

  869.        DG0BDTB=2.*AKSB/DZB(1)
    870. 

  870.        DG0BDTG=0.
    871. 

  871.        DG0GDTG=2.*AKSG/DZG(1)
    872. 

  872.        DG0GDTB=0.
    873. 

  873. 

  874.        F = SB + RB - HB - ELEB - G0B
    875. 

  875.        FX = DRBDTB - DHBDTB - DELEBDTB - DG0BDTB 
    876. 

  876.        FY = DRBDTG - DHBDTG - DELEBDTG - DG0BDTG
    877. 

  877. 

  878.        GF = SG + RG - HG - ELEG - G0G
    879. 

  879.        GX = DRGDTB - DHGDTB - DELEGDTB - DG0GDTB
    880. 

  880.        GY = DRGDTG - DHGDTG - DELEGDTG - DG0GDTG 
    881. 

  881. 

  882.        DTB =  (GF*FY-F*GY)/(FX*GY-GX*FY)
    883. 

  883.        DTG = -(GF+GX*DTB)/GY
    884. 

  884. 

  885.        TB = TBP + DTB
    886. 

  886.        TG = TGP + DTG
    887. 

  887. 

  888.        TBP = TB
    889. 

  889.        TGP = TG
    890. 

  890. 

  891.        TC1=RW*ALPHAC+RW*ALPHAG+W*ALPHAB
    892. 

  892.        TC2=RW*ALPHAC*TA+RW*ALPHAG*TGP+W*ALPHAB*TBP
    893. 

  893.        TC=TC2/TC1
    894. 

  894. 

  895.        QC1=RW*ALPHAC+RW*ALPHAG*BETG+W*ALPHAB*BETB
    896. 

  896.        QC2=RW*ALPHAC*QA+RW*ALPHAG*BETG*QS0G+W*ALPHAB*BETB*QS0B
    897. 

  897.        QC=QC2/QC1
    898. 

  898. 

  899.        DTC=TCP - TC
    900. 

  900.        TCP=TC
    901. 

  901.        QCP=QC
    902. 

  902. 

  903.        IF( ABS(F) < 0.000001 .AND. ABS(DTB) < 0.000001 &
    904. 

  904.         .AND. ABS(GF) < 0.000001 .AND. ABS(DTG) < 0.000001 &
    905. 

  905.         .AND. ABS(DTC) < 0.000001) EXIT
    906. 

  906. 

  907.      END DO
    908. 

  908. 

  909.      CALL multi_layer(num_wall_layers,BOUNDB,G0B,CAPB,AKSB,TBL,DZB,DELT,TBLEND)
    910. 

  910. 

  911.      CALL multi_layer(num_road_layers,BOUNDG,G0G,CAPG,AKSG,TGL,DZG,DELT,TGLEND)
    912. 

  912. 

  913.    ELSE
    914. 

  914. 

  915. !-------------------------------------------------------------------------------
    916. 

  916. !  TB, TG  by Force-Restore Method 
    917. 

  917. !-------------------------------------------------------------------------------
    918. 

  918. 

  919.        ES=6.11*EXP((2.5*10.**6./461.51)*(TBP-273.15)/(273.15*TBP) )
    920. 

  920.        QS0B=0.622*ES/(PS-0.378*ES)       
    921. 

  921. 

  922.        ES=6.11*EXP((2.5*10.**6./461.51)*(TGP-273.15)/(273.15*TGP) )
    923. 

  923.        QS0G=0.622*ES/(PS-0.378*ES)        
    924. 

  924. 

  925.        RG1=EPSG*( RX*VFGS             &
    926. 

  926.        +EPSB*VFGW*SIG*TBP**4./60.     &
    927. 

  927.        -SIG*TGP**4./60. )
    928. 

  928. 

  929.        RB1=EPSB*( RX*VFWS &
    930. 

  930.        +EPSG*VFWG*SIG*TGP**4./60. &
    931. 

  931.        +EPSB*VFWW*SIG*TBP**4./60. &
    932. 

  932.        -SIG*TBP**4./60. )
    933. 

  933. 

  934.        RG2=EPSG*( (1.-EPSB)*(1.-SVF)*VFWS*RX                   &
    935. 

  935.        +(1.-EPSB)*(1.-SVF)*VFWG*EPSG*SIG*TGP**4./60.           &
    936. 

  936.        +EPSB*(1.-EPSB)*(1.-SVF)*(1.-2.*VFWS)*SIG*TBP**4./60. )
    937. 

  937. 

  938.        RB2=EPSB*( (1.-EPSG)*VFWG*VFGS*RX                          &
    939. 

  939.        +(1.-EPSG)*EPSB*VFGW*VFWG*SIG*(TBP**4.)/60.                &
    940. 

  940.        +(1.-EPSB)*VFWS*(1.-2.*VFWS)*RX                            &
    941. 

  941.        +(1.-EPSB)*VFWG*(1.-2.*VFWS)*EPSG*SIG*EPSG*TGP**4./60.     &
    942. 

  942.        +EPSB*(1.-EPSB)*(1.-2.*VFWS)*(1.-2.*VFWS)*SIG*TBP**4./60. )
    943. 

  943. 

  944.        RG=RG1+RG2
    945. 

  945.        RB=RB1+RB2
    946. 

  946. 

  947.        HB=RHO*CP*CHB*UC*(TBP-TCP)*100.
    948. 

  948.        ELEB=RHO*EL*CHB*UC*BETB*(QS0B-QCP)*100.
    949. 

  949.        G0B=SB+RB-HB-ELEB
    950. 

  950. 

  951.        HG=RHO*CP*CHG*UC*(TGP-TCP)*100.
    952. 

  952.        ELEG=RHO*EL*CHG*UC*BETG*(QS0G-QCP)*100.
    953. 

  953.        G0G=SG+RG-HG-ELEG
    954. 

  954. 

  955.        CALL force_restore(CAPB,AKSB,DELT,SB,RB,HB,ELEB,TBLEND,TBP,TB)
    956. 

  956.        CALL force_restore(CAPG,AKSG,DELT,SG,RG,HG,ELEG,TGLEND,TGP,TG)
    957. 

  957. 

  958.        TBP=TB
    959. 

  959.        TGP=TG
    960. 

  960. 

  961.        TC1=RW*ALPHAC+RW*ALPHAG+W*ALPHAB
    962. 

  962.        TC2=RW*ALPHAC*TA+RW*ALPHAG*TGP+W*ALPHAB*TBP
    963. 

  963.        TC=TC2/TC1
    964. 

  964. 

  965.        QC1=RW*ALPHAC+RW*ALPHAG*BETG+W*ALPHAB*BETB
    966. 

  966.        QC2=RW*ALPHAC*QA+RW*ALPHAG*BETG*QS0G+W*ALPHAB*BETB*QS0B
    967. 

  967.        QC=QC2/QC1
    968. 

  968. 

  969.        TCP=TC
    970. 

  970.        QCP=QC
    971. 

  971. 

  972.    END IF
    973. 

  973. 

  974. 

  975.    FLXTHB=HB/RHO/CP/100.
    976. 

  976.    FLXHUMB=ELEB/RHO/EL/100.
    977. 

  977.    FLXTHG=HG/RHO/CP/100.
    978. 

  978.    FLXHUMG=ELEG/RHO/EL/100.
    979. 

  979. 

  980. !-------------------------------------------------------------------------------
    981. 

  981. ! Total Fluxes from Urban Canopy
    982. 

  982. !-------------------------------------------------------------------------------
    983. 

  983. 

  984.    FLXUV  = ( R*CDR + RW*CDC )*UA*UA
    985. 

  985. ! Miao, 2007/01/17, cal. ah
    986. 

  986.    if(ahoption==1) then
    987. 

  987.      FLXTH  = ( R*FLXTHR  + W*FLXTHB  + RW*FLXTHG ) + AH/RHOO/CPP
    988. 

  988.    else
    989. 

  989.      FLXTH  = ( R*FLXTHR  + W*FLXTHB  + RW*FLXTHG )
    990. 

  990.    endif
    991. 

  991.    FLXHUM = ( R*FLXHUMR + W*FLXHUMB + RW*FLXHUMG )
    992. 

  992.    FLXG =   ( R*G0R + W*G0B + RW*G0G )
    993. 

  993.    LNET =     R*RR + W*RB + RW*RG
    994. 

  994. 

  995. !----------------------------------------------------------------------------
    996. 

  996. ! Convert Unit: FLUXES and u* T* q*  --> WRF
    997. 

  997. !----------------------------------------------------------------------------
    998. 

  998. 

  999.    SH    = FLXTH  * RHOO * CPP    ! Sensible heat flux          [W/m/m]
    1000. 

  1000.    LH    = FLXHUM * RHOO * ELL    ! Latent heat flux            [W/m/m]
    1001. 

  1001.    LH_KINEMATIC = FLXHUM * RHOO   ! Latent heat, Kinematic      [kg/m/m/s]
    1002. 

  1002.    LW    = LLG - (LNET*697.7*60.) ! Upward longwave radiation   [W/m/m]
    1003. 

  1003.    SW    = SSG - (SNET*697.7*60.) ! Upward shortwave radiation  [W/m/m]
    1004. 

  1004.    ALB   = 0.
    1005. 

  1005.    IF( ABS(SSG) > 0.0001) ALB = SW/SSG ! Effective albedo [-] 
    1006. 

  1006.    G = -FLXG*697.7*60.            ! [W/m/m]
    1007. 

  1007.    RN = (SNET+LNET)*697.7*60.     ! Net radiation [W/m/m]
    1008. 

  1008. 

  1009.    UST = SQRT(FLXUV)              ! u* [m/s]
    1010. 

  1010.    TST = -FLXTH/UST               ! T* [K]
    1011. 

  1011.    QST = -FLXHUM/UST              ! q* [-]
    1012. 

  1012. 

  1013. !------------------------------------------------------
    1014. 

  1014. !  diagnostic GRID AVERAGED  PSIM  PSIH  TS QS --> WRF
    1015. 

  1015. !------------------------------------------------------
    1016. 

  1016. 

  1017.    Z0 = Z0C 
    1018. 

  1018.    Z0H = Z0HC
    1019. 

  1019.    Z = ZA - ZDC
    1020. 

  1020. 

  1021.    XXX = 0.4*9.81*Z*TST/TA/UST/UST
    1022. 

  1022. 

  1023.    IF ( XXX >= 1. ) XXX = 1.
    1024. 

  1024.    IF ( XXX <= -5. ) XXX = -5.
    1025. 

  1025. 

  1026.    IF ( XXX > 0 ) THEN
    1027. 

  1027.      PSIM = -5. * XXX
    1028. 

  1028.      PSIH = -5. * XXX
    1029. 

  1029.    ELSE
    1030. 

  1030.      X = (1.-16.*XXX)**0.25
    1031. 

  1031.      PSIM = 2.*ALOG((1.+X)/2.) + ALOG((1.+X*X)/2.) - 2.*ATAN(X) + PI/2.
    1032. 

  1032.      PSIH = 2.*ALOG((1.+X*X)/2.)
    1033. 

  1033.    END IF
    1034. 

  1034. 

  1035.    GZ1OZ0 = ALOG(Z/Z0)
    1036. 

  1036.    CD = 0.4**2./(ALOG(Z/Z0)-PSIM)**2.
    1037. 

  1037. !
    1038. 

  1038. !m   CH = 0.4**2./(ALOG(Z/Z0)-PSIM)/(ALOG(Z/Z0H)-PSIH)
    1039. 

  1039. !m   CHS = 0.4*UST/(ALOG(Z/Z0H)-PSIH)
    1040. 

  1040. !m   TS = TA + FLXTH/CH/UA    ! surface potential temp (flux temp)
    1041. 

  1041. !m   QS = QA + FLXHUM/CH/UA   ! surface humidity
    1042. 

  1042. !
    1043. 

  1043.    TS = TA + FLXTH/CHS    ! surface potential temp (flux temp)
    1044. 

  1044.    QS = QA + FLXHUM/CHS   ! surface humidity
    1045. 

  1045. 

  1046. !-------------------------------------------------------
    1047. 

  1047. !  diagnostic  GRID AVERAGED  U10  V10  TH2  Q2 --> WRF
    1048. 

  1048. !-------------------------------------------------------
    1049. 

  1049. 

  1050.    XXX2 = (2./Z)*XXX
    1051. 

  1051.    IF ( XXX2 >= 1. ) XXX2 = 1.
    1052. 

  1052.    IF ( XXX2 <= -5. ) XXX2 = -5.
    1053. 

  1053. 

  1054.    IF ( XXX2 > 0 ) THEN
    1055. 

  1055.       PSIM2 = -5. * XXX2
    1056. 

  1056.       PSIH2 = -5. * XXX2
    1057. 

  1057.    ELSE
    1058. 

  1058.       X = (1.-16.*XXX2)**0.25
    1059. 

  1059.       PSIM2 = 2.*ALOG((1.+X)/2.) + ALOG((1.+X*X)/2.) - 2.*ATAN(X) + 2.*ATAN(1.)
    1060. 

  1060.       PSIH2 = 2.*ALOG((1.+X*X)/2.)
    1061. 

  1061.    END IF
    1062. 

  1062. !
    1063. 

  1063. !m   CHS2 = 0.4*UST/(ALOG(2./Z0H)-PSIH2)
    1064. 

  1064. !
    1065. 

  1065. 

  1066.    XXX10 = (10./Z)*XXX
    1067. 

  1067.    IF ( XXX10 >= 1. ) XXX10 = 1.
    1068. 

  1068.    IF ( XXX10 <= -5. ) XXX10 = -5.
    1069. 

  1069. 

  1070.    IF ( XXX10 > 0 ) THEN
    1071. 

  1071.       PSIM10 = -5. * XXX10
    1072. 

  1072.       PSIH10 = -5. * XXX10
    1073. 

  1073.    ELSE
    1074. 

  1074.       X = (1.-16.*XXX10)**0.25
    1075. 

  1075.       PSIM10 = 2.*ALOG((1.+X)/2.) + ALOG((1.+X*X)/2.) - 2.*ATAN(X) + 2.*ATAN(1.)
    1076. 

  1076.       PSIH10 = 2.*ALOG((1.+X*X)/2.)
    1077. 

  1077.    END IF
    1078. 

  1078. 

  1079.    PSIX = ALOG(Z/Z0) - PSIM
    1080. 

  1080.    PSIT = ALOG(Z/Z0H) - PSIH
    1081. 

  1081. 

  1082.    PSIX2 = ALOG(2./Z0) - PSIM2
    1083. 

  1083.    PSIT2 = ALOG(2./Z0H) - PSIH2
    1084. 

  1084. 

  1085.    PSIX10 = ALOG(10./Z0) - PSIM10
    1086. 

  1086.    PSIT10 = ALOG(10./Z0H) - PSIH10
    1087. 

  1087. 

  1088.    U10 = U1 * (PSIX10/PSIX)       ! u at 10 m [m/s]
    1089. 

  1089.    V10 = V1 * (PSIX10/PSIX)       ! v at 10 m [m/s]
    1090. 

  1090. 

  1091. !   TH2 = TS + (TA-TS)*(PSIT2/PSIT)   ! potential temp at 2 m [K]
    1092. 

  1092. !  TH2 = TS + (TA-TS)*(PSIT2/PSIT)   ! Fei: this seems to be temp (not potential) at 2 m [K]
    1093. 

  1093. !Fei: consistant with M-O theory
    1094. 

  1094.    TH2 = TS + (TA-TS) *(CHS/CHS2) 
    1095. 

  1095. 

  1096.    Q2 = QS + (QA-QS)*(PSIT2/PSIT)    ! humidity at 2 m       [-]
    1097. 

  1097. 

  1098. !   TS = (LW/SIG_SI/0.88)**0.25       ! Radiative temperature [K]
    1099. 

  1099. 

  1100.    END SUBROUTINE urban
    1101. 

  1101. !===============================================================================
    1102. 

  1102. !
    1103. 

  1103. !  mos
    1104. 

  1104. !
    1105. 

  1105. !===============================================================================
    1106. 

  1106.    SUBROUTINE mos(XXX,ALPHA,CD,B1,RIB,Z,Z0,UA,TA,TSF,RHO)
    1107. 

  1107. 

  1108. !  XXX:   z/L (requires iteration by Newton-Rapson method)
    1109. 

  1109. !  B1:    Stanton number
    1110. 

  1110. !  PSIM:  = PSIX of LSM
    1111. 

  1111. !  PSIH:  = PSIT of LSM
    1112. 

  1112. 

  1113.    IMPLICIT NONE
    1114. 

  1114. 

  1115.    REAL, PARAMETER     :: CP=0.24
    1116. 

  1116.    REAL, INTENT(IN)    :: B1, Z, Z0, UA, TA, TSF, RHO
    1117. 

  1117.    REAL, INTENT(OUT)   :: ALPHA, CD
    1118. 

  1118.    REAL, INTENT(INOUT) :: XXX, RIB
    1119. 

  1119.    REAL                :: XXX0, X, X0, FAIH, DPSIM, DPSIH
    1120. 

  1120.    REAL                :: F, DF, XXXP, US, TS, AL, XKB, DD, PSIM, PSIH
    1121. 

  1121.    INTEGER             :: NEWT
    1122. 

  1122.    INTEGER, PARAMETER  :: NEWT_END=10
    1123. 

  1123. 

  1124.    IF(RIB <= -15.) RIB=-15. 
    1125. 

  1125. 

  1126.    IF(RIB < 0.) THEN
    1127. 

  1127. 

  1128.       DO NEWT=1,NEWT_END
    1129. 

  1129. 

  1130.         IF(XXX >= 0.) XXX=-1.E-3
    1131. 

  1131. 

  1132.         XXX0=XXX*Z0/(Z+Z0)
    1133. 

  1133. 

  1134.         X=(1.-16.*XXX)**0.25
    1135. 

  1135.         X0=(1.-16.*XXX0)**0.25
    1136. 

  1136. 

  1137.         PSIM=ALOG((Z+Z0)/Z0) &
    1138. 

  1138.             -ALOG((X+1.)**2.*(X**2.+1.)) &
    1139. 

  1139.             +2.*ATAN(X) &
    1140. 

  1140.             +ALOG((X+1.)**2.*(X0**2.+1.)) &
    1141. 

  1141.             -2.*ATAN(X0)
    1142. 

  1142.         FAIH=1./SQRT(1.-16.*XXX)
    1143. 

  1143.         PSIH=ALOG((Z+Z0)/Z0)+0.4*B1 &
    1144. 

  1144.             -2.*ALOG(SQRT(1.-16.*XXX)+1.) &
    1145. 

  1145.             +2.*ALOG(SQRT(1.-16.*XXX0)+1.)
    1146. 

  1146. 

  1147.         DPSIM=(1.-16.*XXX)**(-0.25)/XXX &
    1148. 

  1148.              -(1.-16.*XXX0)**(-0.25)/XXX
    1149. 

  1149.         DPSIH=1./SQRT(1.-16.*XXX)/XXX &
    1150. 

  1150.              -1./SQRT(1.-16.*XXX0)/XXX
    1151. 

  1151. 

  1152.         F=RIB*PSIM**2./PSIH-XXX
    1153. 

  1153. 

  1154.         DF=RIB*(2.*DPSIM*PSIM*PSIH-DPSIH*PSIM**2.) &
    1155. 

  1155.           /PSIH**2.-1.
    1156. 

  1156. 

  1157.         XXXP=XXX
    1158. 

  1158.         XXX=XXXP-F/DF
    1159. 

  1159.         IF(XXX <= -10.) XXX=-10.
    1160. 

  1160. 

  1161.       END DO
    1162. 

  1162. 

  1163.    ELSE IF(RIB >= 0.142857) THEN
    1164. 

  1164. 

  1165.       XXX=0.714
    1166. 

  1166.       PSIM=ALOG((Z+Z0)/Z0)+7.*XXX
    1167. 

  1167.       PSIH=PSIM+0.4*B1
    1168. 

  1168. 

  1169.    ELSE
    1170. 

  1170. 

  1171.       AL=ALOG((Z+Z0)/Z0)
    1172. 

  1172.       XKB=0.4*B1
    1173. 

  1173.       DD=-4.*RIB*7.*XKB*AL+(AL+XKB)**2.
    1174. 

  1174.       IF(DD <= 0.) DD=0.
    1175. 

  1175.       XXX=(AL+XKB-2.*RIB*7.*AL-SQRT(DD))/(2.*(RIB*7.**2-7.))
    1176. 

  1176.       PSIM=ALOG((Z+Z0)/Z0)+7.*MIN(XXX,0.714)
    1177. 

  1177.       PSIH=PSIM+0.4*B1
    1178. 

  1178. 

  1179.    END IF
    1180. 

  1180. 

  1181.    US=0.4*UA/PSIM             ! u*
    1182. 

  1182.    IF(US <= 0.01) US=0.01
    1183. 

  1183.    TS=0.4*(TA-TSF)/PSIH       ! T*
    1184. 

  1184. 

  1185.    CD=US*US/UA**2.            ! CD
    1186. 

  1186.    ALPHA=RHO*CP*0.4*US/PSIH   ! RHO*CP*CH*U
    1187. 

  1187. 

  1188.    RETURN 
    1189. 

  1189.    END SUBROUTINE mos
    1190. 

  1190. !===============================================================================
    1191. 

  1191. !
    1192. 

  1192. !  louis79
    1193. 

  1193. !
    1194. 

  1194. !===============================================================================
    1195. 

  1195.    SUBROUTINE louis79(ALPHA,CD,RIB,Z,Z0,UA,RHO)
    1196. 

  1196. 

  1197.    IMPLICIT NONE
    1198. 

  1198. 

  1199.    REAL, PARAMETER     :: CP=0.24
    1200. 

  1200.    REAL, INTENT(IN)    :: Z, Z0, UA, RHO
    1201. 

  1201.    REAL, INTENT(OUT)   :: ALPHA, CD
    1202. 

  1202.    REAL, INTENT(INOUT) :: RIB
    1203. 

  1203.    REAL                :: A2, XX, CH, CMB, CHB
    1204. 

  1204. 

  1205.    A2=(0.4/ALOG(Z/Z0))**2.
    1206. 

  1206. 

  1207.    IF(RIB <= -15.) RIB=-15.
    1208. 

  1208. 

  1209.    IF(RIB >= 0.0) THEN
    1210. 

  1210.       IF(RIB >= 0.142857) THEN 
    1211. 

  1211.          XX=0.714
    1212. 

  1212.       ELSE 
    1213. 

  1213.          XX=RIB*LOG(Z/Z0)/(1.-7.*RIB)
    1214. 

  1214.       END IF 
    1215. 

  1215.       CH=0.16/0.74/(LOG(Z/Z0)+7.*MIN(XX,0.714))**2.
    1216. 

  1216.       CD=0.16/(LOG(Z/Z0)+7.*MIN(XX,0.714))**2.
    1217. 

  1217.    ELSE 
    1218. 

  1218.       CMB=7.4*A2*9.4*SQRT(Z/Z0)
    1219. 

  1219.       CHB=5.3*A2*9.4*SQRT(Z/Z0)
    1220. 

  1220.       CH=A2/0.74*(1.-9.4*RIB/(1.+CHB*SQRT(-RIB)))
    1221. 

  1221.       CD=A2*(1.-9.4*RIB/(1.+CHB*SQRT(-RIB)))
    1222. 

  1222.    END IF
    1223. 

  1223. 

  1224.    ALPHA=RHO*CP*CH*UA
    1225. 

  1225. 

  1226.    RETURN 
    1227. 

  1227.    END SUBROUTINE louis79
    1228. 

  1228. !===============================================================================
    1229. 

  1229. !
    1230. 

  1230. !  louis82
    1231. 

  1231. !
    1232. 

  1232. !===============================================================================
    1233. 

  1233.    SUBROUTINE louis82(ALPHA,CD,RIB,Z,Z0,UA,RHO)
    1234. 

  1234. 

  1235.    IMPLICIT NONE
    1236. 

  1236. 

  1237.    REAL, PARAMETER     :: CP=0.24
    1238. 

  1238.    REAL, INTENT(IN)    :: Z, Z0, UA, RHO
    1239. 

  1239.    REAL, INTENT(OUT)   :: ALPHA, CD
    1240. 

  1240.    REAL, INTENT(INOUT) :: RIB 
    1241. 

  1241.    REAL                :: A2, FM, FH, CH, CHH 
    1242. 

  1242. 

  1243.    A2=(0.4/ALOG(Z/Z0))**2.
    1244. 

  1244. 

  1245.    IF(RIB <= -15.) RIB=-15.
    1246. 

  1246. 

  1247.    IF(RIB >= 0.0) THEN 
    1248. 

  1248.       FM=1./((1.+(2.*5.*RIB)/SQRT(1.+5.*RIB)))
    1249. 

  1249.       FH=1./(1.+(3.*5.*RIB)*SQRT(1.+5.*RIB))
    1250. 

  1250.       CH=A2*FH
    1251. 

  1251.       CD=A2*FM
    1252. 

  1252.    ELSE 
    1253. 

  1253.       CHH=5.*3.*5.*A2*SQRT(Z/Z0)
    1254. 

  1254.       FM=1.-(2.*5.*RIB)/(1.+3.*5.*5.*A2*SQRT(Z/Z0+1.)*(-RIB))
    1255. 

  1255.       FH=1.-(3.*5.*RIB)/(1.+CHH*SQRT(-RIB))
    1256. 

  1256.       CH=A2*FH
    1257. 

  1257.       CD=A2*FM
    1258. 

  1258.    END IF
    1259. 

  1259. 

  1260.    ALPHA=RHO*CP*CH*UA
    1261. 

  1261. 

  1262.    RETURN
    1263. 

  1263.    END SUBROUTINE louis82
    1264. 

  1264. !===============================================================================
    1265. 

  1265. !
    1266. 

  1266. !  multi_layer
    1267. 

  1267. !
    1268. 

  1268. !===============================================================================
    1269. 

  1269.    SUBROUTINE multi_layer(KM,BOUND,G0,CAP,AKS,TSL,DZ,DELT,TSLEND)
    1270. 

  1270. 

  1271.    IMPLICIT NONE
    1272. 

  1272. 

  1273.    REAL, INTENT(IN)                   :: G0
    1274. 

  1274. 

  1275.    REAL, INTENT(IN)                   :: CAP
    1276. 

  1276. 

  1277.    REAL, INTENT(IN)                   :: AKS
    1278. 

  1278. 

  1279.    REAL, INTENT(IN)                   :: DELT      ! Time step [ s ]
    1280. 

  1280. 

  1281.    REAL, INTENT(IN)                   :: TSLEND
    1282. 

  1282. 

  1283.    INTEGER, INTENT(IN)                :: KM
    1284. 

  1284. 

  1285.    INTEGER, INTENT(IN)                :: BOUND
    1286. 

  1286. 

  1287.    REAL, DIMENSION(KM), INTENT(IN)    :: DZ
    1288. 

  1288. 

  1289.    REAL, DIMENSION(KM), INTENT(INOUT) :: TSL
    1290. 

  1290. 

  1291.    REAL, DIMENSION(KM)                :: A, B, C, D, X, P, Q
    1292. 

  1292. 

  1293.    REAL                               :: DZEND
    1294. 

  1294. 

  1295.    INTEGER                            :: K
    1296. 

  1296. 

  1297.    DZEND=DZ(KM)
    1298. 

  1298. 

  1299.    A(1) = 0.0
    1300. 

  1300. 

  1301.    B(1) = CAP*DZ(1)/DELT &
    1302. 

  1302.           +2.*AKS/(DZ(1)+DZ(2))
    1303. 

  1303.    C(1) = -2.*AKS/(DZ(1)+DZ(2))
    1304. 

  1304.    D(1) = CAP*DZ(1)/DELT*TSL(1) + G0
    1305. 

  1305. 

  1306.    DO K=2,KM-1
    1307. 

  1307.       A(K) = -2.*AKS/(DZ(K-1)+DZ(K))
    1308. 

  1308.       B(K) = CAP*DZ(K)/DELT + 2.*AKS/(DZ(K-1)+DZ(K)) + 2.*AKS/(DZ(K)+DZ(K+1)) 
    1309. 

  1309.       C(K) = -2.*AKS/(DZ(K)+DZ(K+1))
    1310. 

  1310.       D(K) = CAP*DZ(K)/DELT*TSL(K)
    1311. 

  1311.    END DO 
    1312. 

  1312. 

  1313.    IF(BOUND == 1) THEN                  ! Flux=0
    1314. 

  1314.       A(KM) = -2.*AKS/(DZ(KM-1)+DZ(KM))
    1315. 

  1315.       B(KM) = CAP*DZ(KM)/DELT + 2.*AKS/(DZ(KM-1)+DZ(KM))  
    1316. 

  1316.       C(KM) = 0.0
    1317. 

  1317.       D(KM) = CAP*DZ(KM)/DELT*TSL(KM)
    1318. 

  1318.    ELSE                                 ! T=constant
    1319. 

  1319.       A(KM) = -2.*AKS/(DZ(KM-1)+DZ(KM))
    1320. 

  1320.       B(KM) = CAP*DZ(KM)/DELT + 2.*AKS/(DZ(KM-1)+DZ(KM)) + 2.*AKS/(DZ(KM)+DZEND) 
    1321. 

  1321.       C(KM) = 0.0
    1322. 

  1322.       D(KM) = CAP*DZ(KM)/DELT*TSL(KM) + 2.*AKS*TSLEND/(DZ(KM)+DZEND)
    1323. 

  1323.    END IF 
    1324. 

  1324. 

  1325.    P(1) = -C(1)/B(1)
    1326. 

  1326.    Q(1) =  D(1)/B(1)
    1327. 

  1327. 

  1328.    DO K=2,KM
    1329. 

  1329.       P(K) = -C(K)/(A(K)*P(K-1)+B(K))
    1330. 

  1330.       Q(K) = (-A(K)*Q(K-1)+D(K))/(A(K)*P(K-1)+B(K))
    1331. 

  1331.    END DO 
    1332. 

  1332. 

  1333.    X(KM) = Q(KM)
    1334. 

  1334. 

  1335.    DO K=KM-1,1,-1
    1336. 

  1336.       X(K) = P(K)*X(K+1)+Q(K)
    1337. 

  1337.    END DO 
    1338. 

  1338. 

  1339.    DO K=1,KM
    1340. 

  1340.       TSL(K) = X(K)
    1341. 

  1341.    END DO 
    1342. 

  1342. 

  1343.    RETURN 
    1344. 

  1344.    END SUBROUTINE multi_layer
    1345. 

  1345. !==========…
    1346.

Large files files are truncated, but you can click here to view the full file