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,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, limit to 0 C)
606 SFCTSNO=MIN(SFCTMP,273.15)
607 E2SAT=611.2*EXP(6174.*(1./273.15 - 1./SFCTSNO))
608 Q2SAT=0.622*E2SAT/(SFCPRS-E2SAT)
609 DQSDT2=Q2SAT*6174./(SFCTSNO**2)
610 ! for snow cover fraction at 0 C, ground temp will not change, so DQSDT2 effectively zero
611 IF(T1 .GT. 273. .AND. SNOWC(I,J) .GT. 0.)DQSDT2=DQSDT2*(1.-SNOWC(I,J))
612 ENDIF
613
614 IF(ICE.EQ.0)THEN
615 TBOT=TMN(I,J)
616 ELSE
617 TBOT=271.16
618 ENDIF
619 IF(VEGTYP.EQ.25) SHDFAC=0.0000
620 IF(VEGTYP.EQ.26) SHDFAC=0.0000
621 IF(VEGTYP.EQ.27) SHDFAC=0.0000
622 IF(SOILTYP.EQ.14.AND.XICE(I,J).EQ.0.)THEN
623 IF(IPRINT)PRINT*,' SOIL TYPE FOUND TO BE WATER AT A LAND-POINT'
624 IF(IPRINT)PRINT*,i,j,'RESET SOIL in surfce.F'
625 SOILTYP=7
626 ENDIF
627 CMC=CANWAT(I,J)
628
629 !-------------------------------------------
630 !*** convert snow depth from mm to meter
631 !
632 ! IF(RDMAXALB) THEN
633 ! SNOALB=ALBMAX(I,J)*0.01
634 ! ELSE
635 ! SNOALB=MAXALB(IVGTPK)*0.01
636 ! ENDIF
637 ! IF(RDBRDALB) THEN
638 ! ALBBRD=ALBEDO(I,J)*0.01
639 ! ELSE
640 ! ALBBRD=ALBD(IVGTPK,ISN)*0.01
641 ! ENDIF
642
643 ! SNOALB1=0.80
644 ! SHMIN=0.00
645 ALBBRD=ALBBCK(I,J)
646 Z0BRD=Z0(I,J)
647 !FEI: temporaray arrays above need to be changed later by using SI
648
649 DO 70 NS=1,NSOIL
650 SMC(NS)=SMOIS(I,NS,J)
651 STC(NS)=TSLB(I,NS,J) !STEMP
652 SWC(NS)=SH2O(I,NS,J)
653 70 CONTINUE
654 !
655 if ( (SNEQV.ne.0..AND.SNOWHK.eq.0.).or.(SNOWHK.le.SNEQV) )THEN
656 SNOWHK= 5.*SNEQV
657 endif
658 !
659
660 !Fei: urban. for urban surface, if calling UCM, redefine urban as 5: Cropland/Grassland Mosaic
661
662 IF(UCMCALL == 1 ) THEN
663 IF( IVGTYP(I,J) == 1 .or. IVGTYP(I,J) == 31 .or. &
664 IVGTYP(I,J) == 32 .or. IVGTYP(I,J) == 33) THEN
665 VEGTYP = 5
666 SHDFAC = 0.8
667 ALBEDOK =0.2
668 ALBBRD =0.2
669 IF ( FRC_URB2D(I,J) < 0.99 ) THEN
670 T1= ( TSK(I,J) -FRC_URB2D(I,J) * TS_URB2D (I,J) )/ (1-FRC_URB2D(I,J))
671 ELSE
672 T1 = TSK(I,J)
673 ENDIF
674 ENDIF
675 ELSE
676 IF( IVGTYP(I,J) == 1 .or. IVGTYP(I,J) == 31 .or. &
677 IVGTYP(I,J) == 32 .or. IVGTYP(I,J) == 33) THEN
678 VEGTYP = 1
679 ENDIF
680 ENDIF
681
682 IF(IPRINT) THEN
683 !
684 print*, 'BEFORE SFLX, in Noahlsm_driver'
685 print*, 'ICE', ICE, 'DT',DT, 'ZLVL',ZLVL, 'NSOIL', NSOIL, &
686 'SLDPTH', SLDPTH, 'LOCAL',LOCAL, 'LUTYPE',&
687 LUTYPE, 'SLTYPE',SLTYPE, 'LWDN',LWDN, 'SOLDN',SOLDN, &
688 'SFCPRS',SFCPRS, 'PRCP',PRCP,'SFCTMP',SFCTMP,'Q2K',Q2K, &
689 'TH2',TH2,'Q2SAT',Q2SAT,'DQSDT2',DQSDT2,'VEGTYP', VEGTYP,&
690 'SOILTYP',SOILTYP, 'SLOPETYP',SLOPETYP, 'SHDFAC',SHDFAC,&
691 'SHMIN',SHMIN, 'ALBBRD',ALBBRD,'SNOALB1',SNOALB1,'TBOT',&
692 TBOT, 'Z0BRD',Z0BRD, 'Z0K',Z0K, 'CMC',CMC, 'T1',T1,'STC',&
693 STC, 'SMC',SMC, 'SWC',SWC,'SNOWHK',SNOWHK,'SNEQV',SNEQV,&
694 'ALBEDO',ALBEDO,'CHK',CHK,'ETA',ETA,'SHEAT',SHEAT, &
695 'ETA_KINEMATIC',ETA_KINEMATIC, 'FDOWN',FDOWN,'EC',EC, &
696 'EDIR',EDIR,'ET',ET,'ETT',ETT,'ESNOW',ESNOW,'DRIP',DRIP,&
697 'DEW',DEW,'BETA',BETA,'ETP',ETP,'SSOIL',SSOIL,'FLX1',FLX1,&
698 'FLX2',FLX2,'FLX3',FLX3,'SNOMLT',SNOMLT,'SNCOVR',SNCOVR,&
699 'RUNOFF1',RUNOFF1,'RUNOFF2',RUNOFF2,'RUNOFF3',RUNOFF3, &
700 'RC',RC, 'PC',PC,'RSMIN',RSMIN,'XLAI',XLAI,'RCS',RCS, &
701 'RCT',RCT,'RCQ',RCQ,'RCSOIL',RCSOIL,'SOILW',SOILW, &
702 'SOILM',SOILM,'Q1',Q1,'SMCWLT',SMCWLT,'SMCDRY',SMCDRY,&
703 'SMCREF',SMCREF,'SMCMAX',SMCMAX,'NROOT',NROOT
704 endif
705
706
707 CALL SFLX (FFROZP, ICE,DT,ZLVL,NSOIL,SLDPTH, & !C
708 LOCAL, & !L
709 LUTYPE, SLTYPE, & !CL
710 LWDN,SOLDN,SOLNET,SFCPRS,PRCP,SFCTMP,Q2K,DUMMY, & !F
711 DUMMY,DUMMY, DUMMY, & !F PRCPRAIN not used
712 TH2,Q2SAT,DQSDT2, & !I
713 VEGTYP,SOILTYP,SLOPETYP,SHDFAC,SHMIN,DUMMY, & !I
714 ALBBRD, SNOALB1,TBOT, Z0BRD, Z0K, EMISSI, & !S
715 CMC,T1,STC,SMC,SWC,SNOWHK,SNEQV,ALBEDOK,CHK,dummy,& !H
716 ETA,SHEAT, ETA_KINEMATIC,FDOWN, & !O
717 EC,EDIR,ET,ETT,ESNOW,DRIP,DEW, & !O
718 BETA,ETP,SSOIL, & !O
719 FLX1,FLX2,FLX3, & !O
720 SNOMLT,SNCOVR, & !O
721 RUNOFF1,RUNOFF2,RUNOFF3, & !O
722 RC,PC,RSMIN,XLAI,RCS,RCT,RCQ,RCSOIL, & !O
723 SOILW,SOILM,Q1, & !D
724 SMCWLT,SMCDRY,SMCREF,SMCMAX,NROOT)
725
726
727 IF(IPRINT) THEN
728
729 print*, 'AFTER SFLX, in Noahlsm_driver'
730 print*, 'ICE', ICE, 'DT',DT, 'ZLVL',ZLVL, 'NSOIL', NSOIL, &
731 'SLDPTH', SLDPTH, 'LOCAL',LOCAL, 'LUTYPE',&
732 LUTYPE, 'SLTYPE',SLTYPE, 'LWDN',LWDN, 'SOLDN',SOLDN, &
733 'SFCPRS',SFCPRS, 'PRCP',PRCP,'SFCTMP',SFCTMP,'Q2K',Q2K, &
734 'TH2',TH2,'Q2SAT',Q2SAT,'DQSDT2',DQSDT2,'VEGTYP', VEGTYP,&
735 'SOILTYP',SOILTYP, 'SLOPETYP',SLOPETYP, 'SHDFAC',SHDFAC,&
736 'SHDMIN',SHMIN, 'ALBBRD',ALBBRD,'SNOALB',SNOALB1,'TBOT',&
737 TBOT, 'Z0BRD',Z0BRD, 'Z0K',Z0K, 'CMC',CMC, 'T1',T1,'STC',&
738 STC, 'SMC',SMC, 'SWc',SWC,'SNOWHK',SNOWHK,'SNEQV',SNEQV,&
739 'ALBEDO',ALBEDOK,'CHK',CHK,'ETA',ETA,'SHEAT',SHEAT, &
740 'ETA_KINEMATIC',ETA_KINEMATIC, 'FDOWN',FDOWN,'EC',EC, &
741 'EDIR',EDIR,'ET',ET,'ETT',ETT,'ESNOW',ESNOW,'DRIP',DRIP,&
742 'DEW',DEW,'BETA',BETA,'ETP',ETP,'SSOIL',SSOIL,'FLX1',FLX1,&
743 'FLX2',FLX2,'FLX3',FLX3,'SNOMLT',SNOMLT,'SNCOVR',SNCOVR,&
744 'RUNOFF1',RUNOFF1,'RUNOFF2',RUNOFF2,'RUNOFF3',RUNOFF3, &
745 'RC',RC, 'PC',PC,'RSMIN',RSMIN,'XLAI',XLAI,'RCS',RCS, &
746 'RCT',RCT,'RCQ',RCQ,'RCSOIL',RCSOIL,'SOILW',SOILW, &
747 'SOILM',SOILM,'Q1',Q1,'SMCWLT',SMCWLT,'SMCDRY',SMCDRY,&
748 'SMCREF',SMCREF,'SMCMAX',SMCMAX,'NROOT',NROOT
749 endif
750
751 !*** UPDATE STATE VARIABLES
752 CANWAT(I,J)=CMC
753 SNOW(I,J)=SNEQV*1000.
754 ! SNOWH(I,J)=SNOWHK*1000.
755 SNOWH(I,J)=SNOWHK ! SNOWHK in meters
756 ALBEDO(I,J)=ALBEDOK
757 ! MEK Nov2006 turn off
758 ! ZNT(I,J)=Z0K
759 TSK(I,J)=T1
760 HFX(I,J)=SHEAT
761 ! MEk Jul07 add potential evap accum
762 POTEVP(I,J)=POTEVP(I,J)+ETP*FDTW
763 QFX(I,J)=ETA_KINEMATIC
764
765 LH(I,J)=ETA
766 GRDFLX(I,J)=SSOIL
767 SNOWC(I,J)=SNCOVR
768 CHS2(I,J)=CQS2(I,J)
769 ! prevent diagnostic ground q (q1) from being greater than qsat(tsk)
770 ! as happens over snow cover where the cqs2 value also becomes irrelevant
771 ! by setting cqs2=chs in this situation the 2m q should become just qv(k=1)
772 IF (Q1 .GT. QSFC(I,J)) THEN
773 CQS2(I,J) = CHS(I,J)
774 ENDIF
775 ! QSFC(I,J)=Q1
776 ! Convert QSFC back to mixing ratio
777 QSFC(I,J)= Q1/(1.0-Q1)
778 !
779 DO 80 NS=1,NSOIL
780 SMOIS(I,NS,J)=SMC(NS)
781 TSLB(I,NS,J)=STC(NS) ! STEMP
782 SH2O(I,NS,J)=SWC(NS)
783 80 CONTINUE
784 ! ENDIF
785
786 IF (UCMCALL == 1 ) THEN ! Beginning of UCM CALL if block
787 !--------------------------------------
788 ! URBAN CANOPY MODEL START - urban
789 !--------------------------------------
790 ! Input variables lsm --> urban
791
792
793 IF( IVGTYP(I,J) == 1 .or. IVGTYP(I,J) == 31 .or. &
794 IVGTYP(I,J) == 32 .or. IVGTYP(I,J) == 33 ) THEN
795
796 ! Call urban
797
798 !
799 UTYPE_URB = UTYPE_URB2D(I,J) !urban type (low, high or industrial)
800
801 TA_URB = SFCTMP ! [K]
802 QA_URB = Q2K ! [kg/kg]
803 UA_URB = SQRT(U_PHY(I,1,J)**2.+V_PHY(I,1,J)**2.)
804 U1_URB = U_PHY(I,1,J)
805 V1_URB = V_PHY(I,1,J)
806 IF(UA_URB < 1.) UA_URB=1. ! [m/s]
807 SSG_URB = SOLDN ! [W/m/m]
808 SSGD_URB = 0.8*SOLDN ! [W/m/m]
809 SSGQ_URB = SSG_URB-SSGD_URB ! [W/m/m]
810 LLG_URB = LWDN ! [W/m/m]
811 RAIN_URB = RAINBL(I,J) ! [mm]
812 RHOO_URB = SFCPRS / (287.04 * SFCTMP * (1.0+ 0.61 * Q2K)) ![kg/m/m/m]
813 ZA_URB = ZLVL ! [m]
814 DELT_URB = DT ! [sec]
815 XLAT_URB = XLAT_URB2D(I,J) ! [deg]
816 COSZ_URB = COSZ_URB2D(I,J) !
817 OMG_URB = OMG_URB2D(I,J) !
818 ZNT_URB = ZNT(I,J)
819
820 LSOLAR_URB = .FALSE.
821
822 TR_URB = TR_URB2D(I,J)
823 TB_URB = TB_URB2D(I,J)
824 TG_URB = TG_URB2D(I,J)
825 TC_URB = TC_URB2D(I,J)
826 QC_URB = QC_URB2D(I,J)
827 UC_URB = UC_URB2D(I,J)
828
829 DO K = 1,num_roof_layers
830 TRL_URB(K) = TRL_URB3D(I,K,J)
831 END DO
832 DO K = 1,num_wall_layers
833 TBL_URB(K) = TBL_URB3D(I,K,J)
834 END DO
835 DO K = 1,num_road_layers
836 TGL_URB(K) = TGL_URB3D(I,K,J)
837 END DO
838
839 XXXR_URB = XXXR_URB2D(I,J)
840 XXXB_URB = XXXB_URB2D(I,J)
841 XXXG_URB = XXXG_URB2D(I,J)
842 XXXC_URB = XXXC_URB2D(I,J)
843 !
844 CHS_URB = CHS(I,J)
845 CHS2_URB = CHS2(I,J)
846 !
847
848 ! Call urban
849
850
851 CALL urban(LSOLAR_URB, & ! I
852 num_roof_layers,num_wall_layers,num_road_layers, & ! C
853 DZR,DZB,DZG, & ! C
854 UTYPE_URB,TA_URB,QA_URB,UA_URB,U1_URB,V1_URB,SSG_URB, & ! I
855 SSGD_URB,SSGQ_URB,LLG_URB,RAIN_URB,RHOO_URB, & ! I
856 ZA_URB,DECLIN_URB,COSZ_URB,OMG_URB, & ! I
857 XLAT_URB,DELT_URB,ZNT_URB, & ! I
858 CHS_URB, CHS2_URB, & ! I
859 TR_URB, TB_URB, TG_URB, TC_URB, QC_URB,UC_URB, & ! H
860 TRL_URB,TBL_URB,TGL_URB, & ! H
861 XXXR_URB, XXXB_URB, XXXG_URB, XXXC_URB, & ! H
862 TS_URB,QS_URB,SH_URB,LH_URB,LH_KINEMATIC_URB, & ! O
863 SW_URB,ALB_URB,LW_URB,G_URB,RN_URB,PSIM_URB,PSIH_URB, & ! O
864 GZ1OZ0_URB, & !O
865 U10_URB, V10_URB, TH2_URB, Q2_URB, & ! O
866 UST_URB) !O
867
868
869 IF(IPRINT) THEN
870
871 print*, 'AFTER CALL URBAN'
872 print*,'num_roof_layers',num_roof_layers, 'num_wall_layers', &
873 num_wall_layers, &
874 'DZR',DZR,'DZB',DZB,'DZG',DZG,'UTYPE_URB',UTYPE_URB,'TA_URB', &
875 TA_URB, &
876 'QA_URB',QA_URB,'UA_URB',UA_URB,'U1_URB',U1_URB,'V1_URB', &
877 V1_URB, &
878 'SSG_URB',SSG_URB,'SSGD_URB',SSGD_URB,'SSGQ_URB',SSGQ_URB, &
879 'LLG_URB',LLG_URB,'RAIN_URB',RAIN_URB,'RHOO_URB',RHOO_URB, &
880 'ZA_URB',ZA_URB, 'DECLIN_URB',DECLIN_URB,'COSZ_URB',COSZ_URB,&
881 'OMG_URB',OMG_URB,'XLAT_URB',XLAT_URB,'DELT_URB',DELT_URB, &
882 'ZNT_URB',ZNT_URB,'TR_URB',TR_URB, 'TB_URB',TB_URB,'TG_URB',&
883 TG_URB,'TC_URB',TC_URB,'QC_URB',QC_URB,'TRL_URB',TRL_URB, &
884 'TBL_URB',TBL_URB,'TGL_URB',TGL_URB,'XXXR_URB',XXXR_URB, &
885 'XXXB_URB',XXXB_URB,'XXXG_URB',XXXG_URB,'XXXC_URB',XXXC_URB,&
886 'TS_URB',TS_URB,'QS_URB',QS_URB,'SH_URB',SH_URB,'LH_URB', &
887 LH_URB, 'LH_KINEMATIC_URB',LH_KINEMATIC_URB,'SW_URB',SW_URB,&
888 'ALB_URB',ALB_URB,'LW_URB',LW_URB,'G_URB',G_URB,'RN_URB', &
889 RN_URB, 'PSIM_URB',PSIM_URB,'PSIH_URB',PSIH_URB, &
890 'U10_URB',U10_URB,'V10_URB',V10_URB,'TH2_URB',TH2_URB, &
891 'Q2_URB',Q2_URB,'CHS_URB',CHS_URB,'CHS2_URB',CHS2_URB
892 endif
893
894 TS_URB2D(I,J) = TS_URB
895
896 ALBEDO(I,J) = FRC_URB2D(I,J)*ALB_URB+(1-FRC_URB2D(I,J))*ALBEDOK ![-]
897 HFX(I,J) = FRC_URB2D(I,J)*SH_URB+(1-FRC_URB2D(I,J))*SHEAT ![W/m/m]
898 QFX(I,J) = FRC_URB2D(I,J)*LH_KINEMATIC_URB &
899 + (1-FRC_URB2D(I,J))*ETA_KINEMATIC ![kg/m/m/s]
900 LH(I,J) = FRC_URB2D(I,J)*LH_URB+(1-FRC_URB2D(I,J))*ETA ![W/m/m]
901 GRDFLX(I,J) = FRC_URB2D(I,J)*G_URB+(1-FRC_URB2D(I,J))*SSOIL ![W/m/m]
902 TSK(I,J) = FRC_URB2D(I,J)*TS_URB+(1-FRC_URB2D(I,J))*T1 ![K]
903 QSFC(I,J)= FRC_URB2D(I,J)*QS_URB+(1-FRC_URB2D(I,J))*Q1 ![-]
904
905 IF(IPRINT)THEN
906
907 print*, ' FRC_URB2D', FRC_URB2D, &
908 'ALB_URB',ALB_URB, 'ALBEDOK',ALBEDOK, &
909 'ALBEDO(I,J)', ALBEDO(I,J), &
910 'SH_URB',SH_URB,'SHEAT',SHEAT, 'HFX(I,J)',HFX(I,J), &
911 'LH_KINEMATIC_URB',LH_KINEMATIC_URB,'ETA_KINEMATIC', &
912 ETA_KINEMATIC, 'QFX(I,J)',QFX(I,J), &
913 'LH_URB',LH_URB, 'ETA',ETA, 'LH(I,J)',LH(I,J), &
914 'G_URB',G_URB,'SSOIL',SSOIL,'GRDFLX(I,J)', GRDFLX(I,J),&
915 'TS_URB',TS_URB,'T1',T1,'TSK(I,J)',TSK(I,J), &
916 'QS_URB',QS_URB,'Q1',Q1,'QSFC(I,J)',QSFC(I,J)
917 endif
918
919
920
921
922 ! Renew Urban State Varialbes
923
924 TR_URB2D(I,J) = TR_URB
925 TB_URB2D(I,J) = TB_URB
926 TG_URB2D(I,J) = TG_URB
927 TC_URB2D(I,J) = TC_URB
928 QC_URB2D(I,J) = QC_URB
929 UC_URB2D(I,J) = UC_URB
930
931 DO K = 1,num_roof_layers
932 TRL_URB3D(I,K,J) = TRL_URB(K)
933 END DO
934 DO K = 1,num_wall_layers
935 TBL_URB3D(I,K,J) = TBL_URB(K)
936 END DO
937 DO K = 1,num_road_layers
938 TGL_URB3D(I,K,J) = TGL_URB(K)
939 END DO
940 XXXR_URB2D(I,J) = XXXR_URB
941 XXXB_URB2D(I,J) = XXXB_URB
942 XXXG_URB2D(I,J) = XXXG_URB
943 XXXC_URB2D(I,J) = XXXC_URB
944
945 SH_URB2D(I,J) = SH_URB
946 LH_URB2D(I,J) = LH_URB
947 G_URB2D(I,J) = G_URB
948 RN_URB2D(I,J) = RN_URB
949 PSIM_URB2D(I,J) = PSIM_URB
950 PSIH_URB2D(I,J) = PSIH_URB
951 GZ1OZ0_URB2D(I,J)= GZ1OZ0_URB
952 U10_URB2D(I,J) = U10_URB
953 V10_URB2D(I,J) = V10_URB
954 TH2_URB2D(I,J) = TH2_URB
955 Q2_URB2D(I,J) = Q2_URB
956 UST_URB2D(I,J) = UST_URB
957 AKMS_URB2D(I,J) = KARMAN * UST_URB2D(I,J)/(GZ1OZ0_URB2D(I,J)-PSIM_URB2D(I,J))
958
959 END IF
960
961 ENDIF ! end of UCM CALL if block
962 !--------------------------------------
963 ! Urban Part End - urban
964 !--------------------------------------
965
966 !*** DIAGNOSTICS
967 SMSTAV(I,J)=SOILW
968 SMSTOT(I,J)=SOILM*1000.
969 ! Convert the water unit into mm
970 SFCRUNOFF(I,J)=SFCRUNOFF(I,J)+RUNOFF1*DT*1000.0
971 UDRUNOFF(I,J)=UDRUNOFF(I,J)+(RUNOFF2+RUNOFF3)*DT*1000.0
972 ! snow defined when fraction of frozen precip (FFROZP) > 0.5,
973 IF(FFROZP.GT.0.5)THEN
974 ACSNOW(I,J)=ACSNOW(I,J)+PRCP*DT
975 ENDIF
976 IF(SNOW(I,J).GT.0.)THEN
977 ACSNOM(I,J)=ACSNOM(I,J)+SNOMLT*1000.
978 ! accumulated snow-melt energy
979 SNOPCX(I,J)=SNOPCX(I,J)-SNOMLT/FDTLIW
980 ENDIF
981
982 ENDIF ! endif of land-sea test
983
984 100 CONTINUE ! of I loop
985
986 ENDDO ! of J loop
987 !------------------------------------------------------
988 END SUBROUTINE lsm
989 !------------------------------------------------------
990
991 SUBROUTINE LSMINIT(VEGFRA,SNOW,SNOWC,SNOWH,CANWAT,SMSTAV, &
992 SMSTOT, SFCRUNOFF,UDRUNOFF,ACSNOW, &
993 ACSNOM,IVGTYP,ISLTYP,TSLB,SMOIS,SH2O,ZS,DZS, &
994 FNDSOILW, FNDSNOWH, &
995 num_soil_layers, restart, &
996 allowed_to_read , &
997 ids,ide, jds,jde, kds,kde, &
998 ims,ime, jms,jme, kms,kme, &
999 its,ite, jts,jte, kts,kte )
1000
1001 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
1002 ims,ime, jms,jme, kms,kme, &
1003 its,ite, jts,jte, kts,kte
1004
1005 INTEGER, INTENT(IN) :: num_soil_layers
1006
1007 LOGICAL , INTENT(IN) :: restart , allowed_to_read
1008
1009 REAL, DIMENSION( num_soil_layers), INTENT(INOUT) :: ZS, DZS
1010
1011 REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , &
1012 INTENT(INOUT) :: SMOIS, & !Total soil moisture
1013 SH2O, & !liquid soil moisture
1014 TSLB !STEMP
1015
1016 REAL, DIMENSION( ims:ime, jms:jme ) , &
1017 INTENT(INOUT) :: SNOW, &
1018 SNOWH, &
1019 SNOWC, &
1020 CANWAT, &
1021 SMSTAV, &
1022 SMSTOT, &
1023 SFCRUNOFF, &
1024 UDRUNOFF, &
1025 ACSNOW, &
1026 VEGFRA, &
1027 ACSNOM
1028
1029 INTEGER, DIMENSION( ims:ime, jms:jme ) , &
1030 INTENT(IN) :: IVGTYP, &
1031 ISLTYP
1032
1033 LOGICAL, INTENT(IN) :: FNDSOILW , &
1034 FNDSNOWH
1035
1036 INTEGER :: L
1037 REAL :: BX, SMCMAX, PSISAT, FREE
1038 REAL, PARAMETER :: BLIM = 5.5, HLICE = 3.335E5, &
1039 GRAV = 9.81, T0 = 273.15
1040 INTEGER :: errflag
1041
1042 !
1043
1044
1045 ! initialize three Noah LSM related tables
1046 IF ( allowed_to_read ) THEN
1047 CALL wrf_message( 'INITIALIZE THREE Noah LSM RELATED TABLES' )
1048 CALL LSM_PARM_INIT
1049 ENDIF
1050
1051 IF(.not.restart)THEN
1052
1053 itf=min0(ite,ide-1)
1054 jtf=min0(jte,jde-1)
1055
1056 errflag = 0
1057 DO j = jts,jtf
1058 DO i = its,itf
1059 IF ( ISLTYP( i,j ) .LT. 1 ) THEN
1060 errflag = 1
1061 WRITE(err_message,*)"module_sf_noahlsm.F: lsminit: out of range ISLTYP ",i,j,ISLTYP( i,j )
1062 CALL wrf_message(err_message)
1063 ENDIF
1064 ENDDO
1065 ENDDO
1066 IF ( errflag .EQ. 1 ) THEN
1067 CALL wrf_error_fatal( "module_sf_noahlsm.F: lsminit: out of range value "// &
1068 "of ISLTYP. Is this field in the input?" )
1069 ENDIF
1070
1071 ! initialize soil liquid water content SH2O
1072
1073 ! IF(.NOT.FNDSOILW) THEN
1074
1075 ! If no SWC, do the following
1076 ! PRINT *,'SOIL WATER NOT FOUND - VALUE SET IN LSMINIT'
1077 DO J = jts,jtf
1078 DO I = its,itf
1079 BX = BB(ISLTYP(I,J))
1080 SMCMAX = MAXSMC(ISLTYP(I,J))
1081 PSISAT = SATPSI(ISLTYP(I,J))
1082 if ((bx > 0.0).and.(smcmax > 0.0).and.(psisat > 0.0)) then
1083 DO NS=1, num_soil_layers
1084 ! ----------------------------------------------------------------------
1085 !SH2O <= SMOIS for T < 273.149K (-0.001C)
1086 IF (TSLB(I,NS,J) < 273.149) THEN
1087 ! ----------------------------------------------------------------------
1088 ! first guess following explicit solution for Flerchinger Eqn from Koren
1089 ! et al, JGR, 1999, Eqn 17 (KCOUNT=0 in FUNCTION FRH2O).
1090 ! ISLTPK is soil type
1091 BX = BB(ISLTYP(I,J))
1092 SMCMAX = MAXSMC(ISLTYP(I,J))
1093 PSISAT = SATPSI(ISLTYP(I,J))
1094 IF ( BX > BLIM ) BX = BLIM
1095 FK=(( (HLICE/(GRAV*(-PSISAT))) * &
1096 ((TSLB(I,NS,J)-T0)/TSLB(I,NS,J)) )**(-1/BX) )*SMCMAX
1097 IF (FK < 0.02) FK = 0.02
1098 SH2O(I,NS,J) = MIN( FK, SMOIS(I,NS,J) )
1099 ! ----------------------------------------------------------------------
1100 ! now use iterative solution for liquid soil water content using
1101 ! FUNCTION FRH2O with the initial guess for SH2O from above explicit
1102 ! first guess.
1103 CALL FRH2O (FREE,TSLB(I,NS,J),SMOIS(I,NS,J),SH2O(I,NS,J), &
1104 SMCMAX,BX,PSISAT)
1105 SH2O(I,NS,J) = FREE
1106 ELSE ! of IF (TSLB(I,NS,J)
1107 ! ----------------------------------------------------------------------
1108 ! SH2O = SMOIS ( for T => 273.149K (-0.001C)
1109 SH2O(I,NS,J)=SMOIS(I,NS,J)
1110 ! ----------------------------------------------------------------------
1111 ENDIF ! of IF (TSLB(I,NS,J)
1112 END DO ! of DO NS=1, num_soil_layers
1113 else ! of if ((bx > 0.0)
1114 DO NS=1, num_soil_layers
1115 SH2O(I,NS,J)=SMOIS(I,NS,J)
1116 END DO
1117 endif ! of if ((bx > 0.0)
1118 ENDDO ! DO I = its,itf
1119 ENDDO ! DO J = jts,jtf
1120 ! ENDIF ! of IF(.NOT.FNDSOILW)THEN
1121
1122 ! initialize physical snow height SNOWH
1123
1124 IF(.NOT.FNDSNOWH)THEN
1125 ! If no SNOWH do the following
1126 CALL wrf_message( 'SNOW HEIGHT NOT FOUND - VALUE DEFINED IN LSMINIT' )
1127 DO J = jts,jtf
1128 DO I = its,itf
1129 SNOWH(I,J)=SNOW(I,J)*0.005 ! SNOW in mm and SNOWH in m
1130 ENDDO
1131 ENDDO
1132 ENDIF
1133
1134 ! initialize canopy water to ZERO
1135
1136 ! GO TO 110
1137 ! print*,'Note that canopy water content (CANWAT) is set to ZERO in LSMINIT'
1138 DO J = jts,jtf
1139 DO I = its,itf
1140 CANWAT(I,J)=0.0
1141 ENDDO
1142 ENDDO
1143 110 CONTINUE
1144
1145 ENDIF
1146 !------------------------------------------------------------------------------
1147 END SUBROUTINE lsminit
1148 !------------------------------------------------------------------------------
1149
1150
1151
1152 !
1153 !-----------------------------------------------------------------
1154 SUBROUTINE LSM_PARM_INIT
1155 !-----------------------------------------------------------------
1156
1157 character*4 :: MMINLU, MMINSL
1158
1159 MMINLU='USGS'
1160 MMINSL='STAS'
1161 call SOIL_VEG_GEN_PARM( MMINLU, MMINSL)
1162
1163 !-----------------------------------------------------------------
1164 END SUBROUTINE LSM_PARM_INIT
1165 !-----------------------------------------------------------------
1166
1167 !-----------------------------------------------------------------
1168 SUBROUTINE SOIL_VEG_GEN_PARM( MMINLU, MMINSL)
1169 !-----------------------------------------------------------------
1170
1171 IMPLICIT NONE
1172
1173 integer :: LUMATCH, IINDEX, LC, NUM_SLOPE
1174 integer :: ierr
1175 INTEGER , PARAMETER :: OPEN_OK = 0
1176
1177 character*4 :: MMINLU, MMINSL
1178 character*128 :: mess , message
1179 logical, external :: wrf_dm_on_monitor
1180
1181
1182 !-----SPECIFY VEGETATION RELATED CHARACTERISTICS :
1183 ! ALBBCK: SFC albedo (in percentage)
1184 ! Z0: Roughness length (m)
1185 ! SHDFAC: Green vegetation fraction (in percentage)
1186 ! Note: The ALBEDO, Z0, and SHDFAC values read from the following table
1187 ! ALBEDO, amd Z0 are specified in LAND-USE TABLE; and SHDFAC is
1188 ! the monthly green vegetation data
1189 ! CMXTBL: MAX CNPY Capacity (m)
1190 ! NROTBL: Rooting depth (layer)
1191 ! RSMIN: Mimimum stomatal resistance (s m-1)
1192 ! RSMAX: Max. stomatal resistance (s m-1)
1193 ! RGL: Parameters used in radiation stress function
1194 ! HS: Parameter used in vapor pressure deficit functio
1195 ! TOPT: Optimum transpiration air temperature. (K)
1196 ! CMCMAX: Maximum canopy water capacity
1197 ! CFACTR: Parameter used in the canopy inteception calculati
1198 ! SNUP: Threshold snow depth (in water equivalent m) that
1199 ! implies 100% snow cover
1200 ! LAI: Leaf area index (dimensionless)
1201 ! MAXALB: Upper bound on maximum albedo over deep snow
1202 !
1203 !-----READ IN VEGETAION PROPERTIES FROM VEGPARM.TBL
1204 !
1205
1206 IF ( wrf_dm_on_monitor() ) THEN
1207
1208 OPEN(19, FILE='VEGPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
1209 IF(ierr .NE. OPEN_OK ) THEN
1210 WRITE(message,FMT='(A)') &
1211 'module_sf_noahlsm.F: soil_veg_gen_parm: failure opening VEGPARM.TBL'
1212 CALL wrf_error_fatal ( message )
1213 END IF
1214
1215 WRITE ( mess, * ) 'INPUT LANDUSE = ',MMINLU
1216 CALL wrf_message( mess )
1217
1218 LUMATCH=0
1219
1220 READ (19,*)
1221 READ (19,2000,END=2002)LUTYPE
1222 READ (19,*)LUCATS,IINDEX
1223 2000 FORMAT (A4)
1224
1225 IF(LUTYPE.EQ.MMINLU)THEN
1226 WRITE( mess , * ) 'LANDUSE TYPE = ',LUTYPE,' FOUND', &
1227 LUCATS,' CATEGORIES'
1228 CALL wrf_message( mess )
1229 LUMATCH=1
1230 ENDIF
1231
1232 IF(LUTYPE.EQ.MMINLU)THEN
1233 DO LC=1,LUCATS
1234 READ (19,*)IINDEX,ALBTBL(LC),Z0TBL(LC),SHDTBL(LC), &
1235 NROTBL(LC),RSTBL(LC),RGLTBL(LC),HSTBL(LC), &
1236 SNUPTBL(LC),LAITBL(LC),MAXALB(LC)
1237 ENDDO
1238 !
1239 READ (19,*)
1240 READ (19,*)TOPT_DATA
1241 READ (19,*)
1242 READ (19,*)CMCMAX_DATA
1243 READ (19,*)
1244 READ (19,*)CFACTR_DATA
1245 READ (19,*)
1246 READ (19,*)RSMAX_DATA
1247 READ (19,*)
1248 READ (19,*)BARE
1249 ENDIF
1250 !
1251 2002 CONTINUE
1252
1253 CLOSE (19)
1254 ENDIF
1255
1256 CALL wrf_dm_bcast_string ( LUTYPE , 4 )
1257 CALL wrf_dm_bcast_integer ( LUCATS , 1 )
1258 CALL wrf_dm_bcast_integer ( IINDEX , 1 )
1259 CALL wrf_dm_bcast_integer ( LUMATCH , 1 )
1260 CALL wrf_dm_bcast_real ( ALBTBL , NLUS )
1261 CALL wrf_dm_bcast_real ( Z0TBL , NLUS )
1262 CALL wrf_dm_bcast_real ( SHDTBL , NLUS )
1263 CALL wrf_dm_bcast_real ( NROTBL , NLUS )
1264 CALL wrf_dm_bcast_real ( RSTBL , NLUS )
1265 CALL wrf_dm_bcast_real ( RGLTBL , NLUS )
1266 CALL wrf_dm_bcast_real ( HSTBL , NLUS )
1267 CALL wrf_dm_bcast_real ( SNUPTBL , NLUS )
1268 CALL wrf_dm_bcast_real ( LAITBL , NLUS )
1269 CALL wrf_dm_bcast_real ( MAXALB , NLUS )
1270 CALL wrf_dm_bcast_real ( TOPT_DATA , 1 )
1271 CALL wrf_dm_bcast_real ( CMCMAX_DATA , 1 )
1272 CALL wrf_dm_bcast_real ( CFACTR_DATA , 1 )
1273 CALL wrf_dm_bcast_real ( RSMAX_DATA , 1 )
1274 CALL wrf_dm_bcast_integer ( BARE , 1 )
1275
1276 !
1277 !-----READ IN SOIL PROPERTIES FROM SOILPARM.TBL
1278 !
1279 IF ( wrf_dm_on_monitor() ) THEN
1280 OPEN(19, FILE='SOILPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
1281 IF(ierr .NE. OPEN_OK ) THEN
1282 WRITE(message,FMT='(A)') &
1283 'module_sf_noahlsm.F: soil_veg_gen_parm: failure opening SOILPARM.TBL'
1284 CALL wrf_error_fatal ( message )
1285 END IF
1286
1287 MMINSL='STAS' !oct2
1288 WRITE(mess,*) 'INPUT SOIL TEXTURE CLASSIFICAION = ',MMINSL
1289 CALL wrf_message( mess )
1290
1291 LUMATCH=0
1292
1293 READ (19,*)
1294 READ (19,2000,END=2003)SLTYPE
1295 READ (19,*)SLCATS,IINDEX
1296 IF(SLTYPE.EQ.MMINSL)THEN
1297 WRITE( mess , * ) 'SOIL TEXTURE CLASSIFICATION = ',SLTYPE,' FOUND', &
1298 SLCATS,' CATEGORIES'
1299 CALL wrf_message ( mess )
1300 LUMATCH=1
1301 ENDIF
1302 IF(SLTYPE.EQ.MMINSL)THEN
1303 DO LC=1,SLCATS
1304 READ (19,*) IINDEX,BB(LC),DRYSMC(LC),F11(LC),MAXSMC(LC),&
1305 REFSMC(LC),SATPSI(LC),SATDK(LC), SATDW(LC), &
1306 WLTSMC(LC), QTZ(LC)
1307 ENDDO
1308 ENDIF
1309
1310 2003 CONTINUE
1311
1312 CLOSE (19)
1313 ENDIF
1314
1315 CALL wrf_dm_bcast_integer ( LUMATCH , 1 )
1316 CALL wrf_dm_bcast_string ( SLTYPE , 4 )
1317 CALL wrf_dm_bcast_string ( MMINSL , 4 ) ! since this is reset above, see oct2 ^
1318 CALL wrf_dm_bcast_integer ( SLCATS , 1 )
1319 CALL wrf_dm_bcast_integer ( IINDEX , 1 )
1320 CALL wrf_dm_bcast_real ( BB , NSLTYPE )
1321 CALL wrf_dm_bcast_real ( DRYSMC , NSLTYPE )
1322 CALL wrf_dm_bcast_real ( F11 , NSLTYPE )
1323 CALL wrf_dm_bcast_real ( MAXSMC , NSLTYPE )
1324 CALL wrf_dm_bcast_real ( REFSMC , NSLTYPE )
1325 CALL wrf_dm_bcast_real ( SATPSI , NSLTYPE )
1326 CALL wrf_dm_bcast_real ( SATDK , NSLTYPE )
1327 CALL wrf_dm_bcast_real ( SATDW , NSLTYPE )
1328 CALL wrf_dm_bcast_real ( WLTSMC , NSLTYPE )
1329 CALL wrf_dm_bcast_real ( QTZ , NSLTYPE )
1330
1331 IF(LUMATCH.EQ.0)THEN
1332 CALL wrf_message( 'SOIl TEXTURE IN INPUT FILE DOES NOT ' )
1333 CALL wrf_message( 'MATCH SOILPARM TABLE' )
1334 CALL wrf_error_fatal ( 'INCONSISTENT OR MISSING SOILPARM FILE' )
1335 ENDIF
1336
1337 !
1338 !-----READ IN GENERAL PARAMETERS FROM GENPARM.TBL
1339 !
1340 IF ( wrf_dm_on_monitor() ) THEN
1341 OPEN(19, FILE='GENPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
1342 IF(ierr .NE. OPEN_OK ) THEN
1343 WRITE(message,FMT='(A)') &
1344 'module_sf_noahlsm.F: soil_veg_gen_parm: failure opening GENPARM.TBL'
1345 CALL wrf_error_fatal ( message )
1346 END IF
1347
1348 READ (19,*)
1349 READ (19,*)
1350 READ (19,*) NUM_SLOPE
1351
1352 SLPCATS=NUM_SLOPE
1353
1354 DO LC=1,SLPCATS
1355 READ (19,*)SLOPE_DATA(LC)
1356 ENDDO
1357
1358 READ (19,*)
1359 READ (19,*)SBETA_DATA
1360 READ (19,*)
1361 READ (19,*)FXEXP_DATA
1362 READ (19,*)
1363 READ (19,*)CSOIL_DATA
1364 READ (19,*)
1365 READ (19,*)SALP_DATA
1366 READ (19,*)
1367 READ (19,*)REFDK_DATA
1368 READ (19,*)
1369 READ (19,*)REFKDT_DATA
1370 READ (19,*)
1371 READ (19,*)FRZK_DATA
1372 READ (19,*)
1373 READ (19,*)ZBOT_DATA
1374 READ (19,*)
1375 READ (19,*)CZIL_DATA
1376 READ (19,*)
1377 READ (19,*)SMLOW_DATA
1378 READ (19,*)
1379 READ (19,*)SMHIGH_DATA
1380 CLOSE (19)
1381 ENDIF
1382
1383 CALL wrf_dm_bcast_integer ( NUM_SLOPE , 1 )
1384 CALL wrf_dm_bcast_integer ( SLPCATS , 1 )
1385 CALL wrf_dm_bcast_real ( SLOPE_DATA , NSLOPE )
1386 CALL wrf_dm_bcast_real ( SBETA_DATA , 1 )
1387 CALL wrf_dm_bcast_real ( FXEXP_DATA , 1 )
1388 CALL wrf_dm_bcast_real ( CSOIL_DATA , 1 )
1389 CALL wrf_dm_bcast_real ( SALP_DATA , 1 )
1390 CALL wrf_dm_bcast_real ( REFDK_DATA , 1 )
1391 CALL wrf_dm_bcast_real ( REFKDT_DATA , 1 )
1392 CALL wrf_dm_bcast_real ( FRZK_DATA , 1 )
1393 CALL wrf_dm_bcast_real ( ZBOT_DATA , 1 )
1394 CALL wrf_dm_bcast_real ( CZIL_DATA , 1 )
1395 CALL wrf_dm_bcast_real ( SMLOW_DATA , 1 )
1396 CALL wrf_dm_bcast_real ( SMHIGH_DATA , 1 )
1397
1398
1399 !-----------------------------------------------------------------
1400 END SUBROUTINE SOIL_VEG_GEN_PARM
1401 !-----------------------------------------------------------------
1402
1403 SUBROUTINE SFLX (FFROZP,ICE,DT,ZLVL,NSOIL,SLDPTH, & !C
1404 LOCAL, & !L
1405 LLANDUSE, LSOIL, & !CL
1406 LWDN,SOLDN,SOLNET,SFCPRS,PRCP,SFCTMP,Q2,SFCSPD, & !F
1407 COSZ,PRCPRAIN, SOLARDIRECT, & !F
1408 TH2,Q2SAT,DQSDT2, & !I
1409 VEGTYP,SOILTYP,SLOPETYP,SHDFAC,SHDMIN,SHDMAX, & !I
1410 ALB, SNOALB,TBOT, Z0BRD, Z0, EMISSI, & !S
1411 CMC,T1,STC,SMC,SH2O,SNOWH,SNEQV,ALBEDO,CH,CM, & !H
1412 ! ----------------------------------------------------------------------
1413 ! OUTPUTS, DIAGNOSTICS, PARAMETERS BELOW GENERALLY NOT NECESSARY WHEN
1414 ! COUPLED WITH E.G. A NWP MODEL (SUCH AS THE NOAA/NWS/NCEP MESOSCALE ETA
1415 ! MODEL). OTHER APPLICATIONS MAY REQUIRE DIFFERENT OUTPUT VARIABLES.
1416 ! ----------------------------------------------------------------------
1417 ETA,SHEAT, ETA_KINEMATIC,FDOWN, & !O
1418 EC,EDIR,ET,ETT,ESNOW,DRIP,DEW, & !O
1419 BETA,ETP,SSOIL, & !O
1420 FLX1,FLX2,FLX3, & !O
1421 SNOMLT,SNCOVR, & !O
1422 RUNOFF1,RUNOFF2,RUNOFF3, & !O
1423 RC,PC,RSMIN,XLAI,RCS,RCT,RCQ,RCSOIL, & !O
1424 SOILW,SOILM,Q1, & !D
1425 SMCWLT,SMCDRY,SMCREF,SMCMAX,NROOT) !P
1426 ! ----------------------------------------------------------------------
1427 ! SUBROUTINE SFLX - UNIFIED NOAHLSM VERSION 1.0 JULY 2007
1428 ! ----------------------------------------------------------------------
1429 ! SUB-DRIVER FOR "Noah LSM" FAMILY OF PHYSICS SUBROUTINES FOR A
1430 ! SOIL/VEG/SNOWPACK LAND-SURFACE MODEL TO UPDATE SOIL MOISTURE, SOIL
1431 ! ICE, SOIL TEMPERATURE, SKIN TEMPERATURE, SNOWPACK WATER CONTENT,
1432 ! SNOWDEPTH, AND ALL TERMS OF THE SURFACE ENERGY BALANCE AND SURFACE
1433 ! WATER BALANCE (EXCLUDING INPUT ATMOSPHERIC FORCINGS OF DOWNWARD
1434 ! RADIATION AND PRECIP)
1435 ! ----------------------------------------------------------------------
1436 ! SFLX ARGUMENT LIST KEY:
1437 ! ----------------------------------------------------------------------
1438 ! C CONFIGURATION INFORMATION
1439 ! L LOGICAL
1440 ! CL 4-string character bearing logical meaning
1441 ! F FORCING DATA
1442 ! I OTHER (INPUT) FORCING DATA
1443 ! S SURFACE CHARACTERISTICS
1444 ! H HISTORY (STATE) VARIABLES
1445 ! O OUTPUT VARIABLES
1446 ! D DIAGNOSTIC OUTPUT
1447 ! P Parameters
1448 ! Msic Miscellaneous terms passed from gridded driver
1449 ! ----------------------------------------------------------------------
1450 ! 1. CONFIGURATION INFORMATION (C):
1451 ! ----------------------------------------------------------------------
1452 ! ICE SEA-ICE FLAG (=1: SEA-ICE, =0: LAND)
1453 ! DT TIMESTEP (SEC) (DT SHOULD NOT EXCEED 3600 SECS, RECOMMEND
1454 ! 1800 SECS OR LESS)
1455 ! ZLVL HEIGHT (M) ABOVE GROUND OF ATMOSPHERIC FORCING VARIABLES
1456 ! NSOIL NUMBER OF SOIL LAYERS (AT LEAST 2, AND NOT GREATER THAN
1457 ! PARAMETER NSOLD SET BELOW)
1458 ! SLDPTH THE THICKNESS OF EACH SOIL LAYER (M)
1459 ! ----------------------------------------------------------------------
1460 ! 2. LOGICAL:
1461 ! ----------------------------------------------------------------------
1462 ! LCH Exchange coefficient (Ch) calculation flag (false: using
1463 ! ch-routine SFCDIF; true: Ch is brought in)
1464 ! LOCAL Flag for local-site simulation (where there is no
1465 ! maps for albedo, veg fraction, and roughness
1466 ! true: all LSM parameters (inluding albedo, veg fraction and
1467 ! roughness length) will be defined by three tables
1468 ! LLANDUSE (=USGS, using USGS landuse classification)
1469 ! LSOIL (=STAS, using FAO/STATSGO soil texture classification)
1470 ! ----------------------------------------------------------------------
1471 ! 3. FORCING DATA (F):
1472 ! ----------------------------------------------------------------------
1473 ! LWDN LW DOWNWARD RADIATION (W M-2; POSITIVE, NOT NET LONGWAVE)
1474 ! SOLDN SOLAR DOWNWARD RADIATION (W M-2; POSITIVE, NOT NET SOLAR)
1475 ! SFCPRS PRESSURE AT HEIGHT ZLVL ABOVE GROUND (PASCALS)
1476 ! PRCP PRECIP RATE (KG M-2 S-1) (NOTE, THIS IS A RATE)
1477 ! SFCTMP AIR TEMPERATURE (K) AT HEIGHT ZLVL ABOVE GROUND
1478 ! TH2 AIR POTENTIAL TEMPERATURE (K) AT HEIGHT ZLVL ABOVE GROUND
1479 ! Q2 MIXING RATIO AT HEIGHT ZLVL ABOVE GROUND (KG KG-1)
1480 ! COSZ Solar zenith angle (not used for now)
1481 ! PRCPRAIN Liquid-precipitation rate (KG M-2 S-1) (not used)
1482 ! SOLARDIRECT Direct component of downward solar radiation (W M-2) (not used)
1483 ! ----------------------------------------------------------------------
1484 ! 4. OTHER FORCING (INPUT) DATA (I):
1485 ! ----------------------------------------------------------------------
1486 ! SFCSPD WIND SPEED (M S-1) AT HEIGHT ZLVL ABOVE GROUND
1487 ! Q2SAT SAT SPECIFIC HUMIDITY AT HEIGHT ZLVL ABOVE GROUND (KG KG-1)
1488 ! DQSDT2 SLOPE OF SAT SPECIFIC HUMIDITY CURVE AT T=SFCTMP
1489 ! (KG KG-1 K-1)
1490 ! ----------------------------------------------------------------------
1491 ! 5. CANOPY/SOIL CHARACTERISTICS (S):
1492 ! ----------------------------------------------------------------------
1493 ! VEGTYP VEGETATION TYPE (INTEGER INDEX)
1494 ! SOILTYP SOIL TYPE (INTEGER INDEX)
1495 ! SLOPETYP CLASS OF SFC SLOPE (INTEGER INDEX)
1496 ! SHDFAC AREAL FRACTIONAL COVERAGE OF GREEN VEGETATION
1497 ! (FRACTION= 0.0-1.0)
1498 ! SHDMIN MINIMUM AREAL FRACTIONAL COVERAGE OF GREEN VEGETATION
1499 ! (FRACTION= 0.0-1.0) <= SHDFAC
1500 ! PTU PHOTO THERMAL UNIT (PLANT PHENOLOGY FOR ANNUALS/CROPS)
1501 ! (NOT YET USED, BUT PASSED TO REDPRM FOR FUTURE USE IN
1502 ! VEG PARMS)
1503 ! ALB BACKROUND SNOW-FREE SURFACE ALBEDO (FRACTION), FOR JULIAN
1504 ! DAY OF YEAR (USUALLY FROM TEMPORAL INTERPOLATION OF
1505 ! MONTHLY MEAN VALUES' CALLING PROG MAY OR MAY NOT
1506 ! INCLUDE DIURNAL SUN ANGLE EFFECT)
1507 ! SNOALB UPPER BOUND ON MAXIMUM ALBEDO OVER DEEP SNOW (E.G. FROM
1508 ! ROBINSON AND KUKLA, 1985, J. CLIM. & APPL. METEOR.)
1509 ! TBOT BOTTOM SOIL TEMPERATURE (LOCAL YEARLY-MEAN SFC AIR
1510 ! TEMPERATURE)
1511 ! Z0BRD Background fixed roughness length (M)
1512 ! Z0 Time varying roughness length (M) as function of snow depth
1513 ! ----------------------------------------------------------------------
1514 ! 6. HISTORY (STATE) VARIABLES (H):
1515 ! ----------------------------------------------------------------------
1516 ! CMC CANOPY MOISTURE CONTENT (M)
1517 ! T1 GROUND/CANOPY/SNOWPACK) EFFECTIVE SKIN TEMPERATURE (K)
1518 ! STC(NSOIL) SOIL TEMP (K)
1519 ! SMC(NSOIL) TOTAL SOIL MOISTURE CONTENT (VOLUMETRIC FRACTION)
1520 ! SH2O(NSOIL) UNFROZEN SOIL MOISTURE CONTENT (VOLUMETRIC FRACTION)
1521 ! NOTE: FROZEN SOIL MOISTURE = SMC - SH2O
1522 ! SNOWH ACTUAL SNOW DEPTH (M)
1523 ! SNEQV LIQUID WATER-EQUIVALENT SNOW DEPTH (M)
1524 ! NOTE: SNOW DENSITY = SNEQV/SNOWH
1525 ! ALBEDO SURFACE ALBEDO INCLUDING SNOW EFFECT (UNITLESS FRACTION)
1526 ! =SNOW-FREE ALBEDO (ALB) WHEN SNEQV=0, OR
1527 ! =FCT(MSNOALB,ALB,VEGTYP,SHDFAC,SHDMIN) WHEN SNEQV>0
1528 ! CH SURFACE EXCHANGE COEFFICIENT FOR HEAT AND MOISTURE
1529 ! (M S-1); NOTE: CH IS TECHNICALLY A CONDUCTANCE SINCE
1530 ! IT HAS BEEN MULTIPLIED BY WIND SPEED.
1531 ! CM SURFACE EXCHANGE COEFFICIENT FOR MOMENTUM (M S-1); NOTE:
1532 ! CM IS TECHNICALLY A CONDUCTANCE SINCE IT HAS BEEN
1533 ! MULTIPLIED BY WIND SPEED.
1534 ! ----------------------------------------------------------------------
1535 ! 7. OUTPUT (O):
1536 ! ----------------------------------------------------------------------
1537 ! OUTPUT VARIABLES NECESSARY FOR A COUPLED NUMERICAL WEATHER PREDICTION
1538 ! MODEL, E.G. NOAA/NWS/NCEP MESOSCALE ETA MODEL. FOR THIS APPLICATION,
1539 ! THE REMAINING OUTPUT/DIAGNOSTIC/PARAMETER BLOCKS BELOW ARE NOT
1540 ! NECESSARY. OTHER APPLICATIONS MAY REQUIRE DIFFERENT OUTPUT VARIABLES.
1541 ! ETA ACTUAL LATENT HEAT FLUX (W m-2: NEGATIVE, IF UP FROM
1542 ! SURFACE)
1543 ! ETA_KINEMATIC atctual latent heat flux in Kg m-2 s-1
1544 ! SHEAT SENSIBLE HEAT FLUX (W M-2: NEGATIVE, IF UPWARD FROM
1545 ! SURFACE)
1546 ! FDOWN Radiation forcing at the surface (W m-2) = SOLDN*(1-alb)+LWDN
1547 ! ----------------------------------------------------------------------
1548 ! EC CANOPY WATER EVAPORATION (W m-2)
1549 ! EDIR DIRECT SOIL EVAPORATION (W m-2)
1550 ! ET(NSOIL) PLANT TRANSPIRATION FROM A PARTICULAR ROOT (SOIL) LAYER
1551 ! (W m-2)
1552 ! ETT TOTAL PLANT TRANSPIRATION (W m-2)
1553 ! ESNOW SUBLIMATION FROM (OR DEPOSITION TO IF <0) SNOWPACK
1554 ! (W m-2)
1555 ! DRIP THROUGH-FALL OF PRECIP AND/OR DEW IN EXCESS OF CANOPY
1556 ! WATER-HOLDING CAPACITY (M)
1557 ! DEW DEWFALL (OR FROSTFALL FOR T<273.15) (M)
1558 ! ----------------------------------------------------------------------
1559 ! BETA RATIO OF ACTUAL/POTENTIAL EVAP (DIMENSIONLESS)
1560 ! ETP POTENTIAL EVAPORATION (W m-2)
1561 ! SSOIL SOIL HEAT FLUX (W M-2: NEGATIVE IF DOWNWARD FROM SURFACE)
1562 ! ----------------------------------------------------------------------
1563 ! FLX1 PRECIP-SNOW SFC (W M-2)
1564 ! FLX2 FREEZING RAIN LATENT HEAT FLUX (W M-2)
1565 ! FLX3 PHASE-CHANGE HEAT FLUX FROM SNOWMELT (W M-2)
1566 ! ----------------------------------------------------------------------
1567 ! SNOMLT SNOW MELT (M) (WATER EQUIVALENT)
1568 ! SNCOVR FRACTIONAL SNOW COVER (UNITLESS FRACTION, 0-1)
1569 ! ----------------------------------------------------------------------
1570 ! RUNOFF1 SURFACE RUNOFF (M S-1), NOT INFILTRATING THE SURFACE
1571 ! RUNOFF2 SUBSURFACE RUNOFF (M S-1), DRAINAGE OUT BOTTOM OF LAST
1572 ! SOIL LAYER (BASEFLOW)
1573 ! RUNOFF3 NUMERICAL TRUNCTATION IN EXCESS OF POROSITY (SMCMAX)
1574 ! FOR A GIVEN SOIL LAYER AT THE END OF A TIME STEP (M S-1).
1575 ! Note: the above RUNOFF2 is actually the sum of RUNOFF2 and RUNOFF3
1576 ! ----------------------------------------------------------------------
1577 ! RC CANOPY RESISTANCE (S M-1)
1578 ! PC PLANT COEFFICIENT (UNITLESS FRACTION, 0-1) WHERE PC*ETP
1579 ! = ACTUAL TRANSP
1580 ! XLAI LEAF AREA INDEX (DIMENSIONLESS)
1581 ! RSMIN MINIMUM CANOPY RESISTANCE (S M-1)
1582 ! RCS INCOMING SOLAR RC FACTOR (DIMENSIONLESS)
1583 ! RCT AIR TEMPERATURE RC FACTOR (DIMENSIONLESS)
1584 ! RCQ ATMOS VAPOR PRESSURE DEFICIT RC FACTOR (DIMENSIONLESS)
1585 ! RCSOIL SOIL MOISTURE RC FACTOR (DIMENSIONLESS)
1586 ! ----------------------------------------------------------------------
1587 ! 8. DIAGNOSTIC OUTPUT (D):
1588 ! ----------------------------------------------------------------------
1589 ! SOILW AVAILABLE SOIL MOISTURE IN ROOT ZONE (UNITLESS FRACTION
1590 ! BETWEEN SMCWLT AND SMCMAX)
1591 ! SOILM TOTAL SOIL COLUMN MOISTURE CONTENT (FROZEN+UNFROZEN) (M)
1592 ! Q1 Effective mixing ratio at surface (kg kg-1), used for
1593 ! diagnosing the mixing ratio at 2 meter for coupled model
1594 ! ----------------------------------------------------------------------
1595 ! 9. PARAMETERS (P):
1596 ! ----------------------------------------------------------------------
1597 ! SMCWLT WILTING POINT (VOLUMETRIC)
1598 ! SMCDRY DRY SOIL MOISTURE THRESHOLD WHERE DIRECT EVAP FRM TOP
1599 ! LAYER ENDS (VOLUMETRIC)
1600 ! SMCREF SOIL MOISTURE THRESHOLD WHERE TRANSPIRATION BEGINS TO
1601 ! STRESS (VOLUMETRIC)
1602 ! SMCMAX POROSITY, I.E. SATURATED VALUE OF SOIL MOISTURE
1603 ! (VOLUMETRIC)
1604 ! NROOT NUMBER OF ROOT LAYERS, A FUNCTION OF VEG TYPE, DETERMINED
1605 ! IN SUBROUTINE REDPRM.
1606 ! ----------------------------------------------------------------------
1607
1608
1609 IMPLICIT NONE
1610 ! ----------------------------------------------------------------------
1611
1612 ! DECLARATIONS - LOGICAL AND CHARACTERS
1613 ! ----------------------------------------------------------------------
1614 LOGICAL, INTENT(IN):: LOCAL
1615 LOGICAL :: FRZGRA, SNOWNG
1616 CHARACTER (LEN=4), INTENT(IN):: LLANDUSE, LSOIL
1617
1618 ! ----------------------------------------------------------------------
1619 ! DECLARATIONS - INTEGER
1620 ! ----------------------------------------------------------------------
1621 INTEGER,INTENT(IN) :: ICE,NSOIL,SLOPETYP,SOILTYP,VEGTYP
1622 INTEGER,INTENT(OUT):: NROOT
1623 INTEGER KZ, K, iout
1624
1625 ! ----------------------------------------------------------------------
1626 ! DECLARATIONS - REAL
1627 ! ----------------------------------------------------------------------
1628 REAL, INTENT(IN) :: SHDMIN,SHDMAX,DT,DQSDT2,LWDN,PRCP,PRCPRAIN, &
1629 Q2,Q2SAT,SFCPRS,SFCSPD,SFCTMP, SNOALB, &
1630 SOLDN,SOLNET,TBOT,TH2,ZLVL, &
1631 EMISSI, FFROZP
1632 REAL, INTENT(INOUT):: COSZ, SOLARDIRECT,ALBEDO,CH,CM, &
1633 CMC,SNEQV,SNCOVR,SNOWH,T1,XLAI,SHDFAC,Z0BRD,ALB
1634 REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SLDPTH
1635 REAL, DIMENSION(1:NSOIL), INTENT(OUT):: ET
1636 REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: SH2O, SMC, STC
1637 REAL,DIMENSION(1:NSOIL):: RTDIS, ZSOIL
1638
1639 REAL,INTENT(OUT) :: ETA_KINEMATIC,BETA,DEW,DRIP,EC,EDIR,ESNOW,ETA, &
1640 ETP,FLX1,FLX2,FLX3,SHEAT,PC,RUNOFF1,RUNOFF2, &
1641 RUNOFF3,RC,RSMIN,RCQ,RCS,RCSOIL,RCT,SSOIL, &
1642 SMCDRY,SMCMAX,SMCREF,SMCWLT,SNOMLT, SOILM, &
1643 SOILW,FDOWN,Q1
1644 REAL :: BEXP,CFACTR,CMCMAX,CSOIL,CZIL,DF1,DF1H,DF1A,DKSAT,DWSAT, &
1645 DSOIL,DTOT,ETT,FRCSNO,FRCSOI,EPSCA,F1,FXEXP,FRZX,HS, &
1646 KDT,LVH2O,PRCP1,PSISAT,QUARTZ,R,RCH,REFKDT,RR,RGL, &
1647 RSMAX, &
1648 RSNOW,SNDENS,SNCOND,SBETA,SN_NEW,SLOPE,SNUP,SALP,SOILWM, &
1649 SOILWW,T1V,T24,T2V,TH2V,TOPT,TFREEZ,TSNOW,ZBOT,Z0,PRCPF, &
1650 ETNS,PTU,LSUBS
1651
1652 ! ----------------------------------------------------------------------
1653 ! DECLARATIONS - PARAMETERS
1654 ! ----------------------------------------------------------------------
1655 PARAMETER (TFREEZ = 273.15)
1656 PARAMETER (LVH2O = 2.501E+6)
1657 PARAMETER (LSUBS = 2.83E+6)
1658 PARAMETER (R = 287.04)
1659 ! ----------------------------------------------------------------------
1660 ! INITIALIZATION
1661 ! ----------------------------------------------------------------------
1662 RUNOFF1 = 0.0
1663 RUNOFF2 = 0.0
1664 RUNOFF3 = 0.0
1665 SNOMLT = 0.0
1666
1667 ! ----------------------------------------------------------------------
1668 ! THE VARIABLE "ICE" IS A FLAG DENOTING SEA-ICE CASE
1669 ! ----------------------------------------------------------------------
1670 ! SEA-ICE LAYERS ARE EQUAL THICKNESS AND SUM TO 3 METERS
1671 ! ----------------------------------------------------------------------
1672 IF (ICE == 1) THEN
1673 DO KZ = 1,NSOIL
1674 ZSOIL (KZ) = -3.* FLOAT (KZ)/ FLOAT (NSOIL)
1675 END DO
1676
1677 ! ----------------------------------------------------------------------
1678 ! CALCULATE DEPTH (NEGATIVE) BELOW GROUND FROM TOP SKIN SFC TO BOTTOM OF
1679 ! EACH SOIL LAYER. NOTE: SIGN OF ZSOIL IS NEGATIVE (DENOTING BELOW
1680 ! GROUND)
1681 ! ----------------------------------------------------------------------
1682 ELSE
1683 ZSOIL (1) = - SLDPTH (1)
1684 DO KZ = 2,NSOIL
1685 ZSOIL (KZ) = - SLDPTH (KZ) + ZSOIL (KZ -1)
1686 END DO
1687 ! ----------------------------------------------------------------------
1688 ! NEXT IS CRUCIAL CALL TO SET THE LAND-SURFACE PARAMETERS, INCLUDING
1689 ! SOIL-TYPE AND VEG-TYPE DEPENDENT PARAMETERS.
1690 ! ----------------------------------------------------------------------
1691 CALL REDPRM (VEGTYP,SOILTYP,SLOPETYP,CFACTR,CMCMAX,RSMAX,TOPT, &
1692 REFKDT,KDT,SBETA, SHDFAC,RSMIN,RGL,HS,ZBOT,FRZX, &
1693 PSISAT,SLOPE,SNUP,SALP,BEXP,DKSAT,DWSAT, &
1694 SMCMAX,SMCWLT,SMCREF,SMCDRY,F1,QUARTZ,FXEXP, &
1695 RTDIS,SLDPTH,ZSOIL,NROOT,NSOIL,Z0BRD,CZIL,XLAI, &
1696 CSOIL,ALB,PTU,LLANDUSE,LSOIL,LOCAL)
1697 END IF
1698
1699 !urban
1700 IF(VEGTYP==1)THEN
1701 SHDFAC=0.05
1702 RSMIN=400.0
1703 SMCMAX = 0.45
1704 SMCREF = 0.42
1705 SMCWLT = 0.40
1706 ENDIF
1707
1708 ! ----------------------------------------------------------------------
1709 ! INITIALIZE PRECIPITATION LOGICALS.
1710 ! ----------------------------------------------------------------------
1711 SNOWNG = .FALSE.
1712 FRZGRA = .FALSE.
1713
1714 ! ----------------------------------------------------------------------
1715 ! IF SEA-ICE CASE, ASSIGN DEFAULT WATER-EQUIV SNOW ON TOP
1716 ! ----------------------------------------------------------------------
1717 IF (ICE == 1) THEN
1718 SNEQV = 0.01
1719 SNOWH = 0.05
1720 ! ----------------------------------------------------------------------
1721 ! IF INPUT SNOWPACK IS NONZERO, THEN COMPUTE SNOW DENSITY "SNDENS" AND
1722 ! SNOW THERMAL CONDUCTIVITY "SNCOND" (NOTE THAT CSNOW IS A FUNCTION
1723 ! SUBROUTINE)
1724 ! ----------------------------------------------------------------------
1725 END IF
1726 IF (SNEQV == 0.0) THEN
1727 SNDENS = 0.0
1728 SNOWH = 0.0
1729 SNCOND = 1.0
1730 ELSE
1731 SNDENS = SNEQV / SNOWH
1732 IF(SNDENS > 1.0) THEN
1733 CALL wrf_error_fatal ( 'Physical snow depth is less than snow water equiv.' )
1734 ENDIF
1735 CALL CSNOW (SNCOND,SNDENS)
1736 END IF
1737 ! ----------------------------------------------------------------------
1738 ! DETERMINE IF IT'S PRECIPITATING AND WHAT KIND OF PRECIP IT IS.
1739 ! IF IT'S PRCPING AND THE AIR TEMP IS COLDER THAN 0 C, IT'S SNOWING!
1740 ! IF IT'S PRCPING AND THE AIR TEMP IS WARMER THAN 0 C, BUT THE GRND
1741 ! TEMP IS COLDER THAN 0 C, FREEZING RAIN IS PRESUMED TO BE FALLING.
1742 ! ----------------------------------------------------------------------
1743 IF (PRCP > 0.0) THEN
1744 ! snow defined when fraction of frozen precip (FFROZP) > 0.5,
1745 ! passed in from model microphysics.
1746 IF (FFROZP .GT. 0.5) THEN
1747 SNOWNG = .TRUE.
1748 ELSE
1749 IF (T1 <= TFREEZ) FRZGRA = .TRUE.
1750 END IF
1751 END IF
1752 ! ----------------------------------------------------------------------
1753 ! IF EITHER PRCP FLAG IS SET, DETERMINE NEW SNOWFALL (CONVERTING PRCP
1754 ! RATE FROM KG M-2 S-1 TO A LIQUID EQUIV SNOW DEPTH IN METERS) AND ADD
1755 ! IT TO THE EXISTING SNOWPACK.
1756 ! NOTE THAT SINCE ALL PRECIP IS ADDED TO SNOWPACK, NO PRECIP INFILTRATES
1757 ! INTO THE SOIL SO THAT PRCP1 IS SET TO ZERO.
1758 ! ----------------------------------------------------------------------
1759 IF ( (SNOWNG) .OR. (FRZGRA) ) THEN
1760 SN_NEW = PRCP * DT * 0.001
1761 SNEQV = SNEQV + SN_NEW
1762 PRCPF = 0.0
1763
1764 ! ----------------------------------------------------------------------
1765 ! UPDATE SNOW DENSITY BASED ON NEW SNOWFALL, USING OLD AND NEW SNOW.
1766 ! UPDATE SNOW THERMAL CONDUCTIVITY
1767 ! ----------------------------------------------------------------------
1768 CALL SNOW_NEW (SFCTMP,SN_NEW,SNOWH,SNDENS)
1769 CALL CSNOW (SNCOND,SNDENS)
1770
1771 ! ----------------------------------------------------------------------
1772 ! PRECIP IS LIQUID (RAIN), HENCE SAVE IN THE PRECIP VARIABLE THAT
1773 ! LATER CAN WHOLELY OR PARTIALLY INFILTRATE THE SOIL (ALONG WITH
1774 ! ANY CANOPY "DRIP" ADDED TO THIS LATER)
1775 ! ----------------------------------------------------------------------
1776 ELSE
1777 PRCPF = PRCP
1778 END IF
1779 ! ----------------------------------------------------------------------
1780 ! DETERMINE SNOWCOVER AND ALBEDO OVER LAND.
1781 ! ----------------------------------------------------------------------
1782 ! ----------------------------------------------------------------------
1783 ! IF SNOW DEPTH=0, SET SNOW FRACTION=0, ALBEDO=SNOW FREE ALBEDO.
1784 ! ----------------------------------------------------------------------
1785 IF (ICE == 0) THEN
1786 IF (SNEQV == 0.0) THEN
1787 SNCOVR = 0.0
1788 ALBEDO = ALB
1789 ELSE
1790 ! ----------------------------------------------------------------------
1791 ! DETERMINE SNOW FRACTIONAL COVERAGE.
1792 ! DETERMINE SURFACE ALBEDO MODIFICATION DUE TO SNOWDEPTH STATE.
1793 ! ----------------------------------------------------------------------
1794 CALL SNFRAC (SNEQV,SNUP,SALP,SNOWH,SNCOVR)
1795 ! limit snow cover fraction to maximum of 0.98
1796 SNCOVR = MIN(SNCOVR,0.98)
1797 CALL ALCALC (ALB,SNOALB,SHDFAC,SHDMIN,SNCOVR,TSNOW,ALBEDO)
1798 END IF
1799 ! ----------------------------------------------------------------------
1800 ! SNOW COVER, ALBEDO OVER SEA-ICE
1801 ! ----------------------------------------------------------------------
1802 ELSE
1803 SNCOVR = 1.0
1804 ALBEDO = 0.65
1805 END IF
1806 ! ----------------------------------------------------------------------
1807 ! THERMAL CONDUCTIVITY FOR SEA-ICE CASE
1808 ! ----------------------------------------------------------------------
1809 IF (ICE == 1) THEN
1810 DF1 = 2.2
1811 ELSE
1812 ! ----------------------------------------------------------------------
1813 ! NEXT CALCULATE THE SUBSURFACE HEAT FLUX, WHICH FIRST REQUIRES
1814 ! CALCULATION OF THE THERMAL DIFFUSIVITY. TREATMENT OF THE
1815 ! LATTER FOLLOWS THAT ON PAGES 148-149 FROM "HEAT TRANSFER IN
1816 ! COLD CLIMATES", BY V. J. LUNARDINI (PUBLISHED IN 1981
1817 ! BY VAN NOSTRAND REINHOLD CO.) I.E. TREATMENT OF TWO CONTIGUOUS
1818 ! "PLANE PARALLEL" MEDIUMS (NAMELY HERE THE FIRST SOIL LAYER
1819 ! AND THE SNOWPACK LAYER, IF ANY). THIS DIFFUSIVITY TREATMENT
1820 ! BEHAVES WELL FOR BOTH ZERO AND NONZERO SNOWPACK, INCLUDING THE
1821 ! LIMIT OF VERY THIN SNOWPACK. THIS TREATMENT ALSO ELIMINATES
1822 ! THE NEED TO IMPOSE AN ARBITRARY UPPER BOUND ON SUBSURFACE
1823 ! HEAT FLUX WHEN THE SNOWPACK BECOMES EXTREMELY THIN.
1824 ! ----------------------------------------------------------------------
1825 ! FIRST CALCULATE THERMAL DIFFUSIVITY OF TOP SOIL LAYER, USING
1826 ! BOTH THE FROZEN AND LIQUID SOIL MOISTURE, FOLLOWING THE
1827 ! SOIL THERMAL DIFFUSIVITY FUNCTION OF PETERS-LIDARD ET AL.
1828 ! (1998,JAS, VOL 55, 1209-1224), WHICH REQUIRES THE SPECIFYING
1829 ! THE QUARTZ CONTENT OF THE GIVEN SOIL CLASS (SEE ROUTINE REDPRM)
1830 ! ----------------------------------------------------------------------
1831 ! ----------------------------------------------------------------------
1832 ! NEXT ADD SUBSURFACE HEAT FLUX REDUCTION EFFECT FROM THE
1833 ! OVERLYING GREEN CANOPY, ADAPTED FROM SECTION 2.1.2 OF
1834 ! PETERS-LIDARD ET AL. (1997, JGR, VOL 102(D4))
1835 ! ----------------------------------------------------------------------
1836 CALL TDFCND (DF1,SMC (1),QUARTZ,SMCMAX,SH2O (1))
1837
1838 !urban
1839 IF (VEGTYP==1) DF1=3.24
1840
1841 DF1 = DF1 * EXP (SBETA * SHDFAC)
1842 ! ----------------------------------------------------------------------
1843 ! FINALLY "PLANE PARALLEL" SNOWPACK EFFECT FOLLOWING
1844 ! V.J. LINARDINI REFERENCE CITED ABOVE. NOTE THAT DTOT IS
1845 ! COMBINED DEPTH OF SNOWDEPTH AND THICKNESS OF FIRST SOIL LAYER
1846 ! ----------------------------------------------------------------------
1847 END IF
1848
1849 DSOIL = - (0.5 * ZSOIL (1))
1850 IF (SNEQV == 0.) THEN
1851 SSOIL = DF1 * (T1- STC (1) ) / DSOIL
1852 ELSE
1853 DTOT = SNOWH + DSOIL
1854 FRCSNO = SNOWH / DTOT
1855
1856 ! 1. HARMONIC MEAN (SERIES FLOW)
1857 ! DF1 = (SNCOND*DF1)/(FRCSOI*SNCOND+FRCSNO*DF1)
1858 FRCSOI = DSOIL / DTOT
1859 ! 2. ARITHMETIC MEAN (PARALLEL FLOW)
1860 ! DF1 = FRCSNO*SNCOND + FRCSOI*DF1
1861 DF1H = (SNCOND * DF1)/ (FRCSOI * SNCOND+ FRCSNO * DF1)
1862
1863 ! 3. GEOMETRIC MEAN (INTERMEDIATE BETWEEN HARMONIC AND ARITHMETIC MEAN)
1864 ! DF1 = (SNCOND**FRCSNO)*(DF1**FRCSOI)
1865 ! weigh DF by snow fraction
1866 ! DF1 = DF1H*SNCOVR + DF1A*(1.0-SNCOVR)
1867 ! DF1 = DF1H*SNCOVR + DF1*(1.0-SNCOVR)
1868 DF1A = FRCSNO * SNCOND+ FRCSOI * DF1
1869
1870 ! ----------------------------------------------------------------------
1871 ! CALCULATE SUBSURFACE HEAT FLUX, SSOIL, FROM FINAL THERMAL DIFFUSIVITY
1872 ! OF SURFACE MEDIUMS, DF1 ABOVE, AND SKIN TEMPERATURE AND TOP
1873 ! MID-LAYER SOIL TEMPERATURE
1874 ! ----------------------------------------------------------------------
1875 DF1 = DF1A * SNCOVR + DF1* (1.0- SNCOVR)
1876 SSOIL = DF1 * (T1- STC (1) ) / DTOT
1877 END IF
1878 ! ----------------------------------------------------------------------
1879 ! DETERMINE SURFACE ROUGHNESS OVER SNOWPACK USING SNOW CONDITION FROM
1880 ! THE PREVIOUS TIMESTEP.
1881 ! ----------------------------------------------------------------------
1882 IF (SNCOVR > 0. ) THEN
1883 CALL SNOWZ0 (SNCOVR,Z0,Z0BRD)
1884 ELSE
1885 Z0=Z0BRD
1886 END IF
1887 ! ----------------------------------------------------------------------
1888 ! NEXT CALL ROUTINE SFCDIF TO CALCULATE THE SFC EXCHANGE COEF (CH) FOR
1889 ! HEAT AND MOISTURE.
1890
1891 ! NOTE !!!
1892 ! DO NOT CALL SFCDIF UNTIL AFTER ABOVE CALL TO REDPRM, IN CASE
1893 ! ALTERNATIVE VALUES OF ROUGHNESS LENGTH (Z0) AND ZILINTINKEVICH COEF
1894 ! (CZIL) ARE SET THERE VIA NAMELIST I/O.
1895
1896 ! NOTE !!!
1897 ! ROUTINE SFCDIF RETURNS A CH THAT REPRESENTS THE WIND SPD TIMES THE
1898 ! "ORIGINAL" NONDIMENSIONAL "Ch" TYPICAL IN LITERATURE. HENCE THE CH
1899 ! RETURNED FROM SFCDIF HAS UNITS OF M/S. THE IMPORTANT COMPANION
1900 ! COEFFICIENT OF CH, CARRIED HERE AS "RCH", IS THE CH FROM SFCDIF TIMES
1901 ! AIR DENSITY AND PARAMETER "CP". "RCH" IS COMPUTED IN "CALL PENMAN".
1902 ! RCH RATHER THAN CH IS THE COEFF USUALLY INVOKED LATER IN EQNS.
1903
1904 ! NOTE !!!
1905 ! ----------------------------------------------------------------------
1906 ! SFCDIF ALSO RETURNS THE SURFACE EXCHANGE COEFFICIENT FOR MOMENTUM, CM,
1907 ! ALSO KNOWN AS THE SURFACE DRAGE COEFFICIENT. Needed as a state variable
1908 ! for iterative/implicit solution of CH in SFCDIF
1909 ! ----------------------------------------------------------------------
1910 ! IF(.NOT.LCH) THEN
1911 ! T1V = T1 * (1.0+ 0.61 * Q2)
1912 ! TH2V = TH2 * (1.0+ 0.61 * Q2)
1913 ! CALL SFCDIF_off (ZLVL,Z0,T1V,TH2V,SFCSPD,CZIL,CM,CH)
1914 ! ENDIF
1915
1916 ! ----------------------------------------------------------------------
1917 ! CALL PENMAN SUBROUTINE TO CALCULATE POTENTIAL EVAPORATION (ETP), AND
1918 ! OTHER PARTIAL PRODUCTS AND SUMS SAVE IN COMMON/RITE FOR LATER
1919 ! CALCULATIONS.
1920 ! ----------------------------------------------------------------------
1921 ! ----------------------------------------------------------------------
1922 ! CALCULATE TOTAL DOWNWARD RADIATION (SOLAR PLUS LONGWAVE) NEEDED IN
1923 ! PENMAN EP SUBROUTINE THAT FOLLOWS
1924 ! ----------------------------------------------------------------------
1925 ! FDOWN = SOLDN * (1.0- ALBEDO) + LWDN
1926 FDOWN = SOLNET + LWDN
1927 ! ----------------------------------------------------------------------
1928 ! CALC VIRTUAL TEMPS AND VIRTUAL POTENTIAL TEMPS NEEDED BY SUBROUTINES
1929 ! PENMAN.
1930 T2V = SFCTMP * (1.0+ 0.61 * Q2 )
1931
1932 iout=0
1933 if(iout.eq.1) then
1934 print*,'before penman'
1935 print*,' SFCTMP',SFCTMP,'SFCPRS',SFCPRS,'CH',CH,'T2V',T2V, &
1936 'TH2',TH2,'PRCP',PRCP,'FDOWN',FDOWN,'T24',T24,'SSOIL',SSOIL, &
1937 'Q2',Q2,'Q2SAT',Q2SAT,'ETP',ETP,'RCH',RCH, &
1938 'EPSCA',EPSCA,'RR',RR ,'SNOWNG',SNOWNG,'FRZGRA',FRZGRA, &
1939 'DQSDT2',DQSDT2,'FLX2',FLX2,'SNOWH',SNOWH,'SNEQV',SNEQV, &
1940 ' DSOIL',DSOIL,' FRCSNO',FRCSNO,' SNCOVR',SNCOVR,' DTOT',DTOT, &
1941 ' ZSOIL (1)',ZSOIL(1),' DF1',DF1,'T1',T1,' STC1',STC(1), &
1942 'ALBEDO',ALBEDO,'SMC',SMC,'STC',STC,'SH2O',SH2O
1943 endif
1944
1945 CALL PENMAN (SFCTMP,SFCPRS,CH,T2V,TH2,PRCP,FDOWN,T24,SSOIL, &
1946 Q2,Q2SAT,ETP,RCH,EPSCA,RR,SNOWNG,FRZGRA, &
1947 DQSDT2,FLX2,EMISSI,SNEQV)
1948
1949 ! ----------------------------------------------------------------------
1950 ! CALL CANRES TO CALCULATE THE CANOPY RESISTANCE AND CONVERT IT INTO PC
1951 ! IF NONZERO GREENNESS FRACTION
1952 ! ----------------------------------------------------------------------
1953
1954 ! ----------------------------------------------------------------------
1955 ! FROZEN GROUND EXTENSION: TOTAL SOIL WATER "SMC" WAS REPLACED
1956 ! BY UNFROZEN SOIL WATER "SH2O" IN CALL TO CANRES BELOW
1957 ! ----------------------------------------------------------------------
1958 IF (SHDFAC > 0.) THEN
1959 CALL CANRES (SOLDN,CH,SFCTMP,Q2,SFCPRS,SH2O,ZSOIL,NSOIL, &
1960 SMCWLT,SMCREF,RSMIN,RC,PC,NROOT,Q2SAT,DQSDT2, &
1961 TOPT,RSMAX,RGL,HS,XLAI, &
1962 RCS,RCT,RCQ,RCSOIL,EMISSI)
1963 END IF
1964 ! ----------------------------------------------------------------------
1965 ! NOW DECIDE MAJOR PATHWAY BRANCH TO TAKE DEPENDING ON WHETHER SNOWPACK
1966 ! EXISTS OR NOT:
1967 ! ----------------------------------------------------------------------
1968 ESNOW = 0.0
1969 IF (SNEQV == 0.0) THEN
1970 CALL NOPAC (ETP,ETA,PRCP,SMC,SMCMAX,SMCWLT, &
1971 SMCREF,SMCDRY,CMC,CMCMAX,NSOIL,DT, &
1972 SHDFAC, &
1973 SBETA,Q2,T1,SFCTMP,T24,TH2,FDOWN,F1,EMISSI, &
1974 SSOIL, &
1975 STC,EPSCA,BEXP,PC,RCH,RR,CFACTR, &
1976 SH2O,SLOPE,KDT,FRZX,PSISAT,ZSOIL, &
1977 DKSAT,DWSAT,TBOT,ZBOT,RUNOFF1,RUNOFF2, &
1978 RUNOFF3,EDIR,EC,ET,ETT,NROOT,ICE,RTDIS, &
1979 QUARTZ,FXEXP,CSOIL, &
1980 BETA,DRIP,DEW,FLX1,FLX3,VEGTYP)
1981 ETA_KINEMATIC = ETA
1982 ELSE
1983 CALL SNOPAC (ETP,ETA,PRCP,PRCPF,SNOWNG,SMC,SMCMAX,SMCWLT, &
1984 SMCREF,SMCDRY,CMC,CMCMAX,NSOIL,DT, &
1985 SBETA,DF1, &
1986 Q2,T1,SFCTMP,T24,TH2,FDOWN,F1,SSOIL,STC,EPSCA, &
1987 SFCPRS,BEXP,PC,RCH,RR,CFACTR,SNCOVR,SNEQV,SNDENS,&
1988 SNOWH,SH2O,SLOPE,KDT,FRZX,PSISAT, &
1989 ZSOIL,DWSAT,DKSAT,TBOT,ZBOT,SHDFAC,RUNOFF1, &
1990 RUNOFF2,RUNOFF3,EDIR,EC,ET,ETT,NROOT,SNOMLT, &
1991 ICE,RTDIS,QUARTZ,FXEXP,CSOIL, &
1992 BETA,DRIP,DEW,FLX1,FLX2,FLX3,ESNOW,ETNS,EMISSI, &
1993 VEGTYP)
1994 ETA_KINEMATIC = ESNOW + ETNS
1995 END IF
1996
1997 ! Calculate effective mixing ratio at grnd level (skin)
1998 !
1999 ! Q1=Q2+ETA*CP/RCH
2000 Q1=Q2+ETA_KINEMATIC*CP/RCH
2001 !
2002 ! ----------------------------------------------------------------------
2003 ! DETERMINE SENSIBLE HEAT (H) IN ENERGY UNITS (W M-2)
2004 ! ----------------------------------------------------------------------
2005 SHEAT = - (CH * CP * SFCPRS)/ (R * T2V) * ( TH2- T1 )
2006
2007 ! ----------------------------------------------------------------------
2008 ! CONVERT EVAP TERMS FROM KINEMATIC (KG M-2 S-1) TO ENERGY UNITS (W M-2)
2009 ! ----------------------------------------------------------------------
2010 EDIR = EDIR * LVH2O
2011 EC = EC * LVH2O
2012 DO K=1,4
2013 ET(K) = ET(K) * LVH2O
2014 ENDDO
2015 ETT = ETT * LVH2O
2016 ESNOW = ESNOW * LSUBS
2017 ETP = ETP*((1.-SNCOVR)*LVH2O + SNCOVR*LSUBS)
2018 IF (ETP .GT. 0.) THEN
2019 ETA = EDIR + EC + ETT + ESNOW
2020 ELSE
2021 ETA = ETP
2022 ENDIF
2023 ! ----------------------------------------------------------------------
2024 ! DETERMINE BETA (RATIO OF ACTUAL TO POTENTIAL EVAP)
2025 ! ----------------------------------------------------------------------
2026 IF (ETP == 0.0) THEN
2027 BETA = 0.0
2028 ELSE
2029 BETA = ETA/ETP
2030 ENDIF
2031
2032 ! ----------------------------------------------------------------------
2033 ! CONVERT THE SIGN OF SOIL HEAT FLUX SO THAT:
2034 ! SSOIL>0: WARM THE SURFACE (NIGHT TIME)
2035 ! SSOIL<0: COOL THE SURFACE (DAY TIME)
2036 ! ----------------------------------------------------------------------
2037 SSOIL = -1.0* SSOIL
2038
2039 ! ----------------------------------------------------------------------
2040 ! CONVERT RUNOFF3 (INTERNAL LAYER RUNOFF FROM SUPERSAT) FROM M TO M S-1
2041 ! AND ADD TO SUBSURFACE RUNOFF/DRAINAGE/BASEFLOW
2042 ! ----------------------------------------------------------------------
2043 RUNOFF3 = RUNOFF3/ DT
2044 RUNOFF2 = RUNOFF2+ RUNOFF3
2045
2046 ! ----------------------------------------------------------------------
2047 ! TOTAL COLUMN SOIL MOISTURE IN METERS (SOILM) AND ROOT-ZONE
2048 ! SOIL MOISTURE AVAILABILITY (FRACTION) RELATIVE TO POROSITY/SATURATION
2049 ! ----------------------------------------------------------------------
2050 IF (ICE == 0) THEN
2051 SOILM = -1.0* SMC (1)* ZSOIL (1)
2052 DO K = 2,NSOIL
2053 SOILM = SOILM + SMC (K)* (ZSOIL (K -1) - ZSOIL (K))
2054 END DO
2055 SOILWM = -1.0* (SMCMAX - SMCWLT)* ZSOIL (1)
2056 SOILWW = -1.0* (SMC (1) - SMCWLT)* ZSOIL (1)
2057
2058 IF (NROOT >= 2) THEN
2059 DO K = 2,NROOT
2060 SOILWM = SOILWM + (SMCMAX - SMCWLT)* (ZSOIL (K -1) - ZSOIL (K))
2061 SOILWW = SOILWW + (SMC(K) - SMCWLT)* (ZSOIL (K -1) - ZSOIL (K))
2062 END DO
2063 END IF
2064 IF (SOILWM .LT. 1.E-6) THEN
2065 SOILWM = 0.0
2066 SOILW = 0.0
2067 SOILM = 0.0
2068 ELSE
2069 SOILW = SOILWW / SOILWM
2070 END IF
2071 ELSE
2072 SOILWM = 0.0
2073 SOILW = 0.0
2074 SOILM = 0.0
2075 END IF
2076
2077 ! ----------------------------------------------------------------------
2078 END SUBROUTINE SFLX
2079 ! ----------------------------------------------------------------------
2080
2081 SUBROUTINE ALCALC (ALB,SNOALB,SHDFAC,SHDMIN,SNCOVR,TSNOW,ALBEDO)
2082
2083 ! ----------------------------------------------------------------------
2084 ! CALCULATE ALBEDO INCLUDING SNOW EFFECT (0 -> 1)
2085 ! ALB SNOWFREE ALBEDO
2086 ! SNOALB MAXIMUM (DEEP) SNOW ALBEDO
2087 ! SHDFAC AREAL FRACTIONAL COVERAGE OF GREEN VEGETATION
2088 ! SHDMIN MINIMUM AREAL FRACTIONAL COVERAGE OF GREEN VEGETATION
2089 ! SNCOVR FRACTIONAL SNOW COVER
2090 ! ALBEDO SURFACE ALBEDO INCLUDING SNOW EFFECT
2091 ! TSNOW SNOW SURFACE TEMPERATURE (K)
2092 ! ----------------------------------------------------------------------
2093 IMPLICIT NONE
2094
2095 ! ----------------------------------------------------------------------
2096 ! SNOALB IS ARGUMENT REPRESENTING MAXIMUM ALBEDO OVER DEEP SNOW,
2097 ! AS PASSED INTO SFLX, AND ADAPTED FROM THE SATELLITE-BASED MAXIMUM
2098 ! SNOW ALBEDO FIELDS PROVIDED BY D. ROBINSON AND G. KUKLA
2099 ! (1985, JCAM, VOL 24, 402-411)
2100 ! ----------------------------------------------------------------------
2101 REAL, INTENT(IN) :: ALB, SNOALB, SHDFAC, SHDMIN, SNCOVR, TSNOW
2102 REAL, INTENT(OUT) :: ALBEDO
2103 ! turn of vegetation effect
2104 ! ALBEDO = ALB + (1.0- (SHDFAC - SHDMIN))* SNCOVR * (SNOALB - ALB)
2105 ! ALBEDO = (1.0-SNCOVR)*ALB + SNCOVR*SNOALB !this is equivalent to below
2106 ALBEDO = ALB + SNCOVR*(SNOALB-ALB)
2107
2108 ! BASE FORMULATION (DICKINSON ET AL., 1986, COGLEY ET AL., 1990)
2109 ! IF (TSNOW.LE.263.16) THEN
2110 ! ALBEDO=SNOALB
2111 ! ELSE
2112 ! IF (TSNOW.LT.273.16) THEN
2113 ! TM=0.1*(TSNOW-263.16)
2114 ! ALBEDO=0.5*((0.9-0.2*(TM**3))+(0.8-0.16*(TM**3)))
2115 ! ELSE
2116 ! ALBEDO=0.67
2117 ! ENDIF
2118 ! ENDIF
2119
2120 ! ISBA FORMULATION (VERSEGHY, 1991; BAKER ET AL., 1990)
2121 ! IF (TSNOW.LT.273.16) THEN
2122 ! ALBEDO=SNOALB-0.008*DT/86400
2123 ! ELSE
2124 ! ALBEDO=(SNOALB-0.5)*EXP(-0.24*DT/86400)+0.5
2125 ! ENDIF
2126
2127 IF (ALBEDO > SNOALB) ALBEDO = SNOALB
2128
2129 ! ----------------------------------------------------------------------
2130 END SUBROUTINE ALCALC
2131 ! ----------------------------------------------------------------------
2132
2133 SUBROUTINE CANRES (SOLAR,CH,SFCTMP,Q2,SFCPRS,SMC,ZSOIL,NSOIL, &
2134 SMCWLT,SMCREF,RSMIN,RC,PC,NROOT,Q2SAT,DQSDT2, &
2135 TOPT,RSMAX,RGL,HS,XLAI, &
2136 RCS,RCT,RCQ,RCSOIL,EMISSI)
2137
2138 ! ----------------------------------------------------------------------
2139 ! SUBROUTINE CANRES
2140 ! ----------------------------------------------------------------------
2141 ! CALCULATE CANOPY RESISTANCE WHICH DEPENDS ON INCOMING SOLAR RADIATION,
2142 ! AIR TEMPERATURE, ATMOSPHERIC WATER VAPOR PRESSURE DEFICIT AT THE
2143 ! LOWEST MODEL LEVEL, AND SOIL MOISTURE (PREFERABLY UNFROZEN SOIL
2144 ! MOISTURE RATHER THAN TOTAL)
2145 ! ----------------------------------------------------------------------
2146 ! SOURCE: JARVIS (1976), NOILHAN AND PLANTON (1989, MWR), JACQUEMIN AND
2147 ! NOILHAN (1990, BLM)
2148 ! SEE ALSO: CHEN ET AL (1996, JGR, VOL 101(D3), 7251-7268), EQNS 12-14
2149 ! AND TABLE 2 OF SEC. 3.1.2
2150 ! ----------------------------------------------------------------------
2151 ! INPUT:
2152 ! SOLAR INCOMING SOLAR RADIATION
2153 ! CH SURFACE EXCHANGE COEFFICIENT FOR HEAT AND MOISTURE
2154 ! SFCTMP AIR TEMPERATURE AT 1ST LEVEL ABOVE GROUND
2155 ! Q2 AIR HUMIDITY AT 1ST LEVEL ABOVE GROUND
2156 ! Q2SAT SATURATION AIR HUMIDITY AT 1ST LEVEL ABOVE GROUND
2157 ! DQSDT2 SLOPE OF SATURATION HUMIDITY FUNCTION WRT TEMP
2158 ! SFCPRS SURFACE PRESSURE
2159 ! SMC VOLUMETRIC SOIL MOISTURE
2160 ! ZSOIL SOIL DEPTH (NEGATIVE SIGN, AS IT IS BELOW GROUND)
2161 ! NSOIL NO. OF SOIL LAYERS
2162 ! NROOT NO. OF SOIL LAYERS IN ROOT ZONE (1.LE.NROOT.LE.NSOIL)
2163 ! XLAI LEAF AREA INDEX
2164 ! SMCWLT WILTING POINT
2165 ! SMCREF REFERENCE SOIL MOISTURE (WHERE SOIL WATER DEFICIT STRESS
2166 ! SETS IN)
2167 ! RSMIN, RSMAX, TOPT, RGL, HS ARE CANOPY STRESS PARAMETERS SET IN
2168 ! SURBOUTINE REDPRM
2169 ! OUTPUT:
2170 ! PC PLANT COEFFICIENT
2171 ! RC CANOPY RESISTANCE
2172 ! ----------------------------------------------------------------------
2173
2174 IMPLICIT NONE
2175 INTEGER, INTENT(IN) :: NROOT,NSOIL
2176 INTEGER K
2177 REAL, INTENT(IN) :: CH,DQSDT2,HS,Q2,Q2SAT,RSMIN,RGL,RSMAX, &
2178 SFCPRS,SFCTMP,SMCREF,SMCWLT, SOLAR,TOPT,XLAI, &
2179 EMISSI
2180 REAL,DIMENSION(1:NSOIL), INTENT(IN) :: SMC,ZSOIL
2181 REAL, INTENT(OUT):: PC,RC,RCQ,RCS,RCSOIL,RCT
2182 REAL :: DELTA,FF,GX,P,RR
2183 REAL, DIMENSION(1:NSOIL) :: PART
2184 REAL, PARAMETER :: SLV = 2.501000E6
2185
2186
2187 ! ----------------------------------------------------------------------
2188 ! INITIALIZE CANOPY RESISTANCE MULTIPLIER TERMS.
2189 ! ----------------------------------------------------------------------
2190 RCS = 0.0
2191 RCT = 0.0
2192 RCQ = 0.0
2193 RCSOIL = 0.0
2194
2195 ! ----------------------------------------------------------------------
2196 ! CONTRIBUTION DUE TO INCOMING SOLAR RADIATION
2197 ! ----------------------------------------------------------------------
2198 RC = 0.0
2199 FF = 0.55*2.0* SOLAR / (RGL * XLAI)
2200 RCS = (FF + RSMIN / RSMAX) / (1.0+ FF)
2201
2202 ! ----------------------------------------------------------------------
2203 ! CONTRIBUTION DUE TO AIR TEMPERATURE AT FIRST MODEL LEVEL ABOVE GROUND
2204 ! RCT EXPRESSION FROM NOILHAN AND PLANTON (1989, MWR).
2205 ! ----------------------------------------------------------------------
2206 RCS = MAX (RCS,0.0001)
2207 RCT = 1.0- 0.0016* ( (TOPT - SFCTMP)**2.0)
2208
2209 ! ----------------------------------------------------------------------
2210 ! CONTRIBUTION DUE TO VAPOR PRESSURE DEFICIT AT FIRST MODEL LEVEL.
2211 ! RCQ EXPRESSION FROM SSIB
2212 ! ----------------------------------------------------------------------
2213 RCT = MAX (RCT,0.0001)
2214 RCQ = 1.0/ (1.0+ HS * (Q2SAT - Q2))
2215
2216 ! ----------------------------------------------------------------------
2217 ! CONTRIBUTION DUE TO SOIL MOISTURE AVAILABILITY.
2218 ! DETERMINE CONTRIBUTION FROM EACH SOIL LAYER, THEN ADD THEM UP.
2219 ! ----------------------------------------------------------------------
2220 RCQ = MAX (RCQ,0.01)
2221 GX = (SMC (1) - SMCWLT) / (SMCREF - SMCWLT)
2222 IF (GX > 1.) GX = 1.
2223 IF (GX < 0.) GX = 0.
2224
2225 ! ----------------------------------------------------------------------
2226 ! USE SOIL DEPTH AS WEIGHTING FACTOR
2227 ! ----------------------------------------------------------------------
2228 ! ----------------------------------------------------------------------
2229 ! USE ROOT DISTRIBUTION AS WEIGHTING FACTOR
2230 ! PART(1) = RTDIS(1) * GX
2231 ! ----------------------------------------------------------------------
2232 PART (1) = (ZSOIL (1)/ ZSOIL (NROOT)) * GX
2233 DO K = 2,NROOT
2234 GX = (SMC (K) - SMCWLT) / (SMCREF - SMCWLT)
2235 IF (GX > 1.) GX = 1.
2236 IF (GX < 0.) GX = 0.
2237 ! ----------------------------------------------------------------------
2238 ! USE SOIL DEPTH AS WEIGHTING FACTOR
2239 ! ----------------------------------------------------------------------
2240 ! ----------------------------------------------------------------------
2241 ! USE ROOT DISTRIBUTION AS WEIGHTING FACTOR
2242 ! PART(K) = RTDIS(K) * GX
2243 ! ----------------------------------------------------------------------
2244 PART (K) = ( (ZSOIL (K) - ZSOIL (K -1))/ ZSOIL (NROOT)) * GX
2245 END DO
2246 DO K = 1,NROOT
2247 RCSOIL = RCSOIL + PART (K)
2248 END DO
2249
2250 ! ----------------------------------------------------------------------
2251 ! DETERMINE CANOPY RESISTANCE DUE TO ALL FACTORS. CONVERT CANOPY
2252 ! RESISTANCE (RC) TO PLANT COEFFICIENT (PC) TO BE USED WITH POTENTIAL
2253 ! EVAP IN DETERMINING ACTUAL EVAP. PC IS DETERMINED BY:
2254 ! PC * LINERIZED PENMAN POTENTIAL EVAP =
2255 ! PENMAN-MONTEITH ACTUAL EVAPORATION (CONTAINING RC TERM).
2256 ! ----------------------------------------------------------------------
2257 RCSOIL = MAX (RCSOIL,0.0001)
2258
2259 RC = RSMIN / (XLAI * RCS * RCT * RCQ * RCSOIL)
2260 ! RR = (4.* SIGMA * RD / CP)* (SFCTMP **4.)/ (SFCPRS * CH) + 1.0
2261 RR = (4.* EMISSI *SIGMA * RD / CP)* (SFCTMP **4.)/ (SFCPRS * CH) &
2262 + 1.0
2263
2264 DELTA = (SLV / CP)* DQSDT2
2265
2266 PC = (RR + DELTA)/ (RR * (1. + RC * CH) + DELTA)
2267
2268 ! ----------------------------------------------------------------------
2269 END SUBROUTINE CANRES
2270 ! ----------------------------------------------------------------------
2271
2272 SUBROUTINE CSNOW (SNCOND,DSNOW)
2273
2274 ! ----------------------------------------------------------------------
2275 ! SUBROUTINE CSNOW
2276 ! FUNCTION CSNOW
2277 ! ----------------------------------------------------------------------
2278 ! CALCULATE SNOW TERMAL CONDUCTIVITY
2279 ! ----------------------------------------------------------------------
2280 IMPLICIT NONE
2281 REAL, INTENT(IN) :: DSNOW
2282 REAL, INTENT(OUT):: SNCOND
2283 REAL :: C
2284 REAL, PARAMETER :: UNIT = 0.11631
2285
2286 ! ----------------------------------------------------------------------
2287 ! SNCOND IN UNITS OF CAL/(CM*HR*C), RETURNED IN W/(M*C)
2288 ! CSNOW IN UNITS OF CAL/(CM*HR*C), RETURNED IN W/(M*C)
2289 ! BASIC VERSION IS DYACHKOVA EQUATION (1960), FOR RANGE 0.1-0.4
2290 ! ----------------------------------------------------------------------
2291 C = 0.328*10** (2.25* DSNOW)
2292 ! CSNOW=UNIT*C
2293
2294 ! ----------------------------------------------------------------------
2295 ! DE VAUX EQUATION (1933), IN RANGE 0.1-0.6
2296 ! ----------------------------------------------------------------------
2297 ! SNCOND=0.0293*(1.+100.*DSNOW**2)
2298 ! CSNOW=0.0293*(1.+100.*DSNOW**2)
2299
2300 ! ----------------------------------------------------------------------
2301 ! E. ANDERSEN FROM FLERCHINGER
2302 ! ----------------------------------------------------------------------
2303 ! SNCOND=0.021+2.51*DSNOW**2
2304 ! CSNOW=0.021+2.51*DSNOW**2
2305
2306 ! SNCOND = UNIT * C
2307 ! double snow thermal conductivity
2308 SNCOND = 2.0 * UNIT * C
2309
2310 ! ----------------------------------------------------------------------
2311 END SUBROUTINE CSNOW
2312 ! ----------------------------------------------------------------------
2313
2314 SUBROUTINE DEVAP (EDIR,ETP1,SMC,ZSOIL,SHDFAC,SMCMAX,BEXP, &
2315 DKSAT,DWSAT,SMCDRY,SMCREF,SMCWLT,FXEXP)
2316
2317 ! ----------------------------------------------------------------------
2318 ! SUBROUTINE DEVAP
2319 ! FUNCTION DEVAP
2320 ! ----------------------------------------------------------------------
2321 ! CALCULATE DIRECT SOIL EVAPORATION
2322 ! ----------------------------------------------------------------------
2323 IMPLICIT NONE
2324 REAL, INTENT(IN) :: ETP1,SMC,BEXP,DKSAT,DWSAT,FXEXP, &
2325 SHDFAC,SMCDRY,SMCMAX,ZSOIL,SMCREF,SMCWLT
2326 REAL, INTENT(OUT):: EDIR
2327 REAL :: FX, SRATIO
2328
2329
2330 ! ----------------------------------------------------------------------
2331 ! DIRECT EVAP A FUNCTION OF RELATIVE SOIL MOISTURE AVAILABILITY, LINEAR
2332 ! WHEN FXEXP=1.
2333 ! ----------------------------------------------------------------------
2334 ! ----------------------------------------------------------------------
2335 ! FX > 1 REPRESENTS DEMAND CONTROL
2336 ! FX < 1 REPRESENTS FLUX CONTROL
2337 ! ----------------------------------------------------------------------
2338
2339 SRATIO = (SMC - SMCDRY) / (SMCMAX - SMCDRY)
2340 IF (SRATIO > 0.) THEN
2341 FX = SRATIO**FXEXP
2342 FX = MAX ( MIN ( FX, 1. ) ,0. )
2343 ELSE
2344 FX = 0.
2345 ENDIF
2346
2347 ! ----------------------------------------------------------------------
2348 ! ALLOW FOR THE DIRECT-EVAP-REDUCING EFFECT OF SHADE
2349 ! ----------------------------------------------------------------------
2350 EDIR = FX * ( 1.0- SHDFAC ) * ETP1
2351
2352 ! ----------------------------------------------------------------------
2353 END SUBROUTINE DEVAP
2354 ! ----------------------------------------------------------------------
2355
2356 SUBROUTINE EVAPO (ETA1,SMC,NSOIL,CMC,ETP1,DT,ZSOIL, &
2357 SH2O, &
2358 SMCMAX,BEXP,PC,SMCWLT,DKSAT,DWSAT, &
2359 SMCREF,SHDFAC,CMCMAX, &
2360 SMCDRY,CFACTR, &
2361 EDIR,EC,ET,ETT,SFCTMP,Q2,NROOT,RTDIS,FXEXP)
2362
2363 ! ----------------------------------------------------------------------
2364 ! SUBROUTINE EVAPO
2365 ! ----------------------------------------------------------------------
2366 ! CALCULATE SOIL MOISTURE FLUX. THE SOIL MOISTURE CONTENT (SMC - A PER
2367 ! UNIT VOLUME MEASUREMENT) IS A DEPENDENT VARIABLE THAT IS UPDATED WITH
2368 ! PROGNOSTIC EQNS. THE CANOPY MOISTURE CONTENT (CMC) IS ALSO UPDATED.
2369 ! FROZEN GROUND VERSION: NEW STATES ADDED: SH2O, AND FROZEN GROUND
2370 ! CORRECTION FACTOR, FRZFACT AND PARAMETER SLOPE.
2371 ! ----------------------------------------------------------------------
2372 IMPLICIT NONE
2373 INTEGER, INTENT(IN) :: NSOIL, NROOT
2374 INTEGER :: I,K
2375 REAL, INTENT(IN) :: BEXP, CFACTR,CMC,CMCMAX,DKSAT, &
2376 DT,DWSAT,ETP1,FXEXP,PC,Q2,SFCTMP, &
2377 SHDFAC,SMCDRY,SMCMAX,SMCREF,SMCWLT
2378 REAL, INTENT(OUT) :: EC,EDIR,ETA1,ETT
2379 REAL :: CMC2MS
2380 REAL,DIMENSION(1:NSOIL), INTENT(IN) :: RTDIS, SMC, SH2O, ZSOIL
2381 REAL,DIMENSION(1:NSOIL), INTENT(OUT) :: ET
2382
2383 ! ----------------------------------------------------------------------
2384 ! EXECUTABLE CODE BEGINS HERE IF THE POTENTIAL EVAPOTRANSPIRATION IS
2385 ! GREATER THAN ZERO.
2386 ! ----------------------------------------------------------------------
2387 EDIR = 0.
2388 EC = 0.
2389 ETT = 0.
2390 DO K = 1,NSOIL
2391 ET (K) = 0.
2392 END DO
2393
2394 ! ----------------------------------------------------------------------
2395 ! RETRIEVE DIRECT EVAPORATION FROM SOIL SURFACE. CALL THIS FUNCTION
2396 ! ONLY IF VEG COVER NOT COMPLETE.
2397 ! FROZEN GROUND VERSION: SH2O STATES REPLACE SMC STATES.
2398 ! ----------------------------------------------------------------------
2399 IF (ETP1 > 0.0) THEN
2400 IF (SHDFAC < 1.) THEN
2401 CALL DEVAP (EDIR,ETP1,SMC (1),ZSOIL (1),SHDFAC,SMCMAX, &
2402 BEXP,DKSAT,DWSAT,SMCDRY,SMCREF,SMCWLT,FXEXP)
2403 END IF
2404 ! ----------------------------------------------------------------------
2405 ! INITIALIZE PLANT TOTAL TRANSPIRATION, RETRIEVE PLANT TRANSPIRATION,
2406 ! AND ACCUMULATE IT FOR ALL SOIL LAYERS.
2407 ! ----------------------------------------------------------------------
2408
2409 IF (SHDFAC > 0.0) THEN
2410 CALL TRANSP (ET,NSOIL,ETP1,SH2O,CMC,ZSOIL,SHDFAC,SMCWLT, &
2411 CMCMAX,PC,CFACTR,SMCREF,SFCTMP,Q2,NROOT,RTDIS)
2412 DO K = 1,NSOIL
2413 ETT = ETT + ET ( K )
2414 END DO
2415 ! ----------------------------------------------------------------------
2416 ! CALCULATE CANOPY EVAPORATION.
2417 ! IF STATEMENTS TO AVOID TANGENT LINEAR PROBLEMS NEAR CMC=0.0.
2418 ! ----------------------------------------------------------------------
2419 IF (CMC > 0.0) THEN
2420 EC = SHDFAC * ( ( CMC / CMCMAX ) ** CFACTR ) * ETP1
2421 ELSE
2422 EC = 0.0
2423 END IF
2424 ! ----------------------------------------------------------------------
2425 ! EC SHOULD BE LIMITED BY THE TOTAL AMOUNT OF AVAILABLE WATER ON THE
2426 ! CANOPY. -F.CHEN, 18-OCT-1994
2427 ! ----------------------------------------------------------------------
2428 CMC2MS = CMC / DT
2429 EC = MIN ( CMC2MS, EC )
2430 END IF
2431 END IF
2432 ! ----------------------------------------------------------------------
2433 ! TOTAL UP EVAP AND TRANSP TYPES TO OBTAIN ACTUAL EVAPOTRANSP
2434 ! ----------------------------------------------------------------------
2435 ETA1 = EDIR + ETT + EC
2436
2437 ! ----------------------------------------------------------------------
2438 END SUBROUTINE EVAPO
2439 ! ----------------------------------------------------------------------
2440
2441 SUBROUTINE FRH2O (FREE,TKELV,SMC,SH2O,SMCMAX,BEXP,PSIS)
2442
2443 ! ----------------------------------------------------------------------
2444 ! SUBROUTINE FRH2O
2445 ! ----------------------------------------------------------------------
2446 ! CALCULATE AMOUNT OF SUPERCOOLED LIQUID SOIL WATER CONTENT IF
2447 ! TEMPERATURE IS BELOW 273.15K (T0). REQUIRES NEWTON-TYPE ITERATION TO
2448 ! SOLVE THE NONLINEAR IMPLICIT EQUATION GIVEN IN EQN 17 OF KOREN ET AL
2449 ! (1999, JGR, VOL 104(D16), 19569-19585).
2450 ! ----------------------------------------------------------------------
2451 ! NEW VERSION (JUNE 2001): MUCH FASTER AND MORE ACCURATE NEWTON
2452 ! ITERATION ACHIEVED BY FIRST TAKING LOG OF EQN CITED ABOVE -- LESS THAN
2453 ! 4 (TYPICALLY 1 OR 2) ITERATIONS ACHIEVES CONVERGENCE. ALSO, EXPLICIT
2454 ! 1-STEP SOLUTION OPTION FOR SPECIAL CASE OF PARAMETER CK=0, WHICH
2455 ! REDUCES THE ORIGINAL IMPLICIT EQUATION TO A SIMPLER EXPLICIT FORM,
2456 ! KNOWN AS THE "FLERCHINGER EQN". IMPROVED HANDLING OF SOLUTION IN THE
2457 ! LIMIT OF FREEZING POINT TEMPERATURE T0.
2458 ! ----------------------------------------------------------------------
2459 ! INPUT:
2460
2461 ! TKELV.........TEMPERATURE (Kelvin)
2462 ! SMC...........TOTAL SOIL MOISTURE CONTENT (VOLUMETRIC)
2463 ! SH2O..........LIQUID SOIL MOISTURE CONTENT (VOLUMETRIC)
2464 ! SMCMAX........SATURATION SOIL MOISTURE CONTENT (FROM REDPRM)
2465 ! B.............SOIL TYPE "B" PARAMETER (FROM REDPRM)
2466 ! PSIS..........SATURATED SOIL MATRIC POTENTIAL (FROM REDPRM)
2467
2468 ! OUTPUT:
2469 ! FRH2O.........SUPERCOOLED LIQUID WATER CONTENT
2470 ! FREE..........SUPERCOOLED LIQUID WATER CONTENT
2471 ! ----------------------------------------------------------------------
2472 IMPLICIT NONE
2473 REAL, INTENT(IN) :: BEXP,PSIS,SH2O,SMC,SMCMAX,TKELV
2474 REAL, INTENT(OUT) :: FREE
2475 REAL :: BX,DENOM,DF,DSWL,FK,SWL,SWLK
2476 INTEGER :: NLOG,KCOUNT
2477 ! PARAMETER(CK = 0.0)
2478 REAL, PARAMETER :: CK = 8.0, BLIM = 5.5, ERROR = 0.005, &
2479 HLICE = 3.335E5, GS = 9.81,DICE = 920.0, &
2480 DH2O = 1000.0, T0 = 273.15
2481
2482 ! ----------------------------------------------------------------------
2483 ! LIMITS ON PARAMETER B: B < 5.5 (use parameter BLIM)
2484 ! SIMULATIONS SHOWED IF B > 5.5 UNFROZEN WATER CONTENT IS
2485 ! NON-REALISTICALLY HIGH AT VERY LOW TEMPERATURES.
2486 ! ----------------------------------------------------------------------
2487 BX = BEXP
2488
2489 ! ----------------------------------------------------------------------
2490 ! INITIALIZING ITERATIONS COUNTER AND ITERATIVE SOLUTION FLAG.
2491 ! ----------------------------------------------------------------------
2492 IF (BEXP > BLIM) BX = BLIM
2493 NLOG = 0
2494
2495 ! ----------------------------------------------------------------------
2496 ! IF TEMPERATURE NOT SIGNIFICANTLY BELOW FREEZING (T0), SH2O = SMC
2497 ! ----------------------------------------------------------------------
2498 KCOUNT = 0
2499 ! FRH2O = SMC
2500 IF (TKELV > (T0- 1.E-3)) THEN
2501 FREE = SMC
2502 ELSE
2503
2504 ! ----------------------------------------------------------------------
2505 ! OPTION 1: ITERATED SOLUTION FOR NONZERO CK
2506 ! IN KOREN ET AL, JGR, 1999, EQN 17
2507 ! ----------------------------------------------------------------------
2508 ! INITIAL GUESS FOR SWL (frozen content)
2509 ! ----------------------------------------------------------------------
2510 IF (CK /= 0.0) THEN
2511 SWL = SMC - SH2O
2512 ! ----------------------------------------------------------------------
2513 ! KEEP WITHIN BOUNDS.
2514 ! ----------------------------------------------------------------------
2515 IF (SWL > (SMC -0.02)) SWL = SMC -0.02
2516
2517 ! ----------------------------------------------------------------------
2518 ! START OF ITERATIONS
2519 ! ----------------------------------------------------------------------
2520 IF (SWL < 0.) SWL = 0.
2521 1001 Continue
2522 IF (.NOT.( (NLOG < 10) .AND. (KCOUNT == 0))) goto 1002
2523 NLOG = NLOG +1
2524 DF = ALOG ( ( PSIS * GS / HLICE ) * ( ( 1. + CK * SWL )**2.) * &
2525 ( SMCMAX / (SMC - SWL) )** BX) - ALOG ( - ( &
2526 TKELV - T0)/ TKELV)
2527 DENOM = 2. * CK / ( 1. + CK * SWL ) + BX / ( SMC - SWL )
2528 SWLK = SWL - DF / DENOM
2529 ! ----------------------------------------------------------------------
2530 ! BOUNDS USEFUL FOR MATHEMATICAL SOLUTION.
2531 ! ----------------------------------------------------------------------
2532 IF (SWLK > (SMC -0.02)) SWLK = SMC - 0.02
2533 IF (SWLK < 0.) SWLK = 0.
2534
2535 ! ----------------------------------------------------------------------
2536 ! MATHEMATICAL SOLUTION BOUNDS APPLIED.
2537 ! ----------------------------------------------------------------------
2538 DSWL = ABS (SWLK - SWL)
2539
2540 ! ----------------------------------------------------------------------
2541 ! IF MORE THAN 10 ITERATIONS, USE EXPLICIT METHOD (CK=0 APPROX.)
2542 ! WHEN DSWL LESS OR EQ. ERROR, NO MORE ITERATIONS REQUIRED.
2543 ! ----------------------------------------------------------------------
2544 SWL = SWLK
2545 IF ( DSWL <= ERROR ) THEN
2546 KCOUNT = KCOUNT +1
2547 END IF
2548 ! ----------------------------------------------------------------------
2549 ! END OF ITERATIONS
2550 ! ----------------------------------------------------------------------
2551 ! BOUNDS APPLIED WITHIN DO-BLOCK ARE VALID FOR PHYSICAL SOLUTION.
2552 ! ----------------------------------------------------------------------
2553 ! FRH2O = SMC - SWL
2554 goto 1001
2555 1002 continue
2556 FREE = SMC - SWL
2557 END IF
2558 ! ----------------------------------------------------------------------
2559 ! END OPTION 1
2560 ! ----------------------------------------------------------------------
2561 ! ----------------------------------------------------------------------
2562 ! OPTION 2: EXPLICIT SOLUTION FOR FLERCHINGER EQ. i.e. CK=0
2563 ! IN KOREN ET AL., JGR, 1999, EQN 17
2564 ! APPLY PHYSICAL BOUNDS TO FLERCHINGER SOLUTION
2565 ! ----------------------------------------------------------------------
2566 IF (KCOUNT == 0) THEN
2567 PRINT *,'Flerchinger USEd in NEW version. Iterations=',NLOG
2568 FK = ( ( (HLICE / (GS * ( - PSIS)))* &
2569 ( (TKELV - T0)/ TKELV))** ( -1/ BX))* SMCMAX
2570 ! FRH2O = MIN (FK, SMC)
2571 IF (FK < 0.02) FK = 0.02
2572 FREE = MIN (FK, SMC)
2573 ! ----------------------------------------------------------------------
2574 ! END OPTION 2
2575 ! ----------------------------------------------------------------------
2576 END IF
2577 END IF
2578 ! ----------------------------------------------------------------------
2579 END SUBROUTINE FRH2O
2580 ! ----------------------------------------------------------------------
2581
2582 SUBROUTINE HRT (RHSTS,STC,SMC,SMCMAX,NSOIL,ZSOIL,YY,ZZ1, &
2583 TBOT,ZBOT,PSISAT,SH2O,DT,BEXP, &
2584 F1,DF1,QUARTZ,CSOIL,AI,BI,CI,VEGTYP)
2585
2586 ! ----------------------------------------------------------------------
2587 ! SUBROUTINE HRT
2588 ! ----------------------------------------------------------------------
2589 ! CALCULATE THE RIGHT HAND SIDE OF THE TIME TENDENCY TERM OF THE SOIL
2590 ! THERMAL DIFFUSION EQUATION. ALSO TO COMPUTE ( PREPARE ) THE MATRIX
2591 ! COEFFICIENTS FOR THE TRI-DIAGONAL MATRIX OF THE IMPLICIT TIME SCHEME.
2592 ! ----------------------------------------------------------------------
2593 IMPLICIT NONE
2594 LOGICAL :: ITAVG
2595 INTEGER, INTENT(IN) :: NSOIL, VEGTYP
2596 INTEGER :: I, K
2597
2598 REAL, INTENT(IN) :: BEXP, CSOIL, DF1, DT,F1,PSISAT,QUARTZ, &
2599 SMCMAX ,TBOT,YY,ZZ1, ZBOT
2600 REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SMC,STC,ZSOIL
2601 REAL, DIMENSION(1:NSOIL), INTENT(INOUT):: SH2O
2602 REAL, DIMENSION(1:NSOIL), INTENT(OUT) :: RHSTS
2603 REAL, DIMENSION(1:NSOIL), INTENT(OUT) :: AI, BI,CI
2604 REAL :: DDZ, DDZ2, DENOM, DF1N, DF1K, DTSDZ, &
2605 DTSDZ2,HCPCT,QTOT,SSOIL,SICE,TAVG,TBK, &
2606 TBK1,TSNSR,TSURF,CSOIL_LOC
2607 REAL, PARAMETER :: T0 = 273.15, CAIR = 1004.0, CICE = 2.106E6,&
2608 CH2O = 4.2E6
2609
2610
2611 !urban
2612 IF(VEGTYP==1) then
2613 CSOIL_LOC=3.0E6
2614 ELSE
2615 CSOIL_LOC=CSOIL
2616 ENDIF
2617
2618 ! ----------------------------------------------------------------------
2619 ! INITIALIZE LOGICAL FOR SOIL LAYER TEMPERATURE AVERAGING.
2620 ! ----------------------------------------------------------------------
2621 ITAVG = .TRUE.
2622 ! ----------------------------------------------------------------------
2623 ! BEGIN SECTION FOR TOP SOIL LAYER
2624 ! ----------------------------------------------------------------------
2625 ! CALC THE HEAT CAPACITY OF THE TOP SOIL LAYER
2626 ! ----------------------------------------------------------------------
2627 HCPCT = SH2O (1)* CH2O + (1.0- SMCMAX)* CSOIL_LOC + (SMCMAX - SMC (1))&
2628 * CAIR &
2629 + ( SMC (1) - SH2O (1) )* CICE
2630
2631 ! ----------------------------------------------------------------------
2632 ! CALC THE MATRIX COEFFICIENTS AI, BI, AND CI FOR THE TOP LAYER
2633 ! ----------------------------------------------------------------------
2634 DDZ = 1.0 / ( -0.5 * ZSOIL (2) )
2635 AI (1) = 0.0
2636 CI (1) = (DF1 * DDZ) / (ZSOIL (1) * HCPCT)
2637
2638 ! ----------------------------------------------------------------------
2639 ! CALCULATE THE VERTICAL SOIL TEMP GRADIENT BTWN THE 1ST AND 2ND SOIL
2640 ! LAYERS. THEN CALCULATE THE SUBSURFACE HEAT FLUX. USE THE TEMP
2641 ! GRADIENT AND SUBSFC HEAT FLUX TO CALC "RIGHT-HAND SIDE TENDENCY
2642 ! TERMS", OR "RHSTS", FOR TOP SOIL LAYER.
2643 ! ----------------------------------------------------------------------
2644 BI (1) = - CI (1) + DF1 / (0.5 * ZSOIL (1) * ZSOIL (1)* HCPCT * &
2645 ZZ1)
2646 DTSDZ = (STC (1) - STC (2)) / ( -0.5 * ZSOIL (2))
2647 SSOIL = DF1 * (STC (1) - YY) / (0.5 * ZSOIL (1) * ZZ1)
2648 ! RHSTS(1) = (DF1 * DTSDZ - SSOIL) / (ZSOIL(1) * HCPCT)
2649 DENOM = (ZSOIL (1) * HCPCT)
2650
2651 ! ----------------------------------------------------------------------
2652 ! NEXT CAPTURE THE VERTICAL DIFFERENCE OF THE HEAT FLUX AT TOP AND
2653 ! BOTTOM OF FIRST SOIL LAYER FOR USE IN HEAT FLUX CONSTRAINT APPLIED TO
2654 ! POTENTIAL SOIL FREEZING/THAWING IN ROUTINE SNKSRC.
2655 ! ----------------------------------------------------------------------
2656 ! QTOT = SSOIL - DF1*DTSDZ
2657 RHSTS (1) = (DF1 * DTSDZ - SSOIL) / DENOM
2658
2659 ! ----------------------------------------------------------------------
2660 ! CALCULATE FROZEN WATER CONTENT IN 1ST SOIL LAYER.
2661 ! ----------------------------------------------------------------------
2662 QTOT = -1.0* RHSTS (1)* DENOM
2663
2664 ! ----------------------------------------------------------------------
2665 ! IF TEMPERATURE AVERAGING INVOKED (ITAVG=TRUE; ELSE SKIP):
2666 ! SET TEMP "TSURF" AT TOP OF SOIL COLUMN (FOR USE IN FREEZING SOIL
2667 ! PHYSICS LATER IN FUNCTION SUBROUTINE SNKSRC). IF SNOWPACK CONTENT IS
2668 ! ZERO, THEN TSURF EXPRESSION BELOW GIVES TSURF = SKIN TEMP. IF
2669 ! SNOWPACK IS NONZERO (HENCE ARGUMENT ZZ1=1), THEN TSURF EXPRESSION
2670 ! BELOW YIELDS SOIL COLUMN TOP TEMPERATURE UNDER SNOWPACK. THEN
2671 ! CALCULATE TEMPERATURE AT BOTTOM INTERFACE OF 1ST SOIL LAYER FOR USE
2672 ! LATER IN FUNCTION SUBROUTINE SNKSRC
2673 ! ----------------------------------------------------------------------
2674 SICE = SMC (1) - SH2O (1)
2675 IF (ITAVG) THEN
2676 TSURF = (YY + (ZZ1-1) * STC (1)) / ZZ1
2677 ! ----------------------------------------------------------------------
2678 ! IF FROZEN WATER PRESENT OR ANY OF LAYER-1 MID-POINT OR BOUNDING
2679 ! INTERFACE TEMPERATURES BELOW FREEZING, THEN CALL SNKSRC TO
2680 ! COMPUTE HEAT SOURCE/SINK (AND CHANGE IN FROZEN WATER CONTENT)
2681 ! DUE TO POSSIBLE SOIL WATER PHASE CHANGE
2682 ! ----------------------------------------------------------------------
2683 CALL TBND (STC (1),STC (2),ZSOIL,ZBOT,1,NSOIL,TBK)
2684 IF ( (SICE > 0.) .OR. (STC (1) < T0) .OR. &
2685 (TSURF < T0) .OR. (TBK < T0) ) THEN
2686 ! TSNSR = SNKSRC (TAVG,SMC(1),SH2O(1),
2687 CALL TMPAVG (TAVG,TSURF,STC (1),TBK,ZSOIL,NSOIL,1)
2688 CALL SNKSRC (TSNSR,TAVG,SMC (1),SH2O (1), &
2689 ZSOIL,NSOIL,SMCMAX,PSISAT,BEXP,DT,1,QTOT)
2690 ! RHSTS(1) = RHSTS(1) - TSNSR / ( ZSOIL(1) * HCPCT )
2691 RHSTS (1) = RHSTS (1) - TSNSR / DENOM
2692 END IF
2693 ELSE
2694 ! TSNSR = SNKSRC (STC(1),SMC(1),SH2O(1),
2695 IF ( (SICE > 0.) .OR. (STC (1) < T0) ) THEN
2696 CALL SNKSRC (TSNSR,STC (1),SMC (1),SH2O (1), &
2697 ZSOIL,NSOIL,SMCMAX,PSISAT,BEXP,DT,1,QTOT)
2698 ! RHSTS(1) = RHSTS(1) - TSNSR / ( ZSOIL(1) * HCPCT )
2699 RHSTS (1) = RHSTS (1) - TSNSR / DENOM
2700 END IF
2701 ! ----------------------------------------------------------------------
2702 ! THIS ENDS SECTION FOR TOP SOIL LAYER.
2703 ! ----------------------------------------------------------------------
2704 END IF
2705
2706 ! INITIALIZE DDZ2
2707 ! ----------------------------------------------------------------------
2708
2709 DDZ2 = 0.0
2710 DF1K = DF1
2711
2712 ! ----------------------------------------------------------------------
2713 ! LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABOVE PROCESS
2714 ! (EXCEPT SUBSFC OR "GROUND" HEAT FLUX NOT REPEATED IN LOWER LAYERS)
2715 ! ----------------------------------------------------------------------
2716 ! CALCULATE HEAT CAPACITY FOR THIS SOIL LAYER.
2717 ! ----------------------------------------------------------------------
2718 DO K = 2,NSOIL
2719 HCPCT = SH2O (K)* CH2O + (1.0- SMCMAX)* CSOIL_LOC + (SMCMAX - SMC ( &
2720 K))* CAIR + ( SMC (K) - SH2O (K) )* CICE
2721 ! ----------------------------------------------------------------------
2722 ! THIS SECTION FOR LAYER 2 OR GREATER, BUT NOT LAST LAYER.
2723 ! ----------------------------------------------------------------------
2724 ! CALCULATE THERMAL DIFFUSIVITY FOR THIS LAYER.
2725 ! ----------------------------------------------------------------------
2726 IF (K /= NSOIL) THEN
2727
2728 ! ----------------------------------------------------------------------
2729 ! CALC THE VERTICAL SOIL TEMP GRADIENT THRU THIS LAYER
2730 ! ----------------------------------------------------------------------
2731 CALL TDFCND (DF1N,SMC (K),QUARTZ,SMCMAX,SH2O (K))
2732
2733 !urban
2734 IF(VEGTYP==1) DF1N = 3.24
2735
2736 DENOM = 0.5 * ( ZSOIL (K -1) - ZSOIL (K +1) )
2737
2738 ! ----------------------------------------------------------------------
2739 ! CALC THE MATRIX COEF, CI, AFTER CALC'NG ITS PARTIAL PRODUCT
2740 ! ----------------------------------------------------------------------
2741 DTSDZ2 = ( STC (K) - STC (K +1) ) / DENOM
2742 DDZ2 = 2. / (ZSOIL (K -1) - ZSOIL (K +1))
2743
2744 ! ----------------------------------------------------------------------
2745 ! IF TEMPERATURE AVERAGING INVOKED (ITAVG=TRUE; ELSE SKIP): CALCULATE
2746 ! TEMP AT BOTTOM OF LAYER.
2747 ! ----------------------------------------------------------------------
2748 CI (K) = - DF1N * DDZ2 / ( (ZSOIL (K -1) - ZSOIL (K)) * &
2749 HCPCT)
2750 IF (ITAVG) THEN
2751 CALL TBND (STC (K),STC (K +1),ZSOIL,ZBOT,K,NSOIL,TBK1)
2752 END IF
2753
2754 ELSE
2755 ! ----------------------------------------------------------------------
2756 ! SPECIAL CASE OF BOTTOM SOIL LAYER: CALCULATE THERMAL DIFFUSIVITY FOR
2757 ! BOTTOM LAYER.
2758 ! ----------------------------------------------------------------------
2759
2760 ! ----------------------------------------------------------------------
2761 ! CALC THE VERTICAL SOIL TEMP GRADIENT THRU BOTTOM LAYER.
2762 ! ----------------------------------------------------------------------
2763 CALL TDFCND (DF1N,SMC (K),QUARTZ,SMCMAX,SH2O (K))
2764
2765
2766 !urban
2767 IF(VEGTYP==1) DF1N = 3.24
2768
2769 DENOM = .5 * (ZSOIL (K -1) + ZSOIL (K)) - ZBOT
2770
2771 ! ----------------------------------------------------------------------
2772 ! SET MATRIX COEF, CI TO ZERO IF BOTTOM LAYER.
2773 ! ----------------------------------------------------------------------
2774 DTSDZ2 = (STC (K) - TBOT) / DENOM
2775
2776 ! ----------------------------------------------------------------------
2777 ! IF TEMPERATURE AVERAGING INVOKED (ITAVG=TRUE; ELSE SKIP): CALCULATE
2778 ! TEMP AT BOTTOM OF LAST LAYER.
2779 ! ----------------------------------------------------------------------
2780 CI (K) = 0.
2781 IF (ITAVG) THEN
2782 CALL TBND (STC (K),TBOT,ZSOIL,ZBOT,K,NSOIL,TBK1)
2783 END IF
2784 ! ----------------------------------------------------------------------
2785 ! THIS ENDS SPECIAL LOOP FOR BOTTOM LAYER.
2786 END IF
2787 ! ----------------------------------------------------------------------
2788 ! CALCULATE RHSTS FOR THIS LAYER AFTER CALC'NG A PARTIAL PRODUCT.
2789 ! ----------------------------------------------------------------------
2790 DENOM = ( ZSOIL (K) - ZSOIL (K -1) ) * HCPCT
2791 RHSTS (K) = ( DF1N * DTSDZ2- DF1K * DTSDZ ) / DENOM
2792 QTOT = -1.0* DENOM * RHSTS (K)
2793
2794 SICE = SMC (K) - SH2O (K)
2795 IF (ITAVG) THEN
2796 CALL TMPAVG (TAVG,TBK,STC (K),TBK1,ZSOIL,NSOIL,K)
2797 ! TSNSR = SNKSRC(TAVG,SMC(K),SH2O(K),ZSOIL,NSOIL,
2798 IF ( (SICE > 0.) .OR. (STC (K) < T0) .OR. &
2799 (TBK .lt. T0) .OR. (TBK1 .lt. T0) ) THEN
2800 CALL SNKSRC (TSNSR,TAVG,SMC (K),SH2O (K),ZSOIL,NSOIL, &
2801 SMCMAX,PSISAT,BEXP,DT,K,QTOT)
2802 RHSTS (K) = RHSTS (K) - TSNSR / DENOM
2803 END IF
2804 ELSE
2805 ! TSNSR = SNKSRC(STC(K),SMC(K),SH2O(K),ZSOIL,NSOIL,
2806 IF ( (SICE > 0.) .OR. (STC (K) < T0) ) THEN
2807 CALL SNKSRC (TSNSR,STC (K),SMC (K),SH2O (K),ZSOIL,NSOIL, &
2808 SMCMAX,PSISAT,BEXP,DT,K,QTOT)
2809 RHSTS (K) = RHSTS (K) - TSNSR / DENOM
2810 END IF
2811 END IF
2812
2813 ! ----------------------------------------------------------------------
2814 ! CALC MATRIX COEFS, AI, AND BI FOR THIS LAYER.
2815 ! ----------------------------------------------------------------------
2816 AI (K) = - DF1K * DDZ / ( (ZSOIL (K -1) - ZSOIL (K)) * HCPCT)
2817
2818 ! ----------------------------------------------------------------------
2819 ! RESET VALUES OF DF1, DTSDZ, DDZ, AND TBK FOR LOOP TO NEXT SOIL LAYER.
2820 ! ----------------------------------------------------------------------
2821 BI (K) = - (AI (K) + CI (K))
2822 TBK = TBK1
2823 DF1K = DF1N
2824 DTSDZ = DTSDZ2
2825 DDZ = DDZ2
2826 END DO
2827 ! ----------------------------------------------------------------------
2828 END SUBROUTINE HRT
2829 ! ----------------------------------------------------------------------
2830
2831 SUBROUTINE HRTICE (RHSTS,STC,NSOIL,ZSOIL,YY,ZZ1,DF1,AI,BI,CI)
2832
2833 ! ----------------------------------------------------------------------
2834 ! SUBROUTINE HRTICE
2835 ! ----------------------------------------------------------------------
2836 ! CALCULATE THE RIGHT HAND SIDE OF THE TIME TENDENCY TERM OF THE SOIL
2837 ! THERMAL DIFFUSION EQUATION IN THE CASE OF SEA-ICE PACK. ALSO TO
2838 ! COMPUTE (PREPARE) THE MATRIX COEFFICIENTS FOR THE TRI-DIAGONAL MATRIX
2839 ! OF THE IMPLICIT TIME SCHEME.
2840 ! ----------------------------------------------------------------------
2841 IMPLICIT NONE
2842
2843
2844 INTEGER, INTENT(IN) :: NSOIL
2845 INTEGER :: K
2846
2847 REAL, INTENT(IN) :: DF1,YY,ZZ1
2848 REAL, DIMENSION(1:NSOIL), INTENT(OUT):: AI, BI,CI
2849 REAL, DIMENSION(1:NSOIL), INTENT(IN) :: STC, ZSOIL
2850 REAL, DIMENSION(1:NSOIL), INTENT(OUT):: RHSTS
2851 REAL :: DDZ,DDZ2,DENOM,DTSDZ,DTSDZ2,SSOIL, &
2852 ZBOT,TBOT
2853 REAL, PARAMETER :: HCPCT = 1.72396E+6
2854
2855 ! ----------------------------------------------------------------------
2856 ! SET A NOMINAL UNIVERSAL VALUE OF THE SEA-ICE SPECIFIC HEAT CAPACITY,
2857 ! HCPCT = 1880.0*917.0.
2858 ! ----------------------------------------------------------------------
2859 DATA TBOT /271.16/
2860
2861 ! ----------------------------------------------------------------------
2862 ! THE INPUT ARGUMENT DF1 IS A UNIVERSALLY CONSTANT VALUE OF SEA-ICE
2863 ! THERMAL DIFFUSIVITY, SET IN ROUTINE SNOPAC AS DF1 = 2.2.
2864 ! ----------------------------------------------------------------------
2865 ! SET ICE PACK DEPTH. USE TBOT AS ICE PACK LOWER BOUNDARY TEMPERATURE
2866 ! (THAT OF UNFROZEN SEA WATER AT BOTTOM OF SEA ICE PACK). ASSUME ICE
2867 ! PACK IS OF N=NSOIL LAYERS SPANNING A UNIFORM CONSTANT ICE PACK
2868 ! THICKNESS AS DEFINED BY ZSOIL(NSOIL) IN ROUTINE SFLX.
2869 ! ----------------------------------------------------------------------
2870 ! ----------------------------------------------------------------------
2871 ! CALC THE MATRIX COEFFICIENTS AI, BI, AND CI FOR THE TOP LAYER
2872 ! ----------------------------------------------------------------------
2873 ZBOT = ZSOIL (NSOIL)
2874 DDZ = 1.0 / ( -0.5 * ZSOIL (2) )
2875 AI (1) = 0.0
2876 CI (1) = (DF1 * DDZ) / (ZSOIL (1) * HCPCT)
2877
2878 ! ----------------------------------------------------------------------
2879 ! CALC THE VERTICAL SOIL TEMP GRADIENT BTWN THE TOP AND 2ND SOIL LAYERS.
2880 ! RECALC/ADJUST THE SOIL HEAT FLUX. USE THE GRADIENT AND FLUX TO CALC
2881 ! RHSTS FOR THE TOP SOIL LAYER.
2882 ! ----------------------------------------------------------------------
2883 BI (1) = - CI (1) + DF1/ (0.5 * ZSOIL (1) * ZSOIL (1) * HCPCT * &
2884 ZZ1)
2885 DTSDZ = ( STC (1) - STC (2) ) / ( -0.5 * ZSOIL (2) )
2886 SSOIL = DF1 * ( STC (1) - YY ) / ( 0.5 * ZSOIL (1) * ZZ1 )
2887
2888 ! ----------------------------------------------------------------------
2889 ! INITIALIZE DDZ2
2890 ! ----------------------------------------------------------------------
2891 RHSTS (1) = ( DF1 * DTSDZ - SSOIL ) / ( ZSOIL (1) * HCPCT )
2892
2893 ! ----------------------------------------------------------------------
2894 ! LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABOVE PROCESS
2895 ! ----------------------------------------------------------------------
2896 DDZ2 = 0.0
2897 DO K = 2,NSOIL
2898
2899 ! ----------------------------------------------------------------------
2900 ! CALC THE VERTICAL SOIL TEMP GRADIENT THRU THIS LAYER.
2901 ! ----------------------------------------------------------------------
2902 IF (K /= NSOIL) THEN
2903 DENOM = 0.5 * ( ZSOIL (K -1) - ZSOIL (K +1) )
2904
2905 ! ----------------------------------------------------------------------
2906 ! CALC THE MATRIX COEF, CI, AFTER CALC'NG ITS PARTIAL PRODUCT.
2907 ! ----------------------------------------------------------------------
2908 DTSDZ2 = ( STC (K) - STC (K +1) ) / DENOM
2909 DDZ2 = 2. / (ZSOIL (K -1) - ZSOIL (K +1))
2910 CI (K) = - DF1 * DDZ2 / ( (ZSOIL (K -1) - ZSOIL (K))*HCPCT)
2911
2912 ! ----------------------------------------------------------------------
2913 ! CALC THE VERTICAL SOIL TEMP GRADIENT THRU THE LOWEST LAYER.
2914 ! ----------------------------------------------------------------------
2915 ELSE
2916
2917 ! ----------------------------------------------------------------------
2918 ! SET MATRIX COEF, CI TO ZERO.
2919 ! ----------------------------------------------------------------------
2920 DTSDZ2 = (STC (K) - TBOT)/ (.5 * (ZSOIL (K -1) + ZSOIL (K)) &
2921 - ZBOT)
2922 CI (K) = 0.
2923 ! ----------------------------------------------------------------------
2924 ! CALC RHSTS FOR THIS LAYER AFTER CALC'NG A PARTIAL PRODUCT.
2925 ! ----------------------------------------------------------------------
2926 END IF
2927 DENOM = ( ZSOIL (K) - ZSOIL (K -1) ) * HCPCT
2928
2929 ! ----------------------------------------------------------------------
2930 ! CALC MATRIX COEFS, AI, AND BI FOR THIS LAYER.
2931 ! ----------------------------------------------------------------------
2932 RHSTS (K) = ( DF1 * DTSDZ2- DF1 * DTSDZ ) / DENOM
2933 AI (K) = - DF1 * DDZ / ( (ZSOIL (K -1) - ZSOIL (K)) * HCPCT)
2934
2935 ! ----------------------------------------------------------------------
2936 ! RESET VALUES OF DTSDZ AND DDZ FOR LOOP TO NEXT SOIL LYR.
2937 ! ----------------------------------------------------------------------
2938 BI (K) = - (AI (K) + CI (K))
2939 DTSDZ = DTSDZ2
2940 DDZ = DDZ2
2941 END DO
2942 ! ----------------------------------------------------------------------
2943 END SUBROUTINE HRTICE
2944 ! ----------------------------------------------------------------------
2945
2946 SUBROUTINE HSTEP (STCOUT,STCIN,RHSTS,DT,NSOIL,AI,BI,CI)
2947
2948 ! ----------------------------------------------------------------------
2949 ! SUBROUTINE HSTEP
2950 ! ----------------------------------------------------------------------
2951 ! CALCULATE/UPDATE THE SOIL TEMPERATURE FIELD.
2952 ! ----------------------------------------------------------------------
2953 IMPLICIT NONE
2954 INTEGER, INTENT(IN) :: NSOIL
2955 INTEGER :: K
2956
2957 REAL, DIMENSION(1:NSOIL), INTENT(IN):: STCIN
2958 REAL, DIMENSION(1:NSOIL), INTENT(OUT):: STCOUT
2959 REAL, DIMENSION(1:NSOIL), INTENT(INOUT):: RHSTS
2960 REAL, DIMENSION(1:NSOIL), INTENT(INOUT):: AI,BI,CI
2961 REAL, DIMENSION(1:NSOIL) :: RHSTSin
2962 REAL, DIMENSION(1:NSOIL) :: CIin
2963 REAL :: DT
2964
2965 ! ----------------------------------------------------------------------
2966 ! CREATE FINITE DIFFERENCE VALUES FOR USE IN ROSR12 ROUTINE
2967 ! ----------------------------------------------------------------------
2968 DO K = 1,NSOIL
2969 RHSTS (K) = RHSTS (K) * DT
2970 AI (K) = AI (K) * DT
2971 BI (K) = 1. + BI (K) * DT
2972 CI (K) = CI (K) * DT
2973 END DO
2974 ! ----------------------------------------------------------------------
2975 ! COPY VALUES FOR INPUT VARIABLES BEFORE CALL TO ROSR12
2976 ! ----------------------------------------------------------------------
2977 DO K = 1,NSOIL
2978 RHSTSin (K) = RHSTS (K)
2979 END DO
2980 DO K = 1,NSOIL
2981 CIin (K) = CI (K)
2982 END DO
2983 ! ----------------------------------------------------------------------
2984 ! SOLVE THE TRI-DIAGONAL MATRIX EQUATION
2985 ! ----------------------------------------------------------------------
2986 CALL ROSR12 (CI,AI,BI,CIin,RHSTSin,RHSTS,NSOIL)
2987 ! ----------------------------------------------------------------------
2988 ! CALC/UPDATE THE SOIL TEMPS USING MATRIX SOLUTION
2989 ! ----------------------------------------------------------------------
2990 DO K = 1,NSOIL
2991 STCOUT (K) = STCIN (K) + CI (K)
2992 END DO
2993 ! ----------------------------------------------------------------------
2994 END SUBROUTINE HSTEP
2995 ! ----------------------------------------------------------------------
2996
2997 SUBROUTINE NOPAC (ETP,ETA,PRCP,SMC,SMCMAX,SMCWLT, &
2998 SMCREF,SMCDRY,CMC,CMCMAX,NSOIL,DT,SHDFAC, &
2999 SBETA,Q2,T1,SFCTMP,T24,TH2,FDOWN,F1,EMISSI, &
3000 SSOIL, &
3001 STC,EPSCA,BEXP,PC,RCH,RR,CFACTR, &
3002 SH2O,SLOPE,KDT,FRZFACT,PSISAT,ZSOIL, &
3003 DKSAT,DWSAT,TBOT,ZBOT,RUNOFF1,RUNOFF2, &
3004 RUNOFF3,EDIR,EC,ET,ETT,NROOT,ICE,RTDIS, &
3005 QUARTZ,FXEXP,CSOIL, &
3006 BETA,DRIP,DEW,FLX1,FLX3,VEGTYP)
3007
3008 ! ----------------------------------------------------------------------
3009 ! SUBROUTINE NOPAC
3010 ! ----------------------------------------------------------------------
3011 ! CALCULATE SOIL MOISTURE AND HEAT FLUX VALUES AND UPDATE SOIL MOISTURE
3012 ! CONTENT AND SOIL HEAT CONTENT VALUES FOR THE CASE WHEN NO SNOW PACK IS
3013 ! PRESENT.
3014 ! ----------------------------------------------------------------------
3015 IMPLICIT NONE
3016
3017 INTEGER, INTENT(IN) :: ICE, NROOT,NSOIL,VEGTYP
3018 INTEGER :: K
3019
3020 REAL, INTENT(IN) :: BEXP,CFACTR, CMCMAX,CSOIL,DKSAT,DT,DWSAT, &
3021 EPSCA,ETP,FDOWN,F1,FXEXP,FRZFACT,KDT,PC, &
3022 PRCP,PSISAT,Q2,QUARTZ,RCH,RR,SBETA,SFCTMP,&
3023 SHDFAC,SLOPE,SMCDRY,SMCMAX,SMCREF,SMCWLT, &
3024 T24,TBOT,TH2,ZBOT,EMISSI
3025 REAL, INTENT(INOUT) :: CMC,BETA,T1
3026 REAL, INTENT(OUT) :: DEW,DRIP,EC,EDIR,ETA,ETT,FLX1,FLX3, &
3027 RUNOFF1,RUNOFF2,RUNOFF3,SSOIL
3028 REAL, DIMENSION(1:NSOIL),INTENT(IN) :: RTDIS,ZSOIL
3029 REAL, DIMENSION(1:NSOIL),INTENT(OUT) :: ET
3030 REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: SMC,SH2O,STC
3031 REAL, DIMENSION(1:NSOIL) :: ET1
3032 REAL :: EC1,EDIR1,ETT1,DF1,ETA1,ETP1,PRCP1,YY, &
3033 YYNUM,ZZ1
3034
3035 ! ----------------------------------------------------------------------
3036 ! EXECUTABLE CODE BEGINS HERE:
3037 ! CONVERT ETP Fnd PRCP FROM KG M-2 S-1 TO M S-1 AND INITIALIZE DEW.
3038 ! ----------------------------------------------------------------------
3039 PRCP1 = PRCP * 0.001
3040 ETP1 = ETP * 0.001
3041 DEW = 0.0
3042 ! ----------------------------------------------------------------------
3043 ! INITIALIZE EVAP TERMS.
3044 ! ----------------------------------------------------------------------
3045 EDIR = 0.
3046 EDIR1 = 0.
3047 EC1 = 0.
3048 EC = 0.
3049 DO K = 1,NSOIL
3050 ET(K) = 0.
3051 ET1(K) = 0.
3052 END DO
3053 ETT = 0.
3054 ETT1 = 0.
3055
3056 IF (ETP > 0.0) THEN
3057 CALL EVAPO (ETA1,SMC,NSOIL,CMC,ETP1,DT,ZSOIL, &
3058 SH2O, &
3059 SMCMAX,BEXP,PC,SMCWLT,DKSAT,DWSAT, &
3060 SMCREF,SHDFAC,CMCMAX, &
3061 SMCDRY,CFACTR, &
3062 EDIR1,EC1,ET1,ETT1,SFCTMP,Q2,NROOT,RTDIS,FXEXP)
3063 CALL SMFLX (SMC,NSOIL,CMC,DT,PRCP1,ZSOIL, &
3064 SH2O,SLOPE,KDT,FRZFACT, &
3065 SMCMAX,BEXP,SMCWLT,DKSAT,DWSAT, &
3066 SHDFAC,CMCMAX, &
3067 RUNOFF1,RUNOFF2,RUNOFF3, &
3068 EDIR1,EC1,ET1, &
3069 DRIP)
3070
3071 ! ----------------------------------------------------------------------
3072 ! CONVERT MODELED EVAPOTRANSPIRATION FROM M S-1 TO KG M-2 S-1.
3073 ! ----------------------------------------------------------------------
3074
3075 ETA = ETA1 * 1000.0
3076
3077 ! ----------------------------------------------------------------------
3078 ! IF ETP < 0, ASSUME DEW FORMS (TRANSFORM ETP1 INTO DEW AND REINITIALIZE
3079 ! ETP1 TO ZERO).
3080 ! ----------------------------------------------------------------------
3081 ELSE
3082 DEW = - ETP1
3083
3084 ! ----------------------------------------------------------------------
3085 ! CONVERT PRCP FROM 'KG M-2 S-1' TO 'M S-1' AND ADD DEW AMOUNT.
3086 ! ----------------------------------------------------------------------
3087
3088 PRCP1 = PRCP1+ DEW
3089 CALL SMFLX (SMC,NSOIL,CMC,DT,PRCP1,ZSOIL, &
3090 SH2O,SLOPE,KDT,FRZFACT, &
3091 SMCMAX,BEXP,SMCWLT,DKSAT,DWSAT, &
3092 SHDFAC,CMCMAX, &
3093 RUNOFF1,RUNOFF2,RUNOFF3, &
3094 EDIR1,EC1,ET1, &
3095 DRIP)
3096
3097 ! ----------------------------------------------------------------------
3098 ! CONVERT MODELED EVAPOTRANSPIRATION FROM 'M S-1' TO 'KG M-2 S-1'.
3099 ! ----------------------------------------------------------------------
3100 ! ETA = ETA1 * 1000.0
3101 END IF
3102
3103 ! ----------------------------------------------------------------------
3104 ! BASED ON ETP AND E VALUES, DETERMINE BETA
3105 ! ----------------------------------------------------------------------
3106
3107 IF ( ETP <= 0.0 ) THEN
3108 BETA = 0.0
3109 IF ( ETP < 0.0 ) THEN
3110 BETA = 1.0
3111 ETA = ETP
3112 END IF
3113 ELSE
3114 BETA = ETA / ETP
3115 END IF
3116
3117 ! ----------------------------------------------------------------------
3118 ! CONVERT MODELED EVAPOTRANSPIRATION COMPONENTS 'M S-1' TO 'KG M-2 S-1'.
3119 ! ----------------------------------------------------------------------
3120 EDIR = EDIR1*1000.
3121 EC = EC1*1000.
3122 DO K = 1,NSOIL
3123 ET(K) = ET1(K)*1000.
3124 END DO
3125 ETT = ETT1*1000.
3126
3127 ! ----------------------------------------------------------------------
3128 ! GET SOIL THERMAL DIFFUXIVITY/CONDUCTIVITY FOR TOP SOIL LYR,
3129 ! CALC. ADJUSTED TOP LYR SOIL TEMP AND ADJUSTED SOIL FLUX, THEN
3130 ! CALL SHFLX TO COMPUTE/UPDATE SOIL HEAT FLUX AND SOIL TEMPS.
3131 ! ----------------------------------------------------------------------
3132
3133 CALL TDFCND (DF1,SMC (1),QUARTZ,SMCMAX,SH2O (1))
3134
3135 !urban
3136 IF (VEGTYP==1) DF1=3.24
3137 !
3138
3139 ! ----------------------------------------------------------------------
3140 ! VEGETATION GREENNESS FRACTION REDUCTION IN SUBSURFACE HEAT FLUX
3141 ! VIA REDUCTION FACTOR, WHICH IS CONVENIENT TO APPLY HERE TO THERMAL
3142 ! DIFFUSIVITY THAT IS LATER USED IN HRT TO COMPUTE SUB SFC HEAT FLUX
3143 ! (SEE ADDITIONAL COMMENTS ON VEG EFFECT SUB-SFC HEAT FLX IN
3144 ! ROUTINE SFLX)
3145 ! ----------------------------------------------------------------------
3146 DF1 = DF1 * EXP (SBETA * SHDFAC)
3147 ! ----------------------------------------------------------------------
3148 ! COMPUTE INTERMEDIATE TERMS PASSED TO ROUTINE HRT (VIA ROUTINE
3149 ! SHFLX BELOW) FOR USE IN COMPUTING SUBSURFACE HEAT FLUX IN HRT
3150 ! ----------------------------------------------------------------------
3151 YYNUM = FDOWN - EMISSI*SIGMA * T24
3152 YY = SFCTMP + (YYNUM / RCH + TH2- SFCTMP - BETA * EPSCA) / RR
3153
3154 ZZ1 = DF1 / ( -0.5 * ZSOIL (1) * RCH * RR ) + 1.0
3155
3156 !urban
3157 CALL SHFLX (SSOIL,STC,SMC,SMCMAX,NSOIL,T1,DT,YY,ZZ1,ZSOIL, &
3158 TBOT,ZBOT,SMCWLT,PSISAT,SH2O,BEXP,F1,DF1,ICE, &
3159 QUARTZ,CSOIL,VEGTYP)
3160
3161 ! ----------------------------------------------------------------------
3162 ! SET FLX1 AND FLX3 (SNOPACK PHASE CHANGE HEAT FLUXES) TO ZERO SINCE
3163 ! THEY ARE NOT USED HERE IN SNOPAC. FLX2 (FREEZING RAIN HEAT FLUX) WAS
3164 ! SIMILARLY INITIALIZED IN THE PENMAN ROUTINE.
3165 ! ----------------------------------------------------------------------
3166 FLX1 = 0.0
3167 FLX3 = 0.0
3168
3169 ! ----------------------------------------------------------------------
3170 END SUBROUTINE NOPAC
3171 ! ----------------------------------------------------------------------
3172
3173 SUBROUTINE PENMAN (SFCTMP,SFCPRS,CH,T2V,TH2,PRCP,FDOWN,T24,SSOIL, &
3174 & Q2,Q2SAT,ETP,RCH,EPSCA,RR,SNOWNG,FRZGRA, &
3175 & DQSDT2,FLX2,EMISSI_IN,SNEQV)
3176
3177 ! ----------------------------------------------------------------------
3178 ! SUBROUTINE PENMAN
3179 ! ----------------------------------------------------------------------
3180 ! CALCULATE POTENTIAL EVAPORATION FOR THE CURRENT POINT. VARIOUS
3181 ! PARTIAL SUMS/PRODUCTS ARE ALSO CALCULATED AND PASSED BACK TO THE
3182 ! CALLING ROUTINE FOR LATER USE.
3183 ! ----------------------------------------------------------------------
3184 IMPLICIT NONE
3185 LOGICAL, INTENT(IN) :: SNOWNG, FRZGRA
3186 REAL, INTENT(IN) :: CH, DQSDT2,FDOWN,PRCP, &
3187 Q2, Q2SAT,SSOIL, SFCPRS, SFCTMP, &
3188 T2V, TH2,EMISSI_IN,SNEQV
3189 REAL, INTENT(OUT) :: EPSCA,ETP,FLX2,RCH,RR,T24
3190 REAL :: A, DELTA, FNET,RAD,RHO,EMISSI,ELCP1,LVS
3191
3192 REAL, PARAMETER :: ELCP = 2.4888E+3, LSUBC = 2.501000E+6,CP = 1004.6
3193 REAL, PARAMETER :: LSUBS = 2.83E+6
3194
3195 ! ----------------------------------------------------------------------
3196 ! EXECUTABLE CODE BEGINS HERE:
3197 ! ----------------------------------------------------------------------
3198 ! ----------------------------------------------------------------------
3199 ! PREPARE PARTIAL QUANTITIES FOR PENMAN EQUATION.
3200 ! ----------------------------------------------------------------------
3201 IF(SNEQV == 0.)THEN
3202 EMISSI=EMISSI_IN
3203 ELCP1=ELCP
3204 LVS=LSUBC
3205 ELSE
3206 EMISSI=0.95 ! FOR SNOW
3207 ELCP1=ELCP*LSUBS/LSUBC
3208 LVS=LSUBS
3209 ENDIF
3210
3211 FLX2 = 0.0
3212 ! DELTA = ELCP * DQSDT2
3213 DELTA = ELCP1 * DQSDT2
3214 T24 = SFCTMP * SFCTMP * SFCTMP * SFCTMP
3215 ! RR = T24 * 6.48E-8 / (SFCPRS * CH) + 1.0
3216 RR = EMISSI*T24 * 6.48E-8 / (SFCPRS * CH) + 1.0
3217 RHO = SFCPRS / (RD * T2V)
3218
3219 ! ----------------------------------------------------------------------
3220 ! ADJUST THE PARTIAL SUMS / PRODUCTS WITH THE LATENT HEAT
3221 ! EFFECTS CAUSED BY FALLING PRECIPITATION.
3222 ! ----------------------------------------------------------------------
3223 RCH = RHO * CP * CH
3224 IF (.NOT. SNOWNG) THEN
3225 IF (PRCP > 0.0) RR = RR + CPH2O * PRCP / RCH
3226 ELSE
3227 RR = RR + CPICE * PRCP / RCH
3228 END IF
3229
3230 ! ----------------------------------------------------------------------
3231 ! INCLUDE THE LATENT HEAT EFFECTS OF FRZNG RAIN CONVERTING TO ICE ON
3232 ! IMPACT IN THE CALCULATION OF FLX2 AND FNET.
3233 ! ----------------------------------------------------------------------
3234 ! FNET = FDOWN - SIGMA * T24- SSOIL
3235 FNET = FDOWN - EMISSI*SIGMA * T24- SSOIL
3236 IF (FRZGRA) THEN
3237 FLX2 = - LSUBF * PRCP
3238 FNET = FNET - FLX2
3239 ! ----------------------------------------------------------------------
3240 ! FINISH PENMAN EQUATION CALCULATIONS.
3241 ! ----------------------------------------------------------------------
3242 END IF
3243 RAD = FNET / RCH + TH2- SFCTMP
3244 ! A = ELCP * (Q2SAT - Q2)
3245 A = ELCP1 * (Q2SAT - Q2)
3246 EPSCA = (A * RR + RAD * DELTA) / (DELTA + RR)
3247 ! ETP = EPSCA * RCH / LSUBC
3248 ETP = EPSCA * RCH / LVS
3249
3250 ! ----------------------------------------------------------------------
3251 END SUBROUTINE PENMAN
3252 ! ----------------------------------------------------------------------
3253
3254 SUBROUTINE REDPRM (VEGTYP,SOILTYP,SLOPETYP,CFACTR,CMCMAX,RSMAX, &
3255 TOPT, &
3256 REFKDT,KDT,SBETA, SHDFAC,RSMIN,RGL,HS,ZBOT,FRZX, &
3257 PSISAT,SLOPE,SNUP,SALP,BEXP,DKSAT,DWSAT, &
3258 SMCMAX,SMCWLT,SMCREF,SMCDRY,F1,QUARTZ,FXEXP, &
3259 RTDIS,SLDPTH,ZSOIL, NROOT,NSOIL,Z0BRD,CZIL,LAI, &
3260 CSOIL,ALBBRD,PTU,LLANDUSE,LSOIL,LOCAL)
3261
3262 IMPLICIT NONE
3263 ! ----------------------------------------------------------------------
3264 ! Internally set (default valuess)
3265 ! all soil and vegetation parameters required for the execusion oF
3266 ! the Noah lsm are defined in VEGPARM.TBL, SOILPARM.TB, and GENPARM.TBL.
3267 ! ----------------------------------------------------------------------
3268 ! Vegetation parameters:
3269 ! ALBBRD: SFC background snow-free albedo
3270 ! CMXTBL: MAX CNPY Capacity
3271 ! Z0BRD: Background roughness length
3272 ! SHDFAC: Green vegetation fraction
3273 ! NROOT: Rooting depth
3274 ! RSMIN: Mimimum stomatal resistance
3275 ! RSMAX: Max. stomatal resistance
3276 ! RGL: Parameters used in radiation stress function
3277 ! HS: Parameter used in vapor pressure deficit functio
3278 ! TOPT: Optimum transpiration air temperature.
3279 ! CMCMAX: Maximum canopy water capacity
3280 ! CFACTR: Parameter used in the canopy inteception calculation
3281 ! SNUP: Threshold snow depth (in water equivalent m) that
3282 ! implies 100 percent snow cover
3283 ! LAI: Leaf area index
3284 !
3285 ! ----------------------------------------------------------------------
3286 ! Soil parameters:
3287 ! SMCMAX: MAX soil moisture content (porosity)
3288 ! SMCREF: Reference soil moisture (field capacity)
3289 ! SMCWLT: Wilting point soil moisture
3290 ! SMCWLT: Air dry soil moist content limits
3291 ! SSATPSI: SAT (saturation) soil potential
3292 ! DKSAT: SAT soil conductivity
3293 ! BEXP: B parameter
3294 ! SSATDW: SAT soil diffusivity
3295 ! F1: Soil thermal diffusivity/conductivity coef.
3296 ! QUARTZ: Soil quartz content
3297 ! Modified by F. Chen (12/22/97) to use the STATSGO soil map
3298 ! Modified By F. Chen (01/22/00) to include PLaya, Lava, and White San
3299 ! Modified By F. Chen (08/05/02) to include additional parameters for the Noah
3300 ! NOTE: SATDW = BB*SATDK*(SATPSI/MAXSMC)
3301 ! F11 = ALOG10(SATPSI) + BB*ALOG10(MAXSMC) + 2.0
3302 ! REFSMC1=MAXSMC*(5.79E-9/SATDK)**(1/(2*BB+3)) 5.79E-9 m/s= 0.5 mm
3303 ! REFSMC=REFSMC1+1./3.(MAXSMC-REFSMC1)
3304 ! WLTSMC1=MAXSMC*(200./SATPSI)**(-1./BB) (Wetzel and Chang, 198
3305 ! WLTSMC=WLTSMC1-0.5*WLTSMC1
3306 ! Note: the values for playa is set for it to have a thermal conductivit
3307 ! as sand and to have a hydrulic conductivity as clay
3308 !
3309 ! ----------------------------------------------------------------------
3310 ! Class parameter 'SLOPETYP' was included to estimate linear reservoir
3311 ! coefficient 'SLOPE' to the baseflow runoff out of the bottom layer.
3312 ! lowest class (slopetyp=0) means highest slope parameter = 1.
3313 ! definition of slopetyp from 'zobler' slope type:
3314 ! slope class percent slope
3315 ! 1 0-8
3316 ! 2 8-30
3317 ! 3 > 30
3318 ! 4 0-30
3319 ! 5 0-8 & > 30
3320 ! 6 8-30 & > 30
3321 ! 7 0-8, 8-30, > 30
3322 ! 9 GLACIAL ICE
3323 ! BLANK OCEAN/SEA
3324 ! SLOPE_DATA: linear reservoir coefficient
3325 ! SBETA_DATA: parameter used to caluculate vegetation effect on soil heat
3326 ! FXEXP_DAT: soil evaporation exponent used in DEVAP
3327 ! CSOIL_DATA: soil heat capacity [J M-3 K-1]
3328 ! SALP_DATA: shape parameter of distribution function of snow cover
3329 ! REFDK_DATA and REFKDT_DATA: parameters in the surface runoff parameteriz
3330 ! FRZK_DATA: frozen ground parameter
3331 ! ZBOT_DATA: depth[M] of lower boundary soil temperature
3332 ! CZIL_DATA: calculate roughness length of heat
3333 ! SMLOW_DATA and MHIGH_DATA: two soil moisture wilt, soil moisture referen
3334 ! parameters
3335 ! Set maximum number of soil-, veg-, and slopetyp in data statement.
3336 ! ----------------------------------------------------------------------
3337 INTEGER, PARAMETER :: MAX_SLOPETYP=30,MAX_SOILTYP=30,MAX_VEGTYP=30
3338 LOGICAL :: LOCAL
3339 CHARACTER (LEN=4), INTENT(IN):: LLANDUSE, LSOIL
3340
3341 ! Veg parameters
3342 INTEGER, INTENT(IN) :: VEGTYP
3343 INTEGER, INTENT(OUT) :: NROOT
3344 REAL, INTENT(OUT) :: HS,LAI,RSMIN,RGL,SHDFAC,SNUP,Z0BRD, &
3345 CMCMAX,RSMAX,TOPT,ALBBRD
3346 ! Soil parameters
3347 INTEGER, INTENT(IN) :: SOILTYP
3348 REAL, INTENT(OUT) :: BEXP,DKSAT,DWSAT,F1,QUARTZ,SMCDRY, &
3349 SMCMAX,SMCREF,SMCWLT,PSISAT
3350 ! General parameters
3351 INTEGER, INTENT(IN) :: SLOPETYP,NSOIL
3352 INTEGER :: I
3353
3354 REAL, INTENT(OUT) :: SLOPE,CZIL,SBETA,FXEXP, &
3355 CSOIL,SALP,FRZX,KDT,CFACTR, &
3356 ZBOT,REFKDT,PTU
3357 REAL,DIMENSION(1:NSOIL),INTENT(IN) :: SLDPTH,ZSOIL
3358 REAL,DIMENSION(1:NSOIL),INTENT(OUT):: RTDIS
3359 REAL :: FRZFACT,FRZK,REFDK
3360
3361 ! SAVE
3362 ! ----------------------------------------------------------------------
3363 !
3364 IF (SOILTYP .gt. SLCATS) THEN
3365 CALL wrf_error_fatal ( 'Warning: too many input soil types' )
3366 END IF
3367 IF (VEGTYP .gt. LUCATS) THEN
3368 CALL wrf_error_fatal ( 'Warning: too many input landuse types' )
3369 END IF
3370 IF (SLOPETYP .gt. SLPCATS) THEN
3371 CALL wrf_error_fatal ( 'Warning: too many input slope types' )
3372 END IF
3373
3374 ! ----------------------------------------------------------------------
3375 ! SET-UP SOIL PARAMETERS
3376 ! ----------------------------------------------------------------------
3377 CSOIL = CSOIL_DATA
3378 BEXP = BB (SOILTYP)
3379 DKSAT = SATDK (SOILTYP)
3380 DWSAT = SATDW (SOILTYP)
3381 F1 = F11 (SOILTYP)
3382 PSISAT = SATPSI (SOILTYP)
3383 QUARTZ = QTZ (SOILTYP)
3384 SMCDRY = DRYSMC (SOILTYP)
3385 SMCMAX = MAXSMC (SOILTYP)
3386 SMCREF = REFSMC (SOILTYP)
3387 SMCWLT = WLTSMC (SOILTYP)
3388 ! ----------------------------------------------------------------------
3389 ! Set-up universal parameters (not dependent on SOILTYP, VEGTYP or
3390 ! SLOPETYP)
3391 ! ----------------------------------------------------------------------
3392 ZBOT = ZBOT_DATA
3393 SALP = SALP_DATA
3394 SBETA = SBETA_DATA
3395 REFDK = REFDK_DATA
3396 FRZK = FRZK_DATA
3397 FXEXP = FXEXP_DATA
3398 REFKDT = REFKDT_DATA
3399 PTU = 0. ! (not used yet) to satisify intent(out)
3400 KDT = REFKDT * DKSAT / REFDK
3401 CZIL = CZIL_DATA
3402 SLOPE = SLOPE_DATA (SLOPETYP)
3403
3404 ! ----------------------------------------------------------------------
3405 ! TO ADJUST FRZK PARAMETER TO ACTUAL SOIL TYPE: FRZK * FRZFACT
3406 ! ----------------------------------------------------------------------
3407 FRZFACT = (SMCMAX / SMCREF) * (0.412 / 0.468)
3408 FRZX = FRZK * FRZFACT
3409
3410 ! ----------------------------------------------------------------------
3411 ! SET-UP VEGETATION PARAMETERS
3412 ! ----------------------------------------------------------------------
3413 TOPT = TOPT_DATA
3414 CMCMAX = CMCMAX_DATA
3415 CFACTR = CFACTR_DATA
3416 RSMAX = RSMAX_DATA
3417 NROOT = NROTBL (VEGTYP)
3418 SNUP = SNUPTBL (VEGTYP)
3419 RSMIN = RSTBL (VEGTYP)
3420 RGL = RGLTBL (VEGTYP)
3421 HS = HSTBL (VEGTYP)
3422 LAI = LAITBL (VEGTYP)
3423 IF(LOCAL) THEN
3424 ALBBRD = ALBTBL(VEGTYP)
3425 Z0BRD = Z0TBL(VEGTYP)
3426 SHDFAC = SHDTBL(VEGTYP)
3427 ENDIF
3428
3429 IF (VEGTYP .eq. BARE) SHDFAC = 0.0
3430 IF (NROOT .gt. NSOIL) THEN
3431 WRITE (err_message,*) 'Error: too many root layers ', &
3432 NSOIL,NROOT
3433 CALL wrf_error_fatal ( err_message )
3434 ! ----------------------------------------------------------------------
3435 ! CALCULATE ROOT DISTRIBUTION. PRESENT VERSION ASSUMES UNIFORM
3436 ! DISTRIBUTION BASED ON SOIL LAYER DEPTHS.
3437 ! ----------------------------------------------------------------------
3438 END IF
3439 DO I = 1,NROOT
3440 RTDIS (I) = - SLDPTH (I)/ ZSOIL (NROOT)
3441 ! ----------------------------------------------------------------------
3442 ! SET-UP SLOPE PARAMETER
3443 ! ----------------------------------------------------------------------
3444 END DO
3445
3446 ! print*,'end of PRMRED'
3447 ! print*,'VEGTYP',VEGTYP,'SOILTYP',SOILTYP,'SLOPETYP',SLOPETYP, &
3448 ! & 'CFACTR',CFACTR,'CMCMAX',CMCMAX,'RSMAX',RSMAX,'TOPT',TOPT, &
3449 ! & 'REFKDT',REFKDT,'KDT',KDT,'SBETA',SBETA, 'SHDFAC',SHDFAC, &
3450 ! & 'RSMIN',RSMIN,'RGL',RGL,'HS',HS,'ZBOT',ZBOT,'FRZX',FRZX, &
3451 ! & 'PSISAT',PSISAT,'SLOPE',SLOPE,'SNUP',SNUP,'SALP',SALP,'BEXP', &
3452 ! & BEXP, &
3453 ! & 'DKSAT',DKSAT,'DWSAT',DWSAT, &
3454 ! & 'SMCMAX',SMCMAX,'SMCWLT',SMCWLT,'SMCREF',SMCREF,'SMCDRY',SMCDRY, &
3455 ! & 'F1',F1,'QUARTZ',QUARTZ,'FXEXP',FXEXP, &
3456 ! & 'RTDIS',RTDIS,'SLDPTH',SLDPTH,'ZSOIL',ZSOIL, 'NROOT',NROOT, &
3457 ! & 'NSOIL',NSOIL,'Z0',Z0,'CZIL',CZIL,'LAI',LAI, &
3458 ! & 'CSOIL',CSOIL,'PTU',PTU, &
3459 ! & 'LOCAL', LOCAL
3460
3461 END SUBROUTINE REDPRM
3462
3463 SUBROUTINE ROSR12 (P,A,B,C,D,DELTA,NSOIL)
3464
3465 ! ----------------------------------------------------------------------
3466 ! SUBROUTINE ROSR12
3467 ! ----------------------------------------------------------------------
3468 ! INVERT (SOLVE) THE TRI-DIAGONAL MATRIX PROBLEM SHOWN BELOW:
3469 ! ### ### ### ### ### ###
3470 ! #B(1), C(1), 0 , 0 , 0 , . . . , 0 # # # # #
3471 ! #A(2), B(2), C(2), 0 , 0 , . . . , 0 # # # # #
3472 ! # 0 , A(3), B(3), C(3), 0 , . . . , 0 # # # # D(3) #
3473 ! # 0 , 0 , A(4), B(4), C(4), . . . , 0 # # P(4) # # D(4) #
3474 ! # 0 , 0 , 0 , A(5), B(5), . . . , 0 # # P(5) # # D(5) #
3475 ! # . . # # . # = # . #
3476 ! # . . # # . # # . #
3477 ! # . . # # . # # . #
3478 ! # 0 , . . . , 0 , A(M-2), B(M-2), C(M-2), 0 # #P(M-2)# #D(M-2)#
3479 ! # 0 , . . . , 0 , 0 , A(M-1), B(M-1), C(M-1)# #P(M-1)# #D(M-1)#
3480 ! # 0 , . . . , 0 , 0 , 0 , A(M) , B(M) # # P(M) # # D(M) #
3481 ! ### ### ### ### ### ###
3482 ! ----------------------------------------------------------------------
3483 IMPLICIT NONE
3484
3485 INTEGER, INTENT(IN) :: NSOIL
3486 INTEGER :: K, KK
3487
3488 REAL, DIMENSION(1:NSOIL), INTENT(IN):: A, B, D
3489 REAL, DIMENSION(1:NSOIL),INTENT(INOUT):: C,P,DELTA
3490
3491 ! ----------------------------------------------------------------------
3492 ! INITIALIZE EQN COEF C FOR THE LOWEST SOIL LAYER
3493 ! ----------------------------------------------------------------------
3494 C (NSOIL) = 0.0
3495 P (1) = - C (1) / B (1)
3496 ! ----------------------------------------------------------------------
3497 ! SOLVE THE COEFS FOR THE 1ST SOIL LAYER
3498 ! ----------------------------------------------------------------------
3499
3500 ! ----------------------------------------------------------------------
3501 ! SOLVE THE COEFS FOR SOIL LAYERS 2 THRU NSOIL
3502 ! ----------------------------------------------------------------------
3503 DELTA (1) = D (1) / B (1)
3504 DO K = 2,NSOIL
3505 P (K) = - C (K) * ( 1.0 / (B (K) + A (K) * P (K -1)) )
3506 DELTA (K) = (D (K) - A (K)* DELTA (K -1))* (1.0/ (B (K) + A (K)&
3507 * P (K -1)))
3508 END DO
3509 ! ----------------------------------------------------------------------
3510 ! SET P TO DELTA FOR LOWEST SOIL LAYER
3511 ! ----------------------------------------------------------------------
3512 P (NSOIL) = DELTA (NSOIL)
3513
3514 ! ----------------------------------------------------------------------
3515 ! ADJUST P FOR SOIL LAYERS 2 THRU NSOIL
3516 ! ----------------------------------------------------------------------
3517 DO K = 2,NSOIL
3518 KK = NSOIL - K + 1
3519 P (KK) = P (KK) * P (KK +1) + DELTA (KK)
3520 END DO
3521 ! ----------------------------------------------------------------------
3522 END SUBROUTINE ROSR12
3523 ! ----------------------------------------------------------------------
3524
3525
3526 SUBROUTINE SHFLX (SSOIL,STC,SMC,SMCMAX,NSOIL,T1,DT,YY,ZZ1,ZSOIL, &
3527 TBOT,ZBOT,SMCWLT,PSISAT,SH2O,BEXP,F1,DF1,ICE, &
3528 QUARTZ,CSOIL,VEGTYP)
3529
3530 ! ----------------------------------------------------------------------
3531 ! SUBROUTINE SHFLX
3532 ! ----------------------------------------------------------------------
3533 ! UPDATE THE TEMPERATURE STATE OF THE SOIL COLUMN BASED ON THE THERMAL
3534 ! DIFFUSION EQUATION AND UPDATE THE FROZEN SOIL MOISTURE CONTENT BASED
3535 ! ON THE TEMPERATURE.
3536 ! ----------------------------------------------------------------------
3537 IMPLICIT NONE
3538
3539 INTEGER, INTENT(IN) :: ICE, NSOIL, VEGTYP
3540 INTEGER :: I
3541
3542 REAL, INTENT(IN) :: BEXP,CSOIL,DF1,DT,F1,PSISAT,QUARTZ, &
3543 SMCMAX, SMCWLT, TBOT,YY, ZBOT,ZZ1
3544 REAL, INTENT(INOUT) :: T1
3545 REAL, INTENT(OUT) :: SSOIL
3546 REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SMC,ZSOIL
3547 REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: SH2O
3548 REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: STC
3549 REAL, DIMENSION(1:NSOIL) :: AI, BI, CI, STCF,RHSTS
3550 REAL, PARAMETER :: T0 = 273.15
3551
3552 ! ----------------------------------------------------------------------
3553 ! HRT ROUTINE CALCS THE RIGHT HAND SIDE OF THE SOIL TEMP DIF EQN
3554 ! ----------------------------------------------------------------------
3555
3556 ! ----------------------------------------------------------------------
3557 ! SEA-ICE CASE
3558 ! ----------------------------------------------------------------------
3559 IF (ICE == 1) THEN
3560
3561 CALL HRTICE (RHSTS,STC,NSOIL,ZSOIL,YY,ZZ1,DF1,AI,BI,CI)
3562
3563 CALL HSTEP (STCF,STC,RHSTS,DT,NSOIL,AI,BI,CI)
3564
3565 ! ----------------------------------------------------------------------
3566 ! LAND-MASS CASE
3567 ! ----------------------------------------------------------------------
3568 ELSE
3569 CALL HRT (RHSTS,STC,SMC,SMCMAX,NSOIL,ZSOIL,YY,ZZ1,TBOT, &
3570 ZBOT,PSISAT,SH2O,DT, &
3571 BEXP,F1,DF1,QUARTZ,CSOIL,AI,BI,CI,VEGTYP)
3572
3573 CALL HSTEP (STCF,STC,RHSTS,DT,NSOIL,AI,BI,CI)
3574 END IF
3575 DO I = 1,NSOIL
3576 STC (I) = STCF (I)
3577 END DO
3578
3579 ! ----------------------------------------------------------------------
3580 ! IN THE NO SNOWPACK CASE (VIA ROUTINE NOPAC BRANCH,) UPDATE THE GRND
3581 ! (SKIN) TEMPERATURE HERE IN RESPONSE TO THE UPDATED SOIL TEMPERATURE
3582 ! PROFILE ABOVE. (NOTE: INSPECTION OF ROUTINE SNOPAC SHOWS THAT T1
3583 ! BELOW IS A DUMMY VARIABLE ONLY, AS SKIN TEMPERATURE IS UPDATED
3584 ! DIFFERENTLY IN ROUTINE SNOPAC)
3585 ! ----------------------------------------------------------------------
3586 ! ----------------------------------------------------------------------
3587 ! CALCULATE SURFACE SOIL HEAT FLUX
3588 ! ----------------------------------------------------------------------
3589 T1 = (YY + (ZZ1- 1.0) * STC (1)) / ZZ1
3590 SSOIL = DF1 * (STC (1) - T1) / (0.5 * ZSOIL (1))
3591
3592 ! ----------------------------------------------------------------------
3593 END SUBROUTINE SHFLX
3594 ! ----------------------------------------------------------------------
3595
3596 SUBROUTINE SMFLX (SMC,NSOIL,CMC,DT,PRCP1,ZSOIL, &
3597 & SH2O,SLOPE,KDT,FRZFACT, &
3598 & SMCMAX,BEXP,SMCWLT,DKSAT,DWSAT, &
3599 & SHDFAC,CMCMAX, &
3600 & RUNOFF1,RUNOFF2,RUNOFF3, &
3601 & EDIR,EC,ET, &
3602 & DRIP)
3603
3604 ! ----------------------------------------------------------------------
3605 ! SUBROUTINE SMFLX
3606 ! ----------------------------------------------------------------------
3607 ! CALCULATE SOIL MOISTURE FLUX. THE SOIL MOISTURE CONTENT (SMC - A PER
3608 ! UNIT VOLUME MEASUREMENT) IS A DEPENDENT VARIABLE THAT IS UPDATED WITH
3609 ! PROGNOSTIC EQNS. THE CANOPY MOISTURE CONTENT (CMC) IS ALSO UPDATED.
3610 ! FROZEN GROUND VERSION: NEW STATES ADDED: SH2O, AND FROZEN GROUND
3611 ! CORRECTION FACTOR, FRZFACT AND PARAMETER SLOPE.
3612 ! ----------------------------------------------------------------------
3613 IMPLICIT NONE
3614
3615 INTEGER, INTENT(IN) :: NSOIL
3616 INTEGER :: I,K
3617
3618 REAL, INTENT(IN) :: BEXP, CMCMAX, DKSAT,DWSAT, DT, EC, EDIR, &
3619 KDT, PRCP1, SHDFAC, SLOPE, SMCMAX, SMCWLT
3620 REAL, INTENT(OUT) :: DRIP, RUNOFF1, RUNOFF2, RUNOFF3
3621 REAL, INTENT(INOUT) :: CMC
3622 REAL, DIMENSION(1:NSOIL), INTENT(IN) :: ET,ZSOIL
3623 REAL, DIMENSION(1:NSOIL), INTENT(INOUT):: SMC, SH2O
3624 REAL, DIMENSION(1:NSOIL) :: AI, BI, CI, STCF,RHSTS, RHSTT, &
3625 SICE, SH2OA, SH2OFG
3626 REAL :: DUMMY, EXCESS,FRZFACT,PCPDRP,RHSCT,TRHSCT
3627
3628 ! ----------------------------------------------------------------------
3629 ! EXECUTABLE CODE BEGINS HERE.
3630 ! ----------------------------------------------------------------------
3631 ! ----------------------------------------------------------------------
3632 ! COMPUTE THE RIGHT HAND SIDE OF THE CANOPY EQN TERM ( RHSCT )
3633 ! ----------------------------------------------------------------------
3634 DUMMY = 0.
3635
3636 ! ----------------------------------------------------------------------
3637 ! CONVERT RHSCT (A RATE) TO TRHSCT (AN AMOUNT) AND ADD IT TO EXISTING
3638 ! CMC. IF RESULTING AMT EXCEEDS MAX CAPACITY, IT BECOMES DRIP AND WILL
3639 ! FALL TO THE GRND.
3640 ! ----------------------------------------------------------------------
3641 RHSCT = SHDFAC * PRCP1- EC
3642 DRIP = 0.
3643 TRHSCT = DT * RHSCT
3644 EXCESS = CMC + TRHSCT
3645
3646 ! ----------------------------------------------------------------------
3647 ! PCPDRP IS THE COMBINED PRCP1 AND DRIP (FROM CMC) THAT GOES INTO THE
3648 ! SOIL
3649 ! ----------------------------------------------------------------------
3650 IF (EXCESS > CMCMAX) DRIP = EXCESS - CMCMAX
3651 PCPDRP = (1. - SHDFAC) * PRCP1+ DRIP / DT
3652
3653 ! ----------------------------------------------------------------------
3654 ! STORE ICE CONTENT AT EACH SOIL LAYER BEFORE CALLING SRT and SSTEP
3655 !
3656 DO I = 1,NSOIL
3657 SICE (I) = SMC (I) - SH2O (I)
3658 END DO
3659 ! ----------------------------------------------------------------------
3660 ! CALL SUBROUTINES SRT AND SSTEP TO SOLVE THE SOIL MOISTURE
3661 ! TENDENCY EQUATIONS.
3662 ! IF THE INFILTRATING PRECIP RATE IS NONTRIVIAL,
3663 ! (WE CONSIDER NONTRIVIAL TO BE A PRECIP TOTAL OVER THE TIME STEP
3664 ! EXCEEDING ONE ONE-THOUSANDTH OF THE WATER HOLDING CAPACITY OF
3665 ! THE FIRST SOIL LAYER)
3666 ! THEN CALL THE SRT/SSTEP SUBROUTINE PAIR TWICE IN THE MANNER OF
3667 ! TIME SCHEME "F" (IMPLICIT STATE, AVERAGED COEFFICIENT)
3668 ! OF SECTION 2 OF KALNAY AND KANAMITSU (1988, MWR, VOL 116,
3669 ! PAGES 1945-1958)TO MINIMIZE 2-DELTA-T OSCILLATIONS IN THE
3670 ! SOIL MOISTURE VALUE OF THE TOP SOIL LAYER THAT CAN ARISE BECAUSE
3671 ! OF THE EXTREME NONLINEAR DEPENDENCE OF THE SOIL HYDRAULIC
3672 ! DIFFUSIVITY COEFFICIENT AND THE HYDRAULIC CONDUCTIVITY ON THE
3673 ! SOIL MOISTURE STATE
3674 ! OTHERWISE CALL THE SRT/SSTEP SUBROUTINE PAIR ONCE IN THE MANNER OF
3675 ! TIME SCHEME "D" (IMPLICIT STATE, EXPLICIT COEFFICIENT)
3676 ! OF SECTION 2 OF KALNAY AND KANAMITSU
3677 ! PCPDRP IS UNITS OF KG/M**2/S OR MM/S, ZSOIL IS NEGATIVE DEPTH IN M
3678 ! ----------------------------------------------------------------------
3679 ! IF ( PCPDRP .GT. 0.0 ) THEN
3680
3681 ! ----------------------------------------------------------------------
3682 ! FROZEN GROUND VERSION:
3683 ! SMC STATES REPLACED BY SH2O STATES IN SRT SUBR. SH2O & SICE STATES
3684 ! INC&UDED IN SSTEP SUBR. FROZEN GROUND CORRECTION FACTOR, FRZFACT
3685 ! ADDED. ALL WATER BALANCE CALCULATIONS USING UNFROZEN WATER
3686 ! ----------------------------------------------------------------------
3687 IF ( (PCPDRP * DT) > (0.001*1000.0* (- ZSOIL (1))* SMCMAX) ) THEN
3688 CALL SRT (RHSTT,EDIR,ET,SH2O,SH2O,NSOIL,PCPDRP,ZSOIL, &
3689 DWSAT,DKSAT,SMCMAX,BEXP,RUNOFF1, &
3690 RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZFACT,SICE,AI,BI,CI)
3691 CALL SSTEP (SH2OFG,SH2O,DUMMY,RHSTT,RHSCT,DT,NSOIL,SMCMAX, &
3692 CMCMAX,RUNOFF3,ZSOIL,SMC,SICE,AI,BI,CI)
3693 DO K = 1,NSOIL
3694 SH2OA (K) = (SH2O (K) + SH2OFG (K)) * 0.5
3695 END DO
3696 CALL SRT (RHSTT,EDIR,ET,SH2O,SH2OA,NSOIL,PCPDRP,ZSOIL, &
3697 DWSAT,DKSAT,SMCMAX,BEXP,RUNOFF1, &
3698 RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZFACT,SICE,AI,BI,CI)
3699 CALL SSTEP (SH2O,SH2O,CMC,RHSTT,RHSCT,DT,NSOIL,SMCMAX, &
3700 CMCMAX,RUNOFF3,ZSOIL,SMC,SICE,AI,BI,CI)
3701
3702 ELSE
3703 CALL SRT (RHSTT,EDIR,ET,SH2O,SH2O,NSOIL,PCPDRP,ZSOIL, &
3704 DWSAT,DKSAT,SMCMAX,BEXP,RUNOFF1, &
3705 RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZFACT,SICE,AI,BI,CI)
3706 CALL SSTEP (SH2O,SH2O,CMC,RHSTT,RHSCT,DT,NSOIL,SMCMAX, &
3707 CMCMAX,RUNOFF3,ZSOIL,SMC,SICE,AI,BI,CI)
3708 ! RUNOF = RUNOFF
3709
3710 END IF
3711
3712 ! ----------------------------------------------------------------------
3713 END SUBROUTINE SMFLX
3714 ! ----------------------------------------------------------------------
3715
3716
3717 SUBROUTINE SNFRAC (SNEQV,SNUP,SALP,SNOWH,SNCOVR)
3718
3719 ! ----------------------------------------------------------------------
3720 ! SUBROUTINE SNFRAC
3721 ! ----------------------------------------------------------------------
3722 ! CALCULATE SNOW FRACTION (0 -> 1)
3723 ! SNEQV SNOW WATER EQUIVALENT (M)
3724 ! SNUP THRESHOLD SNEQV DEPTH ABOVE WHICH SNCOVR=1
3725 ! SALP TUNING PARAMETER
3726 ! SNCOVR FRACTIONAL SNOW COVER
3727 ! ----------------------------------------------------------------------
3728 IMPLICIT NONE
3729
3730 REAL, INTENT(IN) :: SNEQV,SNUP,SALP,SNOWH
3731 REAL, INTENT(OUT) :: SNCOVR
3732 REAL :: RSNOW, Z0N
3733
3734 ! ----------------------------------------------------------------------
3735 ! SNUP IS VEG-CLASS DEPENDENT SNOWDEPTH THRESHHOLD (SET IN ROUTINE
3736 ! REDPRM) ABOVE WHICH SNOCVR=1.
3737 ! ----------------------------------------------------------------------
3738 IF (SNEQV < SNUP) THEN
3739 RSNOW = SNEQV / SNUP
3740 SNCOVR = 1. - ( EXP ( - SALP * RSNOW) - RSNOW * EXP ( - SALP))
3741 ELSE
3742 SNCOVR = 1.0
3743 END IF
3744
3745 ! FORMULATION OF DICKINSON ET AL. 1986
3746 ! Z0N = 0.035
3747
3748 ! SNCOVR=SNOWH/(SNOWH + 5*Z0N)
3749
3750 ! FORMULATION OF MARSHALL ET AL. 1994
3751 ! SNCOVR=SNEQV/(SNEQV + 2*Z0N)
3752
3753 ! ----------------------------------------------------------------------
3754 END SUBROUTINE SNFRAC
3755 ! ----------------------------------------------------------------------
3756
3757 SUBROUTINE SNKSRC (TSNSR,TAVG,SMC,SH2O,ZSOIL,NSOIL, &
3758 & SMCMAX,PSISAT,BEXP,DT,K,QTOT)
3759
3760 ! ----------------------------------------------------------------------
3761 ! SUBROUTINE SNKSRC
3762 ! ----------------------------------------------------------------------
3763 ! CALCULATE SINK/SOURCE TERM OF THE TERMAL DIFFUSION EQUATION. (SH2O) IS
3764 ! AVAILABLE LIQUED WATER.
3765 ! ----------------------------------------------------------------------
3766 IMPLICIT NONE
3767
3768 INTEGER, INTENT(IN) :: K,NSOIL
3769
3770 REAL, INTENT(IN) :: BEXP, DT, PSISAT, QTOT, SMC, SMCMAX, &
3771 TAVG
3772 REAL, INTENT(INOUT) :: SH2O
3773
3774 REAL, DIMENSION(1:NSOIL), INTENT(IN):: ZSOIL
3775
3776 REAL :: DF, DZ, DZH, FREE, TSNSR, &
3777 TDN, TM, TUP, TZ, X0, XDN, XH2O, XUP
3778
3779 REAL, PARAMETER :: DH2O = 1.0000E3, HLICE = 3.3350E5, &
3780 T0 = 2.7315E2
3781
3782 IF (K == 1) THEN
3783 DZ = - ZSOIL (1)
3784 ELSE
3785 DZ = ZSOIL (K -1) - ZSOIL (K)
3786 END IF
3787 ! ----------------------------------------------------------------------
3788 ! VIA FUNCTION FRH2O, COMPUTE POTENTIAL OR 'EQUILIBRIUM' UNFROZEN
3789 ! SUPERCOOLED FREE WATER FOR GIVEN SOIL TYPE AND SOIL LAYER TEMPERATURE.
3790 ! FUNCTION FRH20 INVOKES EQN (17) FROM V. KOREN ET AL (1999, JGR, VOL.
3791 ! 104, PG 19573). (ASIDE: LATTER EQN IN JOURNAL IN CENTIGRADE UNITS.
3792 ! ROUTINE FRH2O USE FORM OF EQN IN KELVIN UNITS.)
3793 ! ----------------------------------------------------------------------
3794 ! FREE = FRH2O(TAVG,SMC,SH2O,SMCMAX,BEXP,PSISAT)
3795
3796 ! ----------------------------------------------------------------------
3797 ! IN NEXT BLOCK OF CODE, INVOKE EQN 18 OF V. KOREN ET AL (1999, JGR,
3798 ! VOL. 104, PG 19573.) THAT IS, FIRST ESTIMATE THE NEW AMOUNTOF LIQUID
3799 ! WATER, 'XH2O', IMPLIED BY THE SUM OF (1) THE LIQUID WATER AT THE BEGIN
3800 ! OF CURRENT TIME STEP, AND (2) THE FREEZE OF THAW CHANGE IN LIQUID
3801 ! WATER IMPLIED BY THE HEAT FLUX 'QTOT' PASSED IN FROM ROUTINE HRT.
3802 ! SECOND, DETERMINE IF XH2O NEEDS TO BE BOUNDED BY 'FREE' (EQUIL AMT) OR
3803 ! IF 'FREE' NEEDS TO BE BOUNDED BY XH2O.
3804 ! ----------------------------------------------------------------------
3805 CALL FRH2O (FREE,TAVG,SMC,SH2O,SMCMAX,BEXP,PSISAT)
3806
3807 ! ----------------------------------------------------------------------
3808 ! FIRST, IF FREEZING AND REMAINING LIQUID LESS THAN LOWER BOUND, THEN
3809 ! REDUCE EXTENT OF FREEZING, THEREBY LETTING SOME OR ALL OF HEAT FLUX
3810 ! QTOT COOL THE SOIL TEMP LATER IN ROUTINE HRT.
3811 ! ----------------------------------------------------------------------
3812 XH2O = SH2O + QTOT * DT / (DH2O * HLICE * DZ)
3813 IF ( XH2O < SH2O .AND. XH2O < FREE) THEN
3814 IF ( FREE > SH2O ) THEN
3815 XH2O = SH2O
3816 ELSE
3817 XH2O = FREE
3818 END IF
3819 END IF
3820 ! ----------------------------------------------------------------------
3821 ! SECOND, IF THAWING AND THE INCREASE IN LIQUID WATER GREATER THAN UPPER
3822 ! BOUND, THEN REDUCE EXTENT OF THAW, THEREBY LETTING SOME OR ALL OF HEAT
3823 ! FLUX QTOT WARM THE SOIL TEMP LATER IN ROUTINE HRT.
3824 ! ----------------------------------------------------------------------
3825 IF ( XH2O > SH2O .AND. XH2O > FREE ) THEN
3826 IF ( FREE < SH2O ) THEN
3827 XH2O = SH2O
3828 ELSE
3829 XH2O = FREE
3830 END IF
3831 END IF
3832
3833 ! ----------------------------------------------------------------------
3834 ! CALCULATE PHASE-CHANGE HEAT SOURCE/SINK TERM FOR USE IN ROUTINE HRT
3835 ! AND UPDATE LIQUID WATER TO REFLCET FINAL FREEZE/THAW INCREMENT.
3836 ! ----------------------------------------------------------------------
3837 ! SNKSRC = -DH2O*HLICE*DZ*(XH2O-SH2O)/DT
3838 IF (XH2O < 0.) XH2O = 0.
3839 IF (XH2O > SMC) XH2O = SMC
3840 TSNSR = - DH2O * HLICE * DZ * (XH2O - SH2O)/ DT
3841 SH2O = XH2O
3842
3843 ! ----------------------------------------------------------------------
3844 END SUBROUTINE SNKSRC
3845 ! ----------------------------------------------------------------------
3846
3847 SUBROUTINE SNOPAC (ETP,ETA,PRCP,PRCPF,SNOWNG,SMC,SMCMAX,SMCWLT, &
3848 SMCREF,SMCDRY,CMC,CMCMAX,NSOIL,DT, &
3849 SBETA,DF1, &
3850 Q2,T1,SFCTMP,T24,TH2,FDOWN,F1,SSOIL,STC,EPSCA,&
3851 SFCPRS,BEXP,PC,RCH,RR,CFACTR,SNCOVR,ESD,SNDENS,&
3852 SNOWH,SH2O,SLOPE,KDT,FRZFACT,PSISAT, &
3853 ZSOIL,DWSAT,DKSAT,TBOT,ZBOT,SHDFAC,RUNOFF1, &
3854 RUNOFF2,RUNOFF3,EDIR,EC,ET,ETT,NROOT,SNOMLT, &
3855 ICE,RTDIS,QUARTZ,FXEXP,CSOIL, &
3856 BETA,DRIP,DEW,FLX1,FLX2,FLX3,ESNOW,ETNS,EMISSI,&
3857 VEGTYP)
3858
3859 ! ----------------------------------------------------------------------
3860 ! SUBROUTINE SNOPAC
3861 ! ----------------------------------------------------------------------
3862 ! CALCULATE SOIL MOISTURE AND HEAT FLUX VALUES & UPDATE SOIL MOISTURE
3863 ! CONTENT AND SOIL HEAT CONTENT VALUES FOR THE CASE WHEN A SNOW PACK IS
3864 ! PRESENT.
3865 ! ----------------------------------------------------------------------
3866 IMPLICIT NONE
3867
3868 INTEGER, INTENT(IN) :: ICE, NROOT, NSOIL,VEGTYP
3869 INTEGER :: K
3870 LOGICAL, INTENT(IN) :: SNOWNG
3871 REAL, INTENT(IN) :: BEXP,CFACTR, CMCMAX,CSOIL,DF1,DKSAT, &
3872 DT,DWSAT, EPSCA,ETP,FDOWN,F1,FXEXP, &
3873 FRZFACT,KDT,PC, PRCP,PSISAT,Q2,QUARTZ, &
3874 RCH,RR,SBETA,SFCPRS, SFCTMP, SHDFAC, &
3875 SLOPE,SMCDRY,SMCMAX,SMCREF,SMCWLT, T24, &
3876 TBOT,TH2,ZBOT,EMISSI
3877 REAL, INTENT(INOUT) :: CMC, BETA, ESD,FLX2,PRCPF,SNOWH,SNCOVR, &
3878 SNDENS, T1
3879 REAL, INTENT(OUT) :: DEW,DRIP,EC,EDIR, ETNS, ESNOW,ETT, &
3880 FLX1,FLX3, RUNOFF1,RUNOFF2,RUNOFF3, &
3881 SSOIL,SNOMLT
3882 REAL, DIMENSION(1:NSOIL),INTENT(IN) :: RTDIS,ZSOIL
3883 REAL, DIMENSION(1:NSOIL),INTENT(OUT) :: ET
3884 REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: SMC,SH2O,STC
3885 REAL, DIMENSION(1:NSOIL) :: ET1
3886 REAL :: DENOM,DSOIL,DTOT,EC1,EDIR1,ESDFLX,ETA, &
3887 ETT1, ESNOW1, ESNOW2, ETA1,ETP1,ETP2, &
3888 ETP3, ETNS1, ETANRG, ETAX, EX, FLX3X, &
3889 FRCSNO,FRCSOI, PRCP1, QSAT,RSNOW, SEH, &
3890 SNCOND,SSOIL1, T11,T12, T12A, T12AX, &
3891 T12B, T14, YY, ZZ1, EMISSI_S
3892
3893 REAL, PARAMETER :: ESDMIN = 1.E-6, LSUBC = 2.501000E+6, &
3894 LSUBS = 2.83E+6, TFREEZ = 273.15, &
3895 SNOEXP = 2.0
3896
3897 ! ----------------------------------------------------------------------
3898 ! EXECUTABLE CODE BEGINS HERE:
3899 ! INITIALIZE EVAP TERMS.
3900 ! ----------------------------------------------------------------------
3901 ! conversions:
3902 ! ESNOW [KG M-2 S-1]
3903 ! ESDFLX [KG M-2 S-1] .le. ESNOW
3904 ! ESNOW1 [M S-1]
3905 ! ESNOW2 [M]
3906 ! ETP [KG M-2 S-1]
3907 ! ETP1 [M S-1]
3908 ! ETP2 [M]
3909 ! ----------------------------------------------------------------------
3910 DEW = 0.
3911 EDIR = 0.
3912 EDIR1 = 0.
3913 EC1 = 0.
3914 EC = 0.
3915 EMISSI_S=0.95 ! For snow
3916
3917 DO K = 1,NSOIL
3918 ET (K) = 0.
3919 ET1 (K) = 0.
3920 END DO
3921 ETT = 0.
3922 ETT1 = 0.
3923 ETNS = 0.
3924 ETNS1 = 0.
3925 ESNOW = 0.
3926 ESNOW1 = 0.
3927 ESNOW2 = 0.
3928
3929 ! ----------------------------------------------------------------------
3930 ! CONVERT POTENTIAL EVAP (ETP) FROM KG M-2 S-1 TO ETP1 IN M S-1
3931 ! ----------------------------------------------------------------------
3932 PRCP1 = PRCPF *0.001
3933 ETP1 = ETP * 0.001
3934 ! ----------------------------------------------------------------------
3935 ! IF ETP<0 (DOWNWARD) THEN DEWFALL (=FROSTFALL IN THIS CASE).
3936 ! ----------------------------------------------------------------------
3937 BETA = 1.0
3938 IF (ETP <= 0.0) THEN
3939 IF(ETP == 0.) BETA = 0.0
3940 DEW = -ETP1
3941 ESNOW2 = ETP1*DT
3942 ETANRG = ETP*((1.-SNCOVR)*LSUBC + SNCOVR*LSUBS)
3943 ELSE
3944 IF (ICE /= 1) THEN
3945 IF (SNCOVR < 1.) THEN
3946 CALL EVAPO (ETNS1,SMC,NSOIL,CMC,ETP1,DT,ZSOIL, &
3947 SH2O, &
3948 SMCMAX,BEXP,PC,SMCWLT,DKSAT,DWSAT, &
3949 SMCREF,SHDFAC,CMCMAX, &
3950 SMCDRY,CFACTR, &
3951 EDIR1,EC1,ET1,ETT1,SFCTMP,Q2,NROOT,RTDIS, &
3952 FXEXP)
3953 ! ----------------------------------------------------------------------------
3954 EDIR1 = EDIR1* (1. - SNCOVR)
3955 EC1 = EC1* (1. - SNCOVR)
3956 DO K = 1,NSOIL
3957 ET1 (K) = ET1 (K)* (1. - SNCOVR)
3958 END DO
3959 ETT1 = ETT1*(1.-SNCOVR)
3960 ! ETNS1 = EDIR1+ EC1+ ETT1
3961 ETNS1 = ETNS1*(1.-SNCOVR)
3962 ! ----------------------------------------------------------------------------
3963 EDIR = EDIR1*1000.
3964 EC = EC1*1000.
3965 DO K = 1,NSOIL
3966 ET (K) = ET1 (K)*1000.
3967 END DO
3968 ETT = ETT1*1000.
3969 ETNS = ETNS1*1000.
3970 ETNS = ETNS1*1000.
3971 ! ----------------------------------------------------------------------
3972
3973 ! end IF (SNCOVR .lt. 1.)
3974 END IF
3975 ! end IF (ICE .ne. 1)
3976 END IF
3977 ESNOW = ETP*SNCOVR
3978 ESNOW1 = ESNOW*0.001
3979 ESNOW2 = ESNOW1*DT
3980 ETANRG = ESNOW*LSUBS + ETNS*LSUBC
3981 ! end IF (ETP .le. 0.0)
3982 END IF
3983
3984 ! ----------------------------------------------------------------------
3985 ! IF PRECIP IS FALLING, CALCULATE HEAT FLUX FROM SNOW SFC TO NEWLY
3986 ! ACCUMULATING PRECIP. NOTE THAT THIS REFLECTS THE FLUX APPROPRIATE FOR
3987 ! THE NOT-YET-UPDATED SKIN TEMPERATURE (T1). ASSUMES TEMPERATURE OF THE
3988 ! SNOWFALL STRIKING THE GROUND IS =SFCTMP (LOWEST MODEL LEVEL AIR TEMP).
3989 ! ----------------------------------------------------------------------
3990 FLX1 = 0.0
3991 IF (SNOWNG) THEN
3992 FLX1 = CPICE * PRCP * (T1- SFCTMP)
3993 ELSE
3994 IF (PRCP > 0.0) FLX1 = CPH2O * PRCP * (T1- SFCTMP)
3995 ! ----------------------------------------------------------------------
3996 ! CALCULATE AN 'EFFECTIVE SNOW-GRND SFC TEMP' (T12) BASED ON HEAT FLUXES
3997 ! BETWEEN THE SNOW PACK AND THE SOIL AND ON NET RADIATION.
3998 ! INCLUDE FLX1 (PRECIP-SNOW SFC) AND FLX2 (FREEZING RAIN LATENT HEAT)
3999 ! FLUXES. FLX1 FROM ABOVE, FLX2 BROUGHT IN VIA COMMOM BLOCK RITE.
4000 ! FLX2 REFLECTS FREEZING RAIN LATENT HEAT FLUX USING T1 CALCULATED IN
4001 ! PENMAN.
4002 ! ----------------------------------------------------------------------
4003 END IF
4004 DSOIL = - (0.5 * ZSOIL (1))
4005 DTOT = SNOWH + DSOIL
4006 DENOM = 1.0+ DF1 / (DTOT * RR * RCH)
4007 ! surface emissivity weighted by snow cover fraction
4008 T12A = ( (FDOWN - FLX1 - FLX2 - &
4009 & ((SNCOVR*EMISSI_S)+EMISSI*(1.0-SNCOVR))*SIGMA *T24)/RCH &
4010 & + TH2 - SFCTMP - ETANRG/RCH ) / RR
4011
4012 T12B = DF1 * STC (1) / (DTOT * RR * RCH)
4013
4014 ! ----------------------------------------------------------------------
4015 ! IF THE 'EFFECTIVE SNOW-GRND SFC TEMP' IS AT OR BELOW FREEZING, NO SNOW
4016 ! MELT WILL OCCUR. SET THE SKIN TEMP TO THIS EFFECTIVE TEMP. REDUCE
4017 ! (BY SUBLIMINATION ) OR INCREASE (BY FROST) THE DEPTH OF THE SNOWPACK,
4018 ! DEPENDING ON SIGN OF ETP.
4019 ! UPDATE SOIL HEAT FLUX (SSOIL) USING NEW SKIN TEMPERATURE (T1)
4020 ! SINCE NO SNOWMELT, SET ACCUMULATED SNOWMELT TO ZERO, SET 'EFFECTIVE'
4021 ! PRECIP FROM SNOWMELT TO ZERO, SET PHASE-CHANGE HEAT FLUX FROM SNOWMELT
4022 ! TO ZERO.
4023 ! ----------------------------------------------------------------------
4024 ! SUB-FREEZING BLOCK
4025 ! ----------------------------------------------------------------------
4026 T12 = (SFCTMP + T12A + T12B) / DENOM
4027 IF (T12 <= TFREEZ) THEN
4028 T1 = T12
4029 SSOIL = DF1 * (T1- STC (1)) / DTOT
4030 ! ESD = MAX (0.0, ESD- ETP2)
4031 ESD = MAX(0.0, ESD-ESNOW2)
4032 FLX3 = 0.0
4033 EX = 0.0
4034
4035 SNOMLT = 0.0
4036 ! ----------------------------------------------------------------------
4037 ! IF THE 'EFFECTIVE SNOW-GRND SFC TEMP' IS ABOVE FREEZING, SNOW MELT
4038 ! WILL OCCUR. CALL THE SNOW MELT RATE,EX AND AMT, SNOMLT. REVISE THE
4039 ! EFFECTIVE SNOW DEPTH. REVISE THE SKIN TEMP BECAUSE IT WOULD HAVE CHGD
4040 ! DUE TO THE LATENT HEAT RELEASED BY THE MELTING. CALC THE LATENT HEAT
4041 ! RELEASED, FLX3. SET THE EFFECTIVE PRECIP, PRCP1 TO THE SNOW MELT RATE,
4042 ! EX FOR USE IN SMFLX. ADJUSTMENT TO T1 TO ACCOUNT FOR SNOW PATCHES.
4043 ! CALCULATE QSAT VALID AT FREEZING POINT. NOTE THAT ESAT (SATURATION
4044 ! VAPOR PRESSURE) VALUE OF 6.11E+2 USED HERE IS THAT VALID AT FRZZING
4045 ! POINT. NOTE THAT ETP FROM CALL PENMAN IN SFLX IS IGNORED HERE IN
4046 ! FAVOR OF BULK ETP OVER 'OPEN WATER' AT FREEZING TEMP.
4047 ! UPDATE SOIL HEAT FLUX (S) USING NEW SKIN TEMPERATURE (T1)
4048 ! ----------------------------------------------------------------------
4049 ! ABOVE FREEZING BLOCK
4050 ! ----------------------------------------------------------------------
4051 ELSE
4052 T1 = TFREEZ * SNCOVR ** SNOEXP + T12 * (1.0- SNCOVR ** SNOEXP)
4053 BETA = 1.0
4054
4055 ! ----------------------------------------------------------------------
4056 ! IF POTENTIAL EVAP (SUBLIMATION) GREATER THAN DEPTH OF SNOWPACK.
4057 ! BETA<1
4058 ! SNOWPACK HAS SUBLIMATED AWAY, SET DEPTH TO ZERO.
4059 ! ----------------------------------------------------------------------
4060 SSOIL = DF1 * (T1- STC (1)) / DTOT
4061 IF (ESD-ESNOW2 <= ESDMIN) THEN
4062 ESD = 0.0
4063 EX = 0.0
4064 SNOMLT = 0.0
4065 ! ----------------------------------------------------------------------
4066 ! SUBLIMATION LESS THAN DEPTH OF SNOWPACK
4067 ! SNOWPACK (ESD) REDUCED BY ESNOW2 (DEPTH OF SUBLIMATED SNOW)
4068 ! ----------------------------------------------------------------------
4069 ELSE
4070 ESD = ESD-ESNOW2
4071 ETP3 = ETP * LSUBC
4072 SEH = RCH * (T1- TH2)
4073 T14 = T1* T1
4074 T14 = T14* T14
4075 FLX3 = FDOWN - FLX1 - FLX2 - &
4076 ((SNCOVR*EMISSI_S)+EMISSI*(1-SNCOVR))*SIGMA*T14 - &
4077 SSOIL - SEH - ETANRG
4078 IF (FLX3 <= 0.0) FLX3 = 0.0
4079 ! ----------------------------------------------------------------------
4080 ! SNOWMELT REDUCTION DEPENDING ON SNOW COVER
4081 ! ----------------------------------------------------------------------
4082 EX = FLX3*0.001/ LSUBF
4083
4084 ! ----------------------------------------------------------------------
4085 ! ESDMIN REPRESENTS A SNOWPACK DEPTH THRESHOLD VALUE BELOW WHICH WE
4086 ! CHOOSE NOT TO RETAIN ANY SNOWPACK, AND INSTEAD INCLUDE IT IN SNOWMELT.
4087 ! ----------------------------------------------------------------------
4088 SNOMLT = EX * DT
4089 IF (ESD- SNOMLT >= ESDMIN) THEN
4090 ESD = ESD- SNOMLT
4091 ! ----------------------------------------------------------------------
4092 ! SNOWMELT EXCEEDS SNOW DEPTH
4093 ! ----------------------------------------------------------------------
4094 ELSE
4095 EX = ESD / DT
4096 FLX3 = EX *1000.0* LSUBF
4097 SNOMLT = ESD
4098
4099 ESD = 0.0
4100 ! ----------------------------------------------------------------------
4101 ! END OF 'ESD .LE. ETP2' IF-BLOCK
4102 ! ----------------------------------------------------------------------
4103 END IF
4104 END IF
4105
4106 ! ----------------------------------------------------------------------
4107 ! END OF 'T12 .LE. TFREEZ' IF-BLOCK
4108 ! ----------------------------------------------------------------------
4109 PRCP1 = PRCP1+ EX
4110
4111 ! ----------------------------------------------------------------------
4112 ! SET THE EFFECTIVE POTNL EVAPOTRANSP (ETP1) TO ZERO SINCE THIS IS SNOW
4113 ! CASE, SO SURFACE EVAP NOT CALCULATED FROM EDIR, EC, OR ETT IN SMFLX
4114 ! (BELOW).
4115 ! IF SEAICE (ICE=1) SKIP CALL TO SMFLX.
4116 ! SMFLX RETURNS UPDATED SOIL MOISTURE VALUES. IN THIS, THE SNOW PACK
4117 ! CASE, ETA1 IS NOT USED IN CALCULATION OF EVAP.
4118 ! ----------------------------------------------------------------------
4119 END IF
4120 IF (ICE /= 1) THEN
4121 CALL SMFLX (SMC,NSOIL,CMC,DT,PRCP1,ZSOIL, &
4122 SH2O,SLOPE,KDT,FRZFACT, &
4123 SMCMAX,BEXP,SMCWLT,DKSAT,DWSAT, &
4124 SHDFAC,CMCMAX, &
4125 RUNOFF1,RUNOFF2,RUNOFF3, &
4126 EDIR1,EC1,ET1, &
4127 DRIP)
4128 ! ----------------------------------------------------------------------
4129 ! BEFORE CALL SHFLX IN THIS SNOWPACK CASE, SET ZZ1 AND YY ARGUMENTS TO
4130 ! SPECIAL VALUES THAT ENSURE THAT GROUND HEAT FLUX CALCULATED IN SHFLX
4131 ! MATCHES THAT ALREADY COMPUTER FOR BELOW THE SNOWPACK, THUS THE SFC
4132 ! HEAT FLUX TO BE COMPUTED IN SHFLX WILL EFFECTIVELY BE THE FLUX AT THE
4133 ! SNOW TOP SURFACE. T11 IS A DUMMY ARGUEMENT SO WE WILL NOT USE THE
4134 ! SKIN TEMP VALUE AS REVISED BY SHFLX.
4135 ! ----------------------------------------------------------------------
4136 END IF
4137 ZZ1 = 1.0
4138 YY = STC (1) -0.5* SSOIL * ZSOIL (1)* ZZ1/ DF1
4139
4140 ! ----------------------------------------------------------------------
4141 ! SHFLX WILL CALC/UPDATE THE SOIL TEMPS. NOTE: THE SUB-SFC HEAT FLUX
4142 ! (SSOIL1) AND THE SKIN TEMP (T11) OUTPUT FROM THIS SHFLX CALL ARE NOT
4143 ! USED IN ANY SUBSEQUENT CALCULATIONS. RATHER, THEY ARE DUMMY VARIABLES
4144 ! HERE IN THE SNOPAC CASE, SINCE THE SKIN TEMP AND SUB-SFC HEAT FLUX ARE
4145 ! UPDATED INSTEAD NEAR THE BEGINNING OF THE CALL TO SNOPAC.
4146 ! ----------------------------------------------------------------------
4147 T11 = T1
4148 CALL SHFLX (SSOIL1,STC,SMC,SMCMAX,NSOIL,T11,DT,YY,ZZ1,ZSOIL, &
4149 TBOT,ZBOT,SMCWLT,PSISAT,SH2O,BEXP,F1,DF1,ICE, &
4150 QUARTZ,CSOIL,VEGTYP)
4151
4152 ! ----------------------------------------------------------------------
4153 ! SNOW DEPTH AND DENSITY ADJUSTMENT BASED ON SNOW COMPACTION. YY IS
4154 ! ASSUMED TO BE THE SOIL TEMPERTURE AT THE TOP OF THE SOIL COLUMN.
4155 ! ----------------------------------------------------------------------
4156 IF (ESD > 0.) THEN
4157 CALL SNOWPACK (ESD,DT,SNOWH,SNDENS,T1,YY)
4158 ELSE
4159 ESD = 0.
4160 SNOWH = 0.
4161 SNDENS = 0.
4162 SNCOND = 1.
4163 SNCOVR = 0.
4164 END IF
4165 ! ----------------------------------------------------------------------
4166 END SUBROUTINE SNOPAC
4167 ! ----------------------------------------------------------------------
4168
4169
4170 SUBROUTINE SNOWPACK (ESD,DTSEC,SNOWH,SNDENS,TSNOW,TSOIL)
4171
4172 ! ----------------------------------------------------------------------
4173 ! SUBROUTINE SNOWPACK
4174 ! ----------------------------------------------------------------------
4175 ! CALCULATE COMPACTION OF SNOWPACK UNDER CONDITIONS OF INCREASING SNOW
4176 ! DENSITY, AS OBTAINED FROM AN APPROXIMATE SOLUTION OF E. ANDERSON'S
4177 ! DIFFERENTIAL EQUATION (3.29), NOAA TECHNICAL REPORT NWS 19, BY VICTOR
4178 ! KOREN, 03/25/95.
4179 ! ----------------------------------------------------------------------
4180 ! ESD WATER EQUIVALENT OF SNOW (M)
4181 ! DTSEC TIME STEP (SEC)
4182 ! SNOWH SNOW DEPTH (M)
4183 ! SNDENS SNOW DENSITY (G/CM3=DIMENSIONLESS FRACTION OF H2O DENSITY)
4184 ! TSNOW SNOW SURFACE TEMPERATURE (K)
4185 ! TSOIL SOIL SURFACE TEMPERATURE (K)
4186
4187 ! SUBROUTINE WILL RETURN NEW VALUES OF SNOWH AND SNDENS
4188 ! ----------------------------------------------------------------------
4189 IMPLICIT NONE
4190
4191 INTEGER :: IPOL, J
4192 REAL, INTENT(IN) :: ESD, DTSEC,TSNOW,TSOIL
4193 REAL, INTENT(INOUT) :: SNOWH, SNDENS
4194 REAL :: BFAC,DSX,DTHR,DW,SNOWHC,PEXP, &
4195 TAVGC,TSNOWC,TSOILC,ESDC,ESDCX
4196 REAL, PARAMETER :: C1 = 0.01, C2 = 21.0, G = 9.81, &
4197 KN = 4000.0
4198 ! ----------------------------------------------------------------------
4199 ! CONVERSION INTO SIMULATION UNITS
4200 ! ----------------------------------------------------------------------
4201 SNOWHC = SNOWH *100.
4202 ESDC = ESD *100.
4203 DTHR = DTSEC /3600.
4204 TSNOWC = TSNOW -273.15
4205 TSOILC = TSOIL -273.15
4206
4207 ! ----------------------------------------------------------------------
4208 ! CALCULATING OF AVERAGE TEMPERATURE OF SNOW PACK
4209 ! ----------------------------------------------------------------------
4210 ! ----------------------------------------------------------------------
4211 ! CALCULATING OF SNOW DEPTH AND DENSITY AS A RESULT OF COMPACTION
4212 ! SNDENS=DS0*(EXP(BFAC*ESD)-1.)/(BFAC*ESD)
4213 ! BFAC=DTHR*C1*EXP(0.08*TAVGC-C2*DS0)
4214 ! NOTE: BFAC*ESD IN SNDENS EQN ABOVE HAS TO BE CAREFULLY TREATED
4215 ! NUMERICALLY BELOW:
4216 ! C1 IS THE FRACTIONAL INCREASE IN DENSITY (1/(CM*HR))
4217 ! C2 IS A CONSTANT (CM3/G) KOJIMA ESTIMATED AS 21 CMS/G
4218 ! ----------------------------------------------------------------------
4219 TAVGC = 0.5* (TSNOWC + TSOILC)
4220 IF (ESDC > 1.E-2) THEN
4221 ESDCX = ESDC
4222 ELSE
4223 ESDCX = 1.E-2
4224 END IF
4225
4226 ! DSX = SNDENS*((DEXP(BFAC*ESDC)-1.)/(BFAC*ESDC))
4227 ! ----------------------------------------------------------------------
4228 ! THE FUNCTION OF THE FORM (e**x-1)/x IMBEDDED IN ABOVE EXPRESSION
4229 ! FOR DSX WAS CAUSING NUMERICAL DIFFICULTIES WHEN THE DENOMINATOR "x"
4230 ! (I.E. BFAC*ESDC) BECAME ZERO OR APPROACHED ZERO (DESPITE THE FACT THAT
4231 ! THE ANALYTICAL FUNCTION (e**x-1)/x HAS A WELL DEFINED LIMIT AS
4232 ! "x" APPROACHES ZERO), HENCE BELOW WE REPLACE THE (e**x-1)/x
4233 ! EXPRESSION WITH AN EQUIVALENT, NUMERICALLY WELL-BEHAVED
4234 ! POLYNOMIAL EXPANSION.
4235
4236 ! NUMBER OF TERMS OF POLYNOMIAL EXPANSION, AND HENCE ITS ACCURACY,
4237 ! IS GOVERNED BY ITERATION LIMIT "IPOL".
4238 ! IPOL GREATER THAN 9 ONLY MAKES A DIFFERENCE ON DOUBLE
4239 ! PRECISION (RELATIVE ERRORS GIVEN IN PERCENT %).
4240 ! IPOL=9, FOR REL.ERROR <~ 1.6 E-6 % (8 SIGNIFICANT DIGITS)
4241 ! IPOL=8, FOR REL.ERROR <~ 1.8 E-5 % (7 SIGNIFICANT DIGITS)
4242 ! IPOL=7, FOR REL.ERROR <~ 1.8 E-4 % ...
4243 ! ----------------------------------------------------------------------
4244 BFAC = DTHR * C1* EXP (0.08* TAVGC - C2* SNDENS)
4245 IPOL = 4
4246 PEXP = 0.
4247 ! PEXP = (1. + PEXP)*BFAC*ESDC/REAL(J+1)
4248 DO J = IPOL,1, -1
4249 PEXP = (1. + PEXP)* BFAC * ESDCX / REAL (J +1)
4250 END DO
4251
4252 PEXP = PEXP + 1.
4253 ! ----------------------------------------------------------------------
4254 ! ABOVE LINE ENDS POLYNOMIAL SUBSTITUTION
4255 ! ----------------------------------------------------------------------
4256 ! END OF KOREAN FORMULATION
4257
4258 ! BASE FORMULATION (COGLEY ET AL., 1990)
4259 ! CONVERT DENSITY FROM G/CM3 TO KG/M3
4260 ! DSM=SNDENS*1000.0
4261
4262 ! DSX=DSM+DTSEC*0.5*DSM*G*ESD/
4263 ! & (1E7*EXP(-0.02*DSM+KN/(TAVGC+273.16)-14.643))
4264
4265 ! & CONVERT DENSITY FROM KG/M3 TO G/CM3
4266 ! DSX=DSX/1000.0
4267
4268 ! END OF COGLEY ET AL. FORMULATION
4269
4270 ! ----------------------------------------------------------------------
4271 ! SET UPPER/LOWER LIMIT ON SNOW DENSITY
4272 ! ----------------------------------------------------------------------
4273 DSX = SNDENS * (PEXP)
4274 IF (DSX > 0.40) DSX = 0.40
4275 IF (DSX < 0.05) DSX = 0.05
4276 ! ----------------------------------------------------------------------
4277 ! UPDATE OF SNOW DEPTH AND DENSITY DEPENDING ON LIQUID WATER DURING
4278 ! SNOWMELT. ASSUMED THAT 13% OF LIQUID WATER CAN BE STORED IN SNOW PER
4279 ! DAY DURING SNOWMELT TILL SNOW DENSITY 0.40.
4280 ! ----------------------------------------------------------------------
4281 SNDENS = DSX
4282 IF (TSNOWC >= 0.) THEN
4283 DW = 0.13* DTHR /24.
4284 SNDENS = SNDENS * (1. - DW) + DW
4285 IF (SNDENS >= 0.40) SNDENS = 0.40
4286 ! ----------------------------------------------------------------------
4287 ! CALCULATE SNOW DEPTH (CM) FROM SNOW WATER EQUIVALENT AND SNOW DENSITY.
4288 ! CHANGE SNOW DEPTH UNITS TO METERS
4289 ! ----------------------------------------------------------------------
4290 END IF
4291 SNOWHC = ESDC / SNDENS
4292 SNOWH = SNOWHC *0.01
4293
4294 ! ----------------------------------------------------------------------
4295 END SUBROUTINE SNOWPACK
4296 ! ----------------------------------------------------------------------
4297
4298 SUBROUTINE SNOWZ0 (SNCOVR,Z0, Z0BRD)
4299
4300 ! ----------------------------------------------------------------------
4301 ! SUBROUTINE SNOWZ0
4302 ! ----------------------------------------------------------------------
4303 ! CALCULATE TOTAL ROUGHNESS LENGTH OVER SNOW
4304 ! SNCOVR FRACTIONAL SNOW COVER
4305 ! Z0 ROUGHNESS LENGTH (m)
4306 ! Z0S SNOW ROUGHNESS LENGTH:=0.001 (m)
4307 ! ----------------------------------------------------------------------
4308 IMPLICIT NONE
4309 REAL, INTENT(IN) :: SNCOVR, Z0BRD
4310 REAL, INTENT(OUT) :: Z0
4311 REAL, PARAMETER :: Z0S=0.001
4312
4313 !m Z0 = (1.- SNCOVR)* Z0BRD + SNCOVR * Z0S
4314 Z0 = Z0BRD
4315 ! ----------------------------------------------------------------------
4316 END SUBROUTINE SNOWZ0
4317 ! ----------------------------------------------------------------------
4318
4319
4320 SUBROUTINE SNOW_NEW (TEMP,NEWSN,SNOWH,SNDENS)
4321
4322 ! ----------------------------------------------------------------------
4323 ! SUBROUTINE SNOW_NEW
4324 ! ----------------------------------------------------------------------
4325 ! CALCULATE SNOW DEPTH AND DENSITITY TO ACCOUNT FOR THE NEW SNOWFALL.
4326 ! NEW VALUES OF SNOW DEPTH & DENSITY RETURNED.
4327
4328 ! TEMP AIR TEMPERATURE (K)
4329 ! NEWSN NEW SNOWFALL (M)
4330 ! SNOWH SNOW DEPTH (M)
4331 ! SNDENS SNOW DENSITY (G/CM3=DIMENSIONLESS FRACTION OF H2O DENSITY)
4332 ! ----------------------------------------------------------------------
4333 IMPLICIT NONE
4334 REAL, INTENT(IN) :: NEWSN, TEMP
4335 REAL, INTENT(INOUT) :: SNDENS, SNOWH
4336 REAL :: DSNEW, HNEWC, SNOWHC,NEWSNC,TEMPC
4337
4338 ! ----------------------------------------------------------------------
4339 ! CONVERSION INTO SIMULATION UNITS
4340 ! ----------------------------------------------------------------------
4341 SNOWHC = SNOWH *100.
4342 NEWSNC = NEWSN *100.
4343
4344 ! ----------------------------------------------------------------------
4345 ! CALCULATING NEW SNOWFALL DENSITY DEPENDING ON TEMPERATURE
4346 ! EQUATION FROM GOTTLIB L. 'A GENERAL RUNOFF MODEL FOR SNOWCOVERED
4347 ! AND GLACIERIZED BASIN', 6TH NORDIC HYDROLOGICAL CONFERENCE,
4348 ! VEMADOLEN, SWEDEN, 1980, 172-177PP.
4349 !-----------------------------------------------------------------------
4350 TEMPC = TEMP -273.15
4351 IF (TEMPC <= -15.) THEN
4352 DSNEW = 0.05
4353 ELSE
4354 DSNEW = 0.05+0.0017* (TEMPC +15.)**1.5
4355 END IF
4356 ! ----------------------------------------------------------------------
4357 ! ADJUSTMENT OF SNOW DENSITY DEPENDING ON NEW SNOWFALL
4358 ! ----------------------------------------------------------------------
4359 HNEWC = NEWSNC / DSNEW
4360 IF (SNOWHC + HNEWC .LT. 1.0E-3) THEN
4361 SNDENS = MAX(DSNEW,SNDENS)
4362 ELSE
4363 SNDENS = (SNOWHC * SNDENS + HNEWC * DSNEW)/ (SNOWHC + HNEWC)
4364 ENDIF
4365 SNOWHC = SNOWHC + HNEWC
4366 SNOWH = SNOWHC *0.01
4367
4368 ! ----------------------------------------------------------------------
4369 END SUBROUTINE SNOW_NEW
4370 ! ----------------------------------------------------------------------
4371
4372 SUBROUTINE SRT (RHSTT,EDIR,ET,SH2O,SH2OA,NSOIL,PCPDRP, &
4373 ZSOIL,DWSAT,DKSAT,SMCMAX,BEXP,RUNOFF1, &
4374 RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZX,SICE,AI,BI,CI)
4375
4376 ! ----------------------------------------------------------------------
4377 ! SUBROUTINE SRT
4378 ! ----------------------------------------------------------------------
4379 ! CALCULATE THE RIGHT HAND SIDE OF THE TIME TENDENCY TERM OF THE SOIL
4380 ! WATER DIFFUSION EQUATION. ALSO TO COMPUTE ( PREPARE ) THE MATRIX
4381 ! COEFFICIENTS FOR THE TRI-DIAGONAL MATRIX OF THE IMPLICIT TIME SCHEME.
4382 ! ----------------------------------------------------------------------
4383 IMPLICIT NONE
4384 INTEGER, INTENT(IN) :: NSOIL
4385 INTEGER :: IALP1, IOHINF, J, JJ, K, KS
4386 REAL, INTENT(IN) :: BEXP, DKSAT, DT, DWSAT, EDIR, FRZX, &
4387 KDT, PCPDRP, SLOPE, SMCMAX, SMCWLT
4388 REAL, INTENT(OUT) :: RUNOFF1, RUNOFF2
4389 REAL, DIMENSION(1:NSOIL), INTENT(IN) :: ET, SH2O, SH2OA, SICE, &
4390 ZSOIL
4391 REAL, DIMENSION(1:NSOIL), INTENT(OUT) :: RHSTT
4392 REAL, DIMENSION(1:NSOIL), INTENT(OUT) :: AI, BI, CI
4393 REAL, DIMENSION(1:NSOIL) :: DMAX
4394 REAL :: ACRT, DD, DDT, DDZ, DDZ2, DENOM, &
4395 DENOM2,DICE, DSMDZ, DSMDZ2, DT1, &
4396 FCR,INFMAX,MXSMC,MXSMC2,NUMER,PDDUM, &
4397 PX, SICEMAX,SLOPX, SMCAV, SSTT, &
4398 SUM, VAL, WCND, WCND2, WDF, WDF2
4399 INTEGER, PARAMETER :: CVFRZ = 3
4400
4401 ! ----------------------------------------------------------------------
4402 ! FROZEN GROUND VERSION:
4403 ! REFERENCE FROZEN GROUND PARAMETER, CVFRZ, IS A SHAPE PARAMETER OF
4404 ! AREAL DISTRIBUTION FUNCTION OF SOIL ICE CONTENT WHICH EQUALS 1/CV.
4405 ! CV IS A COEFFICIENT OF SPATIAL VARIATION OF SOIL ICE CONTENT. BASED
4406 ! ON FIELD DATA CV DEPENDS ON AREAL MEAN OF FROZEN DEPTH, AND IT CLOSE
4407 ! TO CONSTANT = 0.6 IF AREAL MEAN FROZEN DEPTH IS ABOVE 20 CM. THAT IS
4408 ! WHY PARAMETER CVFRZ = 3 (INT{1/0.6*0.6}).
4409 ! CURRENT LOGIC DOESN'T ALLOW CVFRZ BE BIGGER THAN 3
4410 ! ----------------------------------------------------------------------
4411
4412 ! ----------------------------------------------------------------------
4413 ! DETERMINE RAINFALL INFILTRATION RATE AND RUNOFF. INCLUDE THE
4414 ! INFILTRATION FORMULE FROM SCHAAKE AND KOREN MODEL.
4415 ! MODIFIED BY Q DUAN
4416 ! ----------------------------------------------------------------------
4417 ! ----------------------------------------------------------------------
4418 ! LET SICEMAX BE THE GREATEST, IF ANY, FROZEN WATER CONTENT WITHIN SOIL
4419 ! LAYERS.
4420 ! ----------------------------------------------------------------------
4421 IOHINF = 1
4422 SICEMAX = 0.0
4423 DO KS = 1,NSOIL
4424 IF (SICE (KS) > SICEMAX) SICEMAX = SICE (KS)
4425 ! ----------------------------------------------------------------------
4426 ! DETERMINE RAINFALL INFILTRATION RATE AND RUNOFF
4427 ! ----------------------------------------------------------------------
4428 END DO
4429 PDDUM = PCPDRP
4430 RUNOFF1 = 0.0
4431
4432 ! ----------------------------------------------------------------------
4433 ! MODIFIED BY Q. DUAN, 5/16/94
4434 ! ----------------------------------------------------------------------
4435 ! IF (IOHINF == 1) THEN
4436
4437 IF (PCPDRP /= 0.0) THEN
4438 DT1 = DT /86400.
4439 SMCAV = SMCMAX - SMCWLT
4440
4441 ! ----------------------------------------------------------------------
4442 ! FROZEN GROUND VERSION:
4443 ! ----------------------------------------------------------------------
4444 DMAX (1)= - ZSOIL (1)* SMCAV
4445
4446 DICE = - ZSOIL (1) * SICE (1)
4447 DMAX (1)= DMAX (1)* (1.0- (SH2OA (1) + SICE (1) - SMCWLT)/ &
4448 SMCAV)
4449
4450 DD = DMAX (1)
4451
4452 ! ----------------------------------------------------------------------
4453 ! FROZEN GROUND VERSION:
4454 ! ----------------------------------------------------------------------
4455 DO KS = 2,NSOIL
4456
4457 DICE = DICE+ ( ZSOIL (KS -1) - ZSOIL (KS) ) * SICE (KS)
4458 DMAX (KS) = (ZSOIL (KS -1) - ZSOIL (KS))* SMCAV
4459 DMAX (KS) = DMAX (KS)* (1.0- (SH2OA (KS) + SICE (KS) &
4460 - SMCWLT)/ SMCAV)
4461 DD = DD+ DMAX (KS)
4462 ! ----------------------------------------------------------------------
4463 ! VAL = (1.-EXP(-KDT*SQRT(DT1)))
4464 ! IN BELOW, REMOVE THE SQRT IN ABOVE
4465 ! ----------------------------------------------------------------------
4466 END DO
4467 VAL = (1. - EXP ( - KDT * DT1))
4468 DDT = DD * VAL
4469 PX = PCPDRP * DT
4470 IF (PX < 0.0) PX = 0.0
4471
4472 ! ----------------------------------------------------------------------
4473 ! FROZEN GROUND VERSION:
4474 ! REDUCTION OF INFILTRATION BASED ON FROZEN GROUND PARAMETERS
4475 ! ----------------------------------------------------------------------
4476 INFMAX = (PX * (DDT / (PX + DDT)))/ DT
4477 FCR = 1.
4478 IF (DICE > 1.E-2) THEN
4479 ACRT = CVFRZ * FRZX / DICE
4480 SUM = 1.
4481 IALP1 = CVFRZ - 1
4482 DO J = 1,IALP1
4483 K = 1
4484 DO JJ = J +1,IALP1
4485 K = K * JJ
4486 END DO
4487 SUM = SUM + (ACRT ** ( CVFRZ - J)) / FLOAT (K)
4488 END DO
4489 FCR = 1. - EXP ( - ACRT) * SUM
4490 END IF
4491
4492 ! ----------------------------------------------------------------------
4493 ! CORRECTION OF INFILTRATION LIMITATION:
4494 ! IF INFMAX .LE. HYDROLIC CONDUCTIVITY ASSIGN INFMAX THE VALUE OF
4495 ! HYDROLIC CONDUCTIVITY
4496 ! ----------------------------------------------------------------------
4497 ! MXSMC = MAX ( SH2OA(1), SH2OA(2) )
4498 INFMAX = INFMAX * FCR
4499
4500 MXSMC = SH2OA (1)
4501 CALL WDFCND (WDF,WCND,MXSMC,SMCMAX,BEXP,DKSAT,DWSAT, &
4502 SICEMAX)
4503 INFMAX = MAX (INFMAX,WCND)
4504
4505 INFMAX = MIN (INFMAX,PX)
4506 IF (PCPDRP > INFMAX) THEN
4507 RUNOFF1 = PCPDRP - INFMAX
4508 PDDUM = INFMAX
4509 END IF
4510 ! ----------------------------------------------------------------------
4511 ! TO AVOID SPURIOUS DRAINAGE BEHAVIOR, 'UPSTREAM DIFFERENCING' IN LINE
4512 ! BELOW REPLACED WITH NEW APPROACH IN 2ND LINE:
4513 ! 'MXSMC = MAX(SH2OA(1), SH2OA(2))'
4514 ! ----------------------------------------------------------------------
4515 END IF
4516
4517 MXSMC = SH2OA (1)
4518 CALL WDFCND (WDF,WCND,MXSMC,SMCMAX,BEXP,DKSAT,DWSAT, &
4519 SICEMAX)
4520 ! ----------------------------------------------------------------------
4521 ! CALC THE MATRIX COEFFICIENTS AI, BI, AND CI FOR THE TOP LAYER
4522 ! ----------------------------------------------------------------------
4523 DDZ = 1. / ( - .5 * ZSOIL (2) )
4524 AI (1) = 0.0
4525 BI (1) = WDF * DDZ / ( - ZSOIL (1) )
4526
4527 ! ----------------------------------------------------------------------
4528 ! CALC RHSTT FOR THE TOP LAYER AFTER CALC'NG THE VERTICAL SOIL MOISTURE
4529 ! GRADIENT BTWN THE TOP AND NEXT TO TOP LAYERS.
4530 ! ----------------------------------------------------------------------
4531 CI (1) = - BI (1)
4532 DSMDZ = ( SH2O (1) - SH2O (2) ) / ( - .5 * ZSOIL (2) )
4533 RHSTT (1) = (WDF * DSMDZ + WCND- PDDUM + EDIR + ET (1))/ ZSOIL (1)
4534
4535 ! ----------------------------------------------------------------------
4536 ! INITIALIZE DDZ2
4537 ! ----------------------------------------------------------------------
4538 SSTT = WDF * DSMDZ + WCND+ EDIR + ET (1)
4539
4540 ! ----------------------------------------------------------------------
4541 ! LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABV PROCESS
4542 ! ----------------------------------------------------------------------
4543 DDZ2 = 0.0
4544 DO K = 2,NSOIL
4545 DENOM2 = (ZSOIL (K -1) - ZSOIL (K))
4546 IF (K /= NSOIL) THEN
4547
4548 ! ----------------------------------------------------------------------
4549 ! AGAIN, TO AVOID SPURIOUS DRAINAGE BEHAVIOR, 'UPSTREAM DIFFERENCING' IN
4550 ! LINE BELOW REPLACED WITH NEW APPROACH IN 2ND LINE:
4551 ! 'MXSMC2 = MAX (SH2OA(K), SH2OA(K+1))'
4552 ! ----------------------------------------------------------------------
4553 SLOPX = 1.
4554
4555 MXSMC2 = SH2OA (K)
4556 CALL WDFCND (WDF2,WCND2,MXSMC2,SMCMAX,BEXP,DKSAT,DWSAT, &
4557 SICEMAX)
4558 ! -----------------------------------------------------------------------
4559 ! CALC SOME PARTIAL PRODUCTS FOR LATER USE IN CALC'NG RHSTT
4560 ! ----------------------------------------------------------------------
4561 DENOM = (ZSOIL (K -1) - ZSOIL (K +1))
4562
4563 ! ----------------------------------------------------------------------
4564 ! CALC THE MATRIX COEF, CI, AFTER CALC'NG ITS PARTIAL PRODUCT
4565 ! ----------------------------------------------------------------------
4566 DSMDZ2 = (SH2O (K) - SH2O (K +1)) / (DENOM * 0.5)
4567 DDZ2 = 2.0 / DENOM
4568 CI (K) = - WDF2 * DDZ2 / DENOM2
4569
4570 ELSE
4571 ! ----------------------------------------------------------------------
4572 ! SLOPE OF BOTTOM LAYER IS INTRODUCED
4573 ! ----------------------------------------------------------------------
4574
4575 ! ----------------------------------------------------------------------
4576 ! RETRIEVE THE SOIL WATER DIFFUSIVITY AND HYDRAULIC CONDUCTIVITY FOR
4577 ! THIS LAYER
4578 ! ----------------------------------------------------------------------
4579 SLOPX = SLOPE
4580 CALL WDFCND (WDF2,WCND2,SH2OA (NSOIL),SMCMAX,BEXP,DKSAT,DWSAT, &
4581 SICEMAX)
4582
4583 ! ----------------------------------------------------------------------
4584 ! CALC A PARTIAL PRODUCT FOR LATER USE IN CALC'NG RHSTT
4585 ! ----------------------------------------------------------------------
4586
4587 ! ----------------------------------------------------------------------
4588 ! SET MATRIX COEF CI TO ZERO
4589 ! ----------------------------------------------------------------------
4590 DSMDZ2 = 0.0
4591 CI (K) = 0.0
4592 ! ----------------------------------------------------------------------
4593 ! CALC RHSTT FOR THIS LAYER AFTER CALC'NG ITS NUMERATOR
4594 ! ----------------------------------------------------------------------
4595 END IF
4596 NUMER = (WDF2 * DSMDZ2) + SLOPX * WCND2- (WDF * DSMDZ) &
4597 - WCND+ ET (K)
4598
4599 ! ----------------------------------------------------------------------
4600 ! CALC MATRIX COEFS, AI, AND BI FOR THIS LAYER
4601 ! ----------------------------------------------------------------------
4602 RHSTT (K) = NUMER / ( - DENOM2)
4603 AI (K) = - WDF * DDZ / DENOM2
4604
4605 ! ----------------------------------------------------------------------
4606 ! RESET VALUES OF WDF, WCND, DSMDZ, AND DDZ FOR LOOP TO NEXT LYR
4607 ! RUNOFF2: SUB-SURFACE OR BASEFLOW RUNOFF
4608 ! ----------------------------------------------------------------------
4609 BI (K) = - ( AI (K) + CI (K) )
4610 IF (K .eq. NSOIL) THEN
4611 RUNOFF2 = SLOPX * WCND2
4612 END IF
4613 IF (K .ne. NSOIL) THEN
4614 WDF = WDF2
4615 WCND = WCND2
4616 DSMDZ = DSMDZ2
4617 DDZ = DDZ2
4618 END IF
4619 END DO
4620 ! ----------------------------------------------------------------------
4621 END SUBROUTINE SRT
4622 ! ----------------------------------------------------------------------
4623
4624 SUBROUTINE SSTEP (SH2OOUT,SH2OIN,CMC,RHSTT,RHSCT,DT, &
4625 NSOIL,SMCMAX,CMCMAX,RUNOFF3,ZSOIL,SMC,SICE, &
4626 AI,BI,CI)
4627
4628 ! ----------------------------------------------------------------------
4629 ! SUBROUTINE SSTEP
4630 ! ----------------------------------------------------------------------
4631 ! CALCULATE/UPDATE SOIL MOISTURE CONTENT VALUES AND CANOPY MOISTURE
4632 ! CONTENT VALUES.
4633 ! ----------------------------------------------------------------------
4634 IMPLICIT NONE
4635 INTEGER, INTENT(IN) :: NSOIL
4636 INTEGER :: I, K, KK11
4637
4638 REAL, INTENT(IN) :: CMCMAX, DT, SMCMAX
4639 REAL, INTENT(OUT) :: RUNOFF3
4640 REAL, INTENT(INOUT) :: CMC
4641 REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SH2OIN, SICE, ZSOIL
4642 REAL, DIMENSION(1:NSOIL), INTENT(OUT) :: SH2OOUT
4643 REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: RHSTT, SMC
4644 REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: AI, BI, CI
4645 REAL, DIMENSION(1:NSOIL) :: RHSTTin
4646 REAL, DIMENSION(1:NSOIL) :: CIin
4647 REAL :: DDZ, RHSCT, STOT, WPLUS
4648
4649 ! ----------------------------------------------------------------------
4650 ! CREATE 'AMOUNT' VALUES OF VARIABLES TO BE INPUT TO THE
4651 ! TRI-DIAGONAL MATRIX ROUTINE.
4652 ! ----------------------------------------------------------------------
4653 DO K = 1,NSOIL
4654 RHSTT (K) = RHSTT (K) * DT
4655 AI (K) = AI (K) * DT
4656 BI (K) = 1. + BI (K) * DT
4657 CI (K) = CI (K) * DT
4658 END DO
4659 ! ----------------------------------------------------------------------
4660 ! COPY VALUES FOR INPUT VARIABLES BEFORE CALL TO ROSR12
4661 ! ----------------------------------------------------------------------
4662 DO K = 1,NSOIL
4663 RHSTTin (K) = RHSTT (K)
4664 END DO
4665 DO K = 1,NSOIL
4666 CIin (K) = CI (K)
4667 END DO
4668 ! ----------------------------------------------------------------------
4669 ! CALL ROSR12 TO SOLVE THE TRI-DIAGONAL MATRIX
4670 ! ----------------------------------------------------------------------
4671 CALL ROSR12 (CI,AI,BI,CIin,RHSTTin,RHSTT,NSOIL)
4672 ! ----------------------------------------------------------------------
4673 ! SUM THE PREVIOUS SMC VALUE AND THE MATRIX SOLUTION TO GET A
4674 ! NEW VALUE. MIN ALLOWABLE VALUE OF SMC WILL BE 0.02.
4675 ! RUNOFF3: RUNOFF WITHIN SOIL LAYERS
4676 ! ----------------------------------------------------------------------
4677 WPLUS = 0.0
4678 RUNOFF3 = 0.
4679
4680 DDZ = - ZSOIL (1)
4681 DO K = 1,NSOIL
4682 IF (K /= 1) DDZ = ZSOIL (K - 1) - ZSOIL (K)
4683 SH2OOUT (K) = SH2OIN (K) + CI (K) + WPLUS / DDZ
4684 STOT = SH2OOUT (K) + SICE (K)
4685 IF (STOT > SMCMAX) THEN
4686 IF (K .eq. 1) THEN
4687 DDZ = - ZSOIL (1)
4688 ELSE
4689 KK11 = K - 1
4690 DDZ = - ZSOIL (K) + ZSOIL (KK11)
4691 END IF
4692 WPLUS = (STOT - SMCMAX) * DDZ
4693 ELSE
4694 WPLUS = 0.
4695 END IF
4696 SMC (K) = MAX ( MIN (STOT,SMCMAX),0.02 )
4697 SH2OOUT (K) = MAX ( (SMC (K) - SICE (K)),0.0)
4698 END DO
4699
4700 ! ----------------------------------------------------------------------
4701 ! UPDATE CANOPY WATER CONTENT/INTERCEPTION (CMC). CONVERT RHSCT TO
4702 ! AN 'AMOUNT' VALUE AND ADD TO PREVIOUS CMC VALUE TO GET NEW CMC.
4703 ! ----------------------------------------------------------------------
4704 RUNOFF3 = WPLUS
4705 CMC = CMC + DT * RHSCT
4706 IF (CMC < 1.E-20) CMC = 0.0
4707 CMC = MIN (CMC,CMCMAX)
4708
4709 ! ----------------------------------------------------------------------
4710 END SUBROUTINE SSTEP
4711 ! ----------------------------------------------------------------------
4712
4713 SUBROUTINE TBND (TU,TB,ZSOIL,ZBOT,K,NSOIL,TBND1)
4714
4715 ! ----------------------------------------------------------------------
4716 ! SUBROUTINE TBND
4717 ! ----------------------------------------------------------------------
4718 ! CALCULATE TEMPERATURE ON THE BOUNDARY OF THE LAYER BY INTERPOLATION OF
4719 ! THE MIDDLE LAYER TEMPERATURES
4720 ! ----------------------------------------------------------------------
4721 IMPLICIT NONE
4722 INTEGER, INTENT(IN) :: NSOIL
4723 INTEGER :: K
4724 REAL, INTENT(IN) :: TB, TU, ZBOT
4725 REAL, INTENT(OUT) :: TBND1
4726 REAL, DIMENSION(1:NSOIL), INTENT(IN) :: ZSOIL
4727 REAL :: ZB, ZUP
4728 REAL, PARAMETER :: T0 = 273.15
4729
4730 ! ----------------------------------------------------------------------
4731 ! USE SURFACE TEMPERATURE ON THE TOP OF THE FIRST LAYER
4732 ! ----------------------------------------------------------------------
4733 IF (K == 1) THEN
4734 ZUP = 0.
4735 ELSE
4736 ZUP = ZSOIL (K -1)
4737 END IF
4738 ! ----------------------------------------------------------------------
4739 ! USE DEPTH OF THE CONSTANT BOTTOM TEMPERATURE WHEN INTERPOLATE
4740 ! TEMPERATURE INTO THE LAST LAYER BOUNDARY
4741 ! ----------------------------------------------------------------------
4742 IF (K == NSOIL) THEN
4743 ZB = 2.* ZBOT - ZSOIL (K)
4744 ELSE
4745 ZB = ZSOIL (K +1)
4746 END IF
4747 ! ----------------------------------------------------------------------
4748 ! LINEAR INTERPOLATION BETWEEN THE AVERAGE LAYER TEMPERATURES
4749 ! ----------------------------------------------------------------------
4750
4751 TBND1 = TU + (TB - TU)* (ZUP - ZSOIL (K))/ (ZUP - ZB)
4752 ! ----------------------------------------------------------------------
4753 END SUBROUTINE TBND
4754 ! ----------------------------------------------------------------------
4755
4756
4757 SUBROUTINE TDFCND ( DF, SMC, QZ, SMCMAX, SH2O)
4758
4759 ! ----------------------------------------------------------------------
4760 ! SUBROUTINE TDFCND
4761 ! ----------------------------------------------------------------------
4762 ! CALCULATE THERMAL DIFFUSIVITY AND CONDUCTIVITY OF THE SOIL FOR A GIVEN
4763 ! POINT AND TIME.
4764 ! ----------------------------------------------------------------------
4765 ! PETERS-LIDARD APPROACH (PETERS-LIDARD et al., 1998)
4766 ! June 2001 CHANGES: FROZEN SOIL CONDITION.
4767 ! ----------------------------------------------------------------------
4768 IMPLICIT NONE
4769 REAL, INTENT(IN) :: QZ, SMC, SMCMAX, SH2O
4770 REAL, INTENT(OUT) :: DF
4771 REAL :: AKE, GAMMD, THKDRY, THKICE, THKO, &
4772 THKQTZ,THKSAT,THKS,THKW,SATRATIO,XU, &
4773 XUNFROZ
4774
4775 ! ----------------------------------------------------------------------
4776 ! WE NOW GET QUARTZ AS AN INPUT ARGUMENT (SET IN ROUTINE REDPRM):
4777 ! DATA QUARTZ /0.82, 0.10, 0.25, 0.60, 0.52,
4778 ! & 0.35, 0.60, 0.40, 0.82/
4779 ! ----------------------------------------------------------------------
4780 ! IF THE SOIL HAS ANY MOISTURE CONTENT COMPUTE A PARTIAL SUM/PRODUCT
4781 ! OTHERWISE USE A CONSTANT VALUE WHICH WORKS WELL WITH MOST SOILS
4782 ! ----------------------------------------------------------------------
4783 ! THKW ......WATER THERMAL CONDUCTIVITY
4784 ! THKQTZ ....THERMAL CONDUCTIVITY FOR QUARTZ
4785 ! THKO ......THERMAL CONDUCTIVITY FOR OTHER SOIL COMPONENTS
4786 ! THKS ......THERMAL CONDUCTIVITY FOR THE SOLIDS COMBINED(QUARTZ+OTHER)
4787 ! THKICE ....ICE THERMAL CONDUCTIVITY
4788 ! SMCMAX ....POROSITY (= SMCMAX)
4789 ! QZ .........QUARTZ CONTENT (SOIL TYPE DEPENDENT)
4790 ! ----------------------------------------------------------------------
4791 ! USE AS IN PETERS-LIDARD, 1998 (MODIF. FROM JOHANSEN, 1975).
4792
4793 ! PABLO GRUNMANN, 08/17/98
4794 ! REFS.:
4795 ! FAROUKI, O.T.,1986: THERMAL PROPERTIES OF SOILS. SERIES ON ROCK
4796 ! AND SOIL MECHANICS, VOL. 11, TRANS TECH, 136 PP.
4797 ! JOHANSEN, O., 1975: THERMAL CONDUCTIVITY OF SOILS. PH.D. THESIS,
4798 ! UNIVERSITY OF TRONDHEIM,
4799 ! PETERS-LIDARD, C. D., ET AL., 1998: THE EFFECT OF SOIL THERMAL
4800 ! CONDUCTIVITY PARAMETERIZATION ON SURFACE ENERGY FLUXES
4801 ! AND TEMPERATURES. JOURNAL OF THE ATMOSPHERIC SCIENCES,
4802 ! VOL. 55, PP. 1209-1224.
4803 ! ----------------------------------------------------------------------
4804 ! NEEDS PARAMETERS
4805 ! POROSITY(SOIL TYPE):
4806 ! POROS = SMCMAX
4807 ! SATURATION RATIO:
4808 ! PARAMETERS W/(M.K)
4809 SATRATIO = SMC / SMCMAX
4810 ! ICE CONDUCTIVITY:
4811 THKICE = 2.2
4812 ! WATER CONDUCTIVITY:
4813 THKW = 0.57
4814 ! THERMAL CONDUCTIVITY OF "OTHER" SOIL COMPONENTS
4815 ! IF (QZ .LE. 0.2) THKO = 3.0
4816 THKO = 2.0
4817 ! QUARTZ' CONDUCTIVITY
4818 THKQTZ = 7.7
4819 ! SOLIDS' CONDUCTIVITY
4820 THKS = (THKQTZ ** QZ)* (THKO ** (1. - QZ))
4821
4822 ! UNFROZEN FRACTION (FROM 1., i.e., 100%LIQUID, TO 0. (100% FROZEN))
4823 XUNFROZ = SH2O / SMC
4824 ! UNFROZEN VOLUME FOR SATURATION (POROSITY*XUNFROZ)
4825 XU = XUNFROZ * SMCMAX
4826
4827 ! SATURATED THERMAL CONDUCTIVITY
4828 THKSAT = THKS ** (1. - SMCMAX)* THKICE ** (SMCMAX - XU)* THKW ** &
4829 (XU)
4830
4831 ! DRY DENSITY IN KG/M3
4832 GAMMD = (1. - SMCMAX)*2700.
4833
4834 ! DRY THERMAL CONDUCTIVITY IN W.M-1.K-1
4835 THKDRY = (0.135* GAMMD+ 64.7)/ (2700. - 0.947* GAMMD)
4836 ! FROZEN
4837 IF ( (SH2O + 0.0005) < SMC ) THEN
4838 AKE = SATRATIO
4839 ! UNFROZEN
4840 ! RANGE OF VALIDITY FOR THE KERSTEN NUMBER (AKE)
4841 ELSE
4842
4843 ! KERSTEN NUMBER (USING "FINE" FORMULA, VALID FOR SOILS CONTAINING AT
4844 ! LEAST 5% OF PARTICLES WITH DIAMETER LESS THAN 2.E-6 METERS.)
4845 ! (FOR "COARSE" FORMULA, SEE PETERS-LIDARD ET AL., 1998).
4846
4847 IF ( SATRATIO > 0.1 ) THEN
4848
4849 AKE = LOG10 (SATRATIO) + 1.0
4850
4851 ! USE K = KDRY
4852 ELSE
4853
4854 AKE = 0.0
4855 END IF
4856 ! THERMAL CONDUCTIVITY
4857
4858 END IF
4859
4860 DF = AKE * (THKSAT - THKDRY) + THKDRY
4861 ! ----------------------------------------------------------------------
4862 END SUBROUTINE TDFCND
4863 ! ----------------------------------------------------------------------
4864
4865 SUBROUTINE TMPAVG (TAVG,TUP,TM,TDN,ZSOIL,NSOIL,K)
4866
4867 ! ----------------------------------------------------------------------
4868 ! SUBROUTINE TMPAVG
4869 ! ----------------------------------------------------------------------
4870 ! CALCULATE SOIL LAYER AVERAGE TEMPERATURE (TAVG) IN FREEZING/THAWING
4871 ! LAYER USING UP, DOWN, AND MIDDLE LAYER TEMPERATURES (TUP, TDN, TM),
4872 ! WHERE TUP IS AT TOP BOUNDARY OF LAYER, TDN IS AT BOTTOM BOUNDARY OF
4873 ! LAYER. TM IS LAYER PROGNOSTIC STATE TEMPERATURE.
4874 ! ----------------------------------------------------------------------
4875 IMPLICIT NONE
4876 INTEGER K
4877
4878 INTEGER NSOIL
4879 REAL DZ
4880 REAL DZH
4881 REAL T0
4882 REAL TAVG
4883 REAL TDN
4884 REAL TM
4885 REAL TUP
4886 REAL X0
4887 REAL XDN
4888 REAL XUP
4889
4890 REAL ZSOIL (NSOIL)
4891
4892 ! ----------------------------------------------------------------------
4893 PARAMETER (T0 = 2.7315E2)
4894 IF (K .eq. 1) THEN
4895 DZ = - ZSOIL (1)
4896 ELSE
4897 DZ = ZSOIL (K -1) - ZSOIL (K)
4898 END IF
4899
4900 DZH = DZ *0.5
4901 IF (TUP .lt. T0) THEN
4902 IF (TM .lt. T0) THEN
4903 ! ----------------------------------------------------------------------
4904 ! TUP, TM, TDN < T0
4905 ! ----------------------------------------------------------------------
4906 IF (TDN .lt. T0) THEN
4907 TAVG = (TUP + 2.0* TM + TDN)/ 4.0
4908 ! ----------------------------------------------------------------------
4909 ! TUP & TM < T0, TDN .ge. T0
4910 ! ----------------------------------------------------------------------
4911 ELSE
4912 X0 = (T0- TM) * DZH / (TDN - TM)
4913 TAVG = 0.5 * (TUP * DZH + TM * (DZH + X0) + T0* ( &
4914 & 2.* DZH - X0)) / DZ
4915 END IF
4916 ELSE
4917 ! ----------------------------------------------------------------------
4918 ! TUP < T0, TM .ge. T0, TDN < T0
4919 ! ----------------------------------------------------------------------
4920 IF (TDN .lt. T0) THEN
4921 XUP = (T0- TUP) * DZH / (TM - TUP)
4922 XDN = DZH - (T0- TM) * DZH / (TDN - TM)
4923 TAVG = 0.5 * (TUP * XUP + T0* (2.* DZ - XUP - XDN) &
4924 & + TDN * XDN) / DZ
4925 ! ----------------------------------------------------------------------
4926 ! TUP < T0, TM .ge. T0, TDN .ge. T0
4927 ! ----------------------------------------------------------------------
4928 ELSE
4929 XUP = (T0- TUP) * DZH / (TM - TUP)
4930 TAVG = 0.5 * (TUP * XUP + T0* (2.* DZ - XUP)) / DZ
4931 END IF
4932 END IF
4933 ELSE
4934 IF (TM .lt. T0) THEN
4935 ! ----------------------------------------------------------------------
4936 ! TUP .ge. T0, TM < T0, TDN < T0
4937 ! ----------------------------------------------------------------------
4938 IF (TDN .lt. T0) THEN
4939 XUP = DZH - (T0- TUP) * DZH / (TM - TUP)
4940 TAVG = 0.5 * (T0* (DZ - XUP) + TM * (DZH + XUP) &
4941 & + TDN * DZH) / DZ
4942 ! ----------------------------------------------------------------------
4943 ! TUP .ge. T0, TM < T0, TDN .ge. T0
4944 ! ----------------------------------------------------------------------
4945 ELSE
4946 XUP = DZH - (T0- TUP) * DZH / (TM - TUP)
4947 XDN = (T0- TM) * DZH / (TDN - TM)
4948 TAVG = 0.5 * (T0* (2.* DZ - XUP - XDN) + TM * &
4949 & (XUP + XDN)) / DZ
4950 END IF
4951 ELSE
4952 ! ----------------------------------------------------------------------
4953 ! TUP .ge. T0, TM .ge. T0, TDN < T0
4954 ! ----------------------------------------------------------------------
4955 IF (TDN .lt. T0) THEN
4956 XDN = DZH - (T0- TM) * DZH / (TDN - TM)
4957 TAVG = (T0* (DZ - XDN) +0.5* (T0+ TDN)* XDN) / DZ
4958 ! ----------------------------------------------------------------------
4959 ! TUP .ge. T0, TM .ge. T0, TDN .ge. T0
4960 ! ----------------------------------------------------------------------
4961 ELSE
4962 TAVG = (TUP + 2.0* TM + TDN) / 4.0
4963 END IF
4964 END IF
4965 END IF
4966 ! ----------------------------------------------------------------------
4967 END SUBROUTINE TMPAVG
4968 ! ----------------------------------------------------------------------
4969
4970 SUBROUTINE TRANSP (ET,NSOIL,ETP1,SMC,CMC,ZSOIL,SHDFAC,SMCWLT, &
4971 & CMCMAX,PC,CFACTR,SMCREF,SFCTMP,Q2,NROOT, &
4972 & RTDIS)
4973
4974 ! ----------------------------------------------------------------------
4975 ! SUBROUTINE TRANSP
4976 ! ----------------------------------------------------------------------
4977 ! CALCULATE TRANSPIRATION FOR THE VEG CLASS.
4978 ! ----------------------------------------------------------------------
4979 IMPLICIT NONE
4980 INTEGER I
4981 INTEGER K
4982 INTEGER NSOIL
4983
4984 INTEGER NROOT
4985 REAL CFACTR
4986 REAL CMC
4987 REAL CMCMAX
4988 REAL DENOM
4989 REAL ET (NSOIL)
4990 REAL ETP1
4991 REAL ETP1A
4992 !.....REAL PART(NSOIL)
4993 REAL GX (7)
4994 REAL PC
4995 REAL Q2
4996 REAL RTDIS (NSOIL)
4997 REAL RTX
4998 REAL SFCTMP
4999 REAL SGX
5000 REAL SHDFAC
5001 REAL SMC (NSOIL)
5002 REAL SMCREF
5003 REAL SMCWLT
5004
5005 ! ----------------------------------------------------------------------
5006 ! INITIALIZE PLANT TRANSP TO ZERO FOR ALL SOIL LAYERS.
5007 ! ----------------------------------------------------------------------
5008 REAL ZSOIL (NSOIL)
5009 DO K = 1,NSOIL
5010 ET (K) = 0.
5011 ! ----------------------------------------------------------------------
5012 ! CALCULATE AN 'ADJUSTED' POTENTIAL TRANSPIRATION
5013 ! IF STATEMENT BELOW TO AVOID TANGENT LINEAR PROBLEMS NEAR ZERO
5014 ! NOTE: GX AND OTHER TERMS BELOW REDISTRIBUTE TRANSPIRATION BY LAYER,
5015 ! ET(K), AS A FUNCTION OF SOIL MOISTURE AVAILABILITY, WHILE PRESERVING
5016 ! TOTAL ETP1A.
5017 ! ----------------------------------------------------------------------
5018 END DO
5019 IF (CMC .ne. 0.0) THEN
5020 ETP1A = SHDFAC * PC * ETP1 * (1.0- (CMC / CMCMAX) ** CFACTR)
5021 ELSE
5022 ETP1A = SHDFAC * PC * ETP1
5023 END IF
5024 SGX = 0.0
5025 DO I = 1,NROOT
5026 GX (I) = ( SMC (I) - SMCWLT ) / ( SMCREF - SMCWLT )
5027 GX (I) = MAX ( MIN ( GX (I), 1. ), 0. )
5028 SGX = SGX + GX (I)
5029 END DO
5030
5031 SGX = SGX / NROOT
5032 DENOM = 0.
5033 DO I = 1,NROOT
5034 RTX = RTDIS (I) + GX (I) - SGX
5035 GX (I) = GX (I) * MAX ( RTX, 0. )
5036 DENOM = DENOM + GX (I)
5037 END DO
5038
5039 IF (DENOM .le. 0.0) DENOM = 1.
5040 DO I = 1,NROOT
5041 ET (I) = ETP1A * GX (I) / DENOM
5042 ! ----------------------------------------------------------------------
5043 ! ABOVE CODE ASSUMES A VERTICALLY UNIFORM ROOT DISTRIBUTION
5044 ! CODE BELOW TESTS A VARIABLE ROOT DISTRIBUTION
5045 ! ----------------------------------------------------------------------
5046 ! ET(1) = ( ZSOIL(1) / ZSOIL(NROOT) ) * GX * ETP1A
5047 ! ET(1) = ( ZSOIL(1) / ZSOIL(NROOT) ) * ETP1A
5048 ! ----------------------------------------------------------------------
5049 ! USING ROOT DISTRIBUTION AS WEIGHTING FACTOR
5050 ! ----------------------------------------------------------------------
5051 ! ET(1) = RTDIS(1) * ETP1A
5052 ! ET(1) = ETP1A * PART(1)
5053 ! ----------------------------------------------------------------------
5054 ! LOOP DOWN THRU THE SOIL LAYERS REPEATING THE OPERATION ABOVE,
5055 ! BUT USING THE THICKNESS OF THE SOIL LAYER (RATHER THAN THE
5056 ! ABSOLUTE DEPTH OF EACH LAYER) IN THE FINAL CALCULATION.
5057 ! ----------------------------------------------------------------------
5058 ! DO K = 2,NROOT
5059 ! GX = ( SMC(K) - SMCWLT ) / ( SMCREF - SMCWLT )
5060 ! GX = MAX ( MIN ( GX, 1. ), 0. )
5061 ! TEST CANOPY RESISTANCE
5062 ! GX = 1.0
5063 ! ET(K) = ((ZSOIL(K)-ZSOIL(K-1))/ZSOIL(NROOT))*GX*ETP1A
5064 ! ET(K) = ((ZSOIL(K)-ZSOIL(K-1))/ZSOIL(NROOT))*ETP1A
5065 ! ----------------------------------------------------------------------
5066 ! USING ROOT DISTRIBUTION AS WEIGHTING FACTOR
5067 ! ----------------------------------------------------------------------
5068 ! ET(K) = RTDIS(K) * ETP1A
5069 ! ET(K) = ETP1A*PART(K)
5070 ! END DO
5071 END DO
5072 ! ----------------------------------------------------------------------
5073 END SUBROUTINE TRANSP
5074 ! ----------------------------------------------------------------------
5075
5076 SUBROUTINE WDFCND (WDF,WCND,SMC,SMCMAX,BEXP,DKSAT,DWSAT, &
5077 & SICEMAX)
5078
5079 ! ----------------------------------------------------------------------
5080 ! SUBROUTINE WDFCND
5081 ! ----------------------------------------------------------------------
5082 ! CALCULATE SOIL WATER DIFFUSIVITY AND SOIL HYDRAULIC CONDUCTIVITY.
5083 ! ----------------------------------------------------------------------
5084 IMPLICIT NONE
5085 REAL BEXP
5086 REAL DKSAT
5087 REAL DWSAT
5088 REAL EXPON
5089 REAL FACTR1
5090 REAL FACTR2
5091 REAL SICEMAX
5092 REAL SMC
5093 REAL SMCMAX
5094 REAL VKwgt
5095 REAL WCND
5096
5097 ! ----------------------------------------------------------------------
5098 ! CALC THE RATIO OF THE ACTUAL TO THE MAX PSBL SOIL H2O CONTENT
5099 ! ----------------------------------------------------------------------
5100 REAL WDF
5101 FACTR1 = 0.2 / SMCMAX
5102
5103 ! ----------------------------------------------------------------------
5104 ! PREP AN EXPNTL COEF AND CALC THE SOIL WATER DIFFUSIVITY
5105 ! ----------------------------------------------------------------------
5106 FACTR2 = SMC / SMCMAX
5107 EXPON = BEXP + 2.0
5108
5109 ! ----------------------------------------------------------------------
5110 ! FROZEN SOIL HYDRAULIC DIFFUSIVITY. VERY SENSITIVE TO THE VERTICAL
5111 ! GRADIENT OF UNFROZEN WATER. THE LATTER GRADIENT CAN BECOME VERY
5112 ! EXTREME IN FREEZING/THAWING SITUATIONS, AND GIVEN THE RELATIVELY
5113 ! FEW AND THICK SOIL LAYERS, THIS GRADIENT SUFFERES SERIOUS
5114 ! TRUNCTION ERRORS YIELDING ERRONEOUSLY HIGH VERTICAL TRANSPORTS OF
5115 ! UNFROZEN WATER IN BOTH DIRECTIONS FROM HUGE HYDRAULIC DIFFUSIVITY.
5116 ! THEREFORE, WE FOUND WE HAD TO ARBITRARILY CONSTRAIN WDF
5117 ! --
5118 ! VERSION D_10CM: ........ FACTR1 = 0.2/SMCMAX
5119 ! WEIGHTED APPROACH...................... PABLO GRUNMANN, 28_SEP_1999.
5120 ! ----------------------------------------------------------------------
5121 WDF = DWSAT * FACTR2 ** EXPON
5122 IF (SICEMAX .gt. 0.0) THEN
5123 VKWGT = 1./ (1. + (500.* SICEMAX)**3.)
5124 WDF = VKWGT * WDF + (1. - VKWGT)* DWSAT * FACTR1** EXPON
5125 ! ----------------------------------------------------------------------
5126 ! RESET THE EXPNTL COEF AND CALC THE HYDRAULIC CONDUCTIVITY
5127 ! ----------------------------------------------------------------------
5128 END IF
5129 EXPON = (2.0 * BEXP) + 3.0
5130 WCND = DKSAT * FACTR2 ** EXPON
5131
5132 ! ----------------------------------------------------------------------
5133 END SUBROUTINE WDFCND
5134 ! ----------------------------------------------------------------------
5135
5136 SUBROUTINE SFCDIF_off (ZLM,Z0,THZ0,THLM,SFCSPD,CZIL,AKMS,AKHS)
5137
5138 ! ----------------------------------------------------------------------
5139 ! SUBROUTINE SFCDIF (renamed SFCDIF_off to avoid clash with Eta PBL)
5140 ! ----------------------------------------------------------------------
5141 ! CALCULATE SURFACE LAYER EXCHANGE COEFFICIENTS VIA ITERATIVE PROCESS.
5142 ! SEE CHEN ET AL (1997, BLM)
5143 ! ----------------------------------------------------------------------
5144
5145 IMPLICIT NONE
5146 REAL WWST, WWST2, G, VKRM, EXCM, BETA, BTG, ELFC, WOLD, WNEW
5147 REAL PIHF, EPSU2, EPSUST, EPSIT, EPSA, ZTMIN, ZTMAX, HPBL, &
5148 & SQVISC
5149 REAL RIC, RRIC, FHNEU, RFC, RFAC, ZZ, PSLMU, PSLMS, PSLHU, &
5150 & PSLHS
5151 REAL XX, PSPMU, YY, PSPMS, PSPHU, PSPHS, ZLM, Z0, THZ0, THLM
5152 REAL SFCSPD, CZIL, AKMS, AKHS, ZILFC, ZU, ZT, RDZ, CXCH
5153 REAL DTHV, DU2, BTGH, WSTAR2, USTAR, ZSLU, ZSLT, RLOGU, RLOGT
5154 REAL RLMO, ZETALT, ZETALU, ZETAU, ZETAT, XLU4, XLT4, XU4, XT4
5155 !CC ......REAL ZTFC
5156
5157 REAL XLU, XLT, XU, XT, PSMZ, SIMM, PSHZ, SIMH, USTARK, RLMN, &
5158 & RLMA
5159
5160 INTEGER ITRMX, ILECH, ITR
5161 PARAMETER &
5162 & (WWST = 1.2,WWST2 = WWST * WWST,G = 9.8,VKRM = 0.40, &
5163 & EXCM = 0.001 &
5164 & ,BETA = 1./270.,BTG = BETA * G,ELFC = VKRM * BTG &
5165 & ,WOLD =.15,WNEW = 1. - WOLD,ITRMX = 05, &
5166 & PIHF = 3.14159265/2.)
5167 PARAMETER &
5168 & (EPSU2 = 1.E-4,EPSUST = 0.07,EPSIT = 1.E-4,EPSA = 1.E-8 &
5169 & ,ZTMIN = -5.,ZTMAX = 1.,HPBL = 1000.0 &
5170 & ,SQVISC = 258.2)
5171 PARAMETER &
5172 & (RIC = 0.183,RRIC = 1.0/ RIC,FHNEU = 0.8,RFC = 0.191 &
5173 & ,RFAC = RIC / (FHNEU * RFC * RFC))
5174
5175 ! ----------------------------------------------------------------------
5176 ! NOTE: THE TWO CODE BLOCKS BELOW DEFINE FUNCTIONS
5177 ! ----------------------------------------------------------------------
5178 ! LECH'S SURFACE FUNCTIONS
5179 ! ----------------------------------------------------------------------
5180 PSLMU (ZZ)= -0.96* log (1.0-4.5* ZZ)
5181 PSLMS (ZZ)= ZZ * RRIC -2.076* (1. -1./ (ZZ +1.))
5182 PSLHU (ZZ)= -0.96* log (1.0-4.5* ZZ)
5183
5184 ! ----------------------------------------------------------------------
5185 ! PAULSON'S SURFACE FUNCTIONS
5186 ! ----------------------------------------------------------------------
5187 PSLHS (ZZ)= ZZ * RFAC -2.076* (1. -1./ (ZZ +1.))
5188 PSPMU (XX)= -2.* log ( (XX +1.)*0.5) - log ( (XX * XX +1.)*0.5) &
5189 & +2.* ATAN (XX) &
5190 &- PIHF
5191 PSPMS (YY)= 5.* YY
5192 PSPHU (XX)= -2.* log ( (XX * XX +1.)*0.5)
5193
5194 ! ----------------------------------------------------------------------
5195 ! THIS ROUTINE SFCDIF CAN HANDLE BOTH OVER OPEN WATER (SEA, OCEAN) AND
5196 ! OVER SOLID SURFACE (LAND, SEA-ICE).
5197 ! ----------------------------------------------------------------------
5198 PSPHS (YY)= 5.* YY
5199
5200 ! ----------------------------------------------------------------------
5201 ! ZTFC: RATIO OF ZOH/ZOM LESS OR EQUAL THAN 1
5202 ! C......ZTFC=0.1
5203 ! CZIL: CONSTANT C IN Zilitinkevich, S. S.1995,:NOTE ABOUT ZT
5204 ! ----------------------------------------------------------------------
5205 ILECH = 0
5206
5207 ! ----------------------------------------------------------------------
5208 ZILFC = - CZIL * VKRM * SQVISC
5209 ! C.......ZT=Z0*ZTFC
5210 ZU = Z0
5211 RDZ = 1./ ZLM
5212 CXCH = EXCM * RDZ
5213 DTHV = THLM - THZ0
5214
5215 ! ----------------------------------------------------------------------
5216 ! BELJARS CORRECTION OF USTAR
5217 ! ----------------------------------------------------------------------
5218 DU2 = MAX (SFCSPD * SFCSPD,EPSU2)
5219 !cc If statements to avoid TANGENT LINEAR problems near zero
5220 BTGH = BTG * HPBL
5221 IF (BTGH * AKHS * DTHV .ne. 0.0) THEN
5222 WSTAR2 = WWST2* ABS (BTGH * AKHS * DTHV)** (2./3.)
5223 ELSE
5224 WSTAR2 = 0.0
5225 END IF
5226
5227 ! ----------------------------------------------------------------------
5228 ! ZILITINKEVITCH APPROACH FOR ZT
5229 ! ----------------------------------------------------------------------
5230 USTAR = MAX (SQRT (AKMS * SQRT (DU2+ WSTAR2)),EPSUST)
5231
5232 ! ----------------------------------------------------------------------
5233 ZT = EXP (ZILFC * SQRT (USTAR * Z0))* Z0
5234 ZSLU = ZLM + ZU
5235 ! PRINT*,'ZSLT=',ZSLT
5236 ! PRINT*,'ZLM=',ZLM
5237 ! PRINT*,'ZT=',ZT
5238
5239 ZSLT = ZLM + ZT
5240 RLOGU = log (ZSLU / ZU)
5241
5242 RLOGT = log (ZSLT / ZT)
5243 ! PRINT*,'RLMO=',RLMO
5244 ! PRINT*,'ELFC=',ELFC
5245 ! PRINT*,'AKHS=',AKHS
5246 ! PRINT*,'DTHV=',DTHV
5247 ! PRINT*,'USTAR=',USTAR
5248
5249 RLMO = ELFC * AKHS * DTHV / USTAR **3
5250 ! ----------------------------------------------------------------------
5251 ! 1./MONIN-OBUKKHOV LENGTH-SCALE
5252 ! ----------------------------------------------------------------------
5253 DO ITR = 1,ITRMX
5254 ZETALT = MAX (ZSLT * RLMO,ZTMIN)
5255 RLMO = ZETALT / ZSLT
5256 ZETALU = ZSLU * RLMO
5257 ZETAU = ZU * RLMO
5258
5259 ZETAT = ZT * RLMO
5260 IF (ILECH .eq. 0) THEN
5261 IF (RLMO .lt. 0.)THEN
5262 XLU4 = 1. -16.* ZETALU
5263 XLT4 = 1. -16.* ZETALT
5264 XU4 = 1. -16.* ZETAU
5265
5266 XT4 = 1. -16.* ZETAT
5267 XLU = SQRT (SQRT (XLU4))
5268 XLT = SQRT (SQRT (XLT4))
5269 XU = SQRT (SQRT (XU4))
5270
5271 XT = SQRT (SQRT (XT4))
5272 ! PRINT*,'-----------1------------'
5273 ! PRINT*,'PSMZ=',PSMZ
5274 ! PRINT*,'PSPMU(ZETAU)=',PSPMU(ZETAU)
5275 ! PRINT*,'XU=',XU
5276 ! PRINT*,'------------------------'
5277 PSMZ = PSPMU (XU)
5278 SIMM = PSPMU (XLU) - PSMZ + RLOGU
5279 PSHZ = PSPHU (XT)
5280 SIMH = PSPHU (XLT) - PSHZ + RLOGT
5281 ELSE
5282 ZETALU = MIN (ZETALU,ZTMAX)
5283 ZETALT = MIN (ZETALT,ZTMAX)
5284 ! PRINT*,'-----------2------------'
5285 ! PRINT*,'PSMZ=',PSMZ
5286 ! PRINT*,'PSPMS(ZETAU)=',PSPMS(ZETAU)
5287 ! PRINT*,'ZETAU=',ZETAU
5288 ! PRINT*,'------------------------'
5289 PSMZ = PSPMS (ZETAU)
5290 SIMM = PSPMS (ZETALU) - PSMZ + RLOGU
5291 PSHZ = PSPHS (ZETAT)
5292 SIMH = PSPHS (ZETALT) - PSHZ + RLOGT
5293 END IF
5294 ! ----------------------------------------------------------------------
5295 ! LECH'S FUNCTIONS
5296 ! ----------------------------------------------------------------------
5297 ELSE
5298 IF (RLMO .lt. 0.)THEN
5299 ! PRINT*,'-----------3------------'
5300 ! PRINT*,'PSMZ=',PSMZ
5301 ! PRINT*,'PSLMU(ZETAU)=',PSLMU(ZETAU)
5302 ! PRINT*,'ZETAU=',ZETAU
5303 ! PRINT*,'------------------------'
5304 PSMZ = PSLMU (ZETAU)
5305 SIMM = PSLMU (ZETALU) - PSMZ + RLOGU
5306 PSHZ = PSLHU (ZETAT)
5307 SIMH = PSLHU (ZETALT) - PSHZ + RLOGT
5308 ELSE
5309 ZETALU = MIN (ZETALU,ZTMAX)
5310
5311 ZETALT = MIN (ZETALT,ZTMAX)
5312 ! PRINT*,'-----------4------------'
5313 ! PRINT*,'PSMZ=',PSMZ
5314 ! PRINT*,'PSLMS(ZETAU)=',PSLMS(ZETAU)
5315 ! PRINT*,'ZETAU=',ZETAU
5316 ! PRINT*,'------------------------'
5317 PSMZ = PSLMS (ZETAU)
5318 SIMM = PSLMS (ZETALU) - PSMZ + RLOGU
5319 PSHZ = PSLHS (ZETAT)
5320 SIMH = PSLHS (ZETALT) - PSHZ + RLOGT
5321 END IF
5322 ! ----------------------------------------------------------------------
5323 ! BELJAARS CORRECTION FOR USTAR
5324 ! ----------------------------------------------------------------------
5325 END IF
5326
5327 ! ----------------------------------------------------------------------
5328 ! ZILITINKEVITCH FIX FOR ZT
5329 ! ----------------------------------------------------------------------
5330 USTAR = MAX (SQRT (AKMS * SQRT (DU2+ WSTAR2)),EPSUST)
5331
5332 ZT = EXP (ZILFC * SQRT (USTAR * Z0))* Z0
5333 ZSLT = ZLM + ZT
5334 !-----------------------------------------------------------------------
5335 RLOGT = log (ZSLT / ZT)
5336 USTARK = USTAR * VKRM
5337 AKMS = MAX (USTARK / SIMM,CXCH)
5338 !-----------------------------------------------------------------------
5339 ! IF STATEMENTS TO AVOID TANGENT LINEAR PROBLEMS NEAR ZERO
5340 !-----------------------------------------------------------------------
5341 AKHS = MAX (USTARK / SIMH,CXCH)
5342 IF (BTGH * AKHS * DTHV .ne. 0.0) THEN
5343 WSTAR2 = WWST2* ABS (BTGH * AKHS * DTHV)** (2./3.)
5344 ELSE
5345 WSTAR2 = 0.0
5346 END IF
5347 !-----------------------------------------------------------------------
5348 RLMN = ELFC * AKHS * DTHV / USTAR **3
5349 !-----------------------------------------------------------------------
5350 ! IF(ABS((RLMN-RLMO)/RLMA).LT.EPSIT) GO TO 110
5351 !-----------------------------------------------------------------------
5352 RLMA = RLMO * WOLD+ RLMN * WNEW
5353 !-----------------------------------------------------------------------
5354 RLMO = RLMA
5355 ! PRINT*,'----------------------------'
5356 ! PRINT*,'SFCDIF OUTPUT ! ! ! ! ! ! ! ! ! ! ! !'
5357
5358 ! PRINT*,'ZLM=',ZLM
5359 ! PRINT*,'Z0=',Z0
5360 ! PRINT*,'THZ0=',THZ0
5361 ! PRINT*,'THLM=',THLM
5362 ! PRINT*,'SFCSPD=',SFCSPD
5363 ! PRINT*,'CZIL=',CZIL
5364 ! PRINT*,'AKMS=',AKMS
5365 ! PRINT*,'AKHS=',AKHS
5366 ! PRINT*,'----------------------------'
5367
5368 END DO
5369 ! ----------------------------------------------------------------------
5370 END SUBROUTINE SFCDIF_off
5371 ! ----------------------------------------------------------------------
5372
5373 END MODULE module_sf_noahlsm