Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

module_sf_noahdrv.F @ 13

Last change on this file since 13 was 1, checked in by lfita, 11 years ago

-- --- Opening of the WRF+LMDZ coupling repository --- -- -

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 81.2 KB

Rev	Line
[1]	1	MODULE module_sf_noahdrv
	2
	3	!-------------------------------
	4	USE module_sf_noahlsm
	5	USE module_sf_urban
	6	USE module_sf_bep
	7	USE module_sf_bep_bem
	8	#ifdef WRF_CHEM
	9	USE module_data_gocart_dust
	10	#endif
	11	!-------------------------------
	12
	13	!
	14	CONTAINS
	15	!
	16	!----------------------------------------------------------------
	17	! Urban related variable are added to arguments - urban
	18	!----------------------------------------------------------------
	19	SUBROUTINE lsm(DZ8W,QV3D,P8W3D,T3D,TSK, &
	20	HFX,QFX,LH,GRDFLX, QGH,GSW,SWDOWN,GLW,SMSTAV,SMSTOT, &
	21	SFCRUNOFF, UDRUNOFF,IVGTYP,ISLTYP,ISURBAN,ISICE,VEGFRA, &
	22	ALBEDO,ALBBCK,ZNT,Z0,TMN,XLAND,XICE,EMISS,EMBCK, &
	23	SNOWC,QSFC,RAINBL,MMINLU, &
	24	num_soil_layers,DT,DZS,ITIMESTEP, &
	25	SMOIS,TSLB,SNOW,CANWAT, &
	26	CHS,CHS2,CQS2,CPM,ROVCP,SR,chklowq,lai,qz0, & !H
	27	myj,frpcpn, &
	28	SH2O,SNOWH, & !H
	29	U_PHY,V_PHY, & !I
	30	SNOALB,SHDMIN,SHDMAX, & !I
	31	SNOTIME, & !?
	32	ACSNOM,ACSNOW, & !O
	33	SNOPCX, & !O
	34	POTEVP, & !O
	35	SMCREL, & !O
	36	XICE_THRESHOLD, &
	37	RDLAI2D,USEMONALB, &
	38	RIB, & !?
	39	NOAHRES, &
	40	ids,ide, jds,jde, kds,kde, &
	41	ims,ime, jms,jme, kms,kme, &
	42	its,ite, jts,jte, kts,kte, &
	43	sf_urban_physics, &
	44	CMR_SFCDIF,CHR_SFCDIF,CMC_SFCDIF,CHC_SFCDIF, &
	45	!Optional Urban
	46	TR_URB2D,TB_URB2D,TG_URB2D,TC_URB2D,QC_URB2D, & !H urban
	47	UC_URB2D, & !H urban
	48	XXXR_URB2D,XXXB_URB2D,XXXG_URB2D,XXXC_URB2D, & !H urban
	49	TRL_URB3D,TBL_URB3D,TGL_URB3D, & !H urban
	50	SH_URB2D,LH_URB2D,G_URB2D,RN_URB2D,TS_URB2D, & !H urban
	51	PSIM_URB2D,PSIH_URB2D,U10_URB2D,V10_URB2D, & !O urban
	52	GZ1OZ0_URB2D, AKMS_URB2D, & !O urban
	53	TH2_URB2D,Q2_URB2D, UST_URB2D, & !O urban
	54	DECLIN_URB,COSZ_URB2D,OMG_URB2D, & !I urban
	55	XLAT_URB2D, & !I urban
	56	num_roof_layers, num_wall_layers, & !I urban
	57	num_road_layers, DZR, DZB, DZG, & !I urban
	58	FRC_URB2D,UTYPE_URB2D, & !O
	59	num_urban_layers, & !I multi-layer urban
	60	trb_urb4d,tw1_urb4d,tw2_urb4d,tgb_urb4d, & !H multi-layer urban
	61	tlev_urb3d,qlev_urb3d, & !H multi-layer urban
	62	tw1lev_urb3d,tw2lev_urb3d, & !H multi-layer urban
	63	tglev_urb3d,tflev_urb3d, & !H multi-layer urban
	64	sf_ac_urb3d,lf_ac_urb3d,cm_ac_urb3d, & !H multi-layer urban
	65	sfvent_urb3d,lfvent_urb3d, & !H multi-layer urban
	66	sfwin1_urb3d,sfwin2_urb3d, & !H multi-layer urban
	67	sfw1_urb3d,sfw2_urb3d,sfr_urb3d,sfg_urb3d, & !H multi-layer urban
	68	th_phy,rho,p_phy,ust, & !I multi-layer urban
	69	gmt,julday,xlong,xlat, & !I multi-layer urban
	70	a_u_bep,a_v_bep,a_t_bep,a_q_bep, & !O multi-layer urban
	71	a_e_bep,b_u_bep,b_v_bep, & !O multi-layer urban
	72	b_t_bep,b_q_bep,b_e_bep,dlg_bep, & !O multi-layer urban
	73	dl_u_bep,sf_bep,vl_bep ) !O multi-layer urban
	74
	75	!----------------------------------------------------------------
	76	IMPLICIT NONE
	77	!----------------------------------------------------------------
	78	!----------------------------------------------------------------
	79	! --- atmospheric (WRF generic) variables
	80	!-- DT time step (seconds)
	81	!-- DZ8W thickness of layers (m)
	82	!-- T3D temperature (K)
	83	!-- QV3D 3D water vapor mixing ratio (Kg/Kg)
	84	!-- P3D 3D pressure (Pa)
	85	!-- FLHC exchange coefficient for heat (m/s)
	86	!-- FLQC exchange coefficient for moisture (m/s)
	87	!-- PSFC surface pressure (Pa)
	88	!-- XLAND land mask (1 for land, 2 for water)
	89	!-- QGH saturated mixing ratio at 2 meter
	90	!-- GSW downward short wave flux at ground surface (W/m^2)
	91	!-- GLW downward long wave flux at ground surface (W/m^2)
	92	!-- History variables
	93	!-- CANWAT canopy moisture content (mm)
	94	!-- TSK surface temperature (K)
	95	!-- TSLB soil temp (k)
	96	!-- SMOIS total soil moisture content (volumetric fraction)
	97	!-- SH2O unfrozen soil moisture content (volumetric fraction)
	98	! note: frozen soil moisture (i.e., soil ice) = SMOIS - SH2O
	99	!-- SNOWH actual snow depth (m)
	100	!-- SNOW liquid water-equivalent snow depth (m)
	101	!-- ALBEDO time-varying surface albedo including snow effect (unitless fraction)
	102	!-- ALBBCK background surface albedo (unitless fraction)
	103	!-- CHS surface exchange coefficient for heat and moisture (m s-1);
	104	!-- CHS2 2m surface exchange coefficient for heat (m s-1);
	105	!-- CQS2 2m surface exchange coefficient for moisture (m s-1);
	106	! --- soil variables
	107	!-- num_soil_layers the number of soil layers
	108	!-- ZS depths of centers of soil layers (m)
	109	!-- DZS thicknesses of soil layers (m)
	110	!-- SLDPTH thickness of each soil layer (m, same as DZS)
	111	!-- TMN soil temperature at lower boundary (K)
	112	!-- SMCWLT wilting point (volumetric)
	113	!-- SMCDRY dry soil moisture threshold where direct evap from
	114	! top soil layer ends (volumetric)
	115	!-- SMCREF soil moisture threshold below which transpiration begins to
	116	! stress (volumetric)
	117	!-- SMCMAX porosity, i.e. saturated value of soil moisture (volumetric)
	118	!-- NROOT number of root layers, a function of veg type, determined
	119	! in subroutine redprm.
	120	!-- SMSTAV Soil moisture availability for evapotranspiration (
	121	! fraction between SMCWLT and SMCMXA)
	122	!-- SMSTOT Total soil moisture content frozen+unfrozen) in the soil column (mm)
	123	! --- snow variables
	124	!-- SNOWC fraction snow coverage (0-1.0)
	125	! --- vegetation variables
	126	!-- SNOALB upper bound on maximum albedo over deep snow
	127	!-- SHDMIN minimum areal fractional coverage of annual green vegetation
	128	!-- SHDMAX maximum areal fractional coverage of annual green vegetation
	129	!-- XLAI leaf area index (dimensionless)
	130	!-- Z0BRD Background fixed roughness length (M)
	131	!-- Z0 Background vroughness length (M) as function
	132	!-- ZNT Time varying roughness length (M) as function
	133	!-- ALBD(IVGTPK,ISN) background albedo reading from a table
	134	! --- LSM output
	135	!-- HFX upward heat flux at the surface (W/m^2)
	136	!-- QFX upward moisture flux at the surface (kg/m^2/s)
	137	!-- LH upward moisture flux at the surface (W m-2)
	138	!-- GRDFLX(I,J) ground heat flux (W m-2)
	139	!-- FDOWN radiation forcing at the surface (W m-2) = SOLDN*(1-alb)+LWDN
	140	!----------------------------------------------------------------------------
	141	!-- EC canopy water evaporation ((W m-2)
	142	!-- EDIR direct soil evaporation (W m-2)
	143	!-- ET plant transpiration from a particular root layer (W m-2)
	144	!-- ETT total plant transpiration (W m-2)
	145	!-- ESNOW sublimation from (or deposition to if <0) snowpack (W m-2)
	146	!-- DRIP through-fall of precip and/or dew in excess of canopy
	147	! water-holding capacity (m)
	148	!-- DEW dewfall (or frostfall for t<273.15) (M)
	149	!-- SMAV Soil Moisture Availability for each layer, as a fraction
	150	! between SMCWLT and SMCMAX (dimensionless fraction)
	151	! ----------------------------------------------------------------------
	152	!-- BETA ratio of actual/potential evap (dimensionless)
	153	!-- ETP potential evaporation (W m-2)
	154	! ----------------------------------------------------------------------
	155	!-- FLX1 precip-snow sfc (W m-2)
	156	!-- FLX2 freezing rain latent heat flux (W m-2)
	157	!-- FLX3 phase-change heat flux from snowmelt (W m-2)
	158	! ----------------------------------------------------------------------
	159	!-- ACSNOM snow melt (mm) (water equivalent)
	160	!-- ACSNOW accumulated snow fall (mm) (water equivalent)
	161	!-- SNOPCX snow phase change heat flux (W/m^2)
	162	!-- POTEVP accumulated potential evaporation (W/m^2)
	163	!-- RIB Documentation needed!!!
	164	! ----------------------------------------------------------------------
	165	!-- RUNOFF1 surface runoff (m s-1), not infiltrating the surface
	166	!-- RUNOFF2 subsurface runoff (m s-1), drainage out bottom of last
	167	! soil layer (baseflow)
	168	! important note: here RUNOFF2 is actually the sum of RUNOFF2 and RUNOFF3
	169	!-- RUNOFF3 numerical trunctation in excess of porosity (smcmax)
	170	! for a given soil layer at the end of a time step (m s-1).
	171	!SFCRUNOFF Surface Runoff (mm)
	172	!UDRUNOFF Total Underground Runoff (mm), which is the sum of RUNOFF2 and RUNOFF3
	173	! ----------------------------------------------------------------------
	174	!-- RC canopy resistance (s m-1)
	175	!-- PC plant coefficient (unitless fraction, 0-1) where PC*ETP = actual transp
	176	!-- RSMIN minimum canopy resistance (s m-1)
	177	!-- RCS incoming solar rc factor (dimensionless)
	178	!-- RCT air temperature rc factor (dimensionless)
	179	!-- RCQ atmos vapor pressure deficit rc factor (dimensionless)
	180	!-- RCSOIL soil moisture rc factor (dimensionless)
	181
	182	!-- EMISS surface emissivity (between 0 and 1)
	183	!-- EMBCK Background surface emissivity (between 0 and 1)
	184
	185	!-- ROVCP R/CP
	186	! (R_d/R_v) (dimensionless)
	187	!-- ids start index for i in domain
	188	!-- ide end index for i in domain
	189	!-- jds start index for j in domain
	190	!-- jde end index for j in domain
	191	!-- kds start index for k in domain
	192	!-- kde end index for k in domain
	193	!-- ims start index for i in memory
	194	!-- ime end index for i in memory
	195	!-- jms start index for j in memory
	196	!-- jme end index for j in memory
	197	!-- kms start index for k in memory
	198	!-- kme end index for k in memory
	199	!-- its start index for i in tile
	200	!-- ite end index for i in tile
	201	!-- jts start index for j in tile
	202	!-- jte end index for j in tile
	203	!-- kts start index for k in tile
	204	!-- kte end index for k in tile
	205	!
	206	!-- SR fraction of frozen precip (0.0 to 1.0)
	207	!----------------------------------------------------------------
	208
	209	! IN only
	210
	211	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
	212	ims,ime, jms,jme, kms,kme, &
	213	its,ite, jts,jte, kts,kte
	214
	215	INTEGER, INTENT(IN ) :: sf_urban_physics !urban
	216	INTEGER, INTENT(IN ) :: isurban
	217	INTEGER, INTENT(IN ) :: isice
	218
	219	REAL, DIMENSION( ims:ime, jms:jme ) , &
	220	INTENT(IN ) :: TMN, &
	221	XLAND, &
	222	XICE, &
	223	VEGFRA, &
	224	SHDMIN, &
	225	SHDMAX, &
	226	SNOALB, &
	227	GSW, &
	228	SWDOWN, & !added 10 jan 2007
	229	GLW, &
	230	RAINBL, &
	231	EMBCK, &
	232	SR
	233
	234	REAL, DIMENSION( ims:ime, jms:jme ) , &
	235	INTENT(INOUT) :: ALBBCK, &
	236	Z0
	237	CHARACTER(LEN=*), INTENT(IN ) :: MMINLU
	238
	239	REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
	240	INTENT(IN ) :: QV3D, &
	241	p8w3D, &
	242	DZ8W, &
	243	T3D
	244	REAL, DIMENSION( ims:ime, jms:jme ) , &
	245	INTENT(IN ) :: QGH, &
	246	CPM
	247
	248	INTEGER, DIMENSION( ims:ime, jms:jme ) , &
	249	INTENT(IN ) :: IVGTYP, &
	250	ISLTYP
	251
	252	INTEGER, INTENT(IN) :: num_soil_layers,ITIMESTEP
	253
	254	REAL, INTENT(IN ) :: DT,ROVCP
	255
	256	REAL, DIMENSION(1:num_soil_layers), INTENT(IN)::DZS
	257
	258	! IN and OUT
	259
	260	REAL, DIMENSION( ims:ime , 1:num_soil_layers, jms:jme ), &
	261	INTENT(INOUT) :: SMOIS, & ! total soil moisture
	262	SH2O, & ! new soil liquid
	263	TSLB ! TSLB STEMP
	264
	265	REAL, DIMENSION( ims:ime , 1:num_soil_layers, jms:jme ), &
	266	INTENT(OUT) :: SMCREL
	267
	268	REAL, DIMENSION( ims:ime, jms:jme ) , &
	269	INTENT(INOUT) :: TSK, & !was TGB (temperature)
	270	HFX, &
	271	QFX, &
	272	LH, &
	273	GRDFLX, &
	274	QSFC,&
	275	CQS2,&
	276	CHS, &
	277	CHS2,&
	278	SNOW, &
	279	SNOWC, &
	280	SNOWH, & !new
	281	CANWAT, &
	282	SMSTAV, &
	283	SMSTOT, &
	284	SFCRUNOFF, &
	285	UDRUNOFF, &
	286	ACSNOM, &
	287	ACSNOW, &
	288	SNOTIME, &
	289	SNOPCX, &
	290	EMISS, &
	291	RIB, &
	292	POTEVP, &
	293	ALBEDO, &
	294	ZNT
	295	REAL, DIMENSION( ims:ime, jms:jme ) , &
	296	INTENT(OUT) :: NOAHRES
	297
	298	REAL, DIMENSION( ims:ime, jms:jme ) , &
	299	INTENT(OUT) :: CHKLOWQ
	300	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: LAI
	301	REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: QZ0
	302
	303	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CMR_SFCDIF
	304	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CHR_SFCDIF
	305	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CMC_SFCDIF
	306	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CHC_SFCDIF
	307	! Local variables (moved here from driver to make routine thread safe, 20031007 jm)
	308
	309	REAL, DIMENSION(1:num_soil_layers) :: ET
	310
	311	REAL, DIMENSION(1:num_soil_layers) :: SMAV
	312
	313	REAL :: BETA, ETP, SSOIL,EC, EDIR, ESNOW, ETT, &
	314	FLX1,FLX2,FLX3, DRIP,DEW,FDOWN,RC,PC,RSMIN,XLAI, &
	315	! RCS,RCT,RCQ,RCSOIL
	316	RCS,RCT,RCQ,RCSOIL,FFROZP
	317
	318	LOGICAL, INTENT(IN ) :: myj,frpcpn
	319
	320	! DECLARATIONS - LOGICAL
	321	! ----------------------------------------------------------------------
	322	LOGICAL, PARAMETER :: LOCAL=.false.
	323	LOGICAL :: FRZGRA, SNOWNG
	324
	325	LOGICAL :: IPRINT
	326
	327	! ----------------------------------------------------------------------
	328	! DECLARATIONS - INTEGER
	329	! ----------------------------------------------------------------------
	330	INTEGER :: I,J, ICE,NSOIL,SLOPETYP,SOILTYP,VEGTYP
	331	INTEGER :: NROOT
	332	INTEGER :: KZ ,K
	333	INTEGER :: NS
	334	! ----------------------------------------------------------------------
	335	! DECLARATIONS - REAL
	336	! ----------------------------------------------------------------------
	337
	338	REAL :: SHMIN,SHMAX,DQSDT2,LWDN,PRCP,PRCPRAIN, &
	339	Q2SAT,Q2SATI,SFCPRS,SFCSPD,SFCTMP,SHDFAC,SNOALB1, &
	340	SOLDN,TBOT,ZLVL, Q2K,ALBBRD, ALBEDOK, ETA, ETA_KINEMATIC, &
	341	EMBRD, &
	342	Z0K,RUNOFF1,RUNOFF2,RUNOFF3,SHEAT,SOLNET,E2SAT,SFCTSNO, &
	343	! mek, WRF testing, expanded diagnostics
	344	SOLUP,LWUP,RNET,RES,Q1SFC,TAIRV,SATFLG
	345	! MEK MAY 2007
	346	REAL :: FDTLIW
	347	! MEK JUL2007 for pot. evap.
	348	REAL :: RIBB
	349	REAL :: FDTW
	350
	351	REAL :: EMISSI
	352
	353	REAL :: SNCOVR,SNEQV,SNOWHK,CMC, CHK,TH2
	354
	355	REAL :: SMCDRY,SMCMAX,SMCREF,SMCWLT,SNOMLT,SOILM,SOILW,Q1,T1
	356	REAL :: SNOTIME1 ! LSTSNW1 INITIAL NUMBER OF TIMESTEPS SINCE LAST SNOWFALL
	357
	358	REAL :: DUMMY,Z0BRD
	359	!
	360	REAL :: COSZ, SOLARDIRECT
	361	!
	362	REAL, DIMENSION(1:num_soil_layers):: SLDPTH, STC,SMC,SWC
	363	!
	364	REAL, DIMENSION(1:num_soil_layers) :: ZSOIL, RTDIS
	365	REAL, PARAMETER :: TRESH=.95E0, A2=17.67,A3=273.15,A4=29.65, &
	366	T0=273.16E0, ELWV=2.50E6, A23M4=A2*(A3-A4)
	367	! MEK MAY 2007
	368	REAL, PARAMETER :: ROW=1.E3,ELIW=XLF,ROWLIW=ROW*ELIW
	369
	370	! ----------------------------------------------------------------------
	371	! DECLARATIONS START - urban
	372	! ----------------------------------------------------------------------
	373
	374	! input variables surface_driver --> lsm
	375	INTEGER, INTENT(IN) :: num_roof_layers
	376	INTEGER, INTENT(IN) :: num_wall_layers
	377	INTEGER, INTENT(IN) :: num_road_layers
	378	REAL, OPTIONAL, DIMENSION(1:num_roof_layers), INTENT(IN) :: DZR
	379	REAL, OPTIONAL, DIMENSION(1:num_wall_layers), INTENT(IN) :: DZB
	380	REAL, OPTIONAL, DIMENSION(1:num_road_layers), INTENT(IN) :: DZG
	381	REAL, OPTIONAL, INTENT(IN) :: DECLIN_URB
	382	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: COSZ_URB2D
	383	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: OMG_URB2D
	384	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: XLAT_URB2D
	385	REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: U_PHY
	386	REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: V_PHY
	387	REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: TH_PHY
	388	REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: P_PHY
	389	REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: RHO
	390	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: UST
	391
	392	LOGICAL, intent(in) :: rdlai2d
	393	LOGICAL, intent(in) :: USEMONALB
	394
	395	! input variables lsm --> urban
	396	INTEGER :: UTYPE_URB ! urban type [urban=1, suburban=2, rural=3]
	397	REAL :: TA_URB ! potential temp at 1st atmospheric level [K]
	398	REAL :: QA_URB ! mixing ratio at 1st atmospheric level [kg/kg]
	399	REAL :: UA_URB ! wind speed at 1st atmospheric level [m/s]
	400	REAL :: U1_URB ! u at 1st atmospheric level [m/s]
	401	REAL :: V1_URB ! v at 1st atmospheric level [m/s]
	402	REAL :: SSG_URB ! downward total short wave radiation [W/m/m]
	403	REAL :: LLG_URB ! downward long wave radiation [W/m/m]
	404	REAL :: RAIN_URB ! precipitation [mm/h]
	405	REAL :: RHOO_URB ! air density [kg/m^3]
	406	REAL :: ZA_URB ! first atmospheric level [m]
	407	REAL :: DELT_URB ! time step [s]
	408	REAL :: SSGD_URB ! downward direct short wave radiation [W/m/m]
	409	REAL :: SSGQ_URB ! downward diffuse short wave radiation [W/m/m]
	410	REAL :: XLAT_URB ! latitude [deg]
	411	REAL :: COSZ_URB ! cosz
	412	REAL :: OMG_URB ! hour angle
	413	REAL :: ZNT_URB ! roughness length [m]
	414	REAL :: TR_URB
	415	REAL :: TB_URB
	416	REAL :: TG_URB
	417	REAL :: TC_URB
	418	REAL :: QC_URB
	419	REAL :: UC_URB
	420	REAL :: XXXR_URB
	421	REAL :: XXXB_URB
	422	REAL :: XXXG_URB
	423	REAL :: XXXC_URB
	424	REAL, DIMENSION(1:num_roof_layers) :: TRL_URB ! roof layer temp [K]
	425	REAL, DIMENSION(1:num_wall_layers) :: TBL_URB ! wall layer temp [K]
	426	REAL, DIMENSION(1:num_road_layers) :: TGL_URB ! road layer temp [K]
	427	LOGICAL :: LSOLAR_URB
	428	! state variable surface_driver <--> lsm <--> urban
	429	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TR_URB2D
	430	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TB_URB2D
	431	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TG_URB2D
	432	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TC_URB2D
	433	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: QC_URB2D
	434	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: UC_URB2D
	435	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXR_URB2D
	436	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXB_URB2D
	437	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXG_URB2D
	438	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXC_URB2D
	439	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: SH_URB2D
	440	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: LH_URB2D
	441	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: G_URB2D
	442	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: RN_URB2D
	443	!
	444	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TS_URB2D
	445
	446	REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_roof_layers, jms:jme ), INTENT(INOUT) :: TRL_URB3D
	447	REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_wall_layers, jms:jme ), INTENT(INOUT) :: TBL_URB3D
	448	REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_road_layers, jms:jme ), INTENT(INOUT) :: TGL_URB3D
	449
	450	! output variable lsm --> surface_driver
	451	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: PSIM_URB2D
	452	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: PSIH_URB2D
	453	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: GZ1OZ0_URB2D
	454	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: U10_URB2D
	455	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: V10_URB2D
	456	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: TH2_URB2D
	457	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: Q2_URB2D
	458	!
	459	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: AKMS_URB2D
	460	!
	461	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: UST_URB2D
	462	REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: FRC_URB2D
	463	INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: UTYPE_URB2D
	464
	465
	466	! output variables urban --> lsm
	467	REAL :: TS_URB ! surface radiative temperature [K]
	468	REAL :: QS_URB ! surface humidity [-]
	469	REAL :: SH_URB ! sensible heat flux [W/m/m]
	470	REAL :: LH_URB ! latent heat flux [W/m/m]
	471	REAL :: LH_KINEMATIC_URB ! latent heat flux, kinetic [kg/m/m/s]
	472	REAL :: SW_URB ! upward short wave radiation flux [W/m/m]
	473	REAL :: ALB_URB ! time-varying albedo [fraction]
	474	REAL :: LW_URB ! upward long wave radiation flux [W/m/m]
	475	REAL :: G_URB ! heat flux into the ground [W/m/m]
	476	REAL :: RN_URB ! net radiation [W/m/m]
	477	REAL :: PSIM_URB ! shear f for momentum [-]
	478	REAL :: PSIH_URB ! shear f for heat [-]
	479	REAL :: GZ1OZ0_URB ! shear f for heat [-]
	480	REAL :: U10_URB ! wind u component at 10 m [m/s]
	481	REAL :: V10_URB ! wind v component at 10 m [m/s]
	482	REAL :: TH2_URB ! potential temperature at 2 m [K]
	483	REAL :: Q2_URB ! humidity at 2 m [-]
	484	REAL :: CHS_URB
	485	REAL :: CHS2_URB
	486	REAL :: UST_URB
	487	! Variables for multi-layer UCM (Martilli et al. 2002)
	488	REAL, OPTIONAL, INTENT(IN ) :: GMT
	489	INTEGER, OPTIONAL, INTENT(IN ) :: JULDAY
	490	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) ::XLAT, XLONG
	491	INTEGER, INTENT(IN ) :: NUM_URBAN_LAYERS
	492	REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: trb_urb4d
	493	REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tw1_urb4d
	494	REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tw2_urb4d
	495	REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tgb_urb4d
	496	REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tlev_urb3d
	497	REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: qlev_urb3d
	498	REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tw1lev_urb3d
	499	REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tw2lev_urb3d
	500	REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tglev_urb3d
	501	REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: tflev_urb3d
	502	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: lf_ac_urb3d
	503	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: sf_ac_urb3d
	504	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: cm_ac_urb3d
	505	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: sfvent_urb3d
	506	REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: lfvent_urb3d
	507	REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfwin1_urb3d
	508	REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfwin2_urb3d
	509	REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfw1_urb3d
	510	REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfw2_urb3d
	511	REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfr_urb3d
	512	REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_layers, jms:jme ), INTENT(INOUT) :: sfg_urb3d
	513	REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_u_bep !Implicit momemtum component X-direction
	514	REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_v_bep !Implicit momemtum component Y-direction
	515	REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_t_bep !Implicit component pot. temperature
	516	REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_q_bep !Implicit momemtum component X-direction
	517	REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_e_bep !Implicit component TKE
	518	REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::b_u_bep !Explicit momentum component X-direction
	519	REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::b_v_bep !Explicit momentum component Y-direction
	520	REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::b_t_bep !Explicit component pot. temperature
	521	REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::b_q_bep !Implicit momemtum component Y-direction
	522	REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::b_e_bep !Explicit component TKE
	523	REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::vl_bep !Fraction air volume in grid cell
	524	REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::dlg_bep !Height above ground
	525	REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::sf_bep !Fraction air at the face of grid cell
	526	REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::dl_u_bep !Length scale
	527
	528	! Local variables for multi-layer UCM (Martilli et al. 2002)
	529	REAL, DIMENSION( ims:ime, jms:jme ) :: HFX_RURAL,LH_RURAL,GRDFLX_RURAL,RN_RURAL
	530	REAL, DIMENSION( ims:ime, jms:jme ) :: QFX_RURAL,QSFC_RURAL,UMOM_RURAL,VMOM_RURAL
	531	REAL, DIMENSION( ims:ime, jms:jme ) :: ALB_RURAL,EMISS_RURAL,UST_RURAL,TSK_RURAL
	532	! REAL, DIMENSION( ims:ime, jms:jme ) :: GRDFLX_URB
	533	! REAL, DIMENSION( ims:ime, jms:jme ) :: QFX_URB,QSFC_URB,UMOM_URB,VMOM_URB
	534	REAL, DIMENSION( ims:ime, jms:jme ) :: HFX_URB,UMOM_URB,VMOM_URB
	535	REAL, DIMENSION( ims:ime, jms:jme ) :: QFX_URB
	536	! REAL, DIMENSION( ims:ime, jms:jme ) :: ALBEDO_URB,EMISS_URB,UMOM,VMOM,UST
	537	REAL, DIMENSION(ims:ime,jms:jme) ::EMISS_URB
	538	REAL, DIMENSION(ims:ime,jms:jme) :: RL_UP_URB
	539	REAL, DIMENSION(ims:ime,jms:jme) ::RS_ABS_URB
	540	REAL, DIMENSION(ims:ime,jms:jme) ::GRDFLX_URB
	541	REAL :: SIGMA_SB,RL_UP_RURAL,RL_UP_TOT,RS_ABS_TOT,UMOM,VMOM
	542	REAL :: r1,r2,r3
	543	REAL :: CMR_URB, CHR_URB, CMC_URB, CHC_URB
	544	! ----------------------------------------------------------------------
	545	! DECLARATIONS END - urban
	546	! ----------------------------------------------------------------------
	547
	548	REAL, PARAMETER :: CAPA=R_D/CP
	549	REAL :: APELM,APES,SFCTH2,PSFC
	550	real, intent(in) :: xice_threshold
	551	character(len=80) :: message_text
	552
	553	! MEK MAY 2007
	554	FDTLIW=DT/ROWLIW
	555	! MEK JUL2007
	556	FDTW=DT/(XLV*RHOWATER)
	557	! debug printout
	558	IPRINT=.false.
	559
	560	! SLOPETYP=2
	561	SLOPETYP=1
	562	! SHDMIN=0.00
	563
	564
	565	NSOIL=num_soil_layers
	566
	567	DO NS=1,NSOIL
	568	SLDPTH(NS)=DZS(NS)
	569	ENDDO
	570
	571	DO J=jts,jte
	572
	573	IF(ITIMESTEP.EQ.1)THEN
	574	DO 50 I=its,ite
	575	!*** initialize soil conditions for IHOP 31 May case
	576	! IF((XLAND(I,J)-1.5) < 0.)THEN
	577	! if (I==108.and.j==85) then
	578	! DO NS=1,NSOIL
	579	! SMOIS(I,NS,J)=0.10
	580	! SH2O(I,NS,J)=0.10
	581	! enddo
	582	! endif
	583	! ENDIF
	584
	585	!*** SET ZERO-VALUE FOR SOME OUTPUT DIAGNOSTIC ARRAYS
	586	IF((XLAND(I,J)-1.5).GE.0.)THEN
	587	! check sea-ice point
	588	#if 0
	589	IF( XICE(I,J).GE. XICE_THRESHOLD .and. IPRINT ) PRINT*, ' sea-ice at water point, I=',I,'J=',J
	590	#endif
	591	!*** Open Water Case
	592	SMSTAV(I,J)=1.0
	593	SMSTOT(I,J)=1.0
	594	DO NS=1,NSOIL
	595	SMOIS(I,NS,J)=1.0
	596	TSLB(I,NS,J)=273.16 !STEMP
	597	SMCREL(I,NS,J)=1.0
	598	ENDDO
	599	ELSE
	600	IF ( XICE(I,J) .GE. XICE_THRESHOLD ) THEN
	601	!*** SEA-ICE CASE
	602	SMSTAV(I,J)=1.0
	603	SMSTOT(I,J)=1.0
	604	DO NS=1,NSOIL
	605	SMOIS(I,NS,J)=1.0
	606	SMCREL(I,NS,J)=1.0
	607	ENDDO
	608	ENDIF
	609	ENDIF
	610	!
	611	50 CONTINUE
	612	ENDIF ! end of initialization over ocean
	613
	614	!-----------------------------------------------------------------------
	615	DO 100 I=its,ite
	616	! surface pressure
	617	PSFC=P8w3D(i,1,j)
	618	! pressure in middle of lowest layer
	619	SFCPRS=(P8W3D(I,KTS+1,j)+P8W3D(i,KTS,j))*0.5
	620	! convert from mixing ratio to specific humidity
	621	Q2K=QV3D(i,1,j)/(1.0+QV3D(i,1,j))
	622	!
	623	! Q2SAT=QGH(I,j)
	624	Q2SAT=QGH(I,J)/(1.0+QGH(I,J)) ! Q2SAT is sp humidity
	625	! add check on myj=.true.
	626	! IF((Q2K.GE.Q2SAT*TRESH).AND.Q2K.LT.QZ0(I,J))THEN
	627	IF((myj).AND.(Q2K.GE.Q2SAT*TRESH).AND.Q2K.LT.QZ0(I,J))THEN
	628	SATFLG=0.
	629	CHKLOWQ(I,J)=0.
	630	ELSE
	631	SATFLG=1.0
	632	CHKLOWQ(I,J)=1.
	633	ENDIF
	634
	635	SFCTMP=T3D(i,1,j)
	636	ZLVL=0.5*DZ8W(i,1,j)
	637
	638	! TH2=SFCTMP+(0.0097545*ZLVL)
	639	! calculate SFCTH2 via Exner function vs lapse-rate (above)
	640	APES=(1.E5/PSFC)**CAPA
	641	APELM=(1.E5/SFCPRS)**CAPA
	642	SFCTH2=SFCTMP*APELM
	643	TH2=SFCTH2/APES
	644	!
	645	EMISSI = EMISS(I,J)
	646	LWDN=GLW(I,J)*EMISSI
	647	! SOLDN is total incoming solar
	648	SOLDN=SWDOWN(I,J)
	649	! GSW is net downward solar
	650	! SOLNET=GSW(I,J)
	651	! use mid-day albedo to determine net downward solar (no solar zenith angle correction)
	652	SOLNET=SOLDN*(1.-ALBEDO(I,J))
	653	PRCP=RAINBL(i,j)/DT
	654	VEGTYP=IVGTYP(I,J)
	655	SOILTYP=ISLTYP(I,J)
	656	SHDFAC=VEGFRA(I,J)/100.
	657	T1=TSK(I,J)
	658	CHK=CHS(I,J)
	659	SHMIN=SHDMIN(I,J)/100. !NEW
	660	SHMAX=SHDMAX(I,J)/100. !NEW
	661	! convert snow water equivalent from mm to meter
	662	SNEQV=SNOW(I,J)*0.001
	663	! snow depth in meters
	664	SNOWHK=SNOWH(I,J)
	665	SNCOVR=SNOWC(I,J)
	666
	667	! if "SR" present, set frac of frozen precip ("FFROZP") = snow-ratio ("SR", range:0-1)
	668	! SR from e.g. Ferrier microphysics
	669	! otherwise define from 1st atmos level temperature
	670	IF(FRPCPN) THEN
	671	FFROZP=SR(I,J)
	672	ELSE
	673	IF (SFCTMP <= 273.15) THEN
	674	FFROZP = 1.0
	675	ELSE
	676	FFROZP = 0.0
	677	ENDIF
	678	ENDIF
	679	!***
	680	IF((XLAND(I,J)-1.5).GE.0.)THEN ! begining of land/sea if block
	681	! Open water points
	682	TSK_RURAL(I,J)=TSK(I,J)
	683	HFX_RURAL(I,J)=HFX(I,J)
	684	QFX_RURAL(I,J)=QFX(I,J)
	685	LH_RURAL(I,J)=LH(I,J)
	686	EMISS_RURAL(I,J)=EMISS(I,J)
	687	GRDFLX_RURAL(I,J)=GRDFLX(I,J)
	688	ELSE
	689	! Land or sea-ice case
	690
	691	IF (XICE(I,J) >= XICE_THRESHOLD) THEN
	692	! Sea-ice point
	693	ICE = 1
	694	ELSE IF ( VEGTYP == ISICE ) THEN
	695	! Land-ice point
	696	ICE = -1
	697	ELSE
	698	! Neither sea ice or land ice.
	699	ICE=0
	700	ENDIF
	701	DQSDT2=Q2SATA23M4/(SFCTMP-A4)*2
	702
	703	IF(SNOW(I,J).GT.0.0)THEN
	704	! snow on surface (use ice saturation properties)
	705	SFCTSNO=SFCTMP
	706	E2SAT=611.2EXP(6174.(1./273.15 - 1./SFCTSNO))
	707	Q2SATI=0.622*E2SAT/(SFCPRS-E2SAT)
	708	Q2SATI=Q2SATI/(1.0+Q2SATI) ! spec. hum.
	709	IF (T1 .GT. 273.14) THEN
	710	! warm ground temps, weight the saturation between ice and water according to SNOWC
	711	Q2SAT=Q2SAT(1.-SNOWC(I,J)) + Q2SATISNOWC(I,J)
	712	DQSDT2=DQSDT2(1.-SNOWC(I,J)) + Q2SATI6174./(SFCTSNO*2)SNOWC(I,J)
	713	ELSE
	714	! cold ground temps, use ice saturation only
	715	Q2SAT=Q2SATI
	716	DQSDT2=Q2SATI6174./(SFCTSNO*2)
	717	ENDIF
	718	! for snow cover fraction at 0 C, ground temp will not change, so DQSDT2 effectively zero
	719	IF(T1 .GT. 273. .AND. SNOWC(I,J) .GT. 0.)DQSDT2=DQSDT2*(1.-SNOWC(I,J))
	720	ENDIF
	721
	722	IF(ICE.EQ.1)THEN
	723	! Sea-ice point has deep-level temperature of -2 C
	724	TBOT=271.16
	725	ELSE
	726	! Land-ice or land points have the usual deep-soil temperature.
	727	TBOT=TMN(I,J)
	728	ENDIF
	729	IF(VEGTYP.EQ.25) SHDFAC=0.0000
	730	IF(VEGTYP.EQ.26) SHDFAC=0.0000
	731	IF(VEGTYP.EQ.27) SHDFAC=0.0000
	732	IF(SOILTYP.EQ.14.AND.XICE(I,J).EQ.0.)THEN
	733	#if 0
	734	IF(IPRINT)PRINT*,' SOIL TYPE FOUND TO BE WATER AT A LAND-POINT'
	735	IF(IPRINT)PRINT*,i,j,'RESET SOIL in surfce.F'
	736	#endif
	737	SOILTYP=7
	738	ENDIF
	739	SNOALB1 = SNOALB(I,J)
	740	CMC=CANWAT(I,J)
	741
	742	!-------------------------------------------
	743	!*** convert snow depth from mm to meter
	744	!
	745	! IF(RDMAXALB) THEN
	746	! SNOALB=ALBMAX(I,J)*0.01
	747	! ELSE
	748	! SNOALB=MAXALB(IVGTPK)*0.01
	749	! ENDIF
	750
	751	! SNOALB1=0.80
	752	! SHMIN=0.00
	753	ALBBRD=ALBBCK(I,J)
	754	Z0BRD=Z0(I,J)
	755	EMBRD=EMBCK(I,J)
	756	SNOTIME1 = SNOTIME(I,J)
	757	RIBB=RIB(I,J)
	758	!FEI: temporaray arrays above need to be changed later by using SI
	759
	760	DO 70 NS=1,NSOIL
	761	SMC(NS)=SMOIS(I,NS,J)
	762	STC(NS)=TSLB(I,NS,J) !STEMP
	763	SWC(NS)=SH2O(I,NS,J)
	764	70 CONTINUE
	765	!
	766	if ( (SNEQV.ne.0..AND.SNOWHK.eq.0.).or.(SNOWHK.le.SNEQV) )THEN
	767	SNOWHK= 5.*SNEQV
	768	endif
	769	!
	770
	771	!Fei: urban. for urban surface, if calling UCM, redefine the natural surface in cities as
	772	! the "NATURAL" category in the VEGPARM.TBL
	773
	774	IF(SF_URBAN_PHYSICS == 1.OR. SF_URBAN_PHYSICS==2.OR.SF_URBAN_PHYSICS==3 ) THEN
	775	IF( IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == 31 .or. &
	776	IVGTYP(I,J) == 32 .or. IVGTYP(I,J) == 33) THEN
	777	VEGTYP = NATURAL
	778	SHDFAC = SHDTBL(NATURAL)
	779	ALBEDOK =0.2 ! 0.2
	780	ALBBRD =0.2 !0.2
	781	EMISSI = 0.98 !for VEGTYP=5
	782	IF ( FRC_URB2D(I,J) < 0.99 ) THEN
	783	if(sf_urban_physics.eq.1)then
	784	T1= ( TSK(I,J) -FRC_URB2D(I,J) * TS_URB2D (I,J) )/ (1-FRC_URB2D(I,J))
	785	elseif((sf_urban_physics.eq.2).OR.(sf_urban_physics.eq.3))then
	786	r1= (tsk(i,j)**4.)
	787	r2= frc_urb2d(i,j)(ts_urb2d(i,j)*4.)
	788	r3= (1.-frc_urb2d(i,j))
	789	t1= ((r1-r2)/r3)**.25
	790	endif
	791	ELSE
	792	T1 = TSK(I,J)
	793	ENDIF
	794	ENDIF
	795	ELSE
	796	IF( IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == 31 .or. &
	797	IVGTYP(I,J) == 32 .or. IVGTYP(I,J) == 33) THEN
	798	VEGTYP = ISURBAN
	799	ENDIF
	800	ENDIF
	801
	802	#if 0
	803	IF(IPRINT) THEN
	804	!
	805	print*, 'BEFORE SFLX, in Noahlsm_driver'
	806	print*, 'ICE', ICE, 'DT',DT, 'ZLVL',ZLVL, 'NSOIL', NSOIL, &
	807	'SLDPTH', SLDPTH, 'LOCAL',LOCAL, 'LUTYPE',&
	808	LUTYPE, 'SLTYPE',SLTYPE, 'LWDN',LWDN, 'SOLDN',SOLDN, &
	809	'SFCPRS',SFCPRS, 'PRCP',PRCP,'SFCTMP',SFCTMP,'Q2K',Q2K, &
	810	'TH2',TH2,'Q2SAT',Q2SAT,'DQSDT2',DQSDT2,'VEGTYP', VEGTYP,&
	811	'SOILTYP',SOILTYP, 'SLOPETYP',SLOPETYP, 'SHDFAC',SHDFAC,&
	812	'SHMIN',SHMIN, 'ALBBRD',ALBBRD,'SNOALB1',SNOALB1,'TBOT',&
	813	TBOT, 'Z0BRD',Z0BRD, 'Z0K',Z0K, 'CMC',CMC, 'T1',T1,'STC',&
	814	STC, 'SMC',SMC, 'SWC',SWC,'SNOWHK',SNOWHK,'SNEQV',SNEQV,&
	815	'ALBEDOK',ALBEDOK,'CHK',CHK,'ETA',ETA,'SHEAT',SHEAT, &
	816	'ETA_KINEMATIC',ETA_KINEMATIC, 'FDOWN',FDOWN,'EC',EC, &
	817	'EDIR',EDIR,'ET',ET,'ETT',ETT,'ESNOW',ESNOW,'DRIP',DRIP,&
	818	'DEW',DEW,'BETA',BETA,'ETP',ETP,'SSOIL',SSOIL,'FLX1',FLX1,&
	819	'FLX2',FLX2,'FLX3',FLX3,'SNOMLT',SNOMLT,'SNCOVR',SNCOVR,&
	820	'RUNOFF1',RUNOFF1,'RUNOFF2',RUNOFF2,'RUNOFF3',RUNOFF3, &
	821	'RC',RC, 'PC',PC,'RSMIN',RSMIN,'XLAI',XLAI,'RCS',RCS, &
	822	'RCT',RCT,'RCQ',RCQ,'RCSOIL',RCSOIL,'SOILW',SOILW, &
	823	'SOILM',SOILM,'Q1',Q1,'SMCWLT',SMCWLT,'SMCDRY',SMCDRY,&
	824	'SMCREF',SMCREF,'SMCMAX',SMCMAX,'NROOT',NROOT
	825	endif
	826	#endif
	827
	828
	829	IF (rdlai2d) THEN
	830	xlai = lai(i,j)
	831	endif
	832
	833	CALL SFLX (FFROZP, ICE, ISURBAN, DT,ZLVL,NSOIL,SLDPTH, & !C
	834	LOCAL, & !L
	835	LUTYPE, SLTYPE, & !CL
	836	LWDN,SOLDN,SOLNET,SFCPRS,PRCP,SFCTMP,Q2K,DUMMY, & !F
	837	DUMMY,DUMMY, DUMMY, & !F PRCPRAIN not used
	838	TH2,Q2SAT,DQSDT2, & !I
	839	VEGTYP,SOILTYP,SLOPETYP,SHDFAC,SHMIN,SHMAX, & !I
	840	ALBBRD, SNOALB1,TBOT, Z0BRD, Z0K, EMISSI, EMBRD, & !S
	841	CMC,T1,STC,SMC,SWC,SNOWHK,SNEQV,ALBEDOK,CHK,dummy,& !H
	842	ETA,SHEAT, ETA_KINEMATIC,FDOWN, & !O
	843	EC,EDIR,ET,ETT,ESNOW,DRIP,DEW, & !O
	844	BETA,ETP,SSOIL, & !O
	845	FLX1,FLX2,FLX3, & !O
	846	SNOMLT,SNCOVR, & !O
	847	RUNOFF1,RUNOFF2,RUNOFF3, & !O
	848	RC,PC,RSMIN,XLAI,RCS,RCT,RCQ,RCSOIL, & !O
	849	SOILW,SOILM,Q1,SMAV, & !D
	850	RDLAI2D,USEMONALB, &
	851	SNOTIME1, &
	852	RIBB, &
	853	SMCWLT,SMCDRY,SMCREF,SMCMAX,NROOT)
	854
	855
	856	lai(i,j) = xlai
	857
	858	#if 0
	859	IF(IPRINT) THEN
	860
	861	print*, 'AFTER SFLX, in Noahlsm_driver'
	862	print*, 'ICE', ICE, 'DT',DT, 'ZLVL',ZLVL, 'NSOIL', NSOIL, &
	863	'SLDPTH', SLDPTH, 'LOCAL',LOCAL, 'LUTYPE',&
	864	LUTYPE, 'SLTYPE',SLTYPE, 'LWDN',LWDN, 'SOLDN',SOLDN, &
	865	'SFCPRS',SFCPRS, 'PRCP',PRCP,'SFCTMP',SFCTMP,'Q2K',Q2K, &
	866	'TH2',TH2,'Q2SAT',Q2SAT,'DQSDT2',DQSDT2,'VEGTYP', VEGTYP,&
	867	'SOILTYP',SOILTYP, 'SLOPETYP',SLOPETYP, 'SHDFAC',SHDFAC,&
	868	'SHDMIN',SHMIN, 'ALBBRD',ALBBRD,'SNOALB',SNOALB1,'TBOT',&
	869	TBOT, 'Z0BRD',Z0BRD, 'Z0K',Z0K, 'CMC',CMC, 'T1',T1,'STC',&
	870	STC, 'SMC',SMC, 'SWc',SWC,'SNOWHK',SNOWHK,'SNEQV',SNEQV,&
	871	'ALBEDOK',ALBEDOK,'CHK',CHK,'ETA',ETA,'SHEAT',SHEAT, &
	872	'ETA_KINEMATIC',ETA_KINEMATIC, 'FDOWN',FDOWN,'EC',EC, &
	873	'EDIR',EDIR,'ET',ET,'ETT',ETT,'ESNOW',ESNOW,'DRIP',DRIP,&
	874	'DEW',DEW,'BETA',BETA,'ETP',ETP,'SSOIL',SSOIL,'FLX1',FLX1,&
	875	'FLX2',FLX2,'FLX3',FLX3,'SNOMLT',SNOMLT,'SNCOVR',SNCOVR,&
	876	'RUNOFF1',RUNOFF1,'RUNOFF2',RUNOFF2,'RUNOFF3',RUNOFF3, &
	877	'RC',RC, 'PC',PC,'RSMIN',RSMIN,'XLAI',XLAI,'RCS',RCS, &
	878	'RCT',RCT,'RCQ',RCQ,'RCSOIL',RCSOIL,'SOILW',SOILW, &
	879	'SOILM',SOILM,'Q1',Q1,'SMCWLT',SMCWLT,'SMCDRY',SMCDRY,&
	880	'SMCREF',SMCREF,'SMCMAX',SMCMAX,'NROOT',NROOT
	881	endif
	882	#endif
	883
	884	!*** UPDATE STATE VARIABLES
	885	CANWAT(I,J)=CMC
	886	SNOW(I,J)=SNEQV*1000.
	887	! SNOWH(I,J)=SNOWHK*1000.
	888	SNOWH(I,J)=SNOWHK ! SNOWHK in meters
	889	ALBEDO(I,J)=ALBEDOK
	890	ALB_RURAL(I,J)=ALBEDOK
	891	ALBBCK(I,J)=ALBBRD
	892	Z0(I,J)=Z0BRD
	893	EMISS(I,J) = EMISSI
	894	EMISS_RURAL(I,J) = EMISSI
	895	! Noah: activate time-varying roughness length (V3.3 Feb 2011)
	896	ZNT(I,J)=Z0K
	897	TSK(I,J)=T1
	898	TSK_RURAL(I,J)=T1
	899	HFX(I,J)=SHEAT
	900	HFX_RURAL(I,J)=SHEAT
	901	! MEk Jul07 add potential evap accum
	902	POTEVP(I,J)=POTEVP(I,J)+ETP*FDTW
	903	QFX(I,J)=ETA_KINEMATIC
	904	QFX_RURAL(I,J)=ETA_KINEMATIC
	905	LH(I,J)=ETA
	906	LH_RURAL(I,J)=ETA
	907	GRDFLX(I,J)=SSOIL
	908	GRDFLX_RURAL(I,J)=SSOIL
	909	SNOWC(I,J)=SNCOVR
	910	CHS2(I,J)=CQS2(I,J)
	911	SNOTIME(I,J) = SNOTIME1
	912	! prevent diagnostic ground q (q1) from being greater than qsat(tsk)
	913	! as happens over snow cover where the cqs2 value also becomes irrelevant
	914	! by setting cqs2=chs in this situation the 2m q should become just qv(k=1)
	915	IF (Q1 .GT. QSFC(I,J)) THEN
	916	CQS2(I,J) = CHS(I,J)
	917	ENDIF
	918	! QSFC(I,J)=Q1
	919	! Convert QSFC back to mixing ratio
	920	QSFC(I,J)= Q1/(1.0-Q1)
	921	!
	922	QSFC_RURAL(I,J)= Q1/(1.0-Q1)
	923	! Calculate momentum flux from rural surface for use with multi-layer UCM (Martilli et al. 2002)
	924
	925	DO 80 NS=1,NSOIL
	926	SMOIS(I,NS,J)=SMC(NS)
	927	TSLB(I,NS,J)=STC(NS) ! STEMP
	928	SH2O(I,NS,J)=SWC(NS)
	929	80 CONTINUE
	930	! ENDIF
	931
	932	!
	933	! Residual of surface energy balance equation terms
	934	!
	935	noahres(i,j) = ( solnet + lwdn ) - sheat + ssoil - eta - ( emissi * STBOLT * (t1**4) ) - flx1 - flx2 - flx3
	936
	937
	938	IF (SF_URBAN_PHYSICS == 1 ) THEN ! Beginning of UCM CALL if block
	939	!--------------------------------------
	940	! URBAN CANOPY MODEL START - urban
	941	!--------------------------------------
	942	! Input variables lsm --> urban
	943
	944
	945	IF( IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == 31 .or. &
	946	IVGTYP(I,J) == 32 .or. IVGTYP(I,J) == 33 ) THEN
	947
	948	! Call urban
	949
	950	!
	951	UTYPE_URB = UTYPE_URB2D(I,J) !urban type (low, high or industrial)
	952
	953	TA_URB = SFCTMP ! [K]
	954	QA_URB = Q2K ! [kg/kg]
	955	UA_URB = SQRT(U_PHY(I,1,J)2.+V_PHY(I,1,J)2.)
	956	U1_URB = U_PHY(I,1,J)
	957	V1_URB = V_PHY(I,1,J)
	958	IF(UA_URB < 1.) UA_URB=1. ! [m/s]
	959	SSG_URB = SOLDN ! [W/m/m]
	960	SSGD_URB = 0.8*SOLDN ! [W/m/m]
	961	SSGQ_URB = SSG_URB-SSGD_URB ! [W/m/m]
	962	LLG_URB = GLW(I,J) ! [W/m/m]
	963	RAIN_URB = RAINBL(I,J) ! [mm]
	964	RHOO_URB = SFCPRS / (287.04 * SFCTMP * (1.0+ 0.61 * Q2K)) ![kg/m/m/m]
	965	ZA_URB = ZLVL ! [m]
	966	DELT_URB = DT ! [sec]
	967	XLAT_URB = XLAT_URB2D(I,J) ! [deg]
	968	COSZ_URB = COSZ_URB2D(I,J) !
	969	OMG_URB = OMG_URB2D(I,J) !
	970	ZNT_URB = ZNT(I,J)
	971
	972	LSOLAR_URB = .FALSE.
	973
	974	TR_URB = TR_URB2D(I,J)
	975	TB_URB = TB_URB2D(I,J)
	976	TG_URB = TG_URB2D(I,J)
	977	TC_URB = TC_URB2D(I,J)
	978	QC_URB = QC_URB2D(I,J)
	979	UC_URB = UC_URB2D(I,J)
	980
	981	DO K = 1,num_roof_layers
	982	TRL_URB(K) = TRL_URB3D(I,K,J)
	983	END DO
	984	DO K = 1,num_wall_layers
	985	TBL_URB(K) = TBL_URB3D(I,K,J)
	986	END DO
	987	DO K = 1,num_road_layers
	988	TGL_URB(K) = TGL_URB3D(I,K,J)
	989	END DO
	990
	991	XXXR_URB = XXXR_URB2D(I,J)
	992	XXXB_URB = XXXB_URB2D(I,J)
	993	XXXG_URB = XXXG_URB2D(I,J)
	994	XXXC_URB = XXXC_URB2D(I,J)
	995	!
	996	!
	997	! Limits to avoid dividing by small number
	998	if (CHS(I,J) < 1.0E-02) then
	999	CHS(I,J) = 1.0E-02
	1000	endif
	1001	if (CHS2(I,J) < 1.0E-02) then
	1002	CHS2(I,J) = 1.0E-02
	1003	endif
	1004	if (CQS2(I,J) < 1.0E-02) then
	1005	CQS2(I,J) = 1.0E-02
	1006	endif
	1007	!
	1008	CHS_URB = CHS(I,J)
	1009	CHS2_URB = CHS2(I,J)
	1010	IF (PRESENT(CMR_SFCDIF)) THEN
	1011	CMR_URB = CMR_SFCDIF(I,J)
	1012	CHR_URB = CHR_SFCDIF(I,J)
	1013	CMC_URB = CMC_SFCDIF(I,J)
	1014	CHC_URB = CHC_SFCDIF(I,J)
	1015	ENDIF
	1016	!
	1017	! Call urban
	1018
	1019	CALL urban(LSOLAR_URB, & ! I
	1020	num_roof_layers,num_wall_layers,num_road_layers, & ! C
	1021	DZR,DZB,DZG, & ! C
	1022	UTYPE_URB,TA_URB,QA_URB,UA_URB,U1_URB,V1_URB,SSG_URB, & ! I
	1023	SSGD_URB,SSGQ_URB,LLG_URB,RAIN_URB,RHOO_URB, & ! I
	1024	ZA_URB,DECLIN_URB,COSZ_URB,OMG_URB, & ! I
	1025	XLAT_URB,DELT_URB,ZNT_URB, & ! I
	1026	CHS_URB, CHS2_URB, & ! I
	1027	TR_URB, TB_URB, TG_URB, TC_URB, QC_URB,UC_URB, & ! H
	1028	TRL_URB,TBL_URB,TGL_URB, & ! H
	1029	XXXR_URB, XXXB_URB, XXXG_URB, XXXC_URB, & ! H
	1030	TS_URB,QS_URB,SH_URB,LH_URB,LH_KINEMATIC_URB, & ! O
	1031	SW_URB,ALB_URB,LW_URB,G_URB,RN_URB,PSIM_URB,PSIH_URB, & ! O
	1032	GZ1OZ0_URB, & !O
	1033	CMR_URB, CHR_URB, CMC_URB, CHC_URB, &
	1034	U10_URB, V10_URB, TH2_URB, Q2_URB, & ! O
	1035	UST_URB) !O
	1036
	1037	#if 0
	1038	IF(IPRINT) THEN
	1039
	1040	print*, 'AFTER CALL URBAN'
	1041	print*,'num_roof_layers',num_roof_layers, 'num_wall_layers', &
	1042	num_wall_layers, &
	1043	'DZR',DZR,'DZB',DZB,'DZG',DZG,'UTYPE_URB',UTYPE_URB,'TA_URB', &
	1044	TA_URB, &
	1045	'QA_URB',QA_URB,'UA_URB',UA_URB,'U1_URB',U1_URB,'V1_URB', &
	1046	V1_URB, &
	1047	'SSG_URB',SSG_URB,'SSGD_URB',SSGD_URB,'SSGQ_URB',SSGQ_URB, &
	1048	'LLG_URB',LLG_URB,'RAIN_URB',RAIN_URB,'RHOO_URB',RHOO_URB, &
	1049	'ZA_URB',ZA_URB, 'DECLIN_URB',DECLIN_URB,'COSZ_URB',COSZ_URB,&
	1050	'OMG_URB',OMG_URB,'XLAT_URB',XLAT_URB,'DELT_URB',DELT_URB, &
	1051	'ZNT_URB',ZNT_URB,'TR_URB',TR_URB, 'TB_URB',TB_URB,'TG_URB',&
	1052	TG_URB,'TC_URB',TC_URB,'QC_URB',QC_URB,'TRL_URB',TRL_URB, &
	1053	'TBL_URB',TBL_URB,'TGL_URB',TGL_URB,'XXXR_URB',XXXR_URB, &
	1054	'XXXB_URB',XXXB_URB,'XXXG_URB',XXXG_URB,'XXXC_URB',XXXC_URB,&
	1055	'TS_URB',TS_URB,'QS_URB',QS_URB,'SH_URB',SH_URB,'LH_URB', &
	1056	LH_URB, 'LH_KINEMATIC_URB',LH_KINEMATIC_URB,'SW_URB',SW_URB,&
	1057	'ALB_URB',ALB_URB,'LW_URB',LW_URB,'G_URB',G_URB,'RN_URB', &
	1058	RN_URB, 'PSIM_URB',PSIM_URB,'PSIH_URB',PSIH_URB, &
	1059	'U10_URB',U10_URB,'V10_URB',V10_URB,'TH2_URB',TH2_URB, &
	1060	'Q2_URB',Q2_URB,'CHS_URB',CHS_URB,'CHS2_URB',CHS2_URB
	1061	endif
	1062	#endif
	1063
	1064	TS_URB2D(I,J) = TS_URB
	1065
	1066	ALBEDO(I,J) = FRC_URB2D(I,J)ALB_URB+(1-FRC_URB2D(I,J))ALBEDOK ![-]
	1067	HFX(I,J) = FRC_URB2D(I,J)SH_URB+(1-FRC_URB2D(I,J))SHEAT ![W/m/m]
	1068	QFX(I,J) = FRC_URB2D(I,J)*LH_KINEMATIC_URB &
	1069	+ (1-FRC_URB2D(I,J))*ETA_KINEMATIC ![kg/m/m/s]
	1070	LH(I,J) = FRC_URB2D(I,J)LH_URB+(1-FRC_URB2D(I,J))ETA ![W/m/m]
	1071	GRDFLX(I,J) = FRC_URB2D(I,J)G_URB+(1-FRC_URB2D(I,J))SSOIL ![W/m/m]
	1072	TSK(I,J) = FRC_URB2D(I,J)TS_URB+(1-FRC_URB2D(I,J))T1 ![K]
	1073	Q1 = FRC_URB2D(I,J)QS_URB+(1-FRC_URB2D(I,J))Q1 ![-]
	1074	! Convert QSFC back to mixing ratio
	1075	QSFC(I,J)= Q1/(1.0-Q1)
	1076	UST(I,J)= FRC_URB2D(I,J)UST_URB+(1-FRC_URB2D(I,J))UST(I,J) ![m/s]
	1077
	1078	#if 0
	1079	IF(IPRINT)THEN
	1080
	1081	print*, ' FRC_URB2D', FRC_URB2D, &
	1082	'ALB_URB',ALB_URB, 'ALBEDOK',ALBEDOK, &
	1083	'ALBEDO(I,J)', ALBEDO(I,J), &
	1084	'SH_URB',SH_URB,'SHEAT',SHEAT, 'HFX(I,J)',HFX(I,J), &
	1085	'LH_KINEMATIC_URB',LH_KINEMATIC_URB,'ETA_KINEMATIC', &
	1086	ETA_KINEMATIC, 'QFX(I,J)',QFX(I,J), &
	1087	'LH_URB',LH_URB, 'ETA',ETA, 'LH(I,J)',LH(I,J), &
	1088	'G_URB',G_URB,'SSOIL',SSOIL,'GRDFLX(I,J)', GRDFLX(I,J),&
	1089	'TS_URB',TS_URB,'T1',T1,'TSK(I,J)',TSK(I,J), &
	1090	'QS_URB',QS_URB,'Q1',Q1,'QSFC(I,J)',QSFC(I,J)
	1091	endif
	1092	#endif
	1093
	1094
	1095
	1096	! Renew Urban State Varialbes
	1097
	1098	TR_URB2D(I,J) = TR_URB
	1099	TB_URB2D(I,J) = TB_URB
	1100	TG_URB2D(I,J) = TG_URB
	1101	TC_URB2D(I,J) = TC_URB
	1102	QC_URB2D(I,J) = QC_URB
	1103	UC_URB2D(I,J) = UC_URB
	1104
	1105	DO K = 1,num_roof_layers
	1106	TRL_URB3D(I,K,J) = TRL_URB(K)
	1107	END DO
	1108	DO K = 1,num_wall_layers
	1109	TBL_URB3D(I,K,J) = TBL_URB(K)
	1110	END DO
	1111	DO K = 1,num_road_layers
	1112	TGL_URB3D(I,K,J) = TGL_URB(K)
	1113	END DO
	1114	XXXR_URB2D(I,J) = XXXR_URB
	1115	XXXB_URB2D(I,J) = XXXB_URB
	1116	XXXG_URB2D(I,J) = XXXG_URB
	1117	XXXC_URB2D(I,J) = XXXC_URB
	1118
	1119	SH_URB2D(I,J) = SH_URB
	1120	LH_URB2D(I,J) = LH_URB
	1121	G_URB2D(I,J) = G_URB
	1122	RN_URB2D(I,J) = RN_URB
	1123	PSIM_URB2D(I,J) = PSIM_URB
	1124	PSIH_URB2D(I,J) = PSIH_URB
	1125	GZ1OZ0_URB2D(I,J)= GZ1OZ0_URB
	1126	U10_URB2D(I,J) = U10_URB
	1127	V10_URB2D(I,J) = V10_URB
	1128	TH2_URB2D(I,J) = TH2_URB
	1129	Q2_URB2D(I,J) = Q2_URB
	1130	UST_URB2D(I,J) = UST_URB
	1131	AKMS_URB2D(I,J) = KARMAN * UST_URB2D(I,J)/(GZ1OZ0_URB2D(I,J)-PSIM_URB2D(I,J))
	1132	IF (PRESENT(CMR_SFCDIF)) THEN
	1133	CMR_SFCDIF(I,J) = CMR_URB
	1134	CHR_SFCDIF(I,J) = CHR_URB
	1135	CMC_SFCDIF(I,J) = CMC_URB
	1136	CHC_SFCDIF(I,J) = CHC_URB
	1137	ENDIF
	1138	END IF
	1139
	1140	ENDIF ! end of UCM CALL if block
	1141	!--------------------------------------
	1142	! Urban Part End - urban
	1143	!--------------------------------------
	1144
	1145	!*** DIAGNOSTICS
	1146	SMSTAV(I,J)=SOILW
	1147	SMSTOT(I,J)=SOILM*1000.
	1148	DO NS=1,NSOIL
	1149	SMCREL(I,NS,J)=SMAV(NS)
	1150	ENDDO
	1151	! Convert the water unit into mm
	1152	SFCRUNOFF(I,J)=SFCRUNOFF(I,J)+RUNOFF1DT1000.0
	1153	UDRUNOFF(I,J)=UDRUNOFF(I,J)+RUNOFF2DT1000.0
	1154	! snow defined when fraction of frozen precip (FFROZP) > 0.5,
	1155	IF(FFROZP.GT.0.5)THEN
	1156	ACSNOW(I,J)=ACSNOW(I,J)+PRCP*DT
	1157	ENDIF
	1158	IF(SNOW(I,J).GT.0.)THEN
	1159	ACSNOM(I,J)=ACSNOM(I,J)+SNOMLT*1000.
	1160	! accumulated snow-melt energy
	1161	SNOPCX(I,J)=SNOPCX(I,J)-SNOMLT/FDTLIW
	1162	ENDIF
	1163
	1164	ENDIF ! endif of land-sea test
	1165
	1166	100 CONTINUE ! of I loop
	1167
	1168	ENDDO ! of J loop
	1169
	1170	IF (SF_URBAN_PHYSICS == 2) THEN
	1171
	1172
	1173	do j=jts,jte
	1174	do i=its,ite
	1175	EMISS_URB(i,j)=0.
	1176	RL_UP_URB(i,j)=0.
	1177	RS_ABS_URB(i,j)=0.
	1178	GRDFLX_URB(i,j)=0.
	1179	b_q_bep(i,kts:kte,j)=0.
	1180	end do
	1181	end do
	1182	CALL BEP(frc_urb2d,utype_urb2d,itimestep,dz8w,dt,u_phy,v_phy, &
	1183	th_phy,rho,p_phy,swdown,glw, &
	1184	gmt,julday,xlong,xlat,declin_urb,cosz_urb2d,omg_urb2d, &
	1185	num_urban_layers, &
	1186	trb_urb4d,tw1_urb4d,tw2_urb4d,tgb_urb4d, &
	1187	sfw1_urb3d,sfw2_urb3d,sfr_urb3d,sfg_urb3d, &
	1188	a_u_bep,a_v_bep,a_t_bep, &
	1189	a_e_bep,b_u_bep,b_v_bep, &
	1190	b_t_bep,b_e_bep,dlg_bep, &
	1191	dl_u_bep,sf_bep,vl_bep, &
	1192	rl_up_urb,rs_abs_urb,emiss_urb,grdflx_urb, &
	1193	ids,ide, jds,jde, kds,kde, &
	1194	ims,ime, jms,jme, kms,kme, &
	1195	its,ite, jts,jte, kts,kte )
	1196
	1197	ENDIF
	1198
	1199
	1200	IF (SF_URBAN_PHYSICS == 3) THEN
	1201
	1202
	1203	do j=jts,jte
	1204	do i=its,ite
	1205	EMISS_URB(i,j)=0.
	1206	RL_UP_URB(i,j)=0.
	1207	RS_ABS_URB(i,j)=0.
	1208	GRDFLX_URB(i,j)=0.
	1209	b_q_bep(i,kts:kte,j)=0.
	1210	end do
	1211	end do
	1212
	1213	CALL BEP_BEM(frc_urb2d,utype_urb2d,itimestep,dz8w,dt,u_phy,v_phy, &
	1214	th_phy,rho,p_phy,swdown,glw, &
	1215	gmt,julday,xlong,xlat,declin_urb,cosz_urb2d,omg_urb2d, &
	1216	num_urban_layers, &
	1217	trb_urb4d,tw1_urb4d,tw2_urb4d,tgb_urb4d, &
	1218	tlev_urb3d,qlev_urb3d,tw1lev_urb3d,tw2lev_urb3d, &
	1219	tglev_urb3d,tflev_urb3d,sf_ac_urb3d,lf_ac_urb3d, &
	1220	cm_ac_urb3d,sfvent_urb3d,lfvent_urb3d, &
	1221	sfwin1_urb3d,sfwin2_urb3d, &
	1222	sfw1_urb3d,sfw2_urb3d,sfr_urb3d,sfg_urb3d, &
	1223	a_u_bep,a_v_bep,a_t_bep, &
	1224	a_e_bep,b_u_bep,b_v_bep, &
	1225	b_t_bep,b_e_bep,b_q_bep,dlg_bep, &
	1226	dl_u_bep,sf_bep,vl_bep, &
	1227	rl_up_urb,rs_abs_urb,emiss_urb,grdflx_urb,qv3d, &
	1228	ids,ide, jds,jde, kds,kde, &
	1229	ims,ime, jms,jme, kms,kme, &
	1230	its,ite, jts,jte, kts,kte )
	1231
	1232	ENDIF
	1233
	1234	if((sf_urban_physics.eq.2).OR.(sf_urban_physics.eq.3))then !Bep begin
	1235	! fix the value of the Stefan-Boltzmann constant
	1236	sigma_sb=5.67e-08
	1237	do j=jts,jte
	1238	do i=its,ite
	1239	UMOM_URB(I,J)=0.
	1240	VMOM_URB(I,J)=0.
	1241	HFX_URB(I,J)=0.
	1242	QFX_URB(I,J)=0.
	1243	do k=kts,kte
	1244	a_u_bep(i,k,j)=a_u_bep(i,k,j)*frc_urb2d(i,j)
	1245	a_v_bep(i,k,j)=a_v_bep(i,k,j)*frc_urb2d(i,j)
	1246	a_t_bep(i,k,j)=a_t_bep(i,k,j)*frc_urb2d(i,j)
	1247	a_q_bep(i,k,j)=0.
	1248	a_e_bep(i,k,j)=0.
	1249	b_u_bep(i,k,j)=b_u_bep(i,k,j)*frc_urb2d(i,j)
	1250	b_v_bep(i,k,j)=b_v_bep(i,k,j)*frc_urb2d(i,j)
	1251	b_t_bep(i,k,j)=b_t_bep(i,k,j)*frc_urb2d(i,j)
	1252	b_q_bep(i,k,j)=b_q_bep(i,k,j)*frc_urb2d(i,j)
	1253	b_e_bep(i,k,j)=b_e_bep(i,k,j)*frc_urb2d(i,j)
	1254	HFX_URB(I,J)=HFX_URB(I,J)+B_T_BEP(I,K,J)RHO(I,K,J)CP* &
	1255	DZ8W(I,K,J)*VL_BEP(I,K,J)
	1256	QFX_URB(I,J)=QFX_URB(I,J)+B_Q_BEP(I,K,J)* &
	1257	DZ8W(I,K,J)*VL_BEP(I,K,J)
	1258	UMOM_URB(I,J)=UMOM_URB(I,J)+ (A_U_BEP(I,K,J)*U_PHY(I,K,J)+ &
	1259	B_U_BEP(I,K,J))DZ8W(I,K,J)VL_BEP(I,K,J)
	1260	VMOM_URB(I,J)=VMOM_URB(I,J)+ (A_V_BEP(I,K,J)*V_PHY(I,K,J)+ &
	1261	B_V_BEP(I,K,J))DZ8W(I,K,J)VL_BEP(I,K,J)
	1262	vl_bep(i,k,j)=(1.-frc_urb2d(i,j))+vl_bep(i,k,j)*frc_urb2d(i,j)
	1263	sf_bep(i,k,j)=(1.-frc_urb2d(i,j))+sf_bep(i,k,j)*frc_urb2d(i,j)
	1264	end do
	1265	a_u_bep(i,1,j)=(1.-frc_urb2d(i,j))(-ust(I,J)ust(I,J))/dz8w(i,1,j)/ &
	1266	((u_phy(i,1,j)2+v_phy(i,1,j)2.)**.5)+a_u_bep(i,1,j)
	1267	a_v_bep(i,1,j)=(1.-frc_urb2d(i,j))(-ust(I,J)ust(I,J))/dz8w(i,1,j)/ &
	1268	((u_phy(i,1,j)2+v_phy(i,1,j)2.)**.5)+a_v_bep(i,1,j)
	1269	b_t_bep(i,1,j)=(1.-frc_urb2d(i,j))*hfx_rural(i,j)/dz8w(i,1,j)/rho(i,1,j)/CP+ &
	1270	b_t_bep(i,1,j)
	1271	b_q_bep(i,1,j)=(1.-frc_urb2d(i,j))*qfx_rural(i,j)/dz8w(i,1,j)/rho(i,1,j)+b_q_bep(i,1,j)
	1272	umom=(1.-frc_urb2d(i,j))ust(i,j)ust(i,j)*u_phy(i,1,j)/ &
	1273	((u_phy(i,1,j)2+v_phy(i,1,j)2.)**.5)+umom_urb(i,j)
	1274	vmom=(1.-frc_urb2d(i,j))ust(i,j)ust(i,j)*v_phy(i,1,j)/ &
	1275	((u_phy(i,1,j)2+v_phy(i,1,j)2.)**.5)+vmom_urb(i,j)
	1276	sf_bep(i,1,j)=1.
	1277
	1278	! compute upward longwave radiation from the rural part and total
	1279	! rl_up_rural=-emiss_rural(i,j)sigma_sb(tsk_rural(i,j)*4.)-(1.-emiss_rural(i,j))glw(i,j)
	1280	! rl_up_tot=(1.-frc_urb2d(i,j))rl_up_rural+frc_urb2d(i,j)rl_up_urb(i,j)
	1281	! emiss(i,j)=(1.-frc_urb2d(i,j))emiss_rural(i,j)+frc_urb2d(i,j)emiss_urb(i,j)
	1282	! using the emissivity and the total longwave upward radiation estimate the averaged skin temperature
	1283	IF (FRC_URB2D(I,J).GT.0.) THEN
	1284	rl_up_rural=-emiss_rural(i,j)sigma_sb(tsk_rural(i,j)*4.)-(1.-emiss_rural(i,j))glw(i,j)
	1285	rl_up_tot=(1.-frc_urb2d(i,j))rl_up_rural+frc_urb2d(i,j)rl_up_urb(i,j)
	1286	emiss(i,j)=(1.-frc_urb2d(i,j))emiss_rural(i,j)+frc_urb2d(i,j)emiss_urb(i,j)
	1287	ts_urb2d(i,j)=(max(0.,(-rl_up_urb(i,j)-(1.-emiss_urb(i,j))glw(i,j))/emiss_urb(i,j)/sigma_sb))*0.25
	1288	tsk(i,j)=(max(0., (-1.rl_up_tot-(1.-emiss(i,j))glw(i,j) )/emiss(i,j)/sigma_sb))**.25
	1289	rs_abs_tot=(1.-frc_urb2d(i,j))swdown(i,j)(1.-albedo(i,j))+frc_urb2d(i,j)*rs_abs_urb(i,j)
	1290	if(swdown(i,j).gt.0.)then
	1291	albedo(i,j)=1.-rs_abs_tot/swdown(i,j)
	1292	else
	1293	albedo(i,j)=alb_rural(i,j)
	1294	endif
	1295	! rename *_urb to sh_urb2d,lh_urb2d,g_urb2d,rn_urb2d
	1296	grdflx(i,j)= (1.-frc_urb2d(i,j))grdflx_rural(i,j)+frc_urb2d(i,j)grdflx_urb(i,j)
	1297	qfx(i,j)=(1.-frc_urb2d(i,j))*qfx_rural(i,j)+qfx_urb(i,j)
	1298	! lh(i,j)=(1.-frc_urb2d(i,j))qfx_rural(i,j)xlv
	1299	lh(i,j)=qfx(i,j)*xlv
	1300	HFX(I,J) = HFX_URB(I,J)+(1-FRC_URB2D(I,J))*HFX_RURAL(I,J) ![W/m/m]
	1301	SH_URB2D(I,J) = HFX_URB(I,J)/FRC_URB2D(I,J)
	1302	LH_URB2D(I,J) = qfx_urb(i,j)*xlv
	1303	G_URB2D(I,J) = grdflx_urb(i,j)
	1304	RN_URB2D(I,J) = rs_abs_urb(i,j)+emiss_urb(i,j)*glw(i,j)-rl_up_urb(i,j)
	1305	ust(i,j)=(umom2.+vmom2.)**.25
	1306	! if(tsk(i,j).gt.350)write(,)'tsk too big!',i,j,tsk(i,j)
	1307	! if(tsk(i,j).lt.260)write(,)'tsk too small!',i,j,tsk(i,j),rl_up_tot,rl_up_urb(i,j),rl_up_rural
	1308	! print*,'ivgtyp,i,j,sigma_sb',ivgtyp(i,j),i,j,sigma_sb
	1309	! print*,'hfx,lh,qfx,grdflx,ts_urb2d',hfx(i,j),lh(i,j),qfx(i,j),grdflx(i,j),ts_urb2d(i,j)
	1310	! print*,'tsk,albedo,emiss',tsk(i,j),albedo(i,j),emiss(i,j)
	1311	! if(i.eq.56.and.j.eq.29)then
	1312	! print*,'ivgtyp, qfx, hfx',ivgtyp(i,j),hfx_rural(i,j),qfx_rural(i,j)
	1313	! print*,'emiss_rural,emiss_urb',emiss_rural(i,j),emiss_urb(i,j)
	1314	! print*,'rl_up_rural,rl_up_urb(i,j)',rl_up_rural,rl_up_urb(i,j)
	1315	! print*,'tsk_rural,ts_urb2d(i,j),tsk',tsk_rural(i,j),ts_urb2d(i,j),tsk(i,j)
	1316	! print,'reconstruction fei',((emiss(i,j)tsk(i,j)*4.-frc_urb2d(i,j)emiss_urb(i,j)ts_urb2d(i,j)4.)/(emiss_rural(i,j)(1.-frc_urb2d(i,j))))**.25
	1317	! print*,'ivgtyp,hfx,hfx_urb,hfx_rural',hfx(i,j),hfx_urb(i,j),hfx_rural(i,j)
	1318	! print*,'lh,lh_rural',lh(i,j),lh_rural(i,j)
	1319	! print*,'qfx',qfx(i,j)
	1320	! print*,'ts_urb2d',ts_urb2d(i,j)
	1321	! print*,'ust',ust(i,j)
	1322	! print*,'swdown,glw',swdown(i,j),glw(i,j)
	1323	! endif
	1324	else
	1325	SH_URB2D(I,J) = 0.
	1326	LH_URB2D(I,J) = 0.
	1327	G_URB2D(I,J) = 0.
	1328	RN_URB2D(I,J) = 0.
	1329	endif
	1330	! IF( IVGTYP(I,J) == 1 .or. IVGTYP(I,J) == 31 .or. &
	1331	! IVGTYP(I,J) == 32 .or. IVGTYP(I,J) == 33) THEN
	1332	! print*,'ivgtyp, qfx, hfx',ivgtyp(i,j),hfx_rural(i,j),qfx_rural(i,j)
	1333	! print*,'ivgtyp,hfx,hfx_urb,hfx_rural',hfx(i,j),hfx_urb(i,j),hfx_rural(i,j)
	1334	! print*,'lh,lh_rural',lh(i,j),lh_rural(i,j)
	1335	! print*,'qfx',qfx(i,j)
	1336	! print*,'ts_urb2d',ts_urb2d(i,j)
	1337	! print*,'ust',ust(i,j)
	1338	! endif
	1339	enddo
	1340	enddo
	1341
	1342
	1343	endif !Bep end
	1344
	1345	!------------------------------------------------------
	1346	END SUBROUTINE lsm
	1347	!------------------------------------------------------
	1348
	1349	SUBROUTINE LSMINIT(VEGFRA,SNOW,SNOWC,SNOWH,CANWAT,SMSTAV, &
	1350	SMSTOT, SFCRUNOFF,UDRUNOFF,ACSNOW, &
	1351	ACSNOM,IVGTYP,ISLTYP,TSLB,SMOIS,SH2O,ZS,DZS, &
	1352	MMINLU, &
	1353	SNOALB, FNDSOILW, FNDSNOWH, RDMAXALB, &
	1354	num_soil_layers, restart, &
	1355	allowed_to_read , &
	1356	ids,ide, jds,jde, kds,kde, &
	1357	ims,ime, jms,jme, kms,kme, &
	1358	its,ite, jts,jte, kts,kte )
	1359
	1360	INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
	1361	ims,ime, jms,jme, kms,kme, &
	1362	its,ite, jts,jte, kts,kte
	1363
	1364	INTEGER, INTENT(IN) :: num_soil_layers
	1365
	1366	LOGICAL , INTENT(IN) :: restart , allowed_to_read
	1367
	1368	REAL, DIMENSION( num_soil_layers), INTENT(INOUT) :: ZS, DZS
	1369
	1370	REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , &
	1371	INTENT(INOUT) :: SMOIS, & !Total soil moisture
	1372	SH2O, & !liquid soil moisture
	1373	TSLB !STEMP
	1374
	1375	REAL, DIMENSION( ims:ime, jms:jme ) , &
	1376	INTENT(INOUT) :: SNOW, &
	1377	SNOWH, &
	1378	SNOWC, &
	1379	SNOALB, &
	1380	CANWAT, &
	1381	SMSTAV, &
	1382	SMSTOT, &
	1383	SFCRUNOFF, &
	1384	UDRUNOFF, &
	1385	ACSNOW, &
	1386	VEGFRA, &
	1387	ACSNOM
	1388
	1389	INTEGER, DIMENSION( ims:ime, jms:jme ) , &
	1390	INTENT(IN) :: IVGTYP, &
	1391	ISLTYP
	1392	CHARACTER(LEN=*), INTENT(IN) :: MMINLU
	1393
	1394	LOGICAL, INTENT(IN) :: FNDSOILW , &
	1395	FNDSNOWH
	1396	LOGICAL, INTENT(IN) :: RDMAXALB
	1397
	1398
	1399	INTEGER :: L
	1400	REAL :: BX, SMCMAX, PSISAT, FREE
	1401	REAL, PARAMETER :: BLIM = 5.5, HLICE = 3.335E5, &
	1402	GRAV = 9.81, T0 = 273.15
	1403	INTEGER :: errflag
	1404
	1405	character*256 :: MMINSL
	1406	MMINSL='STAS'
	1407	!
	1408
	1409	! initialize three Noah LSM related tables
	1410	IF ( allowed_to_read ) THEN
	1411	CALL wrf_message( 'INITIALIZE THREE Noah LSM RELATED TABLES' )
	1412	CALL SOIL_VEG_GEN_PARM( MMINLU, MMINSL )
	1413	ENDIF
	1414
	1415	#ifdef WRF_CHEM
	1416	!
	1417	! need this parameter for dust parameterization in wrf/chem
	1418	!
	1419	do I=1,NSLTYPE
	1420	porosity(i)=maxsmc(i)
	1421	enddo
	1422	#endif
	1423
	1424	IF(.not.restart)THEN
	1425
	1426	itf=min0(ite,ide-1)
	1427	jtf=min0(jte,jde-1)
	1428
	1429	errflag = 0
	1430	DO j = jts,jtf
	1431	DO i = its,itf
	1432	IF ( ISLTYP( i,j ) .LT. 1 ) THEN
	1433	errflag = 1
	1434	WRITE(err_message,*)"module_sf_noahlsm.F: lsminit: out of range ISLTYP ",i,j,ISLTYP( i,j )
	1435	CALL wrf_message(err_message)
	1436	ENDIF
	1437	IF(.not.RDMAXALB) THEN
	1438	SNOALB(i,j)=MAXALB(IVGTYP(i,j))*0.01
	1439	ENDIF
	1440	ENDDO
	1441	ENDDO
	1442	IF ( errflag .EQ. 1 ) THEN
	1443	CALL wrf_error_fatal( "module_sf_noahlsm.F: lsminit: out of range value "// &
	1444	"of ISLTYP. Is this field in the input?" )
	1445	ENDIF
	1446
	1447	! initialize soil liquid water content SH2O
	1448
	1449	! IF(.NOT.FNDSOILW) THEN
	1450
	1451	! If no SWC, do the following
	1452	! PRINT *,'SOIL WATER NOT FOUND - VALUE SET IN LSMINIT'
	1453	DO J = jts,jtf
	1454	DO I = its,itf
	1455	BX = BB(ISLTYP(I,J))
	1456	SMCMAX = MAXSMC(ISLTYP(I,J))
	1457	PSISAT = SATPSI(ISLTYP(I,J))
	1458	if ((bx > 0.0).and.(smcmax > 0.0).and.(psisat > 0.0)) then
	1459	DO NS=1, num_soil_layers
	1460	! ----------------------------------------------------------------------
	1461	!SH2O <= SMOIS for T < 273.149K (-0.001C)
	1462	IF (TSLB(I,NS,J) < 273.149) THEN
	1463	! ----------------------------------------------------------------------
	1464	! first guess following explicit solution for Flerchinger Eqn from Koren
	1465	! et al, JGR, 1999, Eqn 17 (KCOUNT=0 in FUNCTION FRH2O).
	1466	! ISLTPK is soil type
	1467	BX = BB(ISLTYP(I,J))
	1468	SMCMAX = MAXSMC(ISLTYP(I,J))
	1469	PSISAT = SATPSI(ISLTYP(I,J))
	1470	IF ( BX > BLIM ) BX = BLIM
	1471	FK=(( (HLICE/(GRAV(-PSISAT))) &
	1472	((TSLB(I,NS,J)-T0)/TSLB(I,NS,J)) )*(-1/BX) )SMCMAX
	1473	IF (FK < 0.02) FK = 0.02
	1474	SH2O(I,NS,J) = MIN( FK, SMOIS(I,NS,J) )
	1475	! ----------------------------------------------------------------------
	1476	! now use iterative solution for liquid soil water content using
	1477	! FUNCTION FRH2O with the initial guess for SH2O from above explicit
	1478	! first guess.
	1479	CALL FRH2O (FREE,TSLB(I,NS,J),SMOIS(I,NS,J),SH2O(I,NS,J), &
	1480	SMCMAX,BX,PSISAT)
	1481	SH2O(I,NS,J) = FREE
	1482	ELSE ! of IF (TSLB(I,NS,J)
	1483	! ----------------------------------------------------------------------
	1484	! SH2O = SMOIS ( for T => 273.149K (-0.001C)
	1485	SH2O(I,NS,J)=SMOIS(I,NS,J)
	1486	! ----------------------------------------------------------------------
	1487	ENDIF ! of IF (TSLB(I,NS,J)
	1488	END DO ! of DO NS=1, num_soil_layers
	1489	else ! of if ((bx > 0.0)
	1490	DO NS=1, num_soil_layers
	1491	SH2O(I,NS,J)=SMOIS(I,NS,J)
	1492	END DO
	1493	endif ! of if ((bx > 0.0)
	1494	ENDDO ! DO I = its,itf
	1495	ENDDO ! DO J = jts,jtf
	1496	! ENDIF ! of IF(.NOT.FNDSOILW)THEN
	1497
	1498	! initialize physical snow height SNOWH
	1499
	1500	IF(.NOT.FNDSNOWH)THEN
	1501	! If no SNOWH do the following
	1502	CALL wrf_message( 'SNOW HEIGHT NOT FOUND - VALUE DEFINED IN LSMINIT' )
	1503	DO J = jts,jtf
	1504	DO I = its,itf
	1505	SNOWH(I,J)=SNOW(I,J)*0.005 ! SNOW in mm and SNOWH in m
	1506	ENDDO
	1507	ENDDO
	1508	ENDIF
	1509
	1510	! initialize canopy water to ZERO
	1511
	1512	! GO TO 110
	1513	! print*,'Note that canopy water content (CANWAT) is set to ZERO in LSMINIT'
	1514	DO J = jts,jtf
	1515	DO I = its,itf
	1516	CANWAT(I,J)=0.0
	1517	ENDDO
	1518	ENDDO
	1519	110 CONTINUE
	1520
	1521	ENDIF
	1522	!------------------------------------------------------------------------------
	1523	END SUBROUTINE lsminit
	1524	!------------------------------------------------------------------------------
	1525
	1526
	1527
	1528	!-----------------------------------------------------------------
	1529	SUBROUTINE SOIL_VEG_GEN_PARM( MMINLU, MMINSL)
	1530	!-----------------------------------------------------------------
	1531
	1532	USE module_wrf_error
	1533	IMPLICIT NONE
	1534
	1535	CHARACTER(LEN=*), INTENT(IN) :: MMINLU, MMINSL
	1536	integer :: LUMATCH, IINDEX, LC, NUM_SLOPE
	1537	integer :: ierr
	1538	INTEGER , PARAMETER :: OPEN_OK = 0
	1539
	1540	character*128 :: mess , message
	1541	logical, external :: wrf_dm_on_monitor
	1542
	1543
	1544	!-----SPECIFY VEGETATION RELATED CHARACTERISTICS :
	1545	! ALBBCK: SFC albedo (in percentage)
	1546	! Z0: Roughness length (m)
	1547	! SHDFAC: Green vegetation fraction (in percentage)
	1548	! Note: The ALBEDO, Z0, and SHDFAC values read from the following table
	1549	! ALBEDO, amd Z0 are specified in LAND-USE TABLE; and SHDFAC is
	1550	! the monthly green vegetation data
	1551	! CMXTBL: MAX CNPY Capacity (m)
	1552	! NROTBL: Rooting depth (layer)
	1553	! RSMIN: Mimimum stomatal resistance (s m-1)
	1554	! RSMAX: Max. stomatal resistance (s m-1)
	1555	! RGL: Parameters used in radiation stress function
	1556	! HS: Parameter used in vapor pressure deficit functio
	1557	! TOPT: Optimum transpiration air temperature. (K)
	1558	! CMCMAX: Maximum canopy water capacity
	1559	! CFACTR: Parameter used in the canopy inteception calculati
	1560	! SNUP: Threshold snow depth (in water equivalent m) that
	1561	! implies 100% snow cover
	1562	! LAI: Leaf area index (dimensionless)
	1563	! MAXALB: Upper bound on maximum albedo over deep snow
	1564	!
	1565	!-----READ IN VEGETAION PROPERTIES FROM VEGPARM.TBL
	1566	!
	1567
	1568	IF ( wrf_dm_on_monitor() ) THEN
	1569
	1570	OPEN(19, FILE='VEGPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
	1571	IF(ierr .NE. OPEN_OK ) THEN
	1572	WRITE(message,FMT='(A)') &
	1573	'module_sf_noahlsm.F: soil_veg_gen_parm: failure opening VEGPARM.TBL'
	1574	CALL wrf_error_fatal ( message )
	1575	END IF
	1576
	1577
	1578	LUMATCH=0
	1579
	1580	FIND_LUTYPE : DO WHILE (LUMATCH == 0)
	1581	READ (19,*,END=2002)
	1582	READ (19,*,END=2002)LUTYPE
	1583	READ (19,*)LUCATS,IINDEX
	1584
	1585	IF(LUTYPE.EQ.MMINLU)THEN
	1586	WRITE( mess , * ) 'LANDUSE TYPE = ' // TRIM ( LUTYPE ) // ' FOUND', LUCATS,' CATEGORIES'
	1587	CALL wrf_message( mess )
	1588	LUMATCH=1
	1589	ELSE
	1590	call wrf_message ( "Skipping over LUTYPE = " // TRIM ( LUTYPE ) )
	1591	DO LC = 1, LUCATS+12
	1592	read(19,*)
	1593	ENDDO
	1594	ENDIF
	1595	ENDDO FIND_LUTYPE
	1596	! prevent possible array overwrite, Bill Bovermann, IBM, May 6, 2008
	1597	IF ( SIZE(SHDTBL) < LUCATS .OR. &
	1598	SIZE(NROTBL) < LUCATS .OR. &
	1599	SIZE(RSTBL) < LUCATS .OR. &
	1600	SIZE(RGLTBL) < LUCATS .OR. &
	1601	SIZE(HSTBL) < LUCATS .OR. &
	1602	SIZE(SNUPTBL) < LUCATS .OR. &
	1603	SIZE(MAXALB) < LUCATS .OR. &
	1604	SIZE(LAIMINTBL) < LUCATS .OR. &
	1605	SIZE(LAIMAXTBL) < LUCATS .OR. &
	1606	SIZE(Z0MINTBL) < LUCATS .OR. &
	1607	SIZE(Z0MAXTBL) < LUCATS .OR. &
	1608	SIZE(ALBEDOMINTBL) < LUCATS .OR. &
	1609	SIZE(ALBEDOMAXTBL) < LUCATS .OR. &
	1610	SIZE(EMISSMINTBL ) < LUCATS .OR. &
	1611	SIZE(EMISSMAXTBL ) < LUCATS ) THEN
	1612	CALL wrf_error_fatal('Table sizes too small for value of LUCATS in module_sf_noahdrv.F')
	1613	ENDIF
	1614
	1615	IF(LUTYPE.EQ.MMINLU)THEN
	1616	DO LC=1,LUCATS
	1617	READ (19,*)IINDEX,SHDTBL(LC), &
	1618	NROTBL(LC),RSTBL(LC),RGLTBL(LC),HSTBL(LC), &
	1619	SNUPTBL(LC),MAXALB(LC), LAIMINTBL(LC), &
	1620	LAIMAXTBL(LC),EMISSMINTBL(LC), &
	1621	EMISSMAXTBL(LC), ALBEDOMINTBL(LC), &
	1622	ALBEDOMAXTBL(LC), Z0MINTBL(LC), Z0MAXTBL(LC)
	1623	ENDDO
	1624	!
	1625	READ (19,*)
	1626	READ (19,*)TOPT_DATA
	1627	READ (19,*)
	1628	READ (19,*)CMCMAX_DATA
	1629	READ (19,*)
	1630	READ (19,*)CFACTR_DATA
	1631	READ (19,*)
	1632	READ (19,*)RSMAX_DATA
	1633	READ (19,*)
	1634	READ (19,*)BARE
	1635	READ (19,*)
	1636	READ (19,*)NATURAL
	1637	ENDIF
	1638	!
	1639	2002 CONTINUE
	1640
	1641	CLOSE (19)
	1642	IF (LUMATCH == 0) then
	1643	CALL wrf_error_fatal ("Land Use Dataset '"//MMINLU//"' not found in VEGPARM.TBL.")
	1644	ENDIF
	1645	ENDIF
	1646
	1647	CALL wrf_dm_bcast_string ( LUTYPE , 4 )
	1648	CALL wrf_dm_bcast_integer ( LUCATS , 1 )
	1649	CALL wrf_dm_bcast_integer ( IINDEX , 1 )
	1650	CALL wrf_dm_bcast_integer ( LUMATCH , 1 )
	1651	CALL wrf_dm_bcast_real ( SHDTBL , NLUS )
	1652	CALL wrf_dm_bcast_real ( NROTBL , NLUS )
	1653	CALL wrf_dm_bcast_real ( RSTBL , NLUS )
	1654	CALL wrf_dm_bcast_real ( RGLTBL , NLUS )
	1655	CALL wrf_dm_bcast_real ( HSTBL , NLUS )
	1656	CALL wrf_dm_bcast_real ( SNUPTBL , NLUS )
	1657	CALL wrf_dm_bcast_real ( LAIMINTBL , NLUS )
	1658	CALL wrf_dm_bcast_real ( LAIMAXTBL , NLUS )
	1659	CALL wrf_dm_bcast_real ( Z0MINTBL , NLUS )
	1660	CALL wrf_dm_bcast_real ( Z0MAXTBL , NLUS )
	1661	CALL wrf_dm_bcast_real ( EMISSMINTBL , NLUS )
	1662	CALL wrf_dm_bcast_real ( EMISSMAXTBL , NLUS )
	1663	CALL wrf_dm_bcast_real ( ALBEDOMINTBL , NLUS )
	1664	CALL wrf_dm_bcast_real ( ALBEDOMAXTBL , NLUS )
	1665	CALL wrf_dm_bcast_real ( MAXALB , NLUS )
	1666	CALL wrf_dm_bcast_real ( TOPT_DATA , 1 )
	1667	CALL wrf_dm_bcast_real ( CMCMAX_DATA , 1 )
	1668	CALL wrf_dm_bcast_real ( CFACTR_DATA , 1 )
	1669	CALL wrf_dm_bcast_real ( RSMAX_DATA , 1 )
	1670	CALL wrf_dm_bcast_integer ( BARE , 1 )
	1671	CALL wrf_dm_bcast_integer ( NATURAL , 1 )
	1672
	1673	!
	1674	!-----READ IN SOIL PROPERTIES FROM SOILPARM.TBL
	1675	!
	1676	IF ( wrf_dm_on_monitor() ) THEN
	1677	OPEN(19, FILE='SOILPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
	1678	IF(ierr .NE. OPEN_OK ) THEN
	1679	WRITE(message,FMT='(A)') &
	1680	'module_sf_noahlsm.F: soil_veg_gen_parm: failure opening SOILPARM.TBL'
	1681	CALL wrf_error_fatal ( message )
	1682	END IF
	1683
	1684	WRITE(mess,*) 'INPUT SOIL TEXTURE CLASSIFICATION = ', TRIM ( MMINSL )
	1685	CALL wrf_message( mess )
	1686
	1687	LUMATCH=0
	1688
	1689	READ (19,*)
	1690	READ (19,2000,END=2003)SLTYPE
	1691	2000 FORMAT (A4)
	1692	READ (19,*)SLCATS,IINDEX
	1693	IF(SLTYPE.EQ.MMINSL)THEN
	1694	WRITE( mess , * ) 'SOIL TEXTURE CLASSIFICATION = ', TRIM ( SLTYPE ) , ' FOUND', &
	1695	SLCATS,' CATEGORIES'
	1696	CALL wrf_message ( mess )
	1697	LUMATCH=1
	1698	ENDIF
	1699	! prevent possible array overwrite, Bill Bovermann, IBM, May 6, 2008
	1700	IF ( SIZE(BB ) < SLCATS .OR. &
	1701	SIZE(DRYSMC) < SLCATS .OR. &
	1702	SIZE(F11 ) < SLCATS .OR. &
	1703	SIZE(MAXSMC) < SLCATS .OR. &
	1704	SIZE(REFSMC) < SLCATS .OR. &
	1705	SIZE(SATPSI) < SLCATS .OR. &
	1706	SIZE(SATDK ) < SLCATS .OR. &
	1707	SIZE(SATDW ) < SLCATS .OR. &
	1708	SIZE(WLTSMC) < SLCATS .OR. &
	1709	SIZE(QTZ ) < SLCATS ) THEN
	1710	CALL wrf_error_fatal('Table sizes too small for value of SLCATS in module_sf_noahdrv.F')
	1711	ENDIF
	1712	IF(SLTYPE.EQ.MMINSL)THEN
	1713	DO LC=1,SLCATS
	1714	READ (19,*) IINDEX,BB(LC),DRYSMC(LC),F11(LC),MAXSMC(LC),&
	1715	REFSMC(LC),SATPSI(LC),SATDK(LC), SATDW(LC), &
	1716	WLTSMC(LC), QTZ(LC)
	1717	ENDDO
	1718	ENDIF
	1719
	1720	2003 CONTINUE
	1721
	1722	CLOSE (19)
	1723	ENDIF
	1724
	1725	CALL wrf_dm_bcast_integer ( LUMATCH , 1 )
	1726	CALL wrf_dm_bcast_string ( SLTYPE , 4 )
	1727	CALL wrf_dm_bcast_string ( MMINSL , 4 ) ! since this is reset above, see oct2 ^
	1728	CALL wrf_dm_bcast_integer ( SLCATS , 1 )
	1729	CALL wrf_dm_bcast_integer ( IINDEX , 1 )
	1730	CALL wrf_dm_bcast_real ( BB , NSLTYPE )
	1731	CALL wrf_dm_bcast_real ( DRYSMC , NSLTYPE )
	1732	CALL wrf_dm_bcast_real ( F11 , NSLTYPE )
	1733	CALL wrf_dm_bcast_real ( MAXSMC , NSLTYPE )
	1734	CALL wrf_dm_bcast_real ( REFSMC , NSLTYPE )
	1735	CALL wrf_dm_bcast_real ( SATPSI , NSLTYPE )
	1736	CALL wrf_dm_bcast_real ( SATDK , NSLTYPE )
	1737	CALL wrf_dm_bcast_real ( SATDW , NSLTYPE )
	1738	CALL wrf_dm_bcast_real ( WLTSMC , NSLTYPE )
	1739	CALL wrf_dm_bcast_real ( QTZ , NSLTYPE )
	1740
	1741	IF(LUMATCH.EQ.0)THEN
	1742	CALL wrf_message( 'SOIl TEXTURE IN INPUT FILE DOES NOT ' )
	1743	CALL wrf_message( 'MATCH SOILPARM TABLE' )
	1744	CALL wrf_error_fatal ( 'INCONSISTENT OR MISSING SOILPARM FILE' )
	1745	ENDIF
	1746
	1747	!
	1748	!-----READ IN GENERAL PARAMETERS FROM GENPARM.TBL
	1749	!
	1750	IF ( wrf_dm_on_monitor() ) THEN
	1751	OPEN(19, FILE='GENPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
	1752	IF(ierr .NE. OPEN_OK ) THEN
	1753	WRITE(message,FMT='(A)') &
	1754	'module_sf_noahlsm.F: soil_veg_gen_parm: failure opening GENPARM.TBL'
	1755	CALL wrf_error_fatal ( message )
	1756	END IF
	1757
	1758	READ (19,*)
	1759	READ (19,*)
	1760	READ (19,*) NUM_SLOPE
	1761
	1762	SLPCATS=NUM_SLOPE
	1763	! prevent possible array overwrite, Bill Bovermann, IBM, May 6, 2008
	1764	IF ( SIZE(slope_data) < NUM_SLOPE ) THEN
	1765	CALL wrf_error_fatal('NUM_SLOPE too large for slope_data array in module_sf_noahdrv')
	1766	ENDIF
	1767
	1768	DO LC=1,SLPCATS
	1769	READ (19,*)SLOPE_DATA(LC)
	1770	ENDDO
	1771
	1772	READ (19,*)
	1773	READ (19,*)SBETA_DATA
	1774	READ (19,*)
	1775	READ (19,*)FXEXP_DATA
	1776	READ (19,*)
	1777	READ (19,*)CSOIL_DATA
	1778	READ (19,*)
	1779	READ (19,*)SALP_DATA
	1780	READ (19,*)
	1781	READ (19,*)REFDK_DATA
	1782	READ (19,*)
	1783	READ (19,*)REFKDT_DATA
	1784	READ (19,*)
	1785	READ (19,*)FRZK_DATA
	1786	READ (19,*)
	1787	READ (19,*)ZBOT_DATA
	1788	READ (19,*)
	1789	READ (19,*)CZIL_DATA
	1790	READ (19,*)
	1791	READ (19,*)SMLOW_DATA
	1792	READ (19,*)
	1793	READ (19,*)SMHIGH_DATA
	1794	READ (19,*)
	1795	READ (19,*)LVCOEF_DATA
	1796	CLOSE (19)
	1797	ENDIF
	1798
	1799	CALL wrf_dm_bcast_integer ( NUM_SLOPE , 1 )
	1800	CALL wrf_dm_bcast_integer ( SLPCATS , 1 )
	1801	CALL wrf_dm_bcast_real ( SLOPE_DATA , NSLOPE )
	1802	CALL wrf_dm_bcast_real ( SBETA_DATA , 1 )
	1803	CALL wrf_dm_bcast_real ( FXEXP_DATA , 1 )
	1804	CALL wrf_dm_bcast_real ( CSOIL_DATA , 1 )
	1805	CALL wrf_dm_bcast_real ( SALP_DATA , 1 )
	1806	CALL wrf_dm_bcast_real ( REFDK_DATA , 1 )
	1807	CALL wrf_dm_bcast_real ( REFKDT_DATA , 1 )
	1808	CALL wrf_dm_bcast_real ( FRZK_DATA , 1 )
	1809	CALL wrf_dm_bcast_real ( ZBOT_DATA , 1 )
	1810	CALL wrf_dm_bcast_real ( CZIL_DATA , 1 )
	1811	CALL wrf_dm_bcast_real ( SMLOW_DATA , 1 )
	1812	CALL wrf_dm_bcast_real ( SMHIGH_DATA , 1 )
	1813	CALL wrf_dm_bcast_real ( LVCOEF_DATA , 1 )
	1814
	1815
	1816	!-----------------------------------------------------------------
	1817	END SUBROUTINE SOIL_VEG_GEN_PARM
	1818	!-----------------------------------------------------------------
	1819
	1820	END MODULE module_sf_noahdrv

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

Download in other formats:

Original Format