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 lf_urb2d_s, z0_urb2d, & !SLUCM
114 th_phy,rho,p_phy,ust, & !I multi-layer urban
115 gmt,julday,xlong,xlat, & !I multi-layer urban
116 a_u_bep,a_v_bep,a_t_bep,a_q_bep, & !O multi-layer urban
117 a_e_bep,b_u_bep,b_v_bep, & !O multi-layer urban
118 b_t_bep,b_q_bep,b_e_bep,dlg_bep, & !O multi-layer urban
119 dl_u_bep,sf_bep,vl_bep &
121 ,sfcheadrt,INFXSRT,soldrain & !O multi-layer urban
123 ,SDA_HFX, SDA_QFX, HFX_BOTH, QFX_BOTH, QNORM, fasdas & !fasdas
128 !----------------------------------------------------------------
130 !----------------------------------------------------------------
131 !----------------------------------------------------------------
132 ! --- atmospheric (WRF generic) variables
133 !-- DT time step (seconds)
134 !-- DZ8W thickness of layers (m)
135 !-- T3D temperature (K)
136 !-- QV3D 3D water vapor mixing ratio (Kg/Kg)
137 !-- P3D 3D pressure (Pa)
138 !-- FLHC exchange coefficient for heat (m/s)
139 !-- FLQC exchange coefficient for moisture (m/s)
140 !-- PSFC surface pressure (Pa)
141 !-- XLAND land mask (1 for land, 2 for water)
142 !-- QGH saturated mixing ratio at 2 meter
143 !-- GSW downward short wave flux at ground surface (W/m^2)
144 !-- GLW downward long wave flux at ground surface (W/m^2)
145 !-- History variables
146 !-- CANWAT canopy moisture content (mm)
147 !-- TSK surface temperature (K)
148 !-- TSLB soil temp (k)
149 !-- SMOIS total soil moisture content (volumetric fraction)
150 !-- SH2O unfrozen soil moisture content (volumetric fraction)
151 ! note: frozen soil moisture (i.e., soil ice) = SMOIS - SH2O
152 !-- SNOWH actual snow depth (m)
153 !-- SNOW liquid water-equivalent snow depth (m)
154 !-- ALBEDO time-varying surface albedo including snow effect (unitless fraction)
155 !-- ALBBCK background surface albedo (unitless fraction)
156 !-- CHS surface exchange coefficient for heat and moisture (m s-1);
157 !-- CHS2 2m surface exchange coefficient for heat (m s-1);
158 !-- CQS2 2m surface exchange coefficient for moisture (m s-1);
160 !-- num_soil_layers the number of soil layers
161 !-- ZS depths of centers of soil layers (m)
162 !-- DZS thicknesses of soil layers (m)
163 !-- SLDPTH thickness of each soil layer (m, same as DZS)
164 !-- TMN soil temperature at lower boundary (K)
165 !-- SMCWLT wilting point (volumetric)
166 !-- SMCDRY dry soil moisture threshold where direct evap from
167 ! top soil layer ends (volumetric)
168 !-- SMCREF soil moisture threshold below which transpiration begins to
169 ! stress (volumetric)
170 !-- SMCMAX porosity, i.e. saturated value of soil moisture (volumetric)
171 !-- NROOT number of root layers, a function of veg type, determined
172 ! in subroutine redprm.
173 !-- SMSTAV Soil moisture availability for evapotranspiration (
174 ! fraction between SMCWLT and SMCMXA)
175 !-- SMSTOT Total soil moisture content frozen+unfrozen) in the soil column (mm)
177 !-- SNOWC fraction snow coverage (0-1.0)
178 ! --- vegetation variables
179 !-- SNOALB upper bound on maximum albedo over deep snow
180 !-- SHDMIN minimum areal fractional coverage of annual green vegetation
181 !-- SHDMAX maximum areal fractional coverage of annual green vegetation
182 !-- XLAI leaf area index (dimensionless)
183 !-- XLAI2 leaf area index (same as XLAI) passed to output (dimensionless)
184 !-- Z0BRD Background fixed roughness length (M)
185 !-- Z0 Background vroughness length (M) as function
186 !-- ZNT Time varying roughness length (M) as function
187 !-- ALBD(IVGTPK,ISN) background albedo reading from a table
189 !-- HFX upward heat flux at the surface (W/m^2)
190 !-- QFX upward moisture flux at the surface (kg/m^2/s)
191 !-- LH upward moisture flux at the surface (W m-2)
192 !-- GRDFLX(I,J) ground heat flux (W m-2)
193 !-- FDOWN radiation forcing at the surface (W m-2) = SOLDN*(1-alb)+LWDN
194 !----------------------------------------------------------------------------
195 !-- EC canopy water evaporation ((W m-2)
196 !-- EDIR direct soil evaporation (W m-2)
197 !-- ET plant transpiration from a particular root layer (W m-2)
198 !-- ETT total plant transpiration (W m-2)
199 !-- ESNOW sublimation from (or deposition to if <0) snowpack (W m-2)
200 !-- DRIP through-fall of precip and/or dew in excess of canopy
201 ! water-holding capacity (m)
202 !-- DEW dewfall (or frostfall for t<273.15) (M)
203 !-- SMAV Soil Moisture Availability for each layer, as a fraction
204 ! between SMCWLT and SMCMAX (dimensionless fraction)
205 ! ----------------------------------------------------------------------
206 !-- BETA ratio of actual/potential evap (dimensionless)
207 !-- ETP potential evaporation (W m-2)
208 ! ----------------------------------------------------------------------
209 !-- FLX1 precip-snow sfc (W m-2)
210 !-- FLX2 freezing rain latent heat flux (W m-2)
211 !-- FLX3 phase-change heat flux from snowmelt (W m-2)
212 ! ----------------------------------------------------------------------
213 !-- ACSNOM snow melt (mm) (water equivalent)
214 !-- ACSNOW accumulated snow fall (mm) (water equivalent)
215 !-- SNOPCX snow phase change heat flux (W/m^2)
216 !-- POTEVP accumulated potential evaporation (m)
217 !-- RIB Documentation needed!!!
218 ! ----------------------------------------------------------------------
219 !-- RUNOFF1 surface runoff (m s-1), not infiltrating the surface
220 !-- RUNOFF2 subsurface runoff (m s-1), drainage out bottom of last
221 ! soil layer (baseflow)
222 ! important note: here RUNOFF2 is actually the sum of RUNOFF2 and RUNOFF3
223 !-- RUNOFF3 numerical trunctation in excess of porosity (smcmax)
224 ! for a given soil layer at the end of a time step (m s-1).
225 !SFCRUNOFF Surface Runoff (mm)
226 !UDRUNOFF Total Underground Runoff (mm), which is the sum of RUNOFF2 and RUNOFF3
227 ! ----------------------------------------------------------------------
228 !-- RC canopy resistance (s m-1)
229 !-- RC2 canopy resistance (same as RC) passed to output
230 !-- PC plant coefficient (unitless fraction, 0-1) where PC*ETP = actual transp
231 !-- RSMIN minimum canopy resistance (s m-1)
232 !-- RCS incoming solar rc factor (dimensionless)
233 !-- RCT air temperature rc factor (dimensionless)
234 !-- RCQ atmos vapor pressure deficit rc factor (dimensionless)
235 !-- RCSOIL soil moisture rc factor (dimensionless)
237 !-- EMISS surface emissivity (between 0 and 1)
238 !-- EMBCK Background surface emissivity (between 0 and 1)
241 ! (R_d/R_v) (dimensionless)
242 !-- ids start index for i in domain
243 !-- ide end index for i in domain
244 !-- jds start index for j in domain
245 !-- jde end index for j in domain
246 !-- kds start index for k in domain
247 !-- kde end index for k in domain
248 !-- ims start index for i in memory
249 !-- ime end index for i in memory
250 !-- jms start index for j in memory
251 !-- jme end index for j in memory
252 !-- kms start index for k in memory
253 !-- kme end index for k in memory
254 !-- its start index for i in tile
255 !-- ite end index for i in tile
256 !-- jts start index for j in tile
257 !-- jte end index for j in tile
258 !-- kts start index for k in tile
259 !-- kte end index for k in tile
261 !-- SR fraction of frozen precip (0.0 to 1.0)
262 !----------------------------------------------------------------
266 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
267 ims,ime, jms,jme, kms,kme, &
268 its,ite, jts,jte, kts,kte
270 INTEGER, INTENT(IN ) :: sf_urban_physics !urban
271 INTEGER, INTENT(IN ) :: isurban
272 INTEGER, INTENT(IN ) :: isice
273 INTEGER, INTENT(IN ) :: julian, julyr !urban
275 !added by Wei Yu for routing
277 REAL, DIMENSION( ims:ime, jms:jme ) , &
278 INTENT(INOUT) :: sfcheadrt,INFXSRT,soldrain
283 ! new local vars for hydro
284 REAL :: etpnd1_hydro,sfcheadrt_hydro,infxsrt_hydro
286 REAL, DIMENSION( ims:ime, jms:jme ) , &
287 INTENT(IN ) :: TMN, &
295 SWDOWN, & !added 10 jan 2007
303 REAL, DIMENSION( ims:ime, jms:jme ) , &
304 INTENT(INOUT) :: ALBBCK, &
306 CHARACTER(LEN=*), INTENT(IN ) :: MMINLU
308 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
309 INTENT(IN ) :: QV3D, &
313 REAL, DIMENSION( ims:ime, jms:jme ) , &
314 INTENT(IN ) :: QGH, &
317 INTEGER, DIMENSION( ims:ime, jms:jme ) , &
318 INTENT(IN ) :: IVGTYP, &
321 INTEGER, INTENT(IN) :: num_soil_layers,ITIMESTEP
323 REAL, INTENT(IN ) :: DT,ROVCP
325 REAL, DIMENSION(1:num_soil_layers), INTENT(IN)::DZS
329 REAL, DIMENSION( ims:ime , 1:num_soil_layers, jms:jme ), &
330 INTENT(INOUT) :: SMOIS, & ! total soil moisture
331 SH2O, & ! new soil liquid
334 REAL, DIMENSION( ims:ime , 1:num_soil_layers, jms:jme ), &
335 INTENT(OUT) :: SMCREL
337 REAL, DIMENSION( ims:ime, jms:jme ) , &
338 INTENT(INOUT) :: TSK, & !was TGB (temperature)
364 REAL, DIMENSION( ims:ime, jms:jme ) , &
365 INTENT(OUT) :: NOAHRES
366 INTEGER, INTENT(IN) :: OPT_THCND
369 LOGICAL, INTENT(IN) :: UA_PHYS
370 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: FLX4_2D,FVB_2D,FBUR_2D,FGSN_2D
371 REAL :: FLX4,FVB,FBUR,FGSN
373 REAL, DIMENSION( ims:ime, jms:jme ) , &
374 INTENT(OUT) :: CHKLOWQ
375 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: LAI
376 REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: QZ0
377 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: RC2, XLAI2
379 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CMR_SFCDIF
380 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CHR_SFCDIF
381 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CMGR_SFCDIF
382 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CHGR_SFCDIF
383 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CMC_SFCDIF
384 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CHC_SFCDIF
385 ! Local variables (moved here from driver to make routine thread safe, 20031007 jm)
387 REAL, DIMENSION(1:num_soil_layers) :: ET
389 REAL, DIMENSION(1:num_soil_layers) :: SMAV
391 REAL :: BETA, ETP, SSOIL,EC, EDIR, ESNOW, ETT, &
392 FLX1,FLX2,FLX3, DRIP,DEW,FDOWN,RC,PC,RSMIN,XLAI, &
394 RCS,RCT,RCQ,RCSOIL,FFROZP
396 LOGICAL, INTENT(IN ) :: myj,frpcpn
398 ! DECLARATIONS - LOGICAL
399 ! ----------------------------------------------------------------------
400 LOGICAL, PARAMETER :: LOCAL=.false.
401 LOGICAL :: FRZGRA, SNOWNG
405 ! ----------------------------------------------------------------------
406 ! DECLARATIONS - INTEGER
407 ! ----------------------------------------------------------------------
408 INTEGER :: I,J, ICE,NSOIL,SLOPETYP,SOILTYP,VEGTYP
412 ! ----------------------------------------------------------------------
413 ! DECLARATIONS - REAL
414 ! ----------------------------------------------------------------------
416 REAL :: SHMIN,SHMAX,DQSDT2,LWDN,PRCP,PRCPRAIN, &
417 Q2SAT,Q2SATI,SFCPRS,SFCSPD,SFCTMP,SHDFAC,SNOALB1, &
418 SOLDN,TBOT,ZLVL, Q2K,ALBBRD, ALBEDOK, ETA, ETA_KINEMATIC, &
420 Z0K,RUNOFF1,RUNOFF2,RUNOFF3,SHEAT,SOLNET,E2SAT,SFCTSNO, &
421 ! mek, WRF testing, expanded diagnostics
422 SOLUP,LWUP,RNET,RES,Q1SFC,TAIRV,SATFLG
425 ! MEK JUL2007 for pot. evap.
431 REAL :: SNCOVR,SNEQV,SNOWHK,CMC, CHK,TH2
433 REAL :: SMCDRY,SMCMAX,SMCREF,SMCWLT,SNOMLT,SOILM,SOILW,Q1,T1
434 REAL :: SNOTIME1 ! LSTSNW1 INITIAL NUMBER OF TIMESTEPS SINCE LAST SNOWFALL
438 REAL :: COSZ, SOLARDIRECT
440 REAL, DIMENSION(1:num_soil_layers):: SLDPTH, STC,SMC,SWC
442 REAL, DIMENSION(1:num_soil_layers) :: ZSOIL, RTDIS
443 REAL, PARAMETER :: TRESH=.95E0, A2=17.67,A3=273.15,A4=29.65, &
444 T0=273.16E0, ELWV=2.50E6, A23M4=A2*(A3-A4)
446 REAL, PARAMETER :: ROW=1.E3,ELIW=XLF,ROWLIW=ROW*ELIW
448 ! ----------------------------------------------------------------------
449 ! DECLARATIONS START - urban
450 ! ----------------------------------------------------------------------
452 ! input variables surface_driver --> lsm
453 INTEGER, INTENT(IN) :: num_roof_layers
454 INTEGER, INTENT(IN) :: num_wall_layers
455 INTEGER, INTENT(IN) :: num_road_layers
456 REAL, OPTIONAL, DIMENSION(1:num_roof_layers), INTENT(IN) :: DZR
457 REAL, OPTIONAL, DIMENSION(1:num_wall_layers), INTENT(IN) :: DZB
458 REAL, OPTIONAL, DIMENSION(1:num_road_layers), INTENT(IN) :: DZG
459 REAL, OPTIONAL, INTENT(IN) :: DECLIN_URB
460 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: COSZ_URB2D
461 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: OMG_URB2D
462 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: XLAT_URB2D
463 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: U_PHY
464 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: V_PHY
465 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: TH_PHY
466 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: P_PHY
467 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: RHO
468 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: UST
470 LOGICAL, intent(in) :: rdlai2d
471 LOGICAL, intent(in) :: USEMONALB
473 ! input variables lsm --> urban
474 INTEGER :: UTYPE_URB ! urban type [urban=1, suburban=2, rural=3]
475 REAL :: TA_URB ! potential temp at 1st atmospheric level [K]
476 REAL :: QA_URB ! mixing ratio at 1st atmospheric level [kg/kg]
477 REAL :: UA_URB ! wind speed at 1st atmospheric level [m/s]
478 REAL :: U1_URB ! u at 1st atmospheric level [m/s]
479 REAL :: V1_URB ! v at 1st atmospheric level [m/s]
480 REAL :: SSG_URB ! downward total short wave radiation [W/m/m]
481 REAL :: LLG_URB ! downward long wave radiation [W/m/m]
482 REAL :: RAIN_URB ! precipitation [mm/h]
483 REAL :: RHOO_URB ! air density [kg/m^3]
484 REAL :: ZA_URB ! first atmospheric level [m]
485 REAL :: DELT_URB ! time step [s]
486 REAL :: SSGD_URB ! downward direct short wave radiation [W/m/m]
487 REAL :: SSGQ_URB ! downward diffuse short wave radiation [W/m/m]
488 REAL :: XLAT_URB ! latitude [deg]
489 REAL :: COSZ_URB ! cosz
490 REAL :: OMG_URB ! hour angle
491 REAL :: ZNT_URB ! roughness length [m]
502 REAL, DIMENSION(1:num_roof_layers) :: TRL_URB ! roof layer temp [K]
503 REAL, DIMENSION(1:num_wall_layers) :: TBL_URB ! wall layer temp [K]
504 REAL, DIMENSION(1:num_road_layers) :: TGL_URB ! road layer temp [K]
505 LOGICAL :: LSOLAR_URB
507 !===Yang,2014/10/08,hydrological variable for single layer UCM===
508 INTEGER :: jmonth, jday, tloc
509 INTEGER :: IRIOPTION, USOIL, DSOIL
519 REAL, DIMENSION(1:num_roof_layers) :: SMR_URB ! green roof layer moisture
520 REAL, DIMENSION(1:num_roof_layers) :: TGRL_URB ! green roof layer temp [K]
522 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: DRELR_URB2D
523 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: DRELB_URB2D
524 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: DRELG_URB2D
525 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: FLXHUMR_URB2D
526 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: FLXHUMB_URB2D
527 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: FLXHUMG_URB2D
528 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CMCR_URB2D
529 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TGR_URB2D
530 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_roof_layers, jms:jme ), INTENT(INOUT) :: TGRL_URB3D
531 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_roof_layers, jms:jme ), INTENT(INOUT) :: SMR_URB3D
534 ! state variable surface_driver <--> lsm <--> urban
535 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TR_URB2D
536 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TB_URB2D
537 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TG_URB2D
538 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TC_URB2D
539 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: QC_URB2D
540 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: UC_URB2D
541 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXR_URB2D
542 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXB_URB2D
543 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXG_URB2D
544 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXC_URB2D
545 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: SH_URB2D
546 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: LH_URB2D
547 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: G_URB2D
548 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: RN_URB2D
550 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TS_URB2D
552 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_roof_layers, jms:jme ), INTENT(INOUT) :: TRL_URB3D
553 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_wall_layers, jms:jme ), INTENT(INOUT) :: TBL_URB3D
554 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_road_layers, jms:jme ), INTENT(INOUT) :: TGL_URB3D
556 ! output variable lsm --> surface_driver
557 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: PSIM_URB2D
558 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: PSIH_URB2D
559 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: GZ1OZ0_URB2D
560 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: U10_URB2D
561 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: V10_URB2D
562 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: TH2_URB2D
563 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: Q2_URB2D
565 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: AKMS_URB2D
567 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: UST_URB2D
568 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: FRC_URB2D
569 INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: UTYPE_URB2D
572 ! output variables urban --> lsm
573 REAL :: TS_URB ! surface radiative temperature [K]
574 REAL :: QS_URB ! surface humidity [-]
575 REAL :: SH_URB ! sensible heat flux [W/m/m]
576 REAL :: LH_URB ! latent heat flux [W/m/m]
577 REAL :: LH_KINEMATIC_URB ! latent heat flux, kinetic [kg/m/m/s]
578 REAL :: SW_URB ! upward short wave radiation flux [W/m/m]
579 REAL :: ALB_URB ! time-varying albedo [fraction]
580 REAL :: LW_URB ! upward long wave radiation flux [W/m/m]
581 REAL :: G_URB ! heat flux into the ground [W/m/m]
582 REAL :: RN_URB ! net radiation [W/m/m]
583 REAL :: PSIM_URB ! shear f for momentum [-]
584 REAL :: PSIH_URB ! shear f for heat [-]
585 REAL :: GZ1OZ0_URB ! shear f for heat [-]
586 REAL :: U10_URB ! wind u component at 10 m [m/s]
587 REAL :: V10_URB ! wind v component at 10 m [m/s]
588 REAL :: TH2_URB ! potential temperature at 2 m [K]
589 REAL :: Q2_URB ! humidity at 2 m [-]
593 ! NUDAPT Parameters urban --> lam
598 REAL, DIMENSION(4) :: lf_urb
599 ! Distributed aerodynamics parameters
603 ! Variables for multi-layer UCM (Martilli et al. 2002)
604 REAL, OPTIONAL, INTENT(IN ) :: GMT
605 INTEGER, OPTIONAL, INTENT(IN ) :: JULDAY
606 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) ::XLAT, XLONG
607 INTEGER, INTENT(IN ) :: num_urban_ndm
608 INTEGER, INTENT(IN ) :: urban_map_zrd
609 INTEGER, INTENT(IN ) :: urban_map_zwd
610 INTEGER, INTENT(IN ) :: urban_map_gd
611 INTEGER, INTENT(IN ) :: urban_map_zd
612 INTEGER, INTENT(IN ) :: urban_map_zdf
613 INTEGER, INTENT(IN ) :: urban_map_bd
614 INTEGER, INTENT(IN ) :: urban_map_wd
615 INTEGER, INTENT(IN ) :: urban_map_gbd
616 INTEGER, INTENT(IN ) :: urban_map_fbd
617 INTEGER, INTENT(IN ) :: urban_map_zgrd
618 INTEGER, INTENT(IN ) :: NUM_URBAN_HI
619 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: tsk_rural_bep
620 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zrd, jms:jme ), INTENT(INOUT) :: trb_urb4d
621 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zwd, jms:jme ), INTENT(INOUT) :: tw1_urb4d
622 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zwd, jms:jme ), INTENT(INOUT) :: tw2_urb4d
623 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_gd , jms:jme ), INTENT(INOUT) :: tgb_urb4d
624 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_bd , jms:jme ), INTENT(INOUT) :: tlev_urb3d
625 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_bd , jms:jme ), INTENT(INOUT) :: qlev_urb3d
626 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_wd , jms:jme ), INTENT(INOUT) :: tw1lev_urb3d
627 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_wd , jms:jme ), INTENT(INOUT) :: tw2lev_urb3d
628 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_gbd, jms:jme ), INTENT(INOUT) :: tglev_urb3d
629 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_fbd, jms:jme ), INTENT(INOUT) :: tflev_urb3d
630 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: lf_ac_urb3d
631 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: sf_ac_urb3d
632 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: cm_ac_urb3d
633 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: sfvent_urb3d
634 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: lfvent_urb3d
635 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_wd , jms:jme ), INTENT(INOUT) :: sfwin1_urb3d
636 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_wd , jms:jme ), INTENT(INOUT) :: sfwin2_urb3d
637 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zd , jms:jme ), INTENT(INOUT) :: sfw1_urb3d
638 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zd , jms:jme ), INTENT(INOUT) :: sfw2_urb3d
639 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme ), INTENT(INOUT) :: sfr_urb3d
640 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_ndm, jms:jme ), INTENT(INOUT) :: sfg_urb3d
641 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: ep_pv_urb3d !GRZ
642 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme ), INTENT(INOUT) :: t_pv_urb3d !GRZ
643 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zgrd, jms:jme ),INTENT(INOUT) :: trv_urb4d !GRZ
644 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zgrd, jms:jme ),INTENT(INOUT) :: qr_urb4d !GRZ
645 REAL, OPTIONAL, DIMENSION( ims:ime,jms:jme ), INTENT(INOUT) :: qgr_urb3d !GRZ
646 REAL, OPTIONAL, DIMENSION( ims:ime,jms:jme ), INTENT(INOUT) :: tgr_urb3d !GRZ
647 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme ),INTENT(INOUT) :: drain_urb4d !GRZ
648 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: draingr_urb3d !GRZ
649 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme ),INTENT(INOUT) :: sfrv_urb3d !GRZ
650 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme ),INTENT(INOUT) :: lfrv_urb3d !GRZ
651 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme ),INTENT(INOUT) :: dgr_urb3d !GRZ
652 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_ndm, jms:jme ),INTENT(INOUT) :: dg_urb3d !GRZ
653 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme ),INTENT(INOUT) :: lfr_urb3d !GRZ
654 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_ndm, jms:jme ),INTENT(INOUT) :: lfg_urb3d !GRZ
656 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_hi, jms:jme ), INTENT(IN) :: hi_urb2d
657 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: lp_urb2d
658 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: lb_urb2d
659 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: hgt_urb2d
660 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: mh_urb2d
661 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: stdh_urb2d
662 REAL, OPTIONAL, DIMENSION( ims:ime, 4, jms:jme ), INTENT(IN) :: lf_urb2d
663 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: lf_urb2d_s
664 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: z0_urb2d
665 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_u_bep !Implicit momemtum component X-direction
666 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_v_bep !Implicit momemtum component Y-direction
667 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_t_bep !Implicit component pot. temperature
668 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_q_bep !Implicit momemtum component X-direction
669 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_e_bep !Implicit component TKE
670 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::b_u_bep !Explicit momentum component X-direction
671 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::b_v_bep !Explicit momentum component Y-direction
672 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::b_t_bep !Explicit component pot. temperature
673 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::b_q_bep !Implicit momemtum component Y-direction
674 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::b_e_bep !Explicit component TKE
675 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::vl_bep !Fraction air volume in grid cell
676 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::dlg_bep !Height above ground
677 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::sf_bep !Fraction air at the face of grid cell
678 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::dl_u_bep !Length scale
680 ! Local variables for multi-layer UCM (Martilli et al. 2002)
681 REAL, DIMENSION( its:ite, jts:jte ) :: HFX_RURAL,LH_RURAL,GRDFLX_RURAL ! ,RN_RURAL
682 REAL, DIMENSION( its:ite, jts:jte ) :: QFX_RURAL ! ,QSFC_RURAL,UMOM_RURAL,VMOM_RURAL
683 REAL, DIMENSION( its:ite, jts:jte ) :: ALB_RURAL,EMISS_RURAL,TSK_RURAL ! ,UST_RURAL
684 ! REAL, DIMENSION( ims:ime, jms:jme ) :: QSFC_URB
685 REAL, DIMENSION( its:ite, jts:jte ) :: HFX_URB,UMOM_URB,VMOM_URB
686 REAL, DIMENSION( its:ite, jts:jte ) :: QFX_URB
687 ! REAL, DIMENSION( ims:ime, jms:jme ) :: ALBEDO_URB,EMISS_URB,UMOM,VMOM,UST
688 REAL, DIMENSION(its:ite,jts:jte) ::EMISS_URB
689 REAL, DIMENSION(its:ite,jts:jte) :: RL_UP_URB
690 REAL, DIMENSION(its:ite,jts:jte) ::RS_ABS_URB
691 REAL, DIMENSION(its:ite,jts:jte) ::GRDFLX_URB
692 REAL :: SIGMA_SB,RL_UP_RURAL,RL_UP_TOT,RS_ABS_TOT,UMOM,VMOM
693 REAL :: CMR_URB, CHR_URB, CMC_URB, CHC_URB, CMGR_URB, CHGR_URB
694 REAL :: frc_urb,lb_urb
696 ! ----------------------------------------------------------------------
697 ! DECLARATIONS END - urban
698 ! ----------------------------------------------------------------------
700 REAL, PARAMETER :: CAPA=R_D/CP
701 REAL :: APELM,APES,SFCTH2,PSFC
702 real, intent(in) :: xice_threshold
703 character(len=80) :: message_text
707 REAL, DIMENSION( ims:ime, jms:jme ), OPTIONAL, &
708 INTENT(INOUT) :: SDA_HFX, SDA_QFX, HFX_BOTH, QFX_BOTH, QNORM
709 INTEGER, INTENT(IN ) :: fasdas
711 REAL :: XSDA_HFX, XSDA_QFX, XQNORM
712 REAL :: HFX_PHY, QFX_PHY
716 REAL,OPTIONAL,INTENT(IN),DIMENSION( ims:ime, jms:jme ) :: IRR_CHAN
717 REAL :: IRRIGATION_CHANNEL
718 IRRIGATION_CHANNEL =0.0
719 HFX_PHY = 0.0 ! initialize
732 FDTW=DT/(XLV*RHOWATER)
741 NSOIL=num_soil_layers
749 IF(ITIMESTEP.EQ.1)THEN
751 !*** initialize soil conditions for IHOP 31 May case
752 ! IF((XLAND(I,J)-1.5) < 0.)THEN
753 ! if (I==108.and.j==85) then
761 !*** SET ZERO-VALUE FOR SOME OUTPUT DIAGNOSTIC ARRAYS
762 IF((XLAND(I,J)-1.5).GE.0.)THEN
763 ! check sea-ice point
765 IF( XICE(I,J).GE. XICE_THRESHOLD .and. IPRINT ) PRINT*, ' sea-ice at water point, I=',I,'J=',J
772 TSLB(I,NS,J)=273.16 !STEMP
776 IF ( XICE(I,J) .GE. XICE_THRESHOLD ) THEN
788 ENDIF ! end of initialization over ocean
790 !-----------------------------------------------------------------------
794 ! pressure in middle of lowest layer
795 SFCPRS=(P8W3D(I,KTS+1,j)+P8W3D(i,KTS,j))*0.5
796 ! convert from mixing ratio to specific humidity
797 Q2K=QV3D(i,1,j)/(1.0+QV3D(i,1,j))
798 ! initializing local variables
806 Q2SAT=QGH(I,J)/(1.0+QGH(I,J)) ! Q2SAT is sp humidity
807 ! add check on myj=.true.
808 ! IF((Q2K.GE.Q2SAT*TRESH).AND.Q2K.LT.QZ0(I,J))THEN
809 IF((myj).AND.(Q2K.GE.Q2SAT*TRESH).AND.Q2K.LT.QZ0(I,J))THEN
820 ! TH2=SFCTMP+(0.0097545*ZLVL)
821 ! calculate SFCTH2 via Exner function vs lapse-rate (above)
822 APES=(1.E5/PSFC)**CAPA
823 APELM=(1.E5/SFCPRS)**CAPA
829 ! SOLDN is total incoming solar
831 ! GSW is net downward solar
833 ! use mid-day albedo to determine net downward solar (no solar zenith angle correction)
834 SOLNET=SOLDN*(1.-ALBEDO(I,J))
836 IF(PRESENT(IRR_CHAN)) THEN
837 IF(IRR_CHAN(i,j).NE.0) THEN
838 IRRIGATION_CHANNEL=IRR_CHAN(i,j)/DT
840 IRRIGATION_CHANNEL=0.
845 SHDFAC=VEGFRA(I,J)/100.
848 SHMIN=SHDMIN(I,J)/100. !NEW
849 SHMAX=SHDMAX(I,J)/100. !NEW
850 ! convert snow water equivalent from mm to meter
851 SNEQV=SNOW(I,J)*0.001
852 ! snow depth in meters
856 ! if "SR" present, set frac of frozen precip ("FFROZP") = snow-ratio ("SR", range:0-1)
857 ! SR from e.g. Ferrier microphysics
858 ! otherwise define from 1st atmos level temperature
862 IF (SFCTMP <= 273.15) THEN
869 IF((XLAND(I,J)-1.5).GE.0.)THEN ! begining of land/sea if block
871 TSK_RURAL(I,J)=TSK(I,J)
872 HFX_RURAL(I,J)=HFX(I,J)
873 QFX_RURAL(I,J)=QFX(I,J)
874 LH_RURAL(I,J)=LH(I,J)
875 EMISS_RURAL(I,J)=EMISS(I,J)
876 GRDFLX_RURAL(I,J)=GRDFLX(I,J)
879 ! Land or sea-ice case
881 IF (XICE(I,J) >= XICE_THRESHOLD) THEN
884 ELSE IF ( VEGTYP == ISICE ) THEN
888 ! Neither sea ice or land ice.
891 DQSDT2=Q2SAT*A23M4/(SFCTMP-A4)**2
893 IF(SNOW(I,J).GT.0.0)THEN
894 ! snow on surface (use ice saturation properties)
896 E2SAT=611.2*EXP(6174.*(1./273.15 - 1./SFCTSNO))
897 Q2SATI=0.622*E2SAT/(SFCPRS-E2SAT)
898 Q2SATI=Q2SATI/(1.0+Q2SATI) ! spec. hum.
899 IF (T1 .GT. 273.14) THEN
900 ! warm ground temps, weight the saturation between ice and water according to SNOWC
901 Q2SAT=Q2SAT*(1.-SNOWC(I,J)) + Q2SATI*SNOWC(I,J)
902 DQSDT2=DQSDT2*(1.-SNOWC(I,J)) + Q2SATI*6174./(SFCTSNO**2)*SNOWC(I,J)
904 ! cold ground temps, use ice saturation only
906 DQSDT2=Q2SATI*6174./(SFCTSNO**2)
908 ! for snow cover fraction at 0 C, ground temp will not change, so DQSDT2 effectively zero
909 ! V3.8 add condition for SWDOWN to restrict condition to positive forcing (JD)
910 IF(T1 .GT. 273. .AND. SNOWC(I,J) .GT. 0. .AND. SWDOWN(I,J) .GT. 10.)DQSDT2=DQSDT2*(1.-SNOWC(I,J))
913 ! Land-ice or land points use the usual deep-soil temperature.
916 IF(ISURBAN.EQ.1) THEN
917 ! assumes these only need to be set for USGS land data
918 IF(VEGTYP.EQ.25) SHDFAC=0.0000
919 IF(VEGTYP.EQ.26) SHDFAC=0.0000
920 IF(VEGTYP.EQ.27) SHDFAC=0.0000
922 IF(SOILTYP.EQ.14.AND.XICE(I,J).EQ.0.)THEN
924 IF(IPRINT)PRINT*,' SOIL TYPE FOUND TO BE WATER AT A LAND-POINT'
925 IF(IPRINT)PRINT*,i,j,'RESET SOIL in surfce.F'
929 SNOALB1 = SNOALB(I,J)
930 ! converts canwat in mm to CMC in meters
931 CMC=CANWAT(I,J)/1000.
933 !-------------------------------------------
934 !*** convert snow depth from mm to meter
937 ! SNOALB=ALBMAX(I,J)*0.01
939 ! SNOALB=MAXALB(IVGTPK)*0.01
947 SNOTIME1 = SNOTIME(I,J)
949 !FEI: temporaray arrays above need to be changed later by using SI
952 SMC(NS)=SMOIS(I,NS,J)
953 STC(NS)=TSLB(I,NS,J) !STEMP
957 if ( (SNEQV.ne.0..AND.SNOWHK.eq.0.).or.(SNOWHK.le.SNEQV) )THEN
962 !Fei: urban. for urban surface, if calling UCM, redefine the natural surface in cities as
963 ! the "NATURAL" category in the VEGPARM.TBL
964 IF(SF_URBAN_PHYSICS == 1.OR. SF_URBAN_PHYSICS==2.OR.SF_URBAN_PHYSICS==3 ) THEN
965 IF( IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == LCZ_1 .or. IVGTYP(I,J) == LCZ_2 .or. &
966 IVGTYP(I,J) == LCZ_3 .or. IVGTYP(I,J) == LCZ_4 .or. IVGTYP(I,J) == LCZ_5 .or. &
967 IVGTYP(I,J) == LCZ_6 .or. IVGTYP(I,J) == LCZ_7 .or. IVGTYP(I,J) == LCZ_8 .or. &
968 IVGTYP(I,J) == LCZ_9 .or. IVGTYP(I,J) == LCZ_10 .or. IVGTYP(I,J) == LCZ_11 ) THEN
971 SHDFAC = SHDTBL(NATURAL)
974 EMISSI = 0.98 !for VEGTYP=5
975 LWDN = GLW(I,J) * EMISSI
976 SOLNET = SOLDN * (1.0 - ALBEDOK)
978 IF ( FRC_URB2D(I,J) < 0.99 ) THEN
979 if(sf_urban_physics.eq.1)then
980 T1= ( TSK(I,J) -FRC_URB2D(I,J) * TS_URB2D (I,J) )/ (1-FRC_URB2D(I,J))
981 elseif((sf_urban_physics.eq.2).OR.(sf_urban_physics.eq.3))then
982 T1=tsk_rural_bep(i,j)
989 IF( IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == LCZ_1 .or. IVGTYP(I,J) == LCZ_2 .or. &
990 IVGTYP(I,J) == LCZ_3 .or. IVGTYP(I,J) == LCZ_4 .or. IVGTYP(I,J) == LCZ_5 .or. &
991 IVGTYP(I,J) == LCZ_6 .or. IVGTYP(I,J) == LCZ_7 .or. IVGTYP(I,J) == LCZ_8 .or. &
992 IVGTYP(I,J) == LCZ_9 .or. IVGTYP(I,J) == LCZ_10 .or. IVGTYP(I,J) == LCZ_11 ) THEN
998 !===Yang, 2014/10/08, hydrological processes for urban vegetation in single layer UCM===
1002 IRIOPTION=IRI_SCHEME
1003 IF(SF_URBAN_PHYSICS == 1) THEN
1005 tloc=mod(int(OMG/3.14159*180./15.+12.+0.5 ),24)
1006 if (tloc.lt.0) tloc=tloc+24
1007 if (tloc==0) tloc=24
1008 CALL cal_mon_day(julian,julyr,jmonth,jday)
1009 IF( IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == LCZ_1 .or. IVGTYP(I,J) == LCZ_2 .or. &
1010 IVGTYP(I,J) == LCZ_3 .or. IVGTYP(I,J) == LCZ_4 .or. IVGTYP(I,J) == LCZ_5 .or. &
1011 IVGTYP(I,J) == LCZ_6 .or. IVGTYP(I,J) == LCZ_7 .or. IVGTYP(I,J) == LCZ_8 .or. &
1012 IVGTYP(I,J) == LCZ_9 .or. IVGTYP(I,J) == LCZ_10 .or. IVGTYP(I,J) == LCZ_11 ) THEN
1013 AOASIS = oasis ! urban oasis effect
1014 IF (IRIOPTION ==1) THEN
1015 IF (tloc==21 .or. tloc==22) THEN !irrigation on vegetaion in urban area, MAY-SEP, 9-10pm
1016 IF (jmonth==5 .or. jmonth==6 .or. jmonth==7 .or. jmonth==8 .or. jmonth==9) THEN
1017 ! IF (SMC(USOIL) .LT. SMCREF) SMC(USOIL)= SMCREF
1018 ! IF (SMC(DSOIL) .LT. SMCREF) SMC(DSOIL)= SMCREF
1019 IF (SMC(USOIL) .LT. SMCREF) SMC(USOIL)= REFSMC(ISLTYP(I,J))
1020 IF (SMC(DSOIL) .LT. SMCREF) SMC(DSOIL)= REFSMC(ISLTYP(I,J))
1027 IF(SF_URBAN_PHYSICS == 2 .or. SF_URBAN_PHYSICS == 3) THEN
1028 IF(AOASIS > 1.0) THEN
1029 FATAL_ERROR('Urban oasis option is for SF_URBAN_PHYSICS == 1 only')
1031 IF(IRIOPTION == 1) THEN
1032 FATAL_ERROR('Urban irrigation option is for SF_URBAN_PHYSICS == 1 only')
1039 print*, 'BEFORE SFLX, in Noahlsm_driver'
1040 print*, 'ICE', ICE, 'DT',DT, 'ZLVL',ZLVL, 'NSOIL', NSOIL, &
1041 'SLDPTH', SLDPTH, 'LOCAL',LOCAL, 'LUTYPE',&
1042 LUTYPE, 'SLTYPE',SLTYPE, 'LWDN',LWDN, 'SOLDN',SOLDN, &
1043 'SFCPRS',SFCPRS, 'PRCP',PRCP,'SFCTMP',SFCTMP,'Q2K',Q2K, &
1044 'TH2',TH2,'Q2SAT',Q2SAT,'DQSDT2',DQSDT2,'VEGTYP', VEGTYP,&
1045 'SOILTYP',SOILTYP, 'SLOPETYP',SLOPETYP, 'SHDFAC',SHDFAC,&
1046 'SHMIN',SHMIN, 'ALBBRD',ALBBRD,'SNOALB1',SNOALB1,'TBOT',&
1047 TBOT, 'Z0BRD',Z0BRD, 'Z0K',Z0K, 'CMC',CMC, 'T1',T1,'STC',&
1048 STC, 'SMC',SMC, 'SWC',SWC,'SNOWHK',SNOWHK,'SNEQV',SNEQV,&
1049 'ALBEDOK',ALBEDOK,'CHK',CHK,'ETA',ETA,'SHEAT',SHEAT, &
1050 'ETA_KINEMATIC',ETA_KINEMATIC, 'FDOWN',FDOWN,'EC',EC, &
1051 'EDIR',EDIR,'ET',ET,'ETT',ETT,'ESNOW',ESNOW,'DRIP',DRIP,&
1052 'DEW',DEW,'BETA',BETA,'ETP',ETP,'SSOIL',SSOIL,'FLX1',FLX1,&
1053 'FLX2',FLX2,'FLX3',FLX3,'SNOMLT',SNOMLT,'SNCOVR',SNCOVR,&
1054 'RUNOFF1',RUNOFF1,'RUNOFF2',RUNOFF2,'RUNOFF3',RUNOFF3, &
1055 'RC',RC, 'PC',PC,'RSMIN',RSMIN,'XLAI',XLAI,'RCS',RCS, &
1056 'RCT',RCT,'RCQ',RCQ,'RCSOIL',RCSOIL,'SOILW',SOILW, &
1057 'SOILM',SOILM,'Q1',Q1,'SMCWLT',SMCWLT,'SMCDRY',SMCDRY,&
1058 'SMCREF',SMCREF,'SMCMAX',SMCMAX,'NROOT',NROOT
1063 IF (SHDFAC > 0.0 .AND. LAI(I,J) <= 0.0) LAI(I,J) = 0.01
1067 IF ( ICE == 1 ) THEN
1078 ELSEIF (ICE == 0) THEN
1084 IF( fasdas == 1 ) THEN
1087 XSDA_HFX= SDA_HFX(I,J)*RHO(I,1,J)*CPM(I,J)*DZQ ! W/m^2
1088 ! TWG2015 Bugfix remove factor of 1000.0 for correct units
1089 XSDA_QFX= SDA_QFX(I,J)*RHO(I,1,J)*DZQ ! Kg/m2/s of water
1098 sfcheadrt_hydro = sfcheadrt(i,j)
1099 infxsrt_hydro = infxsrt(i,j)
1102 sfcheadrt_hydro = 0.
1105 CALL SFLX (I,J,FFROZP, ISURBAN, DT,ZLVL,NSOIL,SLDPTH, & !C
1107 LUTYPE, SLTYPE, & !CL
1108 LWDN,SOLDN,SOLNET,SFCPRS,PRCP,SFCTMP,Q2K,DUMMY, & !F
1109 DUMMY,DUMMY, DUMMY, & !F PRCPRAIN not used
1110 TH2,Q2SAT,DQSDT2, & !I
1111 VEGTYP,SOILTYP,SLOPETYP,SHDFAC,SHMIN,SHMAX, & !I
1112 ALBBRD, SNOALB1,TBOT, Z0BRD, Z0K, EMISSI, EMBRD, & !S
1113 CMC,T1,STC,SMC,SWC,SNOWHK,SNEQV,ALBEDOK,CHK,dummy,& !H
1114 ETA,SHEAT, ETA_KINEMATIC,FDOWN, & !O
1115 EC,EDIR,ET,ETT,ESNOW,DRIP,DEW, & !O
1116 BETA,ETP,SSOIL, & !O
1117 FLX1,FLX2,FLX3, & !O
1118 FLX4,FVB,FBUR,FGSN,UA_PHYS, & !UA
1120 RUNOFF1,RUNOFF2,RUNOFF3, & !O
1121 RC,PC,RSMIN,XLAI,RCS,RCT,RCQ,RCSOIL, & !O
1122 SOILW,SOILM,Q1,SMAV, & !D
1123 RDLAI2D,USEMONALB, &
1126 SMCWLT,SMCDRY,SMCREF,SMCMAX,NROOT, &
1128 sfcheadrt_hydro, & !I
1129 INFXSRT_hydro,ETPND1_hydro & !O
1130 ,OPT_THCND,AOASIS & !O
1131 ,XSDA_QFX, HFX_PHY, QFX_PHY, XQNORM, fasdas, HCPCT_FASDAS & ! fasdas
1132 ,IRRIGATION_CHANNEL)
1135 soldrain(i,j) = RUNOFF2*DT*1000.0
1136 sfcheadrt(i,j) = sfcheadrt_hydro
1137 infxsrt(i,j) = INFXSRT_hydro
1138 etpnd1 = etpnd1_hydro
1140 ELSEIF (ICE == -1) THEN
1143 ! Set values that the LSM is expected to update,
1144 ! but don't get updated for glacial points.
1146 SOILM = 0.0 !BSINGH(PNNL)- SOILM is undefined for this case, it is used for diagnostics so setting it to zero
1147 XLAI = 0.01 ! KWM Should this be Zero over land ice? Does this value matter?
1158 IF( fasdas == 1 ) THEN
1161 XSDA_HFX= SDA_HFX(I,J)*RHO(I,1,J)*CPM(I,J)*DZQ ! W/m^2
1162 XSDA_QFX= 0.0 ! Kg/m2/s of water
1170 CALL SFLX_GLACIAL(I,J,ISICE,FFROZP,DT,ZLVL,NSOIL,SLDPTH, & !C
1171 & LWDN,SOLNET,SFCPRS,PRCP,SFCTMP,Q2K, & !F
1172 & TH2,Q2SAT,DQSDT2, & !I
1173 & ALBBRD, SNOALB1,TBOT, Z0BRD, Z0K, EMISSI, EMBRD, & !S
1174 & T1,STC(1:NSOIL),SNOWHK,SNEQV,ALBEDOK,CHK, & !H
1175 & ETA,SHEAT,ETA_KINEMATIC,FDOWN, & !O
1178 & FLX1,FLX2,FLX3, & !O
1179 & SNOMLT,SNCOVR, & !O
1188 if (present(rc2) .and. present(xlai2)) then
1189 rc2(I,J) = RC ! for output
1196 print*, 'AFTER SFLX, in Noahlsm_driver'
1197 print*, 'ICE', ICE, 'DT',DT, 'ZLVL',ZLVL, 'NSOIL', NSOIL, &
1198 'SLDPTH', SLDPTH, 'LOCAL',LOCAL, 'LUTYPE',&
1199 LUTYPE, 'SLTYPE',SLTYPE, 'LWDN',LWDN, 'SOLDN',SOLDN, &
1200 'SFCPRS',SFCPRS, 'PRCP',PRCP,'SFCTMP',SFCTMP,'Q2K',Q2K, &
1201 'TH2',TH2,'Q2SAT',Q2SAT,'DQSDT2',DQSDT2,'VEGTYP', VEGTYP,&
1202 'SOILTYP',SOILTYP, 'SLOPETYP',SLOPETYP, 'SHDFAC',SHDFAC,&
1203 'SHDMIN',SHMIN, 'ALBBRD',ALBBRD,'SNOALB',SNOALB1,'TBOT',&
1204 TBOT, 'Z0BRD',Z0BRD, 'Z0K',Z0K, 'CMC',CMC, 'T1',T1,'STC',&
1205 STC, 'SMC',SMC, 'SWc',SWC,'SNOWHK',SNOWHK,'SNEQV',SNEQV,&
1206 'ALBEDOK',ALBEDOK,'CHK',CHK,'ETA',ETA,'SHEAT',SHEAT, &
1207 'ETA_KINEMATIC',ETA_KINEMATIC, 'FDOWN',FDOWN,'EC',EC, &
1208 'EDIR',EDIR,'ET',ET,'ETT',ETT,'ESNOW',ESNOW,'DRIP',DRIP,&
1209 'DEW',DEW,'BETA',BETA,'ETP',ETP,'SSOIL',SSOIL,'FLX1',FLX1,&
1210 'FLX2',FLX2,'FLX3',FLX3,'SNOMLT',SNOMLT,'SNCOVR',SNCOVR,&
1211 'RUNOFF1',RUNOFF1,'RUNOFF2',RUNOFF2,'RUNOFF3',RUNOFF3, &
1212 'RC',RC, 'PC',PC,'RSMIN',RSMIN,'XLAI',XLAI,'RCS',RCS, &
1213 'RCT',RCT,'RCQ',RCQ,'RCSOIL',RCSOIL,'SOILW',SOILW, &
1214 'SOILM',SOILM,'Q1',Q1,'SMCWLT',SMCWLT,'SMCDRY',SMCDRY,&
1215 'SMCREF',SMCREF,'SMCMAX',SMCMAX,'NROOT',NROOT
1219 !*** UPDATE STATE VARIABLES
1220 CANWAT(I,J)=CMC*1000.
1221 SNOW(I,J)=SNEQV*1000.
1222 ! SNOWH(I,J)=SNOWHK*1000.
1223 SNOWH(I,J)=SNOWHK ! SNOWHK in meters
1225 ALB_RURAL(I,J)=ALBEDOK
1229 EMISS_RURAL(I,J) = EMISSI
1230 ! Noah: activate time-varying roughness length (V3.3 Feb 2011)
1235 ! Update Skin Temperature
1236 IF( fasdas == 1 ) THEN
1237 XSDA_QFX= XSDA_QFX*ELWV*XQNORM
1239 !TWG2015 Bugfix to multiply Heat Capacity by Soil Depth for correct
1242 T1 = T1 + (XSDA_HFX-XSDA_QFX)*DT/(HCPCT_FASDAS*DZS(1))
1250 if (present(tsk_rural_bep)) then
1251 IF(SF_URBAN_PHYSICS == 2 .or. SF_URBAN_PHYSICS == 3) THEN
1252 TSK_RURAL_BEP(I,J)=T1
1256 HFX_RURAL(I,J)=SHEAT
1257 ! MEk Jul07 add potential evap accum
1258 POTEVP(I,J)=POTEVP(I,J)+ETP*FDTW
1259 QFX(I,J)=ETA_KINEMATIC
1260 QFX_RURAL(I,J)=ETA_KINEMATIC
1264 ! QFX(I,J) = QFX(I,J) + ETPND1
1265 ! ETA = ETA + ETPND1/2.501E6*dt
1266 !end added by Wei Yu
1273 GRDFLX_RURAL(I,J)=SSOIL
1276 SNOTIME(I,J) = SNOTIME1
1277 ! prevent diagnostic ground q (q1) from being greater than qsat(tsk)
1278 ! as happens over snow cover where the cqs2 value also becomes irrelevant
1279 ! by setting cqs2=chs in this situation the 2m q should become just qv(k=1)
1280 ! ww: comment out this change to avoid Q2 drop due to change of radiative flux
1281 ! IF (Q1 .GT. QSFC(I,J)) THEN
1282 ! CQS2(I,J) = CHS(I,J)
1285 ! Convert QSFC back to mixing ratio
1286 QSFC(I,J)= Q1/(1.0-Q1)
1288 ! QSFC_RURAL(I,J)= Q1/(1.0-Q1)
1289 ! Calculate momentum flux from rural surface for use with multi-layer UCM (Martilli et al. 2002)
1292 SMOIS(I,NS,J)=SMC(NS)
1293 TSLB(I,NS,J)=STC(NS) ! STEMP
1294 SH2O(I,NS,J)=SWC(NS)
1304 ! Residual of surface energy balance equation terms
1308 noahres(i,j) = ( solnet + lwdn ) - sheat + ssoil - eta &
1309 - ( emissi * STBOLT * (t1**4) ) - flx1 - flx2 - flx3 - flx4
1312 noahres(i,j) = ( solnet + lwdn ) - sheat + ssoil - eta &
1313 - ( emissi * STBOLT * (t1**4) ) - flx1 - flx2 - flx3
1317 IF (SF_URBAN_PHYSICS == 1 ) THEN ! Beginning of UCM CALL if block
1318 !--------------------------------------
1319 ! URBAN CANOPY MODEL START - urban
1320 !--------------------------------------
1321 ! Input variables lsm --> urban
1323 IF( IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == LCZ_1 .or. IVGTYP(I,J) == LCZ_2 .or. &
1324 IVGTYP(I,J) == LCZ_3 .or. IVGTYP(I,J) == LCZ_4 .or. IVGTYP(I,J) == LCZ_5 .or. &
1325 IVGTYP(I,J) == LCZ_6 .or. IVGTYP(I,J) == LCZ_7 .or. IVGTYP(I,J) == LCZ_8 .or. &
1326 IVGTYP(I,J) == LCZ_9 .or. IVGTYP(I,J) == LCZ_10 .or. IVGTYP(I,J) == LCZ_11 ) THEN
1332 UTYPE_URB = UTYPE_URB2D(I,J) !urban type (low, high or industrial)
1334 TA_URB = SFCTMP ! [K]
1335 QA_URB = Q2K ! [kg/kg]
1336 UA_URB = SQRT(U_PHY(I,1,J)**2.+V_PHY(I,1,J)**2.)
1337 U1_URB = U_PHY(I,1,J)
1338 V1_URB = V_PHY(I,1,J)
1339 IF(UA_URB < 1.) UA_URB=1. ! [m/s]
1340 SSG_URB = SOLDN ! [W/m/m]
1341 SSGD_URB = 0.8*SOLDN ! [W/m/m]
1342 SSGQ_URB = SSG_URB-SSGD_URB ! [W/m/m]
1343 LLG_URB = GLW(I,J) ! [W/m/m]
1344 RAIN_URB = RAINBL(I,J) / DT * 3600.0 ! [mm/hr]
1345 RHOO_URB = SFCPRS / (287.04 * SFCTMP * (1.0+ 0.61 * Q2K)) ![kg/m/m/m]
1347 DELT_URB = DT ! [sec]
1348 XLAT_URB = XLAT_URB2D(I,J) ! [deg]
1349 COSZ_URB = COSZ_URB2D(I,J) !
1350 OMG_URB = OMG_URB2D(I,J) !
1353 LSOLAR_URB = .FALSE.
1355 TR_URB = TR_URB2D(I,J)
1356 TB_URB = TB_URB2D(I,J)
1357 TG_URB = TG_URB2D(I,J)
1358 TC_URB = TC_URB2D(I,J)
1359 QC_URB = QC_URB2D(I,J)
1360 UC_URB = UC_URB2D(I,J)
1362 TGR_URB = TGR_URB2D(I,J)
1363 CMCR_URB = CMCR_URB2D(I,J)
1364 FLXHUMR_URB = FLXHUMR_URB2D(I,J)
1365 FLXHUMB_URB = FLXHUMB_URB2D(I,J)
1366 FLXHUMG_URB = FLXHUMG_URB2D(I,J)
1367 DRELR_URB = DRELR_URB2D(I,J)
1368 DRELB_URB = DRELB_URB2D(I,J)
1369 DRELG_URB = DRELG_URB2D(I,J)
1371 DO K = 1,num_roof_layers
1372 TRL_URB(K) = TRL_URB3D(I,K,J)
1373 SMR_URB(K) = SMR_URB3D(I,K,J)
1374 TGRL_URB(K)= TGRL_URB3D(I,K,J)
1376 DO K = 1,num_wall_layers
1377 TBL_URB(K) = TBL_URB3D(I,K,J)
1379 DO K = 1,num_road_layers
1380 TGL_URB(K) = TGL_URB3D(I,K,J)
1383 XXXR_URB = XXXR_URB2D(I,J)
1384 XXXB_URB = XXXB_URB2D(I,J)
1385 XXXG_URB = XXXG_URB2D(I,J)
1386 XXXC_URB = XXXC_URB2D(I,J)
1389 ! Limits to avoid dividing by small number
1390 if (CHS(I,J) < 1.0E-02) then
1393 if (CHS2(I,J) < 1.0E-02) then
1396 if (CQS2(I,J) < 1.0E-02) then
1401 CHS2_URB = CHS2(I,J)
1402 IF (PRESENT(CMR_SFCDIF)) THEN
1403 CMR_URB = CMR_SFCDIF(I,J)
1404 CHR_URB = CHR_SFCDIF(I,J)
1405 CMGR_URB = CMGR_SFCDIF(I,J)
1406 CHGR_URB = CHGR_SFCDIF(I,J)
1407 CMC_URB = CMC_SFCDIF(I,J)
1408 CHC_URB = CHC_SFCDIF(I,J)
1412 mh_urb = mh_urb2d(I,J)
1413 stdh_urb = stdh_urb2d(I,J)
1414 lp_urb = lp_urb2d(I,J)
1415 hgt_urb = hgt_urb2d(I,J)
1418 lf_urb(K)=lf_urb2d(I,K,J)
1420 frc_urb = frc_urb2d(I,J)
1421 lb_urb = lb_urb2d(I,J)
1423 if (I.eq.73.and.J.eq.125)THEN
1426 ! Distributed aerodynamics
1427 lf_urb_s = lf_urb2d_s(I, J)
1428 z0_urb = z0_urb2d(I, J)
1429 vegfrac = vegfra(I, J) / 100
1432 CALL cal_mon_day(julian,julyr,jmonth,jday)
1433 CALL urban(LSOLAR_URB, & ! I
1434 num_roof_layers,num_wall_layers,num_road_layers, & ! C
1436 UTYPE_URB,TA_URB,QA_URB,UA_URB,U1_URB,V1_URB,SSG_URB, & ! I
1437 SSGD_URB,SSGQ_URB,LLG_URB,RAIN_URB,RHOO_URB, & ! I
1438 ZA_URB,DECLIN_URB,COSZ_URB,OMG_URB, & ! I
1439 XLAT_URB,DELT_URB,ZNT_URB, & ! I
1440 CHS_URB, CHS2_URB, & ! I
1441 TR_URB, TB_URB, TG_URB, TC_URB, QC_URB,UC_URB, & ! H
1442 TRL_URB,TBL_URB,TGL_URB, & ! H
1443 XXXR_URB, XXXB_URB, XXXG_URB, XXXC_URB, & ! H
1444 TS_URB,QS_URB,SH_URB,LH_URB,LH_KINEMATIC_URB, & ! O
1445 SW_URB,ALB_URB,LW_URB,G_URB,RN_URB,PSIM_URB,PSIH_URB, & ! O
1447 CMR_URB, CHR_URB, CMC_URB, CHC_URB, &
1448 U10_URB, V10_URB, TH2_URB, Q2_URB, & ! O
1449 UST_URB,mh_urb, stdh_urb, lf_urb, lp_urb, & ! 0
1450 hgt_urb,frc_urb,lb_urb, check,CMCR_URB,TGR_URB, & ! H
1451 TGRL_URB,SMR_URB,CMGR_URB,CHGR_URB,jmonth, & ! H
1452 DRELR_URB,DRELB_URB, & ! H
1453 DRELG_URB,FLXHUMR_URB,FLXHUMB_URB,FLXHUMG_URB, &
1454 lf_urb_s, z0_urb, vegfrac)
1459 print*, 'AFTER CALL URBAN'
1460 print*,'num_roof_layers',num_roof_layers, 'num_wall_layers', &
1462 'DZR',DZR,'DZB',DZB,'DZG',DZG,'UTYPE_URB',UTYPE_URB,'TA_URB', &
1464 'QA_URB',QA_URB,'UA_URB',UA_URB,'U1_URB',U1_URB,'V1_URB', &
1466 'SSG_URB',SSG_URB,'SSGD_URB',SSGD_URB,'SSGQ_URB',SSGQ_URB, &
1467 'LLG_URB',LLG_URB,'RAIN_URB',RAIN_URB,'RHOO_URB',RHOO_URB, &
1468 'ZA_URB',ZA_URB, 'DECLIN_URB',DECLIN_URB,'COSZ_URB',COSZ_URB,&
1469 'OMG_URB',OMG_URB,'XLAT_URB',XLAT_URB,'DELT_URB',DELT_URB, &
1470 'ZNT_URB',ZNT_URB,'TR_URB',TR_URB, 'TB_URB',TB_URB,'TG_URB',&
1471 TG_URB,'TC_URB',TC_URB,'QC_URB',QC_URB,'TRL_URB',TRL_URB, &
1472 'TBL_URB',TBL_URB,'TGL_URB',TGL_URB,'XXXR_URB',XXXR_URB, &
1473 'XXXB_URB',XXXB_URB,'XXXG_URB',XXXG_URB,'XXXC_URB',XXXC_URB,&
1474 'TS_URB',TS_URB,'QS_URB',QS_URB,'SH_URB',SH_URB,'LH_URB', &
1475 LH_URB, 'LH_KINEMATIC_URB',LH_KINEMATIC_URB,'SW_URB',SW_URB,&
1476 'ALB_URB',ALB_URB,'LW_URB',LW_URB,'G_URB',G_URB,'RN_URB', &
1477 RN_URB, 'PSIM_URB',PSIM_URB,'PSIH_URB',PSIH_URB, &
1478 'U10_URB',U10_URB,'V10_URB',V10_URB,'TH2_URB',TH2_URB, &
1479 'Q2_URB',Q2_URB,'CHS_URB',CHS_URB,'CHS2_URB',CHS2_URB
1483 TS_URB2D(I,J) = TS_URB
1485 ALBEDO(I,J) = FRC_URB2D(I,J)*ALB_URB+(1-FRC_URB2D(I,J))*ALBEDOK ![-]
1486 HFX(I,J) = FRC_URB2D(I,J)*SH_URB+(1-FRC_URB2D(I,J))*SHEAT ![W/m/m]
1487 QFX(I,J) = FRC_URB2D(I,J)*LH_KINEMATIC_URB &
1488 + (1-FRC_URB2D(I,J))*ETA_KINEMATIC ![kg/m/m/s]
1489 LH(I,J) = FRC_URB2D(I,J)*LH_URB+(1-FRC_URB2D(I,J))*ETA ![W/m/m]
1490 GRDFLX(I,J) = FRC_URB2D(I,J)*G_URB+(1-FRC_URB2D(I,J))*SSOIL ![W/m/m]
1491 TSK(I,J) = FRC_URB2D(I,J)*TS_URB+(1-FRC_URB2D(I,J))*T1 ![K]
1492 Q1 = FRC_URB2D(I,J)*QS_URB+(1-FRC_URB2D(I,J))*Q1 ![-]
1493 ! Convert QSFC back to mixing ratio
1494 QSFC(I,J)= Q1/(1.0-Q1)
1495 UST(I,J)= FRC_URB2D(I,J)*UST_URB+(1-FRC_URB2D(I,J))*UST(I,J) ![m/s]
1500 print*, ' FRC_URB2D', FRC_URB2D, &
1501 'ALB_URB',ALB_URB, 'ALBEDOK',ALBEDOK, &
1502 'ALBEDO(I,J)', ALBEDO(I,J), &
1503 'SH_URB',SH_URB,'SHEAT',SHEAT, 'HFX(I,J)',HFX(I,J), &
1504 'LH_KINEMATIC_URB',LH_KINEMATIC_URB,'ETA_KINEMATIC', &
1505 ETA_KINEMATIC, 'QFX(I,J)',QFX(I,J), &
1506 'LH_URB',LH_URB, 'ETA',ETA, 'LH(I,J)',LH(I,J), &
1507 'G_URB',G_URB,'SSOIL',SSOIL,'GRDFLX(I,J)', GRDFLX(I,J),&
1508 'TS_URB',TS_URB,'T1',T1,'TSK(I,J)',TSK(I,J), &
1509 'QS_URB',QS_URB,'Q1',Q1,'QSFC(I,J)',QSFC(I,J)
1515 ! Renew Urban State Varialbes
1517 TR_URB2D(I,J) = TR_URB
1518 TB_URB2D(I,J) = TB_URB
1519 TG_URB2D(I,J) = TG_URB
1520 TC_URB2D(I,J) = TC_URB
1521 QC_URB2D(I,J) = QC_URB
1522 UC_URB2D(I,J) = UC_URB
1524 TGR_URB2D(I,J) =TGR_URB
1525 CMCR_URB2D(I,J)=CMCR_URB
1526 FLXHUMR_URB2D(I,J)=FLXHUMR_URB
1527 FLXHUMB_URB2D(I,J)=FLXHUMB_URB
1528 FLXHUMG_URB2D(I,J)=FLXHUMG_URB
1529 DRELR_URB2D(I,J) = DRELR_URB
1530 DRELB_URB2D(I,J) = DRELB_URB
1531 DRELG_URB2D(I,J) = DRELG_URB
1533 DO K = 1,num_roof_layers
1534 TRL_URB3D(I,K,J) = TRL_URB(K)
1535 SMR_URB3D(I,K,J) = SMR_URB(K)
1536 TGRL_URB3D(I,K,J)= TGRL_URB(K)
1538 DO K = 1,num_wall_layers
1539 TBL_URB3D(I,K,J) = TBL_URB(K)
1541 DO K = 1,num_road_layers
1542 TGL_URB3D(I,K,J) = TGL_URB(K)
1544 XXXR_URB2D(I,J) = XXXR_URB
1545 XXXB_URB2D(I,J) = XXXB_URB
1546 XXXG_URB2D(I,J) = XXXG_URB
1547 XXXC_URB2D(I,J) = XXXC_URB
1549 SH_URB2D(I,J) = SH_URB
1550 LH_URB2D(I,J) = LH_URB
1551 G_URB2D(I,J) = G_URB
1552 RN_URB2D(I,J) = RN_URB
1553 PSIM_URB2D(I,J) = PSIM_URB
1554 PSIH_URB2D(I,J) = PSIH_URB
1555 GZ1OZ0_URB2D(I,J)= GZ1OZ0_URB
1556 U10_URB2D(I,J) = U10_URB
1557 V10_URB2D(I,J) = V10_URB
1558 TH2_URB2D(I,J) = TH2_URB
1559 Q2_URB2D(I,J) = Q2_URB
1560 UST_URB2D(I,J) = UST_URB
1561 AKMS_URB2D(I,J) = KARMAN * UST_URB2D(I,J)/(GZ1OZ0_URB2D(I,J)-PSIM_URB2D(I,J))
1562 IF (PRESENT(CMR_SFCDIF)) THEN
1563 CMR_SFCDIF(I,J) = CMR_URB
1564 CHR_SFCDIF(I,J) = CHR_URB
1565 CMGR_SFCDIF(I,J) = CMGR_URB
1566 CHGR_SFCDIF(I,J) = CHGR_URB
1567 CMC_SFCDIF(I,J) = CMC_URB
1568 CHC_SFCDIF(I,J) = CHC_URB
1572 ENDIF ! end of UCM CALL if block
1573 !--------------------------------------
1574 ! Urban Part End - urban
1575 !--------------------------------------
1579 SMSTOT(I,J)=SOILM*1000.
1581 SMCREL(I,NS,J)=SMAV(NS)
1585 ! Convert the water unit into mm
1586 SFCRUNOFF(I,J)=SFCRUNOFF(I,J)+RUNOFF1*DT*1000.0
1587 UDRUNOFF(I,J)=UDRUNOFF(I,J)+RUNOFF2*DT*1000.0
1588 ! snow defined when fraction of frozen precip (FFROZP) > 0.5,
1589 IF(FFROZP.GT.0.5)THEN
1590 ACSNOW(I,J)=ACSNOW(I,J)+PRCP*DT
1592 IF(SNOW(I,J).GT.0.)THEN
1593 ACSNOM(I,J)=ACSNOM(I,J)+SNOMLT*1000.
1594 ! accumulated snow-melt energy
1595 SNOPCX(I,J)=SNOPCX(I,J)-SNOMLT/FDTLIW
1598 ENDIF ! endif of land-sea test
1600 ENDDO ILOOP ! of I loop
1601 ENDDO JLOOP ! of J loop
1603 IF (SF_URBAN_PHYSICS == 2) THEN
1612 b_q_bep(i,kts:kte,j)=0.
1615 CALL BEP(frc_urb2d,utype_urb2d,itimestep,dz8w,dt,u_phy,v_phy, &
1616 th_phy,rho,p_phy,swdown,glw, &
1617 gmt,julday,xlong,xlat,declin_urb,cosz_urb2d,omg_urb2d, &
1618 num_urban_ndm, urban_map_zrd, urban_map_zwd, urban_map_gd, &
1619 urban_map_zd, urban_map_zdf, urban_map_bd, urban_map_wd, &
1620 urban_map_gbd, urban_map_fbd, num_urban_hi, &
1621 trb_urb4d,tw1_urb4d,tw2_urb4d,tgb_urb4d, &
1622 sfw1_urb3d,sfw2_urb3d,sfr_urb3d,sfg_urb3d, &
1623 lp_urb2d,hi_urb2d,lb_urb2d,hgt_urb2d, &
1624 a_u_bep,a_v_bep,a_t_bep, &
1625 a_e_bep,b_u_bep,b_v_bep, &
1626 b_t_bep,b_e_bep,b_q_bep,dlg_bep, &
1627 dl_u_bep,sf_bep,vl_bep, &
1628 rl_up_urb,rs_abs_urb,emiss_urb,grdflx_urb, &
1629 ids,ide, jds,jde, kds,kde, &
1630 ims,ime, jms,jme, kms,kme, &
1631 its,ite, jts,jte, kts,kte )
1636 IF (SF_URBAN_PHYSICS == 3) THEN
1645 b_q_bep(i,kts:kte,j)=0.
1650 CALL BEP_BEM(frc_urb2d,utype_urb2d,itimestep,dz8w,dt,u_phy,v_phy, &
1651 th_phy,rho,p_phy,swdown,glw, &
1652 gmt,julday,xlong,xlat,declin_urb,cosz_urb2d,omg_urb2d, &
1653 num_urban_ndm, urban_map_zrd, urban_map_zwd, urban_map_gd, &
1654 urban_map_zd, urban_map_zdf, urban_map_bd, urban_map_wd, &
1655 urban_map_gbd, urban_map_fbd, urban_map_zgrd, num_urban_hi, &
1656 trb_urb4d,tw1_urb4d,tw2_urb4d,tgb_urb4d, &
1657 tlev_urb3d,qlev_urb3d,tw1lev_urb3d,tw2lev_urb3d, &
1658 tglev_urb3d,tflev_urb3d,sf_ac_urb3d,lf_ac_urb3d, &
1659 cm_ac_urb3d,sfvent_urb3d,lfvent_urb3d, &
1660 sfwin1_urb3d,sfwin2_urb3d, &
1661 sfw1_urb3d,sfw2_urb3d,sfr_urb3d,sfg_urb3d, &
1662 ep_pv_urb3d,t_pv_urb3d, & !RMS
1663 trv_urb4d,qr_urb4d,qgr_urb3d,tgr_urb3d, & !RMS
1664 drain_urb4d,draingr_urb3d,sfrv_urb3d, & !RMS
1665 lfrv_urb3d,dgr_urb3d,dg_urb3d,lfr_urb3d,lfg_urb3d, & !RMS
1666 rainbl,swddir,swddif, & !RMS
1667 lp_urb2d,hi_urb2d,lb_urb2d,hgt_urb2d, &
1668 a_u_bep,a_v_bep,a_t_bep, &
1669 a_e_bep,b_u_bep,b_v_bep, &
1670 b_t_bep,b_e_bep,b_q_bep,dlg_bep, &
1671 dl_u_bep,sf_bep,vl_bep, &
1672 rl_up_urb,rs_abs_urb,emiss_urb,grdflx_urb,qv3d, &
1673 ids,ide, jds,jde, kds,kde, &
1674 ims,ime, jms,jme, kms,kme, &
1675 its,ite, jts,jte, kts,kte )
1679 if((sf_urban_physics.eq.2).OR.(sf_urban_physics.eq.3))then !Bep begin
1687 a_u_bep(i,k,j)=a_u_bep(i,k,j)*frc_urb2d(i,j)
1688 a_v_bep(i,k,j)=a_v_bep(i,k,j)*frc_urb2d(i,j)
1689 a_t_bep(i,k,j)=a_t_bep(i,k,j)*frc_urb2d(i,j)
1692 b_u_bep(i,k,j)=b_u_bep(i,k,j)*frc_urb2d(i,j)
1693 b_v_bep(i,k,j)=b_v_bep(i,k,j)*frc_urb2d(i,j)
1694 b_t_bep(i,k,j)=b_t_bep(i,k,j)*frc_urb2d(i,j)
1695 b_q_bep(i,k,j)=b_q_bep(i,k,j)*frc_urb2d(i,j)
1696 b_e_bep(i,k,j)=b_e_bep(i,k,j)*frc_urb2d(i,j)
1697 HFX_URB(I,J)=HFX_URB(I,J)+B_T_BEP(I,K,J)*RHO(I,K,J)*CP* &
1698 DZ8W(I,K,J)*VL_BEP(I,K,J)
1699 QFX_URB(I,J)=QFX_URB(I,J)+B_Q_BEP(I,K,J)* &
1700 DZ8W(I,K,J)*VL_BEP(I,K,J)
1701 UMOM_URB(I,J)=UMOM_URB(I,J)+ (A_U_BEP(I,K,J)*U_PHY(I,K,J)+ &
1702 B_U_BEP(I,K,J))*DZ8W(I,K,J)*VL_BEP(I,K,J)
1703 VMOM_URB(I,J)=VMOM_URB(I,J)+ (A_V_BEP(I,K,J)*V_PHY(I,K,J)+ &
1704 B_V_BEP(I,K,J))*DZ8W(I,K,J)*VL_BEP(I,K,J)
1705 vl_bep(i,k,j)=(1.-frc_urb2d(i,j))+vl_bep(i,k,j)*frc_urb2d(i,j)
1706 sf_bep(i,k,j)=(1.-frc_urb2d(i,j))+sf_bep(i,k,j)*frc_urb2d(i,j)
1708 a_u_bep(i,1,j)=(1.-frc_urb2d(i,j))*(-ust(I,J)*ust(I,J))/dz8w(i,1,j)/ &
1709 ((u_phy(i,1,j)**2+v_phy(i,1,j)**2.)**.5)+a_u_bep(i,1,j)
1710 a_v_bep(i,1,j)=(1.-frc_urb2d(i,j))*(-ust(I,J)*ust(I,J))/dz8w(i,1,j)/ &
1711 ((u_phy(i,1,j)**2+v_phy(i,1,j)**2.)**.5)+a_v_bep(i,1,j)
1712 b_t_bep(i,1,j)=(1.-frc_urb2d(i,j))*hfx_rural(i,j)/dz8w(i,1,j)/rho(i,1,j)/CP+ &
1714 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)
1715 umom=(1.-frc_urb2d(i,j))*ust(i,j)*ust(i,j)*u_phy(i,1,j)/ &
1716 ((u_phy(i,1,j)**2+v_phy(i,1,j)**2.)**.5)+umom_urb(i,j)
1717 vmom=(1.-frc_urb2d(i,j))*ust(i,j)*ust(i,j)*v_phy(i,1,j)/ &
1718 ((u_phy(i,1,j)**2+v_phy(i,1,j)**2.)**.5)+vmom_urb(i,j)
1721 ! compute upward longwave radiation from the rural part and total
1722 ! rl_up_rural=-emiss_rural(i,j)*sigma_sb*(tsk_rural(i,j)**4.)-(1.-emiss_rural(i,j))*glw(i,j)
1723 ! rl_up_tot=(1.-frc_urb2d(i,j))*rl_up_rural+frc_urb2d(i,j)*rl_up_urb(i,j)
1724 ! emiss(i,j)=(1.-frc_urb2d(i,j))*emiss_rural(i,j)+frc_urb2d(i,j)*emiss_urb(i,j)
1725 ! using the emissivity and the total longwave upward radiation estimate the averaged skin temperature
1726 IF (FRC_URB2D(I,J).GT.0.) THEN
1727 rl_up_rural=-emiss_rural(i,j)*sigma_sb*(tsk_rural(i,j)**4.)-(1.-emiss_rural(i,j))*glw(i,j)
1728 rl_up_tot=(1.-frc_urb2d(i,j))*rl_up_rural+frc_urb2d(i,j)*rl_up_urb(i,j)
1729 emiss(i,j)=(1.-frc_urb2d(i,j))*emiss_rural(i,j)+frc_urb2d(i,j)*emiss_urb(i,j)
1730 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
1731 tsk(i,j)=(max(0., (-1.*rl_up_tot-(1.-emiss(i,j))*glw(i,j) )/emiss(i,j)/sigma_sb))**.25
1732 rs_abs_tot=(1.-frc_urb2d(i,j))*swdown(i,j)*(1.-albedo(i,j))+frc_urb2d(i,j)*rs_abs_urb(i,j)
1733 if(swdown(i,j).gt.0.)then
1734 albedo(i,j)=1.-rs_abs_tot/swdown(i,j)
1736 albedo(i,j)=alb_rural(i,j)
1738 ! rename *_urb to sh_urb2d,lh_urb2d,g_urb2d,rn_urb2d
1739 grdflx(i,j)= (1.-frc_urb2d(i,j))*grdflx_rural(i,j)+frc_urb2d(i,j)*grdflx_urb(i,j)
1740 qfx(i,j)=(1.-frc_urb2d(i,j))*qfx_rural(i,j)+qfx_urb(i,j)
1741 ! lh(i,j)=(1.-frc_urb2d(i,j))*qfx_rural(i,j)*xlv
1742 lh(i,j)=qfx(i,j)*xlv
1743 HFX(I,J) = HFX_URB(I,J)+(1-FRC_URB2D(I,J))*HFX_RURAL(I,J) ![W/m/m]
1744 SH_URB2D(I,J) = HFX_URB(I,J)/FRC_URB2D(I,J)
1745 LH_URB2D(I,J) = qfx_urb(i,j)*xlv
1746 G_URB2D(I,J) = grdflx_urb(i,j)
1747 RN_URB2D(I,J) = rs_abs_urb(i,j)+emiss_urb(i,j)*glw(i,j)-rl_up_urb(i,j)
1748 ust(i,j)=(umom**2.+vmom**2.)**.25
1749 ! if(tsk(i,j).gt.350)write(*,*)'tsk too big!',i,j,tsk(i,j)
1750 ! 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
1751 ! print*,'ivgtyp,i,j,sigma_sb',ivgtyp(i,j),i,j,sigma_sb
1752 ! print*,'hfx,lh,qfx,grdflx,ts_urb2d',hfx(i,j),lh(i,j),qfx(i,j),grdflx(i,j),ts_urb2d(i,j)
1753 ! print*,'tsk,albedo,emiss',tsk(i,j),albedo(i,j),emiss(i,j)
1754 ! if(i.eq.56.and.j.eq.29)then
1755 ! print*,'ivgtyp, qfx, hfx',ivgtyp(i,j),hfx_rural(i,j),qfx_rural(i,j)
1756 ! print*,'emiss_rural,emiss_urb',emiss_rural(i,j),emiss_urb(i,j)
1757 ! print*,'rl_up_rural,rl_up_urb(i,j)',rl_up_rural,rl_up_urb(i,j)
1758 ! print*,'tsk_rural,ts_urb2d(i,j),tsk',tsk_rural(i,j),ts_urb2d(i,j),tsk(i,j)
1759 ! 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
1760 ! print*,'ivgtyp,hfx,hfx_urb,hfx_rural',hfx(i,j),hfx_urb(i,j),hfx_rural(i,j)
1761 ! print*,'lh,lh_rural',lh(i,j),lh_rural(i,j)
1762 ! print*,'qfx',qfx(i,j)
1763 ! print*,'ts_urb2d',ts_urb2d(i,j)
1764 ! print*,'ust',ust(i,j)
1765 ! print*,'swdown,glw',swdown(i,j),glw(i,j)
1778 !------------------------------------------------------
1780 !------------------------------------------------------
1782 !For MPAS, the below section of the module is moved to module_physics_lsm_noahinit.F in
1783 !./../core_physics to accomodate differences in the mpi calls between WRF and MPAS.I thought
1784 !that it would be cleaner to do this instead of adding a lot of #ifdef statements throughout
1785 !the initialization subroutine.
1787 #if defined(wrfmodel)
1788 SUBROUTINE LSMINIT(VEGFRA,SNOW,SNOWC,SNOWH,CANWAT,SMSTAV, &
1789 SMSTOT, SFCRUNOFF,UDRUNOFF,ACSNOW, &
1790 ACSNOM,IVGTYP,ISLTYP,TSLB,SMOIS,SH2O,ZS,DZS, &
1792 SNOALB, FNDSOILW, FNDSNOWH, RDMAXALB, &
1793 num_soil_layers, restart, &
1795 irr_rand_field,irr_ph,irr_freq, &
1796 ids,ide, jds,jde, kds,kde, &
1797 ims,ime, jms,jme, kms,kme, &
1798 its,ite, jts,jte, kts,kte &
1801 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
1802 ims,ime, jms,jme, kms,kme, &
1803 its,ite, jts,jte, kts,kte
1805 INTEGER, INTENT(IN) :: num_soil_layers
1807 LOGICAL , INTENT(IN) :: restart , allowed_to_read
1809 REAL, DIMENSION( num_soil_layers), INTENT(INOUT) :: ZS, DZS
1811 REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , &
1812 INTENT(INOUT) :: SMOIS, & !Total soil moisture
1813 SH2O, & !liquid soil moisture
1816 REAL, DIMENSION( ims:ime, jms:jme ) , &
1817 INTENT(INOUT) :: SNOW, &
1830 INTEGER, DIMENSION( ims:ime, jms:jme ) , &
1831 INTENT(IN) :: IVGTYP, &
1833 CHARACTER(LEN=*), INTENT(IN) :: MMINLU
1835 LOGICAL, INTENT(IN) :: FNDSOILW , &
1837 LOGICAL, INTENT(IN) :: RDMAXALB
1841 REAL :: BX, SMCMAX, PSISAT, FREE
1842 REAL, PARAMETER :: BLIM = 5.5, HLICE = 3.335E5, &
1843 GRAV = 9.81, T0 = 273.15
1845 CHARACTER(LEN=80) :: err_message
1846 INTEGER,DIMENSION(ims:ime, jms:jme ),INTENT(INOUT):: irr_rand_field
1847 INTEGER , DIMENSION(jds:jde) :: my_seeds
1848 INTEGER :: irr_ph,irr_freq
1849 REAL,DIMENSION(ims:ime, jms:jme ) :: rand_tmp
1850 character*256 :: MMINSL
1854 ! initialize three Noah LSM related tables
1855 IF ( allowed_to_read ) THEN
1856 CALL wrf_message( 'INITIALIZE THREE Noah LSM RELATED TABLES' )
1857 CALL SOIL_VEG_GEN_PARM( MMINLU, MMINSL )
1861 ! 20130219 - No longer need these - see module_data_gocart_dust
1862 !#if ( WRF_CHEM == 1 )
1864 !! need this parameter for dust parameterization in wrf/chem
1867 ! porosity(i)=maxsmc(i)
1868 ! drypoint(i)=drysmc(i)
1873 IF(.not.restart)THEN
1876 IF (irr_ph.NE.0)THEN
1879 my_seeds(j) =sqrt(ide*(real(j-1)**2))+sqrt(real(jde*i))
1880 ! PRINT*,'myseed', my_seeds(j),j,jts,jds
1882 CALL RANDOM_SEED ( PUT = my_seeds )
1883 CALL RANDOM_NUMBER ( rand_tmp(i,:) )
1884 CALL RANDOM_SEED ( GET = my_seeds )
1885 CALL RANDOM_NUMBER ( rand_tmp(i,:) )
1886 irr_rand_field(i,:)=int(modulo(rand_tmp(i,:)*100,real(irr_freq)))
1897 IF ( ISLTYP( i,j ) .LT. 1 ) THEN
1899 WRITE(err_message,*)"module_sf_noahdrv.F: lsminit: out of range ISLTYP ",i,j,ISLTYP( i,j )
1900 CALL wrf_message(err_message)
1902 IF(.not.RDMAXALB) THEN
1903 SNOALB(i,j)=MAXALB(IVGTYP(i,j))*0.01
1907 IF ( errflag .EQ. 1 ) THEN
1909 CALL wrf_message( "WARNING: message only; was fatal. module_sf_noahdrv.F: lsminit: out of range value "// &
1910 "of ISLTYP. Is this field in the input?" )
1912 CALL wrf_error_fatal( "module_sf_noahdrv.F: lsminit: out of range value "// &
1913 "of ISLTYP. Is this field in the input?" )
1917 ! initialize soil liquid water content SH2O
1919 ! IF(.NOT.FNDSOILW) THEN
1921 ! If no SWC, do the following
1922 ! PRINT *,'SOIL WATER NOT FOUND - VALUE SET IN LSMINIT'
1925 BX = BB(ISLTYP(I,J))
1926 SMCMAX = MAXSMC(ISLTYP(I,J))
1927 PSISAT = SATPSI(ISLTYP(I,J))
1928 if ((bx > 0.0).and.(smcmax > 0.0).and.(psisat > 0.0)) then
1929 DO NS=1, num_soil_layers
1930 ! ----------------------------------------------------------------------
1931 !SH2O <= SMOIS for T < 273.149K (-0.001C)
1932 IF (TSLB(I,NS,J) < 273.149) THEN
1933 ! ----------------------------------------------------------------------
1934 ! first guess following explicit solution for Flerchinger Eqn from Koren
1935 ! et al, JGR, 1999, Eqn 17 (KCOUNT=0 in FUNCTION FRH2O).
1936 ! ISLTPK is soil type
1937 BX = BB(ISLTYP(I,J))
1938 SMCMAX = MAXSMC(ISLTYP(I,J))
1939 PSISAT = SATPSI(ISLTYP(I,J))
1940 IF ( BX > BLIM ) BX = BLIM
1941 FK=(( (HLICE/(GRAV*(-PSISAT))) * &
1942 ((TSLB(I,NS,J)-T0)/TSLB(I,NS,J)) )**(-1/BX) )*SMCMAX
1943 IF (FK < 0.02) FK = 0.02
1944 SH2O(I,NS,J) = MIN( FK, SMOIS(I,NS,J) )
1945 ! ----------------------------------------------------------------------
1946 ! now use iterative solution for liquid soil water content using
1947 ! FUNCTION FRH2O with the initial guess for SH2O from above explicit
1949 CALL FRH2O (FREE,TSLB(I,NS,J),SMOIS(I,NS,J),SH2O(I,NS,J), &
1952 ELSE ! of IF (TSLB(I,NS,J)
1953 ! ----------------------------------------------------------------------
1954 ! SH2O = SMOIS ( for T => 273.149K (-0.001C)
1955 SH2O(I,NS,J)=SMOIS(I,NS,J)
1956 ! ----------------------------------------------------------------------
1957 ENDIF ! of IF (TSLB(I,NS,J)
1958 END DO ! of DO NS=1, num_soil_layers
1959 else ! of if ((bx > 0.0)
1960 DO NS=1, num_soil_layers
1961 SH2O(I,NS,J)=SMOIS(I,NS,J)
1963 endif ! of if ((bx > 0.0)
1964 ENDDO ! DO I = its,itf
1965 ENDDO ! DO J = jts,jtf
1966 ! ENDIF ! of IF(.NOT.FNDSOILW)THEN
1968 ! initialize physical snow height SNOWH
1970 IF(.NOT.FNDSNOWH)THEN
1971 ! If no SNOWH do the following
1972 CALL wrf_message( 'SNOW HEIGHT NOT FOUND - VALUE DEFINED IN LSMINIT' )
1975 SNOWH(I,J)=SNOW(I,J)*0.005 ! SNOW in mm and SNOWH in m
1980 ! initialize canopy water to ZERO
1983 ! print*,'Note that canopy water content (CANWAT) is set to ZERO in LSMINIT'
1992 !------------------------------------------------------------------------------
1993 END SUBROUTINE lsminit
1994 !------------------------------------------------------------------------------
1998 !-----------------------------------------------------------------
1999 SUBROUTINE SOIL_VEG_GEN_PARM( MMINLU, MMINSL)
2000 !-----------------------------------------------------------------
2002 USE module_wrf_error
2005 CHARACTER(LEN=*), INTENT(IN) :: MMINLU, MMINSL
2006 integer :: LUMATCH, IINDEX, LC, NUM_SLOPE
2008 INTEGER , PARAMETER :: OPEN_OK = 0
2010 character*128 :: mess , message
2011 character*256 :: a_string
2012 logical, external :: wrf_dm_on_monitor
2013 integer , parameter :: loop_max = 10
2014 integer :: loop_count
2017 !-----SPECIFY VEGETATION RELATED CHARACTERISTICS :
2018 ! ALBBCK: SFC albedo (in percentage)
2019 ! Z0: Roughness length (m)
2020 ! SHDFAC: Green vegetation fraction (in percentage)
2021 ! Note: The ALBEDO, Z0, and SHDFAC values read from the following table
2022 ! ALBEDO, amd Z0 are specified in LAND-USE TABLE; and SHDFAC is
2023 ! the monthly green vegetation data
2024 ! CMXTBL: MAX CNPY Capacity (m)
2025 ! NROTBL: Rooting depth (layer)
2026 ! RSMIN: Mimimum stomatal resistance (s m-1)
2027 ! RSMAX: Max. stomatal resistance (s m-1)
2028 ! RGL: Parameters used in radiation stress function
2029 ! HS: Parameter used in vapor pressure deficit functio
2030 ! TOPT: Optimum transpiration air temperature. (K)
2031 ! CMCMAX: Maximum canopy water capacity
2032 ! CFACTR: Parameter used in the canopy inteception calculati
2033 ! SNUP: Threshold snow depth (in water equivalent m) that
2034 ! implies 100% snow cover
2035 ! LAI: Leaf area index (dimensionless)
2036 ! MAXALB: Upper bound on maximum albedo over deep snow
2038 !-----READ IN VEGETAION PROPERTIES FROM VEGPARM.TBL
2041 IF ( wrf_dm_on_monitor() ) THEN
2043 OPEN(19, FILE='VEGPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
2044 IF(ierr .NE. OPEN_OK ) THEN
2045 WRITE(message,FMT='(A)') &
2046 'module_sf_noahlsm.F: soil_veg_gen_parm: failure opening VEGPARM.TBL'
2047 CALL wrf_error_fatal ( message )
2054 READ (19,FMT='(A)',END=2002) a_string
2055 FIND_LUTYPE : DO WHILE (LUMATCH == 0)
2056 READ (19,*,END=2002)LUTYPE
2057 READ (19,*)LUCATS,IINDEX
2059 IF(LUTYPE.EQ.MMINLU)THEN
2060 WRITE( mess , * ) 'LANDUSE TYPE = ' // TRIM ( LUTYPE ) // ' FOUND', LUCATS,' CATEGORIES'
2061 CALL wrf_message( mess )
2064 loop_count = loop_count+1
2065 call wrf_message ( "Skipping over LUTYPE = " // TRIM ( LUTYPE ) )
2066 FIND_VEGETATION_PARAMETER_FLAG : DO
2067 READ (19,FMT='(A)', END=2002) a_string
2068 IF ( a_string(1:21) .EQ. 'Vegetation Parameters' ) THEN
2069 EXIT FIND_VEGETATION_PARAMETER_FLAG
2070 ELSE IF ( loop_count .GE. loop_max ) THEN
2071 CALL wrf_error_fatal ( 'Too many loops in VEGPARM.TBL')
2073 ENDDO FIND_VEGETATION_PARAMETER_FLAG
2076 ! prevent possible array overwrite, Bill Bovermann, IBM, May 6, 2008
2077 IF ( SIZE(SHDTBL) < LUCATS .OR. &
2078 SIZE(NROTBL) < LUCATS .OR. &
2079 SIZE(RSTBL) < LUCATS .OR. &
2080 SIZE(RGLTBL) < LUCATS .OR. &
2081 SIZE(HSTBL) < LUCATS .OR. &
2082 SIZE(SNUPTBL) < LUCATS .OR. &
2083 SIZE(MAXALB) < LUCATS .OR. &
2084 SIZE(LAIMINTBL) < LUCATS .OR. &
2085 SIZE(LAIMAXTBL) < LUCATS .OR. &
2086 SIZE(Z0MINTBL) < LUCATS .OR. &
2087 SIZE(Z0MAXTBL) < LUCATS .OR. &
2088 SIZE(ALBEDOMINTBL) < LUCATS .OR. &
2089 SIZE(ALBEDOMAXTBL) < LUCATS .OR. &
2090 SIZE(ZTOPVTBL) < LUCATS .OR. &
2091 SIZE(ZBOTVTBL) < LUCATS .OR. &
2092 SIZE(EMISSMINTBL ) < LUCATS .OR. &
2093 SIZE(EMISSMAXTBL ) < LUCATS ) THEN
2094 CALL wrf_error_fatal('Table sizes too small for value of LUCATS in module_sf_noahdrv.F')
2097 IF(LUTYPE.EQ.MMINLU)THEN
2099 READ (19,*)IINDEX,SHDTBL(LC), &
2100 NROTBL(LC),RSTBL(LC),RGLTBL(LC),HSTBL(LC), &
2101 SNUPTBL(LC),MAXALB(LC), LAIMINTBL(LC), &
2102 LAIMAXTBL(LC),EMISSMINTBL(LC), &
2103 EMISSMAXTBL(LC), ALBEDOMINTBL(LC), &
2104 ALBEDOMAXTBL(LC), Z0MINTBL(LC), Z0MAXTBL(LC),&
2105 ZTOPVTBL(LC), ZBOTVTBL(LC)
2109 READ (19,*)TOPT_DATA
2111 READ (19,*)CMCMAX_DATA
2113 READ (19,*)CFACTR_DATA
2115 READ (19,*)RSMAX_DATA
2122 READ (19,FMT='(A)') a_string
2123 IF ( a_string(1:21) .EQ. 'Vegetation Parameters' ) THEN
2124 CALL wrf_message ("Expected low and high density residential, and high density industrial information in VEGPARM.TBL")
2125 CALL wrf_error_fatal ("This could be caused by using an older VEGPARM.TBL file with a newer WRF source code.")
2152 IF (LUMATCH == 0) then
2153 CALL wrf_error_fatal ("Land Use Dataset '"//MMINLU//"' not found in VEGPARM.TBL.")
2157 CALL wrf_dm_bcast_string ( LUTYPE , 4 )
2158 CALL wrf_dm_bcast_integer ( LUCATS , 1 )
2159 CALL wrf_dm_bcast_integer ( IINDEX , 1 )
2160 CALL wrf_dm_bcast_integer ( LUMATCH , 1 )
2161 CALL wrf_dm_bcast_real ( SHDTBL , NLUS )
2162 CALL wrf_dm_bcast_real ( NROTBL , NLUS )
2163 CALL wrf_dm_bcast_real ( RSTBL , NLUS )
2164 CALL wrf_dm_bcast_real ( RGLTBL , NLUS )
2165 CALL wrf_dm_bcast_real ( HSTBL , NLUS )
2166 CALL wrf_dm_bcast_real ( SNUPTBL , NLUS )
2167 CALL wrf_dm_bcast_real ( LAIMINTBL , NLUS )
2168 CALL wrf_dm_bcast_real ( LAIMAXTBL , NLUS )
2169 CALL wrf_dm_bcast_real ( Z0MINTBL , NLUS )
2170 CALL wrf_dm_bcast_real ( Z0MAXTBL , NLUS )
2171 CALL wrf_dm_bcast_real ( EMISSMINTBL , NLUS )
2172 CALL wrf_dm_bcast_real ( EMISSMAXTBL , NLUS )
2173 CALL wrf_dm_bcast_real ( ALBEDOMINTBL , NLUS )
2174 CALL wrf_dm_bcast_real ( ALBEDOMAXTBL , NLUS )
2175 CALL wrf_dm_bcast_real ( ZTOPVTBL , NLUS )
2176 CALL wrf_dm_bcast_real ( ZBOTVTBL , NLUS )
2177 CALL wrf_dm_bcast_real ( MAXALB , NLUS )
2178 CALL wrf_dm_bcast_real ( TOPT_DATA , 1 )
2179 CALL wrf_dm_bcast_real ( CMCMAX_DATA , 1 )
2180 CALL wrf_dm_bcast_real ( CFACTR_DATA , 1 )
2181 CALL wrf_dm_bcast_real ( RSMAX_DATA , 1 )
2182 CALL wrf_dm_bcast_integer ( BARE , 1 )
2183 CALL wrf_dm_bcast_integer ( NATURAL , 1 )
2184 CALL wrf_dm_bcast_integer ( LCZ_1 , 1 )
2185 CALL wrf_dm_bcast_integer ( LCZ_2 , 1 )
2186 CALL wrf_dm_bcast_integer ( LCZ_3 , 1 )
2187 CALL wrf_dm_bcast_integer ( LCZ_4 , 1 )
2188 CALL wrf_dm_bcast_integer ( LCZ_5 , 1 )
2189 CALL wrf_dm_bcast_integer ( LCZ_6 , 1 )
2190 CALL wrf_dm_bcast_integer ( LCZ_7 , 1 )
2191 CALL wrf_dm_bcast_integer ( LCZ_8 , 1 )
2192 CALL wrf_dm_bcast_integer ( LCZ_9 , 1 )
2193 CALL wrf_dm_bcast_integer ( LCZ_10 , 1 )
2194 CALL wrf_dm_bcast_integer ( LCZ_11 , 1 )
2197 !-----READ IN SOIL PROPERTIES FROM SOILPARM.TBL
2199 IF ( wrf_dm_on_monitor() ) THEN
2200 OPEN(19, FILE='SOILPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
2201 IF(ierr .NE. OPEN_OK ) THEN
2202 WRITE(message,FMT='(A)') &
2203 'module_sf_noahlsm.F: soil_veg_gen_parm: failure opening SOILPARM.TBL'
2204 CALL wrf_error_fatal ( message )
2207 WRITE(mess,*) 'INPUT SOIL TEXTURE CLASSIFICATION = ', TRIM ( MMINSL )
2208 CALL wrf_message( mess )
2213 READ (19,2000,END=2003)SLTYPE
2215 READ (19,*)SLCATS,IINDEX
2216 IF(SLTYPE.EQ.MMINSL)THEN
2217 WRITE( mess , * ) 'SOIL TEXTURE CLASSIFICATION = ', TRIM ( SLTYPE ) , ' FOUND', &
2218 SLCATS,' CATEGORIES'
2219 CALL wrf_message ( mess )
2222 ! prevent possible array overwrite, Bill Bovermann, IBM, May 6, 2008
2223 IF ( SIZE(BB ) < SLCATS .OR. &
2224 SIZE(DRYSMC) < SLCATS .OR. &
2225 SIZE(F11 ) < SLCATS .OR. &
2226 SIZE(MAXSMC) < SLCATS .OR. &
2227 SIZE(REFSMC) < SLCATS .OR. &
2228 SIZE(SATPSI) < SLCATS .OR. &
2229 SIZE(SATDK ) < SLCATS .OR. &
2230 SIZE(SATDW ) < SLCATS .OR. &
2231 SIZE(WLTSMC) < SLCATS .OR. &
2232 SIZE(QTZ ) < SLCATS ) THEN
2233 CALL wrf_error_fatal('Table sizes too small for value of SLCATS in module_sf_noahdrv.F')
2235 IF(SLTYPE.EQ.MMINSL)THEN
2237 READ (19,*) IINDEX,BB(LC),DRYSMC(LC),F11(LC),MAXSMC(LC),&
2238 REFSMC(LC),SATPSI(LC),SATDK(LC), SATDW(LC), &
2248 CALL wrf_dm_bcast_integer ( LUMATCH , 1 )
2249 CALL wrf_dm_bcast_string ( SLTYPE , 4 )
2250 CALL wrf_dm_bcast_string ( MMINSL , 4 ) ! since this is reset above, see oct2 ^
2251 CALL wrf_dm_bcast_integer ( SLCATS , 1 )
2252 CALL wrf_dm_bcast_integer ( IINDEX , 1 )
2253 CALL wrf_dm_bcast_real ( BB , NSLTYPE )
2254 CALL wrf_dm_bcast_real ( DRYSMC , NSLTYPE )
2255 CALL wrf_dm_bcast_real ( F11 , NSLTYPE )
2256 CALL wrf_dm_bcast_real ( MAXSMC , NSLTYPE )
2257 CALL wrf_dm_bcast_real ( REFSMC , NSLTYPE )
2258 CALL wrf_dm_bcast_real ( SATPSI , NSLTYPE )
2259 CALL wrf_dm_bcast_real ( SATDK , NSLTYPE )
2260 CALL wrf_dm_bcast_real ( SATDW , NSLTYPE )
2261 CALL wrf_dm_bcast_real ( WLTSMC , NSLTYPE )
2262 CALL wrf_dm_bcast_real ( QTZ , NSLTYPE )
2264 IF(LUMATCH.EQ.0)THEN
2265 CALL wrf_message( 'SOIl TEXTURE IN INPUT FILE DOES NOT ' )
2266 CALL wrf_message( 'MATCH SOILPARM TABLE' )
2267 CALL wrf_error_fatal ( 'INCONSISTENT OR MISSING SOILPARM FILE' )
2271 !-----READ IN GENERAL PARAMETERS FROM GENPARM.TBL
2273 IF ( wrf_dm_on_monitor() ) THEN
2274 OPEN(19, FILE='GENPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
2275 IF(ierr .NE. OPEN_OK ) THEN
2276 WRITE(message,FMT='(A)') &
2277 'module_sf_noahlsm.F: soil_veg_gen_parm: failure opening GENPARM.TBL'
2278 CALL wrf_error_fatal ( message )
2283 READ (19,*) NUM_SLOPE
2286 ! prevent possible array overwrite, Bill Bovermann, IBM, May 6, 2008
2287 IF ( SIZE(slope_data) < NUM_SLOPE ) THEN
2288 CALL wrf_error_fatal('NUM_SLOPE too large for slope_data array in module_sf_noahdrv')
2292 READ (19,*)SLOPE_DATA(LC)
2296 READ (19,*)SBETA_DATA
2298 READ (19,*)FXEXP_DATA
2300 READ (19,*)CSOIL_DATA
2302 READ (19,*)SALP_DATA
2304 READ (19,*)REFDK_DATA
2306 READ (19,*)REFKDT_DATA
2308 READ (19,*)FRZK_DATA
2310 READ (19,*)ZBOT_DATA
2312 READ (19,*)CZIL_DATA
2314 READ (19,*)SMLOW_DATA
2316 READ (19,*)SMHIGH_DATA
2318 READ (19,*)LVCOEF_DATA
2322 CALL wrf_dm_bcast_integer ( NUM_SLOPE , 1 )
2323 CALL wrf_dm_bcast_integer ( SLPCATS , 1 )
2324 CALL wrf_dm_bcast_real ( SLOPE_DATA , NSLOPE )
2325 CALL wrf_dm_bcast_real ( SBETA_DATA , 1 )
2326 CALL wrf_dm_bcast_real ( FXEXP_DATA , 1 )
2327 CALL wrf_dm_bcast_real ( CSOIL_DATA , 1 )
2328 CALL wrf_dm_bcast_real ( SALP_DATA , 1 )
2329 CALL wrf_dm_bcast_real ( REFDK_DATA , 1 )
2330 CALL wrf_dm_bcast_real ( REFKDT_DATA , 1 )
2331 CALL wrf_dm_bcast_real ( FRZK_DATA , 1 )
2332 CALL wrf_dm_bcast_real ( ZBOT_DATA , 1 )
2333 CALL wrf_dm_bcast_real ( CZIL_DATA , 1 )
2334 CALL wrf_dm_bcast_real ( SMLOW_DATA , 1 )
2335 CALL wrf_dm_bcast_real ( SMHIGH_DATA , 1 )
2336 CALL wrf_dm_bcast_real ( LVCOEF_DATA , 1 )
2339 !-----------------------------------------------------------------
2340 END SUBROUTINE SOIL_VEG_GEN_PARM
2341 !-----------------------------------------------------------------
2343 !===========================================================================
2345 ! subroutine lsm_mosaic: a tiling approach for Noah LSM
2347 !===========================================================================
2349 SUBROUTINE lsm_mosaic(DZ8W,QV3D,P8W3D,T3D,TSK, &
2350 HFX,QFX,LH,GRDFLX, QGH,GSW,SWDOWN,GLW,SMSTAV,SMSTOT, &
2351 SFCRUNOFF, UDRUNOFF,IVGTYP,ISLTYP,ISURBAN,ISICE,VEGFRA, &
2352 ALBEDO,ALBBCK,ZNT,Z0,TMN,XLAND,XICE,EMISS,EMBCK, &
2353 SNOWC,QSFC,RAINBL,MMINLU, &
2354 num_soil_layers,DT,DZS,ITIMESTEP, &
2355 SMOIS,TSLB,SNOW,CANWAT, &
2356 CHS,CHS2,CQS2,CPM,ROVCP,SR,chklowq,lai,qz0, & !H
2360 SNOALB,SHDMIN,SHDMAX,SHDAVG, & !I
2367 RDLAI2D,USEMONALB, &
2369 NOAHRES,OPT_THCND, &
2370 NLCAT,landusef,landusef2, & ! danli mosaic
2371 sf_surface_mosaic,mosaic_cat,mosaic_cat_index, & ! danli mosaic
2372 TSK_mosaic,QSFC_mosaic, & ! danli mosaic
2373 TSLB_mosaic,SMOIS_mosaic,SH2O_mosaic, & ! danli mosaic
2374 CANWAT_mosaic,SNOW_mosaic, & ! danli mosaic
2375 SNOWH_mosaic,SNOWC_mosaic, & ! danli mosaic
2376 ALBEDO_mosaic,ALBBCK_mosaic, & ! danli mosaic
2377 EMISS_mosaic, EMBCK_mosaic, & ! danli mosaic
2378 ZNT_mosaic, Z0_mosaic, & ! danli mosaic
2379 HFX_mosaic,QFX_mosaic, & ! danli mosaic
2380 LH_mosaic, GRDFLX_mosaic, SNOTIME_mosaic, & ! danli mosaic
2381 RC_mosaic, LAI_mosaic, &
2383 ua_phys,flx4_2d,fvb_2d,fbur_2d,fgsn_2d, &
2384 ids,ide, jds,jde, kds,kde, &
2385 ims,ime, jms,jme, kms,kme, &
2386 its,ite, jts,jte, kts,kte, &
2388 CMR_SFCDIF,CHR_SFCDIF,CMC_SFCDIF,CHC_SFCDIF, &
2389 CMGR_SFCDIF,CHGR_SFCDIF, &
2391 TR_URB2D,TB_URB2D,TG_URB2D,TC_URB2D,QC_URB2D, & !H urban
2392 UC_URB2D, & !H urban
2393 XXXR_URB2D,XXXB_URB2D,XXXG_URB2D,XXXC_URB2D, & !H urban
2394 TRL_URB3D,TBL_URB3D,TGL_URB3D, & !H urban
2395 SH_URB2D,LH_URB2D,G_URB2D,RN_URB2D,TS_URB2D, & !H urban
2396 TR_URB2D_mosaic,TB_URB2D_mosaic, & !H urban danli mosaic
2397 TG_URB2D_mosaic,TC_URB2D_mosaic, & !H urban danli mosaic
2398 QC_URB2D_mosaic,UC_URB2D_mosaic, & !H urban danli mosaic
2399 TRL_URB3D_mosaic,TBL_URB3D_mosaic, & !H urban danli mosaic
2400 TGL_URB3D_mosaic, & !H urban danli mosaic
2401 SH_URB2D_mosaic,LH_URB2D_mosaic, & !H urban danli mosaic
2402 G_URB2D_mosaic,RN_URB2D_mosaic, & !H urban danli mosaic
2403 TS_URB2D_mosaic, & !H urban danli mosaic
2404 TS_RUL2D_mosaic, & !H urban danli mosaic
2405 PSIM_URB2D,PSIH_URB2D,U10_URB2D,V10_URB2D, & !O urban
2406 GZ1OZ0_URB2D, AKMS_URB2D, & !O urban
2407 TH2_URB2D,Q2_URB2D, UST_URB2D, & !O urban
2408 DECLIN_URB,COSZ_URB2D,OMG_URB2D, & !I urban
2409 XLAT_URB2D, & !I urban
2410 num_roof_layers, num_wall_layers, & !I urban
2411 num_road_layers, DZR, DZB, DZG, & !I urban
2412 CMCR_URB2D,TGR_URB2D,TGRL_URB3D,SMR_URB3D, & !H urban
2413 julian,julyr, & !H urban
2414 DRELR_URB2D,DRELB_URB2D,DRELG_URB2D, & !H urban
2415 FLXHUMR_URB2D,FLXHUMB_URB2D,FLXHUMG_URB2D, & !H urban
2416 FRC_URB2D,UTYPE_URB2D, & !O
2417 num_urban_ndm, & !I multi-layer urban
2418 urban_map_zrd, & !I multi-layer urban
2419 urban_map_zwd, & !I multi-layer urban
2420 urban_map_gd, & !I multi-layer urban
2421 urban_map_zd, & !I multi-layer urban
2422 urban_map_zdf, & !I multi-layer urban
2423 urban_map_bd, & !I multi-layer urban
2424 urban_map_wd, & !I multi-layer urban
2425 urban_map_gbd, & !I multi-layer urban
2426 urban_map_fbd, & !I multi-layer urban
2427 urban_map_zgrd, & !I multi-layer urban
2428 num_urban_hi, & !I multi-layer urban
2429 use_wudapt_lcz, & !I wudapt
2430 slucm_distributed_drag, & !I slucm
2431 tsk_rural_bep, & !H multi-layer urban
2432 trb_urb4d,tw1_urb4d,tw2_urb4d,tgb_urb4d, & !H multi-layer urban
2433 tlev_urb3d,qlev_urb3d, & !H multi-layer urban
2434 tw1lev_urb3d,tw2lev_urb3d, & !H multi-layer urban
2435 tglev_urb3d,tflev_urb3d, & !H multi-layer urban
2436 sf_ac_urb3d,lf_ac_urb3d,cm_ac_urb3d, & !H multi-layer urban
2437 sfvent_urb3d,lfvent_urb3d, & !H multi-layer urban
2438 sfwin1_urb3d,sfwin2_urb3d, & !H multi-layer urban
2439 sfw1_urb3d,sfw2_urb3d,sfr_urb3d,sfg_urb3d, & !H multi-layer urban
2440 ep_pv_urb3d,t_pv_urb3d, & !RMS
2441 trv_urb4d,qr_urb4d,qgr_urb3d,tgr_urb3d, & !RMS
2442 drain_urb4d,draingr_urb3d,sfrv_urb3d, & !RMS
2443 lfrv_urb3d,dgr_urb3d,dg_urb3d,lfr_urb3d,lfg_urb3d,& !RMS
2444 lp_urb2d,hi_urb2d,lb_urb2d,hgt_urb2d, & !H multi-layer urban
2445 mh_urb2d,stdh_urb2d,lf_urb2d, & !SLUCM
2446 lf_urb2d_s, z0_urb2d, & !SLUCM
2447 th_phy,rho,p_phy,ust, & !I multi-layer urban
2448 gmt,julday,xlong,xlat, & !I multi-layer urban
2449 a_u_bep,a_v_bep,a_t_bep,a_q_bep, & !O multi-layer urban
2450 a_e_bep,b_u_bep,b_v_bep, & !O multi-layer urban
2451 b_t_bep,b_q_bep,b_e_bep,dlg_bep, & !O multi-layer urban
2452 dl_u_bep,sf_bep,vl_bep & !O multi-layer urban
2454 ,sfcheadrt,INFXSRT, soldrain & !hydro
2456 ,SDA_HFX, SDA_QFX, HFX_BOTH, QFX_BOTH, QNORM, fasdas & !fasdas
2461 !----------------------------------------------------------------
2463 !----------------------------------------------------------------
2464 !----------------------------------------------------------------
2465 ! --- atmospheric (WRF generic) variables
2466 !-- DT time step (seconds)
2467 !-- DZ8W thickness of layers (m)
2468 !-- T3D temperature (K)
2469 !-- QV3D 3D water vapor mixing ratio (Kg/Kg)
2470 !-- P3D 3D pressure (Pa)
2471 !-- FLHC exchange coefficient for heat (m/s)
2472 !-- FLQC exchange coefficient for moisture (m/s)
2473 !-- PSFC surface pressure (Pa)
2474 !-- XLAND land mask (1 for land, 2 for water)
2475 !-- QGH saturated mixing ratio at 2 meter
2476 !-- GSW downward short wave flux at ground surface (W/m^2)
2477 !-- GLW downward long wave flux at ground surface (W/m^2)
2478 !-- History variables
2479 !-- CANWAT canopy moisture content (mm)
2480 !-- TSK surface temperature (K)
2481 !-- TSLB soil temp (k)
2482 !-- SMOIS total soil moisture content (volumetric fraction)
2483 !-- SH2O unfrozen soil moisture content (volumetric fraction)
2484 ! note: frozen soil moisture (i.e., soil ice) = SMOIS - SH2O
2485 !-- SNOWH actual snow depth (m)
2486 !-- SNOW liquid water-equivalent snow depth (m)
2487 !-- ALBEDO time-varying surface albedo including snow effect (unitless fraction)
2488 !-- ALBBCK background surface albedo (unitless fraction)
2489 !-- CHS surface exchange coefficient for heat and moisture (m s-1);
2490 !-- CHS2 2m surface exchange coefficient for heat (m s-1);
2491 !-- CQS2 2m surface exchange coefficient for moisture (m s-1);
2492 ! --- soil variables
2493 !-- num_soil_layers the number of soil layers
2494 !-- ZS depths of centers of soil layers (m)
2495 !-- DZS thicknesses of soil layers (m)
2496 !-- SLDPTH thickness of each soil layer (m, same as DZS)
2497 !-- TMN soil temperature at lower boundary (K)
2498 !-- SMCWLT wilting point (volumetric)
2499 !-- SMCDRY dry soil moisture threshold where direct evap from
2500 ! top soil layer ends (volumetric)
2501 !-- SMCREF soil moisture threshold below which transpiration begins to
2502 ! stress (volumetric)
2503 !-- SMCMAX porosity, i.e. saturated value of soil moisture (volumetric)
2504 !-- NROOT number of root layers, a function of veg type, determined
2505 ! in subroutine redprm.
2506 !-- SMSTAV Soil moisture availability for evapotranspiration (
2507 ! fraction between SMCWLT and SMCMXA)
2508 !-- SMSTOT Total soil moisture content frozen+unfrozen) in the soil column (mm)
2509 ! --- snow variables
2510 !-- SNOWC fraction snow coverage (0-1.0)
2511 ! --- vegetation variables
2512 !-- SNOALB upper bound on maximum albedo over deep snow
2513 !-- SHDMIN minimum areal fractional coverage of annual green vegetation
2514 !-- SHDMAX maximum areal fractional coverage of annual green vegetation
2515 !-- XLAI leaf area index (dimensionless)
2516 !-- Z0BRD Background fixed roughness length (M)
2517 !-- Z0 Background vroughness length (M) as function
2518 !-- ZNT Time varying roughness length (M) as function
2519 !-- ALBD(IVGTPK,ISN) background albedo reading from a table
2521 !-- HFX upward heat flux at the surface (W/m^2)
2522 !-- QFX upward moisture flux at the surface (kg/m^2/s)
2523 !-- LH upward moisture flux at the surface (W m-2)
2524 !-- GRDFLX(I,J) ground heat flux (W m-2)
2525 !-- FDOWN radiation forcing at the surface (W m-2) = SOLDN*(1-alb)+LWDN
2526 !----------------------------------------------------------------------------
2527 !-- EC canopy water evaporation ((W m-2)
2528 !-- EDIR direct soil evaporation (W m-2)
2529 !-- ET plant transpiration from a particular root layer (W m-2)
2530 !-- ETT total plant transpiration (W m-2)
2531 !-- ESNOW sublimation from (or deposition to if <0) snowpack (W m-2)
2532 !-- DRIP through-fall of precip and/or dew in excess of canopy
2533 ! water-holding capacity (m)
2534 !-- DEW dewfall (or frostfall for t<273.15) (M)
2535 !-- SMAV Soil Moisture Availability for each layer, as a fraction
2536 ! between SMCWLT and SMCMAX (dimensionless fraction)
2537 ! ----------------------------------------------------------------------
2538 !-- BETA ratio of actual/potential evap (dimensionless)
2539 !-- ETP potential evaporation (W m-2)
2540 ! ----------------------------------------------------------------------
2541 !-- FLX1 precip-snow sfc (W m-2)
2542 !-- FLX2 freezing rain latent heat flux (W m-2)
2543 !-- FLX3 phase-change heat flux from snowmelt (W m-2)
2544 ! ----------------------------------------------------------------------
2545 !-- ACSNOM snow melt (mm) (water equivalent)
2546 !-- ACSNOW accumulated snow fall (mm) (water equivalent)
2547 !-- SNOPCX snow phase change heat flux (W/m^2)
2548 !-- POTEVP accumulated potential evaporation (m)
2549 !-- RIB Documentation needed!!!
2550 ! ----------------------------------------------------------------------
2551 !-- RUNOFF1 surface runoff (m s-1), not infiltrating the surface
2552 !-- RUNOFF2 subsurface runoff (m s-1), drainage out bottom of last
2553 ! soil layer (baseflow)
2554 ! important note: here RUNOFF2 is actually the sum of RUNOFF2 and RUNOFF3
2555 !-- RUNOFF3 numerical trunctation in excess of porosity (smcmax)
2556 ! for a given soil layer at the end of a time step (m s-1).
2557 !SFCRUNOFF Surface Runoff (mm)
2558 !UDRUNOFF Total Underground Runoff (mm), which is the sum of RUNOFF2 and RUNOFF3
2559 ! ----------------------------------------------------------------------
2560 !-- RC canopy resistance (s m-1)
2561 !-- PC plant coefficient (unitless fraction, 0-1) where PC*ETP = actual transp
2562 !-- RSMIN minimum canopy resistance (s m-1)
2563 !-- RCS incoming solar rc factor (dimensionless)
2564 !-- RCT air temperature rc factor (dimensionless)
2565 !-- RCQ atmos vapor pressure deficit rc factor (dimensionless)
2566 !-- RCSOIL soil moisture rc factor (dimensionless)
2568 !-- EMISS surface emissivity (between 0 and 1)
2569 !-- EMBCK Background surface emissivity (between 0 and 1)
2572 ! (R_d/R_v) (dimensionless)
2573 !-- ids start index for i in domain
2574 !-- ide end index for i in domain
2575 !-- jds start index for j in domain
2576 !-- jde end index for j in domain
2577 !-- kds start index for k in domain
2578 !-- kde end index for k in domain
2579 !-- ims start index for i in memory
2580 !-- ime end index for i in memory
2581 !-- jms start index for j in memory
2582 !-- jme end index for j in memory
2583 !-- kms start index for k in memory
2584 !-- kme end index for k in memory
2585 !-- its start index for i in tile
2586 !-- ite end index for i in tile
2587 !-- jts start index for j in tile
2588 !-- jte end index for j in tile
2589 !-- kts start index for k in tile
2590 !-- kte end index for k in tile
2592 !-- SR fraction of frozen precip (0.0 to 1.0)
2593 !----------------------------------------------------------------
2597 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
2598 ims,ime, jms,jme, kms,kme, &
2599 its,ite, jts,jte, kts,kte
2601 INTEGER, INTENT(IN ) :: sf_urban_physics !urban
2602 INTEGER, INTENT(IN ) :: isurban
2603 INTEGER, INTENT(IN ) :: isice
2604 INTEGER, INTENT(IN ) :: julian,julyr
2606 !added by Wei Yu for routing
2608 REAL, DIMENSION( ims:ime, jms:jme ) , &
2609 INTENT(INOUT) :: sfcheadrt,INFXSRT,soldrain
2614 ! new local vars for hydro
2615 REAL :: etpnd1_hydro,sfcheadrt_hydro,infxsrt_hydro
2617 REAL, DIMENSION( ims:ime, jms:jme ) , &
2618 INTENT(IN ) :: TMN, &
2627 SWDOWN, & !added 10 jan 2007
2631 REAL, DIMENSION( ims:ime, jms:jme ) , &
2632 INTENT(INOUT) :: ALBBCK, &
2634 EMBCK ! danli mosaic
2636 CHARACTER(LEN=*), INTENT(IN ) :: MMINLU
2638 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
2639 INTENT(IN ) :: QV3D, &
2643 REAL, DIMENSION( ims:ime, jms:jme ) , &
2644 INTENT(IN ) :: QGH, &
2647 INTEGER, DIMENSION( ims:ime, jms:jme ) , &
2648 INTENT(IN ) :: ISLTYP
2650 INTEGER, DIMENSION( ims:ime, jms:jme ) , &
2651 INTENT(INOUT ) :: IVGTYP ! for mosaic danli
2653 INTEGER, INTENT(IN) :: num_soil_layers,ITIMESTEP
2655 REAL, INTENT(IN ) :: DT,ROVCP
2657 REAL, DIMENSION(1:num_soil_layers), INTENT(IN)::DZS
2661 REAL, DIMENSION( ims:ime , 1:num_soil_layers, jms:jme ), &
2662 INTENT(INOUT) :: SMOIS, & ! total soil moisture
2663 SH2O, & ! new soil liquid
2666 REAL, DIMENSION( ims:ime , 1:num_soil_layers, jms:jme ), &
2667 INTENT(OUT) :: SMCREL
2669 REAL, DIMENSION( ims:ime, jms:jme ) , &
2670 INTENT(INOUT) :: TSK, & !was TGB (temperature)
2696 REAL, DIMENSION( ims:ime, jms:jme ) , &
2697 INTENT(OUT) :: NOAHRES
2698 INTEGER, INTENT(IN) :: OPT_THCND
2701 LOGICAL, INTENT(IN) :: UA_PHYS
2702 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: FLX4_2D,FVB_2D,FBUR_2D,FGSN_2D
2703 REAL :: FLX4,FVB,FBUR,FGSN
2705 REAL, DIMENSION( ims:ime, jms:jme ) , &
2706 INTENT(OUT) :: CHKLOWQ
2707 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: LAI
2708 REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: QZ0
2709 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: RC2, XLAI2
2711 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CMR_SFCDIF
2712 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CHR_SFCDIF
2713 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CMGR_SFCDIF
2714 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CHGR_SFCDIF
2715 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CMC_SFCDIF
2716 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CHC_SFCDIF
2717 ! Local variables (moved here from driver to make routine thread safe, 20031007 jm)
2719 REAL, DIMENSION(1:num_soil_layers) :: ET
2721 REAL, DIMENSION(1:num_soil_layers) :: SMAV
2723 REAL :: BETA, ETP, SSOIL,EC, EDIR, ESNOW, ETT, &
2724 FLX1,FLX2,FLX3, DRIP,DEW,FDOWN,RC,PC,RSMIN,XLAI, &
2725 ! RCS,RCT,RCQ,RCSOIL
2726 RCS,RCT,RCQ,RCSOIL,FFROZP
2728 LOGICAL, INTENT(IN ) :: myj,frpcpn
2730 ! DECLARATIONS - LOGICAL
2731 ! ----------------------------------------------------------------------
2732 LOGICAL, PARAMETER :: LOCAL=.false.
2733 LOGICAL :: FRZGRA, SNOWNG
2737 ! ----------------------------------------------------------------------
2738 ! DECLARATIONS - INTEGER
2739 ! ----------------------------------------------------------------------
2740 INTEGER :: I,J, ICE,NSOIL,SLOPETYP,SOILTYP,VEGTYP
2744 ! ----------------------------------------------------------------------
2745 ! DECLARATIONS - REAL
2746 ! ----------------------------------------------------------------------
2748 REAL :: SHMIN,SHMAX,DQSDT2,LWDN,PRCP,PRCPRAIN, &
2749 Q2SAT,Q2SATI,SFCPRS,SFCSPD,SFCTMP,SHDFAC,SNOALB1, &
2750 SOLDN,TBOT,ZLVL, Q2K,ALBBRD, ALBEDOK, ETA, ETA_KINEMATIC, &
2752 Z0K,RUNOFF1,RUNOFF2,RUNOFF3,SHEAT,SOLNET,E2SAT,SFCTSNO, &
2753 ! mek, WRF testing, expanded diagnostics
2754 SOLUP,LWUP,RNET,RES,Q1SFC,TAIRV,SATFLG
2757 ! MEK JUL2007 for pot. evap.
2763 REAL :: SNCOVR,SNEQV,SNOWHK,CMC, CHK,TH2
2765 REAL :: SMCDRY,SMCMAX,SMCREF,SMCWLT,SNOMLT,SOILM,SOILW,Q1,T1
2766 REAL :: SNOTIME1 ! LSTSNW1 INITIAL NUMBER OF TIMESTEPS SINCE LAST SNOWFALL
2770 REAL :: COSZ, SOLARDIRECT
2772 REAL, DIMENSION(1:num_soil_layers):: SLDPTH, STC,SMC,SWC
2774 REAL, DIMENSION(1:num_soil_layers) :: ZSOIL, RTDIS
2775 REAL, PARAMETER :: TRESH=.95E0, A2=17.67,A3=273.15,A4=29.65, &
2776 T0=273.16E0, ELWV=2.50E6, A23M4=A2*(A3-A4)
2778 REAL, PARAMETER :: ROW=1.E3,ELIW=XLF,ROWLIW=ROW*ELIW
2780 ! ----------------------------------------------------------------------
2781 ! DECLARATIONS START - urban
2782 ! ----------------------------------------------------------------------
2784 ! input variables surface_driver --> lsm
2785 INTEGER, INTENT(IN) :: num_roof_layers
2786 INTEGER, INTENT(IN) :: num_wall_layers
2787 INTEGER, INTENT(IN) :: num_road_layers
2788 REAL, OPTIONAL, DIMENSION(1:num_roof_layers), INTENT(IN) :: DZR
2789 REAL, OPTIONAL, DIMENSION(1:num_wall_layers), INTENT(IN) :: DZB
2790 REAL, OPTIONAL, DIMENSION(1:num_road_layers), INTENT(IN) :: DZG
2791 REAL, OPTIONAL, INTENT(IN) :: DECLIN_URB
2792 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: COSZ_URB2D
2793 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: OMG_URB2D
2794 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: XLAT_URB2D
2795 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: U_PHY
2796 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: V_PHY
2797 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: TH_PHY
2798 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: P_PHY
2799 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: RHO
2800 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: UST
2802 LOGICAL, intent(in) :: rdlai2d
2803 LOGICAL, intent(in) :: USEMONALB
2805 ! input variables lsm --> urban
2806 INTEGER :: UTYPE_URB ! urban type [urban=1, suburban=2, rural=3]
2807 REAL :: TA_URB ! potential temp at 1st atmospheric level [K]
2808 REAL :: QA_URB ! mixing ratio at 1st atmospheric level [kg/kg]
2809 REAL :: UA_URB ! wind speed at 1st atmospheric level [m/s]
2810 REAL :: U1_URB ! u at 1st atmospheric level [m/s]
2811 REAL :: V1_URB ! v at 1st atmospheric level [m/s]
2812 REAL :: SSG_URB ! downward total short wave radiation [W/m/m]
2813 REAL :: LLG_URB ! downward long wave radiation [W/m/m]
2814 REAL :: RAIN_URB ! precipitation [mm/h]
2815 REAL :: RHOO_URB ! air density [kg/m^3]
2816 REAL :: ZA_URB ! first atmospheric level [m]
2817 REAL :: DELT_URB ! time step [s]
2818 REAL :: SSGD_URB ! downward direct short wave radiation [W/m/m]
2819 REAL :: SSGQ_URB ! downward diffuse short wave radiation [W/m/m]
2820 REAL :: XLAT_URB ! latitude [deg]
2821 REAL :: COSZ_URB ! cosz
2822 REAL :: OMG_URB ! hour angle
2823 REAL :: ZNT_URB ! roughness length [m]
2834 REAL, DIMENSION(1:num_roof_layers) :: TRL_URB ! roof layer temp [K]
2835 REAL, DIMENSION(1:num_wall_layers) :: TBL_URB ! wall layer temp [K]
2836 REAL, DIMENSION(1:num_road_layers) :: TGL_URB ! road layer temp [K]
2837 LOGICAL :: LSOLAR_URB
2839 !===Yang,2014/10/08,hydrological variable for single layer UCM===
2840 INTEGER :: jmonth, jday, tloc
2841 INTEGER :: IRIOPTION, USOIL, DSOIL
2851 REAL, DIMENSION(1:num_roof_layers) :: SMR_URB ! green roof layer moisture
2852 REAL, DIMENSION(1:num_roof_layers) :: TGRL_URB ! green roof layer temp [K]
2854 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: DRELR_URB2D
2855 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: DRELB_URB2D
2856 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: DRELG_URB2D
2857 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: FLXHUMR_URB2D
2858 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: FLXHUMB_URB2D
2859 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: FLXHUMG_URB2D
2860 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: CMCR_URB2D
2861 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TGR_URB2D
2862 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_roof_layers, jms:jme ), INTENT(INOUT) :: TGRL_URB3D
2863 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_roof_layers, jms:jme ), INTENT(INOUT) :: SMR_URB3D
2865 ! state variable surface_driver <--> lsm <--> urban
2866 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TR_URB2D
2867 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TB_URB2D
2868 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TG_URB2D
2869 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TC_URB2D
2870 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: QC_URB2D
2871 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: UC_URB2D
2872 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXR_URB2D
2873 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXB_URB2D
2874 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXG_URB2D
2875 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXC_URB2D
2876 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: SH_URB2D
2877 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: LH_URB2D
2878 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: G_URB2D
2879 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: RN_URB2D
2881 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TS_URB2D
2883 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_roof_layers, jms:jme ), INTENT(INOUT) :: TRL_URB3D
2884 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_wall_layers, jms:jme ), INTENT(INOUT) :: TBL_URB3D
2885 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_road_layers, jms:jme ), INTENT(INOUT) :: TGL_URB3D
2887 ! output variable lsm --> surface_driver
2888 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: PSIM_URB2D
2889 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: PSIH_URB2D
2890 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: GZ1OZ0_URB2D
2891 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: U10_URB2D
2892 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: V10_URB2D
2893 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: TH2_URB2D
2894 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: Q2_URB2D
2896 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: AKMS_URB2D
2898 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: UST_URB2D
2899 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: FRC_URB2D
2900 INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: UTYPE_URB2D
2902 ! output variables urban --> lsm
2903 REAL :: TS_URB ! surface radiative temperature [K]
2904 REAL :: QS_URB ! surface humidity [-]
2905 REAL :: SH_URB ! sensible heat flux [W/m/m]
2906 REAL :: LH_URB ! latent heat flux [W/m/m]
2907 REAL :: LH_KINEMATIC_URB ! latent heat flux, kinetic [kg/m/m/s]
2908 REAL :: SW_URB ! upward short wave radiation flux [W/m/m]
2909 REAL :: ALB_URB ! time-varying albedo [fraction]
2910 REAL :: LW_URB ! upward long wave radiation flux [W/m/m]
2911 REAL :: G_URB ! heat flux into the ground [W/m/m]
2912 REAL :: RN_URB ! net radiation [W/m/m]
2913 REAL :: PSIM_URB ! shear f for momentum [-]
2914 REAL :: PSIH_URB ! shear f for heat [-]
2915 REAL :: GZ1OZ0_URB ! shear f for heat [-]
2916 REAL :: U10_URB ! wind u component at 10 m [m/s]
2917 REAL :: V10_URB ! wind v component at 10 m [m/s]
2918 REAL :: TH2_URB ! potential temperature at 2 m [K]
2919 REAL :: Q2_URB ! humidity at 2 m [-]
2923 ! NUDAPT Parameters urban --> lam
2928 REAL, DIMENSION(4) :: lf_urb
2929 ! Distributed aerodynamics parameters
2933 ! Variables for multi-layer UCM (Martilli et al. 2002)
2934 REAL, OPTIONAL, INTENT(IN ) :: GMT
2935 INTEGER, OPTIONAL, INTENT(IN ) :: JULDAY
2936 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) ::XLAT, XLONG
2937 INTEGER, INTENT(IN ) :: num_urban_ndm
2938 INTEGER, INTENT(IN ) :: urban_map_zrd
2939 INTEGER, INTENT(IN ) :: urban_map_zwd
2940 INTEGER, INTENT(IN ) :: urban_map_gd
2941 INTEGER, INTENT(IN ) :: urban_map_zd
2942 INTEGER, INTENT(IN ) :: urban_map_zdf
2943 INTEGER, INTENT(IN ) :: urban_map_bd
2944 INTEGER, INTENT(IN ) :: urban_map_wd
2945 INTEGER, INTENT(IN ) :: urban_map_gbd
2946 INTEGER, INTENT(IN ) :: urban_map_fbd
2947 INTEGER, INTENT(IN ) :: urban_map_zgrd
2948 INTEGER, INTENT(IN ) :: NUM_URBAN_HI
2949 INTEGER, INTENT(IN ) :: use_wudapt_lcz
2950 LOGICAL, INTENT(IN ) :: slucm_distributed_drag
2951 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: tsk_rural_bep
2952 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zrd, jms:jme ), INTENT(INOUT) :: trb_urb4d
2953 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zwd, jms:jme ), INTENT(INOUT) :: tw1_urb4d
2954 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zwd, jms:jme ), INTENT(INOUT) :: tw2_urb4d
2955 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_gd , jms:jme ), INTENT(INOUT) :: tgb_urb4d
2956 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_bd , jms:jme ), INTENT(INOUT) :: tlev_urb3d
2957 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_bd , jms:jme ), INTENT(INOUT) :: qlev_urb3d
2958 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_wd , jms:jme ), INTENT(INOUT) :: tw1lev_urb3d
2959 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_wd , jms:jme ), INTENT(INOUT) :: tw2lev_urb3d
2960 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_gbd, jms:jme ), INTENT(INOUT) :: tglev_urb3d
2961 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_fbd, jms:jme ), INTENT(INOUT) :: tflev_urb3d
2962 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: lf_ac_urb3d
2963 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: sf_ac_urb3d
2964 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: cm_ac_urb3d
2965 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: sfvent_urb3d
2966 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: lfvent_urb3d
2967 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_wd , jms:jme ), INTENT(INOUT) :: sfwin1_urb3d
2968 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_wd , jms:jme ), INTENT(INOUT) :: sfwin2_urb3d
2969 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zd , jms:jme ), INTENT(INOUT) :: sfw1_urb3d
2970 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zd , jms:jme ), INTENT(INOUT) :: sfw2_urb3d
2971 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme ), INTENT(INOUT) :: sfr_urb3d
2972 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_ndm, jms:jme ), INTENT(INOUT) :: sfg_urb3d
2973 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: ep_pv_urb3d !GRZ
2974 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme ), INTENT(INOUT) :: t_pv_urb3d !GRZ
2975 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zgrd, jms:jme ),INTENT(INOUT) :: trv_urb4d !GRZ
2976 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zgrd, jms:jme ),INTENT(INOUT) :: qr_urb4d !GRZ
2977 REAL, OPTIONAL, DIMENSION( ims:ime,jms:jme ), INTENT(INOUT) :: qgr_urb3d !GRZ
2978 REAL, OPTIONAL, DIMENSION( ims:ime,jms:jme ), INTENT(INOUT) :: tgr_urb3d !GRZ
2979 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme ),INTENT(INOUT) :: drain_urb4d !GRZ
2980 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: draingr_urb3d !GRZ
2981 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme ),INTENT(INOUT) :: sfrv_urb3d !GRZ
2982 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme ),INTENT(INOUT) :: lfrv_urb3d !GRZ
2983 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme ),INTENT(INOUT) :: dgr_urb3d !GRZ
2984 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_ndm, jms:jme ),INTENT(INOUT) :: dg_urb3d !GRZ
2985 REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zdf, jms:jme ),INTENT(INOUT) :: lfr_urb3d !GRZ
2986 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_ndm, jms:jme ),INTENT(INOUT) :: lfg_urb3d !GRZ
2987 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_urban_hi, jms:jme ), INTENT(IN) :: hi_urb2d
2988 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: lp_urb2d
2989 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: lb_urb2d
2990 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: hgt_urb2d
2991 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: mh_urb2d
2992 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: stdh_urb2d
2993 REAL, OPTIONAL, DIMENSION( ims:ime, 4, jms:jme ), INTENT(IN) :: lf_urb2d
2994 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: lf_urb2d_s
2995 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: z0_urb2d
2996 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_u_bep !Implicit momemtum component X-direction
2997 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_v_bep !Implicit momemtum component Y-direction
2998 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_t_bep !Implicit component pot. temperature
2999 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_q_bep !Implicit momemtum component X-direction
3000 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::a_e_bep !Implicit component TKE
3001 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::b_u_bep !Explicit momentum component X-direction
3002 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::b_v_bep !Explicit momentum component Y-direction
3003 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::b_t_bep !Explicit component pot. temperature
3004 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::b_q_bep !Implicit momemtum component Y-direction
3005 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::b_e_bep !Explicit component TKE
3006 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::vl_bep !Fraction air volume in grid cell
3007 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::dlg_bep !Height above ground
3008 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::sf_bep !Fraction air at the face of grid cell
3009 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) ::dl_u_bep !Length scale
3011 ! Local variables for multi-layer UCM (Martilli et al. 2002)
3012 REAL, DIMENSION( its:ite, jts:jte ) :: HFX_RURAL,LH_RURAL,GRDFLX_RURAL ! ,RN_RURAL
3013 REAL, DIMENSION( its:ite, jts:jte ) :: QFX_RURAL ! ,QSFC_RURAL,UMOM_RURAL,VMOM_RURAL
3014 REAL, DIMENSION( its:ite, jts:jte ) :: ALB_RURAL,EMISS_RURAL,TSK_RURAL ! ,UST_RURAL
3015 ! REAL, DIMENSION( ims:ime, jms:jme ) :: QSFC_URB
3016 REAL, DIMENSION( its:ite, jts:jte ) :: HFX_URB,UMOM_URB,VMOM_URB
3017 REAL, DIMENSION( its:ite, jts:jte ) :: QFX_URB
3018 ! REAL, DIMENSION( ims:ime, jms:jme ) :: ALBEDO_URB,EMISS_URB,UMOM,VMOM,UST
3019 REAL, DIMENSION(its:ite,jts:jte) ::EMISS_URB
3020 REAL, DIMENSION(its:ite,jts:jte) :: RL_UP_URB
3021 REAL, DIMENSION(its:ite,jts:jte) ::RS_ABS_URB
3022 REAL, DIMENSION(its:ite,jts:jte) ::GRDFLX_URB
3023 REAL :: SIGMA_SB,RL_UP_RURAL,RL_UP_TOT,RS_ABS_TOT,UMOM,VMOM
3024 REAL :: CMR_URB, CHR_URB, CMC_URB, CHC_URB, CMGR_URB, CHGR_URB
3025 REAL :: frc_urb,lb_urb
3027 ! ----------------------------------------------------------------------
3028 ! DECLARATIONS END - urban
3029 ! ----------------------------------------------------------------------
3030 !-------------------------------------------------
3031 ! Noah-mosaic related variables are added to declaration (danli)
3032 !-------------------------------------------------
3034 INTEGER, INTENT(IN) :: sf_surface_mosaic
3035 INTEGER, INTENT(IN) :: mosaic_cat, NLCAT
3036 REAL, DIMENSION( ims:ime, NLCAT, jms:jme ), INTENT(IN) :: landusef
3037 REAL, DIMENSION( ims:ime, NLCAT, jms:jme ), INTENT(INOUT) ::landusef2
3038 INTEGER, DIMENSION( ims:ime, NLCAT, jms:jme ), INTENT(INOUT) :: mosaic_cat_index
3040 REAL, DIMENSION( ims:ime, 1:mosaic_cat, jms:jme ) , OPTIONAL, INTENT(INOUT):: &
3041 TSK_mosaic, QSFC_mosaic, CANWAT_mosaic, SNOW_mosaic,SNOWH_mosaic, SNOWC_mosaic
3042 REAL, DIMENSION( ims:ime, 1:mosaic_cat, jms:jme ) , OPTIONAL, INTENT(INOUT):: &
3043 ALBEDO_mosaic,ALBBCK_mosaic, EMISS_mosaic, EMBCK_mosaic, ZNT_mosaic, Z0_mosaic, &
3044 HFX_mosaic,QFX_mosaic, LH_mosaic, GRDFLX_mosaic,SNOTIME_mosaic
3045 REAL, DIMENSION( ims:ime, 1:num_soil_layers*mosaic_cat, jms:jme ), OPTIONAL, INTENT(INOUT):: &
3046 TSLB_mosaic,SMOIS_mosaic,SH2O_mosaic
3047 REAL, DIMENSION( ims:ime, 1:mosaic_cat, jms:jme ), OPTIONAL, INTENT(INOUT):: LAI_mosaic, RC_mosaic
3049 REAL, DIMENSION( ims:ime, jms:jme ) :: TSK_mosaic_avg, QSFC_mosaic_avg, CANWAT_mosaic_avg,SNOW_mosaic_avg,SNOWH_mosaic_avg, &
3050 SNOWC_mosaic_avg, HFX_mosaic_avg, QFX_mosaic_avg, LH_mosaic_avg, GRDFLX_mosaic_avg, &
3051 ALBEDO_mosaic_avg, ALBBCK_mosaic_avg, EMISS_mosaic_avg, EMBCK_mosaic_avg, &
3052 ZNT_mosaic_avg, Z0_mosaic_avg, LAI_mosaic_avg, RC_mosaic_avg, SNOTIME_mosaic_avg, &
3054 REAL, DIMENSION( ims:ime, 1:num_soil_layers, jms:jme ) :: &
3055 TSLB_mosaic_avg,SMOIS_mosaic_avg,SH2O_mosaic_avg
3057 REAL, DIMENSION( ims:ime, 1:mosaic_cat, jms:jme ) , OPTIONAL, INTENT(INOUT):: &
3058 TR_URB2D_mosaic, TB_URB2D_mosaic, TG_URB2D_mosaic, TC_URB2D_mosaic,QC_URB2D_mosaic, UC_URB2D_mosaic, &
3059 SH_URB2D_mosaic,LH_URB2D_mosaic,G_URB2D_mosaic,RN_URB2D_mosaic,TS_URB2D_mosaic, TS_RUL2D_mosaic
3061 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_roof_layers*mosaic_cat, jms:jme ), INTENT(INOUT) :: TRL_URB3D_mosaic
3062 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_wall_layers*mosaic_cat, jms:jme ), INTENT(INOUT) :: TBL_URB3D_mosaic
3063 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_road_layers*mosaic_cat, jms:jme ), INTENT(INOUT) :: TGL_URB3D_mosaic
3065 INTEGER, DIMENSION( ims:ime, jms:jme ) :: IVGTYP_dominant
3066 INTEGER :: mosaic_i, URBAN_METHOD, zo_avg_option
3068 LOGICAL :: IPRINT_mosaic, Noah_call
3069 !-------------------------------------------------
3070 ! Noah-mosaic related variables declaration end (danli)
3071 !-------------------------------------------------
3073 REAL, PARAMETER :: CAPA=R_D/CP
3074 REAL :: APELM,APES,SFCTH2,PSFC
3075 real, intent(in) :: xice_threshold
3076 character(len=80) :: message_text
3078 ! FASDAS: it doesn't work for mosaic, but we need the variables to call sflx
3080 REAL, DIMENSION( ims:ime, jms:jme ) , &
3081 INTENT(INOUT) :: SDA_HFX, SDA_QFX, HFX_BOTH, QFX_BOTH, QNORM
3082 INTEGER, INTENT(IN ) :: fasdas
3083 REAL :: XSDA_HFX, XSDA_QFX, XQNORM
3084 REAL :: HFX_PHY, QFX_PHY
3086 REAL :: HCPCT_FASDAS
3087 REAL,OPTIONAL,DIMENSION( ims:ime, jms:jme ) :: IRR_CHAN
3088 REAL :: IRRIGATION_CHANNEL
3089 IRRIGATION_CHANNEL=0.0
3090 HFX_PHY = 0.0 ! initialize
3101 FDTW=DT/(XLV*RHOWATER)
3104 IPRINT_mosaic=.false.
3110 NSOIL=num_soil_layers
3116 JLOOP : DO J=jts,jte
3118 IF(ITIMESTEP.EQ.1)THEN
3120 !*** initialize soil conditions for IHOP 31 May case
3121 ! IF((XLAND(I,J)-1.5) < 0.)THEN
3122 ! if (I==108.and.j==85) then
3124 ! SMOIS(I,NS,J)=0.10
3130 !*** SET ZERO-VALUE FOR SOME OUTPUT DIAGNOSTIC ARRAYS
3131 IF((XLAND(I,J)-1.5).GE.0.)THEN
3132 ! check sea-ice point
3134 IF( XICE(I,J).GE. XICE_THRESHOLD .and. IPRINT ) PRINT*, ' sea-ice at water point, I=',I,'J=',J
3136 !*** Open Water Case
3141 TSLB(I,NS,J)=273.16 !STEMP
3145 IF ( XICE(I,J) .GE. XICE_THRESHOLD ) THEN
3157 ENDIF ! end of initialization over ocean
3159 !-----------------------------------------------------------------------
3160 ILOOP : DO I=its,ite
3162 ! initializing local variables
3169 IF (((XLAND(I,J)-1.5).LT.0.) .AND. (XICE(I,J) < XICE_THRESHOLD) ) THEN
3171 IVGTYP_dominant(I,J)=IVGTYP(I,J) ! save this
3173 ! INITIALIZE THE AREA-AVERAGED FLUXES
3175 TSK_mosaic_avg(i,j)= 0.0 ! from 3D to 2D
3176 QSFC_mosaic_avg(i,j)= 0.0
3177 CANWAT_mosaic_avg(i,j)= 0.0
3178 SNOW_mosaic_avg(i,j)= 0.0
3179 SNOWH_mosaic_avg(i,j)= 0.0
3180 SNOWC_mosaic_avg(i,j)= 0.0
3184 TSLB_mosaic_avg(i,NS,j)=0.0
3185 SMOIS_mosaic_avg(i,NS,j)=0.0
3186 SH2O_mosaic_avg(i,NS,j)=0.0
3190 HFX_mosaic_avg(i,j)= 0.0
3191 QFX_mosaic_avg(i,j)= 0.0
3192 LH_mosaic_avg(i,j)= 0.0
3193 GRDFLX_mosaic_avg(i,j)= 0.0
3194 ALBEDO_mosaic_avg(i,j)=0.0
3195 ALBBCK_mosaic_avg(i,j)=0.0
3196 EMISS_mosaic_avg(i,j)=0.0
3197 EMBCK_mosaic_avg(i,j)=0.0
3198 ZNT_mosaic_avg(i,j)=0.0
3199 Z0_mosaic_avg(i,j)=0.0
3200 LAI_mosaic_avg(i,j)=0.0
3201 RC_mosaic_avg(i,j)=0.0
3202 FAREA_mosaic_avg(i,j)=0.0
3204 ! add a new loop for the mosaic_cat
3206 DO mosaic_i = mosaic_cat, 1, -1
3208 ! if (mosaic_cat_index(I,mosaic_i,J) .EQ. 16 ) then
3209 ! PRINT*, 'you still have water tiles at','i=',i,'j=',j, 'mosaic_i',mosaic_i
3210 ! PRINT*, 'xland',xland(i,j),'xice',xice(i,j)
3213 IVGTYP(I,J)=mosaic_cat_index(I,mosaic_i,J) ! replace it with the mosaic one
3214 TSK(I,J)=TSK_mosaic(I,mosaic_i,J) ! from 3D to 2D
3215 QSFC(i,j)=QSFC_mosaic(I,mosaic_i,J)
3216 CANWAT(i,j)=CANWAT_mosaic(i,mosaic_i,j)
3217 SNOW(i,j)=SNOW_mosaic(i,mosaic_i,j)
3218 SNOWH(i,j)=SNOWH_mosaic(i,mosaic_i,j)
3219 SNOWC(i,j)=SNOWC_mosaic(i,mosaic_i,j)
3221 ALBEDO(i,j) = ALBEDO_mosaic(i,mosaic_i,j)
3222 ALBBCK(i,j)= ALBBCK_mosaic(i,mosaic_i,j)
3223 EMISS(i,j)= EMISS_mosaic(i,mosaic_i,j)
3224 EMBCK(i,j)= EMBCK_mosaic(i,mosaic_i,j)
3225 ZNT(i,j)= ZNT_mosaic(i,mosaic_i,j)
3226 Z0(i,j)= Z0_mosaic(i,mosaic_i,j)
3228 SNOTIME(i,j)= SNOTIME_mosaic(i,mosaic_i,j)
3232 TSLB(i,NS,j)=TSLB_mosaic(i,NSOIL*(mosaic_i-1)+NS,j)
3233 SMOIS(i,NS,j)=SMOIS_mosaic(i,NSOIL*(mosaic_i-1)+NS,j)
3234 SH2O(i,NS,j)=SH2O_mosaic(i,NSOIL*(mosaic_i-1)+NS,j)
3238 IF(IPRINT_mosaic) THEN
3240 print*, 'BEFORE SFLX, in Noahdrv.F'
3241 print*, 'mosaic_cat', mosaic_cat, 'IVGTYP',IVGTYP(i,j), 'TSK',TSK(i,j),'HFX',HFX(i,j), 'QSFC', QSFC(i,j), &
3242 '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)
3246 !-----------------------------------------------------------------------
3247 ! insert the NOAH model here for the non-urban one and the urban one DANLI
3248 !-----------------------------------------------------------------------
3252 ! pressure in middle of lowest layer
3253 SFCPRS=(P8W3D(I,KTS+1,j)+P8W3D(i,KTS,j))*0.5
3254 ! convert from mixing ratio to specific humidity
3255 Q2K=QV3D(i,1,j)/(1.0+QV3D(i,1,j))
3258 Q2SAT=QGH(I,J)/(1.0+QGH(I,J)) ! Q2SAT is sp humidity
3259 ! add check on myj=.true.
3260 ! IF((Q2K.GE.Q2SAT*TRESH).AND.Q2K.LT.QZ0(I,J))THEN
3261 IF((myj).AND.(Q2K.GE.Q2SAT*TRESH).AND.Q2K.LT.QZ0(I,J))THEN
3270 ZLVL=0.5*DZ8W(i,1,j)
3272 ! TH2=SFCTMP+(0.0097545*ZLVL)
3273 ! calculate SFCTH2 via Exner function vs lapse-rate (above)
3274 APES=(1.E5/PSFC)**CAPA
3275 APELM=(1.E5/SFCPRS)**CAPA
3280 LWDN=GLW(I,J)*EMISSI
3281 ! SOLDN is total incoming solar
3283 ! GSW is net downward solar
3285 ! use mid-day albedo to determine net downward solar (no solar zenith angle correction)
3286 SOLNET=SOLDN*(1.-ALBEDO(I,J))
3288 IF(PRESENT(IRR_CHAN)) THEN
3289 IF(IRR_CHAN(i,j).NE.0) THEN
3290 IRRIGATION_CHANNEL=IRR_CHAN(i,j)/DT
3292 IRRIGATION_CHANNEL=0.
3297 SHDFAC=VEGFRA(I,J)/100.
3300 SHMIN=SHDMIN(I,J)/100. !NEW
3301 SHMAX=SHDMAX(I,J)/100. !NEW
3302 ! convert snow water equivalent from mm to meter
3303 SNEQV=SNOW(I,J)*0.001
3304 ! snow depth in meters
3308 ! if "SR" present, set frac of frozen precip ("FFROZP") = snow-ratio ("SR", range:0-1)
3309 ! SR from e.g. Ferrier microphysics
3310 ! otherwise define from 1st atmos level temperature
3314 IF (SFCTMP <= 273.15) THEN
3321 IF((XLAND(I,J)-1.5).GE.0.)THEN ! begining of land/sea if block
3323 TSK_RURAL(I,J)=TSK(I,J)
3324 HFX_RURAL(I,J)=HFX(I,J)
3325 QFX_RURAL(I,J)=QFX(I,J)
3326 LH_RURAL(I,J)=LH(I,J)
3327 EMISS_RURAL(I,J)=EMISS(I,J)
3328 GRDFLX_RURAL(I,J)=GRDFLX(I,J)
3330 ! Land or sea-ice case
3332 IF (XICE(I,J) >= XICE_THRESHOLD) THEN
3335 ELSE IF ( VEGTYP == ISICE ) THEN
3339 ! Neither sea ice or land ice.
3342 DQSDT2=Q2SAT*A23M4/(SFCTMP-A4)**2
3344 IF(SNOW(I,J).GT.0.0)THEN
3345 ! snow on surface (use ice saturation properties)
3347 E2SAT=611.2*EXP(6174.*(1./273.15 - 1./SFCTSNO))
3348 Q2SATI=0.622*E2SAT/(SFCPRS-E2SAT)
3349 Q2SATI=Q2SATI/(1.0+Q2SATI) ! spec. hum.
3350 IF (T1 .GT. 273.14) THEN
3351 ! warm ground temps, weight the saturation between ice and water according to SNOWC
3352 Q2SAT=Q2SAT*(1.-SNOWC(I,J)) + Q2SATI*SNOWC(I,J)
3353 DQSDT2=DQSDT2*(1.-SNOWC(I,J)) + Q2SATI*6174./(SFCTSNO**2)*SNOWC(I,J)
3355 ! cold ground temps, use ice saturation only
3357 DQSDT2=Q2SATI*6174./(SFCTSNO**2)
3359 ! for snow cover fraction at 0 C, ground temp will not change, so DQSDT2 effectively zero
3360 IF(T1 .GT. 273. .AND. SNOWC(I,J) .GT. 0.)DQSDT2=DQSDT2*(1.-SNOWC(I,J))
3363 ! Land-ice or land points use the usual deep-soil temperature.
3366 IF(VEGTYP.EQ.25) SHDFAC=0.0000
3367 IF(VEGTYP.EQ.26) SHDFAC=0.0000
3368 IF(VEGTYP.EQ.27) SHDFAC=0.0000
3369 IF(SOILTYP.EQ.14.AND.XICE(I,J).EQ.0.)THEN
3371 IF(IPRINT)PRINT*,' SOIL TYPE FOUND TO BE WATER AT A LAND-POINT'
3372 IF(IPRINT)PRINT*,i,j,'RESET SOIL in surfce.F'
3376 SNOALB1 = SNOALB(I,J)
3377 ! converts canwat in mm to CMC in meters
3378 CMC=CANWAT(I,J)/1000.
3380 !-------------------------------------------
3381 !*** convert snow depth from mm to meter
3384 ! SNOALB=ALBMAX(I,J)*0.01
3386 ! SNOALB=MAXALB(IVGTPK)*0.01
3394 SNOTIME1 = SNOTIME(I,J)
3396 !FEI: temporaray arrays above need to be changed later by using SI
3399 SMC(NS)=SMOIS(I,NS,J)
3400 STC(NS)=TSLB(I,NS,J) !STEMP
3401 SWC(NS)=SH2O(I,NS,J)
3404 if ( (SNEQV.ne.0..AND.SNOWHK.eq.0.).or.(SNOWHK.le.SNEQV) )THEN
3409 !Fei: urban. for urban surface, if calling UCM, redefine the natural surface in cities as
3410 ! the "NATURAL" category in the VEGPARM.TBL
3412 ! IF(SF_URBAN_PHYSICS == 1.OR. SF_URBAN_PHYSICS==2.OR.SF_URBAN_PHYSICS==3 ) THEN
3414 ! IF( IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == LOW_DENSITY_RESIDENTIAL .or. &
3415 ! IVGTYP(I,J) == HIGH_DENSITY_RESIDENTIAL .or. IVGTYP(I,J) == HIGH_INTENSITY_INDUSTRIAL) THEN
3417 ! SHDFAC = SHDTBL(NATURAL)
3418 ! ALBEDOK =0.2 ! 0.2
3420 ! EMISSI = 0.98 !for VEGTYP=5
3421 ! IF ( FRC_URB2D(I,J) < 0.99 ) THEN
3422 ! if(sf_urban_physics.eq.1)then
3423 ! T1= ( TSK(I,J) -FRC_URB2D(I,J) * TS_URB2D (I,J) )/ (1-FRC_URB2D(I,J))
3424 ! elseif((sf_urban_physics.eq.2).OR.(sf_urban_physics.eq.3))then
3425 ! r1= (tsk(i,j)**4.)
3426 ! r2= frc_urb2d(i,j)*(ts_urb2d(i,j)**4.)
3427 ! r3= (1.-frc_urb2d(i,j))
3428 ! t1= ((r1-r2)/r3)**.25
3435 ! IF( IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == LOW_DENSITY_RESIDENTIAL .or. &
3436 ! IVGTYP(I,J) == HIGH_DENSITY_RESIDENTIAL .or. IVGTYP(I,J) == HIGH_INTENSITY_INDUSTRIAL) THEN
3443 If ( SF_URBAN_PHYSICS == 0 ) THEN ! ONLY NOAH
3444 IF( IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == LCZ_1 .or. IVGTYP(I,J) == LCZ_2 .or. &
3445 IVGTYP(I,J) == LCZ_3 .or. IVGTYP(I,J) == LCZ_4 .or. IVGTYP(I,J) == LCZ_5 .or. &
3446 IVGTYP(I,J) == LCZ_6 .or. IVGTYP(I,J) == LCZ_7 .or. IVGTYP(I,J) == LCZ_8 .or. &
3447 IVGTYP(I,J) == LCZ_9 .or. IVGTYP(I,J) == LCZ_10 .or. IVGTYP(I,J) == LCZ_11 ) THEN
3455 IF(SF_URBAN_PHYSICS == 1) THEN
3457 IF( IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == LCZ_1 .or. IVGTYP(I,J) == LCZ_2 .or. &
3458 IVGTYP(I,J) == LCZ_3 .or. IVGTYP(I,J) == LCZ_4 .or. IVGTYP(I,J) == LCZ_5 .or. &
3459 IVGTYP(I,J) == LCZ_6 .or. IVGTYP(I,J) == LCZ_7 .or. IVGTYP(I,J) == LCZ_8 .or. &
3460 IVGTYP(I,J) == LCZ_9 .or. IVGTYP(I,J) == LCZ_10 .or. IVGTYP(I,J) == LCZ_11 ) THEN
3463 SHDFAC = SHDTBL(NATURAL)
3466 EMISSI = 0.98 ! for VEGTYP=5
3467 LWDN = GLW(I,J) * EMISSI
3468 SOLNET = SOLDN * (1.0 - ALBEDOK)
3470 T1= TS_RUL2D_mosaic(I,mosaic_i,J)
3476 !===Yang, 2014/10/08, hydrological processes for urban vegetation in single layer UCM===
3480 IRIOPTION=IRI_SCHEME
3482 tloc=mod(int(OMG/3.14159*180./15.+12.+0.5 ),24)
3483 if (tloc.lt.0) tloc=tloc+24
3484 if (tloc==0) tloc=24
3485 CALL cal_mon_day(julian,julyr,jmonth,jday)
3486 IF(SF_URBAN_PHYSICS == 1) THEN
3487 IF( IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == LCZ_1 .or. IVGTYP(I,J) == LCZ_2 .or. &
3488 IVGTYP(I,J) == LCZ_3 .or. IVGTYP(I,J) == LCZ_4 .or. IVGTYP(I,J) == LCZ_5 .or. &
3489 IVGTYP(I,J) == LCZ_6 .or. IVGTYP(I,J) == LCZ_7 .or. IVGTYP(I,J) == LCZ_8 .or. &
3490 IVGTYP(I,J) == LCZ_9 .or. IVGTYP(I,J) == LCZ_10 .or. IVGTYP(I,J) == LCZ_11 ) THEN
3492 AOASIS = oasis ! urban oasis effect
3493 IF (IRIOPTION ==1) THEN
3494 IF (tloc==21 .or. tloc==22) THEN !irrigation on vegetaion in urban area, MAY-SEP, 9-10pm
3495 IF (jmonth==5 .or. jmonth==6 .or. jmonth==7 .or. jmonth==8 .or. jmonth==9) THEN
3496 IF (SMC(USOIL) .LT. SMCREF) SMC(USOIL)= REFSMC(ISLTYP(I,J))
3497 IF (SMC(DSOIL) .LT. SMCREF) SMC(DSOIL)= REFSMC(ISLTYP(I,J))
3504 IF(SF_URBAN_PHYSICS == 2 .or. SF_URBAN_PHYSICS == 3) THEN
3505 IF(AOASIS > 1.0) THEN
3506 CALL wrf_error_fatal('Urban oasis option is for SF_URBAN_PHYSICS == 1 only')
3508 IF(IRIOPTION == 1) THEN
3509 CALL wrf_error_fatal('Urban irrigation option is for SF_URBAN_PHYSICS == 1 only')
3513 IF( SF_URBAN_PHYSICS==2.OR.SF_URBAN_PHYSICS==3 ) THEN
3514 ! print*, 'MOSAIC is not designed to work with SF_URBAN_PHYSICS=2 or SF_URBAN_PHYSICS=3'
3521 print*, 'BEFORE SFLX, in Noahlsm_driver'
3522 print*, 'ICE', ICE, 'DT',DT, 'ZLVL',ZLVL, 'NSOIL', NSOIL, &
3523 'SLDPTH', SLDPTH, 'LOCAL',LOCAL, 'LUTYPE',&
3524 LUTYPE, 'SLTYPE',SLTYPE, 'LWDN',LWDN, 'SOLDN',SOLDN, &
3525 'SFCPRS',SFCPRS, 'PRCP',PRCP,'SFCTMP',SFCTMP,'Q2K',Q2K, &
3526 'TH2',TH2,'Q2SAT',Q2SAT,'DQSDT2',DQSDT2,'VEGTYP', VEGTYP,&
3527 'SOILTYP',SOILTYP, 'SLOPETYP',SLOPETYP, 'SHDFAC',SHDFAC,&
3528 'SHMIN',SHMIN, 'ALBBRD',ALBBRD,'SNOALB1',SNOALB1,'TBOT',&
3529 TBOT, 'Z0BRD',Z0BRD, 'Z0K',Z0K, 'CMC',CMC, 'T1',T1,'STC',&
3530 STC, 'SMC',SMC, 'SWC',SWC,'SNOWHK',SNOWHK,'SNEQV',SNEQV,&
3531 'ALBEDOK',ALBEDOK,'CHK',CHK,'ETA',ETA,'SHEAT',SHEAT, &
3532 'ETA_KINEMATIC',ETA_KINEMATIC, 'FDOWN',FDOWN,'EC',EC, &
3533 'EDIR',EDIR,'ET',ET,'ETT',ETT,'ESNOW',ESNOW,'DRIP',DRIP,&
3534 'DEW',DEW,'BETA',BETA,'ETP',ETP,'SSOIL',SSOIL,'FLX1',FLX1,&
3535 'FLX2',FLX2,'FLX3',FLX3,'SNOMLT',SNOMLT,'SNCOVR',SNCOVR,&
3536 'RUNOFF1',RUNOFF1,'RUNOFF2',RUNOFF2,'RUNOFF3',RUNOFF3, &
3537 'RC',RC, 'PC',PC,'RSMIN',RSMIN,'XLAI',XLAI,'RCS',RCS, &
3538 'RCT',RCT,'RCQ',RCQ,'RCSOIL',RCSOIL,'SOILW',SOILW, &
3539 'SOILM',SOILM,'Q1',Q1,'SMCWLT',SMCWLT,'SMCDRY',SMCDRY,&
3540 'SMCREF',SMCREF,'SMCMAX',SMCMAX,'NROOT',NROOT
3545 IF (SHDFAC > 0.0 .AND. LAI(I,J) <= 0.0) LAI(I,J) = 0.01
3549 IF ( ICE == 1 ) THEN
3560 ELSEIF (ICE == 0) THEN
3564 CALL SFLX (I,J,FFROZP, ISURBAN, DT,ZLVL,NSOIL,SLDPTH, & !C
3566 LUTYPE, SLTYPE, & !CL
3567 LWDN,SOLDN,SOLNET,SFCPRS,PRCP,SFCTMP,Q2K,DUMMY, & !F
3568 DUMMY,DUMMY, DUMMY, & !F PRCPRAIN not used
3569 TH2,Q2SAT,DQSDT2, & !I
3570 VEGTYP,SOILTYP,SLOPETYP,SHDFAC,SHMIN,SHMAX, & !I
3571 ALBBRD, SNOALB1,TBOT, Z0BRD, Z0K, EMISSI, EMBRD, & !S
3572 CMC,T1,STC,SMC,SWC,SNOWHK,SNEQV,ALBEDOK,CHK,dummy,& !H
3573 ETA,SHEAT, ETA_KINEMATIC,FDOWN, & !O
3574 EC,EDIR,ET,ETT,ESNOW,DRIP,DEW, & !O
3575 BETA,ETP,SSOIL, & !O
3576 FLX1,FLX2,FLX3, & !O
3577 FLX4,FVB,FBUR,FGSN,UA_PHYS, & !UA
3579 RUNOFF1,RUNOFF2,RUNOFF3, & !O
3580 RC,PC,RSMIN,XLAI,RCS,RCT,RCQ,RCSOIL, & !O
3581 SOILW,SOILM,Q1,SMAV, & !D
3582 RDLAI2D,USEMONALB, &
3585 SMCWLT,SMCDRY,SMCREF,SMCMAX,NROOT, &
3587 sfcheadrt_hydro, & !I
3588 INFXSRT_hydro,ETPND1_hydro & !O
3589 ,OPT_THCND,AOASIS & !O
3590 ,XSDA_QFX, HFX_PHY, QFX_PHY, XQNORM, fasdas, HCPCT_FASDAS & ! fasdas vars
3591 ,IRRIGATION_CHANNEL )
3594 soldrain(i,j) = RUNOFF2*DT*1000.0
3595 sfcheadrt(i,j) = sfcheadrt_hydro
3596 infxsrt(i,j) = INFXSRT_hydro
3597 etpnd1 = etpnd1_hydro
3599 ELSEIF (ICE == -1) THEN
3602 ! Set values that the LSM is expected to update,
3603 ! but don't get updated for glacial points.
3605 SOILM = 0.0 !BSINGH(PNNL)- SOILM is undefined for this case, it is used for diagnostics so setting it to zero
3606 XLAI = 0.01 ! KWM Should this be Zero over land ice? Does this value matter?
3614 CALL SFLX_GLACIAL(I,J,ISICE,FFROZP,DT,ZLVL,NSOIL,SLDPTH, & !C
3615 & LWDN,SOLNET,SFCPRS,PRCP,SFCTMP,Q2K, & !F
3616 & TH2,Q2SAT,DQSDT2, & !I
3617 & ALBBRD, SNOALB1,TBOT, Z0BRD, Z0K, EMISSI, EMBRD, & !S
3618 & T1,STC(1:NSOIL),SNOWHK,SNEQV,ALBEDOK,CHK, & !H
3619 & ETA,SHEAT,ETA_KINEMATIC,FDOWN, & !O
3622 & FLX1,FLX2,FLX3, & !O
3623 & SNOMLT,SNCOVR, & !O
3634 print*, 'AFTER SFLX, in Noahlsm_driver'
3635 print*, 'ICE', ICE, 'DT',DT, 'ZLVL',ZLVL, 'NSOIL', NSOIL, &
3636 'SLDPTH', SLDPTH, 'LOCAL',LOCAL, 'LUTYPE',&
3637 LUTYPE, 'SLTYPE',SLTYPE, 'LWDN',LWDN, 'SOLDN',SOLDN, &
3638 'SFCPRS',SFCPRS, 'PRCP',PRCP,'SFCTMP',SFCTMP,'Q2K',Q2K, &
3639 'TH2',TH2,'Q2SAT',Q2SAT,'DQSDT2',DQSDT2,'VEGTYP', VEGTYP,&
3640 'SOILTYP',SOILTYP, 'SLOPETYP',SLOPETYP, 'SHDFAC',SHDFAC,&
3641 'SHDMIN',SHMIN, 'ALBBRD',ALBBRD,'SNOALB',SNOALB1,'TBOT',&
3642 TBOT, 'Z0BRD',Z0BRD, 'Z0K',Z0K, 'CMC',CMC, 'T1',T1,'STC',&
3643 STC, 'SMC',SMC, 'SWc',SWC,'SNOWHK',SNOWHK,'SNEQV',SNEQV,&
3644 'ALBEDOK',ALBEDOK,'CHK',CHK,'ETA',ETA,'SHEAT',SHEAT, &
3645 'ETA_KINEMATIC',ETA_KINEMATIC, 'FDOWN',FDOWN,'EC',EC, &
3646 'EDIR',EDIR,'ET',ET,'ETT',ETT,'ESNOW',ESNOW,'DRIP',DRIP,&
3647 'DEW',DEW,'BETA',BETA,'ETP',ETP,'SSOIL',SSOIL,'FLX1',FLX1,&
3648 'FLX2',FLX2,'FLX3',FLX3,'SNOMLT',SNOMLT,'SNCOVR',SNCOVR,&
3649 'RUNOFF1',RUNOFF1,'RUNOFF2',RUNOFF2,'RUNOFF3',RUNOFF3, &
3650 'RC',RC, 'PC',PC,'RSMIN',RSMIN,'XLAI',XLAI,'RCS',RCS, &
3651 'RCT',RCT,'RCQ',RCQ,'RCSOIL',RCSOIL,'SOILW',SOILW, &
3652 'SOILM',SOILM,'Q1',Q1,'SMCWLT',SMCWLT,'SMCDRY',SMCDRY,&
3653 'SMCREF',SMCREF,'SMCMAX',SMCMAX,'NROOT',NROOT
3657 !*** UPDATE STATE VARIABLES
3658 CANWAT(I,J)=CMC*1000.
3659 SNOW(I,J)=SNEQV*1000.
3660 ! SNOWH(I,J)=SNOWHK*1000.
3661 SNOWH(I,J)=SNOWHK ! SNOWHK in meters
3663 ALB_RURAL(I,J)=ALBEDOK
3667 EMISS_RURAL(I,J) = EMISSI
3668 ! Noah: activate time-varying roughness length (V3.3 Feb 2011)
3673 HFX_RURAL(I,J)=SHEAT
3674 ! MEk Jul07 add potential evap accum
3675 POTEVP(I,J)=POTEVP(I,J)+ETP*FDTW
3676 QFX(I,J)=ETA_KINEMATIC
3677 QFX_RURAL(I,J)=ETA_KINEMATIC
3681 ! QFX(I,J) = QFX(I,J) + ETPND1
3682 ! ETA = ETA + ETPND1/2.501E6*dt
3683 !end added by Wei Yu
3689 GRDFLX_RURAL(I,J)=SSOIL
3692 SNOTIME(I,J) = SNOTIME1
3693 ! prevent diagnostic ground q (q1) from being greater than qsat(tsk)
3694 ! as happens over snow cover where the cqs2 value also becomes irrelevant
3695 ! by setting cqs2=chs in this situation the 2m q should become just qv(k=1)
3696 IF (Q1 .GT. QSFC(I,J)) THEN
3697 CQS2(I,J) = CHS(I,J)
3700 ! Convert QSFC back to mixing ratio
3701 QSFC(I,J)= Q1/(1.0-Q1)
3703 ! QSFC_RURAL(I,J)= Q1/(1.0-Q1)
3704 ! Calculate momentum flux from rural surface for use with multi-layer UCM (Martilli et al. 2002)
3707 SMOIS(I,NS,J)=SMC(NS)
3708 TSLB(I,NS,J)=STC(NS) ! STEMP
3709 SH2O(I,NS,J)=SWC(NS)
3719 ! Residual of surface energy balance equation terms
3723 noahres(i,j) = ( solnet + lwdn ) - sheat + ssoil - eta &
3724 - ( emissi * STBOLT * (t1**4) ) - flx1 - flx2 - flx3 - flx4
3727 noahres(i,j) = ( solnet + lwdn ) - sheat + ssoil - eta &
3728 - ( emissi * STBOLT * (t1**4) ) - flx1 - flx2 - flx3
3731 ENDIF !ENDIF FOR Noah_call
3733 IF (SF_URBAN_PHYSICS == 1 ) THEN ! Beginning of UCM CALL if block
3734 !--------------------------------------
3735 ! URBAN CANOPY MODEL START - urban
3736 !--------------------------------------
3737 ! Input variables lsm --> urban
3738 IF( IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == LCZ_1 .or. IVGTYP(I,J) == LCZ_2 .or. &
3739 IVGTYP(I,J) == LCZ_3 .or. IVGTYP(I,J) == LCZ_4 .or. IVGTYP(I,J) == LCZ_5 .or. &
3740 IVGTYP(I,J) == LCZ_6 .or. IVGTYP(I,J) == LCZ_7 .or. IVGTYP(I,J) == LCZ_8 .or. &
3741 IVGTYP(I,J) == LCZ_9 .or. IVGTYP(I,J) == LCZ_10 .or. IVGTYP(I,J) == LCZ_11 ) THEN
3743 ! UTYPE_URB = UTYPE_URB2D(I,J) !urban type (low, high or industrial)
3744 ! this need to be changed in the mosaic danli
3745 IF (slucm_distributed_drag) THEN
3746 IF (IVGTYP(I, J) == ISURBAN) THEN
3749 ELSE IF (use_wudapt_lcz == 1) THEN
3750 IF(IVGTYP(I,J)==ISURBAN) UTYPE_URB=5
3751 IF(IVGTYP(I,J)==LCZ_1) UTYPE_URB=1
3752 IF(IVGTYP(I,J)==LCZ_2) UTYPE_URB=2
3753 IF(IVGTYP(I,J)==LCZ_3) UTYPE_URB=3
3754 IF(IVGTYP(I,J)==LCZ_4) UTYPE_URB=4
3755 IF(IVGTYP(I,J)==LCZ_5) UTYPE_URB=5
3756 IF(IVGTYP(I,J)==LCZ_6) UTYPE_URB=6
3757 IF(IVGTYP(I,J)==LCZ_7) UTYPE_URB=7
3758 IF(IVGTYP(I,J)==LCZ_8) UTYPE_URB=8
3759 IF(IVGTYP(I,J)==LCZ_9) UTYPE_URB=9
3760 IF(IVGTYP(I,J)==LCZ_10) UTYPE_URB=10
3761 IF(IVGTYP(I,J)==LCZ_11) UTYPE_URB=11
3763 IF(IVGTYP(I,J)==ISURBAN) UTYPE_URB=2
3764 IF(IVGTYP(I,J)==LCZ_1) UTYPE_URB=1 ! LOW_DENSITY_RESIDENTIAL
3765 IF(IVGTYP(I,J)==LCZ_2) UTYPE_URB=2 ! HIGH_DENSITY_RESIDENTIAL
3766 IF(IVGTYP(I,J)==LCZ_3) UTYPE_URB=3 ! HIGH_INTENSITY_INDUSTRIAL
3769 TA_URB = SFCTMP ! [K]
3770 QA_URB = Q2K ! [kg/kg]
3771 UA_URB = SQRT(U_PHY(I,1,J)**2.+V_PHY(I,1,J)**2.)
3772 U1_URB = U_PHY(I,1,J)
3773 V1_URB = V_PHY(I,1,J)
3774 IF(UA_URB < 1.) UA_URB=1. ! [m/s]
3775 SSG_URB = SOLDN ! [W/m/m]
3776 SSGD_URB = 0.8*SOLDN ! [W/m/m]
3777 SSGQ_URB = SSG_URB-SSGD_URB ! [W/m/m]
3778 LLG_URB = GLW(I,J) ! [W/m/m]
3779 RAIN_URB = RAINBL(I,J) / DT * 3600.0 ! [mm/hr]
3780 RHOO_URB = SFCPRS / (287.04 * SFCTMP * (1.0+ 0.61 * Q2K)) ![kg/m/m/m]
3782 DELT_URB = DT ! [sec]
3783 XLAT_URB = XLAT_URB2D(I,J) ! [deg]
3784 COSZ_URB = COSZ_URB2D(I,J) !
3785 OMG_URB = OMG_URB2D(I,J) !
3788 LSOLAR_URB = .FALSE.
3792 TR_URB2D(I,J)=TR_URB2D_mosaic(I,mosaic_i,J) ! replace it with the mosaic one
3793 TB_URB2D(I,J)=TB_URB2D_mosaic(I,mosaic_i,J) ! replace it with the mosaic one
3794 TG_URB2D(I,J)=TG_URB2D_mosaic(I,mosaic_i,J) ! replace it with the mosaic one
3795 TC_URB2D(I,J)=TC_URB2D_mosaic(I,mosaic_i,J) ! replace it with the mosaic one
3796 QC_URB2D(I,J)=QC_URB2D_mosaic(I,mosaic_i,J) ! replace it with the mosaic one
3797 UC_URB2D(I,J)=UC_URB2D_mosaic(I,mosaic_i,J) ! replace it with the mosaic one
3798 TS_URB2D(I,J)=TS_URB2D_mosaic(I,mosaic_i,J) ! replace it with the mosaic one
3800 DO K = 1,num_roof_layers
3801 TRL_URB3D(I,K,J) = TRL_URB3D_mosaic(I,K+(mosaic_i-1)*num_roof_layers,J)
3803 DO K = 1,num_wall_layers
3804 TBL_URB3D(I,K,J) = TBL_URB3D_mosaic(I,K+(mosaic_i-1)*num_roof_layers,J)
3806 DO K = 1,num_road_layers
3807 TGL_URB3D(I,K,J) = TGL_URB3D_mosaic(I,K+(mosaic_i-1)*num_roof_layers,J)
3812 TR_URB = TR_URB2D(I,J)
3813 TB_URB = TB_URB2D(I,J)
3814 TG_URB = TG_URB2D(I,J)
3815 TC_URB = TC_URB2D(I,J)
3816 QC_URB = QC_URB2D(I,J)
3817 UC_URB = UC_URB2D(I,J)
3819 DO K = 1,num_roof_layers
3820 TRL_URB(K) = TRL_URB3D(I,K,J)
3821 SMR_URB(K) = SMR_URB3D(I,K,J)
3822 TGRL_URB(K)= TGRL_URB3D(I,K,J)
3824 DO K = 1,num_wall_layers
3825 TBL_URB(K) = TBL_URB3D(I,K,J)
3827 DO K = 1,num_road_layers
3828 TGL_URB(K) = TGL_URB3D(I,K,J)
3831 TGR_URB = TGR_URB2D(I,J)
3832 CMCR_URB = CMCR_URB2D(I,J)
3833 FLXHUMR_URB = FLXHUMR_URB2D(I,J)
3834 FLXHUMB_URB = FLXHUMB_URB2D(I,J)
3835 FLXHUMG_URB = FLXHUMG_URB2D(I,J)
3836 DRELR_URB = DRELR_URB2D(I,J)
3837 DRELB_URB = DRELB_URB2D(I,J)
3838 DRELG_URB = DRELG_URB2D(I,J)
3840 XXXR_URB = XXXR_URB2D(I,J)
3841 XXXB_URB = XXXB_URB2D(I,J)
3842 XXXG_URB = XXXG_URB2D(I,J)
3843 XXXC_URB = XXXC_URB2D(I,J)
3845 ! Limits to avoid dividing by small number
3846 if (CHS(I,J) < 1.0E-02) then
3849 if (CHS2(I,J) < 1.0E-02) then
3852 if (CQS2(I,J) < 1.0E-02) then
3857 CHS2_URB = CHS2(I,J)
3858 IF (PRESENT(CMR_SFCDIF)) THEN
3859 CMR_URB = CMR_SFCDIF(I,J)
3860 CHR_URB = CHR_SFCDIF(I,J)
3861 CMGR_URB = CMGR_SFCDIF(I,J)
3862 CHGR_URB = CHGR_SFCDIF(I,J)
3863 CMC_URB = CMC_SFCDIF(I,J)
3864 CHC_URB = CHC_SFCDIF(I,J)
3868 mh_urb = mh_urb2d(I,J)
3869 stdh_urb = stdh_urb2d(I,J)
3870 lp_urb = lp_urb2d(I,J)
3871 hgt_urb = hgt_urb2d(I,J)
3874 lf_urb(K)=lf_urb2d(I,K,J)
3876 frc_urb = frc_urb2d(I,J)
3877 lb_urb = lb_urb2d(I,J)
3879 if (I.eq.73.and.J.eq.125)THEN
3882 ! Distributed aerodynamics
3883 lf_urb_s = lf_urb2d_s(I, J)
3884 z0_urb = z0_urb2d(I, J)
3885 vegfrac = vegfra(I, J) / 100.
3888 CALL cal_mon_day(julian,julyr,jmonth,jday)
3889 CALL urban(LSOLAR_URB, & ! I
3890 num_roof_layers,num_wall_layers,num_road_layers, & ! C
3892 UTYPE_URB,TA_URB,QA_URB,UA_URB,U1_URB,V1_URB,SSG_URB, & ! I
3893 SSGD_URB,SSGQ_URB,LLG_URB,RAIN_URB,RHOO_URB, & ! I
3894 ZA_URB,DECLIN_URB,COSZ_URB,OMG_URB, & ! I
3895 XLAT_URB,DELT_URB,ZNT_URB, & ! I
3896 CHS_URB, CHS2_URB, & ! I
3897 TR_URB, TB_URB, TG_URB, TC_URB, QC_URB,UC_URB, & ! H
3898 TRL_URB,TBL_URB,TGL_URB, & ! H
3899 XXXR_URB, XXXB_URB, XXXG_URB, XXXC_URB, & ! H
3900 TS_URB,QS_URB,SH_URB,LH_URB,LH_KINEMATIC_URB, & ! O
3901 SW_URB,ALB_URB,LW_URB,G_URB,RN_URB,PSIM_URB,PSIH_URB, & ! O
3903 CMR_URB, CHR_URB, CMC_URB, CHC_URB, &
3904 U10_URB, V10_URB, TH2_URB, Q2_URB, & ! O
3905 UST_URB,mh_urb, stdh_urb, lf_urb, lp_urb, & ! 0
3906 hgt_urb,frc_urb,lb_urb, check,CMCR_URB,TGR_URB, & ! H
3907 TGRL_URB,SMR_URB,CMGR_URB,CHGR_URB,jmonth, & ! H
3908 DRELR_URB,DRELB_URB, & ! H
3909 DRELG_URB,FLXHUMR_URB,FLXHUMB_URB,FLXHUMG_URB, &
3910 lf_urb_s, z0_urb, vegfrac)
3915 print*, 'AFTER CALL URBAN'
3916 print*,'num_roof_layers',num_roof_layers, 'num_wall_layers', &
3918 'DZR',DZR,'DZB',DZB,'DZG',DZG,'UTYPE_URB',UTYPE_URB,'TA_URB', &
3920 'QA_URB',QA_URB,'UA_URB',UA_URB,'U1_URB',U1_URB,'V1_URB', &
3922 'SSG_URB',SSG_URB,'SSGD_URB',SSGD_URB,'SSGQ_URB',SSGQ_URB, &
3923 'LLG_URB',LLG_URB,'RAIN_URB',RAIN_URB,'RHOO_URB',RHOO_URB, &
3924 'ZA_URB',ZA_URB, 'DECLIN_URB',DECLIN_URB,'COSZ_URB',COSZ_URB,&
3925 'OMG_URB',OMG_URB,'XLAT_URB',XLAT_URB,'DELT_URB',DELT_URB, &
3926 'ZNT_URB',ZNT_URB,'TR_URB',TR_URB, 'TB_URB',TB_URB,'TG_URB',&
3927 TG_URB,'TC_URB',TC_URB,'QC_URB',QC_URB,'TRL_URB',TRL_URB, &
3928 'TBL_URB',TBL_URB,'TGL_URB',TGL_URB,'XXXR_URB',XXXR_URB, &
3929 'XXXB_URB',XXXB_URB,'XXXG_URB',XXXG_URB,'XXXC_URB',XXXC_URB,&
3930 'TS_URB',TS_URB,'QS_URB',QS_URB,'SH_URB',SH_URB,'LH_URB', &
3931 LH_URB, 'LH_KINEMATIC_URB',LH_KINEMATIC_URB,'SW_URB',SW_URB,&
3932 'ALB_URB',ALB_URB,'LW_URB',LW_URB,'G_URB',G_URB,'RN_URB', &
3933 RN_URB, 'PSIM_URB',PSIM_URB,'PSIH_URB',PSIH_URB, &
3934 'U10_URB',U10_URB,'V10_URB',V10_URB,'TH2_URB',TH2_URB, &
3935 'Q2_URB',Q2_URB,'CHS_URB',CHS_URB,'CHS2_URB',CHS2_URB
3939 TS_URB2D(I,J) = TS_URB
3941 ALBEDO(I,J) = FRC_URB2D(I,J)*ALB_URB+(1-FRC_URB2D(I,J))*ALBEDOK ![-]
3942 HFX(I,J) = FRC_URB2D(I,J)*SH_URB+(1-FRC_URB2D(I,J))*SHEAT ![W/m/m]
3943 QFX(I,J) = FRC_URB2D(I,J)*LH_KINEMATIC_URB &
3944 + (1-FRC_URB2D(I,J))*ETA_KINEMATIC ![kg/m/m/s]
3945 LH(I,J) = FRC_URB2D(I,J)*LH_URB+(1-FRC_URB2D(I,J))*ETA ![W/m/m]
3946 GRDFLX(I,J) = FRC_URB2D(I,J)*G_URB+(1-FRC_URB2D(I,J))*SSOIL ![W/m/m]
3947 TSK(I,J) = FRC_URB2D(I,J)*TS_URB+(1-FRC_URB2D(I,J))*T1 ![K]
3948 Q1 = FRC_URB2D(I,J)*QS_URB+(1-FRC_URB2D(I,J))*Q1 ![-]
3949 ! Convert QSFC back to mixing ratio
3950 QSFC(I,J)= Q1/(1.0-Q1)
3951 UST(I,J)= FRC_URB2D(I,J)*UST_URB+(1-FRC_URB2D(I,J))*UST(I,J) ![m/s]
3952 IF (.not. slucm_distributed_drag) THEN
3953 ZNT(I,J)= EXP(FRC_URB2D(I,J)*ALOG(ZNT_URB)+(1-FRC_URB2D(I,J))* ALOG(ZNT(I,J))) ! ADD BY DAN
3959 print*, ' FRC_URB2D', FRC_URB2D, &
3960 'ALB_URB',ALB_URB, 'ALBEDOK',ALBEDOK, &
3961 'ALBEDO(I,J)', ALBEDO(I,J), &
3962 'SH_URB',SH_URB,'SHEAT',SHEAT, 'HFX(I,J)',HFX(I,J), &
3963 'LH_KINEMATIC_URB',LH_KINEMATIC_URB,'ETA_KINEMATIC', &
3964 ETA_KINEMATIC, 'QFX(I,J)',QFX(I,J), &
3965 'LH_URB',LH_URB, 'ETA',ETA, 'LH(I,J)',LH(I,J), &
3966 'G_URB',G_URB,'SSOIL',SSOIL,'GRDFLX(I,J)', GRDFLX(I,J),&
3967 'TS_URB',TS_URB,'T1',T1,'TSK(I,J)',TSK(I,J), &
3968 'QS_URB',QS_URB,'Q1',Q1,'QSFC(I,J)',QSFC(I,J)
3972 ! Renew Urban State Varialbes
3974 TR_URB2D(I,J) = TR_URB
3975 TB_URB2D(I,J) = TB_URB
3976 TG_URB2D(I,J) = TG_URB
3977 TC_URB2D(I,J) = TC_URB
3978 QC_URB2D(I,J) = QC_URB
3979 UC_URB2D(I,J) = UC_URB
3981 DO K = 1,num_roof_layers
3982 TRL_URB3D(I,K,J) = TRL_URB(K)
3983 SMR_URB3D(I,K,J) = SMR_URB(K)
3984 TGRL_URB3D(I,K,J)= TGRL_URB(K)
3986 DO K = 1,num_wall_layers
3987 TBL_URB3D(I,K,J) = TBL_URB(K)
3989 DO K = 1,num_road_layers
3990 TGL_URB3D(I,K,J) = TGL_URB(K)
3993 TGR_URB2D(I,J) =TGR_URB
3994 CMCR_URB2D(I,J)=CMCR_URB
3995 FLXHUMR_URB2D(I,J)=FLXHUMR_URB
3996 FLXHUMB_URB2D(I,J)=FLXHUMB_URB
3997 FLXHUMG_URB2D(I,J)=FLXHUMG_URB
3998 DRELR_URB2D(I,J) = DRELR_URB
3999 DRELB_URB2D(I,J) = DRELB_URB
4000 DRELG_URB2D(I,J) = DRELG_URB
4002 XXXR_URB2D(I,J) = XXXR_URB
4003 XXXB_URB2D(I,J) = XXXB_URB
4004 XXXG_URB2D(I,J) = XXXG_URB
4005 XXXC_URB2D(I,J) = XXXC_URB
4007 SH_URB2D(I,J) = SH_URB
4008 LH_URB2D(I,J) = LH_URB
4009 G_URB2D(I,J) = G_URB
4010 RN_URB2D(I,J) = RN_URB
4011 PSIM_URB2D(I,J) = PSIM_URB
4012 PSIH_URB2D(I,J) = PSIH_URB
4013 GZ1OZ0_URB2D(I,J)= GZ1OZ0_URB
4014 U10_URB2D(I,J) = U10_URB
4015 V10_URB2D(I,J) = V10_URB
4016 TH2_URB2D(I,J) = TH2_URB
4017 Q2_URB2D(I,J) = Q2_URB
4018 UST_URB2D(I,J) = UST_URB
4019 AKMS_URB2D(I,J) = KARMAN * UST_URB2D(I,J)/(GZ1OZ0_URB2D(I,J)-PSIM_URB2D(I,J))
4020 IF (PRESENT(CMR_SFCDIF)) THEN
4021 CMR_SFCDIF(I,J) = CMR_URB
4022 CHR_SFCDIF(I,J) = CHR_URB
4023 CMGR_SFCDIF(I,J) = CMGR_URB
4024 CHGR_SFCDIF(I,J) = CHGR_URB
4025 CMC_SFCDIF(I,J) = CMC_URB
4026 CHC_SFCDIF(I,J) = CHC_URB
4029 ! 2D to 3D mosaic danli
4031 TR_URB2D_mosaic(I,mosaic_i,J)=TR_URB2D(I,J)
4032 TB_URB2D_mosaic(I,mosaic_i,J)=TB_URB2D(I,J)
4033 TG_URB2D_mosaic(I,mosaic_i,J)=TG_URB2D(I,J)
4034 TC_URB2D_mosaic(I,mosaic_i,J)=TC_URB2D(I,J)
4035 QC_URB2D_mosaic(I,mosaic_i,J)=QC_URB2D(I,J)
4036 UC_URB2D_mosaic(I,mosaic_i,J)=UC_URB2D(I,J)
4037 TS_URB2D_mosaic(I,mosaic_i,J)=TS_URB2D(I,J)
4038 TS_RUL2D_mosaic(I,mosaic_i,J)=T1
4040 DO K = 1,num_roof_layers
4041 TRL_URB3D_mosaic(I,K+(mosaic_i-1)*num_roof_layers,J)=TRL_URB3D(I,K,J)
4043 DO K = 1,num_wall_layers
4044 TBL_URB3D_mosaic(I,K+(mosaic_i-1)*num_roof_layers,J)=TBL_URB3D(I,K,J)
4046 DO K = 1,num_road_layers
4047 TGL_URB3D_mosaic(I,K+(mosaic_i-1)*num_roof_layers,J)=TGL_URB3D(I,K,J)
4050 SH_URB2D_mosaic(I,mosaic_i,J) = SH_URB2D(I,J)
4051 LH_URB2D_mosaic(I,mosaic_i,J) = LH_URB2D(I,J)
4052 G_URB2D_mosaic(I,mosaic_i,J) = G_URB2D(I,J)
4053 RN_URB2D_mosaic(I,mosaic_i,J) = RN_URB2D(I,J)
4057 ENDIF ! end of UCM CALL if block
4058 !--------------------------------------
4059 ! Urban Part End - urban
4060 !--------------------------------------
4064 SMSTOT(I,J)=SOILM*1000.
4066 SMCREL(I,NS,J)=SMAV(NS)
4069 ! Convert the water unit into mm
4070 SFCRUNOFF(I,J)=SFCRUNOFF(I,J)+RUNOFF1*DT*1000.0
4071 UDRUNOFF(I,J)=UDRUNOFF(I,J)+RUNOFF2*DT*1000.0
4072 ! snow defined when fraction of frozen precip (FFROZP) > 0.5,
4073 IF(FFROZP.GT.0.5)THEN
4074 ACSNOW(I,J)=ACSNOW(I,J)+PRCP*DT
4076 IF(SNOW(I,J).GT.0.)THEN
4077 ACSNOM(I,J)=ACSNOM(I,J)+SNOMLT*1000.
4078 ! accumulated snow-melt energy
4079 SNOPCX(I,J)=SNOPCX(I,J)-SNOMLT/FDTLIW
4082 ENDIF ! endif of land-sea test
4084 !-----------------------------------------------------------------------
4085 ! Done with the Noah-UCM MOSAIC DANLI
4086 !-----------------------------------------------------------------------
4088 TSK_mosaic(i,mosaic_i,j)=TSK(i,j) ! from 2D to 3D
4089 QSFC_mosaic(i,mosaic_i,j)=QSFC(i,j)
4090 CANWAT_mosaic(i,mosaic_i,j)=CANWAT(i,j)
4091 SNOW_mosaic(i,mosaic_i,j)=SNOW(i,j)
4092 SNOWH_mosaic(i,mosaic_i,j)=SNOWH(i,j)
4093 SNOWC_mosaic(i,mosaic_i,j)=SNOWC(i,j)
4095 ALBEDO_mosaic(i,mosaic_i,j)=ALBEDO(i,j)
4096 ALBBCK_mosaic(i,mosaic_i,j)=ALBBCK(i,j)
4097 EMISS_mosaic(i,mosaic_i,j)=EMISS(i,j)
4098 EMBCK_mosaic(i,mosaic_i,j)=EMBCK(i,j)
4099 ZNT_mosaic(i,mosaic_i,j)=ZNT(i,j)
4100 Z0_mosaic(i,mosaic_i,j)=Z0(i,j)
4101 LAI_mosaic(i,mosaic_i,j)=XLAI
4102 RC_mosaic(i,mosaic_i,j)=RC
4104 HFX_mosaic(i,mosaic_i,j)=HFX(i,j)
4105 QFX_mosaic(i,mosaic_i,j)=QFX(i,j)
4106 LH_mosaic(i,mosaic_i,j)=LH(i,j)
4107 GRDFLX_mosaic(i,mosaic_i,j)=GRDFLX(i,j)
4108 SNOTIME_mosaic(i,mosaic_i,j)=SNOTIME(i,j)
4112 TSLB_mosaic(i,NSOIL*(mosaic_i-1)+NS,j)=TSLB(i,NS,j)
4113 SMOIS_mosaic(i,NSOIL*(mosaic_i-1)+NS,j)=SMOIS(i,NS,j)
4114 SH2O_mosaic(i,NSOIL*(mosaic_i-1)+NS,j)=SH2O(i,NS,j)
4119 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
4120 print*, 'I', I, 'J', J, 'MOSAIC_I', MOSAIC_I
4121 print*, 'mosaic_cat_index',mosaic_cat_index(I,mosaic_i,J), 'landusef2',landusef2(i,mosaic_i,j)
4122 print*, 'TSK_mosaic', TSK_mosaic(i,mosaic_i,j), 'HFX_mosaic', HFX_mosaic(i,mosaic_i,j), &
4123 'LH_mosaic',LH_mosaic(i,mosaic_i,j),'GRDFLX_mosaic',GRDFLX_mosaic(i,mosaic_i,j)
4124 print*, 'ZNT_mosaic', ZNT_mosaic(i, mosaic_i,j), 'Z0_mosaic', Z0_mosaic(i,mosaic_i,j)
4125 print*, 'LAI_mosaic', LAI_mosaic(i, mosaic_i,j)
4126 print*, 'FRC_URB2D',FRC_URB2D(I,J)
4127 print*, 'TS_URB',TS_URB2D(I,J),'T1',T1
4128 print*, 'SH_URB2D',SH_URB2D(I,J),'SHEAT',SHEAT
4129 print*, 'LH_URB',LH_URB2D(I,J),'ETA',ETA
4130 print*, 'TS_RUL2D',TS_RUL2D_mosaic(I,mosaic_i,J)
4135 !-----------------------------------------------------------------------
4136 ! Now let's do the grid-averaging
4137 !-----------------------------------------------------------------------
4139 FAREA = landusef2(i,mosaic_i,j)
4141 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
4143 QSFC_mosaic_avg(i,j) = QSFC_mosaic_avg(i,j) + QSFC_mosaic(i,mosaic_i,j)*FAREA
4144 CANWAT_mosaic_avg(i,j) = CANWAT_mosaic_avg(i,j) + CANWAT_mosaic(i,mosaic_i,j)*FAREA
4145 SNOW_mosaic_avg(i,j) = SNOW_mosaic_avg(i,j) + SNOW_mosaic(i,mosaic_i,j)*FAREA
4146 SNOWH_mosaic_avg(i,j) = SNOWH_mosaic_avg(i,j) + SNOWH_mosaic(i,mosaic_i,j)*FAREA
4147 SNOWC_mosaic_avg(i,j) = SNOWC_mosaic_avg(i,j) + SNOWC_mosaic(i,mosaic_i,j)*FAREA
4151 TSLB_mosaic_avg(i,NS,j)=TSLB_mosaic_avg(i,NS,j) + TSLB_mosaic(i,NS*mosaic_i,j)*FAREA
4152 SMOIS_mosaic_avg(i,NS,j)=SMOIS_mosaic_avg(i,NS,j) + SMOIS_mosaic(i,NS*mosaic_i,j)*FAREA
4153 SH2O_mosaic_avg(i,NS,j)=SH2O_mosaic_avg(i,NS,j) + SH2O_mosaic(i,NS*mosaic_i,j)*FAREA
4157 FAREA_mosaic_avg(i,j)=FAREA_mosaic_avg(i,j)+FAREA
4158 HFX_mosaic_avg(i,j) = HFX_mosaic_avg(i,j) + HFX_mosaic(i,mosaic_i,j)*FAREA
4159 QFX_mosaic_avg(i,j) = QFX_mosaic_avg(i,j) + QFX_mosaic(i,mosaic_i,j)*FAREA
4160 LH_mosaic_avg(i,j) = LH_mosaic_avg(i,j) + LH_mosaic(i,mosaic_i,j)*FAREA
4161 GRDFLX_mosaic_avg(i,j)=GRDFLX_mosaic_avg(i,j)+GRDFLX_mosaic(i,mosaic_i,j)*FAREA
4163 ALBEDO_mosaic_avg(i,j)=ALBEDO_mosaic_avg(i,j)+ALBEDO_mosaic(i,mosaic_i,j)*FAREA
4164 ALBBCK_mosaic_avg(i,j)=ALBBCK_mosaic_avg(i,j)+ALBBCK_mosaic(i,mosaic_i,j)*FAREA
4165 EMISS_mosaic_avg(i,j)=EMISS_mosaic_avg(i,j)+EMISS_mosaic(i,mosaic_i,j)*FAREA
4166 EMBCK_mosaic_avg(i,j)=EMBCK_mosaic_avg(i,j)+EMBCK_mosaic(i,mosaic_i,j)*FAREA
4167 ZNT_mosaic_avg(i,j)=ZNT_mosaic_avg(i,j)+ALOG(ZNT_mosaic(i,mosaic_i,j))*FAREA
4168 Z0_mosaic_avg(i,j)=Z0_mosaic_avg(i,j)+ALOG(Z0_mosaic(i,mosaic_i,j))*FAREA
4169 LAI_mosaic_avg(i,j)=LAI_mosaic_avg(i,j)+LAI_mosaic(i,mosaic_i,j)*FAREA
4170 if(RC_mosaic(i,mosaic_i,j) .Gt. 0.0) Then
4171 RC_mosaic_avg(i,j) = RC_mosaic_avg(i,j)+1.0/RC_mosaic(i,mosaic_i,j)*FAREA
4173 RC_mosaic_avg(i,j) = RC_mosaic_avg(i,j) + RC_mosaic(i,mosaic_i,j)*FAREA
4175 ENDDO ! ENDDO FOR mosaic_i = 1, mosaic_cat
4177 !-----------------------------------------------------------------------
4178 ! Now let's send the 3D values to the 2D variables that might be needed in other routines
4179 !-----------------------------------------------------------------------
4181 IVGTYP(I,J)=IVGTYP_dominant(I,J) ! the dominant vege category
4182 ALBEDO(i,j)=ALBEDO_mosaic_avg(i,j)
4183 ALBBCK(i,j)=ALBBCK_mosaic_avg(i,j)
4184 EMISS(i,j)= EMISS_mosaic_avg(i,j)
4185 EMBCK(i,j)= EMBCK_mosaic_avg(i,j)
4186 ZNT(i,j)= EXP(ZNT_mosaic_avg(i,j)/FAREA_mosaic_avg(i,j))
4187 Z0(i,j)= EXP(Z0_mosaic_avg(i,j)/FAREA_mosaic_avg(i,j))
4188 XLAI2(i,j)= LAI_mosaic_avg(i,j)
4189 IF (RC_mosaic_avg(i,j) .Gt. 0.0) THEN
4190 rc2(i,j) = 1.0/(RC_mosaic_avg(i,j))
4192 !RC_mosaic_avg was zero for all tiles (cell over water), thus RC2 set to zero to avoid infinity
4193 rc2(i,j) = RC_mosaic_avg(i,j)
4195 TSK(i,j)=(TSK_mosaic_avg(I,J)/EMISS_mosaic_avg(I,J))**(0.25) ! from 3D to 2D
4196 QSFC(i,j)=QSFC_mosaic_avg(I,J)
4197 CANWAT(i,j) = CANWAT_mosaic_avg(i,j)
4198 SNOW(i,j) = SNOW_mosaic_avg(i,j)
4199 SNOWH(i,j) = SNOWH_mosaic_avg(i,j)
4200 SNOWC(i,j) = SNOWC_mosaic_avg(i,j)
4202 HFX(i,j) = HFX_mosaic_avg(i,j)
4203 QFX(i,j) = QFX_mosaic_avg(i,j)
4204 LH(i,j) = LH_mosaic_avg(i,j)
4205 GRDFLX(i,j)=GRDFLX_mosaic_avg(i,j)
4209 TSLB(i,NS,j)=TSLB_mosaic_avg(i,NS,j)
4210 SMOIS(i,NS,j)=SMOIS_mosaic_avg(i,NS,j)
4211 SH2O(i,NS,j)=SH2O_mosaic_avg(i,NS,j)
4215 ELSE ! This corresponds to IF ((sf_surface_mosaic == 1) .AND. ((XLAND(I,J)-1.5).LT.0.) .AND. (XICE(I,J) < XICE_THRESHOLD) ) THEN
4219 ! pressure in middle of lowest layer
4220 SFCPRS=(P8W3D(I,KTS+1,j)+P8W3D(i,KTS,j))*0.5
4221 ! convert from mixing ratio to specific humidity
4222 Q2K=QV3D(i,1,j)/(1.0+QV3D(i,1,j))
4225 Q2SAT=QGH(I,J)/(1.0+QGH(I,J)) ! Q2SAT is sp humidity
4226 ! add check on myj=.true.
4227 ! IF((Q2K.GE.Q2SAT*TRESH).AND.Q2K.LT.QZ0(I,J))THEN
4228 IF((myj).AND.(Q2K.GE.Q2SAT*TRESH).AND.Q2K.LT.QZ0(I,J))THEN
4237 ZLVL=0.5*DZ8W(i,1,j)
4239 ! TH2=SFCTMP+(0.0097545*ZLVL)
4240 ! calculate SFCTH2 via Exner function vs lapse-rate (above)
4241 APES=(1.E5/PSFC)**CAPA
4242 APELM=(1.E5/SFCPRS)**CAPA
4247 LWDN=GLW(I,J)*EMISSI
4248 ! SOLDN is total incoming solar
4250 ! GSW is net downward solar
4252 ! use mid-day albedo to determine net downward solar (no solar zenith angle correction)
4253 SOLNET=SOLDN*(1.-ALBEDO(I,J))
4257 SHDFAC=VEGFRA(I,J)/100.
4260 SHMIN=SHDMIN(I,J)/100. !NEW
4261 SHMAX=SHDMAX(I,J)/100. !NEW
4262 ! convert snow water equivalent from mm to meter
4263 SNEQV=SNOW(I,J)*0.001
4264 ! snow depth in meters
4268 ! if "SR" present, set frac of frozen precip ("FFROZP") = snow-ratio ("SR", range:0-1)
4269 ! SR from e.g. Ferrier microphysics
4270 ! otherwise define from 1st atmos level temperature
4274 IF (SFCTMP <= 273.15) THEN
4281 IF((XLAND(I,J)-1.5).GE.0.)THEN ! begining of land/sea if block
4283 TSK_RURAL(I,J)=TSK(I,J)
4284 HFX_RURAL(I,J)=HFX(I,J)
4285 QFX_RURAL(I,J)=QFX(I,J)
4286 LH_RURAL(I,J)=LH(I,J)
4287 EMISS_RURAL(I,J)=EMISS(I,J)
4288 GRDFLX_RURAL(I,J)=GRDFLX(I,J)
4290 ! Land or sea-ice case
4292 IF (XICE(I,J) >= XICE_THRESHOLD) THEN
4295 ELSE IF ( VEGTYP == ISICE ) THEN
4299 ! Neither sea ice or land ice.
4302 DQSDT2=Q2SAT*A23M4/(SFCTMP-A4)**2
4304 IF(SNOW(I,J).GT.0.0)THEN
4305 ! snow on surface (use ice saturation properties)
4307 E2SAT=611.2*EXP(6174.*(1./273.15 - 1./SFCTSNO))
4308 Q2SATI=0.622*E2SAT/(SFCPRS-E2SAT)
4309 Q2SATI=Q2SATI/(1.0+Q2SATI) ! spec. hum.
4310 IF (T1 .GT. 273.14) THEN
4311 ! warm ground temps, weight the saturation between ice and water according to SNOWC
4312 Q2SAT=Q2SAT*(1.-SNOWC(I,J)) + Q2SATI*SNOWC(I,J)
4313 DQSDT2=DQSDT2*(1.-SNOWC(I,J)) + Q2SATI*6174./(SFCTSNO**2)*SNOWC(I,J)
4315 ! cold ground temps, use ice saturation only
4317 DQSDT2=Q2SATI*6174./(SFCTSNO**2)
4319 ! for snow cover fraction at 0 C, ground temp will not change, so DQSDT2 effectively zero
4320 IF(T1 .GT. 273. .AND. SNOWC(I,J) .GT. 0.)DQSDT2=DQSDT2*(1.-SNOWC(I,J))
4323 ! Land-ice or land points use the usual deep-soil temperature.
4326 IF(VEGTYP.EQ.25) SHDFAC=0.0000
4327 IF(VEGTYP.EQ.26) SHDFAC=0.0000
4328 IF(VEGTYP.EQ.27) SHDFAC=0.0000
4329 IF(SOILTYP.EQ.14.AND.XICE(I,J).EQ.0.)THEN
4331 IF(IPRINT)PRINT*,' SOIL TYPE FOUND TO BE WATER AT A LAND-POINT'
4332 IF(IPRINT)PRINT*,i,j,'RESET SOIL in surfce.F'
4336 SNOALB1 = SNOALB(I,J)
4337 CMC=CANWAT(I,J)/1000.
4339 !-------------------------------------------
4340 !*** convert snow depth from mm to meter
4343 ! SNOALB=ALBMAX(I,J)*0.01
4345 ! SNOALB=MAXALB(IVGTPK)*0.01
4353 SNOTIME1 = SNOTIME(I,J)
4355 !FEI: temporaray arrays above need to be changed later by using SI
4358 SMC(NS)=SMOIS(I,NS,J)
4359 STC(NS)=TSLB(I,NS,J) !STEMP
4360 SWC(NS)=SH2O(I,NS,J)
4363 if ( (SNEQV.ne.0..AND.SNOWHK.eq.0.).or.(SNOWHK.le.SNEQV) )THEN
4368 !Fei: urban. for urban surface, if calling UCM, redefine the natural surface in cities as
4369 ! the "NATURAL" category in the VEGPARM.TBL
4370 IF(SF_URBAN_PHYSICS == 1.OR. SF_URBAN_PHYSICS==2.OR.SF_URBAN_PHYSICS==3 ) THEN
4372 IF( IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == LCZ_1 .or. IVGTYP(I,J) == LCZ_2 .or. &
4373 IVGTYP(I,J) == LCZ_3 .or. IVGTYP(I,J) == LCZ_4 .or. IVGTYP(I,J) == LCZ_5 .or. &
4374 IVGTYP(I,J) == LCZ_6 .or. IVGTYP(I,J) == LCZ_7 .or. IVGTYP(I,J) == LCZ_8 .or. &
4375 IVGTYP(I,J) == LCZ_9 .or. IVGTYP(I,J) == LCZ_10 .or. IVGTYP(I,J) == LCZ_11 ) THEN
4378 SHDFAC = SHDTBL(NATURAL)
4381 EMISSI = 0.98 !for VEGTYP=5
4382 LWDN = GLW(I,J) * EMISSI
4383 SOLNET = SOLDN * (1.0 - ALBEDOK)
4385 IF ( FRC_URB2D(I,J) < 0.99 ) THEN
4386 if(sf_urban_physics.eq.1)then
4387 T1= ( TSK(I,J) -FRC_URB2D(I,J) * TS_URB2D (I,J) )/ (1-FRC_URB2D(I,J))
4388 elseif((sf_urban_physics.eq.2).OR.(sf_urban_physics.eq.3))then
4389 T1=tsk_rural_bep(i,j)
4396 IF( IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == LCZ_1 .or. IVGTYP(I,J) == LCZ_2 .or. &
4397 IVGTYP(I,J) == LCZ_3 .or. IVGTYP(I,J) == LCZ_4 .or. IVGTYP(I,J) == LCZ_5 .or. &
4398 IVGTYP(I,J) == LCZ_6 .or. IVGTYP(I,J) == LCZ_7 .or. IVGTYP(I,J) == LCZ_8 .or. &
4399 IVGTYP(I,J) == LCZ_9 .or. IVGTYP(I,J) == LCZ_10 .or. IVGTYP(I,J) == LCZ_11 ) THEN
4405 !===Yang, 2014/10/08, hydrological processes for urban vegetation in single layer UCM===
4409 IRIOPTION=IRI_SCHEME
4411 tloc=mod(int(OMG/3.14159*180./15.+12.+0.5 ),24)
4412 if (tloc.lt.0) tloc=tloc+24
4413 if (tloc==0) tloc=24
4414 CALL cal_mon_day(julian,julyr,jmonth,jday)
4415 IF(SF_URBAN_PHYSICS == 1) THEN
4416 IF( IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == LCZ_1 .or. IVGTYP(I,J) == LCZ_2 .or. &
4417 IVGTYP(I,J) == LCZ_3 .or. IVGTYP(I,J) == LCZ_4 .or. IVGTYP(I,J) == LCZ_5 .or. &
4418 IVGTYP(I,J) == LCZ_6 .or. IVGTYP(I,J) == LCZ_7 .or. IVGTYP(I,J) == LCZ_8 .or. &
4419 IVGTYP(I,J) == LCZ_9 .or. IVGTYP(I,J) == LCZ_10 .or. IVGTYP(I,J) == LCZ_11 ) THEN
4421 AOASIS = oasis ! urban oasis effect
4422 IF (IRIOPTION ==1) THEN
4423 IF (tloc==21 .or. tloc==22) THEN !irrigation on vegetaion in urban area, MAY-SEP, 9-10pm
4424 IF (jmonth==5 .or. jmonth==6 .or. jmonth==7 .or. jmonth==8 .or. jmonth==9) THEN
4425 IF (SMC(USOIL) .LT. SMCREF) SMC(USOIL)= REFSMC(ISLTYP(I,J))
4426 IF (SMC(DSOIL) .LT. SMCREF) SMC(DSOIL)= REFSMC(ISLTYP(I,J))
4433 IF(SF_URBAN_PHYSICS == 2 .or. SF_URBAN_PHYSICS == 3) THEN
4434 IF(AOASIS > 1.0) THEN
4435 CALL wrf_error_fatal('Urban oasis option is for SF_URBAN_PHYSICS == 1 only')
4437 IF(IRIOPTION == 1) THEN
4438 CALL wrf_error_fatal('Urban irrigation option is for SF_URBAN_PHYSICS == 1 only')
4445 print*, 'BEFORE SFLX, in Noahlsm_driver'
4446 print*, 'ICE', ICE, 'DT',DT, 'ZLVL',ZLVL, 'NSOIL', NSOIL, &
4447 'SLDPTH', SLDPTH, 'LOCAL',LOCAL, 'LUTYPE',&
4448 LUTYPE, 'SLTYPE',SLTYPE, 'LWDN',LWDN, 'SOLDN',SOLDN, &
4449 'SFCPRS',SFCPRS, 'PRCP',PRCP,'SFCTMP',SFCTMP,'Q2K',Q2K, &
4450 'TH2',TH2,'Q2SAT',Q2SAT,'DQSDT2',DQSDT2,'VEGTYP', VEGTYP,&
4451 'SOILTYP',SOILTYP, 'SLOPETYP',SLOPETYP, 'SHDFAC',SHDFAC,&
4452 'SHMIN',SHMIN, 'ALBBRD',ALBBRD,'SNOALB1',SNOALB1,'TBOT',&
4453 TBOT, 'Z0BRD',Z0BRD, 'Z0K',Z0K, 'CMC',CMC, 'T1',T1,'STC',&
4454 STC, 'SMC',SMC, 'SWC',SWC,'SNOWHK',SNOWHK,'SNEQV',SNEQV,&
4455 'ALBEDOK',ALBEDOK,'CHK',CHK,'ETA',ETA,'SHEAT',SHEAT, &
4456 'ETA_KINEMATIC',ETA_KINEMATIC, 'FDOWN',FDOWN,'EC',EC, &
4457 'EDIR',EDIR,'ET',ET,'ETT',ETT,'ESNOW',ESNOW,'DRIP',DRIP,&
4458 'DEW',DEW,'BETA',BETA,'ETP',ETP,'SSOIL',SSOIL,'FLX1',FLX1,&
4459 'FLX2',FLX2,'FLX3',FLX3,'SNOMLT',SNOMLT,'SNCOVR',SNCOVR,&
4460 'RUNOFF1',RUNOFF1,'RUNOFF2',RUNOFF2,'RUNOFF3',RUNOFF3, &
4461 'RC',RC, 'PC',PC,'RSMIN',RSMIN,'XLAI',XLAI,'RCS',RCS, &
4462 'RCT',RCT,'RCQ',RCQ,'RCSOIL',RCSOIL,'SOILW',SOILW, &
4463 'SOILM',SOILM,'Q1',Q1,'SMCWLT',SMCWLT,'SMCDRY',SMCDRY,&
4464 'SMCREF',SMCREF,'SMCMAX',SMCMAX,'NROOT',NROOT
4472 IF ( ICE == 1 ) THEN
4483 ELSEIF (ICE == 0) THEN
4487 CALL SFLX (I,J,FFROZP, ISURBAN, DT,ZLVL,NSOIL,SLDPTH, & !C
4489 LUTYPE, SLTYPE, & !CL
4490 LWDN,SOLDN,SOLNET,SFCPRS,PRCP,SFCTMP,Q2K,DUMMY, & !F
4491 DUMMY,DUMMY, DUMMY, & !F PRCPRAIN not used
4492 TH2,Q2SAT,DQSDT2, & !I
4493 VEGTYP,SOILTYP,SLOPETYP,SHDFAC,SHMIN,SHMAX, & !I
4494 ALBBRD, SNOALB1,TBOT, Z0BRD, Z0K, EMISSI, EMBRD, & !S
4495 CMC,T1,STC,SMC,SWC,SNOWHK,SNEQV,ALBEDOK,CHK,dummy,& !H
4496 ETA,SHEAT, ETA_KINEMATIC,FDOWN, & !O
4497 EC,EDIR,ET,ETT,ESNOW,DRIP,DEW, & !O
4498 BETA,ETP,SSOIL, & !O
4499 FLX1,FLX2,FLX3, & !O
4500 FLX4,FVB,FBUR,FGSN,UA_PHYS, & !UA
4502 RUNOFF1,RUNOFF2,RUNOFF3, & !O
4503 RC,PC,RSMIN,XLAI,RCS,RCT,RCQ,RCSOIL, & !O
4504 SOILW,SOILM,Q1,SMAV, & !D
4505 RDLAI2D,USEMONALB, &
4508 SMCWLT,SMCDRY,SMCREF,SMCMAX,NROOT, &
4510 sfcheadrt_hydro, & !I
4511 INFXSRT_hydro,ETPND1_hydro & !O
4512 ,OPT_THCND,AOASIS & !O
4513 ,XSDA_QFX, HFX_PHY, QFX_PHY, XQNORM, fasdas, HCPCT_FASDAS & ! fasdas vars
4514 ,IRRIGATION_CHANNEL )
4517 soldrain(i,j) = RUNOFF2*DT*1000.0
4518 sfcheadrt(i,j) = sfcheadrt_hydro
4519 infxsrt(i,j) = INFXSRT_hydro
4520 etpnd1 = etpnd1_hydro
4522 ELSEIF (ICE == -1) THEN
4525 ! Set values that the LSM is expected to update,
4526 ! but don't get updated for glacial points.
4528 SOILM = 0.0 !BSINGH(PNNL)- SOILM is undefined for this case, it is used for diagnostics so setting it to zero
4529 XLAI = 0.01 ! KWM Should this be Zero over land ice? Does this value matter?
4537 CALL SFLX_GLACIAL(I,J,ISICE,FFROZP,DT,ZLVL,NSOIL,SLDPTH, & !C
4538 & LWDN,SOLNET,SFCPRS,PRCP,SFCTMP,Q2K, & !F
4539 & TH2,Q2SAT,DQSDT2, & !I
4540 & ALBBRD, SNOALB1,TBOT, Z0BRD, Z0K, EMISSI, EMBRD, & !S
4541 & T1,STC(1:NSOIL),SNOWHK,SNEQV,ALBEDOK,CHK, & !H
4542 & ETA,SHEAT,ETA_KINEMATIC,FDOWN, & !O
4545 & FLX1,FLX2,FLX3, & !O
4546 & SNOMLT,SNCOVR, & !O
4559 print*, 'AFTER SFLX, in Noahlsm_driver'
4560 print*, 'ICE', ICE, 'DT',DT, 'ZLVL',ZLVL, 'NSOIL', NSOIL, &
4561 'SLDPTH', SLDPTH, 'LOCAL',LOCAL, 'LUTYPE',&
4562 LUTYPE, 'SLTYPE',SLTYPE, 'LWDN',LWDN, 'SOLDN',SOLDN, &
4563 'SFCPRS',SFCPRS, 'PRCP',PRCP,'SFCTMP',SFCTMP,'Q2K',Q2K, &
4564 'TH2',TH2,'Q2SAT',Q2SAT,'DQSDT2',DQSDT2,'VEGTYP', VEGTYP,&
4565 'SOILTYP',SOILTYP, 'SLOPETYP',SLOPETYP, 'SHDFAC',SHDFAC,&
4566 'SHDMIN',SHMIN, 'ALBBRD',ALBBRD,'SNOALB',SNOALB1,'TBOT',&
4567 TBOT, 'Z0BRD',Z0BRD, 'Z0K',Z0K, 'CMC',CMC, 'T1',T1,'STC',&
4568 STC, 'SMC',SMC, 'SWc',SWC,'SNOWHK',SNOWHK,'SNEQV',SNEQV,&
4569 'ALBEDOK',ALBEDOK,'CHK',CHK,'ETA',ETA,'SHEAT',SHEAT, &
4570 'ETA_KINEMATIC',ETA_KINEMATIC, 'FDOWN',FDOWN,'EC',EC, &
4571 'EDIR',EDIR,'ET',ET,'ETT',ETT,'ESNOW',ESNOW,'DRIP',DRIP,&
4572 'DEW',DEW,'BETA',BETA,'ETP',ETP,'SSOIL',SSOIL,'FLX1',FLX1,&
4573 'FLX2',FLX2,'FLX3',FLX3,'SNOMLT',SNOMLT,'SNCOVR',SNCOVR,&
4574 'RUNOFF1',RUNOFF1,'RUNOFF2',RUNOFF2,'RUNOFF3',RUNOFF3, &
4575 'RC',RC, 'PC',PC,'RSMIN',RSMIN,'XLAI',XLAI,'RCS',RCS, &
4576 'RCT',RCT,'RCQ',RCQ,'RCSOIL',RCSOIL,'SOILW',SOILW, &
4577 'SOILM',SOILM,'Q1',Q1,'SMCWLT',SMCWLT,'SMCDRY',SMCDRY,&
4578 'SMCREF',SMCREF,'SMCMAX',SMCMAX,'NROOT',NROOT
4582 !*** UPDATE STATE VARIABLES
4583 CANWAT(I,J)=CMC*1000.
4584 SNOW(I,J)=SNEQV*1000.
4585 ! SNOWH(I,J)=SNOWHK*1000.
4586 SNOWH(I,J)=SNOWHK ! SNOWHK in meters
4588 ALB_RURAL(I,J)=ALBEDOK
4592 EMISS_RURAL(I,J) = EMISSI
4593 ! Noah: activate time-varying roughness length (V3.3 Feb 2011)
4598 HFX_RURAL(I,J)=SHEAT
4599 ! MEk Jul07 add potential evap accum
4600 POTEVP(I,J)=POTEVP(I,J)+ETP*FDTW
4601 QFX(I,J)=ETA_KINEMATIC
4602 QFX_RURAL(I,J)=ETA_KINEMATIC
4606 ! QFX(I,J) = QFX(I,J) + ETPND1
4607 ! ETA = ETA + ETPND1/2.501E6*dt
4608 !end added by Wei Yu
4614 GRDFLX_RURAL(I,J)=SSOIL
4617 SNOTIME(I,J) = SNOTIME1
4618 ! prevent diagnostic ground q (q1) from being greater than qsat(tsk)
4619 ! as happens over snow cover where the cqs2 value also becomes irrelevant
4620 ! by setting cqs2=chs in this situation the 2m q should become just qv(k=1)
4621 IF (Q1 .GT. QSFC(I,J)) THEN
4622 CQS2(I,J) = CHS(I,J)
4625 ! Convert QSFC back to mixing ratio
4626 QSFC(I,J)= Q1/(1.0-Q1)
4628 ! QSFC_RURAL(I,J)= Q1/(1.0-Q1)
4629 ! Calculate momentum flux from rural surface for use with multi-layer UCM (Martilli et al. 2002)
4632 SMOIS(I,NS,J)=SMC(NS)
4633 TSLB(I,NS,J)=STC(NS) ! STEMP
4634 SH2O(I,NS,J)=SWC(NS)
4643 ! Residual of surface energy balance equation terms
4647 noahres(i,j) = ( solnet + lwdn ) - sheat + ssoil - eta &
4648 - ( emissi * STBOLT * (t1**4) ) - flx1 - flx2 - flx3 - flx4
4651 noahres(i,j) = ( solnet + lwdn ) - sheat + ssoil - eta &
4652 - ( emissi * STBOLT * (t1**4) ) - flx1 - flx2 - flx3
4655 IF (SF_URBAN_PHYSICS == 1 ) THEN ! Beginning of UCM CALL if block
4656 !--------------------------------------
4657 ! URBAN CANOPY MODEL START - urban
4658 !--------------------------------------
4659 ! Input variables lsm --> urban
4660 IF( IVGTYP(I,J) == ISURBAN .or. IVGTYP(I,J) == LCZ_1 .or. IVGTYP(I,J) == LCZ_2 .or. &
4661 IVGTYP(I,J) == LCZ_3 .or. IVGTYP(I,J) == LCZ_4 .or. IVGTYP(I,J) == LCZ_5 .or. &
4662 IVGTYP(I,J) == LCZ_6 .or. IVGTYP(I,J) == LCZ_7 .or. IVGTYP(I,J) == LCZ_8 .or. &
4663 IVGTYP(I,J) == LCZ_9 .or. IVGTYP(I,J) == LCZ_10 .or. IVGTYP(I,J) == LCZ_11 ) THEN
4668 UTYPE_URB = UTYPE_URB2D(I,J) !urban type (low, high or industrial)
4670 TA_URB = SFCTMP ! [K]
4671 QA_URB = Q2K ! [kg/kg]
4672 UA_URB = SQRT(U_PHY(I,1,J)**2.+V_PHY(I,1,J)**2.)
4673 U1_URB = U_PHY(I,1,J)
4674 V1_URB = V_PHY(I,1,J)
4675 IF(UA_URB < 1.) UA_URB=1. ! [m/s]
4676 SSG_URB = SOLDN ! [W/m/m]
4677 SSGD_URB = 0.8*SOLDN ! [W/m/m]
4678 SSGQ_URB = SSG_URB-SSGD_URB ! [W/m/m]
4679 LLG_URB = GLW(I,J) ! [W/m/m]
4680 RAIN_URB = RAINBL(I,J) / DT * 3600.0 ! [mm/hr]
4681 RHOO_URB = SFCPRS / (287.04 * SFCTMP * (1.0+ 0.61 * Q2K)) ![kg/m/m/m]
4683 DELT_URB = DT ! [sec]
4684 XLAT_URB = XLAT_URB2D(I,J) ! [deg]
4685 COSZ_URB = COSZ_URB2D(I,J) !
4686 OMG_URB = OMG_URB2D(I,J) !
4689 LSOLAR_URB = .FALSE.
4691 TR_URB = TR_URB2D(I,J)
4692 TB_URB = TB_URB2D(I,J)
4693 TG_URB = TG_URB2D(I,J)
4694 TC_URB = TC_URB2D(I,J)
4695 QC_URB = QC_URB2D(I,J)
4696 UC_URB = UC_URB2D(I,J)
4698 DO K = 1,num_roof_layers
4699 TRL_URB(K) = TRL_URB3D(I,K,J)
4700 SMR_URB(K) = SMR_URB3D(I,K,J)
4701 TGRL_URB(K)= TGRL_URB3D(I,K,J)
4703 DO K = 1,num_wall_layers
4704 TBL_URB(K) = TBL_URB3D(I,K,J)
4706 DO K = 1,num_road_layers
4707 TGL_URB(K) = TGL_URB3D(I,K,J)
4710 TGR_URB = TGR_URB2D(I,J)
4711 CMCR_URB = CMCR_URB2D(I,J)
4712 FLXHUMR_URB = FLXHUMR_URB2D(I,J)
4713 FLXHUMB_URB = FLXHUMB_URB2D(I,J)
4714 FLXHUMG_URB = FLXHUMG_URB2D(I,J)
4715 DRELR_URB = DRELR_URB2D(I,J)
4716 DRELB_URB = DRELB_URB2D(I,J)
4717 DRELG_URB = DRELG_URB2D(I,J)
4719 XXXR_URB = XXXR_URB2D(I,J)
4720 XXXB_URB = XXXB_URB2D(I,J)
4721 XXXG_URB = XXXG_URB2D(I,J)
4722 XXXC_URB = XXXC_URB2D(I,J)
4724 ! Limits to avoid dividing by small number
4725 if (CHS(I,J) < 1.0E-02) then
4728 if (CHS2(I,J) < 1.0E-02) then
4731 if (CQS2(I,J) < 1.0E-02) then
4736 CHS2_URB = CHS2(I,J)
4737 IF (PRESENT(CMR_SFCDIF)) THEN
4738 CMR_URB = CMR_SFCDIF(I,J)
4739 CHR_URB = CHR_SFCDIF(I,J)
4740 CMGR_URB = CMGR_SFCDIF(I,J)
4741 CHGR_URB = CHGR_SFCDIF(I,J)
4742 CMC_URB = CMC_SFCDIF(I,J)
4743 CHC_URB = CHC_SFCDIF(I,J)
4747 mh_urb = mh_urb2d(I,J)
4748 stdh_urb = stdh_urb2d(I,J)
4749 lp_urb = lp_urb2d(I,J)
4750 hgt_urb = hgt_urb2d(I,J)
4753 lf_urb(K)=lf_urb2d(I,K,J)
4755 frc_urb = frc_urb2d(I,J)
4756 lb_urb = lb_urb2d(I,J)
4758 if (I.eq.73.and.J.eq.125)THEN
4761 ! Distributed aerodynamics
4762 lf_urb_s = lf_urb2d_s(I, J)
4763 z0_urb = z0_urb2d(I, J)
4764 vegfrac = vegfra(I, J) / 100.0
4768 CALL cal_mon_day(julian,julyr,jmonth,jday)
4769 CALL urban(LSOLAR_URB, & ! I
4770 num_roof_layers,num_wall_layers,num_road_layers, & ! C
4772 UTYPE_URB,TA_URB,QA_URB,UA_URB,U1_URB,V1_URB,SSG_URB, & ! I
4773 SSGD_URB,SSGQ_URB,LLG_URB,RAIN_URB,RHOO_URB, & ! I
4774 ZA_URB,DECLIN_URB,COSZ_URB,OMG_URB, & ! I
4775 XLAT_URB,DELT_URB,ZNT_URB, & ! I
4776 CHS_URB, CHS2_URB, & ! I
4777 TR_URB, TB_URB, TG_URB, TC_URB, QC_URB,UC_URB, & ! H
4778 TRL_URB,TBL_URB,TGL_URB, & ! H
4779 XXXR_URB, XXXB_URB, XXXG_URB, XXXC_URB, & ! H
4780 TS_URB,QS_URB,SH_URB,LH_URB,LH_KINEMATIC_URB, & ! O
4781 SW_URB,ALB_URB,LW_URB,G_URB,RN_URB,PSIM_URB,PSIH_URB, & ! O
4783 CMR_URB, CHR_URB, CMC_URB, CHC_URB, &
4784 U10_URB, V10_URB, TH2_URB, Q2_URB, & ! O
4785 UST_URB,mh_urb, stdh_urb, lf_urb, lp_urb, & ! 0
4786 hgt_urb,frc_urb,lb_urb, check,CMCR_URB,TGR_URB, & ! H
4787 TGRL_URB,SMR_URB,CMGR_URB,CHGR_URB,jmonth, & ! H
4788 DRELR_URB,DRELB_URB, & ! H
4789 DRELG_URB,FLXHUMR_URB,FLXHUMB_URB,FLXHUMG_URB, &
4790 lf_urb_s, z0_urb, vegfrac)
4795 print*, 'AFTER CALL URBAN'
4796 print*,'num_roof_layers',num_roof_layers, 'num_wall_layers', &
4798 'DZR',DZR,'DZB',DZB,'DZG',DZG,'UTYPE_URB',UTYPE_URB,'TA_URB', &
4800 'QA_URB',QA_URB,'UA_URB',UA_URB,'U1_URB',U1_URB,'V1_URB', &
4802 'SSG_URB',SSG_URB,'SSGD_URB',SSGD_URB,'SSGQ_URB',SSGQ_URB, &
4803 'LLG_URB',LLG_URB,'RAIN_URB',RAIN_URB,'RHOO_URB',RHOO_URB, &
4804 'ZA_URB',ZA_URB, 'DECLIN_URB',DECLIN_URB,'COSZ_URB',COSZ_URB,&
4805 'OMG_URB',OMG_URB,'XLAT_URB',XLAT_URB,'DELT_URB',DELT_URB, &
4806 'ZNT_URB',ZNT_URB,'TR_URB',TR_URB, 'TB_URB',TB_URB,'TG_URB',&
4807 TG_URB,'TC_URB',TC_URB,'QC_URB',QC_URB,'TRL_URB',TRL_URB, &
4808 'TBL_URB',TBL_URB,'TGL_URB',TGL_URB,'XXXR_URB',XXXR_URB, &
4809 'XXXB_URB',XXXB_URB,'XXXG_URB',XXXG_URB,'XXXC_URB',XXXC_URB,&
4810 'TS_URB',TS_URB,'QS_URB',QS_URB,'SH_URB',SH_URB,'LH_URB', &
4811 LH_URB, 'LH_KINEMATIC_URB',LH_KINEMATIC_URB,'SW_URB',SW_URB,&
4812 'ALB_URB',ALB_URB,'LW_URB',LW_URB,'G_URB',G_URB,'RN_URB', &
4813 RN_URB, 'PSIM_URB',PSIM_URB,'PSIH_URB',PSIH_URB, &
4814 'U10_URB',U10_URB,'V10_URB',V10_URB,'TH2_URB',TH2_URB, &
4815 'Q2_URB',Q2_URB,'CHS_URB',CHS_URB,'CHS2_URB',CHS2_URB
4819 TS_URB2D(I,J) = TS_URB
4821 ALBEDO(I,J) = FRC_URB2D(I,J)*ALB_URB+(1-FRC_URB2D(I,J))*ALBEDOK ![-]
4822 HFX(I,J) = FRC_URB2D(I,J)*SH_URB+(1-FRC_URB2D(I,J))*SHEAT ![W/m/m]
4823 QFX(I,J) = FRC_URB2D(I,J)*LH_KINEMATIC_URB &
4824 + (1-FRC_URB2D(I,J))*ETA_KINEMATIC ![kg/m/m/s]
4825 LH(I,J) = FRC_URB2D(I,J)*LH_URB+(1-FRC_URB2D(I,J))*ETA ![W/m/m]
4826 GRDFLX(I,J) = FRC_URB2D(I,J)*G_URB+(1-FRC_URB2D(I,J))*SSOIL ![W/m/m]
4827 TSK(I,J) = FRC_URB2D(I,J)*TS_URB+(1-FRC_URB2D(I,J))*T1 ![K]
4828 Q1 = FRC_URB2D(I,J)*QS_URB+(1-FRC_URB2D(I,J))*Q1 ![-]
4829 ! Convert QSFC back to mixing ratio
4830 QSFC(I,J)= Q1/(1.0-Q1)
4831 UST(I,J)= FRC_URB2D(I,J)*UST_URB+(1-FRC_URB2D(I,J))*UST(I,J) ![m/s]
4836 print*, ' FRC_URB2D', FRC_URB2D, &
4837 'ALB_URB',ALB_URB, 'ALBEDOK',ALBEDOK, &
4838 'ALBEDO(I,J)', ALBEDO(I,J), &
4839 'SH_URB',SH_URB,'SHEAT',SHEAT, 'HFX(I,J)',HFX(I,J), &
4840 'LH_KINEMATIC_URB',LH_KINEMATIC_URB,'ETA_KINEMATIC', &
4841 ETA_KINEMATIC, 'QFX(I,J)',QFX(I,J), &
4842 'LH_URB',LH_URB, 'ETA',ETA, 'LH(I,J)',LH(I,J), &
4843 'G_URB',G_URB,'SSOIL',SSOIL,'GRDFLX(I,J)', GRDFLX(I,J),&
4844 'TS_URB',TS_URB,'T1',T1,'TSK(I,J)',TSK(I,J), &
4845 'QS_URB',QS_URB,'Q1',Q1,'QSFC(I,J)',QSFC(I,J)
4849 ! Renew Urban State Varialbes
4851 TR_URB2D(I,J) = TR_URB
4852 TB_URB2D(I,J) = TB_URB
4853 TG_URB2D(I,J) = TG_URB
4854 TC_URB2D(I,J) = TC_URB
4855 QC_URB2D(I,J) = QC_URB
4856 UC_URB2D(I,J) = UC_URB
4858 DO K = 1,num_roof_layers
4859 TRL_URB3D(I,K,J) = TRL_URB(K)
4860 SMR_URB3D(I,K,J) = SMR_URB(K)
4861 TGRL_URB3D(I,K,J)= TGRL_URB(K)
4863 DO K = 1,num_wall_layers
4864 TBL_URB3D(I,K,J) = TBL_URB(K)
4866 DO K = 1,num_road_layers
4867 TGL_URB3D(I,K,J) = TGL_URB(K)
4870 TGR_URB2D(I,J) =TGR_URB
4871 CMCR_URB2D(I,J)=CMCR_URB
4872 FLXHUMR_URB2D(I,J)=FLXHUMR_URB
4873 FLXHUMB_URB2D(I,J)=FLXHUMB_URB
4874 FLXHUMG_URB2D(I,J)=FLXHUMG_URB
4875 DRELR_URB2D(I,J) = DRELR_URB
4876 DRELB_URB2D(I,J) = DRELB_URB
4877 DRELG_URB2D(I,J) = DRELG_URB
4879 XXXR_URB2D(I,J) = XXXR_URB
4880 XXXB_URB2D(I,J) = XXXB_URB
4881 XXXG_URB2D(I,J) = XXXG_URB
4882 XXXC_URB2D(I,J) = XXXC_URB
4884 SH_URB2D(I,J) = SH_URB
4885 LH_URB2D(I,J) = LH_URB
4886 G_URB2D(I,J) = G_URB
4887 RN_URB2D(I,J) = RN_URB
4888 PSIM_URB2D(I,J) = PSIM_URB
4889 PSIH_URB2D(I,J) = PSIH_URB
4890 GZ1OZ0_URB2D(I,J)= GZ1OZ0_URB
4891 U10_URB2D(I,J) = U10_URB
4892 V10_URB2D(I,J) = V10_URB
4893 TH2_URB2D(I,J) = TH2_URB
4894 Q2_URB2D(I,J) = Q2_URB
4895 UST_URB2D(I,J) = UST_URB
4896 AKMS_URB2D(I,J) = KARMAN * UST_URB2D(I,J)/(GZ1OZ0_URB2D(I,J)-PSIM_URB2D(I,J))
4897 IF (PRESENT(CMR_SFCDIF)) THEN
4898 CMR_SFCDIF(I,J) = CMR_URB
4899 CHR_SFCDIF(I,J) = CHR_URB
4900 CMGR_SFCDIF(I,J) = CMGR_URB
4901 CHGR_SFCDIF(I,J) = CHGR_URB
4902 CMC_SFCDIF(I,J) = CMC_URB
4903 CHC_SFCDIF(I,J) = CHC_URB
4907 ENDIF ! end of UCM CALL if block
4908 !--------------------------------------
4909 ! Urban Part End - urban
4910 !--------------------------------------
4914 SMSTOT(I,J)=SOILM*1000.
4916 SMCREL(I,NS,J)=SMAV(NS)
4919 ! Convert the water unit into mm
4920 SFCRUNOFF(I,J)=SFCRUNOFF(I,J)+RUNOFF1*DT*1000.0
4921 UDRUNOFF(I,J)=UDRUNOFF(I,J)+RUNOFF2*DT*1000.0
4922 ! snow defined when fraction of frozen precip (FFROZP) > 0.5,
4923 IF(FFROZP.GT.0.5)THEN
4924 ACSNOW(I,J)=ACSNOW(I,J)+PRCP*DT
4926 IF(SNOW(I,J).GT.0.)THEN
4927 ACSNOM(I,J)=ACSNOM(I,J)+SNOMLT*1000.
4928 ! accumulated snow-melt energy
4929 SNOPCX(I,J)=SNOPCX(I,J)-SNOMLT/FDTLIW
4932 ENDIF ! endif of land-sea test
4934 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
4936 ENDDO ILOOP ! of I loop
4937 ENDDO JLOOP ! of J loop
4939 !------------------------------------------------------
4940 END SUBROUTINE lsm_mosaic
4941 !------------------------------------------------------
4942 !===========================================================================
4944 ! subroutine lsm_mosaic_init: initialization of mosaic state variables
4946 !===========================================================================
4948 SUBROUTINE lsm_mosaic_init(IVGTYP,ISWATER,ISURBAN,ISICE, XLAND, XICE,fractional_seaice, &
4949 TSK,TSLB,SMOIS,SH2O,SNOW,SNOWC,SNOWH,CANWAT, &
4950 ids,ide, jds,jde, kds,kde, &
4951 ims,ime, jms,jme, kms,kme, &
4952 its,ite, jts,jte, kts,kte, restart, &
4953 landusef,landusef2,NLCAT,num_soil_layers &
4954 ,sf_surface_mosaic, mosaic_cat &
4956 ,TSK_mosaic,TSLB_mosaic &
4957 ,SMOIS_mosaic,SH2O_mosaic &
4958 ,CANWAT_mosaic,SNOW_mosaic &
4959 ,SNOWH_mosaic,SNOWC_mosaic &
4960 ,ALBEDO,ALBBCK, EMISS, EMBCK,Z0 & !danli
4961 ,ALBEDO_mosaic,ALBBCK_mosaic, EMISS_mosaic & !danli
4962 ,EMBCK_mosaic, ZNT_mosaic, Z0_mosaic & !danli
4963 ,TR_URB2D_mosaic,TB_URB2D_mosaic & !danli mosaic
4964 ,TG_URB2D_mosaic,TC_URB2D_mosaic & !danli mosaic
4965 ,QC_URB2D_mosaic & !danli mosaic
4966 ,TRL_URB3D_mosaic,TBL_URB3D_mosaic & !danli mosaic
4967 ,TGL_URB3D_mosaic & !danli mosaic
4968 ,SH_URB2D_mosaic,LH_URB2D_mosaic & !danli mosaic
4969 ,G_URB2D_mosaic,RN_URB2D_mosaic & !danli mosaic
4970 ,TS_URB2D_mosaic & !danli mosaic
4971 ,TS_RUL2D_mosaic & !danli mosaic
4974 INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, &
4975 ims,ime, jms,jme, kms,kme, &
4976 its,ite, jts,jte, kts,kte
4978 INTEGER, INTENT(IN) :: NLCAT, num_soil_layers, ISWATER,ISURBAN, ISICE, fractional_seaice
4980 LOGICAL , INTENT(IN) :: restart
4982 ! REAL, DIMENSION( num_soil_layers), INTENT(INOUT) :: ZS, DZS
4984 REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , &
4985 INTENT(IN) :: SMOIS, & !Total soil moisture
4986 SH2O, & !liquid soil moisture
4989 REAL, DIMENSION( ims:ime, jms:jme ) , &
4990 INTENT(IN) :: SNOW, &
4996 INTEGER, INTENT(IN) :: sf_surface_mosaic
4997 INTEGER, INTENT(IN) :: mosaic_cat
4998 INTEGER, DIMENSION( ims:ime, jms:jme ),INTENT(IN) :: IVGTYP
4999 REAL, DIMENSION( ims:ime, NLCAT, jms:jme ) , INTENT(IN):: LANDUSEF
5000 REAL, DIMENSION( ims:ime, NLCAT, jms:jme ) , INTENT(INOUT):: LANDUSEF2
5002 INTEGER, DIMENSION( ims:ime, NLCAT, jms:jme ), INTENT(INOUT) :: mosaic_cat_index
5004 REAL, DIMENSION( ims:ime, 1:mosaic_cat, jms:jme ) , OPTIONAL, INTENT(INOUT):: &
5005 TSK_mosaic, CANWAT_mosaic, SNOW_mosaic,SNOWH_mosaic, SNOWC_mosaic
5006 REAL, DIMENSION( ims:ime, 1:num_soil_layers*mosaic_cat, jms:jme ), OPTIONAL, INTENT(INOUT):: &
5007 TSLB_mosaic,SMOIS_mosaic,SH2O_mosaic
5009 REAL, DIMENSION( ims:ime, jms:jme ) , INTENT(IN):: ALBEDO, ALBBCK, EMISS, EMBCK, Z0
5010 REAL, DIMENSION( ims:ime, 1:mosaic_cat, jms:jme ) , OPTIONAL, INTENT(INOUT):: &
5011 ALBEDO_mosaic,ALBBCK_mosaic, EMISS_mosaic, EMBCK_mosaic, ZNT_mosaic, Z0_mosaic
5013 REAL, DIMENSION( ims:ime, 1:mosaic_cat, jms:jme ) , OPTIONAL, INTENT(INOUT):: &
5014 TR_URB2D_mosaic, TB_URB2D_mosaic, TG_URB2D_mosaic, TC_URB2D_mosaic,QC_URB2D_mosaic, &
5015 SH_URB2D_mosaic,LH_URB2D_mosaic,G_URB2D_mosaic,RN_URB2D_mosaic,TS_URB2D_mosaic, TS_RUL2D_mosaic
5017 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_soil_layers*mosaic_cat, jms:jme ), INTENT(INOUT) :: TRL_URB3D_mosaic
5018 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_soil_layers*mosaic_cat, jms:jme ), INTENT(INOUT) :: TBL_URB3D_mosaic
5019 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_soil_layers*mosaic_cat, jms:jme ), INTENT(INOUT) :: TGL_URB3D_mosaic
5021 INTEGER :: ij,i,j,mosaic_i,LastSwap,NumPairs,soil_k, Temp2,Temp5,Temp7, ICE,temp_index
5022 REAL :: Temp, Temp3,Temp4,Temp6,xice_threshold
5024 CHARACTER(len=256) :: message_text
5028 if ( fractional_seaice == 0 ) then
5029 xice_threshold = 0.5
5030 else if ( fractional_seaice == 1 ) then
5031 xice_threshold = 0.02
5034 IF(.not.restart)THEN
5035 !===========================================================================
5037 !===========================================================================
5046 IF ((xland(i,j).LT. 1.5 ) .AND. (IVGTYP(i,j) .EQ. ISWATER)) THEN
5047 PRINT*, 'BEFORE MOSAIC_INIT'
5048 CALL wrf_message("BEFORE MOSAIC_INIT")
5049 WRITE(message_text,fmt='(a,2I6,2F8.2,2I6)') 'I,J,xland,xice,mosaic_cat_index,ivgtyp = ', &
5050 I,J,xland(i,j),xice(i,j),mosaic_cat_index(I,1,J),IVGTYP(i,j)
5051 CALL wrf_message(message_text)
5059 LANDUSEF2(i,mosaic_i,j)=LANDUSEF(i,mosaic_i,j)
5060 mosaic_cat_index(i,mosaic_i,j)=mosaic_i
5071 IF (NumPairs == 0) EXIT
5073 DO mosaic_i=1, NumPairs
5074 IF(LANDUSEF2(i,mosaic_i, j) < LANDUSEF2(i,mosaic_i+1, j) ) THEN
5075 Temp = LANDUSEF2(i,mosaic_i, j)
5076 LANDUSEF2(i,mosaic_i, j)=LANDUSEF2(i,mosaic_i+1, j)
5077 LANDUSEF2(i,mosaic_i+1, j)=Temp
5080 Temp2 = mosaic_cat_index(i,mosaic_i,j)
5081 mosaic_cat_index(i,mosaic_i,j)=mosaic_cat_index(i,mosaic_i+1,j)
5082 mosaic_cat_index(i,mosaic_i+1,j)=Temp2
5085 NumPairs = LastSwap - 1
5091 !===========================================================================
5092 ! For non-seaice grids, eliminate the seaice-tiles
5093 !===========================================================================
5098 IF (XLAND(I,J).LT.1.5) THEN
5101 IF( XICE(I,J).GE. XICE_THRESHOLD ) THEN
5102 WRITE (message_text,fmt='(a,2I5)') 'sea-ice at point, I and J = ', i,j
5103 CALL wrf_message(message_text)
5107 IF (ICE == 1) Then ! sea-ice case , eliminate sea-ice if they are not the dominant ones
5109 IF (IVGTYP(i,j) == isice) THEN ! if this grid cell is dominanted by ice, then do nothing
5113 DO mosaic_i=2,mosaic_cat
5114 IF (mosaic_cat_index(i,mosaic_i,j) == isice ) THEN
5115 Temp4=LANDUSEF2(i,mosaic_i,j)
5116 Temp5=mosaic_cat_index(i,mosaic_i,j)
5118 LANDUSEF2(i,mosaic_i:NLCAT-1,j)=LANDUSEF2(i,mosaic_i+1:NLCAT,j)
5119 mosaic_cat_index(i,mosaic_i:NLCAT-1,j)=mosaic_cat_index(i,mosaic_i+1:NLCAT,j)
5121 LANDUSEF2(i,NLCAT,j)=Temp4
5122 mosaic_cat_index(i,NLCAT,j)=Temp5
5126 ENDIF ! for (IVGTYP(i,j) == isice )
5128 ELSEIF (ICE ==0) THEN
5130 IF ((mosaic_cat_index(I,1,J) .EQ. ISWATER)) THEN
5132 ! xland < 1.5 but the dominant land use category based on our calculation is water
5134 IF (IVGTYP(i,j) .EQ. ISWATER) THEN
5136 ! xland < 1.5 but the dominant land use category based on the geogrid calculation is water, this must be wrong
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 "In addition, the dominant land use category based on the geogrid calculation is water, this must be wrong")
5144 ENDIF ! for (IVGTYP(i,j) .EQ. ISWATER)
5146 IF (IVGTYP(i,j) .NE. ISWATER) THEN
5148 ! 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
5150 Temp4=LANDUSEF2(i,1,j)
5151 Temp5=mosaic_cat_index(i,1,j)
5153 LANDUSEF2(i,1:NLCAT-1,j)=LANDUSEF2(i,2:NLCAT,j)
5154 mosaic_cat_index(i,1:NLCAT-1,j)=mosaic_cat_index(i,2:NLCAT,j)
5156 LANDUSEF2(i,NLCAT,j)=Temp4
5157 mosaic_cat_index(i,NLCAT,j)=Temp5
5159 CALL wrf_message("IN MOSAIC_INIT")
5160 WRITE(message_text,fmt='(a,3I6,2F8.2)') 'I,J,IVGTYP,XLAND,XICE = ',I,J,IVGTYP(I,J),xland(i,j),xice(i,j)
5161 CALL wrf_message(message_text)
5162 CALL wrf_message("xland < 1.5 but the dominant land use category based on our calculation is water."//&
5163 "this is fine as long as we change our calculation so that the dominant land use category is"//&
5164 "stwiched back to not water.")
5165 WRITE(message_text,fmt='(a,2I6)') 'land use category has been switched, before and after values are ', &
5166 temp5,mosaic_cat_index(i,1,j)
5167 CALL wrf_message(message_text)
5168 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)
5169 CALL wrf_message(message_text)
5171 ENDIF ! for (IVGTYP(i,j) .NE. ISWATER)
5173 ELSE ! for (mosaic_cat_index(I,1,J) .EQ. ISWATER)
5175 DO mosaic_i=2,mosaic_cat
5176 IF (mosaic_cat_index(i,mosaic_i,j) == iswater ) THEN
5177 Temp4=LANDUSEF2(i,mosaic_i,j)
5178 Temp5=mosaic_cat_index(i,mosaic_i,j)
5180 LANDUSEF2(i,mosaic_i:NLCAT-1,j)=LANDUSEF2(i,mosaic_i+1:NLCAT,j)
5181 mosaic_cat_index(i,mosaic_i:NLCAT-1,j)=mosaic_cat_index(i,mosaic_i+1:NLCAT,j)
5183 LANDUSEF2(i,NLCAT,j)=Temp4
5184 mosaic_cat_index(i,NLCAT,j)=Temp5
5188 ENDIF ! for (mosaic_cat_index(I,1,J) .EQ. ISWATER)
5190 ENDIF ! for ICE == 1
5192 ELSE ! FOR (XLAND(I,J).LT.1.5)
5196 IF( XICE(I,J).GE. XICE_THRESHOLD ) THEN
5197 WRITE (message_text,fmt='(a,2I6)') 'sea-ice at water point, I and J = ', i,j
5198 CALL wrf_message(message_text)
5202 IF ((mosaic_cat_index(I,1,J) .NE. ISWATER)) THEN
5204 ! xland > 1.5 and the dominant land use category based on our calculation is not water
5206 IF (IVGTYP(i,j) .NE. ISWATER) THEN
5208 ! xland > 1.5 but the dominant land use category based on the geogrid calculation is not water, this must be wrong
5209 CALL wrf_message("IN MOSAIC_INIT")
5210 WRITE(message_text,fmt='(a,3I6,2F8.2)') 'I,J,IVGTYP,XLAND,XICE = ',I,J,IVGTYP(I,J),xland(i,j),xice(i,j)
5211 CALL wrf_message(message_text)
5212 CALL wrf_message("xland > 1.5 but the dominant land use category based on our calculation is not water."// &
5213 "in addition, the dominant land use category based on the geogrid calculation is not water,"// &
5214 "this must be wrong.")
5215 ENDIF ! for (IVGTYP(i,j) .NE. ISWATER)
5217 IF (IVGTYP(i,j) .EQ. ISWATER) THEN
5219 ! 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
5221 CALL wrf_message("IN MOSAIC_INIT")
5222 WRITE(message_text,fmt='(a,3I6,2F8.2)') 'I,J,IVGTYP,XLAND,XICE = ',I,J,IVGTYP(I,J),xland(i,j),xice(i,j)
5223 CALL wrf_message(message_text)
5224 CALL wrf_message("xland > 1.5 but the dominant land use category based on our calculation is not water."// &
5225 "however, the dominant land use category based on the geogrid calculation is water")
5226 CALL wrf_message("This is fine. We do not need to do anyting because in the noaddrv, "//&
5227 "we use xland as a criterion for whether using"// &
5228 "mosaic or not when xland > 1.5, no mosaic will be used anyway")
5230 ENDIF ! for (IVGTYP(i,j) .NE. ISWATER)
5232 ENDIF ! for (mosaic_cat_index(I,1,J) .NE. ISWATER)
5234 ENDIF ! FOR (XLAND(I,J).LT.1.5)
5239 !===========================================================================
5241 !===========================================================================
5248 DO mosaic_i=1,mosaic_cat
5249 Temp6=Temp6+LANDUSEF2(i,mosaic_i,j)
5252 if (Temp6 .LT. 1e-5) then
5255 WRITE (message_text,fmt='(a,e8.1)') 'the total land surface fraction is less than ', temp6
5256 CALL wrf_message(message_text)
5257 WRITE (message_text,fmt='(a,2I6,4F8.2)') 'some landusef values at i,j are ', &
5258 i,j,landusef2(i,1,j),landusef2(i,2,j),landusef2(i,3,j),landusef2(i,4,j)
5259 CALL wrf_message(message_text)
5260 WRITE (message_text,fmt='(a,2I6,3I6)') 'some mosaic cat values at i,j are ', &
5261 i,j,mosaic_cat_index(i,1,j),mosaic_cat_index(i,2,j),mosaic_cat_index(i,3,j)
5262 CALL wrf_message(message_text)
5266 LANDUSEF2(i,1:mosaic_cat, j)=LANDUSEF2(i,1:mosaic_cat,j)*(1/Temp6)
5271 !===========================================================================
5272 ! initilize the variables
5273 !===========================================================================
5278 DO mosaic_i=1,mosaic_cat
5280 TSK_mosaic(i,mosaic_i,j)=TSK(i,j)
5281 CANWAT_mosaic(i,mosaic_i,j)=CANWAT(i,j)
5282 SNOW_mosaic(i,mosaic_i,j)=SNOW(i,j)
5283 SNOWH_mosaic(i,mosaic_i,j)=SNOWH(i,j)
5284 SNOWC_mosaic(i,mosaic_i,j)=SNOWC(i,j)
5286 ALBEDO_mosaic(i,mosaic_i,j)=ALBEDO(i,j)
5287 ALBBCK_mosaic(i,mosaic_i,j)=ALBBCK(i,j)
5288 EMISS_mosaic(i,mosaic_i,j)=EMISS(i,j)
5289 EMBCK_mosaic(i,mosaic_i,j)=EMBCK(i,j)
5290 ZNT_mosaic(i,mosaic_i,j)=Z0(i,j)
5291 Z0_mosaic(i,mosaic_i,j)=Z0(i,j)
5293 DO soil_k=1,num_soil_layers
5295 TSLB_mosaic(i,num_soil_layers*(mosaic_i-1)+soil_k,j)=TSLB(i,soil_k,j)
5296 SMOIS_mosaic(i,num_soil_layers*(mosaic_i-1)+soil_k,j)=SMOIS(i,soil_k,j)
5297 SH2O_mosaic(i,num_soil_layers*(mosaic_i-1)+soil_k,j)=SH2O(i,soil_k,j)
5301 TR_URB2D_mosaic(i,mosaic_i,j)=TSK(i,j)
5302 TB_URB2D_mosaic(i,mosaic_i,j)=TSK(i,j)
5303 TG_URB2D_mosaic(i,mosaic_i,j)=TSK(i,j)
5304 TC_URB2D_mosaic(i,mosaic_i,j)=TSK(i,j)
5305 TS_URB2D_mosaic(i,mosaic_i,j)=TSK(i,j)
5306 TS_RUL2D_mosaic(i,mosaic_i,j)=TSK(i,j)
5307 QC_URB2D_mosaic(i,mosaic_i,j)=0.01
5308 SH_URB2D_mosaic(i,mosaic_i,j)=0
5309 LH_URB2D_mosaic(i,mosaic_i,j)=0
5310 G_URB2D_mosaic(i,mosaic_i,j)=0
5311 RN_URB2D_mosaic(i,mosaic_i,j)=0
5313 TRL_URB3D_mosaic(I,4*(mosaic_i-1)+1,J)=TSLB(I,1,J)+0.
5314 TRL_URB3D_mosaic(I,4*(mosaic_i-1)+2,J)=0.5*(TSLB(I,1,J)+TSLB(I,2,J))
5315 TRL_URB3D_mosaic(I,4*(mosaic_i-1)+3,J)=TSLB(I,2,J)+0.
5316 TRL_URB3D_mosaic(I,4*(mosaic_i-1)+4,J)=TSLB(I,2,J)+(TSLB(I,3,J)-TSLB(I,2,J))*0.29
5318 TBL_URB3D_mosaic(I,4*(mosaic_i-1)+1,J)=TSLB(I,1,J)+0.
5319 TBL_URB3D_mosaic(I,4*(mosaic_i-1)+2,J)=0.5*(TSLB(I,1,J)+TSLB(I,2,J))
5320 TBL_URB3D_mosaic(I,4*(mosaic_i-1)+3,J)=TSLB(I,2,J)+0.
5321 TBL_URB3D_mosaic(I,4*(mosaic_i-1)+4,J)=TSLB(I,2,J)+(TSLB(I,3,J)-TSLB(I,2,J))*0.29
5323 TGL_URB3D_mosaic(I,4*(mosaic_i-1)+1,J)=TSLB(I,1,J)
5324 TGL_URB3D_mosaic(I,4*(mosaic_i-1)+2,J)=TSLB(I,2,J)
5325 TGL_URB3D_mosaic(I,4*(mosaic_i-1)+3,J)=TSLB(I,3,J)
5326 TGL_URB3D_mosaic(I,4*(mosaic_i-1)+4,J)=TSLB(I,4,J)
5337 IF ((xland(i,j).LT. 1.5 ) .AND. (mosaic_cat_index(I,1,J) .EQ. ISWATER)) THEN
5338 CALL wrf_message("After MOSAIC_INIT")
5339 WRITE (message_text,fmt='(a,2I6,2F8.2,2I6)') 'weird xland,xice,mosaic_cat_index and ivgtyp at I,J = ', &
5340 i,j,xland(i,j),xice(i,j),mosaic_cat_index(I,1,J),IVGTYP(i,j)
5341 CALL wrf_message(message_text)
5347 ENDIF ! for not restart
5349 !--------------------------------
5350 END SUBROUTINE lsm_mosaic_init
5351 !--------------------------------
5354 END MODULE module_sf_noahdrv