| blob | 84e53ec930854ddf368c56c8fbcd86784e4adf57 |
1 MODULE module_sf_urban
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 !===============================================================================
9 CHARACTER(LEN=4) :: LU_DATA_TYPE
11 INTEGER :: ICATE
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
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
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
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
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
70 INTEGER :: allocate_status
72 ! INTEGER :: num_roof_layers
73 ! INTEGER :: num_wall_layers
74 ! INTEGER :: num_road_layers
76 CONTAINS
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 !===============================================================================
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 )
280 IMPLICIT NONE
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 [-]
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
296 !-------------------------------------------------------------------------------
297 ! C: configuration variables
298 !-------------------------------------------------------------------------------
300 LOGICAL, INTENT(IN) :: LSOLAR ! logical [true=both, false=SSG only]
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.
306 INTEGER, INTENT(IN) :: num_roof_layers
307 INTEGER, INTENT(IN) :: num_wall_layers
308 INTEGER, INTENT(IN) :: num_road_layers
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]
315 !-------------------------------------------------------------------------------
316 ! I: input variables from LSM to Urban
317 !-------------------------------------------------------------------------------
319 INTEGER, INTENT(IN) :: UTYPE ! urban type [1=Commercial/Industrial, 2=High-intensity residential,
320 ! 3=low-intensity residential]
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]
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]
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
348 !-------------------------------------------------------------------------------
349 ! O: output variables from Urban to LSM
350 !-------------------------------------------------------------------------------
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]
373 !-------------------------------------------------------------------------------
374 ! H: Historical (state) variables of Urban : LSM <--> Urban
375 !-------------------------------------------------------------------------------
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)
391 REAL, INTENT(INOUT):: TR, TB, TG, TC, QC, UC
392 REAL, INTENT(INOUT):: XXXR, XXXB, XXXG, XXXC
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
398 !-------------------------------------------------------------------------------
399 ! L: Local variables from read_param
400 !-------------------------------------------------------------------------------
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
411 REAL :: TH2X !m
413 INTEGER :: BOUNDR, BOUNDB, BOUNDG
414 INTEGER :: CH_SCHEME, TS_SCHEME
416 LOGICAL :: SHADOW ! [true=consider svf and shadow effects, false=consider svf effect only]
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 !-------------------------------------------------------------------------------
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
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
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]
471 REAL :: XXX, X, Z0, Z0H, CD, CH
472 REAL :: XXX2, PSIM2, PSIH2, XXX10, PSIM10, PSIH10
473 REAL :: PSIX, PSIT, PSIX2, PSIT2, PSIX10, PSIT10
475 REAL :: TRP, TBP, TGP, TCP, QCP, TST, QST
477 INTEGER :: iteration, K
478 INTEGER :: tloc
480 !-------------------------------------------------------------------------------
481 ! Set parameters
482 !-------------------------------------------------------------------------------
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
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)
502 ! Miao, 2007/01/17, cal. ah
503 if(ahoption==1) AH=AH*ahdiuprf(tloc)
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
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
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
524 !-------------------------------------------------------------------------------
525 ! Convert unit from MKS to cgs
526 ! Renew surface and layer temperatures
527 !-------------------------------------------------------------------------------
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
535 TRP=TR
536 TBP=TB
537 TGP=TG
538 TCP=TC
539 QCP=QC
541 TAV=TA*(1.+0.61*QA)
542 PS=RHOO*287.*TAV/100. ![hPa]
544 !-------------------------------------------------------------------------------
545 ! Canopy wind
546 !-------------------------------------------------------------------------------
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
561 !-------------------------------------------------------------------------------
562 ! Net Short Wave Radiation at roof, wall, and road
563 !-------------------------------------------------------------------------------
565 SHADOW = .false.
566 ! SHADOW = .true.
568 IF (SSG > 0.0) THEN
570 IF(.NOT.SHADOW) THEN ! no shadow effects model
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)
578 ELSE ! shadow effects model
580 FAI=XLAT*PI/180.
582 THEATAS=ABS(ASIN(COSZ))
583 THEATAZ=ABS(ACOS(COSZ))
585 SNT=COS(DECLIN)*SIN(OMG)/COS(THEATAS)
586 CNT=(COSZ*SIN(FAI)-SIN(DECLIN))/COS(THEATAS)/COS(FAI)
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.))
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)
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
615 SLX=(SLX1+SLX2+SLX3+SLX4+SLX5+SLX6+SLX7+SLX8)/8.
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)
623 END IF
625 SR=SR1
626 SG=SG1+SG2
627 SB=SB1+SB2
629 SNET=R*SR+W*SB+RW*SG
631 ELSE
633 SR=0.
634 SG=0.
635 SB=0.
636 SNET=0.
638 END IF
640 !-------------------------------------------------------------------------------
641 ! Roof
642 !-------------------------------------------------------------------------------
644 !-------------------------------------------------------------------------------
645 ! CHR, CDR, BETR
646 !-------------------------------------------------------------------------------
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)
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)
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
664 IF(RAIN > 1.) BETR=0.7
666 IF (TS_SCHEME == 1) THEN
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
683 ! TRP=350.
685 DO ITERATION=1,20
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.)
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.)
697 F = SR + RR - HR - ELER - G0R
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)
704 DFDT = DRRDTR - DHRDTR - DELERDTR - DG0RDTR
705 DTR = F/DFDT
707 TR = TRP - DTR
708 TRP = TR
710 IF( ABS(F) < 0.000001 .AND. ABS(DTR) < 0.000001 ) EXIT
712 END DO
714 ! multi-layer heat equation model
716 CALL multi_layer(num_roof_layers,BOUNDR,G0R,CAPR,AKSR,TRL,DZR,DELT,TRLEND)
718 ELSE
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)
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
728 CALL force_restore(CAPR,AKSR,DELT,SR,RR,HR,ELER,TRLEND,TRP,TR)
730 TRP=TR
732 END IF
734 FLXTHR=HR/RHO/CP/100.
735 FLXHUMR=ELER/RHO/EL/100.
737 !-------------------------------------------------------------------------------
738 ! Wall and Road
739 !-------------------------------------------------------------------------------
741 !-------------------------------------------------------------------------------
742 ! CHC, CHB, CDB, BETB, CHG, CDG, BETG
743 !-------------------------------------------------------------------------------
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
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
757 IF (CH_SCHEME == 1) THEN
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)
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)
768 ELSE
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.
775 END IF
777 CHC=ALPHAC/RHO/CP/UA
778 CHB=ALPHAB/RHO/CP/UC
779 CHG=ALPHAG/RHO/CP/UC
781 BETB=0.0
782 IF(RAIN > 1.) BETG=0.7
784 IF (TS_SCHEME == 1) THEN
786 !-------------------------------------------------------------------------------
787 ! TB, TG Solving Non-Linear Simultaneous Equation by Newton-Rapson
788 ! TBL,TGL Solving Heat Equation by Tri Diagonal Matrix Algorithm
789 !-------------------------------------------------------------------------------
791 ! TBP=350.
792 ! TGP=350.
794 DO ITERATION=1,20
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.)
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.)
806 RG1=EPSG*( RX*VFGS &
807 +EPSB*VFGW*SIG*TBP**4./60. &
808 -SIG*TGP**4./60. )
810 RB1=EPSB*( RX*VFWS &
811 +EPSG*VFWG*SIG*TGP**4./60. &
812 +EPSB*VFWW*SIG*TBP**4./60. &
813 -SIG*TBP**4./60. )
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. )
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. )
825 RG=RG1+RG2
826 RB=RB1+RB2
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.
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.
839 DRBDTB=DRBDTB1+DRBDTB2
840 DRBDTG=DRBDTG1+DRBDTG2
841 DRGDTB=DRGDTB1+DRGDTB2
842 DRGDTG=DRGDTG1+DRGDTG2
844 HB=RHO*CP*CHB*UC*(TBP-TCP)*100.
845 HG=RHO*CP*CHG*UC*(TGP-TCP)*100.
847 DTCDTB=W*ALPHAB/(RW*ALPHAC+RW*ALPHAG+W*ALPHAB)
848 DTCDTG=RW*ALPHAG/(RW*ALPHAC+RW*ALPHAG+W*ALPHAB)
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.
855 ELEB=RHO*EL*CHB*UC*BETB*(QS0B-QCP)*100.
856 ELEG=RHO*EL*CHG*UC*BETG*(QS0G-QCP)*100.
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)
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.
866 G0B=AKSB*(TBP-TBL(1))/(DZB(1)/2.)
867 G0G=AKSG*(TGP-TGL(1))/(DZG(1)/2.)
869 DG0BDTB=2.*AKSB/DZB(1)
870 DG0BDTG=0.
871 DG0GDTG=2.*AKSG/DZG(1)
872 DG0GDTB=0.
874 F = SB + RB - HB - ELEB - G0B
875 FX = DRBDTB - DHBDTB - DELEBDTB - DG0BDTB
876 FY = DRBDTG - DHBDTG - DELEBDTG - DG0BDTG
878 GF = SG + RG - HG - ELEG - G0G
879 GX = DRGDTB - DHGDTB - DELEGDTB - DG0GDTB
880 GY = DRGDTG - DHGDTG - DELEGDTG - DG0GDTG
882 DTB = (GF*FY-F*GY)/(FX*GY-GX*FY)
883 DTG = -(GF+GX*DTB)/GY
885 TB = TBP + DTB
886 TG = TGP + DTG
888 TBP = TB
889 TGP = TG
891 TC1=RW*ALPHAC+RW*ALPHAG+W*ALPHAB
892 TC2=RW*ALPHAC*TA+RW*ALPHAG*TGP+W*ALPHAB*TBP
893 TC=TC2/TC1
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
899 DTC=TCP - TC
900 TCP=TC
901 QCP=QC
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
907 END DO
909 CALL multi_layer(num_wall_layers,BOUNDB,G0B,CAPB,AKSB,TBL,DZB,DELT,TBLEND)
911 CALL multi_layer(num_road_layers,BOUNDG,G0G,CAPG,AKSG,TGL,DZG,DELT,TGLEND)
913 ELSE
915 !-------------------------------------------------------------------------------
916 ! TB, TG by Force-Restore Method
917 !-------------------------------------------------------------------------------
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)
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)
925 RG1=EPSG*( RX*VFGS &
926 +EPSB*VFGW*SIG*TBP**4./60. &
927 -SIG*TGP**4./60. )
929 RB1=EPSB*( RX*VFWS &
930 +EPSG*VFWG*SIG*TGP**4./60. &
931 +EPSB*VFWW*SIG*TBP**4./60. &
932 -SIG*TBP**4./60. )
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. )
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. )
944 RG=RG1+RG2
945 RB=RB1+RB2
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
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
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)
958 TBP=TB
959 TGP=TG
961 TC1=RW*ALPHAC+RW*ALPHAG+W*ALPHAB
962 TC2=RW*ALPHAC*TA+RW*ALPHAG*TGP+W*ALPHAB*TBP
963 TC=TC2/TC1
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
969 TCP=TC
970 QCP=QC
972 END IF
975 FLXTHB=HB/RHO/CP/100.
976 FLXHUMB=ELEB/RHO/EL/100.
977 FLXTHG=HG/RHO/CP/100.
978 FLXHUMG=ELEG/RHO/EL/100.
980 !-------------------------------------------------------------------------------
981 ! Total Fluxes from Urban Canopy
982 !-------------------------------------------------------------------------------
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
995 !----------------------------------------------------------------------------
996 ! Convert Unit: FLUXES and u* T* q* --> WRF
997 !----------------------------------------------------------------------------
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]
1009 UST = SQRT(FLXUV) ! u* [m/s]
1010 TST = -FLXTH/UST ! T* [K]
1011 QST = -FLXHUM/UST ! q* [-]
1013 !------------------------------------------------------
1014 ! diagnostic GRID AVERAGED PSIM PSIH TS QS --> WRF
1015 !------------------------------------------------------
1017 Z0 = Z0C
1018 Z0H = Z0HC
1019 Z = ZA - ZDC
1021 XXX = 0.4*9.81*Z*TST/TA/UST/UST
1023 IF ( XXX >= 1. ) XXX = 1.
1024 IF ( XXX <= -5. ) XXX = -5.
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
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
1046 !-------------------------------------------------------
1047 ! diagnostic GRID AVERAGED U10 V10 TH2 Q2 --> WRF
1048 !-------------------------------------------------------
1050 XXX2 = (2./Z)*XXX
1051 IF ( XXX2 >= 1. ) XXX2 = 1.
1052 IF ( XXX2 <= -5. ) XXX2 = -5.
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 !
1066 XXX10 = (10./Z)*XXX
1067 IF ( XXX10 >= 1. ) XXX10 = 1.
1068 IF ( XXX10 <= -5. ) XXX10 = -5.
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
1079 PSIX = ALOG(Z/Z0) - PSIM
1080 PSIT = ALOG(Z/Z0H) - PSIH
1082 PSIX2 = ALOG(2./Z0) - PSIM2
1083 PSIT2 = ALOG(2./Z0H) - PSIH2
1085 PSIX10 = ALOG(10./Z0) - PSIM10
1086 PSIT10 = ALOG(10./Z0H) - PSIH10
1088 U10 = U1 * (PSIX10/PSIX) ! u at 10 m [m/s]
1089 V10 = V1 * (PSIX10/PSIX) ! v at 10 m [m/s]
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)
1096 Q2 = QS + (QA-QS)*(PSIT2/PSIT) ! humidity at 2 m [-]
1098 ! TS = (LW/SIG_SI/0.88)**0.25 ! Radiative temperature [K]
1100 END SUBROUTINE urban
1101 !===============================================================================
1102 !
1103 ! mos
1104 !
1105 !===============================================================================
1106 SUBROUTINE mos(XXX,ALPHA,CD,B1,RIB,Z,Z0,UA,TA,TSF,RHO)
1108 ! XXX: z/L (requires iteration by Newton-Rapson method)
1109 ! B1: Stanton number
1110 ! PSIM: = PSIX of LSM
1111 ! PSIH: = PSIT of LSM
1113 IMPLICIT NONE
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
1124 IF(RIB <= -15.) RIB=-15.
1126 IF(RIB < 0.) THEN
1128 DO NEWT=1,NEWT_END
1130 IF(XXX >= 0.) XXX=-1.E-3
1132 XXX0=XXX*Z0/(Z+Z0)
1134 X=(1.-16.*XXX)**0.25
1135 X0=(1.-16.*XXX0)**0.25
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.)
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
1152 F=RIB*PSIM**2./PSIH-XXX
1154 DF=RIB*(2.*DPSIM*PSIM*PSIH-DPSIH*PSIM**2.) &
1155 /PSIH**2.-1.
1157 XXXP=XXX
1158 XXX=XXXP-F/DF
1159 IF(XXX <= -10.) XXX=-10.
1161 END DO
1163 ELSE IF(RIB >= 0.142857) THEN
1165 XXX=0.714
1166 PSIM=ALOG((Z+Z0)/Z0)+7.*XXX
1167 PSIH=PSIM+0.4*B1
1169 ELSE
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
1179 END IF
1181 US=0.4*UA/PSIM ! u*
1182 IF(US <= 0.01) US=0.01
1183 TS=0.4*(TA-TSF)/PSIH ! T*
1185 CD=US*US/UA**2. ! CD
1186 ALPHA=RHO*CP*0.4*US/PSIH ! RHO*CP*CH*U
1188 RETURN
1189 END SUBROUTINE mos
1190 !===============================================================================
1191 !
1192 ! louis79
1193 !
1194 !===============================================================================
1195 SUBROUTINE louis79(ALPHA,CD,RIB,Z,Z0,UA,RHO)
1197 IMPLICIT NONE
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
1205 A2=(0.4/ALOG(Z/Z0))**2.
1207 IF(RIB <= -15.) RIB=-15.
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
1224 ALPHA=RHO*CP*CH*UA
1226 RETURN
1227 END SUBROUTINE louis79
1228 !===============================================================================
1229 !
1230 ! louis82
1231 !
1232 !===============================================================================
1233 SUBROUTINE louis82(ALPHA,CD,RIB,Z,Z0,UA,RHO)
1235 IMPLICIT NONE
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
1243 A2=(0.4/ALOG(Z/Z0))**2.
1245 IF(RIB <= -15.) RIB=-15.
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
1260 ALPHA=RHO*CP*CH*UA
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)
1271 IMPLICIT NONE
1273 REAL, INTENT(IN) :: G0
1275 REAL, INTENT(IN) :: CAP
1277 REAL, INTENT(IN) :: AKS
1279 REAL, INTENT(IN) :: DELT ! Time step [ s ]
1281 REAL, INTENT(IN) :: TSLEND
1283 INTEGER, INTENT(IN) :: KM
1285 INTEGER, INTENT(IN) :: BOUND
1287 REAL, DIMENSION(KM), INTENT(IN) :: DZ
1289 REAL, DIMENSION(KM), INTENT(INOUT) :: TSL
1291 REAL, DIMENSION(KM) :: A, B, C, D, X, P, Q
1293 REAL :: DZEND
1295 INTEGER :: K
1297 DZEND=DZ(KM)
1299 A(1) = 0.0
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
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
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
1325 P(1) = -C(1)/B(1)
1326 Q(1) = D(1)/B(1)
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
1333 X(KM) = Q(KM)
1335 DO K=KM-1,1,-1
1336 X(K) = P(K)*X(K+1)+Q(K)
1337 END DO
1339 DO K=1,KM
1340 TSL(K) = X(K)
1341 END DO
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
1363 INTEGER, INTENT(IN) :: UTYPE
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
1380 !end BEP
1382 INTEGER, INTENT(OUT) :: BOUNDR,BOUNDB,BOUNDG,CH_SCHEME,TS_SCHEME
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)
1398 !m FRC_URB= FRC_URB_TBL(UTYPE)
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)
1427 !for BEP
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
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 )
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)
1454 IMPLICIT NONE
1456 INTEGER, INTENT(IN) :: num_soil_layers
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
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
1474 REAL, allocatable, dimension(:) :: ROOF_WIDTH
1475 REAL, allocatable, dimension(:) :: ROAD_WIDTH
1477 character(len=512) :: string
1478 character(len=128) :: name
1479 integer :: indx
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
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
1500 ICATE=0
1502 OPEN (UNIT=11, &
1503 FILE='URBPARM.TBL', &
1504 ACCESS='SEQUENTIAL', &
1505 STATUS='OLD', &
1506 ACTION='READ', &
1507 POSITION='REWIND', &
1508 IOSTAT=IOSTATUS)
1510 IF (IOSTATUS > 0) THEN
1511 CALL wrf_error_fatal('ERROR OPEN URBPARM.TBL')
1512 ENDIF
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)))
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
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')
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
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
1756 else if (name == "BUILDING HEIGHTS") then
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
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
1810 CLOSE(11)
1812 ! Assign a few table values that do not need to come from URBPARM.TBL
1814 Z0HB_TBL = 0.1 * Z0B_TBL
1815 Z0HG_TBL = 0.1 * Z0G_TBL
1817 DO LC = 1, ICATE
1819 ! HGT: Normalized height
1820 HGT_TBL(LC) = ZR_TBL(LC) / ( ROAD_WIDTH(LC) + ROOF_WIDTH(LC) )
1822 ! R: Normalized Roof Width (a.k.a. "building coverage ratio")
1823 R_TBL(LC) = ROOF_WIDTH(LC) / ( ROAD_WIDTH(LC) + ROOF_WIDTH(LC) )
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
1830 ! The following urban canyon geometry parameters are following Macdonald's (1998) formulations
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 )
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
1845 ZDC_TBL(LC) = ZR_TBL(LC) * ( 1.0 + ( alpha_macd ** ( -Lambda_P ) ) * ( Lambda_P - 1.0 ) )
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))
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
1859 !
1860 ! Z0HC still one-tenth of Z0C, as before ?
1861 !
1863 Z0HC_TBL(LC) = 0.1 * Z0C_TBL(LC)
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
1877 VFWS=VFWS/99.
1878 VFWS=VFWS*2.
1880 VFGS=1.-2.*VFWS*HGT_TBL(LC)/RW_TBL(LC)
1881 SVF_TBL(LC)=VFGS
1882 END DO
1884 deallocate(roof_width)
1885 deallocate(road_width)
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
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
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
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
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
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
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
1993 INTEGER :: I,J,K
1995 DO I=ims,ime
1996 DO J=jms,jme
1998 ! XXXR_URB2D(I,J)=0.
1999 ! XXXB_URB2D(I,J)=0.
2000 ! XXXG_URB2D(I,J)=0.
2001 ! XXXC_URB2D(I,J)=0.
2003 SH_URB2D(I,J)=0.
2004 LH_URB2D(I,J)=0.
2005 G_URB2D(I,J)=0.
2006 RN_URB2D(I,J)=0.
2008 !m
2009 FRC_URB2D(I,J)=0.
2010 UTYPE_URB2D(I,J)=0
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
2034 QC_URB2D(I,J)=0.01
2036 IF (.not.restart) THEN
2038 XXXR_URB2D(I,J)=0.
2039 XXXB_URB2D(I,J)=0.
2040 XXXG_URB2D(I,J)=0.
2041 XXXC_URB2D(I,J)=0.
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.
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.
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
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.
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
2076 ! END DO
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
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
2110 ENDIF
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.
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
2135 endif
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)
2168 REAL, INTENT(IN) :: CAP,AKS,DELT,S,R,H,LE,TSLEND,TSP
2169 REAL, INTENT(OUT) :: TS
2170 REAL :: C1,C2
2172 C2=24.*3600./2./3.14159
2173 C1=SQRT(0.5*C2*CAP*AKS)
2175 TS = TSP + DELT*( (S+R-H-LE)/C1 -(TSP-TSLEND)/C2 )
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)
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
2189 TS1 = TSP - 5.
2190 TS2 = TSP + 5.
2192 DO ITERATION = 1,22
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
2202 TS=0.5*(TS1+TS2)
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
2212 IF (F1*F > 0.0) THEN
2213 TS1=TS
2214 ELSE
2215 TS2=TS
2216 END IF
2218 END DO
2220 RETURN
2221 END SUBROUTINE bisection
2222 !===========================================================================
2224 SUBROUTINE SFCDIF_URB (ZLM,Z0,THZ0,THLM,SFCSPD,AKANDA,AKMS,AKHS,RLMO,CD)
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 ! ----------------------------------------------------------------------
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
2245 REAL XLU, XLT, XU, XT, PSMZ, SIMM, PSHZ, SIMH, USTARK, RLMN, &
2246 & RLMA
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))
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)
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)
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
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
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
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
2316 ! ----------------------------------------------------------------------
2317 ! ZILITINKEVITCH APPROACH FOR ZT
2318 ! ----------------------------------------------------------------------
2319 USTAR = MAX (SQRT (AKMS * SQRT (DU2+ WSTAR2)),EPSUST)
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
2326 ZSLU = ZLM + ZU
2328 ZSLT = ZLM + ZT
2329 RLOGU = log (ZSLU / ZU)
2331 RLOGT = log (ZSLT / ZT)
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
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
2350 XT4 = 1. -16.* ZETAT
2351 XLU = SQRT (SQRT (XLU4))
2352 XLT = SQRT (SQRT (XLT4))
2353 XU = SQRT (SQRT (XU4))
2355 XT = SQRT (SQRT (XT4))
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
2414 END DO
2416 CD = USTAR*USTAR/SFCSPD**2
2417 ! ----------------------------------------------------------------------
2418 END SUBROUTINE SFCDIF_URB
2419 ! ----------------------------------------------------------------------
2420 !===========================================================================
2421 END MODULE module_sf_urban