Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_sf_urban.F @ 3094

Last change on this file since 3094 was 2759, checked in by aslmd, 2 years ago
adding unmodified code from WRFV3.0.1.1, expurged from useless data +1M size
File size: 59.8 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format