| blob | 5035c586d28ceb4efa76765ba82747747689c8fc |
1 MODULE module_sf_noahdrv
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 !-------------------------------
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
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)
182 !-- EMISS surface emissivity (between 0 and 1)
183 !-- EMBCK Background surface emissivity (between 0 and 1)
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 !----------------------------------------------------------------
209 ! IN only
211 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
212 ims,ime, jms,jme, kms,kme, &
213 its,ite, jts,jte, kts,kte
215 INTEGER, INTENT(IN ) :: sf_urban_physics !urban
216 INTEGER, INTENT(IN ) :: isurban
217 INTEGER, INTENT(IN ) :: isice
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
234 REAL, DIMENSION( ims:ime, jms:jme ) , &
235 INTENT(INOUT) :: ALBBCK, &
236 Z0
237 CHARACTER(LEN=*), INTENT(IN ) :: MMINLU
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
248 INTEGER, DIMENSION( ims:ime, jms:jme ) , &
249 INTENT(IN ) :: IVGTYP, &
250 ISLTYP
252 INTEGER, INTENT(IN) :: num_soil_layers,ITIMESTEP
254 REAL, INTENT(IN ) :: DT,ROVCP
256 REAL, DIMENSION(1:num_soil_layers), INTENT(IN)::DZS
258 ! IN and OUT
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
265 REAL, DIMENSION( ims:ime , 1:num_soil_layers, jms:jme ), &
266 INTENT(OUT) :: SMCREL
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
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
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)
309 REAL, DIMENSION(1:num_soil_layers) :: ET
311 REAL, DIMENSION(1:num_soil_layers) :: SMAV
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
318 LOGICAL, INTENT(IN ) :: myj,frpcpn
320 ! DECLARATIONS - LOGICAL
321 ! ----------------------------------------------------------------------
322 LOGICAL, PARAMETER :: LOCAL=.false.
323 LOGICAL :: FRZGRA, SNOWNG
325 LOGICAL :: IPRINT
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 ! ----------------------------------------------------------------------
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
351 REAL :: EMISSI
353 REAL :: SNCOVR,SNEQV,SNOWHK,CMC, CHK,TH2
355 REAL :: SMCDRY,SMCMAX,SMCREF,SMCWLT,SNOMLT,SOILM,SOILW,Q1,T1
356 REAL :: SNOTIME1 ! LSTSNW1 INITIAL NUMBER OF TIMESTEPS SINCE LAST SNOWFALL
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
370 ! ----------------------------------------------------------------------
371 ! DECLARATIONS START - urban
372 ! ----------------------------------------------------------------------
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
392 LOGICAL, intent(in) :: rdlai2d
393 LOGICAL, intent(in) :: USEMONALB
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
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
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
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
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 ! ----------------------------------------------------------------------
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
553 ! MEK MAY 2007
554 FDTLIW=DT/ROWLIW
555 ! MEK JUL2007
556 FDTW=DT/(XLV*RHOWATER)
557 ! debug printout
558 IPRINT=.false.
560 ! SLOPETYP=2
561 SLOPETYP=1
562 ! SHDMIN=0.00
565 NSOIL=num_soil_layers
567 DO NS=1,NSOIL
568 SLDPTH(NS)=DZS(NS)
569 ENDDO
571 JLOOP : DO J=jts,jte
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
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
614 !-----------------------------------------------------------------------
615 ILOOP : DO 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
635 SFCTMP=T3D(i,1,j)
636 ZLVL=0.5*DZ8W(i,1,j)
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)
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
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=Q2SAT*A23M4/(SFCTMP-A4)**2
703 IF(SNOW(I,J).GT.0.0)THEN
704 ! snow on surface (use ice saturation properties)
705 SFCTSNO=SFCTMP
706 E2SAT=611.2*EXP(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)) + Q2SATI*SNOWC(I,J)
712 DQSDT2=DQSDT2*(1.-SNOWC(I,J)) + Q2SATI*6174./(SFCTSNO**2)*SNOWC(I,J)
713 ELSE
714 ! cold ground temps, use ice saturation only
715 Q2SAT=Q2SATI
716 DQSDT2=Q2SATI*6174./(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
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)
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
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
760 DO 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 ENDDO
765 !
766 if ( (SNEQV.ne.0..AND.SNOWHK.eq.0.).or.(SNOWHK.le.SNEQV) )THEN
767 SNOWHK= 5.*SNEQV
768 endif
769 !
771 !Fei: urban. for urban surface, if calling UCM, redefine the natural surface in cities as
772 ! the "NATURAL" category in the VEGPARM.TBL
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
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
828 IF (rdlai2d) THEN
829 xlai = lai(i,j)
830 endif
832 IF ( XICE(I,J) >= XICE_THRESHOLD ) THEN
834 ! Sea-ice case
836 DO NS = 1, NSOIL
837 SH2O(I,NS,J) = 1.0
838 ENDDO
839 LAI(I,J) = 0.01
841 CYCLE ILOOP
843 ELSE
845 ! Land and glacial land.
847 CALL SFLX (I,J,ICE,FFROZP, ISURBAN, DT,ZLVL,NSOIL,SLDPTH, & !C
848 LOCAL, & !L
849 LUTYPE, SLTYPE, & !CL
850 LWDN,SOLDN,SOLNET,SFCPRS,PRCP,SFCTMP,Q2K,DUMMY, & !F
851 DUMMY,DUMMY, DUMMY, & !F PRCPRAIN not used
852 TH2,Q2SAT,DQSDT2, & !I
853 VEGTYP,SOILTYP,SLOPETYP,SHDFAC,SHMIN,SHMAX, & !I
854 ALBBRD, SNOALB1,TBOT, Z0BRD, Z0K, EMISSI, EMBRD, & !S
855 CMC,T1,STC,SMC,SWC,SNOWHK,SNEQV,ALBEDOK,CHK,dummy,& !H
856 ETA,SHEAT, ETA_KINEMATIC,FDOWN, & !O
857 EC,EDIR,ET,ETT,ESNOW,DRIP,DEW, & !O
858 BETA,ETP,SSOIL, & !O
859 FLX1,FLX2,FLX3, & !O
860 SNOMLT,SNCOVR, & !O
861 RUNOFF1,RUNOFF2,RUNOFF3, & !O
862 RC,PC,RSMIN,XLAI,RCS,RCT,RCQ,RCSOIL, & !O
863 SOILW,SOILM,Q1,SMAV, & !D
864 RDLAI2D,USEMONALB, &
865 SNOTIME1, &
866 RIBB, &
867 SMCWLT,SMCDRY,SMCREF,SMCMAX,NROOT)
868 ENDIF
870 lai(i,j) = xlai
872 #if 0
873 IF(IPRINT) THEN
875 print*, 'AFTER SFLX, in Noahlsm_driver'
876 print*, 'ICE', ICE, 'DT',DT, 'ZLVL',ZLVL, 'NSOIL', NSOIL, &
877 'SLDPTH', SLDPTH, 'LOCAL',LOCAL, 'LUTYPE',&
878 LUTYPE, 'SLTYPE',SLTYPE, 'LWDN',LWDN, 'SOLDN',SOLDN, &
879 'SFCPRS',SFCPRS, 'PRCP',PRCP,'SFCTMP',SFCTMP,'Q2K',Q2K, &
880 'TH2',TH2,'Q2SAT',Q2SAT,'DQSDT2',DQSDT2,'VEGTYP', VEGTYP,&
881 'SOILTYP',SOILTYP, 'SLOPETYP',SLOPETYP, 'SHDFAC',SHDFAC,&
882 'SHDMIN',SHMIN, 'ALBBRD',ALBBRD,'SNOALB',SNOALB1,'TBOT',&
883 TBOT, 'Z0BRD',Z0BRD, 'Z0K',Z0K, 'CMC',CMC, 'T1',T1,'STC',&
884 STC, 'SMC',SMC, 'SWc',SWC,'SNOWHK',SNOWHK,'SNEQV',SNEQV,&
885 'ALBEDOK',ALBEDOK,'CHK',CHK,'ETA',ETA,'SHEAT',SHEAT, &
886 'ETA_KINEMATIC',ETA_KINEMATIC, 'FDOWN',FDOWN,'EC',EC, &
887 'EDIR',EDIR,'ET',ET,'ETT',ETT,'ESNOW',ESNOW,'DRIP',DRIP,&
888 'DEW',DEW,'BETA',BETA,'ETP',ETP,'SSOIL',SSOIL,'FLX1',FLX1,&
889 'FLX2',FLX2,'FLX3',FLX3,'SNOMLT',SNOMLT,'SNCOVR',SNCOVR,&
890 'RUNOFF1',RUNOFF1,'RUNOFF2',RUNOFF2,'RUNOFF3',RUNOFF3, &
891 'RC',RC, 'PC',PC,'RSMIN',RSMIN,'XLAI',XLAI,'RCS',RCS, &
892 'RCT',RCT,'RCQ',RCQ,'RCSOIL',RCSOIL,'SOILW',SOILW, &
893 'SOILM',SOILM,'Q1',Q1,'SMCWLT',SMCWLT,'SMCDRY',SMCDRY,&
894 'SMCREF',SMCREF,'SMCMAX',SMCMAX,'NROOT',NROOT
895 endif
896 #endif
898 !*** UPDATE STATE VARIABLES
899 CANWAT(I,J)=CMC
900 SNOW(I,J)=SNEQV*1000.
901 ! SNOWH(I,J)=SNOWHK*1000.
902 SNOWH(I,J)=SNOWHK ! SNOWHK in meters
903 ALBEDO(I,J)=ALBEDOK
904 ALB_RURAL(I,J)=ALBEDOK
905 ALBBCK(I,J)=ALBBRD
906 Z0(I,J)=Z0BRD
907 EMISS(I,J) = EMISSI
908 EMISS_RURAL(I,J) = EMISSI
909 ! Noah: activate time-varying roughness length (V3.3 Feb 2011)
910 ZNT(I,J)=Z0K
911 TSK(I,J)=T1
912 TSK_RURAL(I,J)=T1
913 HFX(I,J)=SHEAT
914 HFX_RURAL(I,J)=SHEAT
915 ! MEk Jul07 add potential evap accum
916 POTEVP(I,J)=POTEVP(I,J)+ETP*FDTW
917 QFX(I,J)=ETA_KINEMATIC
918 QFX_RURAL(I,J)=ETA_KINEMATIC
919 LH(I,J)=ETA
920 LH_RURAL(I,J)=ETA
921 GRDFLX(I,J)=SSOIL
922 GRDFLX_RURAL(I,J)=SSOIL
923 SNOWC(I,J)=SNCOVR
924 CHS2(I,J)=CQS2(I,J)
925 SNOTIME(I,J) = SNOTIME1
926 ! prevent diagnostic ground q (q1) from being greater than qsat(tsk)
927 ! as happens over snow cover where the cqs2 value also becomes irrelevant
928 ! by setting cqs2=chs in this situation the 2m q should become just qv(k=1)
929 IF (Q1 .GT. QSFC(I,J)) THEN
930 CQS2(I,J) = CHS(I,J)
931 ENDIF
932 ! QSFC(I,J)=Q1
933 ! Convert QSFC back to mixing ratio
934 QSFC(I,J)= Q1/(1.0-Q1)
935 !
936 QSFC_RURAL(I,J)= Q1/(1.0-Q1)
937 ! Calculate momentum flux from rural surface for use with multi-layer UCM (Martilli et al. 2002)
939 DO 80 NS=1,NSOIL
940 SMOIS(I,NS,J)=SMC(NS)
941 TSLB(I,NS,J)=STC(NS) ! STEMP
942 SH2O(I,NS,J)=SWC(NS)
943 80 CONTINUE
944 ! ENDIF
946 !
947 ! Residual of surface energy balance equation terms
948 !
949 noahres(i,j) = ( solnet + lwdn ) - sheat + ssoil - eta - ( emissi * STBOLT * (t1**4) ) - flx1 - flx2 - flx3
952 IF (SF_URBAN_PHYSICS == 1 ) THEN ! Beginning of UCM CALL if block
953 !--------------------------------------
954 ! URBAN CANOPY MODEL START - urban
955 !--------------------------------------
956 ! Input variables lsm --> urban
959 IF( IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == 31 .or. &
960 IVGTYP(I,J) == 32 .or. IVGTYP(I,J) == 33 ) THEN
962 ! Call urban
964 !
965 UTYPE_URB = UTYPE_URB2D(I,J) !urban type (low, high or industrial)
967 TA_URB = SFCTMP ! [K]
968 QA_URB = Q2K ! [kg/kg]
969 UA_URB = SQRT(U_PHY(I,1,J)**2.+V_PHY(I,1,J)**2.)
970 U1_URB = U_PHY(I,1,J)
971 V1_URB = V_PHY(I,1,J)
972 IF(UA_URB < 1.) UA_URB=1. ! [m/s]
973 SSG_URB = SOLDN ! [W/m/m]
974 SSGD_URB = 0.8*SOLDN ! [W/m/m]
975 SSGQ_URB = SSG_URB-SSGD_URB ! [W/m/m]
976 LLG_URB = GLW(I,J) ! [W/m/m]
977 RAIN_URB = RAINBL(I,J) ! [mm]
978 RHOO_URB = SFCPRS / (287.04 * SFCTMP * (1.0+ 0.61 * Q2K)) ![kg/m/m/m]
979 ZA_URB = ZLVL ! [m]
980 DELT_URB = DT ! [sec]
981 XLAT_URB = XLAT_URB2D(I,J) ! [deg]
982 COSZ_URB = COSZ_URB2D(I,J) !
983 OMG_URB = OMG_URB2D(I,J) !
984 ZNT_URB = ZNT(I,J)
986 LSOLAR_URB = .FALSE.
988 TR_URB = TR_URB2D(I,J)
989 TB_URB = TB_URB2D(I,J)
990 TG_URB = TG_URB2D(I,J)
991 TC_URB = TC_URB2D(I,J)
992 QC_URB = QC_URB2D(I,J)
993 UC_URB = UC_URB2D(I,J)
995 DO K = 1,num_roof_layers
996 TRL_URB(K) = TRL_URB3D(I,K,J)
997 END DO
998 DO K = 1,num_wall_layers
999 TBL_URB(K) = TBL_URB3D(I,K,J)
1000 END DO
1001 DO K = 1,num_road_layers
1002 TGL_URB(K) = TGL_URB3D(I,K,J)
1003 END DO
1005 XXXR_URB = XXXR_URB2D(I,J)
1006 XXXB_URB = XXXB_URB2D(I,J)
1007 XXXG_URB = XXXG_URB2D(I,J)
1008 XXXC_URB = XXXC_URB2D(I,J)
1009 !
1010 !
1011 ! Limits to avoid dividing by small number
1012 if (CHS(I,J) < 1.0E-02) then
1013 CHS(I,J) = 1.0E-02
1014 endif
1015 if (CHS2(I,J) < 1.0E-02) then
1016 CHS2(I,J) = 1.0E-02
1017 endif
1018 if (CQS2(I,J) < 1.0E-02) then
1019 CQS2(I,J) = 1.0E-02
1020 endif
1021 !
1022 CHS_URB = CHS(I,J)
1023 CHS2_URB = CHS2(I,J)
1024 IF (PRESENT(CMR_SFCDIF)) THEN
1025 CMR_URB = CMR_SFCDIF(I,J)
1026 CHR_URB = CHR_SFCDIF(I,J)
1027 CMC_URB = CMC_SFCDIF(I,J)
1028 CHC_URB = CHC_SFCDIF(I,J)
1029 ENDIF
1030 !
1031 ! Call urban
1033 CALL urban(LSOLAR_URB, & ! I
1034 num_roof_layers,num_wall_layers,num_road_layers, & ! C
1035 DZR,DZB,DZG, & ! C
1036 UTYPE_URB,TA_URB,QA_URB,UA_URB,U1_URB,V1_URB,SSG_URB, & ! I
1037 SSGD_URB,SSGQ_URB,LLG_URB,RAIN_URB,RHOO_URB, & ! I
1038 ZA_URB,DECLIN_URB,COSZ_URB,OMG_URB, & ! I
1039 XLAT_URB,DELT_URB,ZNT_URB, & ! I
1040 CHS_URB, CHS2_URB, & ! I
1041 TR_URB, TB_URB, TG_URB, TC_URB, QC_URB,UC_URB, & ! H
1042 TRL_URB,TBL_URB,TGL_URB, & ! H
1043 XXXR_URB, XXXB_URB, XXXG_URB, XXXC_URB, & ! H
1044 TS_URB,QS_URB,SH_URB,LH_URB,LH_KINEMATIC_URB, & ! O
1045 SW_URB,ALB_URB,LW_URB,G_URB,RN_URB,PSIM_URB,PSIH_URB, & ! O
1046 GZ1OZ0_URB, & !O
1047 CMR_URB, CHR_URB, CMC_URB, CHC_URB, &
1048 U10_URB, V10_URB, TH2_URB, Q2_URB, & ! O
1049 UST_URB) !O
1051 #if 0
1052 IF(IPRINT) THEN
1054 print*, 'AFTER CALL URBAN'
1055 print*,'num_roof_layers',num_roof_layers, 'num_wall_layers', &
1056 num_wall_layers, &
1057 'DZR',DZR,'DZB',DZB,'DZG',DZG,'UTYPE_URB',UTYPE_URB,'TA_URB', &
1058 TA_URB, &
1059 'QA_URB',QA_URB,'UA_URB',UA_URB,'U1_URB',U1_URB,'V1_URB', &
1060 V1_URB, &
1061 'SSG_URB',SSG_URB,'SSGD_URB',SSGD_URB,'SSGQ_URB',SSGQ_URB, &
1062 'LLG_URB',LLG_URB,'RAIN_URB',RAIN_URB,'RHOO_URB',RHOO_URB, &
1063 'ZA_URB',ZA_URB, 'DECLIN_URB',DECLIN_URB,'COSZ_URB',COSZ_URB,&
1064 'OMG_URB',OMG_URB,'XLAT_URB',XLAT_URB,'DELT_URB',DELT_URB, &
1065 'ZNT_URB',ZNT_URB,'TR_URB',TR_URB, 'TB_URB',TB_URB,'TG_URB',&
1066 TG_URB,'TC_URB',TC_URB,'QC_URB',QC_URB,'TRL_URB',TRL_URB, &
1067 'TBL_URB',TBL_URB,'TGL_URB',TGL_URB,'XXXR_URB',XXXR_URB, &
1068 'XXXB_URB',XXXB_URB,'XXXG_URB',XXXG_URB,'XXXC_URB',XXXC_URB,&
1069 'TS_URB',TS_URB,'QS_URB',QS_URB,'SH_URB',SH_URB,'LH_URB', &
1070 LH_URB, 'LH_KINEMATIC_URB',LH_KINEMATIC_URB,'SW_URB',SW_URB,&
1071 'ALB_URB',ALB_URB,'LW_URB',LW_URB,'G_URB',G_URB,'RN_URB', &
1072 RN_URB, 'PSIM_URB',PSIM_URB,'PSIH_URB',PSIH_URB, &
1073 'U10_URB',U10_URB,'V10_URB',V10_URB,'TH2_URB',TH2_URB, &
1074 'Q2_URB',Q2_URB,'CHS_URB',CHS_URB,'CHS2_URB',CHS2_URB
1075 endif
1076 #endif
1078 TS_URB2D(I,J) = TS_URB
1080 ALBEDO(I,J) = FRC_URB2D(I,J)*ALB_URB+(1-FRC_URB2D(I,J))*ALBEDOK ![-]
1081 HFX(I,J) = FRC_URB2D(I,J)*SH_URB+(1-FRC_URB2D(I,J))*SHEAT ![W/m/m]
1082 QFX(I,J) = FRC_URB2D(I,J)*LH_KINEMATIC_URB &
1083 + (1-FRC_URB2D(I,J))*ETA_KINEMATIC ![kg/m/m/s]
1084 LH(I,J) = FRC_URB2D(I,J)*LH_URB+(1-FRC_URB2D(I,J))*ETA ![W/m/m]
1085 GRDFLX(I,J) = FRC_URB2D(I,J)*G_URB+(1-FRC_URB2D(I,J))*SSOIL ![W/m/m]
1086 TSK(I,J) = FRC_URB2D(I,J)*TS_URB+(1-FRC_URB2D(I,J))*T1 ![K]
1087 Q1 = FRC_URB2D(I,J)*QS_URB+(1-FRC_URB2D(I,J))*Q1 ![-]
1088 ! Convert QSFC back to mixing ratio
1089 QSFC(I,J)= Q1/(1.0-Q1)
1090 UST(I,J)= FRC_URB2D(I,J)*UST_URB+(1-FRC_URB2D(I,J))*UST(I,J) ![m/s]
1092 #if 0
1093 IF(IPRINT)THEN
1095 print*, ' FRC_URB2D', FRC_URB2D, &
1096 'ALB_URB',ALB_URB, 'ALBEDOK',ALBEDOK, &
1097 'ALBEDO(I,J)', ALBEDO(I,J), &
1098 'SH_URB',SH_URB,'SHEAT',SHEAT, 'HFX(I,J)',HFX(I,J), &
1099 'LH_KINEMATIC_URB',LH_KINEMATIC_URB,'ETA_KINEMATIC', &
1100 ETA_KINEMATIC, 'QFX(I,J)',QFX(I,J), &
1101 'LH_URB',LH_URB, 'ETA',ETA, 'LH(I,J)',LH(I,J), &
1102 'G_URB',G_URB,'SSOIL',SSOIL,'GRDFLX(I,J)', GRDFLX(I,J),&
1103 'TS_URB',TS_URB,'T1',T1,'TSK(I,J)',TSK(I,J), &
1104 'QS_URB',QS_URB,'Q1',Q1,'QSFC(I,J)',QSFC(I,J)
1105 endif
1106 #endif
1110 ! Renew Urban State Varialbes
1112 TR_URB2D(I,J) = TR_URB
1113 TB_URB2D(I,J) = TB_URB
1114 TG_URB2D(I,J) = TG_URB
1115 TC_URB2D(I,J) = TC_URB
1116 QC_URB2D(I,J) = QC_URB
1117 UC_URB2D(I,J) = UC_URB
1119 DO K = 1,num_roof_layers
1120 TRL_URB3D(I,K,J) = TRL_URB(K)
1121 END DO
1122 DO K = 1,num_wall_layers
1123 TBL_URB3D(I,K,J) = TBL_URB(K)
1124 END DO
1125 DO K = 1,num_road_layers
1126 TGL_URB3D(I,K,J) = TGL_URB(K)
1127 END DO
1128 XXXR_URB2D(I,J) = XXXR_URB
1129 XXXB_URB2D(I,J) = XXXB_URB
1130 XXXG_URB2D(I,J) = XXXG_URB
1131 XXXC_URB2D(I,J) = XXXC_URB
1133 SH_URB2D(I,J) = SH_URB
1134 LH_URB2D(I,J) = LH_URB
1135 G_URB2D(I,J) = G_URB
1136 RN_URB2D(I,J) = RN_URB
1137 PSIM_URB2D(I,J) = PSIM_URB
1138 PSIH_URB2D(I,J) = PSIH_URB
1139 GZ1OZ0_URB2D(I,J)= GZ1OZ0_URB
1140 U10_URB2D(I,J) = U10_URB
1141 V10_URB2D(I,J) = V10_URB
1142 TH2_URB2D(I,J) = TH2_URB
1143 Q2_URB2D(I,J) = Q2_URB
1144 UST_URB2D(I,J) = UST_URB
1145 AKMS_URB2D(I,J) = KARMAN * UST_URB2D(I,J)/(GZ1OZ0_URB2D(I,J)-PSIM_URB2D(I,J))
1146 IF (PRESENT(CMR_SFCDIF)) THEN
1147 CMR_SFCDIF(I,J) = CMR_URB
1148 CHR_SFCDIF(I,J) = CHR_URB
1149 CMC_SFCDIF(I,J) = CMC_URB
1150 CHC_SFCDIF(I,J) = CHC_URB
1151 ENDIF
1152 END IF
1154 ENDIF ! end of UCM CALL if block
1155 !--------------------------------------
1156 ! Urban Part End - urban
1157 !--------------------------------------
1159 !*** DIAGNOSTICS
1160 SMSTAV(I,J)=SOILW
1161 SMSTOT(I,J)=SOILM*1000.
1162 DO NS=1,NSOIL
1163 SMCREL(I,NS,J)=SMAV(NS)
1164 ENDDO
1165 ! Convert the water unit into mm
1166 SFCRUNOFF(I,J)=SFCRUNOFF(I,J)+RUNOFF1*DT*1000.0
1167 UDRUNOFF(I,J)=UDRUNOFF(I,J)+RUNOFF2*DT*1000.0
1168 ! snow defined when fraction of frozen precip (FFROZP) > 0.5,
1169 IF(FFROZP.GT.0.5)THEN
1170 ACSNOW(I,J)=ACSNOW(I,J)+PRCP*DT
1171 ENDIF
1172 IF(SNOW(I,J).GT.0.)THEN
1173 ACSNOM(I,J)=ACSNOM(I,J)+SNOMLT*1000.
1174 ! accumulated snow-melt energy
1175 SNOPCX(I,J)=SNOPCX(I,J)-SNOMLT/FDTLIW
1176 ENDIF
1178 ENDIF ! endif of land-sea test
1180 ENDDO ILOOP ! of I loop
1181 ENDDO JLOOP ! of J loop
1183 IF (SF_URBAN_PHYSICS == 2) THEN
1186 do j=jts,jte
1187 do i=its,ite
1188 EMISS_URB(i,j)=0.
1189 RL_UP_URB(i,j)=0.
1190 RS_ABS_URB(i,j)=0.
1191 GRDFLX_URB(i,j)=0.
1192 b_q_bep(i,kts:kte,j)=0.
1193 end do
1194 end do
1195 CALL BEP(frc_urb2d,utype_urb2d,itimestep,dz8w,dt,u_phy,v_phy, &
1196 th_phy,rho,p_phy,swdown,glw, &
1197 gmt,julday,xlong,xlat,declin_urb,cosz_urb2d,omg_urb2d, &
1198 num_urban_layers, &
1199 trb_urb4d,tw1_urb4d,tw2_urb4d,tgb_urb4d, &
1200 sfw1_urb3d,sfw2_urb3d,sfr_urb3d,sfg_urb3d, &
1201 a_u_bep,a_v_bep,a_t_bep, &
1202 a_e_bep,b_u_bep,b_v_bep, &
1203 b_t_bep,b_e_bep,dlg_bep, &
1204 dl_u_bep,sf_bep,vl_bep, &
1205 rl_up_urb,rs_abs_urb,emiss_urb,grdflx_urb, &
1206 ids,ide, jds,jde, kds,kde, &
1207 ims,ime, jms,jme, kms,kme, &
1208 its,ite, jts,jte, kts,kte )
1210 ENDIF
1213 IF (SF_URBAN_PHYSICS == 3) THEN
1216 do j=jts,jte
1217 do i=its,ite
1218 EMISS_URB(i,j)=0.
1219 RL_UP_URB(i,j)=0.
1220 RS_ABS_URB(i,j)=0.
1221 GRDFLX_URB(i,j)=0.
1222 b_q_bep(i,kts:kte,j)=0.
1223 end do
1224 end do
1226 CALL BEP_BEM(frc_urb2d,utype_urb2d,itimestep,dz8w,dt,u_phy,v_phy, &
1227 th_phy,rho,p_phy,swdown,glw, &
1228 gmt,julday,xlong,xlat,declin_urb,cosz_urb2d,omg_urb2d, &
1229 num_urban_layers, &
1230 trb_urb4d,tw1_urb4d,tw2_urb4d,tgb_urb4d, &
1231 tlev_urb3d,qlev_urb3d,tw1lev_urb3d,tw2lev_urb3d, &
1232 tglev_urb3d,tflev_urb3d,sf_ac_urb3d,lf_ac_urb3d, &
1233 cm_ac_urb3d,sfvent_urb3d,lfvent_urb3d, &
1234 sfwin1_urb3d,sfwin2_urb3d, &
1235 sfw1_urb3d,sfw2_urb3d,sfr_urb3d,sfg_urb3d, &
1236 a_u_bep,a_v_bep,a_t_bep, &
1237 a_e_bep,b_u_bep,b_v_bep, &
1238 b_t_bep,b_e_bep,b_q_bep,dlg_bep, &
1239 dl_u_bep,sf_bep,vl_bep, &
1240 rl_up_urb,rs_abs_urb,emiss_urb,grdflx_urb,qv3d, &
1241 ids,ide, jds,jde, kds,kde, &
1242 ims,ime, jms,jme, kms,kme, &
1243 its,ite, jts,jte, kts,kte )
1245 ENDIF
1247 if((sf_urban_physics.eq.2).OR.(sf_urban_physics.eq.3))then !Bep begin
1248 ! fix the value of the Stefan-Boltzmann constant
1249 sigma_sb=5.67e-08
1250 do j=jts,jte
1251 do i=its,ite
1252 UMOM_URB(I,J)=0.
1253 VMOM_URB(I,J)=0.
1254 HFX_URB(I,J)=0.
1255 QFX_URB(I,J)=0.
1256 do k=kts,kte
1257 a_u_bep(i,k,j)=a_u_bep(i,k,j)*frc_urb2d(i,j)
1258 a_v_bep(i,k,j)=a_v_bep(i,k,j)*frc_urb2d(i,j)
1259 a_t_bep(i,k,j)=a_t_bep(i,k,j)*frc_urb2d(i,j)
1260 a_q_bep(i,k,j)=0.
1261 a_e_bep(i,k,j)=0.
1262 b_u_bep(i,k,j)=b_u_bep(i,k,j)*frc_urb2d(i,j)
1263 b_v_bep(i,k,j)=b_v_bep(i,k,j)*frc_urb2d(i,j)
1264 b_t_bep(i,k,j)=b_t_bep(i,k,j)*frc_urb2d(i,j)
1265 b_q_bep(i,k,j)=b_q_bep(i,k,j)*frc_urb2d(i,j)
1266 b_e_bep(i,k,j)=b_e_bep(i,k,j)*frc_urb2d(i,j)
1267 HFX_URB(I,J)=HFX_URB(I,J)+B_T_BEP(I,K,J)*RHO(I,K,J)*CP* &
1268 DZ8W(I,K,J)*VL_BEP(I,K,J)
1269 QFX_URB(I,J)=QFX_URB(I,J)+B_Q_BEP(I,K,J)* &
1270 DZ8W(I,K,J)*VL_BEP(I,K,J)
1271 UMOM_URB(I,J)=UMOM_URB(I,J)+ (A_U_BEP(I,K,J)*U_PHY(I,K,J)+ &
1272 B_U_BEP(I,K,J))*DZ8W(I,K,J)*VL_BEP(I,K,J)
1273 VMOM_URB(I,J)=VMOM_URB(I,J)+ (A_V_BEP(I,K,J)*V_PHY(I,K,J)+ &
1274 B_V_BEP(I,K,J))*DZ8W(I,K,J)*VL_BEP(I,K,J)
1275 vl_bep(i,k,j)=(1.-frc_urb2d(i,j))+vl_bep(i,k,j)*frc_urb2d(i,j)
1276 sf_bep(i,k,j)=(1.-frc_urb2d(i,j))+sf_bep(i,k,j)*frc_urb2d(i,j)
1277 end do
1278 a_u_bep(i,1,j)=(1.-frc_urb2d(i,j))*(-ust(I,J)*ust(I,J))/dz8w(i,1,j)/ &
1279 ((u_phy(i,1,j)**2+v_phy(i,1,j)**2.)**.5)+a_u_bep(i,1,j)
1280 a_v_bep(i,1,j)=(1.-frc_urb2d(i,j))*(-ust(I,J)*ust(I,J))/dz8w(i,1,j)/ &
1281 ((u_phy(i,1,j)**2+v_phy(i,1,j)**2.)**.5)+a_v_bep(i,1,j)
1282 b_t_bep(i,1,j)=(1.-frc_urb2d(i,j))*hfx_rural(i,j)/dz8w(i,1,j)/rho(i,1,j)/CP+ &
1283 b_t_bep(i,1,j)
1284 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)
1285 umom=(1.-frc_urb2d(i,j))*ust(i,j)*ust(i,j)*u_phy(i,1,j)/ &
1286 ((u_phy(i,1,j)**2+v_phy(i,1,j)**2.)**.5)+umom_urb(i,j)
1287 vmom=(1.-frc_urb2d(i,j))*ust(i,j)*ust(i,j)*v_phy(i,1,j)/ &
1288 ((u_phy(i,1,j)**2+v_phy(i,1,j)**2.)**.5)+vmom_urb(i,j)
1289 sf_bep(i,1,j)=1.
1291 ! compute upward longwave radiation from the rural part and total
1292 ! rl_up_rural=-emiss_rural(i,j)*sigma_sb*(tsk_rural(i,j)**4.)-(1.-emiss_rural(i,j))*glw(i,j)
1293 ! rl_up_tot=(1.-frc_urb2d(i,j))*rl_up_rural+frc_urb2d(i,j)*rl_up_urb(i,j)
1294 ! emiss(i,j)=(1.-frc_urb2d(i,j))*emiss_rural(i,j)+frc_urb2d(i,j)*emiss_urb(i,j)
1295 ! using the emissivity and the total longwave upward radiation estimate the averaged skin temperature
1296 IF (FRC_URB2D(I,J).GT.0.) THEN
1297 rl_up_rural=-emiss_rural(i,j)*sigma_sb*(tsk_rural(i,j)**4.)-(1.-emiss_rural(i,j))*glw(i,j)
1298 rl_up_tot=(1.-frc_urb2d(i,j))*rl_up_rural+frc_urb2d(i,j)*rl_up_urb(i,j)
1299 emiss(i,j)=(1.-frc_urb2d(i,j))*emiss_rural(i,j)+frc_urb2d(i,j)*emiss_urb(i,j)
1300 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
1301 tsk(i,j)=(max(0., (-1.*rl_up_tot-(1.-emiss(i,j))*glw(i,j) )/emiss(i,j)/sigma_sb))**.25
1302 rs_abs_tot=(1.-frc_urb2d(i,j))*swdown(i,j)*(1.-albedo(i,j))+frc_urb2d(i,j)*rs_abs_urb(i,j)
1303 if(swdown(i,j).gt.0.)then
1304 albedo(i,j)=1.-rs_abs_tot/swdown(i,j)
1305 else
1306 albedo(i,j)=alb_rural(i,j)
1307 endif
1308 ! rename *_urb to sh_urb2d,lh_urb2d,g_urb2d,rn_urb2d
1309 grdflx(i,j)= (1.-frc_urb2d(i,j))*grdflx_rural(i,j)+frc_urb2d(i,j)*grdflx_urb(i,j)
1310 qfx(i,j)=(1.-frc_urb2d(i,j))*qfx_rural(i,j)+qfx_urb(i,j)
1311 ! lh(i,j)=(1.-frc_urb2d(i,j))*qfx_rural(i,j)*xlv
1312 lh(i,j)=qfx(i,j)*xlv
1313 HFX(I,J) = HFX_URB(I,J)+(1-FRC_URB2D(I,J))*HFX_RURAL(I,J) ![W/m/m]
1314 SH_URB2D(I,J) = HFX_URB(I,J)/FRC_URB2D(I,J)
1315 LH_URB2D(I,J) = qfx_urb(i,j)*xlv
1316 G_URB2D(I,J) = grdflx_urb(i,j)
1317 RN_URB2D(I,J) = rs_abs_urb(i,j)+emiss_urb(i,j)*glw(i,j)-rl_up_urb(i,j)
1318 ust(i,j)=(umom**2.+vmom**2.)**.25
1319 ! if(tsk(i,j).gt.350)write(*,*)'tsk too big!',i,j,tsk(i,j)
1320 ! 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
1321 ! print*,'ivgtyp,i,j,sigma_sb',ivgtyp(i,j),i,j,sigma_sb
1322 ! print*,'hfx,lh,qfx,grdflx,ts_urb2d',hfx(i,j),lh(i,j),qfx(i,j),grdflx(i,j),ts_urb2d(i,j)
1323 ! print*,'tsk,albedo,emiss',tsk(i,j),albedo(i,j),emiss(i,j)
1324 ! if(i.eq.56.and.j.eq.29)then
1325 ! print*,'ivgtyp, qfx, hfx',ivgtyp(i,j),hfx_rural(i,j),qfx_rural(i,j)
1326 ! print*,'emiss_rural,emiss_urb',emiss_rural(i,j),emiss_urb(i,j)
1327 ! print*,'rl_up_rural,rl_up_urb(i,j)',rl_up_rural,rl_up_urb(i,j)
1328 ! print*,'tsk_rural,ts_urb2d(i,j),tsk',tsk_rural(i,j),ts_urb2d(i,j),tsk(i,j)
1329 ! 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
1330 ! print*,'ivgtyp,hfx,hfx_urb,hfx_rural',hfx(i,j),hfx_urb(i,j),hfx_rural(i,j)
1331 ! print*,'lh,lh_rural',lh(i,j),lh_rural(i,j)
1332 ! print*,'qfx',qfx(i,j)
1333 ! print*,'ts_urb2d',ts_urb2d(i,j)
1334 ! print*,'ust',ust(i,j)
1335 ! print*,'swdown,glw',swdown(i,j),glw(i,j)
1336 ! endif
1337 else
1338 SH_URB2D(I,J) = 0.
1339 LH_URB2D(I,J) = 0.
1340 G_URB2D(I,J) = 0.
1341 RN_URB2D(I,J) = 0.
1342 endif
1343 ! IF( IVGTYP(I,J) == 1 .or. IVGTYP(I,J) == 31 .or. &
1344 ! IVGTYP(I,J) == 32 .or. IVGTYP(I,J) == 33) THEN
1345 ! print*,'ivgtyp, qfx, hfx',ivgtyp(i,j),hfx_rural(i,j),qfx_rural(i,j)
1346 ! print*,'ivgtyp,hfx,hfx_urb,hfx_rural',hfx(i,j),hfx_urb(i,j),hfx_rural(i,j)
1347 ! print*,'lh,lh_rural',lh(i,j),lh_rural(i,j)
1348 ! print*,'qfx',qfx(i,j)
1349 ! print*,'ts_urb2d',ts_urb2d(i,j)
1350 ! print*,'ust',ust(i,j)
1351 ! endif
1352 enddo
1353 enddo
1356 endif !Bep end
1358 !------------------------------------------------------
1359 END SUBROUTINE lsm
1360 !------------------------------------------------------
1362 SUBROUTINE LSMINIT(VEGFRA,SNOW,SNOWC,SNOWH,CANWAT,SMSTAV, &
1363 SMSTOT, SFCRUNOFF,UDRUNOFF,ACSNOW, &
1364 ACSNOM,IVGTYP,ISLTYP,TSLB,SMOIS,SH2O,ZS,DZS, &
1365 MMINLU, &
1366 SNOALB, FNDSOILW, FNDSNOWH, RDMAXALB, &
1367 num_soil_layers, restart, &
1368 allowed_to_read , &
1369 ids,ide, jds,jde, kds,kde, &
1370 ims,ime, jms,jme, kms,kme, &
1371 its,ite, jts,jte, kts,kte )
1373 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
1374 ims,ime, jms,jme, kms,kme, &
1375 its,ite, jts,jte, kts,kte
1377 INTEGER, INTENT(IN) :: num_soil_layers
1379 LOGICAL , INTENT(IN) :: restart , allowed_to_read
1381 REAL, DIMENSION( num_soil_layers), INTENT(INOUT) :: ZS, DZS
1383 REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , &
1384 INTENT(INOUT) :: SMOIS, & !Total soil moisture
1385 SH2O, & !liquid soil moisture
1386 TSLB !STEMP
1388 REAL, DIMENSION( ims:ime, jms:jme ) , &
1389 INTENT(INOUT) :: SNOW, &
1390 SNOWH, &
1391 SNOWC, &
1392 SNOALB, &
1393 CANWAT, &
1394 SMSTAV, &
1395 SMSTOT, &
1396 SFCRUNOFF, &
1397 UDRUNOFF, &
1398 ACSNOW, &
1399 VEGFRA, &
1400 ACSNOM
1402 INTEGER, DIMENSION( ims:ime, jms:jme ) , &
1403 INTENT(IN) :: IVGTYP, &
1404 ISLTYP
1405 CHARACTER(LEN=*), INTENT(IN) :: MMINLU
1407 LOGICAL, INTENT(IN) :: FNDSOILW , &
1408 FNDSNOWH
1409 LOGICAL, INTENT(IN) :: RDMAXALB
1412 INTEGER :: L
1413 REAL :: BX, SMCMAX, PSISAT, FREE
1414 REAL, PARAMETER :: BLIM = 5.5, HLICE = 3.335E5, &
1415 GRAV = 9.81, T0 = 273.15
1416 INTEGER :: errflag
1418 character*256 :: MMINSL
1419 MMINSL='STAS'
1420 !
1422 ! initialize three Noah LSM related tables
1423 IF ( allowed_to_read ) THEN
1424 CALL wrf_message( 'INITIALIZE THREE Noah LSM RELATED TABLES' )
1425 CALL SOIL_VEG_GEN_PARM( MMINLU, MMINSL )
1426 ENDIF
1428 #ifdef WRF_CHEM
1429 !
1430 ! need this parameter for dust parameterization in wrf/chem
1431 !
1432 do I=1,NSLTYPE
1433 porosity(i)=maxsmc(i)
1434 drypoint(i)=drysmc(i)
1435 enddo
1436 #endif
1438 IF(.not.restart)THEN
1440 itf=min0(ite,ide-1)
1441 jtf=min0(jte,jde-1)
1443 errflag = 0
1444 DO j = jts,jtf
1445 DO i = its,itf
1446 IF ( ISLTYP( i,j ) .LT. 1 ) THEN
1447 errflag = 1
1448 WRITE(err_message,*)"module_sf_noahlsm.F: lsminit: out of range ISLTYP ",i,j,ISLTYP( i,j )
1449 CALL wrf_message(err_message)
1450 ENDIF
1451 IF(.not.RDMAXALB) THEN
1452 SNOALB(i,j)=MAXALB(IVGTYP(i,j))*0.01
1453 ENDIF
1454 ENDDO
1455 ENDDO
1456 IF ( errflag .EQ. 1 ) THEN
1457 CALL wrf_error_fatal( "module_sf_noahlsm.F: lsminit: out of range value "// &
1458 "of ISLTYP. Is this field in the input?" )
1459 ENDIF
1461 ! initialize soil liquid water content SH2O
1463 ! IF(.NOT.FNDSOILW) THEN
1465 ! If no SWC, do the following
1466 ! PRINT *,'SOIL WATER NOT FOUND - VALUE SET IN LSMINIT'
1467 DO J = jts,jtf
1468 DO I = its,itf
1469 BX = BB(ISLTYP(I,J))
1470 SMCMAX = MAXSMC(ISLTYP(I,J))
1471 PSISAT = SATPSI(ISLTYP(I,J))
1472 if ((bx > 0.0).and.(smcmax > 0.0).and.(psisat > 0.0)) then
1473 DO NS=1, num_soil_layers
1474 ! ----------------------------------------------------------------------
1475 !SH2O <= SMOIS for T < 273.149K (-0.001C)
1476 IF (TSLB(I,NS,J) < 273.149) THEN
1477 ! ----------------------------------------------------------------------
1478 ! first guess following explicit solution for Flerchinger Eqn from Koren
1479 ! et al, JGR, 1999, Eqn 17 (KCOUNT=0 in FUNCTION FRH2O).
1480 ! ISLTPK is soil type
1481 BX = BB(ISLTYP(I,J))
1482 SMCMAX = MAXSMC(ISLTYP(I,J))
1483 PSISAT = SATPSI(ISLTYP(I,J))
1484 IF ( BX > BLIM ) BX = BLIM
1485 FK=(( (HLICE/(GRAV*(-PSISAT))) * &
1486 ((TSLB(I,NS,J)-T0)/TSLB(I,NS,J)) )**(-1/BX) )*SMCMAX
1487 IF (FK < 0.02) FK = 0.02
1488 SH2O(I,NS,J) = MIN( FK, SMOIS(I,NS,J) )
1489 ! ----------------------------------------------------------------------
1490 ! now use iterative solution for liquid soil water content using
1491 ! FUNCTION FRH2O with the initial guess for SH2O from above explicit
1492 ! first guess.
1493 CALL FRH2O (FREE,TSLB(I,NS,J),SMOIS(I,NS,J),SH2O(I,NS,J), &
1494 SMCMAX,BX,PSISAT)
1495 SH2O(I,NS,J) = FREE
1496 ELSE ! of IF (TSLB(I,NS,J)
1497 ! ----------------------------------------------------------------------
1498 ! SH2O = SMOIS ( for T => 273.149K (-0.001C)
1499 SH2O(I,NS,J)=SMOIS(I,NS,J)
1500 ! ----------------------------------------------------------------------
1501 ENDIF ! of IF (TSLB(I,NS,J)
1502 END DO ! of DO NS=1, num_soil_layers
1503 else ! of if ((bx > 0.0)
1504 DO NS=1, num_soil_layers
1505 SH2O(I,NS,J)=SMOIS(I,NS,J)
1506 END DO
1507 endif ! of if ((bx > 0.0)
1508 ENDDO ! DO I = its,itf
1509 ENDDO ! DO J = jts,jtf
1510 ! ENDIF ! of IF(.NOT.FNDSOILW)THEN
1512 ! initialize physical snow height SNOWH
1514 IF(.NOT.FNDSNOWH)THEN
1515 ! If no SNOWH do the following
1516 CALL wrf_message( 'SNOW HEIGHT NOT FOUND - VALUE DEFINED IN LSMINIT' )
1517 DO J = jts,jtf
1518 DO I = its,itf
1519 SNOWH(I,J)=SNOW(I,J)*0.005 ! SNOW in mm and SNOWH in m
1520 ENDDO
1521 ENDDO
1522 ENDIF
1524 ! initialize canopy water to ZERO
1526 ! GO TO 110
1527 ! print*,'Note that canopy water content (CANWAT) is set to ZERO in LSMINIT'
1528 DO J = jts,jtf
1529 DO I = its,itf
1530 CANWAT(I,J)=0.0
1531 ENDDO
1532 ENDDO
1533 110 CONTINUE
1535 ENDIF
1536 !------------------------------------------------------------------------------
1537 END SUBROUTINE lsminit
1538 !------------------------------------------------------------------------------
1542 !-----------------------------------------------------------------
1543 SUBROUTINE SOIL_VEG_GEN_PARM( MMINLU, MMINSL)
1544 !-----------------------------------------------------------------
1546 USE module_wrf_error
1547 IMPLICIT NONE
1549 CHARACTER(LEN=*), INTENT(IN) :: MMINLU, MMINSL
1550 integer :: LUMATCH, IINDEX, LC, NUM_SLOPE
1551 integer :: ierr
1552 INTEGER , PARAMETER :: OPEN_OK = 0
1554 character*128 :: mess , message
1555 logical, external :: wrf_dm_on_monitor
1558 !-----SPECIFY VEGETATION RELATED CHARACTERISTICS :
1559 ! ALBBCK: SFC albedo (in percentage)
1560 ! Z0: Roughness length (m)
1561 ! SHDFAC: Green vegetation fraction (in percentage)
1562 ! Note: The ALBEDO, Z0, and SHDFAC values read from the following table
1563 ! ALBEDO, amd Z0 are specified in LAND-USE TABLE; and SHDFAC is
1564 ! the monthly green vegetation data
1565 ! CMXTBL: MAX CNPY Capacity (m)
1566 ! NROTBL: Rooting depth (layer)
1567 ! RSMIN: Mimimum stomatal resistance (s m-1)
1568 ! RSMAX: Max. stomatal resistance (s m-1)
1569 ! RGL: Parameters used in radiation stress function
1570 ! HS: Parameter used in vapor pressure deficit functio
1571 ! TOPT: Optimum transpiration air temperature. (K)
1572 ! CMCMAX: Maximum canopy water capacity
1573 ! CFACTR: Parameter used in the canopy inteception calculati
1574 ! SNUP: Threshold snow depth (in water equivalent m) that
1575 ! implies 100% snow cover
1576 ! LAI: Leaf area index (dimensionless)
1577 ! MAXALB: Upper bound on maximum albedo over deep snow
1578 !
1579 !-----READ IN VEGETAION PROPERTIES FROM VEGPARM.TBL
1580 !
1582 IF ( wrf_dm_on_monitor() ) THEN
1584 OPEN(19, FILE='VEGPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
1585 IF(ierr .NE. OPEN_OK ) THEN
1586 WRITE(message,FMT='(A)') &
1587 'module_sf_noahlsm.F: soil_veg_gen_parm: failure opening VEGPARM.TBL'
1588 CALL wrf_error_fatal ( message )
1589 END IF
1592 LUMATCH=0
1594 FIND_LUTYPE : DO WHILE (LUMATCH == 0)
1595 READ (19,*,END=2002)
1596 READ (19,*,END=2002)LUTYPE
1597 READ (19,*)LUCATS,IINDEX
1599 IF(LUTYPE.EQ.MMINLU)THEN
1600 WRITE( mess , * ) 'LANDUSE TYPE = ' // TRIM ( LUTYPE ) // ' FOUND', LUCATS,' CATEGORIES'
1601 CALL wrf_message( mess )
1602 LUMATCH=1
1603 ELSE
1604 call wrf_message ( "Skipping over LUTYPE = " // TRIM ( LUTYPE ) )
1605 DO LC = 1, LUCATS+12
1606 read(19,*)
1607 ENDDO
1608 ENDIF
1609 ENDDO FIND_LUTYPE
1610 ! prevent possible array overwrite, Bill Bovermann, IBM, May 6, 2008
1611 IF ( SIZE(SHDTBL) < LUCATS .OR. &
1612 SIZE(NROTBL) < LUCATS .OR. &
1613 SIZE(RSTBL) < LUCATS .OR. &
1614 SIZE(RGLTBL) < LUCATS .OR. &
1615 SIZE(HSTBL) < LUCATS .OR. &
1616 SIZE(SNUPTBL) < LUCATS .OR. &
1617 SIZE(MAXALB) < LUCATS .OR. &
1618 SIZE(LAIMINTBL) < LUCATS .OR. &
1619 SIZE(LAIMAXTBL) < LUCATS .OR. &
1620 SIZE(Z0MINTBL) < LUCATS .OR. &
1621 SIZE(Z0MAXTBL) < LUCATS .OR. &
1622 SIZE(ALBEDOMINTBL) < LUCATS .OR. &
1623 SIZE(ALBEDOMAXTBL) < LUCATS .OR. &
1624 SIZE(EMISSMINTBL ) < LUCATS .OR. &
1625 SIZE(EMISSMAXTBL ) < LUCATS ) THEN
1626 CALL wrf_error_fatal('Table sizes too small for value of LUCATS in module_sf_noahdrv.F')
1627 ENDIF
1629 IF(LUTYPE.EQ.MMINLU)THEN
1630 DO LC=1,LUCATS
1631 READ (19,*)IINDEX,SHDTBL(LC), &
1632 NROTBL(LC),RSTBL(LC),RGLTBL(LC),HSTBL(LC), &
1633 SNUPTBL(LC),MAXALB(LC), LAIMINTBL(LC), &
1634 LAIMAXTBL(LC),EMISSMINTBL(LC), &
1635 EMISSMAXTBL(LC), ALBEDOMINTBL(LC), &
1636 ALBEDOMAXTBL(LC), Z0MINTBL(LC), Z0MAXTBL(LC)
1637 ENDDO
1638 !
1639 READ (19,*)
1640 READ (19,*)TOPT_DATA
1641 READ (19,*)
1642 READ (19,*)CMCMAX_DATA
1643 READ (19,*)
1644 READ (19,*)CFACTR_DATA
1645 READ (19,*)
1646 READ (19,*)RSMAX_DATA
1647 READ (19,*)
1648 READ (19,*)BARE
1649 READ (19,*)
1650 READ (19,*)NATURAL
1651 ENDIF
1652 !
1653 2002 CONTINUE
1655 CLOSE (19)
1656 IF (LUMATCH == 0) then
1657 CALL wrf_error_fatal ("Land Use Dataset '"//MMINLU//"' not found in VEGPARM.TBL.")
1658 ENDIF
1659 ENDIF
1661 CALL wrf_dm_bcast_string ( LUTYPE , 4 )
1662 CALL wrf_dm_bcast_integer ( LUCATS , 1 )
1663 CALL wrf_dm_bcast_integer ( IINDEX , 1 )
1664 CALL wrf_dm_bcast_integer ( LUMATCH , 1 )
1665 CALL wrf_dm_bcast_real ( SHDTBL , NLUS )
1666 CALL wrf_dm_bcast_real ( NROTBL , NLUS )
1667 CALL wrf_dm_bcast_real ( RSTBL , NLUS )
1668 CALL wrf_dm_bcast_real ( RGLTBL , NLUS )
1669 CALL wrf_dm_bcast_real ( HSTBL , NLUS )
1670 CALL wrf_dm_bcast_real ( SNUPTBL , NLUS )
1671 CALL wrf_dm_bcast_real ( LAIMINTBL , NLUS )
1672 CALL wrf_dm_bcast_real ( LAIMAXTBL , NLUS )
1673 CALL wrf_dm_bcast_real ( Z0MINTBL , NLUS )
1674 CALL wrf_dm_bcast_real ( Z0MAXTBL , NLUS )
1675 CALL wrf_dm_bcast_real ( EMISSMINTBL , NLUS )
1676 CALL wrf_dm_bcast_real ( EMISSMAXTBL , NLUS )
1677 CALL wrf_dm_bcast_real ( ALBEDOMINTBL , NLUS )
1678 CALL wrf_dm_bcast_real ( ALBEDOMAXTBL , NLUS )
1679 CALL wrf_dm_bcast_real ( MAXALB , NLUS )
1680 CALL wrf_dm_bcast_real ( TOPT_DATA , 1 )
1681 CALL wrf_dm_bcast_real ( CMCMAX_DATA , 1 )
1682 CALL wrf_dm_bcast_real ( CFACTR_DATA , 1 )
1683 CALL wrf_dm_bcast_real ( RSMAX_DATA , 1 )
1684 CALL wrf_dm_bcast_integer ( BARE , 1 )
1685 CALL wrf_dm_bcast_integer ( NATURAL , 1 )
1687 !
1688 !-----READ IN SOIL PROPERTIES FROM SOILPARM.TBL
1689 !
1690 IF ( wrf_dm_on_monitor() ) THEN
1691 OPEN(19, FILE='SOILPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
1692 IF(ierr .NE. OPEN_OK ) THEN
1693 WRITE(message,FMT='(A)') &
1694 'module_sf_noahlsm.F: soil_veg_gen_parm: failure opening SOILPARM.TBL'
1695 CALL wrf_error_fatal ( message )
1696 END IF
1698 WRITE(mess,*) 'INPUT SOIL TEXTURE CLASSIFICATION = ', TRIM ( MMINSL )
1699 CALL wrf_message( mess )
1701 LUMATCH=0
1703 READ (19,*)
1704 READ (19,2000,END=2003)SLTYPE
1705 2000 FORMAT (A4)
1706 READ (19,*)SLCATS,IINDEX
1707 IF(SLTYPE.EQ.MMINSL)THEN
1708 WRITE( mess , * ) 'SOIL TEXTURE CLASSIFICATION = ', TRIM ( SLTYPE ) , ' FOUND', &
1709 SLCATS,' CATEGORIES'
1710 CALL wrf_message ( mess )
1711 LUMATCH=1
1712 ENDIF
1713 ! prevent possible array overwrite, Bill Bovermann, IBM, May 6, 2008
1714 IF ( SIZE(BB ) < SLCATS .OR. &
1715 SIZE(DRYSMC) < SLCATS .OR. &
1716 SIZE(F11 ) < SLCATS .OR. &
1717 SIZE(MAXSMC) < SLCATS .OR. &
1718 SIZE(REFSMC) < SLCATS .OR. &
1719 SIZE(SATPSI) < SLCATS .OR. &
1720 SIZE(SATDK ) < SLCATS .OR. &
1721 SIZE(SATDW ) < SLCATS .OR. &
1722 SIZE(WLTSMC) < SLCATS .OR. &
1723 SIZE(QTZ ) < SLCATS ) THEN
1724 CALL wrf_error_fatal('Table sizes too small for value of SLCATS in module_sf_noahdrv.F')
1725 ENDIF
1726 IF(SLTYPE.EQ.MMINSL)THEN
1727 DO LC=1,SLCATS
1728 READ (19,*) IINDEX,BB(LC),DRYSMC(LC),F11(LC),MAXSMC(LC),&
1729 REFSMC(LC),SATPSI(LC),SATDK(LC), SATDW(LC), &
1730 WLTSMC(LC), QTZ(LC)
1731 ENDDO
1732 ENDIF
1734 2003 CONTINUE
1736 CLOSE (19)
1737 ENDIF
1739 CALL wrf_dm_bcast_integer ( LUMATCH , 1 )
1740 CALL wrf_dm_bcast_string ( SLTYPE , 4 )
1741 CALL wrf_dm_bcast_string ( MMINSL , 4 ) ! since this is reset above, see oct2 ^
1742 CALL wrf_dm_bcast_integer ( SLCATS , 1 )
1743 CALL wrf_dm_bcast_integer ( IINDEX , 1 )
1744 CALL wrf_dm_bcast_real ( BB , NSLTYPE )
1745 CALL wrf_dm_bcast_real ( DRYSMC , NSLTYPE )
1746 CALL wrf_dm_bcast_real ( F11 , NSLTYPE )
1747 CALL wrf_dm_bcast_real ( MAXSMC , NSLTYPE )
1748 CALL wrf_dm_bcast_real ( REFSMC , NSLTYPE )
1749 CALL wrf_dm_bcast_real ( SATPSI , NSLTYPE )
1750 CALL wrf_dm_bcast_real ( SATDK , NSLTYPE )
1751 CALL wrf_dm_bcast_real ( SATDW , NSLTYPE )
1752 CALL wrf_dm_bcast_real ( WLTSMC , NSLTYPE )
1753 CALL wrf_dm_bcast_real ( QTZ , NSLTYPE )
1755 IF(LUMATCH.EQ.0)THEN
1756 CALL wrf_message( 'SOIl TEXTURE IN INPUT FILE DOES NOT ' )
1757 CALL wrf_message( 'MATCH SOILPARM TABLE' )
1758 CALL wrf_error_fatal ( 'INCONSISTENT OR MISSING SOILPARM FILE' )
1759 ENDIF
1761 !
1762 !-----READ IN GENERAL PARAMETERS FROM GENPARM.TBL
1763 !
1764 IF ( wrf_dm_on_monitor() ) THEN
1765 OPEN(19, FILE='GENPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
1766 IF(ierr .NE. OPEN_OK ) THEN
1767 WRITE(message,FMT='(A)') &
1768 'module_sf_noahlsm.F: soil_veg_gen_parm: failure opening GENPARM.TBL'
1769 CALL wrf_error_fatal ( message )
1770 END IF
1772 READ (19,*)
1773 READ (19,*)
1774 READ (19,*) NUM_SLOPE
1776 SLPCATS=NUM_SLOPE
1777 ! prevent possible array overwrite, Bill Bovermann, IBM, May 6, 2008
1778 IF ( SIZE(slope_data) < NUM_SLOPE ) THEN
1779 CALL wrf_error_fatal('NUM_SLOPE too large for slope_data array in module_sf_noahdrv')
1780 ENDIF
1782 DO LC=1,SLPCATS
1783 READ (19,*)SLOPE_DATA(LC)
1784 ENDDO
1786 READ (19,*)
1787 READ (19,*)SBETA_DATA
1788 READ (19,*)
1789 READ (19,*)FXEXP_DATA
1790 READ (19,*)
1791 READ (19,*)CSOIL_DATA
1792 READ (19,*)
1793 READ (19,*)SALP_DATA
1794 READ (19,*)
1795 READ (19,*)REFDK_DATA
1796 READ (19,*)
1797 READ (19,*)REFKDT_DATA
1798 READ (19,*)
1799 READ (19,*)FRZK_DATA
1800 READ (19,*)
1801 READ (19,*)ZBOT_DATA
1802 READ (19,*)
1803 READ (19,*)CZIL_DATA
1804 READ (19,*)
1805 READ (19,*)SMLOW_DATA
1806 READ (19,*)
1807 READ (19,*)SMHIGH_DATA
1808 READ (19,*)
1809 READ (19,*)LVCOEF_DATA
1810 CLOSE (19)
1811 ENDIF
1813 CALL wrf_dm_bcast_integer ( NUM_SLOPE , 1 )
1814 CALL wrf_dm_bcast_integer ( SLPCATS , 1 )
1815 CALL wrf_dm_bcast_real ( SLOPE_DATA , NSLOPE )
1816 CALL wrf_dm_bcast_real ( SBETA_DATA , 1 )
1817 CALL wrf_dm_bcast_real ( FXEXP_DATA , 1 )
1818 CALL wrf_dm_bcast_real ( CSOIL_DATA , 1 )
1819 CALL wrf_dm_bcast_real ( SALP_DATA , 1 )
1820 CALL wrf_dm_bcast_real ( REFDK_DATA , 1 )
1821 CALL wrf_dm_bcast_real ( REFKDT_DATA , 1 )
1822 CALL wrf_dm_bcast_real ( FRZK_DATA , 1 )
1823 CALL wrf_dm_bcast_real ( ZBOT_DATA , 1 )
1824 CALL wrf_dm_bcast_real ( CZIL_DATA , 1 )
1825 CALL wrf_dm_bcast_real ( SMLOW_DATA , 1 )
1826 CALL wrf_dm_bcast_real ( SMHIGH_DATA , 1 )
1827 CALL wrf_dm_bcast_real ( LVCOEF_DATA , 1 )
1830 !-----------------------------------------------------------------
1831 END SUBROUTINE SOIL_VEG_GEN_PARM
1832 !-----------------------------------------------------------------
1834 END MODULE module_sf_noahdrv