Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

module_sf_noahdrv.F @ 1

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

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

WRF: version v3.3
LMDZ: version v1818

More details in:

File size: 81.2 KB

Line
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