module_sf_noahlsm.F
References to this file elsewhere.
1 MODULE module_sf_noahlsm
2
3 USE module_model_constants
4
5 !-------------------------------
6 USE module_sf_urban
7 !-------------------------------
8
9 ! REAL, PARAMETER :: CP = 1004.5
10 REAL, PARAMETER :: RD = 287.04, SIGMA = 5.67E-8, &
11 CPH2O = 4.218E+3,CPICE = 2.106E+3, &
12 LSUBF = 3.335E+5
13
14 ! VEGETATION PARAMETERS
15 INTEGER :: LUCATS , BARE
16 integer, PARAMETER :: NLUS=50
17 CHARACTER*4 LUTYPE
18 INTEGER, DIMENSION(1:NLUS) :: NROTBL
19 real, dimension(1:NLUS) :: SNUPTBL, RSTBL, RGLTBL, HSTBL, LAITBL, &
20 ALBTBL, Z0TBL, SHDTBL, MAXALB
21 REAL :: TOPT_DATA,CMCMAX_DATA,CFACTR_DATA,RSMAX_DATA
22
23 ! SOIL PARAMETERS
24 INTEGER :: SLCATS
25 INTEGER, PARAMETER :: NSLTYPE=30
26 CHARACTER*4 SLTYPE
27 REAL, DIMENSION (1:NSLTYPE) :: BB,DRYSMC,F11, &
28 MAXSMC, REFSMC,SATPSI,SATDK,SATDW, WLTSMC,QTZ
29
30 ! LSM GENERAL PARAMETERS
31 INTEGER :: SLPCATS
32 INTEGER, PARAMETER :: NSLOPE=30
33 REAL, DIMENSION (1:NSLOPE) :: SLOPE_DATA
34 REAL :: SBETA_DATA,FXEXP_DATA,CSOIL_DATA,SALP_DATA,REFDK_DATA, &
35 REFKDT_DATA,FRZK_DATA,ZBOT_DATA, SMLOW_DATA,SMHIGH_DATA, &
36 CZIL_DATA
37
38 CHARACTER*256 :: err_message
39
40 !
41 CONTAINS
42 !
43 !----------------------------------------------------------------
44 ! Urban related variable are added to arguments - urban
45 !----------------------------------------------------------------
46 SUBROUTINE lsm(DZ8W,QV3D,P8W3D,T3D,TSK, &
47 HFX,QFX,LH,GRDFLX, QGH,GSW,SWDOWN,GLW,SMSTAV,SMSTOT, &
48 SFCRUNOFF, UDRUNOFF,IVGTYP,ISLTYP,VEGFRA, &
49 ALBEDO,ALBBCK,ZNT,Z0,TMN,XLAND,XICE, EMISS, &
50 SNOWC,QSFC,RAINBL, &
51 num_soil_layers,DT,DZS,ITIMESTEP, &
52 SMOIS,TSLB,SNOW,CANWAT, &
53 CHS,CHS2,CQS2,CPM,ROVCP,SR,chklowq,qz0, & !H
54 myj,frpcpn, &
55 SH2O,SNOWH, & !H
56 U_PHY,V_PHY, & !I
57 SNOALB,SHDMIN,SHDMAX, & !I
58 ACSNOM,ACSNOW, & !O
59 SNOPCX, & !O
60 ! MEK JUL2007
61 POTEVP, & !O
62 ids,ide, jds,jde, kds,kde, &
63 ims,ime, jms,jme, kms,kme, &
64 its,ite, jts,jte, kts,kte, &
65 ucmcall, &
66 !Optional Urban
67 TR_URB2D,TB_URB2D,TG_URB2D,TC_URB2D,QC_URB2D, & !H urban
68 UC_URB2D, & !H urban
69 XXXR_URB2D,XXXB_URB2D,XXXG_URB2D,XXXC_URB2D, & !H urban
70 TRL_URB3D,TBL_URB3D,TGL_URB3D, & !H urban
71 SH_URB2D,LH_URB2D,G_URB2D,RN_URB2D,TS_URB2D, & !H urban
72 PSIM_URB2D,PSIH_URB2D,U10_URB2D,V10_URB2D, & !O urban
73 GZ1OZ0_URB2D, AKMS_URB2D, & !O urban
74 TH2_URB2D,Q2_URB2D, UST_URB2D, & !O urban
75 DECLIN_URB,COSZ_URB2D,OMG_URB2D, & !I urban
76 XLAT_URB2D, & !I urban
77 num_roof_layers, num_wall_layers, & !I urban
78 num_road_layers, DZR, DZB, DZG, & !I urban
79 FRC_URB2D,UTYPE_URB2D) !O
80 !----------------------------------------------------------------
81 IMPLICIT NONE
82 !----------------------------------------------------------------
83 !----------------------------------------------------------------
84 ! --- atmospheric (WRF generic) variables
85 !-- DT time step (seconds)
86 !-- DZ8W thickness of layers (m)
87 !-- T3D temperature (K)
88 !-- QV3D 3D water vapor mixing ratio (Kg/Kg)
89 !-- P3D 3D pressure (Pa)
90 !-- FLHC exchange coefficient for heat (m/s)
91 !-- FLQC exchange coefficient for moisture (m/s)
92 !-- PSFC surface pressure (Pa)
93 !-- XLAND land mask (1 for land, 2 for water)
94 !-- QGH saturated mixing ratio at 2 meter
95 !-- GSW downward short wave flux at ground surface (W/m^2)
96 !-- GLW downward long wave flux at ground surface (W/m^2)
97 !-- History variables
98 !-- CANWAT canopy moisture content (mm)
99 !-- TSK surface temperature (K)
100 !-- TSLB soil temp (k)
101 !-- SMOIS total soil moisture content (volumetric fraction)
102 !-- SH2O unfrozen soil moisture content (volumetric fraction)
103 ! note: frozen soil moisture (i.e., soil ice) = SMOIS - SH2O
104 !-- SNOWH actual snow depth (m)
105 !-- SNOW liquid water-equivalent snow depth (m)
106 !-- ALBEDO time-varying surface albedo including snow effect (unitless fraction)
107 !-- ALBBCK background surface albedo (unitless fraction)
108 !-- CHS surface exchange coefficient for heat and moisture (m s-1);
109 !-- CHS2 2m surface exchange coefficient for heat (m s-1);
110 !-- CQS2 2m surface exchange coefficient for moisture (m s-1);
111 ! --- soil variables
112 !-- num_soil_layers the number of soil layers
113 !-- ZS depths of centers of soil layers (m)
114 !-- DZS thicknesses of soil layers (m)
115 !-- SLDPTH thickness of each soil layer (m, same as DZS)
116 !-- TMN soil temperature at lower boundary (K)
117 !-- SMCWLT wilting point (volumetric)
118 !-- SMCDRY dry soil moisture threshold where direct evap from
119 ! top soil layer ends (volumetric)
120 !-- SMCREF soil moisture threshold below which transpiration begins to
121 ! stress (volumetric)
122 !-- SMCMAX porosity, i.e. saturated value of soil moisture (volumetric)
123 !-- NROOT number of root layers, a function of veg type, determined
124 ! in subroutine redprm.
125 !-- SMSTAV Soil moisture availability for evapotranspiration (
126 ! fraction between SMCWLT and SMCMXA)
127 !-- SMSTOT Total soil moisture content frozen+unfrozen) in the soil column (mm)
128 ! --- snow variables
129 !-- SNOWC fraction snow coverage (0-1.0)
130 ! --- vegetation variables
131 !-- SNOALB upper bound on maximum albedo over deep snow
132 !-- SHDMIN minimum areal fractional coverage of annual green vegetation
133 !-- SHDMAX maximum areal fractional coverage of annual green vegetation
134 !-- XLAI leaf area index (dimensionless)
135 !-- Z0BRD Background fixed roughness length (M)
136 !-- Z0 Background vroughness length (M) as function
137 !-- ZNT Time varying roughness length (M) as function
138 !-- ALBD(IVGTPK,ISN) background albedo reading from a table
139 ! --- LSM output
140 !-- HFX upward heat flux at the surface (W/m^2)
141 !-- QFX upward moisture flux at the surface (kg/m^2/s)
142 !-- LH upward moisture flux at the surface (W m-2)
143 !-- GRDFLX(I,J) ground heat flux (W m-2)
144 !-- FDOWN radiation forcing at the surface (W m-2) = SOLDN*(1-alb)+LWDN
145 !----------------------------------------------------------------------------
146 !-- EC canopy water evaporation ((W m-2)
147 !-- EDIR direct soil evaporation (W m-2)
148 !-- ET plant transpiration from a particular root layer (W m-2)
149 !-- ETT total plant transpiration (W m-2)
150 !-- ESNOW sublimation from (or deposition to if <0) snowpack (W m-2)
151 !-- DRIP through-fall of precip and/or dew in excess of canopy
152 ! water-holding capacity (m)
153 !-- DEW dewfall (or frostfall for t<273.15) (M)
154 ! ----------------------------------------------------------------------
155 !-- BETA ratio of actual/potential evap (dimensionless)
156 !-- ETP potential evaporation (W m-2)
157 ! ----------------------------------------------------------------------
158 !-- FLX1 precip-snow sfc (W m-2)
159 !-- FLX2 freezing rain latent heat flux (W m-2)
160 !-- FLX3 phase-change heat flux from snowmelt (W m-2)
161 ! ----------------------------------------------------------------------
162 !-- ACSNOM snow melt (mm) (water equivalent)
163 !-- ACSNOW accumulated snow fall (mm) (water equivalent)
164 !-- SNOPCX snow phase change heat flux (W/m^2)
165 !-- POTEVP accumulated potential evaporation (W/m^2)
166 ! ----------------------------------------------------------------------
167 !-- RUNOFF1 surface runoff (m s-1), not infiltrating the surface
168 !-- RUNOFF2 subsurface runoff (m s-1), drainage out bottom of last
169 ! soil layer (baseflow)
170 ! important note: here RUNOFF2 is actually the sum of RUNOFF2 and RUNOFF3
171 !-- RUNOFF3 numerical trunctation in excess of porosity (smcmax)
172 ! for a given soil layer at the end of a time step (m s-1).
173 ! ----------------------------------------------------------------------
174 !-- RC canopy resistance (s m-1)
175 !-- PC plant coefficient (unitless fraction, 0-1) where PC*ETP = actual transp
176 !-- RSMIN minimum canopy resistance (s m-1)
177 !-- RCS incoming solar rc factor (dimensionless)
178 !-- RCT air temperature rc factor (dimensionless)
179 !-- RCQ atmos vapor pressure deficit rc factor (dimensionless)
180 !-- RCSOIL soil moisture rc factor (dimensionless)
181
182 !-- EMISS surface emissivity (between 0 and 1)
183
184 !-- ROVCP R/CP
185 ! (R_d/R_v) (dimensionless)
186 !-- ids start index for i in domain
187 !-- ide end index for i in domain
188 !-- jds start index for j in domain
189 !-- jde end index for j in domain
190 !-- kds start index for k in domain
191 !-- kde end index for k in domain
192 !-- ims start index for i in memory
193 !-- ime end index for i in memory
194 !-- jms start index for j in memory
195 !-- jme end index for j in memory
196 !-- kms start index for k in memory
197 !-- kme end index for k in memory
198 !-- its start index for i in tile
199 !-- ite end index for i in tile
200 !-- jts start index for j in tile
201 !-- jte end index for j in tile
202 !-- kts start index for k in tile
203 !-- kte end index for k in tile
204 !
205 !-- SR fraction of frozen precip (0.0 to 1.0)
206 !----------------------------------------------------------------
207
208 ! IN only
209
210 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
211 ims,ime, jms,jme, kms,kme, &
212 its,ite, jts,jte, kts,kte
213
214 INTEGER, INTENT(IN ) :: ucmcall !urban
215
216 REAL, DIMENSION( ims:ime, jms:jme ) , &
217 INTENT(IN ) :: TMN, &
218 XLAND, &
219 XICE, &
220 VEGFRA, &
221 SHDMIN, &
222 SHDMAX, &
223 SNOALB, &
224 GSW, &
225 SWDOWN, & !added 10 jan 2007
226 GLW, &
227 Z0, &
228 ALBBCK, &
229 RAINBL, &
230 EMISS, &
231 SR
232
233
234 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
235 INTENT(IN ) :: QV3D, &
236 p8w3D, &
237 DZ8W, &
238 T3D
239 REAL, DIMENSION( ims:ime, jms:jme ) , &
240 INTENT(IN ) :: QGH, &
241 CHS, &
242 CPM
243
244 INTEGER, DIMENSION( ims:ime, jms:jme ) , &
245 INTENT(IN ) :: IVGTYP, &
246 ISLTYP
247
248 INTEGER, INTENT(IN) :: num_soil_layers,ITIMESTEP
249
250 REAL, INTENT(IN ) :: DT,ROVCP
251
252 REAL, DIMENSION(1:num_soil_layers), INTENT(IN)::DZS
253
254 ! IN and OUT
255
256 REAL, DIMENSION( ims:ime , 1:num_soil_layers, jms:jme ), &
257 INTENT(INOUT) :: SMOIS, & ! total soil moisture
258 SH2O, & ! new soil liquid
259 TSLB ! TSLB STEMP
260
261 REAL, DIMENSION( ims:ime, jms:jme ) , &
262 INTENT(INOUT) :: TSK, & !was TGB (temperature)
263 HFX, &
264 QFX, &
265 LH, &
266 GRDFLX, &
267 QSFC,&
268 CQS2,&
269 CHS2,&
270 SNOW, &
271 SNOWC, &
272 SNOWH, & !new
273 CANWAT, &
274 SMSTAV, &
275 SMSTOT, &
276 SFCRUNOFF, &
277 UDRUNOFF, &
278 ACSNOM, &
279 ACSNOW, &
280 SNOPCX, &
281 ! MEK JUL2007
282 POTEVP, &
283 ALBEDO, &
284 ZNT
285
286 REAL, DIMENSION( ims:ime, jms:jme ) , &
287 INTENT(OUT) :: CHKLOWQ
288 REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: QZ0
289
290 ! Local variables (moved here from driver to make routine thread safe, 20031007 jm)
291
292 REAL, DIMENSION(1:num_soil_layers) :: ET
293
294 REAL :: BETA, ETP, SSOIL,EC, EDIR, ESNOW, ETT, &
295 FLX1,FLX2,FLX3, DRIP,DEW,FDOWN,RC,PC,RSMIN,XLAI, &
296 ! RCS,RCT,RCQ,RCSOIL
297 RCS,RCT,RCQ,RCSOIL,FFROZP
298
299 LOGICAL, INTENT(IN ) :: myj,frpcpn
300
301 ! DECLARATIONS - LOGICAL
302 ! ----------------------------------------------------------------------
303 LOGICAL, PARAMETER :: LOCAL=.false.
304 LOGICAL :: FRZGRA, SNOWNG
305
306 LOGICAL :: IPRINT
307
308 ! ----------------------------------------------------------------------
309 ! DECLARATIONS - INTEGER
310 ! ----------------------------------------------------------------------
311 INTEGER :: I,J, ICE,NSOIL,SLOPETYP,SOILTYP,VEGTYP
312 INTEGER :: NROOT
313 INTEGER :: KZ ,K
314 INTEGER :: NS
315 ! ----------------------------------------------------------------------
316 ! DECLARATIONS - REAL
317 ! ----------------------------------------------------------------------
318
319 REAL :: SHMIN,SHMAX,DQSDT2,LWDN,PRCP,PRCPRAIN, &
320 Q2SAT,Q2SATI,SFCPRS,SFCSPD,SFCTMP,SHDFAC,SNOALB1, &
321 SOLDN,TBOT,ZLVL, Q2K,ALBBRD, ALBEDOK, ETA, ETA_KINEMATIC, &
322 Z0K,RUNOFF1,RUNOFF2,RUNOFF3,SHEAT,SOLNET,E2SAT,SFCTSNO, &
323 ! mek, WRF testing, expanded diagnostics
324 SOLUP,LWUP,RNET,RES,Q1SFC,TAIRV,RHO,SATFLG
325 ! MEK MAY 2007
326 REAL :: FDTLIW
327 ! MEK JUL2007 for pot. evap.
328 REAL :: FDTW
329
330 REAL :: EMISSI
331
332 REAL :: SNCOVR,SNEQV,SNOWHK,CMC, CHK,TH2
333
334 REAL :: SMCDRY,SMCMAX,SMCREF,SMCWLT,SNOMLT,SOILM,SOILW,Q1,T1
335
336 REAL :: DUMMY,Z0BRD
337 !
338 REAL :: COSZ, SOLARDIRECT
339 !
340 REAL, DIMENSION(1:num_soil_layers):: SLDPTH, STC,SMC,SWC
341 !
342 REAL, DIMENSION(1:num_soil_layers) :: ZSOIL, RTDIS
343 REAL, PARAMETER :: TRESH=.95E0, A2=17.67,A3=273.15,A4=29.65, &
344 T0=273.16E0, ELWV=2.50E6, A23M4=A2*(A3-A4)
345 ! MEK MAY 2007
346 REAL, PARAMETER :: ROW=1.E3,ELIW=XLF,ROWLIW=ROW*ELIW
347
348 ! ----------------------------------------------------------------------
349 ! DECLARATIONS START - urban
350 ! ----------------------------------------------------------------------
351
352 ! input variables surface_driver --> lsm
353 INTEGER, INTENT(IN) :: num_roof_layers
354 INTEGER, INTENT(IN) :: num_wall_layers
355 INTEGER, INTENT(IN) :: num_road_layers
356 REAL, OPTIONAL, DIMENSION(1:num_roof_layers), INTENT(IN) :: DZR
357 REAL, OPTIONAL, DIMENSION(1:num_wall_layers), INTENT(IN) :: DZB
358 REAL, OPTIONAL, DIMENSION(1:num_road_layers), INTENT(IN) :: DZG
359 REAL, OPTIONAL, INTENT(IN) :: DECLIN_URB
360 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: COSZ_URB2D
361 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: OMG_URB2D
362 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: XLAT_URB2D
363 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: U_PHY
364 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: V_PHY
365
366 ! input variables lsm --> urban
367 INTEGER :: UTYPE_URB ! urban type [urban=1, suburban=2, rural=3]
368 REAL :: TA_URB ! potential temp at 1st atmospheric level [K]
369 REAL :: QA_URB ! mixing ratio at 1st atmospheric level [kg/kg]
370 REAL :: UA_URB ! wind speed at 1st atmospheric level [m/s]
371 REAL :: U1_URB ! u at 1st atmospheric level [m/s]
372 REAL :: V1_URB ! v at 1st atmospheric level [m/s]
373 REAL :: SSG_URB ! downward total short wave radiation [W/m/m]
374 REAL :: LLG_URB ! downward long wave radiation [W/m/m]
375 REAL :: RAIN_URB ! precipitation [mm/h]
376 REAL :: RHOO_URB ! air density [kg/m^3]
377 REAL :: ZA_URB ! first atmospheric level [m]
378 REAL :: DELT_URB ! time step [s]
379 REAL :: SSGD_URB ! downward direct short wave radiation [W/m/m]
380 REAL :: SSGQ_URB ! downward diffuse short wave radiation [W/m/m]
381 REAL :: XLAT_URB ! latitude [deg]
382 REAL :: COSZ_URB ! cosz
383 REAL :: OMG_URB ! hour angle
384 REAL :: ZNT_URB ! roughness length [m]
385 REAL :: TR_URB
386 REAL :: TB_URB
387 REAL :: TG_URB
388 REAL :: TC_URB
389 REAL :: QC_URB
390 REAL :: UC_URB
391 REAL :: XXXR_URB
392 REAL :: XXXB_URB
393 REAL :: XXXG_URB
394 REAL :: XXXC_URB
395 REAL, DIMENSION(1:num_roof_layers) :: TRL_URB ! roof layer temp [K]
396 REAL, DIMENSION(1:num_wall_layers) :: TBL_URB ! wall layer temp [K]
397 REAL, DIMENSION(1:num_road_layers) :: TGL_URB ! road layer temp [K]
398 LOGICAL :: LSOLAR_URB
399 ! state variable surface_driver <--> lsm <--> urban
400 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TR_URB2D
401 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TB_URB2D
402 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TG_URB2D
403 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TC_URB2D
404 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: QC_URB2D
405 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: UC_URB2D
406 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXR_URB2D
407 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXB_URB2D
408 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXG_URB2D
409 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXC_URB2D
410 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: SH_URB2D
411 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: LH_URB2D
412 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: G_URB2D
413 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: RN_URB2D
414 !
415 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TS_URB2D
416
417 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_roof_layers, jms:jme ), INTENT(INOUT) :: TRL_URB3D
418 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_wall_layers, jms:jme ), INTENT(INOUT) :: TBL_URB3D
419 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_road_layers, jms:jme ), INTENT(INOUT) :: TGL_URB3D
420
421 ! output variable lsm --> surface_driver
422 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: PSIM_URB2D
423 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: PSIH_URB2D
424 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: GZ1OZ0_URB2D
425 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: U10_URB2D
426 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: V10_URB2D
427 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: TH2_URB2D
428 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: Q2_URB2D
429 !
430 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: AKMS_URB2D
431 !
432 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: UST_URB2D
433 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: FRC_URB2D
434 INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: UTYPE_URB2D
435
436
437 ! output variables urban --> lsm
438 REAL :: TS_URB ! surface radiative temperature [K]
439 REAL :: QS_URB ! surface humidity [-]
440 REAL :: SH_URB ! sensible heat flux [W/m/m]
441 REAL :: LH_URB ! latent heat flux [W/m/m]
442 REAL :: LH_KINEMATIC_URB ! latent heat flux, kinetic [kg/m/m/s]
443 REAL :: SW_URB ! upward short wave radiation flux [W/m/m]
444 REAL :: ALB_URB ! time-varying albedo [fraction]
445 REAL :: LW_URB ! upward long wave radiation flux [W/m/m]
446 REAL :: G_URB ! heat flux into the ground [W/m/m]
447 REAL :: RN_URB ! net radiation [W/m/m]
448 REAL :: PSIM_URB ! shear f for momentum [-]
449 REAL :: PSIH_URB ! shear f for heat [-]
450 REAL :: GZ1OZ0_URB ! shear f for heat [-]
451 REAL :: U10_URB ! wind u component at 10 m [m/s]
452 REAL :: V10_URB ! wind v component at 10 m [m/s]
453 REAL :: TH2_URB ! potential temperature at 2 m [K]
454 REAL :: Q2_URB ! humidity at 2 m [-]
455 REAL :: CHS_URB
456 REAL :: CHS2_URB
457 REAL :: UST_URB
458
459 ! ----------------------------------------------------------------------
460 ! DECLARATIONS END - urban
461 ! ----------------------------------------------------------------------
462
463 REAL, PARAMETER :: CAPA=R_D/CP
464 REAL :: APELM,APES,SFCTH2,PSFC
465
466 ! PRINT *,'THIS IS UNIFIED NOAH LSM'
467
468 ! MEK MAY 2007
469 FDTLIW=DT/ROWLIW
470 ! MEK JUL2007
471 FDTW=DT/(XLV*RHOWATER)
472 ! debug printout
473 IPRINT=.false.
474
475 ! SLOPETYP=2
476 SLOPETYP=1
477 ! SHDMIN=0.00
478
479
480 NSOIL=num_soil_layers
481
482 DO NS=1,NSOIL
483 SLDPTH(NS)=DZS(NS)
484 ENDDO
485
486 DO J=jts,jte
487
488 IF(ITIMESTEP.EQ.1)THEN
489 DO 50 I=its,ite
490 !*** initialize soil conditions for IHOP 31 May case
491 ! IF((XLAND(I,J)-1.5) < 0.)THEN
492 ! if (I==108.and.j==85) then
493 ! DO NS=1,NSOIL
494 ! SMOIS(I,NS,J)=0.10
495 ! SH2O(I,NS,J)=0.10
496 ! enddo
497 ! endif
498 ! ENDIF
499
500 !*** SET ZERO-VALUE FOR SOME OUTPUT DIAGNOSTIC ARRAYS
501 IF((XLAND(I,J)-1.5).GE.0.)THEN
502 ! check sea-ice point
503 IF(XICE(I,J).EQ.1..and.IPRINT)PRINT*,' sea-ice at water point, I=',I, &
504 'J=',J
505 !*** Open Water Case
506 SMSTAV(I,J)=1.0
507 SMSTOT(I,J)=1.0
508 DO NS=1,NSOIL
509 SMOIS(I,NS,J)=1.0
510 TSLB(I,NS,J)=273.16 !STEMP
511 ENDDO
512 ELSE
513 IF(XICE(I,J).EQ.1.)THEN
514 !*** SEA-ICE CASE
515 SMSTAV(I,J)=1.0
516 SMSTOT(I,J)=1.0
517 DO NS=1,NSOIL
518 SMOIS(I,NS,J)=1.0
519 ENDDO
520 ENDIF
521 ENDIF
522 !
523 50 CONTINUE
524 ENDIF ! end of initialization over ocean
525
526 !-----------------------------------------------------------------------
527 DO 100 I=its,ite
528 ! surface pressure
529 PSFC=P8w3D(i,1,j)
530 ! pressure in middle of lowest layer
531 SFCPRS=(P8W3D(I,KTS+1,j)+P8W3D(i,KTS,j))*0.5
532 ! convert from mixing ratio to specific humidity
533 Q2K=QV3D(i,1,j)/(1.0+QV3D(i,1,j))
534 !
535 ! Q2SAT=QGH(I,j)
536 Q2SAT=QGH(I,J)/(1.0+QGH(I,J)) ! Q2SAT is sp humidity
537 ! add check on myj=.true.
538 ! IF((Q2K.GE.Q2SAT*TRESH).AND.Q2K.LT.QZ0(I,J))THEN
539 IF((myj).AND.(Q2K.GE.Q2SAT*TRESH).AND.Q2K.LT.QZ0(I,J))THEN
540 SATFLG=0.
541 CHKLOWQ(I,J)=0.
542 ELSE
543 SATFLG=1.0
544 CHKLOWQ(I,J)=1.
545 ENDIF
546
547 SFCTMP=T3D(i,1,j)
548 ZLVL=0.5*DZ8W(i,1,j)
549
550 ! TH2=SFCTMP+(0.0097545*ZLVL)
551 ! calculate SFCTH2 via Exner function vs lapse-rate (above)
552 APES=(1.E5/PSFC)**CAPA
553 APELM=(1.E5/SFCPRS)**CAPA
554 SFCTH2=SFCTMP*APELM
555 TH2=SFCTH2/APES
556 !
557 EMISSI = EMISS(I,J)
558 LWDN=GLW(I,J)*EMISSI
559 ! SOLDN is total incoming solar
560 SOLDN=SWDOWN(I,J)
561 ! GSW is net downward solar
562 ! SOLNET=GSW(I,J)
563 ! use mid-day albedo to determine net downward solar (no solar zenith angle correction)
564 SOLNET=SOLDN*(1.-ALBEDO(I,J))
565 PRCP=RAINBL(i,j)/DT
566 VEGTYP=IVGTYP(I,J)
567 SOILTYP=ISLTYP(I,J)
568 SHDFAC=VEGFRA(I,J)/100.
569 T1=TSK(I,J)
570 CHK=CHS(I,J)
571 SHMIN=SHDMIN(I,J)/100. !NEW
572 SHMAX=SHDMAX(I,J)/100. !NEW
573 SNOALB1=SNOALB(I,J) !NEW
574 ! convert snow water equivalent from mm to meter
575 SNEQV=SNOW(I,J)*0.001
576 ! snow depth in meters
577 SNOWHK=SNOWH(I,J)
578
579 ! if "SR" present, set frac of frozen precip ("FFROZP") = snow-ratio ("SR", range:0-1)
580 ! SR from e.g. Ferrier microphysics
581 ! otherwise define from 1st atmos level temperature
582 IF(FRPCPN) THEN
583 FFROZP=SR(I,J)
584 ELSE
585 IF (SFCTMP <= 273.15) THEN
586 FFROZP = 1.0
587 ELSE
588 FFROZP = 0.0
589 ENDIF
590 ENDIF
591 !***
592 IF((XLAND(I,J)-1.5).GE.0.)THEN ! begining of land/sea if block
593 ! Open water points
594 ELSE
595 ! Land or sea-ice case
596
597 IF (XICE(I,J) .GT. 0.5) THEN
598 ICE=1
599 ELSE
600 ICE=0
601 ENDIF
602 DQSDT2=Q2SAT*A23M4/(SFCTMP-A4)**2
603
604 IF(SNOW(I,J).GT.0.0)THEN
605 ! snow on surface (use ice saturation properties)
606 SFCTSNO=SFCTMP
607 E2SAT=611.2*EXP(6174.*(1./273.15 - 1./SFCTSNO))
608 Q2SATI=0.622*E2SAT/(SFCPRS-E2SAT)
609 Q2SATI=Q2SATI/(1.0+Q2SATI) ! spec. hum.
610 IF(T1 .GT. 273.15)THEN
611 ! warm ground temps, weight the saturation between ice and water according to SNOWC
612 Q2SAT=Q2SAT*(1.-SNOWC(I,J)) + Q2SATI*SNOWC(I,J)
613 DQSDT2=DQSDT2*(1.-SNOWC(I,J)) + Q2SATI*6174./(SFCTSNO**2)*SNOWC(I,J)
614 ELSE
615 ! cold ground temps, use ice saturation only
616 Q2SAT=Q2SATI
617 DQSDT2=Q2SATI*6174./(SFCTSNO**2)
618 ENDIF
619 ! for snow cover fraction at 0 C, ground temp will not change, so DQSDT2 effectively zero
620 IF(T1 .GT. 273. .AND. SNOWC(I,J) .GT. 0.)DQSDT2=DQSDT2*(1.-SNOWC(I,J))
621 ENDIF
622
623 IF(ICE.EQ.0)THEN
624 TBOT=TMN(I,J)
625 ELSE
626 TBOT=271.16
627 ENDIF
628 IF(VEGTYP.EQ.25) SHDFAC=0.0000
629 IF(VEGTYP.EQ.26) SHDFAC=0.0000
630 IF(VEGTYP.EQ.27) SHDFAC=0.0000
631 IF(SOILTYP.EQ.14.AND.XICE(I,J).EQ.0.)THEN
632 IF(IPRINT)PRINT*,' SOIL TYPE FOUND TO BE WATER AT A LAND-POINT'
633 IF(IPRINT)PRINT*,i,j,'RESET SOIL in surfce.F'
634 SOILTYP=7
635 ENDIF
636 CMC=CANWAT(I,J)
637
638 !-------------------------------------------
639 !*** convert snow depth from mm to meter
640 !
641 ! IF(RDMAXALB) THEN
642 ! SNOALB=ALBMAX(I,J)*0.01
643 ! ELSE
644 ! SNOALB=MAXALB(IVGTPK)*0.01
645 ! ENDIF
646 ! IF(RDBRDALB) THEN
647 ! ALBBRD=ALBEDO(I,J)*0.01
648 ! ELSE
649 ! ALBBRD=ALBD(IVGTPK,ISN)*0.01
650 ! ENDIF
651
652 ! SNOALB1=0.80
653 ! SHMIN=0.00
654 ALBBRD=ALBBCK(I,J)
655 Z0BRD=Z0(I,J)
656 !FEI: temporaray arrays above need to be changed later by using SI
657
658 DO 70 NS=1,NSOIL
659 SMC(NS)=SMOIS(I,NS,J)
660 STC(NS)=TSLB(I,NS,J) !STEMP
661 SWC(NS)=SH2O(I,NS,J)
662 70 CONTINUE
663 !
664 if ( (SNEQV.ne.0..AND.SNOWHK.eq.0.).or.(SNOWHK.le.SNEQV) )THEN
665 SNOWHK= 5.*SNEQV
666 endif
667 !
668
669 !Fei: urban. for urban surface, if calling UCM, redefine urban as 5: Cropland/Grassland Mosaic
670
671 IF(UCMCALL == 1 ) THEN
672 IF( IVGTYP(I,J) == 1 .or. IVGTYP(I,J) == 31 .or. &
673 IVGTYP(I,J) == 32 .or. IVGTYP(I,J) == 33) THEN
674 VEGTYP = 5
675 SHDFAC = 0.8
676 ALBEDOK =0.2
677 ALBBRD =0.2
678 IF ( FRC_URB2D(I,J) < 0.99 ) THEN
679 T1= ( TSK(I,J) -FRC_URB2D(I,J) * TS_URB2D (I,J) )/ (1-FRC_URB2D(I,J))
680 ELSE
681 T1 = TSK(I,J)
682 ENDIF
683 ENDIF
684 ELSE
685 IF( IVGTYP(I,J) == 1 .or. IVGTYP(I,J) == 31 .or. &
686 IVGTYP(I,J) == 32 .or. IVGTYP(I,J) == 33) THEN
687 VEGTYP = 1
688 ENDIF
689 ENDIF
690
691 IF(IPRINT) THEN
692 !
693 print*, 'BEFORE SFLX, in Noahlsm_driver'
694 print*, 'ICE', ICE, 'DT',DT, 'ZLVL',ZLVL, 'NSOIL', NSOIL, &
695 'SLDPTH', SLDPTH, 'LOCAL',LOCAL, 'LUTYPE',&
696 LUTYPE, 'SLTYPE',SLTYPE, 'LWDN',LWDN, 'SOLDN',SOLDN, &
697 'SFCPRS',SFCPRS, 'PRCP',PRCP,'SFCTMP',SFCTMP,'Q2K',Q2K, &
698 'TH2',TH2,'Q2SAT',Q2SAT,'DQSDT2',DQSDT2,'VEGTYP', VEGTYP,&
699 'SOILTYP',SOILTYP, 'SLOPETYP',SLOPETYP, 'SHDFAC',SHDFAC,&
700 'SHMIN',SHMIN, 'ALBBRD',ALBBRD,'SNOALB1',SNOALB1,'TBOT',&
701 TBOT, 'Z0BRD',Z0BRD, 'Z0K',Z0K, 'CMC',CMC, 'T1',T1,'STC',&
702 STC, 'SMC',SMC, 'SWC',SWC,'SNOWHK',SNOWHK,'SNEQV',SNEQV,&
703 'ALBEDO',ALBEDO,'CHK',CHK,'ETA',ETA,'SHEAT',SHEAT, &
704 'ETA_KINEMATIC',ETA_KINEMATIC, 'FDOWN',FDOWN,'EC',EC, &
705 'EDIR',EDIR,'ET',ET,'ETT',ETT,'ESNOW',ESNOW,'DRIP',DRIP,&
706 'DEW',DEW,'BETA',BETA,'ETP',ETP,'SSOIL',SSOIL,'FLX1',FLX1,&
707 'FLX2',FLX2,'FLX3',FLX3,'SNOMLT',SNOMLT,'SNCOVR',SNCOVR,&
708 'RUNOFF1',RUNOFF1,'RUNOFF2',RUNOFF2,'RUNOFF3',RUNOFF3, &
709 'RC',RC, 'PC',PC,'RSMIN',RSMIN,'XLAI',XLAI,'RCS',RCS, &
710 'RCT',RCT,'RCQ',RCQ,'RCSOIL',RCSOIL,'SOILW',SOILW, &
711 'SOILM',SOILM,'Q1',Q1,'SMCWLT',SMCWLT,'SMCDRY',SMCDRY,&
712 'SMCREF',SMCREF,'SMCMAX',SMCMAX,'NROOT',NROOT
713 endif
714
715
716 CALL SFLX (FFROZP, ICE,DT,ZLVL,NSOIL,SLDPTH, & !C
717 LOCAL, & !L
718 LUTYPE, SLTYPE, & !CL
719 LWDN,SOLDN,SOLNET,SFCPRS,PRCP,SFCTMP,Q2K,DUMMY, & !F
720 DUMMY,DUMMY, DUMMY, & !F PRCPRAIN not used
721 TH2,Q2SAT,DQSDT2, & !I
722 VEGTYP,SOILTYP,SLOPETYP,SHDFAC,SHMIN,DUMMY, & !I
723 ALBBRD, SNOALB1,TBOT, Z0BRD, Z0K, EMISSI, & !S
724 CMC,T1,STC,SMC,SWC,SNOWHK,SNEQV,ALBEDOK,CHK,dummy,& !H
725 ETA,SHEAT, ETA_KINEMATIC,FDOWN, & !O
726 EC,EDIR,ET,ETT,ESNOW,DRIP,DEW, & !O
727 BETA,ETP,SSOIL, & !O
728 FLX1,FLX2,FLX3, & !O
729 SNOMLT,SNCOVR, & !O
730 RUNOFF1,RUNOFF2,RUNOFF3, & !O
731 RC,PC,RSMIN,XLAI,RCS,RCT,RCQ,RCSOIL, & !O
732 SOILW,SOILM,Q1, & !D
733 SMCWLT,SMCDRY,SMCREF,SMCMAX,NROOT)
734
735
736 IF(IPRINT) THEN
737
738 print*, 'AFTER SFLX, in Noahlsm_driver'
739 print*, 'ICE', ICE, 'DT',DT, 'ZLVL',ZLVL, 'NSOIL', NSOIL, &
740 'SLDPTH', SLDPTH, 'LOCAL',LOCAL, 'LUTYPE',&
741 LUTYPE, 'SLTYPE',SLTYPE, 'LWDN',LWDN, 'SOLDN',SOLDN, &
742 'SFCPRS',SFCPRS, 'PRCP',PRCP,'SFCTMP',SFCTMP,'Q2K',Q2K, &
743 'TH2',TH2,'Q2SAT',Q2SAT,'DQSDT2',DQSDT2,'VEGTYP', VEGTYP,&
744 'SOILTYP',SOILTYP, 'SLOPETYP',SLOPETYP, 'SHDFAC',SHDFAC,&
745 'SHDMIN',SHMIN, 'ALBBRD',ALBBRD,'SNOALB',SNOALB1,'TBOT',&
746 TBOT, 'Z0BRD',Z0BRD, 'Z0K',Z0K, 'CMC',CMC, 'T1',T1,'STC',&
747 STC, 'SMC',SMC, 'SWc',SWC,'SNOWHK',SNOWHK,'SNEQV',SNEQV,&
748 'ALBEDO',ALBEDOK,'CHK',CHK,'ETA',ETA,'SHEAT',SHEAT, &
749 'ETA_KINEMATIC',ETA_KINEMATIC, 'FDOWN',FDOWN,'EC',EC, &
750 'EDIR',EDIR,'ET',ET,'ETT',ETT,'ESNOW',ESNOW,'DRIP',DRIP,&
751 'DEW',DEW,'BETA',BETA,'ETP',ETP,'SSOIL',SSOIL,'FLX1',FLX1,&
752 'FLX2',FLX2,'FLX3',FLX3,'SNOMLT',SNOMLT,'SNCOVR',SNCOVR,&
753 'RUNOFF1',RUNOFF1,'RUNOFF2',RUNOFF2,'RUNOFF3',RUNOFF3, &
754 'RC',RC, 'PC',PC,'RSMIN',RSMIN,'XLAI',XLAI,'RCS',RCS, &
755 'RCT',RCT,'RCQ',RCQ,'RCSOIL',RCSOIL,'SOILW',SOILW, &
756 'SOILM',SOILM,'Q1',Q1,'SMCWLT',SMCWLT,'SMCDRY',SMCDRY,&
757 'SMCREF',SMCREF,'SMCMAX',SMCMAX,'NROOT',NROOT
758 endif
759
760 !*** UPDATE STATE VARIABLES
761 CANWAT(I,J)=CMC
762 SNOW(I,J)=SNEQV*1000.
763 ! SNOWH(I,J)=SNOWHK*1000.
764 SNOWH(I,J)=SNOWHK ! SNOWHK in meters
765 ALBEDO(I,J)=ALBEDOK
766 ! MEK Nov2006 turn off
767 ! ZNT(I,J)=Z0K
768 TSK(I,J)=T1
769 HFX(I,J)=SHEAT
770 ! MEk Jul07 add potential evap accum
771 POTEVP(I,J)=POTEVP(I,J)+ETP*FDTW
772 QFX(I,J)=ETA_KINEMATIC
773
774 LH(I,J)=ETA
775 GRDFLX(I,J)=SSOIL
776 SNOWC(I,J)=SNCOVR
777 CHS2(I,J)=CQS2(I,J)
778 ! prevent diagnostic ground q (q1) from being greater than qsat(tsk)
779 ! as happens over snow cover where the cqs2 value also becomes irrelevant
780 ! by setting cqs2=chs in this situation the 2m q should become just qv(k=1)
781 IF (Q1 .GT. QSFC(I,J)) THEN
782 CQS2(I,J) = CHS(I,J)
783 ENDIF
784 ! QSFC(I,J)=Q1
785 ! Convert QSFC back to mixing ratio
786 QSFC(I,J)= Q1/(1.0-Q1)
787 !
788 DO 80 NS=1,NSOIL
789 SMOIS(I,NS,J)=SMC(NS)
790 TSLB(I,NS,J)=STC(NS) ! STEMP
791 SH2O(I,NS,J)=SWC(NS)
792 80 CONTINUE
793 ! ENDIF
794
795 IF (UCMCALL == 1 ) THEN ! Beginning of UCM CALL if block
796 !--------------------------------------
797 ! URBAN CANOPY MODEL START - urban
798 !--------------------------------------
799 ! Input variables lsm --> urban
800
801
802 IF( IVGTYP(I,J) == 1 .or. IVGTYP(I,J) == 31 .or. &
803 IVGTYP(I,J) == 32 .or. IVGTYP(I,J) == 33 ) THEN
804
805 ! Call urban
806
807 !
808 UTYPE_URB = UTYPE_URB2D(I,J) !urban type (low, high or industrial)
809
810 TA_URB = SFCTMP ! [K]
811 QA_URB = Q2K ! [kg/kg]
812 UA_URB = SQRT(U_PHY(I,1,J)**2.+V_PHY(I,1,J)**2.)
813 U1_URB = U_PHY(I,1,J)
814 V1_URB = V_PHY(I,1,J)
815 IF(UA_URB < 1.) UA_URB=1. ! [m/s]
816 SSG_URB = SOLDN ! [W/m/m]
817 SSGD_URB = 0.8*SOLDN ! [W/m/m]
818 SSGQ_URB = SSG_URB-SSGD_URB ! [W/m/m]
819 LLG_URB = LWDN ! [W/m/m]
820 RAIN_URB = RAINBL(I,J) ! [mm]
821 RHOO_URB = SFCPRS / (287.04 * SFCTMP * (1.0+ 0.61 * Q2K)) ![kg/m/m/m]
822 ZA_URB = ZLVL ! [m]
823 DELT_URB = DT ! [sec]
824 XLAT_URB = XLAT_URB2D(I,J) ! [deg]
825 COSZ_URB = COSZ_URB2D(I,J) !
826 OMG_URB = OMG_URB2D(I,J) !
827 ZNT_URB = ZNT(I,J)
828
829 LSOLAR_URB = .FALSE.
830
831 TR_URB = TR_URB2D(I,J)
832 TB_URB = TB_URB2D(I,J)
833 TG_URB = TG_URB2D(I,J)
834 TC_URB = TC_URB2D(I,J)
835 QC_URB = QC_URB2D(I,J)
836 UC_URB = UC_URB2D(I,J)
837
838 DO K = 1,num_roof_layers
839 TRL_URB(K) = TRL_URB3D(I,K,J)
840 END DO
841 DO K = 1,num_wall_layers
842 TBL_URB(K) = TBL_URB3D(I,K,J)
843 END DO
844 DO K = 1,num_road_layers
845 TGL_URB(K) = TGL_URB3D(I,K,J)
846 END DO
847
848 XXXR_URB = XXXR_URB2D(I,J)
849 XXXB_URB = XXXB_URB2D(I,J)
850 XXXG_URB = XXXG_URB2D(I,J)
851 XXXC_URB = XXXC_URB2D(I,J)
852 !
853 CHS_URB = CHS(I,J)
854 CHS2_URB = CHS2(I,J)
855 !
856
857 ! Call urban
858
859
860 CALL urban(LSOLAR_URB, & ! I
861 num_roof_layers,num_wall_layers,num_road_layers, & ! C
862 DZR,DZB,DZG, & ! C
863 UTYPE_URB,TA_URB,QA_URB,UA_URB,U1_URB,V1_URB,SSG_URB, & ! I
864 SSGD_URB,SSGQ_URB,LLG_URB,RAIN_URB,RHOO_URB, & ! I
865 ZA_URB,DECLIN_URB,COSZ_URB,OMG_URB, & ! I
866 XLAT_URB,DELT_URB,ZNT_URB, & ! I
867 CHS_URB, CHS2_URB, & ! I
868 TR_URB, TB_URB, TG_URB, TC_URB, QC_URB,UC_URB, & ! H
869 TRL_URB,TBL_URB,TGL_URB, & ! H
870 XXXR_URB, XXXB_URB, XXXG_URB, XXXC_URB, & ! H
871 TS_URB,QS_URB,SH_URB,LH_URB,LH_KINEMATIC_URB, & ! O
872 SW_URB,ALB_URB,LW_URB,G_URB,RN_URB,PSIM_URB,PSIH_URB, & ! O
873 GZ1OZ0_URB, & !O
874 U10_URB, V10_URB, TH2_URB, Q2_URB, & ! O
875 UST_URB) !O
876
877
878 IF(IPRINT) THEN
879
880 print*, 'AFTER CALL URBAN'
881 print*,'num_roof_layers',num_roof_layers, 'num_wall_layers', &
882 num_wall_layers, &
883 'DZR',DZR,'DZB',DZB,'DZG',DZG,'UTYPE_URB',UTYPE_URB,'TA_URB', &
884 TA_URB, &
885 'QA_URB',QA_URB,'UA_URB',UA_URB,'U1_URB',U1_URB,'V1_URB', &
886 V1_URB, &
887 'SSG_URB',SSG_URB,'SSGD_URB',SSGD_URB,'SSGQ_URB',SSGQ_URB, &
888 'LLG_URB',LLG_URB,'RAIN_URB',RAIN_URB,'RHOO_URB',RHOO_URB, &
889 'ZA_URB',ZA_URB, 'DECLIN_URB',DECLIN_URB,'COSZ_URB',COSZ_URB,&
890 'OMG_URB',OMG_URB,'XLAT_URB',XLAT_URB,'DELT_URB',DELT_URB, &
891 'ZNT_URB',ZNT_URB,'TR_URB',TR_URB, 'TB_URB',TB_URB,'TG_URB',&
892 TG_URB,'TC_URB',TC_URB,'QC_URB',QC_URB,'TRL_URB',TRL_URB, &
893 'TBL_URB',TBL_URB,'TGL_URB',TGL_URB,'XXXR_URB',XXXR_URB, &
894 'XXXB_URB',XXXB_URB,'XXXG_URB',XXXG_URB,'XXXC_URB',XXXC_URB,&
895 'TS_URB',TS_URB,'QS_URB',QS_URB,'SH_URB',SH_URB,'LH_URB', &
896 LH_URB, 'LH_KINEMATIC_URB',LH_KINEMATIC_URB,'SW_URB',SW_URB,&
897 'ALB_URB',ALB_URB,'LW_URB',LW_URB,'G_URB',G_URB,'RN_URB', &
898 RN_URB, 'PSIM_URB',PSIM_URB,'PSIH_URB',PSIH_URB, &
899 'U10_URB',U10_URB,'V10_URB',V10_URB,'TH2_URB',TH2_URB, &
900 'Q2_URB',Q2_URB,'CHS_URB',CHS_URB,'CHS2_URB',CHS2_URB
901 endif
902
903 TS_URB2D(I,J) = TS_URB
904
905 ALBEDO(I,J) = FRC_URB2D(I,J)*ALB_URB+(1-FRC_URB2D(I,J))*ALBEDOK ![-]
906 HFX(I,J) = FRC_URB2D(I,J)*SH_URB+(1-FRC_URB2D(I,J))*SHEAT ![W/m/m]
907 QFX(I,J) = FRC_URB2D(I,J)*LH_KINEMATIC_URB &
908 + (1-FRC_URB2D(I,J))*ETA_KINEMATIC ![kg/m/m/s]
909 LH(I,J) = FRC_URB2D(I,J)*LH_URB+(1-FRC_URB2D(I,J))*ETA ![W/m/m]
910 GRDFLX(I,J) = FRC_URB2D(I,J)*G_URB+(1-FRC_URB2D(I,J))*SSOIL ![W/m/m]
911 TSK(I,J) = FRC_URB2D(I,J)*TS_URB+(1-FRC_URB2D(I,J))*T1 ![K]
912 QSFC(I,J)= FRC_URB2D(I,J)*QS_URB+(1-FRC_URB2D(I,J))*Q1 ![-]
913
914 IF(IPRINT)THEN
915
916 print*, ' FRC_URB2D', FRC_URB2D, &
917 'ALB_URB',ALB_URB, 'ALBEDOK',ALBEDOK, &
918 'ALBEDO(I,J)', ALBEDO(I,J), &
919 'SH_URB',SH_URB,'SHEAT',SHEAT, 'HFX(I,J)',HFX(I,J), &
920 'LH_KINEMATIC_URB',LH_KINEMATIC_URB,'ETA_KINEMATIC', &
921 ETA_KINEMATIC, 'QFX(I,J)',QFX(I,J), &
922 'LH_URB',LH_URB, 'ETA',ETA, 'LH(I,J)',LH(I,J), &
923 'G_URB',G_URB,'SSOIL',SSOIL,'GRDFLX(I,J)', GRDFLX(I,J),&
924 'TS_URB',TS_URB,'T1',T1,'TSK(I,J)',TSK(I,J), &
925 'QS_URB',QS_URB,'Q1',Q1,'QSFC(I,J)',QSFC(I,J)
926 endif
927
928
929
930
931 ! Renew Urban State Varialbes
932
933 TR_URB2D(I,J) = TR_URB
934 TB_URB2D(I,J) = TB_URB
935 TG_URB2D(I,J) = TG_URB
936 TC_URB2D(I,J) = TC_URB
937 QC_URB2D(I,J) = QC_URB
938 UC_URB2D(I,J) = UC_URB
939
940 DO K = 1,num_roof_layers
941 TRL_URB3D(I,K,J) = TRL_URB(K)
942 END DO
943 DO K = 1,num_wall_layers
944 TBL_URB3D(I,K,J) = TBL_URB(K)
945 END DO
946 DO K = 1,num_road_layers
947 TGL_URB3D(I,K,J) = TGL_URB(K)
948 END DO
949 XXXR_URB2D(I,J) = XXXR_URB
950 XXXB_URB2D(I,J) = XXXB_URB
951 XXXG_URB2D(I,J) = XXXG_URB
952 XXXC_URB2D(I,J) = XXXC_URB
953
954 SH_URB2D(I,J) = SH_URB
955 LH_URB2D(I,J) = LH_URB
956 G_URB2D(I,J) = G_URB
957 RN_URB2D(I,J) = RN_URB
958 PSIM_URB2D(I,J) = PSIM_URB
959 PSIH_URB2D(I,J) = PSIH_URB
960 GZ1OZ0_URB2D(I,J)= GZ1OZ0_URB
961 U10_URB2D(I,J) = U10_URB
962 V10_URB2D(I,J) = V10_URB
963 TH2_URB2D(I,J) = TH2_URB
964 Q2_URB2D(I,J) = Q2_URB
965 UST_URB2D(I,J) = UST_URB
966 AKMS_URB2D(I,J) = KARMAN * UST_URB2D(I,J)/(GZ1OZ0_URB2D(I,J)-PSIM_URB2D(I,J))
967
968 END IF
969
970 ENDIF ! end of UCM CALL if block
971 !--------------------------------------
972 ! Urban Part End - urban
973 !--------------------------------------
974
975 !*** DIAGNOSTICS
976 SMSTAV(I,J)=SOILW
977 SMSTOT(I,J)=SOILM*1000.
978 ! Convert the water unit into mm
979 SFCRUNOFF(I,J)=SFCRUNOFF(I,J)+RUNOFF1*DT*1000.0
980 UDRUNOFF(I,J)=UDRUNOFF(I,J)+(RUNOFF2+RUNOFF3)*DT*1000.0
981 ! snow defined when fraction of frozen precip (FFROZP) > 0.5,
982 IF(FFROZP.GT.0.5)THEN
983 ACSNOW(I,J)=ACSNOW(I,J)+PRCP*DT
984 ENDIF
985 IF(SNOW(I,J).GT.0.)THEN
986 ACSNOM(I,J)=ACSNOM(I,J)+SNOMLT*1000.
987 ! accumulated snow-melt energy
988 SNOPCX(I,J)=SNOPCX(I,J)-SNOMLT/FDTLIW
989 ENDIF
990
991 ENDIF ! endif of land-sea test
992
993 100 CONTINUE ! of I loop
994
995 ENDDO ! of J loop
996 !------------------------------------------------------
997 END SUBROUTINE lsm
998 !------------------------------------------------------
999
1000 SUBROUTINE LSMINIT(VEGFRA,SNOW,SNOWC,SNOWH,CANWAT,SMSTAV, &
1001 SMSTOT, SFCRUNOFF,UDRUNOFF,ACSNOW, &
1002 ACSNOM,IVGTYP,ISLTYP,TSLB,SMOIS,SH2O,ZS,DZS, &
1003 FNDSOILW, FNDSNOWH, &
1004 num_soil_layers, restart, &
1005 allowed_to_read , &
1006 ids,ide, jds,jde, kds,kde, &
1007 ims,ime, jms,jme, kms,kme, &
1008 its,ite, jts,jte, kts,kte )
1009
1010 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
1011 ims,ime, jms,jme, kms,kme, &
1012 its,ite, jts,jte, kts,kte
1013
1014 INTEGER, INTENT(IN) :: num_soil_layers
1015
1016 LOGICAL , INTENT(IN) :: restart , allowed_to_read
1017
1018 REAL, DIMENSION( num_soil_layers), INTENT(INOUT) :: ZS, DZS
1019
1020 REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , &
1021 INTENT(INOUT) :: SMOIS, & !Total soil moisture
1022 SH2O, & !liquid soil moisture
1023 TSLB !STEMP
1024
1025 REAL, DIMENSION( ims:ime, jms:jme ) , &
1026 INTENT(INOUT) :: SNOW, &
1027 SNOWH, &
1028 SNOWC, &
1029 CANWAT, &
1030 SMSTAV, &
1031 SMSTOT, &
1032 SFCRUNOFF, &
1033 UDRUNOFF, &
1034 ACSNOW, &
1035 VEGFRA, &
1036 ACSNOM
1037
1038 INTEGER, DIMENSION( ims:ime, jms:jme ) , &
1039 INTENT(IN) :: IVGTYP, &
1040 ISLTYP
1041
1042 LOGICAL, INTENT(IN) :: FNDSOILW , &
1043 FNDSNOWH
1044
1045 INTEGER :: L
1046 REAL :: BX, SMCMAX, PSISAT, FREE
1047 REAL, PARAMETER :: BLIM = 5.5, HLICE = 3.335E5, &
1048 GRAV = 9.81, T0 = 273.15
1049 INTEGER :: errflag
1050
1051 !
1052
1053
1054 ! initialize three Noah LSM related tables
1055 IF ( allowed_to_read ) THEN
1056 CALL wrf_message( 'INITIALIZE THREE Noah LSM RELATED TABLES' )
1057 CALL LSM_PARM_INIT
1058 ENDIF
1059
1060 IF(.not.restart)THEN
1061
1062 itf=min0(ite,ide-1)
1063 jtf=min0(jte,jde-1)
1064
1065 errflag = 0
1066 DO j = jts,jtf
1067 DO i = its,itf
1068 IF ( ISLTYP( i,j ) .LT. 1 ) THEN
1069 errflag = 1
1070 WRITE(err_message,*)"module_sf_noahlsm.F: lsminit: out of range ISLTYP ",i,j,ISLTYP( i,j )
1071 CALL wrf_message(err_message)
1072 ENDIF
1073 ENDDO
1074 ENDDO
1075 IF ( errflag .EQ. 1 ) THEN
1076 CALL wrf_error_fatal( "module_sf_noahlsm.F: lsminit: out of range value "// &
1077 "of ISLTYP. Is this field in the input?" )
1078 ENDIF
1079
1080 ! initialize soil liquid water content SH2O
1081
1082 ! IF(.NOT.FNDSOILW) THEN
1083
1084 ! If no SWC, do the following
1085 ! PRINT *,'SOIL WATER NOT FOUND - VALUE SET IN LSMINIT'
1086 DO J = jts,jtf
1087 DO I = its,itf
1088 BX = BB(ISLTYP(I,J))
1089 SMCMAX = MAXSMC(ISLTYP(I,J))
1090 PSISAT = SATPSI(ISLTYP(I,J))
1091 if ((bx > 0.0).and.(smcmax > 0.0).and.(psisat > 0.0)) then
1092 DO NS=1, num_soil_layers
1093 ! ----------------------------------------------------------------------
1094 !SH2O <= SMOIS for T < 273.149K (-0.001C)
1095 IF (TSLB(I,NS,J) < 273.149) THEN
1096 ! ----------------------------------------------------------------------
1097 ! first guess following explicit solution for Flerchinger Eqn from Koren
1098 ! et al, JGR, 1999, Eqn 17 (KCOUNT=0 in FUNCTION FRH2O).
1099 ! ISLTPK is soil type
1100 BX = BB(ISLTYP(I,J))
1101 SMCMAX = MAXSMC(ISLTYP(I,J))
1102 PSISAT = SATPSI(ISLTYP(I,J))
1103 IF ( BX > BLIM ) BX = BLIM
1104 FK=(( (HLICE/(GRAV*(-PSISAT))) * &
1105 ((TSLB(I,NS,J)-T0)/TSLB(I,NS,J)) )**(-1/BX) )*SMCMAX
1106 IF (FK < 0.02) FK = 0.02
1107 SH2O(I,NS,J) = MIN( FK, SMOIS(I,NS,J) )
1108 ! ----------------------------------------------------------------------
1109 ! now use iterative solution for liquid soil water content using
1110 ! FUNCTION FRH2O with the initial guess for SH2O from above explicit
1111 ! first guess.
1112 CALL FRH2O (FREE,TSLB(I,NS,J),SMOIS(I,NS,J),SH2O(I,NS,J), &
1113 SMCMAX,BX,PSISAT)
1114 SH2O(I,NS,J) = FREE
1115 ELSE ! of IF (TSLB(I,NS,J)
1116 ! ----------------------------------------------------------------------
1117 ! SH2O = SMOIS ( for T => 273.149K (-0.001C)
1118 SH2O(I,NS,J)=SMOIS(I,NS,J)
1119 ! ----------------------------------------------------------------------
1120 ENDIF ! of IF (TSLB(I,NS,J)
1121 END DO ! of DO NS=1, num_soil_layers
1122 else ! of if ((bx > 0.0)
1123 DO NS=1, num_soil_layers
1124 SH2O(I,NS,J)=SMOIS(I,NS,J)
1125 END DO
1126 endif ! of if ((bx > 0.0)
1127 ENDDO ! DO I = its,itf
1128 ENDDO ! DO J = jts,jtf
1129 ! ENDIF ! of IF(.NOT.FNDSOILW)THEN
1130
1131 ! initialize physical snow height SNOWH
1132
1133 IF(.NOT.FNDSNOWH)THEN
1134 ! If no SNOWH do the following
1135 CALL wrf_message( 'SNOW HEIGHT NOT FOUND - VALUE DEFINED IN LSMINIT' )
1136 DO J = jts,jtf
1137 DO I = its,itf
1138 SNOWH(I,J)=SNOW(I,J)*0.005 ! SNOW in mm and SNOWH in m
1139 ENDDO
1140 ENDDO
1141 ENDIF
1142
1143 ! initialize canopy water to ZERO
1144
1145 ! GO TO 110
1146 ! print*,'Note that canopy water content (CANWAT) is set to ZERO in LSMINIT'
1147 DO J = jts,jtf
1148 DO I = its,itf
1149 CANWAT(I,J)=0.0
1150 ENDDO
1151 ENDDO
1152 110 CONTINUE
1153
1154 ENDIF
1155 !------------------------------------------------------------------------------
1156 END SUBROUTINE lsminit
1157 !------------------------------------------------------------------------------
1158
1159
1160
1161 !
1162 !-----------------------------------------------------------------
1163 SUBROUTINE LSM_PARM_INIT
1164 !-----------------------------------------------------------------
1165
1166 character*4 :: MMINLU, MMINSL
1167
1168 MMINLU='USGS'
1169 MMINSL='STAS'
1170 call SOIL_VEG_GEN_PARM( MMINLU, MMINSL)
1171
1172 !-----------------------------------------------------------------
1173 END SUBROUTINE LSM_PARM_INIT
1174 !-----------------------------------------------------------------
1175
1176 !-----------------------------------------------------------------
1177 SUBROUTINE SOIL_VEG_GEN_PARM( MMINLU, MMINSL)
1178 !-----------------------------------------------------------------
1179
1180 IMPLICIT NONE
1181
1182 integer :: LUMATCH, IINDEX, LC, NUM_SLOPE
1183 integer :: ierr
1184 INTEGER , PARAMETER :: OPEN_OK = 0
1185
1186 character*4 :: MMINLU, MMINSL
1187 character*128 :: mess , message
1188 logical, external :: wrf_dm_on_monitor
1189
1190
1191 !-----SPECIFY VEGETATION RELATED CHARACTERISTICS :
1192 ! ALBBCK: SFC albedo (in percentage)
1193 ! Z0: Roughness length (m)
1194 ! SHDFAC: Green vegetation fraction (in percentage)
1195 ! Note: The ALBEDO, Z0, and SHDFAC values read from the following table
1196 ! ALBEDO, amd Z0 are specified in LAND-USE TABLE; and SHDFAC is
1197 ! the monthly green vegetation data
1198 ! CMXTBL: MAX CNPY Capacity (m)
1199 ! NROTBL: Rooting depth (layer)
1200 ! RSMIN: Mimimum stomatal resistance (s m-1)
1201 ! RSMAX: Max. stomatal resistance (s m-1)
1202 ! RGL: Parameters used in radiation stress function
1203 ! HS: Parameter used in vapor pressure deficit functio
1204 ! TOPT: Optimum transpiration air temperature. (K)
1205 ! CMCMAX: Maximum canopy water capacity
1206 ! CFACTR: Parameter used in the canopy inteception calculati
1207 ! SNUP: Threshold snow depth (in water equivalent m) that
1208 ! implies 100% snow cover
1209 ! LAI: Leaf area index (dimensionless)
1210 ! MAXALB: Upper bound on maximum albedo over deep snow
1211 !
1212 !-----READ IN VEGETAION PROPERTIES FROM VEGPARM.TBL
1213 !
1214
1215 IF ( wrf_dm_on_monitor() ) THEN
1216
1217 OPEN(19, FILE='VEGPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
1218 IF(ierr .NE. OPEN_OK ) THEN
1219 WRITE(message,FMT='(A)') &
1220 'module_sf_noahlsm.F: soil_veg_gen_parm: failure opening VEGPARM.TBL'
1221 CALL wrf_error_fatal ( message )
1222 END IF
1223
1224 WRITE ( mess, * ) 'INPUT LANDUSE = ',MMINLU
1225 CALL wrf_message( mess )
1226
1227 LUMATCH=0
1228
1229 READ (19,*)
1230 READ (19,2000,END=2002)LUTYPE
1231 READ (19,*)LUCATS,IINDEX
1232 2000 FORMAT (A4)
1233
1234 IF(LUTYPE.EQ.MMINLU)THEN
1235 WRITE( mess , * ) 'LANDUSE TYPE = ',LUTYPE,' FOUND', &
1236 LUCATS,' CATEGORIES'
1237 CALL wrf_message( mess )
1238 LUMATCH=1
1239 ENDIF
1240
1241 IF(LUTYPE.EQ.MMINLU)THEN
1242 DO LC=1,LUCATS
1243 READ (19,*)IINDEX,ALBTBL(LC),Z0TBL(LC),SHDTBL(LC), &
1244 NROTBL(LC),RSTBL(LC),RGLTBL(LC),HSTBL(LC), &
1245 SNUPTBL(LC),LAITBL(LC),MAXALB(LC)
1246 ENDDO
1247 !
1248 READ (19,*)
1249 READ (19,*)TOPT_DATA
1250 READ (19,*)
1251 READ (19,*)CMCMAX_DATA
1252 READ (19,*)
1253 READ (19,*)CFACTR_DATA
1254 READ (19,*)
1255 READ (19,*)RSMAX_DATA
1256 READ (19,*)
1257 READ (19,*)BARE
1258 ENDIF
1259 !
1260 2002 CONTINUE
1261
1262 CLOSE (19)
1263 ENDIF
1264
1265 CALL wrf_dm_bcast_string ( LUTYPE , 4 )
1266 CALL wrf_dm_bcast_integer ( LUCATS , 1 )
1267 CALL wrf_dm_bcast_integer ( IINDEX , 1 )
1268 CALL wrf_dm_bcast_integer ( LUMATCH , 1 )
1269 CALL wrf_dm_bcast_real ( ALBTBL , NLUS )
1270 CALL wrf_dm_bcast_real ( Z0TBL , NLUS )
1271 CALL wrf_dm_bcast_real ( SHDTBL , NLUS )
1272 CALL wrf_dm_bcast_real ( NROTBL , NLUS )
1273 CALL wrf_dm_bcast_real ( RSTBL , NLUS )
1274 CALL wrf_dm_bcast_real ( RGLTBL , NLUS )
1275 CALL wrf_dm_bcast_real ( HSTBL , NLUS )
1276 CALL wrf_dm_bcast_real ( SNUPTBL , NLUS )
1277 CALL wrf_dm_bcast_real ( LAITBL , NLUS )
1278 CALL wrf_dm_bcast_real ( MAXALB , NLUS )
1279 CALL wrf_dm_bcast_real ( TOPT_DATA , 1 )
1280 CALL wrf_dm_bcast_real ( CMCMAX_DATA , 1 )
1281 CALL wrf_dm_bcast_real ( CFACTR_DATA , 1 )
1282 CALL wrf_dm_bcast_real ( RSMAX_DATA , 1 )
1283 CALL wrf_dm_bcast_integer ( BARE , 1 )
1284
1285 !
1286 !-----READ IN SOIL PROPERTIES FROM SOILPARM.TBL
1287 !
1288 IF ( wrf_dm_on_monitor() ) THEN
1289 OPEN(19, FILE='SOILPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
1290 IF(ierr .NE. OPEN_OK ) THEN
1291 WRITE(message,FMT='(A)') &
1292 'module_sf_noahlsm.F: soil_veg_gen_parm: failure opening SOILPARM.TBL'
1293 CALL wrf_error_fatal ( message )
1294 END IF
1295
1296 MMINSL='STAS' !oct2
1297 WRITE(mess,*) 'INPUT SOIL TEXTURE CLASSIFICAION = ',MMINSL
1298 CALL wrf_message( mess )
1299
1300 LUMATCH=0
1301
1302 READ (19,*)
1303 READ (19,2000,END=2003)SLTYPE
1304 READ (19,*)SLCATS,IINDEX
1305 IF(SLTYPE.EQ.MMINSL)THEN
1306 WRITE( mess , * ) 'SOIL TEXTURE CLASSIFICATION = ',SLTYPE,' FOUND', &
1307 SLCATS,' CATEGORIES'
1308 CALL wrf_message ( mess )
1309 LUMATCH=1
1310 ENDIF
1311 IF(SLTYPE.EQ.MMINSL)THEN
1312 DO LC=1,SLCATS
1313 READ (19,*) IINDEX,BB(LC),DRYSMC(LC),F11(LC),MAXSMC(LC),&
1314 REFSMC(LC),SATPSI(LC),SATDK(LC), SATDW(LC), &
1315 WLTSMC(LC), QTZ(LC)
1316 ENDDO
1317 ENDIF
1318
1319 2003 CONTINUE
1320
1321 CLOSE (19)
1322 ENDIF
1323
1324 CALL wrf_dm_bcast_integer ( LUMATCH , 1 )
1325 CALL wrf_dm_bcast_string ( SLTYPE , 4 )
1326 CALL wrf_dm_bcast_string ( MMINSL , 4 ) ! since this is reset above, see oct2 ^
1327 CALL wrf_dm_bcast_integer ( SLCATS , 1 )
1328 CALL wrf_dm_bcast_integer ( IINDEX , 1 )
1329 CALL wrf_dm_bcast_real ( BB , NSLTYPE )
1330 CALL wrf_dm_bcast_real ( DRYSMC , NSLTYPE )
1331 CALL wrf_dm_bcast_real ( F11 , NSLTYPE )
1332 CALL wrf_dm_bcast_real ( MAXSMC , NSLTYPE )
1333 CALL wrf_dm_bcast_real ( REFSMC , NSLTYPE )
1334 CALL wrf_dm_bcast_real ( SATPSI , NSLTYPE )
1335 CALL wrf_dm_bcast_real ( SATDK , NSLTYPE )
1336 CALL wrf_dm_bcast_real ( SATDW , NSLTYPE )
1337 CALL wrf_dm_bcast_real ( WLTSMC , NSLTYPE )
1338 CALL wrf_dm_bcast_real ( QTZ , NSLTYPE )
1339
1340 IF(LUMATCH.EQ.0)THEN
1341 CALL wrf_message( 'SOIl TEXTURE IN INPUT FILE DOES NOT ' )
1342 CALL wrf_message( 'MATCH SOILPARM TABLE' )
1343 CALL wrf_error_fatal ( 'INCONSISTENT OR MISSING SOILPARM FILE' )
1344 ENDIF
1345
1346 !
1347 !-----READ IN GENERAL PARAMETERS FROM GENPARM.TBL
1348 !
1349 IF ( wrf_dm_on_monitor() ) THEN
1350 OPEN(19, FILE='GENPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
1351 IF(ierr .NE. OPEN_OK ) THEN
1352 WRITE(message,FMT='(A)') &
1353 'module_sf_noahlsm.F: soil_veg_gen_parm: failure opening GENPARM.TBL'
1354 CALL wrf_error_fatal ( message )
1355 END IF
1356
1357 READ (19,*)
1358 READ (19,*)
1359 READ (19,*) NUM_SLOPE
1360
1361 SLPCATS=NUM_SLOPE
1362
1363 DO LC=1,SLPCATS
1364 READ (19,*)SLOPE_DATA(LC)
1365 ENDDO
1366
1367 READ (19,*)
1368 READ (19,*)SBETA_DATA
1369 READ (19,*)
1370 READ (19,*)FXEXP_DATA
1371 READ (19,*)
1372 READ (19,*)CSOIL_DATA
1373 READ (19,*)
1374 READ (19,*)SALP_DATA
1375 READ (19,*)
1376 READ (19,*)REFDK_DATA
1377 READ (19,*)
1378 READ (19,*)REFKDT_DATA
1379 READ (19,*)
1380 READ (19,*)FRZK_DATA
1381 READ (19,*)
1382 READ (19,*)ZBOT_DATA
1383 READ (19,*)
1384 READ (19,*)CZIL_DATA
1385 READ (19,*)
1386 READ (19,*)SMLOW_DATA
1387 READ (19,*)
1388 READ (19,*)SMHIGH_DATA
1389 CLOSE (19)
1390 ENDIF
1391
1392 CALL wrf_dm_bcast_integer ( NUM_SLOPE , 1 )
1393 CALL wrf_dm_bcast_integer ( SLPCATS , 1 )
1394 CALL wrf_dm_bcast_real ( SLOPE_DATA , NSLOPE )
1395 CALL wrf_dm_bcast_real ( SBETA_DATA , 1 )
1396 CALL wrf_dm_bcast_real ( FXEXP_DATA , 1 )
1397 CALL wrf_dm_bcast_real ( CSOIL_DATA , 1 )
1398 CALL wrf_dm_bcast_real ( SALP_DATA , 1 )
1399 CALL wrf_dm_bcast_real ( REFDK_DATA , 1 )
1400 CALL wrf_dm_bcast_real ( REFKDT_DATA , 1 )
1401 CALL wrf_dm_bcast_real ( FRZK_DATA , 1 )
1402 CALL wrf_dm_bcast_real ( ZBOT_DATA , 1 )
1403 CALL wrf_dm_bcast_real ( CZIL_DATA , 1 )
1404 CALL wrf_dm_bcast_real ( SMLOW_DATA , 1 )
1405 CALL wrf_dm_bcast_real ( SMHIGH_DATA , 1 )
1406
1407
1408 !-----------------------------------------------------------------
1409 END SUBROUTINE SOIL_VEG_GEN_PARM
1410 !-----------------------------------------------------------------
1411
1412 SUBROUTINE SFLX (FFROZP,ICE,DT,ZLVL,NSOIL,SLDPTH, & !C
1413 LOCAL, & !L
1414 LLANDUSE, LSOIL, & !CL
1415 LWDN,SOLDN,SOLNET,SFCPRS,PRCP,SFCTMP,Q2,SFCSPD, & !F
1416 COSZ,PRCPRAIN, SOLARDIRECT, & !F
1417 TH2,Q2SAT,DQSDT2, & !I
1418 VEGTYP,SOILTYP,SLOPETYP,SHDFAC,SHDMIN,SHDMAX, & !I
1419 ALB, SNOALB,TBOT, Z0BRD, Z0, EMISSI, & !S
1420 CMC,T1,STC,SMC,SH2O,SNOWH,SNEQV,ALBEDO,CH,CM, & !H
1421 ! ----------------------------------------------------------------------
1422 ! OUTPUTS, DIAGNOSTICS, PARAMETERS BELOW GENERALLY NOT NECESSARY WHEN
1423 ! COUPLED WITH E.G. A NWP MODEL (SUCH AS THE NOAA/NWS/NCEP MESOSCALE ETA
1424 ! MODEL). OTHER APPLICATIONS MAY REQUIRE DIFFERENT OUTPUT VARIABLES.
1425 ! ----------------------------------------------------------------------
1426 ETA,SHEAT, ETA_KINEMATIC,FDOWN, & !O
1427 EC,EDIR,ET,ETT,ESNOW,DRIP,DEW, & !O
1428 BETA,ETP,SSOIL, & !O
1429 FLX1,FLX2,FLX3, & !O
1430 SNOMLT,SNCOVR, & !O
1431 RUNOFF1,RUNOFF2,RUNOFF3, & !O
1432 RC,PC,RSMIN,XLAI,RCS,RCT,RCQ,RCSOIL, & !O
1433 SOILW,SOILM,Q1, & !D
1434 SMCWLT,SMCDRY,SMCREF,SMCMAX,NROOT) !P
1435 ! ----------------------------------------------------------------------
1436 ! SUBROUTINE SFLX - UNIFIED NOAHLSM VERSION 1.0 JULY 2007
1437 ! ----------------------------------------------------------------------
1438 ! SUB-DRIVER FOR "Noah LSM" FAMILY OF PHYSICS SUBROUTINES FOR A
1439 ! SOIL/VEG/SNOWPACK LAND-SURFACE MODEL TO UPDATE SOIL MOISTURE, SOIL
1440 ! ICE, SOIL TEMPERATURE, SKIN TEMPERATURE, SNOWPACK WATER CONTENT,
1441 ! SNOWDEPTH, AND ALL TERMS OF THE SURFACE ENERGY BALANCE AND SURFACE
1442 ! WATER BALANCE (EXCLUDING INPUT ATMOSPHERIC FORCINGS OF DOWNWARD
1443 ! RADIATION AND PRECIP)
1444 ! ----------------------------------------------------------------------
1445 ! SFLX ARGUMENT LIST KEY:
1446 ! ----------------------------------------------------------------------
1447 ! C CONFIGURATION INFORMATION
1448 ! L LOGICAL
1449 ! CL 4-string character bearing logical meaning
1450 ! F FORCING DATA
1451 ! I OTHER (INPUT) FORCING DATA
1452 ! S SURFACE CHARACTERISTICS
1453 ! H HISTORY (STATE) VARIABLES
1454 ! O OUTPUT VARIABLES
1455 ! D DIAGNOSTIC OUTPUT
1456 ! P Parameters
1457 ! Msic Miscellaneous terms passed from gridded driver
1458 ! ----------------------------------------------------------------------
1459 ! 1. CONFIGURATION INFORMATION (C):
1460 ! ----------------------------------------------------------------------
1461 ! ICE SEA-ICE FLAG (=1: SEA-ICE, =0: LAND)
1462 ! DT TIMESTEP (SEC) (DT SHOULD NOT EXCEED 3600 SECS, RECOMMEND
1463 ! 1800 SECS OR LESS)
1464 ! ZLVL HEIGHT (M) ABOVE GROUND OF ATMOSPHERIC FORCING VARIABLES
1465 ! NSOIL NUMBER OF SOIL LAYERS (AT LEAST 2, AND NOT GREATER THAN
1466 ! PARAMETER NSOLD SET BELOW)
1467 ! SLDPTH THE THICKNESS OF EACH SOIL LAYER (M)
1468 ! ----------------------------------------------------------------------
1469 ! 2. LOGICAL:
1470 ! ----------------------------------------------------------------------
1471 ! LCH Exchange coefficient (Ch) calculation flag (false: using
1472 ! ch-routine SFCDIF; true: Ch is brought in)
1473 ! LOCAL Flag for local-site simulation (where there is no
1474 ! maps for albedo, veg fraction, and roughness
1475 ! true: all LSM parameters (inluding albedo, veg fraction and
1476 ! roughness length) will be defined by three tables
1477 ! LLANDUSE (=USGS, using USGS landuse classification)
1478 ! LSOIL (=STAS, using FAO/STATSGO soil texture classification)
1479 ! ----------------------------------------------------------------------
1480 ! 3. FORCING DATA (F):
1481 ! ----------------------------------------------------------------------
1482 ! LWDN LW DOWNWARD RADIATION (W M-2; POSITIVE, NOT NET LONGWAVE)
1483 ! SOLDN SOLAR DOWNWARD RADIATION (W M-2; POSITIVE, NOT NET SOLAR)
1484 ! SFCPRS PRESSURE AT HEIGHT ZLVL ABOVE GROUND (PASCALS)
1485 ! PRCP PRECIP RATE (KG M-2 S-1) (NOTE, THIS IS A RATE)
1486 ! SFCTMP AIR TEMPERATURE (K) AT HEIGHT ZLVL ABOVE GROUND
1487 ! TH2 AIR POTENTIAL TEMPERATURE (K) AT HEIGHT ZLVL ABOVE GROUND
1488 ! Q2 MIXING RATIO AT HEIGHT ZLVL ABOVE GROUND (KG KG-1)
1489 ! COSZ Solar zenith angle (not used for now)
1490 ! PRCPRAIN Liquid-precipitation rate (KG M-2 S-1) (not used)
1491 ! SOLARDIRECT Direct component of downward solar radiation (W M-2) (not used)
1492 ! ----------------------------------------------------------------------
1493 ! 4. OTHER FORCING (INPUT) DATA (I):
1494 ! ----------------------------------------------------------------------
1495 ! SFCSPD WIND SPEED (M S-1) AT HEIGHT ZLVL ABOVE GROUND
1496 ! Q2SAT SAT SPECIFIC HUMIDITY AT HEIGHT ZLVL ABOVE GROUND (KG KG-1)
1497 ! DQSDT2 SLOPE OF SAT SPECIFIC HUMIDITY CURVE AT T=SFCTMP
1498 ! (KG KG-1 K-1)
1499 ! ----------------------------------------------------------------------
1500 ! 5. CANOPY/SOIL CHARACTERISTICS (S):
1501 ! ----------------------------------------------------------------------
1502 ! VEGTYP VEGETATION TYPE (INTEGER INDEX)
1503 ! SOILTYP SOIL TYPE (INTEGER INDEX)
1504 ! SLOPETYP CLASS OF SFC SLOPE (INTEGER INDEX)
1505 ! SHDFAC AREAL FRACTIONAL COVERAGE OF GREEN VEGETATION
1506 ! (FRACTION= 0.0-1.0)
1507 ! SHDMIN MINIMUM AREAL FRACTIONAL COVERAGE OF GREEN VEGETATION
1508 ! (FRACTION= 0.0-1.0) <= SHDFAC
1509 ! PTU PHOTO THERMAL UNIT (PLANT PHENOLOGY FOR ANNUALS/CROPS)
1510 ! (NOT YET USED, BUT PASSED TO REDPRM FOR FUTURE USE IN
1511 ! VEG PARMS)
1512 ! ALB BACKROUND SNOW-FREE SURFACE ALBEDO (FRACTION), FOR JULIAN
1513 ! DAY OF YEAR (USUALLY FROM TEMPORAL INTERPOLATION OF
1514 ! MONTHLY MEAN VALUES' CALLING PROG MAY OR MAY NOT
1515 ! INCLUDE DIURNAL SUN ANGLE EFFECT)
1516 ! SNOALB UPPER BOUND ON MAXIMUM ALBEDO OVER DEEP SNOW (E.G. FROM
1517 ! ROBINSON AND KUKLA, 1985, J. CLIM. & APPL. METEOR.)
1518 ! TBOT BOTTOM SOIL TEMPERATURE (LOCAL YEARLY-MEAN SFC AIR
1519 ! TEMPERATURE)
1520 ! Z0BRD Background fixed roughness length (M)
1521 ! Z0 Time varying roughness length (M) as function of snow depth
1522 ! ----------------------------------------------------------------------
1523 ! 6. HISTORY (STATE) VARIABLES (H):
1524 ! ----------------------------------------------------------------------
1525 ! CMC CANOPY MOISTURE CONTENT (M)
1526 ! T1 GROUND/CANOPY/SNOWPACK) EFFECTIVE SKIN TEMPERATURE (K)
1527 ! STC(NSOIL) SOIL TEMP (K)
1528 ! SMC(NSOIL) TOTAL SOIL MOISTURE CONTENT (VOLUMETRIC FRACTION)
1529 ! SH2O(NSOIL) UNFROZEN SOIL MOISTURE CONTENT (VOLUMETRIC FRACTION)
1530 ! NOTE: FROZEN SOIL MOISTURE = SMC - SH2O
1531 ! SNOWH ACTUAL SNOW DEPTH (M)
1532 ! SNEQV LIQUID WATER-EQUIVALENT SNOW DEPTH (M)
1533 ! NOTE: SNOW DENSITY = SNEQV/SNOWH
1534 ! ALBEDO SURFACE ALBEDO INCLUDING SNOW EFFECT (UNITLESS FRACTION)
1535 ! =SNOW-FREE ALBEDO (ALB) WHEN SNEQV=0, OR
1536 ! =FCT(MSNOALB,ALB,VEGTYP,SHDFAC,SHDMIN) WHEN SNEQV>0
1537 ! CH SURFACE EXCHANGE COEFFICIENT FOR HEAT AND MOISTURE
1538 ! (M S-1); NOTE: CH IS TECHNICALLY A CONDUCTANCE SINCE
1539 ! IT HAS BEEN MULTIPLIED BY WIND SPEED.
1540 ! CM SURFACE EXCHANGE COEFFICIENT FOR MOMENTUM (M S-1); NOTE:
1541 ! CM IS TECHNICALLY A CONDUCTANCE SINCE IT HAS BEEN
1542 ! MULTIPLIED BY WIND SPEED.
1543 ! ----------------------------------------------------------------------
1544 ! 7. OUTPUT (O):
1545 ! ----------------------------------------------------------------------
1546 ! OUTPUT VARIABLES NECESSARY FOR A COUPLED NUMERICAL WEATHER PREDICTION
1547 ! MODEL, E.G. NOAA/NWS/NCEP MESOSCALE ETA MODEL. FOR THIS APPLICATION,
1548 ! THE REMAINING OUTPUT/DIAGNOSTIC/PARAMETER BLOCKS BELOW ARE NOT
1549 ! NECESSARY. OTHER APPLICATIONS MAY REQUIRE DIFFERENT OUTPUT VARIABLES.
1550 ! ETA ACTUAL LATENT HEAT FLUX (W m-2: NEGATIVE, IF UP FROM
1551 ! SURFACE)
1552 ! ETA_KINEMATIC atctual latent heat flux in Kg m-2 s-1
1553 ! SHEAT SENSIBLE HEAT FLUX (W M-2: NEGATIVE, IF UPWARD FROM
1554 ! SURFACE)
1555 ! FDOWN Radiation forcing at the surface (W m-2) = SOLDN*(1-alb)+LWDN
1556 ! ----------------------------------------------------------------------
1557 ! EC CANOPY WATER EVAPORATION (W m-2)
1558 ! EDIR DIRECT SOIL EVAPORATION (W m-2)
1559 ! ET(NSOIL) PLANT TRANSPIRATION FROM A PARTICULAR ROOT (SOIL) LAYER
1560 ! (W m-2)
1561 ! ETT TOTAL PLANT TRANSPIRATION (W m-2)
1562 ! ESNOW SUBLIMATION FROM (OR DEPOSITION TO IF <0) SNOWPACK
1563 ! (W m-2)
1564 ! DRIP THROUGH-FALL OF PRECIP AND/OR DEW IN EXCESS OF CANOPY
1565 ! WATER-HOLDING CAPACITY (M)
1566 ! DEW DEWFALL (OR FROSTFALL FOR T<273.15) (M)
1567 ! ----------------------------------------------------------------------
1568 ! BETA RATIO OF ACTUAL/POTENTIAL EVAP (DIMENSIONLESS)
1569 ! ETP POTENTIAL EVAPORATION (W m-2)
1570 ! SSOIL SOIL HEAT FLUX (W M-2: NEGATIVE IF DOWNWARD FROM SURFACE)
1571 ! ----------------------------------------------------------------------
1572 ! FLX1 PRECIP-SNOW SFC (W M-2)
1573 ! FLX2 FREEZING RAIN LATENT HEAT FLUX (W M-2)
1574 ! FLX3 PHASE-CHANGE HEAT FLUX FROM SNOWMELT (W M-2)
1575 ! ----------------------------------------------------------------------
1576 ! SNOMLT SNOW MELT (M) (WATER EQUIVALENT)
1577 ! SNCOVR FRACTIONAL SNOW COVER (UNITLESS FRACTION, 0-1)
1578 ! ----------------------------------------------------------------------
1579 ! RUNOFF1 SURFACE RUNOFF (M S-1), NOT INFILTRATING THE SURFACE
1580 ! RUNOFF2 SUBSURFACE RUNOFF (M S-1), DRAINAGE OUT BOTTOM OF LAST
1581 ! SOIL LAYER (BASEFLOW)
1582 ! RUNOFF3 NUMERICAL TRUNCTATION IN EXCESS OF POROSITY (SMCMAX)
1583 ! FOR A GIVEN SOIL LAYER AT THE END OF A TIME STEP (M S-1).
1584 ! Note: the above RUNOFF2 is actually the sum of RUNOFF2 and RUNOFF3
1585 ! ----------------------------------------------------------------------
1586 ! RC CANOPY RESISTANCE (S M-1)
1587 ! PC PLANT COEFFICIENT (UNITLESS FRACTION, 0-1) WHERE PC*ETP
1588 ! = ACTUAL TRANSP
1589 ! XLAI LEAF AREA INDEX (DIMENSIONLESS)
1590 ! RSMIN MINIMUM CANOPY RESISTANCE (S M-1)
1591 ! RCS INCOMING SOLAR RC FACTOR (DIMENSIONLESS)
1592 ! RCT AIR TEMPERATURE RC FACTOR (DIMENSIONLESS)
1593 ! RCQ ATMOS VAPOR PRESSURE DEFICIT RC FACTOR (DIMENSIONLESS)
1594 ! RCSOIL SOIL MOISTURE RC FACTOR (DIMENSIONLESS)
1595 ! ----------------------------------------------------------------------
1596 ! 8. DIAGNOSTIC OUTPUT (D):
1597 ! ----------------------------------------------------------------------
1598 ! SOILW AVAILABLE SOIL MOISTURE IN ROOT ZONE (UNITLESS FRACTION
1599 ! BETWEEN SMCWLT AND SMCMAX)
1600 ! SOILM TOTAL SOIL COLUMN MOISTURE CONTENT (FROZEN+UNFROZEN) (M)
1601 ! Q1 Effective mixing ratio at surface (kg kg-1), used for
1602 ! diagnosing the mixing ratio at 2 meter for coupled model
1603 ! ----------------------------------------------------------------------
1604 ! 9. PARAMETERS (P):
1605 ! ----------------------------------------------------------------------
1606 ! SMCWLT WILTING POINT (VOLUMETRIC)
1607 ! SMCDRY DRY SOIL MOISTURE THRESHOLD WHERE DIRECT EVAP FRM TOP
1608 ! LAYER ENDS (VOLUMETRIC)
1609 ! SMCREF SOIL MOISTURE THRESHOLD WHERE TRANSPIRATION BEGINS TO
1610 ! STRESS (VOLUMETRIC)
1611 ! SMCMAX POROSITY, I.E. SATURATED VALUE OF SOIL MOISTURE
1612 ! (VOLUMETRIC)
1613 ! NROOT NUMBER OF ROOT LAYERS, A FUNCTION OF VEG TYPE, DETERMINED
1614 ! IN SUBROUTINE REDPRM.
1615 ! ----------------------------------------------------------------------
1616
1617
1618 IMPLICIT NONE
1619 ! ----------------------------------------------------------------------
1620
1621 ! DECLARATIONS - LOGICAL AND CHARACTERS
1622 ! ----------------------------------------------------------------------
1623 LOGICAL, INTENT(IN):: LOCAL
1624 LOGICAL :: FRZGRA, SNOWNG
1625 CHARACTER (LEN=4), INTENT(IN):: LLANDUSE, LSOIL
1626
1627 ! ----------------------------------------------------------------------
1628 ! DECLARATIONS - INTEGER
1629 ! ----------------------------------------------------------------------
1630 INTEGER,INTENT(IN) :: ICE,NSOIL,SLOPETYP,SOILTYP,VEGTYP
1631 INTEGER,INTENT(OUT):: NROOT
1632 INTEGER KZ, K, iout
1633
1634 ! ----------------------------------------------------------------------
1635 ! DECLARATIONS - REAL
1636 ! ----------------------------------------------------------------------
1637 REAL, INTENT(IN) :: SHDMIN,SHDMAX,DT,DQSDT2,LWDN,PRCP,PRCPRAIN, &
1638 Q2,Q2SAT,SFCPRS,SFCSPD,SFCTMP, SNOALB, &
1639 SOLDN,SOLNET,TBOT,TH2,ZLVL, &
1640 EMISSI, FFROZP
1641 REAL, INTENT(INOUT):: COSZ, SOLARDIRECT,ALBEDO,CH,CM, &
1642 CMC,SNEQV,SNCOVR,SNOWH,T1,XLAI,SHDFAC,Z0BRD,ALB
1643 REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SLDPTH
1644 REAL, DIMENSION(1:NSOIL), INTENT(OUT):: ET
1645 REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: SH2O, SMC, STC
1646 REAL,DIMENSION(1:NSOIL):: RTDIS, ZSOIL
1647
1648 REAL,INTENT(OUT) :: ETA_KINEMATIC,BETA,DEW,DRIP,EC,EDIR,ESNOW,ETA, &
1649 ETP,FLX1,FLX2,FLX3,SHEAT,PC,RUNOFF1,RUNOFF2, &
1650 RUNOFF3,RC,RSMIN,RCQ,RCS,RCSOIL,RCT,SSOIL, &
1651 SMCDRY,SMCMAX,SMCREF,SMCWLT,SNOMLT, SOILM, &
1652 SOILW,FDOWN,Q1
1653 REAL :: BEXP,CFACTR,CMCMAX,CSOIL,CZIL,DF1,DF1H,DF1A,DKSAT,DWSAT, &
1654 DSOIL,DTOT,ETT,FRCSNO,FRCSOI,EPSCA,F1,FXEXP,FRZX,HS, &
1655 KDT,LVH2O,PRCP1,PSISAT,QUARTZ,R,RCH,REFKDT,RR,RGL, &
1656 RSMAX, &
1657 RSNOW,SNDENS,SNCOND,SBETA,SN_NEW,SLOPE,SNUP,SALP,SOILWM, &
1658 SOILWW,T1V,T24,T2V,TH2V,TOPT,TFREEZ,TSNOW,ZBOT,Z0,PRCPF, &
1659 ETNS,PTU,LSUBS
1660
1661 ! ----------------------------------------------------------------------
1662 ! DECLARATIONS - PARAMETERS
1663 ! ----------------------------------------------------------------------
1664 PARAMETER (TFREEZ = 273.15)
1665 PARAMETER (LVH2O = 2.501E+6)
1666 PARAMETER (LSUBS = 2.83E+6)
1667 PARAMETER (R = 287.04)
1668 ! ----------------------------------------------------------------------
1669 ! INITIALIZATION
1670 ! ----------------------------------------------------------------------
1671 RUNOFF1 = 0.0
1672 RUNOFF2 = 0.0
1673 RUNOFF3 = 0.0
1674 SNOMLT = 0.0
1675
1676 ! ----------------------------------------------------------------------
1677 ! THE VARIABLE "ICE" IS A FLAG DENOTING SEA-ICE CASE
1678 ! ----------------------------------------------------------------------
1679 ! SEA-ICE LAYERS ARE EQUAL THICKNESS AND SUM TO 3 METERS
1680 ! ----------------------------------------------------------------------
1681 IF (ICE == 1) THEN
1682 DO KZ = 1,NSOIL
1683 ZSOIL (KZ) = -3.* FLOAT (KZ)/ FLOAT (NSOIL)
1684 END DO
1685
1686 ! ----------------------------------------------------------------------
1687 ! CALCULATE DEPTH (NEGATIVE) BELOW GROUND FROM TOP SKIN SFC TO BOTTOM OF
1688 ! EACH SOIL LAYER. NOTE: SIGN OF ZSOIL IS NEGATIVE (DENOTING BELOW
1689 ! GROUND)
1690 ! ----------------------------------------------------------------------
1691 ELSE
1692 ZSOIL (1) = - SLDPTH (1)
1693 DO KZ = 2,NSOIL
1694 ZSOIL (KZ) = - SLDPTH (KZ) + ZSOIL (KZ -1)
1695 END DO
1696 ! ----------------------------------------------------------------------
1697 ! NEXT IS CRUCIAL CALL TO SET THE LAND-SURFACE PARAMETERS, INCLUDING
1698 ! SOIL-TYPE AND VEG-TYPE DEPENDENT PARAMETERS.
1699 ! ----------------------------------------------------------------------
1700 CALL REDPRM (VEGTYP,SOILTYP,SLOPETYP,CFACTR,CMCMAX,RSMAX,TOPT, &
1701 REFKDT,KDT,SBETA, SHDFAC,RSMIN,RGL,HS,ZBOT,FRZX, &
1702 PSISAT,SLOPE,SNUP,SALP,BEXP,DKSAT,DWSAT, &
1703 SMCMAX,SMCWLT,SMCREF,SMCDRY,F1,QUARTZ,FXEXP, &
1704 RTDIS,SLDPTH,ZSOIL,NROOT,NSOIL,Z0BRD,CZIL,XLAI, &
1705 CSOIL,ALB,PTU,LLANDUSE,LSOIL,LOCAL)
1706 END IF
1707
1708 !urban
1709 IF(VEGTYP==1)THEN
1710 SHDFAC=0.05
1711 RSMIN=400.0
1712 SMCMAX = 0.45
1713 SMCREF = 0.42
1714 SMCWLT = 0.40
1715 ENDIF
1716
1717 ! ----------------------------------------------------------------------
1718 ! INITIALIZE PRECIPITATION LOGICALS.
1719 ! ----------------------------------------------------------------------
1720 SNOWNG = .FALSE.
1721 FRZGRA = .FALSE.
1722
1723 ! ----------------------------------------------------------------------
1724 ! IF SEA-ICE CASE, ASSIGN DEFAULT WATER-EQUIV SNOW ON TOP
1725 ! ----------------------------------------------------------------------
1726 IF (ICE == 1) THEN
1727 SNEQV = 0.01
1728 SNOWH = 0.05
1729 ! ----------------------------------------------------------------------
1730 ! IF INPUT SNOWPACK IS NONZERO, THEN COMPUTE SNOW DENSITY "SNDENS" AND
1731 ! SNOW THERMAL CONDUCTIVITY "SNCOND" (NOTE THAT CSNOW IS A FUNCTION
1732 ! SUBROUTINE)
1733 ! ----------------------------------------------------------------------
1734 END IF
1735 IF (SNEQV == 0.0) THEN
1736 SNDENS = 0.0
1737 SNOWH = 0.0
1738 SNCOND = 1.0
1739 ELSE
1740 SNDENS = SNEQV / SNOWH
1741 IF(SNDENS > 1.0) THEN
1742 CALL wrf_error_fatal ( 'Physical snow depth is less than snow water equiv.' )
1743 ENDIF
1744 CALL CSNOW (SNCOND,SNDENS)
1745 END IF
1746 ! ----------------------------------------------------------------------
1747 ! DETERMINE IF IT'S PRECIPITATING AND WHAT KIND OF PRECIP IT IS.
1748 ! IF IT'S PRCPING AND THE AIR TEMP IS COLDER THAN 0 C, IT'S SNOWING!
1749 ! IF IT'S PRCPING AND THE AIR TEMP IS WARMER THAN 0 C, BUT THE GRND
1750 ! TEMP IS COLDER THAN 0 C, FREEZING RAIN IS PRESUMED TO BE FALLING.
1751 ! ----------------------------------------------------------------------
1752 IF (PRCP > 0.0) THEN
1753 ! snow defined when fraction of frozen precip (FFROZP) > 0.5,
1754 ! passed in from model microphysics.
1755 IF (FFROZP .GT. 0.5) THEN
1756 SNOWNG = .TRUE.
1757 ELSE
1758 IF (T1 <= TFREEZ) FRZGRA = .TRUE.
1759 END IF
1760 END IF
1761 ! ----------------------------------------------------------------------
1762 ! IF EITHER PRCP FLAG IS SET, DETERMINE NEW SNOWFALL (CONVERTING PRCP
1763 ! RATE FROM KG M-2 S-1 TO A LIQUID EQUIV SNOW DEPTH IN METERS) AND ADD
1764 ! IT TO THE EXISTING SNOWPACK.
1765 ! NOTE THAT SINCE ALL PRECIP IS ADDED TO SNOWPACK, NO PRECIP INFILTRATES
1766 ! INTO THE SOIL SO THAT PRCP1 IS SET TO ZERO.
1767 ! ----------------------------------------------------------------------
1768 IF ( (SNOWNG) .OR. (FRZGRA) ) THEN
1769 SN_NEW = PRCP * DT * 0.001
1770 SNEQV = SNEQV + SN_NEW
1771 PRCPF = 0.0
1772
1773 ! ----------------------------------------------------------------------
1774 ! UPDATE SNOW DENSITY BASED ON NEW SNOWFALL, USING OLD AND NEW SNOW.
1775 ! UPDATE SNOW THERMAL CONDUCTIVITY
1776 ! ----------------------------------------------------------------------
1777 CALL SNOW_NEW (SFCTMP,SN_NEW,SNOWH,SNDENS)
1778 CALL CSNOW (SNCOND,SNDENS)
1779
1780 ! ----------------------------------------------------------------------
1781 ! PRECIP IS LIQUID (RAIN), HENCE SAVE IN THE PRECIP VARIABLE THAT
1782 ! LATER CAN WHOLELY OR PARTIALLY INFILTRATE THE SOIL (ALONG WITH
1783 ! ANY CANOPY "DRIP" ADDED TO THIS LATER)
1784 ! ----------------------------------------------------------------------
1785 ELSE
1786 PRCPF = PRCP
1787 END IF
1788 ! ----------------------------------------------------------------------
1789 ! DETERMINE SNOWCOVER AND ALBEDO OVER LAND.
1790 ! ----------------------------------------------------------------------
1791 ! ----------------------------------------------------------------------
1792 ! IF SNOW DEPTH=0, SET SNOW FRACTION=0, ALBEDO=SNOW FREE ALBEDO.
1793 ! ----------------------------------------------------------------------
1794 IF (ICE == 0) THEN
1795 IF (SNEQV == 0.0) THEN
1796 SNCOVR = 0.0
1797 ALBEDO = ALB
1798 ELSE
1799 ! ----------------------------------------------------------------------
1800 ! DETERMINE SNOW FRACTIONAL COVERAGE.
1801 ! DETERMINE SURFACE ALBEDO MODIFICATION DUE TO SNOWDEPTH STATE.
1802 ! ----------------------------------------------------------------------
1803 CALL SNFRAC (SNEQV,SNUP,SALP,SNOWH,SNCOVR)
1804 ! limit snow cover fraction to maximum of 0.98
1805 SNCOVR = MIN(SNCOVR,0.98)
1806 CALL ALCALC (ALB,SNOALB,SHDFAC,SHDMIN,SNCOVR,TSNOW,ALBEDO)
1807 END IF
1808 ! ----------------------------------------------------------------------
1809 ! SNOW COVER, ALBEDO OVER SEA-ICE
1810 ! ----------------------------------------------------------------------
1811 ELSE
1812 SNCOVR = 1.0
1813 ALBEDO = 0.65
1814 END IF
1815 ! ----------------------------------------------------------------------
1816 ! THERMAL CONDUCTIVITY FOR SEA-ICE CASE
1817 ! ----------------------------------------------------------------------
1818 IF (ICE == 1) THEN
1819 DF1 = 2.2
1820 ELSE
1821 ! ----------------------------------------------------------------------
1822 ! NEXT CALCULATE THE SUBSURFACE HEAT FLUX, WHICH FIRST REQUIRES
1823 ! CALCULATION OF THE THERMAL DIFFUSIVITY. TREATMENT OF THE
1824 ! LATTER FOLLOWS THAT ON PAGES 148-149 FROM "HEAT TRANSFER IN
1825 ! COLD CLIMATES", BY V. J. LUNARDINI (PUBLISHED IN 1981
1826 ! BY VAN NOSTRAND REINHOLD CO.) I.E. TREATMENT OF TWO CONTIGUOUS
1827 ! "PLANE PARALLEL" MEDIUMS (NAMELY HERE THE FIRST SOIL LAYER
1828 ! AND THE SNOWPACK LAYER, IF ANY). THIS DIFFUSIVITY TREATMENT
1829 ! BEHAVES WELL FOR BOTH ZERO AND NONZERO SNOWPACK, INCLUDING THE
1830 ! LIMIT OF VERY THIN SNOWPACK. THIS TREATMENT ALSO ELIMINATES
1831 ! THE NEED TO IMPOSE AN ARBITRARY UPPER BOUND ON SUBSURFACE
1832 ! HEAT FLUX WHEN THE SNOWPACK BECOMES EXTREMELY THIN.
1833 ! ----------------------------------------------------------------------
1834 ! FIRST CALCULATE THERMAL DIFFUSIVITY OF TOP SOIL LAYER, USING
1835 ! BOTH THE FROZEN AND LIQUID SOIL MOISTURE, FOLLOWING THE
1836 ! SOIL THERMAL DIFFUSIVITY FUNCTION OF PETERS-LIDARD ET AL.
1837 ! (1998,JAS, VOL 55, 1209-1224), WHICH REQUIRES THE SPECIFYING
1838 ! THE QUARTZ CONTENT OF THE GIVEN SOIL CLASS (SEE ROUTINE REDPRM)
1839 ! ----------------------------------------------------------------------
1840 ! ----------------------------------------------------------------------
1841 ! NEXT ADD SUBSURFACE HEAT FLUX REDUCTION EFFECT FROM THE
1842 ! OVERLYING GREEN CANOPY, ADAPTED FROM SECTION 2.1.2 OF
1843 ! PETERS-LIDARD ET AL. (1997, JGR, VOL 102(D4))
1844 ! ----------------------------------------------------------------------
1845 CALL TDFCND (DF1,SMC (1),QUARTZ,SMCMAX,SH2O (1))
1846
1847 !urban
1848 IF (VEGTYP==1) DF1=3.24
1849
1850 DF1 = DF1 * EXP (SBETA * SHDFAC)
1851 ! ----------------------------------------------------------------------
1852 ! FINALLY "PLANE PARALLEL" SNOWPACK EFFECT FOLLOWING
1853 ! V.J. LINARDINI REFERENCE CITED ABOVE. NOTE THAT DTOT IS
1854 ! COMBINED DEPTH OF SNOWDEPTH AND THICKNESS OF FIRST SOIL LAYER
1855 ! ----------------------------------------------------------------------
1856 END IF
1857
1858 DSOIL = - (0.5 * ZSOIL (1))
1859 IF (SNEQV == 0.) THEN
1860 SSOIL = DF1 * (T1- STC (1) ) / DSOIL
1861 ELSE
1862 DTOT = SNOWH + DSOIL
1863 FRCSNO = SNOWH / DTOT
1864
1865 ! 1. HARMONIC MEAN (SERIES FLOW)
1866 ! DF1 = (SNCOND*DF1)/(FRCSOI*SNCOND+FRCSNO*DF1)
1867 FRCSOI = DSOIL / DTOT
1868 ! 2. ARITHMETIC MEAN (PARALLEL FLOW)
1869 ! DF1 = FRCSNO*SNCOND + FRCSOI*DF1
1870 DF1H = (SNCOND * DF1)/ (FRCSOI * SNCOND+ FRCSNO * DF1)
1871
1872 ! 3. GEOMETRIC MEAN (INTERMEDIATE BETWEEN HARMONIC AND ARITHMETIC MEAN)
1873 ! DF1 = (SNCOND**FRCSNO)*(DF1**FRCSOI)
1874 ! weigh DF by snow fraction
1875 ! DF1 = DF1H*SNCOVR + DF1A*(1.0-SNCOVR)
1876 ! DF1 = DF1H*SNCOVR + DF1*(1.0-SNCOVR)
1877 DF1A = FRCSNO * SNCOND+ FRCSOI * DF1
1878
1879 ! ----------------------------------------------------------------------
1880 ! CALCULATE SUBSURFACE HEAT FLUX, SSOIL, FROM FINAL THERMAL DIFFUSIVITY
1881 ! OF SURFACE MEDIUMS, DF1 ABOVE, AND SKIN TEMPERATURE AND TOP
1882 ! MID-LAYER SOIL TEMPERATURE
1883 ! ----------------------------------------------------------------------
1884 DF1 = DF1A * SNCOVR + DF1* (1.0- SNCOVR)
1885 SSOIL = DF1 * (T1- STC (1) ) / DTOT
1886 END IF
1887 ! ----------------------------------------------------------------------
1888 ! DETERMINE SURFACE ROUGHNESS OVER SNOWPACK USING SNOW CONDITION FROM
1889 ! THE PREVIOUS TIMESTEP.
1890 ! ----------------------------------------------------------------------
1891 IF (SNCOVR > 0. ) THEN
1892 CALL SNOWZ0 (SNCOVR,Z0,Z0BRD)
1893 ELSE
1894 Z0=Z0BRD
1895 END IF
1896 ! ----------------------------------------------------------------------
1897 ! NEXT CALL ROUTINE SFCDIF TO CALCULATE THE SFC EXCHANGE COEF (CH) FOR
1898 ! HEAT AND MOISTURE.
1899
1900 ! NOTE !!!
1901 ! DO NOT CALL SFCDIF UNTIL AFTER ABOVE CALL TO REDPRM, IN CASE
1902 ! ALTERNATIVE VALUES OF ROUGHNESS LENGTH (Z0) AND ZILINTINKEVICH COEF
1903 ! (CZIL) ARE SET THERE VIA NAMELIST I/O.
1904
1905 ! NOTE !!!
1906 ! ROUTINE SFCDIF RETURNS A CH THAT REPRESENTS THE WIND SPD TIMES THE
1907 ! "ORIGINAL" NONDIMENSIONAL "Ch" TYPICAL IN LITERATURE. HENCE THE CH
1908 ! RETURNED FROM SFCDIF HAS UNITS OF M/S. THE IMPORTANT COMPANION
1909 ! COEFFICIENT OF CH, CARRIED HERE AS "RCH", IS THE CH FROM SFCDIF TIMES
1910 ! AIR DENSITY AND PARAMETER "CP". "RCH" IS COMPUTED IN "CALL PENMAN".
1911 ! RCH RATHER THAN CH IS THE COEFF USUALLY INVOKED LATER IN EQNS.
1912
1913 ! NOTE !!!
1914 ! ----------------------------------------------------------------------
1915 ! SFCDIF ALSO RETURNS THE SURFACE EXCHANGE COEFFICIENT FOR MOMENTUM, CM,
1916 ! ALSO KNOWN AS THE SURFACE DRAGE COEFFICIENT. Needed as a state variable
1917 ! for iterative/implicit solution of CH in SFCDIF
1918 ! ----------------------------------------------------------------------
1919 ! IF(.NOT.LCH) THEN
1920 ! T1V = T1 * (1.0+ 0.61 * Q2)
1921 ! TH2V = TH2 * (1.0+ 0.61 * Q2)
1922 ! CALL SFCDIF_off (ZLVL,Z0,T1V,TH2V,SFCSPD,CZIL,CM,CH)
1923 ! ENDIF
1924
1925 ! ----------------------------------------------------------------------
1926 ! CALL PENMAN SUBROUTINE TO CALCULATE POTENTIAL EVAPORATION (ETP), AND
1927 ! OTHER PARTIAL PRODUCTS AND SUMS SAVE IN COMMON/RITE FOR LATER
1928 ! CALCULATIONS.
1929 ! ----------------------------------------------------------------------
1930 ! ----------------------------------------------------------------------
1931 ! CALCULATE TOTAL DOWNWARD RADIATION (SOLAR PLUS LONGWAVE) NEEDED IN
1932 ! PENMAN EP SUBROUTINE THAT FOLLOWS
1933 ! ----------------------------------------------------------------------
1934 ! FDOWN = SOLDN * (1.0- ALBEDO) + LWDN
1935 FDOWN = SOLNET + LWDN
1936 ! ----------------------------------------------------------------------
1937 ! CALC VIRTUAL TEMPS AND VIRTUAL POTENTIAL TEMPS NEEDED BY SUBROUTINES
1938 ! PENMAN.
1939 T2V = SFCTMP * (1.0+ 0.61 * Q2 )
1940
1941 iout=0
1942 if(iout.eq.1) then
1943 print*,'before penman'
1944 print*,' SFCTMP',SFCTMP,'SFCPRS',SFCPRS,'CH',CH,'T2V',T2V, &
1945 'TH2',TH2,'PRCP',PRCP,'FDOWN',FDOWN,'T24',T24,'SSOIL',SSOIL, &
1946 'Q2',Q2,'Q2SAT',Q2SAT,'ETP',ETP,'RCH',RCH, &
1947 'EPSCA',EPSCA,'RR',RR ,'SNOWNG',SNOWNG,'FRZGRA',FRZGRA, &
1948 'DQSDT2',DQSDT2,'FLX2',FLX2,'SNOWH',SNOWH,'SNEQV',SNEQV, &
1949 ' DSOIL',DSOIL,' FRCSNO',FRCSNO,' SNCOVR',SNCOVR,' DTOT',DTOT, &
1950 ' ZSOIL (1)',ZSOIL(1),' DF1',DF1,'T1',T1,' STC1',STC(1), &
1951 'ALBEDO',ALBEDO,'SMC',SMC,'STC',STC,'SH2O',SH2O
1952 endif
1953
1954 CALL PENMAN (SFCTMP,SFCPRS,CH,T2V,TH2,PRCP,FDOWN,T24,SSOIL, &
1955 Q2,Q2SAT,ETP,RCH,EPSCA,RR,SNOWNG,FRZGRA, &
1956 DQSDT2,FLX2,EMISSI,SNEQV)
1957
1958 ! ----------------------------------------------------------------------
1959 ! CALL CANRES TO CALCULATE THE CANOPY RESISTANCE AND CONVERT IT INTO PC
1960 ! IF NONZERO GREENNESS FRACTION
1961 ! ----------------------------------------------------------------------
1962
1963 ! ----------------------------------------------------------------------
1964 ! FROZEN GROUND EXTENSION: TOTAL SOIL WATER "SMC" WAS REPLACED
1965 ! BY UNFROZEN SOIL WATER "SH2O" IN CALL TO CANRES BELOW
1966 ! ----------------------------------------------------------------------
1967 IF (SHDFAC > 0.) THEN
1968 CALL CANRES (SOLDN,CH,SFCTMP,Q2,SFCPRS,SH2O,ZSOIL,NSOIL, &
1969 SMCWLT,SMCREF,RSMIN,RC,PC,NROOT,Q2SAT,DQSDT2, &
1970 TOPT,RSMAX,RGL,HS,XLAI, &
1971 RCS,RCT,RCQ,RCSOIL,EMISSI)
1972 END IF
1973 ! ----------------------------------------------------------------------
1974 ! NOW DECIDE MAJOR PATHWAY BRANCH TO TAKE DEPENDING ON WHETHER SNOWPACK
1975 ! EXISTS OR NOT:
1976 ! ----------------------------------------------------------------------
1977 ESNOW = 0.0
1978 IF (SNEQV == 0.0) THEN
1979 CALL NOPAC (ETP,ETA,PRCP,SMC,SMCMAX,SMCWLT, &
1980 SMCREF,SMCDRY,CMC,CMCMAX,NSOIL,DT, &
1981 SHDFAC, &
1982 SBETA,Q2,T1,SFCTMP,T24,TH2,FDOWN,F1,EMISSI, &
1983 SSOIL, &
1984 STC,EPSCA,BEXP,PC,RCH,RR,CFACTR, &
1985 SH2O,SLOPE,KDT,FRZX,PSISAT,ZSOIL, &
1986 DKSAT,DWSAT,TBOT,ZBOT,RUNOFF1,RUNOFF2, &
1987 RUNOFF3,EDIR,EC,ET,ETT,NROOT,ICE,RTDIS, &
1988 QUARTZ,FXEXP,CSOIL, &
1989 BETA,DRIP,DEW,FLX1,FLX3,VEGTYP)
1990 ETA_KINEMATIC = ETA
1991 ELSE
1992 CALL SNOPAC (ETP,ETA,PRCP,PRCPF,SNOWNG,SMC,SMCMAX,SMCWLT, &
1993 SMCREF,SMCDRY,CMC,CMCMAX,NSOIL,DT, &
1994 SBETA,DF1, &
1995 Q2,T1,SFCTMP,T24,TH2,FDOWN,F1,SSOIL,STC,EPSCA, &
1996 SFCPRS,BEXP,PC,RCH,RR,CFACTR,SNCOVR,SNEQV,SNDENS,&
1997 SNOWH,SH2O,SLOPE,KDT,FRZX,PSISAT, &
1998 ZSOIL,DWSAT,DKSAT,TBOT,ZBOT,SHDFAC,RUNOFF1, &
1999 RUNOFF2,RUNOFF3,EDIR,EC,ET,ETT,NROOT,SNOMLT, &
2000 ICE,RTDIS,QUARTZ,FXEXP,CSOIL, &
2001 BETA,DRIP,DEW,FLX1,FLX2,FLX3,ESNOW,ETNS,EMISSI, &
2002 VEGTYP)
2003 ETA_KINEMATIC = ESNOW + ETNS
2004 END IF
2005
2006 ! Calculate effective mixing ratio at grnd level (skin)
2007 !
2008 ! Q1=Q2+ETA*CP/RCH
2009 Q1=Q2+ETA_KINEMATIC*CP/RCH
2010 !
2011 ! ----------------------------------------------------------------------
2012 ! DETERMINE SENSIBLE HEAT (H) IN ENERGY UNITS (W M-2)
2013 ! ----------------------------------------------------------------------
2014 SHEAT = - (CH * CP * SFCPRS)/ (R * T2V) * ( TH2- T1 )
2015
2016 ! ----------------------------------------------------------------------
2017 ! CONVERT EVAP TERMS FROM KINEMATIC (KG M-2 S-1) TO ENERGY UNITS (W M-2)
2018 ! ----------------------------------------------------------------------
2019 EDIR = EDIR * LVH2O
2020 EC = EC * LVH2O
2021 DO K=1,4
2022 ET(K) = ET(K) * LVH2O
2023 ENDDO
2024 ETT = ETT * LVH2O
2025 ESNOW = ESNOW * LSUBS
2026 ETP = ETP*((1.-SNCOVR)*LVH2O + SNCOVR*LSUBS)
2027 IF (ETP .GT. 0.) THEN
2028 ETA = EDIR + EC + ETT + ESNOW
2029 ELSE
2030 ETA = ETP
2031 ENDIF
2032 ! ----------------------------------------------------------------------
2033 ! DETERMINE BETA (RATIO OF ACTUAL TO POTENTIAL EVAP)
2034 ! ----------------------------------------------------------------------
2035 IF (ETP == 0.0) THEN
2036 BETA = 0.0
2037 ELSE
2038 BETA = ETA/ETP
2039 ENDIF
2040
2041 ! ----------------------------------------------------------------------
2042 ! CONVERT THE SIGN OF SOIL HEAT FLUX SO THAT:
2043 ! SSOIL>0: WARM THE SURFACE (NIGHT TIME)
2044 ! SSOIL<0: COOL THE SURFACE (DAY TIME)
2045 ! ----------------------------------------------------------------------
2046 SSOIL = -1.0* SSOIL
2047
2048 ! ----------------------------------------------------------------------
2049 ! CONVERT RUNOFF3 (INTERNAL LAYER RUNOFF FROM SUPERSAT) FROM M TO M S-1
2050 ! AND ADD TO SUBSURFACE RUNOFF/DRAINAGE/BASEFLOW
2051 ! ----------------------------------------------------------------------
2052 RUNOFF3 = RUNOFF3/ DT
2053 RUNOFF2 = RUNOFF2+ RUNOFF3
2054
2055 ! ----------------------------------------------------------------------
2056 ! TOTAL COLUMN SOIL MOISTURE IN METERS (SOILM) AND ROOT-ZONE
2057 ! SOIL MOISTURE AVAILABILITY (FRACTION) RELATIVE TO POROSITY/SATURATION
2058 ! ----------------------------------------------------------------------
2059 IF (ICE == 0) THEN
2060 SOILM = -1.0* SMC (1)* ZSOIL (1)
2061 DO K = 2,NSOIL
2062 SOILM = SOILM + SMC (K)* (ZSOIL (K -1) - ZSOIL (K))
2063 END DO
2064 SOILWM = -1.0* (SMCMAX - SMCWLT)* ZSOIL (1)
2065 SOILWW = -1.0* (SMC (1) - SMCWLT)* ZSOIL (1)
2066
2067 IF (NROOT >= 2) THEN
2068 DO K = 2,NROOT
2069 SOILWM = SOILWM + (SMCMAX - SMCWLT)* (ZSOIL (K -1) - ZSOIL (K))
2070 SOILWW = SOILWW + (SMC(K) - SMCWLT)* (ZSOIL (K -1) - ZSOIL (K))
2071 END DO
2072 END IF
2073 IF (SOILWM .LT. 1.E-6) THEN
2074 SOILWM = 0.0
2075 SOILW = 0.0
2076 SOILM = 0.0
2077 ELSE
2078 SOILW = SOILWW / SOILWM
2079 END IF
2080 ELSE
2081 SOILWM = 0.0
2082 SOILW = 0.0
2083 SOILM = 0.0
2084 END IF
2085
2086 ! ----------------------------------------------------------------------
2087 END SUBROUTINE SFLX
2088 ! ----------------------------------------------------------------------
2089
2090 SUBROUTINE ALCALC (ALB,SNOALB,SHDFAC,SHDMIN,SNCOVR,TSNOW,ALBEDO)
2091
2092 ! ----------------------------------------------------------------------
2093 ! CALCULATE ALBEDO INCLUDING SNOW EFFECT (0 -> 1)
2094 ! ALB SNOWFREE ALBEDO
2095 ! SNOALB MAXIMUM (DEEP) SNOW ALBEDO
2096 ! SHDFAC AREAL FRACTIONAL COVERAGE OF GREEN VEGETATION
2097 ! SHDMIN MINIMUM AREAL FRACTIONAL COVERAGE OF GREEN VEGETATION
2098 ! SNCOVR FRACTIONAL SNOW COVER
2099 ! ALBEDO SURFACE ALBEDO INCLUDING SNOW EFFECT
2100 ! TSNOW SNOW SURFACE TEMPERATURE (K)
2101 ! ----------------------------------------------------------------------
2102 IMPLICIT NONE
2103
2104 ! ----------------------------------------------------------------------
2105 ! SNOALB IS ARGUMENT REPRESENTING MAXIMUM ALBEDO OVER DEEP SNOW,
2106 ! AS PASSED INTO SFLX, AND ADAPTED FROM THE SATELLITE-BASED MAXIMUM
2107 ! SNOW ALBEDO FIELDS PROVIDED BY D. ROBINSON AND G. KUKLA
2108 ! (1985, JCAM, VOL 24, 402-411)
2109 ! ----------------------------------------------------------------------
2110 REAL, INTENT(IN) :: ALB, SNOALB, SHDFAC, SHDMIN, SNCOVR, TSNOW
2111 REAL, INTENT(OUT) :: ALBEDO
2112 ! turn of vegetation effect
2113 ! ALBEDO = ALB + (1.0- (SHDFAC - SHDMIN))* SNCOVR * (SNOALB - ALB)
2114 ! ALBEDO = (1.0-SNCOVR)*ALB + SNCOVR*SNOALB !this is equivalent to below
2115 ALBEDO = ALB + SNCOVR*(SNOALB-ALB)
2116
2117 ! BASE FORMULATION (DICKINSON ET AL., 1986, COGLEY ET AL., 1990)
2118 ! IF (TSNOW.LE.263.16) THEN
2119 ! ALBEDO=SNOALB
2120 ! ELSE
2121 ! IF (TSNOW.LT.273.16) THEN
2122 ! TM=0.1*(TSNOW-263.16)
2123 ! ALBEDO=0.5*((0.9-0.2*(TM**3))+(0.8-0.16*(TM**3)))
2124 ! ELSE
2125 ! ALBEDO=0.67
2126 ! ENDIF
2127 ! ENDIF
2128
2129 ! ISBA FORMULATION (VERSEGHY, 1991; BAKER ET AL., 1990)
2130 ! IF (TSNOW.LT.273.16) THEN
2131 ! ALBEDO=SNOALB-0.008*DT/86400
2132 ! ELSE
2133 ! ALBEDO=(SNOALB-0.5)*EXP(-0.24*DT/86400)+0.5
2134 ! ENDIF
2135
2136 IF (ALBEDO > SNOALB) ALBEDO = SNOALB
2137
2138 ! ----------------------------------------------------------------------
2139 END SUBROUTINE ALCALC
2140 ! ----------------------------------------------------------------------
2141
2142 SUBROUTINE CANRES (SOLAR,CH,SFCTMP,Q2,SFCPRS,SMC,ZSOIL,NSOIL, &
2143 SMCWLT,SMCREF,RSMIN,RC,PC,NROOT,Q2SAT,DQSDT2, &
2144 TOPT,RSMAX,RGL,HS,XLAI, &
2145 RCS,RCT,RCQ,RCSOIL,EMISSI)
2146
2147 ! ----------------------------------------------------------------------
2148 ! SUBROUTINE CANRES
2149 ! ----------------------------------------------------------------------
2150 ! CALCULATE CANOPY RESISTANCE WHICH DEPENDS ON INCOMING SOLAR RADIATION,
2151 ! AIR TEMPERATURE, ATMOSPHERIC WATER VAPOR PRESSURE DEFICIT AT THE
2152 ! LOWEST MODEL LEVEL, AND SOIL MOISTURE (PREFERABLY UNFROZEN SOIL
2153 ! MOISTURE RATHER THAN TOTAL)
2154 ! ----------------------------------------------------------------------
2155 ! SOURCE: JARVIS (1976), NOILHAN AND PLANTON (1989, MWR), JACQUEMIN AND
2156 ! NOILHAN (1990, BLM)
2157 ! SEE ALSO: CHEN ET AL (1996, JGR, VOL 101(D3), 7251-7268), EQNS 12-14
2158 ! AND TABLE 2 OF SEC. 3.1.2
2159 ! ----------------------------------------------------------------------
2160 ! INPUT:
2161 ! SOLAR INCOMING SOLAR RADIATION
2162 ! CH SURFACE EXCHANGE COEFFICIENT FOR HEAT AND MOISTURE
2163 ! SFCTMP AIR TEMPERATURE AT 1ST LEVEL ABOVE GROUND
2164 ! Q2 AIR HUMIDITY AT 1ST LEVEL ABOVE GROUND
2165 ! Q2SAT SATURATION AIR HUMIDITY AT 1ST LEVEL ABOVE GROUND
2166 ! DQSDT2 SLOPE OF SATURATION HUMIDITY FUNCTION WRT TEMP
2167 ! SFCPRS SURFACE PRESSURE
2168 ! SMC VOLUMETRIC SOIL MOISTURE
2169 ! ZSOIL SOIL DEPTH (NEGATIVE SIGN, AS IT IS BELOW GROUND)
2170 ! NSOIL NO. OF SOIL LAYERS
2171 ! NROOT NO. OF SOIL LAYERS IN ROOT ZONE (1.LE.NROOT.LE.NSOIL)
2172 ! XLAI LEAF AREA INDEX
2173 ! SMCWLT WILTING POINT
2174 ! SMCREF REFERENCE SOIL MOISTURE (WHERE SOIL WATER DEFICIT STRESS
2175 ! SETS IN)
2176 ! RSMIN, RSMAX, TOPT, RGL, HS ARE CANOPY STRESS PARAMETERS SET IN
2177 ! SURBOUTINE REDPRM
2178 ! OUTPUT:
2179 ! PC PLANT COEFFICIENT
2180 ! RC CANOPY RESISTANCE
2181 ! ----------------------------------------------------------------------
2182
2183 IMPLICIT NONE
2184 INTEGER, INTENT(IN) :: NROOT,NSOIL
2185 INTEGER K
2186 REAL, INTENT(IN) :: CH,DQSDT2,HS,Q2,Q2SAT,RSMIN,RGL,RSMAX, &
2187 SFCPRS,SFCTMP,SMCREF,SMCWLT, SOLAR,TOPT,XLAI, &
2188 EMISSI
2189 REAL,DIMENSION(1:NSOIL), INTENT(IN) :: SMC,ZSOIL
2190 REAL, INTENT(OUT):: PC,RC,RCQ,RCS,RCSOIL,RCT
2191 REAL :: DELTA,FF,GX,P,RR
2192 REAL, DIMENSION(1:NSOIL) :: PART
2193 REAL, PARAMETER :: SLV = 2.501000E6
2194
2195
2196 ! ----------------------------------------------------------------------
2197 ! INITIALIZE CANOPY RESISTANCE MULTIPLIER TERMS.
2198 ! ----------------------------------------------------------------------
2199 RCS = 0.0
2200 RCT = 0.0
2201 RCQ = 0.0
2202 RCSOIL = 0.0
2203
2204 ! ----------------------------------------------------------------------
2205 ! CONTRIBUTION DUE TO INCOMING SOLAR RADIATION
2206 ! ----------------------------------------------------------------------
2207 RC = 0.0
2208 FF = 0.55*2.0* SOLAR / (RGL * XLAI)
2209 RCS = (FF + RSMIN / RSMAX) / (1.0+ FF)
2210
2211 ! ----------------------------------------------------------------------
2212 ! CONTRIBUTION DUE TO AIR TEMPERATURE AT FIRST MODEL LEVEL ABOVE GROUND
2213 ! RCT EXPRESSION FROM NOILHAN AND PLANTON (1989, MWR).
2214 ! ----------------------------------------------------------------------
2215 RCS = MAX (RCS,0.0001)
2216 RCT = 1.0- 0.0016* ( (TOPT - SFCTMP)**2.0)
2217
2218 ! ----------------------------------------------------------------------
2219 ! CONTRIBUTION DUE TO VAPOR PRESSURE DEFICIT AT FIRST MODEL LEVEL.
2220 ! RCQ EXPRESSION FROM SSIB
2221 ! ----------------------------------------------------------------------
2222 RCT = MAX (RCT,0.0001)
2223 RCQ = 1.0/ (1.0+ HS * (Q2SAT - Q2))
2224
2225 ! ----------------------------------------------------------------------
2226 ! CONTRIBUTION DUE TO SOIL MOISTURE AVAILABILITY.
2227 ! DETERMINE CONTRIBUTION FROM EACH SOIL LAYER, THEN ADD THEM UP.
2228 ! ----------------------------------------------------------------------
2229 RCQ = MAX (RCQ,0.01)
2230 GX = (SMC (1) - SMCWLT) / (SMCREF - SMCWLT)
2231 IF (GX > 1.) GX = 1.
2232 IF (GX < 0.) GX = 0.
2233
2234 ! ----------------------------------------------------------------------
2235 ! USE SOIL DEPTH AS WEIGHTING FACTOR
2236 ! ----------------------------------------------------------------------
2237 ! ----------------------------------------------------------------------
2238 ! USE ROOT DISTRIBUTION AS WEIGHTING FACTOR
2239 ! PART(1) = RTDIS(1) * GX
2240 ! ----------------------------------------------------------------------
2241 PART (1) = (ZSOIL (1)/ ZSOIL (NROOT)) * GX
2242 DO K = 2,NROOT
2243 GX = (SMC (K) - SMCWLT) / (SMCREF - SMCWLT)
2244 IF (GX > 1.) GX = 1.
2245 IF (GX < 0.) GX = 0.
2246 ! ----------------------------------------------------------------------
2247 ! USE SOIL DEPTH AS WEIGHTING FACTOR
2248 ! ----------------------------------------------------------------------
2249 ! ----------------------------------------------------------------------
2250 ! USE ROOT DISTRIBUTION AS WEIGHTING FACTOR
2251 ! PART(K) = RTDIS(K) * GX
2252 ! ----------------------------------------------------------------------
2253 PART (K) = ( (ZSOIL (K) - ZSOIL (K -1))/ ZSOIL (NROOT)) * GX
2254 END DO
2255 DO K = 1,NROOT
2256 RCSOIL = RCSOIL + PART (K)
2257 END DO
2258
2259 ! ----------------------------------------------------------------------
2260 ! DETERMINE CANOPY RESISTANCE DUE TO ALL FACTORS. CONVERT CANOPY
2261 ! RESISTANCE (RC) TO PLANT COEFFICIENT (PC) TO BE USED WITH POTENTIAL
2262 ! EVAP IN DETERMINING ACTUAL EVAP. PC IS DETERMINED BY:
2263 ! PC * LINERIZED PENMAN POTENTIAL EVAP =
2264 ! PENMAN-MONTEITH ACTUAL EVAPORATION (CONTAINING RC TERM).
2265 ! ----------------------------------------------------------------------
2266 RCSOIL = MAX (RCSOIL,0.0001)
2267
2268 RC = RSMIN / (XLAI * RCS * RCT * RCQ * RCSOIL)
2269 ! RR = (4.* SIGMA * RD / CP)* (SFCTMP **4.)/ (SFCPRS * CH) + 1.0
2270 RR = (4.* EMISSI *SIGMA * RD / CP)* (SFCTMP **4.)/ (SFCPRS * CH) &
2271 + 1.0
2272
2273 DELTA = (SLV / CP)* DQSDT2
2274
2275 PC = (RR + DELTA)/ (RR * (1. + RC * CH) + DELTA)
2276
2277 ! ----------------------------------------------------------------------
2278 END SUBROUTINE CANRES
2279 ! ----------------------------------------------------------------------
2280
2281 SUBROUTINE CSNOW (SNCOND,DSNOW)
2282
2283 ! ----------------------------------------------------------------------
2284 ! SUBROUTINE CSNOW
2285 ! FUNCTION CSNOW
2286 ! ----------------------------------------------------------------------
2287 ! CALCULATE SNOW TERMAL CONDUCTIVITY
2288 ! ----------------------------------------------------------------------
2289 IMPLICIT NONE
2290 REAL, INTENT(IN) :: DSNOW
2291 REAL, INTENT(OUT):: SNCOND
2292 REAL :: C
2293 REAL, PARAMETER :: UNIT = 0.11631
2294
2295 ! ----------------------------------------------------------------------
2296 ! SNCOND IN UNITS OF CAL/(CM*HR*C), RETURNED IN W/(M*C)
2297 ! CSNOW IN UNITS OF CAL/(CM*HR*C), RETURNED IN W/(M*C)
2298 ! BASIC VERSION IS DYACHKOVA EQUATION (1960), FOR RANGE 0.1-0.4
2299 ! ----------------------------------------------------------------------
2300 C = 0.328*10** (2.25* DSNOW)
2301 ! CSNOW=UNIT*C
2302
2303 ! ----------------------------------------------------------------------
2304 ! DE VAUX EQUATION (1933), IN RANGE 0.1-0.6
2305 ! ----------------------------------------------------------------------
2306 ! SNCOND=0.0293*(1.+100.*DSNOW**2)
2307 ! CSNOW=0.0293*(1.+100.*DSNOW**2)
2308
2309 ! ----------------------------------------------------------------------
2310 ! E. ANDERSEN FROM FLERCHINGER
2311 ! ----------------------------------------------------------------------
2312 ! SNCOND=0.021+2.51*DSNOW**2
2313 ! CSNOW=0.021+2.51*DSNOW**2
2314
2315 ! SNCOND = UNIT * C
2316 ! double snow thermal conductivity
2317 SNCOND = 2.0 * UNIT * C
2318
2319 ! ----------------------------------------------------------------------
2320 END SUBROUTINE CSNOW
2321 ! ----------------------------------------------------------------------
2322
2323 SUBROUTINE DEVAP (EDIR,ETP1,SMC,ZSOIL,SHDFAC,SMCMAX,BEXP, &
2324 DKSAT,DWSAT,SMCDRY,SMCREF,SMCWLT,FXEXP)
2325
2326 ! ----------------------------------------------------------------------
2327 ! SUBROUTINE DEVAP
2328 ! FUNCTION DEVAP
2329 ! ----------------------------------------------------------------------
2330 ! CALCULATE DIRECT SOIL EVAPORATION
2331 ! ----------------------------------------------------------------------
2332 IMPLICIT NONE
2333 REAL, INTENT(IN) :: ETP1,SMC,BEXP,DKSAT,DWSAT,FXEXP, &
2334 SHDFAC,SMCDRY,SMCMAX,ZSOIL,SMCREF,SMCWLT
2335 REAL, INTENT(OUT):: EDIR
2336 REAL :: FX, SRATIO
2337
2338
2339 ! ----------------------------------------------------------------------
2340 ! DIRECT EVAP A FUNCTION OF RELATIVE SOIL MOISTURE AVAILABILITY, LINEAR
2341 ! WHEN FXEXP=1.
2342 ! ----------------------------------------------------------------------
2343 ! ----------------------------------------------------------------------
2344 ! FX > 1 REPRESENTS DEMAND CONTROL
2345 ! FX < 1 REPRESENTS FLUX CONTROL
2346 ! ----------------------------------------------------------------------
2347
2348 SRATIO = (SMC - SMCDRY) / (SMCMAX - SMCDRY)
2349 IF (SRATIO > 0.) THEN
2350 FX = SRATIO**FXEXP
2351 FX = MAX ( MIN ( FX, 1. ) ,0. )
2352 ELSE
2353 FX = 0.
2354 ENDIF
2355
2356 ! ----------------------------------------------------------------------
2357 ! ALLOW FOR THE DIRECT-EVAP-REDUCING EFFECT OF SHADE
2358 ! ----------------------------------------------------------------------
2359 EDIR = FX * ( 1.0- SHDFAC ) * ETP1
2360
2361 ! ----------------------------------------------------------------------
2362 END SUBROUTINE DEVAP
2363 ! ----------------------------------------------------------------------
2364
2365 SUBROUTINE EVAPO (ETA1,SMC,NSOIL,CMC,ETP1,DT,ZSOIL, &
2366 SH2O, &
2367 SMCMAX,BEXP,PC,SMCWLT,DKSAT,DWSAT, &
2368 SMCREF,SHDFAC,CMCMAX, &
2369 SMCDRY,CFACTR, &
2370 EDIR,EC,ET,ETT,SFCTMP,Q2,NROOT,RTDIS,FXEXP)
2371
2372 ! ----------------------------------------------------------------------
2373 ! SUBROUTINE EVAPO
2374 ! ----------------------------------------------------------------------
2375 ! CALCULATE SOIL MOISTURE FLUX. THE SOIL MOISTURE CONTENT (SMC - A PER
2376 ! UNIT VOLUME MEASUREMENT) IS A DEPENDENT VARIABLE THAT IS UPDATED WITH
2377 ! PROGNOSTIC EQNS. THE CANOPY MOISTURE CONTENT (CMC) IS ALSO UPDATED.
2378 ! FROZEN GROUND VERSION: NEW STATES ADDED: SH2O, AND FROZEN GROUND
2379 ! CORRECTION FACTOR, FRZFACT AND PARAMETER SLOPE.
2380 ! ----------------------------------------------------------------------
2381 IMPLICIT NONE
2382 INTEGER, INTENT(IN) :: NSOIL, NROOT
2383 INTEGER :: I,K
2384 REAL, INTENT(IN) :: BEXP, CFACTR,CMC,CMCMAX,DKSAT, &
2385 DT,DWSAT,ETP1,FXEXP,PC,Q2,SFCTMP, &
2386 SHDFAC,SMCDRY,SMCMAX,SMCREF,SMCWLT
2387 REAL, INTENT(OUT) :: EC,EDIR,ETA1,ETT
2388 REAL :: CMC2MS
2389 REAL,DIMENSION(1:NSOIL), INTENT(IN) :: RTDIS, SMC, SH2O, ZSOIL
2390 REAL,DIMENSION(1:NSOIL), INTENT(OUT) :: ET
2391
2392 ! ----------------------------------------------------------------------
2393 ! EXECUTABLE CODE BEGINS HERE IF THE POTENTIAL EVAPOTRANSPIRATION IS
2394 ! GREATER THAN ZERO.
2395 ! ----------------------------------------------------------------------
2396 EDIR = 0.
2397 EC = 0.
2398 ETT = 0.
2399 DO K = 1,NSOIL
2400 ET (K) = 0.
2401 END DO
2402
2403 ! ----------------------------------------------------------------------
2404 ! RETRIEVE DIRECT EVAPORATION FROM SOIL SURFACE. CALL THIS FUNCTION
2405 ! ONLY IF VEG COVER NOT COMPLETE.
2406 ! FROZEN GROUND VERSION: SH2O STATES REPLACE SMC STATES.
2407 ! ----------------------------------------------------------------------
2408 IF (ETP1 > 0.0) THEN
2409 IF (SHDFAC < 1.) THEN
2410 CALL DEVAP (EDIR,ETP1,SMC (1),ZSOIL (1),SHDFAC,SMCMAX, &
2411 BEXP,DKSAT,DWSAT,SMCDRY,SMCREF,SMCWLT,FXEXP)
2412 END IF
2413 ! ----------------------------------------------------------------------
2414 ! INITIALIZE PLANT TOTAL TRANSPIRATION, RETRIEVE PLANT TRANSPIRATION,
2415 ! AND ACCUMULATE IT FOR ALL SOIL LAYERS.
2416 ! ----------------------------------------------------------------------
2417
2418 IF (SHDFAC > 0.0) THEN
2419 CALL TRANSP (ET,NSOIL,ETP1,SH2O,CMC,ZSOIL,SHDFAC,SMCWLT, &
2420 CMCMAX,PC,CFACTR,SMCREF,SFCTMP,Q2,NROOT,RTDIS)
2421 DO K = 1,NSOIL
2422 ETT = ETT + ET ( K )
2423 END DO
2424 ! ----------------------------------------------------------------------
2425 ! CALCULATE CANOPY EVAPORATION.
2426 ! IF STATEMENTS TO AVOID TANGENT LINEAR PROBLEMS NEAR CMC=0.0.
2427 ! ----------------------------------------------------------------------
2428 IF (CMC > 0.0) THEN
2429 EC = SHDFAC * ( ( CMC / CMCMAX ) ** CFACTR ) * ETP1
2430 ELSE
2431 EC = 0.0
2432 END IF
2433 ! ----------------------------------------------------------------------
2434 ! EC SHOULD BE LIMITED BY THE TOTAL AMOUNT OF AVAILABLE WATER ON THE
2435 ! CANOPY. -F.CHEN, 18-OCT-1994
2436 ! ----------------------------------------------------------------------
2437 CMC2MS = CMC / DT
2438 EC = MIN ( CMC2MS, EC )
2439 END IF
2440 END IF
2441 ! ----------------------------------------------------------------------
2442 ! TOTAL UP EVAP AND TRANSP TYPES TO OBTAIN ACTUAL EVAPOTRANSP
2443 ! ----------------------------------------------------------------------
2444 ETA1 = EDIR + ETT + EC
2445
2446 ! ----------------------------------------------------------------------
2447 END SUBROUTINE EVAPO
2448 ! ----------------------------------------------------------------------
2449
2450 SUBROUTINE FRH2O (FREE,TKELV,SMC,SH2O,SMCMAX,BEXP,PSIS)
2451
2452 ! ----------------------------------------------------------------------
2453 ! SUBROUTINE FRH2O
2454 ! ----------------------------------------------------------------------
2455 ! CALCULATE AMOUNT OF SUPERCOOLED LIQUID SOIL WATER CONTENT IF
2456 ! TEMPERATURE IS BELOW 273.15K (T0). REQUIRES NEWTON-TYPE ITERATION TO
2457 ! SOLVE THE NONLINEAR IMPLICIT EQUATION GIVEN IN EQN 17 OF KOREN ET AL
2458 ! (1999, JGR, VOL 104(D16), 19569-19585).
2459 ! ----------------------------------------------------------------------
2460 ! NEW VERSION (JUNE 2001): MUCH FASTER AND MORE ACCURATE NEWTON
2461 ! ITERATION ACHIEVED BY FIRST TAKING LOG OF EQN CITED ABOVE -- LESS THAN
2462 ! 4 (TYPICALLY 1 OR 2) ITERATIONS ACHIEVES CONVERGENCE. ALSO, EXPLICIT
2463 ! 1-STEP SOLUTION OPTION FOR SPECIAL CASE OF PARAMETER CK=0, WHICH
2464 ! REDUCES THE ORIGINAL IMPLICIT EQUATION TO A SIMPLER EXPLICIT FORM,
2465 ! KNOWN AS THE "FLERCHINGER EQN". IMPROVED HANDLING OF SOLUTION IN THE
2466 ! LIMIT OF FREEZING POINT TEMPERATURE T0.
2467 ! ----------------------------------------------------------------------
2468 ! INPUT:
2469
2470 ! TKELV.........TEMPERATURE (Kelvin)
2471 ! SMC...........TOTAL SOIL MOISTURE CONTENT (VOLUMETRIC)
2472 ! SH2O..........LIQUID SOIL MOISTURE CONTENT (VOLUMETRIC)
2473 ! SMCMAX........SATURATION SOIL MOISTURE CONTENT (FROM REDPRM)
2474 ! B.............SOIL TYPE "B" PARAMETER (FROM REDPRM)
2475 ! PSIS..........SATURATED SOIL MATRIC POTENTIAL (FROM REDPRM)
2476
2477 ! OUTPUT:
2478 ! FRH2O.........SUPERCOOLED LIQUID WATER CONTENT
2479 ! FREE..........SUPERCOOLED LIQUID WATER CONTENT
2480 ! ----------------------------------------------------------------------
2481 IMPLICIT NONE
2482 REAL, INTENT(IN) :: BEXP,PSIS,SH2O,SMC,SMCMAX,TKELV
2483 REAL, INTENT(OUT) :: FREE
2484 REAL :: BX,DENOM,DF,DSWL,FK,SWL,SWLK
2485 INTEGER :: NLOG,KCOUNT
2486 ! PARAMETER(CK = 0.0)
2487 REAL, PARAMETER :: CK = 8.0, BLIM = 5.5, ERROR = 0.005, &
2488 HLICE = 3.335E5, GS = 9.81,DICE = 920.0, &
2489 DH2O = 1000.0, T0 = 273.15
2490
2491 ! ----------------------------------------------------------------------
2492 ! LIMITS ON PARAMETER B: B < 5.5 (use parameter BLIM)
2493 ! SIMULATIONS SHOWED IF B > 5.5 UNFROZEN WATER CONTENT IS
2494 ! NON-REALISTICALLY HIGH AT VERY LOW TEMPERATURES.
2495 ! ----------------------------------------------------------------------
2496 BX = BEXP
2497
2498 ! ----------------------------------------------------------------------
2499 ! INITIALIZING ITERATIONS COUNTER AND ITERATIVE SOLUTION FLAG.
2500 ! ----------------------------------------------------------------------
2501 IF (BEXP > BLIM) BX = BLIM
2502 NLOG = 0
2503
2504 ! ----------------------------------------------------------------------
2505 ! IF TEMPERATURE NOT SIGNIFICANTLY BELOW FREEZING (T0), SH2O = SMC
2506 ! ----------------------------------------------------------------------
2507 KCOUNT = 0
2508 ! FRH2O = SMC
2509 IF (TKELV > (T0- 1.E-3)) THEN
2510 FREE = SMC
2511 ELSE
2512
2513 ! ----------------------------------------------------------------------
2514 ! OPTION 1: ITERATED SOLUTION FOR NONZERO CK
2515 ! IN KOREN ET AL, JGR, 1999, EQN 17
2516 ! ----------------------------------------------------------------------
2517 ! INITIAL GUESS FOR SWL (frozen content)
2518 ! ----------------------------------------------------------------------
2519 IF (CK /= 0.0) THEN
2520 SWL = SMC - SH2O
2521 ! ----------------------------------------------------------------------
2522 ! KEEP WITHIN BOUNDS.
2523 ! ----------------------------------------------------------------------
2524 IF (SWL > (SMC -0.02)) SWL = SMC -0.02
2525
2526 ! ----------------------------------------------------------------------
2527 ! START OF ITERATIONS
2528 ! ----------------------------------------------------------------------
2529 IF (SWL < 0.) SWL = 0.
2530 1001 Continue
2531 IF (.NOT.( (NLOG < 10) .AND. (KCOUNT == 0))) goto 1002
2532 NLOG = NLOG +1
2533 DF = ALOG ( ( PSIS * GS / HLICE ) * ( ( 1. + CK * SWL )**2.) * &
2534 ( SMCMAX / (SMC - SWL) )** BX) - ALOG ( - ( &
2535 TKELV - T0)/ TKELV)
2536 DENOM = 2. * CK / ( 1. + CK * SWL ) + BX / ( SMC - SWL )
2537 SWLK = SWL - DF / DENOM
2538 ! ----------------------------------------------------------------------
2539 ! BOUNDS USEFUL FOR MATHEMATICAL SOLUTION.
2540 ! ----------------------------------------------------------------------
2541 IF (SWLK > (SMC -0.02)) SWLK = SMC - 0.02
2542 IF (SWLK < 0.) SWLK = 0.
2543
2544 ! ----------------------------------------------------------------------
2545 ! MATHEMATICAL SOLUTION BOUNDS APPLIED.
2546 ! ----------------------------------------------------------------------
2547 DSWL = ABS (SWLK - SWL)
2548
2549 ! ----------------------------------------------------------------------
2550 ! IF MORE THAN 10 ITERATIONS, USE EXPLICIT METHOD (CK=0 APPROX.)
2551 ! WHEN DSWL LESS OR EQ. ERROR, NO MORE ITERATIONS REQUIRED.
2552 ! ----------------------------------------------------------------------
2553 SWL = SWLK
2554 IF ( DSWL <= ERROR ) THEN
2555 KCOUNT = KCOUNT +1
2556 END IF
2557 ! ----------------------------------------------------------------------
2558 ! END OF ITERATIONS
2559 ! ----------------------------------------------------------------------
2560 ! BOUNDS APPLIED WITHIN DO-BLOCK ARE VALID FOR PHYSICAL SOLUTION.
2561 ! ----------------------------------------------------------------------
2562 ! FRH2O = SMC - SWL
2563 goto 1001
2564 1002 continue
2565 FREE = SMC - SWL
2566 END IF
2567 ! ----------------------------------------------------------------------
2568 ! END OPTION 1
2569 ! ----------------------------------------------------------------------
2570 ! ----------------------------------------------------------------------
2571 ! OPTION 2: EXPLICIT SOLUTION FOR FLERCHINGER EQ. i.e. CK=0
2572 ! IN KOREN ET AL., JGR, 1999, EQN 17
2573 ! APPLY PHYSICAL BOUNDS TO FLERCHINGER SOLUTION
2574 ! ----------------------------------------------------------------------
2575 IF (KCOUNT == 0) THEN
2576 PRINT *,'Flerchinger USEd in NEW version. Iterations=',NLOG
2577 FK = ( ( (HLICE / (GS * ( - PSIS)))* &
2578 ( (TKELV - T0)/ TKELV))** ( -1/ BX))* SMCMAX
2579 ! FRH2O = MIN (FK, SMC)
2580 IF (FK < 0.02) FK = 0.02
2581 FREE = MIN (FK, SMC)
2582 ! ----------------------------------------------------------------------
2583 ! END OPTION 2
2584 ! ----------------------------------------------------------------------
2585 END IF
2586 END IF
2587 ! ----------------------------------------------------------------------
2588 END SUBROUTINE FRH2O
2589 ! ----------------------------------------------------------------------
2590
2591 SUBROUTINE HRT (RHSTS,STC,SMC,SMCMAX,NSOIL,ZSOIL,YY,ZZ1, &
2592 TBOT,ZBOT,PSISAT,SH2O,DT,BEXP, &
2593 F1,DF1,QUARTZ,CSOIL,AI,BI,CI,VEGTYP)
2594
2595 ! ----------------------------------------------------------------------
2596 ! SUBROUTINE HRT
2597 ! ----------------------------------------------------------------------
2598 ! CALCULATE THE RIGHT HAND SIDE OF THE TIME TENDENCY TERM OF THE SOIL
2599 ! THERMAL DIFFUSION EQUATION. ALSO TO COMPUTE ( PREPARE ) THE MATRIX
2600 ! COEFFICIENTS FOR THE TRI-DIAGONAL MATRIX OF THE IMPLICIT TIME SCHEME.
2601 ! ----------------------------------------------------------------------
2602 IMPLICIT NONE
2603 LOGICAL :: ITAVG
2604 INTEGER, INTENT(IN) :: NSOIL, VEGTYP
2605 INTEGER :: I, K
2606
2607 REAL, INTENT(IN) :: BEXP, CSOIL, DF1, DT,F1,PSISAT,QUARTZ, &
2608 SMCMAX ,TBOT,YY,ZZ1, ZBOT
2609 REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SMC,STC,ZSOIL
2610 REAL, DIMENSION(1:NSOIL), INTENT(INOUT):: SH2O
2611 REAL, DIMENSION(1:NSOIL), INTENT(OUT) :: RHSTS
2612 REAL, DIMENSION(1:NSOIL), INTENT(OUT) :: AI, BI,CI
2613 REAL :: DDZ, DDZ2, DENOM, DF1N, DF1K, DTSDZ, &
2614 DTSDZ2,HCPCT,QTOT,SSOIL,SICE,TAVG,TBK, &
2615 TBK1,TSNSR,TSURF,CSOIL_LOC
2616 REAL, PARAMETER :: T0 = 273.15, CAIR = 1004.0, CICE = 2.106E6,&
2617 CH2O = 4.2E6
2618
2619
2620 !urban
2621 IF(VEGTYP==1) then
2622 CSOIL_LOC=3.0E6
2623 ELSE
2624 CSOIL_LOC=CSOIL
2625 ENDIF
2626
2627 ! ----------------------------------------------------------------------
2628 ! INITIALIZE LOGICAL FOR SOIL LAYER TEMPERATURE AVERAGING.
2629 ! ----------------------------------------------------------------------
2630 ITAVG = .TRUE.
2631 ! ----------------------------------------------------------------------
2632 ! BEGIN SECTION FOR TOP SOIL LAYER
2633 ! ----------------------------------------------------------------------
2634 ! CALC THE HEAT CAPACITY OF THE TOP SOIL LAYER
2635 ! ----------------------------------------------------------------------
2636 HCPCT = SH2O (1)* CH2O + (1.0- SMCMAX)* CSOIL_LOC + (SMCMAX - SMC (1))&
2637 * CAIR &
2638 + ( SMC (1) - SH2O (1) )* CICE
2639
2640 ! ----------------------------------------------------------------------
2641 ! CALC THE MATRIX COEFFICIENTS AI, BI, AND CI FOR THE TOP LAYER
2642 ! ----------------------------------------------------------------------
2643 DDZ = 1.0 / ( -0.5 * ZSOIL (2) )
2644 AI (1) = 0.0
2645 CI (1) = (DF1 * DDZ) / (ZSOIL (1) * HCPCT)
2646
2647 ! ----------------------------------------------------------------------
2648 ! CALCULATE THE VERTICAL SOIL TEMP GRADIENT BTWN THE 1ST AND 2ND SOIL
2649 ! LAYERS. THEN CALCULATE THE SUBSURFACE HEAT FLUX. USE THE TEMP
2650 ! GRADIENT AND SUBSFC HEAT FLUX TO CALC "RIGHT-HAND SIDE TENDENCY
2651 ! TERMS", OR "RHSTS", FOR TOP SOIL LAYER.
2652 ! ----------------------------------------------------------------------
2653 BI (1) = - CI (1) + DF1 / (0.5 * ZSOIL (1) * ZSOIL (1)* HCPCT * &
2654 ZZ1)
2655 DTSDZ = (STC (1) - STC (2)) / ( -0.5 * ZSOIL (2))
2656 SSOIL = DF1 * (STC (1) - YY) / (0.5 * ZSOIL (1) * ZZ1)
2657 ! RHSTS(1) = (DF1 * DTSDZ - SSOIL) / (ZSOIL(1) * HCPCT)
2658 DENOM = (ZSOIL (1) * HCPCT)
2659
2660 ! ----------------------------------------------------------------------
2661 ! NEXT CAPTURE THE VERTICAL DIFFERENCE OF THE HEAT FLUX AT TOP AND
2662 ! BOTTOM OF FIRST SOIL LAYER FOR USE IN HEAT FLUX CONSTRAINT APPLIED TO
2663 ! POTENTIAL SOIL FREEZING/THAWING IN ROUTINE SNKSRC.
2664 ! ----------------------------------------------------------------------
2665 ! QTOT = SSOIL - DF1*DTSDZ
2666 RHSTS (1) = (DF1 * DTSDZ - SSOIL) / DENOM
2667
2668 ! ----------------------------------------------------------------------
2669 ! CALCULATE FROZEN WATER CONTENT IN 1ST SOIL LAYER.
2670 ! ----------------------------------------------------------------------
2671 QTOT = -1.0* RHSTS (1)* DENOM
2672
2673 ! ----------------------------------------------------------------------
2674 ! IF TEMPERATURE AVERAGING INVOKED (ITAVG=TRUE; ELSE SKIP):
2675 ! SET TEMP "TSURF" AT TOP OF SOIL COLUMN (FOR USE IN FREEZING SOIL
2676 ! PHYSICS LATER IN FUNCTION SUBROUTINE SNKSRC). IF SNOWPACK CONTENT IS
2677 ! ZERO, THEN TSURF EXPRESSION BELOW GIVES TSURF = SKIN TEMP. IF
2678 ! SNOWPACK IS NONZERO (HENCE ARGUMENT ZZ1=1), THEN TSURF EXPRESSION
2679 ! BELOW YIELDS SOIL COLUMN TOP TEMPERATURE UNDER SNOWPACK. THEN
2680 ! CALCULATE TEMPERATURE AT BOTTOM INTERFACE OF 1ST SOIL LAYER FOR USE
2681 ! LATER IN FUNCTION SUBROUTINE SNKSRC
2682 ! ----------------------------------------------------------------------
2683 SICE = SMC (1) - SH2O (1)
2684 IF (ITAVG) THEN
2685 TSURF = (YY + (ZZ1-1) * STC (1)) / ZZ1
2686 ! ----------------------------------------------------------------------
2687 ! IF FROZEN WATER PRESENT OR ANY OF LAYER-1 MID-POINT OR BOUNDING
2688 ! INTERFACE TEMPERATURES BELOW FREEZING, THEN CALL SNKSRC TO
2689 ! COMPUTE HEAT SOURCE/SINK (AND CHANGE IN FROZEN WATER CONTENT)
2690 ! DUE TO POSSIBLE SOIL WATER PHASE CHANGE
2691 ! ----------------------------------------------------------------------
2692 CALL TBND (STC (1),STC (2),ZSOIL,ZBOT,1,NSOIL,TBK)
2693 IF ( (SICE > 0.) .OR. (STC (1) < T0) .OR. &
2694 (TSURF < T0) .OR. (TBK < T0) ) THEN
2695 ! TSNSR = SNKSRC (TAVG,SMC(1),SH2O(1),
2696 CALL TMPAVG (TAVG,TSURF,STC (1),TBK,ZSOIL,NSOIL,1)
2697 CALL SNKSRC (TSNSR,TAVG,SMC (1),SH2O (1), &
2698 ZSOIL,NSOIL,SMCMAX,PSISAT,BEXP,DT,1,QTOT)
2699 ! RHSTS(1) = RHSTS(1) - TSNSR / ( ZSOIL(1) * HCPCT )
2700 RHSTS (1) = RHSTS (1) - TSNSR / DENOM
2701 END IF
2702 ELSE
2703 ! TSNSR = SNKSRC (STC(1),SMC(1),SH2O(1),
2704 IF ( (SICE > 0.) .OR. (STC (1) < T0) ) THEN
2705 CALL SNKSRC (TSNSR,STC (1),SMC (1),SH2O (1), &
2706 ZSOIL,NSOIL,SMCMAX,PSISAT,BEXP,DT,1,QTOT)
2707 ! RHSTS(1) = RHSTS(1) - TSNSR / ( ZSOIL(1) * HCPCT )
2708 RHSTS (1) = RHSTS (1) - TSNSR / DENOM
2709 END IF
2710 ! ----------------------------------------------------------------------
2711 ! THIS ENDS SECTION FOR TOP SOIL LAYER.
2712 ! ----------------------------------------------------------------------
2713 END IF
2714
2715 ! INITIALIZE DDZ2
2716 ! ----------------------------------------------------------------------
2717
2718 DDZ2 = 0.0
2719 DF1K = DF1
2720
2721 ! ----------------------------------------------------------------------
2722 ! LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABOVE PROCESS
2723 ! (EXCEPT SUBSFC OR "GROUND" HEAT FLUX NOT REPEATED IN LOWER LAYERS)
2724 ! ----------------------------------------------------------------------
2725 ! CALCULATE HEAT CAPACITY FOR THIS SOIL LAYER.
2726 ! ----------------------------------------------------------------------
2727 DO K = 2,NSOIL
2728 HCPCT = SH2O (K)* CH2O + (1.0- SMCMAX)* CSOIL_LOC + (SMCMAX - SMC ( &
2729 K))* CAIR + ( SMC (K) - SH2O (K) )* CICE
2730 ! ----------------------------------------------------------------------
2731 ! THIS SECTION FOR LAYER 2 OR GREATER, BUT NOT LAST LAYER.
2732 ! ----------------------------------------------------------------------
2733 ! CALCULATE THERMAL DIFFUSIVITY FOR THIS LAYER.
2734 ! ----------------------------------------------------------------------
2735 IF (K /= NSOIL) THEN
2736
2737 ! ----------------------------------------------------------------------
2738 ! CALC THE VERTICAL SOIL TEMP GRADIENT THRU THIS LAYER
2739 ! ----------------------------------------------------------------------
2740 CALL TDFCND (DF1N,SMC (K),QUARTZ,SMCMAX,SH2O (K))
2741
2742 !urban
2743 IF(VEGTYP==1) DF1N = 3.24
2744
2745 DENOM = 0.5 * ( ZSOIL (K -1) - ZSOIL (K +1) )
2746
2747 ! ----------------------------------------------------------------------
2748 ! CALC THE MATRIX COEF, CI, AFTER CALC'NG ITS PARTIAL PRODUCT
2749 ! ----------------------------------------------------------------------
2750 DTSDZ2 = ( STC (K) - STC (K +1) ) / DENOM
2751 DDZ2 = 2. / (ZSOIL (K -1) - ZSOIL (K +1))
2752
2753 ! ----------------------------------------------------------------------
2754 ! IF TEMPERATURE AVERAGING INVOKED (ITAVG=TRUE; ELSE SKIP): CALCULATE
2755 ! TEMP AT BOTTOM OF LAYER.
2756 ! ----------------------------------------------------------------------
2757 CI (K) = - DF1N * DDZ2 / ( (ZSOIL (K -1) - ZSOIL (K)) * &
2758 HCPCT)
2759 IF (ITAVG) THEN
2760 CALL TBND (STC (K),STC (K +1),ZSOIL,ZBOT,K,NSOIL,TBK1)
2761 END IF
2762
2763 ELSE
2764 ! ----------------------------------------------------------------------
2765 ! SPECIAL CASE OF BOTTOM SOIL LAYER: CALCULATE THERMAL DIFFUSIVITY FOR
2766 ! BOTTOM LAYER.
2767 ! ----------------------------------------------------------------------
2768
2769 ! ----------------------------------------------------------------------
2770 ! CALC THE VERTICAL SOIL TEMP GRADIENT THRU BOTTOM LAYER.
2771 ! ----------------------------------------------------------------------
2772 CALL TDFCND (DF1N,SMC (K),QUARTZ,SMCMAX,SH2O (K))
2773
2774
2775 !urban
2776 IF(VEGTYP==1) DF1N = 3.24
2777
2778 DENOM = .5 * (ZSOIL (K -1) + ZSOIL (K)) - ZBOT
2779
2780 ! ----------------------------------------------------------------------
2781 ! SET MATRIX COEF, CI TO ZERO IF BOTTOM LAYER.
2782 ! ----------------------------------------------------------------------
2783 DTSDZ2 = (STC (K) - TBOT) / DENOM
2784
2785 ! ----------------------------------------------------------------------
2786 ! IF TEMPERATURE AVERAGING INVOKED (ITAVG=TRUE; ELSE SKIP): CALCULATE
2787 ! TEMP AT BOTTOM OF LAST LAYER.
2788 ! ----------------------------------------------------------------------
2789 CI (K) = 0.
2790 IF (ITAVG) THEN
2791 CALL TBND (STC (K),TBOT,ZSOIL,ZBOT,K,NSOIL,TBK1)
2792 END IF
2793 ! ----------------------------------------------------------------------
2794 ! THIS ENDS SPECIAL LOOP FOR BOTTOM LAYER.
2795 END IF
2796 ! ----------------------------------------------------------------------
2797 ! CALCULATE RHSTS FOR THIS LAYER AFTER CALC'NG A PARTIAL PRODUCT.
2798 ! ----------------------------------------------------------------------
2799 DENOM = ( ZSOIL (K) - ZSOIL (K -1) ) * HCPCT
2800 RHSTS (K) = ( DF1N * DTSDZ2- DF1K * DTSDZ ) / DENOM
2801 QTOT = -1.0* DENOM * RHSTS (K)
2802
2803 SICE = SMC (K) - SH2O (K)
2804 IF (ITAVG) THEN
2805 CALL TMPAVG (TAVG,TBK,STC (K),TBK1,ZSOIL,NSOIL,K)
2806 ! TSNSR = SNKSRC(TAVG,SMC(K),SH2O(K),ZSOIL,NSOIL,
2807 IF ( (SICE > 0.) .OR. (STC (K) < T0) .OR. &
2808 (TBK .lt. T0) .OR. (TBK1 .lt. T0) ) THEN
2809 CALL SNKSRC (TSNSR,TAVG,SMC (K),SH2O (K),ZSOIL,NSOIL, &
2810 SMCMAX,PSISAT,BEXP,DT,K,QTOT)
2811 RHSTS (K) = RHSTS (K) - TSNSR / DENOM
2812 END IF
2813 ELSE
2814 ! TSNSR = SNKSRC(STC(K),SMC(K),SH2O(K),ZSOIL,NSOIL,
2815 IF ( (SICE > 0.) .OR. (STC (K) < T0) ) THEN
2816 CALL SNKSRC (TSNSR,STC (K),SMC (K),SH2O (K),ZSOIL,NSOIL, &
2817 SMCMAX,PSISAT,BEXP,DT,K,QTOT)
2818 RHSTS (K) = RHSTS (K) - TSNSR / DENOM
2819 END IF
2820 END IF
2821
2822 ! ----------------------------------------------------------------------
2823 ! CALC MATRIX COEFS, AI, AND BI FOR THIS LAYER.
2824 ! ----------------------------------------------------------------------
2825 AI (K) = - DF1K * DDZ / ( (ZSOIL (K -1) - ZSOIL (K)) * HCPCT)
2826
2827 ! ----------------------------------------------------------------------
2828 ! RESET VALUES OF DF1, DTSDZ, DDZ, AND TBK FOR LOOP TO NEXT SOIL LAYER.
2829 ! ----------------------------------------------------------------------
2830 BI (K) = - (AI (K) + CI (K))
2831 TBK = TBK1
2832 DF1K = DF1N
2833 DTSDZ = DTSDZ2
2834 DDZ = DDZ2
2835 END DO
2836 ! ----------------------------------------------------------------------
2837 END SUBROUTINE HRT
2838 ! ----------------------------------------------------------------------
2839
2840 SUBROUTINE HRTICE (RHSTS,STC,NSOIL,ZSOIL,YY,ZZ1,DF1,AI,BI,CI)
2841
2842 ! ----------------------------------------------------------------------
2843 ! SUBROUTINE HRTICE
2844 ! ----------------------------------------------------------------------
2845 ! CALCULATE THE RIGHT HAND SIDE OF THE TIME TENDENCY TERM OF THE SOIL
2846 ! THERMAL DIFFUSION EQUATION IN THE CASE OF SEA-ICE PACK. ALSO TO
2847 ! COMPUTE (PREPARE) THE MATRIX COEFFICIENTS FOR THE TRI-DIAGONAL MATRIX
2848 ! OF THE IMPLICIT TIME SCHEME.
2849 ! ----------------------------------------------------------------------
2850 IMPLICIT NONE
2851
2852
2853 INTEGER, INTENT(IN) :: NSOIL
2854 INTEGER :: K
2855
2856 REAL, INTENT(IN) :: DF1,YY,ZZ1
2857 REAL, DIMENSION(1:NSOIL), INTENT(OUT):: AI, BI,CI
2858 REAL, DIMENSION(1:NSOIL), INTENT(IN) :: STC, ZSOIL
2859 REAL, DIMENSION(1:NSOIL), INTENT(OUT):: RHSTS
2860 REAL :: DDZ,DDZ2,DENOM,DTSDZ,DTSDZ2,SSOIL, &
2861 ZBOT,TBOT
2862 REAL, PARAMETER :: HCPCT = 1.72396E+6
2863
2864 ! ----------------------------------------------------------------------
2865 ! SET A NOMINAL UNIVERSAL VALUE OF THE SEA-ICE SPECIFIC HEAT CAPACITY,
2866 ! HCPCT = 1880.0*917.0.
2867 ! ----------------------------------------------------------------------
2868 DATA TBOT /271.16/
2869
2870 ! ----------------------------------------------------------------------
2871 ! THE INPUT ARGUMENT DF1 IS A UNIVERSALLY CONSTANT VALUE OF SEA-ICE
2872 ! THERMAL DIFFUSIVITY, SET IN ROUTINE SNOPAC AS DF1 = 2.2.
2873 ! ----------------------------------------------------------------------
2874 ! SET ICE PACK DEPTH. USE TBOT AS ICE PACK LOWER BOUNDARY TEMPERATURE
2875 ! (THAT OF UNFROZEN SEA WATER AT BOTTOM OF SEA ICE PACK). ASSUME ICE
2876 ! PACK IS OF N=NSOIL LAYERS SPANNING A UNIFORM CONSTANT ICE PACK
2877 ! THICKNESS AS DEFINED BY ZSOIL(NSOIL) IN ROUTINE SFLX.
2878 ! ----------------------------------------------------------------------
2879 ! ----------------------------------------------------------------------
2880 ! CALC THE MATRIX COEFFICIENTS AI, BI, AND CI FOR THE TOP LAYER
2881 ! ----------------------------------------------------------------------
2882 ZBOT = ZSOIL (NSOIL)
2883 DDZ = 1.0 / ( -0.5 * ZSOIL (2) )
2884 AI (1) = 0.0
2885 CI (1) = (DF1 * DDZ) / (ZSOIL (1) * HCPCT)
2886
2887 ! ----------------------------------------------------------------------
2888 ! CALC THE VERTICAL SOIL TEMP GRADIENT BTWN THE TOP AND 2ND SOIL LAYERS.
2889 ! RECALC/ADJUST THE SOIL HEAT FLUX. USE THE GRADIENT AND FLUX TO CALC
2890 ! RHSTS FOR THE TOP SOIL LAYER.
2891 ! ----------------------------------------------------------------------
2892 BI (1) = - CI (1) + DF1/ (0.5 * ZSOIL (1) * ZSOIL (1) * HCPCT * &
2893 ZZ1)
2894 DTSDZ = ( STC (1) - STC (2) ) / ( -0.5 * ZSOIL (2) )
2895 SSOIL = DF1 * ( STC (1) - YY ) / ( 0.5 * ZSOIL (1) * ZZ1 )
2896
2897 ! ----------------------------------------------------------------------
2898 ! INITIALIZE DDZ2
2899 ! ----------------------------------------------------------------------
2900 RHSTS (1) = ( DF1 * DTSDZ - SSOIL ) / ( ZSOIL (1) * HCPCT )
2901
2902 ! ----------------------------------------------------------------------
2903 ! LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABOVE PROCESS
2904 ! ----------------------------------------------------------------------
2905 DDZ2 = 0.0
2906 DO K = 2,NSOIL
2907
2908 ! ----------------------------------------------------------------------
2909 ! CALC THE VERTICAL SOIL TEMP GRADIENT THRU THIS LAYER.
2910 ! ----------------------------------------------------------------------
2911 IF (K /= NSOIL) THEN
2912 DENOM = 0.5 * ( ZSOIL (K -1) - ZSOIL (K +1) )
2913
2914 ! ----------------------------------------------------------------------
2915 ! CALC THE MATRIX COEF, CI, AFTER CALC'NG ITS PARTIAL PRODUCT.
2916 ! ----------------------------------------------------------------------
2917 DTSDZ2 = ( STC (K) - STC (K +1) ) / DENOM
2918 DDZ2 = 2. / (ZSOIL (K -1) - ZSOIL (K +1))
2919 CI (K) = - DF1 * DDZ2 / ( (ZSOIL (K -1) - ZSOIL (K))*HCPCT)
2920
2921 ! ----------------------------------------------------------------------
2922 ! CALC THE VERTICAL SOIL TEMP GRADIENT THRU THE LOWEST LAYER.
2923 ! ----------------------------------------------------------------------
2924 ELSE
2925
2926 ! ----------------------------------------------------------------------
2927 ! SET MATRIX COEF, CI TO ZERO.
2928 ! ----------------------------------------------------------------------
2929 DTSDZ2 = (STC (K) - TBOT)/ (.5 * (ZSOIL (K -1) + ZSOIL (K)) &
2930 - ZBOT)
2931 CI (K) = 0.
2932 ! ----------------------------------------------------------------------
2933 ! CALC RHSTS FOR THIS LAYER AFTER CALC'NG A PARTIAL PRODUCT.
2934 ! ----------------------------------------------------------------------
2935 END IF
2936 DENOM = ( ZSOIL (K) - ZSOIL (K -1) ) * HCPCT
2937
2938 ! ----------------------------------------------------------------------
2939 ! CALC MATRIX COEFS, AI, AND BI FOR THIS LAYER.
2940 ! ----------------------------------------------------------------------
2941 RHSTS (K) = ( DF1 * DTSDZ2- DF1 * DTSDZ ) / DENOM
2942 AI (K) = - DF1 * DDZ / ( (ZSOIL (K -1) - ZSOIL (K)) * HCPCT)
2943
2944 ! ----------------------------------------------------------------------
2945 ! RESET VALUES OF DTSDZ AND DDZ FOR LOOP TO NEXT SOIL LYR.
2946 ! ----------------------------------------------------------------------
2947 BI (K) = - (AI (K) + CI (K))
2948 DTSDZ = DTSDZ2
2949 DDZ = DDZ2
2950 END DO
2951 ! ----------------------------------------------------------------------
2952 END SUBROUTINE HRTICE
2953 ! ----------------------------------------------------------------------
2954
2955 SUBROUTINE HSTEP (STCOUT,STCIN,RHSTS,DT,NSOIL,AI,BI,CI)
2956
2957 ! ----------------------------------------------------------------------
2958 ! SUBROUTINE HSTEP
2959 ! ----------------------------------------------------------------------
2960 ! CALCULATE/UPDATE THE SOIL TEMPERATURE FIELD.
2961 ! ----------------------------------------------------------------------
2962 IMPLICIT NONE
2963 INTEGER, INTENT(IN) :: NSOIL
2964 INTEGER :: K
2965
2966 REAL, DIMENSION(1:NSOIL), INTENT(IN):: STCIN
2967 REAL, DIMENSION(1:NSOIL), INTENT(OUT):: STCOUT
2968 REAL, DIMENSION(1:NSOIL), INTENT(INOUT):: RHSTS
2969 REAL, DIMENSION(1:NSOIL), INTENT(INOUT):: AI,BI,CI
2970 REAL, DIMENSION(1:NSOIL) :: RHSTSin
2971 REAL, DIMENSION(1:NSOIL) :: CIin
2972 REAL :: DT
2973
2974 ! ----------------------------------------------------------------------
2975 ! CREATE FINITE DIFFERENCE VALUES FOR USE IN ROSR12 ROUTINE
2976 ! ----------------------------------------------------------------------
2977 DO K = 1,NSOIL
2978 RHSTS (K) = RHSTS (K) * DT
2979 AI (K) = AI (K) * DT
2980 BI (K) = 1. + BI (K) * DT
2981 CI (K) = CI (K) * DT
2982 END DO
2983 ! ----------------------------------------------------------------------
2984 ! COPY VALUES FOR INPUT VARIABLES BEFORE CALL TO ROSR12
2985 ! ----------------------------------------------------------------------
2986 DO K = 1,NSOIL
2987 RHSTSin (K) = RHSTS (K)
2988 END DO
2989 DO K = 1,NSOIL
2990 CIin (K) = CI (K)
2991 END DO
2992 ! ----------------------------------------------------------------------
2993 ! SOLVE THE TRI-DIAGONAL MATRIX EQUATION
2994 ! ----------------------------------------------------------------------
2995 CALL ROSR12 (CI,AI,BI,CIin,RHSTSin,RHSTS,NSOIL)
2996 ! ----------------------------------------------------------------------
2997 ! CALC/UPDATE THE SOIL TEMPS USING MATRIX SOLUTION
2998 ! ----------------------------------------------------------------------
2999 DO K = 1,NSOIL
3000 STCOUT (K) = STCIN (K) + CI (K)
3001 END DO
3002 ! ----------------------------------------------------------------------
3003 END SUBROUTINE HSTEP
3004 ! ----------------------------------------------------------------------
3005
3006 SUBROUTINE NOPAC (ETP,ETA,PRCP,SMC,SMCMAX,SMCWLT, &
3007 SMCREF,SMCDRY,CMC,CMCMAX,NSOIL,DT,SHDFAC, &
3008 SBETA,Q2,T1,SFCTMP,T24,TH2,FDOWN,F1,EMISSI, &
3009 SSOIL, &
3010 STC,EPSCA,BEXP,PC,RCH,RR,CFACTR, &
3011 SH2O,SLOPE,KDT,FRZFACT,PSISAT,ZSOIL, &
3012 DKSAT,DWSAT,TBOT,ZBOT,RUNOFF1,RUNOFF2, &
3013 RUNOFF3,EDIR,EC,ET,ETT,NROOT,ICE,RTDIS, &
3014 QUARTZ,FXEXP,CSOIL, &
3015 BETA,DRIP,DEW,FLX1,FLX3,VEGTYP)
3016
3017 ! ----------------------------------------------------------------------
3018 ! SUBROUTINE NOPAC
3019 ! ----------------------------------------------------------------------
3020 ! CALCULATE SOIL MOISTURE AND HEAT FLUX VALUES AND UPDATE SOIL MOISTURE
3021 ! CONTENT AND SOIL HEAT CONTENT VALUES FOR THE CASE WHEN NO SNOW PACK IS
3022 ! PRESENT.
3023 ! ----------------------------------------------------------------------
3024 IMPLICIT NONE
3025
3026 INTEGER, INTENT(IN) :: ICE, NROOT,NSOIL,VEGTYP
3027 INTEGER :: K
3028
3029 REAL, INTENT(IN) :: BEXP,CFACTR, CMCMAX,CSOIL,DKSAT,DT,DWSAT, &
3030 EPSCA,ETP,FDOWN,F1,FXEXP,FRZFACT,KDT,PC, &
3031 PRCP,PSISAT,Q2,QUARTZ,RCH,RR,SBETA,SFCTMP,&
3032 SHDFAC,SLOPE,SMCDRY,SMCMAX,SMCREF,SMCWLT, &
3033 T24,TBOT,TH2,ZBOT,EMISSI
3034 REAL, INTENT(INOUT) :: CMC,BETA,T1
3035 REAL, INTENT(OUT) :: DEW,DRIP,EC,EDIR,ETA,ETT,FLX1,FLX3, &
3036 RUNOFF1,RUNOFF2,RUNOFF3,SSOIL
3037 REAL, DIMENSION(1:NSOIL),INTENT(IN) :: RTDIS,ZSOIL
3038 REAL, DIMENSION(1:NSOIL),INTENT(OUT) :: ET
3039 REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: SMC,SH2O,STC
3040 REAL, DIMENSION(1:NSOIL) :: ET1
3041 REAL :: EC1,EDIR1,ETT1,DF1,ETA1,ETP1,PRCP1,YY, &
3042 YYNUM,ZZ1
3043
3044 ! ----------------------------------------------------------------------
3045 ! EXECUTABLE CODE BEGINS HERE:
3046 ! CONVERT ETP Fnd PRCP FROM KG M-2 S-1 TO M S-1 AND INITIALIZE DEW.
3047 ! ----------------------------------------------------------------------
3048 PRCP1 = PRCP * 0.001
3049 ETP1 = ETP * 0.001
3050 DEW = 0.0
3051 ! ----------------------------------------------------------------------
3052 ! INITIALIZE EVAP TERMS.
3053 ! ----------------------------------------------------------------------
3054 EDIR = 0.
3055 EDIR1 = 0.
3056 EC1 = 0.
3057 EC = 0.
3058 DO K = 1,NSOIL
3059 ET(K) = 0.
3060 ET1(K) = 0.
3061 END DO
3062 ETT = 0.
3063 ETT1 = 0.
3064
3065 IF (ETP > 0.0) THEN
3066 CALL EVAPO (ETA1,SMC,NSOIL,CMC,ETP1,DT,ZSOIL, &
3067 SH2O, &
3068 SMCMAX,BEXP,PC,SMCWLT,DKSAT,DWSAT, &
3069 SMCREF,SHDFAC,CMCMAX, &
3070 SMCDRY,CFACTR, &
3071 EDIR1,EC1,ET1,ETT1,SFCTMP,Q2,NROOT,RTDIS,FXEXP)
3072 CALL SMFLX (SMC,NSOIL,CMC,DT,PRCP1,ZSOIL, &
3073 SH2O,SLOPE,KDT,FRZFACT, &
3074 SMCMAX,BEXP,SMCWLT,DKSAT,DWSAT, &
3075 SHDFAC,CMCMAX, &
3076 RUNOFF1,RUNOFF2,RUNOFF3, &
3077 EDIR1,EC1,ET1, &
3078 DRIP)
3079
3080 ! ----------------------------------------------------------------------
3081 ! CONVERT MODELED EVAPOTRANSPIRATION FROM M S-1 TO KG M-2 S-1.
3082 ! ----------------------------------------------------------------------
3083
3084 ETA = ETA1 * 1000.0
3085
3086 ! ----------------------------------------------------------------------
3087 ! IF ETP < 0, ASSUME DEW FORMS (TRANSFORM ETP1 INTO DEW AND REINITIALIZE
3088 ! ETP1 TO ZERO).
3089 ! ----------------------------------------------------------------------
3090 ELSE
3091 DEW = - ETP1
3092
3093 ! ----------------------------------------------------------------------
3094 ! CONVERT PRCP FROM 'KG M-2 S-1' TO 'M S-1' AND ADD DEW AMOUNT.
3095 ! ----------------------------------------------------------------------
3096
3097 PRCP1 = PRCP1+ DEW
3098 CALL SMFLX (SMC,NSOIL,CMC,DT,PRCP1,ZSOIL, &
3099 SH2O,SLOPE,KDT,FRZFACT, &
3100 SMCMAX,BEXP,SMCWLT,DKSAT,DWSAT, &
3101 SHDFAC,CMCMAX, &
3102 RUNOFF1,RUNOFF2,RUNOFF3, &
3103 EDIR1,EC1,ET1, &
3104 DRIP)
3105
3106 ! ----------------------------------------------------------------------
3107 ! CONVERT MODELED EVAPOTRANSPIRATION FROM 'M S-1' TO 'KG M-2 S-1'.
3108 ! ----------------------------------------------------------------------
3109 ! ETA = ETA1 * 1000.0
3110 END IF
3111
3112 ! ----------------------------------------------------------------------
3113 ! BASED ON ETP AND E VALUES, DETERMINE BETA
3114 ! ----------------------------------------------------------------------
3115
3116 IF ( ETP <= 0.0 ) THEN
3117 BETA = 0.0
3118 IF ( ETP < 0.0 ) THEN
3119 BETA = 1.0
3120 ETA = ETP
3121 END IF
3122 ELSE
3123 BETA = ETA / ETP
3124 END IF
3125
3126 ! ----------------------------------------------------------------------
3127 ! CONVERT MODELED EVAPOTRANSPIRATION COMPONENTS 'M S-1' TO 'KG M-2 S-1'.
3128 ! ----------------------------------------------------------------------
3129 EDIR = EDIR1*1000.
3130 EC = EC1*1000.
3131 DO K = 1,NSOIL
3132 ET(K) = ET1(K)*1000.
3133 END DO
3134 ETT = ETT1*1000.
3135
3136 ! ----------------------------------------------------------------------
3137 ! GET SOIL THERMAL DIFFUXIVITY/CONDUCTIVITY FOR TOP SOIL LYR,
3138 ! CALC. ADJUSTED TOP LYR SOIL TEMP AND ADJUSTED SOIL FLUX, THEN
3139 ! CALL SHFLX TO COMPUTE/UPDATE SOIL HEAT FLUX AND SOIL TEMPS.
3140 ! ----------------------------------------------------------------------
3141
3142 CALL TDFCND (DF1,SMC (1),QUARTZ,SMCMAX,SH2O (1))
3143
3144 !urban
3145 IF (VEGTYP==1) DF1=3.24
3146 !
3147
3148 ! ----------------------------------------------------------------------
3149 ! VEGETATION GREENNESS FRACTION REDUCTION IN SUBSURFACE HEAT FLUX
3150 ! VIA REDUCTION FACTOR, WHICH IS CONVENIENT TO APPLY HERE TO THERMAL
3151 ! DIFFUSIVITY THAT IS LATER USED IN HRT TO COMPUTE SUB SFC HEAT FLUX
3152 ! (SEE ADDITIONAL COMMENTS ON VEG EFFECT SUB-SFC HEAT FLX IN
3153 ! ROUTINE SFLX)
3154 ! ----------------------------------------------------------------------
3155 DF1 = DF1 * EXP (SBETA * SHDFAC)
3156 ! ----------------------------------------------------------------------
3157 ! COMPUTE INTERMEDIATE TERMS PASSED TO ROUTINE HRT (VIA ROUTINE
3158 ! SHFLX BELOW) FOR USE IN COMPUTING SUBSURFACE HEAT FLUX IN HRT
3159 ! ----------------------------------------------------------------------
3160 YYNUM = FDOWN - EMISSI*SIGMA * T24
3161 YY = SFCTMP + (YYNUM / RCH + TH2- SFCTMP - BETA * EPSCA) / RR
3162
3163 ZZ1 = DF1 / ( -0.5 * ZSOIL (1) * RCH * RR ) + 1.0
3164
3165 !urban
3166 CALL SHFLX (SSOIL,STC,SMC,SMCMAX,NSOIL,T1,DT,YY,ZZ1,ZSOIL, &
3167 TBOT,ZBOT,SMCWLT,PSISAT,SH2O,BEXP,F1,DF1,ICE, &
3168 QUARTZ,CSOIL,VEGTYP)
3169
3170 ! ----------------------------------------------------------------------
3171 ! SET FLX1 AND FLX3 (SNOPACK PHASE CHANGE HEAT FLUXES) TO ZERO SINCE
3172 ! THEY ARE NOT USED HERE IN SNOPAC. FLX2 (FREEZING RAIN HEAT FLUX) WAS
3173 ! SIMILARLY INITIALIZED IN THE PENMAN ROUTINE.
3174 ! ----------------------------------------------------------------------
3175 FLX1 = 0.0
3176 FLX3 = 0.0
3177
3178 ! ----------------------------------------------------------------------
3179 END SUBROUTINE NOPAC
3180 ! ----------------------------------------------------------------------
3181
3182 SUBROUTINE PENMAN (SFCTMP,SFCPRS,CH,T2V,TH2,PRCP,FDOWN,T24,SSOIL, &
3183 & Q2,Q2SAT,ETP,RCH,EPSCA,RR,SNOWNG,FRZGRA, &
3184 & DQSDT2,FLX2,EMISSI_IN,SNEQV)
3185
3186 ! ----------------------------------------------------------------------
3187 ! SUBROUTINE PENMAN
3188 ! ----------------------------------------------------------------------
3189 ! CALCULATE POTENTIAL EVAPORATION FOR THE CURRENT POINT. VARIOUS
3190 ! PARTIAL SUMS/PRODUCTS ARE ALSO CALCULATED AND PASSED BACK TO THE
3191 ! CALLING ROUTINE FOR LATER USE.
3192 ! ----------------------------------------------------------------------
3193 IMPLICIT NONE
3194 LOGICAL, INTENT(IN) :: SNOWNG, FRZGRA
3195 REAL, INTENT(IN) :: CH, DQSDT2,FDOWN,PRCP, &
3196 Q2, Q2SAT,SSOIL, SFCPRS, SFCTMP, &
3197 T2V, TH2,EMISSI_IN,SNEQV
3198 REAL, INTENT(OUT) :: EPSCA,ETP,FLX2,RCH,RR,T24
3199 REAL :: A, DELTA, FNET,RAD,RHO,EMISSI,ELCP1,LVS
3200
3201 REAL, PARAMETER :: ELCP = 2.4888E+3, LSUBC = 2.501000E+6,CP = 1004.6
3202 REAL, PARAMETER :: LSUBS = 2.83E+6
3203
3204 ! ----------------------------------------------------------------------
3205 ! EXECUTABLE CODE BEGINS HERE:
3206 ! ----------------------------------------------------------------------
3207 ! ----------------------------------------------------------------------
3208 ! PREPARE PARTIAL QUANTITIES FOR PENMAN EQUATION.
3209 ! ----------------------------------------------------------------------
3210 IF(SNEQV == 0.)THEN
3211 EMISSI=EMISSI_IN
3212 ELCP1=ELCP
3213 LVS=LSUBC
3214 ELSE
3215 EMISSI=0.95 ! FOR SNOW
3216 ELCP1=ELCP*LSUBS/LSUBC
3217 LVS=LSUBS
3218 ENDIF
3219
3220 FLX2 = 0.0
3221 ! DELTA = ELCP * DQSDT2
3222 DELTA = ELCP1 * DQSDT2
3223 T24 = SFCTMP * SFCTMP * SFCTMP * SFCTMP
3224 ! RR = T24 * 6.48E-8 / (SFCPRS * CH) + 1.0
3225 RR = EMISSI*T24 * 6.48E-8 / (SFCPRS * CH) + 1.0
3226 RHO = SFCPRS / (RD * T2V)
3227
3228 ! ----------------------------------------------------------------------
3229 ! ADJUST THE PARTIAL SUMS / PRODUCTS WITH THE LATENT HEAT
3230 ! EFFECTS CAUSED BY FALLING PRECIPITATION.
3231 ! ----------------------------------------------------------------------
3232 RCH = RHO * CP * CH
3233 IF (.NOT. SNOWNG) THEN
3234 IF (PRCP > 0.0) RR = RR + CPH2O * PRCP / RCH
3235 ELSE
3236 RR = RR + CPICE * PRCP / RCH
3237 END IF
3238
3239 ! ----------------------------------------------------------------------
3240 ! INCLUDE THE LATENT HEAT EFFECTS OF FRZNG RAIN CONVERTING TO ICE ON
3241 ! IMPACT IN THE CALCULATION OF FLX2 AND FNET.
3242 ! ----------------------------------------------------------------------
3243 ! FNET = FDOWN - SIGMA * T24- SSOIL
3244 FNET = FDOWN - EMISSI*SIGMA * T24- SSOIL
3245 IF (FRZGRA) THEN
3246 FLX2 = - LSUBF * PRCP
3247 FNET = FNET - FLX2
3248 ! ----------------------------------------------------------------------
3249 ! FINISH PENMAN EQUATION CALCULATIONS.
3250 ! ----------------------------------------------------------------------
3251 END IF
3252 RAD = FNET / RCH + TH2- SFCTMP
3253 ! A = ELCP * (Q2SAT - Q2)
3254 A = ELCP1 * (Q2SAT - Q2)
3255 EPSCA = (A * RR + RAD * DELTA) / (DELTA + RR)
3256 ! ETP = EPSCA * RCH / LSUBC
3257 ETP = EPSCA * RCH / LVS
3258
3259 ! ----------------------------------------------------------------------
3260 END SUBROUTINE PENMAN
3261 ! ----------------------------------------------------------------------
3262
3263 SUBROUTINE REDPRM (VEGTYP,SOILTYP,SLOPETYP,CFACTR,CMCMAX,RSMAX, &
3264 TOPT, &
3265 REFKDT,KDT,SBETA, SHDFAC,RSMIN,RGL,HS,ZBOT,FRZX, &
3266 PSISAT,SLOPE,SNUP,SALP,BEXP,DKSAT,DWSAT, &
3267 SMCMAX,SMCWLT,SMCREF,SMCDRY,F1,QUARTZ,FXEXP, &
3268 RTDIS,SLDPTH,ZSOIL, NROOT,NSOIL,Z0BRD,CZIL,LAI, &
3269 CSOIL,ALBBRD,PTU,LLANDUSE,LSOIL,LOCAL)
3270
3271 IMPLICIT NONE
3272 ! ----------------------------------------------------------------------
3273 ! Internally set (default valuess)
3274 ! all soil and vegetation parameters required for the execusion oF
3275 ! the Noah lsm are defined in VEGPARM.TBL, SOILPARM.TB, and GENPARM.TBL.
3276 ! ----------------------------------------------------------------------
3277 ! Vegetation parameters:
3278 ! ALBBRD: SFC background snow-free albedo
3279 ! CMXTBL: MAX CNPY Capacity
3280 ! Z0BRD: Background roughness length
3281 ! SHDFAC: Green vegetation fraction
3282 ! NROOT: Rooting depth
3283 ! RSMIN: Mimimum stomatal resistance
3284 ! RSMAX: Max. stomatal resistance
3285 ! RGL: Parameters used in radiation stress function
3286 ! HS: Parameter used in vapor pressure deficit functio
3287 ! TOPT: Optimum transpiration air temperature.
3288 ! CMCMAX: Maximum canopy water capacity
3289 ! CFACTR: Parameter used in the canopy inteception calculation
3290 ! SNUP: Threshold snow depth (in water equivalent m) that
3291 ! implies 100 percent snow cover
3292 ! LAI: Leaf area index
3293 !
3294 ! ----------------------------------------------------------------------
3295 ! Soil parameters:
3296 ! SMCMAX: MAX soil moisture content (porosity)
3297 ! SMCREF: Reference soil moisture (field capacity)
3298 ! SMCWLT: Wilting point soil moisture
3299 ! SMCWLT: Air dry soil moist content limits
3300 ! SSATPSI: SAT (saturation) soil potential
3301 ! DKSAT: SAT soil conductivity
3302 ! BEXP: B parameter
3303 ! SSATDW: SAT soil diffusivity
3304 ! F1: Soil thermal diffusivity/conductivity coef.
3305 ! QUARTZ: Soil quartz content
3306 ! Modified by F. Chen (12/22/97) to use the STATSGO soil map
3307 ! Modified By F. Chen (01/22/00) to include PLaya, Lava, and White San
3308 ! Modified By F. Chen (08/05/02) to include additional parameters for the Noah
3309 ! NOTE: SATDW = BB*SATDK*(SATPSI/MAXSMC)
3310 ! F11 = ALOG10(SATPSI) + BB*ALOG10(MAXSMC) + 2.0
3311 ! REFSMC1=MAXSMC*(5.79E-9/SATDK)**(1/(2*BB+3)) 5.79E-9 m/s= 0.5 mm
3312 ! REFSMC=REFSMC1+1./3.(MAXSMC-REFSMC1)
3313 ! WLTSMC1=MAXSMC*(200./SATPSI)**(-1./BB) (Wetzel and Chang, 198
3314 ! WLTSMC=WLTSMC1-0.5*WLTSMC1
3315 ! Note: the values for playa is set for it to have a thermal conductivit
3316 ! as sand and to have a hydrulic conductivity as clay
3317 !
3318 ! ----------------------------------------------------------------------
3319 ! Class parameter 'SLOPETYP' was included to estimate linear reservoir
3320 ! coefficient 'SLOPE' to the baseflow runoff out of the bottom layer.
3321 ! lowest class (slopetyp=0) means highest slope parameter = 1.
3322 ! definition of slopetyp from 'zobler' slope type:
3323 ! slope class percent slope
3324 ! 1 0-8
3325 ! 2 8-30
3326 ! 3 > 30
3327 ! 4 0-30
3328 ! 5 0-8 & > 30
3329 ! 6 8-30 & > 30
3330 ! 7 0-8, 8-30, > 30
3331 ! 9 GLACIAL ICE
3332 ! BLANK OCEAN/SEA
3333 ! SLOPE_DATA: linear reservoir coefficient
3334 ! SBETA_DATA: parameter used to caluculate vegetation effect on soil heat
3335 ! FXEXP_DAT: soil evaporation exponent used in DEVAP
3336 ! CSOIL_DATA: soil heat capacity [J M-3 K-1]
3337 ! SALP_DATA: shape parameter of distribution function of snow cover
3338 ! REFDK_DATA and REFKDT_DATA: parameters in the surface runoff parameteriz
3339 ! FRZK_DATA: frozen ground parameter
3340 ! ZBOT_DATA: depth[M] of lower boundary soil temperature
3341 ! CZIL_DATA: calculate roughness length of heat
3342 ! SMLOW_DATA and MHIGH_DATA: two soil moisture wilt, soil moisture referen
3343 ! parameters
3344 ! Set maximum number of soil-, veg-, and slopetyp in data statement.
3345 ! ----------------------------------------------------------------------
3346 INTEGER, PARAMETER :: MAX_SLOPETYP=30,MAX_SOILTYP=30,MAX_VEGTYP=30
3347 LOGICAL :: LOCAL
3348 CHARACTER (LEN=4), INTENT(IN):: LLANDUSE, LSOIL
3349
3350 ! Veg parameters
3351 INTEGER, INTENT(IN) :: VEGTYP
3352 INTEGER, INTENT(OUT) :: NROOT
3353 REAL, INTENT(OUT) :: HS,LAI,RSMIN,RGL,SHDFAC,SNUP,Z0BRD, &
3354 CMCMAX,RSMAX,TOPT,ALBBRD
3355 ! Soil parameters
3356 INTEGER, INTENT(IN) :: SOILTYP
3357 REAL, INTENT(OUT) :: BEXP,DKSAT,DWSAT,F1,QUARTZ,SMCDRY, &
3358 SMCMAX,SMCREF,SMCWLT,PSISAT
3359 ! General parameters
3360 INTEGER, INTENT(IN) :: SLOPETYP,NSOIL
3361 INTEGER :: I
3362
3363 REAL, INTENT(OUT) :: SLOPE,CZIL,SBETA,FXEXP, &
3364 CSOIL,SALP,FRZX,KDT,CFACTR, &
3365 ZBOT,REFKDT,PTU
3366 REAL,DIMENSION(1:NSOIL),INTENT(IN) :: SLDPTH,ZSOIL
3367 REAL,DIMENSION(1:NSOIL),INTENT(OUT):: RTDIS
3368 REAL :: FRZFACT,FRZK,REFDK
3369
3370 ! SAVE
3371 ! ----------------------------------------------------------------------
3372 !
3373 IF (SOILTYP .gt. SLCATS) THEN
3374 CALL wrf_error_fatal ( 'Warning: too many input soil types' )
3375 END IF
3376 IF (VEGTYP .gt. LUCATS) THEN
3377 CALL wrf_error_fatal ( 'Warning: too many input landuse types' )
3378 END IF
3379 IF (SLOPETYP .gt. SLPCATS) THEN
3380 CALL wrf_error_fatal ( 'Warning: too many input slope types' )
3381 END IF
3382
3383 ! ----------------------------------------------------------------------
3384 ! SET-UP SOIL PARAMETERS
3385 ! ----------------------------------------------------------------------
3386 CSOIL = CSOIL_DATA
3387 BEXP = BB (SOILTYP)
3388 DKSAT = SATDK (SOILTYP)
3389 DWSAT = SATDW (SOILTYP)
3390 F1 = F11 (SOILTYP)
3391 PSISAT = SATPSI (SOILTYP)
3392 QUARTZ = QTZ (SOILTYP)
3393 SMCDRY = DRYSMC (SOILTYP)
3394 SMCMAX = MAXSMC (SOILTYP)
3395 SMCREF = REFSMC (SOILTYP)
3396 SMCWLT = WLTSMC (SOILTYP)
3397 ! ----------------------------------------------------------------------
3398 ! Set-up universal parameters (not dependent on SOILTYP, VEGTYP or
3399 ! SLOPETYP)
3400 ! ----------------------------------------------------------------------
3401 ZBOT = ZBOT_DATA
3402 SALP = SALP_DATA
3403 SBETA = SBETA_DATA
3404 REFDK = REFDK_DATA
3405 FRZK = FRZK_DATA
3406 FXEXP = FXEXP_DATA
3407 REFKDT = REFKDT_DATA
3408 PTU = 0. ! (not used yet) to satisify intent(out)
3409 KDT = REFKDT * DKSAT / REFDK
3410 CZIL = CZIL_DATA
3411 SLOPE = SLOPE_DATA (SLOPETYP)
3412
3413 ! ----------------------------------------------------------------------
3414 ! TO ADJUST FRZK PARAMETER TO ACTUAL SOIL TYPE: FRZK * FRZFACT
3415 ! ----------------------------------------------------------------------
3416 FRZFACT = (SMCMAX / SMCREF) * (0.412 / 0.468)
3417 FRZX = FRZK * FRZFACT
3418
3419 ! ----------------------------------------------------------------------
3420 ! SET-UP VEGETATION PARAMETERS
3421 ! ----------------------------------------------------------------------
3422 TOPT = TOPT_DATA
3423 CMCMAX = CMCMAX_DATA
3424 CFACTR = CFACTR_DATA
3425 RSMAX = RSMAX_DATA
3426 NROOT = NROTBL (VEGTYP)
3427 SNUP = SNUPTBL (VEGTYP)
3428 RSMIN = RSTBL (VEGTYP)
3429 RGL = RGLTBL (VEGTYP)
3430 HS = HSTBL (VEGTYP)
3431 LAI = LAITBL (VEGTYP)
3432 IF(LOCAL) THEN
3433 ALBBRD = ALBTBL(VEGTYP)
3434 Z0BRD = Z0TBL(VEGTYP)
3435 SHDFAC = SHDTBL(VEGTYP)
3436 ENDIF
3437
3438 IF (VEGTYP .eq. BARE) SHDFAC = 0.0
3439 IF (NROOT .gt. NSOIL) THEN
3440 WRITE (err_message,*) 'Error: too many root layers ', &
3441 NSOIL,NROOT
3442 CALL wrf_error_fatal ( err_message )
3443 ! ----------------------------------------------------------------------
3444 ! CALCULATE ROOT DISTRIBUTION. PRESENT VERSION ASSUMES UNIFORM
3445 ! DISTRIBUTION BASED ON SOIL LAYER DEPTHS.
3446 ! ----------------------------------------------------------------------
3447 END IF
3448 DO I = 1,NROOT
3449 RTDIS (I) = - SLDPTH (I)/ ZSOIL (NROOT)
3450 ! ----------------------------------------------------------------------
3451 ! SET-UP SLOPE PARAMETER
3452 ! ----------------------------------------------------------------------
3453 END DO
3454
3455 ! print*,'end of PRMRED'
3456 ! print*,'VEGTYP',VEGTYP,'SOILTYP',SOILTYP,'SLOPETYP',SLOPETYP, &
3457 ! & 'CFACTR',CFACTR,'CMCMAX',CMCMAX,'RSMAX',RSMAX,'TOPT',TOPT, &
3458 ! & 'REFKDT',REFKDT,'KDT',KDT,'SBETA',SBETA, 'SHDFAC',SHDFAC, &
3459 ! & 'RSMIN',RSMIN,'RGL',RGL,'HS',HS,'ZBOT',ZBOT,'FRZX',FRZX, &
3460 ! & 'PSISAT',PSISAT,'SLOPE',SLOPE,'SNUP',SNUP,'SALP',SALP,'BEXP', &
3461 ! & BEXP, &
3462 ! & 'DKSAT',DKSAT,'DWSAT',DWSAT, &
3463 ! & 'SMCMAX',SMCMAX,'SMCWLT',SMCWLT,'SMCREF',SMCREF,'SMCDRY',SMCDRY, &
3464 ! & 'F1',F1,'QUARTZ',QUARTZ,'FXEXP',FXEXP, &
3465 ! & 'RTDIS',RTDIS,'SLDPTH',SLDPTH,'ZSOIL',ZSOIL, 'NROOT',NROOT, &
3466 ! & 'NSOIL',NSOIL,'Z0',Z0,'CZIL',CZIL,'LAI',LAI, &
3467 ! & 'CSOIL',CSOIL,'PTU',PTU, &
3468 ! & 'LOCAL', LOCAL
3469
3470 END SUBROUTINE REDPRM
3471
3472 SUBROUTINE ROSR12 (P,A,B,C,D,DELTA,NSOIL)
3473
3474 ! ----------------------------------------------------------------------
3475 ! SUBROUTINE ROSR12
3476 ! ----------------------------------------------------------------------
3477 ! INVERT (SOLVE) THE TRI-DIAGONAL MATRIX PROBLEM SHOWN BELOW:
3478 ! ### ### ### ### ### ###
3479 ! #B(1), C(1), 0 , 0 , 0 , . . . , 0 # # # # #
3480 ! #A(2), B(2), C(2), 0 , 0 , . . . , 0 # # # # #
3481 ! # 0 , A(3), B(3), C(3), 0 , . . . , 0 # # # # D(3) #
3482 ! # 0 , 0 , A(4), B(4), C(4), . . . , 0 # # P(4) # # D(4) #
3483 ! # 0 , 0 , 0 , A(5), B(5), . . . , 0 # # P(5) # # D(5) #
3484 ! # . . # # . # = # . #
3485 ! # . . # # . # # . #
3486 ! # . . # # . # # . #
3487 ! # 0 , . . . , 0 , A(M-2), B(M-2), C(M-2), 0 # #P(M-2)# #D(M-2)#
3488 ! # 0 , . . . , 0 , 0 , A(M-1), B(M-1), C(M-1)# #P(M-1)# #D(M-1)#
3489 ! # 0 , . . . , 0 , 0 , 0 , A(M) , B(M) # # P(M) # # D(M) #
3490 ! ### ### ### ### ### ###
3491 ! ----------------------------------------------------------------------
3492 IMPLICIT NONE
3493
3494 INTEGER, INTENT(IN) :: NSOIL
3495 INTEGER :: K, KK
3496
3497 REAL, DIMENSION(1:NSOIL), INTENT(IN):: A, B, D
3498 REAL, DIMENSION(1:NSOIL),INTENT(INOUT):: C,P,DELTA
3499
3500 ! ----------------------------------------------------------------------
3501 ! INITIALIZE EQN COEF C FOR THE LOWEST SOIL LAYER
3502 ! ----------------------------------------------------------------------
3503 C (NSOIL) = 0.0
3504 P (1) = - C (1) / B (1)
3505 ! ----------------------------------------------------------------------
3506 ! SOLVE THE COEFS FOR THE 1ST SOIL LAYER
3507 ! ----------------------------------------------------------------------
3508
3509 ! ----------------------------------------------------------------------
3510 ! SOLVE THE COEFS FOR SOIL LAYERS 2 THRU NSOIL
3511 ! ----------------------------------------------------------------------
3512 DELTA (1) = D (1) / B (1)
3513 DO K = 2,NSOIL
3514 P (K) = - C (K) * ( 1.0 / (B (K) + A (K) * P (K -1)) )
3515 DELTA (K) = (D (K) - A (K)* DELTA (K -1))* (1.0/ (B (K) + A (K)&
3516 * P (K -1)))
3517 END DO
3518 ! ----------------------------------------------------------------------
3519 ! SET P TO DELTA FOR LOWEST SOIL LAYER
3520 ! ----------------------------------------------------------------------
3521 P (NSOIL) = DELTA (NSOIL)
3522
3523 ! ----------------------------------------------------------------------
3524 ! ADJUST P FOR SOIL LAYERS 2 THRU NSOIL
3525 ! ----------------------------------------------------------------------
3526 DO K = 2,NSOIL
3527 KK = NSOIL - K + 1
3528 P (KK) = P (KK) * P (KK +1) + DELTA (KK)
3529 END DO
3530 ! ----------------------------------------------------------------------
3531 END SUBROUTINE ROSR12
3532 ! ----------------------------------------------------------------------
3533
3534
3535 SUBROUTINE SHFLX (SSOIL,STC,SMC,SMCMAX,NSOIL,T1,DT,YY,ZZ1,ZSOIL, &
3536 TBOT,ZBOT,SMCWLT,PSISAT,SH2O,BEXP,F1,DF1,ICE, &
3537 QUARTZ,CSOIL,VEGTYP)
3538
3539 ! ----------------------------------------------------------------------
3540 ! SUBROUTINE SHFLX
3541 ! ----------------------------------------------------------------------
3542 ! UPDATE THE TEMPERATURE STATE OF THE SOIL COLUMN BASED ON THE THERMAL
3543 ! DIFFUSION EQUATION AND UPDATE THE FROZEN SOIL MOISTURE CONTENT BASED
3544 ! ON THE TEMPERATURE.
3545 ! ----------------------------------------------------------------------
3546 IMPLICIT NONE
3547
3548 INTEGER, INTENT(IN) :: ICE, NSOIL, VEGTYP
3549 INTEGER :: I
3550
3551 REAL, INTENT(IN) :: BEXP,CSOIL,DF1,DT,F1,PSISAT,QUARTZ, &
3552 SMCMAX, SMCWLT, TBOT,YY, ZBOT,ZZ1
3553 REAL, INTENT(INOUT) :: T1
3554 REAL, INTENT(OUT) :: SSOIL
3555 REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SMC,ZSOIL
3556 REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: SH2O
3557 REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: STC
3558 REAL, DIMENSION(1:NSOIL) :: AI, BI, CI, STCF,RHSTS
3559 REAL, PARAMETER :: T0 = 273.15
3560
3561 ! ----------------------------------------------------------------------
3562 ! HRT ROUTINE CALCS THE RIGHT HAND SIDE OF THE SOIL TEMP DIF EQN
3563 ! ----------------------------------------------------------------------
3564
3565 ! ----------------------------------------------------------------------
3566 ! SEA-ICE CASE
3567 ! ----------------------------------------------------------------------
3568 IF (ICE == 1) THEN
3569
3570 CALL HRTICE (RHSTS,STC,NSOIL,ZSOIL,YY,ZZ1,DF1,AI,BI,CI)
3571
3572 CALL HSTEP (STCF,STC,RHSTS,DT,NSOIL,AI,BI,CI)
3573
3574 ! ----------------------------------------------------------------------
3575 ! LAND-MASS CASE
3576 ! ----------------------------------------------------------------------
3577 ELSE
3578 CALL HRT (RHSTS,STC,SMC,SMCMAX,NSOIL,ZSOIL,YY,ZZ1,TBOT, &
3579 ZBOT,PSISAT,SH2O,DT, &
3580 BEXP,F1,DF1,QUARTZ,CSOIL,AI,BI,CI,VEGTYP)
3581
3582 CALL HSTEP (STCF,STC,RHSTS,DT,NSOIL,AI,BI,CI)
3583 END IF
3584 DO I = 1,NSOIL
3585 STC (I) = STCF (I)
3586 END DO
3587
3588 ! ----------------------------------------------------------------------
3589 ! IN THE NO SNOWPACK CASE (VIA ROUTINE NOPAC BRANCH,) UPDATE THE GRND
3590 ! (SKIN) TEMPERATURE HERE IN RESPONSE TO THE UPDATED SOIL TEMPERATURE
3591 ! PROFILE ABOVE. (NOTE: INSPECTION OF ROUTINE SNOPAC SHOWS THAT T1
3592 ! BELOW IS A DUMMY VARIABLE ONLY, AS SKIN TEMPERATURE IS UPDATED
3593 ! DIFFERENTLY IN ROUTINE SNOPAC)
3594 ! ----------------------------------------------------------------------
3595 ! ----------------------------------------------------------------------
3596 ! CALCULATE SURFACE SOIL HEAT FLUX
3597 ! ----------------------------------------------------------------------
3598 T1 = (YY + (ZZ1- 1.0) * STC (1)) / ZZ1
3599 SSOIL = DF1 * (STC (1) - T1) / (0.5 * ZSOIL (1))
3600
3601 ! ----------------------------------------------------------------------
3602 END SUBROUTINE SHFLX
3603 ! ----------------------------------------------------------------------
3604
3605 SUBROUTINE SMFLX (SMC,NSOIL,CMC,DT,PRCP1,ZSOIL, &
3606 & SH2O,SLOPE,KDT,FRZFACT, &
3607 & SMCMAX,BEXP,SMCWLT,DKSAT,DWSAT, &
3608 & SHDFAC,CMCMAX, &
3609 & RUNOFF1,RUNOFF2,RUNOFF3, &
3610 & EDIR,EC,ET, &
3611 & DRIP)
3612
3613 ! ----------------------------------------------------------------------
3614 ! SUBROUTINE SMFLX
3615 ! ----------------------------------------------------------------------
3616 ! CALCULATE SOIL MOISTURE FLUX. THE SOIL MOISTURE CONTENT (SMC - A PER
3617 ! UNIT VOLUME MEASUREMENT) IS A DEPENDENT VARIABLE THAT IS UPDATED WITH
3618 ! PROGNOSTIC EQNS. THE CANOPY MOISTURE CONTENT (CMC) IS ALSO UPDATED.
3619 ! FROZEN GROUND VERSION: NEW STATES ADDED: SH2O, AND FROZEN GROUND
3620 ! CORRECTION FACTOR, FRZFACT AND PARAMETER SLOPE.
3621 ! ----------------------------------------------------------------------
3622 IMPLICIT NONE
3623
3624 INTEGER, INTENT(IN) :: NSOIL
3625 INTEGER :: I,K
3626
3627 REAL, INTENT(IN) :: BEXP, CMCMAX, DKSAT,DWSAT, DT, EC, EDIR, &
3628 KDT, PRCP1, SHDFAC, SLOPE, SMCMAX, SMCWLT
3629 REAL, INTENT(OUT) :: DRIP, RUNOFF1, RUNOFF2, RUNOFF3
3630 REAL, INTENT(INOUT) :: CMC
3631 REAL, DIMENSION(1:NSOIL), INTENT(IN) :: ET,ZSOIL
3632 REAL, DIMENSION(1:NSOIL), INTENT(INOUT):: SMC, SH2O
3633 REAL, DIMENSION(1:NSOIL) :: AI, BI, CI, STCF,RHSTS, RHSTT, &
3634 SICE, SH2OA, SH2OFG
3635 REAL :: DUMMY, EXCESS,FRZFACT,PCPDRP,RHSCT,TRHSCT
3636
3637 ! ----------------------------------------------------------------------
3638 ! EXECUTABLE CODE BEGINS HERE.
3639 ! ----------------------------------------------------------------------
3640 ! ----------------------------------------------------------------------
3641 ! COMPUTE THE RIGHT HAND SIDE OF THE CANOPY EQN TERM ( RHSCT )
3642 ! ----------------------------------------------------------------------
3643 DUMMY = 0.
3644
3645 ! ----------------------------------------------------------------------
3646 ! CONVERT RHSCT (A RATE) TO TRHSCT (AN AMOUNT) AND ADD IT TO EXISTING
3647 ! CMC. IF RESULTING AMT EXCEEDS MAX CAPACITY, IT BECOMES DRIP AND WILL
3648 ! FALL TO THE GRND.
3649 ! ----------------------------------------------------------------------
3650 RHSCT = SHDFAC * PRCP1- EC
3651 DRIP = 0.
3652 TRHSCT = DT * RHSCT
3653 EXCESS = CMC + TRHSCT
3654
3655 ! ----------------------------------------------------------------------
3656 ! PCPDRP IS THE COMBINED PRCP1 AND DRIP (FROM CMC) THAT GOES INTO THE
3657 ! SOIL
3658 ! ----------------------------------------------------------------------
3659 IF (EXCESS > CMCMAX) DRIP = EXCESS - CMCMAX
3660 PCPDRP = (1. - SHDFAC) * PRCP1+ DRIP / DT
3661
3662 ! ----------------------------------------------------------------------
3663 ! STORE ICE CONTENT AT EACH SOIL LAYER BEFORE CALLING SRT and SSTEP
3664 !
3665 DO I = 1,NSOIL
3666 SICE (I) = SMC (I) - SH2O (I)
3667 END DO
3668 ! ----------------------------------------------------------------------
3669 ! CALL SUBROUTINES SRT AND SSTEP TO SOLVE THE SOIL MOISTURE
3670 ! TENDENCY EQUATIONS.
3671 ! IF THE INFILTRATING PRECIP RATE IS NONTRIVIAL,
3672 ! (WE CONSIDER NONTRIVIAL TO BE A PRECIP TOTAL OVER THE TIME STEP
3673 ! EXCEEDING ONE ONE-THOUSANDTH OF THE WATER HOLDING CAPACITY OF
3674 ! THE FIRST SOIL LAYER)
3675 ! THEN CALL THE SRT/SSTEP SUBROUTINE PAIR TWICE IN THE MANNER OF
3676 ! TIME SCHEME "F" (IMPLICIT STATE, AVERAGED COEFFICIENT)
3677 ! OF SECTION 2 OF KALNAY AND KANAMITSU (1988, MWR, VOL 116,
3678 ! PAGES 1945-1958)TO MINIMIZE 2-DELTA-T OSCILLATIONS IN THE
3679 ! SOIL MOISTURE VALUE OF THE TOP SOIL LAYER THAT CAN ARISE BECAUSE
3680 ! OF THE EXTREME NONLINEAR DEPENDENCE OF THE SOIL HYDRAULIC
3681 ! DIFFUSIVITY COEFFICIENT AND THE HYDRAULIC CONDUCTIVITY ON THE
3682 ! SOIL MOISTURE STATE
3683 ! OTHERWISE CALL THE SRT/SSTEP SUBROUTINE PAIR ONCE IN THE MANNER OF
3684 ! TIME SCHEME "D" (IMPLICIT STATE, EXPLICIT COEFFICIENT)
3685 ! OF SECTION 2 OF KALNAY AND KANAMITSU
3686 ! PCPDRP IS UNITS OF KG/M**2/S OR MM/S, ZSOIL IS NEGATIVE DEPTH IN M
3687 ! ----------------------------------------------------------------------
3688 ! IF ( PCPDRP .GT. 0.0 ) THEN
3689
3690 ! ----------------------------------------------------------------------
3691 ! FROZEN GROUND VERSION:
3692 ! SMC STATES REPLACED BY SH2O STATES IN SRT SUBR. SH2O & SICE STATES
3693 ! INC&UDED IN SSTEP SUBR. FROZEN GROUND CORRECTION FACTOR, FRZFACT
3694 ! ADDED. ALL WATER BALANCE CALCULATIONS USING UNFROZEN WATER
3695 ! ----------------------------------------------------------------------
3696 IF ( (PCPDRP * DT) > (0.001*1000.0* (- ZSOIL (1))* SMCMAX) ) THEN
3697 CALL SRT (RHSTT,EDIR,ET,SH2O,SH2O,NSOIL,PCPDRP,ZSOIL, &
3698 DWSAT,DKSAT,SMCMAX,BEXP,RUNOFF1, &
3699 RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZFACT,SICE,AI,BI,CI)
3700 CALL SSTEP (SH2OFG,SH2O,DUMMY,RHSTT,RHSCT,DT,NSOIL,SMCMAX, &
3701 CMCMAX,RUNOFF3,ZSOIL,SMC,SICE,AI,BI,CI)
3702 DO K = 1,NSOIL
3703 SH2OA (K) = (SH2O (K) + SH2OFG (K)) * 0.5
3704 END DO
3705 CALL SRT (RHSTT,EDIR,ET,SH2O,SH2OA,NSOIL,PCPDRP,ZSOIL, &
3706 DWSAT,DKSAT,SMCMAX,BEXP,RUNOFF1, &
3707 RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZFACT,SICE,AI,BI,CI)
3708 CALL SSTEP (SH2O,SH2O,CMC,RHSTT,RHSCT,DT,NSOIL,SMCMAX, &
3709 CMCMAX,RUNOFF3,ZSOIL,SMC,SICE,AI,BI,CI)
3710
3711 ELSE
3712 CALL SRT (RHSTT,EDIR,ET,SH2O,SH2O,NSOIL,PCPDRP,ZSOIL, &
3713 DWSAT,DKSAT,SMCMAX,BEXP,RUNOFF1, &
3714 RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZFACT,SICE,AI,BI,CI)
3715 CALL SSTEP (SH2O,SH2O,CMC,RHSTT,RHSCT,DT,NSOIL,SMCMAX, &
3716 CMCMAX,RUNOFF3,ZSOIL,SMC,SICE,AI,BI,CI)
3717 ! RUNOF = RUNOFF
3718
3719 END IF
3720
3721 ! ----------------------------------------------------------------------
3722 END SUBROUTINE SMFLX
3723 ! ----------------------------------------------------------------------
3724
3725
3726 SUBROUTINE SNFRAC (SNEQV,SNUP,SALP,SNOWH,SNCOVR)
3727
3728 ! ----------------------------------------------------------------------
3729 ! SUBROUTINE SNFRAC
3730 ! ----------------------------------------------------------------------
3731 ! CALCULATE SNOW FRACTION (0 -> 1)
3732 ! SNEQV SNOW WATER EQUIVALENT (M)
3733 ! SNUP THRESHOLD SNEQV DEPTH ABOVE WHICH SNCOVR=1
3734 ! SALP TUNING PARAMETER
3735 ! SNCOVR FRACTIONAL SNOW COVER
3736 ! ----------------------------------------------------------------------
3737 IMPLICIT NONE
3738
3739 REAL, INTENT(IN) :: SNEQV,SNUP,SALP,SNOWH
3740 REAL, INTENT(OUT) :: SNCOVR
3741 REAL :: RSNOW, Z0N
3742
3743 ! ----------------------------------------------------------------------
3744 ! SNUP IS VEG-CLASS DEPENDENT SNOWDEPTH THRESHHOLD (SET IN ROUTINE
3745 ! REDPRM) ABOVE WHICH SNOCVR=1.
3746 ! ----------------------------------------------------------------------
3747 IF (SNEQV < SNUP) THEN
3748 RSNOW = SNEQV / SNUP
3749 SNCOVR = 1. - ( EXP ( - SALP * RSNOW) - RSNOW * EXP ( - SALP))
3750 ELSE
3751 SNCOVR = 1.0
3752 END IF
3753
3754 ! FORMULATION OF DICKINSON ET AL. 1986
3755 ! Z0N = 0.035
3756
3757 ! SNCOVR=SNOWH/(SNOWH + 5*Z0N)
3758
3759 ! FORMULATION OF MARSHALL ET AL. 1994
3760 ! SNCOVR=SNEQV/(SNEQV + 2*Z0N)
3761
3762 ! ----------------------------------------------------------------------
3763 END SUBROUTINE SNFRAC
3764 ! ----------------------------------------------------------------------
3765
3766 SUBROUTINE SNKSRC (TSNSR,TAVG,SMC,SH2O,ZSOIL,NSOIL, &
3767 & SMCMAX,PSISAT,BEXP,DT,K,QTOT)
3768
3769 ! ----------------------------------------------------------------------
3770 ! SUBROUTINE SNKSRC
3771 ! ----------------------------------------------------------------------
3772 ! CALCULATE SINK/SOURCE TERM OF THE TERMAL DIFFUSION EQUATION. (SH2O) IS
3773 ! AVAILABLE LIQUED WATER.
3774 ! ----------------------------------------------------------------------
3775 IMPLICIT NONE
3776
3777 INTEGER, INTENT(IN) :: K,NSOIL
3778
3779 REAL, INTENT(IN) :: BEXP, DT, PSISAT, QTOT, SMC, SMCMAX, &
3780 TAVG
3781 REAL, INTENT(INOUT) :: SH2O
3782
3783 REAL, DIMENSION(1:NSOIL), INTENT(IN):: ZSOIL
3784
3785 REAL :: DF, DZ, DZH, FREE, TSNSR, &
3786 TDN, TM, TUP, TZ, X0, XDN, XH2O, XUP
3787
3788 REAL, PARAMETER :: DH2O = 1.0000E3, HLICE = 3.3350E5, &
3789 T0 = 2.7315E2
3790
3791 IF (K == 1) THEN
3792 DZ = - ZSOIL (1)
3793 ELSE
3794 DZ = ZSOIL (K -1) - ZSOIL (K)
3795 END IF
3796 ! ----------------------------------------------------------------------
3797 ! VIA FUNCTION FRH2O, COMPUTE POTENTIAL OR 'EQUILIBRIUM' UNFROZEN
3798 ! SUPERCOOLED FREE WATER FOR GIVEN SOIL TYPE AND SOIL LAYER TEMPERATURE.
3799 ! FUNCTION FRH20 INVOKES EQN (17) FROM V. KOREN ET AL (1999, JGR, VOL.
3800 ! 104, PG 19573). (ASIDE: LATTER EQN IN JOURNAL IN CENTIGRADE UNITS.
3801 ! ROUTINE FRH2O USE FORM OF EQN IN KELVIN UNITS.)
3802 ! ----------------------------------------------------------------------
3803 ! FREE = FRH2O(TAVG,SMC,SH2O,SMCMAX,BEXP,PSISAT)
3804
3805 ! ----------------------------------------------------------------------
3806 ! IN NEXT BLOCK OF CODE, INVOKE EQN 18 OF V. KOREN ET AL (1999, JGR,
3807 ! VOL. 104, PG 19573.) THAT IS, FIRST ESTIMATE THE NEW AMOUNTOF LIQUID
3808 ! WATER, 'XH2O', IMPLIED BY THE SUM OF (1) THE LIQUID WATER AT THE BEGIN
3809 ! OF CURRENT TIME STEP, AND (2) THE FREEZE OF THAW CHANGE IN LIQUID
3810 ! WATER IMPLIED BY THE HEAT FLUX 'QTOT' PASSED IN FROM ROUTINE HRT.
3811 ! SECOND, DETERMINE IF XH2O NEEDS TO BE BOUNDED BY 'FREE' (EQUIL AMT) OR
3812 ! IF 'FREE' NEEDS TO BE BOUNDED BY XH2O.
3813 ! ----------------------------------------------------------------------
3814 CALL FRH2O (FREE,TAVG,SMC,SH2O,SMCMAX,BEXP,PSISAT)
3815
3816 ! ----------------------------------------------------------------------
3817 ! FIRST, IF FREEZING AND REMAINING LIQUID LESS THAN LOWER BOUND, THEN
3818 ! REDUCE EXTENT OF FREEZING, THEREBY LETTING SOME OR ALL OF HEAT FLUX
3819 ! QTOT COOL THE SOIL TEMP LATER IN ROUTINE HRT.
3820 ! ----------------------------------------------------------------------
3821 XH2O = SH2O + QTOT * DT / (DH2O * HLICE * DZ)
3822 IF ( XH2O < SH2O .AND. XH2O < FREE) THEN
3823 IF ( FREE > SH2O ) THEN
3824 XH2O = SH2O
3825 ELSE
3826 XH2O = FREE
3827 END IF
3828 END IF
3829 ! ----------------------------------------------------------------------
3830 ! SECOND, IF THAWING AND THE INCREASE IN LIQUID WATER GREATER THAN UPPER
3831 ! BOUND, THEN REDUCE EXTENT OF THAW, THEREBY LETTING SOME OR ALL OF HEAT
3832 ! FLUX QTOT WARM THE SOIL TEMP LATER IN ROUTINE HRT.
3833 ! ----------------------------------------------------------------------
3834 IF ( XH2O > SH2O .AND. XH2O > FREE ) THEN
3835 IF ( FREE < SH2O ) THEN
3836 XH2O = SH2O
3837 ELSE
3838 XH2O = FREE
3839 END IF
3840 END IF
3841
3842 ! ----------------------------------------------------------------------
3843 ! CALCULATE PHASE-CHANGE HEAT SOURCE/SINK TERM FOR USE IN ROUTINE HRT
3844 ! AND UPDATE LIQUID WATER TO REFLCET FINAL FREEZE/THAW INCREMENT.
3845 ! ----------------------------------------------------------------------
3846 ! SNKSRC = -DH2O*HLICE*DZ*(XH2O-SH2O)/DT
3847 IF (XH2O < 0.) XH2O = 0.
3848 IF (XH2O > SMC) XH2O = SMC
3849 TSNSR = - DH2O * HLICE * DZ * (XH2O - SH2O)/ DT
3850 SH2O = XH2O
3851
3852 ! ----------------------------------------------------------------------
3853 END SUBROUTINE SNKSRC
3854 ! ----------------------------------------------------------------------
3855
3856 SUBROUTINE SNOPAC (ETP,ETA,PRCP,PRCPF,SNOWNG,SMC,SMCMAX,SMCWLT, &
3857 SMCREF,SMCDRY,CMC,CMCMAX,NSOIL,DT, &
3858 SBETA,DF1, &
3859 Q2,T1,SFCTMP,T24,TH2,FDOWN,F1,SSOIL,STC,EPSCA,&
3860 SFCPRS,BEXP,PC,RCH,RR,CFACTR,SNCOVR,ESD,SNDENS,&
3861 SNOWH,SH2O,SLOPE,KDT,FRZFACT,PSISAT, &
3862 ZSOIL,DWSAT,DKSAT,TBOT,ZBOT,SHDFAC,RUNOFF1, &
3863 RUNOFF2,RUNOFF3,EDIR,EC,ET,ETT,NROOT,SNOMLT, &
3864 ICE,RTDIS,QUARTZ,FXEXP,CSOIL, &
3865 BETA,DRIP,DEW,FLX1,FLX2,FLX3,ESNOW,ETNS,EMISSI,&
3866 VEGTYP)
3867
3868 ! ----------------------------------------------------------------------
3869 ! SUBROUTINE SNOPAC
3870 ! ----------------------------------------------------------------------
3871 ! CALCULATE SOIL MOISTURE AND HEAT FLUX VALUES & UPDATE SOIL MOISTURE
3872 ! CONTENT AND SOIL HEAT CONTENT VALUES FOR THE CASE WHEN A SNOW PACK IS
3873 ! PRESENT.
3874 ! ----------------------------------------------------------------------
3875 IMPLICIT NONE
3876
3877 INTEGER, INTENT(IN) :: ICE, NROOT, NSOIL,VEGTYP
3878 INTEGER :: K
3879 LOGICAL, INTENT(IN) :: SNOWNG
3880 REAL, INTENT(IN) :: BEXP,CFACTR, CMCMAX,CSOIL,DF1,DKSAT, &
3881 DT,DWSAT, EPSCA,ETP,FDOWN,F1,FXEXP, &
3882 FRZFACT,KDT,PC, PRCP,PSISAT,Q2,QUARTZ, &
3883 RCH,RR,SBETA,SFCPRS, SFCTMP, SHDFAC, &
3884 SLOPE,SMCDRY,SMCMAX,SMCREF,SMCWLT, T24, &
3885 TBOT,TH2,ZBOT,EMISSI
3886 REAL, INTENT(INOUT) :: CMC, BETA, ESD,FLX2,PRCPF,SNOWH,SNCOVR, &
3887 SNDENS, T1
3888 REAL, INTENT(OUT) :: DEW,DRIP,EC,EDIR, ETNS, ESNOW,ETT, &
3889 FLX1,FLX3, RUNOFF1,RUNOFF2,RUNOFF3, &
3890 SSOIL,SNOMLT
3891 REAL, DIMENSION(1:NSOIL),INTENT(IN) :: RTDIS,ZSOIL
3892 REAL, DIMENSION(1:NSOIL),INTENT(OUT) :: ET
3893 REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: SMC,SH2O,STC
3894 REAL, DIMENSION(1:NSOIL) :: ET1
3895 REAL :: DENOM,DSOIL,DTOT,EC1,EDIR1,ESDFLX,ETA, &
3896 ETT1, ESNOW1, ESNOW2, ETA1,ETP1,ETP2, &
3897 ETP3, ETNS1, ETANRG, ETAX, EX, FLX3X, &
3898 FRCSNO,FRCSOI, PRCP1, QSAT,RSNOW, SEH, &
3899 SNCOND,SSOIL1, T11,T12, T12A, T12AX, &
3900 T12B, T14, YY, ZZ1, EMISSI_S
3901
3902 REAL, PARAMETER :: ESDMIN = 1.E-6, LSUBC = 2.501000E+6, &
3903 LSUBS = 2.83E+6, TFREEZ = 273.15, &
3904 SNOEXP = 2.0
3905
3906 ! ----------------------------------------------------------------------
3907 ! EXECUTABLE CODE BEGINS HERE:
3908 ! INITIALIZE EVAP TERMS.
3909 ! ----------------------------------------------------------------------
3910 ! conversions:
3911 ! ESNOW [KG M-2 S-1]
3912 ! ESDFLX [KG M-2 S-1] .le. ESNOW
3913 ! ESNOW1 [M S-1]
3914 ! ESNOW2 [M]
3915 ! ETP [KG M-2 S-1]
3916 ! ETP1 [M S-1]
3917 ! ETP2 [M]
3918 ! ----------------------------------------------------------------------
3919 DEW = 0.
3920 EDIR = 0.
3921 EDIR1 = 0.
3922 EC1 = 0.
3923 EC = 0.
3924 EMISSI_S=0.95 ! For snow
3925
3926 DO K = 1,NSOIL
3927 ET (K) = 0.
3928 ET1 (K) = 0.
3929 END DO
3930 ETT = 0.
3931 ETT1 = 0.
3932 ETNS = 0.
3933 ETNS1 = 0.
3934 ESNOW = 0.
3935 ESNOW1 = 0.
3936 ESNOW2 = 0.
3937
3938 ! ----------------------------------------------------------------------
3939 ! CONVERT POTENTIAL EVAP (ETP) FROM KG M-2 S-1 TO ETP1 IN M S-1
3940 ! ----------------------------------------------------------------------
3941 PRCP1 = PRCPF *0.001
3942 ETP1 = ETP * 0.001
3943 ! ----------------------------------------------------------------------
3944 ! IF ETP<0 (DOWNWARD) THEN DEWFALL (=FROSTFALL IN THIS CASE).
3945 ! ----------------------------------------------------------------------
3946 BETA = 1.0
3947 IF (ETP <= 0.0) THEN
3948 IF(ETP == 0.) BETA = 0.0
3949 DEW = -ETP1
3950 ESNOW2 = ETP1*DT
3951 ETANRG = ETP*((1.-SNCOVR)*LSUBC + SNCOVR*LSUBS)
3952 ELSE
3953 IF (ICE /= 1) THEN
3954 IF (SNCOVR < 1.) THEN
3955 CALL EVAPO (ETNS1,SMC,NSOIL,CMC,ETP1,DT,ZSOIL, &
3956 SH2O, &
3957 SMCMAX,BEXP,PC,SMCWLT,DKSAT,DWSAT, &
3958 SMCREF,SHDFAC,CMCMAX, &
3959 SMCDRY,CFACTR, &
3960 EDIR1,EC1,ET1,ETT1,SFCTMP,Q2,NROOT,RTDIS, &
3961 FXEXP)
3962 ! ----------------------------------------------------------------------------
3963 EDIR1 = EDIR1* (1. - SNCOVR)
3964 EC1 = EC1* (1. - SNCOVR)
3965 DO K = 1,NSOIL
3966 ET1 (K) = ET1 (K)* (1. - SNCOVR)
3967 END DO
3968 ETT1 = ETT1*(1.-SNCOVR)
3969 ! ETNS1 = EDIR1+ EC1+ ETT1
3970 ETNS1 = ETNS1*(1.-SNCOVR)
3971 ! ----------------------------------------------------------------------------
3972 EDIR = EDIR1*1000.
3973 EC = EC1*1000.
3974 DO K = 1,NSOIL
3975 ET (K) = ET1 (K)*1000.
3976 END DO
3977 ETT = ETT1*1000.
3978 ETNS = ETNS1*1000.
3979 ETNS = ETNS1*1000.
3980 ! ----------------------------------------------------------------------
3981
3982 ! end IF (SNCOVR .lt. 1.)
3983 END IF
3984 ! end IF (ICE .ne. 1)
3985 END IF
3986 ESNOW = ETP*SNCOVR
3987 ESNOW1 = ESNOW*0.001
3988 ESNOW2 = ESNOW1*DT
3989 ETANRG = ESNOW*LSUBS + ETNS*LSUBC
3990 ! end IF (ETP .le. 0.0)
3991 END IF
3992
3993 ! ----------------------------------------------------------------------
3994 ! IF PRECIP IS FALLING, CALCULATE HEAT FLUX FROM SNOW SFC TO NEWLY
3995 ! ACCUMULATING PRECIP. NOTE THAT THIS REFLECTS THE FLUX APPROPRIATE FOR
3996 ! THE NOT-YET-UPDATED SKIN TEMPERATURE (T1). ASSUMES TEMPERATURE OF THE
3997 ! SNOWFALL STRIKING THE GROUND IS =SFCTMP (LOWEST MODEL LEVEL AIR TEMP).
3998 ! ----------------------------------------------------------------------
3999 FLX1 = 0.0
4000 IF (SNOWNG) THEN
4001 FLX1 = CPICE * PRCP * (T1- SFCTMP)
4002 ELSE
4003 IF (PRCP > 0.0) FLX1 = CPH2O * PRCP * (T1- SFCTMP)
4004 ! ----------------------------------------------------------------------
4005 ! CALCULATE AN 'EFFECTIVE SNOW-GRND SFC TEMP' (T12) BASED ON HEAT FLUXES
4006 ! BETWEEN THE SNOW PACK AND THE SOIL AND ON NET RADIATION.
4007 ! INCLUDE FLX1 (PRECIP-SNOW SFC) AND FLX2 (FREEZING RAIN LATENT HEAT)
4008 ! FLUXES. FLX1 FROM ABOVE, FLX2 BROUGHT IN VIA COMMOM BLOCK RITE.
4009 ! FLX2 REFLECTS FREEZING RAIN LATENT HEAT FLUX USING T1 CALCULATED IN
4010 ! PENMAN.
4011 ! ----------------------------------------------------------------------
4012 END IF
4013 DSOIL = - (0.5 * ZSOIL (1))
4014 DTOT = SNOWH + DSOIL
4015 DENOM = 1.0+ DF1 / (DTOT * RR * RCH)
4016 ! surface emissivity weighted by snow cover fraction
4017 T12A = ( (FDOWN - FLX1 - FLX2 - &
4018 & ((SNCOVR*EMISSI_S)+EMISSI*(1.0-SNCOVR))*SIGMA *T24)/RCH &
4019 & + TH2 - SFCTMP - ETANRG/RCH ) / RR
4020
4021 T12B = DF1 * STC (1) / (DTOT * RR * RCH)
4022
4023 ! ----------------------------------------------------------------------
4024 ! IF THE 'EFFECTIVE SNOW-GRND SFC TEMP' IS AT OR BELOW FREEZING, NO SNOW
4025 ! MELT WILL OCCUR. SET THE SKIN TEMP TO THIS EFFECTIVE TEMP. REDUCE
4026 ! (BY SUBLIMINATION ) OR INCREASE (BY FROST) THE DEPTH OF THE SNOWPACK,
4027 ! DEPENDING ON SIGN OF ETP.
4028 ! UPDATE SOIL HEAT FLUX (SSOIL) USING NEW SKIN TEMPERATURE (T1)
4029 ! SINCE NO SNOWMELT, SET ACCUMULATED SNOWMELT TO ZERO, SET 'EFFECTIVE'
4030 ! PRECIP FROM SNOWMELT TO ZERO, SET PHASE-CHANGE HEAT FLUX FROM SNOWMELT
4031 ! TO ZERO.
4032 ! ----------------------------------------------------------------------
4033 ! SUB-FREEZING BLOCK
4034 ! ----------------------------------------------------------------------
4035 T12 = (SFCTMP + T12A + T12B) / DENOM
4036 IF (T12 <= TFREEZ) THEN
4037 T1 = T12
4038 SSOIL = DF1 * (T1- STC (1)) / DTOT
4039 ! ESD = MAX (0.0, ESD- ETP2)
4040 ESD = MAX(0.0, ESD-ESNOW2)
4041 FLX3 = 0.0
4042 EX = 0.0
4043
4044 SNOMLT = 0.0
4045 ! ----------------------------------------------------------------------
4046 ! IF THE 'EFFECTIVE SNOW-GRND SFC TEMP' IS ABOVE FREEZING, SNOW MELT
4047 ! WILL OCCUR. CALL THE SNOW MELT RATE,EX AND AMT, SNOMLT. REVISE THE
4048 ! EFFECTIVE SNOW DEPTH. REVISE THE SKIN TEMP BECAUSE IT WOULD HAVE CHGD
4049 ! DUE TO THE LATENT HEAT RELEASED BY THE MELTING. CALC THE LATENT HEAT
4050 ! RELEASED, FLX3. SET THE EFFECTIVE PRECIP, PRCP1 TO THE SNOW MELT RATE,
4051 ! EX FOR USE IN SMFLX. ADJUSTMENT TO T1 TO ACCOUNT FOR SNOW PATCHES.
4052 ! CALCULATE QSAT VALID AT FREEZING POINT. NOTE THAT ESAT (SATURATION
4053 ! VAPOR PRESSURE) VALUE OF 6.11E+2 USED HERE IS THAT VALID AT FRZZING
4054 ! POINT. NOTE THAT ETP FROM CALL PENMAN IN SFLX IS IGNORED HERE IN
4055 ! FAVOR OF BULK ETP OVER 'OPEN WATER' AT FREEZING TEMP.
4056 ! UPDATE SOIL HEAT FLUX (S) USING NEW SKIN TEMPERATURE (T1)
4057 ! ----------------------------------------------------------------------
4058 ! ABOVE FREEZING BLOCK
4059 ! ----------------------------------------------------------------------
4060 ELSE
4061 T1 = TFREEZ * SNCOVR ** SNOEXP + T12 * (1.0- SNCOVR ** SNOEXP)
4062 BETA = 1.0
4063
4064 ! ----------------------------------------------------------------------
4065 ! IF POTENTIAL EVAP (SUBLIMATION) GREATER THAN DEPTH OF SNOWPACK.
4066 ! BETA<1
4067 ! SNOWPACK HAS SUBLIMATED AWAY, SET DEPTH TO ZERO.
4068 ! ----------------------------------------------------------------------
4069 SSOIL = DF1 * (T1- STC (1)) / DTOT
4070 IF (ESD-ESNOW2 <= ESDMIN) THEN
4071 ESD = 0.0
4072 EX = 0.0
4073 SNOMLT = 0.0
4074 ! ----------------------------------------------------------------------
4075 ! SUBLIMATION LESS THAN DEPTH OF SNOWPACK
4076 ! SNOWPACK (ESD) REDUCED BY ESNOW2 (DEPTH OF SUBLIMATED SNOW)
4077 ! ----------------------------------------------------------------------
4078 ELSE
4079 ESD = ESD-ESNOW2
4080 ETP3 = ETP * LSUBC
4081 SEH = RCH * (T1- TH2)
4082 T14 = T1* T1
4083 T14 = T14* T14
4084 FLX3 = FDOWN - FLX1 - FLX2 - &
4085 ((SNCOVR*EMISSI_S)+EMISSI*(1-SNCOVR))*SIGMA*T14 - &
4086 SSOIL - SEH - ETANRG
4087 IF (FLX3 <= 0.0) FLX3 = 0.0
4088 ! ----------------------------------------------------------------------
4089 ! SNOWMELT REDUCTION DEPENDING ON SNOW COVER
4090 ! ----------------------------------------------------------------------
4091 EX = FLX3*0.001/ LSUBF
4092
4093 ! ----------------------------------------------------------------------
4094 ! ESDMIN REPRESENTS A SNOWPACK DEPTH THRESHOLD VALUE BELOW WHICH WE
4095 ! CHOOSE NOT TO RETAIN ANY SNOWPACK, AND INSTEAD INCLUDE IT IN SNOWMELT.
4096 ! ----------------------------------------------------------------------
4097 SNOMLT = EX * DT
4098 IF (ESD- SNOMLT >= ESDMIN) THEN
4099 ESD = ESD- SNOMLT
4100 ! ----------------------------------------------------------------------
4101 ! SNOWMELT EXCEEDS SNOW DEPTH
4102 ! ----------------------------------------------------------------------
4103 ELSE
4104 EX = ESD / DT
4105 FLX3 = EX *1000.0* LSUBF
4106 SNOMLT = ESD
4107
4108 ESD = 0.0
4109 ! ----------------------------------------------------------------------
4110 ! END OF 'ESD .LE. ETP2' IF-BLOCK
4111 ! ----------------------------------------------------------------------
4112 END IF
4113 END IF
4114
4115 ! ----------------------------------------------------------------------
4116 ! END OF 'T12 .LE. TFREEZ' IF-BLOCK
4117 ! ----------------------------------------------------------------------
4118 PRCP1 = PRCP1+ EX
4119
4120 ! ----------------------------------------------------------------------
4121 ! SET THE EFFECTIVE POTNL EVAPOTRANSP (ETP1) TO ZERO SINCE THIS IS SNOW
4122 ! CASE, SO SURFACE EVAP NOT CALCULATED FROM EDIR, EC, OR ETT IN SMFLX
4123 ! (BELOW).
4124 ! IF SEAICE (ICE=1) SKIP CALL TO SMFLX.
4125 ! SMFLX RETURNS UPDATED SOIL MOISTURE VALUES. IN THIS, THE SNOW PACK
4126 ! CASE, ETA1 IS NOT USED IN CALCULATION OF EVAP.
4127 ! ----------------------------------------------------------------------
4128 END IF
4129 IF (ICE /= 1) THEN
4130 CALL SMFLX (SMC,NSOIL,CMC,DT,PRCP1,ZSOIL, &
4131 SH2O,SLOPE,KDT,FRZFACT, &
4132 SMCMAX,BEXP,SMCWLT,DKSAT,DWSAT, &
4133 SHDFAC,CMCMAX, &
4134 RUNOFF1,RUNOFF2,RUNOFF3, &
4135 EDIR1,EC1,ET1, &
4136 DRIP)
4137 ! ----------------------------------------------------------------------
4138 ! BEFORE CALL SHFLX IN THIS SNOWPACK CASE, SET ZZ1 AND YY ARGUMENTS TO
4139 ! SPECIAL VALUES THAT ENSURE THAT GROUND HEAT FLUX CALCULATED IN SHFLX
4140 ! MATCHES THAT ALREADY COMPUTER FOR BELOW THE SNOWPACK, THUS THE SFC
4141 ! HEAT FLUX TO BE COMPUTED IN SHFLX WILL EFFECTIVELY BE THE FLUX AT THE
4142 ! SNOW TOP SURFACE. T11 IS A DUMMY ARGUEMENT SO WE WILL NOT USE THE
4143 ! SKIN TEMP VALUE AS REVISED BY SHFLX.
4144 ! ----------------------------------------------------------------------
4145 END IF
4146 ZZ1 = 1.0
4147 YY = STC (1) -0.5* SSOIL * ZSOIL (1)* ZZ1/ DF1
4148
4149 ! ----------------------------------------------------------------------
4150 ! SHFLX WILL CALC/UPDATE THE SOIL TEMPS. NOTE: THE SUB-SFC HEAT FLUX
4151 ! (SSOIL1) AND THE SKIN TEMP (T11) OUTPUT FROM THIS SHFLX CALL ARE NOT
4152 ! USED IN ANY SUBSEQUENT CALCULATIONS. RATHER, THEY ARE DUMMY VARIABLES
4153 ! HERE IN THE SNOPAC CASE, SINCE THE SKIN TEMP AND SUB-SFC HEAT FLUX ARE
4154 ! UPDATED INSTEAD NEAR THE BEGINNING OF THE CALL TO SNOPAC.
4155 ! ----------------------------------------------------------------------
4156 T11 = T1
4157 CALL SHFLX (SSOIL1,STC,SMC,SMCMAX,NSOIL,T11,DT,YY,ZZ1,ZSOIL, &
4158 TBOT,ZBOT,SMCWLT,PSISAT,SH2O,BEXP,F1,DF1,ICE, &
4159 QUARTZ,CSOIL,VEGTYP)
4160
4161 ! ----------------------------------------------------------------------
4162 ! SNOW DEPTH AND DENSITY ADJUSTMENT BASED ON SNOW COMPACTION. YY IS
4163 ! ASSUMED TO BE THE SOIL TEMPERTURE AT THE TOP OF THE SOIL COLUMN.
4164 ! ----------------------------------------------------------------------
4165 IF (ESD > 0.) THEN
4166 CALL SNOWPACK (ESD,DT,SNOWH,SNDENS,T1,YY)
4167 ELSE
4168 ESD = 0.
4169 SNOWH = 0.
4170 SNDENS = 0.
4171 SNCOND = 1.
4172 SNCOVR = 0.
4173 END IF
4174 ! ----------------------------------------------------------------------
4175 END SUBROUTINE SNOPAC
4176 ! ----------------------------------------------------------------------
4177
4178
4179 SUBROUTINE SNOWPACK (ESD,DTSEC,SNOWH,SNDENS,TSNOW,TSOIL)
4180
4181 ! ----------------------------------------------------------------------
4182 ! SUBROUTINE SNOWPACK
4183 ! ----------------------------------------------------------------------
4184 ! CALCULATE COMPACTION OF SNOWPACK UNDER CONDITIONS OF INCREASING SNOW
4185 ! DENSITY, AS OBTAINED FROM AN APPROXIMATE SOLUTION OF E. ANDERSON'S
4186 ! DIFFERENTIAL EQUATION (3.29), NOAA TECHNICAL REPORT NWS 19, BY VICTOR
4187 ! KOREN, 03/25/95.
4188 ! ----------------------------------------------------------------------
4189 ! ESD WATER EQUIVALENT OF SNOW (M)
4190 ! DTSEC TIME STEP (SEC)
4191 ! SNOWH SNOW DEPTH (M)
4192 ! SNDENS SNOW DENSITY (G/CM3=DIMENSIONLESS FRACTION OF H2O DENSITY)
4193 ! TSNOW SNOW SURFACE TEMPERATURE (K)
4194 ! TSOIL SOIL SURFACE TEMPERATURE (K)
4195
4196 ! SUBROUTINE WILL RETURN NEW VALUES OF SNOWH AND SNDENS
4197 ! ----------------------------------------------------------------------
4198 IMPLICIT NONE
4199
4200 INTEGER :: IPOL, J
4201 REAL, INTENT(IN) :: ESD, DTSEC,TSNOW,TSOIL
4202 REAL, INTENT(INOUT) :: SNOWH, SNDENS
4203 REAL :: BFAC,DSX,DTHR,DW,SNOWHC,PEXP, &
4204 TAVGC,TSNOWC,TSOILC,ESDC,ESDCX
4205 REAL, PARAMETER :: C1 = 0.01, C2 = 21.0, G = 9.81, &
4206 KN = 4000.0
4207 ! ----------------------------------------------------------------------
4208 ! CONVERSION INTO SIMULATION UNITS
4209 ! ----------------------------------------------------------------------
4210 SNOWHC = SNOWH *100.
4211 ESDC = ESD *100.
4212 DTHR = DTSEC /3600.
4213 TSNOWC = TSNOW -273.15
4214 TSOILC = TSOIL -273.15
4215
4216 ! ----------------------------------------------------------------------
4217 ! CALCULATING OF AVERAGE TEMPERATURE OF SNOW PACK
4218 ! ----------------------------------------------------------------------
4219 ! ----------------------------------------------------------------------
4220 ! CALCULATING OF SNOW DEPTH AND DENSITY AS A RESULT OF COMPACTION
4221 ! SNDENS=DS0*(EXP(BFAC*ESD)-1.)/(BFAC*ESD)
4222 ! BFAC=DTHR*C1*EXP(0.08*TAVGC-C2*DS0)
4223 ! NOTE: BFAC*ESD IN SNDENS EQN ABOVE HAS TO BE CAREFULLY TREATED
4224 ! NUMERICALLY BELOW:
4225 ! C1 IS THE FRACTIONAL INCREASE IN DENSITY (1/(CM*HR))
4226 ! C2 IS A CONSTANT (CM3/G) KOJIMA ESTIMATED AS 21 CMS/G
4227 ! ----------------------------------------------------------------------
4228 TAVGC = 0.5* (TSNOWC + TSOILC)
4229 IF (ESDC > 1.E-2) THEN
4230 ESDCX = ESDC
4231 ELSE
4232 ESDCX = 1.E-2
4233 END IF
4234
4235 ! DSX = SNDENS*((DEXP(BFAC*ESDC)-1.)/(BFAC*ESDC))
4236 ! ----------------------------------------------------------------------
4237 ! THE FUNCTION OF THE FORM (e**x-1)/x IMBEDDED IN ABOVE EXPRESSION
4238 ! FOR DSX WAS CAUSING NUMERICAL DIFFICULTIES WHEN THE DENOMINATOR "x"
4239 ! (I.E. BFAC*ESDC) BECAME ZERO OR APPROACHED ZERO (DESPITE THE FACT THAT
4240 ! THE ANALYTICAL FUNCTION (e**x-1)/x HAS A WELL DEFINED LIMIT AS
4241 ! "x" APPROACHES ZERO), HENCE BELOW WE REPLACE THE (e**x-1)/x
4242 ! EXPRESSION WITH AN EQUIVALENT, NUMERICALLY WELL-BEHAVED
4243 ! POLYNOMIAL EXPANSION.
4244
4245 ! NUMBER OF TERMS OF POLYNOMIAL EXPANSION, AND HENCE ITS ACCURACY,
4246 ! IS GOVERNED BY ITERATION LIMIT "IPOL".
4247 ! IPOL GREATER THAN 9 ONLY MAKES A DIFFERENCE ON DOUBLE
4248 ! PRECISION (RELATIVE ERRORS GIVEN IN PERCENT %).
4249 ! IPOL=9, FOR REL.ERROR <~ 1.6 E-6 % (8 SIGNIFICANT DIGITS)
4250 ! IPOL=8, FOR REL.ERROR <~ 1.8 E-5 % (7 SIGNIFICANT DIGITS)
4251 ! IPOL=7, FOR REL.ERROR <~ 1.8 E-4 % ...
4252 ! ----------------------------------------------------------------------
4253 BFAC = DTHR * C1* EXP (0.08* TAVGC - C2* SNDENS)
4254 IPOL = 4
4255 PEXP = 0.
4256 ! PEXP = (1. + PEXP)*BFAC*ESDC/REAL(J+1)
4257 DO J = IPOL,1, -1
4258 PEXP = (1. + PEXP)* BFAC * ESDCX / REAL (J +1)
4259 END DO
4260
4261 PEXP = PEXP + 1.
4262 ! ----------------------------------------------------------------------
4263 ! ABOVE LINE ENDS POLYNOMIAL SUBSTITUTION
4264 ! ----------------------------------------------------------------------
4265 ! END OF KOREAN FORMULATION
4266
4267 ! BASE FORMULATION (COGLEY ET AL., 1990)
4268 ! CONVERT DENSITY FROM G/CM3 TO KG/M3
4269 ! DSM=SNDENS*1000.0
4270
4271 ! DSX=DSM+DTSEC*0.5*DSM*G*ESD/
4272 ! & (1E7*EXP(-0.02*DSM+KN/(TAVGC+273.16)-14.643))
4273
4274 ! & CONVERT DENSITY FROM KG/M3 TO G/CM3
4275 ! DSX=DSX/1000.0
4276
4277 ! END OF COGLEY ET AL. FORMULATION
4278
4279 ! ----------------------------------------------------------------------
4280 ! SET UPPER/LOWER LIMIT ON SNOW DENSITY
4281 ! ----------------------------------------------------------------------
4282 DSX = SNDENS * (PEXP)
4283 IF (DSX > 0.40) DSX = 0.40
4284 IF (DSX < 0.05) DSX = 0.05
4285 ! ----------------------------------------------------------------------
4286 ! UPDATE OF SNOW DEPTH AND DENSITY DEPENDING ON LIQUID WATER DURING
4287 ! SNOWMELT. ASSUMED THAT 13% OF LIQUID WATER CAN BE STORED IN SNOW PER
4288 ! DAY DURING SNOWMELT TILL SNOW DENSITY 0.40.
4289 ! ----------------------------------------------------------------------
4290 SNDENS = DSX
4291 IF (TSNOWC >= 0.) THEN
4292 DW = 0.13* DTHR /24.
4293 SNDENS = SNDENS * (1. - DW) + DW
4294 IF (SNDENS >= 0.40) SNDENS = 0.40
4295 ! ----------------------------------------------------------------------
4296 ! CALCULATE SNOW DEPTH (CM) FROM SNOW WATER EQUIVALENT AND SNOW DENSITY.
4297 ! CHANGE SNOW DEPTH UNITS TO METERS
4298 ! ----------------------------------------------------------------------
4299 END IF
4300 SNOWHC = ESDC / SNDENS
4301 SNOWH = SNOWHC *0.01
4302
4303 ! ----------------------------------------------------------------------
4304 END SUBROUTINE SNOWPACK
4305 ! ----------------------------------------------------------------------
4306
4307 SUBROUTINE SNOWZ0 (SNCOVR,Z0, Z0BRD)
4308
4309 ! ----------------------------------------------------------------------
4310 ! SUBROUTINE SNOWZ0
4311 ! ----------------------------------------------------------------------
4312 ! CALCULATE TOTAL ROUGHNESS LENGTH OVER SNOW
4313 ! SNCOVR FRACTIONAL SNOW COVER
4314 ! Z0 ROUGHNESS LENGTH (m)
4315 ! Z0S SNOW ROUGHNESS LENGTH:=0.001 (m)
4316 ! ----------------------------------------------------------------------
4317 IMPLICIT NONE
4318 REAL, INTENT(IN) :: SNCOVR, Z0BRD
4319 REAL, INTENT(OUT) :: Z0
4320 REAL, PARAMETER :: Z0S=0.001
4321
4322 !m Z0 = (1.- SNCOVR)* Z0BRD + SNCOVR * Z0S
4323 Z0 = Z0BRD
4324 ! ----------------------------------------------------------------------
4325 END SUBROUTINE SNOWZ0
4326 ! ----------------------------------------------------------------------
4327
4328
4329 SUBROUTINE SNOW_NEW (TEMP,NEWSN,SNOWH,SNDENS)
4330
4331 ! ----------------------------------------------------------------------
4332 ! SUBROUTINE SNOW_NEW
4333 ! ----------------------------------------------------------------------
4334 ! CALCULATE SNOW DEPTH AND DENSITITY TO ACCOUNT FOR THE NEW SNOWFALL.
4335 ! NEW VALUES OF SNOW DEPTH & DENSITY RETURNED.
4336
4337 ! TEMP AIR TEMPERATURE (K)
4338 ! NEWSN NEW SNOWFALL (M)
4339 ! SNOWH SNOW DEPTH (M)
4340 ! SNDENS SNOW DENSITY (G/CM3=DIMENSIONLESS FRACTION OF H2O DENSITY)
4341 ! ----------------------------------------------------------------------
4342 IMPLICIT NONE
4343 REAL, INTENT(IN) :: NEWSN, TEMP
4344 REAL, INTENT(INOUT) :: SNDENS, SNOWH
4345 REAL :: DSNEW, HNEWC, SNOWHC,NEWSNC,TEMPC
4346
4347 ! ----------------------------------------------------------------------
4348 ! CONVERSION INTO SIMULATION UNITS
4349 ! ----------------------------------------------------------------------
4350 SNOWHC = SNOWH *100.
4351 NEWSNC = NEWSN *100.
4352
4353 ! ----------------------------------------------------------------------
4354 ! CALCULATING NEW SNOWFALL DENSITY DEPENDING ON TEMPERATURE
4355 ! EQUATION FROM GOTTLIB L. 'A GENERAL RUNOFF MODEL FOR SNOWCOVERED
4356 ! AND GLACIERIZED BASIN', 6TH NORDIC HYDROLOGICAL CONFERENCE,
4357 ! VEMADOLEN, SWEDEN, 1980, 172-177PP.
4358 !-----------------------------------------------------------------------
4359 TEMPC = TEMP -273.15
4360 IF (TEMPC <= -15.) THEN
4361 DSNEW = 0.05
4362 ELSE
4363 DSNEW = 0.05+0.0017* (TEMPC +15.)**1.5
4364 END IF
4365 ! ----------------------------------------------------------------------
4366 ! ADJUSTMENT OF SNOW DENSITY DEPENDING ON NEW SNOWFALL
4367 ! ----------------------------------------------------------------------
4368 HNEWC = NEWSNC / DSNEW
4369 IF (SNOWHC + HNEWC .LT. 1.0E-3) THEN
4370 SNDENS = MAX(DSNEW,SNDENS)
4371 ELSE
4372 SNDENS = (SNOWHC * SNDENS + HNEWC * DSNEW)/ (SNOWHC + HNEWC)
4373 ENDIF
4374 SNOWHC = SNOWHC + HNEWC
4375 SNOWH = SNOWHC *0.01
4376
4377 ! ----------------------------------------------------------------------
4378 END SUBROUTINE SNOW_NEW
4379 ! ----------------------------------------------------------------------
4380
4381 SUBROUTINE SRT (RHSTT,EDIR,ET,SH2O,SH2OA,NSOIL,PCPDRP, &
4382 ZSOIL,DWSAT,DKSAT,SMCMAX,BEXP,RUNOFF1, &
4383 RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZX,SICE,AI,BI,CI)
4384
4385 ! ----------------------------------------------------------------------
4386 ! SUBROUTINE SRT
4387 ! ----------------------------------------------------------------------
4388 ! CALCULATE THE RIGHT HAND SIDE OF THE TIME TENDENCY TERM OF THE SOIL
4389 ! WATER DIFFUSION EQUATION. ALSO TO COMPUTE ( PREPARE ) THE MATRIX
4390 ! COEFFICIENTS FOR THE TRI-DIAGONAL MATRIX OF THE IMPLICIT TIME SCHEME.
4391 ! ----------------------------------------------------------------------
4392 IMPLICIT NONE
4393 INTEGER, INTENT(IN) :: NSOIL
4394 INTEGER :: IALP1, IOHINF, J, JJ, K, KS
4395 REAL, INTENT(IN) :: BEXP, DKSAT, DT, DWSAT, EDIR, FRZX, &
4396 KDT, PCPDRP, SLOPE, SMCMAX, SMCWLT
4397 REAL, INTENT(OUT) :: RUNOFF1, RUNOFF2
4398 REAL, DIMENSION(1:NSOIL), INTENT(IN) :: ET, SH2O, SH2OA, SICE, &
4399 ZSOIL
4400 REAL, DIMENSION(1:NSOIL), INTENT(OUT) :: RHSTT
4401 REAL, DIMENSION(1:NSOIL), INTENT(OUT) :: AI, BI, CI
4402 REAL, DIMENSION(1:NSOIL) :: DMAX
4403 REAL :: ACRT, DD, DDT, DDZ, DDZ2, DENOM, &
4404 DENOM2,DICE, DSMDZ, DSMDZ2, DT1, &
4405 FCR,INFMAX,MXSMC,MXSMC2,NUMER,PDDUM, &
4406 PX, SICEMAX,SLOPX, SMCAV, SSTT, &
4407 SUM, VAL, WCND, WCND2, WDF, WDF2
4408 INTEGER, PARAMETER :: CVFRZ = 3
4409
4410 ! ----------------------------------------------------------------------
4411 ! FROZEN GROUND VERSION:
4412 ! REFERENCE FROZEN GROUND PARAMETER, CVFRZ, IS A SHAPE PARAMETER OF
4413 ! AREAL DISTRIBUTION FUNCTION OF SOIL ICE CONTENT WHICH EQUALS 1/CV.
4414 ! CV IS A COEFFICIENT OF SPATIAL VARIATION OF SOIL ICE CONTENT. BASED
4415 ! ON FIELD DATA CV DEPENDS ON AREAL MEAN OF FROZEN DEPTH, AND IT CLOSE
4416 ! TO CONSTANT = 0.6 IF AREAL MEAN FROZEN DEPTH IS ABOVE 20 CM. THAT IS
4417 ! WHY PARAMETER CVFRZ = 3 (INT{1/0.6*0.6}).
4418 ! CURRENT LOGIC DOESN'T ALLOW CVFRZ BE BIGGER THAN 3
4419 ! ----------------------------------------------------------------------
4420
4421 ! ----------------------------------------------------------------------
4422 ! DETERMINE RAINFALL INFILTRATION RATE AND RUNOFF. INCLUDE THE
4423 ! INFILTRATION FORMULE FROM SCHAAKE AND KOREN MODEL.
4424 ! MODIFIED BY Q DUAN
4425 ! ----------------------------------------------------------------------
4426 ! ----------------------------------------------------------------------
4427 ! LET SICEMAX BE THE GREATEST, IF ANY, FROZEN WATER CONTENT WITHIN SOIL
4428 ! LAYERS.
4429 ! ----------------------------------------------------------------------
4430 IOHINF = 1
4431 SICEMAX = 0.0
4432 DO KS = 1,NSOIL
4433 IF (SICE (KS) > SICEMAX) SICEMAX = SICE (KS)
4434 ! ----------------------------------------------------------------------
4435 ! DETERMINE RAINFALL INFILTRATION RATE AND RUNOFF
4436 ! ----------------------------------------------------------------------
4437 END DO
4438 PDDUM = PCPDRP
4439 RUNOFF1 = 0.0
4440
4441 ! ----------------------------------------------------------------------
4442 ! MODIFIED BY Q. DUAN, 5/16/94
4443 ! ----------------------------------------------------------------------
4444 ! IF (IOHINF == 1) THEN
4445
4446 IF (PCPDRP /= 0.0) THEN
4447 DT1 = DT /86400.
4448 SMCAV = SMCMAX - SMCWLT
4449
4450 ! ----------------------------------------------------------------------
4451 ! FROZEN GROUND VERSION:
4452 ! ----------------------------------------------------------------------
4453 DMAX (1)= - ZSOIL (1)* SMCAV
4454
4455 DICE = - ZSOIL (1) * SICE (1)
4456 DMAX (1)= DMAX (1)* (1.0- (SH2OA (1) + SICE (1) - SMCWLT)/ &
4457 SMCAV)
4458
4459 DD = DMAX (1)
4460
4461 ! ----------------------------------------------------------------------
4462 ! FROZEN GROUND VERSION:
4463 ! ----------------------------------------------------------------------
4464 DO KS = 2,NSOIL
4465
4466 DICE = DICE+ ( ZSOIL (KS -1) - ZSOIL (KS) ) * SICE (KS)
4467 DMAX (KS) = (ZSOIL (KS -1) - ZSOIL (KS))* SMCAV
4468 DMAX (KS) = DMAX (KS)* (1.0- (SH2OA (KS) + SICE (KS) &
4469 - SMCWLT)/ SMCAV)
4470 DD = DD+ DMAX (KS)
4471 ! ----------------------------------------------------------------------
4472 ! VAL = (1.-EXP(-KDT*SQRT(DT1)))
4473 ! IN BELOW, REMOVE THE SQRT IN ABOVE
4474 ! ----------------------------------------------------------------------
4475 END DO
4476 VAL = (1. - EXP ( - KDT * DT1))
4477 DDT = DD * VAL
4478 PX = PCPDRP * DT
4479 IF (PX < 0.0) PX = 0.0
4480
4481 ! ----------------------------------------------------------------------
4482 ! FROZEN GROUND VERSION:
4483 ! REDUCTION OF INFILTRATION BASED ON FROZEN GROUND PARAMETERS
4484 ! ----------------------------------------------------------------------
4485 INFMAX = (PX * (DDT / (PX + DDT)))/ DT
4486 FCR = 1.
4487 IF (DICE > 1.E-2) THEN
4488 ACRT = CVFRZ * FRZX / DICE
4489 SUM = 1.
4490 IALP1 = CVFRZ - 1
4491 DO J = 1,IALP1
4492 K = 1
4493 DO JJ = J +1,IALP1
4494 K = K * JJ
4495 END DO
4496 SUM = SUM + (ACRT ** ( CVFRZ - J)) / FLOAT (K)
4497 END DO
4498 FCR = 1. - EXP ( - ACRT) * SUM
4499 END IF
4500
4501 ! ----------------------------------------------------------------------
4502 ! CORRECTION OF INFILTRATION LIMITATION:
4503 ! IF INFMAX .LE. HYDROLIC CONDUCTIVITY ASSIGN INFMAX THE VALUE OF
4504 ! HYDROLIC CONDUCTIVITY
4505 ! ----------------------------------------------------------------------
4506 ! MXSMC = MAX ( SH2OA(1), SH2OA(2) )
4507 INFMAX = INFMAX * FCR
4508
4509 MXSMC = SH2OA (1)
4510 CALL WDFCND (WDF,WCND,MXSMC,SMCMAX,BEXP,DKSAT,DWSAT, &
4511 SICEMAX)
4512 INFMAX = MAX (INFMAX,WCND)
4513
4514 INFMAX = MIN (INFMAX,PX)
4515 IF (PCPDRP > INFMAX) THEN
4516 RUNOFF1 = PCPDRP - INFMAX
4517 PDDUM = INFMAX
4518 END IF
4519 ! ----------------------------------------------------------------------
4520 ! TO AVOID SPURIOUS DRAINAGE BEHAVIOR, 'UPSTREAM DIFFERENCING' IN LINE
4521 ! BELOW REPLACED WITH NEW APPROACH IN 2ND LINE:
4522 ! 'MXSMC = MAX(SH2OA(1), SH2OA(2))'
4523 ! ----------------------------------------------------------------------
4524 END IF
4525
4526 MXSMC = SH2OA (1)
4527 CALL WDFCND (WDF,WCND,MXSMC,SMCMAX,BEXP,DKSAT,DWSAT, &
4528 SICEMAX)
4529 ! ----------------------------------------------------------------------
4530 ! CALC THE MATRIX COEFFICIENTS AI, BI, AND CI FOR THE TOP LAYER
4531 ! ----------------------------------------------------------------------
4532 DDZ = 1. / ( - .5 * ZSOIL (2) )
4533 AI (1) = 0.0
4534 BI (1) = WDF * DDZ / ( - ZSOIL (1) )
4535
4536 ! ----------------------------------------------------------------------
4537 ! CALC RHSTT FOR THE TOP LAYER AFTER CALC'NG THE VERTICAL SOIL MOISTURE
4538 ! GRADIENT BTWN THE TOP AND NEXT TO TOP LAYERS.
4539 ! ----------------------------------------------------------------------
4540 CI (1) = - BI (1)
4541 DSMDZ = ( SH2O (1) - SH2O (2) ) / ( - .5 * ZSOIL (2) )
4542 RHSTT (1) = (WDF * DSMDZ + WCND- PDDUM + EDIR + ET (1))/ ZSOIL (1)
4543
4544 ! ----------------------------------------------------------------------
4545 ! INITIALIZE DDZ2
4546 ! ----------------------------------------------------------------------
4547 SSTT = WDF * DSMDZ + WCND+ EDIR + ET (1)
4548
4549 ! ----------------------------------------------------------------------
4550 ! LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABV PROCESS
4551 ! ----------------------------------------------------------------------
4552 DDZ2 = 0.0
4553 DO K = 2,NSOIL
4554 DENOM2 = (ZSOIL (K -1) - ZSOIL (K))
4555 IF (K /= NSOIL) THEN
4556
4557 ! ----------------------------------------------------------------------
4558 ! AGAIN, TO AVOID SPURIOUS DRAINAGE BEHAVIOR, 'UPSTREAM DIFFERENCING' IN
4559 ! LINE BELOW REPLACED WITH NEW APPROACH IN 2ND LINE:
4560 ! 'MXSMC2 = MAX (SH2OA(K), SH2OA(K+1))'
4561 ! ----------------------------------------------------------------------
4562 SLOPX = 1.
4563
4564 MXSMC2 = SH2OA (K)
4565 CALL WDFCND (WDF2,WCND2,MXSMC2,SMCMAX,BEXP,DKSAT,DWSAT, &
4566 SICEMAX)
4567 ! -----------------------------------------------------------------------
4568 ! CALC SOME PARTIAL PRODUCTS FOR LATER USE IN CALC'NG RHSTT
4569 ! ----------------------------------------------------------------------
4570 DENOM = (ZSOIL (K -1) - ZSOIL (K +1))
4571
4572 ! ----------------------------------------------------------------------
4573 ! CALC THE MATRIX COEF, CI, AFTER CALC'NG ITS PARTIAL PRODUCT
4574 ! ----------------------------------------------------------------------
4575 DSMDZ2 = (SH2O (K) - SH2O (K +1)) / (DENOM * 0.5)
4576 DDZ2 = 2.0 / DENOM
4577 CI (K) = - WDF2 * DDZ2 / DENOM2
4578
4579 ELSE
4580 ! ----------------------------------------------------------------------
4581 ! SLOPE OF BOTTOM LAYER IS INTRODUCED
4582 ! ----------------------------------------------------------------------
4583
4584 ! ----------------------------------------------------------------------
4585 ! RETRIEVE THE SOIL WATER DIFFUSIVITY AND HYDRAULIC CONDUCTIVITY FOR
4586 ! THIS LAYER
4587 ! ----------------------------------------------------------------------
4588 SLOPX = SLOPE
4589 CALL WDFCND (WDF2,WCND2,SH2OA (NSOIL),SMCMAX,BEXP,DKSAT,DWSAT, &
4590 SICEMAX)
4591
4592 ! ----------------------------------------------------------------------
4593 ! CALC A PARTIAL PRODUCT FOR LATER USE IN CALC'NG RHSTT
4594 ! ----------------------------------------------------------------------
4595
4596 ! ----------------------------------------------------------------------
4597 ! SET MATRIX COEF CI TO ZERO
4598 ! ----------------------------------------------------------------------
4599 DSMDZ2 = 0.0
4600 CI (K) = 0.0
4601 ! ----------------------------------------------------------------------
4602 ! CALC RHSTT FOR THIS LAYER AFTER CALC'NG ITS NUMERATOR
4603 ! ----------------------------------------------------------------------
4604 END IF
4605 NUMER = (WDF2 * DSMDZ2) + SLOPX * WCND2- (WDF * DSMDZ) &
4606 - WCND+ ET (K)
4607
4608 ! ----------------------------------------------------------------------
4609 ! CALC MATRIX COEFS, AI, AND BI FOR THIS LAYER
4610 ! ----------------------------------------------------------------------
4611 RHSTT (K) = NUMER / ( - DENOM2)
4612 AI (K) = - WDF * DDZ / DENOM2
4613
4614 ! ----------------------------------------------------------------------
4615 ! RESET VALUES OF WDF, WCND, DSMDZ, AND DDZ FOR LOOP TO NEXT LYR
4616 ! RUNOFF2: SUB-SURFACE OR BASEFLOW RUNOFF
4617 ! ----------------------------------------------------------------------
4618 BI (K) = - ( AI (K) + CI (K) )
4619 IF (K .eq. NSOIL) THEN
4620 RUNOFF2 = SLOPX * WCND2
4621 END IF
4622 IF (K .ne. NSOIL) THEN
4623 WDF = WDF2
4624 WCND = WCND2
4625 DSMDZ = DSMDZ2
4626 DDZ = DDZ2
4627 END IF
4628 END DO
4629 ! ----------------------------------------------------------------------
4630 END SUBROUTINE SRT
4631 ! ----------------------------------------------------------------------
4632
4633 SUBROUTINE SSTEP (SH2OOUT,SH2OIN,CMC,RHSTT,RHSCT,DT, &
4634 NSOIL,SMCMAX,CMCMAX,RUNOFF3,ZSOIL,SMC,SICE, &
4635 AI,BI,CI)
4636
4637 ! ----------------------------------------------------------------------
4638 ! SUBROUTINE SSTEP
4639 ! ----------------------------------------------------------------------
4640 ! CALCULATE/UPDATE SOIL MOISTURE CONTENT VALUES AND CANOPY MOISTURE
4641 ! CONTENT VALUES.
4642 ! ----------------------------------------------------------------------
4643 IMPLICIT NONE
4644 INTEGER, INTENT(IN) :: NSOIL
4645 INTEGER :: I, K, KK11
4646
4647 REAL, INTENT(IN) :: CMCMAX, DT, SMCMAX
4648 REAL, INTENT(OUT) :: RUNOFF3
4649 REAL, INTENT(INOUT) :: CMC
4650 REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SH2OIN, SICE, ZSOIL
4651 REAL, DIMENSION(1:NSOIL), INTENT(OUT) :: SH2OOUT
4652 REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: RHSTT, SMC
4653 REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: AI, BI, CI
4654 REAL, DIMENSION(1:NSOIL) :: RHSTTin
4655 REAL, DIMENSION(1:NSOIL) :: CIin
4656 REAL :: DDZ, RHSCT, STOT, WPLUS
4657
4658 ! ----------------------------------------------------------------------
4659 ! CREATE 'AMOUNT' VALUES OF VARIABLES TO BE INPUT TO THE
4660 ! TRI-DIAGONAL MATRIX ROUTINE.
4661 ! ----------------------------------------------------------------------
4662 DO K = 1,NSOIL
4663 RHSTT (K) = RHSTT (K) * DT
4664 AI (K) = AI (K) * DT
4665 BI (K) = 1. + BI (K) * DT
4666 CI (K) = CI (K) * DT
4667 END DO
4668 ! ----------------------------------------------------------------------
4669 ! COPY VALUES FOR INPUT VARIABLES BEFORE CALL TO ROSR12
4670 ! ----------------------------------------------------------------------
4671 DO K = 1,NSOIL
4672 RHSTTin (K) = RHSTT (K)
4673 END DO
4674 DO K = 1,NSOIL
4675 CIin (K) = CI (K)
4676 END DO
4677 ! ----------------------------------------------------------------------
4678 ! CALL ROSR12 TO SOLVE THE TRI-DIAGONAL MATRIX
4679 ! ----------------------------------------------------------------------
4680 CALL ROSR12 (CI,AI,BI,CIin,RHSTTin,RHSTT,NSOIL)
4681 ! ----------------------------------------------------------------------
4682 ! SUM THE PREVIOUS SMC VALUE AND THE MATRIX SOLUTION TO GET A
4683 ! NEW VALUE. MIN ALLOWABLE VALUE OF SMC WILL BE 0.02.
4684 ! RUNOFF3: RUNOFF WITHIN SOIL LAYERS
4685 ! ----------------------------------------------------------------------
4686 WPLUS = 0.0
4687 RUNOFF3 = 0.
4688
4689 DDZ = - ZSOIL (1)
4690 DO K = 1,NSOIL
4691 IF (K /= 1) DDZ = ZSOIL (K - 1) - ZSOIL (K)
4692 SH2OOUT (K) = SH2OIN (K) + CI (K) + WPLUS / DDZ
4693 STOT = SH2OOUT (K) + SICE (K)
4694 IF (STOT > SMCMAX) THEN
4695 IF (K .eq. 1) THEN
4696 DDZ = - ZSOIL (1)
4697 ELSE
4698 KK11 = K - 1
4699 DDZ = - ZSOIL (K) + ZSOIL (KK11)
4700 END IF
4701 WPLUS = (STOT - SMCMAX) * DDZ
4702 ELSE
4703 WPLUS = 0.
4704 END IF
4705 SMC (K) = MAX ( MIN (STOT,SMCMAX),0.02 )
4706 SH2OOUT (K) = MAX ( (SMC (K) - SICE (K)),0.0)
4707 END DO
4708
4709 ! ----------------------------------------------------------------------
4710 ! UPDATE CANOPY WATER CONTENT/INTERCEPTION (CMC). CONVERT RHSCT TO
4711 ! AN 'AMOUNT' VALUE AND ADD TO PREVIOUS CMC VALUE TO GET NEW CMC.
4712 ! ----------------------------------------------------------------------
4713 RUNOFF3 = WPLUS
4714 CMC = CMC + DT * RHSCT
4715 IF (CMC < 1.E-20) CMC = 0.0
4716 CMC = MIN (CMC,CMCMAX)
4717
4718 ! ----------------------------------------------------------------------
4719 END SUBROUTINE SSTEP
4720 ! ----------------------------------------------------------------------
4721
4722 SUBROUTINE TBND (TU,TB,ZSOIL,ZBOT,K,NSOIL,TBND1)
4723
4724 ! ----------------------------------------------------------------------
4725 ! SUBROUTINE TBND
4726 ! ----------------------------------------------------------------------
4727 ! CALCULATE TEMPERATURE ON THE BOUNDARY OF THE LAYER BY INTERPOLATION OF
4728 ! THE MIDDLE LAYER TEMPERATURES
4729 ! ----------------------------------------------------------------------
4730 IMPLICIT NONE
4731 INTEGER, INTENT(IN) :: NSOIL
4732 INTEGER :: K
4733 REAL, INTENT(IN) :: TB, TU, ZBOT
4734 REAL, INTENT(OUT) :: TBND1
4735 REAL, DIMENSION(1:NSOIL), INTENT(IN) :: ZSOIL
4736 REAL :: ZB, ZUP
4737 REAL, PARAMETER :: T0 = 273.15
4738
4739 ! ----------------------------------------------------------------------
4740 ! USE SURFACE TEMPERATURE ON THE TOP OF THE FIRST LAYER
4741 ! ----------------------------------------------------------------------
4742 IF (K == 1) THEN
4743 ZUP = 0.
4744 ELSE
4745 ZUP = ZSOIL (K -1)
4746 END IF
4747 ! ----------------------------------------------------------------------
4748 ! USE DEPTH OF THE CONSTANT BOTTOM TEMPERATURE WHEN INTERPOLATE
4749 ! TEMPERATURE INTO THE LAST LAYER BOUNDARY
4750 ! ----------------------------------------------------------------------
4751 IF (K == NSOIL) THEN
4752 ZB = 2.* ZBOT - ZSOIL (K)
4753 ELSE
4754 ZB = ZSOIL (K +1)
4755 END IF
4756 ! ----------------------------------------------------------------------
4757 ! LINEAR INTERPOLATION BETWEEN THE AVERAGE LAYER TEMPERATURES
4758 ! ----------------------------------------------------------------------
4759
4760 TBND1 = TU + (TB - TU)* (ZUP - ZSOIL (K))/ (ZUP - ZB)
4761 ! ----------------------------------------------------------------------
4762 END SUBROUTINE TBND
4763 ! ----------------------------------------------------------------------
4764
4765
4766 SUBROUTINE TDFCND ( DF, SMC, QZ, SMCMAX, SH2O)
4767
4768 ! ----------------------------------------------------------------------
4769 ! SUBROUTINE TDFCND
4770 ! ----------------------------------------------------------------------
4771 ! CALCULATE THERMAL DIFFUSIVITY AND CONDUCTIVITY OF THE SOIL FOR A GIVEN
4772 ! POINT AND TIME.
4773 ! ----------------------------------------------------------------------
4774 ! PETERS-LIDARD APPROACH (PETERS-LIDARD et al., 1998)
4775 ! June 2001 CHANGES: FROZEN SOIL CONDITION.
4776 ! ----------------------------------------------------------------------
4777 IMPLICIT NONE
4778 REAL, INTENT(IN) :: QZ, SMC, SMCMAX, SH2O
4779 REAL, INTENT(OUT) :: DF
4780 REAL :: AKE, GAMMD, THKDRY, THKICE, THKO, &
4781 THKQTZ,THKSAT,THKS,THKW,SATRATIO,XU, &
4782 XUNFROZ
4783
4784 ! ----------------------------------------------------------------------
4785 ! WE NOW GET QUARTZ AS AN INPUT ARGUMENT (SET IN ROUTINE REDPRM):
4786 ! DATA QUARTZ /0.82, 0.10, 0.25, 0.60, 0.52,
4787 ! & 0.35, 0.60, 0.40, 0.82/
4788 ! ----------------------------------------------------------------------
4789 ! IF THE SOIL HAS ANY MOISTURE CONTENT COMPUTE A PARTIAL SUM/PRODUCT
4790 ! OTHERWISE USE A CONSTANT VALUE WHICH WORKS WELL WITH MOST SOILS
4791 ! ----------------------------------------------------------------------
4792 ! THKW ......WATER THERMAL CONDUCTIVITY
4793 ! THKQTZ ....THERMAL CONDUCTIVITY FOR QUARTZ
4794 ! THKO ......THERMAL CONDUCTIVITY FOR OTHER SOIL COMPONENTS
4795 ! THKS ......THERMAL CONDUCTIVITY FOR THE SOLIDS COMBINED(QUARTZ+OTHER)
4796 ! THKICE ....ICE THERMAL CONDUCTIVITY
4797 ! SMCMAX ....POROSITY (= SMCMAX)
4798 ! QZ .........QUARTZ CONTENT (SOIL TYPE DEPENDENT)
4799 ! ----------------------------------------------------------------------
4800 ! USE AS IN PETERS-LIDARD, 1998 (MODIF. FROM JOHANSEN, 1975).
4801
4802 ! PABLO GRUNMANN, 08/17/98
4803 ! REFS.:
4804 ! FAROUKI, O.T.,1986: THERMAL PROPERTIES OF SOILS. SERIES ON ROCK
4805 ! AND SOIL MECHANICS, VOL. 11, TRANS TECH, 136 PP.
4806 ! JOHANSEN, O., 1975: THERMAL CONDUCTIVITY OF SOILS. PH.D. THESIS,
4807 ! UNIVERSITY OF TRONDHEIM,
4808 ! PETERS-LIDARD, C. D., ET AL., 1998: THE EFFECT OF SOIL THERMAL
4809 ! CONDUCTIVITY PARAMETERIZATION ON SURFACE ENERGY FLUXES
4810 ! AND TEMPERATURES. JOURNAL OF THE ATMOSPHERIC SCIENCES,
4811 ! VOL. 55, PP. 1209-1224.
4812 ! ----------------------------------------------------------------------
4813 ! NEEDS PARAMETERS
4814 ! POROSITY(SOIL TYPE):
4815 ! POROS = SMCMAX
4816 ! SATURATION RATIO:
4817 ! PARAMETERS W/(M.K)
4818 SATRATIO = SMC / SMCMAX
4819 ! ICE CONDUCTIVITY:
4820 THKICE = 2.2
4821 ! WATER CONDUCTIVITY:
4822 THKW = 0.57
4823 ! THERMAL CONDUCTIVITY OF "OTHER" SOIL COMPONENTS
4824 ! IF (QZ .LE. 0.2) THKO = 3.0
4825 THKO = 2.0
4826 ! QUARTZ' CONDUCTIVITY
4827 THKQTZ = 7.7
4828 ! SOLIDS' CONDUCTIVITY
4829 THKS = (THKQTZ ** QZ)* (THKO ** (1. - QZ))
4830
4831 ! UNFROZEN FRACTION (FROM 1., i.e., 100%LIQUID, TO 0. (100% FROZEN))
4832 XUNFROZ = SH2O / SMC
4833 ! UNFROZEN VOLUME FOR SATURATION (POROSITY*XUNFROZ)
4834 XU = XUNFROZ * SMCMAX
4835
4836 ! SATURATED THERMAL CONDUCTIVITY
4837 THKSAT = THKS ** (1. - SMCMAX)* THKICE ** (SMCMAX - XU)* THKW ** &
4838 (XU)
4839
4840 ! DRY DENSITY IN KG/M3
4841 GAMMD = (1. - SMCMAX)*2700.
4842
4843 ! DRY THERMAL CONDUCTIVITY IN W.M-1.K-1
4844 THKDRY = (0.135* GAMMD+ 64.7)/ (2700. - 0.947* GAMMD)
4845 ! FROZEN
4846 IF ( (SH2O + 0.0005) < SMC ) THEN
4847 AKE = SATRATIO
4848 ! UNFROZEN
4849 ! RANGE OF VALIDITY FOR THE KERSTEN NUMBER (AKE)
4850 ELSE
4851
4852 ! KERSTEN NUMBER (USING "FINE" FORMULA, VALID FOR SOILS CONTAINING AT
4853 ! LEAST 5% OF PARTICLES WITH DIAMETER LESS THAN 2.E-6 METERS.)
4854 ! (FOR "COARSE" FORMULA, SEE PETERS-LIDARD ET AL., 1998).
4855
4856 IF ( SATRATIO > 0.1 ) THEN
4857
4858 AKE = LOG10 (SATRATIO) + 1.0
4859
4860 ! USE K = KDRY
4861 ELSE
4862
4863 AKE = 0.0
4864 END IF
4865 ! THERMAL CONDUCTIVITY
4866
4867 END IF
4868
4869 DF = AKE * (THKSAT - THKDRY) + THKDRY
4870 ! ----------------------------------------------------------------------
4871 END SUBROUTINE TDFCND
4872 ! ----------------------------------------------------------------------
4873
4874 SUBROUTINE TMPAVG (TAVG,TUP,TM,TDN,ZSOIL,NSOIL,K)
4875
4876 ! ----------------------------------------------------------------------
4877 ! SUBROUTINE TMPAVG
4878 ! ----------------------------------------------------------------------
4879 ! CALCULATE SOIL LAYER AVERAGE TEMPERATURE (TAVG) IN FREEZING/THAWING
4880 ! LAYER USING UP, DOWN, AND MIDDLE LAYER TEMPERATURES (TUP, TDN, TM),
4881 ! WHERE TUP IS AT TOP BOUNDARY OF LAYER, TDN IS AT BOTTOM BOUNDARY OF
4882 ! LAYER. TM IS LAYER PROGNOSTIC STATE TEMPERATURE.
4883 ! ----------------------------------------------------------------------
4884 IMPLICIT NONE
4885 INTEGER K
4886
4887 INTEGER NSOIL
4888 REAL DZ
4889 REAL DZH
4890 REAL T0
4891 REAL TAVG
4892 REAL TDN
4893 REAL TM
4894 REAL TUP
4895 REAL X0
4896 REAL XDN
4897 REAL XUP
4898
4899 REAL ZSOIL (NSOIL)
4900
4901 ! ----------------------------------------------------------------------
4902 PARAMETER (T0 = 2.7315E2)
4903 IF (K .eq. 1) THEN
4904 DZ = - ZSOIL (1)
4905 ELSE
4906 DZ = ZSOIL (K -1) - ZSOIL (K)
4907 END IF
4908
4909 DZH = DZ *0.5
4910 IF (TUP .lt. T0) THEN
4911 IF (TM .lt. T0) THEN
4912 ! ----------------------------------------------------------------------
4913 ! TUP, TM, TDN < T0
4914 ! ----------------------------------------------------------------------
4915 IF (TDN .lt. T0) THEN
4916 TAVG = (TUP + 2.0* TM + TDN)/ 4.0
4917 ! ----------------------------------------------------------------------
4918 ! TUP & TM < T0, TDN .ge. T0
4919 ! ----------------------------------------------------------------------
4920 ELSE
4921 X0 = (T0- TM) * DZH / (TDN - TM)
4922 TAVG = 0.5 * (TUP * DZH + TM * (DZH + X0) + T0* ( &
4923 & 2.* DZH - X0)) / DZ
4924 END IF
4925 ELSE
4926 ! ----------------------------------------------------------------------
4927 ! TUP < T0, TM .ge. T0, TDN < T0
4928 ! ----------------------------------------------------------------------
4929 IF (TDN .lt. T0) THEN
4930 XUP = (T0- TUP) * DZH / (TM - TUP)
4931 XDN = DZH - (T0- TM) * DZH / (TDN - TM)
4932 TAVG = 0.5 * (TUP * XUP + T0* (2.* DZ - XUP - XDN) &
4933 & + TDN * XDN) / DZ
4934 ! ----------------------------------------------------------------------
4935 ! TUP < T0, TM .ge. T0, TDN .ge. T0
4936 ! ----------------------------------------------------------------------
4937 ELSE
4938 XUP = (T0- TUP) * DZH / (TM - TUP)
4939 TAVG = 0.5 * (TUP * XUP + T0* (2.* DZ - XUP)) / DZ
4940 END IF
4941 END IF
4942 ELSE
4943 IF (TM .lt. T0) THEN
4944 ! ----------------------------------------------------------------------
4945 ! TUP .ge. T0, TM < T0, TDN < T0
4946 ! ----------------------------------------------------------------------
4947 IF (TDN .lt. T0) THEN
4948 XUP = DZH - (T0- TUP) * DZH / (TM - TUP)
4949 TAVG = 0.5 * (T0* (DZ - XUP) + TM * (DZH + XUP) &
4950 & + TDN * DZH) / DZ
4951 ! ----------------------------------------------------------------------
4952 ! TUP .ge. T0, TM < T0, TDN .ge. T0
4953 ! ----------------------------------------------------------------------
4954 ELSE
4955 XUP = DZH - (T0- TUP) * DZH / (TM - TUP)
4956 XDN = (T0- TM) * DZH / (TDN - TM)
4957 TAVG = 0.5 * (T0* (2.* DZ - XUP - XDN) + TM * &
4958 & (XUP + XDN)) / DZ
4959 END IF
4960 ELSE
4961 ! ----------------------------------------------------------------------
4962 ! TUP .ge. T0, TM .ge. T0, TDN < T0
4963 ! ----------------------------------------------------------------------
4964 IF (TDN .lt. T0) THEN
4965 XDN = DZH - (T0- TM) * DZH / (TDN - TM)
4966 TAVG = (T0* (DZ - XDN) +0.5* (T0+ TDN)* XDN) / DZ
4967 ! ----------------------------------------------------------------------
4968 ! TUP .ge. T0, TM .ge. T0, TDN .ge. T0
4969 ! ----------------------------------------------------------------------
4970 ELSE
4971 TAVG = (TUP + 2.0* TM + TDN) / 4.0
4972 END IF
4973 END IF
4974 END IF
4975 ! ----------------------------------------------------------------------
4976 END SUBROUTINE TMPAVG
4977 ! ----------------------------------------------------------------------
4978
4979 SUBROUTINE TRANSP (ET,NSOIL,ETP1,SMC,CMC,ZSOIL,SHDFAC,SMCWLT, &
4980 & CMCMAX,PC,CFACTR,SMCREF,SFCTMP,Q2,NROOT, &
4981 & RTDIS)
4982
4983 ! ----------------------------------------------------------------------
4984 ! SUBROUTINE TRANSP
4985 ! ----------------------------------------------------------------------
4986 ! CALCULATE TRANSPIRATION FOR THE VEG CLASS.
4987 ! ----------------------------------------------------------------------
4988 IMPLICIT NONE
4989 INTEGER I
4990 INTEGER K
4991 INTEGER NSOIL
4992
4993 INTEGER NROOT
4994 REAL CFACTR
4995 REAL CMC
4996 REAL CMCMAX
4997 REAL DENOM
4998 REAL ET (NSOIL)
4999 REAL ETP1
5000 REAL ETP1A
5001 !.....REAL PART(NSOIL)
5002 REAL GX (7)
5003 REAL PC
5004 REAL Q2
5005 REAL RTDIS (NSOIL)
5006 REAL RTX
5007 REAL SFCTMP
5008 REAL SGX
5009 REAL SHDFAC
5010 REAL SMC (NSOIL)
5011 REAL SMCREF
5012 REAL SMCWLT
5013
5014 ! ----------------------------------------------------------------------
5015 ! INITIALIZE PLANT TRANSP TO ZERO FOR ALL SOIL LAYERS.
5016 ! ----------------------------------------------------------------------
5017 REAL ZSOIL (NSOIL)
5018 DO K = 1,NSOIL
5019 ET (K) = 0.
5020 ! ----------------------------------------------------------------------
5021 ! CALCULATE AN 'ADJUSTED' POTENTIAL TRANSPIRATION
5022 ! IF STATEMENT BELOW TO AVOID TANGENT LINEAR PROBLEMS NEAR ZERO
5023 ! NOTE: GX AND OTHER TERMS BELOW REDISTRIBUTE TRANSPIRATION BY LAYER,
5024 ! ET(K), AS A FUNCTION OF SOIL MOISTURE AVAILABILITY, WHILE PRESERVING
5025 ! TOTAL ETP1A.
5026 ! ----------------------------------------------------------------------
5027 END DO
5028 IF (CMC .ne. 0.0) THEN
5029 ETP1A = SHDFAC * PC * ETP1 * (1.0- (CMC / CMCMAX) ** CFACTR)
5030 ELSE
5031 ETP1A = SHDFAC * PC * ETP1
5032 END IF
5033 SGX = 0.0
5034 DO I = 1,NROOT
5035 GX (I) = ( SMC (I) - SMCWLT ) / ( SMCREF - SMCWLT )
5036 GX (I) = MAX ( MIN ( GX (I), 1. ), 0. )
5037 SGX = SGX + GX (I)
5038 END DO
5039
5040 SGX = SGX / NROOT
5041 DENOM = 0.
5042 DO I = 1,NROOT
5043 RTX = RTDIS (I) + GX (I) - SGX
5044 GX (I) = GX (I) * MAX ( RTX, 0. )
5045 DENOM = DENOM + GX (I)
5046 END DO
5047
5048 IF (DENOM .le. 0.0) DENOM = 1.
5049 DO I = 1,NROOT
5050 ET (I) = ETP1A * GX (I) / DENOM
5051 ! ----------------------------------------------------------------------
5052 ! ABOVE CODE ASSUMES A VERTICALLY UNIFORM ROOT DISTRIBUTION
5053 ! CODE BELOW TESTS A VARIABLE ROOT DISTRIBUTION
5054 ! ----------------------------------------------------------------------
5055 ! ET(1) = ( ZSOIL(1) / ZSOIL(NROOT) ) * GX * ETP1A
5056 ! ET(1) = ( ZSOIL(1) / ZSOIL(NROOT) ) * ETP1A
5057 ! ----------------------------------------------------------------------
5058 ! USING ROOT DISTRIBUTION AS WEIGHTING FACTOR
5059 ! ----------------------------------------------------------------------
5060 ! ET(1) = RTDIS(1) * ETP1A
5061 ! ET(1) = ETP1A * PART(1)
5062 ! ----------------------------------------------------------------------
5063 ! LOOP DOWN THRU THE SOIL LAYERS REPEATING THE OPERATION ABOVE,
5064 ! BUT USING THE THICKNESS OF THE SOIL LAYER (RATHER THAN THE
5065 ! ABSOLUTE DEPTH OF EACH LAYER) IN THE FINAL CALCULATION.
5066 ! ----------------------------------------------------------------------
5067 ! DO K = 2,NROOT
5068 ! GX = ( SMC(K) - SMCWLT ) / ( SMCREF - SMCWLT )
5069 ! GX = MAX ( MIN ( GX, 1. ), 0. )
5070 ! TEST CANOPY RESISTANCE
5071 ! GX = 1.0
5072 ! ET(K) = ((ZSOIL(K)-ZSOIL(K-1))/ZSOIL(NROOT))*GX*ETP1A
5073 ! ET(K) = ((ZSOIL(K)-ZSOIL(K-1))/ZSOIL(NROOT))*ETP1A
5074 ! ----------------------------------------------------------------------
5075 ! USING ROOT DISTRIBUTION AS WEIGHTING FACTOR
5076 ! ----------------------------------------------------------------------
5077 ! ET(K) = RTDIS(K) * ETP1A
5078 ! ET(K) = ETP1A*PART(K)
5079 ! END DO
5080 END DO
5081 ! ----------------------------------------------------------------------
5082 END SUBROUTINE TRANSP
5083 ! ----------------------------------------------------------------------
5084
5085 SUBROUTINE WDFCND (WDF,WCND,SMC,SMCMAX,BEXP,DKSAT,DWSAT, &
5086 & SICEMAX)
5087
5088 ! ----------------------------------------------------------------------
5089 ! SUBROUTINE WDFCND
5090 ! ----------------------------------------------------------------------
5091 ! CALCULATE SOIL WATER DIFFUSIVITY AND SOIL HYDRAULIC CONDUCTIVITY.
5092 ! ----------------------------------------------------------------------
5093 IMPLICIT NONE
5094 REAL BEXP
5095 REAL DKSAT
5096 REAL DWSAT
5097 REAL EXPON
5098 REAL FACTR1
5099 REAL FACTR2
5100 REAL SICEMAX
5101 REAL SMC
5102 REAL SMCMAX
5103 REAL VKwgt
5104 REAL WCND
5105
5106 ! ----------------------------------------------------------------------
5107 ! CALC THE RATIO OF THE ACTUAL TO THE MAX PSBL SOIL H2O CONTENT
5108 ! ----------------------------------------------------------------------
5109 REAL WDF
5110 FACTR1 = 0.2 / SMCMAX
5111
5112 ! ----------------------------------------------------------------------
5113 ! PREP AN EXPNTL COEF AND CALC THE SOIL WATER DIFFUSIVITY
5114 ! ----------------------------------------------------------------------
5115 FACTR2 = SMC / SMCMAX
5116 EXPON = BEXP + 2.0
5117
5118 ! ----------------------------------------------------------------------
5119 ! FROZEN SOIL HYDRAULIC DIFFUSIVITY. VERY SENSITIVE TO THE VERTICAL
5120 ! GRADIENT OF UNFROZEN WATER. THE LATTER GRADIENT CAN BECOME VERY
5121 ! EXTREME IN FREEZING/THAWING SITUATIONS, AND GIVEN THE RELATIVELY
5122 ! FEW AND THICK SOIL LAYERS, THIS GRADIENT SUFFERES SERIOUS
5123 ! TRUNCTION ERRORS YIELDING ERRONEOUSLY HIGH VERTICAL TRANSPORTS OF
5124 ! UNFROZEN WATER IN BOTH DIRECTIONS FROM HUGE HYDRAULIC DIFFUSIVITY.
5125 ! THEREFORE, WE FOUND WE HAD TO ARBITRARILY CONSTRAIN WDF
5126 ! --
5127 ! VERSION D_10CM: ........ FACTR1 = 0.2/SMCMAX
5128 ! WEIGHTED APPROACH...................... PABLO GRUNMANN, 28_SEP_1999.
5129 ! ----------------------------------------------------------------------
5130 WDF = DWSAT * FACTR2 ** EXPON
5131 IF (SICEMAX .gt. 0.0) THEN
5132 VKWGT = 1./ (1. + (500.* SICEMAX)**3.)
5133 WDF = VKWGT * WDF + (1. - VKWGT)* DWSAT * FACTR1** EXPON
5134 ! ----------------------------------------------------------------------
5135 ! RESET THE EXPNTL COEF AND CALC THE HYDRAULIC CONDUCTIVITY
5136 ! ----------------------------------------------------------------------
5137 END IF
5138 EXPON = (2.0 * BEXP) + 3.0
5139 WCND = DKSAT * FACTR2 ** EXPON
5140
5141 ! ----------------------------------------------------------------------
5142 END SUBROUTINE WDFCND
5143 ! ----------------------------------------------------------------------
5144
5145 SUBROUTINE SFCDIF_off (ZLM,Z0,THZ0,THLM,SFCSPD,CZIL,AKMS,AKHS)
5146
5147 ! ----------------------------------------------------------------------
5148 ! SUBROUTINE SFCDIF (renamed SFCDIF_off to avoid clash with Eta PBL)
5149 ! ----------------------------------------------------------------------
5150 ! CALCULATE SURFACE LAYER EXCHANGE COEFFICIENTS VIA ITERATIVE PROCESS.
5151 ! SEE CHEN ET AL (1997, BLM)
5152 ! ----------------------------------------------------------------------
5153
5154 IMPLICIT NONE
5155 REAL WWST, WWST2, G, VKRM, EXCM, BETA, BTG, ELFC, WOLD, WNEW
5156 REAL PIHF, EPSU2, EPSUST, EPSIT, EPSA, ZTMIN, ZTMAX, HPBL, &
5157 & SQVISC
5158 REAL RIC, RRIC, FHNEU, RFC, RFAC, ZZ, PSLMU, PSLMS, PSLHU, &
5159 & PSLHS
5160 REAL XX, PSPMU, YY, PSPMS, PSPHU, PSPHS, ZLM, Z0, THZ0, THLM
5161 REAL SFCSPD, CZIL, AKMS, AKHS, ZILFC, ZU, ZT, RDZ, CXCH
5162 REAL DTHV, DU2, BTGH, WSTAR2, USTAR, ZSLU, ZSLT, RLOGU, RLOGT
5163 REAL RLMO, ZETALT, ZETALU, ZETAU, ZETAT, XLU4, XLT4, XU4, XT4
5164 !CC ......REAL ZTFC
5165
5166 REAL XLU, XLT, XU, XT, PSMZ, SIMM, PSHZ, SIMH, USTARK, RLMN, &
5167 & RLMA
5168
5169 INTEGER ITRMX, ILECH, ITR
5170 PARAMETER &
5171 & (WWST = 1.2,WWST2 = WWST * WWST,G = 9.8,VKRM = 0.40, &
5172 & EXCM = 0.001 &
5173 & ,BETA = 1./270.,BTG = BETA * G,ELFC = VKRM * BTG &
5174 & ,WOLD =.15,WNEW = 1. - WOLD,ITRMX = 05, &
5175 & PIHF = 3.14159265/2.)
5176 PARAMETER &
5177 & (EPSU2 = 1.E-4,EPSUST = 0.07,EPSIT = 1.E-4,EPSA = 1.E-8 &
5178 & ,ZTMIN = -5.,ZTMAX = 1.,HPBL = 1000.0 &
5179 & ,SQVISC = 258.2)
5180 PARAMETER &
5181 & (RIC = 0.183,RRIC = 1.0/ RIC,FHNEU = 0.8,RFC = 0.191 &
5182 & ,RFAC = RIC / (FHNEU * RFC * RFC))
5183
5184 ! ----------------------------------------------------------------------
5185 ! NOTE: THE TWO CODE BLOCKS BELOW DEFINE FUNCTIONS
5186 ! ----------------------------------------------------------------------
5187 ! LECH'S SURFACE FUNCTIONS
5188 ! ----------------------------------------------------------------------
5189 PSLMU (ZZ)= -0.96* log (1.0-4.5* ZZ)
5190 PSLMS (ZZ)= ZZ * RRIC -2.076* (1. -1./ (ZZ +1.))
5191 PSLHU (ZZ)= -0.96* log (1.0-4.5* ZZ)
5192
5193 ! ----------------------------------------------------------------------
5194 ! PAULSON'S SURFACE FUNCTIONS
5195 ! ----------------------------------------------------------------------
5196 PSLHS (ZZ)= ZZ * RFAC -2.076* (1. -1./ (ZZ +1.))
5197 PSPMU (XX)= -2.* log ( (XX +1.)*0.5) - log ( (XX * XX +1.)*0.5) &
5198 & +2.* ATAN (XX) &
5199 &- PIHF
5200 PSPMS (YY)= 5.* YY
5201 PSPHU (XX)= -2.* log ( (XX * XX +1.)*0.5)
5202
5203 ! ----------------------------------------------------------------------
5204 ! THIS ROUTINE SFCDIF CAN HANDLE BOTH OVER OPEN WATER (SEA, OCEAN) AND
5205 ! OVER SOLID SURFACE (LAND, SEA-ICE).
5206 ! ----------------------------------------------------------------------
5207 PSPHS (YY)= 5.* YY
5208
5209 ! ----------------------------------------------------------------------
5210 ! ZTFC: RATIO OF ZOH/ZOM LESS OR EQUAL THAN 1
5211 ! C......ZTFC=0.1
5212 ! CZIL: CONSTANT C IN Zilitinkevich, S. S.1995,:NOTE ABOUT ZT
5213 ! ----------------------------------------------------------------------
5214 ILECH = 0
5215
5216 ! ----------------------------------------------------------------------
5217 ZILFC = - CZIL * VKRM * SQVISC
5218 ! C.......ZT=Z0*ZTFC
5219 ZU = Z0
5220 RDZ = 1./ ZLM
5221 CXCH = EXCM * RDZ
5222 DTHV = THLM - THZ0
5223
5224 ! ----------------------------------------------------------------------
5225 ! BELJARS CORRECTION OF USTAR
5226 ! ----------------------------------------------------------------------
5227 DU2 = MAX (SFCSPD * SFCSPD,EPSU2)
5228 !cc If statements to avoid TANGENT LINEAR problems near zero
5229 BTGH = BTG * HPBL
5230 IF (BTGH * AKHS * DTHV .ne. 0.0) THEN
5231 WSTAR2 = WWST2* ABS (BTGH * AKHS * DTHV)** (2./3.)
5232 ELSE
5233 WSTAR2 = 0.0
5234 END IF
5235
5236 ! ----------------------------------------------------------------------
5237 ! ZILITINKEVITCH APPROACH FOR ZT
5238 ! ----------------------------------------------------------------------
5239 USTAR = MAX (SQRT (AKMS * SQRT (DU2+ WSTAR2)),EPSUST)
5240
5241 ! ----------------------------------------------------------------------
5242 ZT = EXP (ZILFC * SQRT (USTAR * Z0))* Z0
5243 ZSLU = ZLM + ZU
5244 ! PRINT*,'ZSLT=',ZSLT
5245 ! PRINT*,'ZLM=',ZLM
5246 ! PRINT*,'ZT=',ZT
5247
5248 ZSLT = ZLM + ZT
5249 RLOGU = log (ZSLU / ZU)
5250
5251 RLOGT = log (ZSLT / ZT)
5252 ! PRINT*,'RLMO=',RLMO
5253 ! PRINT*,'ELFC=',ELFC
5254 ! PRINT*,'AKHS=',AKHS
5255 ! PRINT*,'DTHV=',DTHV
5256 ! PRINT*,'USTAR=',USTAR
5257
5258 RLMO = ELFC * AKHS * DTHV / USTAR **3
5259 ! ----------------------------------------------------------------------
5260 ! 1./MONIN-OBUKKHOV LENGTH-SCALE
5261 ! ----------------------------------------------------------------------
5262 DO ITR = 1,ITRMX
5263 ZETALT = MAX (ZSLT * RLMO,ZTMIN)
5264 RLMO = ZETALT / ZSLT
5265 ZETALU = ZSLU * RLMO
5266 ZETAU = ZU * RLMO
5267
5268 ZETAT = ZT * RLMO
5269 IF (ILECH .eq. 0) THEN
5270 IF (RLMO .lt. 0.)THEN
5271 XLU4 = 1. -16.* ZETALU
5272 XLT4 = 1. -16.* ZETALT
5273 XU4 = 1. -16.* ZETAU
5274
5275 XT4 = 1. -16.* ZETAT
5276 XLU = SQRT (SQRT (XLU4))
5277 XLT = SQRT (SQRT (XLT4))
5278 XU = SQRT (SQRT (XU4))
5279
5280 XT = SQRT (SQRT (XT4))
5281 ! PRINT*,'-----------1------------'
5282 ! PRINT*,'PSMZ=',PSMZ
5283 ! PRINT*,'PSPMU(ZETAU)=',PSPMU(ZETAU)
5284 ! PRINT*,'XU=',XU
5285 ! PRINT*,'------------------------'
5286 PSMZ = PSPMU (XU)
5287 SIMM = PSPMU (XLU) - PSMZ + RLOGU
5288 PSHZ = PSPHU (XT)
5289 SIMH = PSPHU (XLT) - PSHZ + RLOGT
5290 ELSE
5291 ZETALU = MIN (ZETALU,ZTMAX)
5292 ZETALT = MIN (ZETALT,ZTMAX)
5293 ! PRINT*,'-----------2------------'
5294 ! PRINT*,'PSMZ=',PSMZ
5295 ! PRINT*,'PSPMS(ZETAU)=',PSPMS(ZETAU)
5296 ! PRINT*,'ZETAU=',ZETAU
5297 ! PRINT*,'------------------------'
5298 PSMZ = PSPMS (ZETAU)
5299 SIMM = PSPMS (ZETALU) - PSMZ + RLOGU
5300 PSHZ = PSPHS (ZETAT)
5301 SIMH = PSPHS (ZETALT) - PSHZ + RLOGT
5302 END IF
5303 ! ----------------------------------------------------------------------
5304 ! LECH'S FUNCTIONS
5305 ! ----------------------------------------------------------------------
5306 ELSE
5307 IF (RLMO .lt. 0.)THEN
5308 ! PRINT*,'-----------3------------'
5309 ! PRINT*,'PSMZ=',PSMZ
5310 ! PRINT*,'PSLMU(ZETAU)=',PSLMU(ZETAU)
5311 ! PRINT*,'ZETAU=',ZETAU
5312 ! PRINT*,'------------------------'
5313 PSMZ = PSLMU (ZETAU)
5314 SIMM = PSLMU (ZETALU) - PSMZ + RLOGU
5315 PSHZ = PSLHU (ZETAT)
5316 SIMH = PSLHU (ZETALT) - PSHZ + RLOGT
5317 ELSE
5318 ZETALU = MIN (ZETALU,ZTMAX)
5319
5320 ZETALT = MIN (ZETALT,ZTMAX)
5321 ! PRINT*,'-----------4------------'
5322 ! PRINT*,'PSMZ=',PSMZ
5323 ! PRINT*,'PSLMS(ZETAU)=',PSLMS(ZETAU)
5324 ! PRINT*,'ZETAU=',ZETAU
5325 ! PRINT*,'------------------------'
5326 PSMZ = PSLMS (ZETAU)
5327 SIMM = PSLMS (ZETALU) - PSMZ + RLOGU
5328 PSHZ = PSLHS (ZETAT)
5329 SIMH = PSLHS (ZETALT) - PSHZ + RLOGT
5330 END IF
5331 ! ----------------------------------------------------------------------
5332 ! BELJAARS CORRECTION FOR USTAR
5333 ! ----------------------------------------------------------------------
5334 END IF
5335
5336 ! ----------------------------------------------------------------------
5337 ! ZILITINKEVITCH FIX FOR ZT
5338 ! ----------------------------------------------------------------------
5339 USTAR = MAX (SQRT (AKMS * SQRT (DU2+ WSTAR2)),EPSUST)
5340
5341 ZT = EXP (ZILFC * SQRT (USTAR * Z0))* Z0
5342 ZSLT = ZLM + ZT
5343 !-----------------------------------------------------------------------
5344 RLOGT = log (ZSLT / ZT)
5345 USTARK = USTAR * VKRM
5346 AKMS = MAX (USTARK / SIMM,CXCH)
5347 !-----------------------------------------------------------------------
5348 ! IF STATEMENTS TO AVOID TANGENT LINEAR PROBLEMS NEAR ZERO
5349 !-----------------------------------------------------------------------
5350 AKHS = MAX (USTARK / SIMH,CXCH)
5351 IF (BTGH * AKHS * DTHV .ne. 0.0) THEN
5352 WSTAR2 = WWST2* ABS (BTGH * AKHS * DTHV)** (2./3.)
5353 ELSE
5354 WSTAR2 = 0.0
5355 END IF
5356 !-----------------------------------------------------------------------
5357 RLMN = ELFC * AKHS * DTHV / USTAR **3
5358 !-----------------------------------------------------------------------
5359 ! IF(ABS((RLMN-RLMO)/RLMA).LT.EPSIT) GO TO 110
5360 !-----------------------------------------------------------------------
5361 RLMA = RLMO * WOLD+ RLMN * WNEW
5362 !-----------------------------------------------------------------------
5363 RLMO = RLMA
5364 ! PRINT*,'----------------------------'
5365 ! PRINT*,'SFCDIF OUTPUT ! ! ! ! ! ! ! ! ! ! ! !'
5366
5367 ! PRINT*,'ZLM=',ZLM
5368 ! PRINT*,'Z0=',Z0
5369 ! PRINT*,'THZ0=',THZ0
5370 ! PRINT*,'THLM=',THLM
5371 ! PRINT*,'SFCSPD=',SFCSPD
5372 ! PRINT*,'CZIL=',CZIL
5373 ! PRINT*,'AKMS=',AKMS
5374 ! PRINT*,'AKHS=',AKHS
5375 ! PRINT*,'----------------------------'
5376
5377 END DO
5378 ! ----------------------------------------------------------------------
5379 END SUBROUTINE SFCDIF_off
5380 ! ----------------------------------------------------------------------
5381
5382 END MODULE module_sf_noahlsm