Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_sf_urban.F @ 3094

Last change on this file since 3094 was 11, checked in by aslmd, 14 years ago
spiga@svn-planeto:ajoute le modele meso-echelle martien
File size: 58.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
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, = FLXTHRHOOCPP
112	! LH [W/m/m] : Latent heat flux, = FLXHUMRHOOELL
113	! LH_INEMATIC [kg/m/m/sec]: Moisture Kinematic flux, = FLXHUM*RHOO
114	! SW [W/m/m] : Upward shortwave radiation flux,
115	! = SSG-SNET697.760. (697.760.=100.100.*4.186)
116	! ALB [-] : Time-varying albedo
117	! LW [W/m/m] : Upward longwave radiation flux,
118	! = LNET697.760.-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] : CHU at ZA, CHU 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.61QA)
451	PS=RHOO287.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(CDSAS/(2.XLB2.))*(1./3.)
463	UC=UREXP(-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=SXVFGS(1.-ALBG)
484	SB1=SXVFWS(1.-ALBB)
485	SG2=SB1ALBB/(1.-ALBB)VFGW*(1.-ALBG)
486	SB2=SG1ALBG/(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=(SNTCOS(PI/8.) -CNTSIN(PI/8.))
499	HOUI2=(SNTCOS(2.PI/8.) -CNTSIN(2.PI/8.))
500	HOUI3=(SNTCOS(3.PI/8.) -CNTSIN(3.PI/8.))
501	HOUI4=(SNTCOS(4.PI/8.) -CNTSIN(4.PI/8.))
502	HOUI5=(SNTCOS(5.PI/8.) -CNTSIN(5.PI/8.))
503	HOUI6=(SNTCOS(6.PI/8.) -CNTSIN(6.PI/8.))
504	HOUI7=(SNTCOS(7.PI/8.) -CNTSIN(7.PI/8.))
505	HOUI8=(SNTCOS(8.PI/8.) -CNTSIN(8.PI/8.))
506
507	SLX1=HGTABS(TAN(THEATAZ))ABS(HOUI1)
508	SLX2=HGTABS(TAN(THEATAZ))ABS(HOUI2)
509	SLX3=HGTABS(TAN(THEATAZ))ABS(HOUI3)
510	SLX4=HGTABS(TAN(THEATAZ))ABS(HOUI4)
511	SLX5=HGTABS(TAN(THEATAZ))ABS(HOUI5)
512	SLX6=HGTABS(TAN(THEATAZ))ABS(HOUI6)
513	SLX7=HGTABS(TAN(THEATAZ))ABS(HOUI7)
514	SLX8=HGTABS(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=RSR+WSB+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.82./(TA+TRP))(TA-TRP)(Z+Z0R)/(UAUA)
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.1110.(TETENATSC/(TETENB+TSC)) ! Tetens
571	! DESDT=( 6.1078(2500.-2.4TSC)/ & ! Tetens
572	! (0.46151(TSC+273.15)2.) )10.*(7.5TSC/(237.3+TSC))
573	! ES=6.11EXP((2.510.*6./461.51)(TRP-273.15)/(273.15*TRP) ) ! Clausius-Clapeyron
574	! DESDT=(2.510.6./461.51)ES/(TRP**2.) ! Clausius-Clapeyron
575	! QS0R=0.622ES/(PS-0.378ES)
576	! DQS0RDTR = DESDT0.622PS/((PS-0.378ES)*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.11EXP( (2.510.*6./461.51)(TRP-273.15)/(273.15*TRP) )
584	DESDT=(2.510.6./461.51)ES/(TRP**2.)
585	QS0R=0.622ES/(PS-0.378ES)
586	DQS0RDTR = DESDT0.622PS/((PS-0.378ES)*2.)
587
588	RR=EPSR(RX-SIG(TRP**4.)/60.)
589	HR=RHOCPCHRUA(TRP-TA)*100.
590	ELER=RHOELCHRUABETR(QS0R-QA)100.
591	G0R=AKSR*(TRP-TRL(1))/(DZR(1)/2.)
592
593	F = SR + RR - HR - ELER - G0R
594
595	DRRDTR = (-4.EPSRSIGTRP*3.)/60.
596	DHRDTR = RHOCPCHRUA100.
597	DELERDTR = RHOELCHRUABETRDQS0RDTR100.
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.11EXP( (2.510.*6./461.51)(TRP-273.15)/(273.15*TRP) )
617	QS0R=0.622ES/(PS-0.378ES)
618
619	RR=EPSR(RX-SIG(TRP**4.)/60.)
620	HR=RHOCPCHRUA(TRP-TA)*100.
621	ELER=RHOELCHRUABETR(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.82./(TA+TCP))(TA-TCP)(Z+Z0C)/(UAUA)
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.82./(TCP+TBP))(TCP-TBP)(Z+Z0B)/(UCUC)
653	RIBG=(9.82./(TCP+TGP))(TCP-TGP)(Z+Z0G)/(UCUC)
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=RHOCP(6.15+4.18*UC)/1200.
661	IF(UC > 5.) ALPHAB=RHOCP(7.51UC*0.78)/1200.
662	ALPHAG=RHOCP(6.15+4.18*UC)/1200.
663	IF(UC > 5.) ALPHAG=RHOCP(7.51UC*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.11EXP( (2.510.*6./461.51)(TBP-273.15)/(273.15*TBP) )
687	DESDT=(2.510.6./461.51)ES/(TBP**2.)
688	QS0B=0.622ES/(PS-0.378ES)
689	DQS0BDTB=DESDT0.622PS/((PS-0.378ES)*2.)
690
691	ES=6.11EXP( (2.510.*6./461.51)(TGP-273.15)/(273.15*TGP) )
692	DESDT=(2.510.6./461.51)ES/(TGP**2.)
693	QS0G=0.622ES/(PS-0.378ES)
694	DQS0GDTG=DESDT0.22PS/((PS-0.378ES)*2.)
695
696	RG1=EPSG( RXVFGS &
697	+EPSBVFGWSIGTBP*4./60. &
698	-SIGTGP*4./60. )
699
700	RB1=EPSB( RXVFWS &
701	+EPSGVFWGSIGTGP*4./60. &
702	+EPSBVFWWSIGTBP*4./60. &
703	-SIGTBP*4./60. )
704
705	RG2=EPSG( (1.-EPSB)(1.-SVF)VFWSRX &
706	+(1.-EPSB)(1.-SVF)VFWGEPSGSIGTGP*4./60. &
707	+EPSB(1.-EPSB)(1.-SVF)(1.-2.VFWS)SIGTBP**4./60. )
708
709	RB2=EPSB( (1.-EPSG)VFWGVFGSRX &
710	+(1.-EPSG)EPSBVFGWVFWGSIG(TBP*4.)/60. &
711	+(1.-EPSB)VFWS(1.-2.VFWS)RX &
712	+(1.-EPSB)VFWG(1.-2.VFWS)EPSGSIGEPSGTGP*4./60. &
713	+EPSB(1.-EPSB)(1.-2.VFWS)(1.-2.VFWS)SIGTBP*4./60. )
714
715	RG=RG1+RG2
716	RB=RB1+RB2
717
718	DRBDTB1=EPSB(4.EPSBSIGTB*3.VFWW-4.SIGTB**3.)/60.
719	DRBDTG1=EPSB(4.EPSGSIGTG*3.VFWG)/60.
720	DRBDTB2=EPSB(4.(1.-EPSG)EPSBSIGTB3.VFGW*VFWG &
721	+4.EPSB(1.-EPSB)SIGTB*3.VFWW*VFWW)/60.
722	DRBDTG2=EPSB(4.(1.-EPSB)EPSGSIGTG3.VFWG*VFWW)/60.
723
724	DRGDTB1=EPSG(4.EPSBSIGTB*3.VFGW)/60.
725	DRGDTG1=EPSG(-4.SIGTG*3.)/60.
726	DRGDTB2=EPSG(4.EPSB(1.-EPSB)SIGTB3.VFWW*VFGW)/60.
727	DRGDTG2=EPSG(4.(1.-EPSB)EPSGSIGTG3.VFWG*VFGW)/60.
728
729	DRBDTB=DRBDTB1+DRBDTB2
730	DRBDTG=DRBDTG1+DRBDTG2
731	DRGDTB=DRGDTB1+DRGDTB2
732	DRGDTG=DRGDTG1+DRGDTG2
733
734	HB=RHOCPCHBUC(TBP-TCP)*100.
735	HG=RHOCPCHGUC(TGP-TCP)*100.
736
737	DTCDTB=WALPHAB/(RWALPHAC+RWALPHAG+WALPHAB)
738	DTCDTG=RWALPHAG/(RWALPHAC+RWALPHAG+WALPHAB)
739
740	DHBDTB=RHOCPCHBUC(1.-DTCDTB)*100.
741	DHBDTG=RHOCPCHBUC(0.-DTCDTG)*100.
742	DHGDTG=RHOCPCHGUC(1.-DTCDTG)*100.
743	DHGDTB=RHOCPCHGUC(0.-DTCDTB)*100.
744
745	ELEB=RHOELCHBUCBETB(QS0B-QCP)100.
746	ELEG=RHOELCHGUCBETG(QS0G-QCP)100.
747
748	DQCDTB=WALPHABBETBDQS0BDTB/(RWALPHAC+RWALPHAGBETG+WALPHABBETB)
749	DQCDTG=RWALPHAGBETGDQS0GDTG/(RWALPHAC+RWALPHAGBETG+WALPHABBETB)
750
751	DELEBDTB=RHOELCHBUCBETB(DQS0BDTB-DQCDTB)100.
752	DELEBDTG=RHOELCHBUCBETB(0.-DQCDTG)100.
753	DELEGDTG=RHOELCHGUCBETG(DQS0GDTG-DQCDTG)100.
754	DELEGDTB=RHOELCHGUCBETG(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 = (GFFY-FGY)/(FXGY-GXFY)
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=RWALPHAC+RWALPHAG+W*ALPHAB
782	TC2=RWALPHACTA+RWALPHAGTGP+WALPHABTBP
783	TC=TC2/TC1
784
785	QC1=RWALPHAC+RWALPHAGBETG+WALPHAB*BETB
786	QC2=RWALPHACQA+RWALPHAGBETGQS0G+WALPHABBETBQS0B
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.11EXP((2.510.*6./461.51)(TBP-273.15)/(273.15*TBP) )
810	QS0B=0.622ES/(PS-0.378ES)
811
812	ES=6.11EXP((2.510.*6./461.51)(TGP-273.15)/(273.15*TGP) )
813	QS0G=0.622ES/(PS-0.378ES)
814
815	RG1=EPSG( RXVFGS &
816	+EPSBVFGWSIGTBP*4./60. &
817	-SIGTGP*4./60. )
818
819	RB1=EPSB( RXVFWS &
820	+EPSGVFWGSIGTGP*4./60. &
821	+EPSBVFWWSIGTBP*4./60. &
822	-SIGTBP*4./60. )
823
824	RG2=EPSG( (1.-EPSB)(1.-SVF)VFWSRX &
825	+(1.-EPSB)(1.-SVF)VFWGEPSGSIGTGP*4./60. &
826	+EPSB(1.-EPSB)(1.-SVF)(1.-2.VFWS)SIGTBP**4./60. )
827
828	RB2=EPSB( (1.-EPSG)VFWGVFGSRX &
829	+(1.-EPSG)EPSBVFGWVFWGSIG(TBP*4.)/60. &
830	+(1.-EPSB)VFWS(1.-2.VFWS)RX &
831	+(1.-EPSB)VFWG(1.-2.VFWS)EPSGSIGEPSGTGP*4./60. &
832	+EPSB(1.-EPSB)(1.-2.VFWS)(1.-2.VFWS)SIGTBP*4./60. )
833
834	RG=RG1+RG2
835	RB=RB1+RB2
836
837	HB=RHOCPCHBUC(TBP-TCP)*100.
838	ELEB=RHOELCHBUCBETB(QS0B-QCP)100.
839	G0B=SB+RB-HB-ELEB
840
841	HG=RHOCPCHGUC(TGP-TCP)*100.
842	ELEG=RHOELCHGUCBETG(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=RWALPHAC+RWALPHAG+W*ALPHAB
852	TC2=RWALPHACTA+RWALPHAGTGP+WALPHABTBP
853	TC=TC2/TC1
854
855	QC1=RWALPHAC+RWALPHAGBETG+WALPHAB*BETB
856	QC2=RWALPHACQA+RWALPHAGBETGQS0G+WALPHABBETBQS0B
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 = ( RCDR + RWCDC )UAUA
874	FLXTH = ( RFLXTHR + WFLXTHB + RW*FLXTHG )
875	FLXHUM = ( RFLXHUMR + WFLXHUMB + RW*FLXHUMG )
876	FLXG = ( RG0R + WG0B + RW*G0G )
877	LNET = RRR + WRB + 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 - (LNET697.760.) ! Upward longwave radiation [W/m/m]
887	SW = SSG - (SNET697.760.) ! Upward shortwave radiation [W/m/m]
888	ALB = 0.
889	IF( ABS(SSG) > 0.0001) ALB = SW/SSG ! Effective albedo [-]
890	G = -FLXG697.760. ! [W/m/m]
891	RN = (SNET+LNET)697.760. ! 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.49.81Z*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.+XX)/2.) - 2.*ATAN(X) + PI/2.
916	PSIH = 2.ALOG((1.+XX)/2.)
917	END IF
918
919	GZ1OZ0 = ALOG(Z/Z0)
920	CD = 0.42./(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.+XX)/2.) - 2.ATAN(X) + 2.ATAN(1.)
944	PSIH2 = 2.ALOG((1.+XX)/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.+XX)/2.) - 2.ATAN(X) + 2.ATAN(1.)
960	PSIH10 = 2.ALOG((1.+XX)/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=RIBPSIM*2./PSIH-XXX
1039
1040	DF=RIB(2.DPSIMPSIMPSIH-DPSIHPSIM*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.RIB7.XKBAL+(AL+XKB)**2.
1060	IF(DD <= 0.) DD=0.
1061	XXX=(AL+XKB-2.RIB7.AL-SQRT(DD))/(2.(RIB7.*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.4UA/PSIM ! u
1068	IF(US <= 0.01) US=0.01
1069	TS=0.4(TA-TSF)/PSIH ! T
1070
1071	CD=USUS/UA*2. ! CD
1072	ALPHA=RHOCP0.4US/PSIH ! RHOCPCHU
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=RIBLOG(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.4A29.4*SQRT(Z/Z0)
1105	CHB=5.3A29.4*SQRT(Z/Z0)
1106	CH=A2/0.74(1.-9.4RIB/(1.+CHB*SQRT(-RIB)))
1107	CD=A2(1.-9.4RIB/(1.+CHB*SQRT(-RIB)))
1108	END IF
1109
1110	ALPHA=RHOCPCH*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.A2SQRT(Z/Z0)
1140	FM=1.-(2.5.RIB)/(1.+3.5.5.A2SQRT(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=RHOCPCH*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) = CAPDZ(1)/DELTTSL(1) + G0
1181
1182	DO K=2,KM-1
1183	A(K) = -2.*AKS/(DZ(K-1)+DZ(K))
1184	B(K) = CAPDZ(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) = CAPDZ(K)/DELTTSL(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) = CAPDZ(KM)/DELT + 2.AKS/(DZ(KM-1)+DZ(KM))
1192	C(KM) = 0.0
1193	D(KM) = CAPDZ(KM)/DELTTSL(KM)
1194	ELSE ! T=constant
1195	A(KM) = -2.*AKS/(DZ(KM-1)+DZ(KM))
1196	B(KM) = CAPDZ(KM)/DELT + 2.AKS/(DZ(KM-1)+DZ(KM)) + 2.*AKS/(DZ(KM)+DZEND)
1197	C(KM) = 0.0
1198	D(KM) = CAPDZ(KM)/DELTTSL(KM) + 2.AKSTSLEND/(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(HGT2.+RW_TBL(LC)*2.))
1488	end do
1489
1490	VFWS=VFWS/99.
1491	VFWS=VFWS*2.
1492
1493	VFGS=1.-2.VFWSHGT_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.5C2CAP*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.11EXP( (2.510.*6./461.51)(TS1-273.15)/(273.15*TS1) )
1667	QS0=0.622ES/(PS-0.378ES)
1668	R=EPS(RX-SIG(TS1**4.)/60.)
1669	H=RHOCPCHUA(TS1-TA)*100.
1670	ELE=RHOELCHUABET(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.11EXP( (2.510.*6./461.51)(TS-273.15)/(273.15*TS) )
1677	QS0=0.622ES/(PS-0.378ES)
1678	R=EPS(RX-SIG(TS**4.)/60.)
1679	H=RHOCPCHUA(TS-TA)*100.
1680	ELE=RHOELCHUABET(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

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

Download in other formats:

Original Format