source: lmdz_wrf/trunk/WRFV3/phys/module_sf_urban.F @ 1371

Last change on this file since 1371 was 1, checked in by lfita, 11 years ago
  • -- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 96.6 KB
Line 
1MODULE module_sf_urban
2
3!===============================================================================
4!     Single-Layer Urban Canopy Model for WRF Noah-LSM
5!     Original Version: 2002/11/06 by Hiroyuki Kusaka
6!     Last Update:      2006/08/24 by Fei Chen and Mukul Tewari (NCAR/RAL) 
7!===============================================================================
8
9   CHARACTER(LEN=4)                :: LU_DATA_TYPE
10
11   INTEGER                         :: ICATE
12
13   REAL, ALLOCATABLE, DIMENSION(:) :: ZR_TBL
14   REAL, ALLOCATABLE, DIMENSION(:) :: Z0C_TBL
15   REAL, ALLOCATABLE, DIMENSION(:) :: Z0HC_TBL
16   REAL, ALLOCATABLE, DIMENSION(:) :: ZDC_TBL
17   REAL, ALLOCATABLE, DIMENSION(:) :: SVF_TBL
18   REAL, ALLOCATABLE, DIMENSION(:) :: R_TBL
19   REAL, ALLOCATABLE, DIMENSION(:) :: RW_TBL
20   REAL, ALLOCATABLE, DIMENSION(:) :: HGT_TBL
21   REAL, ALLOCATABLE, DIMENSION(:) :: AH_TBL
22   REAL, ALLOCATABLE, DIMENSION(:) :: BETR_TBL
23   REAL, ALLOCATABLE, DIMENSION(:) :: BETB_TBL
24   REAL, ALLOCATABLE, DIMENSION(:) :: BETG_TBL
25   REAL, ALLOCATABLE, DIMENSION(:) :: FRC_URB_TBL
26
27   REAL, ALLOCATABLE, DIMENSION(:) :: COP_TBL
28   REAL, ALLOCATABLE, DIMENSION(:) :: PWIN_TBL
29   REAL, ALLOCATABLE, DIMENSION(:) :: BETA_TBL
30   INTEGER, ALLOCATABLE, DIMENSION(:) :: SW_COND_TBL
31   REAL, ALLOCATABLE, DIMENSION(:) :: TIME_ON_TBL
32   REAL, ALLOCATABLE, DIMENSION(:) :: TIME_OFF_TBL
33   REAL, ALLOCATABLE, DIMENSION(:) :: TARGTEMP_TBL
34   REAL, ALLOCATABLE, DIMENSION(:) :: GAPTEMP_TBL
35   REAL, ALLOCATABLE, DIMENSION(:) :: TARGHUM_TBL
36   REAL, ALLOCATABLE, DIMENSION(:) :: GAPHUM_TBL
37   REAL, ALLOCATABLE, DIMENSION(:) :: PERFLO_TBL
38   REAL, ALLOCATABLE, DIMENSION(:) :: HSESF_TBL
39
40   REAL, ALLOCATABLE, DIMENSION(:) :: CAPR_TBL, CAPB_TBL, CAPG_TBL
41   REAL, ALLOCATABLE, DIMENSION(:) :: AKSR_TBL, AKSB_TBL, AKSG_TBL
42   REAL, ALLOCATABLE, DIMENSION(:) :: ALBR_TBL, ALBB_TBL, ALBG_TBL
43   REAL, ALLOCATABLE, DIMENSION(:) :: EPSR_TBL, EPSB_TBL, EPSG_TBL
44   REAL, ALLOCATABLE, DIMENSION(:) :: Z0R_TBL,  Z0B_TBL,  Z0G_TBL
45   REAL, ALLOCATABLE, DIMENSION(:) :: SIGMA_ZED_TBL
46   REAL, ALLOCATABLE, DIMENSION(:) :: Z0HB_TBL, Z0HG_TBL
47   REAL, ALLOCATABLE, DIMENSION(:) :: TRLEND_TBL, TBLEND_TBL, TGLEND_TBL
48   REAL, ALLOCATABLE, DIMENSION(:) :: AKANDA_URBAN_TBL
49!for BEP
50
51   ! MAXDIRS :: The maximum number of street directions we're allowed to define
52   INTEGER, PARAMETER :: MAXDIRS = 3
53   ! MAXHGTS :: The maximum number of building height bins we're allowed to define
54   INTEGER, PARAMETER :: MAXHGTS = 50
55
56   INTEGER, ALLOCATABLE, DIMENSION(:)   :: NUMDIR_TBL
57   REAL,    ALLOCATABLE, DIMENSION(:,:) :: STREET_DIRECTION_TBL
58   REAL,    ALLOCATABLE, DIMENSION(:,:) :: STREET_WIDTH_TBL
59   REAL,    ALLOCATABLE, DIMENSION(:,:) :: BUILDING_WIDTH_TBL
60   INTEGER, ALLOCATABLE, DIMENSION(:)   :: NUMHGT_TBL
61   REAL,    ALLOCATABLE, DIMENSION(:,:) :: HEIGHT_BIN_TBL
62   REAL,    ALLOCATABLE, DIMENSION(:,:) :: HPERCENT_BIN_TBL
63!end BEP
64   INTEGER                         :: BOUNDR_DATA,BOUNDB_DATA,BOUNDG_DATA
65   INTEGER                         :: CH_SCHEME_DATA, TS_SCHEME_DATA
66   INTEGER                         :: ahoption        ! Miao, 2007/01/17, cal. ah
67   REAL, DIMENSION(1:24)           :: ahdiuprf        ! ah diurnal profile, tloc: 1-24
68   REAL, DIMENSION(1:24)           :: hsequip_tbl
69
70   INTEGER                         :: allocate_status
71
72!   INTEGER                         :: num_roof_layers
73!   INTEGER                         :: num_wall_layers
74!   INTEGER                         :: num_road_layers
75
76   CONTAINS
77
78!===============================================================================
79!
80! Author:
81!      Hiroyuki KUSAKA, PhD
82!      University of Tsukuba, JAPAN
83!      (CRIEPI, NCAR/MMM visiting scientist, 2002-2004)
84!      kusaka@ccs.tsukuba.ac.jp
85!
86! Co-Researchers:
87!     Fei CHEN, PhD
88!      NCAR/RAP feichen@ucar.edu
89!     Mukul TEWARI, PhD
90!      NCAR/RAP mukul@ucar.edu
91!
92! Purpose:
93!     Calculate surface temeprature, fluxes, canopy air temperature, and canopy wind
94!
95! Subroutines:
96!     module_sf_urban
97!       |- urban
98!            |- read_param
99!            |- mos or jurges
100!            |- multi_layer or force_restore
101!       |- urban_param_init <-- URBPARM.TBL
102!       |- urban_var_init
103!
104! Input Data from WRF [MKS unit]:
105!
106!     UTYPE  [-]     : Urban type. 1=Commercial/Industrial; 2=High-intensity residential;
107!                    : 3=low-intensity residential
108!     TA     [K]     : Potential temperature at 1st wrf level (absolute temp)
109!     QA     [kg/kg] : Mixing ratio at 1st atmospheric level
110!     UA     [m/s]   : Wind speed at 1st atmospheric level
111!     SSG    [W/m/m] : Short wave downward radiation at a flat surface
112!                      Note this is the total of direct and diffusive solar
113!                      downward radiation. If without two components, the
114!                      single solar downward can be used instead.
115!                      SSG = SSGD + SSGQ
116!     LSOLAR [-]     : Indicating the input type of solar downward radiation
117!                      True: both direct and diffusive solar radiation
118!                      are available
119!                      False: only total downward ridiation is available.
120!     SSGD   [W/m/m] : Direct solar radiation at a flat surface
121!                      if SSGD is not available, one can assume a ratio SRATIO
122!                      (e.g., 0.7), so that SSGD = SRATIO*SSG
123!     SSGQ   [W/m/m] : Diffuse solar radiation at a flat surface
124!                      If SSGQ is not available, SSGQ = SSG - SSGD
125!     LLG    [W/m/m] : Long wave downward radiation at a flat surface
126!     RAIN   [mm/h]  : Precipitation
127!     RHOO   [kg/m/m/m] : Air density
128!     ZA     [m]     : First atmospheric level
129!                      as a lowest boundary condition
130!     DECLIN [rad]   : solar declination
131!     COSZ           : = sin(fai)*sin(del)+cos(fai)*cos(del)*cos(omg)
132!     OMG    [rad]   : solar hour angle
133!     XLAT   [deg]   : latitude
134!     DELT   [sec]   : Time step
135!     ZNT    [m]     : Roughnes length
136!
137! Output Data to WRF [MKS unit]:
138!
139!     TS  [K]            : Surface potential temperature (absolute temp)
140!     QS  [-]            : Surface humidity
141!
142!     SH  [W/m/m/]       : Sensible heat flux, = FLXTH*RHOO*CPP
143!     LH  [W/m/m]        : Latent heat flux, = FLXHUM*RHOO*ELL
144!     LH_INEMATIC [kg/m/m/sec]: Moisture Kinematic flux, = FLXHUM*RHOO
145!     SW  [W/m/m]        : Upward shortwave radiation flux,
146!                          = SSG-SNET*697.7*60. (697.7*60.=100.*100.*4.186)
147!     ALB [-]            : Time-varying albedo
148!     LW  [W/m/m]        : Upward longwave radiation flux,
149!                          = LNET*697.7*60.-LLG
150!     G   [W/m/m]        : Heat Flux into the Ground
151!     RN  [W/m/m]        : Net radiation
152!
153!     PSIM  [-]          : Diagnostic similarity stability function for momentum
154!     PSIH  [-]          : Diagnostic similarity stability function for heat
155!
156!     TC  [K]            : Diagnostic canopy air temperature
157!     QC  [-]            : Diagnostic canopy humidity
158!
159!     TH2 [K]            : Diagnostic potential temperature at 2 m
160!     Q2  [-]            : Diagnostic humidity at 2 m
161!     U10 [m/s]          : Diagnostic u wind component at 10 m
162!     V10 [m/s]          : Diagnostic v wind component at 10 m
163!
164!     CHS, CHS2 [m/s]    : CH*U at ZA, CH*U at 2 m (not used)
165!
166! Important parameters:
167!
168! Morphology of the urban canyon:
169! These parameters assigned in the URBPARM.TBL
170!
171!    ZR  [m]             : roof level (building height)
172!    SIGMA_ZED [m]       : Standard Deviation of roof height
173!    ROOF_WIDTH [m]      : roof (i.e., building) width
174!    ROAD_WIDTH [m]      : road width
175!
176! Parameters derived from the morphological terms above.
177! These parameters are computed in the code.
178!
179!    HGT [-]             : normalized building height
180!    SVF [-]             : sky view factor
181!    R   [-]             : Normalized roof width (a.k.a. "building coverage ratio")
182!    RW  [-]             : = 1 - R
183!    Z0C [m]             : Roughness length above canyon for momentum (1/10 of ZR)
184!    Z0HC [m]            : Roughness length above canyon for heat (1/10 of Z0C)
185!    ZDC [m]             : Zero plane displacement height (1/5 of ZR)
186!
187! Following parameter are assigned in run/URBPARM.TBL
188!
189!  AH  [ W m{-2} ]        : anthropogenic heat  ( W m{-2} in the table, converted internally to cal cm{-2} )
190!  CAPR[ J m{-3} K{-1} ]  : heat capacity of roof ( units converted in code to [ cal cm{-3} deg{-1} ] )
191!  CAPB[ J m{-3} K{-1} ]  : heat capacity of building wall ( units converted in code to [ cal cm{-3} deg{-1} ] )
192!  CAPG[ J m{-3} K{-1} ]  : heat capacity of road ( units converted in code to [ cal cm{-3} deg{-1} ] )
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!  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!  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!  ALBR [-]               : surface albedo of roof
197!  ALBB [-]               : surface albedo of building wall
198!  ALBG [-]               : surface albedo of road
199!  EPSR [-]               : surface emissivity of roof
200!  EPSB [-]               : surface emissivity of building wall
201!  EPSG [-]               : surface emissivity of road
202!  Z0B [m]                : roughness length for momentum of building wall (only for CH_SCHEME = 1)
203!  Z0G [m]                : roughness length for momentum of road (only for CH_SCHEME = 1)
204!  Z0HB [m]               : roughness length for heat of building wall (only for CH_SCHEME = 1)
205!  Z0HG [m]               : roughness length for heat of road
206!  num_roof_layers        : number of layers within roof
207!  num_wall_layers        : number of layers within building walls
208!  num_road_layers        : number of layers within below road surface
209!   NOTE: for now, these layers are defined as same as the number of soil layers in namelist.input
210!  DZR [cm]               : thickness of each roof layer
211!  DZB [cm]               : thickness of each building wall layer
212!  DZG [cm]               : thickness of each ground layer
213!  BOUNDR [integer 1 or 2] : Boundary Condition for Roof Layer Temp [1: Zero-Flux, 2: T = Constant]
214!  BOUNDB [integer 1 or 2] : Boundary Condition for Building Wall Layer Temp [1: Zero-Flux, 2: T = Constant]
215!  BOUNDG [integer 1 or 2] : Boundary Condition for Road Layer Temp [1: Zero-Flux, 2: T = Constant]
216!  TRLEND [K]             : lower boundary condition of roof temperature
217!  TBLEND [K]             : lower boundary condition of building temperature
218!  TGLEND [K]             : lower boundary condition of ground temperature
219!  CH_SCHEME [integer 1 or 2] : Sfc exchange scheme used for building wall and road
220!                         [1: M-O Similarity Theory,  2: Empirical Form (recommend)]
221!  TS_SCHEME [integer 1 or 2] : Scheme for computing surface temperature (for roof, wall, and road)
222!                         [1: 4-layer model,  2: Force-Restore method]
223!
224!for BEP
225!  numdir [ - ]             : Number of street directions defined for a particular urban category
226!  street_direction [ deg ] : Direction of streets for a particular urban category and a particular street direction
227!  street_width [ m ]       : Width of street for a particular urban category and a particular street direction
228!  building_width [ m ]     : Width of buildings for a particular urban category and a particular street direction
229!  numhgt [ - ]             : Number of building height levels defined for a particular urban category
230!  height_bin [ m ]         : Building height bins defined for a particular urban category.
231!  hpercent_bin [ % ]       : Percentage of a particular urban category populated by buildings of particular height bins
232!end BEP
233! Moved from URBPARM.TBL
234!
235!  BETR [-]            : minimum moisture availability of roof
236!  BETB [-]            : minimum moisture availability of building wall
237!  BETG [-]            : minimum moisture availability of road
238!  Z0R [m]                : roughness length for momentum of roof
239!  Z0HB [m]               : roughness length for heat of building wall (only for CH_SCHEME = 1)
240!  Z0HG [m]               : roughness length for heat of road
241!  num_roof_layers        : number of layers within roof
242!  num_wall_layers        : number of layers within building walls
243!  num_road_layers        : number of layers within below road surface
244!   NOTE: for now, these layers are defined as same as the number of soil layers in namelist.input
245!
246! References:
247!     Kusaka and Kimura (2004) J.Appl.Meteor., vol.43, p1899-1910
248!     Kusaka and Kimura (2004) J.Meteor.Soc.Japan, vol.82, p45-65
249!     Kusaka et al. (2001) Bound.-Layer Meteor., vol.101, p329-358
250!
251! History:
252!     2006/06          modified by H. Kusaka (Univ. Tsukuba), M. Tewari
253!     2005/10/26,      modified by Fei Chen, Mukul Tewari
254!     2003/07/21 WRF , modified  by H. Kusaka of CRIEPI (NCAR/MMM)
255!     2001/08/26 PhD , modified  by H. Kusaka of CRIEPI (Univ.Tsukuba)
256!     1999/08/25 LCM , developed by H. Kusaka of CRIEPI (Univ.Tsukuba)
257!
258!===============================================================================
259!
260!  subroutine urban:
261!
262!===============================================================================
263
264   SUBROUTINE urban(LSOLAR,                                           & ! L
265                    num_roof_layers,num_wall_layers,num_road_layers,  & ! I
266                    DZR,DZB,DZG,                                      & ! I
267                    UTYPE,TA,QA,UA,U1,V1,SSG,SSGD,SSGQ,LLG,RAIN,RHOO, & ! I
268                    ZA,DECLIN,COSZ,OMG,XLAT,DELT,ZNT,                 & ! I
269                    CHS, CHS2,                                        & ! I
270                    TR, TB, TG, TC, QC, UC,                           & ! H
271                    TRL,TBL,TGL,                                      & ! H
272                    XXXR, XXXB, XXXG, XXXC,                           & ! H
273                    TS,QS,SH,LH,LH_KINEMATIC,                         & ! O
274                    SW,ALB,LW,G,RN,PSIM,PSIH,                         & ! O
275                    GZ1OZ0,                                           & ! O
276                    CMR_URB,CHR_URB,CMC_URB,CHC_URB,                  & ! I/O
277                    U10,V10,TH2,Q2,UST                                & ! O
278                    )
279
280   IMPLICIT NONE
281
282   REAL, PARAMETER    :: CP=0.24          ! heat capacity of dry air  [cgs unit]
283   REAL, PARAMETER    :: EL=583.          ! latent heat of vaporation [cgs unit]
284   REAL, PARAMETER    :: SIG=8.17E-11     ! stefun bolzman constant   [cgs unit]
285   REAL, PARAMETER    :: SIG_SI=5.67E-8   !                           [MKS unit]
286   REAL, PARAMETER    :: AK=0.4           ! kalman const.                [-]
287   REAL, PARAMETER    :: PI=3.14159       ! pi                           [-]
288   REAL, PARAMETER    :: TETENA=7.5       ! const. of Tetens Equation    [-]
289   REAL, PARAMETER    :: TETENB=237.3     ! const. of Tetens Equation    [-]
290   REAL, PARAMETER    :: SRATIO=0.75      ! ratio between direct/total solar [-]
291
292   REAL, PARAMETER    :: CPP=1004.5       ! heat capacity of dry air    [J/K/kg]
293   REAL, PARAMETER    :: ELL=2.442E+06    ! latent heat of vaporization [J/kg]
294   REAL, PARAMETER    :: XKA=2.4E-5
295
296!-------------------------------------------------------------------------------
297! C: configuration variables
298!-------------------------------------------------------------------------------
299
300   LOGICAL, INTENT(IN) :: LSOLAR  ! logical    [true=both, false=SSG only]
301
302!  The following variables are also model configuration variables, but are
303!  defined in the URBAN.TBL and in the contains statement in the top of
304!  the module_urban_init, so we should not declare them here.
305
306  INTEGER, INTENT(IN) :: num_roof_layers
307  INTEGER, INTENT(IN) :: num_wall_layers
308  INTEGER, INTENT(IN) :: num_road_layers
309
310
311  REAL, INTENT(IN), DIMENSION(1:num_roof_layers) :: DZR ! grid interval of roof layers [cm]
312  REAL, INTENT(IN), DIMENSION(1:num_wall_layers) :: DZB ! grid interval of wall layers [cm]
313  REAL, INTENT(IN), DIMENSION(1:num_road_layers) :: DZG ! grid interval of road layers [cm]
314
315!-------------------------------------------------------------------------------
316! I: input variables from LSM to Urban
317!-------------------------------------------------------------------------------
318
319   INTEGER, INTENT(IN) :: UTYPE ! urban type [1=Commercial/Industrial, 2=High-intensity residential,
320                                ! 3=low-intensity residential]
321
322   REAL, INTENT(IN)    :: TA   ! potential temp at 1st atmospheric level [K]
323   REAL, INTENT(IN)    :: QA   ! mixing ratio at 1st atmospheric level  [kg/kg]
324   REAL, INTENT(IN)    :: UA   ! wind speed at 1st atmospheric level    [m/s]
325   REAL, INTENT(IN)    :: U1   ! u at 1st atmospheric level             [m/s]
326   REAL, INTENT(IN)    :: V1   ! v at 1st atmospheric level             [m/s]
327   REAL, INTENT(IN)    :: SSG  ! downward total short wave radiation    [W/m/m]
328   REAL, INTENT(IN)    :: LLG  ! downward long wave radiation           [W/m/m]
329   REAL, INTENT(IN)    :: RAIN ! precipitation                          [mm/h]
330   REAL, INTENT(IN)    :: RHOO ! air density                            [kg/m^3]
331   REAL, INTENT(IN)    :: ZA   ! first atmospheric level                [m]
332   REAL, INTENT(IN)    :: DECLIN ! solar declination                    [rad]
333   REAL, INTENT(IN)    :: COSZ ! sin(fai)*sin(del)+cos(fai)*cos(del)*cos(omg)
334   REAL, INTENT(IN)    :: OMG  ! solar hour angle                       [rad]
335
336   REAL, INTENT(IN)    :: XLAT ! latitude                               [deg]
337   REAL, INTENT(IN)    :: DELT ! time step                              [s]
338   REAL, INTENT(IN)    :: ZNT  ! roughness length                       [m]
339   REAL, INTENT(IN)    :: CHS,CHS2 ! CH*U at za and 2 m             [m/s]
340
341   REAL, INTENT(INOUT) :: SSGD ! downward direct short wave radiation   [W/m/m]
342   REAL, INTENT(INOUT) :: SSGQ ! downward diffuse short wave radiation  [W/m/m]
343   REAL, INTENT(INOUT) :: CMR_URB
344   REAL, INTENT(INOUT) :: CHR_URB
345   REAL, INTENT(INOUT) :: CMC_URB
346   REAL, INTENT(INOUT) :: CHC_URB
347
348!-------------------------------------------------------------------------------
349! O: output variables from Urban to LSM
350!-------------------------------------------------------------------------------
351
352   REAL, INTENT(OUT) :: TS     ! surface potential temperature    [K]
353   REAL, INTENT(OUT) :: QS     ! surface humidity                 [K]
354   REAL, INTENT(OUT) :: SH     ! sensible heat flux               [W/m/m]
355   REAL, INTENT(OUT) :: LH     ! latent heat flux                 [W/m/m]
356   REAL, INTENT(OUT) :: LH_KINEMATIC ! latent heat, kinetic     [kg/m/m/s]
357   REAL, INTENT(OUT) :: SW     ! upward short wave radiation flux [W/m/m]
358   REAL, INTENT(OUT) :: ALB    ! time-varying albedo            [fraction]
359   REAL, INTENT(OUT) :: LW     ! upward long wave radiation flux  [W/m/m]
360   REAL, INTENT(OUT) :: G      ! heat flux into the ground        [W/m/m]
361   REAL, INTENT(OUT) :: RN     ! net radition                     [W/m/m]
362   REAL, INTENT(OUT) :: PSIM   ! similality stability shear function for momentum
363   REAL, INTENT(OUT) :: PSIH   ! similality stability shear function for heat
364   REAL, INTENT(OUT) :: GZ1OZ0   
365   REAL, INTENT(OUT) :: U10    ! u at 10m                         [m/s]
366   REAL, INTENT(OUT) :: V10    ! u at 10m                         [m/s]
367   REAL, INTENT(OUT) :: TH2    ! potential temperature at 2 m     [K]
368   REAL, INTENT(OUT) :: Q2     ! humidity at 2 m                  [-]
369!m   REAL, INTENT(OUT) :: CHS,CHS2 ! CH*U at za and 2 m             [m/s]
370   REAL, INTENT(OUT) :: UST    ! friction velocity                [m/s]
371
372
373!-------------------------------------------------------------------------------
374! H: Historical (state) variables of Urban : LSM <--> Urban
375!-------------------------------------------------------------------------------
376
377! TR: roof temperature              [K];      TRP: at previous time step [K]
378! TB: building wall temperature     [K];      TBP: at previous time step [K]
379! TG: road temperature              [K];      TGP: at previous time step [K]
380! TC: urban-canopy air temperature  [K];      TCP: at previous time step [K]
381!                                                  (absolute temperature)
382! QC: urban-canopy air mixing ratio [kg/kg];  QCP: at previous time step [kg/kg]
383!
384! XXXR: Monin-Obkhov length for roof          [dimensionless]
385! XXXB: Monin-Obkhov length for building wall [dimensionless]
386! XXXG: Monin-Obkhov length for road          [dimensionless]
387! XXXC: Monin-Obkhov length for urban-canopy  [dimensionless]
388!
389! TRL, TBL, TGL: layer temperature [K] (absolute temperature)
390
391   REAL, INTENT(INOUT):: TR, TB, TG, TC, QC, UC
392   REAL, INTENT(INOUT):: XXXR, XXXB, XXXG, XXXC
393
394   REAL, DIMENSION(1:num_roof_layers), INTENT(INOUT) :: TRL
395   REAL, DIMENSION(1:num_wall_layers), INTENT(INOUT) :: TBL
396   REAL, DIMENSION(1:num_road_layers), INTENT(INOUT) :: TGL
397
398!-------------------------------------------------------------------------------
399! L:  Local variables from read_param
400!-------------------------------------------------------------------------------
401
402   REAL :: ZR, Z0C, Z0HC, ZDC, SVF, R, RW, HGT, AH
403   REAL :: SIGMA_ZED
404   REAL :: CAPR, CAPB, CAPG, AKSR, AKSB, AKSG, ALBR, ALBB, ALBG
405   REAL :: EPSR, EPSB, EPSG, Z0R,  Z0B,  Z0G,  Z0HB, Z0HG
406   REAL :: TRLEND,TBLEND,TGLEND
407   REAL :: T1VR, T1VC,TH2V
408   REAL :: RLMO_URB
409   REAL :: AKANDA_URBAN
410   
411   REAL :: TH2X                                                !m
412   
413   INTEGER :: BOUNDR, BOUNDB, BOUNDG
414   INTEGER :: CH_SCHEME, TS_SCHEME
415
416   LOGICAL :: SHADOW  ! [true=consider svf and shadow effects, false=consider svf effect only]
417
418!for BEP
419   INTEGER                        :: NUMDIR
420   REAL,    DIMENSION ( MAXDIRS ) :: STREET_DIRECTION
421   REAL,    DIMENSION ( MAXDIRS ) :: STREET_WIDTH
422   REAL,    DIMENSION ( MAXDIRS ) :: BUILDING_WIDTH
423   INTEGER                        :: NUMHGT
424   REAL,    DIMENSION ( MAXHGTS ) :: HEIGHT_BIN
425   REAL,    DIMENSION ( MAXHGTS ) :: HPERCENT_BIN
426!end BEP
427!-------------------------------------------------------------------------------
428! L: Local variables
429!-------------------------------------------------------------------------------
430
431   REAL :: BETR, BETB, BETG
432   REAL :: SX, SD, SQ, RX
433   REAL :: UR, ZC, XLB, BB
434   REAL :: Z, RIBB, RIBG, RIBC, BHR, BHB, BHG, BHC
435   REAL :: TSC, LNET, SNET, FLXUV, THG, FLXTH, FLXHUM, FLXG
436   REAL :: W, VFGS, VFGW, VFWG, VFWS, VFWW
437   REAL :: HOUI1, HOUI2, HOUI3, HOUI4, HOUI5, HOUI6, HOUI7, HOUI8
438   REAL :: SLX, SLX1, SLX2, SLX3, SLX4, SLX5, SLX6, SLX7, SLX8
439   REAL :: FLXTHR, FLXTHB, FLXTHG, FLXHUMR, FLXHUMB, FLXHUMG
440   REAL :: SR, SB, SG, RR, RB, RG
441   REAL :: SR1, SR2, SB1, SB2, SG1, SG2, RR1, RR2, RB1, RB2, RG1, RG2
442   REAL :: HR, HB, HG, ELER, ELEB, ELEG, G0R, G0B, G0G
443   REAL :: ALPHAC, ALPHAR, ALPHAB, ALPHAG
444   REAL :: CHC, CHR, CHB, CHG, CDC, CDR, CDB, CDG
445   REAL :: C1R, C1B, C1G, TE, TC1, TC2, QC1, QC2, QS0R, QS0B, QS0G,RHO,ES
446
447   REAL :: DESDT
448   REAL :: F
449   REAL :: DQS0RDTR
450   REAL :: DRRDTR, DHRDTR, DELERDTR, DG0RDTR
451   REAL :: DTR, DFDT
452   REAL :: FX, FY, GF, GX, GY
453   REAL :: DTCDTB, DTCDTG
454   REAL :: DQCDTB, DQCDTG
455   REAL :: DRBDTB1,  DRBDTG1,  DRBDTB2,  DRBDTG2
456   REAL :: DRGDTB1,  DRGDTG1,  DRGDTB2,  DRGDTG2
457   REAL :: DRBDTB,   DRBDTG,   DRGDTB,   DRGDTG
458   REAL :: DHBDTB,   DHBDTG,   DHGDTB,   DHGDTG
459   REAL :: DELEBDTB, DELEBDTG, DELEGDTG, DELEGDTB
460   REAL :: DG0BDTB,  DG0BDTG,  DG0GDTG,  DG0GDTB
461   REAL :: DQS0BDTB, DQS0GDTG
462   REAL :: DTB, DTG, DTC
463
464   REAL :: THEATAZ    ! Solar Zenith Angle [rad]
465   REAL :: THEATAS    ! = PI/2. - THETAZ
466   REAL :: FAI        ! Latitude [rad]
467   REAL :: CNT,SNT
468   REAL :: PS         ! Surface Pressure [hPa]
469   REAL :: TAV        ! Vertial Temperature [K]
470
471   REAL :: XXX, X, Z0, Z0H, CD, CH
472   REAL :: XXX2, PSIM2, PSIH2, XXX10, PSIM10, PSIH10
473   REAL :: PSIX, PSIT, PSIX2, PSIT2, PSIX10, PSIT10
474
475   REAL :: TRP, TBP, TGP, TCP, QCP, TST, QST
476
477   INTEGER :: iteration, K
478   INTEGER :: tloc
479
480!-------------------------------------------------------------------------------
481! Set parameters
482!-------------------------------------------------------------------------------
483
484! Miao, 2007/01/17, cal. ah
485   if(ahoption==1) then
486     tloc=mod(int(OMG/PI*180./15.+12.+0.5 ),24)
487     if(tloc==0) tloc=24
488   endif
489
490   CALL read_param(UTYPE,ZR,SIGMA_ZED,Z0C,Z0HC,ZDC,SVF,R,RW,HGT,  &
491                   AH,CAPR,CAPB,CAPG,AKSR,AKSB,AKSG,ALBR,ALBB,    &
492                   ALBG,EPSR,EPSB,EPSG,Z0R,Z0B,Z0G,Z0HB,Z0HG,     &
493                   BETR,BETB,BETG,TRLEND,TBLEND,TGLEND,           &
494!for BEP
495                   NUMDIR, STREET_DIRECTION, STREET_WIDTH,        &
496                   BUILDING_WIDTH, NUMHGT, HEIGHT_BIN,            &
497                   HPERCENT_BIN,                                  &
498!end BEP
499                   BOUNDR,BOUNDB,BOUNDG,CH_SCHEME,TS_SCHEME,      &
500                   AKANDA_URBAN)
501
502! Miao, 2007/01/17, cal. ah
503   if(ahoption==1) AH=AH*ahdiuprf(tloc)
504
505   IF( ZDC+Z0C+2. >= ZA) THEN
506    CALL wrf_error_fatal ("ZDC + Z0C + 2m is larger than the 1st WRF level "// &
507                           "Stop in subroutine urban - change ZDC and Z0C" )
508   END IF
509
510   IF(.NOT.LSOLAR) THEN
511     SSGD = SRATIO*SSG
512     SSGQ = SSG - SSGD
513   ENDIF
514   SSGD = SRATIO*SSG   ! No radiation scheme has SSGD and SSGQ.
515   SSGQ = SSG - SSGD
516
517   W=2.*1.*HGT
518   VFGS=SVF
519   VFGW=1.-SVF
520   VFWG=(1.-SVF)*(1.-R)/W
521   VFWS=VFWG
522   VFWW=1.-2.*VFWG
523
524!-------------------------------------------------------------------------------
525! Convert unit from MKS to cgs
526! Renew surface and layer temperatures
527!-------------------------------------------------------------------------------
528
529   SX=(SSGD+SSGQ)/697.7/60.  ! downward short wave radition [ly/min]
530   SD=SSGD/697.7/60.         ! downward direct short wave radiation
531   SQ=SSGQ/697.7/60.         ! downward diffuse short wave radiation
532   RX=LLG/697.7/60.          ! downward long wave radiation
533   RHO=RHOO*0.001            ! air density at first atmospheric level
534
535   TRP=TR
536   TBP=TB
537   TGP=TG
538   TCP=TC
539   QCP=QC
540
541   TAV=TA*(1.+0.61*QA)
542   PS=RHOO*287.*TAV/100. ![hPa]
543
544!-------------------------------------------------------------------------------
545! Canopy wind
546!-------------------------------------------------------------------------------
547
548   IF ( ZR + 2. < ZA ) THEN
549     UR=UA*LOG((ZR-ZDC)/Z0C)/LOG((ZA-ZDC)/Z0C)
550     ZC=0.7*ZR
551     XLB=0.4*(ZR-ZDC)
552     ! BB formulation from Inoue (1963)
553     BB = 0.4 * ZR / ( XLB * alog( ( ZR - ZDC ) / Z0C ) )
554     UC=UR*EXP(-BB*(1.-ZC/ZR))
555   ELSE
556     ! PRINT *, 'Warning ZR + 2m  is larger than the 1st WRF level'
557     ZC=ZA/2.
558     UC=UA/2.
559   END IF
560
561!-------------------------------------------------------------------------------
562! Net Short Wave Radiation at roof, wall, and road
563!-------------------------------------------------------------------------------
564
565   SHADOW = .false.
566!   SHADOW = .true.
567
568   IF (SSG > 0.0) THEN
569
570     IF(.NOT.SHADOW) THEN              ! no shadow effects model
571
572       SR1=SX*(1.-ALBR)
573       SG1=SX*VFGS*(1.-ALBG)
574       SB1=SX*VFWS*(1.-ALBB)
575       SG2=SB1*ALBB/(1.-ALBB)*VFGW*(1.-ALBG)
576       SB2=SG1*ALBG/(1.-ALBG)*VFWG*(1.-ALBB)
577
578     ELSE                              ! shadow effects model
579
580       FAI=XLAT*PI/180.
581
582       THEATAS=ABS(ASIN(COSZ))
583       THEATAZ=ABS(ACOS(COSZ))
584
585       SNT=COS(DECLIN)*SIN(OMG)/COS(THEATAS)
586       CNT=(COSZ*SIN(FAI)-SIN(DECLIN))/COS(THEATAS)/COS(FAI)
587
588       HOUI1=(SNT*COS(PI/8.)    -CNT*SIN(PI/8.))
589       HOUI2=(SNT*COS(2.*PI/8.) -CNT*SIN(2.*PI/8.))
590       HOUI3=(SNT*COS(3.*PI/8.) -CNT*SIN(3.*PI/8.))
591       HOUI4=(SNT*COS(4.*PI/8.) -CNT*SIN(4.*PI/8.))
592       HOUI5=(SNT*COS(5.*PI/8.) -CNT*SIN(5.*PI/8.))
593       HOUI6=(SNT*COS(6.*PI/8.) -CNT*SIN(6.*PI/8.))
594       HOUI7=(SNT*COS(7.*PI/8.) -CNT*SIN(7.*PI/8.))
595       HOUI8=(SNT*COS(8.*PI/8.) -CNT*SIN(8.*PI/8.))
596
597       SLX1=HGT*ABS(TAN(THEATAZ))*ABS(HOUI1)
598       SLX2=HGT*ABS(TAN(THEATAZ))*ABS(HOUI2)
599       SLX3=HGT*ABS(TAN(THEATAZ))*ABS(HOUI3)
600       SLX4=HGT*ABS(TAN(THEATAZ))*ABS(HOUI4)
601       SLX5=HGT*ABS(TAN(THEATAZ))*ABS(HOUI5)
602       SLX6=HGT*ABS(TAN(THEATAZ))*ABS(HOUI6)
603       SLX7=HGT*ABS(TAN(THEATAZ))*ABS(HOUI7)
604       SLX8=HGT*ABS(TAN(THEATAZ))*ABS(HOUI8)
605
606       IF(SLX1 > RW) SLX1=RW
607       IF(SLX2 > RW) SLX2=RW
608       IF(SLX3 > RW) SLX3=RW
609       IF(SLX4 > RW) SLX4=RW
610       IF(SLX5 > RW) SLX5=RW
611       IF(SLX6 > RW) SLX6=RW
612       IF(SLX7 > RW) SLX7=RW
613       IF(SLX8 > RW) SLX8=RW
614
615       SLX=(SLX1+SLX2+SLX3+SLX4+SLX5+SLX6+SLX7+SLX8)/8.
616
617       SR1=SD*(1.-ALBR)+SQ*(1.-ALBR)
618       SG1=SD*(RW-SLX)/RW*(1.-ALBG)+SQ*VFGS*(1.-ALBG)
619       SB1=SD*SLX/W*(1.-ALBB)+SQ*VFWS*(1.-ALBB)
620       SG2=SB1*ALBB/(1.-ALBB)*VFGW*(1.-ALBG)
621       SB2=SG1*ALBG/(1.-ALBG)*VFWG*(1.-ALBB)
622
623     END IF
624
625     SR=SR1
626     SG=SG1+SG2
627     SB=SB1+SB2
628
629     SNET=R*SR+W*SB+RW*SG
630
631   ELSE
632
633     SR=0.
634     SG=0.
635     SB=0.
636     SNET=0.
637
638   END IF
639
640!-------------------------------------------------------------------------------
641! Roof
642!-------------------------------------------------------------------------------
643
644!-------------------------------------------------------------------------------
645! CHR, CDR, BETR
646!-------------------------------------------------------------------------------
647
648   ! Z=ZA-ZDC
649   ! BHR=LOG(Z0R/Z0HR)/0.4
650   ! RIBR=(9.8*2./(TA+TRP))*(TA-TRP)*(Z+Z0R)/(UA*UA)
651   ! CALL mos(XXXR,ALPHAR,CDR,BHR,RIBR,Z,Z0R,UA,TA,TRP,RHO)
652
653   ! Alternative option for MOST using SFCDIF routine from Noah
654   ! Virtual temperatures needed by SFCDIF
655   T1VR = TRP* (1.0+ 0.61 * QA)
656   TH2V = (TA + ( 0.0098 * ZA)) * (1.0+ 0.61 * QA)
657 
658   ! note that CHR_URB contains UA (=CHR_MOS*UA)
659   RLMO_URB=0.0
660   CALL SFCDIF_URB (ZA,Z0R,T1VR,TH2V,UA,AKANDA_URBAN,CMR_URB,CHR_URB,RLMO_URB,CDR)
661   ALPHAR =  RHO*CP*CHR_URB
662   CHR=ALPHAR/RHO/CP/UA
663
664   IF(RAIN > 1.) BETR=0.7 
665
666   IF (TS_SCHEME == 1) THEN
667
668!-------------------------------------------------------------------------------
669! TR  Solving Non-Linear Equation by Newton-Rapson
670! TRL Solving Heat Equation by Tri Diagonal Matrix Algorithm 
671!-------------------------------------------------------------------------------
672!      TSC=TRP-273.15                         
673!      ES=EXP(19.482-4303.4/(TSC+243.5))        ! WMO
674!      ES=6.11*10.**(TETENA*TSC/(TETENB+TSC))   ! Tetens
675!      DESDT=( 6.1078*(2500.-2.4*TSC)/ &        ! Tetens
676!            (0.46151*(TSC+273.15)**2.) )*10.**(7.5*TSC/(237.3+TSC))
677!      ES=6.11*EXP((2.5*10.**6./461.51)*(TRP-273.15)/(273.15*TRP) ) ! Clausius-Clapeyron
678!      DESDT=(2.5*10.**6./461.51)*ES/(TRP**2.)                      ! Clausius-Clapeyron
679!      QS0R=0.622*ES/(PS-0.378*ES)
680!      DQS0RDTR = DESDT*0.622*PS/((PS-0.378*ES)**2.)
681!      DQS0RDTR = 17.269*(273.15-35.86)/((TRP-35.86)**2.)*QS0R
682
683!    TRP=350.
684
685     DO ITERATION=1,20
686
687       ES=6.11*EXP( (2.5*10.**6./461.51)*(TRP-273.15)/(273.15*TRP) )
688       DESDT=(2.5*10.**6./461.51)*ES/(TRP**2.)
689       QS0R=0.622*ES/(PS-0.378*ES)
690       DQS0RDTR = DESDT*0.622*PS/((PS-0.378*ES)**2.)
691
692       RR=EPSR*(RX-SIG*(TRP**4.)/60.)
693       HR=RHO*CP*CHR*UA*(TRP-TA)*100.
694       ELER=RHO*EL*CHR*UA*BETR*(QS0R-QA)*100.
695       G0R=AKSR*(TRP-TRL(1))/(DZR(1)/2.)
696
697       F = SR + RR - HR - ELER - G0R
698
699       DRRDTR = (-4.*EPSR*SIG*TRP**3.)/60.
700       DHRDTR = RHO*CP*CHR*UA*100.
701       DELERDTR = RHO*EL*CHR*UA*BETR*DQS0RDTR*100.
702       DG0RDTR =  2.*AKSR/DZR(1)
703
704       DFDT = DRRDTR - DHRDTR - DELERDTR - DG0RDTR
705       DTR = F/DFDT
706
707       TR = TRP - DTR
708       TRP = TR
709
710       IF( ABS(F) < 0.000001 .AND. ABS(DTR) < 0.000001 ) EXIT
711
712     END DO
713
714! multi-layer heat equation model
715
716     CALL multi_layer(num_roof_layers,BOUNDR,G0R,CAPR,AKSR,TRL,DZR,DELT,TRLEND)
717
718   ELSE
719
720     ES=6.11*EXP( (2.5*10.**6./461.51)*(TRP-273.15)/(273.15*TRP) )
721     QS0R=0.622*ES/(PS-0.378*ES)       
722
723     RR=EPSR*(RX-SIG*(TRP**4.)/60.)
724     HR=RHO*CP*CHR*UA*(TRP-TA)*100.
725     ELER=RHO*EL*CHR*UA*BETR*(QS0R-QA)*100.
726     G0R=SR+RR-HR-ELER
727
728     CALL force_restore(CAPR,AKSR,DELT,SR,RR,HR,ELER,TRLEND,TRP,TR)
729
730     TRP=TR
731
732   END IF
733
734   FLXTHR=HR/RHO/CP/100.
735   FLXHUMR=ELER/RHO/EL/100.
736
737!-------------------------------------------------------------------------------
738!  Wall and Road
739!-------------------------------------------------------------------------------
740
741!-------------------------------------------------------------------------------
742!  CHC, CHB, CDB, BETB, CHG, CDG, BETG
743!-------------------------------------------------------------------------------
744
745   ! Z=ZA-ZDC
746   ! BHC=LOG(Z0C/Z0HC)/0.4
747   ! RIBC=(9.8*2./(TA+TCP))*(TA-TCP)*(Z+Z0C)/(UA*UA)
748   !
749   ! CALL mos(XXXC,ALPHAC,CDC,BHC,RIBC,Z,Z0C,UA,TA,TCP,RHO)
750   ! Virtual temperatures needed by SFCDIF routine from Noah
751
752   T1VC = TCP* (1.0+ 0.61 * QA)
753   RLMO_URB=0.0
754   CALL SFCDIF_URB(ZA,Z0C,T1VC,TH2V,UA,AKANDA_URBAN,CMC_URB,CHC_URB,RLMO_URB,CDC)
755   ALPHAC =  RHO*CP*CHC_URB
756
757   IF (CH_SCHEME == 1) THEN
758
759     Z=ZDC
760     BHB=LOG(Z0B/Z0HB)/0.4
761     BHG=LOG(Z0G/Z0HG)/0.4
762     RIBB=(9.8*2./(TCP+TBP))*(TCP-TBP)*(Z+Z0B)/(UC*UC)
763     RIBG=(9.8*2./(TCP+TGP))*(TCP-TGP)*(Z+Z0G)/(UC*UC)
764
765     CALL mos(XXXB,ALPHAB,CDB,BHB,RIBB,Z,Z0B,UC,TCP,TBP,RHO)
766     CALL mos(XXXG,ALPHAG,CDG,BHG,RIBG,Z,Z0G,UC,TCP,TGP,RHO)
767
768   ELSE
769
770     ALPHAB=RHO*CP*(6.15+4.18*UC)/1200.
771     IF(UC > 5.) ALPHAB=RHO*CP*(7.51*UC**0.78)/1200.
772     ALPHAG=RHO*CP*(6.15+4.18*UC)/1200.
773     IF(UC > 5.) ALPHAG=RHO*CP*(7.51*UC**0.78)/1200.
774
775   END IF
776
777   CHC=ALPHAC/RHO/CP/UA
778   CHB=ALPHAB/RHO/CP/UC
779   CHG=ALPHAG/RHO/CP/UC
780
781   BETB=0.0
782   IF(RAIN > 1.) BETG=0.7
783
784   IF (TS_SCHEME == 1) THEN
785
786!-------------------------------------------------------------------------------
787!  TB, TG  Solving Non-Linear Simultaneous Equation by Newton-Rapson
788!  TBL,TGL Solving Heat Equation by Tri Diagonal Matrix Algorithm 
789!-------------------------------------------------------------------------------
790
791!    TBP=350.
792!    TGP=350.
793
794     DO ITERATION=1,20
795
796       ES=6.11*EXP( (2.5*10.**6./461.51)*(TBP-273.15)/(273.15*TBP) )
797       DESDT=(2.5*10.**6./461.51)*ES/(TBP**2.)
798       QS0B=0.622*ES/(PS-0.378*ES)
799       DQS0BDTB=DESDT*0.622*PS/((PS-0.378*ES)**2.)
800
801       ES=6.11*EXP( (2.5*10.**6./461.51)*(TGP-273.15)/(273.15*TGP) )
802       DESDT=(2.5*10.**6./461.51)*ES/(TGP**2.)
803       QS0G=0.622*ES/(PS-0.378*ES)       
804       DQS0GDTG=DESDT*0.22*PS/((PS-0.378*ES)**2.)
805
806       RG1=EPSG*( RX*VFGS          &
807       +EPSB*VFGW*SIG*TBP**4./60.  &
808       -SIG*TGP**4./60. )
809
810       RB1=EPSB*( RX*VFWS         &
811       +EPSG*VFWG*SIG*TGP**4./60. &
812       +EPSB*VFWW*SIG*TBP**4./60. &
813       -SIG*TBP**4./60. )
814
815       RG2=EPSG*( (1.-EPSB)*(1.-SVF)*VFWS*RX                  &
816       +(1.-EPSB)*(1.-SVF)*VFWG*EPSG*SIG*TGP**4./60.          &
817       +EPSB*(1.-EPSB)*(1.-SVF)*(1.-2.*VFWS)*SIG*TBP**4./60. )
818
819       RB2=EPSB*( (1.-EPSG)*VFWG*VFGS*RX                          &
820       +(1.-EPSG)*EPSB*VFGW*VFWG*SIG*(TBP**4.)/60.                &
821       +(1.-EPSB)*VFWS*(1.-2.*VFWS)*RX                            &
822       +(1.-EPSB)*VFWG*(1.-2.*VFWS)*EPSG*SIG*EPSG*TGP**4./60.     &
823       +EPSB*(1.-EPSB)*(1.-2.*VFWS)*(1.-2.*VFWS)*SIG*TBP**4./60. )
824
825       RG=RG1+RG2
826       RB=RB1+RB2
827
828       DRBDTB1=EPSB*(4.*EPSB*SIG*TB**3.*VFWW-4.*SIG*TB**3.)/60.
829       DRBDTG1=EPSB*(4.*EPSG*SIG*TG**3.*VFWG)/60.
830       DRBDTB2=EPSB*(4.*(1.-EPSG)*EPSB*SIG*TB**3.*VFGW*VFWG &
831               +4.*EPSB*(1.-EPSB)*SIG*TB**3.*VFWW*VFWW)/60.
832       DRBDTG2=EPSB*(4.*(1.-EPSB)*EPSG*SIG*TG**3.*VFWG*VFWW)/60.
833
834       DRGDTB1=EPSG*(4.*EPSB*SIG*TB**3.*VFGW)/60.
835       DRGDTG1=EPSG*(-4.*SIG*TG**3.)/60.
836       DRGDTB2=EPSG*(4.*EPSB*(1.-EPSB)*SIG*TB**3.*VFWW*VFGW)/60.
837       DRGDTG2=EPSG*(4.*(1.-EPSB)*EPSG*SIG*TG**3.*VFWG*VFGW)/60.
838
839       DRBDTB=DRBDTB1+DRBDTB2
840       DRBDTG=DRBDTG1+DRBDTG2
841       DRGDTB=DRGDTB1+DRGDTB2
842       DRGDTG=DRGDTG1+DRGDTG2
843
844       HB=RHO*CP*CHB*UC*(TBP-TCP)*100.
845       HG=RHO*CP*CHG*UC*(TGP-TCP)*100.
846
847       DTCDTB=W*ALPHAB/(RW*ALPHAC+RW*ALPHAG+W*ALPHAB)
848       DTCDTG=RW*ALPHAG/(RW*ALPHAC+RW*ALPHAG+W*ALPHAB)
849
850       DHBDTB=RHO*CP*CHB*UC*(1.-DTCDTB)*100.
851       DHBDTG=RHO*CP*CHB*UC*(0.-DTCDTG)*100.
852       DHGDTG=RHO*CP*CHG*UC*(1.-DTCDTG)*100.
853       DHGDTB=RHO*CP*CHG*UC*(0.-DTCDTB)*100.
854
855       ELEB=RHO*EL*CHB*UC*BETB*(QS0B-QCP)*100.
856       ELEG=RHO*EL*CHG*UC*BETG*(QS0G-QCP)*100.
857
858       DQCDTB=W*ALPHAB*BETB*DQS0BDTB/(RW*ALPHAC+RW*ALPHAG*BETG+W*ALPHAB*BETB)
859       DQCDTG=RW*ALPHAG*BETG*DQS0GDTG/(RW*ALPHAC+RW*ALPHAG*BETG+W*ALPHAB*BETB)
860
861       DELEBDTB=RHO*EL*CHB*UC*BETB*(DQS0BDTB-DQCDTB)*100.
862       DELEBDTG=RHO*EL*CHB*UC*BETB*(0.-DQCDTG)*100.
863       DELEGDTG=RHO*EL*CHG*UC*BETG*(DQS0GDTG-DQCDTG)*100.
864       DELEGDTB=RHO*EL*CHG*UC*BETG*(0.-DQCDTB)*100.
865
866       G0B=AKSB*(TBP-TBL(1))/(DZB(1)/2.)
867       G0G=AKSG*(TGP-TGL(1))/(DZG(1)/2.)
868
869       DG0BDTB=2.*AKSB/DZB(1)
870       DG0BDTG=0.
871       DG0GDTG=2.*AKSG/DZG(1)
872       DG0GDTB=0.
873
874       F = SB + RB - HB - ELEB - G0B
875       FX = DRBDTB - DHBDTB - DELEBDTB - DG0BDTB
876       FY = DRBDTG - DHBDTG - DELEBDTG - DG0BDTG
877
878       GF = SG + RG - HG - ELEG - G0G
879       GX = DRGDTB - DHGDTB - DELEGDTB - DG0GDTB
880       GY = DRGDTG - DHGDTG - DELEGDTG - DG0GDTG
881
882       DTB =  (GF*FY-F*GY)/(FX*GY-GX*FY)
883       DTG = -(GF+GX*DTB)/GY
884
885       TB = TBP + DTB
886       TG = TGP + DTG
887
888       TBP = TB
889       TGP = TG
890
891       TC1=RW*ALPHAC+RW*ALPHAG+W*ALPHAB
892       TC2=RW*ALPHAC*TA+RW*ALPHAG*TGP+W*ALPHAB*TBP
893       TC=TC2/TC1
894
895       QC1=RW*ALPHAC+RW*ALPHAG*BETG+W*ALPHAB*BETB
896       QC2=RW*ALPHAC*QA+RW*ALPHAG*BETG*QS0G+W*ALPHAB*BETB*QS0B
897       QC=QC2/QC1
898
899       DTC=TCP - TC
900       TCP=TC
901       QCP=QC
902
903       IF( ABS(F) < 0.000001 .AND. ABS(DTB) < 0.000001 &
904        .AND. ABS(GF) < 0.000001 .AND. ABS(DTG) < 0.000001 &
905        .AND. ABS(DTC) < 0.000001) EXIT
906
907     END DO
908
909     CALL multi_layer(num_wall_layers,BOUNDB,G0B,CAPB,AKSB,TBL,DZB,DELT,TBLEND)
910
911     CALL multi_layer(num_road_layers,BOUNDG,G0G,CAPG,AKSG,TGL,DZG,DELT,TGLEND)
912
913   ELSE
914
915!-------------------------------------------------------------------------------
916!  TB, TG  by Force-Restore Method
917!-------------------------------------------------------------------------------
918
919       ES=6.11*EXP((2.5*10.**6./461.51)*(TBP-273.15)/(273.15*TBP) )
920       QS0B=0.622*ES/(PS-0.378*ES)       
921
922       ES=6.11*EXP((2.5*10.**6./461.51)*(TGP-273.15)/(273.15*TGP) )
923       QS0G=0.622*ES/(PS-0.378*ES)       
924
925       RG1=EPSG*( RX*VFGS             &
926       +EPSB*VFGW*SIG*TBP**4./60.     &
927       -SIG*TGP**4./60. )
928
929       RB1=EPSB*( RX*VFWS &
930       +EPSG*VFWG*SIG*TGP**4./60. &
931       +EPSB*VFWW*SIG*TBP**4./60. &
932       -SIG*TBP**4./60. )
933
934       RG2=EPSG*( (1.-EPSB)*(1.-SVF)*VFWS*RX                   &
935       +(1.-EPSB)*(1.-SVF)*VFWG*EPSG*SIG*TGP**4./60.           &
936       +EPSB*(1.-EPSB)*(1.-SVF)*(1.-2.*VFWS)*SIG*TBP**4./60. )
937
938       RB2=EPSB*( (1.-EPSG)*VFWG*VFGS*RX                          &
939       +(1.-EPSG)*EPSB*VFGW*VFWG*SIG*(TBP**4.)/60.                &
940       +(1.-EPSB)*VFWS*(1.-2.*VFWS)*RX                            &
941       +(1.-EPSB)*VFWG*(1.-2.*VFWS)*EPSG*SIG*EPSG*TGP**4./60.     &
942       +EPSB*(1.-EPSB)*(1.-2.*VFWS)*(1.-2.*VFWS)*SIG*TBP**4./60. )
943
944       RG=RG1+RG2
945       RB=RB1+RB2
946
947       HB=RHO*CP*CHB*UC*(TBP-TCP)*100.
948       ELEB=RHO*EL*CHB*UC*BETB*(QS0B-QCP)*100.
949       G0B=SB+RB-HB-ELEB
950
951       HG=RHO*CP*CHG*UC*(TGP-TCP)*100.
952       ELEG=RHO*EL*CHG*UC*BETG*(QS0G-QCP)*100.
953       G0G=SG+RG-HG-ELEG
954
955       CALL force_restore(CAPB,AKSB,DELT,SB,RB,HB,ELEB,TBLEND,TBP,TB)
956       CALL force_restore(CAPG,AKSG,DELT,SG,RG,HG,ELEG,TGLEND,TGP,TG)
957
958       TBP=TB
959       TGP=TG
960
961       TC1=RW*ALPHAC+RW*ALPHAG+W*ALPHAB
962       TC2=RW*ALPHAC*TA+RW*ALPHAG*TGP+W*ALPHAB*TBP
963       TC=TC2/TC1
964
965       QC1=RW*ALPHAC+RW*ALPHAG*BETG+W*ALPHAB*BETB
966       QC2=RW*ALPHAC*QA+RW*ALPHAG*BETG*QS0G+W*ALPHAB*BETB*QS0B
967       QC=QC2/QC1
968
969       TCP=TC
970       QCP=QC
971
972   END IF
973
974
975   FLXTHB=HB/RHO/CP/100.
976   FLXHUMB=ELEB/RHO/EL/100.
977   FLXTHG=HG/RHO/CP/100.
978   FLXHUMG=ELEG/RHO/EL/100.
979
980!-------------------------------------------------------------------------------
981! Total Fluxes from Urban Canopy
982!-------------------------------------------------------------------------------
983
984   FLXUV  = ( R*CDR + RW*CDC )*UA*UA
985! Miao, 2007/01/17, cal. ah
986   if(ahoption==1) then
987     FLXTH  = ( R*FLXTHR  + W*FLXTHB  + RW*FLXTHG ) + AH/RHOO/CPP
988   else
989     FLXTH  = ( R*FLXTHR  + W*FLXTHB  + RW*FLXTHG )
990   endif
991   FLXHUM = ( R*FLXHUMR + W*FLXHUMB + RW*FLXHUMG )
992   FLXG =   ( R*G0R + W*G0B + RW*G0G )
993   LNET =     R*RR + W*RB + RW*RG
994
995!----------------------------------------------------------------------------
996! Convert Unit: FLUXES and u* T* q*  --> WRF
997!----------------------------------------------------------------------------
998
999   SH    = FLXTH  * RHOO * CPP    ! Sensible heat flux          [W/m/m]
1000   LH    = FLXHUM * RHOO * ELL    ! Latent heat flux            [W/m/m]
1001   LH_KINEMATIC = FLXHUM * RHOO   ! Latent heat, Kinematic      [kg/m/m/s]
1002   LW    = LLG - (LNET*697.7*60.) ! Upward longwave radiation   [W/m/m]
1003   SW    = SSG - (SNET*697.7*60.) ! Upward shortwave radiation  [W/m/m]
1004   ALB   = 0.
1005   IF( ABS(SSG) > 0.0001) ALB = SW/SSG ! Effective albedo [-]
1006   G = -FLXG*697.7*60.            ! [W/m/m]
1007   RN = (SNET+LNET)*697.7*60.     ! Net radiation [W/m/m]
1008
1009   UST = SQRT(FLXUV)              ! u* [m/s]
1010   TST = -FLXTH/UST               ! T* [K]
1011   QST = -FLXHUM/UST              ! q* [-]
1012
1013!------------------------------------------------------
1014!  diagnostic GRID AVERAGED  PSIM  PSIH  TS QS --> WRF
1015!------------------------------------------------------
1016
1017   Z0 = Z0C
1018   Z0H = Z0HC
1019   Z = ZA - ZDC
1020
1021   XXX = 0.4*9.81*Z*TST/TA/UST/UST
1022
1023   IF ( XXX >= 1. ) XXX = 1.
1024   IF ( XXX <= -5. ) XXX = -5.
1025
1026   IF ( XXX > 0 ) THEN
1027     PSIM = -5. * XXX
1028     PSIH = -5. * XXX
1029   ELSE
1030     X = (1.-16.*XXX)**0.25
1031     PSIM = 2.*ALOG((1.+X)/2.) + ALOG((1.+X*X)/2.) - 2.*ATAN(X) + PI/2.
1032     PSIH = 2.*ALOG((1.+X*X)/2.)
1033   END IF
1034
1035   GZ1OZ0 = ALOG(Z/Z0)
1036   CD = 0.4**2./(ALOG(Z/Z0)-PSIM)**2.
1037!
1038!m   CH = 0.4**2./(ALOG(Z/Z0)-PSIM)/(ALOG(Z/Z0H)-PSIH)
1039!m   CHS = 0.4*UST/(ALOG(Z/Z0H)-PSIH)
1040!m   TS = TA + FLXTH/CH/UA    ! surface potential temp (flux temp)
1041!m   QS = QA + FLXHUM/CH/UA   ! surface humidity
1042!
1043   TS = TA + FLXTH/CHS    ! surface potential temp (flux temp)
1044   QS = QA + FLXHUM/CHS   ! surface humidity
1045
1046!-------------------------------------------------------
1047!  diagnostic  GRID AVERAGED  U10  V10  TH2  Q2 --> WRF
1048!-------------------------------------------------------
1049
1050   XXX2 = (2./Z)*XXX
1051   IF ( XXX2 >= 1. ) XXX2 = 1.
1052   IF ( XXX2 <= -5. ) XXX2 = -5.
1053
1054   IF ( XXX2 > 0 ) THEN
1055      PSIM2 = -5. * XXX2
1056      PSIH2 = -5. * XXX2
1057   ELSE
1058      X = (1.-16.*XXX2)**0.25
1059      PSIM2 = 2.*ALOG((1.+X)/2.) + ALOG((1.+X*X)/2.) - 2.*ATAN(X) + 2.*ATAN(1.)
1060      PSIH2 = 2.*ALOG((1.+X*X)/2.)
1061   END IF
1062!
1063!m   CHS2 = 0.4*UST/(ALOG(2./Z0H)-PSIH2)
1064!
1065
1066   XXX10 = (10./Z)*XXX
1067   IF ( XXX10 >= 1. ) XXX10 = 1.
1068   IF ( XXX10 <= -5. ) XXX10 = -5.
1069
1070   IF ( XXX10 > 0 ) THEN
1071      PSIM10 = -5. * XXX10
1072      PSIH10 = -5. * XXX10
1073   ELSE
1074      X = (1.-16.*XXX10)**0.25
1075      PSIM10 = 2.*ALOG((1.+X)/2.) + ALOG((1.+X*X)/2.) - 2.*ATAN(X) + 2.*ATAN(1.)
1076      PSIH10 = 2.*ALOG((1.+X*X)/2.)
1077   END IF
1078
1079   PSIX = ALOG(Z/Z0) - PSIM
1080   PSIT = ALOG(Z/Z0H) - PSIH
1081
1082   PSIX2 = ALOG(2./Z0) - PSIM2
1083   PSIT2 = ALOG(2./Z0H) - PSIH2
1084
1085   PSIX10 = ALOG(10./Z0) - PSIM10
1086   PSIT10 = ALOG(10./Z0H) - PSIH10
1087
1088   U10 = U1 * (PSIX10/PSIX)       ! u at 10 m [m/s]
1089   V10 = V1 * (PSIX10/PSIX)       ! v at 10 m [m/s]
1090
1091!   TH2 = TS + (TA-TS)*(PSIT2/PSIT)   ! potential temp at 2 m [K]
1092!  TH2 = TS + (TA-TS)*(PSIT2/PSIT)   ! Fei: this seems to be temp (not potential) at 2 m [K]
1093!Fei: consistant with M-O theory
1094   TH2 = TS + (TA-TS) *(CHS/CHS2)
1095
1096   Q2 = QS + (QA-QS)*(PSIT2/PSIT)    ! humidity at 2 m       [-]
1097
1098!   TS = (LW/SIG_SI/0.88)**0.25       ! Radiative temperature [K]
1099
1100   END SUBROUTINE urban
1101!===============================================================================
1102!
1103!  mos
1104!
1105!===============================================================================
1106   SUBROUTINE mos(XXX,ALPHA,CD,B1,RIB,Z,Z0,UA,TA,TSF,RHO)
1107
1108!  XXX:   z/L (requires iteration by Newton-Rapson method)
1109!  B1:    Stanton number
1110!  PSIM:  = PSIX of LSM
1111!  PSIH:  = PSIT of LSM
1112
1113   IMPLICIT NONE
1114
1115   REAL, PARAMETER     :: CP=0.24
1116   REAL, INTENT(IN)    :: B1, Z, Z0, UA, TA, TSF, RHO
1117   REAL, INTENT(OUT)   :: ALPHA, CD
1118   REAL, INTENT(INOUT) :: XXX, RIB
1119   REAL                :: XXX0, X, X0, FAIH, DPSIM, DPSIH
1120   REAL                :: F, DF, XXXP, US, TS, AL, XKB, DD, PSIM, PSIH
1121   INTEGER             :: NEWT
1122   INTEGER, PARAMETER  :: NEWT_END=10
1123
1124   IF(RIB <= -15.) RIB=-15.
1125
1126   IF(RIB < 0.) THEN
1127
1128      DO NEWT=1,NEWT_END
1129
1130        IF(XXX >= 0.) XXX=-1.E-3
1131
1132        XXX0=XXX*Z0/(Z+Z0)
1133
1134        X=(1.-16.*XXX)**0.25
1135        X0=(1.-16.*XXX0)**0.25
1136
1137        PSIM=ALOG((Z+Z0)/Z0) &
1138            -ALOG((X+1.)**2.*(X**2.+1.)) &
1139            +2.*ATAN(X) &
1140            +ALOG((X+1.)**2.*(X0**2.+1.)) &
1141            -2.*ATAN(X0)
1142        FAIH=1./SQRT(1.-16.*XXX)
1143        PSIH=ALOG((Z+Z0)/Z0)+0.4*B1 &
1144            -2.*ALOG(SQRT(1.-16.*XXX)+1.) &
1145            +2.*ALOG(SQRT(1.-16.*XXX0)+1.)
1146
1147        DPSIM=(1.-16.*XXX)**(-0.25)/XXX &
1148             -(1.-16.*XXX0)**(-0.25)/XXX
1149        DPSIH=1./SQRT(1.-16.*XXX)/XXX &
1150             -1./SQRT(1.-16.*XXX0)/XXX
1151
1152        F=RIB*PSIM**2./PSIH-XXX
1153
1154        DF=RIB*(2.*DPSIM*PSIM*PSIH-DPSIH*PSIM**2.) &
1155          /PSIH**2.-1.
1156
1157        XXXP=XXX
1158        XXX=XXXP-F/DF
1159        IF(XXX <= -10.) XXX=-10.
1160
1161      END DO
1162
1163   ELSE IF(RIB >= 0.142857) THEN
1164
1165      XXX=0.714
1166      PSIM=ALOG((Z+Z0)/Z0)+7.*XXX
1167      PSIH=PSIM+0.4*B1
1168
1169   ELSE
1170
1171      AL=ALOG((Z+Z0)/Z0)
1172      XKB=0.4*B1
1173      DD=-4.*RIB*7.*XKB*AL+(AL+XKB)**2.
1174      IF(DD <= 0.) DD=0.
1175      XXX=(AL+XKB-2.*RIB*7.*AL-SQRT(DD))/(2.*(RIB*7.**2-7.))
1176      PSIM=ALOG((Z+Z0)/Z0)+7.*MIN(XXX,0.714)
1177      PSIH=PSIM+0.4*B1
1178
1179   END IF
1180
1181   US=0.4*UA/PSIM             ! u*
1182   IF(US <= 0.01) US=0.01
1183   TS=0.4*(TA-TSF)/PSIH       ! T*
1184
1185   CD=US*US/UA**2.            ! CD
1186   ALPHA=RHO*CP*0.4*US/PSIH   ! RHO*CP*CH*U
1187
1188   RETURN
1189   END SUBROUTINE mos
1190!===============================================================================
1191!
1192!  louis79
1193!
1194!===============================================================================
1195   SUBROUTINE louis79(ALPHA,CD,RIB,Z,Z0,UA,RHO)
1196
1197   IMPLICIT NONE
1198
1199   REAL, PARAMETER     :: CP=0.24
1200   REAL, INTENT(IN)    :: Z, Z0, UA, RHO
1201   REAL, INTENT(OUT)   :: ALPHA, CD
1202   REAL, INTENT(INOUT) :: RIB
1203   REAL                :: A2, XX, CH, CMB, CHB
1204
1205   A2=(0.4/ALOG(Z/Z0))**2.
1206
1207   IF(RIB <= -15.) RIB=-15.
1208
1209   IF(RIB >= 0.0) THEN
1210      IF(RIB >= 0.142857) THEN
1211         XX=0.714
1212      ELSE
1213         XX=RIB*LOG(Z/Z0)/(1.-7.*RIB)
1214      END IF
1215      CH=0.16/0.74/(LOG(Z/Z0)+7.*MIN(XX,0.714))**2.
1216      CD=0.16/(LOG(Z/Z0)+7.*MIN(XX,0.714))**2.
1217   ELSE
1218      CMB=7.4*A2*9.4*SQRT(Z/Z0)
1219      CHB=5.3*A2*9.4*SQRT(Z/Z0)
1220      CH=A2/0.74*(1.-9.4*RIB/(1.+CHB*SQRT(-RIB)))
1221      CD=A2*(1.-9.4*RIB/(1.+CHB*SQRT(-RIB)))
1222   END IF
1223
1224   ALPHA=RHO*CP*CH*UA
1225
1226   RETURN
1227   END SUBROUTINE louis79
1228!===============================================================================
1229!
1230!  louis82
1231!
1232!===============================================================================
1233   SUBROUTINE louis82(ALPHA,CD,RIB,Z,Z0,UA,RHO)
1234
1235   IMPLICIT NONE
1236
1237   REAL, PARAMETER     :: CP=0.24
1238   REAL, INTENT(IN)    :: Z, Z0, UA, RHO
1239   REAL, INTENT(OUT)   :: ALPHA, CD
1240   REAL, INTENT(INOUT) :: RIB
1241   REAL                :: A2, FM, FH, CH, CHH
1242
1243   A2=(0.4/ALOG(Z/Z0))**2.
1244
1245   IF(RIB <= -15.) RIB=-15.
1246
1247   IF(RIB >= 0.0) THEN
1248      FM=1./((1.+(2.*5.*RIB)/SQRT(1.+5.*RIB)))
1249      FH=1./(1.+(3.*5.*RIB)*SQRT(1.+5.*RIB))
1250      CH=A2*FH
1251      CD=A2*FM
1252   ELSE
1253      CHH=5.*3.*5.*A2*SQRT(Z/Z0)
1254      FM=1.-(2.*5.*RIB)/(1.+3.*5.*5.*A2*SQRT(Z/Z0+1.)*(-RIB))
1255      FH=1.-(3.*5.*RIB)/(1.+CHH*SQRT(-RIB))
1256      CH=A2*FH
1257      CD=A2*FM
1258   END IF
1259
1260   ALPHA=RHO*CP*CH*UA
1261
1262   RETURN
1263   END SUBROUTINE louis82
1264!===============================================================================
1265!
1266!  multi_layer
1267!
1268!===============================================================================
1269   SUBROUTINE multi_layer(KM,BOUND,G0,CAP,AKS,TSL,DZ,DELT,TSLEND)
1270
1271   IMPLICIT NONE
1272
1273   REAL, INTENT(IN)                   :: G0
1274
1275   REAL, INTENT(IN)                   :: CAP
1276
1277   REAL, INTENT(IN)                   :: AKS
1278
1279   REAL, INTENT(IN)                   :: DELT      ! Time step [ s ]
1280
1281   REAL, INTENT(IN)                   :: TSLEND
1282
1283   INTEGER, INTENT(IN)                :: KM
1284
1285   INTEGER, INTENT(IN)                :: BOUND
1286
1287   REAL, DIMENSION(KM), INTENT(IN)    :: DZ
1288
1289   REAL, DIMENSION(KM), INTENT(INOUT) :: TSL
1290
1291   REAL, DIMENSION(KM)                :: A, B, C, D, X, P, Q
1292
1293   REAL                               :: DZEND
1294
1295   INTEGER                            :: K
1296
1297   DZEND=DZ(KM)
1298
1299   A(1) = 0.0
1300
1301   B(1) = CAP*DZ(1)/DELT &
1302          +2.*AKS/(DZ(1)+DZ(2))
1303   C(1) = -2.*AKS/(DZ(1)+DZ(2))
1304   D(1) = CAP*DZ(1)/DELT*TSL(1) + G0
1305
1306   DO K=2,KM-1
1307      A(K) = -2.*AKS/(DZ(K-1)+DZ(K))
1308      B(K) = CAP*DZ(K)/DELT + 2.*AKS/(DZ(K-1)+DZ(K)) + 2.*AKS/(DZ(K)+DZ(K+1))
1309      C(K) = -2.*AKS/(DZ(K)+DZ(K+1))
1310      D(K) = CAP*DZ(K)/DELT*TSL(K)
1311   END DO
1312
1313   IF(BOUND == 1) THEN                  ! Flux=0
1314      A(KM) = -2.*AKS/(DZ(KM-1)+DZ(KM))
1315      B(KM) = CAP*DZ(KM)/DELT + 2.*AKS/(DZ(KM-1)+DZ(KM)) 
1316      C(KM) = 0.0
1317      D(KM) = CAP*DZ(KM)/DELT*TSL(KM)
1318   ELSE                                 ! T=constant
1319      A(KM) = -2.*AKS/(DZ(KM-1)+DZ(KM))
1320      B(KM) = CAP*DZ(KM)/DELT + 2.*AKS/(DZ(KM-1)+DZ(KM)) + 2.*AKS/(DZ(KM)+DZEND)
1321      C(KM) = 0.0
1322      D(KM) = CAP*DZ(KM)/DELT*TSL(KM) + 2.*AKS*TSLEND/(DZ(KM)+DZEND)
1323   END IF
1324
1325   P(1) = -C(1)/B(1)
1326   Q(1) =  D(1)/B(1)
1327
1328   DO K=2,KM
1329      P(K) = -C(K)/(A(K)*P(K-1)+B(K))
1330      Q(K) = (-A(K)*Q(K-1)+D(K))/(A(K)*P(K-1)+B(K))
1331   END DO
1332
1333   X(KM) = Q(KM)
1334
1335   DO K=KM-1,1,-1
1336      X(K) = P(K)*X(K+1)+Q(K)
1337   END DO
1338
1339   DO K=1,KM
1340      TSL(K) = X(K)
1341   END DO
1342
1343   RETURN
1344   END SUBROUTINE multi_layer
1345!===============================================================================
1346!
1347!  subroutine read_param
1348!
1349!===============================================================================
1350   SUBROUTINE read_param(UTYPE,                                        & ! in
1351                         ZR,SIGMA_ZED,Z0C,Z0HC,ZDC,SVF,R,RW,HGT,AH,    & ! out
1352                         CAPR,CAPB,CAPG,AKSR,AKSB,AKSG,ALBR,ALBB,ALBG, & ! out
1353                         EPSR,EPSB,EPSG,Z0R,Z0B,Z0G,Z0HB,Z0HG,         & ! out
1354                         BETR,BETB,BETG,TRLEND,TBLEND,TGLEND,          & ! out
1355!for BEP
1356                         NUMDIR, STREET_DIRECTION, STREET_WIDTH,       & ! out
1357                         BUILDING_WIDTH, NUMHGT, HEIGHT_BIN,           & ! out
1358                         HPERCENT_BIN,                                 & ! out
1359!end BEP
1360                         BOUNDR,BOUNDB,BOUNDG,CH_SCHEME,TS_SCHEME,     & ! out
1361                         AKANDA_URBAN)                                   ! out
1362
1363   INTEGER, INTENT(IN)  :: UTYPE
1364
1365   REAL, INTENT(OUT)    :: ZR,Z0C,Z0HC,ZDC,SVF,R,RW,HGT,AH,              &
1366                           CAPR,CAPB,CAPG,AKSR,AKSB,AKSG,ALBR,ALBB,ALBG, &
1367                           SIGMA_ZED,                                    &
1368                           EPSR,EPSB,EPSG,Z0R,Z0B,Z0G,Z0HB,Z0HG,         &
1369                           BETR,BETB,BETG,TRLEND,TBLEND,TGLEND
1370   REAL, INTENT(OUT)    :: AKANDA_URBAN
1371!for BEP
1372   INTEGER,                     INTENT(OUT) :: NUMDIR
1373   REAL,    DIMENSION(MAXDIRS), INTENT(OUT) :: STREET_DIRECTION
1374   REAL,    DIMENSION(MAXDIRS), INTENT(OUT) :: STREET_WIDTH
1375   REAL,    DIMENSION(MAXDIRS), INTENT(OUT) :: BUILDING_WIDTH
1376   INTEGER,                     INTENT(OUT) :: NUMHGT
1377   REAL,    DIMENSION(MAXHGTS), INTENT(OUT) :: HEIGHT_BIN
1378   REAL,    DIMENSION(MAXHGTS), INTENT(OUT) :: HPERCENT_BIN
1379   
1380!end BEP
1381
1382   INTEGER, INTENT(OUT) :: BOUNDR,BOUNDB,BOUNDG,CH_SCHEME,TS_SCHEME
1383
1384   ZR =     ZR_TBL(UTYPE)
1385   SIGMA_ZED = SIGMA_ZED_TBL(UTYPE)
1386   Z0C=     Z0C_TBL(UTYPE)
1387   Z0HC=    Z0HC_TBL(UTYPE)
1388   ZDC=     ZDC_TBL(UTYPE)
1389   SVF=     SVF_TBL(UTYPE)
1390   R=       R_TBL(UTYPE)
1391   RW=      RW_TBL(UTYPE)
1392   HGT=     HGT_TBL(UTYPE)
1393   AH=      AH_TBL(UTYPE)
1394   BETR=    BETR_TBL(UTYPE)
1395   BETB=    BETB_TBL(UTYPE)
1396   BETG=    BETG_TBL(UTYPE)
1397
1398!m   FRC_URB= FRC_URB_TBL(UTYPE)
1399
1400   CAPR=    CAPR_TBL(UTYPE)
1401   CAPB=    CAPB_TBL(UTYPE)
1402   CAPG=    CAPG_TBL(UTYPE)
1403   AKSR=    AKSR_TBL(UTYPE)
1404   AKSB=    AKSB_TBL(UTYPE)
1405   AKSG=    AKSG_TBL(UTYPE)
1406   ALBR=    ALBR_TBL(UTYPE)
1407   ALBB=    ALBB_TBL(UTYPE)
1408   ALBG=    ALBG_TBL(UTYPE)
1409   EPSR=    EPSR_TBL(UTYPE)
1410   EPSB=    EPSB_TBL(UTYPE)
1411   EPSG=    EPSG_TBL(UTYPE)
1412   Z0R=     Z0R_TBL(UTYPE)
1413   Z0B=     Z0B_TBL(UTYPE)
1414   Z0G=     Z0G_TBL(UTYPE)
1415   Z0HB=    Z0HB_TBL(UTYPE)
1416   Z0HG=    Z0HG_TBL(UTYPE)
1417   TRLEND=  TRLEND_TBL(UTYPE)
1418   TBLEND=  TBLEND_TBL(UTYPE)
1419   TGLEND=  TGLEND_TBL(UTYPE)
1420   BOUNDR=  BOUNDR_DATA
1421   BOUNDB=  BOUNDB_DATA
1422   BOUNDG=  BOUNDG_DATA
1423   CH_SCHEME = CH_SCHEME_DATA
1424   TS_SCHEME = TS_SCHEME_DATA
1425   AKANDA_URBAN = AKANDA_URBAN_TBL(UTYPE)
1426
1427!for BEP
1428
1429   STREET_DIRECTION = -1.E36
1430   STREET_WIDTH     = -1.E36
1431   BUILDING_WIDTH   = -1.E36
1432   HEIGHT_BIN       = -1.E36
1433   HPERCENT_BIN     = -1.E36
1434
1435   NUMDIR                     = NUMDIR_TBL ( UTYPE )
1436   STREET_DIRECTION(1:NUMDIR) = STREET_DIRECTION_TBL( 1:NUMDIR, UTYPE )
1437   STREET_WIDTH    (1:NUMDIR) = STREET_WIDTH_TBL    ( 1:NUMDIR, UTYPE )
1438   BUILDING_WIDTH  (1:NUMDIR) = BUILDING_WIDTH_TBL  ( 1:NUMDIR, UTYPE )
1439   NUMHGT                     = NUMHGT_TBL ( UTYPE )
1440   HEIGHT_BIN      (1:NUMHGT) = HEIGHT_BIN_TBL   ( 1:NUMHGT , UTYPE )
1441   HPERCENT_BIN    (1:NUMHGT) = HPERCENT_BIN_TBL ( 1:NUMHGT , UTYPE )
1442
1443!end BEP
1444   END SUBROUTINE read_param
1445!===============================================================================
1446!
1447! subroutine urban_param_init: Read parameters from URBPARM.TBL
1448!
1449!===============================================================================
1450   SUBROUTINE urban_param_init(DZR,DZB,DZG,num_soil_layers, &
1451                               sf_urban_physics)
1452!                              num_roof_layers,num_wall_layers,num_road_layers)
1453
1454   IMPLICIT NONE
1455
1456   INTEGER, INTENT(IN) :: num_soil_layers
1457
1458!   REAL, DIMENSION(1:num_roof_layers), INTENT(INOUT) :: DZR
1459!   REAL, DIMENSION(1:num_wall_layers), INTENT(INOUT) :: DZB
1460!   REAL, DIMENSION(1:num_road_layers), INTENT(INOUT) :: DZG
1461   REAL, DIMENSION(1:num_soil_layers), INTENT(INOUT) :: DZR
1462   REAL, DIMENSION(1:num_soil_layers), INTENT(INOUT) :: DZB
1463   REAL, DIMENSION(1:num_soil_layers), INTENT(INOUT) :: DZG
1464   INTEGER,                            INTENT(IN)    :: SF_URBAN_PHYSICS
1465
1466   INTEGER :: LC, K
1467   INTEGER :: IOSTATUS, ALLOCATE_STATUS
1468   INTEGER :: num_roof_layers
1469   INTEGER :: num_wall_layers
1470   INTEGER :: num_road_layers
1471   INTEGER :: dummy
1472   REAL    :: DHGT, HGT, VFWS, VFGS
1473
1474   REAL, allocatable, dimension(:) :: ROOF_WIDTH
1475   REAL, allocatable, dimension(:) :: ROAD_WIDTH
1476
1477   character(len=512) :: string
1478   character(len=128) :: name
1479   integer :: indx
1480
1481   real, parameter :: VonK = 0.4
1482   real :: lambda_p
1483   real :: lambda_f
1484   real :: Cd
1485   real :: alpha_macd
1486   real :: beta_macd
1487   real :: lambda_fr
1488
1489
1490!for BEP
1491   real :: dummy_hgt
1492   real :: dummy_pct
1493   real :: pctsum
1494!end BEP
1495   num_roof_layers = num_soil_layers
1496   num_wall_layers = num_soil_layers
1497   num_road_layers = num_soil_layers
1498
1499
1500   ICATE=0
1501
1502   OPEN (UNIT=11,                &
1503         FILE='URBPARM.TBL', &
1504         ACCESS='SEQUENTIAL',    &
1505         STATUS='OLD',           &
1506         ACTION='READ',          &
1507         POSITION='REWIND',      &
1508         IOSTAT=IOSTATUS)
1509
1510   IF (IOSTATUS > 0) THEN
1511   CALL wrf_error_fatal('ERROR OPEN URBPARM.TBL')
1512   ENDIF
1513
1514   READLOOP : do
1515      read(11,'(A512)', iostat=iostatus) string
1516      if (iostatus /= 0) exit READLOOP
1517      if (string(1:1) == "#") cycle READLOOP
1518      if (trim(string) == "") cycle READLOOP
1519      indx = index(string,":")
1520      if (indx <= 0) cycle READLOOP
1521      name = trim(adjustl(string(1:indx-1)))
1522     
1523      ! Here are the variables we expect to be defined in the URBPARM.TBL:
1524      if (name == "Number of urban categories") then
1525         read(string(indx+1:),*) icate
1526         IF (.not. ALLOCATED(ZR_TBL)) then
1527            ALLOCATE( ZR_TBL(ICATE), stat=allocate_status )
1528            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating ZR_TBL in urban_param_init')
1529            ALLOCATE( SIGMA_ZED_TBL(ICATE), stat=allocate_status )
1530            if(allocate_status /= 0)CALL wrf_error_fatal('Error allocating SIGMA_ZED_TBL in urban_param_init')
1531            ALLOCATE( Z0C_TBL(ICATE), stat=allocate_status )
1532            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating Z0C_TBL in urban_param_init')
1533            ALLOCATE( Z0HC_TBL(ICATE), stat=allocate_status )
1534            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating Z0HC_TBL in urban_param_init')
1535            ALLOCATE( ZDC_TBL(ICATE), stat=allocate_status )
1536            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating ZDC_TBL in urban_param_init')
1537            ALLOCATE( SVF_TBL(ICATE), stat=allocate_status )
1538            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating SVF_TBL in urban_param_init')
1539            ALLOCATE( R_TBL(ICATE), stat=allocate_status )
1540            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating R_TBL in urban_param_init')
1541            ALLOCATE( RW_TBL(ICATE), stat=allocate_status )
1542            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating RW_TBL in urban_param_init')
1543            ALLOCATE( HGT_TBL(ICATE), stat=allocate_status )
1544            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating HGT_TBL in urban_param_init')
1545            ALLOCATE( AH_TBL(ICATE), stat=allocate_status )
1546            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating AH_TBL in urban_param_init')
1547            ALLOCATE( BETR_TBL(ICATE), stat=allocate_status )
1548            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating BETR_TBL in urban_param_init')
1549            ALLOCATE( BETB_TBL(ICATE), stat=allocate_status )
1550            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating BETB_TBL in urban_param_init')
1551            ALLOCATE( BETG_TBL(ICATE), stat=allocate_status )
1552            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating BETG_TBL in urban_param_init')
1553            ALLOCATE( CAPR_TBL(ICATE), stat=allocate_status )
1554            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating CAPR_TBL in urban_param_init')
1555            ALLOCATE( CAPB_TBL(ICATE), stat=allocate_status )
1556            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating CAPB_TBL in urban_param_init')
1557            ALLOCATE( CAPG_TBL(ICATE), stat=allocate_status )
1558            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating CAPG_TBL in urban_param_init')
1559            ALLOCATE( AKSR_TBL(ICATE), stat=allocate_status )
1560            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating AKSR_TBL in urban_param_init')
1561            ALLOCATE( AKSB_TBL(ICATE), stat=allocate_status )
1562            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating AKSB_TBL in urban_param_init')
1563            ALLOCATE( AKSG_TBL(ICATE), stat=allocate_status )
1564            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating AKSG_TBL in urban_param_init')
1565            ALLOCATE( ALBR_TBL(ICATE), stat=allocate_status )
1566            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating ALBR_TBL in urban_param_init')
1567            ALLOCATE( ALBB_TBL(ICATE), stat=allocate_status )
1568            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating ALBB_TBL in urban_param_init')
1569            ALLOCATE( ALBG_TBL(ICATE), stat=allocate_status )
1570            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating ALBG_TBL in urban_param_init')
1571            ALLOCATE( EPSR_TBL(ICATE), stat=allocate_status )
1572            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating EPSR_TBL in urban_param_init')
1573            ALLOCATE( EPSB_TBL(ICATE), stat=allocate_status )
1574            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating EPSB_TBL in urban_param_init')
1575            ALLOCATE( EPSG_TBL(ICATE), stat=allocate_status )
1576            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating EPSG_TBL in urban_param_init')
1577            ALLOCATE( Z0R_TBL(ICATE), stat=allocate_status )
1578            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating Z0R_TBL in urban_param_init')
1579            ALLOCATE( Z0B_TBL(ICATE), stat=allocate_status )
1580            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating Z0B_TBL in urban_param_init')
1581            ALLOCATE( Z0G_TBL(ICATE), stat=allocate_status )
1582            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating Z0G_TBL in urban_param_init')
1583            ALLOCATE( AKANDA_URBAN_TBL(ICATE), stat=allocate_status )
1584            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating AKANDA_URBAN_TBL in urban_param_init')
1585            ALLOCATE( Z0HB_TBL(ICATE), stat=allocate_status )
1586            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating Z0HB_TBL in urban_param_init')
1587            ALLOCATE( Z0HG_TBL(ICATE), stat=allocate_status )
1588            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating Z0HG_TBL in urban_param_init')
1589            ALLOCATE( TRLEND_TBL(ICATE), stat=allocate_status )
1590            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TRLEND_TBL in urban_param_init')
1591            ALLOCATE( TBLEND_TBL(ICATE), stat=allocate_status )
1592            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TBLEND_TBL in urban_param_init')
1593            ALLOCATE( TGLEND_TBL(ICATE), stat=allocate_status )
1594            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TGLEND_TBL in urban_param_init')
1595            ALLOCATE( FRC_URB_TBL(ICATE), stat=allocate_status )
1596            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating FRC_URB_TBL in urban_param_init')
1597            ! ALLOCATE( ROOF_WIDTH(ICATE), stat=allocate_status )
1598            ! if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating ROOF_WIDTH in urban_param_init')
1599            ! ALLOCATE( ROAD_WIDTH(ICATE), stat=allocate_status )
1600            ! if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating ROAD_WIDTH in urban_param_init')
1601            !for BEP
1602            ALLOCATE( NUMDIR_TBL(ICATE), stat=allocate_status )
1603            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating NUMDIR_TBL in urban_param_init')
1604            ALLOCATE( STREET_DIRECTION_TBL(MAXDIRS , ICATE), stat=allocate_status )
1605            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating STREET_DIRECTION_TBL in urban_param_init')
1606            ALLOCATE( STREET_WIDTH_TBL(MAXDIRS , ICATE), stat=allocate_status )
1607            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating STREET_WIDTH_TBL in urban_param_init')
1608            ALLOCATE( BUILDING_WIDTH_TBL(MAXDIRS , ICATE), stat=allocate_status )
1609            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating BUILDING_WIDTH_TBL in urban_param_init')
1610            ALLOCATE( NUMHGT_TBL(ICATE), stat=allocate_status )
1611            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating NUMHGT_TBL in urban_param_init')
1612            ALLOCATE( HEIGHT_BIN_TBL(MAXHGTS , ICATE), stat=allocate_status )
1613            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating HEIGHT_BIN_TBL in urban_param_init')
1614            ALLOCATE( HPERCENT_BIN_TBL(MAXHGTS , ICATE), stat=allocate_status )
1615            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating HPERCENT_BIN_TBL in urban_param_init')
1616            ALLOCATE( COP_TBL(ICATE), stat=allocate_status )
1617            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating COP_TBL in urban_param_init')
1618            ALLOCATE( PWIN_TBL(ICATE), stat=allocate_status )
1619            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating PWIN_TBL in urban_param_init')
1620            ALLOCATE( BETA_TBL(ICATE), stat=allocate_status )
1621            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating BETA_TBL in urban_param_init')
1622            ALLOCATE( SW_COND_TBL(ICATE), stat=allocate_status )
1623            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating SW_COND_TBL in urban_param_init')
1624            ALLOCATE( TIME_ON_TBL(ICATE), stat=allocate_status )
1625            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TIME_ON_TBL in urban_param_init')
1626            ALLOCATE( TIME_OFF_TBL(ICATE), stat=allocate_status )
1627            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TIME_OFF_TBL in urban_param_init')
1628            ALLOCATE( TARGTEMP_TBL(ICATE), stat=allocate_status )
1629            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TARGTEMP_TBL in urban_param_init')
1630            ALLOCATE( GAPTEMP_TBL(ICATE), stat=allocate_status )
1631            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating GAPTEMP_TBL in urban_param_init')
1632            ALLOCATE( TARGHUM_TBL(ICATE), stat=allocate_status )
1633            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating TARGHUM_TBL in urban_param_init')
1634            ALLOCATE( GAPHUM_TBL(ICATE), stat=allocate_status )
1635            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating GAPHUM_TBL in urban_param_init')
1636            ALLOCATE( PERFLO_TBL(ICATE), stat=allocate_status )
1637            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating PERFLO_TBL in urban_param_init')
1638            ALLOCATE( HSESF_TBL(ICATE), stat=allocate_status )
1639            if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating HSESF_TBL in urban_param_init')
1640         endif
1641         numdir_tbl = 0
1642         street_direction_tbl = -1.E36
1643         street_width_tbl = 0
1644         building_width_tbl = 0
1645         numhgt_tbl = 0
1646         height_bin_tbl = -1.E36
1647         hpercent_bin_tbl = -1.E36
1648!end BEP
1649
1650      else if (name == "ZR") then
1651         read(string(indx+1:),*) zr_tbl(1:icate)
1652      else if (name == "SIGMA_ZED") then
1653         read(string(indx+1:),*) sigma_zed_tbl(1:icate)
1654      else if (name == "ROOF_WIDTH") then
1655         ALLOCATE( ROOF_WIDTH(ICATE), stat=allocate_status )
1656         if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating ROOF_WIDTH in urban_param_init')
1657
1658         read(string(indx+1:),*) roof_width(1:icate)
1659      else if (name == "ROAD_WIDTH") then
1660         ALLOCATE( ROAD_WIDTH(ICATE), stat=allocate_status )
1661         if(allocate_status /= 0) CALL wrf_error_fatal('Error allocating ROAD_WIDTH in urban_param_init')
1662         read(string(indx+1:),*) road_width(1:icate)
1663      else if (name == "AH") then
1664         read(string(indx+1:),*) ah_tbl(1:icate)
1665      else if (name == "FRC_URB") then
1666         read(string(indx+1:),*) frc_urb_tbl(1:icate)
1667      else if (name == "CAPR") then
1668         read(string(indx+1:),*) capr_tbl(1:icate)
1669         ! Convert CAPR_TBL from J m{-3} K{-1} to cal cm{-3} deg{-1}
1670         capr_tbl = capr_tbl * ( 1.0 / 4.1868 ) * 1.E-6
1671      else if (name == "CAPB") then
1672         read(string(indx+1:),*) capb_tbl(1:icate)
1673         ! Convert CABR_TBL from J m{-3} K{-1} to cal cm{-3} deg{-1}
1674         capb_tbl = capb_tbl * ( 1.0 / 4.1868 ) * 1.E-6
1675      else if (name == "CAPG") then
1676         read(string(indx+1:),*) capg_tbl(1:icate)
1677         ! Convert CABG_TBL from J m{-3} K{-1} to cal cm{-3} deg{-1}
1678         capg_tbl = capg_tbl * ( 1.0 / 4.1868 ) * 1.E-6
1679      else if (name == "AKSR") then
1680         read(string(indx+1:),*) aksr_tbl(1:icate)
1681         ! Convert AKSR_TBL from J m{-1} s{-1} K{-1} to cal cm{-1} s{-1} deg{-1}
1682         AKSR_TBL = AKSR_TBL * ( 1.0 / 4.1868 ) * 1.E-2
1683      else if (name == "AKSB") then
1684         read(string(indx+1:),*) aksb_tbl(1:icate)
1685         ! Convert AKSB_TBL from J m{-1} s{-1} K{-1} to cal cm{-1} s{-1} deg{-1}
1686         AKSB_TBL = AKSB_TBL * ( 1.0 / 4.1868 ) * 1.E-2
1687      else if (name == "AKSG") then
1688         read(string(indx+1:),*) aksg_tbl(1:icate)
1689         ! Convert AKSG_TBL from J m{-1} s{-1} K{-1} to cal cm{-1} s{-1} deg{-1}
1690         AKSG_TBL = AKSG_TBL * ( 1.0 / 4.1868 ) * 1.E-2
1691      else if (name == "ALBR") then
1692         read(string(indx+1:),*) albr_tbl(1:icate)
1693      else if (name == "ALBB") then
1694         read(string(indx+1:),*) albb_tbl(1:icate)
1695      else if (name == "ALBG") then
1696         read(string(indx+1:),*) albg_tbl(1:icate)
1697      else if (name == "EPSR") then
1698         read(string(indx+1:),*) epsr_tbl(1:icate)
1699      else if (name == "EPSB") then
1700         read(string(indx+1:),*) epsb_tbl(1:icate)
1701      else if (name == "EPSG") then
1702         read(string(indx+1:),*) epsg_tbl(1:icate)
1703      else if (name == "AKANDA_URBAN") then
1704         read(string(indx+1:),*) akanda_urban_tbl(1:icate)
1705      else if (name == "Z0B") then
1706         read(string(indx+1:),*) z0b_tbl(1:icate)
1707      else if (name == "Z0G") then
1708         read(string(indx+1:),*) z0g_tbl(1:icate)
1709      else if (name == "DDZR") then
1710         read(string(indx+1:),*) dzr(1:num_roof_layers)
1711         ! Convert thicknesses from m to cm
1712         dzr = dzr * 100.0
1713      else if (name == "DDZB") then
1714         read(string(indx+1:),*) dzb(1:num_wall_layers)
1715         ! Convert thicknesses from m to cm
1716         dzb = dzb * 100.0
1717      else if (name == "DDZG") then
1718         read(string(indx+1:),*) dzg(1:num_road_layers)
1719         ! Convert thicknesses from m to cm
1720         dzg = dzg * 100.0
1721      else if (name == "BOUNDR") then
1722         read(string(indx+1:),*) boundr_data
1723      else if (name == "BOUNDB") then
1724         read(string(indx+1:),*) boundb_data
1725      else if (name == "BOUNDG") then
1726         read(string(indx+1:),*) boundg_data
1727      else if (name == "TRLEND") then
1728         read(string(indx+1:),*) trlend_tbl(1:icate)
1729      else if (name == "TBLEND") then
1730         read(string(indx+1:),*) tblend_tbl(1:icate)
1731      else if (name == "TGLEND") then
1732         read(string(indx+1:),*) tglend_tbl(1:icate)
1733      else if (name == "CH_SCHEME") then
1734         read(string(indx+1:),*) ch_scheme_data
1735      else if (name == "TS_SCHEME") then
1736         read(string(indx+1:),*) ts_scheme_data
1737      else if (name == "AHOPTION") then
1738         read(string(indx+1:),*) ahoption
1739      else if (name == "AHDIUPRF") then
1740         read(string(indx+1:),*) ahdiuprf(1:24)
1741!for BEP
1742      else if (name == "STREET PARAMETERS") then
1743
1744         STREETLOOP : do
1745            read(11,'(A512)', iostat=iostatus) string
1746            if (string(1:1) == "#") cycle STREETLOOP
1747            if (trim(string) == "") cycle STREETLOOP
1748            if (string == "END STREET PARAMETERS") exit STREETLOOP
1749            read(string, *) k ! , dirst, ws, bs
1750            numdir_tbl(k) = numdir_tbl(k) + 1
1751            read(string, *) k, street_direction_tbl(numdir_tbl(k),k), &
1752                               street_width_tbl(numdir_tbl(k),k), &
1753                               building_width_tbl(numdir_tbl(k),k)
1754         enddo STREETLOOP
1755
1756      else if (name == "BUILDING HEIGHTS") then
1757
1758         read(string(indx+1:),*) k
1759         HEIGHTLOOP : do
1760            read(11,'(A512)', iostat=iostatus) string
1761            if (string(1:1) == "#") cycle HEIGHTLOOP
1762            if (trim(string) == "") cycle HEIGHTLOOP
1763            if (string == "END BUILDING HEIGHTS") exit HEIGHTLOOP
1764            read(string,*) dummy_hgt, dummy_pct
1765            numhgt_tbl(k) = numhgt_tbl(k) + 1
1766            height_bin_tbl(numhgt_tbl(k), k) = dummy_hgt
1767            hpercent_bin_tbl(numhgt_tbl(k),k) = dummy_pct
1768           
1769         enddo HEIGHTLOOP
1770         pctsum = sum ( hpercent_bin_tbl(:,k) , mask=(hpercent_bin_tbl(:,k)>-1.E25 ) )
1771         if ( pctsum /= 100.) then
1772            write (*,'(//,"Building height percentages for category ", I2, " must sum to 100.0")') k
1773            write (*,'("Currently, they sum to ", F6.2,/)') pctsum
1774            CALL wrf_error_fatal('pctsum is not equal to 100.')
1775         endif
1776      else if ( name == "Z0R") then
1777         read(string(indx+1:),*) Z0R_tbl(1:icate)
1778      else if ( name == "COP") then
1779         read(string(indx+1:),*) cop_tbl(1:icate)
1780      else if ( name == "PWIN") then
1781         read(string(indx+1:),*) pwin_tbl(1:icate)
1782      else if ( name == "BETA") then
1783         read(string(indx+1:),*) beta_tbl(1:icate)
1784      else if ( name == "SW_COND") then
1785         read(string(indx+1:),*) sw_cond_tbl(1:icate)
1786      else if ( name == "TIME_ON") then
1787         read(string(indx+1:),*) time_on_tbl(1:icate)
1788      else if ( name == "TIME_OFF") then
1789         read(string(indx+1:),*) time_off_tbl(1:icate)
1790      else if ( name == "TARGTEMP") then
1791         read(string(indx+1:),*) targtemp_tbl(1:icate)
1792      else if ( name == "GAPTEMP") then
1793         read(string(indx+1:),*) gaptemp_tbl(1:icate)
1794      else if ( name == "TARGHUM") then
1795         read(string(indx+1:),*) targhum_tbl(1:icate)
1796      else if ( name == "GAPHUM") then
1797         read(string(indx+1:),*) gaphum_tbl(1:icate)
1798      else if ( name == "PERFLO") then
1799         read(string(indx+1:),*) perflo_tbl(1:icate)
1800      else if (name == "HSEQUIP") then
1801         read(string(indx+1:),*) hsequip_tbl(1:24)
1802      else if (name == "HSEQUIP_SCALE_FACTOR") then
1803         read(string(indx+1:),*) hsesf_tbl(1:icate)
1804!end BEP         
1805      else
1806         CALL wrf_error_fatal('URBPARM.TBL: Unrecognized NAME = "'//trim(name)//'" in Subr URBAN_PARAM_INIT')
1807      endif
1808   enddo READLOOP
1809
1810   CLOSE(11)
1811
1812   ! Assign a few table values that do not need to come from URBPARM.TBL
1813
1814   Z0HB_TBL = 0.1 * Z0B_TBL
1815   Z0HG_TBL = 0.1 * Z0G_TBL
1816
1817   DO LC = 1, ICATE
1818
1819      ! HGT:  Normalized height
1820      HGT_TBL(LC) = ZR_TBL(LC) / ( ROAD_WIDTH(LC) + ROOF_WIDTH(LC) )
1821
1822      ! R:  Normalized Roof Width (a.k.a. "building coverage ratio")
1823      R_TBL(LC)  = ROOF_WIDTH(LC) / ( ROAD_WIDTH(LC) + ROOF_WIDTH(LC) )
1824
1825      RW_TBL(LC) = 1.0 - R_TBL(LC)
1826      BETR_TBL(LC) = 0.0
1827      BETB_TBL(LC) = 0.0
1828      BETG_TBL(LC) = 0.0
1829
1830      ! The following urban canyon geometry parameters are following Macdonald's (1998) formulations
1831     
1832      ! Lambda_P :: Plan areal fraction, which corresponds to R for a 2-d canyon.
1833      ! Lambda_F :: Frontal area index, which corresponds to HGT for a 2-d canyon
1834      ! Cd       :: Drag coefficient ( 1.2 from Grimmond and Oke, 1998 )
1835      ! Alpha_macd :: Emperical coefficient ( 4.43 from Macdonald et al., 1998 )
1836      ! Beta_macd  :: Correction factor for the drag coefficient ( 1.0 from Macdonald et al., 1998 )
1837
1838      Lambda_P = R_TBL(LC)
1839      Lambda_F = HGT_TBL(LC)
1840      Cd         = 1.2
1841      alpha_macd = 4.43
1842      beta_macd  = 1.0
1843
1844
1845      ZDC_TBL(LC) = ZR_TBL(LC) * ( 1.0 + ( alpha_macd ** ( -Lambda_P ) )  * ( Lambda_P - 1.0 ) )
1846
1847      Z0C_TBL(LC) = ZR_TBL(LC) * ( 1.0 - ZDC_TBL(LC)/ZR_TBL(LC) ) * &
1848           exp (-(0.5 * beta_macd * Cd / (VonK**2) * ( 1.0-ZDC_TBL(LC)/ZR_TBL(LC) ) * Lambda_F )**(-0.5))
1849
1850      IF (SF_URBAN_PHYSICS == 1) THEN
1851         ! Include roof height variability in Macdonald
1852         ! to parameterize Z0R as a function of ZR_SD (Standard Deviation)
1853         Lambda_FR  = SIGMA_ZED_TBL(LC) / ( ROAD_WIDTH(LC) + ROOF_WIDTH(LC) )
1854         Z0R_TBL(LC) = ZR_TBL(LC) * ( 1.0 - ZDC_TBL(LC)/ZR_TBL(LC) ) &
1855              * exp ( -(0.5 * beta_macd * Cd / (VonK**2) &
1856              * ( 1.0-ZDC_TBL(LC)/ZR_TBL(LC) ) * Lambda_FR )**(-0.5))
1857      ENDIF
1858
1859      !
1860      ! Z0HC still one-tenth of Z0C, as before ?
1861      !
1862
1863      Z0HC_TBL(LC) = 0.1 * Z0C_TBL(LC)
1864
1865      !
1866      ! Calculate Sky View Factor:
1867      !
1868      DHGT=HGT_TBL(LC)/100.
1869      HGT=0.
1870      VFWS=0.
1871      HGT=HGT_TBL(LC)-DHGT/2.
1872      do k=1,99
1873         HGT=HGT-DHGT
1874         VFWS=VFWS+0.25*(1.-HGT/SQRT(HGT**2.+RW_TBL(LC)**2.))
1875      end do
1876
1877     VFWS=VFWS/99.
1878     VFWS=VFWS*2.
1879
1880     VFGS=1.-2.*VFWS*HGT_TBL(LC)/RW_TBL(LC)
1881     SVF_TBL(LC)=VFGS
1882   END DO
1883
1884   deallocate(roof_width)
1885   deallocate(road_width)
1886
1887   END SUBROUTINE urban_param_init
1888!===========================================================================
1889!
1890! subroutine urban_var_init: initialization of urban state variables
1891!
1892!===========================================================================
1893   SUBROUTINE urban_var_init(ISURBAN, TSURFACE0_URB,TLAYER0_URB,TDEEP0_URB,IVGTYP,  & ! in
1894                             ims,ime,jms,jme,kms,kme,num_soil_layers,         & ! in
1895!                             num_roof_layers,num_wall_layers,num_road_layers, & ! in
1896                             restart,sf_urban_physics,                     & !in
1897                             XXXR_URB2D,XXXB_URB2D,XXXG_URB2D,XXXC_URB2D,  & ! inout
1898                             TR_URB2D,TB_URB2D,TG_URB2D,TC_URB2D,QC_URB2D, & ! inout
1899                             TRL_URB3D,TBL_URB3D,TGL_URB3D,                & ! inout
1900                             SH_URB2D,LH_URB2D,G_URB2D,RN_URB2D,           & ! inout
1901                             TS_URB2D,                                     & ! inout
1902                             num_urban_layers,                             & ! in
1903                             TRB_URB4D,TW1_URB4D,TW2_URB4D,TGB_URB4D,      & ! inout
1904                             TLEV_URB3D,QLEV_URB3D,                        & ! inout
1905                             TW1LEV_URB3D,TW2LEV_URB3D,                    & ! inout
1906                             TGLEV_URB3D,TFLEV_URB3D,                      & ! inout
1907                             SF_AC_URB3D,LF_AC_URB3D,CM_AC_URB3D,          & ! inout
1908                             SFVENT_URB3D,LFVENT_URB3D,                    & ! inout
1909                             SFWIN1_URB3D,SFWIN2_URB3D,                    & ! inout
1910                             SFW1_URB3D,SFW2_URB3D,SFR_URB3D,SFG_URB3D,    & ! inout
1911                             A_U_BEP,A_V_BEP,A_T_BEP,A_Q_BEP,              & ! inout multi-layer urban
1912                             A_E_BEP,B_U_BEP,B_V_BEP,                      & ! inout multi-layer urban
1913                             B_T_BEP,B_Q_BEP,B_E_BEP,DLG_BEP,              & ! inout multi-layer urban
1914                             DL_U_BEP,SF_BEP,VL_BEP,                       & ! inout multi-layer urban
1915                             FRC_URB2D, UTYPE_URB2D)                         ! inout
1916   IMPLICIT NONE
1917
1918   INTEGER, INTENT(IN) :: ISURBAN, sf_urban_physics
1919   INTEGER, INTENT(IN) :: ims,ime,jms,jme,kms,kme,num_soil_layers
1920   INTEGER, INTENT(IN) :: num_urban_layers                            !multi-layer urban
1921!   INTEGER, INTENT(IN) :: num_roof_layers, num_wall_layers, num_road_layers
1922
1923   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN)                    :: TSURFACE0_URB
1924   REAL, DIMENSION( ims:ime, 1:num_soil_layers, jms:jme ), INTENT(IN) :: TLAYER0_URB
1925   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN)                    :: TDEEP0_URB
1926   INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(IN)                 :: IVGTYP
1927   LOGICAL , INTENT(IN) :: restart
1928 
1929   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TR_URB2D
1930   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TB_URB2D
1931   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TG_URB2D
1932   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TC_URB2D
1933   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: QC_URB2D
1934   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXR_URB2D
1935   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXB_URB2D
1936   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXG_URB2D
1937   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXC_URB2D
1938
1939!   REAL, DIMENSION(ims:ime, 1:num_roof_layers, jms:jme), INTENT(INOUT) :: TRL_URB3D
1940!   REAL, DIMENSION(ims:ime, 1:num_wall_layers, jms:jme), INTENT(INOUT) :: TBL_URB3D
1941!   REAL, DIMENSION(ims:ime, 1:num_road_layers, jms:jme), INTENT(INOUT) :: TGL_URB3D
1942   REAL, DIMENSION(ims:ime, 1:num_soil_layers, jms:jme), INTENT(INOUT) :: TRL_URB3D
1943   REAL, DIMENSION(ims:ime, 1:num_soil_layers, jms:jme), INTENT(INOUT) :: TBL_URB3D
1944   REAL, DIMENSION(ims:ime, 1:num_soil_layers, jms:jme), INTENT(INOUT) :: TGL_URB3D
1945
1946   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: SH_URB2D
1947   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: LH_URB2D
1948   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: G_URB2D
1949   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: RN_URB2D
1950   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TS_URB2D
1951
1952! multi-layer UCM variables
1953   REAL, DIMENSION(ims:ime, 1:num_urban_layers, jms:jme), INTENT(INOUT) :: TRB_URB4D
1954   REAL, DIMENSION(ims:ime, 1:num_urban_layers, jms:jme), INTENT(INOUT) :: TW1_URB4D
1955   REAL, DIMENSION(ims:ime, 1:num_urban_layers, jms:jme), INTENT(INOUT) :: TW2_URB4D
1956   REAL, DIMENSION(ims:ime, 1:num_urban_layers, jms:jme), INTENT(INOUT) :: TGB_URB4D
1957   REAL, DIMENSION(ims:ime, 1:num_urban_layers, jms:jme), INTENT(INOUT) :: TLEV_URB3D
1958   REAL, DIMENSION(ims:ime, 1:num_urban_layers, jms:jme), INTENT(INOUT) :: QLEV_URB3D
1959   REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: TW1LEV_URB3D
1960   REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: TW2LEV_URB3D
1961   REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: TGLEV_URB3D
1962   REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: TFLEV_URB3D
1963   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: LF_AC_URB3D
1964   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: SF_AC_URB3D
1965   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CM_AC_URB3D
1966   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: SFVENT_URB3D
1967   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: LFVENT_URB3D
1968   REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: SFWIN1_URB3D
1969   REAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: SFWIN2_URB3D
1970   REAL, DIMENSION(ims:ime, 1:num_urban_layers, jms:jme), INTENT(INOUT) :: SFW1_URB3D
1971   REAL, DIMENSION(ims:ime, 1:num_urban_layers, jms:jme), INTENT(INOUT) :: SFW2_URB3D
1972   REAL, DIMENSION(ims:ime, 1:num_urban_layers, jms:jme), INTENT(INOUT) :: SFR_URB3D
1973   REAL, DIMENSION(ims:ime, 1:num_urban_layers, jms:jme), INTENT(INOUT) :: SFG_URB3D
1974   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: A_U_BEP
1975   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: A_V_BEP
1976   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: A_T_BEP
1977   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: A_Q_BEP
1978   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: A_E_BEP
1979   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: B_U_BEP
1980   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: B_V_BEP
1981   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: B_T_BEP
1982   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: B_Q_BEP
1983   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: B_E_BEP
1984   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: VL_BEP
1985   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: DLG_BEP
1986   REAL, DIMENSION(ims:ime, kms:kme,jms:jme),INTENT(INOUT) :: SF_BEP
1987   REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: DL_U_BEP
1988!
1989   REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: FRC_URB2D
1990   INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: UTYPE_URB2D
1991   INTEGER                                            :: UTYPE_URB
1992
1993   INTEGER :: I,J,K
1994
1995   DO I=ims,ime
1996    DO J=jms,jme
1997
1998!      XXXR_URB2D(I,J)=0.
1999!      XXXB_URB2D(I,J)=0.
2000!      XXXG_URB2D(I,J)=0.
2001!      XXXC_URB2D(I,J)=0.
2002
2003      SH_URB2D(I,J)=0.
2004      LH_URB2D(I,J)=0.
2005      G_URB2D(I,J)=0.
2006      RN_URB2D(I,J)=0.
2007     
2008!m
2009      FRC_URB2D(I,J)=0.
2010      UTYPE_URB2D(I,J)=0
2011
2012            IF( IVGTYP(I,J) == ISURBAN)  THEN
2013               UTYPE_URB2D(I,J) = 2  ! for default. high-intensity
2014               UTYPE_URB = UTYPE_URB2D(I,J)  ! for default. high-intensity
2015               FRC_URB2D(I,J) = FRC_URB_TBL(UTYPE_URB)
2016            ENDIF
2017            IF( IVGTYP(I,J) == 31) THEN
2018               UTYPE_URB2D(I,J) = 3  ! low-intensity residential
2019               UTYPE_URB = UTYPE_URB2D(I,J)  ! low-intensity residential
2020               FRC_URB2D(I,J) = FRC_URB_TBL(UTYPE_URB)
2021            ENDIF
2022            IF( IVGTYP(I,J) == 32) THEN
2023               UTYPE_URB2D(I,J) = 2  ! high-intensity
2024               UTYPE_URB = UTYPE_URB2D(I,J)  ! high-intensity
2025               FRC_URB2D(I,J) = FRC_URB_TBL(UTYPE_URB)
2026            ENDIF
2027            IF( IVGTYP(I,J) == 33) THEN
2028               UTYPE_URB2D(I,J) = 1  !  Commercial/Industrial/Transportation
2029               UTYPE_URB = UTYPE_URB2D(I,J)  !  Commercial/Industrial/Transportation
2030               FRC_URB2D(I,J) = FRC_URB_TBL(UTYPE_URB)
2031            ENDIF
2032
2033
2034      QC_URB2D(I,J)=0.01
2035
2036       IF (.not.restart) THEN
2037
2038      XXXR_URB2D(I,J)=0.
2039      XXXB_URB2D(I,J)=0.
2040      XXXG_URB2D(I,J)=0.
2041      XXXC_URB2D(I,J)=0.
2042
2043
2044      TC_URB2D(I,J)=TSURFACE0_URB(I,J)+0.
2045      TR_URB2D(I,J)=TSURFACE0_URB(I,J)+0.
2046      TB_URB2D(I,J)=TSURFACE0_URB(I,J)+0.
2047      TG_URB2D(I,J)=TSURFACE0_URB(I,J)+0.
2048!
2049      TS_URB2D(I,J)=TSURFACE0_URB(I,J)+0.
2050
2051!      DO K=1,num_roof_layers
2052!     DO K=1,num_soil_layers
2053!         TRL_URB3D(I,1,J)=TLAYER0_URB(I,1,J)+0.
2054!         TRL_URB3D(I,2,J)=TLAYER0_URB(I,2,J)+0.
2055!         TRL_URB3D(I,3,J)=TLAYER0_URB(I,3,J)+0.
2056!         TRL_URB3D(I,4,J)=TLAYER0_URB(I,4,J)+0.
2057
2058          TRL_URB3D(I,1,J)=TLAYER0_URB(I,1,J)+0.
2059          TRL_URB3D(I,2,J)=0.5*(TLAYER0_URB(I,1,J)+TLAYER0_URB(I,2,J))
2060          TRL_URB3D(I,3,J)=TLAYER0_URB(I,2,J)+0.
2061          TRL_URB3D(I,4,J)=TLAYER0_URB(I,2,J)+(TLAYER0_URB(I,3,J)-TLAYER0_URB(I,2,J))*0.29
2062!     END DO
2063
2064!      DO K=1,num_wall_layers
2065!      DO K=1,num_soil_layers
2066!m        TBL_URB3D(I,1,J)=TLAYER0_URB(I,1,J)+0.
2067!m        TBL_URB3D(I,2,J)=TLAYER0_URB(I,2,J)+0.
2068!m        TBL_URB3D(I,3,J)=TLAYER0_URB(I,3,J)+0.
2069!m        TBL_URB3D(I,4,J)=TLAYER0_URB(I,4,J)+0.
2070
2071        TBL_URB3D(I,1,J)=TLAYER0_URB(I,1,J)+0.
2072        TBL_URB3D(I,2,J)=0.5*(TLAYER0_URB(I,1,J)+TLAYER0_URB(I,2,J))
2073        TBL_URB3D(I,3,J)=TLAYER0_URB(I,2,J)+0.
2074        TBL_URB3D(I,4,J)=TLAYER0_URB(I,2,J)+(TLAYER0_URB(I,3,J)-TLAYER0_URB(I,2,J))*0.29
2075
2076!      END DO
2077
2078!      DO K=1,num_road_layers
2079      DO K=1,num_soil_layers
2080        TGL_URB3D(I,K,J)=TLAYER0_URB(I,K,J)+0.
2081      END DO
2082     
2083! multi-layer urban
2084!      IF( sf_urban_physics .EQ. 2)THEN
2085       IF((SF_URBAN_PHYSICS.eq.2).OR.(SF_URBAN_PHYSICS.eq.3)) THEN
2086      DO k=1,num_urban_layers
2087!        TRB_URB4D(I,k,J)=TSURFACE0_URB(I,J)
2088!        TW1_URB4D(I,k,J)=TSURFACE0_URB(I,J)
2089!        TW2_URB4D(I,k,J)=TSURFACE0_URB(I,J)
2090!        TGB_URB4D(I,k,J)=TSURFACE0_URB(I,J)
2091!MT        TRB_URB4D(I,K,J)=tlayer0_urb(I,1,J)
2092!MT        TW1_URB4D(I,K,J)=tlayer0_urb(I,1,J)
2093!MT        TW2_URB4D(I,K,J)=tlayer0_urb(I,1,J)
2094        IF (UTYPE_URB2D(I,J) > 0) THEN
2095           TRB_URB4D(I,K,J)=TBLEND_TBL(UTYPE_URB2D(I,J))
2096           TW1_URB4D(I,K,J)=TBLEND_TBL(UTYPE_URB2D(I,J))
2097           TW2_URB4D(I,K,J)=TBLEND_TBL(UTYPE_URB2D(I,J))
2098        ELSE
2099           TRB_URB4D(I,K,J)=tlayer0_urb(I,1,J)
2100           TW1_URB4D(I,K,J)=tlayer0_urb(I,1,J)
2101           TW2_URB4D(I,K,J)=tlayer0_urb(I,1,J)
2102        ENDIF
2103        TGB_URB4D(I,K,J)=tlayer0_urb(I,1,J)
2104        SFW1_URB3D(I,K,J)=0.
2105        SFW2_URB3D(I,K,J)=0.
2106        SFR_URB3D(I,K,J)=0.
2107        SFG_URB3D(I,K,J)=0.
2108      ENDDO
2109       
2110       ENDIF
2111             
2112      if (SF_URBAN_PHYSICS.EQ.3) then
2113         LF_AC_URB3D(I,J)=0.
2114         SF_AC_URB3D(I,J)=0.
2115         CM_AC_URB3D(I,J)=0.
2116         SFVENT_URB3D(I,J)=0.
2117         LFVENT_URB3D(I,J)=0.
2118
2119         DO K=1,num_urban_layers
2120            TLEV_URB3D(I,K,J)=tlayer0_urb(I,1,J)
2121            TW1LEV_URB3D(I,K,J)=tlayer0_urb(I,1,J)
2122            TW2LEV_URB3D(I,K,J)=tlayer0_urb(I,1,J)
2123            TGLEV_URB3D(I,K,J)=tlayer0_urb(I,1,J)
2124            TFLEV_URB3D(I,K,J)=tlayer0_urb(I,1,J)
2125            QLEV_URB3D(I,K,J)=0.01
2126            SFWIN1_URB3D(I,K,J)=0.
2127            SFWIN2_URB3D(I,K,J)=0.
2128!rm         LF_AC_URB3D(I,J)=0.
2129!rm         SF_AC_URB3D(I,J)=0.
2130!rm         CM_AC_URB3D(I,J)=0.
2131!rm         SFVENT_URB3D(I,J)=0.
2132!rm         LFVENT_URB3D(I,J)=0.
2133         ENDDO
2134
2135      endif
2136
2137!      IF( sf_urban_physics .EQ. 2 )THEN
2138      IF((SF_URBAN_PHYSICS.eq.2).OR.(SF_URBAN_PHYSICS.eq.3)) THEN
2139      DO K= KMS,KME
2140          SF_BEP(I,K,J)=1.
2141          VL_BEP(I,K,J)=1.
2142          A_U_BEP(I,K,J)=0.
2143          A_V_BEP(I,K,J)=0.
2144          A_T_BEP(I,K,J)=0.
2145          A_E_BEP(I,K,J)=0.
2146          A_Q_BEP(I,K,J)=0.
2147          B_U_BEP(I,K,J)=0.
2148          B_V_BEP(I,K,J)=0.
2149          B_T_BEP(I,K,J)=0.
2150          B_E_BEP(I,K,J)=0.
2151          B_Q_BEP(I,K,J)=0.
2152          DLG_BEP(I,K,J)=0.
2153          DL_U_BEP(I,K,J)=0.
2154      END DO
2155      ENDIF       !sf_urban_physics=2
2156     ENDIF        !restart
2157    END DO
2158   END DO
2159   RETURN
2160   END SUBROUTINE urban_var_init
2161!===========================================================================
2162!
2163!  force_restore
2164!
2165!===========================================================================
2166   SUBROUTINE force_restore(CAP,AKS,DELT,S,R,H,LE,TSLEND,TSP,TS)
2167
2168     REAL, INTENT(IN)  :: CAP,AKS,DELT,S,R,H,LE,TSLEND,TSP
2169     REAL, INTENT(OUT) :: TS
2170     REAL              :: C1,C2
2171
2172     C2=24.*3600./2./3.14159
2173     C1=SQRT(0.5*C2*CAP*AKS)
2174
2175     TS = TSP + DELT*( (S+R-H-LE)/C1 -(TSP-TSLEND)/C2 )
2176
2177   END SUBROUTINE force_restore
2178!===========================================================================
2179!
2180!  bisection (not used)
2181!
2182!==============================================================================
2183   SUBROUTINE bisection(TSP,PS,S,EPS,RX,SIG,RHO,CP,CH,UA,QA,TA,EL,BET,AKS,TSL,DZ,TS)
2184
2185     REAL, INTENT(IN) :: TSP,PS,S,EPS,RX,SIG,RHO,CP,CH,UA,QA,TA,EL,BET,AKS,TSL,DZ
2186     REAL, INTENT(OUT) :: TS
2187     REAL :: ES,QS0,R,H,ELE,G0,F1,F
2188
2189     TS1 = TSP - 5.
2190     TS2 = TSP + 5.
2191
2192     DO ITERATION = 1,22
2193
2194       ES=6.11*EXP( (2.5*10.**6./461.51)*(TS1-273.15)/(273.15*TS1) )
2195       QS0=0.622*ES/(PS-0.378*ES)
2196       R=EPS*(RX-SIG*(TS1**4.)/60.)
2197       H=RHO*CP*CH*UA*(TS1-TA)*100.
2198       ELE=RHO*EL*CH*UA*BET*(QS0-QA)*100.
2199       G0=AKS*(TS1-TSL)/(DZ/2.)
2200       F1= S + R - H - ELE - G0
2201
2202       TS=0.5*(TS1+TS2)
2203
2204       ES=6.11*EXP( (2.5*10.**6./461.51)*(TS-273.15)/(273.15*TS) )
2205       QS0=0.622*ES/(PS-0.378*ES)
2206       R=EPS*(RX-SIG*(TS**4.)/60.)
2207       H=RHO*CP*CH*UA*(TS-TA)*100.
2208       ELE=RHO*EL*CH*UA*BET*(QS0-QA)*100.
2209       G0=AKS*(TS-TSL)/(DZ/2.)
2210       F = S + R - H - ELE - G0
2211
2212       IF (F1*F > 0.0) THEN
2213         TS1=TS
2214        ELSE
2215         TS2=TS
2216       END IF
2217
2218     END DO
2219
2220     RETURN
2221END SUBROUTINE bisection
2222!===========================================================================
2223
2224SUBROUTINE SFCDIF_URB (ZLM,Z0,THZ0,THLM,SFCSPD,AKANDA,AKMS,AKHS,RLMO,CD)
2225
2226! ----------------------------------------------------------------------
2227! SUBROUTINE SFCDIF_URB (Urban version of SFCDIF_off)
2228! ----------------------------------------------------------------------
2229! CALCULATE SURFACE LAYER EXCHANGE COEFFICIENTS VIA ITERATIVE PROCESS.
2230! SEE CHEN ET AL (1997, BLM)
2231! ----------------------------------------------------------------------
2232
2233      IMPLICIT NONE
2234      REAL     WWST, WWST2, G, VKRM, EXCM, BETA, BTG, ELFC, WOLD, WNEW
2235      REAL     PIHF, EPSU2, EPSUST, EPSIT, EPSA, ZTMIN, ZTMAX, HPBL,     &
2236     & SQVISC
2237      REAL     RIC, RRIC, FHNEU, RFC,RLMO_THR, RFAC, ZZ, PSLMU, PSLMS, PSLHU,     &
2238     & PSLHS
2239      REAL     XX, PSPMU, YY, PSPMS, PSPHU, PSPHS, ZLM, Z0, THZ0, THLM
2240      REAL     SFCSPD, AKANDA, AKMS, AKHS, ZU, ZT, RDZ, CXCH
2241      REAL     DTHV, DU2, BTGH, WSTAR2, USTAR, ZSLU, ZSLT, RLOGU, RLOGT
2242      REAL     RLMO, ZETALT, ZETALU, ZETAU, ZETAT, XLU4, XLT4, XU4, XT4
2243!CC   ......REAL ZTFC
2244
2245      REAL     XLU, XLT, XU, XT, PSMZ, SIMM, PSHZ, SIMH, USTARK, RLMN,  &
2246     &         RLMA
2247
2248      INTEGER  ITRMX, ILECH, ITR
2249      REAL,    INTENT(OUT) :: CD
2250      PARAMETER                                                         &
2251     &        (WWST = 1.2,WWST2 = WWST * WWST,G = 9.8,VKRM = 0.40,      &
2252     &         EXCM = 0.001                                             &
2253     &        ,BETA = 1./270.,BTG = BETA * G,ELFC = VKRM * BTG          &
2254     &                  ,WOLD =.15,WNEW = 1. - WOLD,ITRMX = 05,         &
2255     &                   PIHF = 3.14159265/2.)
2256      PARAMETER                                                         &
2257     &         (EPSU2 = 1.E-4,EPSUST = 0.07,EPSIT = 1.E-4,EPSA = 1.E-8  &
2258     &         ,ZTMIN = -5.,ZTMAX = 1.,HPBL = 1000.0                    &
2259     &          ,SQVISC = 258.2)
2260      PARAMETER                                                         &
2261     &       (RIC = 0.183,RRIC = 1.0/ RIC,FHNEU = 0.8,RFC = 0.191       &
2262     &        ,RLMO_THR = 0.001,RFAC = RIC / (FHNEU * RFC * RFC))
2263
2264! ----------------------------------------------------------------------
2265! NOTE: THE TWO CODE BLOCKS BELOW DEFINE FUNCTIONS
2266! ----------------------------------------------------------------------
2267! LECH'S SURFACE FUNCTIONS
2268! ----------------------------------------------------------------------
2269      PSLMU (ZZ)= -0.96* log (1.0-4.5* ZZ)
2270      PSLMS (ZZ)= ZZ * RRIC -2.076* (1. -1./ (ZZ +1.))
2271      PSLHU (ZZ)= -0.96* log (1.0-4.5* ZZ)
2272
2273! ----------------------------------------------------------------------
2274! PAULSON'S SURFACE FUNCTIONS
2275! ----------------------------------------------------------------------
2276      PSLHS (ZZ)= ZZ * RFAC -2.076* (1. -1./ (ZZ +1.))
2277      PSPMU (XX)= -2.* log ( (XX +1.)*0.5) - log ( (XX * XX +1.)*0.5)   &
2278     &        +2.* ATAN (XX)                                            &
2279     &- PIHF
2280      PSPMS (YY)= 5.* YY
2281      PSPHU (XX)= -2.* log ( (XX * XX +1.)*0.5)
2282
2283! ----------------------------------------------------------------------
2284! THIS ROUTINE SFCDIF CAN HANDLE BOTH OVER OPEN WATER (SEA, OCEAN) AND
2285! OVER SOLID SURFACE (LAND, SEA-ICE).
2286! ----------------------------------------------------------------------
2287      PSPHS (YY)= 5.* YY
2288
2289! ----------------------------------------------------------------------
2290!     ZTFC: RATIO OF ZOH/ZOM  LESS OR EQUAL THAN 1
2291!     C......ZTFC=0.1
2292!     CZIL: CONSTANT C IN Zilitinkevich, S. S.1995,:NOTE ABOUT ZT
2293! ----------------------------------------------------------------------
2294      ILECH = 0
2295
2296! ----------------------------------------------------------------------
2297!      ZILFC = - CZIL * VKRM * SQVISC
2298!     C.......ZT=Z0*ZTFC
2299      ZU = Z0
2300      RDZ = 1./ ZLM
2301      CXCH = EXCM * RDZ
2302      DTHV = THLM - THZ0
2303
2304! ----------------------------------------------------------------------
2305! BELJARS CORRECTION OF USTAR
2306! ----------------------------------------------------------------------
2307      DU2 = MAX (SFCSPD * SFCSPD,EPSU2)
2308!cc   If statements to avoid TANGENT LINEAR problems near zero
2309      BTGH = BTG * HPBL
2310      IF (BTGH * AKHS * DTHV .ne. 0.0) THEN
2311         WSTAR2 = WWST2* ABS (BTGH * AKHS * DTHV)** (2./3.)
2312      ELSE
2313         WSTAR2 = 0.0
2314      END IF
2315
2316! ----------------------------------------------------------------------
2317! ZILITINKEVITCH APPROACH FOR ZT
2318! ----------------------------------------------------------------------
2319      USTAR = MAX (SQRT (AKMS * SQRT (DU2+ WSTAR2)),EPSUST)
2320
2321! ----------------------------------------------------------------------
2322! KCL/TL Try Kanda approach instead (Kanda et al. 2007, JAMC)
2323!      ZT = EXP (ZILFC * SQRT (USTAR * Z0))* Z0
2324      ZT = EXP (2.0-AKANDA*(SQVISC**2 * USTAR * Z0)**0.25)* Z0
2325     
2326      ZSLU = ZLM + ZU
2327
2328      ZSLT = ZLM + ZT
2329      RLOGU = log (ZSLU / ZU)
2330
2331      RLOGT = log (ZSLT / ZT)
2332
2333      RLMO = ELFC * AKHS * DTHV / USTAR **3
2334! ----------------------------------------------------------------------
2335! 1./MONIN-OBUKKHOV LENGTH-SCALE
2336! ----------------------------------------------------------------------
2337      DO ITR = 1,ITRMX
2338         ZETALT = MAX (ZSLT * RLMO,ZTMIN)
2339         RLMO = ZETALT / ZSLT
2340         ZETALU = ZSLU * RLMO
2341         ZETAU = ZU * RLMO
2342
2343         ZETAT = ZT * RLMO
2344         IF (ILECH .eq. 0) THEN
2345            IF (RLMO .lt. 0.0)THEN
2346               XLU4 = 1. -16.* ZETALU
2347               XLT4 = 1. -16.* ZETALT
2348               XU4 = 1. -16.* ZETAU
2349
2350               XT4 = 1. -16.* ZETAT
2351               XLU = SQRT (SQRT (XLU4))
2352               XLT = SQRT (SQRT (XLT4))
2353               XU = SQRT (SQRT (XU4))
2354
2355               XT = SQRT (SQRT (XT4))
2356
2357               PSMZ = PSPMU (XU)
2358               SIMM = PSPMU (XLU) - PSMZ + RLOGU
2359               PSHZ = PSPHU (XT)
2360               SIMH = PSPHU (XLT) - PSHZ + RLOGT
2361            ELSE
2362               ZETALU = MIN (ZETALU,ZTMAX)
2363               ZETALT = MIN (ZETALT,ZTMAX)
2364               PSMZ = PSPMS (ZETAU)
2365               SIMM = PSPMS (ZETALU) - PSMZ + RLOGU
2366               PSHZ = PSPHS (ZETAT)
2367               SIMH = PSPHS (ZETALT) - PSHZ + RLOGT
2368            END IF
2369! ----------------------------------------------------------------------
2370! LECH'S FUNCTIONS
2371! ----------------------------------------------------------------------
2372         ELSE
2373            IF (RLMO .lt. 0.)THEN
2374               PSMZ = PSLMU (ZETAU)
2375               SIMM = PSLMU (ZETALU) - PSMZ + RLOGU
2376               PSHZ = PSLHU (ZETAT)
2377               SIMH = PSLHU (ZETALT) - PSHZ + RLOGT
2378            ELSE
2379               ZETALU = MIN (ZETALU,ZTMAX)
2380               ZETALT = MIN (ZETALT,ZTMAX)
2381               PSMZ = PSLMS (ZETAU)
2382               SIMM = PSLMS (ZETALU) - PSMZ + RLOGU
2383               PSHZ = PSLHS (ZETAT)
2384               SIMH = PSLHS (ZETALT) - PSHZ + RLOGT
2385            END IF
2386! ----------------------------------------------------------------------
2387! BELJAARS CORRECTION FOR USTAR
2388! ----------------------------------------------------------------------
2389         END IF
2390            USTAR = MAX (SQRT (AKMS * SQRT (DU2+ WSTAR2)),EPSUST)
2391            !KCL/TL
2392            !ZT = EXP (ZILFC * SQRT (USTAR * Z0))* Z0
2393            ZT = EXP (2.0-AKANDA*(SQVISC**2 * USTAR * Z0)**0.25)* Z0
2394            ZSLT = ZLM + ZT
2395            RLOGT = log (ZSLT / ZT)
2396            USTARK = USTAR * VKRM
2397            AKMS = MAX (USTARK / SIMM,CXCH)
2398            AKHS = MAX (USTARK / SIMH,CXCH)
2399!
2400         IF (BTGH * AKHS * DTHV .ne. 0.0) THEN
2401            WSTAR2 = WWST2* ABS (BTGH * AKHS * DTHV)** (2./3.)
2402         ELSE
2403            WSTAR2 = 0.0
2404         END IF
2405!-----------------------------------------------------------------------
2406         RLMN = ELFC * AKHS * DTHV / USTAR **3
2407!-----------------------------------------------------------------------
2408!     IF(ABS((RLMN-RLMO)/RLMA).LT.EPSIT)    GO TO 110
2409!-----------------------------------------------------------------------
2410         RLMA = RLMO * WOLD+ RLMN * WNEW
2411!-----------------------------------------------------------------------
2412         RLMO = RLMA
2413
2414      END DO
2415
2416      CD = USTAR*USTAR/SFCSPD**2
2417! ----------------------------------------------------------------------
2418  END SUBROUTINE SFCDIF_URB
2419! ----------------------------------------------------------------------
2420!===========================================================================
2421END MODULE module_sf_urban
Note: See TracBrowser for help on using the repository browser.