Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

module_sf_urban.F @ 3026

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

Rev	Line
[2759]	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