1 MODULE module_sf_noahdrv
3 !-------------------------------
4 USE module_sf_noahlsm, only: SFLX, XLF, XLV, CP, R_D, RHOWATER, NATURAL, SHDTBL, LUTYPE, SLTYPE, STBOLT, &
5 & KARMAN, LUCATS, NROTBL, RSTBL, RGLTBL, HSTBL, SNUPTBL, MAXALB, LAIMINTBL, &
6 & LAIMAXTBL, Z0MINTBL, Z0MAXTBL, ALBEDOMINTBL, ALBEDOMAXTBL, EMISSMINTBL, &
7 & EMISSMAXTBL, TOPT_DATA, CMCMAX_DATA, CFACTR_DATA, RSMAX_DATA, BARE, NLUS, &
8 & SLCATS, BB, DRYSMC, F11, MAXSMC, REFSMC, SATPSI, SATDK, SATDW, WLTSMC, QTZ, &
9 & NSLTYPE, SLPCATS, SLOPE_DATA, SBETA_DATA, FXEXP_DATA, CSOIL_DATA, &
10 & SALP_DATA, REFDK_DATA, REFKDT_DATA, FRZK_DATA, ZBOT_DATA, CZIL_DATA, &
11 & SMLOW_DATA, SMHIGH_DATA, LVCOEF_DATA, NSLOPE, &
12 & FRH2O,ZTOPVTBL,ZBOTVTBL, &
13 & LCZ_1,LCZ_2,LCZ_3,LCZ_4,LCZ_5,LCZ_6,LCZ_7,LCZ_8,LCZ_9,LCZ_10,LCZ_11
14 USE module_sf_urban, only: urban, oasis, IRI_SCHEME
15 USE module_sf_noahlsm_glacial_only, only: sflx_glacial
16 USE module_sf_bep, only: bep
17 USE module_sf_bep_bem, only: bep_bem
19 use mpas_atmphys_date_time, only: cal_mon_day
20 use mpas_atmphys_utilities, only: physics_error_fatal
21 #define FATAL_ERROR(M) call physics_error_fatal( M )
23 use module_ra_gfdleta, only: cal_mon_day
25 #define FATAL_ERROR(M) call wrf_error_fatal( M )
28 USE module_data_gocart_dust
30 !-------------------------------
35 !----------------------------------------------------------------
36 ! Urban related variable are added to arguments - urban
37 !----------------------------------------------------------------
38 SUBROUTINE lsm(DZ8W,QV3D,P8W3D,T3D,TSK, &
39 HFX,QFX,LH,GRDFLX, QGH,GSW,SWDOWN,SWDDIR,SWDDIF,&
41 SFCRUNOFF, UDRUNOFF,IVGTYP,ISLTYP,ISURBAN,ISICE,VEGFRA, &
42 ALBEDO,ALBBCK,ZNT,Z0,TMN,XLAND,XICE,EMISS,EMBCK, &
43 SNOWC,QSFC,RAINBL,MMINLU, &
44 num_soil_layers,DT,DZS,ITIMESTEP, &
45 SMOIS,TSLB,SNOW,CANWAT, &
46 CHS,CHS2,CQS2,CPM,ROVCP,SR,chklowq,lai,qz0, & !H
50 SNOALB,SHDMIN,SHDMAX, & !I
61 ua_phys,flx4_2d,fvb_2d,fbur_2d,fgsn_2d, &
62 ids,ide, jds,jde, kds,kde, &
63 ims,ime, jms,jme, kms,kme, &
64 its,ite, jts,jte, kts,kte, &
66 CMR_SFCDIF,CHR_SFCDIF,CMC_SFCDIF,CHC_SFCDIF, &
67 CMGR_SFCDIF,CHGR_SFCDIF, &
69 TR_URB2D,TB_URB2D,TG_URB2D,TC_URB2D,QC_URB2D, & !H urban
71 XXXR_URB2D,XXXB_URB2D,XXXG_URB2D,XXXC_URB2D, & !H urban
72 TRL_URB3D,TBL_URB3D,TGL_URB3D, & !H urban
73 SH_URB2D,LH_URB2D,G_URB2D,RN_URB2D,TS_URB2D, & !H urban
74 PSIM_URB2D,PSIH_URB2D,U10_URB2D,V10_URB2D, & !O urban
75 GZ1OZ0_URB2D, AKMS_URB2D, & !O urban
76 TH2_URB2D,Q2_URB2D, UST_URB2D, & !O urban
77 DECLIN_URB,COSZ_URB2D,OMG_URB2D, & !I urban
78 XLAT_URB2D, & !I urban
79 num_roof_layers, num_wall_layers, & !I urban
80 num_road_layers, DZR, DZB, DZG, & !I urban
81 CMCR_URB2D,TGR_URB2D,TGRL_URB3D,SMR_URB3D, & !H urban
82 DRELR_URB2D,DRELB_URB2D,DRELG_URB2D, & !H urban
83 FLXHUMR_URB2D,FLXHUMB_URB2D,FLXHUMG_URB2D, & !H urban
84 julian, julyr, & !H urban
85 FRC_URB2D,UTYPE_URB2D, & !O
86 num_urban_ndm, & !I multi-layer urban
87 urban_map_zrd, & !I multi-layer urban
88 urban_map_zwd, & !I multi-layer urban
89 urban_map_gd, & !I multi-layer urban
90 urban_map_zd, & !I multi-layer urban
91 urban_map_zdf, & !I multi-layer urban
92 urban_map_bd, & !I multi-layer urban
93 urban_map_wd, & !I multi-layer urban
94 urban_map_gbd, & !I multi-layer urban
95 urban_map_fbd, & !I multi-layer urban
96 urban_map_zgrd, & !I multi-layer urban
97 num_urban_hi, & !I multi-layer urban
98 tsk_rural_bep, & !H multi-layer urban
99 trb_urb4d,tw1_urb4d,tw2_urb4d,tgb_urb4d, & !H multi-layer urban
100 tlev_urb3d,qlev_urb3d, & !H multi-layer urban
101 tw1lev_urb3d,tw2lev_urb3d, & !H multi-layer urban
102 tglev_urb3d,tflev_urb3d, & !H multi-layer urban
103 sf_ac_urb3d,lf_ac_urb3d,cm_ac_urb3d, & !H multi-layer urban
104 sfvent_urb3d,lfvent_urb3d, & !H multi-layer urban
105 sfwin1_urb3d,sfwin2_urb3d, & !H multi-layer urban
106 sfw1_urb3d,sfw2_urb3d,sfr_urb3d,sfg_urb3d, & !H multi-layer urban
107 ep_pv_urb3d,t_pv_urb3d, & !RMS
108 trv_urb4d,qr_urb4d,qgr_urb3d,tgr_urb3d, & !RMS
109 drain_urb4d,draingr_urb3d,sfrv_urb3d, & !RMS
110 lfrv_urb3d,dgr_urb3d,dg_urb3d,lfr_urb3d,lfg_urb3d,& !RMS
111 lp_urb2d,hi_urb2d,lb_urb2d,hgt_urb2d, & !H multi-layer urban
112 mh_urb2d,stdh_urb2d,lf_urb2d, & !SLUCM
113 th_phy,rho,p_phy,ust, & !I multi-layer urban
114 gmt,julday,xlong,xlat, & !I multi-layer urban
115 a_u_bep,a_v_bep,a_t_bep,a_q_bep, & !O multi-layer urban
116 a_e_bep,b_u_bep,b_v_bep, & !O multi-layer urban
117 b_t_bep,b_q_bep,b_e_bep,dlg_bep, & !O multi-layer urban
118 dl_u_bep,sf_bep,vl_bep &
120 ,sfcheadrt,INFXSRT,soldrain & !O multi-layer urban
122 ,SDA_HFX, SDA_QFX, HFX_BOTH, QFX_BOTH, QNORM, fasdas & !fasdas
127 !----------------------------------------------------------------
129 !----------------------------------------------------------------
130 !----------------------------------------------------------------
131 ! --- atmospheric (WRF generic) variables
132 !-- DT time step (seconds)
133 !-- DZ8W thickness of layers (m)
134 !-- T3D temperature (K)
135 !-- QV3D 3D water vapor mixing ratio (Kg/Kg)
136 !-- P3D 3D pressure (Pa)
137 !-- FLHC exchange coefficient for heat (m/s)
138 !-- FLQC exchange coefficient for moisture (m/s)
139 !-- PSFC surface pressure (Pa)
140 !-- XLAND land mask (1 for land, 2 for water)
141 !-- QGH saturated mixing ratio at 2 meter
142 !-- GSW downward short wave flux at ground surface (W/m^2)
143 !-- GLW downward long wave flux at ground surface (W/m^2)
144 !-- History variables
145 !-- CANWAT canopy moisture content (mm)
146 !-- TSK surface temperature (K)
147 !-- TSLB soil temp (k)
148 !-- SMOIS total soil moisture content (volumetric fraction)
149 !-- SH2O unfrozen soil moisture content (volumetric fraction)
150 ! note: frozen soil moisture (i.e., soil ice) = SMOIS - SH2O
151 !-- SNOWH actual snow depth (m)
152 !-- SNOW liquid water-equivalent snow depth (m)
153 !-- ALBEDO time-varying surface albedo including snow effect (unitless fraction)
154 !-- ALBBCK background surface albedo (unitless fraction)
155 !-- CHS surface exchange coefficient for heat and moisture (m s-1);
156 !-- CHS2 2m surface exchange coefficient for heat (m s-1);
157 !-- CQS2 2m surface exchange coefficient for moisture (m s-1);
159 !-- num_soil_layers the number of soil layers
160 !-- ZS depths of centers of soil layers (m)
161 !-- DZS thicknesses of soil layers (m)
162 !-- SLDPTH thickness of each soil layer (m, same as DZS)
163 !-- TMN soil temperature at lower boundary (K)
164 !-- SMCWLT wilting point (volumetric)
165 !-- SMCDRY dry soil moisture threshold where direct evap from
166 ! top soil layer ends (volumetric)
167 !-- SMCREF soil moisture threshold below which transpiration begins to
168 ! stress (volumetric)
169 !-- SMCMAX porosity, i.e. saturated value of soil moisture (volumetric)
170 !-- NROOT number of root layers, a function of veg type, determined
171 ! in subroutine redprm.
172 !-- SMSTAV Soil moisture availability for evapotranspiration (
173 ! fraction between SMCWLT and SMCMXA)
174 !-- SMSTOT Total soil moisture content frozen+unfrozen) in the soil column (mm)
176 !-- SNOWC fraction snow coverage (0-1.0)
177 ! --- vegetation variables
178 !-- SNOALB upper bound on maximum albedo over deep snow
179 !-- SHDMIN minimum areal fractional coverage of annual green vegetation
180 !-- SHDMAX maximum areal fractional coverage of annual green vegetation
181 !-- XLAI leaf area index (dimensionless)
182 !-- XLAI2 leaf area index (same as XLAI) passed to output (dimensionless)
183 !-- Z0BRD Background fixed roughness length (M)
184 !-- Z0 Background vroughness length (M) as function
185 !-- ZNT Time varying roughness length (M) as function
186 !-- ALBD(IVGTPK,ISN) background albedo reading from a table
188 !-- HFX upward heat flux at the surface (W/m^2)
189 !-- QFX upward moisture flux at the surface (kg/m^2/s)
190 !-- LH upward moisture flux at the surface (W m-2)
191 !-- GRDFLX(I,J) ground heat flux (W m-2)
192 !-- FDOWN radiation forcing at the surface (W m-2) = SOLDN*(1-alb)+LWDN
193 !----------------------------------------------------------------------------
194 !-- EC canopy water evaporation ((W m-2)
195 !-- EDIR direct soil evaporation (W m-2)
196 !-- ET plant transpiration from a particular root layer (W m-2)
197 !-- ETT total plant transpiration (W m-2)
198 !-- ESNOW sublimation from (or deposition to if <0) snowpack (W m-2)
199 !-- DRIP through-fall of precip and/or dew in excess of canopy
200 ! water-holding capacity (m)
201 !-- DEW dewfall (or frostfall for t<273.15) (M)
202 !-- SMAV Soil Moisture Availability for each layer, as a fraction
203 ! between SMCWLT and SMCMAX (dimensionless fraction)
204 ! ----------------------------------------------------------------------
205 !-- BETA ratio of actual/potential evap (dimensionless)
206 !-- ETP potential evaporation (W m-2)
207 ! ----------------------------------------------------------------------
208 !-- FLX1 precip-snow sfc (W m-2)
209 !-- FLX2 freezing rain latent heat flux (W m-2)
210 !-- FLX3 phase-change heat flux from snowmelt (W m-2)
211 ! ----------------------------------------------------------------------
212 !-- ACSNOM snow melt (mm) (water equivalent)
213 !-- ACSNOW accumulated snow fall (mm) (water equivalent)
214 !-- SNOPCX snow phase change heat flux (W/m^2)
215 !-- POTEVP accumulated potential evaporation (m)
216 !-- RIB Documentation needed!!!
217 ! ----------------------------------------------------------------------
218 !-- RUNOFF1 surface runoff (m s-1), not infiltrating the surface
219 !-- RUNOFF2 subsurface runoff (m s-1), drainage out bottom of last
220 ! soil layer (baseflow)
221 ! important note: here RUNOFF2 is actually the sum of RUNOFF2 and RUNOFF3
222 !-- RUNOFF3 numerical trunctation in excess of porosity (smcmax)
223 ! for a given soil layer at the end of a time step (m s-1).
224 !SFCRUNOFF Surface Runoff (mm)
225 !UDRUNOFF Total Underground Runoff (mm), which is the sum of RUNOFF2 and RUNOFF3
226 ! ----------------------------------------------------------------------
227 !-- RC canopy resistance (s m-1)
228 !-- RC2 canopy resistance (same as RC) passed to output
229 !-- PC plant coefficient (unitless fraction, 0-1) where PC*ETP = actual transp
230 !-- RSMIN minimum canopy resistance (s m-1)
231 !-- RCS incoming solar rc factor (dimensionless)
232 !-- RCT air temperature rc factor (dimensionless)
233 !-- RCQ atmos vapor pressure deficit rc factor (dimensionless)
234 !-- RCSOIL soil moisture rc factor (dimensionless)
236 !-- EMISS surface emissivity (between 0 and 1)
237 !-- EMBCK Background surface emissivity (between 0 and 1)
240 ! (R_d/R_v) (dimensionless)
241 !-- ids start index for i in domain
242 !-- ide end index for i in domain
243 !-- jds start index for j in domain
244 !-- jde end index for j in domain
245 !-- kds start index for k in domain
246 !-- kde end index for k in domain
247 !-- ims start index for i in memory
248 !-- ime end index for i in memory
249 !-- jms start index for j in memory
250 !-- jme end index for j in memory
251 !-- kms start index for k in memory
252 !-- kme end index for k in memory
253 !-- its start index for i in tile
254 !-- ite end index for i in tile
255 !-- jts start index for j in tile
256 !-- jte end index for j in tile
257 !-- kts start index for k in tile
258 !-- kte end index for k in tile
260 !-- SR fraction of frozen precip (0.0 to 1.0)
261 !----------------------------------------------------------------
265 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
266 ims,ime, jms,jme, kms,kme, &
267 its,ite, jts,jte, kts,kte
269 INTEGER, INTENT(IN ) :: sf_urban_physics !urban
270 INTEGER, INTENT(IN ) :: isurban
271 INTEGER, INTENT(IN ) :: isice
272 INTEGER, INTENT(IN ) :: julian, julyr !urban
274 !added by Wei Yu for routing
276 REAL, DIMENSION( ims:ime, jms:jme ) , &
277 INTENT(INOUT) :: sfcheadrt,INFXSRT,soldrain
282 ! new local vars for hydro
283 REAL :: etpnd1_hydro,sfcheadrt_hydro,infxsrt_hydro
285 REAL, DIMENSION( ims:ime, jms:jme ) , &
286 INTENT(IN ) :: TMN, &
294 SWDOWN, & !added 10 jan 2007
302 REAL, DIMENSION( ims:ime, jms:jme ) , &
303 INTENT(INOUT) :: ALBBCK, &
305 CHARACTER(LEN=*), INTENT(IN ) :: MMINLU
307 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
308 INTENT(IN ) :: QV3D, &
312 REAL, DIMENSION( ims:ime, jms:jme ) , &
313 INTENT(IN ) :: QGH, &
316 INTEGER, DIMENSION( ims:ime, jms:jme ) , &
317 INTENT(IN ) :: IVGTYP, &
320 INTEGER, INTENT(IN) :: num_soil_layers,ITIMESTEP
322 REAL, INTENT(IN ) :: DT,ROVCP
324 REAL, DIMENSION(1:num_soil_layers), INTENT(IN)::DZS
328 REAL, DIMENSION( ims:ime , 1:num_soil_layers, jms:jme ), &
329 INTENT(INOUT) :: SMOIS, & ! total soil moisture
330 SH2O, & ! new soil liquid
333 REAL, DIMENSION( ims:ime , 1:num_soil_layers, jms:jme ), &
334 INTENT(OUT) :: SMCREL
336 REAL, DIMENSION( ims:ime, jms:jme ) , &
337 INTENT(INOUT) :: TSK, & !was TGB (temperature)
363 REAL, DIMENSION( ims:ime, jms:jme ) , &
364 INTENT(OUT) :: NOAHRES
365 INTEGER, INTENT(IN) :: OPT_THCND
368 LOGICAL, INTENT(IN) :: UA_PHYS
369 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: FLX4_2D,FVB_2D,FBUR_2D,FGSN_2D
370 REAL :: FLX4,FVB,FBUR,FGSN
372 REAL, DIMENSION( ims:ime, jms:jme ) , &
373 INTENT(OUT) :: CHKLOWQ
374 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: LAI
375 REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: QZ0
376 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: RC2, XLAI2
378 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CMR_SFCDIF
379 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CHR_SFCDIF
380 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CMGR_SFCDIF
381 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CHGR_SFCDIF
382 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CMC_SFCDIF
383 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CHC_SFCDIF
384 ! Local variables (moved here from driver to make routine thread safe, 20031007 jm)
386 REAL, DIMENSION(1:num_soil_layers) :: ET
388 REAL, DIMENSION(1:num_soil_layers) :: SMAV
390 REAL :: BETA, ETP, SSOIL,EC, EDIR, ESNOW, ETT, &
391 FLX1,FLX2,FLX3, DRIP,DEW,FDOWN,RC,PC,RSMIN,XLAI, &
393 RCS,RCT,RCQ,RCSOIL,FFROZP
395 LOGICAL, INTENT(IN ) :: myj,frpcpn
397 ! DECLARATIONS - LOGICAL
398 ! ----------------------------------------------------------------------
399 LOGICAL, PARAMETER :: LOCAL=.false.
400 LOGICAL :: FRZGRA, SNOWNG
404 ! ----------------------------------------------------------------------
405 ! DECLARATIONS - INTEGER
406 ! ----------------------------------------------------------------------
407 INTEGER :: I,J, ICE,NSOIL,SLOPETYP,SOILTYP,VEGTYP
411 ! ----------------------------------------------------------------------
412 ! DECLARATIONS - REAL
413 ! ----------------------------------------------------------------------
415 REAL :: SHMIN,SHMAX,DQSDT2,LWDN,PRCP,PRCPRAIN, &
416 Q2SAT,Q2SATI,SFCPRS,SFCSPD,SFCTMP,SHDFAC,SNOALB1, &
417 SOLDN,TBOT,ZLVL, Q2K,ALBBRD, ALBEDOK, ETA, ETA_KINEMATIC, &
419 Z0K,RUNOFF1,RUNOFF2,RUNOFF3,SHEAT,SOLNET,E2SAT,SFCTSNO, &
420 ! mek, WRF testing, expanded diagnostics
421 SOLUP,LWUP,RNET,RES,Q1SFC,TAIRV,SATFLG
424 ! MEK JUL2007 for pot. evap.
430 REAL :: SNCOVR,SNEQV,SNOWHK,CMC, CHK,TH2
432 REAL :: SMCDRY,SMCMAX,SMCREF,SMCWLT,SNOMLT,SOILM,SOILW,Q1,T1
433 REAL :: SNOTIME1 ! LSTSNW1 INITIAL NUMBER OF TIMESTEPS SINCE LAST SNOWFALL
437 REAL :: COSZ, SOLARDIRECT
439 REAL, DIMENSION(1:num_soil_layers):: SLDPTH, STC,SMC,SWC
441 REAL, DIMENSION(1:num_soil_layers) :: ZSOIL, RTDIS
442 REAL, PARAMETER :: TRESH=.95E0, A2=17.67,A3=273.15,A4=29.65, &
443 T0=273.16E0, ELWV=2.50E6, A23M4=A2*(A3-A4)
445 REAL, PARAMETER :: ROW=1.E3,ELIW=XLF,ROWLIW=ROW*ELIW
447 ! ----------------------------------------------------------------------
448 ! DECLARATIONS START - urban
449 ! ----------------------------------------------------------------------
451 ! input variables surface_driver --> lsm
452 INTEGER, INTENT(IN) :: num_roof_layers
453 INTEGER, INTENT(IN) :: num_wall_layers
454 INTEGER, INTENT(IN) :: num_road_layers
455 REAL, OPTIONAL, DIMENSION(1:num_roof_layers), INTENT(IN) :: DZR
456 REAL, OPTIONAL, DIMENSION(1:num_wall_layers), INTENT(IN) :: DZB
457 REAL, OPTIONAL, DIMENSION(1:num_road_layers), INTENT(IN) :: DZG
458 REAL, OPTIONAL, INTENT(IN) :: DECLIN_URB
459 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: COSZ_URB2D
460 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: OMG_URB2D
461 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: XLAT_URB2D
462 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: U_PHY
463 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: V_PHY
464 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: TH_PHY
465 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: P_PHY
466 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: RHO
467 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: UST
469 LOGICAL, intent(in) :: rdlai2d
470 LOGICAL, intent(in) :: USEMONALB
472 ! input variables lsm --> urban
473 INTEGER :: UTYPE_URB ! urban type [urban=1, suburban=2, rural=3]
474 REAL :: TA_URB ! potential temp at 1st atmospheric level [K]
475 REAL :: QA_URB ! mixing ratio at 1st atmospheric level [kg/kg]
476 REAL :: UA_URB ! wind speed at 1st atmospheric level [m/s]
477 REAL :: U1_URB ! u at 1st atmospheric level [m/s]
478 REAL :: V1_URB ! v at 1st atmospheric level [m/s]
479 REAL :: SSG_URB ! downward total short wave radiation [W/m/m]
480 REAL :: LLG_URB ! downward long wave radiation [W/m/m]
481 REAL :: RAIN_URB ! precipitation [mm/h]
482 REAL :: RHOO_URB ! air density [kg/m^3]
483 REAL :: ZA_URB ! first atmospheric level [m]
484 REAL :: DELT_URB ! time step [s]
485 REAL :: SSGD_URB ! downward direct short wave radiation [W/m/m]
486 REAL :: SSGQ_URB ! downward diffuse short wave radiation [W/m/m]
487 REAL :: XLAT_URB ! latitude [deg]
488 REAL :: COSZ_URB ! cosz
489 REAL :: OMG_URB ! hour angle
490 REAL :: ZNT_URB ! roughness length [m]
501 REAL, DIMENSION(1:num_roof_layers) :: TRL_URB ! roof layer temp [K]
502 REAL, DIMENSION(1:num_wall_layers) :: TBL_URB ! wall layer temp [K]
503 REAL, DIMENSION(1:num_road_layers) :: TGL_URB ! road layer temp [K]
504 LOGICAL :: LSOLAR_URB
506 !===Yang,2014/10/08,hydrological variable for single layer UCM===
507 INTEGER :: jmonth, jday, tloc
508 INTEGER :: IRIOPTION, USOIL, DSOIL
518 REAL, DIMENSION(1:num_roof_layers) :: SMR_URB ! green roof layer moisture
519 REAL, DIMENSION(1:num_roof_layers) :: TGRL_URB ! green roof layer temp [K]
521 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: DRELR_URB2D
522 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: DRELB_URB2D
523 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: DRELG_URB2D
524 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: FLXHUMR_URB2D
525 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: FLXHUMB_URB2D
526 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: FLXHUMG_URB2D
527 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CMCR_URB2D
528 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TGR_URB2D
529 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_roof_layers, jms:jme ), INTENT(INOUT) :: TGRL_URB3D
530 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_roof_layers, jms:jme ), INTENT(INOUT) :: SMR_URB3D
533 ! state variable surface_driver <--> lsm <--> urban
534 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TR_URB2D
535 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TB_URB2D
536 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TG_URB2D
537 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TC_URB2D
538 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: QC_URB2D
539 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: UC_URB2D
540 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXR_URB2D
541 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXB_URB2D
542 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXG_URB2D
543 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXC_URB2D
544 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: SH_URB2D
545 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: LH_URB2D
546 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: G_URB2D
547 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: RN_URB2D
549 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TS_URB2D
551 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_roof_layers, jms:jme ), INTENT(INOUT) :: TRL_URB3D
552 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_wall_layers, jms:jme ), INTENT(INOUT) :: TBL_URB3D
553 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_road_layers, jms:jme ), INTENT(INOUT) :: TGL_URB3D
555 ! output variable lsm --> surface_driver
556 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: PSIM_URB2D
557 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: PSIH_URB2D
558 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: GZ1OZ0_URB2D
559 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: U10_URB2D
560 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: V10_URB2D
561 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: TH2_URB2D
562 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: Q2_URB2D
564 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: AKMS_URB2D
566 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: UST_URB2D
567 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: FRC_URB2D
568 INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: UTYPE_URB2D
571 ! output variables urban --> lsm
572 REAL :: TS_URB ! surface radiative temperature [K]
573 REAL :: QS_URB ! surface humidity [-]
574 REAL :: SH_URB ! sensible heat flux [W/m/m]
575 REAL :: LH_URB ! latent heat flux [W/m/m]
576 REAL :: LH_KINEMATIC_URB ! latent heat flux, kinetic [kg/m/m/s]
577 REAL :: SW_URB ! upward short wave radiation flux [W/m/m]
578 REAL :: ALB_URB ! time-varying albedo [fraction]
579 REAL :: LW_URB ! upward long wave radiation flux [W/m/m]
580 REAL :: G_URB ! heat flux into the ground [W/m/m]
581 REAL :: RN_URB ! net radiation [W/m/m]
582 REAL :: PSIM_URB ! shear f for momentum [-]
583 REAL :: PSIH_URB ! shear f for heat [-]
584 REAL :: GZ1OZ0_URB ! shear f for heat [-]
585 REAL :: U10_URB ! wind u component at 10 m [m/s]
586 REAL :: V10_URB ! wind v component at 10 m [m/s]
587 REAL :: TH2_URB ! potential temperature at 2 m [K]
588 REAL :: Q2_URB ! humidity at 2 m [-]
592 ! NUDAPT Parameters urban --> lam
597 REAL, DIMENSION(4) :: lf_urb
598 ! Variables for multi-layer UCM (Martilli et al. 2002)
599 REAL, OPTIONAL, INTENT(IN ) :: GMT
600 INTEGER, OPTIONAL, INTENT(IN ) :: JULDAY
601 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) ::XLAT, XLONG
602 INTEGER, INTENT(IN ) :: num_urban_ndm
603 INTEGER, INTENT(IN ) :: urban_map_zrd
604 INTEGER, INTENT(IN ) :: urban_map_zwd
605 INTEGER, INTENT(IN ) :: urban_map_gd
606 INTEGER, INTENT(IN ) :: urban_map_zd
607 INTEGER, INTENT(IN ) :: urban_map_zdf
608 INTEGER, INTENT(IN ) :: urban_map_bd
609 INTEGER, INTENT(IN ) :: urban_map_wd
610 INTEGER, INTENT(IN ) :: urban_map_gbd
611 INTEGER, INTENT(IN ) :: urban_map_fbd
612 INTEGER, INTENT(IN ) :: urban_map_zgrd
613 INTEGER, INTENT(IN ) :: NUM_URBAN_HI
614 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: tsk_rural_bep
615 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zrd, jms:jme ), INTENT(INOUT) :: trb_urb4d
616 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zwd, jms:jme ), INTENT(INOUT) :: tw1_urb4d
617 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zwd, jms:jme ), INTENT(INOUT) :: tw2_urb4d
618 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_gd , jms:jme ), INTENT(INOUT) :: tgb_urb4d
619 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_bd , jms:jme ), INTENT(INOUT) :: tlev_urb3d
620 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_bd , jms:jme ), INTENT(INOUT) :: qlev_urb3d
621 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_wd , jms:jme ), INTENT(INOUT) :: tw1lev_urb3d
622 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_wd , jms:jme ), INTENT(INOUT) :: tw2lev_urb3d
623 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_gbd, jms:jme ), INTENT(INOUT) :: tglev_urb3d
624 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_fbd, jms:jme ), INTENT(INOUT) :: tflev_urb3d
625 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: lf_ac_urb3d
626 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: sf_ac_urb3d
627 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: cm_ac_urb3d
628 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: sfvent_urb3d
629 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: lfvent_urb3d
630 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_wd , jms:jme ), INTENT(INOUT) :: sfwin1_urb3d
631 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_wd , jms:jme ), INTENT(INOUT) :: sfwin2_urb3d
632 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zd , jms:jme ), INTENT(INOUT) :: sfw1_urb3d
633 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zd , jms:jme ), INTENT(INOUT) :: sfw2_urb3d
634 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme ), INTENT(INOUT) :: sfr_urb3d
635 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_ndm, jms:jme ), INTENT(INOUT) :: sfg_urb3d
636 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: ep_pv_urb3d !GRZ
637 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme ), INTENT(INOUT) :: t_pv_urb3d !GRZ
638 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zgrd, jms:jme ),INTENT(INOUT) :: trv_urb4d !GRZ
639 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zgrd, jms:jme ),INTENT(INOUT) :: qr_urb4d !GRZ
640 REAL, OPTIONAL, DIMENSION( ims:ime,jms:jme ), INTENT(INOUT) :: qgr_urb3d !GRZ
641 REAL, OPTIONAL, DIMENSION( ims:ime,jms:jme ), INTENT(INOUT) :: tgr_urb3d !GRZ
642 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme ),INTENT(INOUT) :: drain_urb4d !GRZ
643 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: draingr_urb3d !GRZ
644 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme ),INTENT(INOUT) :: sfrv_urb3d !GRZ
645 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme ),INTENT(INOUT) :: lfrv_urb3d !GRZ
646 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme ),INTENT(INOUT) :: dgr_urb3d !GRZ
647 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_ndm, jms:jme ),INTENT(INOUT) :: dg_urb3d !GRZ
648 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme ),INTENT(INOUT) :: lfr_urb3d !GRZ
649 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_ndm, jms:jme ),INTENT(INOUT) :: lfg_urb3d !GRZ
651 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_hi, jms:jme ), INTENT(IN) :: hi_urb2d
652 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: lp_urb2d
653 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: lb_urb2d
654 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: hgt_urb2d
655 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: mh_urb2d
656 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: stdh_urb2d
657 REAL, OPTIONAL, DIMENSION( ims:ime, 4, jms:jme ), INTENT(IN) :: lf_urb2d
658 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_u_bep !Implicit momemtum component X-direction
659 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_v_bep !Implicit momemtum component Y-direction
660 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_t_bep !Implicit component pot. temperature
661 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_q_bep !Implicit momemtum component X-direction
662 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_e_bep !Implicit component TKE
663 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::b_u_bep !Explicit momentum component X-direction
664 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::b_v_bep !Explicit momentum component Y-direction
665 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::b_t_bep !Explicit component pot. temperature
666 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::b_q_bep !Implicit momemtum component Y-direction
667 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::b_e_bep !Explicit component TKE
668 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::vl_bep !Fraction air volume in grid cell
669 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::dlg_bep !Height above ground
670 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::sf_bep !Fraction air at the face of grid cell
671 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::dl_u_bep !Length scale
673 ! Local variables for multi-layer UCM (Martilli et al. 2002)
674 REAL, DIMENSION( its:ite, jts:jte ) :: HFX_RURAL,LH_RURAL,GRDFLX_RURAL ! ,RN_RURAL
675 REAL, DIMENSION( its:ite, jts:jte ) :: QFX_RURAL ! ,QSFC_RURAL,UMOM_RURAL,VMOM_RURAL
676 REAL, DIMENSION( its:ite, jts:jte ) :: ALB_RURAL,EMISS_RURAL,TSK_RURAL ! ,UST_RURAL
677 ! REAL, DIMENSION( ims:ime, jms:jme ) :: QSFC_URB
678 REAL, DIMENSION( its:ite, jts:jte ) :: HFX_URB,UMOM_URB,VMOM_URB
679 REAL, DIMENSION( its:ite, jts:jte ) :: QFX_URB
680 ! REAL, DIMENSION( ims:ime, jms:jme ) :: ALBEDO_URB,EMISS_URB,UMOM,VMOM,UST
681 REAL, DIMENSION(its:ite,jts:jte) ::EMISS_URB
682 REAL, DIMENSION(its:ite,jts:jte) :: RL_UP_URB
683 REAL, DIMENSION(its:ite,jts:jte) ::RS_ABS_URB
684 REAL, DIMENSION(its:ite,jts:jte) ::GRDFLX_URB
685 REAL :: SIGMA_SB,RL_UP_RURAL,RL_UP_TOT,RS_ABS_TOT,UMOM,VMOM
686 REAL :: CMR_URB, CHR_URB, CMC_URB, CHC_URB, CMGR_URB, CHGR_URB
687 REAL :: frc_urb,lb_urb
689 ! ----------------------------------------------------------------------
690 ! DECLARATIONS END - urban
691 ! ----------------------------------------------------------------------
693 REAL, PARAMETER :: CAPA=R_D/CP
694 REAL :: APELM,APES,SFCTH2,PSFC
695 real, intent(in) :: xice_threshold
696 character(len=80) :: message_text
700 REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, &
701 INTENT(INOUT) :: SDA_HFX, SDA_QFX, HFX_BOTH, QFX_BOTH, QNORM
702 INTEGER, INTENT(IN ) :: fasdas
704 REAL :: XSDA_HFX, XSDA_QFX, XQNORM
705 REAL :: HFX_PHY, QFX_PHY
709 REAL,OPTIONAL,INTENT(IN),DIMENSION( ims:ime, jms:jme ) :: IRR_CHAN
710 REAL :: IRRIGATION_CHANNEL
711 IRRIGATION_CHANNEL =0.0
712 HFX_PHY = 0.0 ! initialize
725 FDTW=DT/(XLV*RHOWATER)
734 NSOIL=num_soil_layers
742 IF(ITIMESTEP.EQ.1)THEN
744 !*** initialize soil conditions for IHOP 31 May case
745 ! IF((XLAND(I,J)-1.5) < 0.)THEN
746 ! if (I==108.and.j==85) then
754 !*** SET ZERO-VALUE FOR SOME OUTPUT DIAGNOSTIC ARRAYS
755 IF((XLAND(I,J)-1.5).GE.0.)THEN
756 ! check sea-ice point
758 IF( XICE(I,J).GE. XICE_THRESHOLD .and. IPRINT ) PRINT*, ' sea-ice at water point, I=',I,'J=',J
765 TSLB(I,NS,J)=273.16 !STEMP
769 IF ( XICE(I,J) .GE. XICE_THRESHOLD ) THEN
781 ENDIF ! end of initialization over ocean
783 !-----------------------------------------------------------------------
787 ! pressure in middle of lowest layer
788 SFCPRS=(P8W3D(I,KTS+1,j)+P8W3D(i,KTS,j))*0.5
789 ! convert from mixing ratio to specific humidity
790 Q2K=QV3D(i,1,j)/(1.0+QV3D(i,1,j))
791 ! initializing local variables
799 Q2SAT=QGH(I,J)/(1.0+QGH(I,J)) ! Q2SAT is sp humidity
800 ! add check on myj=.true.
801 ! IF((Q2K.GE.Q2SAT*TRESH).AND.Q2K.LT.QZ0(I,J))THEN
802 IF((myj).AND.(Q2K.GE.Q2SAT*TRESH).AND.Q2K.LT.QZ0(I,J))THEN
813 ! TH2=SFCTMP+(0.0097545*ZLVL)
814 ! calculate SFCTH2 via Exner function vs lapse-rate (above)
815 APES=(1.E5/PSFC)**CAPA
816 APELM=(1.E5/SFCPRS)**CAPA
822 ! SOLDN is total incoming solar
824 ! GSW is net downward solar
826 ! use mid-day albedo to determine net downward solar (no solar zenith angle correction)
827 SOLNET=SOLDN*(1.-ALBEDO(I,J))
829 IF(PRESENT(IRR_CHAN)) THEN
830 IF(IRR_CHAN(i,j).NE.0) THEN
831 IRRIGATION_CHANNEL=IRR_CHAN(i,j)/DT
833 IRRIGATION_CHANNEL=0.
838 SHDFAC=VEGFRA(I,J)/100.
841 SHMIN=SHDMIN(I,J)/100. !NEW
842 SHMAX=SHDMAX(I,J)/100. !NEW
843 ! convert snow water equivalent from mm to meter
844 SNEQV=SNOW(I,J)*0.001
845 ! snow depth in meters
849 ! if "SR" present, set frac of frozen precip ("FFROZP") = snow-ratio ("SR", range:0-1)
850 ! SR from e.g. Ferrier microphysics
851 ! otherwise define from 1st atmos level temperature
855 IF (SFCTMP <= 273.15) THEN
862 IF((XLAND(I,J)-1.5).GE.0.)THEN ! begining of land/sea if block
864 TSK_RURAL(I,J)=TSK(I,J)
865 HFX_RURAL(I,J)=HFX(I,J)
866 QFX_RURAL(I,J)=QFX(I,J)
867 LH_RURAL(I,J)=LH(I,J)
868 EMISS_RURAL(I,J)=EMISS(I,J)
869 GRDFLX_RURAL(I,J)=GRDFLX(I,J)
872 ! Land or sea-ice case
874 IF (XICE(I,J) >= XICE_THRESHOLD) THEN
877 ELSE IF ( VEGTYP == ISICE ) THEN
881 ! Neither sea ice or land ice.
884 DQSDT2=Q2SAT*A23M4/(SFCTMP-A4)**2
886 IF(SNOW(I,J).GT.0.0)THEN
887 ! snow on surface (use ice saturation properties)
889 E2SAT=611.2*EXP(6174.*(1./273.15 - 1./SFCTSNO))
890 Q2SATI=0.622*E2SAT/(SFCPRS-E2SAT)
891 Q2SATI=Q2SATI/(1.0+Q2SATI) ! spec. hum.
892 IF (T1 .GT. 273.14) THEN
893 ! warm ground temps, weight the saturation between ice and water according to SNOWC
894 Q2SAT=Q2SAT*(1.-SNOWC(I,J)) + Q2SATI*SNOWC(I,J)
895 DQSDT2=DQSDT2*(1.-SNOWC(I,J)) + Q2SATI*6174./(SFCTSNO**2)*SNOWC(I,J)
897 ! cold ground temps, use ice saturation only
899 DQSDT2=Q2SATI*6174./(SFCTSNO**2)
901 ! for snow cover fraction at 0 C, ground temp will not change, so DQSDT2 effectively zero
902 ! V3.8 add condition for SWDOWN to restrict condition to positive forcing (JD)
903 IF(T1 .GT. 273. .AND. SNOWC(I,J) .GT. 0. .AND. SWDOWN(I,J) .GT. 10.)DQSDT2=DQSDT2*(1.-SNOWC(I,J))
906 ! Land-ice or land points use the usual deep-soil temperature.
909 IF(ISURBAN.EQ.1) THEN
910 ! assumes these only need to be set for USGS land data
911 IF(VEGTYP.EQ.25) SHDFAC=0.0000
912 IF(VEGTYP.EQ.26) SHDFAC=0.0000
913 IF(VEGTYP.EQ.27) SHDFAC=0.0000
915 IF(SOILTYP.EQ.14.AND.XICE(I,J).EQ.0.)THEN
917 IF(IPRINT)PRINT*,' SOIL TYPE FOUND TO BE WATER AT A LAND-POINT'
918 IF(IPRINT)PRINT*,i,j,'RESET SOIL in surfce.F'
922 SNOALB1 = SNOALB(I,J)
923 ! converts canwat in mm to CMC in meters
924 CMC=CANWAT(I,J)/1000.
926 !-------------------------------------------
927 !*** convert snow depth from mm to meter
930 ! SNOALB=ALBMAX(I,J)*0.01
932 ! SNOALB=MAXALB(IVGTPK)*0.01
940 SNOTIME1 = SNOTIME(I,J)
942 !FEI: temporaray arrays above need to be changed later by using SI
945 SMC(NS)=SMOIS(I,NS,J)
946 STC(NS)=TSLB(I,NS,J) !STEMP
950 if ( (SNEQV.ne.0..AND.SNOWHK.eq.0.).or.(SNOWHK.le.SNEQV) )THEN
955 !Fei: urban. for urban surface, if calling UCM, redefine the natural surface in cities as
956 ! the "NATURAL" category in the VEGPARM.TBL
957 IF(SF_URBAN_PHYSICS == 1.OR. SF_URBAN_PHYSICS==2.OR.SF_URBAN_PHYSICS==3 ) THEN
958 IF( IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == LCZ_1 .or. IVGTYP(I,J) == LCZ_2 .or. &
959 IVGTYP(I,J) == LCZ_3 .or. IVGTYP(I,J) == LCZ_4 .or. IVGTYP(I,J) == LCZ_5 .or. &
960 IVGTYP(I,J) == LCZ_6 .or. IVGTYP(I,J) == LCZ_7 .or. IVGTYP(I,J) == LCZ_8 .or. &
961 IVGTYP(I,J) == LCZ_9 .or. IVGTYP(I,J) == LCZ_10 .or. IVGTYP(I,J) == LCZ_11 ) THEN
964 SHDFAC = SHDTBL(NATURAL)
967 EMISSI = 0.98 !for VEGTYP=5
968 LWDN = GLW(I,J) * EMISSI
969 SOLNET = SOLDN * (1.0 - ALBEDOK)
971 IF ( FRC_URB2D(I,J) < 0.99 ) THEN
972 if(sf_urban_physics.eq.1)then
973 T1= ( TSK(I,J) -FRC_URB2D(I,J) * TS_URB2D (I,J) )/ (1-FRC_URB2D(I,J))
974 elseif((sf_urban_physics.eq.2).OR.(sf_urban_physics.eq.3))then
975 T1=tsk_rural_bep(i,j)
982 IF( IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == LCZ_1 .or. IVGTYP(I,J) == LCZ_2 .or. &
983 IVGTYP(I,J) == LCZ_3 .or. IVGTYP(I,J) == LCZ_4 .or. IVGTYP(I,J) == LCZ_5 .or. &
984 IVGTYP(I,J) == LCZ_6 .or. IVGTYP(I,J) == LCZ_7 .or. IVGTYP(I,J) == LCZ_8 .or. &
985 IVGTYP(I,J) == LCZ_9 .or. IVGTYP(I,J) == LCZ_10 .or. IVGTYP(I,J) == LCZ_11 ) THEN
991 !===Yang, 2014/10/08, hydrological processes for urban vegetation in single layer UCM===
996 IF(SF_URBAN_PHYSICS == 1) THEN
998 tloc=mod(int(OMG/3.14159*180./15.+12.+0.5 ),24)
999 if (tloc.lt.0) tloc=tloc+24
1000 if (tloc==0) tloc=24
1001 CALL cal_mon_day(julian,julyr,jmonth,jday)
1002 IF( IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == LCZ_1 .or. IVGTYP(I,J) == LCZ_2 .or. &
1003 IVGTYP(I,J) == LCZ_3 .or. IVGTYP(I,J) == LCZ_4 .or. IVGTYP(I,J) == LCZ_5 .or. &
1004 IVGTYP(I,J) == LCZ_6 .or. IVGTYP(I,J) == LCZ_7 .or. IVGTYP(I,J) == LCZ_8 .or. &
1005 IVGTYP(I,J) == LCZ_9 .or. IVGTYP(I,J) == LCZ_10 .or. IVGTYP(I,J) == LCZ_11 ) THEN
1006 AOASIS = oasis ! urban oasis effect
1007 IF (IRIOPTION ==1) THEN
1008 IF (tloc==21 .or. tloc==22) THEN !irrigation on vegetaion in urban area, MAY-SEP, 9-10pm
1009 IF (jmonth==5 .or. jmonth==6 .or. jmonth==7 .or. jmonth==8 .or. jmonth==9) THEN
1010 ! IF (SMC(USOIL) .LT. SMCREF) SMC(USOIL)= SMCREF
1011 ! IF (SMC(DSOIL) .LT. SMCREF) SMC(DSOIL)= SMCREF
1012 IF (SMC(USOIL) .LT. SMCREF) SMC(USOIL)= REFSMC(ISLTYP(I,J))
1013 IF (SMC(DSOIL) .LT. SMCREF) SMC(DSOIL)= REFSMC(ISLTYP(I,J))
1020 IF(SF_URBAN_PHYSICS == 2 .or. SF_URBAN_PHYSICS == 3) THEN
1021 IF(AOASIS > 1.0) THEN
1022 FATAL_ERROR('Urban oasis option is for SF_URBAN_PHYSICS == 1 only')
1024 IF(IRIOPTION == 1) THEN
1025 FATAL_ERROR('Urban irrigation option is for SF_URBAN_PHYSICS == 1 only')
1032 print*, 'BEFORE SFLX, in Noahlsm_driver'
1033 print*, 'ICE', ICE, 'DT',DT, 'ZLVL',ZLVL, 'NSOIL', NSOIL, &
1034 'SLDPTH', SLDPTH, 'LOCAL',LOCAL, 'LUTYPE',&
1035 LUTYPE, 'SLTYPE',SLTYPE, 'LWDN',LWDN, 'SOLDN',SOLDN, &
1036 'SFCPRS',SFCPRS, 'PRCP',PRCP,'SFCTMP',SFCTMP,'Q2K',Q2K, &
1037 'TH2',TH2,'Q2SAT',Q2SAT,'DQSDT2',DQSDT2,'VEGTYP', VEGTYP,&
1038 'SOILTYP',SOILTYP, 'SLOPETYP',SLOPETYP, 'SHDFAC',SHDFAC,&
1039 'SHMIN',SHMIN, 'ALBBRD',ALBBRD,'SNOALB1',SNOALB1,'TBOT',&
1040 TBOT, 'Z0BRD',Z0BRD, 'Z0K',Z0K, 'CMC',CMC, 'T1',T1,'STC',&
1041 STC, 'SMC',SMC, 'SWC',SWC,'SNOWHK',SNOWHK,'SNEQV',SNEQV,&
1042 'ALBEDOK',ALBEDOK,'CHK',CHK,'ETA',ETA,'SHEAT',SHEAT, &
1043 'ETA_KINEMATIC',ETA_KINEMATIC, 'FDOWN',FDOWN,'EC',EC, &
1044 'EDIR',EDIR,'ET',ET,'ETT',ETT,'ESNOW',ESNOW,'DRIP',DRIP,&
1045 'DEW',DEW,'BETA',BETA,'ETP',ETP,'SSOIL',SSOIL,'FLX1',FLX1,&
1046 'FLX2',FLX2,'FLX3',FLX3,'SNOMLT',SNOMLT,'SNCOVR',SNCOVR,&
1047 'RUNOFF1',RUNOFF1,'RUNOFF2',RUNOFF2,'RUNOFF3',RUNOFF3, &
1048 'RC',RC, 'PC',PC,'RSMIN',RSMIN,'XLAI',XLAI,'RCS',RCS, &
1049 'RCT',RCT,'RCQ',RCQ,'RCSOIL',RCSOIL,'SOILW',SOILW, &
1050 'SOILM',SOILM,'Q1',Q1,'SMCWLT',SMCWLT,'SMCDRY',SMCDRY,&
1051 'SMCREF',SMCREF,'SMCMAX',SMCMAX,'NROOT',NROOT
1056 IF (SHDFAC > 0.0 .AND. LAI(I,J) <= 0.0) LAI(I,J) = 0.01
1060 IF ( ICE == 1 ) THEN
1071 ELSEIF (ICE == 0) THEN
1077 IF( fasdas == 1 ) THEN
1080 XSDA_HFX= SDA_HFX(I,J)*RHO(I,1,J)*CPM(I,J)*DZQ ! W/m^2
1081 ! TWG2015 Bugfix remove factor of 1000.0 for correct units
1082 XSDA_QFX= SDA_QFX(I,J)*RHO(I,1,J)*DZQ ! Kg/m2/s of water
1091 sfcheadrt_hydro = sfcheadrt(i,j)
1092 infxsrt_hydro = infxsrt(i,j)
1095 sfcheadrt_hydro = 0.
1098 CALL SFLX (I,J,FFROZP, ISURBAN, DT,ZLVL,NSOIL,SLDPTH, & !C
1100 LUTYPE, SLTYPE, & !CL
1101 LWDN,SOLDN,SOLNET,SFCPRS,PRCP,SFCTMP,Q2K,DUMMY, & !F
1102 DUMMY,DUMMY, DUMMY, & !F PRCPRAIN not used
1103 TH2,Q2SAT,DQSDT2, & !I
1104 VEGTYP,SOILTYP,SLOPETYP,SHDFAC,SHMIN,SHMAX, & !I
1105 ALBBRD, SNOALB1,TBOT, Z0BRD, Z0K, EMISSI, EMBRD, & !S
1106 CMC,T1,STC,SMC,SWC,SNOWHK,SNEQV,ALBEDOK,CHK,dummy,& !H
1107 ETA,SHEAT, ETA_KINEMATIC,FDOWN, & !O
1108 EC,EDIR,ET,ETT,ESNOW,DRIP,DEW, & !O
1109 BETA,ETP,SSOIL, & !O
1110 FLX1,FLX2,FLX3, & !O
1111 FLX4,FVB,FBUR,FGSN,UA_PHYS, & !UA
1113 RUNOFF1,RUNOFF2,RUNOFF3, & !O
1114 RC,PC,RSMIN,XLAI,RCS,RCT,RCQ,RCSOIL, & !O
1115 SOILW,SOILM,Q1,SMAV, & !D
1116 RDLAI2D,USEMONALB, &
1119 SMCWLT,SMCDRY,SMCREF,SMCMAX,NROOT, &
1121 sfcheadrt_hydro, & !I
1122 INFXSRT_hydro,ETPND1_hydro & !O
1123 ,OPT_THCND,AOASIS & !O
1124 ,XSDA_QFX, HFX_PHY, QFX_PHY, XQNORM, fasdas, HCPCT_FASDAS & ! fasdas
1125 ,IRRIGATION_CHANNEL)
1128 soldrain(i,j) = RUNOFF2*DT*1000.0
1129 sfcheadrt(i,j) = sfcheadrt_hydro
1130 infxsrt(i,j) = INFXSRT_hydro
1131 etpnd1 = etpnd1_hydro
1133 ELSEIF (ICE == -1) THEN
1136 ! Set values that the LSM is expected to update,
1137 ! but don't get updated for glacial points.
1139 SOILM = 0.0 !BSINGH(PNNL)- SOILM is undefined for this case, it is used for diagnostics so setting it to zero
1140 XLAI = 0.01 ! KWM Should this be Zero over land ice? Does this value matter?
1151 IF( fasdas == 1 ) THEN
1154 XSDA_HFX= SDA_HFX(I,J)*RHO(I,1,J)*CPM(I,J)*DZQ ! W/m^2
1155 XSDA_QFX= 0.0 ! Kg/m2/s of water
1163 CALL SFLX_GLACIAL(I,J,ISICE,FFROZP,DT,ZLVL,NSOIL,SLDPTH, & !C
1164 & LWDN,SOLNET,SFCPRS,PRCP,SFCTMP,Q2K, & !F
1165 & TH2,Q2SAT,DQSDT2, & !I
1166 & ALBBRD, SNOALB1,TBOT, Z0BRD, Z0K, EMISSI, EMBRD, & !S
1167 & T1,STC(1:NSOIL),SNOWHK,SNEQV,ALBEDOK,CHK, & !H
1168 & ETA,SHEAT,ETA_KINEMATIC,FDOWN, & !O
1171 & FLX1,FLX2,FLX3, & !O
1172 & SNOMLT,SNCOVR, & !O
1181 if (present(rc2) .and. present(xlai2)) then
1182 rc2(I,J) = RC ! for output
1189 print*, 'AFTER SFLX, in Noahlsm_driver'
1190 print*, 'ICE', ICE, 'DT',DT, 'ZLVL',ZLVL, 'NSOIL', NSOIL, &
1191 'SLDPTH', SLDPTH, 'LOCAL',LOCAL, 'LUTYPE',&
1192 LUTYPE, 'SLTYPE',SLTYPE, 'LWDN',LWDN, 'SOLDN',SOLDN, &
1193 'SFCPRS',SFCPRS, 'PRCP',PRCP,'SFCTMP',SFCTMP,'Q2K',Q2K, &
1194 'TH2',TH2,'Q2SAT',Q2SAT,'DQSDT2',DQSDT2,'VEGTYP', VEGTYP,&
1195 'SOILTYP',SOILTYP, 'SLOPETYP',SLOPETYP, 'SHDFAC',SHDFAC,&
1196 'SHDMIN',SHMIN, 'ALBBRD',ALBBRD,'SNOALB',SNOALB1,'TBOT',&
1197 TBOT, 'Z0BRD',Z0BRD, 'Z0K',Z0K, 'CMC',CMC, 'T1',T1,'STC',&
1198 STC, 'SMC',SMC, 'SWc',SWC,'SNOWHK',SNOWHK,'SNEQV',SNEQV,&
1199 'ALBEDOK',ALBEDOK,'CHK',CHK,'ETA',ETA,'SHEAT',SHEAT, &
1200 'ETA_KINEMATIC',ETA_KINEMATIC, 'FDOWN',FDOWN,'EC',EC, &
1201 'EDIR',EDIR,'ET',ET,'ETT',ETT,'ESNOW',ESNOW,'DRIP',DRIP,&
1202 'DEW',DEW,'BETA',BETA,'ETP',ETP,'SSOIL',SSOIL,'FLX1',FLX1,&
1203 'FLX2',FLX2,'FLX3',FLX3,'SNOMLT',SNOMLT,'SNCOVR',SNCOVR,&
1204 'RUNOFF1',RUNOFF1,'RUNOFF2',RUNOFF2,'RUNOFF3',RUNOFF3, &
1205 'RC',RC, 'PC',PC,'RSMIN',RSMIN,'XLAI',XLAI,'RCS',RCS, &
1206 'RCT',RCT,'RCQ',RCQ,'RCSOIL',RCSOIL,'SOILW',SOILW, &
1207 'SOILM',SOILM,'Q1',Q1,'SMCWLT',SMCWLT,'SMCDRY',SMCDRY,&
1208 'SMCREF',SMCREF,'SMCMAX',SMCMAX,'NROOT',NROOT
1212 !*** UPDATE STATE VARIABLES
1213 CANWAT(I,J)=CMC*1000.
1214 SNOW(I,J)=SNEQV*1000.
1215 ! SNOWH(I,J)=SNOWHK*1000.
1216 SNOWH(I,J)=SNOWHK ! SNOWHK in meters
1218 ALB_RURAL(I,J)=ALBEDOK
1222 EMISS_RURAL(I,J) = EMISSI
1223 ! Noah: activate time-varying roughness length (V3.3 Feb 2011)
1228 ! Update Skin Temperature
1229 IF( fasdas == 1 ) THEN
1230 XSDA_QFX= XSDA_QFX*ELWV*XQNORM
1232 !TWG2015 Bugfix to multiply Heat Capacity by Soil Depth for correct
1235 T1 = T1 + (XSDA_HFX-XSDA_QFX)*DT/(HCPCT_FASDAS*DZS(1))
1243 if (present(tsk_rural_bep)) then
1244 IF(SF_URBAN_PHYSICS == 2 .or. SF_URBAN_PHYSICS == 3) THEN
1245 TSK_RURAL_BEP(I,J)=T1
1249 HFX_RURAL(I,J)=SHEAT
1250 ! MEk Jul07 add potential evap accum
1251 POTEVP(I,J)=POTEVP(I,J)+ETP*FDTW
1252 QFX(I,J)=ETA_KINEMATIC
1253 QFX_RURAL(I,J)=ETA_KINEMATIC
1257 ! QFX(I,J) = QFX(I,J) + ETPND1
1258 ! ETA = ETA + ETPND1/2.501E6*dt
1259 !end added by Wei Yu
1266 GRDFLX_RURAL(I,J)=SSOIL
1269 SNOTIME(I,J) = SNOTIME1
1270 ! prevent diagnostic ground q (q1) from being greater than qsat(tsk)
1271 ! as happens over snow cover where the cqs2 value also becomes irrelevant
1272 ! by setting cqs2=chs in this situation the 2m q should become just qv(k=1)
1273 ! ww: comment out this change to avoid Q2 drop due to change of radiative flux
1274 ! IF (Q1 .GT. QSFC(I,J)) THEN
1275 ! CQS2(I,J) = CHS(I,J)
1278 ! Convert QSFC back to mixing ratio
1279 QSFC(I,J)= Q1/(1.0-Q1)
1281 ! QSFC_RURAL(I,J)= Q1/(1.0-Q1)
1282 ! Calculate momentum flux from rural surface for use with multi-layer UCM (Martilli et al. 2002)
1285 SMOIS(I,NS,J)=SMC(NS)
1286 TSLB(I,NS,J)=STC(NS) ! STEMP
1287 SH2O(I,NS,J)=SWC(NS)
1297 ! Residual of surface energy balance equation terms
1301 noahres(i,j) = ( solnet + lwdn ) - sheat + ssoil - eta &
1302 - ( emissi * STBOLT * (t1**4) ) - flx1 - flx2 - flx3 - flx4
1305 noahres(i,j) = ( solnet + lwdn ) - sheat + ssoil - eta &
1306 - ( emissi * STBOLT * (t1**4) ) - flx1 - flx2 - flx3
1310 IF (SF_URBAN_PHYSICS == 1 ) THEN ! Beginning of UCM CALL if block
1311 !--------------------------------------
1312 ! URBAN CANOPY MODEL START - urban
1313 !--------------------------------------
1314 ! Input variables lsm --> urban
1316 IF( IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == LCZ_1 .or. IVGTYP(I,J) == LCZ_2 .or. &
1317 IVGTYP(I,J) == LCZ_3 .or. IVGTYP(I,J) == LCZ_4 .or. IVGTYP(I,J) == LCZ_5 .or. &
1318 IVGTYP(I,J) == LCZ_6 .or. IVGTYP(I,J) == LCZ_7 .or. IVGTYP(I,J) == LCZ_8 .or. &
1319 IVGTYP(I,J) == LCZ_9 .or. IVGTYP(I,J) == LCZ_10 .or. IVGTYP(I,J) == LCZ_11 ) THEN
1325 UTYPE_URB = UTYPE_URB2D(I,J) !urban type (low, high or industrial)
1327 TA_URB = SFCTMP ! [K]
1328 QA_URB = Q2K ! [kg/kg]
1329 UA_URB = SQRT(U_PHY(I,1,J)**2.+V_PHY(I,1,J)**2.)
1330 U1_URB = U_PHY(I,1,J)
1331 V1_URB = V_PHY(I,1,J)
1332 IF(UA_URB < 1.) UA_URB=1. ! [m/s]
1333 SSG_URB = SOLDN ! [W/m/m]
1334 SSGD_URB = 0.8*SOLDN ! [W/m/m]
1335 SSGQ_URB = SSG_URB-SSGD_URB ! [W/m/m]
1336 LLG_URB = GLW(I,J) ! [W/m/m]
1337 RAIN_URB = RAINBL(I,J) / DT * 3600.0 ! [mm/hr]
1338 RHOO_URB = SFCPRS / (287.04 * SFCTMP * (1.0+ 0.61 * Q2K)) ![kg/m/m/m]
1340 DELT_URB = DT ! [sec]
1341 XLAT_URB = XLAT_URB2D(I,J) ! [deg]
1342 COSZ_URB = COSZ_URB2D(I,J) !
1343 OMG_URB = OMG_URB2D(I,J) !
1346 LSOLAR_URB = .FALSE.
1348 TR_URB = TR_URB2D(I,J)
1349 TB_URB = TB_URB2D(I,J)
1350 TG_URB = TG_URB2D(I,J)
1351 TC_URB = TC_URB2D(I,J)
1352 QC_URB = QC_URB2D(I,J)
1353 UC_URB = UC_URB2D(I,J)
1355 TGR_URB = TGR_URB2D(I,J)
1356 CMCR_URB = CMCR_URB2D(I,J)
1357 FLXHUMR_URB = FLXHUMR_URB2D(I,J)
1358 FLXHUMB_URB = FLXHUMB_URB2D(I,J)
1359 FLXHUMG_URB = FLXHUMG_URB2D(I,J)
1360 DRELR_URB = DRELR_URB2D(I,J)
1361 DRELB_URB = DRELB_URB2D(I,J)
1362 DRELG_URB = DRELG_URB2D(I,J)
1364 DO K = 1,num_roof_layers
1365 TRL_URB(K) = TRL_URB3D(I,K,J)
1366 SMR_URB(K) = SMR_URB3D(I,K,J)
1367 TGRL_URB(K)= TGRL_URB3D(I,K,J)
1369 DO K = 1,num_wall_layers
1370 TBL_URB(K) = TBL_URB3D(I,K,J)
1372 DO K = 1,num_road_layers
1373 TGL_URB(K) = TGL_URB3D(I,K,J)
1376 XXXR_URB = XXXR_URB2D(I,J)
1377 XXXB_URB = XXXB_URB2D(I,J)
1378 XXXG_URB = XXXG_URB2D(I,J)
1379 XXXC_URB = XXXC_URB2D(I,J)
1382 ! Limits to avoid dividing by small number
1383 if (CHS(I,J) < 1.0E-02) then
1386 if (CHS2(I,J) < 1.0E-02) then
1389 if (CQS2(I,J) < 1.0E-02) then
1394 CHS2_URB = CHS2(I,J)
1395 IF (PRESENT(CMR_SFCDIF)) THEN
1396 CMR_URB = CMR_SFCDIF(I,J)
1397 CHR_URB = CHR_SFCDIF(I,J)
1398 CMGR_URB = CMGR_SFCDIF(I,J)
1399 CHGR_URB = CHGR_SFCDIF(I,J)
1400 CMC_URB = CMC_SFCDIF(I,J)
1401 CHC_URB = CHC_SFCDIF(I,J)
1405 mh_urb = mh_urb2d(I,J)
1406 stdh_urb = stdh_urb2d(I,J)
1407 lp_urb = lp_urb2d(I,J)
1408 hgt_urb = hgt_urb2d(I,J)
1411 lf_urb(K)=lf_urb2d(I,K,J)
1413 frc_urb = frc_urb2d(I,J)
1414 lb_urb = lb_urb2d(I,J)
1416 if (I.eq.73.and.J.eq.125)THEN
1421 CALL cal_mon_day(julian,julyr,jmonth,jday)
1422 CALL urban(LSOLAR_URB, & ! I
1423 num_roof_layers,num_wall_layers,num_road_layers, & ! C
1425 UTYPE_URB,TA_URB,QA_URB,UA_URB,U1_URB,V1_URB,SSG_URB, & ! I
1426 SSGD_URB,SSGQ_URB,LLG_URB,RAIN_URB,RHOO_URB, & ! I
1427 ZA_URB,DECLIN_URB,COSZ_URB,OMG_URB, & ! I
1428 XLAT_URB,DELT_URB,ZNT_URB, & ! I
1429 CHS_URB, CHS2_URB, & ! I
1430 TR_URB, TB_URB, TG_URB, TC_URB, QC_URB,UC_URB, & ! H
1431 TRL_URB,TBL_URB,TGL_URB, & ! H
1432 XXXR_URB, XXXB_URB, XXXG_URB, XXXC_URB, & ! H
1433 TS_URB,QS_URB,SH_URB,LH_URB,LH_KINEMATIC_URB, & ! O
1434 SW_URB,ALB_URB,LW_URB,G_URB,RN_URB,PSIM_URB,PSIH_URB, & ! O
1436 CMR_URB, CHR_URB, CMC_URB, CHC_URB, &
1437 U10_URB, V10_URB, TH2_URB, Q2_URB, & ! O
1438 UST_URB,mh_urb, stdh_urb, lf_urb, lp_urb, & ! 0
1439 hgt_urb,frc_urb,lb_urb, check,CMCR_URB,TGR_URB, & ! H
1440 TGRL_URB,SMR_URB,CMGR_URB,CHGR_URB,jmonth, & ! H
1441 DRELR_URB,DRELB_URB, & ! H
1442 DRELG_URB,FLXHUMR_URB,FLXHUMB_URB,FLXHUMG_URB)
1447 print*, 'AFTER CALL URBAN'
1448 print*,'num_roof_layers',num_roof_layers, 'num_wall_layers', &
1450 'DZR',DZR,'DZB',DZB,'DZG',DZG,'UTYPE_URB',UTYPE_URB,'TA_URB', &
1452 'QA_URB',QA_URB,'UA_URB',UA_URB,'U1_URB',U1_URB,'V1_URB', &
1454 'SSG_URB',SSG_URB,'SSGD_URB',SSGD_URB,'SSGQ_URB',SSGQ_URB, &
1455 'LLG_URB',LLG_URB,'RAIN_URB',RAIN_URB,'RHOO_URB',RHOO_URB, &
1456 'ZA_URB',ZA_URB, 'DECLIN_URB',DECLIN_URB,'COSZ_URB',COSZ_URB,&
1457 'OMG_URB',OMG_URB,'XLAT_URB',XLAT_URB,'DELT_URB',DELT_URB, &
1458 'ZNT_URB',ZNT_URB,'TR_URB',TR_URB, 'TB_URB',TB_URB,'TG_URB',&
1459 TG_URB,'TC_URB',TC_URB,'QC_URB',QC_URB,'TRL_URB',TRL_URB, &
1460 'TBL_URB',TBL_URB,'TGL_URB',TGL_URB,'XXXR_URB',XXXR_URB, &
1461 'XXXB_URB',XXXB_URB,'XXXG_URB',XXXG_URB,'XXXC_URB',XXXC_URB,&
1462 'TS_URB',TS_URB,'QS_URB',QS_URB,'SH_URB',SH_URB,'LH_URB', &
1463 LH_URB, 'LH_KINEMATIC_URB',LH_KINEMATIC_URB,'SW_URB',SW_URB,&
1464 'ALB_URB',ALB_URB,'LW_URB',LW_URB,'G_URB',G_URB,'RN_URB', &
1465 RN_URB, 'PSIM_URB',PSIM_URB,'PSIH_URB',PSIH_URB, &
1466 'U10_URB',U10_URB,'V10_URB',V10_URB,'TH2_URB',TH2_URB, &
1467 'Q2_URB',Q2_URB,'CHS_URB',CHS_URB,'CHS2_URB',CHS2_URB
1471 TS_URB2D(I,J) = TS_URB
1473 ALBEDO(I,J) = FRC_URB2D(I,J)*ALB_URB+(1-FRC_URB2D(I,J))*ALBEDOK ![-]
1474 HFX(I,J) = FRC_URB2D(I,J)*SH_URB+(1-FRC_URB2D(I,J))*SHEAT ![W/m/m]
1475 QFX(I,J) = FRC_URB2D(I,J)*LH_KINEMATIC_URB &
1476 + (1-FRC_URB2D(I,J))*ETA_KINEMATIC ![kg/m/m/s]
1477 LH(I,J) = FRC_URB2D(I,J)*LH_URB+(1-FRC_URB2D(I,J))*ETA ![W/m/m]
1478 GRDFLX(I,J) = FRC_URB2D(I,J)*G_URB+(1-FRC_URB2D(I,J))*SSOIL ![W/m/m]
1479 TSK(I,J) = FRC_URB2D(I,J)*TS_URB+(1-FRC_URB2D(I,J))*T1 ![K]
1480 Q1 = FRC_URB2D(I,J)*QS_URB+(1-FRC_URB2D(I,J))*Q1 ![-]
1481 ! Convert QSFC back to mixing ratio
1482 QSFC(I,J)= Q1/(1.0-Q1)
1483 UST(I,J)= FRC_URB2D(I,J)*UST_URB+(1-FRC_URB2D(I,J))*UST(I,J) ![m/s]
1488 print*, ' FRC_URB2D', FRC_URB2D, &
1489 'ALB_URB',ALB_URB, 'ALBEDOK',ALBEDOK, &
1490 'ALBEDO(I,J)', ALBEDO(I,J), &
1491 'SH_URB',SH_URB,'SHEAT',SHEAT, 'HFX(I,J)',HFX(I,J), &
1492 'LH_KINEMATIC_URB',LH_KINEMATIC_URB,'ETA_KINEMATIC', &
1493 ETA_KINEMATIC, 'QFX(I,J)',QFX(I,J), &
1494 'LH_URB',LH_URB, 'ETA',ETA, 'LH(I,J)',LH(I,J), &
1495 'G_URB',G_URB,'SSOIL',SSOIL,'GRDFLX(I,J)', GRDFLX(I,J),&
1496 'TS_URB',TS_URB,'T1',T1,'TSK(I,J)',TSK(I,J), &
1497 'QS_URB',QS_URB,'Q1',Q1,'QSFC(I,J)',QSFC(I,J)
1503 ! Renew Urban State Varialbes
1505 TR_URB2D(I,J) = TR_URB
1506 TB_URB2D(I,J) = TB_URB
1507 TG_URB2D(I,J) = TG_URB
1508 TC_URB2D(I,J) = TC_URB
1509 QC_URB2D(I,J) = QC_URB
1510 UC_URB2D(I,J) = UC_URB
1512 TGR_URB2D(I,J) =TGR_URB
1513 CMCR_URB2D(I,J)=CMCR_URB
1514 FLXHUMR_URB2D(I,J)=FLXHUMR_URB
1515 FLXHUMB_URB2D(I,J)=FLXHUMB_URB
1516 FLXHUMG_URB2D(I,J)=FLXHUMG_URB
1517 DRELR_URB2D(I,J) = DRELR_URB
1518 DRELB_URB2D(I,J) = DRELB_URB
1519 DRELG_URB2D(I,J) = DRELG_URB
1521 DO K = 1,num_roof_layers
1522 TRL_URB3D(I,K,J) = TRL_URB(K)
1523 SMR_URB3D(I,K,J) = SMR_URB(K)
1524 TGRL_URB3D(I,K,J)= TGRL_URB(K)
1526 DO K = 1,num_wall_layers
1527 TBL_URB3D(I,K,J) = TBL_URB(K)
1529 DO K = 1,num_road_layers
1530 TGL_URB3D(I,K,J) = TGL_URB(K)
1532 XXXR_URB2D(I,J) = XXXR_URB
1533 XXXB_URB2D(I,J) = XXXB_URB
1534 XXXG_URB2D(I,J) = XXXG_URB
1535 XXXC_URB2D(I,J) = XXXC_URB
1537 SH_URB2D(I,J) = SH_URB
1538 LH_URB2D(I,J) = LH_URB
1539 G_URB2D(I,J) = G_URB
1540 RN_URB2D(I,J) = RN_URB
1541 PSIM_URB2D(I,J) = PSIM_URB
1542 PSIH_URB2D(I,J) = PSIH_URB
1543 GZ1OZ0_URB2D(I,J)= GZ1OZ0_URB
1544 U10_URB2D(I,J) = U10_URB
1545 V10_URB2D(I,J) = V10_URB
1546 TH2_URB2D(I,J) = TH2_URB
1547 Q2_URB2D(I,J) = Q2_URB
1548 UST_URB2D(I,J) = UST_URB
1549 AKMS_URB2D(I,J) = KARMAN * UST_URB2D(I,J)/(GZ1OZ0_URB2D(I,J)-PSIM_URB2D(I,J))
1550 IF (PRESENT(CMR_SFCDIF)) THEN
1551 CMR_SFCDIF(I,J) = CMR_URB
1552 CHR_SFCDIF(I,J) = CHR_URB
1553 CMGR_SFCDIF(I,J) = CMGR_URB
1554 CHGR_SFCDIF(I,J) = CHGR_URB
1555 CMC_SFCDIF(I,J) = CMC_URB
1556 CHC_SFCDIF(I,J) = CHC_URB
1560 ENDIF ! end of UCM CALL if block
1561 !--------------------------------------
1562 ! Urban Part End - urban
1563 !--------------------------------------
1567 SMSTOT(I,J)=SOILM*1000.
1569 SMCREL(I,NS,J)=SMAV(NS)
1573 ! Convert the water unit into mm
1574 SFCRUNOFF(I,J)=SFCRUNOFF(I,J)+RUNOFF1*DT*1000.0
1575 UDRUNOFF(I,J)=UDRUNOFF(I,J)+RUNOFF2*DT*1000.0
1576 ! snow defined when fraction of frozen precip (FFROZP) > 0.5,
1577 IF(FFROZP.GT.0.5)THEN
1578 ACSNOW(I,J)=ACSNOW(I,J)+PRCP*DT
1580 IF(SNOW(I,J).GT.0.)THEN
1581 ACSNOM(I,J)=ACSNOM(I,J)+SNOMLT*1000.
1582 ! accumulated snow-melt energy
1583 SNOPCX(I,J)=SNOPCX(I,J)-SNOMLT/FDTLIW
1586 ENDIF ! endif of land-sea test
1588 ENDDO ILOOP ! of I loop
1589 ENDDO JLOOP ! of J loop
1591 IF (SF_URBAN_PHYSICS == 2) THEN
1600 b_q_bep(i,kts:kte,j)=0.
1603 CALL BEP(frc_urb2d,utype_urb2d,itimestep,dz8w,dt,u_phy,v_phy, &
1604 th_phy,rho,p_phy,swdown,glw, &
1605 gmt,julday,xlong,xlat,declin_urb,cosz_urb2d,omg_urb2d, &
1606 num_urban_ndm, urban_map_zrd, urban_map_zwd, urban_map_gd, &
1607 urban_map_zd, urban_map_zdf, urban_map_bd, urban_map_wd, &
1608 urban_map_gbd, urban_map_fbd, num_urban_hi, &
1609 trb_urb4d,tw1_urb4d,tw2_urb4d,tgb_urb4d, &
1610 sfw1_urb3d,sfw2_urb3d,sfr_urb3d,sfg_urb3d, &
1611 lp_urb2d,hi_urb2d,lb_urb2d,hgt_urb2d, &
1612 a_u_bep,a_v_bep,a_t_bep, &
1613 a_e_bep,b_u_bep,b_v_bep, &
1614 b_t_bep,b_e_bep,b_q_bep,dlg_bep, &
1615 dl_u_bep,sf_bep,vl_bep, &
1616 rl_up_urb,rs_abs_urb,emiss_urb,grdflx_urb, &
1617 ids,ide, jds,jde, kds,kde, &
1618 ims,ime, jms,jme, kms,kme, &
1619 its,ite, jts,jte, kts,kte )
1624 IF (SF_URBAN_PHYSICS == 3) THEN
1633 b_q_bep(i,kts:kte,j)=0.
1638 CALL BEP_BEM(frc_urb2d,utype_urb2d,itimestep,dz8w,dt,u_phy,v_phy, &
1639 th_phy,rho,p_phy,swdown,glw, &
1640 gmt,julday,xlong,xlat,declin_urb,cosz_urb2d,omg_urb2d, &
1641 num_urban_ndm, urban_map_zrd, urban_map_zwd, urban_map_gd, &
1642 urban_map_zd, urban_map_zdf, urban_map_bd, urban_map_wd, &
1643 urban_map_gbd, urban_map_fbd, urban_map_zgrd, num_urban_hi, &
1644 trb_urb4d,tw1_urb4d,tw2_urb4d,tgb_urb4d, &
1645 tlev_urb3d,qlev_urb3d,tw1lev_urb3d,tw2lev_urb3d, &
1646 tglev_urb3d,tflev_urb3d,sf_ac_urb3d,lf_ac_urb3d, &
1647 cm_ac_urb3d,sfvent_urb3d,lfvent_urb3d, &
1648 sfwin1_urb3d,sfwin2_urb3d, &
1649 sfw1_urb3d,sfw2_urb3d,sfr_urb3d,sfg_urb3d, &
1650 ep_pv_urb3d,t_pv_urb3d, & !RMS
1651 trv_urb4d,qr_urb4d,qgr_urb3d,tgr_urb3d, & !RMS
1652 drain_urb4d,draingr_urb3d,sfrv_urb3d, & !RMS
1653 lfrv_urb3d,dgr_urb3d,dg_urb3d,lfr_urb3d,lfg_urb3d, & !RMS
1654 rainbl,swddir,swddif, & !RMS
1655 lp_urb2d,hi_urb2d,lb_urb2d,hgt_urb2d, &
1656 a_u_bep,a_v_bep,a_t_bep, &
1657 a_e_bep,b_u_bep,b_v_bep, &
1658 b_t_bep,b_e_bep,b_q_bep,dlg_bep, &
1659 dl_u_bep,sf_bep,vl_bep, &
1660 rl_up_urb,rs_abs_urb,emiss_urb,grdflx_urb,qv3d, &
1661 ids,ide, jds,jde, kds,kde, &
1662 ims,ime, jms,jme, kms,kme, &
1663 its,ite, jts,jte, kts,kte )
1667 if((sf_urban_physics.eq.2).OR.(sf_urban_physics.eq.3))then !Bep begin
1675 a_u_bep(i,k,j)=a_u_bep(i,k,j)*frc_urb2d(i,j)
1676 a_v_bep(i,k,j)=a_v_bep(i,k,j)*frc_urb2d(i,j)
1677 a_t_bep(i,k,j)=a_t_bep(i,k,j)*frc_urb2d(i,j)
1680 b_u_bep(i,k,j)=b_u_bep(i,k,j)*frc_urb2d(i,j)
1681 b_v_bep(i,k,j)=b_v_bep(i,k,j)*frc_urb2d(i,j)
1682 b_t_bep(i,k,j)=b_t_bep(i,k,j)*frc_urb2d(i,j)
1683 b_q_bep(i,k,j)=b_q_bep(i,k,j)*frc_urb2d(i,j)
1684 b_e_bep(i,k,j)=b_e_bep(i,k,j)*frc_urb2d(i,j)
1685 HFX_URB(I,J)=HFX_URB(I,J)+B_T_BEP(I,K,J)*RHO(I,K,J)*CP* &
1686 DZ8W(I,K,J)*VL_BEP(I,K,J)
1687 QFX_URB(I,J)=QFX_URB(I,J)+B_Q_BEP(I,K,J)* &
1688 DZ8W(I,K,J)*VL_BEP(I,K,J)
1689 UMOM_URB(I,J)=UMOM_URB(I,J)+ (A_U_BEP(I,K,J)*U_PHY(I,K,J)+ &
1690 B_U_BEP(I,K,J))*DZ8W(I,K,J)*VL_BEP(I,K,J)
1691 VMOM_URB(I,J)=VMOM_URB(I,J)+ (A_V_BEP(I,K,J)*V_PHY(I,K,J)+ &
1692 B_V_BEP(I,K,J))*DZ8W(I,K,J)*VL_BEP(I,K,J)
1693 vl_bep(i,k,j)=(1.-frc_urb2d(i,j))+vl_bep(i,k,j)*frc_urb2d(i,j)
1694 sf_bep(i,k,j)=(1.-frc_urb2d(i,j))+sf_bep(i,k,j)*frc_urb2d(i,j)
1696 a_u_bep(i,1,j)=(1.-frc_urb2d(i,j))*(-ust(I,J)*ust(I,J))/dz8w(i,1,j)/ &
1697 ((u_phy(i,1,j)**2+v_phy(i,1,j)**2.)**.5)+a_u_bep(i,1,j)
1698 a_v_bep(i,1,j)=(1.-frc_urb2d(i,j))*(-ust(I,J)*ust(I,J))/dz8w(i,1,j)/ &
1699 ((u_phy(i,1,j)**2+v_phy(i,1,j)**2.)**.5)+a_v_bep(i,1,j)
1700 b_t_bep(i,1,j)=(1.-frc_urb2d(i,j))*hfx_rural(i,j)/dz8w(i,1,j)/rho(i,1,j)/CP+ &
1702 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)
1703 umom=(1.-frc_urb2d(i,j))*ust(i,j)*ust(i,j)*u_phy(i,1,j)/ &
1704 ((u_phy(i,1,j)**2+v_phy(i,1,j)**2.)**.5)+umom_urb(i,j)
1705 vmom=(1.-frc_urb2d(i,j))*ust(i,j)*ust(i,j)*v_phy(i,1,j)/ &
1706 ((u_phy(i,1,j)**2+v_phy(i,1,j)**2.)**.5)+vmom_urb(i,j)
1709 ! compute upward longwave radiation from the rural part and total
1710 ! rl_up_rural=-emiss_rural(i,j)*sigma_sb*(tsk_rural(i,j)**4.)-(1.-emiss_rural(i,j))*glw(i,j)
1711 ! rl_up_tot=(1.-frc_urb2d(i,j))*rl_up_rural+frc_urb2d(i,j)*rl_up_urb(i,j)
1712 ! emiss(i,j)=(1.-frc_urb2d(i,j))*emiss_rural(i,j)+frc_urb2d(i,j)*emiss_urb(i,j)
1713 ! using the emissivity and the total longwave upward radiation estimate the averaged skin temperature
1714 IF (FRC_URB2D(I,J).GT.0.) THEN
1715 rl_up_rural=-emiss_rural(i,j)*sigma_sb*(tsk_rural(i,j)**4.)-(1.-emiss_rural(i,j))*glw(i,j)
1716 rl_up_tot=(1.-frc_urb2d(i,j))*rl_up_rural+frc_urb2d(i,j)*rl_up_urb(i,j)
1717 emiss(i,j)=(1.-frc_urb2d(i,j))*emiss_rural(i,j)+frc_urb2d(i,j)*emiss_urb(i,j)
1718 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
1719 tsk(i,j)=(max(0., (-1.*rl_up_tot-(1.-emiss(i,j))*glw(i,j) )/emiss(i,j)/sigma_sb))**.25
1720 rs_abs_tot=(1.-frc_urb2d(i,j))*swdown(i,j)*(1.-albedo(i,j))+frc_urb2d(i,j)*rs_abs_urb(i,j)
1721 if(swdown(i,j).gt.0.)then
1722 albedo(i,j)=1.-rs_abs_tot/swdown(i,j)
1724 albedo(i,j)=alb_rural(i,j)
1726 ! rename *_urb to sh_urb2d,lh_urb2d,g_urb2d,rn_urb2d
1727 grdflx(i,j)= (1.-frc_urb2d(i,j))*grdflx_rural(i,j)+frc_urb2d(i,j)*grdflx_urb(i,j)
1728 qfx(i,j)=(1.-frc_urb2d(i,j))*qfx_rural(i,j)+qfx_urb(i,j)
1729 ! lh(i,j)=(1.-frc_urb2d(i,j))*qfx_rural(i,j)*xlv
1730 lh(i,j)=qfx(i,j)*xlv
1731 HFX(I,J) = HFX_URB(I,J)+(1-FRC_URB2D(I,J))*HFX_RURAL(I,J) ![W/m/m]
1732 SH_URB2D(I,J) = HFX_URB(I,J)/FRC_URB2D(I,J)
1733 LH_URB2D(I,J) = qfx_urb(i,j)*xlv
1734 G_URB2D(I,J) = grdflx_urb(i,j)
1735 RN_URB2D(I,J) = rs_abs_urb(i,j)+emiss_urb(i,j)*glw(i,j)-rl_up_urb(i,j)
1736 ust(i,j)=(umom**2.+vmom**2.)**.25
1737 ! if(tsk(i,j).gt.350)write(*,*)'tsk too big!',i,j,tsk(i,j)
1738 ! 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
1739 ! print*,'ivgtyp,i,j,sigma_sb',ivgtyp(i,j),i,j,sigma_sb
1740 ! print*,'hfx,lh,qfx,grdflx,ts_urb2d',hfx(i,j),lh(i,j),qfx(i,j),grdflx(i,j),ts_urb2d(i,j)
1741 ! print*,'tsk,albedo,emiss',tsk(i,j),albedo(i,j),emiss(i,j)
1742 ! if(i.eq.56.and.j.eq.29)then
1743 ! print*,'ivgtyp, qfx, hfx',ivgtyp(i,j),hfx_rural(i,j),qfx_rural(i,j)
1744 ! print*,'emiss_rural,emiss_urb',emiss_rural(i,j),emiss_urb(i,j)
1745 ! print*,'rl_up_rural,rl_up_urb(i,j)',rl_up_rural,rl_up_urb(i,j)
1746 ! print*,'tsk_rural,ts_urb2d(i,j),tsk',tsk_rural(i,j),ts_urb2d(i,j),tsk(i,j)
1747 ! 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
1748 ! print*,'ivgtyp,hfx,hfx_urb,hfx_rural',hfx(i,j),hfx_urb(i,j),hfx_rural(i,j)
1749 ! print*,'lh,lh_rural',lh(i,j),lh_rural(i,j)
1750 ! print*,'qfx',qfx(i,j)
1751 ! print*,'ts_urb2d',ts_urb2d(i,j)
1752 ! print*,'ust',ust(i,j)
1753 ! print*,'swdown,glw',swdown(i,j),glw(i,j)
1766 !------------------------------------------------------
1768 !------------------------------------------------------
1770 !For MPAS, the below section of the module is moved to module_physics_lsm_noahinit.F in
1771 !./../core_physics to accomodate differences in the mpi calls between WRF and MPAS.I thought
1772 !that it would be cleaner to do this instead of adding a lot of #ifdef statements throughout
1773 !the initialization subroutine.
1775 #if defined(wrfmodel)
1776 SUBROUTINE LSMINIT(VEGFRA,SNOW,SNOWC,SNOWH,CANWAT,SMSTAV, &
1777 SMSTOT, SFCRUNOFF,UDRUNOFF,ACSNOW, &
1778 ACSNOM,IVGTYP,ISLTYP,TSLB,SMOIS,SH2O,ZS,DZS, &
1780 SNOALB, FNDSOILW, FNDSNOWH, RDMAXALB, &
1781 num_soil_layers, restart, &
1783 irr_rand_field,irr_ph,irr_freq, &
1784 ids,ide, jds,jde, kds,kde, &
1785 ims,ime, jms,jme, kms,kme, &
1786 its,ite, jts,jte, kts,kte &
1789 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
1790 ims,ime, jms,jme, kms,kme, &
1791 its,ite, jts,jte, kts,kte
1793 INTEGER, INTENT(IN) :: num_soil_layers
1795 LOGICAL , INTENT(IN) :: restart , allowed_to_read
1797 REAL, DIMENSION( num_soil_layers), INTENT(INOUT) :: ZS, DZS
1799 REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , &
1800 INTENT(INOUT) :: SMOIS, & !Total soil moisture
1801 SH2O, & !liquid soil moisture
1804 REAL, DIMENSION( ims:ime, jms:jme ) , &
1805 INTENT(INOUT) :: SNOW, &
1818 INTEGER, DIMENSION( ims:ime, jms:jme ) , &
1819 INTENT(IN) :: IVGTYP, &
1821 CHARACTER(LEN=*), INTENT(IN) :: MMINLU
1823 LOGICAL, INTENT(IN) :: FNDSOILW , &
1825 LOGICAL, INTENT(IN) :: RDMAXALB
1829 REAL :: BX, SMCMAX, PSISAT, FREE
1830 REAL, PARAMETER :: BLIM = 5.5, HLICE = 3.335E5, &
1831 GRAV = 9.81, T0 = 273.15
1833 CHARACTER(LEN=80) :: err_message
1834 INTEGER,DIMENSION(ims:ime, jms:jme ),INTENT(INOUT):: irr_rand_field
1835 INTEGER , DIMENSION(jds:jde) :: my_seeds
1836 INTEGER :: irr_ph,irr_freq
1837 REAL,DIMENSION(ims:ime, jms:jme ) :: rand_tmp
1838 character*256 :: MMINSL
1842 ! initialize three Noah LSM related tables
1843 IF ( allowed_to_read ) THEN
1844 CALL wrf_message( 'INITIALIZE THREE Noah LSM RELATED TABLES' )
1845 CALL SOIL_VEG_GEN_PARM( MMINLU, MMINSL )
1849 ! 20130219 - No longer need these - see module_data_gocart_dust
1850 !#if ( WRF_CHEM == 1 )
1852 !! need this parameter for dust parameterization in wrf/chem
1855 ! porosity(i)=maxsmc(i)
1856 ! drypoint(i)=drysmc(i)
1861 IF(.not.restart)THEN
1864 IF (irr_ph.NE.0)THEN
1867 my_seeds(j) =sqrt(ide*(real(j-1)**2))+sqrt(real(jde*i))
1868 ! PRINT*,'myseed', my_seeds(j),j,jts,jds
1870 CALL RANDOM_SEED ( PUT = my_seeds )
1871 CALL RANDOM_NUMBER ( rand_tmp(i,:) )
1872 CALL RANDOM_SEED ( GET = my_seeds )
1873 CALL RANDOM_NUMBER ( rand_tmp(i,:) )
1874 irr_rand_field(i,:)=int(modulo(rand_tmp(i,:)*100,real(irr_freq)))
1885 IF ( ISLTYP( i,j ) .LT. 1 ) THEN
1887 WRITE(err_message,*)"module_sf_noahdrv.F: lsminit: out of range ISLTYP ",i,j,ISLTYP( i,j )
1888 CALL wrf_message(err_message)
1890 IF(.not.RDMAXALB) THEN
1891 SNOALB(i,j)=MAXALB(IVGTYP(i,j))*0.01
1895 IF ( errflag .EQ. 1 ) THEN
1897 CALL wrf_message( "WARNING: message only; was fatal. module_sf_noahdrv.F: lsminit: out of range value "// &
1898 "of ISLTYP. Is this field in the input?" )
1900 CALL wrf_error_fatal( "module_sf_noahdrv.F: lsminit: out of range value "// &
1901 "of ISLTYP. Is this field in the input?" )
1905 ! initialize soil liquid water content SH2O
1907 ! IF(.NOT.FNDSOILW) THEN
1909 ! If no SWC, do the following
1910 ! PRINT *,'SOIL WATER NOT FOUND - VALUE SET IN LSMINIT'
1913 BX = BB(ISLTYP(I,J))
1914 SMCMAX = MAXSMC(ISLTYP(I,J))
1915 PSISAT = SATPSI(ISLTYP(I,J))
1916 if ((bx > 0.0).and.(smcmax > 0.0).and.(psisat > 0.0)) then
1917 DO NS=1, num_soil_layers
1918 ! ----------------------------------------------------------------------
1919 !SH2O <= SMOIS for T < 273.149K (-0.001C)
1920 IF (TSLB(I,NS,J) < 273.149) THEN
1921 ! ----------------------------------------------------------------------
1922 ! first guess following explicit solution for Flerchinger Eqn from Koren
1923 ! et al, JGR, 1999, Eqn 17 (KCOUNT=0 in FUNCTION FRH2O).
1924 ! ISLTPK is soil type
1925 BX = BB(ISLTYP(I,J))
1926 SMCMAX = MAXSMC(ISLTYP(I,J))
1927 PSISAT = SATPSI(ISLTYP(I,J))
1928 IF ( BX > BLIM ) BX = BLIM
1929 FK=(( (HLICE/(GRAV*(-PSISAT))) * &
1930 ((TSLB(I,NS,J)-T0)/TSLB(I,NS,J)) )**(-1/BX) )*SMCMAX
1931 IF (FK < 0.02) FK = 0.02
1932 SH2O(I,NS,J) = MIN( FK, SMOIS(I,NS,J) )
1933 ! ----------------------------------------------------------------------
1934 ! now use iterative solution for liquid soil water content using
1935 ! FUNCTION FRH2O with the initial guess for SH2O from above explicit
1937 CALL FRH2O (FREE,TSLB(I,NS,J),SMOIS(I,NS,J),SH2O(I,NS,J), &
1940 ELSE ! of IF (TSLB(I,NS,J)
1941 ! ----------------------------------------------------------------------
1942 ! SH2O = SMOIS ( for T => 273.149K (-0.001C)
1943 SH2O(I,NS,J)=SMOIS(I,NS,J)
1944 ! ----------------------------------------------------------------------
1945 ENDIF ! of IF (TSLB(I,NS,J)
1946 END DO ! of DO NS=1, num_soil_layers
1947 else ! of if ((bx > 0.0)
1948 DO NS=1, num_soil_layers
1949 SH2O(I,NS,J)=SMOIS(I,NS,J)
1951 endif ! of if ((bx > 0.0)
1952 ENDDO ! DO I = its,itf
1953 ENDDO ! DO J = jts,jtf
1954 ! ENDIF ! of IF(.NOT.FNDSOILW)THEN
1956 ! initialize physical snow height SNOWH
1958 IF(.NOT.FNDSNOWH)THEN
1959 ! If no SNOWH do the following
1960 CALL wrf_message( 'SNOW HEIGHT NOT FOUND - VALUE DEFINED IN LSMINIT' )
1963 SNOWH(I,J)=SNOW(I,J)*0.005 ! SNOW in mm and SNOWH in m
1968 ! initialize canopy water to ZERO
1971 ! print*,'Note that canopy water content (CANWAT) is set to ZERO in LSMINIT'
1980 !------------------------------------------------------------------------------
1981 END SUBROUTINE lsminit
1982 !------------------------------------------------------------------------------
1986 !-----------------------------------------------------------------
1987 SUBROUTINE SOIL_VEG_GEN_PARM( MMINLU, MMINSL)
1988 !-----------------------------------------------------------------
1990 USE module_wrf_error
1993 CHARACTER(LEN=*), INTENT(IN) :: MMINLU, MMINSL
1994 integer :: LUMATCH, IINDEX, LC, NUM_SLOPE
1996 INTEGER , PARAMETER :: OPEN_OK = 0
1998 character*128 :: mess , message
1999 character*256 :: a_string
2000 logical, external :: wrf_dm_on_monitor
2001 integer , parameter :: loop_max = 10
2002 integer :: loop_count
2005 !-----SPECIFY VEGETATION RELATED CHARACTERISTICS :
2006 ! ALBBCK: SFC albedo (in percentage)
2007 ! Z0: Roughness length (m)
2008 ! SHDFAC: Green vegetation fraction (in percentage)
2009 ! Note: The ALBEDO, Z0, and SHDFAC values read from the following table
2010 ! ALBEDO, amd Z0 are specified in LAND-USE TABLE; and SHDFAC is
2011 ! the monthly green vegetation data
2012 ! CMXTBL: MAX CNPY Capacity (m)
2013 ! NROTBL: Rooting depth (layer)
2014 ! RSMIN: Mimimum stomatal resistance (s m-1)
2015 ! RSMAX: Max. stomatal resistance (s m-1)
2016 ! RGL: Parameters used in radiation stress function
2017 ! HS: Parameter used in vapor pressure deficit functio
2018 ! TOPT: Optimum transpiration air temperature. (K)
2019 ! CMCMAX: Maximum canopy water capacity
2020 ! CFACTR: Parameter used in the canopy inteception calculati
2021 ! SNUP: Threshold snow depth (in water equivalent m) that
2022 ! implies 100% snow cover
2023 ! LAI: Leaf area index (dimensionless)
2024 ! MAXALB: Upper bound on maximum albedo over deep snow
2026 !-----READ IN VEGETAION PROPERTIES FROM VEGPARM.TBL
2029 IF ( wrf_dm_on_monitor() ) THEN
2031 OPEN(19, FILE='VEGPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
2032 IF(ierr .NE. OPEN_OK ) THEN
2033 WRITE(message,FMT='(A)') &
2034 'module_sf_noahlsm.F: soil_veg_gen_parm: failure opening VEGPARM.TBL'
2035 CALL wrf_error_fatal ( message )
2042 READ (19,FMT='(A)',END=2002) a_string
2043 FIND_LUTYPE : DO WHILE (LUMATCH == 0)
2044 READ (19,*,END=2002)LUTYPE
2045 READ (19,*)LUCATS,IINDEX
2047 IF(LUTYPE.EQ.MMINLU)THEN
2048 WRITE( mess , * ) 'LANDUSE TYPE = ' // TRIM ( LUTYPE ) // ' FOUND', LUCATS,' CATEGORIES'
2049 CALL wrf_message( mess )
2052 loop_count = loop_count+1
2053 call wrf_message ( "Skipping over LUTYPE = " // TRIM ( LUTYPE ) )
2054 FIND_VEGETATION_PARAMETER_FLAG : DO
2055 READ (19,FMT='(A)', END=2002) a_string
2056 IF ( a_string(1:21) .EQ. 'Vegetation Parameters' ) THEN
2057 EXIT FIND_VEGETATION_PARAMETER_FLAG
2058 ELSE IF ( loop_count .GE. loop_max ) THEN
2059 CALL wrf_error_fatal ( 'Too many loops in VEGPARM.TBL')
2061 ENDDO FIND_VEGETATION_PARAMETER_FLAG
2064 ! prevent possible array overwrite, Bill Bovermann, IBM, May 6, 2008
2065 IF ( SIZE(SHDTBL) < LUCATS .OR. &
2066 SIZE(NROTBL) < LUCATS .OR. &
2067 SIZE(RSTBL) < LUCATS .OR. &
2068 SIZE(RGLTBL) < LUCATS .OR. &
2069 SIZE(HSTBL) < LUCATS .OR. &
2070 SIZE(SNUPTBL) < LUCATS .OR. &
2071 SIZE(MAXALB) < LUCATS .OR. &
2072 SIZE(LAIMINTBL) < LUCATS .OR. &
2073 SIZE(LAIMAXTBL) < LUCATS .OR. &
2074 SIZE(Z0MINTBL) < LUCATS .OR. &
2075 SIZE(Z0MAXTBL) < LUCATS .OR. &
2076 SIZE(ALBEDOMINTBL) < LUCATS .OR. &
2077 SIZE(ALBEDOMAXTBL) < LUCATS .OR. &
2078 SIZE(ZTOPVTBL) < LUCATS .OR. &
2079 SIZE(ZBOTVTBL) < LUCATS .OR. &
2080 SIZE(EMISSMINTBL ) < LUCATS .OR. &
2081 SIZE(EMISSMAXTBL ) < LUCATS ) THEN
2082 CALL wrf_error_fatal('Table sizes too small for value of LUCATS in module_sf_noahdrv.F')
2085 IF(LUTYPE.EQ.MMINLU)THEN
2087 READ (19,*)IINDEX,SHDTBL(LC), &
2088 NROTBL(LC),RSTBL(LC),RGLTBL(LC),HSTBL(LC), &
2089 SNUPTBL(LC),MAXALB(LC), LAIMINTBL(LC), &
2090 LAIMAXTBL(LC),EMISSMINTBL(LC), &
2091 EMISSMAXTBL(LC), ALBEDOMINTBL(LC), &
2092 ALBEDOMAXTBL(LC), Z0MINTBL(LC), Z0MAXTBL(LC),&
2093 ZTOPVTBL(LC), ZBOTVTBL(LC)
2097 READ (19,*)TOPT_DATA
2099 READ (19,*)CMCMAX_DATA
2101 READ (19,*)CFACTR_DATA
2103 READ (19,*)RSMAX_DATA
2110 READ (19,FMT='(A)') a_string
2111 IF ( a_string(1:21) .EQ. 'Vegetation Parameters' ) THEN
2112 CALL wrf_message ("Expected low and high density residential, and high density industrial information in VEGPARM.TBL")
2113 CALL wrf_error_fatal ("This could be caused by using an older VEGPARM.TBL file with a newer WRF source code.")
2140 IF (LUMATCH == 0) then
2141 CALL wrf_error_fatal ("Land Use Dataset '"//MMINLU//"' not found in VEGPARM.TBL.")
2145 CALL wrf_dm_bcast_string ( LUTYPE , 4 )
2146 CALL wrf_dm_bcast_integer ( LUCATS , 1 )
2147 CALL wrf_dm_bcast_integer ( IINDEX , 1 )
2148 CALL wrf_dm_bcast_integer ( LUMATCH , 1 )
2149 CALL wrf_dm_bcast_real ( SHDTBL , NLUS )
2150 CALL wrf_dm_bcast_real ( NROTBL , NLUS )
2151 CALL wrf_dm_bcast_real ( RSTBL , NLUS )
2152 CALL wrf_dm_bcast_real ( RGLTBL , NLUS )
2153 CALL wrf_dm_bcast_real ( HSTBL , NLUS )
2154 CALL wrf_dm_bcast_real ( SNUPTBL , NLUS )
2155 CALL wrf_dm_bcast_real ( LAIMINTBL , NLUS )
2156 CALL wrf_dm_bcast_real ( LAIMAXTBL , NLUS )
2157 CALL wrf_dm_bcast_real ( Z0MINTBL , NLUS )
2158 CALL wrf_dm_bcast_real ( Z0MAXTBL , NLUS )
2159 CALL wrf_dm_bcast_real ( EMISSMINTBL , NLUS )
2160 CALL wrf_dm_bcast_real ( EMISSMAXTBL , NLUS )
2161 CALL wrf_dm_bcast_real ( ALBEDOMINTBL , NLUS )
2162 CALL wrf_dm_bcast_real ( ALBEDOMAXTBL , NLUS )
2163 CALL wrf_dm_bcast_real ( ZTOPVTBL , NLUS )
2164 CALL wrf_dm_bcast_real ( ZBOTVTBL , NLUS )
2165 CALL wrf_dm_bcast_real ( MAXALB , NLUS )
2166 CALL wrf_dm_bcast_real ( TOPT_DATA , 1 )
2167 CALL wrf_dm_bcast_real ( CMCMAX_DATA , 1 )
2168 CALL wrf_dm_bcast_real ( CFACTR_DATA , 1 )
2169 CALL wrf_dm_bcast_real ( RSMAX_DATA , 1 )
2170 CALL wrf_dm_bcast_integer ( BARE , 1 )
2171 CALL wrf_dm_bcast_integer ( NATURAL , 1 )
2172 CALL wrf_dm_bcast_integer ( LCZ_1 , 1 )
2173 CALL wrf_dm_bcast_integer ( LCZ_2 , 1 )
2174 CALL wrf_dm_bcast_integer ( LCZ_3 , 1 )
2175 CALL wrf_dm_bcast_integer ( LCZ_4 , 1 )
2176 CALL wrf_dm_bcast_integer ( LCZ_5 , 1 )
2177 CALL wrf_dm_bcast_integer ( LCZ_6 , 1 )
2178 CALL wrf_dm_bcast_integer ( LCZ_7 , 1 )
2179 CALL wrf_dm_bcast_integer ( LCZ_8 , 1 )
2180 CALL wrf_dm_bcast_integer ( LCZ_9 , 1 )
2181 CALL wrf_dm_bcast_integer ( LCZ_10 , 1 )
2182 CALL wrf_dm_bcast_integer ( LCZ_11 , 1 )
2185 !-----READ IN SOIL PROPERTIES FROM SOILPARM.TBL
2187 IF ( wrf_dm_on_monitor() ) THEN
2188 OPEN(19, FILE='SOILPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
2189 IF(ierr .NE. OPEN_OK ) THEN
2190 WRITE(message,FMT='(A)') &
2191 'module_sf_noahlsm.F: soil_veg_gen_parm: failure opening SOILPARM.TBL'
2192 CALL wrf_error_fatal ( message )
2195 WRITE(mess,*) 'INPUT SOIL TEXTURE CLASSIFICATION = ', TRIM ( MMINSL )
2196 CALL wrf_message( mess )
2201 READ (19,2000,END=2003)SLTYPE
2203 READ (19,*)SLCATS,IINDEX
2204 IF(SLTYPE.EQ.MMINSL)THEN
2205 WRITE( mess , * ) 'SOIL TEXTURE CLASSIFICATION = ', TRIM ( SLTYPE ) , ' FOUND', &
2206 SLCATS,' CATEGORIES'
2207 CALL wrf_message ( mess )
2210 ! prevent possible array overwrite, Bill Bovermann, IBM, May 6, 2008
2211 IF ( SIZE(BB ) < SLCATS .OR. &
2212 SIZE(DRYSMC) < SLCATS .OR. &
2213 SIZE(F11 ) < SLCATS .OR. &
2214 SIZE(MAXSMC) < SLCATS .OR. &
2215 SIZE(REFSMC) < SLCATS .OR. &
2216 SIZE(SATPSI) < SLCATS .OR. &
2217 SIZE(SATDK ) < SLCATS .OR. &
2218 SIZE(SATDW ) < SLCATS .OR. &
2219 SIZE(WLTSMC) < SLCATS .OR. &
2220 SIZE(QTZ ) < SLCATS ) THEN
2221 CALL wrf_error_fatal('Table sizes too small for value of SLCATS in module_sf_noahdrv.F')
2223 IF(SLTYPE.EQ.MMINSL)THEN
2225 READ (19,*) IINDEX,BB(LC),DRYSMC(LC),F11(LC),MAXSMC(LC),&
2226 REFSMC(LC),SATPSI(LC),SATDK(LC), SATDW(LC), &
2236 CALL wrf_dm_bcast_integer ( LUMATCH , 1 )
2237 CALL wrf_dm_bcast_string ( SLTYPE , 4 )
2238 CALL wrf_dm_bcast_string ( MMINSL , 4 ) ! since this is reset above, see oct2 ^
2239 CALL wrf_dm_bcast_integer ( SLCATS , 1 )
2240 CALL wrf_dm_bcast_integer ( IINDEX , 1 )
2241 CALL wrf_dm_bcast_real ( BB , NSLTYPE )
2242 CALL wrf_dm_bcast_real ( DRYSMC , NSLTYPE )
2243 CALL wrf_dm_bcast_real ( F11 , NSLTYPE )
2244 CALL wrf_dm_bcast_real ( MAXSMC , NSLTYPE )
2245 CALL wrf_dm_bcast_real ( REFSMC , NSLTYPE )
2246 CALL wrf_dm_bcast_real ( SATPSI , NSLTYPE )
2247 CALL wrf_dm_bcast_real ( SATDK , NSLTYPE )
2248 CALL wrf_dm_bcast_real ( SATDW , NSLTYPE )
2249 CALL wrf_dm_bcast_real ( WLTSMC , NSLTYPE )
2250 CALL wrf_dm_bcast_real ( QTZ , NSLTYPE )
2252 IF(LUMATCH.EQ.0)THEN
2253 CALL wrf_message( 'SOIl TEXTURE IN INPUT FILE DOES NOT ' )
2254 CALL wrf_message( 'MATCH SOILPARM TABLE' )
2255 CALL wrf_error_fatal ( 'INCONSISTENT OR MISSING SOILPARM FILE' )
2259 !-----READ IN GENERAL PARAMETERS FROM GENPARM.TBL
2261 IF ( wrf_dm_on_monitor() ) THEN
2262 OPEN(19, FILE='GENPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
2263 IF(ierr .NE. OPEN_OK ) THEN
2264 WRITE(message,FMT='(A)') &
2265 'module_sf_noahlsm.F: soil_veg_gen_parm: failure opening GENPARM.TBL'
2266 CALL wrf_error_fatal ( message )
2271 READ (19,*) NUM_SLOPE
2274 ! prevent possible array overwrite, Bill Bovermann, IBM, May 6, 2008
2275 IF ( SIZE(slope_data) < NUM_SLOPE ) THEN
2276 CALL wrf_error_fatal('NUM_SLOPE too large for slope_data array in module_sf_noahdrv')
2280 READ (19,*)SLOPE_DATA(LC)
2284 READ (19,*)SBETA_DATA
2286 READ (19,*)FXEXP_DATA
2288 READ (19,*)CSOIL_DATA
2290 READ (19,*)SALP_DATA
2292 READ (19,*)REFDK_DATA
2294 READ (19,*)REFKDT_DATA
2296 READ (19,*)FRZK_DATA
2298 READ (19,*)ZBOT_DATA
2300 READ (19,*)CZIL_DATA
2302 READ (19,*)SMLOW_DATA
2304 READ (19,*)SMHIGH_DATA
2306 READ (19,*)LVCOEF_DATA
2310 CALL wrf_dm_bcast_integer ( NUM_SLOPE , 1 )
2311 CALL wrf_dm_bcast_integer ( SLPCATS , 1 )
2312 CALL wrf_dm_bcast_real ( SLOPE_DATA , NSLOPE )
2313 CALL wrf_dm_bcast_real ( SBETA_DATA , 1 )
2314 CALL wrf_dm_bcast_real ( FXEXP_DATA , 1 )
2315 CALL wrf_dm_bcast_real ( CSOIL_DATA , 1 )
2316 CALL wrf_dm_bcast_real ( SALP_DATA , 1 )
2317 CALL wrf_dm_bcast_real ( REFDK_DATA , 1 )
2318 CALL wrf_dm_bcast_real ( REFKDT_DATA , 1 )
2319 CALL wrf_dm_bcast_real ( FRZK_DATA , 1 )
2320 CALL wrf_dm_bcast_real ( ZBOT_DATA , 1 )
2321 CALL wrf_dm_bcast_real ( CZIL_DATA , 1 )
2322 CALL wrf_dm_bcast_real ( SMLOW_DATA , 1 )
2323 CALL wrf_dm_bcast_real ( SMHIGH_DATA , 1 )
2324 CALL wrf_dm_bcast_real ( LVCOEF_DATA , 1 )
2327 !-----------------------------------------------------------------
2328 END SUBROUTINE SOIL_VEG_GEN_PARM
2329 !-----------------------------------------------------------------
2331 !===========================================================================
2333 ! subroutine lsm_mosaic: a tiling approach for Noah LSM
2335 !===========================================================================
2337 SUBROUTINE lsm_mosaic(DZ8W,QV3D,P8W3D,T3D,TSK, &
2338 HFX,QFX,LH,GRDFLX, QGH,GSW,SWDOWN,GLW,SMSTAV,SMSTOT, &
2339 SFCRUNOFF, UDRUNOFF,IVGTYP,ISLTYP,ISURBAN,ISICE,VEGFRA, &
2340 ALBEDO,ALBBCK,ZNT,Z0,TMN,XLAND,XICE,EMISS,EMBCK, &
2341 SNOWC,QSFC,RAINBL,MMINLU, &
2342 num_soil_layers,DT,DZS,ITIMESTEP, &
2343 SMOIS,TSLB,SNOW,CANWAT, &
2344 CHS,CHS2,CQS2,CPM,ROVCP,SR,chklowq,lai,qz0, & !H
2348 SNOALB,SHDMIN,SHDMAX, & !I
2355 RDLAI2D,USEMONALB, &
2357 NOAHRES,OPT_THCND, &
2358 NLCAT,landusef,landusef2, & ! danli mosaic
2359 sf_surface_mosaic,mosaic_cat,mosaic_cat_index, & ! danli mosaic
2360 TSK_mosaic,QSFC_mosaic, & ! danli mosaic
2361 TSLB_mosaic,SMOIS_mosaic,SH2O_mosaic, & ! danli mosaic
2362 CANWAT_mosaic,SNOW_mosaic, & ! danli mosaic
2363 SNOWH_mosaic,SNOWC_mosaic, & ! danli mosaic
2364 ALBEDO_mosaic,ALBBCK_mosaic, & ! danli mosaic
2365 EMISS_mosaic, EMBCK_mosaic, & ! danli mosaic
2366 ZNT_mosaic, Z0_mosaic, & ! danli mosaic
2367 HFX_mosaic,QFX_mosaic, & ! danli mosaic
2368 LH_mosaic, GRDFLX_mosaic, SNOTIME_mosaic, & ! danli mosaic
2369 RC_mosaic, LAI_mosaic, &
2371 ua_phys,flx4_2d,fvb_2d,fbur_2d,fgsn_2d, &
2372 ids,ide, jds,jde, kds,kde, &
2373 ims,ime, jms,jme, kms,kme, &
2374 its,ite, jts,jte, kts,kte, &
2376 CMR_SFCDIF,CHR_SFCDIF,CMC_SFCDIF,CHC_SFCDIF, &
2377 CMGR_SFCDIF,CHGR_SFCDIF, &
2379 TR_URB2D,TB_URB2D,TG_URB2D,TC_URB2D,QC_URB2D, & !H urban
2380 UC_URB2D, & !H urban
2381 XXXR_URB2D,XXXB_URB2D,XXXG_URB2D,XXXC_URB2D, & !H urban
2382 TRL_URB3D,TBL_URB3D,TGL_URB3D, & !H urban
2383 SH_URB2D,LH_URB2D,G_URB2D,RN_URB2D,TS_URB2D, & !H urban
2384 TR_URB2D_mosaic,TB_URB2D_mosaic, & !H urban danli mosaic
2385 TG_URB2D_mosaic,TC_URB2D_mosaic, & !H urban danli mosaic
2386 QC_URB2D_mosaic,UC_URB2D_mosaic, & !H urban danli mosaic
2387 TRL_URB3D_mosaic,TBL_URB3D_mosaic, & !H urban danli mosaic
2388 TGL_URB3D_mosaic, & !H urban danli mosaic
2389 SH_URB2D_mosaic,LH_URB2D_mosaic, & !H urban danli mosaic
2390 G_URB2D_mosaic,RN_URB2D_mosaic, & !H urban danli mosaic
2391 TS_URB2D_mosaic, & !H urban danli mosaic
2392 TS_RUL2D_mosaic, & !H urban danli mosaic
2393 PSIM_URB2D,PSIH_URB2D,U10_URB2D,V10_URB2D, & !O urban
2394 GZ1OZ0_URB2D, AKMS_URB2D, & !O urban
2395 TH2_URB2D,Q2_URB2D, UST_URB2D, & !O urban
2396 DECLIN_URB,COSZ_URB2D,OMG_URB2D, & !I urban
2397 XLAT_URB2D, & !I urban
2398 num_roof_layers, num_wall_layers, & !I urban
2399 num_road_layers, DZR, DZB, DZG, & !I urban
2400 CMCR_URB2D,TGR_URB2D,TGRL_URB3D,SMR_URB3D, & !H urban
2401 julian,julyr, & !H urban
2402 DRELR_URB2D,DRELB_URB2D,DRELG_URB2D, & !H urban
2403 FLXHUMR_URB2D,FLXHUMB_URB2D,FLXHUMG_URB2D, & !H urban
2404 FRC_URB2D,UTYPE_URB2D, & !O
2405 num_urban_ndm, & !I multi-layer urban
2406 urban_map_zrd, & !I multi-layer urban
2407 urban_map_zwd, & !I multi-layer urban
2408 urban_map_gd, & !I multi-layer urban
2409 urban_map_zd, & !I multi-layer urban
2410 urban_map_zdf, & !I multi-layer urban
2411 urban_map_bd, & !I multi-layer urban
2412 urban_map_wd, & !I multi-layer urban
2413 urban_map_gbd, & !I multi-layer urban
2414 urban_map_fbd, & !I multi-layer urban
2415 urban_map_zgrd, & !I multi-layer urban
2416 num_urban_hi, & !I multi-layer urban
2417 use_wudapt_lcz, & !I wudapt
2418 tsk_rural_bep, & !H multi-layer urban
2419 trb_urb4d,tw1_urb4d,tw2_urb4d,tgb_urb4d, & !H multi-layer urban
2420 tlev_urb3d,qlev_urb3d, & !H multi-layer urban
2421 tw1lev_urb3d,tw2lev_urb3d, & !H multi-layer urban
2422 tglev_urb3d,tflev_urb3d, & !H multi-layer urban
2423 sf_ac_urb3d,lf_ac_urb3d,cm_ac_urb3d, & !H multi-layer urban
2424 sfvent_urb3d,lfvent_urb3d, & !H multi-layer urban
2425 sfwin1_urb3d,sfwin2_urb3d, & !H multi-layer urban
2426 sfw1_urb3d,sfw2_urb3d,sfr_urb3d,sfg_urb3d, & !H multi-layer urban
2427 ep_pv_urb3d,t_pv_urb3d, & !RMS
2428 trv_urb4d,qr_urb4d,qgr_urb3d,tgr_urb3d, & !RMS
2429 drain_urb4d,draingr_urb3d,sfrv_urb3d, & !RMS
2430 lfrv_urb3d,dgr_urb3d,dg_urb3d,lfr_urb3d,lfg_urb3d,& !RMS
2431 lp_urb2d,hi_urb2d,lb_urb2d,hgt_urb2d, & !H multi-layer urban
2432 mh_urb2d,stdh_urb2d,lf_urb2d, & !SLUCM
2433 th_phy,rho,p_phy,ust, & !I multi-layer urban
2434 gmt,julday,xlong,xlat, & !I multi-layer urban
2435 a_u_bep,a_v_bep,a_t_bep,a_q_bep, & !O multi-layer urban
2436 a_e_bep,b_u_bep,b_v_bep, & !O multi-layer urban
2437 b_t_bep,b_q_bep,b_e_bep,dlg_bep, & !O multi-layer urban
2438 dl_u_bep,sf_bep,vl_bep & !O multi-layer urban
2440 ,sfcheadrt,INFXSRT, soldrain & !hydro
2442 ,SDA_HFX, SDA_QFX, HFX_BOTH, QFX_BOTH, QNORM, fasdas & !fasdas
2447 !----------------------------------------------------------------
2449 !----------------------------------------------------------------
2450 !----------------------------------------------------------------
2451 ! --- atmospheric (WRF generic) variables
2452 !-- DT time step (seconds)
2453 !-- DZ8W thickness of layers (m)
2454 !-- T3D temperature (K)
2455 !-- QV3D 3D water vapor mixing ratio (Kg/Kg)
2456 !-- P3D 3D pressure (Pa)
2457 !-- FLHC exchange coefficient for heat (m/s)
2458 !-- FLQC exchange coefficient for moisture (m/s)
2459 !-- PSFC surface pressure (Pa)
2460 !-- XLAND land mask (1 for land, 2 for water)
2461 !-- QGH saturated mixing ratio at 2 meter
2462 !-- GSW downward short wave flux at ground surface (W/m^2)
2463 !-- GLW downward long wave flux at ground surface (W/m^2)
2464 !-- History variables
2465 !-- CANWAT canopy moisture content (mm)
2466 !-- TSK surface temperature (K)
2467 !-- TSLB soil temp (k)
2468 !-- SMOIS total soil moisture content (volumetric fraction)
2469 !-- SH2O unfrozen soil moisture content (volumetric fraction)
2470 ! note: frozen soil moisture (i.e., soil ice) = SMOIS - SH2O
2471 !-- SNOWH actual snow depth (m)
2472 !-- SNOW liquid water-equivalent snow depth (m)
2473 !-- ALBEDO time-varying surface albedo including snow effect (unitless fraction)
2474 !-- ALBBCK background surface albedo (unitless fraction)
2475 !-- CHS surface exchange coefficient for heat and moisture (m s-1);
2476 !-- CHS2 2m surface exchange coefficient for heat (m s-1);
2477 !-- CQS2 2m surface exchange coefficient for moisture (m s-1);
2478 ! --- soil variables
2479 !-- num_soil_layers the number of soil layers
2480 !-- ZS depths of centers of soil layers (m)
2481 !-- DZS thicknesses of soil layers (m)
2482 !-- SLDPTH thickness of each soil layer (m, same as DZS)
2483 !-- TMN soil temperature at lower boundary (K)
2484 !-- SMCWLT wilting point (volumetric)
2485 !-- SMCDRY dry soil moisture threshold where direct evap from
2486 ! top soil layer ends (volumetric)
2487 !-- SMCREF soil moisture threshold below which transpiration begins to
2488 ! stress (volumetric)
2489 !-- SMCMAX porosity, i.e. saturated value of soil moisture (volumetric)
2490 !-- NROOT number of root layers, a function of veg type, determined
2491 ! in subroutine redprm.
2492 !-- SMSTAV Soil moisture availability for evapotranspiration (
2493 ! fraction between SMCWLT and SMCMXA)
2494 !-- SMSTOT Total soil moisture content frozen+unfrozen) in the soil column (mm)
2495 ! --- snow variables
2496 !-- SNOWC fraction snow coverage (0-1.0)
2497 ! --- vegetation variables
2498 !-- SNOALB upper bound on maximum albedo over deep snow
2499 !-- SHDMIN minimum areal fractional coverage of annual green vegetation
2500 !-- SHDMAX maximum areal fractional coverage of annual green vegetation
2501 !-- XLAI leaf area index (dimensionless)
2502 !-- Z0BRD Background fixed roughness length (M)
2503 !-- Z0 Background vroughness length (M) as function
2504 !-- ZNT Time varying roughness length (M) as function
2505 !-- ALBD(IVGTPK,ISN) background albedo reading from a table
2507 !-- HFX upward heat flux at the surface (W/m^2)
2508 !-- QFX upward moisture flux at the surface (kg/m^2/s)
2509 !-- LH upward moisture flux at the surface (W m-2)
2510 !-- GRDFLX(I,J) ground heat flux (W m-2)
2511 !-- FDOWN radiation forcing at the surface (W m-2) = SOLDN*(1-alb)+LWDN
2512 !----------------------------------------------------------------------------
2513 !-- EC canopy water evaporation ((W m-2)
2514 !-- EDIR direct soil evaporation (W m-2)
2515 !-- ET plant transpiration from a particular root layer (W m-2)
2516 !-- ETT total plant transpiration (W m-2)
2517 !-- ESNOW sublimation from (or deposition to if <0) snowpack (W m-2)
2518 !-- DRIP through-fall of precip and/or dew in excess of canopy
2519 ! water-holding capacity (m)
2520 !-- DEW dewfall (or frostfall for t<273.15) (M)
2521 !-- SMAV Soil Moisture Availability for each layer, as a fraction
2522 ! between SMCWLT and SMCMAX (dimensionless fraction)
2523 ! ----------------------------------------------------------------------
2524 !-- BETA ratio of actual/potential evap (dimensionless)
2525 !-- ETP potential evaporation (W m-2)
2526 ! ----------------------------------------------------------------------
2527 !-- FLX1 precip-snow sfc (W m-2)
2528 !-- FLX2 freezing rain latent heat flux (W m-2)
2529 !-- FLX3 phase-change heat flux from snowmelt (W m-2)
2530 ! ----------------------------------------------------------------------
2531 !-- ACSNOM snow melt (mm) (water equivalent)
2532 !-- ACSNOW accumulated snow fall (mm) (water equivalent)
2533 !-- SNOPCX snow phase change heat flux (W/m^2)
2534 !-- POTEVP accumulated potential evaporation (m)
2535 !-- RIB Documentation needed!!!
2536 ! ----------------------------------------------------------------------
2537 !-- RUNOFF1 surface runoff (m s-1), not infiltrating the surface
2538 !-- RUNOFF2 subsurface runoff (m s-1), drainage out bottom of last
2539 ! soil layer (baseflow)
2540 ! important note: here RUNOFF2 is actually the sum of RUNOFF2 and RUNOFF3
2541 !-- RUNOFF3 numerical trunctation in excess of porosity (smcmax)
2542 ! for a given soil layer at the end of a time step (m s-1).
2543 !SFCRUNOFF Surface Runoff (mm)
2544 !UDRUNOFF Total Underground Runoff (mm), which is the sum of RUNOFF2 and RUNOFF3
2545 ! ----------------------------------------------------------------------
2546 !-- RC canopy resistance (s m-1)
2547 !-- PC plant coefficient (unitless fraction, 0-1) where PC*ETP = actual transp
2548 !-- RSMIN minimum canopy resistance (s m-1)
2549 !-- RCS incoming solar rc factor (dimensionless)
2550 !-- RCT air temperature rc factor (dimensionless)
2551 !-- RCQ atmos vapor pressure deficit rc factor (dimensionless)
2552 !-- RCSOIL soil moisture rc factor (dimensionless)
2554 !-- EMISS surface emissivity (between 0 and 1)
2555 !-- EMBCK Background surface emissivity (between 0 and 1)
2558 ! (R_d/R_v) (dimensionless)
2559 !-- ids start index for i in domain
2560 !-- ide end index for i in domain
2561 !-- jds start index for j in domain
2562 !-- jde end index for j in domain
2563 !-- kds start index for k in domain
2564 !-- kde end index for k in domain
2565 !-- ims start index for i in memory
2566 !-- ime end index for i in memory
2567 !-- jms start index for j in memory
2568 !-- jme end index for j in memory
2569 !-- kms start index for k in memory
2570 !-- kme end index for k in memory
2571 !-- its start index for i in tile
2572 !-- ite end index for i in tile
2573 !-- jts start index for j in tile
2574 !-- jte end index for j in tile
2575 !-- kts start index for k in tile
2576 !-- kte end index for k in tile
2578 !-- SR fraction of frozen precip (0.0 to 1.0)
2579 !----------------------------------------------------------------
2583 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
2584 ims,ime, jms,jme, kms,kme, &
2585 its,ite, jts,jte, kts,kte
2587 INTEGER, INTENT(IN ) :: sf_urban_physics !urban
2588 INTEGER, INTENT(IN ) :: isurban
2589 INTEGER, INTENT(IN ) :: isice
2590 INTEGER, INTENT(IN ) :: julian,julyr
2592 !added by Wei Yu for routing
2594 REAL, DIMENSION( ims:ime, jms:jme ) , &
2595 INTENT(INOUT) :: sfcheadrt,INFXSRT,soldrain
2600 ! new local vars for hydro
2601 REAL :: etpnd1_hydro,sfcheadrt_hydro,infxsrt_hydro
2603 REAL, DIMENSION( ims:ime, jms:jme ) , &
2604 INTENT(IN ) :: TMN, &
2612 SWDOWN, & !added 10 jan 2007
2616 REAL, DIMENSION( ims:ime, jms:jme ) , &
2617 INTENT(INOUT) :: ALBBCK, &
2619 EMBCK ! danli mosaic
2621 CHARACTER(LEN=*), INTENT(IN ) :: MMINLU
2623 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
2624 INTENT(IN ) :: QV3D, &
2628 REAL, DIMENSION( ims:ime, jms:jme ) , &
2629 INTENT(IN ) :: QGH, &
2632 INTEGER, DIMENSION( ims:ime, jms:jme ) , &
2633 INTENT(IN ) :: ISLTYP
2635 INTEGER, DIMENSION( ims:ime, jms:jme ) , &
2636 INTENT(INOUT ) :: IVGTYP ! for mosaic danli
2638 INTEGER, INTENT(IN) :: num_soil_layers,ITIMESTEP
2640 REAL, INTENT(IN ) :: DT,ROVCP
2642 REAL, DIMENSION(1:num_soil_layers), INTENT(IN)::DZS
2646 REAL, DIMENSION( ims:ime , 1:num_soil_layers, jms:jme ), &
2647 INTENT(INOUT) :: SMOIS, & ! total soil moisture
2648 SH2O, & ! new soil liquid
2651 REAL, DIMENSION( ims:ime , 1:num_soil_layers, jms:jme ), &
2652 INTENT(OUT) :: SMCREL
2654 REAL, DIMENSION( ims:ime, jms:jme ) , &
2655 INTENT(INOUT) :: TSK, & !was TGB (temperature)
2681 REAL, DIMENSION( ims:ime, jms:jme ) , &
2682 INTENT(OUT) :: NOAHRES
2683 INTEGER, INTENT(IN) :: OPT_THCND
2686 LOGICAL, INTENT(IN) :: UA_PHYS
2687 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: FLX4_2D,FVB_2D,FBUR_2D,FGSN_2D
2688 REAL :: FLX4,FVB,FBUR,FGSN
2690 REAL, DIMENSION( ims:ime, jms:jme ) , &
2691 INTENT(OUT) :: CHKLOWQ
2692 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: LAI
2693 REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: QZ0
2694 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: RC2, XLAI2
2696 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CMR_SFCDIF
2697 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CHR_SFCDIF
2698 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CMGR_SFCDIF
2699 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CHGR_SFCDIF
2700 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CMC_SFCDIF
2701 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CHC_SFCDIF
2702 ! Local variables (moved here from driver to make routine thread safe, 20031007 jm)
2704 REAL, DIMENSION(1:num_soil_layers) :: ET
2706 REAL, DIMENSION(1:num_soil_layers) :: SMAV
2708 REAL :: BETA, ETP, SSOIL,EC, EDIR, ESNOW, ETT, &
2709 FLX1,FLX2,FLX3, DRIP,DEW,FDOWN,RC,PC,RSMIN,XLAI, &
2710 ! RCS,RCT,RCQ,RCSOIL
2711 RCS,RCT,RCQ,RCSOIL,FFROZP
2713 LOGICAL, INTENT(IN ) :: myj,frpcpn
2715 ! DECLARATIONS - LOGICAL
2716 ! ----------------------------------------------------------------------
2717 LOGICAL, PARAMETER :: LOCAL=.false.
2718 LOGICAL :: FRZGRA, SNOWNG
2722 ! ----------------------------------------------------------------------
2723 ! DECLARATIONS - INTEGER
2724 ! ----------------------------------------------------------------------
2725 INTEGER :: I,J, ICE,NSOIL,SLOPETYP,SOILTYP,VEGTYP
2729 ! ----------------------------------------------------------------------
2730 ! DECLARATIONS - REAL
2731 ! ----------------------------------------------------------------------
2733 REAL :: SHMIN,SHMAX,DQSDT2,LWDN,PRCP,PRCPRAIN, &
2734 Q2SAT,Q2SATI,SFCPRS,SFCSPD,SFCTMP,SHDFAC,SNOALB1, &
2735 SOLDN,TBOT,ZLVL, Q2K,ALBBRD, ALBEDOK, ETA, ETA_KINEMATIC, &
2737 Z0K,RUNOFF1,RUNOFF2,RUNOFF3,SHEAT,SOLNET,E2SAT,SFCTSNO, &
2738 ! mek, WRF testing, expanded diagnostics
2739 SOLUP,LWUP,RNET,RES,Q1SFC,TAIRV,SATFLG
2742 ! MEK JUL2007 for pot. evap.
2748 REAL :: SNCOVR,SNEQV,SNOWHK,CMC, CHK,TH2
2750 REAL :: SMCDRY,SMCMAX,SMCREF,SMCWLT,SNOMLT,SOILM,SOILW,Q1,T1
2751 REAL :: SNOTIME1 ! LSTSNW1 INITIAL NUMBER OF TIMESTEPS SINCE LAST SNOWFALL
2755 REAL :: COSZ, SOLARDIRECT
2757 REAL, DIMENSION(1:num_soil_layers):: SLDPTH, STC,SMC,SWC
2759 REAL, DIMENSION(1:num_soil_layers) :: ZSOIL, RTDIS
2760 REAL, PARAMETER :: TRESH=.95E0, A2=17.67,A3=273.15,A4=29.65, &
2761 T0=273.16E0, ELWV=2.50E6, A23M4=A2*(A3-A4)
2763 REAL, PARAMETER :: ROW=1.E3,ELIW=XLF,ROWLIW=ROW*ELIW
2765 ! ----------------------------------------------------------------------
2766 ! DECLARATIONS START - urban
2767 ! ----------------------------------------------------------------------
2769 ! input variables surface_driver --> lsm
2770 INTEGER, INTENT(IN) :: num_roof_layers
2771 INTEGER, INTENT(IN) :: num_wall_layers
2772 INTEGER, INTENT(IN) :: num_road_layers
2773 REAL, OPTIONAL, DIMENSION(1:num_roof_layers), INTENT(IN) :: DZR
2774 REAL, OPTIONAL, DIMENSION(1:num_wall_layers), INTENT(IN) :: DZB
2775 REAL, OPTIONAL, DIMENSION(1:num_road_layers), INTENT(IN) :: DZG
2776 REAL, OPTIONAL, INTENT(IN) :: DECLIN_URB
2777 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: COSZ_URB2D
2778 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: OMG_URB2D
2779 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: XLAT_URB2D
2780 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: U_PHY
2781 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: V_PHY
2782 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: TH_PHY
2783 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: P_PHY
2784 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: RHO
2785 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: UST
2787 LOGICAL, intent(in) :: rdlai2d
2788 LOGICAL, intent(in) :: USEMONALB
2790 ! input variables lsm --> urban
2791 INTEGER :: UTYPE_URB ! urban type [urban=1, suburban=2, rural=3]
2792 REAL :: TA_URB ! potential temp at 1st atmospheric level [K]
2793 REAL :: QA_URB ! mixing ratio at 1st atmospheric level [kg/kg]
2794 REAL :: UA_URB ! wind speed at 1st atmospheric level [m/s]
2795 REAL :: U1_URB ! u at 1st atmospheric level [m/s]
2796 REAL :: V1_URB ! v at 1st atmospheric level [m/s]
2797 REAL :: SSG_URB ! downward total short wave radiation [W/m/m]
2798 REAL :: LLG_URB ! downward long wave radiation [W/m/m]
2799 REAL :: RAIN_URB ! precipitation [mm/h]
2800 REAL :: RHOO_URB ! air density [kg/m^3]
2801 REAL :: ZA_URB ! first atmospheric level [m]
2802 REAL :: DELT_URB ! time step [s]
2803 REAL :: SSGD_URB ! downward direct short wave radiation [W/m/m]
2804 REAL :: SSGQ_URB ! downward diffuse short wave radiation [W/m/m]
2805 REAL :: XLAT_URB ! latitude [deg]
2806 REAL :: COSZ_URB ! cosz
2807 REAL :: OMG_URB ! hour angle
2808 REAL :: ZNT_URB ! roughness length [m]
2819 REAL, DIMENSION(1:num_roof_layers) :: TRL_URB ! roof layer temp [K]
2820 REAL, DIMENSION(1:num_wall_layers) :: TBL_URB ! wall layer temp [K]
2821 REAL, DIMENSION(1:num_road_layers) :: TGL_URB ! road layer temp [K]
2822 LOGICAL :: LSOLAR_URB
2824 !===Yang,2014/10/08,hydrological variable for single layer UCM===
2825 INTEGER :: jmonth, jday, tloc
2826 INTEGER :: IRIOPTION, USOIL, DSOIL
2836 REAL, DIMENSION(1:num_roof_layers) :: SMR_URB ! green roof layer moisture
2837 REAL, DIMENSION(1:num_roof_layers) :: TGRL_URB ! green roof layer temp [K]
2839 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: DRELR_URB2D
2840 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: DRELB_URB2D
2841 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: DRELG_URB2D
2842 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: FLXHUMR_URB2D
2843 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: FLXHUMB_URB2D
2844 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: FLXHUMG_URB2D
2845 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CMCR_URB2D
2846 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TGR_URB2D
2847 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_roof_layers, jms:jme ), INTENT(INOUT) :: TGRL_URB3D
2848 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_roof_layers, jms:jme ), INTENT(INOUT) :: SMR_URB3D
2850 ! state variable surface_driver <--> lsm <--> urban
2851 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TR_URB2D
2852 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TB_URB2D
2853 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TG_URB2D
2854 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TC_URB2D
2855 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: QC_URB2D
2856 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: UC_URB2D
2857 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXR_URB2D
2858 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXB_URB2D
2859 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXG_URB2D
2860 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXC_URB2D
2861 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: SH_URB2D
2862 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: LH_URB2D
2863 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: G_URB2D
2864 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: RN_URB2D
2866 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TS_URB2D
2868 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_roof_layers, jms:jme ), INTENT(INOUT) :: TRL_URB3D
2869 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_wall_layers, jms:jme ), INTENT(INOUT) :: TBL_URB3D
2870 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_road_layers, jms:jme ), INTENT(INOUT) :: TGL_URB3D
2872 ! output variable lsm --> surface_driver
2873 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: PSIM_URB2D
2874 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: PSIH_URB2D
2875 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: GZ1OZ0_URB2D
2876 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: U10_URB2D
2877 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: V10_URB2D
2878 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: TH2_URB2D
2879 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: Q2_URB2D
2881 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: AKMS_URB2D
2883 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: UST_URB2D
2884 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: FRC_URB2D ! change this to inout, danli mosaic
2885 INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: UTYPE_URB2D
2887 ! output variables urban --> lsm
2888 REAL :: TS_URB ! surface radiative temperature [K]
2889 REAL :: QS_URB ! surface humidity [-]
2890 REAL :: SH_URB ! sensible heat flux [W/m/m]
2891 REAL :: LH_URB ! latent heat flux [W/m/m]
2892 REAL :: LH_KINEMATIC_URB ! latent heat flux, kinetic [kg/m/m/s]
2893 REAL :: SW_URB ! upward short wave radiation flux [W/m/m]
2894 REAL :: ALB_URB ! time-varying albedo [fraction]
2895 REAL :: LW_URB ! upward long wave radiation flux [W/m/m]
2896 REAL :: G_URB ! heat flux into the ground [W/m/m]
2897 REAL :: RN_URB ! net radiation [W/m/m]
2898 REAL :: PSIM_URB ! shear f for momentum [-]
2899 REAL :: PSIH_URB ! shear f for heat [-]
2900 REAL :: GZ1OZ0_URB ! shear f for heat [-]
2901 REAL :: U10_URB ! wind u component at 10 m [m/s]
2902 REAL :: V10_URB ! wind v component at 10 m [m/s]
2903 REAL :: TH2_URB ! potential temperature at 2 m [K]
2904 REAL :: Q2_URB ! humidity at 2 m [-]
2908 ! NUDAPT Parameters urban --> lam
2913 REAL, DIMENSION(4) :: lf_urb
2914 ! Variables for multi-layer UCM (Martilli et al. 2002)
2915 REAL, OPTIONAL, INTENT(IN ) :: GMT
2916 INTEGER, OPTIONAL, INTENT(IN ) :: JULDAY
2917 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) ::XLAT, XLONG
2918 INTEGER, INTENT(IN ) :: num_urban_ndm
2919 INTEGER, INTENT(IN ) :: urban_map_zrd
2920 INTEGER, INTENT(IN ) :: urban_map_zwd
2921 INTEGER, INTENT(IN ) :: urban_map_gd
2922 INTEGER, INTENT(IN ) :: urban_map_zd
2923 INTEGER, INTENT(IN ) :: urban_map_zdf
2924 INTEGER, INTENT(IN ) :: urban_map_bd
2925 INTEGER, INTENT(IN ) :: urban_map_wd
2926 INTEGER, INTENT(IN ) :: urban_map_gbd
2927 INTEGER, INTENT(IN ) :: urban_map_fbd
2928 INTEGER, INTENT(IN ) :: urban_map_zgrd
2929 INTEGER, INTENT(IN ) :: NUM_URBAN_HI
2930 INTEGER, INTENT(IN ) :: use_wudapt_lcz
2931 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: tsk_rural_bep
2932 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zrd, jms:jme ), INTENT(INOUT) :: trb_urb4d
2933 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zwd, jms:jme ), INTENT(INOUT) :: tw1_urb4d
2934 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zwd, jms:jme ), INTENT(INOUT) :: tw2_urb4d
2935 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_gd , jms:jme ), INTENT(INOUT) :: tgb_urb4d
2936 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_bd , jms:jme ), INTENT(INOUT) :: tlev_urb3d
2937 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_bd , jms:jme ), INTENT(INOUT) :: qlev_urb3d
2938 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_wd , jms:jme ), INTENT(INOUT) :: tw1lev_urb3d
2939 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_wd , jms:jme ), INTENT(INOUT) :: tw2lev_urb3d
2940 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_gbd, jms:jme ), INTENT(INOUT) :: tglev_urb3d
2941 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_fbd, jms:jme ), INTENT(INOUT) :: tflev_urb3d
2942 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: lf_ac_urb3d
2943 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: sf_ac_urb3d
2944 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: cm_ac_urb3d
2945 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: sfvent_urb3d
2946 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: lfvent_urb3d
2947 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_wd , jms:jme ), INTENT(INOUT) :: sfwin1_urb3d
2948 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_wd , jms:jme ), INTENT(INOUT) :: sfwin2_urb3d
2949 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zd , jms:jme ), INTENT(INOUT) :: sfw1_urb3d
2950 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zd , jms:jme ), INTENT(INOUT) :: sfw2_urb3d
2951 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme ), INTENT(INOUT) :: sfr_urb3d
2952 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_ndm, jms:jme ), INTENT(INOUT) :: sfg_urb3d
2953 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: ep_pv_urb3d !GRZ
2954 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme ), INTENT(INOUT) :: t_pv_urb3d !GRZ
2955 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zgrd, jms:jme ),INTENT(INOUT) :: trv_urb4d !GRZ
2956 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zgrd, jms:jme ),INTENT(INOUT) :: qr_urb4d !GRZ
2957 REAL, OPTIONAL, DIMENSION( ims:ime,jms:jme ), INTENT(INOUT) :: qgr_urb3d !GRZ
2958 REAL, OPTIONAL, DIMENSION( ims:ime,jms:jme ), INTENT(INOUT) :: tgr_urb3d !GRZ
2959 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme ),INTENT(INOUT) :: drain_urb4d !GRZ
2960 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: draingr_urb3d !GRZ
2961 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme ),INTENT(INOUT) :: sfrv_urb3d !GRZ
2962 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme ),INTENT(INOUT) :: lfrv_urb3d !GRZ
2963 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme ),INTENT(INOUT) :: dgr_urb3d !GRZ
2964 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_ndm, jms:jme ),INTENT(INOUT) :: dg_urb3d !GRZ
2965 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme ),INTENT(INOUT) :: lfr_urb3d !GRZ
2966 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_ndm, jms:jme ),INTENT(INOUT) :: lfg_urb3d !GRZ
2967 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_hi, jms:jme ), INTENT(IN) :: hi_urb2d
2968 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: lp_urb2d
2969 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: lb_urb2d
2970 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: hgt_urb2d
2971 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: mh_urb2d
2972 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: stdh_urb2d
2973 REAL, OPTIONAL, DIMENSION( ims:ime, 4, jms:jme ), INTENT(IN) :: lf_urb2d
2974 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_u_bep !Implicit momemtum component X-direction
2975 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_v_bep !Implicit momemtum component Y-direction
2976 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_t_bep !Implicit component pot. temperature
2977 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_q_bep !Implicit momemtum component X-direction
2978 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_e_bep !Implicit component TKE
2979 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::b_u_bep !Explicit momentum component X-direction
2980 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::b_v_bep !Explicit momentum component Y-direction
2981 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::b_t_bep !Explicit component pot. temperature
2982 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::b_q_bep !Implicit momemtum component Y-direction
2983 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::b_e_bep !Explicit component TKE
2984 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::vl_bep !Fraction air volume in grid cell
2985 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::dlg_bep !Height above ground
2986 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::sf_bep !Fraction air at the face of grid cell
2987 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::dl_u_bep !Length scale
2989 ! Local variables for multi-layer UCM (Martilli et al. 2002)
2990 REAL, DIMENSION( its:ite, jts:jte ) :: HFX_RURAL,LH_RURAL,GRDFLX_RURAL ! ,RN_RURAL
2991 REAL, DIMENSION( its:ite, jts:jte ) :: QFX_RURAL ! ,QSFC_RURAL,UMOM_RURAL,VMOM_RURAL
2992 REAL, DIMENSION( its:ite, jts:jte ) :: ALB_RURAL,EMISS_RURAL,TSK_RURAL ! ,UST_RURAL
2993 ! REAL, DIMENSION( ims:ime, jms:jme ) :: QSFC_URB
2994 REAL, DIMENSION( its:ite, jts:jte ) :: HFX_URB,UMOM_URB,VMOM_URB
2995 REAL, DIMENSION( its:ite, jts:jte ) :: QFX_URB
2996 ! REAL, DIMENSION( ims:ime, jms:jme ) :: ALBEDO_URB,EMISS_URB,UMOM,VMOM,UST
2997 REAL, DIMENSION(its:ite,jts:jte) ::EMISS_URB
2998 REAL, DIMENSION(its:ite,jts:jte) :: RL_UP_URB
2999 REAL, DIMENSION(its:ite,jts:jte) ::RS_ABS_URB
3000 REAL, DIMENSION(its:ite,jts:jte) ::GRDFLX_URB
3001 REAL :: SIGMA_SB,RL_UP_RURAL,RL_UP_TOT,RS_ABS_TOT,UMOM,VMOM
3002 REAL :: CMR_URB, CHR_URB, CMC_URB, CHC_URB, CMGR_URB, CHGR_URB
3003 REAL :: frc_urb,lb_urb
3005 ! ----------------------------------------------------------------------
3006 ! DECLARATIONS END - urban
3007 ! ----------------------------------------------------------------------
3008 !-------------------------------------------------
3009 ! Noah-mosaic related variables are added to declaration (danli)
3010 !-------------------------------------------------
3012 INTEGER, INTENT(IN) :: sf_surface_mosaic
3013 INTEGER, INTENT(IN) :: mosaic_cat, NLCAT
3014 REAL, DIMENSION( ims:ime, NLCAT, jms:jme ), INTENT(IN) :: landusef
3015 REAL, DIMENSION( ims:ime, NLCAT, jms:jme ), INTENT(INOUT) ::landusef2
3016 INTEGER, DIMENSION( ims:ime, NLCAT, jms:jme ), INTENT(INOUT) :: mosaic_cat_index
3018 REAL, DIMENSION( ims:ime, 1:mosaic_cat, jms:jme ) , OPTIONAL, INTENT(INOUT):: &
3019 TSK_mosaic, QSFC_mosaic, CANWAT_mosaic, SNOW_mosaic,SNOWH_mosaic, SNOWC_mosaic
3020 REAL, DIMENSION( ims:ime, 1:mosaic_cat, jms:jme ) , OPTIONAL, INTENT(INOUT):: &
3021 ALBEDO_mosaic,ALBBCK_mosaic, EMISS_mosaic, EMBCK_mosaic, ZNT_mosaic, Z0_mosaic, &
3022 HFX_mosaic,QFX_mosaic, LH_mosaic, GRDFLX_mosaic,SNOTIME_mosaic
3023 REAL, DIMENSION( ims:ime, 1:num_soil_layers*mosaic_cat, jms:jme ), OPTIONAL, INTENT(INOUT):: &
3024 TSLB_mosaic,SMOIS_mosaic,SH2O_mosaic
3025 REAL, DIMENSION( ims:ime, 1:mosaic_cat, jms:jme ), OPTIONAL, INTENT(INOUT):: LAI_mosaic, RC_mosaic
3027 REAL, DIMENSION( ims:ime, jms:jme ) :: TSK_mosaic_avg, QSFC_mosaic_avg, CANWAT_mosaic_avg,SNOW_mosaic_avg,SNOWH_mosaic_avg, &
3028 SNOWC_mosaic_avg, HFX_mosaic_avg, QFX_mosaic_avg, LH_mosaic_avg, GRDFLX_mosaic_avg, &
3029 ALBEDO_mosaic_avg, ALBBCK_mosaic_avg, EMISS_mosaic_avg, EMBCK_mosaic_avg, &
3030 ZNT_mosaic_avg, Z0_mosaic_avg, LAI_mosaic_avg, RC_mosaic_avg, SNOTIME_mosaic_avg, &
3032 REAL, DIMENSION( ims:ime, 1:num_soil_layers, jms:jme ) :: &
3033 TSLB_mosaic_avg,SMOIS_mosaic_avg,SH2O_mosaic_avg
3035 REAL, DIMENSION( ims:ime, 1:mosaic_cat, jms:jme ) , OPTIONAL, INTENT(INOUT):: &
3036 TR_URB2D_mosaic, TB_URB2D_mosaic, TG_URB2D_mosaic, TC_URB2D_mosaic,QC_URB2D_mosaic, UC_URB2D_mosaic, &
3037 SH_URB2D_mosaic,LH_URB2D_mosaic,G_URB2D_mosaic,RN_URB2D_mosaic,TS_URB2D_mosaic, TS_RUL2D_mosaic
3039 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_roof_layers*mosaic_cat, jms:jme ), INTENT(INOUT) :: TRL_URB3D_mosaic
3040 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_wall_layers*mosaic_cat, jms:jme ), INTENT(INOUT) :: TBL_URB3D_mosaic
3041 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_road_layers*mosaic_cat, jms:jme ), INTENT(INOUT) :: TGL_URB3D_mosaic
3043 INTEGER, DIMENSION( ims:ime, jms:jme ) :: IVGTYP_dominant
3044 INTEGER :: mosaic_i, URBAN_METHOD, zo_avg_option
3046 LOGICAL :: IPRINT_mosaic, Noah_call
3047 !-------------------------------------------------
3048 ! Noah-mosaic related variables declaration end (danli)
3049 !-------------------------------------------------
3051 REAL, PARAMETER :: CAPA=R_D/CP
3052 REAL :: APELM,APES,SFCTH2,PSFC
3053 real, intent(in) :: xice_threshold
3054 character(len=80) :: message_text
3056 ! FASDAS: it doesn't work for mosaic, but we need the variables to call sflx
3058 REAL, DIMENSION( ims:ime, jms:jme ) , &
3059 INTENT(INOUT) :: SDA_HFX, SDA_QFX, HFX_BOTH, QFX_BOTH, QNORM
3060 INTEGER, INTENT(IN ) :: fasdas
3061 REAL :: XSDA_HFX, XSDA_QFX, XQNORM
3062 REAL :: HFX_PHY, QFX_PHY
3064 REAL :: HCPCT_FASDAS
3065 REAL,OPTIONAL,DIMENSION( ims:ime, jms:jme ) :: IRR_CHAN
3066 REAL :: IRRIGATION_CHANNEL
3067 IRRIGATION_CHANNEL=0.0
3068 HFX_PHY = 0.0 ! initialize
3079 FDTW=DT/(XLV*RHOWATER)
3082 IPRINT_mosaic=.false.
3088 NSOIL=num_soil_layers
3094 JLOOP : DO J=jts,jte
3096 IF(ITIMESTEP.EQ.1)THEN
3098 !*** initialize soil conditions for IHOP 31 May case
3099 ! IF((XLAND(I,J)-1.5) < 0.)THEN
3100 ! if (I==108.and.j==85) then
3102 ! SMOIS(I,NS,J)=0.10
3108 !*** SET ZERO-VALUE FOR SOME OUTPUT DIAGNOSTIC ARRAYS
3109 IF((XLAND(I,J)-1.5).GE.0.)THEN
3110 ! check sea-ice point
3112 IF( XICE(I,J).GE. XICE_THRESHOLD .and. IPRINT ) PRINT*, ' sea-ice at water point, I=',I,'J=',J
3114 !*** Open Water Case
3119 TSLB(I,NS,J)=273.16 !STEMP
3123 IF ( XICE(I,J) .GE. XICE_THRESHOLD ) THEN
3135 ENDIF ! end of initialization over ocean
3137 !-----------------------------------------------------------------------
3138 ILOOP : DO I=its,ite
3140 ! initializing local variables
3147 IF (((XLAND(I,J)-1.5).LT.0.) .AND. (XICE(I,J) < XICE_THRESHOLD) ) THEN
3149 IVGTYP_dominant(I,J)=IVGTYP(I,J) ! save this
3151 ! INITIALIZE THE AREA-AVERAGED FLUXES
3153 TSK_mosaic_avg(i,j)= 0.0 ! from 3D to 2D
3154 QSFC_mosaic_avg(i,j)= 0.0
3155 CANWAT_mosaic_avg(i,j)= 0.0
3156 SNOW_mosaic_avg(i,j)= 0.0
3157 SNOWH_mosaic_avg(i,j)= 0.0
3158 SNOWC_mosaic_avg(i,j)= 0.0
3162 TSLB_mosaic_avg(i,NS,j)=0.0
3163 SMOIS_mosaic_avg(i,NS,j)=0.0
3164 SH2O_mosaic_avg(i,NS,j)=0.0
3168 HFX_mosaic_avg(i,j)= 0.0
3169 QFX_mosaic_avg(i,j)= 0.0
3170 LH_mosaic_avg(i,j)= 0.0
3171 GRDFLX_mosaic_avg(i,j)= 0.0
3172 ALBEDO_mosaic_avg(i,j)=0.0
3173 ALBBCK_mosaic_avg(i,j)=0.0
3174 EMISS_mosaic_avg(i,j)=0.0
3175 EMBCK_mosaic_avg(i,j)=0.0
3176 ZNT_mosaic_avg(i,j)=0.0
3177 Z0_mosaic_avg(i,j)=0.0
3178 LAI_mosaic_avg(i,j)=0.0
3179 RC_mosaic_avg(i,j)=0.0
3180 FAREA_mosaic_avg(i,j)=0.0
3182 ! add a new loop for the mosaic_cat
3184 DO mosaic_i = mosaic_cat, 1, -1
3186 ! if (mosaic_cat_index(I,mosaic_i,J) .EQ. 16 ) then
3187 ! PRINT*, 'you still have water tiles at','i=',i,'j=',j, 'mosaic_i',mosaic_i
3188 ! PRINT*, 'xland',xland(i,j),'xice',xice(i,j)
3191 IVGTYP(I,J)=mosaic_cat_index(I,mosaic_i,J) ! replace it with the mosaic one
3192 TSK(I,J)=TSK_mosaic(I,mosaic_i,J) ! from 3D to 2D
3193 QSFC(i,j)=QSFC_mosaic(I,mosaic_i,J)
3194 CANWAT(i,j)=CANWAT_mosaic(i,mosaic_i,j)
3195 SNOW(i,j)=SNOW_mosaic(i,mosaic_i,j)
3196 SNOWH(i,j)=SNOWH_mosaic(i,mosaic_i,j)
3197 SNOWC(i,j)=SNOWC_mosaic(i,mosaic_i,j)
3199 ALBEDO(i,j) = ALBEDO_mosaic(i,mosaic_i,j)
3200 ALBBCK(i,j)= ALBBCK_mosaic(i,mosaic_i,j)
3201 EMISS(i,j)= EMISS_mosaic(i,mosaic_i,j)
3202 EMBCK(i,j)= EMBCK_mosaic(i,mosaic_i,j)
3203 ZNT(i,j)= ZNT_mosaic(i,mosaic_i,j)
3204 Z0(i,j)= Z0_mosaic(i,mosaic_i,j)
3206 SNOTIME(i,j)= SNOTIME_mosaic(i,mosaic_i,j)
3210 TSLB(i,NS,j)=TSLB_mosaic(i,NSOIL*(mosaic_i-1)+NS,j)
3211 SMOIS(i,NS,j)=SMOIS_mosaic(i,NSOIL*(mosaic_i-1)+NS,j)
3212 SH2O(i,NS,j)=SH2O_mosaic(i,NSOIL*(mosaic_i-1)+NS,j)
3216 IF(IPRINT_mosaic) THEN
3218 print*, 'BEFORE SFLX, in Noahdrv.F'
3219 print*, 'mosaic_cat', mosaic_cat, 'IVGTYP',IVGTYP(i,j), 'TSK',TSK(i,j),'HFX',HFX(i,j), 'QSFC', QSFC(i,j), &
3220 'CANWAT', CANWAT(i,j), 'SNOW',SNOW(i,j), 'ALBEDO',ALBEDO(i,j), 'TSLB',TSLB(i,1,j),'CHS',CHS(i,j),'ZNT',ZNT(i,j)
3224 !-----------------------------------------------------------------------
3225 ! insert the NOAH model here for the non-urban one and the urban one DANLI
3226 !-----------------------------------------------------------------------
3230 ! pressure in middle of lowest layer
3231 SFCPRS=(P8W3D(I,KTS+1,j)+P8W3D(i,KTS,j))*0.5
3232 ! convert from mixing ratio to specific humidity
3233 Q2K=QV3D(i,1,j)/(1.0+QV3D(i,1,j))
3236 Q2SAT=QGH(I,J)/(1.0+QGH(I,J)) ! Q2SAT is sp humidity
3237 ! add check on myj=.true.
3238 ! IF((Q2K.GE.Q2SAT*TRESH).AND.Q2K.LT.QZ0(I,J))THEN
3239 IF((myj).AND.(Q2K.GE.Q2SAT*TRESH).AND.Q2K.LT.QZ0(I,J))THEN
3248 ZLVL=0.5*DZ8W(i,1,j)
3250 ! TH2=SFCTMP+(0.0097545*ZLVL)
3251 ! calculate SFCTH2 via Exner function vs lapse-rate (above)
3252 APES=(1.E5/PSFC)**CAPA
3253 APELM=(1.E5/SFCPRS)**CAPA
3258 LWDN=GLW(I,J)*EMISSI
3259 ! SOLDN is total incoming solar
3261 ! GSW is net downward solar
3263 ! use mid-day albedo to determine net downward solar (no solar zenith angle correction)
3264 SOLNET=SOLDN*(1.-ALBEDO(I,J))
3266 IF(PRESENT(IRR_CHAN)) THEN
3267 IF(IRR_CHAN(i,j).NE.0) THEN
3268 IRRIGATION_CHANNEL=IRR_CHAN(i,j)/DT
3270 IRRIGATION_CHANNEL=0.
3275 SHDFAC=VEGFRA(I,J)/100.
3278 SHMIN=SHDMIN(I,J)/100. !NEW
3279 SHMAX=SHDMAX(I,J)/100. !NEW
3280 ! convert snow water equivalent from mm to meter
3281 SNEQV=SNOW(I,J)*0.001
3282 ! snow depth in meters
3286 ! if "SR" present, set frac of frozen precip ("FFROZP") = snow-ratio ("SR", range:0-1)
3287 ! SR from e.g. Ferrier microphysics
3288 ! otherwise define from 1st atmos level temperature
3292 IF (SFCTMP <= 273.15) THEN
3299 IF((XLAND(I,J)-1.5).GE.0.)THEN ! begining of land/sea if block
3301 TSK_RURAL(I,J)=TSK(I,J)
3302 HFX_RURAL(I,J)=HFX(I,J)
3303 QFX_RURAL(I,J)=QFX(I,J)
3304 LH_RURAL(I,J)=LH(I,J)
3305 EMISS_RURAL(I,J)=EMISS(I,J)
3306 GRDFLX_RURAL(I,J)=GRDFLX(I,J)
3308 ! Land or sea-ice case
3310 IF (XICE(I,J) >= XICE_THRESHOLD) THEN
3313 ELSE IF ( VEGTYP == ISICE ) THEN
3317 ! Neither sea ice or land ice.
3320 DQSDT2=Q2SAT*A23M4/(SFCTMP-A4)**2
3322 IF(SNOW(I,J).GT.0.0)THEN
3323 ! snow on surface (use ice saturation properties)
3325 E2SAT=611.2*EXP(6174.*(1./273.15 - 1./SFCTSNO))
3326 Q2SATI=0.622*E2SAT/(SFCPRS-E2SAT)
3327 Q2SATI=Q2SATI/(1.0+Q2SATI) ! spec. hum.
3328 IF (T1 .GT. 273.14) THEN
3329 ! warm ground temps, weight the saturation between ice and water according to SNOWC
3330 Q2SAT=Q2SAT*(1.-SNOWC(I,J)) + Q2SATI*SNOWC(I,J)
3331 DQSDT2=DQSDT2*(1.-SNOWC(I,J)) + Q2SATI*6174./(SFCTSNO**2)*SNOWC(I,J)
3333 ! cold ground temps, use ice saturation only
3335 DQSDT2=Q2SATI*6174./(SFCTSNO**2)
3337 ! for snow cover fraction at 0 C, ground temp will not change, so DQSDT2 effectively zero
3338 IF(T1 .GT. 273. .AND. SNOWC(I,J) .GT. 0.)DQSDT2=DQSDT2*(1.-SNOWC(I,J))
3341 ! Land-ice or land points use the usual deep-soil temperature.
3344 IF(VEGTYP.EQ.25) SHDFAC=0.0000
3345 IF(VEGTYP.EQ.26) SHDFAC=0.0000
3346 IF(VEGTYP.EQ.27) SHDFAC=0.0000
3347 IF(SOILTYP.EQ.14.AND.XICE(I,J).EQ.0.)THEN
3349 IF(IPRINT)PRINT*,' SOIL TYPE FOUND TO BE WATER AT A LAND-POINT'
3350 IF(IPRINT)PRINT*,i,j,'RESET SOIL in surfce.F'
3354 SNOALB1 = SNOALB(I,J)
3355 ! converts canwat in mm to CMC in meters
3356 CMC=CANWAT(I,J)/1000.
3358 !-------------------------------------------
3359 !*** convert snow depth from mm to meter
3362 ! SNOALB=ALBMAX(I,J)*0.01
3364 ! SNOALB=MAXALB(IVGTPK)*0.01
3372 SNOTIME1 = SNOTIME(I,J)
3374 !FEI: temporaray arrays above need to be changed later by using SI
3377 SMC(NS)=SMOIS(I,NS,J)
3378 STC(NS)=TSLB(I,NS,J) !STEMP
3379 SWC(NS)=SH2O(I,NS,J)
3382 if ( (SNEQV.ne.0..AND.SNOWHK.eq.0.).or.(SNOWHK.le.SNEQV) )THEN
3387 !Fei: urban. for urban surface, if calling UCM, redefine the natural surface in cities as
3388 ! the "NATURAL" category in the VEGPARM.TBL
3390 ! IF(SF_URBAN_PHYSICS == 1.OR. SF_URBAN_PHYSICS==2.OR.SF_URBAN_PHYSICS==3 ) THEN
3392 ! IF( IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == LOW_DENSITY_RESIDENTIAL .or. &
3393 ! IVGTYP(I,J) == HIGH_DENSITY_RESIDENTIAL .or. IVGTYP(I,J) == HIGH_INTENSITY_INDUSTRIAL) THEN
3395 ! SHDFAC = SHDTBL(NATURAL)
3396 ! ALBEDOK =0.2 ! 0.2
3398 ! EMISSI = 0.98 !for VEGTYP=5
3399 ! IF ( FRC_URB2D(I,J) < 0.99 ) THEN
3400 ! if(sf_urban_physics.eq.1)then
3401 ! T1= ( TSK(I,J) -FRC_URB2D(I,J) * TS_URB2D (I,J) )/ (1-FRC_URB2D(I,J))
3402 ! elseif((sf_urban_physics.eq.2).OR.(sf_urban_physics.eq.3))then
3403 ! r1= (tsk(i,j)**4.)
3404 ! r2= frc_urb2d(i,j)*(ts_urb2d(i,j)**4.)
3405 ! r3= (1.-frc_urb2d(i,j))
3406 ! t1= ((r1-r2)/r3)**.25
3413 ! IF( IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == LOW_DENSITY_RESIDENTIAL .or. &
3414 ! IVGTYP(I,J) == HIGH_DENSITY_RESIDENTIAL .or. IVGTYP(I,J) == HIGH_INTENSITY_INDUSTRIAL) THEN
3421 If ( SF_URBAN_PHYSICS == 0 ) THEN ! ONLY NOAH
3422 IF( IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == LCZ_1 .or. IVGTYP(I,J) == LCZ_2 .or. &
3423 IVGTYP(I,J) == LCZ_3 .or. IVGTYP(I,J) == LCZ_4 .or. IVGTYP(I,J) == LCZ_5 .or. &
3424 IVGTYP(I,J) == LCZ_6 .or. IVGTYP(I,J) == LCZ_7 .or. IVGTYP(I,J) == LCZ_8 .or. &
3425 IVGTYP(I,J) == LCZ_9 .or. IVGTYP(I,J) == LCZ_10 .or. IVGTYP(I,J) == LCZ_11 ) THEN
3433 IF(SF_URBAN_PHYSICS == 1) THEN
3435 IF( IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == LCZ_1 .or. IVGTYP(I,J) == LCZ_2 .or. &
3436 IVGTYP(I,J) == LCZ_3 .or. IVGTYP(I,J) == LCZ_4 .or. IVGTYP(I,J) == LCZ_5 .or. &
3437 IVGTYP(I,J) == LCZ_6 .or. IVGTYP(I,J) == LCZ_7 .or. IVGTYP(I,J) == LCZ_8 .or. &
3438 IVGTYP(I,J) == LCZ_9 .or. IVGTYP(I,J) == LCZ_10 .or. IVGTYP(I,J) == LCZ_11 ) THEN
3441 SHDFAC = SHDTBL(NATURAL)
3444 EMISSI = 0.98 ! for VEGTYP=5
3445 LWDN = GLW(I,J) * EMISSI
3446 SOLNET = SOLDN * (1.0 - ALBEDOK)
3448 T1= TS_RUL2D_mosaic(I,mosaic_i,J)
3454 !===Yang, 2014/10/08, hydrological processes for urban vegetation in single layer UCM===
3458 IRIOPTION=IRI_SCHEME
3460 tloc=mod(int(OMG/3.14159*180./15.+12.+0.5 ),24)
3461 if (tloc.lt.0) tloc=tloc+24
3462 if (tloc==0) tloc=24
3463 CALL cal_mon_day(julian,julyr,jmonth,jday)
3464 IF(SF_URBAN_PHYSICS == 1) THEN
3465 IF( IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == LCZ_1 .or. IVGTYP(I,J) == LCZ_2 .or. &
3466 IVGTYP(I,J) == LCZ_3 .or. IVGTYP(I,J) == LCZ_4 .or. IVGTYP(I,J) == LCZ_5 .or. &
3467 IVGTYP(I,J) == LCZ_6 .or. IVGTYP(I,J) == LCZ_7 .or. IVGTYP(I,J) == LCZ_8 .or. &
3468 IVGTYP(I,J) == LCZ_9 .or. IVGTYP(I,J) == LCZ_10 .or. IVGTYP(I,J) == LCZ_11 ) THEN
3470 AOASIS = oasis ! urban oasis effect
3471 IF (IRIOPTION ==1) THEN
3472 IF (tloc==21 .or. tloc==22) THEN !irrigation on vegetaion in urban area, MAY-SEP, 9-10pm
3473 IF (jmonth==5 .or. jmonth==6 .or. jmonth==7 .or. jmonth==8 .or. jmonth==9) THEN
3474 IF (SMC(USOIL) .LT. SMCREF) SMC(USOIL)= REFSMC(ISLTYP(I,J))
3475 IF (SMC(DSOIL) .LT. SMCREF) SMC(DSOIL)= REFSMC(ISLTYP(I,J))
3482 IF(SF_URBAN_PHYSICS == 2 .or. SF_URBAN_PHYSICS == 3) THEN
3483 IF(AOASIS > 1.0) THEN
3484 CALL wrf_error_fatal('Urban oasis option is for SF_URBAN_PHYSICS == 1 only')
3486 IF(IRIOPTION == 1) THEN
3487 CALL wrf_error_fatal('Urban irrigation option is for SF_URBAN_PHYSICS == 1 only')
3491 IF( SF_URBAN_PHYSICS==2.OR.SF_URBAN_PHYSICS==3 ) THEN
3492 ! print*, 'MOSAIC is not designed to work with SF_URBAN_PHYSICS=2 or SF_URBAN_PHYSICS=3'
3499 print*, 'BEFORE SFLX, in Noahlsm_driver'
3500 print*, 'ICE', ICE, 'DT',DT, 'ZLVL',ZLVL, 'NSOIL', NSOIL, &
3501 'SLDPTH', SLDPTH, 'LOCAL',LOCAL, 'LUTYPE',&
3502 LUTYPE, 'SLTYPE',SLTYPE, 'LWDN',LWDN, 'SOLDN',SOLDN, &
3503 'SFCPRS',SFCPRS, 'PRCP',PRCP,'SFCTMP',SFCTMP,'Q2K',Q2K, &
3504 'TH2',TH2,'Q2SAT',Q2SAT,'DQSDT2',DQSDT2,'VEGTYP', VEGTYP,&
3505 'SOILTYP',SOILTYP, 'SLOPETYP',SLOPETYP, 'SHDFAC',SHDFAC,&
3506 'SHMIN',SHMIN, 'ALBBRD',ALBBRD,'SNOALB1',SNOALB1,'TBOT',&
3507 TBOT, 'Z0BRD',Z0BRD, 'Z0K',Z0K, 'CMC',CMC, 'T1',T1,'STC',&
3508 STC, 'SMC',SMC, 'SWC',SWC,'SNOWHK',SNOWHK,'SNEQV',SNEQV,&
3509 'ALBEDOK',ALBEDOK,'CHK',CHK,'ETA',ETA,'SHEAT',SHEAT, &
3510 'ETA_KINEMATIC',ETA_KINEMATIC, 'FDOWN',FDOWN,'EC',EC, &
3511 'EDIR',EDIR,'ET',ET,'ETT',ETT,'ESNOW',ESNOW,'DRIP',DRIP,&
3512 'DEW',DEW,'BETA',BETA,'ETP',ETP,'SSOIL',SSOIL,'FLX1',FLX1,&
3513 'FLX2',FLX2,'FLX3',FLX3,'SNOMLT',SNOMLT,'SNCOVR',SNCOVR,&
3514 'RUNOFF1',RUNOFF1,'RUNOFF2',RUNOFF2,'RUNOFF3',RUNOFF3, &
3515 'RC',RC, 'PC',PC,'RSMIN',RSMIN,'XLAI',XLAI,'RCS',RCS, &
3516 'RCT',RCT,'RCQ',RCQ,'RCSOIL',RCSOIL,'SOILW',SOILW, &
3517 'SOILM',SOILM,'Q1',Q1,'SMCWLT',SMCWLT,'SMCDRY',SMCDRY,&
3518 'SMCREF',SMCREF,'SMCMAX',SMCMAX,'NROOT',NROOT
3523 IF (SHDFAC > 0.0 .AND. LAI(I,J) <= 0.0) LAI(I,J) = 0.01
3527 IF ( ICE == 1 ) THEN
3538 ELSEIF (ICE == 0) THEN
3542 CALL SFLX (I,J,FFROZP, ISURBAN, DT,ZLVL,NSOIL,SLDPTH, & !C
3544 LUTYPE, SLTYPE, & !CL
3545 LWDN,SOLDN,SOLNET,SFCPRS,PRCP,SFCTMP,Q2K,DUMMY, & !F
3546 DUMMY,DUMMY, DUMMY, & !F PRCPRAIN not used
3547 TH2,Q2SAT,DQSDT2, & !I
3548 VEGTYP,SOILTYP,SLOPETYP,SHDFAC,SHMIN,SHMAX, & !I
3549 ALBBRD, SNOALB1,TBOT, Z0BRD, Z0K, EMISSI, EMBRD, & !S
3550 CMC,T1,STC,SMC,SWC,SNOWHK,SNEQV,ALBEDOK,CHK,dummy,& !H
3551 ETA,SHEAT, ETA_KINEMATIC,FDOWN, & !O
3552 EC,EDIR,ET,ETT,ESNOW,DRIP,DEW, & !O
3553 BETA,ETP,SSOIL, & !O
3554 FLX1,FLX2,FLX3, & !O
3555 FLX4,FVB,FBUR,FGSN,UA_PHYS, & !UA
3557 RUNOFF1,RUNOFF2,RUNOFF3, & !O
3558 RC,PC,RSMIN,XLAI,RCS,RCT,RCQ,RCSOIL, & !O
3559 SOILW,SOILM,Q1,SMAV, & !D
3560 RDLAI2D,USEMONALB, &
3563 SMCWLT,SMCDRY,SMCREF,SMCMAX,NROOT, &
3565 sfcheadrt_hydro, & !I
3566 INFXSRT_hydro,ETPND1_hydro & !O
3567 ,OPT_THCND,AOASIS & !O
3568 ,XSDA_QFX, HFX_PHY, QFX_PHY, XQNORM, fasdas, HCPCT_FASDAS & ! fasdas vars
3569 ,IRRIGATION_CHANNEL )
3572 soldrain(i,j) = RUNOFF2*DT*1000.0
3573 sfcheadrt(i,j) = sfcheadrt_hydro
3574 infxsrt(i,j) = INFXSRT_hydro
3575 etpnd1 = etpnd1_hydro
3577 ELSEIF (ICE == -1) THEN
3580 ! Set values that the LSM is expected to update,
3581 ! but don't get updated for glacial points.
3583 SOILM = 0.0 !BSINGH(PNNL)- SOILM is undefined for this case, it is used for diagnostics so setting it to zero
3584 XLAI = 0.01 ! KWM Should this be Zero over land ice? Does this value matter?
3592 CALL SFLX_GLACIAL(I,J,ISICE,FFROZP,DT,ZLVL,NSOIL,SLDPTH, & !C
3593 & LWDN,SOLNET,SFCPRS,PRCP,SFCTMP,Q2K, & !F
3594 & TH2,Q2SAT,DQSDT2, & !I
3595 & ALBBRD, SNOALB1,TBOT, Z0BRD, Z0K, EMISSI, EMBRD, & !S
3596 & T1,STC(1:NSOIL),SNOWHK,SNEQV,ALBEDOK,CHK, & !H
3597 & ETA,SHEAT,ETA_KINEMATIC,FDOWN, & !O
3600 & FLX1,FLX2,FLX3, & !O
3601 & SNOMLT,SNCOVR, & !O
3612 print*, 'AFTER SFLX, in Noahlsm_driver'
3613 print*, 'ICE', ICE, 'DT',DT, 'ZLVL',ZLVL, 'NSOIL', NSOIL, &
3614 'SLDPTH', SLDPTH, 'LOCAL',LOCAL, 'LUTYPE',&
3615 LUTYPE, 'SLTYPE',SLTYPE, 'LWDN',LWDN, 'SOLDN',SOLDN, &
3616 'SFCPRS',SFCPRS, 'PRCP',PRCP,'SFCTMP',SFCTMP,'Q2K',Q2K, &
3617 'TH2',TH2,'Q2SAT',Q2SAT,'DQSDT2',DQSDT2,'VEGTYP', VEGTYP,&
3618 'SOILTYP',SOILTYP, 'SLOPETYP',SLOPETYP, 'SHDFAC',SHDFAC,&
3619 'SHDMIN',SHMIN, 'ALBBRD',ALBBRD,'SNOALB',SNOALB1,'TBOT',&
3620 TBOT, 'Z0BRD',Z0BRD, 'Z0K',Z0K, 'CMC',CMC, 'T1',T1,'STC',&
3621 STC, 'SMC',SMC, 'SWc',SWC,'SNOWHK',SNOWHK,'SNEQV',SNEQV,&
3622 'ALBEDOK',ALBEDOK,'CHK',CHK,'ETA',ETA,'SHEAT',SHEAT, &
3623 'ETA_KINEMATIC',ETA_KINEMATIC, 'FDOWN',FDOWN,'EC',EC, &
3624 'EDIR',EDIR,'ET',ET,'ETT',ETT,'ESNOW',ESNOW,'DRIP',DRIP,&
3625 'DEW',DEW,'BETA',BETA,'ETP',ETP,'SSOIL',SSOIL,'FLX1',FLX1,&
3626 'FLX2',FLX2,'FLX3',FLX3,'SNOMLT',SNOMLT,'SNCOVR',SNCOVR,&
3627 'RUNOFF1',RUNOFF1,'RUNOFF2',RUNOFF2,'RUNOFF3',RUNOFF3, &
3628 'RC',RC, 'PC',PC,'RSMIN',RSMIN,'XLAI',XLAI,'RCS',RCS, &
3629 'RCT',RCT,'RCQ',RCQ,'RCSOIL',RCSOIL,'SOILW',SOILW, &
3630 'SOILM',SOILM,'Q1',Q1,'SMCWLT',SMCWLT,'SMCDRY',SMCDRY,&
3631 'SMCREF',SMCREF,'SMCMAX',SMCMAX,'NROOT',NROOT
3635 !*** UPDATE STATE VARIABLES
3636 CANWAT(I,J)=CMC*1000.
3637 SNOW(I,J)=SNEQV*1000.
3638 ! SNOWH(I,J)=SNOWHK*1000.
3639 SNOWH(I,J)=SNOWHK ! SNOWHK in meters
3641 ALB_RURAL(I,J)=ALBEDOK
3645 EMISS_RURAL(I,J) = EMISSI
3646 ! Noah: activate time-varying roughness length (V3.3 Feb 2011)
3651 HFX_RURAL(I,J)=SHEAT
3652 ! MEk Jul07 add potential evap accum
3653 POTEVP(I,J)=POTEVP(I,J)+ETP*FDTW
3654 QFX(I,J)=ETA_KINEMATIC
3655 QFX_RURAL(I,J)=ETA_KINEMATIC
3659 ! QFX(I,J) = QFX(I,J) + ETPND1
3660 ! ETA = ETA + ETPND1/2.501E6*dt
3661 !end added by Wei Yu
3667 GRDFLX_RURAL(I,J)=SSOIL
3670 SNOTIME(I,J) = SNOTIME1
3671 ! prevent diagnostic ground q (q1) from being greater than qsat(tsk)
3672 ! as happens over snow cover where the cqs2 value also becomes irrelevant
3673 ! by setting cqs2=chs in this situation the 2m q should become just qv(k=1)
3674 IF (Q1 .GT. QSFC(I,J)) THEN
3675 CQS2(I,J) = CHS(I,J)
3678 ! Convert QSFC back to mixing ratio
3679 QSFC(I,J)= Q1/(1.0-Q1)
3681 ! QSFC_RURAL(I,J)= Q1/(1.0-Q1)
3682 ! Calculate momentum flux from rural surface for use with multi-layer UCM (Martilli et al. 2002)
3685 SMOIS(I,NS,J)=SMC(NS)
3686 TSLB(I,NS,J)=STC(NS) ! STEMP
3687 SH2O(I,NS,J)=SWC(NS)
3697 ! Residual of surface energy balance equation terms
3701 noahres(i,j) = ( solnet + lwdn ) - sheat + ssoil - eta &
3702 - ( emissi * STBOLT * (t1**4) ) - flx1 - flx2 - flx3 - flx4
3705 noahres(i,j) = ( solnet + lwdn ) - sheat + ssoil - eta &
3706 - ( emissi * STBOLT * (t1**4) ) - flx1 - flx2 - flx3
3709 ENDIF !ENDIF FOR Noah_call
3711 IF (SF_URBAN_PHYSICS == 1 ) THEN ! Beginning of UCM CALL if block
3712 !--------------------------------------
3713 ! URBAN CANOPY MODEL START - urban
3714 !--------------------------------------
3715 ! Input variables lsm --> urban
3716 IF( IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == LCZ_1 .or. IVGTYP(I,J) == LCZ_2 .or. &
3717 IVGTYP(I,J) == LCZ_3 .or. IVGTYP(I,J) == LCZ_4 .or. IVGTYP(I,J) == LCZ_5 .or. &
3718 IVGTYP(I,J) == LCZ_6 .or. IVGTYP(I,J) == LCZ_7 .or. IVGTYP(I,J) == LCZ_8 .or. &
3719 IVGTYP(I,J) == LCZ_9 .or. IVGTYP(I,J) == LCZ_10 .or. IVGTYP(I,J) == LCZ_11 ) THEN
3722 ! UTYPE_URB = UTYPE_URB2D(I,J) !urban type (low, high or industrial)
3723 ! this need to be changed in the mosaic danli
3724 IF (use_wudapt_lcz == 1) THEN
3725 IF(IVGTYP(I,J)==ISURBAN) UTYPE_URB=5
3726 IF(IVGTYP(I,J)==LCZ_1) UTYPE_URB=1
3727 IF(IVGTYP(I,J)==LCZ_2) UTYPE_URB=2
3728 IF(IVGTYP(I,J)==LCZ_3) UTYPE_URB=3
3729 IF(IVGTYP(I,J)==LCZ_4) UTYPE_URB=4
3730 IF(IVGTYP(I,J)==LCZ_5) UTYPE_URB=5
3731 IF(IVGTYP(I,J)==LCZ_6) UTYPE_URB=6
3732 IF(IVGTYP(I,J)==LCZ_7) UTYPE_URB=7
3733 IF(IVGTYP(I,J)==LCZ_8) UTYPE_URB=8
3734 IF(IVGTYP(I,J)==LCZ_9) UTYPE_URB=9
3735 IF(IVGTYP(I,J)==LCZ_10) UTYPE_URB=10
3736 IF(IVGTYP(I,J)==LCZ_11) UTYPE_URB=11
3739 IF(UTYPE_URB==1) FRC_URB2D(I,J)=1.
3740 IF(UTYPE_URB==2) FRC_URB2D(I,J)=0.99
3741 IF(UTYPE_URB==3) FRC_URB2D(I,J)=1.00
3742 IF(UTYPE_URB==4) FRC_URB2D(I,J)=0.65
3743 IF(UTYPE_URB==5) FRC_URB2D(I,J)=0.7
3744 IF(UTYPE_URB==6) FRC_URB2D(I,J)=0.65
3745 IF(UTYPE_URB==7) FRC_URB2D(I,J)=0.3
3746 IF(UTYPE_URB==8) FRC_URB2D(I,J)=0.85
3747 IF(UTYPE_URB==9) FRC_URB2D(I,J)=0.3
3748 IF(UTYPE_URB==10) FRC_URB2D(I,J)=0.55
3749 IF(UTYPE_URB==11) FRC_URB2D(I,J)=1.
3751 IF(IVGTYP(I,J)==ISURBAN) UTYPE_URB=2
3752 IF(IVGTYP(I,J)==LCZ_1) UTYPE_URB=1 ! LOW_DENSITY_RESIDENTIAL
3753 IF(IVGTYP(I,J)==LCZ_2) UTYPE_URB=2 ! HIGH_DENSITY_RESIDENTIAL
3754 IF(IVGTYP(I,J)==LCZ_3) UTYPE_URB=3 ! HIGH_INTENSITY_INDUSTRIAL
3756 IF(UTYPE_URB==1) FRC_URB2D(I,J)=0.5
3757 IF(UTYPE_URB==2) FRC_URB2D(I,J)=0.9
3758 IF(UTYPE_URB==3) FRC_URB2D(I,J)=0.95
3761 TA_URB = SFCTMP ! [K]
3762 QA_URB = Q2K ! [kg/kg]
3763 UA_URB = SQRT(U_PHY(I,1,J)**2.+V_PHY(I,1,J)**2.)
3764 U1_URB = U_PHY(I,1,J)
3765 V1_URB = V_PHY(I,1,J)
3766 IF(UA_URB < 1.) UA_URB=1. ! [m/s]
3767 SSG_URB = SOLDN ! [W/m/m]
3768 SSGD_URB = 0.8*SOLDN ! [W/m/m]
3769 SSGQ_URB = SSG_URB-SSGD_URB ! [W/m/m]
3770 LLG_URB = GLW(I,J) ! [W/m/m]
3771 RAIN_URB = RAINBL(I,J) / DT * 3600.0 ! [mm/hr]
3772 RHOO_URB = SFCPRS / (287.04 * SFCTMP * (1.0+ 0.61 * Q2K)) ![kg/m/m/m]
3774 DELT_URB = DT ! [sec]
3775 XLAT_URB = XLAT_URB2D(I,J) ! [deg]
3776 COSZ_URB = COSZ_URB2D(I,J) !
3777 OMG_URB = OMG_URB2D(I,J) !
3780 LSOLAR_URB = .FALSE.
3784 TR_URB2D(I,J)=TR_URB2D_mosaic(I,mosaic_i,J) ! replace it with the mosaic one
3785 TB_URB2D(I,J)=TB_URB2D_mosaic(I,mosaic_i,J) ! replace it with the mosaic one
3786 TG_URB2D(I,J)=TG_URB2D_mosaic(I,mosaic_i,J) ! replace it with the mosaic one
3787 TC_URB2D(I,J)=TC_URB2D_mosaic(I,mosaic_i,J) ! replace it with the mosaic one
3788 QC_URB2D(I,J)=QC_URB2D_mosaic(I,mosaic_i,J) ! replace it with the mosaic one
3789 UC_URB2D(I,J)=UC_URB2D_mosaic(I,mosaic_i,J) ! replace it with the mosaic one
3790 TS_URB2D(I,J)=TS_URB2D_mosaic(I,mosaic_i,J) ! replace it with the mosaic one
3792 DO K = 1,num_roof_layers
3793 TRL_URB3D(I,K,J) = TRL_URB3D_mosaic(I,K+(mosaic_i-1)*num_roof_layers,J)
3795 DO K = 1,num_wall_layers
3796 TBL_URB3D(I,K,J) = TBL_URB3D_mosaic(I,K+(mosaic_i-1)*num_roof_layers,J)
3798 DO K = 1,num_road_layers
3799 TGL_URB3D(I,K,J) = TGL_URB3D_mosaic(I,K+(mosaic_i-1)*num_roof_layers,J)
3804 TR_URB = TR_URB2D(I,J)
3805 TB_URB = TB_URB2D(I,J)
3806 TG_URB = TG_URB2D(I,J)
3807 TC_URB = TC_URB2D(I,J)
3808 QC_URB = QC_URB2D(I,J)
3809 UC_URB = UC_URB2D(I,J)
3811 DO K = 1,num_roof_layers
3812 TRL_URB(K) = TRL_URB3D(I,K,J)
3813 SMR_URB(K) = SMR_URB3D(I,K,J)
3814 TGRL_URB(K)= TGRL_URB3D(I,K,J)
3816 DO K = 1,num_wall_layers
3817 TBL_URB(K) = TBL_URB3D(I,K,J)
3819 DO K = 1,num_road_layers
3820 TGL_URB(K) = TGL_URB3D(I,K,J)
3823 TGR_URB = TGR_URB2D(I,J)
3824 CMCR_URB = CMCR_URB2D(I,J)
3825 FLXHUMR_URB = FLXHUMR_URB2D(I,J)
3826 FLXHUMB_URB = FLXHUMB_URB2D(I,J)
3827 FLXHUMG_URB = FLXHUMG_URB2D(I,J)
3828 DRELR_URB = DRELR_URB2D(I,J)
3829 DRELB_URB = DRELB_URB2D(I,J)
3830 DRELG_URB = DRELG_URB2D(I,J)
3832 XXXR_URB = XXXR_URB2D(I,J)
3833 XXXB_URB = XXXB_URB2D(I,J)
3834 XXXG_URB = XXXG_URB2D(I,J)
3835 XXXC_URB = XXXC_URB2D(I,J)
3837 ! Limits to avoid dividing by small number
3838 if (CHS(I,J) < 1.0E-02) then
3841 if (CHS2(I,J) < 1.0E-02) then
3844 if (CQS2(I,J) < 1.0E-02) then
3849 CHS2_URB = CHS2(I,J)
3850 IF (PRESENT(CMR_SFCDIF)) THEN
3851 CMR_URB = CMR_SFCDIF(I,J)
3852 CHR_URB = CHR_SFCDIF(I,J)
3853 CMGR_URB = CMGR_SFCDIF(I,J)
3854 CHGR_URB = CHGR_SFCDIF(I,J)
3855 CMC_URB = CMC_SFCDIF(I,J)
3856 CHC_URB = CHC_SFCDIF(I,J)
3860 mh_urb = mh_urb2d(I,J)
3861 stdh_urb = stdh_urb2d(I,J)
3862 lp_urb = lp_urb2d(I,J)
3863 hgt_urb = hgt_urb2d(I,J)
3866 lf_urb(K)=lf_urb2d(I,K,J)
3868 frc_urb = frc_urb2d(I,J)
3869 lb_urb = lb_urb2d(I,J)
3871 if (I.eq.73.and.J.eq.125)THEN
3876 CALL cal_mon_day(julian,julyr,jmonth,jday)
3877 CALL urban(LSOLAR_URB, & ! I
3878 num_roof_layers,num_wall_layers,num_road_layers, & ! C
3880 UTYPE_URB,TA_URB,QA_URB,UA_URB,U1_URB,V1_URB,SSG_URB, & ! I
3881 SSGD_URB,SSGQ_URB,LLG_URB,RAIN_URB,RHOO_URB, & ! I
3882 ZA_URB,DECLIN_URB,COSZ_URB,OMG_URB, & ! I
3883 XLAT_URB,DELT_URB,ZNT_URB, & ! I
3884 CHS_URB, CHS2_URB, & ! I
3885 TR_URB, TB_URB, TG_URB, TC_URB, QC_URB,UC_URB, & ! H
3886 TRL_URB,TBL_URB,TGL_URB, & ! H
3887 XXXR_URB, XXXB_URB, XXXG_URB, XXXC_URB, & ! H
3888 TS_URB,QS_URB,SH_URB,LH_URB,LH_KINEMATIC_URB, & ! O
3889 SW_URB,ALB_URB,LW_URB,G_URB,RN_URB,PSIM_URB,PSIH_URB, & ! O
3891 CMR_URB, CHR_URB, CMC_URB, CHC_URB, &
3892 U10_URB, V10_URB, TH2_URB, Q2_URB, & ! O
3893 UST_URB,mh_urb, stdh_urb, lf_urb, lp_urb, & ! 0
3894 hgt_urb,frc_urb,lb_urb, check,CMCR_URB,TGR_URB, & ! H
3895 TGRL_URB,SMR_URB,CMGR_URB,CHGR_URB,jmonth, & ! H
3896 DRELR_URB,DRELB_URB, & ! H
3897 DRELG_URB,FLXHUMR_URB,FLXHUMB_URB,FLXHUMG_URB)
3902 print*, 'AFTER CALL URBAN'
3903 print*,'num_roof_layers',num_roof_layers, 'num_wall_layers', &
3905 'DZR',DZR,'DZB',DZB,'DZG',DZG,'UTYPE_URB',UTYPE_URB,'TA_URB', &
3907 'QA_URB',QA_URB,'UA_URB',UA_URB,'U1_URB',U1_URB,'V1_URB', &
3909 'SSG_URB',SSG_URB,'SSGD_URB',SSGD_URB,'SSGQ_URB',SSGQ_URB, &
3910 'LLG_URB',LLG_URB,'RAIN_URB',RAIN_URB,'RHOO_URB',RHOO_URB, &
3911 'ZA_URB',ZA_URB, 'DECLIN_URB',DECLIN_URB,'COSZ_URB',COSZ_URB,&
3912 'OMG_URB',OMG_URB,'XLAT_URB',XLAT_URB,'DELT_URB',DELT_URB, &
3913 'ZNT_URB',ZNT_URB,'TR_URB',TR_URB, 'TB_URB',TB_URB,'TG_URB',&
3914 TG_URB,'TC_URB',TC_URB,'QC_URB',QC_URB,'TRL_URB',TRL_URB, &
3915 'TBL_URB',TBL_URB,'TGL_URB',TGL_URB,'XXXR_URB',XXXR_URB, &
3916 'XXXB_URB',XXXB_URB,'XXXG_URB',XXXG_URB,'XXXC_URB',XXXC_URB,&
3917 'TS_URB',TS_URB,'QS_URB',QS_URB,'SH_URB',SH_URB,'LH_URB', &
3918 LH_URB, 'LH_KINEMATIC_URB',LH_KINEMATIC_URB,'SW_URB',SW_URB,&
3919 'ALB_URB',ALB_URB,'LW_URB',LW_URB,'G_URB',G_URB,'RN_URB', &
3920 RN_URB, 'PSIM_URB',PSIM_URB,'PSIH_URB',PSIH_URB, &
3921 'U10_URB',U10_URB,'V10_URB',V10_URB,'TH2_URB',TH2_URB, &
3922 'Q2_URB',Q2_URB,'CHS_URB',CHS_URB,'CHS2_URB',CHS2_URB
3926 TS_URB2D(I,J) = TS_URB
3928 ALBEDO(I,J) = FRC_URB2D(I,J)*ALB_URB+(1-FRC_URB2D(I,J))*ALBEDOK ![-]
3929 HFX(I,J) = FRC_URB2D(I,J)*SH_URB+(1-FRC_URB2D(I,J))*SHEAT ![W/m/m]
3930 QFX(I,J) = FRC_URB2D(I,J)*LH_KINEMATIC_URB &
3931 + (1-FRC_URB2D(I,J))*ETA_KINEMATIC ![kg/m/m/s]
3932 LH(I,J) = FRC_URB2D(I,J)*LH_URB+(1-FRC_URB2D(I,J))*ETA ![W/m/m]
3933 GRDFLX(I,J) = FRC_URB2D(I,J)*G_URB+(1-FRC_URB2D(I,J))*SSOIL ![W/m/m]
3934 TSK(I,J) = FRC_URB2D(I,J)*TS_URB+(1-FRC_URB2D(I,J))*T1 ![K]
3935 Q1 = FRC_URB2D(I,J)*QS_URB+(1-FRC_URB2D(I,J))*Q1 ![-]
3936 ! Convert QSFC back to mixing ratio
3937 QSFC(I,J)= Q1/(1.0-Q1)
3938 UST(I,J)= FRC_URB2D(I,J)*UST_URB+(1-FRC_URB2D(I,J))*UST(I,J) ![m/s]
3939 ZNT(I,J)= EXP(FRC_URB2D(I,J)*ALOG(ZNT_URB)+(1-FRC_URB2D(I,J))* ALOG(ZNT(I,J))) ! ADD BY DAN
3944 print*, ' FRC_URB2D', FRC_URB2D, &
3945 'ALB_URB',ALB_URB, 'ALBEDOK',ALBEDOK, &
3946 'ALBEDO(I,J)', ALBEDO(I,J), &
3947 'SH_URB',SH_URB,'SHEAT',SHEAT, 'HFX(I,J)',HFX(I,J), &
3948 'LH_KINEMATIC_URB',LH_KINEMATIC_URB,'ETA_KINEMATIC', &
3949 ETA_KINEMATIC, 'QFX(I,J)',QFX(I,J), &
3950 'LH_URB',LH_URB, 'ETA',ETA, 'LH(I,J)',LH(I,J), &
3951 'G_URB',G_URB,'SSOIL',SSOIL,'GRDFLX(I,J)', GRDFLX(I,J),&
3952 'TS_URB',TS_URB,'T1',T1,'TSK(I,J)',TSK(I,J), &
3953 'QS_URB',QS_URB,'Q1',Q1,'QSFC(I,J)',QSFC(I,J)
3957 ! Renew Urban State Varialbes
3959 TR_URB2D(I,J) = TR_URB
3960 TB_URB2D(I,J) = TB_URB
3961 TG_URB2D(I,J) = TG_URB
3962 TC_URB2D(I,J) = TC_URB
3963 QC_URB2D(I,J) = QC_URB
3964 UC_URB2D(I,J) = UC_URB
3966 DO K = 1,num_roof_layers
3967 TRL_URB3D(I,K,J) = TRL_URB(K)
3968 SMR_URB3D(I,K,J) = SMR_URB(K)
3969 TGRL_URB3D(I,K,J)= TGRL_URB(K)
3971 DO K = 1,num_wall_layers
3972 TBL_URB3D(I,K,J) = TBL_URB(K)
3974 DO K = 1,num_road_layers
3975 TGL_URB3D(I,K,J) = TGL_URB(K)
3978 TGR_URB2D(I,J) =TGR_URB
3979 CMCR_URB2D(I,J)=CMCR_URB
3980 FLXHUMR_URB2D(I,J)=FLXHUMR_URB
3981 FLXHUMB_URB2D(I,J)=FLXHUMB_URB
3982 FLXHUMG_URB2D(I,J)=FLXHUMG_URB
3983 DRELR_URB2D(I,J) = DRELR_URB
3984 DRELB_URB2D(I,J) = DRELB_URB
3985 DRELG_URB2D(I,J) = DRELG_URB
3987 XXXR_URB2D(I,J) = XXXR_URB
3988 XXXB_URB2D(I,J) = XXXB_URB
3989 XXXG_URB2D(I,J) = XXXG_URB
3990 XXXC_URB2D(I,J) = XXXC_URB
3992 SH_URB2D(I,J) = SH_URB
3993 LH_URB2D(I,J) = LH_URB
3994 G_URB2D(I,J) = G_URB
3995 RN_URB2D(I,J) = RN_URB
3996 PSIM_URB2D(I,J) = PSIM_URB
3997 PSIH_URB2D(I,J) = PSIH_URB
3998 GZ1OZ0_URB2D(I,J)= GZ1OZ0_URB
3999 U10_URB2D(I,J) = U10_URB
4000 V10_URB2D(I,J) = V10_URB
4001 TH2_URB2D(I,J) = TH2_URB
4002 Q2_URB2D(I,J) = Q2_URB
4003 UST_URB2D(I,J) = UST_URB
4004 AKMS_URB2D(I,J) = KARMAN * UST_URB2D(I,J)/(GZ1OZ0_URB2D(I,J)-PSIM_URB2D(I,J))
4005 IF (PRESENT(CMR_SFCDIF)) THEN
4006 CMR_SFCDIF(I,J) = CMR_URB
4007 CHR_SFCDIF(I,J) = CHR_URB
4008 CMGR_SFCDIF(I,J) = CMGR_URB
4009 CHGR_SFCDIF(I,J) = CHGR_URB
4010 CMC_SFCDIF(I,J) = CMC_URB
4011 CHC_SFCDIF(I,J) = CHC_URB
4014 ! 2D to 3D mosaic danli
4016 TR_URB2D_mosaic(I,mosaic_i,J)=TR_URB2D(I,J)
4017 TB_URB2D_mosaic(I,mosaic_i,J)=TB_URB2D(I,J)
4018 TG_URB2D_mosaic(I,mosaic_i,J)=TG_URB2D(I,J)
4019 TC_URB2D_mosaic(I,mosaic_i,J)=TC_URB2D(I,J)
4020 QC_URB2D_mosaic(I,mosaic_i,J)=QC_URB2D(I,J)
4021 UC_URB2D_mosaic(I,mosaic_i,J)=UC_URB2D(I,J)
4022 TS_URB2D_mosaic(I,mosaic_i,J)=TS_URB2D(I,J)
4023 TS_RUL2D_mosaic(I,mosaic_i,J)=T1
4025 DO K = 1,num_roof_layers
4026 TRL_URB3D_mosaic(I,K+(mosaic_i-1)*num_roof_layers,J)=TRL_URB3D(I,K,J)
4028 DO K = 1,num_wall_layers
4029 TBL_URB3D_mosaic(I,K+(mosaic_i-1)*num_roof_layers,J)=TBL_URB3D(I,K,J)
4031 DO K = 1,num_road_layers
4032 TGL_URB3D_mosaic(I,K+(mosaic_i-1)*num_roof_layers,J)=TGL_URB3D(I,K,J)
4035 SH_URB2D_mosaic(I,mosaic_i,J) = SH_URB2D(I,J)
4036 LH_URB2D_mosaic(I,mosaic_i,J) = LH_URB2D(I,J)
4037 G_URB2D_mosaic(I,mosaic_i,J) = G_URB2D(I,J)
4038 RN_URB2D_mosaic(I,mosaic_i,J) = RN_URB2D(I,J)
4042 ENDIF ! end of UCM CALL if block
4043 !--------------------------------------
4044 ! Urban Part End - urban
4045 !--------------------------------------
4049 SMSTOT(I,J)=SOILM*1000.
4051 SMCREL(I,NS,J)=SMAV(NS)
4054 ! Convert the water unit into mm
4055 SFCRUNOFF(I,J)=SFCRUNOFF(I,J)+RUNOFF1*DT*1000.0
4056 UDRUNOFF(I,J)=UDRUNOFF(I,J)+RUNOFF2*DT*1000.0
4057 ! snow defined when fraction of frozen precip (FFROZP) > 0.5,
4058 IF(FFROZP.GT.0.5)THEN
4059 ACSNOW(I,J)=ACSNOW(I,J)+PRCP*DT
4061 IF(SNOW(I,J).GT.0.)THEN
4062 ACSNOM(I,J)=ACSNOM(I,J)+SNOMLT*1000.
4063 ! accumulated snow-melt energy
4064 SNOPCX(I,J)=SNOPCX(I,J)-SNOMLT/FDTLIW
4067 ENDIF ! endif of land-sea test
4069 !-----------------------------------------------------------------------
4070 ! Done with the Noah-UCM MOSAIC DANLI
4071 !-----------------------------------------------------------------------
4073 TSK_mosaic(i,mosaic_i,j)=TSK(i,j) ! from 2D to 3D
4074 QSFC_mosaic(i,mosaic_i,j)=QSFC(i,j)
4075 CANWAT_mosaic(i,mosaic_i,j)=CANWAT(i,j)
4076 SNOW_mosaic(i,mosaic_i,j)=SNOW(i,j)
4077 SNOWH_mosaic(i,mosaic_i,j)=SNOWH(i,j)
4078 SNOWC_mosaic(i,mosaic_i,j)=SNOWC(i,j)
4080 ALBEDO_mosaic(i,mosaic_i,j)=ALBEDO(i,j)
4081 ALBBCK_mosaic(i,mosaic_i,j)=ALBBCK(i,j)
4082 EMISS_mosaic(i,mosaic_i,j)=EMISS(i,j)
4083 EMBCK_mosaic(i,mosaic_i,j)=EMBCK(i,j)
4084 ZNT_mosaic(i,mosaic_i,j)=ZNT(i,j)
4085 Z0_mosaic(i,mosaic_i,j)=Z0(i,j)
4086 LAI_mosaic(i,mosaic_i,j)=XLAI
4087 RC_mosaic(i,mosaic_i,j)=RC
4089 HFX_mosaic(i,mosaic_i,j)=HFX(i,j)
4090 QFX_mosaic(i,mosaic_i,j)=QFX(i,j)
4091 LH_mosaic(i,mosaic_i,j)=LH(i,j)
4092 GRDFLX_mosaic(i,mosaic_i,j)=GRDFLX(i,j)
4093 SNOTIME_mosaic(i,mosaic_i,j)=SNOTIME(i,j)
4097 TSLB_mosaic(i,NSOIL*(mosaic_i-1)+NS,j)=TSLB(i,NS,j)
4098 SMOIS_mosaic(i,NSOIL*(mosaic_i-1)+NS,j)=SMOIS(i,NS,j)
4099 SH2O_mosaic(i,NSOIL*(mosaic_i-1)+NS,j)=SH2O(i,NS,j)
4104 IF(TSK_mosaic(i,mosaic_i,j) > 350 .OR. TSK_mosaic(i,mosaic_i,j) < 250 .OR. abs(HFX_mosaic(i,mosaic_i,j)) > 700 ) THEN
4105 print*, 'I', I, 'J', J, 'MOSAIC_I', MOSAIC_I
4106 print*, 'mosaic_cat_index',mosaic_cat_index(I,mosaic_i,J), 'landusef2',landusef2(i,mosaic_i,j)
4107 print*, 'TSK_mosaic', TSK_mosaic(i,mosaic_i,j), 'HFX_mosaic', HFX_mosaic(i,mosaic_i,j), &
4108 'LH_mosaic',LH_mosaic(i,mosaic_i,j),'GRDFLX_mosaic',GRDFLX_mosaic(i,mosaic_i,j)
4109 print*, 'ZNT_mosaic', ZNT_mosaic(i, mosaic_i,j), 'Z0_mosaic', Z0_mosaic(i,mosaic_i,j)
4110 print*, 'LAI_mosaic', LAI_mosaic(i, mosaic_i,j)
4111 print*, 'FRC_URB2D',FRC_URB2D(I,J)
4112 print*, 'TS_URB',TS_URB2D(I,J),'T1',T1
4113 print*, 'SH_URB2D',SH_URB2D(I,J),'SHEAT',SHEAT
4114 print*, 'LH_URB',LH_URB2D(I,J),'ETA',ETA
4115 print*, 'TS_RUL2D',TS_RUL2D_mosaic(I,mosaic_i,J)
4120 !-----------------------------------------------------------------------
4121 ! Now let's do the grid-averaging
4122 !-----------------------------------------------------------------------
4124 FAREA = landusef2(i,mosaic_i,j)
4126 TSK_mosaic_avg(i,j) = TSK_mosaic_avg(i,j) + (EMISS_mosaic(i,mosaic_i,j)*TSK_mosaic(i,mosaic_i,j)**4)*FAREA ! conserve the longwave radiation
4128 QSFC_mosaic_avg(i,j) = QSFC_mosaic_avg(i,j) + QSFC_mosaic(i,mosaic_i,j)*FAREA
4129 CANWAT_mosaic_avg(i,j) = CANWAT_mosaic_avg(i,j) + CANWAT_mosaic(i,mosaic_i,j)*FAREA
4130 SNOW_mosaic_avg(i,j) = SNOW_mosaic_avg(i,j) + SNOW_mosaic(i,mosaic_i,j)*FAREA
4131 SNOWH_mosaic_avg(i,j) = SNOWH_mosaic_avg(i,j) + SNOWH_mosaic(i,mosaic_i,j)*FAREA
4132 SNOWC_mosaic_avg(i,j) = SNOWC_mosaic_avg(i,j) + SNOWC_mosaic(i,mosaic_i,j)*FAREA
4136 TSLB_mosaic_avg(i,NS,j)=TSLB_mosaic_avg(i,NS,j) + TSLB_mosaic(i,NS*mosaic_i,j)*FAREA
4137 SMOIS_mosaic_avg(i,NS,j)=SMOIS_mosaic_avg(i,NS,j) + SMOIS_mosaic(i,NS*mosaic_i,j)*FAREA
4138 SH2O_mosaic_avg(i,NS,j)=SH2O_mosaic_avg(i,NS,j) + SH2O_mosaic(i,NS*mosaic_i,j)*FAREA
4142 FAREA_mosaic_avg(i,j)=FAREA_mosaic_avg(i,j)+FAREA
4143 HFX_mosaic_avg(i,j) = HFX_mosaic_avg(i,j) + HFX_mosaic(i,mosaic_i,j)*FAREA
4144 QFX_mosaic_avg(i,j) = QFX_mosaic_avg(i,j) + QFX_mosaic(i,mosaic_i,j)*FAREA
4145 LH_mosaic_avg(i,j) = LH_mosaic_avg(i,j) + LH_mosaic(i,mosaic_i,j)*FAREA
4146 GRDFLX_mosaic_avg(i,j)=GRDFLX_mosaic_avg(i,j)+GRDFLX_mosaic(i,mosaic_i,j)*FAREA
4148 ALBEDO_mosaic_avg(i,j)=ALBEDO_mosaic_avg(i,j)+ALBEDO_mosaic(i,mosaic_i,j)*FAREA
4149 ALBBCK_mosaic_avg(i,j)=ALBBCK_mosaic_avg(i,j)+ALBBCK_mosaic(i,mosaic_i,j)*FAREA
4150 EMISS_mosaic_avg(i,j)=EMISS_mosaic_avg(i,j)+EMISS_mosaic(i,mosaic_i,j)*FAREA
4151 EMBCK_mosaic_avg(i,j)=EMBCK_mosaic_avg(i,j)+EMBCK_mosaic(i,mosaic_i,j)*FAREA
4152 ZNT_mosaic_avg(i,j)=ZNT_mosaic_avg(i,j)+ALOG(ZNT_mosaic(i,mosaic_i,j))*FAREA
4153 Z0_mosaic_avg(i,j)=Z0_mosaic_avg(i,j)+ALOG(Z0_mosaic(i,mosaic_i,j))*FAREA
4154 LAI_mosaic_avg(i,j)=LAI_mosaic_avg(i,j)+LAI_mosaic(i,mosaic_i,j)*FAREA
4155 if(RC_mosaic(i,mosaic_i,j) .Gt. 0.0) Then
4156 RC_mosaic_avg(i,j) = RC_mosaic_avg(i,j)+1.0/RC_mosaic(i,mosaic_i,j)*FAREA
4158 RC_mosaic_avg(i,j) = RC_mosaic_avg(i,j) + RC_mosaic(i,mosaic_i,j)*FAREA
4160 ENDDO ! ENDDO FOR mosaic_i = 1, mosaic_cat
4162 !-----------------------------------------------------------------------
4163 ! Now let's send the 3D values to the 2D variables that might be needed in other routines
4164 !-----------------------------------------------------------------------
4166 IVGTYP(I,J)=IVGTYP_dominant(I,J) ! the dominant vege category
4167 ALBEDO(i,j)=ALBEDO_mosaic_avg(i,j)
4168 ALBBCK(i,j)=ALBBCK_mosaic_avg(i,j)
4169 EMISS(i,j)= EMISS_mosaic_avg(i,j)
4170 EMBCK(i,j)= EMBCK_mosaic_avg(i,j)
4171 ZNT(i,j)= EXP(ZNT_mosaic_avg(i,j)/FAREA_mosaic_avg(i,j))
4172 Z0(i,j)= EXP(Z0_mosaic_avg(i,j)/FAREA_mosaic_avg(i,j))
4173 XLAI2(i,j)= LAI_mosaic_avg(i,j)
4174 IF (RC_mosaic_avg(i,j) .Gt. 0.0) THEN
4175 rc2(i,j) = 1.0/(RC_mosaic_avg(i,j))
4177 !RC_mosaic_avg was zero for all tiles (cell over water), thus RC2 set to zero to avoid infinity
4178 rc2(i,j) = RC_mosaic_avg(i,j)
4180 TSK(i,j)=(TSK_mosaic_avg(I,J)/EMISS_mosaic_avg(I,J))**(0.25) ! from 3D to 2D
4181 QSFC(i,j)=QSFC_mosaic_avg(I,J)
4182 CANWAT(i,j) = CANWAT_mosaic_avg(i,j)
4183 SNOW(i,j) = SNOW_mosaic_avg(i,j)
4184 SNOWH(i,j) = SNOWH_mosaic_avg(i,j)
4185 SNOWC(i,j) = SNOWC_mosaic_avg(i,j)
4187 HFX(i,j) = HFX_mosaic_avg(i,j)
4188 QFX(i,j) = QFX_mosaic_avg(i,j)
4189 LH(i,j) = LH_mosaic_avg(i,j)
4190 GRDFLX(i,j)=GRDFLX_mosaic_avg(i,j)
4194 TSLB(i,NS,j)=TSLB_mosaic_avg(i,NS,j)
4195 SMOIS(i,NS,j)=SMOIS_mosaic_avg(i,NS,j)
4196 SH2O(i,NS,j)=SH2O_mosaic_avg(i,NS,j)
4200 ELSE ! This corresponds to IF ((sf_surface_mosaic == 1) .AND. ((XLAND(I,J)-1.5).LT.0.) .AND. (XICE(I,J) < XICE_THRESHOLD) ) THEN
4204 ! pressure in middle of lowest layer
4205 SFCPRS=(P8W3D(I,KTS+1,j)+P8W3D(i,KTS,j))*0.5
4206 ! convert from mixing ratio to specific humidity
4207 Q2K=QV3D(i,1,j)/(1.0+QV3D(i,1,j))
4210 Q2SAT=QGH(I,J)/(1.0+QGH(I,J)) ! Q2SAT is sp humidity
4211 ! add check on myj=.true.
4212 ! IF((Q2K.GE.Q2SAT*TRESH).AND.Q2K.LT.QZ0(I,J))THEN
4213 IF((myj).AND.(Q2K.GE.Q2SAT*TRESH).AND.Q2K.LT.QZ0(I,J))THEN
4222 ZLVL=0.5*DZ8W(i,1,j)
4224 ! TH2=SFCTMP+(0.0097545*ZLVL)
4225 ! calculate SFCTH2 via Exner function vs lapse-rate (above)
4226 APES=(1.E5/PSFC)**CAPA
4227 APELM=(1.E5/SFCPRS)**CAPA
4232 LWDN=GLW(I,J)*EMISSI
4233 ! SOLDN is total incoming solar
4235 ! GSW is net downward solar
4237 ! use mid-day albedo to determine net downward solar (no solar zenith angle correction)
4238 SOLNET=SOLDN*(1.-ALBEDO(I,J))
4242 SHDFAC=VEGFRA(I,J)/100.
4245 SHMIN=SHDMIN(I,J)/100. !NEW
4246 SHMAX=SHDMAX(I,J)/100. !NEW
4247 ! convert snow water equivalent from mm to meter
4248 SNEQV=SNOW(I,J)*0.001
4249 ! snow depth in meters
4253 ! if "SR" present, set frac of frozen precip ("FFROZP") = snow-ratio ("SR", range:0-1)
4254 ! SR from e.g. Ferrier microphysics
4255 ! otherwise define from 1st atmos level temperature
4259 IF (SFCTMP <= 273.15) THEN
4266 IF((XLAND(I,J)-1.5).GE.0.)THEN ! begining of land/sea if block
4268 TSK_RURAL(I,J)=TSK(I,J)
4269 HFX_RURAL(I,J)=HFX(I,J)
4270 QFX_RURAL(I,J)=QFX(I,J)
4271 LH_RURAL(I,J)=LH(I,J)
4272 EMISS_RURAL(I,J)=EMISS(I,J)
4273 GRDFLX_RURAL(I,J)=GRDFLX(I,J)
4275 ! Land or sea-ice case
4277 IF (XICE(I,J) >= XICE_THRESHOLD) THEN
4280 ELSE IF ( VEGTYP == ISICE ) THEN
4284 ! Neither sea ice or land ice.
4287 DQSDT2=Q2SAT*A23M4/(SFCTMP-A4)**2
4289 IF(SNOW(I,J).GT.0.0)THEN
4290 ! snow on surface (use ice saturation properties)
4292 E2SAT=611.2*EXP(6174.*(1./273.15 - 1./SFCTSNO))
4293 Q2SATI=0.622*E2SAT/(SFCPRS-E2SAT)
4294 Q2SATI=Q2SATI/(1.0+Q2SATI) ! spec. hum.
4295 IF (T1 .GT. 273.14) THEN
4296 ! warm ground temps, weight the saturation between ice and water according to SNOWC
4297 Q2SAT=Q2SAT*(1.-SNOWC(I,J)) + Q2SATI*SNOWC(I,J)
4298 DQSDT2=DQSDT2*(1.-SNOWC(I,J)) + Q2SATI*6174./(SFCTSNO**2)*SNOWC(I,J)
4300 ! cold ground temps, use ice saturation only
4302 DQSDT2=Q2SATI*6174./(SFCTSNO**2)
4304 ! for snow cover fraction at 0 C, ground temp will not change, so DQSDT2 effectively zero
4305 IF(T1 .GT. 273. .AND. SNOWC(I,J) .GT. 0.)DQSDT2=DQSDT2*(1.-SNOWC(I,J))
4308 ! Land-ice or land points use the usual deep-soil temperature.
4311 IF(VEGTYP.EQ.25) SHDFAC=0.0000
4312 IF(VEGTYP.EQ.26) SHDFAC=0.0000
4313 IF(VEGTYP.EQ.27) SHDFAC=0.0000
4314 IF(SOILTYP.EQ.14.AND.XICE(I,J).EQ.0.)THEN
4316 IF(IPRINT)PRINT*,' SOIL TYPE FOUND TO BE WATER AT A LAND-POINT'
4317 IF(IPRINT)PRINT*,i,j,'RESET SOIL in surfce.F'
4321 SNOALB1 = SNOALB(I,J)
4322 CMC=CANWAT(I,J)/1000.
4324 !-------------------------------------------
4325 !*** convert snow depth from mm to meter
4328 ! SNOALB=ALBMAX(I,J)*0.01
4330 ! SNOALB=MAXALB(IVGTPK)*0.01
4338 SNOTIME1 = SNOTIME(I,J)
4340 !FEI: temporaray arrays above need to be changed later by using SI
4343 SMC(NS)=SMOIS(I,NS,J)
4344 STC(NS)=TSLB(I,NS,J) !STEMP
4345 SWC(NS)=SH2O(I,NS,J)
4348 if ( (SNEQV.ne.0..AND.SNOWHK.eq.0.).or.(SNOWHK.le.SNEQV) )THEN
4353 !Fei: urban. for urban surface, if calling UCM, redefine the natural surface in cities as
4354 ! the "NATURAL" category in the VEGPARM.TBL
4355 IF(SF_URBAN_PHYSICS == 1.OR. SF_URBAN_PHYSICS==2.OR.SF_URBAN_PHYSICS==3 ) THEN
4357 IF( IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == LCZ_1 .or. IVGTYP(I,J) == LCZ_2 .or. &
4358 IVGTYP(I,J) == LCZ_3 .or. IVGTYP(I,J) == LCZ_4 .or. IVGTYP(I,J) == LCZ_5 .or. &
4359 IVGTYP(I,J) == LCZ_6 .or. IVGTYP(I,J) == LCZ_7 .or. IVGTYP(I,J) == LCZ_8 .or. &
4360 IVGTYP(I,J) == LCZ_9 .or. IVGTYP(I,J) == LCZ_10 .or. IVGTYP(I,J) == LCZ_11 ) THEN
4363 SHDFAC = SHDTBL(NATURAL)
4366 EMISSI = 0.98 !for VEGTYP=5
4367 LWDN = GLW(I,J) * EMISSI
4368 SOLNET = SOLDN * (1.0 - ALBEDOK)
4370 IF ( FRC_URB2D(I,J) < 0.99 ) THEN
4371 if(sf_urban_physics.eq.1)then
4372 T1= ( TSK(I,J) -FRC_URB2D(I,J) * TS_URB2D (I,J) )/ (1-FRC_URB2D(I,J))
4373 elseif((sf_urban_physics.eq.2).OR.(sf_urban_physics.eq.3))then
4374 T1=tsk_rural_bep(i,j)
4381 IF( IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == LCZ_1 .or. IVGTYP(I,J) == LCZ_2 .or. &
4382 IVGTYP(I,J) == LCZ_3 .or. IVGTYP(I,J) == LCZ_4 .or. IVGTYP(I,J) == LCZ_5 .or. &
4383 IVGTYP(I,J) == LCZ_6 .or. IVGTYP(I,J) == LCZ_7 .or. IVGTYP(I,J) == LCZ_8 .or. &
4384 IVGTYP(I,J) == LCZ_9 .or. IVGTYP(I,J) == LCZ_10 .or. IVGTYP(I,J) == LCZ_11 ) THEN
4390 !===Yang, 2014/10/08, hydrological processes for urban vegetation in single layer UCM===
4394 IRIOPTION=IRI_SCHEME
4396 tloc=mod(int(OMG/3.14159*180./15.+12.+0.5 ),24)
4397 if (tloc.lt.0) tloc=tloc+24
4398 if (tloc==0) tloc=24
4399 CALL cal_mon_day(julian,julyr,jmonth,jday)
4400 IF(SF_URBAN_PHYSICS == 1) THEN
4401 IF( IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == LCZ_1 .or. IVGTYP(I,J) == LCZ_2 .or. &
4402 IVGTYP(I,J) == LCZ_3 .or. IVGTYP(I,J) == LCZ_4 .or. IVGTYP(I,J) == LCZ_5 .or. &
4403 IVGTYP(I,J) == LCZ_6 .or. IVGTYP(I,J) == LCZ_7 .or. IVGTYP(I,J) == LCZ_8 .or. &
4404 IVGTYP(I,J) == LCZ_9 .or. IVGTYP(I,J) == LCZ_10 .or. IVGTYP(I,J) == LCZ_11 ) THEN
4406 AOASIS = oasis ! urban oasis effect
4407 IF (IRIOPTION ==1) THEN
4408 IF (tloc==21 .or. tloc==22) THEN !irrigation on vegetaion in urban area, MAY-SEP, 9-10pm
4409 IF (jmonth==5 .or. jmonth==6 .or. jmonth==7 .or. jmonth==8 .or. jmonth==9) THEN
4410 IF (SMC(USOIL) .LT. SMCREF) SMC(USOIL)= REFSMC(ISLTYP(I,J))
4411 IF (SMC(DSOIL) .LT. SMCREF) SMC(DSOIL)= REFSMC(ISLTYP(I,J))
4418 IF(SF_URBAN_PHYSICS == 2 .or. SF_URBAN_PHYSICS == 3) THEN
4419 IF(AOASIS > 1.0) THEN
4420 CALL wrf_error_fatal('Urban oasis option is for SF_URBAN_PHYSICS == 1 only')
4422 IF(IRIOPTION == 1) THEN
4423 CALL wrf_error_fatal('Urban irrigation option is for SF_URBAN_PHYSICS == 1 only')
4430 print*, 'BEFORE SFLX, in Noahlsm_driver'
4431 print*, 'ICE', ICE, 'DT',DT, 'ZLVL',ZLVL, 'NSOIL', NSOIL, &
4432 'SLDPTH', SLDPTH, 'LOCAL',LOCAL, 'LUTYPE',&
4433 LUTYPE, 'SLTYPE',SLTYPE, 'LWDN',LWDN, 'SOLDN',SOLDN, &
4434 'SFCPRS',SFCPRS, 'PRCP',PRCP,'SFCTMP',SFCTMP,'Q2K',Q2K, &
4435 'TH2',TH2,'Q2SAT',Q2SAT,'DQSDT2',DQSDT2,'VEGTYP', VEGTYP,&
4436 'SOILTYP',SOILTYP, 'SLOPETYP',SLOPETYP, 'SHDFAC',SHDFAC,&
4437 'SHMIN',SHMIN, 'ALBBRD',ALBBRD,'SNOALB1',SNOALB1,'TBOT',&
4438 TBOT, 'Z0BRD',Z0BRD, 'Z0K',Z0K, 'CMC',CMC, 'T1',T1,'STC',&
4439 STC, 'SMC',SMC, 'SWC',SWC,'SNOWHK',SNOWHK,'SNEQV',SNEQV,&
4440 'ALBEDOK',ALBEDOK,'CHK',CHK,'ETA',ETA,'SHEAT',SHEAT, &
4441 'ETA_KINEMATIC',ETA_KINEMATIC, 'FDOWN',FDOWN,'EC',EC, &
4442 'EDIR',EDIR,'ET',ET,'ETT',ETT,'ESNOW',ESNOW,'DRIP',DRIP,&
4443 'DEW',DEW,'BETA',BETA,'ETP',ETP,'SSOIL',SSOIL,'FLX1',FLX1,&
4444 'FLX2',FLX2,'FLX3',FLX3,'SNOMLT',SNOMLT,'SNCOVR',SNCOVR,&
4445 'RUNOFF1',RUNOFF1,'RUNOFF2',RUNOFF2,'RUNOFF3',RUNOFF3, &
4446 'RC',RC, 'PC',PC,'RSMIN',RSMIN,'XLAI',XLAI,'RCS',RCS, &
4447 'RCT',RCT,'RCQ',RCQ,'RCSOIL',RCSOIL,'SOILW',SOILW, &
4448 'SOILM',SOILM,'Q1',Q1,'SMCWLT',SMCWLT,'SMCDRY',SMCDRY,&
4449 'SMCREF',SMCREF,'SMCMAX',SMCMAX,'NROOT',NROOT
4457 IF ( ICE == 1 ) THEN
4468 ELSEIF (ICE == 0) THEN
4472 CALL SFLX (I,J,FFROZP, ISURBAN, DT,ZLVL,NSOIL,SLDPTH, & !C
4474 LUTYPE, SLTYPE, & !CL
4475 LWDN,SOLDN,SOLNET,SFCPRS,PRCP,SFCTMP,Q2K,DUMMY, & !F
4476 DUMMY,DUMMY, DUMMY, & !F PRCPRAIN not used
4477 TH2,Q2SAT,DQSDT2, & !I
4478 VEGTYP,SOILTYP,SLOPETYP,SHDFAC,SHMIN,SHMAX, & !I
4479 ALBBRD, SNOALB1,TBOT, Z0BRD, Z0K, EMISSI, EMBRD, & !S
4480 CMC,T1,STC,SMC,SWC,SNOWHK,SNEQV,ALBEDOK,CHK,dummy,& !H
4481 ETA,SHEAT, ETA_KINEMATIC,FDOWN, & !O
4482 EC,EDIR,ET,ETT,ESNOW,DRIP,DEW, & !O
4483 BETA,ETP,SSOIL, & !O
4484 FLX1,FLX2,FLX3, & !O
4485 FLX4,FVB,FBUR,FGSN,UA_PHYS, & !UA
4487 RUNOFF1,RUNOFF2,RUNOFF3, & !O
4488 RC,PC,RSMIN,XLAI,RCS,RCT,RCQ,RCSOIL, & !O
4489 SOILW,SOILM,Q1,SMAV, & !D
4490 RDLAI2D,USEMONALB, &
4493 SMCWLT,SMCDRY,SMCREF,SMCMAX,NROOT, &
4495 sfcheadrt_hydro, & !I
4496 INFXSRT_hydro,ETPND1_hydro & !O
4497 ,OPT_THCND,AOASIS & !O
4498 ,XSDA_QFX, HFX_PHY, QFX_PHY, XQNORM, fasdas, HCPCT_FASDAS & ! fasdas vars
4499 ,IRRIGATION_CHANNEL )
4502 soldrain(i,j) = RUNOFF2*DT*1000.0
4503 sfcheadrt(i,j) = sfcheadrt_hydro
4504 infxsrt(i,j) = INFXSRT_hydro
4505 etpnd1 = etpnd1_hydro
4507 ELSEIF (ICE == -1) THEN
4510 ! Set values that the LSM is expected to update,
4511 ! but don't get updated for glacial points.
4513 SOILM = 0.0 !BSINGH(PNNL)- SOILM is undefined for this case, it is used for diagnostics so setting it to zero
4514 XLAI = 0.01 ! KWM Should this be Zero over land ice? Does this value matter?
4522 CALL SFLX_GLACIAL(I,J,ISICE,FFROZP,DT,ZLVL,NSOIL,SLDPTH, & !C
4523 & LWDN,SOLNET,SFCPRS,PRCP,SFCTMP,Q2K, & !F
4524 & TH2,Q2SAT,DQSDT2, & !I
4525 & ALBBRD, SNOALB1,TBOT, Z0BRD, Z0K, EMISSI, EMBRD, & !S
4526 & T1,STC(1:NSOIL),SNOWHK,SNEQV,ALBEDOK,CHK, & !H
4527 & ETA,SHEAT,ETA_KINEMATIC,FDOWN, & !O
4530 & FLX1,FLX2,FLX3, & !O
4531 & SNOMLT,SNCOVR, & !O
4544 print*, 'AFTER SFLX, in Noahlsm_driver'
4545 print*, 'ICE', ICE, 'DT',DT, 'ZLVL',ZLVL, 'NSOIL', NSOIL, &
4546 'SLDPTH', SLDPTH, 'LOCAL',LOCAL, 'LUTYPE',&
4547 LUTYPE, 'SLTYPE',SLTYPE, 'LWDN',LWDN, 'SOLDN',SOLDN, &
4548 'SFCPRS',SFCPRS, 'PRCP',PRCP,'SFCTMP',SFCTMP,'Q2K',Q2K, &
4549 'TH2',TH2,'Q2SAT',Q2SAT,'DQSDT2',DQSDT2,'VEGTYP', VEGTYP,&
4550 'SOILTYP',SOILTYP, 'SLOPETYP',SLOPETYP, 'SHDFAC',SHDFAC,&
4551 'SHDMIN',SHMIN, 'ALBBRD',ALBBRD,'SNOALB',SNOALB1,'TBOT',&
4552 TBOT, 'Z0BRD',Z0BRD, 'Z0K',Z0K, 'CMC',CMC, 'T1',T1,'STC',&
4553 STC, 'SMC',SMC, 'SWc',SWC,'SNOWHK',SNOWHK,'SNEQV',SNEQV,&
4554 'ALBEDOK',ALBEDOK,'CHK',CHK,'ETA',ETA,'SHEAT',SHEAT, &
4555 'ETA_KINEMATIC',ETA_KINEMATIC, 'FDOWN',FDOWN,'EC',EC, &
4556 'EDIR',EDIR,'ET',ET,'ETT',ETT,'ESNOW',ESNOW,'DRIP',DRIP,&
4557 'DEW',DEW,'BETA',BETA,'ETP',ETP,'SSOIL',SSOIL,'FLX1',FLX1,&
4558 'FLX2',FLX2,'FLX3',FLX3,'SNOMLT',SNOMLT,'SNCOVR',SNCOVR,&
4559 'RUNOFF1',RUNOFF1,'RUNOFF2',RUNOFF2,'RUNOFF3',RUNOFF3, &
4560 'RC',RC, 'PC',PC,'RSMIN',RSMIN,'XLAI',XLAI,'RCS',RCS, &
4561 'RCT',RCT,'RCQ',RCQ,'RCSOIL',RCSOIL,'SOILW',SOILW, &
4562 'SOILM',SOILM,'Q1',Q1,'SMCWLT',SMCWLT,'SMCDRY',SMCDRY,&
4563 'SMCREF',SMCREF,'SMCMAX',SMCMAX,'NROOT',NROOT
4567 !*** UPDATE STATE VARIABLES
4568 CANWAT(I,J)=CMC*1000.
4569 SNOW(I,J)=SNEQV*1000.
4570 ! SNOWH(I,J)=SNOWHK*1000.
4571 SNOWH(I,J)=SNOWHK ! SNOWHK in meters
4573 ALB_RURAL(I,J)=ALBEDOK
4577 EMISS_RURAL(I,J) = EMISSI
4578 ! Noah: activate time-varying roughness length (V3.3 Feb 2011)
4583 HFX_RURAL(I,J)=SHEAT
4584 ! MEk Jul07 add potential evap accum
4585 POTEVP(I,J)=POTEVP(I,J)+ETP*FDTW
4586 QFX(I,J)=ETA_KINEMATIC
4587 QFX_RURAL(I,J)=ETA_KINEMATIC
4591 ! QFX(I,J) = QFX(I,J) + ETPND1
4592 ! ETA = ETA + ETPND1/2.501E6*dt
4593 !end added by Wei Yu
4599 GRDFLX_RURAL(I,J)=SSOIL
4602 SNOTIME(I,J) = SNOTIME1
4603 ! prevent diagnostic ground q (q1) from being greater than qsat(tsk)
4604 ! as happens over snow cover where the cqs2 value also becomes irrelevant
4605 ! by setting cqs2=chs in this situation the 2m q should become just qv(k=1)
4606 IF (Q1 .GT. QSFC(I,J)) THEN
4607 CQS2(I,J) = CHS(I,J)
4610 ! Convert QSFC back to mixing ratio
4611 QSFC(I,J)= Q1/(1.0-Q1)
4613 ! QSFC_RURAL(I,J)= Q1/(1.0-Q1)
4614 ! Calculate momentum flux from rural surface for use with multi-layer UCM (Martilli et al. 2002)
4617 SMOIS(I,NS,J)=SMC(NS)
4618 TSLB(I,NS,J)=STC(NS) ! STEMP
4619 SH2O(I,NS,J)=SWC(NS)
4628 ! Residual of surface energy balance equation terms
4632 noahres(i,j) = ( solnet + lwdn ) - sheat + ssoil - eta &
4633 - ( emissi * STBOLT * (t1**4) ) - flx1 - flx2 - flx3 - flx4
4636 noahres(i,j) = ( solnet + lwdn ) - sheat + ssoil - eta &
4637 - ( emissi * STBOLT * (t1**4) ) - flx1 - flx2 - flx3
4640 IF (SF_URBAN_PHYSICS == 1 ) THEN ! Beginning of UCM CALL if block
4641 !--------------------------------------
4642 ! URBAN CANOPY MODEL START - urban
4643 !--------------------------------------
4644 ! Input variables lsm --> urban
4645 IF( IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == LCZ_1 .or. IVGTYP(I,J) == LCZ_2 .or. &
4646 IVGTYP(I,J) == LCZ_3 .or. IVGTYP(I,J) == LCZ_4 .or. IVGTYP(I,J) == LCZ_5 .or. &
4647 IVGTYP(I,J) == LCZ_6 .or. IVGTYP(I,J) == LCZ_7 .or. IVGTYP(I,J) == LCZ_8 .or. &
4648 IVGTYP(I,J) == LCZ_9 .or. IVGTYP(I,J) == LCZ_10 .or. IVGTYP(I,J) == LCZ_11 ) THEN
4653 UTYPE_URB = UTYPE_URB2D(I,J) !urban type (low, high or industrial)
4655 TA_URB = SFCTMP ! [K]
4656 QA_URB = Q2K ! [kg/kg]
4657 UA_URB = SQRT(U_PHY(I,1,J)**2.+V_PHY(I,1,J)**2.)
4658 U1_URB = U_PHY(I,1,J)
4659 V1_URB = V_PHY(I,1,J)
4660 IF(UA_URB < 1.) UA_URB=1. ! [m/s]
4661 SSG_URB = SOLDN ! [W/m/m]
4662 SSGD_URB = 0.8*SOLDN ! [W/m/m]
4663 SSGQ_URB = SSG_URB-SSGD_URB ! [W/m/m]
4664 LLG_URB = GLW(I,J) ! [W/m/m]
4665 RAIN_URB = RAINBL(I,J) / DT * 3600.0 ! [mm/hr]
4666 RHOO_URB = SFCPRS / (287.04 * SFCTMP * (1.0+ 0.61 * Q2K)) ![kg/m/m/m]
4668 DELT_URB = DT ! [sec]
4669 XLAT_URB = XLAT_URB2D(I,J) ! [deg]
4670 COSZ_URB = COSZ_URB2D(I,J) !
4671 OMG_URB = OMG_URB2D(I,J) !
4674 LSOLAR_URB = .FALSE.
4676 TR_URB = TR_URB2D(I,J)
4677 TB_URB = TB_URB2D(I,J)
4678 TG_URB = TG_URB2D(I,J)
4679 TC_URB = TC_URB2D(I,J)
4680 QC_URB = QC_URB2D(I,J)
4681 UC_URB = UC_URB2D(I,J)
4683 DO K = 1,num_roof_layers
4684 TRL_URB(K) = TRL_URB3D(I,K,J)
4685 SMR_URB(K) = SMR_URB3D(I,K,J)
4686 TGRL_URB(K)= TGRL_URB3D(I,K,J)
4688 DO K = 1,num_wall_layers
4689 TBL_URB(K) = TBL_URB3D(I,K,J)
4691 DO K = 1,num_road_layers
4692 TGL_URB(K) = TGL_URB3D(I,K,J)
4695 TGR_URB = TGR_URB2D(I,J)
4696 CMCR_URB = CMCR_URB2D(I,J)
4697 FLXHUMR_URB = FLXHUMR_URB2D(I,J)
4698 FLXHUMB_URB = FLXHUMB_URB2D(I,J)
4699 FLXHUMG_URB = FLXHUMG_URB2D(I,J)
4700 DRELR_URB = DRELR_URB2D(I,J)
4701 DRELB_URB = DRELB_URB2D(I,J)
4702 DRELG_URB = DRELG_URB2D(I,J)
4704 XXXR_URB = XXXR_URB2D(I,J)
4705 XXXB_URB = XXXB_URB2D(I,J)
4706 XXXG_URB = XXXG_URB2D(I,J)
4707 XXXC_URB = XXXC_URB2D(I,J)
4709 ! Limits to avoid dividing by small number
4710 if (CHS(I,J) < 1.0E-02) then
4713 if (CHS2(I,J) < 1.0E-02) then
4716 if (CQS2(I,J) < 1.0E-02) then
4721 CHS2_URB = CHS2(I,J)
4722 IF (PRESENT(CMR_SFCDIF)) THEN
4723 CMR_URB = CMR_SFCDIF(I,J)
4724 CHR_URB = CHR_SFCDIF(I,J)
4725 CMGR_URB = CMGR_SFCDIF(I,J)
4726 CHGR_URB = CHGR_SFCDIF(I,J)
4727 CMC_URB = CMC_SFCDIF(I,J)
4728 CHC_URB = CHC_SFCDIF(I,J)
4732 mh_urb = mh_urb2d(I,J)
4733 stdh_urb = stdh_urb2d(I,J)
4734 lp_urb = lp_urb2d(I,J)
4735 hgt_urb = hgt_urb2d(I,J)
4738 lf_urb(K)=lf_urb2d(I,K,J)
4740 frc_urb = frc_urb2d(I,J)
4741 lb_urb = lb_urb2d(I,J)
4743 if (I.eq.73.and.J.eq.125)THEN
4748 CALL cal_mon_day(julian,julyr,jmonth,jday)
4749 CALL urban(LSOLAR_URB, & ! I
4750 num_roof_layers,num_wall_layers,num_road_layers, & ! C
4752 UTYPE_URB,TA_URB,QA_URB,UA_URB,U1_URB,V1_URB,SSG_URB, & ! I
4753 SSGD_URB,SSGQ_URB,LLG_URB,RAIN_URB,RHOO_URB, & ! I
4754 ZA_URB,DECLIN_URB,COSZ_URB,OMG_URB, & ! I
4755 XLAT_URB,DELT_URB,ZNT_URB, & ! I
4756 CHS_URB, CHS2_URB, & ! I
4757 TR_URB, TB_URB, TG_URB, TC_URB, QC_URB,UC_URB, & ! H
4758 TRL_URB,TBL_URB,TGL_URB, & ! H
4759 XXXR_URB, XXXB_URB, XXXG_URB, XXXC_URB, & ! H
4760 TS_URB,QS_URB,SH_URB,LH_URB,LH_KINEMATIC_URB, & ! O
4761 SW_URB,ALB_URB,LW_URB,G_URB,RN_URB,PSIM_URB,PSIH_URB, & ! O
4763 CMR_URB, CHR_URB, CMC_URB, CHC_URB, &
4764 U10_URB, V10_URB, TH2_URB, Q2_URB, & ! O
4765 UST_URB,mh_urb, stdh_urb, lf_urb, lp_urb, & ! 0
4766 hgt_urb,frc_urb,lb_urb, check,CMCR_URB,TGR_URB, & ! H
4767 TGRL_URB,SMR_URB,CMGR_URB,CHGR_URB,jmonth, & ! H
4768 DRELR_URB,DRELB_URB, & ! H
4769 DRELG_URB,FLXHUMR_URB,FLXHUMB_URB,FLXHUMG_URB)
4774 print*, 'AFTER CALL URBAN'
4775 print*,'num_roof_layers',num_roof_layers, 'num_wall_layers', &
4777 'DZR',DZR,'DZB',DZB,'DZG',DZG,'UTYPE_URB',UTYPE_URB,'TA_URB', &
4779 'QA_URB',QA_URB,'UA_URB',UA_URB,'U1_URB',U1_URB,'V1_URB', &
4781 'SSG_URB',SSG_URB,'SSGD_URB',SSGD_URB,'SSGQ_URB',SSGQ_URB, &
4782 'LLG_URB',LLG_URB,'RAIN_URB',RAIN_URB,'RHOO_URB',RHOO_URB, &
4783 'ZA_URB',ZA_URB, 'DECLIN_URB',DECLIN_URB,'COSZ_URB',COSZ_URB,&
4784 'OMG_URB',OMG_URB,'XLAT_URB',XLAT_URB,'DELT_URB',DELT_URB, &
4785 'ZNT_URB',ZNT_URB,'TR_URB',TR_URB, 'TB_URB',TB_URB,'TG_URB',&
4786 TG_URB,'TC_URB',TC_URB,'QC_URB',QC_URB,'TRL_URB',TRL_URB, &
4787 'TBL_URB',TBL_URB,'TGL_URB',TGL_URB,'XXXR_URB',XXXR_URB, &
4788 'XXXB_URB',XXXB_URB,'XXXG_URB',XXXG_URB,'XXXC_URB',XXXC_URB,&
4789 'TS_URB',TS_URB,'QS_URB',QS_URB,'SH_URB',SH_URB,'LH_URB', &
4790 LH_URB, 'LH_KINEMATIC_URB',LH_KINEMATIC_URB,'SW_URB',SW_URB,&
4791 'ALB_URB',ALB_URB,'LW_URB',LW_URB,'G_URB',G_URB,'RN_URB', &
4792 RN_URB, 'PSIM_URB',PSIM_URB,'PSIH_URB',PSIH_URB, &
4793 'U10_URB',U10_URB,'V10_URB',V10_URB,'TH2_URB',TH2_URB, &
4794 'Q2_URB',Q2_URB,'CHS_URB',CHS_URB,'CHS2_URB',CHS2_URB
4798 TS_URB2D(I,J) = TS_URB
4800 ALBEDO(I,J) = FRC_URB2D(I,J)*ALB_URB+(1-FRC_URB2D(I,J))*ALBEDOK ![-]
4801 HFX(I,J) = FRC_URB2D(I,J)*SH_URB+(1-FRC_URB2D(I,J))*SHEAT ![W/m/m]
4802 QFX(I,J) = FRC_URB2D(I,J)*LH_KINEMATIC_URB &
4803 + (1-FRC_URB2D(I,J))*ETA_KINEMATIC ![kg/m/m/s]
4804 LH(I,J) = FRC_URB2D(I,J)*LH_URB+(1-FRC_URB2D(I,J))*ETA ![W/m/m]
4805 GRDFLX(I,J) = FRC_URB2D(I,J)*G_URB+(1-FRC_URB2D(I,J))*SSOIL ![W/m/m]
4806 TSK(I,J) = FRC_URB2D(I,J)*TS_URB+(1-FRC_URB2D(I,J))*T1 ![K]
4807 Q1 = FRC_URB2D(I,J)*QS_URB+(1-FRC_URB2D(I,J))*Q1 ![-]
4808 ! Convert QSFC back to mixing ratio
4809 QSFC(I,J)= Q1/(1.0-Q1)
4810 UST(I,J)= FRC_URB2D(I,J)*UST_URB+(1-FRC_URB2D(I,J))*UST(I,J) ![m/s]
4815 print*, ' FRC_URB2D', FRC_URB2D, &
4816 'ALB_URB',ALB_URB, 'ALBEDOK',ALBEDOK, &
4817 'ALBEDO(I,J)', ALBEDO(I,J), &
4818 'SH_URB',SH_URB,'SHEAT',SHEAT, 'HFX(I,J)',HFX(I,J), &
4819 'LH_KINEMATIC_URB',LH_KINEMATIC_URB,'ETA_KINEMATIC', &
4820 ETA_KINEMATIC, 'QFX(I,J)',QFX(I,J), &
4821 'LH_URB',LH_URB, 'ETA',ETA, 'LH(I,J)',LH(I,J), &
4822 'G_URB',G_URB,'SSOIL',SSOIL,'GRDFLX(I,J)', GRDFLX(I,J),&
4823 'TS_URB',TS_URB,'T1',T1,'TSK(I,J)',TSK(I,J), &
4824 'QS_URB',QS_URB,'Q1',Q1,'QSFC(I,J)',QSFC(I,J)
4828 ! Renew Urban State Varialbes
4830 TR_URB2D(I,J) = TR_URB
4831 TB_URB2D(I,J) = TB_URB
4832 TG_URB2D(I,J) = TG_URB
4833 TC_URB2D(I,J) = TC_URB
4834 QC_URB2D(I,J) = QC_URB
4835 UC_URB2D(I,J) = UC_URB
4837 DO K = 1,num_roof_layers
4838 TRL_URB3D(I,K,J) = TRL_URB(K)
4839 SMR_URB3D(I,K,J) = SMR_URB(K)
4840 TGRL_URB3D(I,K,J)= TGRL_URB(K)
4842 DO K = 1,num_wall_layers
4843 TBL_URB3D(I,K,J) = TBL_URB(K)
4845 DO K = 1,num_road_layers
4846 TGL_URB3D(I,K,J) = TGL_URB(K)
4849 TGR_URB2D(I,J) =TGR_URB
4850 CMCR_URB2D(I,J)=CMCR_URB
4851 FLXHUMR_URB2D(I,J)=FLXHUMR_URB
4852 FLXHUMB_URB2D(I,J)=FLXHUMB_URB
4853 FLXHUMG_URB2D(I,J)=FLXHUMG_URB
4854 DRELR_URB2D(I,J) = DRELR_URB
4855 DRELB_URB2D(I,J) = DRELB_URB
4856 DRELG_URB2D(I,J) = DRELG_URB
4858 XXXR_URB2D(I,J) = XXXR_URB
4859 XXXB_URB2D(I,J) = XXXB_URB
4860 XXXG_URB2D(I,J) = XXXG_URB
4861 XXXC_URB2D(I,J) = XXXC_URB
4863 SH_URB2D(I,J) = SH_URB
4864 LH_URB2D(I,J) = LH_URB
4865 G_URB2D(I,J) = G_URB
4866 RN_URB2D(I,J) = RN_URB
4867 PSIM_URB2D(I,J) = PSIM_URB
4868 PSIH_URB2D(I,J) = PSIH_URB
4869 GZ1OZ0_URB2D(I,J)= GZ1OZ0_URB
4870 U10_URB2D(I,J) = U10_URB
4871 V10_URB2D(I,J) = V10_URB
4872 TH2_URB2D(I,J) = TH2_URB
4873 Q2_URB2D(I,J) = Q2_URB
4874 UST_URB2D(I,J) = UST_URB
4875 AKMS_URB2D(I,J) = KARMAN * UST_URB2D(I,J)/(GZ1OZ0_URB2D(I,J)-PSIM_URB2D(I,J))
4876 IF (PRESENT(CMR_SFCDIF)) THEN
4877 CMR_SFCDIF(I,J) = CMR_URB
4878 CHR_SFCDIF(I,J) = CHR_URB
4879 CMGR_SFCDIF(I,J) = CMGR_URB
4880 CHGR_SFCDIF(I,J) = CHGR_URB
4881 CMC_SFCDIF(I,J) = CMC_URB
4882 CHC_SFCDIF(I,J) = CHC_URB
4886 ENDIF ! end of UCM CALL if block
4887 !--------------------------------------
4888 ! Urban Part End - urban
4889 !--------------------------------------
4893 SMSTOT(I,J)=SOILM*1000.
4895 SMCREL(I,NS,J)=SMAV(NS)
4898 ! Convert the water unit into mm
4899 SFCRUNOFF(I,J)=SFCRUNOFF(I,J)+RUNOFF1*DT*1000.0
4900 UDRUNOFF(I,J)=UDRUNOFF(I,J)+RUNOFF2*DT*1000.0
4901 ! snow defined when fraction of frozen precip (FFROZP) > 0.5,
4902 IF(FFROZP.GT.0.5)THEN
4903 ACSNOW(I,J)=ACSNOW(I,J)+PRCP*DT
4905 IF(SNOW(I,J).GT.0.)THEN
4906 ACSNOM(I,J)=ACSNOM(I,J)+SNOMLT*1000.
4907 ! accumulated snow-melt energy
4908 SNOPCX(I,J)=SNOPCX(I,J)-SNOMLT/FDTLIW
4911 ENDIF ! endif of land-sea test
4913 ENDIF ! ENDIF FOR MOSAIC DANLI ! This corresponds to IF ((sf_surface_mosaic == 1) .AND. ((XLAND(I,J)-1.5).LT.0.) .AND. (XICE(I,J) < XICE_THRESHOLD) ) THEN
4915 ENDDO ILOOP ! of I loop
4916 ENDDO JLOOP ! of J loop
4918 !------------------------------------------------------
4919 END SUBROUTINE lsm_mosaic
4920 !------------------------------------------------------
4921 !===========================================================================
4923 ! subroutine lsm_mosaic_init: initialization of mosaic state variables
4925 !===========================================================================
4927 SUBROUTINE lsm_mosaic_init(IVGTYP,ISWATER,ISURBAN,ISICE, XLAND, XICE,fractional_seaice, &
4928 TSK,TSLB,SMOIS,SH2O,SNOW,SNOWC,SNOWH,CANWAT, &
4929 ids,ide, jds,jde, kds,kde, &
4930 ims,ime, jms,jme, kms,kme, &
4931 its,ite, jts,jte, kts,kte, restart, &
4932 landusef,landusef2,NLCAT,num_soil_layers &
4933 ,sf_surface_mosaic, mosaic_cat &
4935 ,TSK_mosaic,TSLB_mosaic &
4936 ,SMOIS_mosaic,SH2O_mosaic &
4937 ,CANWAT_mosaic,SNOW_mosaic &
4938 ,SNOWH_mosaic,SNOWC_mosaic &
4939 ,ALBEDO,ALBBCK, EMISS, EMBCK,Z0 & !danli
4940 ,ALBEDO_mosaic,ALBBCK_mosaic, EMISS_mosaic & !danli
4941 ,EMBCK_mosaic, ZNT_mosaic, Z0_mosaic & !danli
4942 ,TR_URB2D_mosaic,TB_URB2D_mosaic & !danli mosaic
4943 ,TG_URB2D_mosaic,TC_URB2D_mosaic & !danli mosaic
4944 ,QC_URB2D_mosaic & !danli mosaic
4945 ,TRL_URB3D_mosaic,TBL_URB3D_mosaic & !danli mosaic
4946 ,TGL_URB3D_mosaic & !danli mosaic
4947 ,SH_URB2D_mosaic,LH_URB2D_mosaic & !danli mosaic
4948 ,G_URB2D_mosaic,RN_URB2D_mosaic & !danli mosaic
4949 ,TS_URB2D_mosaic & !danli mosaic
4950 ,TS_RUL2D_mosaic & !danli mosaic
4953 INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
4954 ims,ime, jms,jme, kms,kme, &
4955 its,ite, jts,jte, kts,kte
4957 INTEGER, INTENT(IN) :: NLCAT, num_soil_layers, ISWATER,ISURBAN, ISICE, fractional_seaice
4959 LOGICAL , INTENT(IN) :: restart
4961 ! REAL, DIMENSION( num_soil_layers), INTENT(INOUT) :: ZS, DZS
4963 REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , &
4964 INTENT(IN) :: SMOIS, & !Total soil moisture
4965 SH2O, & !liquid soil moisture
4968 REAL, DIMENSION( ims:ime, jms:jme ) , &
4969 INTENT(IN) :: SNOW, &
4975 INTEGER, INTENT(IN) :: sf_surface_mosaic
4976 INTEGER, INTENT(IN) :: mosaic_cat
4977 INTEGER, DIMENSION( ims:ime, jms:jme ),INTENT(IN) :: IVGTYP
4978 REAL, DIMENSION( ims:ime, NLCAT, jms:jme ) , INTENT(IN):: LANDUSEF
4979 REAL, DIMENSION( ims:ime, NLCAT, jms:jme ) , INTENT(INOUT):: LANDUSEF2
4981 INTEGER, DIMENSION( ims:ime, NLCAT, jms:jme ), INTENT(INOUT) :: mosaic_cat_index
4983 REAL, DIMENSION( ims:ime, 1:mosaic_cat, jms:jme ) , OPTIONAL, INTENT(INOUT):: &
4984 TSK_mosaic, CANWAT_mosaic, SNOW_mosaic,SNOWH_mosaic, SNOWC_mosaic
4985 REAL, DIMENSION( ims:ime, 1:num_soil_layers*mosaic_cat, jms:jme ), OPTIONAL, INTENT(INOUT):: &
4986 TSLB_mosaic,SMOIS_mosaic,SH2O_mosaic
4988 REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(IN):: ALBEDO, ALBBCK, EMISS, EMBCK, Z0
4989 REAL, DIMENSION( ims:ime, 1:mosaic_cat, jms:jme ) , OPTIONAL, INTENT(INOUT):: &
4990 ALBEDO_mosaic,ALBBCK_mosaic, EMISS_mosaic, EMBCK_mosaic, ZNT_mosaic, Z0_mosaic
4992 REAL, DIMENSION( ims:ime, 1:mosaic_cat, jms:jme ) , OPTIONAL, INTENT(INOUT):: &
4993 TR_URB2D_mosaic, TB_URB2D_mosaic, TG_URB2D_mosaic, TC_URB2D_mosaic,QC_URB2D_mosaic, &
4994 SH_URB2D_mosaic,LH_URB2D_mosaic,G_URB2D_mosaic,RN_URB2D_mosaic,TS_URB2D_mosaic, TS_RUL2D_mosaic
4996 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_soil_layers*mosaic_cat, jms:jme ), INTENT(INOUT) :: TRL_URB3D_mosaic
4997 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_soil_layers*mosaic_cat, jms:jme ), INTENT(INOUT) :: TBL_URB3D_mosaic
4998 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_soil_layers*mosaic_cat, jms:jme ), INTENT(INOUT) :: TGL_URB3D_mosaic
5000 INTEGER :: ij,i,j,mosaic_i,LastSwap,NumPairs,soil_k, Temp2,Temp5,Temp7, ICE,temp_index
5001 REAL :: Temp, Temp3,Temp4,Temp6,xice_threshold
5003 CHARACTER(len=256) :: message_text
5007 if ( fractional_seaice == 0 ) then
5008 xice_threshold = 0.5
5009 else if ( fractional_seaice == 1 ) then
5010 xice_threshold = 0.02
5013 IF(.not.restart)THEN
5014 !===========================================================================
5016 !===========================================================================
5025 IF ((xland(i,j).LT. 1.5 ) .AND. (IVGTYP(i,j) .EQ. ISWATER)) THEN
5026 PRINT*, 'BEFORE MOSAIC_INIT'
5027 CALL wrf_message("BEFORE MOSAIC_INIT")
5028 WRITE(message_text,fmt='(a,2I6,2F8.2,2I6)') 'I,J,xland,xice,mosaic_cat_index,ivgtyp = ', &
5029 I,J,xland(i,j),xice(i,j),mosaic_cat_index(I,1,J),IVGTYP(i,j)
5030 CALL wrf_message(message_text)
5038 LANDUSEF2(i,mosaic_i,j)=LANDUSEF(i,mosaic_i,j)
5039 mosaic_cat_index(i,mosaic_i,j)=mosaic_i
5050 IF (NumPairs == 0) EXIT
5052 DO mosaic_i=1, NumPairs
5053 IF(LANDUSEF2(i,mosaic_i, j) < LANDUSEF2(i,mosaic_i+1, j) ) THEN
5054 Temp = LANDUSEF2(i,mosaic_i, j)
5055 LANDUSEF2(i,mosaic_i, j)=LANDUSEF2(i,mosaic_i+1, j)
5056 LANDUSEF2(i,mosaic_i+1, j)=Temp
5059 Temp2 = mosaic_cat_index(i,mosaic_i,j)
5060 mosaic_cat_index(i,mosaic_i,j)=mosaic_cat_index(i,mosaic_i+1,j)
5061 mosaic_cat_index(i,mosaic_i+1,j)=Temp2
5064 NumPairs = LastSwap - 1
5070 !===========================================================================
5071 ! For non-seaice grids, eliminate the seaice-tiles
5072 !===========================================================================
5077 IF (XLAND(I,J).LT.1.5) THEN
5080 IF( XICE(I,J).GE. XICE_THRESHOLD ) THEN
5081 WRITE (message_text,fmt='(a,2I5)') 'sea-ice at point, I and J = ', i,j
5082 CALL wrf_message(message_text)
5086 IF (ICE == 1) Then ! sea-ice case , eliminate sea-ice if they are not the dominant ones
5088 IF (IVGTYP(i,j) == isice) THEN ! if this grid cell is dominanted by ice, then do nothing
5092 DO mosaic_i=2,mosaic_cat
5093 IF (mosaic_cat_index(i,mosaic_i,j) == isice ) THEN
5094 Temp4=LANDUSEF2(i,mosaic_i,j)
5095 Temp5=mosaic_cat_index(i,mosaic_i,j)
5097 LANDUSEF2(i,mosaic_i:NLCAT-1,j)=LANDUSEF2(i,mosaic_i+1:NLCAT,j)
5098 mosaic_cat_index(i,mosaic_i:NLCAT-1,j)=mosaic_cat_index(i,mosaic_i+1:NLCAT,j)
5100 LANDUSEF2(i,NLCAT,j)=Temp4
5101 mosaic_cat_index(i,NLCAT,j)=Temp5
5105 ENDIF ! for (IVGTYP(i,j) == isice )
5107 ELSEIF (ICE ==0) THEN
5109 IF ((mosaic_cat_index(I,1,J) .EQ. ISWATER)) THEN
5111 ! xland < 1.5 but the dominant land use category based on our calculation is water
5113 IF (IVGTYP(i,j) .EQ. ISWATER) THEN
5115 ! xland < 1.5 but the dominant land use category based on the geogrid calculation is water, this must be wrong
5117 CALL wrf_message("IN MOSAIC_INIT")
5118 WRITE(message_text,fmt='(a,3I6,2F8.2)') 'I,J,IVGTYP,XLAND,XICE = ',I,J,IVGTYP(I,J),xland(i,j),xice(i,j)
5119 CALL wrf_message(message_text)
5120 CALL wrf_message("xland < 1.5 but the dominant land use category based on our calculation is water."//&
5121 "In addition, the dominant land use category based on the geogrid calculation is water, this must be wrong")
5123 ENDIF ! for (IVGTYP(i,j) .EQ. ISWATER)
5125 IF (IVGTYP(i,j) .NE. ISWATER) THEN
5127 ! xland < 1.5, the dominant land use category based on our calculation is water, but based on the geogrid calculation is not water, which might be due to the inconsistence between land use data and land-sea mask
5129 Temp4=LANDUSEF2(i,1,j)
5130 Temp5=mosaic_cat_index(i,1,j)
5132 LANDUSEF2(i,1:NLCAT-1,j)=LANDUSEF2(i,2:NLCAT,j)
5133 mosaic_cat_index(i,1:NLCAT-1,j)=mosaic_cat_index(i,2:NLCAT,j)
5135 LANDUSEF2(i,NLCAT,j)=Temp4
5136 mosaic_cat_index(i,NLCAT,j)=Temp5
5138 CALL wrf_message("IN MOSAIC_INIT")
5139 WRITE(message_text,fmt='(a,3I6,2F8.2)') 'I,J,IVGTYP,XLAND,XICE = ',I,J,IVGTYP(I,J),xland(i,j),xice(i,j)
5140 CALL wrf_message(message_text)
5141 CALL wrf_message("xland < 1.5 but the dominant land use category based on our calculation is water."//&
5142 "this is fine as long as we change our calculation so that the dominant land use category is"//&
5143 "stwiched back to not water.")
5144 WRITE(message_text,fmt='(a,2I6)') 'land use category has been switched, before and after values are ', &
5145 temp5,mosaic_cat_index(i,1,j)
5146 CALL wrf_message(message_text)
5147 WRITE(message_text,fmt='(a,2I6)') 'new dominant and second dominant cat are ', mosaic_cat_index(i,1,j),mosaic_cat_index(i,2,j)
5148 CALL wrf_message(message_text)
5150 ENDIF ! for (IVGTYP(i,j) .NE. ISWATER)
5152 ELSE ! for (mosaic_cat_index(I,1,J) .EQ. ISWATER)
5154 DO mosaic_i=2,mosaic_cat
5155 IF (mosaic_cat_index(i,mosaic_i,j) == iswater ) THEN
5156 Temp4=LANDUSEF2(i,mosaic_i,j)
5157 Temp5=mosaic_cat_index(i,mosaic_i,j)
5159 LANDUSEF2(i,mosaic_i:NLCAT-1,j)=LANDUSEF2(i,mosaic_i+1:NLCAT,j)
5160 mosaic_cat_index(i,mosaic_i:NLCAT-1,j)=mosaic_cat_index(i,mosaic_i+1:NLCAT,j)
5162 LANDUSEF2(i,NLCAT,j)=Temp4
5163 mosaic_cat_index(i,NLCAT,j)=Temp5
5167 ENDIF ! for (mosaic_cat_index(I,1,J) .EQ. ISWATER)
5169 ENDIF ! for ICE == 1
5171 ELSE ! FOR (XLAND(I,J).LT.1.5)
5175 IF( XICE(I,J).GE. XICE_THRESHOLD ) THEN
5176 WRITE (message_text,fmt='(a,2I6)') 'sea-ice at water point, I and J = ', i,j
5177 CALL wrf_message(message_text)
5181 IF ((mosaic_cat_index(I,1,J) .NE. ISWATER)) THEN
5183 ! xland > 1.5 and the dominant land use category based on our calculation is not water
5185 IF (IVGTYP(i,j) .NE. ISWATER) THEN
5187 ! xland > 1.5 but the dominant land use category based on the geogrid calculation is not water, this must be wrong
5188 CALL wrf_message("IN MOSAIC_INIT")
5189 WRITE(message_text,fmt='(a,3I6,2F8.2)') 'I,J,IVGTYP,XLAND,XICE = ',I,J,IVGTYP(I,J),xland(i,j),xice(i,j)
5190 CALL wrf_message(message_text)
5191 CALL wrf_message("xland > 1.5 but the dominant land use category based on our calculation is not water."// &
5192 "in addition, the dominant land use category based on the geogrid calculation is not water,"// &
5193 "this must be wrong.")
5194 ENDIF ! for (IVGTYP(i,j) .NE. ISWATER)
5196 IF (IVGTYP(i,j) .EQ. ISWATER) THEN
5198 ! xland > 1.5, the dominant land use category based on our calculation is not water, but based on the geogrid calculation is water, which might be due to the inconsistence between land use data and land-sea mask
5200 CALL wrf_message("IN MOSAIC_INIT")
5201 WRITE(message_text,fmt='(a,3I6,2F8.2)') 'I,J,IVGTYP,XLAND,XICE = ',I,J,IVGTYP(I,J),xland(i,j),xice(i,j)
5202 CALL wrf_message(message_text)
5203 CALL wrf_message("xland > 1.5 but the dominant land use category based on our calculation is not water."// &
5204 "however, the dominant land use category based on the geogrid calculation is water")
5205 CALL wrf_message("This is fine. We do not need to do anyting because in the noaddrv, "//&
5206 "we use xland as a criterion for whether using"// &
5207 "mosaic or not when xland > 1.5, no mosaic will be used anyway")
5209 ENDIF ! for (IVGTYP(i,j) .NE. ISWATER)
5211 ENDIF ! for (mosaic_cat_index(I,1,J) .NE. ISWATER)
5213 ENDIF ! FOR (XLAND(I,J).LT.1.5)
5218 !===========================================================================
5220 !===========================================================================
5227 DO mosaic_i=1,mosaic_cat
5228 Temp6=Temp6+LANDUSEF2(i,mosaic_i,j)
5231 if (Temp6 .LT. 1e-5) then
5234 WRITE (message_text,fmt='(a,e8.1)') 'the total land surface fraction is less than ', temp6
5235 CALL wrf_message(message_text)
5236 WRITE (message_text,fmt='(a,2I6,4F8.2)') 'some landusef values at i,j are ', &
5237 i,j,landusef2(i,1,j),landusef2(i,2,j),landusef2(i,3,j),landusef2(i,4,j)
5238 CALL wrf_message(message_text)
5239 WRITE (message_text,fmt='(a,2I6,3I6)') 'some mosaic cat values at i,j are ', &
5240 i,j,mosaic_cat_index(i,1,j),mosaic_cat_index(i,2,j),mosaic_cat_index(i,3,j)
5241 CALL wrf_message(message_text)
5245 LANDUSEF2(i,1:mosaic_cat, j)=LANDUSEF2(i,1:mosaic_cat,j)*(1/Temp6)
5250 !===========================================================================
5251 ! initilize the variables
5252 !===========================================================================
5257 DO mosaic_i=1,mosaic_cat
5259 TSK_mosaic(i,mosaic_i,j)=TSK(i,j)
5260 CANWAT_mosaic(i,mosaic_i,j)=CANWAT(i,j)
5261 SNOW_mosaic(i,mosaic_i,j)=SNOW(i,j)
5262 SNOWH_mosaic(i,mosaic_i,j)=SNOWH(i,j)
5263 SNOWC_mosaic(i,mosaic_i,j)=SNOWC(i,j)
5265 ALBEDO_mosaic(i,mosaic_i,j)=ALBEDO(i,j)
5266 ALBBCK_mosaic(i,mosaic_i,j)=ALBBCK(i,j)
5267 EMISS_mosaic(i,mosaic_i,j)=EMISS(i,j)
5268 EMBCK_mosaic(i,mosaic_i,j)=EMBCK(i,j)
5269 ZNT_mosaic(i,mosaic_i,j)=Z0(i,j)
5270 Z0_mosaic(i,mosaic_i,j)=Z0(i,j)
5272 DO soil_k=1,num_soil_layers
5274 TSLB_mosaic(i,num_soil_layers*(mosaic_i-1)+soil_k,j)=TSLB(i,soil_k,j)
5275 SMOIS_mosaic(i,num_soil_layers*(mosaic_i-1)+soil_k,j)=SMOIS(i,soil_k,j)
5276 SH2O_mosaic(i,num_soil_layers*(mosaic_i-1)+soil_k,j)=SH2O(i,soil_k,j)
5280 TR_URB2D_mosaic(i,mosaic_i,j)=TSK(i,j)
5281 TB_URB2D_mosaic(i,mosaic_i,j)=TSK(i,j)
5282 TG_URB2D_mosaic(i,mosaic_i,j)=TSK(i,j)
5283 TC_URB2D_mosaic(i,mosaic_i,j)=TSK(i,j)
5284 TS_URB2D_mosaic(i,mosaic_i,j)=TSK(i,j)
5285 TS_RUL2D_mosaic(i,mosaic_i,j)=TSK(i,j)
5286 QC_URB2D_mosaic(i,mosaic_i,j)=0.01
5287 SH_URB2D_mosaic(i,mosaic_i,j)=0
5288 LH_URB2D_mosaic(i,mosaic_i,j)=0
5289 G_URB2D_mosaic(i,mosaic_i,j)=0
5290 RN_URB2D_mosaic(i,mosaic_i,j)=0
5292 TRL_URB3D_mosaic(I,4*(mosaic_i-1)+1,J)=TSLB(I,1,J)+0.
5293 TRL_URB3D_mosaic(I,4*(mosaic_i-1)+2,J)=0.5*(TSLB(I,1,J)+TSLB(I,2,J))
5294 TRL_URB3D_mosaic(I,4*(mosaic_i-1)+3,J)=TSLB(I,2,J)+0.
5295 TRL_URB3D_mosaic(I,4*(mosaic_i-1)+4,J)=TSLB(I,2,J)+(TSLB(I,3,J)-TSLB(I,2,J))*0.29
5297 TBL_URB3D_mosaic(I,4*(mosaic_i-1)+1,J)=TSLB(I,1,J)+0.
5298 TBL_URB3D_mosaic(I,4*(mosaic_i-1)+2,J)=0.5*(TSLB(I,1,J)+TSLB(I,2,J))
5299 TBL_URB3D_mosaic(I,4*(mosaic_i-1)+3,J)=TSLB(I,2,J)+0.
5300 TBL_URB3D_mosaic(I,4*(mosaic_i-1)+4,J)=TSLB(I,2,J)+(TSLB(I,3,J)-TSLB(I,2,J))*0.29
5302 TGL_URB3D_mosaic(I,4*(mosaic_i-1)+1,J)=TSLB(I,1,J)
5303 TGL_URB3D_mosaic(I,4*(mosaic_i-1)+2,J)=TSLB(I,2,J)
5304 TGL_URB3D_mosaic(I,4*(mosaic_i-1)+3,J)=TSLB(I,3,J)
5305 TGL_URB3D_mosaic(I,4*(mosaic_i-1)+4,J)=TSLB(I,4,J)
5316 IF ((xland(i,j).LT. 1.5 ) .AND. (mosaic_cat_index(I,1,J) .EQ. ISWATER)) THEN
5317 CALL wrf_message("After MOSAIC_INIT")
5318 WRITE (message_text,fmt='(a,2I6,2F8.2,2I6)') 'weird xland,xice,mosaic_cat_index and ivgtyp at I,J = ', &
5319 i,j,xland(i,j),xice(i,j),mosaic_cat_index(I,1,J),IVGTYP(i,j)
5320 CALL wrf_message(message_text)
5326 ENDIF ! for not restart
5328 !--------------------------------
5329 END SUBROUTINE lsm_mosaic_init
5330 !--------------------------------
5333 END MODULE module_sf_noahdrv