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