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,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, & !H
54 SH2O,SNOWH, & !H
55 U_PHY,V_PHY, & !I
56 SNOALB,SHDMIN,SHDMAX, & !I
57 ACSNOM,ACSNOW, & !O
58 ids,ide, jds,jde, kds,kde, &
59 ims,ime, jms,jme, kms,kme, &
60 its,ite, jts,jte, kts,kte, &
61 ucmcall, &
62 !Optional Urban
63 TR_URB2D,TB_URB2D,TG_URB2D,TC_URB2D,QC_URB2D, & !H urban
64 UC_URB2D, & !H urban
65 XXXR_URB2D,XXXB_URB2D,XXXG_URB2D,XXXC_URB2D, & !H urban
66 TRL_URB3D,TBL_URB3D,TGL_URB3D, & !H urban
67 SH_URB2D,LH_URB2D,G_URB2D,RN_URB2D,TS_URB2D, & !H urban
68 PSIM_URB2D,PSIH_URB2D,U10_URB2D,V10_URB2D, & !O urban
69 GZ1OZ0_URB2D, AKMS_URB2D, & !O urban
70 TH2_URB2D,Q2_URB2D, UST_URB2D, & !O urban
71 DECLIN_URB,COSZ_URB2D,OMG_URB2D, & !I urban
72 XLAT_URB2D, & !I urban
73 num_roof_layers, num_wall_layers, & !I urban
74 num_road_layers, DZR, DZB, DZG, & !I urban
75 FRC_URB2D,UTYPE_URB2D) !O
76 !----------------------------------------------------------------
77 IMPLICIT NONE
78 !----------------------------------------------------------------
79 !----------------------------------------------------------------
80 ! --- atmospheric (WRF generic) variables
81 !-- DT time step (seconds)
82 !-- DZ8W thickness of layers (m)
83 !-- T3D temperature (K)
84 !-- QV3D 3D water vapor mixing ratio (Kg/Kg)
85 !-- P3D 3D pressure (Pa)
86 !-- FLHC exchange coefficient for heat (m/s)
87 !-- FLQC exchange coefficient for moisture (m/s)
88 !-- PSFC surface pressure (Pa)
89 !-- XLAND land mask (1 for land, 2 for water)
90 !-- QGH saturated mixing ratio at 2 meter
91 !-- GSW downward short wave flux at ground surface (W/m^2)
92 !-- GLW downward long wave flux at ground surface (W/m^2)
93 !-- History variables
94 !-- CANWAT canopy moisture content (mm)
95 !-- TSK surface temperature (K)
96 !-- TSLB soil temp (k)
97 !-- SMOIS total soil moisture content (volumetric fraction)
98 !-- SH2O unfrozen soil moisture content (volumetric fraction)
99 ! note: frozen soil moisture (i.e., soil ice) = SMOIS - SH2O
100 !-- SNOWH actual snow depth (m)
101 !-- SNOW liquid water-equivalent snow depth (m)
102 !-- ALBEDO time-varying surface albedo including snow effect (unitless fraction)
103 !-- ALBBCK background surface albedo (unitless fraction)
104 !-- CHS surface exchange coefficient for heat and moisture (m s-1);
105 !-- CHS2 2m surface exchange coefficient for heat (m s-1);
106 !-- CQS2 2m surface exchange coefficient for moisture (m s-1);
107 ! --- soil variables
108 !-- num_soil_layers the number of soil layers
109 !-- ZS depths of centers of soil layers (m)
110 !-- DZS thicknesses of soil layers (m)
111 !-- SLDPTH thickness of each soil layer (m, same as DZS)
112 !-- TMN soil temperature at lower boundary (K)
113 !-- SMCWLT wilting point (volumetric)
114 !-- SMCDRY dry soil moisture threshold where direct evap from
115 ! top soil layer ends (volumetric)
116 !-- SMCREF soil moisture threshold below which transpiration begins to
117 ! stress (volumetric)
118 !-- SMCMAX porosity, i.e. saturated value of soil moisture (volumetric)
119 !-- NROOT number of root layers, a function of veg type, determined
120 ! in subroutine redprm.
121 !-- SMSTAV Soil moisture availability for evapotranspiration (
122 ! fraction between SMCWLT and SMCMXA)
123 !-- SMSTOT Total soil moisture content frozen+unfrozen) in the soil column (mm)
124 ! --- snow variables
125 !-- SNOWC fraction snow coverage (0-1.0)
126 ! --- vegetation variables
127 !-- SNOALB upper bound on maximum albedo over deep snow
128 !-- SHDMIN minimum areal fractional coverage of annual green vegetation
129 !-- SHDMAX maximum areal fractional coverage of annual green vegetation
130 !-- XLAI leaf area index (dimensionless)
131 !-- Z0BRD Background fixed roughness length (M)
132 !-- Z0 Background vroughness length (M) as function
133 !-- ZNT Time varying roughness length (M) as function
134 !-- ALBD(IVGTPK,ISN) background albedo reading from a table
135 ! --- LSM output
136 !-- HFX upward heat flux at the surface (W/m^2)
137 !-- QFX upward moisture flux at the surface (kg/m^2/s)
138 !-- LH upward moisture flux at the surface (W m-2)
139 !-- GRDFLX(I,J) ground heat flux (W m-2)
140 !-- FDOWN radiation forcing at the surface (W m-2) = SOLDN*(1-alb)+LWDN
141 !----------------------------------------------------------------------------
142 !-- EC canopy water evaporation ((W m-2)
143 !-- EDIR direct soil evaporation (W m-2)
144 !-- ET plant transpiration from a particular root layer (W m-2)
145 !-- ETT total plant transpiration (W m-2)
146 !-- ESNOW sublimation from (or deposition to if <0) snowpack (W m-2)
147 !-- DRIP through-fall of precip and/or dew in excess of canopy
148 ! water-holding capacity (m)
149 !-- DEW dewfall (or frostfall for t<273.15) (M)
150 ! ----------------------------------------------------------------------
151 !-- BETA ratio of actual/potential evap (dimensionless)
152 !-- ETP potential evaporation (W m-2)
153 ! ----------------------------------------------------------------------
154 !-- FLX1 precip-snow sfc (W m-2)
155 !-- FLX2 freezing rain latent heat flux (W m-2)
156 !-- FLX3 phase-change heat flux from snowmelt (W m-2)
157 ! ----------------------------------------------------------------------
158 !-- ACSNOM snow melt (mm) (water equivalent)
159 !-- ACSNOW accumulated snow fall (mm) (water equivalent)
160 ! ----------------------------------------------------------------------
161 !-- RUNOFF1 surface runoff (m s-1), not infiltrating the surface
162 !-- RUNOFF2 subsurface runoff (m s-1), drainage out bottom of last
163 ! soil layer (baseflow)
164 ! important note: here RUNOFF2 is actually the sum of RUNOFF2 and RUNOFF3
165 !-- RUNOFF3 numerical trunctation in excess of porosity (smcmax)
166 ! for a given soil layer at the end of a time step (m s-1).
167 ! ----------------------------------------------------------------------
168 !-- RC canopy resistance (s m-1)
169 !-- PC plant coefficient (unitless fraction, 0-1) where PC*ETP = actual transp
170 !-- RSMIN minimum canopy resistance (s m-1)
171 !-- RCS incoming solar rc factor (dimensionless)
172 !-- RCT air temperature rc factor (dimensionless)
173 !-- RCQ atmos vapor pressure deficit rc factor (dimensionless)
174 !-- RCSOIL soil moisture rc factor (dimensionless)
175
176 !-- EMISS surface emissivity (between 0 and 1)
177
178 !-- ROVCP R/CP
179 ! (R_d/R_v) (dimensionless)
180 !-- ids start index for i in domain
181 !-- ide end index for i in domain
182 !-- jds start index for j in domain
183 !-- jde end index for j in domain
184 !-- kds start index for k in domain
185 !-- kde end index for k in domain
186 !-- ims start index for i in memory
187 !-- ime end index for i in memory
188 !-- jms start index for j in memory
189 !-- jme end index for j in memory
190 !-- kms start index for k in memory
191 !-- kme end index for k in memory
192 !-- its start index for i in tile
193 !-- ite end index for i in tile
194 !-- jts start index for j in tile
195 !-- jte end index for j in tile
196 !-- kts start index for k in tile
197 !-- kte end index for k in tile
198 !----------------------------------------------------------------
199
200 ! IN only
201
202 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
203 ims,ime, jms,jme, kms,kme, &
204 its,ite, jts,jte, kts,kte
205
206 INTEGER, INTENT(IN ) :: ucmcall !urban
207
208 REAL, DIMENSION( ims:ime, jms:jme ) , &
209 INTENT(IN ) :: TMN, &
210 XLAND, &
211 XICE, &
212 VEGFRA, &
213 SHDMIN, &
214 SHDMAX, &
215 SNOALB, &
216 GSW, &
217 GLW, &
218 Z0, &
219 ALBBCK, &
220 RAINBL, &
221 EMISS
222
223 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
224 INTENT(IN ) :: QV3D, &
225 p8w3D, &
226 DZ8W, &
227 T3D
228 REAL, DIMENSION( ims:ime, jms:jme ) , &
229 INTENT(IN ) :: QGH, &
230 CHS, &
231 CPM
232
233 INTEGER, DIMENSION( ims:ime, jms:jme ) , &
234 INTENT(IN ) :: IVGTYP, &
235 ISLTYP
236
237 INTEGER, INTENT(IN) :: num_soil_layers,ITIMESTEP
238
239 REAL, INTENT(IN ) :: DT,ROVCP
240
241 REAL, DIMENSION(1:num_soil_layers), INTENT(IN)::DZS
242
243 ! IN and OUT
244
245 REAL, DIMENSION( ims:ime , 1:num_soil_layers, jms:jme ), &
246 INTENT(INOUT) :: SMOIS, & ! total soil moisture
247 SH2O, & ! new soil liquid
248 TSLB ! TSLB STEMP
249
250 REAL, DIMENSION( ims:ime, jms:jme ) , &
251 INTENT(INOUT) :: TSK, & !was TGB (temperature)
252 HFX, &
253 QFX, &
254 LH, &
255 GRDFLX, &
256 QSFC,&
257 CQS2,&
258 CHS2,&
259 SNOW, &
260 SNOWC, &
261 SNOWH, & !new
262 CANWAT, &
263 SMSTAV, &
264 SMSTOT, &
265 SFCRUNOFF, &
266 UDRUNOFF, &
267 ACSNOM, &
268 ACSNOW, &
269 ALBEDO, &
270 ZNT
271
272
273 ! Local variables (moved here from driver to make routine thread safe, 20031007 jm)
274
275 REAL, DIMENSION(1:num_soil_layers) :: ET
276
277 REAL :: BETA, ETP, SSOIL,EC, EDIR, ESNOW, ETT, &
278 FLX1,FLX2,FLX3, DRIP,DEW,FDOWN,RC,PC,RSMIN,XLAI, &
279 RCS,RCT,RCQ,RCSOIL
280
281
282 ! DECLARATIONS - LOGICAL
283 ! ----------------------------------------------------------------------
284 LOGICAL, PARAMETER :: LOCAL=.false.
285 LOGICAL :: FRZGRA, SNOWNG
286
287 LOGICAL :: IPRINT
288
289 ! ----------------------------------------------------------------------
290 ! DECLARATIONS - INTEGER
291 ! ----------------------------------------------------------------------
292 INTEGER :: I,J, ICE,NSOIL,SLOPETYP,SOILTYP,VEGTYP
293 INTEGER :: NROOT
294 INTEGER :: KZ ,K
295 INTEGER :: NS
296 ! ----------------------------------------------------------------------
297 ! DECLARATIONS - REAL
298 ! ----------------------------------------------------------------------
299
300 REAL :: SHMIN,SHMAX,DQSDT2,LWDN,PRCP,PRCPRAIN, &
301 Q2SAT,SFCPRS,SFCSPD,SFCTMP,SHDFAC,SNOALB1, &
302 SOLDN,TBOT,ZLVL, Q2K,ALBBRD, ALBEDOK, ETA, ETA_KINEMATIC, &
303 Z0K,RUNOFF1,RUNOFF2,RUNOFF3,SHEAT
304
305 REAL :: EMISSI
306
307 REAL :: SNCOVR,SNEQV,SNOWHK,CMC, CHK,TH2
308
309 REAL :: SMCDRY,SMCMAX,SMCREF,SMCWLT,SNOMLT,SOILM,SOILW,Q1,T1
310
311 REAL :: DUMMY,Z0BRD
312 !
313 REAL :: COSZ, SOLARDIRECT
314 !
315 REAL, DIMENSION(1:num_soil_layers):: SLDPTH, STC,SMC,SWC
316 !
317 REAL, DIMENSION(1:num_soil_layers) :: ZSOIL, RTDIS
318 REAL, PARAMETER :: TRESH=.95E0, A2=17.67,A3=273.15,A4=29.65, &
319 T0=273.16E0, ELWV=2.50E6, A23M4=A2*(A3-A4)
320
321 ! ----------------------------------------------------------------------
322 ! DECLARATIONS START - urban
323 ! ----------------------------------------------------------------------
324
325 ! input variables surface_driver --> lsm
326 INTEGER, INTENT(IN) :: num_roof_layers
327 INTEGER, INTENT(IN) :: num_wall_layers
328 INTEGER, INTENT(IN) :: num_road_layers
329 REAL, OPTIONAL, DIMENSION(1:num_roof_layers), INTENT(IN) :: DZR
330 REAL, OPTIONAL, DIMENSION(1:num_wall_layers), INTENT(IN) :: DZB
331 REAL, OPTIONAL, DIMENSION(1:num_road_layers), INTENT(IN) :: DZG
332 REAL, OPTIONAL, INTENT(IN) :: DECLIN_URB
333 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: COSZ_URB2D
334 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: OMG_URB2D
335 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: XLAT_URB2D
336 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: U_PHY
337 REAL, OPTIONAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: V_PHY
338
339 ! input variables lsm --> urban
340 INTEGER :: UTYPE_URB ! urban type [urban=1, suburban=2, rural=3]
341 REAL :: TA_URB ! potential temp at 1st atmospheric level [K]
342 REAL :: QA_URB ! mixing ratio at 1st atmospheric level [kg/kg]
343 REAL :: UA_URB ! wind speed at 1st atmospheric level [m/s]
344 REAL :: U1_URB ! u at 1st atmospheric level [m/s]
345 REAL :: V1_URB ! v at 1st atmospheric level [m/s]
346 REAL :: SSG_URB ! downward total short wave radiation [W/m/m]
347 REAL :: LLG_URB ! downward long wave radiation [W/m/m]
348 REAL :: RAIN_URB ! precipitation [mm/h]
349 REAL :: RHOO_URB ! air density [kg/m^3]
350 REAL :: ZA_URB ! first atmospheric level [m]
351 REAL :: DELT_URB ! time step [s]
352 REAL :: SSGD_URB ! downward direct short wave radiation [W/m/m]
353 REAL :: SSGQ_URB ! downward diffuse short wave radiation [W/m/m]
354 REAL :: XLAT_URB ! latitude [deg]
355 REAL :: COSZ_URB ! cosz
356 REAL :: OMG_URB ! hour angle
357 REAL :: ZNT_URB ! roughness length [m]
358 REAL :: TR_URB
359 REAL :: TB_URB
360 REAL :: TG_URB
361 REAL :: TC_URB
362 REAL :: QC_URB
363 REAL :: UC_URB
364 REAL :: XXXR_URB
365 REAL :: XXXB_URB
366 REAL :: XXXG_URB
367 REAL :: XXXC_URB
368 REAL, DIMENSION(1:num_roof_layers) :: TRL_URB ! roof layer temp [K]
369 REAL, DIMENSION(1:num_wall_layers) :: TBL_URB ! wall layer temp [K]
370 REAL, DIMENSION(1:num_road_layers) :: TGL_URB ! road layer temp [K]
371 LOGICAL :: LSOLAR_URB
372 ! state variable surface_driver <--> lsm <--> urban
373 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TR_URB2D
374 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TB_URB2D
375 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TG_URB2D
376 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TC_URB2D
377 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: QC_URB2D
378 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: UC_URB2D
379 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXR_URB2D
380 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXB_URB2D
381 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXG_URB2D
382 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: XXXC_URB2D
383 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: SH_URB2D
384 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: LH_URB2D
385 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: G_URB2D
386 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: RN_URB2D
387 !
388 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: TS_URB2D
389
390 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_roof_layers, jms:jme ), INTENT(INOUT) :: TRL_URB3D
391 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_wall_layers, jms:jme ), INTENT(INOUT) :: TBL_URB3D
392 REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_road_layers, jms:jme ), INTENT(INOUT) :: TGL_URB3D
393
394 ! output variable lsm --> surface_driver
395 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: PSIM_URB2D
396 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: PSIH_URB2D
397 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: GZ1OZ0_URB2D
398 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: U10_URB2D
399 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: V10_URB2D
400 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: TH2_URB2D
401 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: Q2_URB2D
402 !
403 REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: AKMS_URB2D
404 !
405 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: UST_URB2D
406 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: FRC_URB2D
407 INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: UTYPE_URB2D
408
409
410 ! output variables urban --> lsm
411 REAL :: TS_URB ! surface radiative temperature [K]
412 REAL :: QS_URB ! surface humidity [-]
413 REAL :: SH_URB ! sensible heat flux [W/m/m]
414 REAL :: LH_URB ! latent heat flux [W/m/m]
415 REAL :: LH_KINEMATIC_URB ! latent heat flux, kinetic [kg/m/m/s]
416 REAL :: SW_URB ! upward short wave radiation flux [W/m/m]
417 REAL :: ALB_URB ! time-varying albedo [fraction]
418 REAL :: LW_URB ! upward long wave radiation flux [W/m/m]
419 REAL :: G_URB ! heat flux into the ground [W/m/m]
420 REAL :: RN_URB ! net radiation [W/m/m]
421 REAL :: PSIM_URB ! shear f for momentum [-]
422 REAL :: PSIH_URB ! shear f for heat [-]
423 REAL :: GZ1OZ0_URB ! shear f for heat [-]
424 REAL :: U10_URB ! wind u component at 10 m [m/s]
425 REAL :: V10_URB ! wind v component at 10 m [m/s]
426 REAL :: TH2_URB ! potential temperature at 2 m [K]
427 REAL :: Q2_URB ! humidity at 2 m [-]
428 REAL :: CHS_URB
429 REAL :: CHS2_URB
430 REAL :: UST_URB
431
432 ! ----------------------------------------------------------------------
433 ! DECLARATIONS END - urban
434 ! ----------------------------------------------------------------------
435
436 ! debug printout
437 IPRINT=.false.
438
439 SLOPETYP=2
440 ! SHDMIN=0.00
441
442
443 NSOIL=num_soil_layers
444
445 DO NS=1,NSOIL
446 SLDPTH(NS)=DZS(NS)
447 ENDDO
448
449 DO J=jts,jte
450
451 IF(ITIMESTEP.EQ.1)THEN
452 DO 50 I=its,ite
453 !*** initialize soil conditions for IHOP 31 May case
454 ! IF((XLAND(I,J)-1.5) < 0.)THEN
455 ! if (I==108.and.j==85) then
456 ! DO NS=1,NSOIL
457 ! SMOIS(I,NS,J)=0.10
458 ! SH2O(I,NS,J)=0.10
459 ! enddo
460 ! endif
461 ! ENDIF
462
463 !*** SET ZERO-VALUE FOR SOME OUTPUT DIAGNOSTIC ARRAYS
464 IF((XLAND(I,J)-1.5).GE.0.)THEN
465 ! check sea-ice point
466 IF(XICE(I,J).EQ.1..and.IPRINT)PRINT*,' sea-ice at water point, I=',I, &
467 'J=',J
468 !*** Open Water Case
469 SMSTAV(I,J)=1.0
470 SMSTOT(I,J)=1.0
471 DO NS=1,NSOIL
472 SMOIS(I,NS,J)=1.0
473 TSLB(I,NS,J)=273.16 !STEMP
474 ENDDO
475 ELSE
476 IF(XICE(I,J).EQ.1.)THEN
477 !*** SEA-ICE CASE
478 SMSTAV(I,J)=1.0
479 SMSTOT(I,J)=1.0
480 DO NS=1,NSOIL
481 SMOIS(I,NS,J)=1.0
482 ENDDO
483 ENDIF
484 ENDIF
485 !
486 50 CONTINUE
487 ENDIF ! end of initialization over ocean
488
489 !-----------------------------------------------------------------------
490 DO 100 I=its,ite
491 SFCPRS=P8w3D(i,1,j)
492 Q2K=QV3D(i,1,j)
493 Q2SAT=QGH(I,j)
494 SFCTMP=T3D(i,1,j)
495 ZLVL=0.5*DZ8W(i,1,j)
496 TH2=SFCTMP+(0.0097545*ZLVL)
497 EMISSI = EMISS(I,J)
498 LWDN=GLW(I,J)*EMISSI
499 SOLDN=GSW(I,J)/(1.0-albedo(i,j))
500 PRCP=RAINBL(i,j)/DT
501 VEGTYP=IVGTYP(I,J)
502 SOILTYP=ISLTYP(I,J)
503 SHDFAC=VEGFRA(I,J)/100.
504 T1=TSK(I,J)
505 CHK=CHS(I,J)
506 SHMIN=SHDMIN(I,J)/100. !NEW
507 SHMAX=SHDMAX(I,J)/100. !NEW
508 SNOALB1=SNOALB(I,J) !NEW
509 ! convert snow depth from mm to meter
510 SNEQV=SNOW(I,J)*0.001
511 ! SNOWHK=SNOWH(I,J)*0.001 ! check the unit of snowh which is assumed to be mm for now
512 SNOWHK=SNOWH(I,J) ! check the unit of snowh which is assumed to be m for now
513
514
515 !***
516 IF((XLAND(I,J)-1.5).GE.0.)THEN ! begining of land/sea if block
517 ! Open water points
518 !CC Q2SAT=PQ0/SFCPRS*EXP(A2*(TGDSA(I)-A3)/(TGDSA(I)-A4))
519 ! HFX(I,J)=CHFF*(TGDSA(I)-TH(I))
520 ! QFX(I,J)=RHO(I)*CHS(I)*(Q2SAT-QV(I))
521 ! SFCEVP(I,J)=SFCEVP(I,J)+QFX(I,J)*DT
522
523 ELSE
524 ! Land or sea-ice case
525
526 IF (XICE(I,J) .GT. 0.5) THEN
527 ICE=1
528 ELSE
529 ICE=0
530 ENDIF
531 DQSDT2=Q2SAT*A23M4/(SFCTMP-A4)**2
532
533 IF(ICE.EQ.0)THEN
534 TBOT=TMN(I,J)
535 ELSE
536 TBOT=271.16
537 ENDIF
538 IF(VEGTYP.EQ.25) SHDFAC=0.0000
539 IF(VEGTYP.EQ.26) SHDFAC=0.0000
540 IF(VEGTYP.EQ.27) SHDFAC=0.0000
541 IF(SOILTYP.EQ.14.AND.XICE(I,J).EQ.0.)THEN
542 IF(IPRINT)PRINT*,' SOIL TYPE FOUND TO BE WATER AT A LAND-POINT'
543 IF(IPRINT)PRINT*,i,j,'RESET SOIL in surfce.F'
544 SOILTYP=7
545 ENDIF
546 CMC=CANWAT(I,J)
547
548 !-------------------------------------------
549 !*** convert snow depth from mm to meter
550 !
551 ! IF(RDMAXALB) THEN
552 ! SNOALB=ALBMAX(I,J)*0.01
553 ! ELSE
554 ! SNOALB=MAXALB(IVGTPK)*0.01
555 ! ENDIF
556 ! IF(RDBRDALB) THEN
557 ! ALBBRD=ALBEDO(I,J)*0.01
558 ! ELSE
559 ! ALBBRD=ALBD(IVGTPK,ISN)*0.01
560 ! ENDIF
561
562 ! SNOALB1=0.80
563 ! SHMIN=0.00
564 ALBBRD=ALBBCK(I,J)
565 Z0BRD=Z0(I,J)
566 !FEI: temporaray arrays above need to be changed later by using SI
567
568 DO 70 NS=1,NSOIL
569 SMC(NS)=SMOIS(I,NS,J)
570 STC(NS)=TSLB(I,NS,J) !STEMP
571 SWC(NS)=SH2O(I,NS,J)
572 70 CONTINUE
573 !
574 if ( (SNEQV.ne.0..AND.SNOWHK.eq.0.).or.(SNOWHK.le.SNEQV) )THEN
575 SNOWHK= 5.*SNEQV
576 endif
577 !
578
579 !Fei: urban. for urban surface, if calling UCM, redefine urban as 5: Cropland/Grassland Mosaic
580
581 IF(UCMCALL == 1 ) THEN
582 IF( IVGTYP(I,J) == 1 .or. IVGTYP(I,J) == 31 .or. &
583 IVGTYP(I,J) == 32 .or. IVGTYP(I,J) == 33) THEN
584 VEGTYP = 5
585 SHDFAC = 0.8
586 ALBEDOK =0.2
587 ALBBRD =0.2
588 IF ( FRC_URB2D(I,J) < 0.99 ) THEN
589 T1= ( TSK(I,J) -FRC_URB2D(I,J) * TS_URB2D (I,J) )/ (1-FRC_URB2D(I,J))
590 ELSE
591 T1 = TSK(I,J)
592 ENDIF
593 ENDIF
594 ELSE
595 IF( IVGTYP(I,J) == 1 .or. IVGTYP(I,J) == 31 .or. &
596 IVGTYP(I,J) == 32 .or. IVGTYP(I,J) == 33) THEN
597 VEGTYP = 1
598 ENDIF
599 ENDIF
600
601 IF(IPRINT) THEN
602 !
603 print*, 'BEFORE SFLX, in Noahlsm_driver'
604 print*, 'ICE', ICE, 'DT',DT, 'ZLVL',ZLVL, 'NSOIL', NSOIL, &
605 'SLDPTH', SLDPTH, 'LOCAL',LOCAL, 'LUTYPE',&
606 LUTYPE, 'SLTYPE',SLTYPE, 'LWDN',LWDN, 'SOLDN',SOLDN, &
607 'SFCPRS',SFCPRS, 'PRCP',PRCP,'SFCTMP',SFCTMP,'Q2K',Q2K, &
608 'TH2',TH2,'Q2SAT',Q2SAT,'DQSDT2',DQSDT2,'VEGTYP', VEGTYP,&
609 'SOILTYP',SOILTYP, 'SLOPETYP',SLOPETYP, 'SHDFAC',SHDFAC,&
610 'SHMIN',SHMIN, 'ALBBRD',ALBBRD,'SNOALB1',SNOALB1,'TBOT',&
611 TBOT, 'Z0BRD',Z0BRD, 'Z0K',Z0K, 'CMC',CMC, 'T1',T1,'STC',&
612 STC, 'SMC',SMC, 'SWC',SWC,'SNOWHK',SNOWHK,'SNEQV',SNEQV,&
613 'ALBEDO',ALBEDO,'CHK',CHK,'ETA',ETA,'SHEAT',SHEAT, &
614 'ETA_KINEMATIC',ETA_KINEMATIC, 'FDOWN',FDOWN,'EC',EC, &
615 'EDIR',EDIR,'ET',ET,'ETT',ETT,'ESNOW',ESNOW,'DRIP',DRIP,&
616 'DEW',DEW,'BETA',BETA,'ETP',ETP,'SSOIL',SSOIL,'FLX1',FLX1,&
617 'FLX2',FLX2,'FLX3',FLX3,'SNOMLT',SNOMLT,'SNCOVR',SNCOVR,&
618 'RUNOFF1',RUNOFF1,'RUNOFF2',RUNOFF2,'RUNOFF3',RUNOFF3, &
619 'RC',RC, 'PC',PC,'RSMIN',RSMIN,'XLAI',XLAI,'RCS',RCS, &
620 'RCT',RCT,'RCQ',RCQ,'RCSOIL',RCSOIL,'SOILW',SOILW, &
621 'SOILM',SOILM,'Q1',Q1,'SMCWLT',SMCWLT,'SMCDRY',SMCDRY,&
622 'SMCREF',SMCREF,'SMCMAX',SMCMAX,'NROOT',NROOT
623 endif
624
625
626 CALL SFLX (ICE,DT,ZLVL,NSOIL,SLDPTH, & !C
627 LOCAL, & !L
628 LUTYPE, SLTYPE, & !CL
629 LWDN,SOLDN,SFCPRS,PRCP,SFCTMP,Q2K,DUMMY, & !F
630 DUMMY,DUMMY, DUMMY, & !F PRCPRAIN not used
631 TH2,Q2SAT,DQSDT2, & !I
632 VEGTYP,SOILTYP,SLOPETYP,SHDFAC,SHMIN,DUMMY, & !I
633 ALBBRD, SNOALB1,TBOT, Z0BRD, Z0K, EMISSI, & !S
634 CMC,T1,STC,SMC,SWC,SNOWHK,SNEQV,ALBEDOK,CHK,dummy,& !H
635 ETA,SHEAT, ETA_KINEMATIC,FDOWN, & !O
636 EC,EDIR,ET,ETT,ESNOW,DRIP,DEW, & !O
637 BETA,ETP,SSOIL, & !O
638 FLX1,FLX2,FLX3, & !O
639 SNOMLT,SNCOVR, & !O
640 RUNOFF1,RUNOFF2,RUNOFF3, & !O
641 RC,PC,RSMIN,XLAI,RCS,RCT,RCQ,RCSOIL, & !O
642 SOILW,SOILM,Q1, & !D
643 SMCWLT,SMCDRY,SMCREF,SMCMAX,NROOT)
644
645
646 IF(IPRINT) THEN
647
648 print*, 'AFTER SFLX, in Noahlsm_driver'
649 print*, 'ICE', ICE, 'DT',DT, 'ZLVL',ZLVL, 'NSOIL', NSOIL, &
650 'SLDPTH', SLDPTH, 'LOCAL',LOCAL, 'LUTYPE',&
651 LUTYPE, 'SLTYPE',SLTYPE, 'LWDN',LWDN, 'SOLDN',SOLDN, &
652 'SFCPRS',SFCPRS, 'PRCP',PRCP,'SFCTMP',SFCTMP,'Q2K',Q2K, &
653 'TH2',TH2,'Q2SAT',Q2SAT,'DQSDT2',DQSDT2,'VEGTYP', VEGTYP,&
654 'SOILTYP',SOILTYP, 'SLOPETYP',SLOPETYP, 'SHDFAC',SHDFAC,&
655 'SHDMIN',SHMIN, 'ALBBRD',ALBBRD,'SNOALB',SNOALB1,'TBOT',&
656 TBOT, 'Z0BRD',Z0BRD, 'Z0K',Z0K, 'CMC',CMC, 'T1',T1,'STC',&
657 STC, 'SMC',SMC, 'SWc',SWC,'SNOWHK',SNOWHK,'SNEQV',SNEQV,&
658 'ALBEDO',ALBEDOK,'CHK',CHK,'ETA',ETA,'SHEAT',SHEAT, &
659 'ETA_KINEMATIC',ETA_KINEMATIC, 'FDOWN',FDOWN,'EC',EC, &
660 'EDIR',EDIR,'ET',ET,'ETT',ETT,'ESNOW',ESNOW,'DRIP',DRIP,&
661 'DEW',DEW,'BETA',BETA,'ETP',ETP,'SSOIL',SSOIL,'FLX1',FLX1,&
662 'FLX2',FLX2,'FLX3',FLX3,'SNOMLT',SNOMLT,'SNCOVR',SNCOVR,&
663 'RUNOFF1',RUNOFF1,'RUNOFF2',RUNOFF2,'RUNOFF3',RUNOFF3, &
664 'RC',RC, 'PC',PC,'RSMIN',RSMIN,'XLAI',XLAI,'RCS',RCS, &
665 'RCT',RCT,'RCQ',RCQ,'RCSOIL',RCSOIL,'SOILW',SOILW, &
666 'SOILM',SOILM,'Q1',Q1,'SMCWLT',SMCWLT,'SMCDRY',SMCDRY,&
667 'SMCREF',SMCREF,'SMCMAX',SMCMAX,'NROOT',NROOT
668 endif
669
670 !*** UPDATE STATE VARIABLES
671 CANWAT(I,J)=CMC
672 SNOW(I,J)=SNEQV*1000.
673 ! SNOWH(I,J)=SNOWHK*1000.
674 SNOWH(I,J)=SNOWHK ! SNOWHK is assumed in meters
675 ALBEDO(I,J)=ALBEDOK
676 ZNT(I,J)=Z0K
677 TSK(I,J)=T1
678 HFX(I,J)=SHEAT
679 QFX(I,J)=ETA_KINEMATIC
680 LH(I,J)=ETA
681 GRDFLX(I,J)=SSOIL
682 SNOWC(I,J)=SNCOVR
683 CHS2(I,J)=CQS2(I,J)
684 ! prevent diagnostic ground q (q1) from being greater than qsat(tsk)
685 ! as happens over snow cover where the cqs2 value also becomes irrelevant
686 ! by setting cqs2=chs in this situation the 2m q should become just qv(k=1)
687 IF (Q1 .GT. QSFC(I,J)) THEN
688 CQS2(I,J) = CHS(I,J)
689 ENDIF
690 QSFC(I,J)=Q1
691 !
692 DO 80 NS=1,NSOIL
693 SMOIS(I,NS,J)=SMC(NS)
694 TSLB(I,NS,J)=STC(NS) ! STEMP
695 SH2O(I,NS,J)=SWC(NS)
696 80 CONTINUE
697 ! ENDIF
698
699 IF (UCMCALL == 1 ) THEN ! Beginning of UCM CALL if block
700 !--------------------------------------
701 ! URBAN CANOPY MODEL START - urban
702 !--------------------------------------
703 ! Input variables lsm --> urban
704
705
706 IF( IVGTYP(I,J) == 1 .or. IVGTYP(I,J) == 31 .or. &
707 IVGTYP(I,J) == 32 .or. IVGTYP(I,J) == 33 ) THEN
708
709 ! Call urban
710
711 !
712 UTYPE_URB = UTYPE_URB2D(I,J) !urban type (low, high or industrial)
713
714 TA_URB = SFCTMP ! [K]
715 QA_URB = Q2K ! [kg/kg]
716 UA_URB = SQRT(U_PHY(I,1,J)**2.+V_PHY(I,1,J)**2.)
717 U1_URB = U_PHY(I,1,J)
718 V1_URB = V_PHY(I,1,J)
719 IF(UA_URB < 1.) UA_URB=1. ! [m/s]
720 SSG_URB = SOLDN ! [W/m/m]
721 SSGD_URB = 0.8*SOLDN ! [W/m/m]
722 SSGQ_URB = SSG_URB-SSGD_URB ! [W/m/m]
723 LLG_URB = LWDN ! [W/m/m]
724 RAIN_URB = RAINBL(I,J) ! [mm]
725 RHOO_URB = SFCPRS / (287.04 * SFCTMP * (1.0+ 0.61 * Q2K)) ![kg/m/m/m]
726 ZA_URB = ZLVL ! [m]
727 DELT_URB = DT ! [sec]
728 XLAT_URB = XLAT_URB2D(I,J) ! [deg]
729 COSZ_URB = COSZ_URB2D(I,J) !
730 OMG_URB = OMG_URB2D(I,J) !
731 ZNT_URB = ZNT(I,J)
732
733 LSOLAR_URB = .FALSE.
734
735 TR_URB = TR_URB2D(I,J)
736 TB_URB = TB_URB2D(I,J)
737 TG_URB = TG_URB2D(I,J)
738 TC_URB = TC_URB2D(I,J)
739 QC_URB = QC_URB2D(I,J)
740 UC_URB = UC_URB2D(I,J)
741
742 DO K = 1,num_roof_layers
743 TRL_URB(K) = TRL_URB3D(I,K,J)
744 END DO
745 DO K = 1,num_wall_layers
746 TBL_URB(K) = TBL_URB3D(I,K,J)
747 END DO
748 DO K = 1,num_road_layers
749 TGL_URB(K) = TGL_URB3D(I,K,J)
750 END DO
751
752 XXXR_URB = XXXR_URB2D(I,J)
753 XXXB_URB = XXXB_URB2D(I,J)
754 XXXG_URB = XXXG_URB2D(I,J)
755 XXXC_URB = XXXC_URB2D(I,J)
756 !
757 CHS_URB = CHS(I,J)
758 CHS2_URB = CHS2(I,J)
759 !
760
761 ! Call urban
762
763
764 CALL urban(LSOLAR_URB, & ! I
765 num_roof_layers,num_wall_layers,num_road_layers, & ! C
766 DZR,DZB,DZG, & ! C
767 UTYPE_URB,TA_URB,QA_URB,UA_URB,U1_URB,V1_URB,SSG_URB, & ! I
768 SSGD_URB,SSGQ_URB,LLG_URB,RAIN_URB,RHOO_URB, & ! I
769 ZA_URB,DECLIN_URB,COSZ_URB,OMG_URB, & ! I
770 XLAT_URB,DELT_URB,ZNT_URB, & ! I
771 CHS_URB, CHS2_URB, & ! I
772 TR_URB, TB_URB, TG_URB, TC_URB, QC_URB,UC_URB, & ! H
773 TRL_URB,TBL_URB,TGL_URB, & ! H
774 XXXR_URB, XXXB_URB, XXXG_URB, XXXC_URB, & ! H
775 TS_URB,QS_URB,SH_URB,LH_URB,LH_KINEMATIC_URB, & ! O
776 SW_URB,ALB_URB,LW_URB,G_URB,RN_URB,PSIM_URB,PSIH_URB, & ! O
777 GZ1OZ0_URB, & !O
778 U10_URB, V10_URB, TH2_URB, Q2_URB, & ! O
779 UST_URB) !O
780
781
782 IF(IPRINT) THEN
783
784 print*, 'AFTER CALL URBAN'
785 print*,'num_roof_layers',num_roof_layers, 'num_wall_layers', &
786 num_wall_layers, &
787 'DZR',DZR,'DZB',DZB,'DZG',DZG,'UTYPE_URB',UTYPE_URB,'TA_URB', &
788 TA_URB, &
789 'QA_URB',QA_URB,'UA_URB',UA_URB,'U1_URB',U1_URB,'V1_URB', &
790 V1_URB, &
791 'SSG_URB',SSG_URB,'SSGD_URB',SSGD_URB,'SSGQ_URB',SSGQ_URB, &
792 'LLG_URB',LLG_URB,'RAIN_URB',RAIN_URB,'RHOO_URB',RHOO_URB, &
793 'ZA_URB',ZA_URB, 'DECLIN_URB',DECLIN_URB,'COSZ_URB',COSZ_URB,&
794 'OMG_URB',OMG_URB,'XLAT_URB',XLAT_URB,'DELT_URB',DELT_URB, &
795 'ZNT_URB',ZNT_URB,'TR_URB',TR_URB, 'TB_URB',TB_URB,'TG_URB',&
796 TG_URB,'TC_URB',TC_URB,'QC_URB',QC_URB,'TRL_URB',TRL_URB, &
797 'TBL_URB',TBL_URB,'TGL_URB',TGL_URB,'XXXR_URB',XXXR_URB, &
798 'XXXB_URB',XXXB_URB,'XXXG_URB',XXXG_URB,'XXXC_URB',XXXC_URB,&
799 'TS_URB',TS_URB,'QS_URB',QS_URB,'SH_URB',SH_URB,'LH_URB', &
800 LH_URB, 'LH_KINEMATIC_URB',LH_KINEMATIC_URB,'SW_URB',SW_URB,&
801 'ALB_URB',ALB_URB,'LW_URB',LW_URB,'G_URB',G_URB,'RN_URB', &
802 RN_URB, 'PSIM_URB',PSIM_URB,'PSIH_URB',PSIH_URB, &
803 'U10_URB',U10_URB,'V10_URB',V10_URB,'TH2_URB',TH2_URB, &
804 'Q2_URB',Q2_URB,'CHS_URB',CHS_URB,'CHS2_URB',CHS2_URB
805 endif
806
807 TS_URB2D(I,J) = TS_URB
808
809 ALBEDO(I,J) = FRC_URB2D(I,J)*ALB_URB+(1-FRC_URB2D(I,J))*ALBEDOK ![-]
810 HFX(I,J) = FRC_URB2D(I,J)*SH_URB+(1-FRC_URB2D(I,J))*SHEAT ![W/m/m]
811 QFX(I,J) = FRC_URB2D(I,J)*LH_KINEMATIC_URB &
812 + (1-FRC_URB2D(I,J))*ETA_KINEMATIC ![kg/m/m/s]
813 LH(I,J) = FRC_URB2D(I,J)*LH_URB+(1-FRC_URB2D(I,J))*ETA ![W/m/m]
814 GRDFLX(I,J) = FRC_URB2D(I,J)*G_URB+(1-FRC_URB2D(I,J))*SSOIL ![W/m/m]
815 TSK(I,J) = FRC_URB2D(I,J)*TS_URB+(1-FRC_URB2D(I,J))*T1 ![K]
816 QSFC(I,J)= FRC_URB2D(I,J)*QS_URB+(1-FRC_URB2D(I,J))*Q1 ![-]
817
818 IF(IPRINT)THEN
819
820 print*, ' FRC_URB2D', FRC_URB2D, &
821 'ALB_URB',ALB_URB, 'ALBEDOK',ALBEDOK, &
822 'ALBEDO(I,J)', ALBEDO(I,J), &
823 'SH_URB',SH_URB,'SHEAT',SHEAT, 'HFX(I,J)',HFX(I,J), &
824 'LH_KINEMATIC_URB',LH_KINEMATIC_URB,'ETA_KINEMATIC', &
825 ETA_KINEMATIC, 'QFX(I,J)',QFX(I,J), &
826 'LH_URB',LH_URB, 'ETA',ETA, 'LH(I,J)',LH(I,J), &
827 'G_URB',G_URB,'SSOIL',SSOIL,'GRDFLX(I,J)', GRDFLX(I,J),&
828 'TS_URB',TS_URB,'T1',T1,'TSK(I,J)',TSK(I,J), &
829 'QS_URB',QS_URB,'Q1',Q1,'QSFC(I,J)',QSFC(I,J)
830 endif
831
832
833
834
835 ! Renew Urban State Varialbes
836
837 TR_URB2D(I,J) = TR_URB
838 TB_URB2D(I,J) = TB_URB
839 TG_URB2D(I,J) = TG_URB
840 TC_URB2D(I,J) = TC_URB
841 QC_URB2D(I,J) = QC_URB
842 UC_URB2D(I,J) = UC_URB
843
844 DO K = 1,num_roof_layers
845 TRL_URB3D(I,K,J) = TRL_URB(K)
846 END DO
847 DO K = 1,num_wall_layers
848 TBL_URB3D(I,K,J) = TBL_URB(K)
849 END DO
850 DO K = 1,num_road_layers
851 TGL_URB3D(I,K,J) = TGL_URB(K)
852 END DO
853 XXXR_URB2D(I,J) = XXXR_URB
854 XXXB_URB2D(I,J) = XXXB_URB
855 XXXG_URB2D(I,J) = XXXG_URB
856 XXXC_URB2D(I,J) = XXXC_URB
857
858 SH_URB2D(I,J) = SH_URB
859 LH_URB2D(I,J) = LH_URB
860 G_URB2D(I,J) = G_URB
861 RN_URB2D(I,J) = RN_URB
862 PSIM_URB2D(I,J) = PSIM_URB
863 PSIH_URB2D(I,J) = PSIH_URB
864 GZ1OZ0_URB2D(I,J)= GZ1OZ0_URB
865 U10_URB2D(I,J) = U10_URB
866 V10_URB2D(I,J) = V10_URB
867 TH2_URB2D(I,J) = TH2_URB
868 Q2_URB2D(I,J) = Q2_URB
869 UST_URB2D(I,J) = UST_URB
870 AKMS_URB2D(I,J) = KARMAN * UST_URB2D(I,J)/(GZ1OZ0_URB2D(I,J)-PSIM_URB2D(I,J))
871
872 END IF
873
874 ENDIF ! end of UCM CALL if block
875 !--------------------------------------
876 ! Urban Part End - urban
877 !--------------------------------------
878
879 !*** DIAGNOSTICS
880 SMSTAV(I,J)=SOILW
881 SMSTOT(I,J)=SOILM*1000.
882 ! Convert the water unit into mm
883 SFCRUNOFF(I,J)=SFCRUNOFF(I,J)+RUNOFF1*DT*1000.0
884 UDRUNOFF(I,J)=UDRUNOFF(I,J)+(RUNOFF2+RUNOFF3)*DT*1000.0
885 IF(SFCTMP<=T0)THEN
886 ! ACSNOW(I,J)=ACSNOW(I,J)+PREC(I)*DT
887 ACSNOW(I,J)=ACSNOW(I,J)+PRCP*DT
888 ENDIF
889 IF(SNOW(I,J).GT.0.)THEN
890 ACSNOM(I,J)=ACSNOM(I,J)+SNOMLT*1000.
891 ENDIF
892
893 ENDIF ! endif of land-sea test
894
895 100 CONTINUE ! of I loop
896
897 ENDDO ! of J loop
898 !------------------------------------------------------
899 END SUBROUTINE lsm
900 !------------------------------------------------------
901
902 SUBROUTINE LSMINIT(VEGFRA,SNOW,SNOWC,SNOWH,CANWAT,SMSTAV, &
903 SMSTOT, SFCRUNOFF,UDRUNOFF,ACSNOW, &
904 ACSNOM,IVGTYP,ISLTYP,TSLB,SMOIS,SH2O,ZS,DZS, &
905 FNDSOILW, FNDSNOWH, &
906 num_soil_layers, restart, &
907 allowed_to_read , &
908 ids,ide, jds,jde, kds,kde, &
909 ims,ime, jms,jme, kms,kme, &
910 its,ite, jts,jte, kts,kte )
911
912 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
913 ims,ime, jms,jme, kms,kme, &
914 its,ite, jts,jte, kts,kte
915
916 INTEGER, INTENT(IN) :: num_soil_layers
917
918 LOGICAL , INTENT(IN) :: restart , allowed_to_read
919
920 REAL, DIMENSION( num_soil_layers), INTENT(INOUT) :: ZS, DZS
921
922 REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , &
923 INTENT(INOUT) :: SMOIS, & !Total soil moisture
924 SH2O, & !liquid soil moisture
925 TSLB !STEMP
926
927 REAL, DIMENSION( ims:ime, jms:jme ) , &
928 INTENT(INOUT) :: SNOW, &
929 SNOWH, &
930 SNOWC, &
931 CANWAT, &
932 SMSTAV, &
933 SMSTOT, &
934 SFCRUNOFF, &
935 UDRUNOFF, &
936 ACSNOW, &
937 VEGFRA, &
938 ACSNOM
939
940 INTEGER, DIMENSION( ims:ime, jms:jme ) , &
941 INTENT(IN) :: IVGTYP, &
942 ISLTYP
943
944 LOGICAL, INTENT(IN) :: FNDSOILW , &
945 FNDSNOWH
946
947 INTEGER :: L
948 REAL :: BX, SMCMAX, PSISAT, FREE
949 REAL, PARAMETER :: BLIM = 5.5, HLICE = 3.335E5, &
950 GRAV = 9.81, T0 = 273.15
951 INTEGER :: errflag
952
953 !
954
955
956 ! initialize three Noah LSM related tables
957 IF ( allowed_to_read ) THEN
958 CALL wrf_message( 'INITIALIZE THREE Noah LSM RELATED TABLES' )
959 CALL LSM_PARM_INIT
960 ENDIF
961
962 IF(.not.restart)THEN
963
964 itf=min0(ite,ide-1)
965 jtf=min0(jte,jde-1)
966
967 errflag = 0
968 DO j = jts,jtf
969 DO i = its,itf
970 IF ( ISLTYP( i,j ) .LT. 1 ) THEN
971 errflag = 1
972 WRITE(err_message,*)"module_sf_noahlsm.F: lsminit: out of range ISLTYP ",i,j,ISLTYP( i,j )
973 CALL wrf_message(err_message)
974 ENDIF
975 ENDDO
976 ENDDO
977 IF ( errflag .EQ. 1 ) THEN
978 CALL wrf_error_fatal( "module_sf_noahlsm.F: lsminit: out of range value "// &
979 "of ISLTYP. Is this field in the input?" )
980 ENDIF
981
982 ! initialize soil liquid water content SH2O
983
984 ! IF(.NOT.FNDSOILW) THEN
985
986 ! If no SWC, do the following
987 ! PRINT *,'SOIL WATER NOT FOUND - VALUE SET IN LSMINIT'
988 DO J = jts,jtf
989 DO I = its,itf
990 BX = BB(ISLTYP(I,J))
991 SMCMAX = MAXSMC(ISLTYP(I,J))
992 PSISAT = SATPSI(ISLTYP(I,J))
993 if ((bx > 0.0).and.(smcmax > 0.0).and.(psisat > 0.0)) then
994 DO NS=1, num_soil_layers
995 ! ----------------------------------------------------------------------
996 !SH2O <= SMOIS for T < 273.149K (-0.001C)
997 IF (TSLB(I,NS,J) < 273.149) THEN
998 ! ----------------------------------------------------------------------
999 ! first guess following explicit solution for Flerchinger Eqn from Koren
1000 ! et al, JGR, 1999, Eqn 17 (KCOUNT=0 in FUNCTION FRH2O).
1001 ! ISLTPK is soil type
1002 BX = BB(ISLTYP(I,J))
1003 SMCMAX = MAXSMC(ISLTYP(I,J))
1004 PSISAT = SATPSI(ISLTYP(I,J))
1005 IF ( BX > BLIM ) BX = BLIM
1006 FK=(( (HLICE/(GRAV*(-PSISAT))) * &
1007 ((TSLB(I,NS,J)-T0)/TSLB(I,NS,J)) )**(-1/BX) )*SMCMAX
1008 IF (FK < 0.02) FK = 0.02
1009 SH2O(I,NS,J) = MIN( FK, SMOIS(I,NS,J) )
1010 ! ----------------------------------------------------------------------
1011 ! now use iterative solution for liquid soil water content using
1012 ! FUNCTION FRH2O with the initial guess for SH2O from above explicit
1013 ! first guess.
1014 CALL FRH2O (FREE,TSLB(I,NS,J),SMOIS(I,NS,J),SH2O(I,NS,J), &
1015 SMCMAX,BX,PSISAT)
1016 SH2O(I,NS,J) = FREE
1017 ELSE ! of IF (TSLB(I,NS,J)
1018 ! ----------------------------------------------------------------------
1019 ! SH2O = SMOIS ( for T => 273.149K (-0.001C)
1020 SH2O(I,NS,J)=SMOIS(I,NS,J)
1021 ! ----------------------------------------------------------------------
1022 ENDIF ! of IF (TSLB(I,NS,J)
1023 END DO ! of DO NS=1, num_soil_layers
1024 else ! of if ((bx > 0.0)
1025 DO NS=1, num_soil_layers
1026 SH2O(I,NS,J)=SMOIS(I,NS,J)
1027 END DO
1028 endif ! of if ((bx > 0.0)
1029 ENDDO ! DO I = its,itf
1030 ENDDO ! DO J = jts,jtf
1031 ! ENDIF ! of IF(.NOT.FNDSOILW)THEN
1032
1033 ! initialize physical snow height SNOWH
1034
1035 IF(.NOT.FNDSNOWH)THEN
1036 ! If no SNOWH do the following
1037 CALL wrf_message( 'SNOW HEIGHT NOT FOUND - VALUE DEFINED IN LSMINIT' )
1038 DO J = jts,jtf
1039 DO I = its,itf
1040 SNOWH(I,J)=SNOW(I,J)*0.005 ! SNOW in mm and SNOWH in m
1041 ENDDO
1042 ENDDO
1043 ENDIF
1044
1045 ! initialize canopy water to ZERO
1046
1047 ! GO TO 110
1048 ! print*,'Note that canopy water content (CANWAT) is set to ZERO in LSMINIT'
1049 DO J = jts,jtf
1050 DO I = its,itf
1051 CANWAT(I,J)=0.0
1052 ENDDO
1053 ENDDO
1054 110 CONTINUE
1055
1056 ENDIF
1057 !------------------------------------------------------------------------------
1058 END SUBROUTINE lsminit
1059 !------------------------------------------------------------------------------
1060
1061
1062
1063 !
1064 !-----------------------------------------------------------------
1065 SUBROUTINE LSM_PARM_INIT
1066 !-----------------------------------------------------------------
1067
1068 character*4 :: MMINLU, MMINSL
1069
1070 MMINLU='USGS'
1071 MMINSL='STAS'
1072 call SOIL_VEG_GEN_PARM( MMINLU, MMINSL)
1073
1074 !-----------------------------------------------------------------
1075 END SUBROUTINE LSM_PARM_INIT
1076 !-----------------------------------------------------------------
1077
1078 !-----------------------------------------------------------------
1079 SUBROUTINE SOIL_VEG_GEN_PARM( MMINLU, MMINSL)
1080 !-----------------------------------------------------------------
1081
1082 IMPLICIT NONE
1083
1084 integer :: LUMATCH, IINDEX, LC, NUM_SLOPE
1085 integer :: ierr
1086 INTEGER , PARAMETER :: OPEN_OK = 0
1087
1088 character*4 :: MMINLU, MMINSL
1089 character*128 :: mess , message
1090 logical, external :: wrf_dm_on_monitor
1091
1092
1093 !-----SPECIFY VEGETATION RELATED CHARACTERISTICS :
1094 ! ALBBCK: SFC albedo (in percentage)
1095 ! Z0: Roughness length (m)
1096 ! SHDFAC: Green vegetation fraction (in percentage)
1097 ! Note: The ALBEDO, Z0, and SHDFAC values read from the following table
1098 ! ALBEDO, amd Z0 are specified in LAND-USE TABLE; and SHDFAC is
1099 ! the monthly green vegetation data
1100 ! CMXTBL: MAX CNPY Capacity (m)
1101 ! NROTBL: Rooting depth (layer)
1102 ! RSMIN: Mimimum stomatal resistance (s m-1)
1103 ! RSMAX: Max. stomatal resistance (s m-1)
1104 ! RGL: Parameters used in radiation stress function
1105 ! HS: Parameter used in vapor pressure deficit functio
1106 ! TOPT: Optimum transpiration air temperature. (K)
1107 ! CMCMAX: Maximum canopy water capacity
1108 ! CFACTR: Parameter used in the canopy inteception calculati
1109 ! SNUP: Threshold snow depth (in water equivalent m) that
1110 ! implies 100% snow cover
1111 ! LAI: Leaf area index (dimensionless)
1112 ! MAXALB: Upper bound on maximum albedo over deep snow
1113 !
1114 !-----READ IN VEGETAION PROPERTIES FROM VEGPARM.TBL
1115 !
1116
1117 IF ( wrf_dm_on_monitor() ) THEN
1118
1119 OPEN(19, FILE='VEGPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
1120 IF(ierr .NE. OPEN_OK ) THEN
1121 WRITE(message,FMT='(A)') &
1122 'module_sf_noahlsm.F: soil_veg_gen_parm: failure opening VEGPARM.TBL'
1123 CALL wrf_error_fatal ( message )
1124 END IF
1125
1126 WRITE ( mess, * ) 'INPUT LANDUSE = ',MMINLU
1127 CALL wrf_message( mess )
1128
1129 LUMATCH=0
1130
1131 READ (19,*)
1132 READ (19,2000,END=2002)LUTYPE
1133 READ (19,*)LUCATS,IINDEX
1134 2000 FORMAT (A4)
1135
1136 IF(LUTYPE.EQ.MMINLU)THEN
1137 WRITE( mess , * ) 'LANDUSE TYPE = ',LUTYPE,' FOUND', &
1138 LUCATS,' CATEGORIES'
1139 CALL wrf_message( mess )
1140 LUMATCH=1
1141 ENDIF
1142
1143 IF(LUTYPE.EQ.MMINLU)THEN
1144 DO LC=1,LUCATS
1145 READ (19,*)IINDEX,ALBTBL(LC),Z0TBL(LC),SHDTBL(LC), &
1146 NROTBL(LC),RSTBL(LC),RGLTBL(LC),HSTBL(LC), &
1147 SNUPTBL(LC),LAITBL(LC),MAXALB(LC)
1148 ENDDO
1149 !
1150 READ (19,*)
1151 READ (19,*)TOPT_DATA
1152 READ (19,*)
1153 READ (19,*)CMCMAX_DATA
1154 READ (19,*)
1155 READ (19,*)CFACTR_DATA
1156 READ (19,*)
1157 READ (19,*)RSMAX_DATA
1158 READ (19,*)
1159 READ (19,*)BARE
1160 ENDIF
1161 !
1162 2002 CONTINUE
1163
1164 CLOSE (19)
1165 ENDIF
1166
1167 CALL wrf_dm_bcast_string ( LUTYPE , 4 )
1168 CALL wrf_dm_bcast_integer ( LUCATS , 1 )
1169 CALL wrf_dm_bcast_integer ( IINDEX , 1 )
1170 CALL wrf_dm_bcast_integer ( LUMATCH , 1 )
1171 CALL wrf_dm_bcast_real ( ALBTBL , NLUS )
1172 CALL wrf_dm_bcast_real ( Z0TBL , NLUS )
1173 CALL wrf_dm_bcast_real ( SHDTBL , NLUS )
1174 CALL wrf_dm_bcast_real ( NROTBL , NLUS )
1175 CALL wrf_dm_bcast_real ( RSTBL , NLUS )
1176 CALL wrf_dm_bcast_real ( RGLTBL , NLUS )
1177 CALL wrf_dm_bcast_real ( HSTBL , NLUS )
1178 CALL wrf_dm_bcast_real ( SNUPTBL , NLUS )
1179 CALL wrf_dm_bcast_real ( LAITBL , NLUS )
1180 CALL wrf_dm_bcast_real ( MAXALB , NLUS )
1181 CALL wrf_dm_bcast_real ( TOPT_DATA , 1 )
1182 CALL wrf_dm_bcast_real ( CMCMAX_DATA , 1 )
1183 CALL wrf_dm_bcast_real ( CFACTR_DATA , 1 )
1184 CALL wrf_dm_bcast_real ( RSMAX_DATA , 1 )
1185 CALL wrf_dm_bcast_integer ( BARE , 1 )
1186
1187 !
1188 !-----READ IN SOIL PROPERTIES FROM SOILPARM.TBL
1189 !
1190 IF ( wrf_dm_on_monitor() ) THEN
1191 OPEN(19, FILE='SOILPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
1192 IF(ierr .NE. OPEN_OK ) THEN
1193 WRITE(message,FMT='(A)') &
1194 'module_sf_noahlsm.F: soil_veg_gen_parm: failure opening SOILPARM.TBL'
1195 CALL wrf_error_fatal ( message )
1196 END IF
1197
1198 MMINSL='STAS' !oct2
1199 WRITE(mess,*) 'INPUT SOIL TEXTURE CLASSIFICAION = ',MMINSL
1200 CALL wrf_message( mess )
1201
1202 LUMATCH=0
1203
1204 READ (19,*)
1205 READ (19,2000,END=2003)SLTYPE
1206 READ (19,*)SLCATS,IINDEX
1207 IF(SLTYPE.EQ.MMINSL)THEN
1208 WRITE( mess , * ) 'SOIL TEXTURE CLASSIFICATION = ',SLTYPE,' FOUND', &
1209 SLCATS,' CATEGORIES'
1210 CALL wrf_message ( mess )
1211 LUMATCH=1
1212 ENDIF
1213 IF(SLTYPE.EQ.MMINSL)THEN
1214 DO LC=1,SLCATS
1215 READ (19,*) IINDEX,BB(LC),DRYSMC(LC),F11(LC),MAXSMC(LC),&
1216 REFSMC(LC),SATPSI(LC),SATDK(LC), SATDW(LC), &
1217 WLTSMC(LC), QTZ(LC)
1218 ENDDO
1219 ENDIF
1220
1221 2003 CONTINUE
1222
1223 CLOSE (19)
1224 ENDIF
1225
1226 CALL wrf_dm_bcast_integer ( LUMATCH , 1 )
1227 CALL wrf_dm_bcast_string ( SLTYPE , 4 )
1228 CALL wrf_dm_bcast_string ( MMINSL , 4 ) ! since this is reset above, see oct2 ^
1229 CALL wrf_dm_bcast_integer ( SLCATS , 1 )
1230 CALL wrf_dm_bcast_integer ( IINDEX , 1 )
1231 CALL wrf_dm_bcast_real ( BB , NSLTYPE )
1232 CALL wrf_dm_bcast_real ( DRYSMC , NSLTYPE )
1233 CALL wrf_dm_bcast_real ( F11 , NSLTYPE )
1234 CALL wrf_dm_bcast_real ( MAXSMC , NSLTYPE )
1235 CALL wrf_dm_bcast_real ( REFSMC , NSLTYPE )
1236 CALL wrf_dm_bcast_real ( SATPSI , NSLTYPE )
1237 CALL wrf_dm_bcast_real ( SATDK , NSLTYPE )
1238 CALL wrf_dm_bcast_real ( SATDW , NSLTYPE )
1239 CALL wrf_dm_bcast_real ( WLTSMC , NSLTYPE )
1240 CALL wrf_dm_bcast_real ( QTZ , NSLTYPE )
1241
1242 IF(LUMATCH.EQ.0)THEN
1243 CALL wrf_message( 'SOIl TEXTURE IN INPUT FILE DOES NOT ' )
1244 CALL wrf_message( 'MATCH SOILPARM TABLE' )
1245 CALL wrf_error_fatal ( 'INCONSISTENT OR MISSING SOILPARM FILE' )
1246 ENDIF
1247
1248 !
1249 !-----READ IN GENERAL PARAMETERS FROM GENPARM.TBL
1250 !
1251 IF ( wrf_dm_on_monitor() ) THEN
1252 OPEN(19, FILE='GENPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
1253 IF(ierr .NE. OPEN_OK ) THEN
1254 WRITE(message,FMT='(A)') &
1255 'module_sf_noahlsm.F: soil_veg_gen_parm: failure opening GENPARM.TBL'
1256 CALL wrf_error_fatal ( message )
1257 END IF
1258
1259 READ (19,*)
1260 READ (19,*)
1261 READ (19,*) NUM_SLOPE
1262
1263 SLPCATS=NUM_SLOPE
1264
1265 DO LC=1,SLPCATS
1266 READ (19,*)SLOPE_DATA(LC)
1267 ENDDO
1268
1269 READ (19,*)
1270 READ (19,*)SBETA_DATA
1271 READ (19,*)
1272 READ (19,*)FXEXP_DATA
1273 READ (19,*)
1274 READ (19,*)CSOIL_DATA
1275 READ (19,*)
1276 READ (19,*)SALP_DATA
1277 READ (19,*)
1278 READ (19,*)REFDK_DATA
1279 READ (19,*)
1280 READ (19,*)REFKDT_DATA
1281 READ (19,*)
1282 READ (19,*)FRZK_DATA
1283 READ (19,*)
1284 READ (19,*)ZBOT_DATA
1285 READ (19,*)
1286 READ (19,*)CZIL_DATA
1287 READ (19,*)
1288 READ (19,*)SMLOW_DATA
1289 READ (19,*)
1290 READ (19,*)SMHIGH_DATA
1291 CLOSE (19)
1292 ENDIF
1293
1294 CALL wrf_dm_bcast_integer ( NUM_SLOPE , 1 )
1295 CALL wrf_dm_bcast_integer ( SLPCATS , 1 )
1296 CALL wrf_dm_bcast_real ( SLOPE_DATA , NSLOPE )
1297 CALL wrf_dm_bcast_real ( SBETA_DATA , 1 )
1298 CALL wrf_dm_bcast_real ( FXEXP_DATA , 1 )
1299 CALL wrf_dm_bcast_real ( CSOIL_DATA , 1 )
1300 CALL wrf_dm_bcast_real ( SALP_DATA , 1 )
1301 CALL wrf_dm_bcast_real ( REFDK_DATA , 1 )
1302 CALL wrf_dm_bcast_real ( REFKDT_DATA , 1 )
1303 CALL wrf_dm_bcast_real ( FRZK_DATA , 1 )
1304 CALL wrf_dm_bcast_real ( ZBOT_DATA , 1 )
1305 CALL wrf_dm_bcast_real ( CZIL_DATA , 1 )
1306 CALL wrf_dm_bcast_real ( SMLOW_DATA , 1 )
1307 CALL wrf_dm_bcast_real ( SMHIGH_DATA , 1 )
1308
1309
1310 !-----------------------------------------------------------------
1311 END SUBROUTINE SOIL_VEG_GEN_PARM
1312 !-----------------------------------------------------------------
1313
1314 SUBROUTINE SFLX (ICE,DT,ZLVL,NSOIL,SLDPTH, & !C
1315 LOCAL, & !L
1316 LLANDUSE, LSOIL, & !CL
1317 LWDN,SOLDN,SFCPRS,PRCP,SFCTMP,Q2,SFCSPD, & !F
1318 COSZ,PRCPRAIN, SOLARDIRECT, & !F
1319 TH2,Q2SAT,DQSDT2, & !I
1320 VEGTYP,SOILTYP,SLOPETYP,SHDFAC,SHDMIN,SHDMAX, & !I
1321 ALB, SNOALB,TBOT, Z0BRD, Z0, EMISSI, & !S
1322 CMC,T1,STC,SMC,SH2O,SNOWH,SNEQV,ALBEDO,CH,CM, & !H
1323 ! ----------------------------------------------------------------------
1324 ! OUTPUTS, DIAGNOSTICS, PARAMETERS BELOW GENERALLY NOT NECESSARY WHEN
1325 ! COUPLED WITH E.G. A NWP MODEL (SUCH AS THE NOAA/NWS/NCEP MESOSCALE ETA
1326 ! MODEL). OTHER APPLICATIONS MAY REQUIRE DIFFERENT OUTPUT VARIABLES.
1327 ! ----------------------------------------------------------------------
1328 ETA,SHEAT, ETA_KINEMATIC,FDOWN, & !O
1329 EC,EDIR,ET,ETT,ESNOW,DRIP,DEW, & !O
1330 BETA,ETP,SSOIL, & !O
1331 FLX1,FLX2,FLX3, & !O
1332 SNOMLT,SNCOVR, & !O
1333 RUNOFF1,RUNOFF2,RUNOFF3, & !O
1334 RC,PC,RSMIN,XLAI,RCS,RCT,RCQ,RCSOIL, & !O
1335 SOILW,SOILM,Q1, & !D
1336 SMCWLT,SMCDRY,SMCREF,SMCMAX,NROOT) !P
1337 ! ----------------------------------------------------------------------
1338 ! SUBROUTINE SFLX - VERSION 3.X - October 2002
1339 ! ----------------------------------------------------------------------
1340 ! SUB-DRIVER FOR "Noah LSM" FAMILY OF PHYSICS SUBROUTINES FOR A
1341 ! SOIL/VEG/SNOWPACK LAND-SURFACE MODEL TO UPDATE SOIL MOISTURE, SOIL
1342 ! ICE, SOIL TEMPERATURE, SKIN TEMPERATURE, SNOWPACK WATER CONTENT,
1343 ! SNOWDEPTH, AND ALL TERMS OF THE SURFACE ENERGY BALANCE AND SURFACE
1344 ! WATER BALANCE (EXCLUDING INPUT ATMOSPHERIC FORCINGS OF DOWNWARD
1345 ! RADIATION AND PRECIP)
1346 ! ----------------------------------------------------------------------
1347 ! SFLX ARGUMENT LIST KEY:
1348 ! ----------------------------------------------------------------------
1349 ! C CONFIGURATION INFORMATION
1350 ! L LOGICAL
1351 ! CL 4-string character bearing logical meaning
1352 ! F FORCING DATA
1353 ! I OTHER (INPUT) FORCING DATA
1354 ! S SURFACE CHARACTERISTICS
1355 ! H HISTORY (STATE) VARIABLES
1356 ! O OUTPUT VARIABLES
1357 ! D DIAGNOSTIC OUTPUT
1358 ! P Parameters
1359 ! Msic Miscellaneous terms passed from gridded driver
1360 ! ----------------------------------------------------------------------
1361 ! 1. CONFIGURATION INFORMATION (C):
1362 ! ----------------------------------------------------------------------
1363 ! ICE SEA-ICE FLAG (=1: SEA-ICE, =0: LAND)
1364 ! DT TIMESTEP (SEC) (DT SHOULD NOT EXCEED 3600 SECS, RECOMMEND
1365 ! 1800 SECS OR LESS)
1366 ! ZLVL HEIGHT (M) ABOVE GROUND OF ATMOSPHERIC FORCING VARIABLES
1367 ! NSOIL NUMBER OF SOIL LAYERS (AT LEAST 2, AND NOT GREATER THAN
1368 ! PARAMETER NSOLD SET BELOW)
1369 ! SLDPTH THE THICKNESS OF EACH SOIL LAYER (M)
1370 ! ----------------------------------------------------------------------
1371 ! 2. LOGICAL:
1372 ! ----------------------------------------------------------------------
1373 ! LCH Exchange coefficient (Ch) calculation flag (false: using
1374 ! ch-routine SFCDIF; true: Ch is brought in)
1375 ! LOCAL Flag for local-site simulation (where there is no
1376 ! maps for albedo, veg fraction, and roughness
1377 ! true: all LSM parameters (inluding albedo, veg fraction and
1378 ! roughness length) will be defined by three tables
1379 ! LLANDUSE (=USGS, using USGS landuse classification)
1380 ! LSOIL (=STAS, using FAO/STATSGO soil texture classification)
1381 ! ----------------------------------------------------------------------
1382 ! 3. FORCING DATA (F):
1383 ! ----------------------------------------------------------------------
1384 ! LWDN LW DOWNWARD RADIATION (W M-2; POSITIVE, NOT NET LONGWAVE)
1385 ! SOLDN SOLAR DOWNWARD RADIATION (W M-2; POSITIVE, NOT NET SOLAR)
1386 ! SFCPRS PRESSURE AT HEIGHT ZLVL ABOVE GROUND (PASCALS)
1387 ! PRCP PRECIP RATE (KG M-2 S-1) (NOTE, THIS IS A RATE)
1388 ! SFCTMP AIR TEMPERATURE (K) AT HEIGHT ZLVL ABOVE GROUND
1389 ! TH2 AIR POTENTIAL TEMPERATURE (K) AT HEIGHT ZLVL ABOVE GROUND
1390 ! Q2 MIXING RATIO AT HEIGHT ZLVL ABOVE GROUND (KG KG-1)
1391 ! COSZ Solar zenith angle (not used for now)
1392 ! PRCPRAIN Liquid-precipitation rate (KG M-2 S-1) (not used)
1393 ! SOLARDIRECT Direct component of downward solar radiation (W M-2) (not used)
1394 ! ----------------------------------------------------------------------
1395 ! 4. OTHER FORCING (INPUT) DATA (I):
1396 ! ----------------------------------------------------------------------
1397 ! SFCSPD WIND SPEED (M S-1) AT HEIGHT ZLVL ABOVE GROUND
1398 ! Q2SAT SAT MIXING RATIO AT HEIGHT ZLVL ABOVE GROUND (KG KG-1)
1399 ! DQSDT2 SLOPE OF SAT SPECIFIC HUMIDITY CURVE AT T=SFCTMP
1400 ! (KG KG-1 K-1)
1401 ! ----------------------------------------------------------------------
1402 ! 5. CANOPY/SOIL CHARACTERISTICS (S):
1403 ! ----------------------------------------------------------------------
1404 ! VEGTYP VEGETATION TYPE (INTEGER INDEX)
1405 ! SOILTYP SOIL TYPE (INTEGER INDEX)
1406 ! SLOPETYP CLASS OF SFC SLOPE (INTEGER INDEX)
1407 ! SHDFAC AREAL FRACTIONAL COVERAGE OF GREEN VEGETATION
1408 ! (FRACTION= 0.0-1.0)
1409 ! SHDMIN MINIMUM AREAL FRACTIONAL COVERAGE OF GREEN VEGETATION
1410 ! (FRACTION= 0.0-1.0) <= SHDFAC
1411 ! PTU PHOTO THERMAL UNIT (PLANT PHENOLOGY FOR ANNUALS/CROPS)
1412 ! (NOT YET USED, BUT PASSED TO REDPRM FOR FUTURE USE IN
1413 ! VEG PARMS)
1414 ! ALB BACKROUND SNOW-FREE SURFACE ALBEDO (FRACTION), FOR JULIAN
1415 ! DAY OF YEAR (USUALLY FROM TEMPORAL INTERPOLATION OF
1416 ! MONTHLY MEAN VALUES' CALLING PROG MAY OR MAY NOT
1417 ! INCLUDE DIURNAL SUN ANGLE EFFECT)
1418 ! SNOALB UPPER BOUND ON MAXIMUM ALBEDO OVER DEEP SNOW (E.G. FROM
1419 ! ROBINSON AND KUKLA, 1985, J. CLIM. & APPL. METEOR.)
1420 ! TBOT BOTTOM SOIL TEMPERATURE (LOCAL YEARLY-MEAN SFC AIR
1421 ! TEMPERATURE)
1422 ! Z0BRD Background fixed roughness length (M)
1423 ! Z0 Time varying roughness length (M) as function of snow depth
1424 ! ----------------------------------------------------------------------
1425 ! 6. HISTORY (STATE) VARIABLES (H):
1426 ! ----------------------------------------------------------------------
1427 ! CMC CANOPY MOISTURE CONTENT (M)
1428 ! T1 GROUND/CANOPY/SNOWPACK) EFFECTIVE SKIN TEMPERATURE (K)
1429 ! STC(NSOIL) SOIL TEMP (K)
1430 ! SMC(NSOIL) TOTAL SOIL MOISTURE CONTENT (VOLUMETRIC FRACTION)
1431 ! SH2O(NSOIL) UNFROZEN SOIL MOISTURE CONTENT (VOLUMETRIC FRACTION)
1432 ! NOTE: FROZEN SOIL MOISTURE = SMC - SH2O
1433 ! SNOWH ACTUAL SNOW DEPTH (M)
1434 ! SNEQV LIQUID WATER-EQUIVALENT SNOW DEPTH (M)
1435 ! NOTE: SNOW DENSITY = SNEQV/SNOWH
1436 ! ALBEDO SURFACE ALBEDO INCLUDING SNOW EFFECT (UNITLESS FRACTION)
1437 ! =SNOW-FREE ALBEDO (ALB) WHEN SNEQV=0, OR
1438 ! =FCT(MSNOALB,ALB,VEGTYP,SHDFAC,SHDMIN) WHEN SNEQV>0
1439 ! CH SURFACE EXCHANGE COEFFICIENT FOR HEAT AND MOISTURE
1440 ! (M S-1); NOTE: CH IS TECHNICALLY A CONDUCTANCE SINCE
1441 ! IT HAS BEEN MULTIPLIED BY WIND SPEED.
1442 ! CM SURFACE EXCHANGE COEFFICIENT FOR MOMENTUM (M S-1); NOTE:
1443 ! CM IS TECHNICALLY A CONDUCTANCE SINCE IT HAS BEEN
1444 ! MULTIPLIED BY WIND SPEED.
1445 ! ----------------------------------------------------------------------
1446 ! 7. OUTPUT (O):
1447 ! ----------------------------------------------------------------------
1448 ! OUTPUT VARIABLES NECESSARY FOR A COUPLED NUMERICAL WEATHER PREDICTION
1449 ! MODEL, E.G. NOAA/NWS/NCEP MESOSCALE ETA MODEL. FOR THIS APPLICATION,
1450 ! THE REMAINING OUTPUT/DIAGNOSTIC/PARAMETER BLOCKS BELOW ARE NOT
1451 ! NECESSARY. OTHER APPLICATIONS MAY REQUIRE DIFFERENT OUTPUT VARIABLES.
1452 ! ETA ACTUAL LATENT HEAT FLUX (W m-2: NEGATIVE, IF UP FROM
1453 ! SURFACE)
1454 ! ETA_KINEMATIC atctual latent heat flux in Kg m-2 s-1
1455 ! SHEAT SENSIBLE HEAT FLUX (W M-2: NEGATIVE, IF UPWARD FROM
1456 ! SURFACE)
1457 ! FDOWN Radiation forcing at the surface (W m-2) = SOLDN*(1-alb)+LWDN
1458 ! ----------------------------------------------------------------------
1459 ! EC CANOPY WATER EVAPORATION (W m-2)
1460 ! EDIR DIRECT SOIL EVAPORATION (W m-2)
1461 ! ET(NSOIL) PLANT TRANSPIRATION FROM A PARTICULAR ROOT (SOIL) LAYER
1462 ! (W m-2)
1463 ! ETT TOTAL PLANT TRANSPIRATION (W m-2)
1464 ! ESNOW SUBLIMATION FROM (OR DEPOSITION TO IF <0) SNOWPACK
1465 ! (W m-2)
1466 ! DRIP THROUGH-FALL OF PRECIP AND/OR DEW IN EXCESS OF CANOPY
1467 ! WATER-HOLDING CAPACITY (M)
1468 ! DEW DEWFALL (OR FROSTFALL FOR T<273.15) (M)
1469 ! ----------------------------------------------------------------------
1470 ! BETA RATIO OF ACTUAL/POTENTIAL EVAP (DIMENSIONLESS)
1471 ! ETP POTENTIAL EVAPORATION (W m-2)
1472 ! SSOIL SOIL HEAT FLUX (W M-2: NEGATIVE IF DOWNWARD FROM SURFACE)
1473 ! ----------------------------------------------------------------------
1474 ! FLX1 PRECIP-SNOW SFC (W M-2)
1475 ! FLX2 FREEZING RAIN LATENT HEAT FLUX (W M-2)
1476 ! FLX3 PHASE-CHANGE HEAT FLUX FROM SNOWMELT (W M-2)
1477 ! ----------------------------------------------------------------------
1478 ! SNOMLT SNOW MELT (M) (WATER EQUIVALENT)
1479 ! SNCOVR FRACTIONAL SNOW COVER (UNITLESS FRACTION, 0-1)
1480 ! ----------------------------------------------------------------------
1481 ! RUNOFF1 SURFACE RUNOFF (M S-1), NOT INFILTRATING THE SURFACE
1482 ! RUNOFF2 SUBSURFACE RUNOFF (M S-1), DRAINAGE OUT BOTTOM OF LAST
1483 ! SOIL LAYER (BASEFLOW)
1484 ! RUNOFF3 NUMERICAL TRUNCTATION IN EXCESS OF POROSITY (SMCMAX)
1485 ! FOR A GIVEN SOIL LAYER AT THE END OF A TIME STEP (M S-1).
1486 ! Note: the above RUNOFF2 is actually the sum of RUNOFF2 and RUNOFF3
1487 ! ----------------------------------------------------------------------
1488 ! RC CANOPY RESISTANCE (S M-1)
1489 ! PC PLANT COEFFICIENT (UNITLESS FRACTION, 0-1) WHERE PC*ETP
1490 ! = ACTUAL TRANSP
1491 ! XLAI LEAF AREA INDEX (DIMENSIONLESS)
1492 ! RSMIN MINIMUM CANOPY RESISTANCE (S M-1)
1493 ! RCS INCOMING SOLAR RC FACTOR (DIMENSIONLESS)
1494 ! RCT AIR TEMPERATURE RC FACTOR (DIMENSIONLESS)
1495 ! RCQ ATMOS VAPOR PRESSURE DEFICIT RC FACTOR (DIMENSIONLESS)
1496 ! RCSOIL SOIL MOISTURE RC FACTOR (DIMENSIONLESS)
1497 ! ----------------------------------------------------------------------
1498 ! 8. DIAGNOSTIC OUTPUT (D):
1499 ! ----------------------------------------------------------------------
1500 ! SOILW AVAILABLE SOIL MOISTURE IN ROOT ZONE (UNITLESS FRACTION
1501 ! BETWEEN SMCWLT AND SMCMAX)
1502 ! SOILM TOTAL SOIL COLUMN MOISTURE CONTENT (FROZEN+UNFROZEN) (M)
1503 ! Q1 Effective mixing ratio at surface (kg kg-1), used for
1504 ! diagnosing the mixing ratio at 2 meter for coupled model
1505 ! ----------------------------------------------------------------------
1506 ! 9. PARAMETERS (P):
1507 ! ----------------------------------------------------------------------
1508 ! SMCWLT WILTING POINT (VOLUMETRIC)
1509 ! SMCDRY DRY SOIL MOISTURE THRESHOLD WHERE DIRECT EVAP FRM TOP
1510 ! LAYER ENDS (VOLUMETRIC)
1511 ! SMCREF SOIL MOISTURE THRESHOLD WHERE TRANSPIRATION BEGINS TO
1512 ! STRESS (VOLUMETRIC)
1513 ! SMCMAX POROSITY, I.E. SATURATED VALUE OF SOIL MOISTURE
1514 ! (VOLUMETRIC)
1515 ! NROOT NUMBER OF ROOT LAYERS, A FUNCTION OF VEG TYPE, DETERMINED
1516 ! IN SUBROUTINE REDPRM.
1517 ! ----------------------------------------------------------------------
1518
1519
1520 IMPLICIT NONE
1521 ! ----------------------------------------------------------------------
1522
1523 ! DECLARATIONS - LOGICAL AND CHARACTERS
1524 ! ----------------------------------------------------------------------
1525 LOGICAL, INTENT(IN):: LOCAL
1526 LOGICAL :: FRZGRA, SNOWNG
1527 CHARACTER (LEN=4), INTENT(IN):: LLANDUSE, LSOIL
1528
1529 ! ----------------------------------------------------------------------
1530 ! DECLARATIONS - INTEGER
1531 ! ----------------------------------------------------------------------
1532 INTEGER,INTENT(IN) :: ICE,NSOIL,SLOPETYP,SOILTYP,VEGTYP
1533 INTEGER,INTENT(OUT):: NROOT
1534 INTEGER KZ, K, iout
1535
1536 ! ----------------------------------------------------------------------
1537 ! DECLARATIONS - REAL
1538 ! ----------------------------------------------------------------------
1539 REAL, INTENT(IN) :: SHDMIN,SHDMAX,DT,DQSDT2,LWDN,PRCP,PRCPRAIN, &
1540 Q2,Q2SAT,SFCPRS,SFCSPD,SFCTMP, SNOALB, &
1541 SOLDN,TBOT,TH2,ZLVL, &
1542 EMISSI
1543 REAL, INTENT(INOUT):: COSZ, SOLARDIRECT,ALBEDO,CH,CM, &
1544 CMC,SNEQV,SNCOVR,SNOWH,T1,XLAI,SHDFAC,Z0BRD,ALB
1545 REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SLDPTH
1546 REAL, DIMENSION(1:NSOIL), INTENT(OUT):: ET
1547 REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: SH2O, SMC, STC
1548 REAL,DIMENSION(1:NSOIL):: RTDIS, ZSOIL
1549
1550 REAL,INTENT(OUT) :: ETA_KINEMATIC,BETA,DEW,DRIP,EC,EDIR,ESNOW,ETA, &
1551 ETP,FLX1,FLX2,FLX3,SHEAT,PC,RUNOFF1,RUNOFF2, &
1552 RUNOFF3,RC,RSMIN,RCQ,RCS,RCSOIL,RCT,SSOIL, &
1553 SMCDRY,SMCMAX,SMCREF,SMCWLT,SNOMLT, SOILM, &
1554 SOILW,FDOWN,Q1
1555 REAL :: BEXP,CFACTR,CMCMAX,CSOIL,CZIL,DF1,DF1H,DF1A,DKSAT,DWSAT, &
1556 DSOIL,DTOT,ETT,FRCSNO,FRCSOI,EPSCA,F1,FXEXP,FRZX,HS, &
1557 KDT,LVH2O,PRCP1,PSISAT,QUARTZ,R,RCH,REFKDT,RR,RGL, &
1558 RSMAX, &
1559 RSNOW,SNDENS,SNCOND,SBETA,SN_NEW,SLOPE,SNUP,SALP,SOILWM, &
1560 SOILWW,T1V,T24,T2V,TH2V,TOPT,TFREEZ,TSNOW,ZBOT,Z0,PRCPF, &
1561 ETNS,PTU,LSUBS
1562
1563 ! ----------------------------------------------------------------------
1564 ! DECLARATIONS - PARAMETERS
1565 ! ----------------------------------------------------------------------
1566 PARAMETER (TFREEZ = 273.15)
1567 PARAMETER (LVH2O = 2.501E+6)
1568 PARAMETER (LSUBS = 2.83E+6)
1569 PARAMETER (R = 287.04)
1570 ! ----------------------------------------------------------------------
1571 ! INITIALIZATION
1572 ! ----------------------------------------------------------------------
1573 RUNOFF1 = 0.0
1574 RUNOFF2 = 0.0
1575 RUNOFF3 = 0.0
1576 SNOMLT = 0.0
1577
1578 ! ----------------------------------------------------------------------
1579 ! THE VARIABLE "ICE" IS A FLAG DENOTING SEA-ICE CASE
1580 ! ----------------------------------------------------------------------
1581 ! SEA-ICE LAYERS ARE EQUAL THICKNESS AND SUM TO 3 METERS
1582 ! ----------------------------------------------------------------------
1583 IF (ICE == 1) THEN
1584 DO KZ = 1,NSOIL
1585 ZSOIL (KZ) = -3.* FLOAT (KZ)/ FLOAT (NSOIL)
1586 END DO
1587
1588 ! ----------------------------------------------------------------------
1589 ! CALCULATE DEPTH (NEGATIVE) BELOW GROUND FROM TOP SKIN SFC TO BOTTOM OF
1590 ! EACH SOIL LAYER. NOTE: SIGN OF ZSOIL IS NEGATIVE (DENOTING BELOW
1591 ! GROUND)
1592 ! ----------------------------------------------------------------------
1593 ELSE
1594 ZSOIL (1) = - SLDPTH (1)
1595 DO KZ = 2,NSOIL
1596 ZSOIL (KZ) = - SLDPTH (KZ) + ZSOIL (KZ -1)
1597 END DO
1598 ! ----------------------------------------------------------------------
1599 ! NEXT IS CRUCIAL CALL TO SET THE LAND-SURFACE PARAMETERS, INCLUDING
1600 ! SOIL-TYPE AND VEG-TYPE DEPENDENT PARAMETERS.
1601 ! ----------------------------------------------------------------------
1602 CALL REDPRM (VEGTYP,SOILTYP,SLOPETYP,CFACTR,CMCMAX,RSMAX,TOPT, &
1603 REFKDT,KDT,SBETA, SHDFAC,RSMIN,RGL,HS,ZBOT,FRZX, &
1604 PSISAT,SLOPE,SNUP,SALP,BEXP,DKSAT,DWSAT, &
1605 SMCMAX,SMCWLT,SMCREF,SMCDRY,F1,QUARTZ,FXEXP, &
1606 RTDIS,SLDPTH,ZSOIL,NROOT,NSOIL,Z0BRD,CZIL,XLAI, &
1607 CSOIL,ALB,PTU,LLANDUSE,LSOIL,LOCAL)
1608 END IF
1609
1610 !urban change
1611 IF(VEGTYP==1)THEN
1612 SHDFAC=0.05
1613 RSMIN=400.0
1614 SMCMAX = 0.45
1615 SMCREF = 0.42
1616 SMCWLT = 0.40
1617 ENDIF
1618
1619 ! ----------------------------------------------------------------------
1620 ! INITIALIZE PRECIPITATION LOGICALS.
1621 ! ----------------------------------------------------------------------
1622 SNOWNG = .FALSE.
1623 FRZGRA = .FALSE.
1624
1625 ! ----------------------------------------------------------------------
1626 ! IF SEA-ICE CASE, ASSIGN DEFAULT WATER-EQUIV SNOW ON TOP
1627 ! ----------------------------------------------------------------------
1628 IF (ICE == 1) THEN
1629 SNEQV = 0.01
1630 SNOWH = 0.05
1631 ! ----------------------------------------------------------------------
1632 ! IF INPUT SNOWPACK IS NONZERO, THEN COMPUTE SNOW DENSITY "SNDENS" AND
1633 ! SNOW THERMAL CONDUCTIVITY "SNCOND" (NOTE THAT CSNOW IS A FUNCTION
1634 ! SUBROUTINE)
1635 ! ----------------------------------------------------------------------
1636 END IF
1637 IF (SNEQV == 0.0) THEN
1638 SNDENS = 0.0
1639 SNOWH = 0.0
1640 SNCOND = 1.0
1641 ELSE
1642 SNDENS = SNEQV / SNOWH
1643 IF(SNDENS > 1.0) THEN
1644 CALL wrf_error_fatal ( 'Physical snow depth is less than snow water equiv.' )
1645 ENDIF
1646 CALL CSNOW (SNCOND,SNDENS)
1647 END IF
1648 ! ----------------------------------------------------------------------
1649 ! DETERMINE IF IT'S PRECIPITATING AND WHAT KIND OF PRECIP IT IS.
1650 ! IF IT'S PRCPING AND THE AIR TEMP IS COLDER THAN 0 C, IT'S SNOWING!
1651 ! IF IT'S PRCPING AND THE AIR TEMP IS WARMER THAN 0 C, BUT THE GRND
1652 ! TEMP IS COLDER THAN 0 C, FREEZING RAIN IS PRESUMED TO BE FALLING.
1653 ! ----------------------------------------------------------------------
1654 IF (PRCP > 0.0) THEN
1655 IF (SFCTMP <= TFREEZ) THEN
1656 SNOWNG = .TRUE.
1657 ELSE
1658 IF (T1 <= TFREEZ) FRZGRA = .TRUE.
1659 END IF
1660 END IF
1661 ! ----------------------------------------------------------------------
1662 ! IF EITHER PRCP FLAG IS SET, DETERMINE NEW SNOWFALL (CONVERTING PRCP
1663 ! RATE FROM KG M-2 S-1 TO A LIQUID EQUIV SNOW DEPTH IN METERS) AND ADD
1664 ! IT TO THE EXISTING SNOWPACK.
1665 ! NOTE THAT SINCE ALL PRECIP IS ADDED TO SNOWPACK, NO PRECIP INFILTRATES
1666 ! INTO THE SOIL SO THAT PRCP1 IS SET TO ZERO.
1667 ! ----------------------------------------------------------------------
1668 IF ( (SNOWNG) .OR. (FRZGRA) ) THEN
1669 SN_NEW = PRCP * DT * 0.001
1670 ! PRCP1 = 0.0
1671 ! change name of PRCP1 to PRCPF
1672 SNEQV = SNEQV + SN_NEW
1673 PRCPF = 0.0
1674
1675 ! ----------------------------------------------------------------------
1676 ! UPDATE SNOW DENSITY BASED ON NEW SNOWFALL, USING OLD AND NEW SNOW.
1677 ! UPDATE SNOW THERMAL CONDUCTIVITY
1678 ! ----------------------------------------------------------------------
1679 CALL SNOW_NEW (SFCTMP,SN_NEW,SNOWH,SNDENS)
1680 CALL CSNOW (SNCOND,SNDENS)
1681
1682 ! ----------------------------------------------------------------------
1683 ! PRECIP IS LIQUID (RAIN), HENCE SAVE IN THE PRECIP VARIABLE THAT
1684 ! LATER CAN WHOLELY OR PARTIALLY INFILTRATE THE SOIL (ALONG WITH
1685 ! ANY CANOPY "DRIP" ADDED TO THIS LATER)
1686 ! ----------------------------------------------------------------------
1687 ! PRCP1 = PRCP
1688 ELSE
1689 ! change name of PRCP1 to PRCPF
1690
1691 PRCPF = PRCP
1692 END IF
1693 ! ----------------------------------------------------------------------
1694 ! DETERMINE SNOWCOVER AND ALBEDO OVER LAND.
1695 ! ----------------------------------------------------------------------
1696 ! ----------------------------------------------------------------------
1697 ! IF SNOW DEPTH=0, SET SNOW FRACTION=0, ALBEDO=SNOW FREE ALBEDO.
1698 ! ----------------------------------------------------------------------
1699 IF (ICE == 0) THEN
1700 IF (SNEQV == 0.0) THEN
1701 SNCOVR = 0.0
1702 ALBEDO = ALB
1703 ELSE
1704 ! ----------------------------------------------------------------------
1705 ! DETERMINE SNOW FRACTIONAL COVERAGE.
1706 ! DETERMINE SURFACE ALBEDO MODIFICATION DUE TO SNOWDEPTH STATE.
1707 ! ----------------------------------------------------------------------
1708 CALL SNFRAC (SNEQV,SNUP,SALP,SNOWH,SNCOVR)
1709 CALL ALCALC (ALB,SNOALB,SHDFAC,SHDMIN,SNCOVR,TSNOW,ALBEDO)
1710 END IF
1711 ! ----------------------------------------------------------------------
1712 ! SNOW COVER, ALBEDO OVER SEA-ICE
1713 ! ----------------------------------------------------------------------
1714 ELSE
1715 SNCOVR = 1.0
1716 ALBEDO = 0.60
1717 END IF
1718 ! ----------------------------------------------------------------------
1719 ! THERMAL CONDUCTIVITY FOR SEA-ICE CASE
1720 ! ----------------------------------------------------------------------
1721 IF (ICE == 1) THEN
1722 DF1 = 2.2
1723 ELSE
1724 ! ----------------------------------------------------------------------
1725 ! NEXT CALCULATE THE SUBSURFACE HEAT FLUX, WHICH FIRST REQUIRES
1726 ! CALCULATION OF THE THERMAL DIFFUSIVITY. TREATMENT OF THE
1727 ! LATTER FOLLOWS THAT ON PAGES 148-149 FROM "HEAT TRANSFER IN
1728 ! COLD CLIMATES", BY V. J. LUNARDINI (PUBLISHED IN 1981
1729 ! BY VAN NOSTRAND REINHOLD CO.) I.E. TREATMENT OF TWO CONTIGUOUS
1730 ! "PLANE PARALLEL" MEDIUMS (NAMELY HERE THE FIRST SOIL LAYER
1731 ! AND THE SNOWPACK LAYER, IF ANY). THIS DIFFUSIVITY TREATMENT
1732 ! BEHAVES WELL FOR BOTH ZERO AND NONZERO SNOWPACK, INCLUDING THE
1733 ! LIMIT OF VERY THIN SNOWPACK. THIS TREATMENT ALSO ELIMINATES
1734 ! THE NEED TO IMPOSE AN ARBITRARY UPPER BOUND ON SUBSURFACE
1735 ! HEAT FLUX WHEN THE SNOWPACK BECOMES EXTREMELY THIN.
1736 ! ----------------------------------------------------------------------
1737 ! FIRST CALCULATE THERMAL DIFFUSIVITY OF TOP SOIL LAYER, USING
1738 ! BOTH THE FROZEN AND LIQUID SOIL MOISTURE, FOLLOWING THE
1739 ! SOIL THERMAL DIFFUSIVITY FUNCTION OF PETERS-LIDARD ET AL.
1740 ! (1998,JAS, VOL 55, 1209-1224), WHICH REQUIRES THE SPECIFYING
1741 ! THE QUARTZ CONTENT OF THE GIVEN SOIL CLASS (SEE ROUTINE REDPRM)
1742 ! ----------------------------------------------------------------------
1743 ! ----------------------------------------------------------------------
1744 ! NEXT ADD SUBSURFACE HEAT FLUX REDUCTION EFFECT FROM THE
1745 ! OVERLYING GREEN CANOPY, ADAPTED FROM SECTION 2.1.2 OF
1746 ! PETERS-LIDARD ET AL. (1997, JGR, VOL 102(D4))
1747 ! ----------------------------------------------------------------------
1748 CALL TDFCND (DF1,SMC (1),QUARTZ,SMCMAX,SH2O (1))
1749
1750 !urban change
1751 IF (VEGTYP==1) DF1=3.24
1752
1753 DF1 = DF1 * EXP (SBETA * SHDFAC)
1754 ! ----------------------------------------------------------------------
1755 ! FINALLY "PLANE PARALLEL" SNOWPACK EFFECT FOLLOWING
1756 ! V.J. LINARDINI REFERENCE CITED ABOVE. NOTE THAT DTOT IS
1757 ! COMBINED DEPTH OF SNOWDEPTH AND THICKNESS OF FIRST SOIL LAYER
1758 ! ----------------------------------------------------------------------
1759 END IF
1760
1761 DSOIL = - (0.5 * ZSOIL (1))
1762 IF (SNEQV == 0.) THEN
1763 SSOIL = DF1 * (T1- STC (1) ) / DSOIL
1764 ELSE
1765 DTOT = SNOWH + DSOIL
1766 FRCSNO = SNOWH / DTOT
1767
1768 ! 1. HARMONIC MEAN (SERIES FLOW)
1769 ! DF1 = (SNCOND*DF1)/(FRCSOI*SNCOND+FRCSNO*DF1)
1770 FRCSOI = DSOIL / DTOT
1771 ! 2. ARITHMETIC MEAN (PARALLEL FLOW)
1772 ! DF1 = FRCSNO*SNCOND + FRCSOI*DF1
1773 DF1H = (SNCOND * DF1)/ (FRCSOI * SNCOND+ FRCSNO * DF1)
1774
1775 ! 3. GEOMETRIC MEAN (INTERMEDIATE BETWEEN HARMONIC AND ARITHMETIC MEAN)
1776 ! DF1 = (SNCOND**FRCSNO)*(DF1**FRCSOI)
1777 ! TEST - MBEK, 10 Jan 2002
1778 ! weigh DF by snow fraction
1779 ! DF1 = DF1H*SNCOVR + DF1A*(1.0-SNCOVR)
1780 ! DF1 = DF1H*SNCOVR + DF1*(1.0-SNCOVR)
1781 DF1A = FRCSNO * SNCOND+ FRCSOI * DF1
1782
1783 ! ----------------------------------------------------------------------
1784 ! CALCULATE SUBSURFACE HEAT FLUX, SSOIL, FROM FINAL THERMAL DIFFUSIVITY
1785 ! OF SURFACE MEDIUMS, DF1 ABOVE, AND SKIN TEMPERATURE AND TOP
1786 ! MID-LAYER SOIL TEMPERATURE
1787 ! ----------------------------------------------------------------------
1788 DF1 = DF1A * SNCOVR + DF1* (1.0- SNCOVR)
1789 SSOIL = DF1 * (T1- STC (1) ) / DTOT
1790 END IF
1791 ! ----------------------------------------------------------------------
1792 ! DETERMINE SURFACE ROUGHNESS OVER SNOWPACK USING SNOW CONDITION FROM
1793 ! THE PREVIOUS TIMESTEP.
1794 ! ----------------------------------------------------------------------
1795 IF (SNCOVR > 0. ) THEN
1796 CALL SNOWZ0 (SNCOVR,Z0,Z0BRD)
1797 ELSE
1798 Z0=Z0BRD
1799 END IF
1800 ! ----------------------------------------------------------------------
1801 ! NEXT CALL ROUTINE SFCDIF TO CALCULATE THE SFC EXCHANGE COEF (CH) FOR
1802 ! HEAT AND MOISTURE.
1803
1804 ! NOTE !!!
1805 ! DO NOT CALL SFCDIF UNTIL AFTER ABOVE CALL TO REDPRM, IN CASE
1806 ! ALTERNATIVE VALUES OF ROUGHNESS LENGTH (Z0) AND ZILINTINKEVICH COEF
1807 ! (CZIL) ARE SET THERE VIA NAMELIST I/O.
1808
1809 ! NOTE !!!
1810 ! ROUTINE SFCDIF RETURNS A CH THAT REPRESENTS THE WIND SPD TIMES THE
1811 ! "ORIGINAL" NONDIMENSIONAL "Ch" TYPICAL IN LITERATURE. HENCE THE CH
1812 ! RETURNED FROM SFCDIF HAS UNITS OF M/S. THE IMPORTANT COMPANION
1813 ! COEFFICIENT OF CH, CARRIED HERE AS "RCH", IS THE CH FROM SFCDIF TIMES
1814 ! AIR DENSITY AND PARAMETER "CP". "RCH" IS COMPUTED IN "CALL PENMAN".
1815 ! RCH RATHER THAN CH IS THE COEFF USUALLY INVOKED LATER IN EQNS.
1816
1817 ! NOTE !!!
1818 ! ----------------------------------------------------------------------
1819 ! SFCDIF ALSO RETURNS THE SURFACE EXCHANGE COEFFICIENT FOR MOMENTUM, CM,
1820 ! ALSO KNOWN AS THE SURFACE DRAGE COEFFICIENT. Needed as a state variable
1821 ! for iterative/implicit solution of CH in SFCDIF
1822 ! ----------------------------------------------------------------------
1823 ! IF(.NOT.LCH) THEN
1824 ! T1V = T1 * (1.0+ 0.61 * Q2)
1825 ! TH2V = TH2 * (1.0+ 0.61 * Q2)
1826 ! CALL SFCDIF_off (ZLVL,Z0,T1V,TH2V,SFCSPD,CZIL,CM,CH)
1827 ! ENDIF
1828
1829 ! ----------------------------------------------------------------------
1830 ! CALL PENMAN SUBROUTINE TO CALCULATE POTENTIAL EVAPORATION (ETP), AND
1831 ! OTHER PARTIAL PRODUCTS AND SUMS SAVE IN COMMON/RITE FOR LATER
1832 ! CALCULATIONS.
1833 ! ----------------------------------------------------------------------
1834 ! ----------------------------------------------------------------------
1835 ! CALCULATE TOTAL DOWNWARD RADIATION (SOLAR PLUS LONGWAVE) NEEDED IN
1836 ! PENMAN EP SUBROUTINE THAT FOLLOWS
1837 ! ----------------------------------------------------------------------
1838 FDOWN = SOLDN * (1.0- ALBEDO) + LWDN
1839 ! ----------------------------------------------------------------------
1840 ! CALC VIRTUAL TEMPS AND VIRTUAL POTENTIAL TEMPS NEEDED BY SUBROUTINES
1841 ! PENMAN.
1842 T2V = SFCTMP * (1.0+ 0.61 * Q2 )
1843
1844 iout=0
1845 if(iout.eq.1) then
1846 print*,'before penman'
1847 print*,' SFCTMP',SFCTMP,'SFCPRS',SFCPRS,'CH',CH,'T2V',T2V, &
1848 'TH2',TH2,'PRCP',PRCP,'FDOWN',FDOWN,'T24',T24,'SSOIL',SSOIL, &
1849 'Q2',Q2,'Q2SAT',Q2SAT,'ETP',ETP,'RCH',RCH, &
1850 'EPSCA',EPSCA,'RR',RR ,'SNOWNG',SNOWNG,'FRZGRA',FRZGRA, &
1851 'DQSDT2',DQSDT2,'FLX2',FLX2,'SNOWH',SNOWH,'SNEQV',SNEQV, &
1852 ' DSOIL',DSOIL,' FRCSNO',FRCSNO,' SNCOVR',SNCOVR,' DTOT',DTOT, &
1853 ' ZSOIL (1)',ZSOIL(1),' DF1',DF1,'T1',T1,' STC1',STC(1), &
1854 'ALBEDO',ALBEDO,'SMC',SMC,'STC',STC,'SH2O',SH2O
1855 endif
1856
1857 CALL PENMAN (SFCTMP,SFCPRS,CH,T2V,TH2,PRCP,FDOWN,T24,SSOIL, &
1858 Q2,Q2SAT,ETP,RCH,EPSCA,RR,SNOWNG,FRZGRA, &
1859 DQSDT2,FLX2,EMISSI)
1860
1861 ! ----------------------------------------------------------------------
1862 ! CALL CANRES TO CALCULATE THE CANOPY RESISTANCE AND CONVERT IT INTO PC
1863 ! IF NONZERO GREENNESS FRACTION
1864 ! ----------------------------------------------------------------------
1865
1866 ! print*,'after penman ETP',ETP
1867 ! ----------------------------------------------------------------------
1868 ! FROZEN GROUND EXTENSION: TOTAL SOIL WATER "SMC" WAS REPLACED
1869 ! BY UNFROZEN SOIL WATER "SH2O" IN CALL TO CANRES BELOW
1870 ! ----------------------------------------------------------------------
1871 IF (SHDFAC > 0.) THEN
1872 CALL CANRES (SOLDN,CH,SFCTMP,Q2,SFCPRS,SH2O,ZSOIL,NSOIL, &
1873 SMCWLT,SMCREF,RSMIN,RC,PC,NROOT,Q2SAT,DQSDT2, &
1874 TOPT,RSMAX,RGL,HS,XLAI, &
1875 RCS,RCT,RCQ,RCSOIL,EMISSI)
1876 END IF
1877 ! ----------------------------------------------------------------------
1878 ! NOW DECIDE MAJOR PATHWAY BRANCH TO TAKE DEPENDING ON WHETHER SNOWPACK
1879 ! EXISTS OR NOT:
1880 ! ----------------------------------------------------------------------
1881 ESNOW = 0.0
1882 IF (SNEQV == 0.0) THEN
1883 CALL NOPAC (ETP,ETA,PRCP,SMC,SMCMAX,SMCWLT, &
1884 SMCREF,SMCDRY,CMC,CMCMAX,NSOIL,DT, &
1885 SHDFAC, &
1886 SBETA,Q2,T1,SFCTMP,T24,TH2,FDOWN,F1,EMISSI, &
1887 SSOIL, &
1888 STC,EPSCA,BEXP,PC,RCH,RR,CFACTR, &
1889 SH2O,SLOPE,KDT,FRZX,PSISAT,ZSOIL, &
1890 DKSAT,DWSAT,TBOT,ZBOT,RUNOFF1,RUNOFF2, &
1891 RUNOFF3,EDIR,EC,ET,ETT,NROOT,ICE,RTDIS, &
1892 QUARTZ,FXEXP,CSOIL, &
1893 BETA,DRIP,DEW,FLX1,FLX3,VEGTYP)
1894 ETA_KINEMATIC = ETA
1895 ETA = ETA * LVH2O
1896 ETP = ETP*LVH2O
1897 ! BETA = ETA/ETP
1898 ELSE
1899 CALL SNOPAC (ETP,ETA,PRCP,PRCPF,SNOWNG,SMC,SMCMAX,SMCWLT, &
1900 SMCREF,SMCDRY,CMC,CMCMAX,NSOIL,DT, &
1901 SBETA,DF1, &
1902 Q2,T1,SFCTMP,T24,TH2,FDOWN,F1,SSOIL,STC,EPSCA, &
1903 SFCPRS,BEXP,PC,RCH,RR,CFACTR,SNCOVR,SNEQV,SNDENS,&
1904 SNOWH,SH2O,SLOPE,KDT,FRZX,PSISAT, &
1905 ZSOIL,DWSAT,DKSAT,TBOT,ZBOT,SHDFAC,RUNOFF1, &
1906 RUNOFF2,RUNOFF3,EDIR,EC,ET,ETT,NROOT,SNOMLT, &
1907 ICE,RTDIS,QUARTZ,FXEXP,CSOIL, &
1908 BETA,DRIP,DEW,FLX1,FLX2,FLX3,ESNOW,ETNS,EMISSI, &
1909 VEGTYP)
1910 ETA_KINEMATIC = ESNOW + ETNS
1911 ESNOW = ESNOW * LSUBS
1912 ! ETA = ESNOW * LSUBS + ETNS * LVH2O
1913 ETA = ESNOW + ETNS * LVH2O
1914 ETP = ETP*LSUBS
1915 ! ETP = ETP* ( (SNCOVR*LSUBS+(1.-SNCOVR)*LVH2O )
1916 ! BETA = ETA/ETP
1917 END IF
1918
1919 ! Calculate effective mixing ratio at grnd level (skin)
1920 !
1921 ! Q1=Q2+ETA*CP/RCH
1922 Q1=Q2+ETA_KINEMATIC*CP/RCH
1923 !
1924 ! ----------------------------------------------------------------------
1925 ! PREPARE SENSIBLE HEAT (H) FOR RETURN TO PARENT MODEL
1926 ! ----------------------------------------------------------------------
1927 ! ----------------------------------------------------------------------
1928 ! CONVERT UNITS AND/OR SIGN OF TOTAL EVAP (ETA), POTENTIAL EVAP (ETP),
1929 ! SUBSURFACE HEAT FLUX (S), AND RUNOFFS FOR WHAT PARENT MODEL EXPECTS
1930 ! CONVERT ETA FROM KG M-2 S-1 TO W M-2
1931 ! ----------------------------------------------------------------------
1932 !ek ETA = ETA*LVH2O
1933 SHEAT = - (CH * CP * SFCPRS)/ (R * T2V) * ( TH2- T1 )
1934
1935 ! ----------------------------------------------------------------------
1936 ! CONVERT OTHER EVAP TERMS FROM M S-1 TO KG M-2 S-1, THEN TO W M-2
1937 ! ----------------------------------------------------------------------
1938 ! EDIR = EDIR*1000.
1939 ! EC = EC*1000.
1940 ! DO K = 1,NSOIL
1941 ! ET(K) = ET(K)*1000.
1942 ! END DO
1943 ! ETT = ETT*1000.
1944
1945 ! ETP = ETP * LVH2O
1946 EC = EC * LVH2O
1947 EDIR = EDIR * LVH2O
1948 DO K = 1,NSOIL
1949 ET (K) = ET (K)* LVH2O
1950 END DO
1951 ETT = ETT * LVH2O
1952
1953 ! ----------------------------------------------------------------------
1954 ! CONVERT THE SIGN OF SOIL HEAT FLUX SO THAT:
1955 ! SSOIL>0: WARM THE SURFACE (NIGHT TIME)
1956 ! SSOIL<0: COOL THE SURFACE (DAY TIME)
1957 ! ----------------------------------------------------------------------
1958 ! ESNOW = ESNOW * LVH2O
1959
1960 ! ----------------------------------------------------------------------
1961 ! CONVERT RUNOFF3 (INTERNAL LAYER RUNOFF FROM SUPERSAT) FROM M TO M S-1
1962 ! AND ADD TO SUBSURFACE RUNOFF/DRAINAGE/BASEFLOW
1963 ! ----------------------------------------------------------------------
1964 SSOIL = -1.0* SSOIL
1965 RUNOFF3 = RUNOFF3/ DT
1966
1967 ! ----------------------------------------------------------------------
1968 ! TOTAL COLUMN SOIL MOISTURE IN METERS (SOILM) AND ROOT-ZONE
1969 ! SOIL MOISTURE AVAILABILITY (FRACTION) RELATIVE TO POROSITY/SATURATION
1970 ! ----------------------------------------------------------------------
1971 RUNOFF2 = RUNOFF2+ RUNOFF3
1972 IF (ICE == 0) THEN
1973 SOILM = -1.0* SMC (1)* ZSOIL (1)
1974 DO K = 2,NSOIL
1975 SOILM = SOILM + SMC (K)* (ZSOIL (K -1) - ZSOIL (K))
1976 END DO
1977 SOILWM = -1.0* (SMCMAX - SMCWLT)* ZSOIL (1)
1978 SOILWW = -1.0* (SMC (1) - SMCWLT)* ZSOIL (1)
1979
1980 IF (NROOT >= 2) THEN
1981 DO K = 2,NROOT
1982 SOILWM = SOILWM + (SMCMAX - SMCWLT)* (ZSOIL (K -1) - ZSOIL (K))
1983 SOILWW = SOILWW + (SMC(K) - SMCWLT)* (ZSOIL (K -1) - ZSOIL (K))
1984 END DO
1985 END IF
1986 IF (SOILWM .LT. 1.E-6) THEN
1987 SOILWM = 0.0
1988 SOILW = 0.0
1989 SOILM = 0.0
1990 ELSE
1991 SOILW = SOILWW / SOILWM
1992 END IF
1993 ELSE
1994 SOILWM = 0.0
1995 SOILW = 0.0
1996 SOILM = 0.0
1997 END IF
1998
1999 ! ----------------------------------------------------------------------
2000 END SUBROUTINE SFLX
2001 ! ----------------------------------------------------------------------
2002
2003 SUBROUTINE ALCALC (ALB,SNOALB,SHDFAC,SHDMIN,SNCOVR,TSNOW,ALBEDO)
2004
2005 ! ----------------------------------------------------------------------
2006 ! CALCULATE ALBEDO INCLUDING SNOW EFFECT (0 -> 1)
2007 ! ALB SNOWFREE ALBEDO
2008 ! SNOALB MAXIMUM (DEEP) SNOW ALBEDO
2009 ! SHDFAC AREAL FRACTIONAL COVERAGE OF GREEN VEGETATION
2010 ! SHDMIN MINIMUM AREAL FRACTIONAL COVERAGE OF GREEN VEGETATION
2011 ! SNCOVR FRACTIONAL SNOW COVER
2012 ! ALBEDO SURFACE ALBEDO INCLUDING SNOW EFFECT
2013 ! TSNOW SNOW SURFACE TEMPERATURE (K)
2014 ! ----------------------------------------------------------------------
2015 IMPLICIT NONE
2016
2017 ! ----------------------------------------------------------------------
2018 ! SNOALB IS ARGUMENT REPRESENTING MAXIMUM ALBEDO OVER DEEP SNOW,
2019 ! AS PASSED INTO SFLX, AND ADAPTED FROM THE SATELLITE-BASED MAXIMUM
2020 ! SNOW ALBEDO FIELDS PROVIDED BY D. ROBINSON AND G. KUKLA
2021 ! (1985, JCAM, VOL 24, 402-411)
2022 ! ----------------------------------------------------------------------
2023 REAL, INTENT(IN) :: ALB, SNOALB, SHDFAC, SHDMIN, SNCOVR, TSNOW
2024 REAL, INTENT(OUT) :: ALBEDO
2025 ALBEDO = ALB + (1.0- (SHDFAC - SHDMIN))* SNCOVR * (SNOALB - ALB)
2026
2027 ! BASE FORMULATION (DICKINSON ET AL., 1986, COGLEY ET AL., 1990)
2028 ! IF (TSNOW.LE.263.16) THEN
2029 ! ALBEDO=SNOALB
2030 ! ELSE
2031 ! IF (TSNOW.LT.273.16) THEN
2032 ! TM=0.1*(TSNOW-263.16)
2033 ! ALBEDO=0.5*((0.9-0.2*(TM**3))+(0.8-0.16*(TM**3)))
2034 ! ELSE
2035 ! ALBEDO=0.67
2036 ! ENDIF
2037 ! ENDIF
2038
2039 ! ISBA FORMULATION (VERSEGHY, 1991; BAKER ET AL., 1990)
2040 ! IF (TSNOW.LT.273.16) THEN
2041 ! ALBEDO=SNOALB-0.008*DT/86400
2042 ! ELSE
2043 ! ALBEDO=(SNOALB-0.5)*EXP(-0.24*DT/86400)+0.5
2044 ! ENDIF
2045
2046 IF (ALBEDO > SNOALB) ALBEDO = SNOALB
2047
2048 ! ----------------------------------------------------------------------
2049 END SUBROUTINE ALCALC
2050 ! ----------------------------------------------------------------------
2051
2052 SUBROUTINE CANRES (SOLAR,CH,SFCTMP,Q2,SFCPRS,SMC,ZSOIL,NSOIL, &
2053 SMCWLT,SMCREF,RSMIN,RC,PC,NROOT,Q2SAT,DQSDT2, &
2054 TOPT,RSMAX,RGL,HS,XLAI, &
2055 RCS,RCT,RCQ,RCSOIL,EMISSI)
2056
2057 ! ----------------------------------------------------------------------
2058 ! SUBROUTINE CANRES
2059 ! ----------------------------------------------------------------------
2060 ! CALCULATE CANOPY RESISTANCE WHICH DEPENDS ON INCOMING SOLAR RADIATION,
2061 ! AIR TEMPERATURE, ATMOSPHERIC WATER VAPOR PRESSURE DEFICIT AT THE
2062 ! LOWEST MODEL LEVEL, AND SOIL MOISTURE (PREFERABLY UNFROZEN SOIL
2063 ! MOISTURE RATHER THAN TOTAL)
2064 ! ----------------------------------------------------------------------
2065 ! SOURCE: JARVIS (1976), NOILHAN AND PLANTON (1989, MWR), JACQUEMIN AND
2066 ! NOILHAN (1990, BLM)
2067 ! SEE ALSO: CHEN ET AL (1996, JGR, VOL 101(D3), 7251-7268), EQNS 12-14
2068 ! AND TABLE 2 OF SEC. 3.1.2
2069 ! ----------------------------------------------------------------------
2070 ! INPUT:
2071 ! SOLAR INCOMING SOLAR RADIATION
2072 ! CH SURFACE EXCHANGE COEFFICIENT FOR HEAT AND MOISTURE
2073 ! SFCTMP AIR TEMPERATURE AT 1ST LEVEL ABOVE GROUND
2074 ! Q2 AIR HUMIDITY AT 1ST LEVEL ABOVE GROUND
2075 ! Q2SAT SATURATION AIR HUMIDITY AT 1ST LEVEL ABOVE GROUND
2076 ! DQSDT2 SLOPE OF SATURATION HUMIDITY FUNCTION WRT TEMP
2077 ! SFCPRS SURFACE PRESSURE
2078 ! SMC VOLUMETRIC SOIL MOISTURE
2079 ! ZSOIL SOIL DEPTH (NEGATIVE SIGN, AS IT IS BELOW GROUND)
2080 ! NSOIL NO. OF SOIL LAYERS
2081 ! NROOT NO. OF SOIL LAYERS IN ROOT ZONE (1.LE.NROOT.LE.NSOIL)
2082 ! XLAI LEAF AREA INDEX
2083 ! SMCWLT WILTING POINT
2084 ! SMCREF REFERENCE SOIL MOISTURE (WHERE SOIL WATER DEFICIT STRESS
2085 ! SETS IN)
2086 ! RSMIN, RSMAX, TOPT, RGL, HS ARE CANOPY STRESS PARAMETERS SET IN
2087 ! SURBOUTINE REDPRM
2088 ! OUTPUT:
2089 ! PC PLANT COEFFICIENT
2090 ! RC CANOPY RESISTANCE
2091 ! ----------------------------------------------------------------------
2092
2093 IMPLICIT NONE
2094 INTEGER, INTENT(IN) :: NROOT,NSOIL
2095 INTEGER K
2096 REAL, INTENT(IN) :: CH,DQSDT2,HS,Q2,Q2SAT,RSMIN,RGL,RSMAX, &
2097 SFCPRS,SFCTMP,SMCREF,SMCWLT, SOLAR,TOPT,XLAI, &
2098 EMISSI
2099 REAL,DIMENSION(1:NSOIL), INTENT(IN) :: SMC,ZSOIL
2100 REAL, INTENT(OUT):: PC,RC,RCQ,RCS,RCSOIL,RCT
2101 REAL :: DELTA,FF,GX,P,RR
2102 REAL, DIMENSION(1:NSOIL) :: PART
2103 REAL, PARAMETER :: SLV = 2.501000E6
2104
2105
2106 ! ----------------------------------------------------------------------
2107 ! INITIALIZE CANOPY RESISTANCE MULTIPLIER TERMS.
2108 ! ----------------------------------------------------------------------
2109 RCS = 0.0
2110 RCT = 0.0
2111 RCQ = 0.0
2112 RCSOIL = 0.0
2113
2114 ! ----------------------------------------------------------------------
2115 ! CONTRIBUTION DUE TO INCOMING SOLAR RADIATION
2116 ! ----------------------------------------------------------------------
2117 RC = 0.0
2118 FF = 0.55*2.0* SOLAR / (RGL * XLAI)
2119 RCS = (FF + RSMIN / RSMAX) / (1.0+ FF)
2120
2121 ! ----------------------------------------------------------------------
2122 ! CONTRIBUTION DUE TO AIR TEMPERATURE AT FIRST MODEL LEVEL ABOVE GROUND
2123 ! RCT EXPRESSION FROM NOILHAN AND PLANTON (1989, MWR).
2124 ! ----------------------------------------------------------------------
2125 RCS = MAX (RCS,0.0001)
2126 RCT = 1.0- 0.0016* ( (TOPT - SFCTMP)**2.0)
2127
2128 ! ----------------------------------------------------------------------
2129 ! CONTRIBUTION DUE TO VAPOR PRESSURE DEFICIT AT FIRST MODEL LEVEL.
2130 ! RCQ EXPRESSION FROM SSIB
2131 ! ----------------------------------------------------------------------
2132 RCT = MAX (RCT,0.0001)
2133 RCQ = 1.0/ (1.0+ HS * (Q2SAT - Q2))
2134
2135 ! ----------------------------------------------------------------------
2136 ! CONTRIBUTION DUE TO SOIL MOISTURE AVAILABILITY.
2137 ! DETERMINE CONTRIBUTION FROM EACH SOIL LAYER, THEN ADD THEM UP.
2138 ! ----------------------------------------------------------------------
2139 RCQ = MAX (RCQ,0.01)
2140 GX = (SMC (1) - SMCWLT) / (SMCREF - SMCWLT)
2141 IF (GX > 1.) GX = 1.
2142 IF (GX < 0.) GX = 0.
2143
2144 ! ----------------------------------------------------------------------
2145 ! USE SOIL DEPTH AS WEIGHTING FACTOR
2146 ! ----------------------------------------------------------------------
2147 ! ----------------------------------------------------------------------
2148 ! USE ROOT DISTRIBUTION AS WEIGHTING FACTOR
2149 ! PART(1) = RTDIS(1) * GX
2150 ! ----------------------------------------------------------------------
2151 PART (1) = (ZSOIL (1)/ ZSOIL (NROOT)) * GX
2152 DO K = 2,NROOT
2153 GX = (SMC (K) - SMCWLT) / (SMCREF - SMCWLT)
2154 IF (GX > 1.) GX = 1.
2155 IF (GX < 0.) GX = 0.
2156 ! ----------------------------------------------------------------------
2157 ! USE SOIL DEPTH AS WEIGHTING FACTOR
2158 ! ----------------------------------------------------------------------
2159 ! ----------------------------------------------------------------------
2160 ! USE ROOT DISTRIBUTION AS WEIGHTING FACTOR
2161 ! PART(K) = RTDIS(K) * GX
2162 ! ----------------------------------------------------------------------
2163 PART (K) = ( (ZSOIL (K) - ZSOIL (K -1))/ ZSOIL (NROOT)) * GX
2164 END DO
2165 DO K = 1,NROOT
2166 RCSOIL = RCSOIL + PART (K)
2167 END DO
2168
2169 ! ----------------------------------------------------------------------
2170 ! DETERMINE CANOPY RESISTANCE DUE TO ALL FACTORS. CONVERT CANOPY
2171 ! RESISTANCE (RC) TO PLANT COEFFICIENT (PC) TO BE USED WITH POTENTIAL
2172 ! EVAP IN DETERMINING ACTUAL EVAP. PC IS DETERMINED BY:
2173 ! PC * LINERIZED PENMAN POTENTIAL EVAP =
2174 ! PENMAN-MONTEITH ACTUAL EVAPORATION (CONTAINING RC TERM).
2175 ! ----------------------------------------------------------------------
2176 RCSOIL = MAX (RCSOIL,0.0001)
2177
2178 RC = RSMIN / (XLAI * RCS * RCT * RCQ * RCSOIL)
2179 ! RR = (4.* SIGMA * RD / CP)* (SFCTMP **4.)/ (SFCPRS * CH) + 1.0
2180 RR = (4.* EMISSI *SIGMA * RD / CP)* (SFCTMP **4.)/ (SFCPRS * CH) &
2181 + 1.0
2182
2183 DELTA = (SLV / CP)* DQSDT2
2184
2185 PC = (RR + DELTA)/ (RR * (1. + RC * CH) + DELTA)
2186
2187 ! ----------------------------------------------------------------------
2188 END SUBROUTINE CANRES
2189 ! ----------------------------------------------------------------------
2190
2191 SUBROUTINE CSNOW (SNCOND,DSNOW)
2192
2193 ! ----------------------------------------------------------------------
2194 ! SUBROUTINE CSNOW
2195 ! FUNCTION CSNOW
2196 ! ----------------------------------------------------------------------
2197 ! CALCULATE SNOW TERMAL CONDUCTIVITY
2198 ! ----------------------------------------------------------------------
2199 IMPLICIT NONE
2200 REAL, INTENT(IN) :: DSNOW
2201 REAL, INTENT(OUT):: SNCOND
2202 REAL :: C
2203 REAL, PARAMETER :: UNIT = 0.11631
2204
2205 ! ----------------------------------------------------------------------
2206 ! SNCOND IN UNITS OF CAL/(CM*HR*C), RETURNED IN W/(M*C)
2207 ! CSNOW IN UNITS OF CAL/(CM*HR*C), RETURNED IN W/(M*C)
2208 ! BASIC VERSION IS DYACHKOVA EQUATION (1960), FOR RANGE 0.1-0.4
2209 ! ----------------------------------------------------------------------
2210 C = 0.328*10** (2.25* DSNOW)
2211 ! CSNOW=UNIT*C
2212
2213 ! ----------------------------------------------------------------------
2214 ! DE VAUX EQUATION (1933), IN RANGE 0.1-0.6
2215 ! ----------------------------------------------------------------------
2216 ! SNCOND=0.0293*(1.+100.*DSNOW**2)
2217 ! CSNOW=0.0293*(1.+100.*DSNOW**2)
2218
2219 ! ----------------------------------------------------------------------
2220 ! E. ANDERSEN FROM FLERCHINGER
2221 ! ----------------------------------------------------------------------
2222 ! SNCOND=0.021+2.51*DSNOW**2
2223 ! CSNOW=0.021+2.51*DSNOW**2
2224
2225 SNCOND = UNIT * C
2226
2227 ! ----------------------------------------------------------------------
2228 END SUBROUTINE CSNOW
2229 ! ----------------------------------------------------------------------
2230
2231 SUBROUTINE DEVAP (EDIR,ETP1,SMC,ZSOIL,SHDFAC,SMCMAX,BEXP, &
2232 DKSAT,DWSAT,SMCDRY,SMCREF,SMCWLT,FXEXP)
2233
2234 ! ----------------------------------------------------------------------
2235 ! SUBROUTINE DEVAP
2236 ! FUNCTION DEVAP
2237 ! ----------------------------------------------------------------------
2238 ! CALCULATE DIRECT SOIL EVAPORATION
2239 ! ----------------------------------------------------------------------
2240 IMPLICIT NONE
2241 REAL, INTENT(IN) :: ETP1,SMC,BEXP,DKSAT,DWSAT,FXEXP, &
2242 SHDFAC,SMCDRY,SMCMAX,ZSOIL,SMCREF,SMCWLT
2243 REAL, INTENT(OUT):: EDIR
2244 REAL :: FX, SRATIO
2245
2246
2247 ! ----------------------------------------------------------------------
2248 ! DIRECT EVAP A FUNCTION OF RELATIVE SOIL MOISTURE AVAILABILITY, LINEAR
2249 ! WHEN FXEXP=1.
2250 ! ----------------------------------------------------------------------
2251 ! ----------------------------------------------------------------------
2252 ! FX > 1 REPRESENTS DEMAND CONTROL
2253 ! FX < 1 REPRESENTS FLUX CONTROL
2254 ! ----------------------------------------------------------------------
2255
2256 SRATIO = (SMC - SMCDRY) / (SMCMAX - SMCDRY)
2257 IF (SRATIO > 0.) THEN
2258 FX = SRATIO**FXEXP
2259 FX = MAX ( MIN ( FX, 1. ) ,0. )
2260 ELSE
2261 FX = 0.
2262 ENDIF
2263
2264 ! ----------------------------------------------------------------------
2265 ! ALLOW FOR THE DIRECT-EVAP-REDUCING EFFECT OF SHADE
2266 ! ----------------------------------------------------------------------
2267 EDIR = FX * ( 1.0- SHDFAC ) * ETP1
2268
2269 ! ----------------------------------------------------------------------
2270 END SUBROUTINE DEVAP
2271 ! ----------------------------------------------------------------------
2272
2273 SUBROUTINE EVAPO (ETA1,SMC,NSOIL,CMC,ETP1,DT,ZSOIL, &
2274 SH2O, &
2275 SMCMAX,BEXP,PC,SMCWLT,DKSAT,DWSAT, &
2276 SMCREF,SHDFAC,CMCMAX, &
2277 SMCDRY,CFACTR, &
2278 EDIR,EC,ET,ETT,SFCTMP,Q2,NROOT,RTDIS,FXEXP)
2279
2280 ! ----------------------------------------------------------------------
2281 ! SUBROUTINE EVAPO
2282 ! ----------------------------------------------------------------------
2283 ! CALCULATE SOIL MOISTURE FLUX. THE SOIL MOISTURE CONTENT (SMC - A PER
2284 ! UNIT VOLUME MEASUREMENT) IS A DEPENDENT VARIABLE THAT IS UPDATED WITH
2285 ! PROGNOSTIC EQNS. THE CANOPY MOISTURE CONTENT (CMC) IS ALSO UPDATED.
2286 ! FROZEN GROUND VERSION: NEW STATES ADDED: SH2O, AND FROZEN GROUND
2287 ! CORRECTION FACTOR, FRZFACT AND PARAMETER SLOPE.
2288 ! ----------------------------------------------------------------------
2289 IMPLICIT NONE
2290 INTEGER, INTENT(IN) :: NSOIL, NROOT
2291 INTEGER :: I,K
2292 REAL, INTENT(IN) :: BEXP, CFACTR,CMC,CMCMAX,DKSAT, &
2293 DT,DWSAT,ETP1,FXEXP,PC,Q2,SFCTMP, &
2294 SHDFAC,SMCDRY,SMCMAX,SMCREF,SMCWLT
2295 REAL, INTENT(OUT) :: EC,EDIR,ETA1,ETT
2296 REAL :: CMC2MS
2297 REAL,DIMENSION(1:NSOIL), INTENT(IN) :: RTDIS, SMC, SH2O, ZSOIL
2298 REAL,DIMENSION(1:NSOIL), INTENT(OUT) :: ET
2299
2300 ! ----------------------------------------------------------------------
2301 ! EXECUTABLE CODE BEGINS HERE IF THE POTENTIAL EVAPOTRANSPIRATION IS
2302 ! GREATER THAN ZERO.
2303 ! ----------------------------------------------------------------------
2304 EDIR = 0.
2305 EC = 0.
2306 ETT = 0.
2307 DO K = 1,NSOIL
2308 ET (K) = 0.
2309 END DO
2310
2311 ! print*,'SHDFAC',SHDFAC,'SMCMAX', SMCMAX,'BEXP',BEXP,
2312 ! & 'DKSAT',DKSAT,'DWSAT',DWSAT,'DRY',SMCDRY,'REF',SMCREF,
2313 ! & 'WLT',SMCWLT,'FXEX',FXEXP,'SMC',SMC
2314 ! ----------------------------------------------------------------------
2315 ! RETRIEVE DIRECT EVAPORATION FROM SOIL SURFACE. CALL THIS FUNCTION
2316 ! ONLY IF VEG COVER NOT COMPLETE.
2317 ! FROZEN GROUND VERSION: SH2O STATES REPLACE SMC STATES.
2318 ! ----------------------------------------------------------------------
2319 IF (ETP1 > 0.0) THEN
2320 IF (SHDFAC < 1.) THEN
2321 ! CALL DEVAP (EDIR,ETP1,SH2O (1),ZSOIL (1),SHDFAC,SMCMAX,
2322 CALL DEVAP (EDIR,ETP1,SMC (1),ZSOIL (1),SHDFAC,SMCMAX, &
2323 BEXP,DKSAT,DWSAT,SMCDRY,SMCREF,SMCWLT,FXEXP)
2324 END IF
2325 ! ----------------------------------------------------------------------
2326 ! INITIALIZE PLANT TOTAL TRANSPIRATION, RETRIEVE PLANT TRANSPIRATION,
2327 ! AND ACCUMULATE IT FOR ALL SOIL LAYERS.
2328 ! ----------------------------------------------------------------------
2329
2330 IF (SHDFAC > 0.0) THEN
2331 CALL TRANSP (ET,NSOIL,ETP1,SH2O,CMC,ZSOIL,SHDFAC,SMCWLT, &
2332 CMCMAX,PC,CFACTR,SMCREF,SFCTMP,Q2,NROOT,RTDIS)
2333 DO K = 1,NSOIL
2334 ETT = ETT + ET ( K )
2335 END DO
2336 ! ----------------------------------------------------------------------
2337 ! CALCULATE CANOPY EVAPORATION.
2338 ! IF STATEMENTS TO AVOID TANGENT LINEAR PROBLEMS NEAR CMC=0.0.
2339 ! ----------------------------------------------------------------------
2340 IF (CMC > 0.0) THEN
2341 EC = SHDFAC * ( ( CMC / CMCMAX ) ** CFACTR ) * ETP1
2342 ELSE
2343 EC = 0.0
2344 END IF
2345 ! ----------------------------------------------------------------------
2346 ! EC SHOULD BE LIMITED BY THE TOTAL AMOUNT OF AVAILABLE WATER ON THE
2347 ! CANOPY. -F.CHEN, 18-OCT-1994
2348 ! ----------------------------------------------------------------------
2349 CMC2MS = CMC / DT
2350 EC = MIN ( CMC2MS, EC )
2351 END IF
2352 END IF
2353 ! ----------------------------------------------------------------------
2354 ! TOTAL UP EVAP AND TRANSP TYPES TO OBTAIN ACTUAL EVAPOTRANSP
2355 ! ----------------------------------------------------------------------
2356 ETA1 = EDIR + ETT + EC
2357
2358 ! ----------------------------------------------------------------------
2359 END SUBROUTINE EVAPO
2360 ! ----------------------------------------------------------------------
2361
2362 SUBROUTINE FRH2O (FREE,TKELV,SMC,SH2O,SMCMAX,BEXP,PSIS)
2363
2364 ! ----------------------------------------------------------------------
2365 ! SUBROUTINE FRH2O
2366 ! ----------------------------------------------------------------------
2367 ! CALCULATE AMOUNT OF SUPERCOOLED LIQUID SOIL WATER CONTENT IF
2368 ! TEMPERATURE IS BELOW 273.15K (T0). REQUIRES NEWTON-TYPE ITERATION TO
2369 ! SOLVE THE NONLINEAR IMPLICIT EQUATION GIVEN IN EQN 17 OF KOREN ET AL
2370 ! (1999, JGR, VOL 104(D16), 19569-19585).
2371 ! ----------------------------------------------------------------------
2372 ! NEW VERSION (JUNE 2001): MUCH FASTER AND MORE ACCURATE NEWTON
2373 ! ITERATION ACHIEVED BY FIRST TAKING LOG OF EQN CITED ABOVE -- LESS THAN
2374 ! 4 (TYPICALLY 1 OR 2) ITERATIONS ACHIEVES CONVERGENCE. ALSO, EXPLICIT
2375 ! 1-STEP SOLUTION OPTION FOR SPECIAL CASE OF PARAMETER CK=0, WHICH
2376 ! REDUCES THE ORIGINAL IMPLICIT EQUATION TO A SIMPLER EXPLICIT FORM,
2377 ! KNOWN AS THE "FLERCHINGER EQN". IMPROVED HANDLING OF SOLUTION IN THE
2378 ! LIMIT OF FREEZING POINT TEMPERATURE T0.
2379 ! ----------------------------------------------------------------------
2380 ! INPUT:
2381
2382 ! TKELV.........TEMPERATURE (Kelvin)
2383 ! SMC...........TOTAL SOIL MOISTURE CONTENT (VOLUMETRIC)
2384 ! SH2O..........LIQUID SOIL MOISTURE CONTENT (VOLUMETRIC)
2385 ! SMCMAX........SATURATION SOIL MOISTURE CONTENT (FROM REDPRM)
2386 ! B.............SOIL TYPE "B" PARAMETER (FROM REDPRM)
2387 ! PSIS..........SATURATED SOIL MATRIC POTENTIAL (FROM REDPRM)
2388
2389 ! OUTPUT:
2390 ! FRH2O.........SUPERCOOLED LIQUID WATER CONTENT
2391 ! FREE..........SUPERCOOLED LIQUID WATER CONTENT
2392 ! ----------------------------------------------------------------------
2393 IMPLICIT NONE
2394 REAL, INTENT(IN) :: BEXP,PSIS,SH2O,SMC,SMCMAX,TKELV
2395 REAL, INTENT(OUT) :: FREE
2396 REAL :: BX,DENOM,DF,DSWL,FK,SWL,SWLK
2397 INTEGER :: NLOG,KCOUNT
2398 ! PARAMETER(CK = 0.0)
2399 REAL, PARAMETER :: CK = 8.0, BLIM = 5.5, ERROR = 0.005, &
2400 HLICE = 3.335E5, GS = 9.81,DICE = 920.0, &
2401 DH2O = 1000.0, T0 = 273.15
2402
2403 ! ----------------------------------------------------------------------
2404 ! LIMITS ON PARAMETER B: B < 5.5 (use parameter BLIM)
2405 ! SIMULATIONS SHOWED IF B > 5.5 UNFROZEN WATER CONTENT IS
2406 ! NON-REALISTICALLY HIGH AT VERY LOW TEMPERATURES.
2407 ! ----------------------------------------------------------------------
2408 BX = BEXP
2409
2410 ! ----------------------------------------------------------------------
2411 ! INITIALIZING ITERATIONS COUNTER AND ITERATIVE SOLUTION FLAG.
2412 ! ----------------------------------------------------------------------
2413 IF (BEXP > BLIM) BX = BLIM
2414 NLOG = 0
2415
2416 ! ----------------------------------------------------------------------
2417 ! IF TEMPERATURE NOT SIGNIFICANTLY BELOW FREEZING (T0), SH2O = SMC
2418 ! ----------------------------------------------------------------------
2419 KCOUNT = 0
2420 ! FRH2O = SMC
2421 IF (TKELV > (T0- 1.E-3)) THEN
2422 FREE = SMC
2423 ELSE
2424
2425 ! ----------------------------------------------------------------------
2426 ! OPTION 1: ITERATED SOLUTION FOR NONZERO CK
2427 ! IN KOREN ET AL, JGR, 1999, EQN 17
2428 ! ----------------------------------------------------------------------
2429 ! INITIAL GUESS FOR SWL (frozen content)
2430 ! ----------------------------------------------------------------------
2431 IF (CK /= 0.0) THEN
2432 SWL = SMC - SH2O
2433 ! ----------------------------------------------------------------------
2434 ! KEEP WITHIN BOUNDS.
2435 ! ----------------------------------------------------------------------
2436 IF (SWL > (SMC -0.02)) SWL = SMC -0.02
2437
2438 ! ----------------------------------------------------------------------
2439 ! START OF ITERATIONS
2440 ! ----------------------------------------------------------------------
2441 IF (SWL < 0.) SWL = 0.
2442 1001 Continue
2443 IF (.NOT.( (NLOG < 10) .AND. (KCOUNT == 0))) goto 1002
2444 NLOG = NLOG +1
2445 DF = ALOG ( ( PSIS * GS / HLICE ) * ( ( 1. + CK * SWL )**2.) * &
2446 ( SMCMAX / (SMC - SWL) )** BX) - ALOG ( - ( &
2447 TKELV - T0)/ TKELV)
2448 DENOM = 2. * CK / ( 1. + CK * SWL ) + BX / ( SMC - SWL )
2449 SWLK = SWL - DF / DENOM
2450 ! ----------------------------------------------------------------------
2451 ! BOUNDS USEFUL FOR MATHEMATICAL SOLUTION.
2452 ! ----------------------------------------------------------------------
2453 IF (SWLK > (SMC -0.02)) SWLK = SMC - 0.02
2454 IF (SWLK < 0.) SWLK = 0.
2455
2456 ! ----------------------------------------------------------------------
2457 ! MATHEMATICAL SOLUTION BOUNDS APPLIED.
2458 ! ----------------------------------------------------------------------
2459 DSWL = ABS (SWLK - SWL)
2460
2461 ! ----------------------------------------------------------------------
2462 ! IF MORE THAN 10 ITERATIONS, USE EXPLICIT METHOD (CK=0 APPROX.)
2463 ! WHEN DSWL LESS OR EQ. ERROR, NO MORE ITERATIONS REQUIRED.
2464 ! ----------------------------------------------------------------------
2465 SWL = SWLK
2466 IF ( DSWL <= ERROR ) THEN
2467 KCOUNT = KCOUNT +1
2468 END IF
2469 ! ----------------------------------------------------------------------
2470 ! END OF ITERATIONS
2471 ! ----------------------------------------------------------------------
2472 ! BOUNDS APPLIED WITHIN DO-BLOCK ARE VALID FOR PHYSICAL SOLUTION.
2473 ! ----------------------------------------------------------------------
2474 ! FRH2O = SMC - SWL
2475 goto 1001
2476 1002 continue
2477 FREE = SMC - SWL
2478 END IF
2479 ! ----------------------------------------------------------------------
2480 ! END OPTION 1
2481 ! ----------------------------------------------------------------------
2482 ! ----------------------------------------------------------------------
2483 ! OPTION 2: EXPLICIT SOLUTION FOR FLERCHINGER EQ. i.e. CK=0
2484 ! IN KOREN ET AL., JGR, 1999, EQN 17
2485 ! APPLY PHYSICAL BOUNDS TO FLERCHINGER SOLUTION
2486 ! ----------------------------------------------------------------------
2487 IF (KCOUNT == 0) THEN
2488 PRINT *,'Flerchinger USEd in NEW version. Iterations=',NLOG
2489 FK = ( ( (HLICE / (GS * ( - PSIS)))* &
2490 ( (TKELV - T0)/ TKELV))** ( -1/ BX))* SMCMAX
2491 ! FRH2O = MIN (FK, SMC)
2492 IF (FK < 0.02) FK = 0.02
2493 FREE = MIN (FK, SMC)
2494 ! ----------------------------------------------------------------------
2495 ! END OPTION 2
2496 ! ----------------------------------------------------------------------
2497 END IF
2498 END IF
2499 ! ----------------------------------------------------------------------
2500 END SUBROUTINE FRH2O
2501 ! ----------------------------------------------------------------------
2502
2503 SUBROUTINE HRT (RHSTS,STC,SMC,SMCMAX,NSOIL,ZSOIL,YY,ZZ1, &
2504 TBOT,ZBOT,PSISAT,SH2O,DT,BEXP, &
2505 F1,DF1,QUARTZ,CSOIL,AI,BI,CI,VEGTYP)
2506
2507 ! ----------------------------------------------------------------------
2508 ! SUBROUTINE HRT
2509 ! ----------------------------------------------------------------------
2510 ! CALCULATE THE RIGHT HAND SIDE OF THE TIME TENDENCY TERM OF THE SOIL
2511 ! THERMAL DIFFUSION EQUATION. ALSO TO COMPUTE ( PREPARE ) THE MATRIX
2512 ! COEFFICIENTS FOR THE TRI-DIAGONAL MATRIX OF THE IMPLICIT TIME SCHEME.
2513 ! ----------------------------------------------------------------------
2514 IMPLICIT NONE
2515 LOGICAL :: ITAVG
2516 INTEGER, INTENT(IN) :: NSOIL, VEGTYP
2517 INTEGER :: I, K
2518
2519 REAL, INTENT(IN) :: BEXP, CSOIL, DF1, DT,F1,PSISAT,QUARTZ, &
2520 SMCMAX ,TBOT,YY,ZZ1, ZBOT
2521 REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SMC,STC,ZSOIL
2522 REAL, DIMENSION(1:NSOIL), INTENT(INOUT):: SH2O
2523 REAL, DIMENSION(1:NSOIL), INTENT(OUT) :: RHSTS
2524 REAL, DIMENSION(1:NSOIL), INTENT(OUT) :: AI, BI,CI
2525 REAL :: DDZ, DDZ2, DENOM, DF1N, DF1K, DTSDZ, &
2526 DTSDZ2,HCPCT,QTOT,SSOIL,SICE,TAVG,TBK, &
2527 TBK1,TSNSR,TSURF,CSOIL_LOC
2528 REAL, PARAMETER :: T0 = 273.15, CAIR = 1004.0, CICE = 2.106E6,&
2529 CH2O = 4.2E6
2530
2531
2532 !urban
2533 IF(VEGTYP==1) then
2534 CSOIL_LOC=3.0E6
2535 ELSE
2536 CSOIL_LOC=CSOIL
2537 ENDIF
2538
2539 ! ----------------------------------------------------------------------
2540 ! INITIALIZE LOGICAL FOR SOIL LAYER TEMPERATURE AVERAGING.
2541 ! ----------------------------------------------------------------------
2542 ITAVG = .TRUE.
2543 ! ----------------------------------------------------------------------
2544 ! BEGIN SECTION FOR TOP SOIL LAYER
2545 ! ----------------------------------------------------------------------
2546 ! CALC THE HEAT CAPACITY OF THE TOP SOIL LAYER
2547 ! ----------------------------------------------------------------------
2548 HCPCT = SH2O (1)* CH2O + (1.0- SMCMAX)* CSOIL_LOC + (SMCMAX - SMC (1))&
2549 * CAIR &
2550 + ( SMC (1) - SH2O (1) )* CICE
2551
2552 ! ----------------------------------------------------------------------
2553 ! CALC THE MATRIX COEFFICIENTS AI, BI, AND CI FOR THE TOP LAYER
2554 ! ----------------------------------------------------------------------
2555 DDZ = 1.0 / ( -0.5 * ZSOIL (2) )
2556 AI (1) = 0.0
2557 CI (1) = (DF1 * DDZ) / (ZSOIL (1) * HCPCT)
2558
2559 ! ----------------------------------------------------------------------
2560 ! CALCULATE THE VERTICAL SOIL TEMP GRADIENT BTWN THE 1ST AND 2ND SOIL
2561 ! LAYERS. THEN CALCULATE THE SUBSURFACE HEAT FLUX. USE THE TEMP
2562 ! GRADIENT AND SUBSFC HEAT FLUX TO CALC "RIGHT-HAND SIDE TENDENCY
2563 ! TERMS", OR "RHSTS", FOR TOP SOIL LAYER.
2564 ! ----------------------------------------------------------------------
2565 BI (1) = - CI (1) + DF1 / (0.5 * ZSOIL (1) * ZSOIL (1)* HCPCT * &
2566 ZZ1)
2567 DTSDZ = (STC (1) - STC (2)) / ( -0.5 * ZSOIL (2))
2568 SSOIL = DF1 * (STC (1) - YY) / (0.5 * ZSOIL (1) * ZZ1)
2569 ! RHSTS(1) = (DF1 * DTSDZ - SSOIL) / (ZSOIL(1) * HCPCT)
2570 DENOM = (ZSOIL (1) * HCPCT)
2571
2572 ! ----------------------------------------------------------------------
2573 ! NEXT CAPTURE THE VERTICAL DIFFERENCE OF THE HEAT FLUX AT TOP AND
2574 ! BOTTOM OF FIRST SOIL LAYER FOR USE IN HEAT FLUX CONSTRAINT APPLIED TO
2575 ! POTENTIAL SOIL FREEZING/THAWING IN ROUTINE SNKSRC.
2576 ! ----------------------------------------------------------------------
2577 ! QTOT = SSOIL - DF1*DTSDZ
2578 RHSTS (1) = (DF1 * DTSDZ - SSOIL) / DENOM
2579
2580 ! ----------------------------------------------------------------------
2581 ! CALCULATE FROZEN WATER CONTENT IN 1ST SOIL LAYER.
2582 ! ----------------------------------------------------------------------
2583 QTOT = -1.0* RHSTS (1)* DENOM
2584
2585 ! ----------------------------------------------------------------------
2586 ! IF TEMPERATURE AVERAGING INVOKED (ITAVG=TRUE; ELSE SKIP):
2587 ! SET TEMP "TSURF" AT TOP OF SOIL COLUMN (FOR USE IN FREEZING SOIL
2588 ! PHYSICS LATER IN FUNCTION SUBROUTINE SNKSRC). IF SNOWPACK CONTENT IS
2589 ! ZERO, THEN TSURF EXPRESSION BELOW GIVES TSURF = SKIN TEMP. IF
2590 ! SNOWPACK IS NONZERO (HENCE ARGUMENT ZZ1=1), THEN TSURF EXPRESSION
2591 ! BELOW YIELDS SOIL COLUMN TOP TEMPERATURE UNDER SNOWPACK. THEN
2592 ! CALCULATE TEMPERATURE AT BOTTOM INTERFACE OF 1ST SOIL LAYER FOR USE
2593 ! LATER IN FUNCTION SUBROUTINE SNKSRC
2594 ! ----------------------------------------------------------------------
2595 SICE = SMC (1) - SH2O (1)
2596 IF (ITAVG) THEN
2597 TSURF = (YY + (ZZ1-1) * STC (1)) / ZZ1
2598 ! ----------------------------------------------------------------------
2599 ! IF FROZEN WATER PRESENT OR ANY OF LAYER-1 MID-POINT OR BOUNDING
2600 ! INTERFACE TEMPERATURES BELOW FREEZING, THEN CALL SNKSRC TO
2601 ! COMPUTE HEAT SOURCE/SINK (AND CHANGE IN FROZEN WATER CONTENT)
2602 ! DUE TO POSSIBLE SOIL WATER PHASE CHANGE
2603 ! ----------------------------------------------------------------------
2604 CALL TBND (STC (1),STC (2),ZSOIL,ZBOT,1,NSOIL,TBK)
2605 IF ( (SICE > 0.) .OR. (STC (1) < T0) .OR. &
2606 (TSURF < T0) .OR. (TBK < T0) ) THEN
2607 ! TSNSR = SNKSRC (TAVG,SMC(1),SH2O(1),
2608 CALL TMPAVG (TAVG,TSURF,STC (1),TBK,ZSOIL,NSOIL,1)
2609 CALL SNKSRC (TSNSR,TAVG,SMC (1),SH2O (1), &
2610 ZSOIL,NSOIL,SMCMAX,PSISAT,BEXP,DT,1,QTOT)
2611 ! RHSTS(1) = RHSTS(1) - TSNSR / ( ZSOIL(1) * HCPCT )
2612 RHSTS (1) = RHSTS (1) - TSNSR / DENOM
2613 END IF
2614 ELSE
2615 ! TSNSR = SNKSRC (STC(1),SMC(1),SH2O(1),
2616 IF ( (SICE > 0.) .OR. (STC (1) < T0) ) THEN
2617 CALL SNKSRC (TSNSR,STC (1),SMC (1),SH2O (1), &
2618 ZSOIL,NSOIL,SMCMAX,PSISAT,BEXP,DT,1,QTOT)
2619 ! RHSTS(1) = RHSTS(1) - TSNSR / ( ZSOIL(1) * HCPCT )
2620 RHSTS (1) = RHSTS (1) - TSNSR / DENOM
2621 END IF
2622 ! ----------------------------------------------------------------------
2623 ! THIS ENDS SECTION FOR TOP SOIL LAYER.
2624 ! ----------------------------------------------------------------------
2625 END IF
2626
2627 ! INITIALIZE DDZ2
2628 ! ----------------------------------------------------------------------
2629
2630 DDZ2 = 0.0
2631 DF1K = DF1
2632
2633 ! ----------------------------------------------------------------------
2634 ! LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABOVE PROCESS
2635 ! (EXCEPT SUBSFC OR "GROUND" HEAT FLUX NOT REPEATED IN LOWER LAYERS)
2636 ! ----------------------------------------------------------------------
2637 ! CALCULATE HEAT CAPACITY FOR THIS SOIL LAYER.
2638 ! ----------------------------------------------------------------------
2639 DO K = 2,NSOIL
2640 HCPCT = SH2O (K)* CH2O + (1.0- SMCMAX)* CSOIL_LOC + (SMCMAX - SMC ( &
2641 K))* CAIR + ( SMC (K) - SH2O (K) )* CICE
2642 ! ----------------------------------------------------------------------
2643 ! THIS SECTION FOR LAYER 2 OR GREATER, BUT NOT LAST LAYER.
2644 ! ----------------------------------------------------------------------
2645 ! CALCULATE THERMAL DIFFUSIVITY FOR THIS LAYER.
2646 ! ----------------------------------------------------------------------
2647 IF (K /= NSOIL) THEN
2648
2649 ! ----------------------------------------------------------------------
2650 ! CALC THE VERTICAL SOIL TEMP GRADIENT THRU THIS LAYER
2651 ! ----------------------------------------------------------------------
2652 CALL TDFCND (DF1N,SMC (K),QUARTZ,SMCMAX,SH2O (K))
2653
2654 !urban
2655 IF(VEGTYP==1) DF1N = 3.24
2656
2657 DENOM = 0.5 * ( ZSOIL (K -1) - ZSOIL (K +1) )
2658
2659 ! ----------------------------------------------------------------------
2660 ! CALC THE MATRIX COEF, CI, AFTER CALC'NG ITS PARTIAL PRODUCT
2661 ! ----------------------------------------------------------------------
2662 DTSDZ2 = ( STC (K) - STC (K +1) ) / DENOM
2663 DDZ2 = 2. / (ZSOIL (K -1) - ZSOIL (K +1))
2664
2665 ! ----------------------------------------------------------------------
2666 ! IF TEMPERATURE AVERAGING INVOKED (ITAVG=TRUE; ELSE SKIP): CALCULATE
2667 ! TEMP AT BOTTOM OF LAYER.
2668 ! ----------------------------------------------------------------------
2669 CI (K) = - DF1N * DDZ2 / ( (ZSOIL (K -1) - ZSOIL (K)) * &
2670 HCPCT)
2671 IF (ITAVG) THEN
2672 CALL TBND (STC (K),STC (K +1),ZSOIL,ZBOT,K,NSOIL,TBK1)
2673 END IF
2674
2675 ELSE
2676 ! ----------------------------------------------------------------------
2677 ! SPECIAL CASE OF BOTTOM SOIL LAYER: CALCULATE THERMAL DIFFUSIVITY FOR
2678 ! BOTTOM LAYER.
2679 ! ----------------------------------------------------------------------
2680
2681 ! ----------------------------------------------------------------------
2682 ! CALC THE VERTICAL SOIL TEMP GRADIENT THRU BOTTOM LAYER.
2683 ! ----------------------------------------------------------------------
2684 CALL TDFCND (DF1N,SMC (K),QUARTZ,SMCMAX,SH2O (K))
2685
2686
2687 !urban
2688 IF(VEGTYP==1) DF1N = 3.24
2689
2690 DENOM = .5 * (ZSOIL (K -1) + ZSOIL (K)) - ZBOT
2691
2692 ! ----------------------------------------------------------------------
2693 ! SET MATRIX COEF, CI TO ZERO IF BOTTOM LAYER.
2694 ! ----------------------------------------------------------------------
2695 DTSDZ2 = (STC (K) - TBOT) / DENOM
2696
2697 ! ----------------------------------------------------------------------
2698 ! IF TEMPERATURE AVERAGING INVOKED (ITAVG=TRUE; ELSE SKIP): CALCULATE
2699 ! TEMP AT BOTTOM OF LAST LAYER.
2700 ! ----------------------------------------------------------------------
2701 CI (K) = 0.
2702 IF (ITAVG) THEN
2703 CALL TBND (STC (K),TBOT,ZSOIL,ZBOT,K,NSOIL,TBK1)
2704 END IF
2705 ! ----------------------------------------------------------------------
2706 ! THIS ENDS SPECIAL LOOP FOR BOTTOM LAYER.
2707 END IF
2708 ! ----------------------------------------------------------------------
2709 ! CALCULATE RHSTS FOR THIS LAYER AFTER CALC'NG A PARTIAL PRODUCT.
2710 ! ----------------------------------------------------------------------
2711 DENOM = ( ZSOIL (K) - ZSOIL (K -1) ) * HCPCT
2712 RHSTS (K) = ( DF1N * DTSDZ2- DF1K * DTSDZ ) / DENOM
2713 QTOT = -1.0* DENOM * RHSTS (K)
2714
2715 SICE = SMC (K) - SH2O (K)
2716 IF (ITAVG) THEN
2717 CALL TMPAVG (TAVG,TBK,STC (K),TBK1,ZSOIL,NSOIL,K)
2718 ! TSNSR = SNKSRC(TAVG,SMC(K),SH2O(K),ZSOIL,NSOIL,
2719 IF ( (SICE > 0.) .OR. (STC (K) < T0) .OR. &
2720 (TBK .lt. T0) .OR. (TBK1 .lt. T0) ) THEN
2721 CALL SNKSRC (TSNSR,TAVG,SMC (K),SH2O (K),ZSOIL,NSOIL, &
2722 SMCMAX,PSISAT,BEXP,DT,K,QTOT)
2723 RHSTS (K) = RHSTS (K) - TSNSR / DENOM
2724 END IF
2725 ELSE
2726 ! TSNSR = SNKSRC(STC(K),SMC(K),SH2O(K),ZSOIL,NSOIL,
2727 IF ( (SICE > 0.) .OR. (STC (K) < T0) ) THEN
2728 CALL SNKSRC (TSNSR,STC (K),SMC (K),SH2O (K),ZSOIL,NSOIL, &
2729 SMCMAX,PSISAT,BEXP,DT,K,QTOT)
2730 RHSTS (K) = RHSTS (K) - TSNSR / DENOM
2731 END IF
2732 END IF
2733
2734 ! ----------------------------------------------------------------------
2735 ! CALC MATRIX COEFS, AI, AND BI FOR THIS LAYER.
2736 ! ----------------------------------------------------------------------
2737 AI (K) = - DF1 * DDZ / ( (ZSOIL (K -1) - ZSOIL (K)) * HCPCT)
2738
2739 ! ----------------------------------------------------------------------
2740 ! RESET VALUES OF DF1, DTSDZ, DDZ, AND TBK FOR LOOP TO NEXT SOIL LAYER.
2741 ! ----------------------------------------------------------------------
2742 BI (K) = - (AI (K) + CI (K))
2743 TBK = TBK1
2744 DF1K = DF1N
2745 DTSDZ = DTSDZ2
2746 DDZ = DDZ2
2747 END DO
2748 ! ----------------------------------------------------------------------
2749 END SUBROUTINE HRT
2750 ! ----------------------------------------------------------------------
2751
2752 SUBROUTINE HRTICE (RHSTS,STC,NSOIL,ZSOIL,YY,ZZ1,DF1,AI,BI,CI)
2753
2754 ! ----------------------------------------------------------------------
2755 ! SUBROUTINE HRTICE
2756 ! ----------------------------------------------------------------------
2757 ! CALCULATE THE RIGHT HAND SIDE OF THE TIME TENDENCY TERM OF THE SOIL
2758 ! THERMAL DIFFUSION EQUATION IN THE CASE OF SEA-ICE PACK. ALSO TO
2759 ! COMPUTE (PREPARE) THE MATRIX COEFFICIENTS FOR THE TRI-DIAGONAL MATRIX
2760 ! OF THE IMPLICIT TIME SCHEME.
2761 ! ----------------------------------------------------------------------
2762 IMPLICIT NONE
2763
2764
2765 INTEGER, INTENT(IN) :: NSOIL
2766 INTEGER :: K
2767
2768 REAL, INTENT(IN) :: DF1,YY,ZZ1
2769 REAL, DIMENSION(1:NSOIL), INTENT(OUT):: AI, BI,CI
2770 REAL, DIMENSION(1:NSOIL), INTENT(IN) :: STC, ZSOIL
2771 REAL, DIMENSION(1:NSOIL), INTENT(OUT):: RHSTS
2772 REAL :: DDZ,DDZ2,DENOM,DTSDZ,DTSDZ2,SSOIL, &
2773 ZBOT,TBOT
2774 REAL, PARAMETER :: HCPCT = 1.72396E+6
2775
2776 ! ----------------------------------------------------------------------
2777 ! SET A NOMINAL UNIVERSAL VALUE OF THE SEA-ICE SPECIFIC HEAT CAPACITY,
2778 ! HCPCT = 1880.0*917.0.
2779 ! ----------------------------------------------------------------------
2780 DATA TBOT /271.16/
2781
2782 ! ----------------------------------------------------------------------
2783 ! THE INPUT ARGUMENT DF1 IS A UNIVERSALLY CONSTANT VALUE OF SEA-ICE
2784 ! THERMAL DIFFUSIVITY, SET IN ROUTINE SNOPAC AS DF1 = 2.2.
2785 ! ----------------------------------------------------------------------
2786 ! SET ICE PACK DEPTH. USE TBOT AS ICE PACK LOWER BOUNDARY TEMPERATURE
2787 ! (THAT OF UNFROZEN SEA WATER AT BOTTOM OF SEA ICE PACK). ASSUME ICE
2788 ! PACK IS OF N=NSOIL LAYERS SPANNING A UNIFORM CONSTANT ICE PACK
2789 ! THICKNESS AS DEFINED BY ZSOIL(NSOIL) IN ROUTINE SFLX.
2790 ! ----------------------------------------------------------------------
2791 ! ----------------------------------------------------------------------
2792 ! CALC THE MATRIX COEFFICIENTS AI, BI, AND CI FOR THE TOP LAYER
2793 ! ----------------------------------------------------------------------
2794 ZBOT = ZSOIL (NSOIL)
2795 DDZ = 1.0 / ( -0.5 * ZSOIL (2) )
2796 AI (1) = 0.0
2797 CI (1) = (DF1 * DDZ) / (ZSOIL (1) * HCPCT)
2798
2799 ! ----------------------------------------------------------------------
2800 ! CALC THE VERTICAL SOIL TEMP GRADIENT BTWN THE TOP AND 2ND SOIL LAYERS.
2801 ! RECALC/ADJUST THE SOIL HEAT FLUX. USE THE GRADIENT AND FLUX TO CALC
2802 ! RHSTS FOR THE TOP SOIL LAYER.
2803 ! ----------------------------------------------------------------------
2804 BI (1) = - CI (1) + DF1/ (0.5 * ZSOIL (1) * ZSOIL (1) * HCPCT * &
2805 ZZ1)
2806 DTSDZ = ( STC (1) - STC (2) ) / ( -0.5 * ZSOIL (2) )
2807 SSOIL = DF1 * ( STC (1) - YY ) / ( 0.5 * ZSOIL (1) * ZZ1 )
2808
2809 ! ----------------------------------------------------------------------
2810 ! INITIALIZE DDZ2
2811 ! ----------------------------------------------------------------------
2812 RHSTS (1) = ( DF1 * DTSDZ - SSOIL ) / ( ZSOIL (1) * HCPCT )
2813
2814 ! ----------------------------------------------------------------------
2815 ! LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABOVE PROCESS
2816 ! ----------------------------------------------------------------------
2817 DDZ2 = 0.0
2818 DO K = 2,NSOIL
2819
2820 ! ----------------------------------------------------------------------
2821 ! CALC THE VERTICAL SOIL TEMP GRADIENT THRU THIS LAYER.
2822 ! ----------------------------------------------------------------------
2823 IF (K /= NSOIL) THEN
2824 DENOM = 0.5 * ( ZSOIL (K -1) - ZSOIL (K +1) )
2825
2826 ! ----------------------------------------------------------------------
2827 ! CALC THE MATRIX COEF, CI, AFTER CALC'NG ITS PARTIAL PRODUCT.
2828 ! ----------------------------------------------------------------------
2829 DTSDZ2 = ( STC (K) - STC (K +1) ) / DENOM
2830 DDZ2 = 2. / (ZSOIL (K -1) - ZSOIL (K +1))
2831 CI (K) = - DF1 * DDZ2 / ( (ZSOIL (K -1) - ZSOIL (K))*HCPCT)
2832
2833 ! ----------------------------------------------------------------------
2834 ! CALC THE VERTICAL SOIL TEMP GRADIENT THRU THE LOWEST LAYER.
2835 ! ----------------------------------------------------------------------
2836 ELSE
2837
2838 ! ----------------------------------------------------------------------
2839 ! SET MATRIX COEF, CI TO ZERO.
2840 ! ----------------------------------------------------------------------
2841 DTSDZ2 = (STC (K) - TBOT)/ (.5 * (ZSOIL (K -1) + ZSOIL (K)) &
2842 - ZBOT)
2843 CI (K) = 0.
2844 ! ----------------------------------------------------------------------
2845 ! CALC RHSTS FOR THIS LAYER AFTER CALC'NG A PARTIAL PRODUCT.
2846 ! ----------------------------------------------------------------------
2847 END IF
2848 DENOM = ( ZSOIL (K) - ZSOIL (K -1) ) * HCPCT
2849
2850 ! ----------------------------------------------------------------------
2851 ! CALC MATRIX COEFS, AI, AND BI FOR THIS LAYER.
2852 ! ----------------------------------------------------------------------
2853 RHSTS (K) = ( DF1 * DTSDZ2- DF1 * DTSDZ ) / DENOM
2854 AI (K) = - DF1 * DDZ / ( (ZSOIL (K -1) - ZSOIL (K)) * HCPCT)
2855
2856 ! ----------------------------------------------------------------------
2857 ! RESET VALUES OF DTSDZ AND DDZ FOR LOOP TO NEXT SOIL LYR.
2858 ! ----------------------------------------------------------------------
2859 BI (K) = - (AI (K) + CI (K))
2860 DTSDZ = DTSDZ2
2861 DDZ = DDZ2
2862 END DO
2863 ! ----------------------------------------------------------------------
2864 END SUBROUTINE HRTICE
2865 ! ----------------------------------------------------------------------
2866
2867 SUBROUTINE HSTEP (STCOUT,STCIN,RHSTS,DT,NSOIL,AI,BI,CI)
2868
2869 ! ----------------------------------------------------------------------
2870 ! SUBROUTINE HSTEP
2871 ! ----------------------------------------------------------------------
2872 ! CALCULATE/UPDATE THE SOIL TEMPERATURE FIELD.
2873 ! ----------------------------------------------------------------------
2874 IMPLICIT NONE
2875 INTEGER, INTENT(IN) :: NSOIL
2876 INTEGER :: K
2877
2878 REAL, DIMENSION(1:NSOIL), INTENT(IN):: STCIN
2879 REAL, DIMENSION(1:NSOIL), INTENT(OUT):: STCOUT
2880 REAL, DIMENSION(1:NSOIL), INTENT(INOUT):: RHSTS
2881 REAL, DIMENSION(1:NSOIL), INTENT(INOUT):: AI,BI,CI
2882 REAL, DIMENSION(1:NSOIL) :: RHSTSin
2883 REAL, DIMENSION(1:NSOIL) :: CIin
2884 REAL :: DT
2885
2886 ! ----------------------------------------------------------------------
2887 ! CREATE FINITE DIFFERENCE VALUES FOR USE IN ROSR12 ROUTINE
2888 ! ----------------------------------------------------------------------
2889 DO K = 1,NSOIL
2890 RHSTS (K) = RHSTS (K) * DT
2891 AI (K) = AI (K) * DT
2892 BI (K) = 1. + BI (K) * DT
2893 CI (K) = CI (K) * DT
2894 END DO
2895 ! ----------------------------------------------------------------------
2896 ! COPY VALUES FOR INPUT VARIABLES BEFORE CALL TO ROSR12
2897 ! ----------------------------------------------------------------------
2898 DO K = 1,NSOIL
2899 RHSTSin (K) = RHSTS (K)
2900 END DO
2901 DO K = 1,NSOIL
2902 CIin (K) = CI (K)
2903 END DO
2904 ! ----------------------------------------------------------------------
2905 ! SOLVE THE TRI-DIAGONAL MATRIX EQUATION
2906 ! ----------------------------------------------------------------------
2907 CALL ROSR12 (CI,AI,BI,CIin,RHSTSin,RHSTS,NSOIL)
2908 ! ----------------------------------------------------------------------
2909 ! CALC/UPDATE THE SOIL TEMPS USING MATRIX SOLUTION
2910 ! ----------------------------------------------------------------------
2911 DO K = 1,NSOIL
2912 STCOUT (K) = STCIN (K) + CI (K)
2913 END DO
2914 ! ----------------------------------------------------------------------
2915 END SUBROUTINE HSTEP
2916 ! ----------------------------------------------------------------------
2917
2918 SUBROUTINE NOPAC (ETP,ETA,PRCP,SMC,SMCMAX,SMCWLT, &
2919 SMCREF,SMCDRY,CMC,CMCMAX,NSOIL,DT,SHDFAC, &
2920 SBETA,Q2,T1,SFCTMP,T24,TH2,FDOWN,F1,EMISSI, &
2921 SSOIL, &
2922 STC,EPSCA,BEXP,PC,RCH,RR,CFACTR, &
2923 SH2O,SLOPE,KDT,FRZFACT,PSISAT,ZSOIL, &
2924 DKSAT,DWSAT,TBOT,ZBOT,RUNOFF1,RUNOFF2, &
2925 RUNOFF3,EDIR,EC,ET,ETT,NROOT,ICE,RTDIS, &
2926 QUARTZ,FXEXP,CSOIL, &
2927 BETA,DRIP,DEW,FLX1,FLX3,VEGTYP)
2928
2929 ! ----------------------------------------------------------------------
2930 ! SUBROUTINE NOPAC
2931 ! ----------------------------------------------------------------------
2932 ! CALCULATE SOIL MOISTURE AND HEAT FLUX VALUES AND UPDATE SOIL MOISTURE
2933 ! CONTENT AND SOIL HEAT CONTENT VALUES FOR THE CASE WHEN NO SNOW PACK IS
2934 ! PRESENT.
2935 ! ----------------------------------------------------------------------
2936 IMPLICIT NONE
2937
2938 INTEGER, INTENT(IN) :: ICE, NROOT,NSOIL,VEGTYP
2939 INTEGER :: K
2940
2941 REAL, INTENT(IN) :: BEXP,CFACTR, CMCMAX,CSOIL,DKSAT,DT,DWSAT, &
2942 EPSCA,ETP,FDOWN,F1,FXEXP,FRZFACT,KDT,PC, &
2943 PRCP,PSISAT,Q2,QUARTZ,RCH,RR,SBETA,SFCTMP,&
2944 SHDFAC,SLOPE,SMCDRY,SMCMAX,SMCREF,SMCWLT, &
2945 T24,TBOT,TH2,ZBOT,EMISSI
2946 REAL, INTENT(INOUT) :: CMC,BETA,T1
2947 REAL, INTENT(OUT) :: DEW,DRIP,EC,EDIR,ETA,ETT,FLX1,FLX3, &
2948 RUNOFF1,RUNOFF2,RUNOFF3,SSOIL
2949 REAL, DIMENSION(1:NSOIL),INTENT(IN) :: RTDIS,ZSOIL
2950 REAL, DIMENSION(1:NSOIL),INTENT(OUT) :: ET
2951 REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: SMC,SH2O,STC
2952 REAL, DIMENSION(1:NSOIL) :: ET1
2953 REAL :: EC1,EDIR1,ETT1,DF1,ETA1,ETP1,PRCP1,YY, &
2954 YYNUM,ZZ1
2955
2956 ! ----------------------------------------------------------------------
2957 ! EXECUTABLE CODE BEGINS HERE:
2958 ! CONVERT ETP Fnd PRCP FROM KG M-2 S-1 TO M S-1 AND INITIALIZE DEW.
2959 ! ----------------------------------------------------------------------
2960 PRCP1 = PRCP * 0.001
2961 ETP1 = ETP * 0.001
2962 DEW = 0.0
2963 ! ----------------------------------------------------------------------
2964 ! INITIALIZE EVAP TERMS.
2965 ! ----------------------------------------------------------------------
2966 EDIR = 0.
2967 EDIR1 = 0.
2968 EC1 = 0.
2969 EC = 0.
2970 DO K = 1,NSOIL
2971 ET(K) = 0.
2972 ET1(K) = 0.
2973 END DO
2974 ETT = 0.
2975 ETT1 = 0.
2976
2977 IF (ETP > 0.0) THEN
2978 CALL EVAPO (ETA1,SMC,NSOIL,CMC,ETP1,DT,ZSOIL, &
2979 SH2O, &
2980 SMCMAX,BEXP,PC,SMCWLT,DKSAT,DWSAT, &
2981 SMCREF,SHDFAC,CMCMAX, &
2982 SMCDRY,CFACTR, &
2983 EDIR1,EC1,ET1,ETT1,SFCTMP,Q2,NROOT,RTDIS,FXEXP)
2984 CALL SMFLX (SMC,NSOIL,CMC,DT,PRCP1,ZSOIL, &
2985 SH2O,SLOPE,KDT,FRZFACT, &
2986 SMCMAX,BEXP,SMCWLT,DKSAT,DWSAT, &
2987 SHDFAC,CMCMAX, &
2988 RUNOFF1,RUNOFF2,RUNOFF3, &
2989 EDIR1,EC1,ET1, &
2990 DRIP)
2991
2992 ! ----------------------------------------------------------------------
2993 ! CONVERT MODELED EVAPOTRANSPIRATION FROM M S-1 TO KG M-2 S-1.
2994 ! ----------------------------------------------------------------------
2995
2996 ETA = ETA1 * 1000.0
2997
2998 ! ----------------------------------------------------------------------
2999 ! IF ETP < 0, ASSUME DEW FORMS (TRANSFORM ETP1 INTO DEW AND REINITIALIZE
3000 ! ETP1 TO ZERO).
3001 ! ----------------------------------------------------------------------
3002 ELSE
3003 DEW = - ETP1
3004
3005 ! ----------------------------------------------------------------------
3006 ! CONVERT PRCP FROM 'KG M-2 S-1' TO 'M S-1' AND ADD DEW AMOUNT.
3007 ! ----------------------------------------------------------------------
3008 ETP1 = 0.0
3009
3010 PRCP1 = PRCP1+ DEW
3011 CALL EVAPO (ETA1,SMC,NSOIL,CMC,ETP1,DT,ZSOIL, &
3012 SH2O, &
3013 SMCMAX,BEXP,PC,SMCWLT,DKSAT,DWSAT, &
3014 SMCREF,SHDFAC,CMCMAX, &
3015 SMCDRY,CFACTR, &
3016 EDIR1,EC1,ET1,ETT1,SFCTMP,Q2,NROOT,RTDIS,FXEXP)
3017 CALL SMFLX (SMC,NSOIL,CMC,DT,PRCP1,ZSOIL, &
3018 SH2O,SLOPE,KDT,FRZFACT, &
3019 SMCMAX,BEXP,SMCWLT,DKSAT,DWSAT, &
3020 SHDFAC,CMCMAX, &
3021 RUNOFF1,RUNOFF2,RUNOFF3, &
3022 EDIR1,EC1,ET1, &
3023 DRIP)
3024
3025 ! ----------------------------------------------------------------------
3026 ! CONVERT MODELED EVAPOTRANSPIRATION FROM 'M S-1' TO 'KG M-2 S-1'.
3027 ! ----------------------------------------------------------------------
3028 ETA = ETA1 * 1000.0
3029 END IF
3030
3031 ! ----------------------------------------------------------------------
3032 ! BASED ON ETP AND E VALUES, DETERMINE BETA
3033 ! ----------------------------------------------------------------------
3034
3035 IF ( ETP <= 0.0 ) THEN
3036 BETA = 0.0
3037 IF ( ETP < 0.0 ) THEN
3038 BETA = 1.0
3039 ETA = ETP
3040 END IF
3041 ELSE
3042 BETA = ETA / ETP
3043 END IF
3044
3045 ! ----------------------------------------------------------------------
3046 ! CONVERT MODELED EVAPOTRANSPIRATION COMPONENTS 'M S-1' TO 'KG M-2 S-1'.
3047 ! ----------------------------------------------------------------------
3048 EDIR = EDIR1*1000.
3049 EC = EC1*1000.
3050 DO K = 1,NSOIL
3051 ET(K) = ET1(K)*1000.
3052 END DO
3053 ETT = ETT1*1000.
3054
3055 ! ----------------------------------------------------------------------
3056 ! GET SOIL THERMAL DIFFUXIVITY/CONDUCTIVITY FOR TOP SOIL LYR,
3057 ! CALC. ADJUSTED TOP LYR SOIL TEMP AND ADJUSTED SOIL FLUX, THEN
3058 ! CALL SHFLX TO COMPUTE/UPDATE SOIL HEAT FLUX AND SOIL TEMPS.
3059 ! ----------------------------------------------------------------------
3060
3061 CALL TDFCND (DF1,SMC (1),QUARTZ,SMCMAX,SH2O (1))
3062
3063 !urban
3064 IF (VEGTYP==1) DF1=3.24
3065 !
3066
3067 ! ----------------------------------------------------------------------
3068 ! VEGETATION GREENNESS FRACTION REDUCTION IN SUBSURFACE HEAT FLUX
3069 ! VIA REDUCTION FACTOR, WHICH IS CONVENIENT TO APPLY HERE TO THERMAL
3070 ! DIFFUSIVITY THAT IS LATER USED IN HRT TO COMPUTE SUB SFC HEAT FLUX
3071 ! (SEE ADDITIONAL COMMENTS ON VEG EFFECT SUB-SFC HEAT FLX IN
3072 ! ROUTINE SFLX)
3073 ! ----------------------------------------------------------------------
3074 DF1 = DF1 * EXP (SBETA * SHDFAC)
3075 ! ----------------------------------------------------------------------
3076 ! COMPUTE INTERMEDIATE TERMS PASSED TO ROUTINE HRT (VIA ROUTINE
3077 ! SHFLX BELOW) FOR USE IN COMPUTING SUBSURFACE HEAT FLUX IN HRT
3078 ! ----------------------------------------------------------------------
3079 ! YYNUM = FDOWN - SIGMA * T24
3080 YYNUM = FDOWN - EMISSI*SIGMA * T24
3081 YY = SFCTMP + (YYNUM / RCH + TH2- SFCTMP - BETA * EPSCA) / RR
3082
3083 ZZ1 = DF1 / ( -0.5 * ZSOIL (1) * RCH * RR ) + 1.0
3084
3085 !urban
3086 CALL SHFLX (SSOIL,STC,SMC,SMCMAX,NSOIL,T1,DT,YY,ZZ1,ZSOIL, &
3087 TBOT,ZBOT,SMCWLT,PSISAT,SH2O,BEXP,F1,DF1,ICE, &
3088 QUARTZ,CSOIL,VEGTYP)
3089
3090 ! ----------------------------------------------------------------------
3091 ! SET FLX1 AND FLX3 (SNOPACK PHASE CHANGE HEAT FLUXES) TO ZERO SINCE
3092 ! THEY ARE NOT USED HERE IN SNOPAC. FLX2 (FREEZING RAIN HEAT FLUX) WAS
3093 ! SIMILARLY INITIALIZED IN THE PENMAN ROUTINE.
3094 ! ----------------------------------------------------------------------
3095 FLX1 = 0.0
3096 FLX3 = 0.0
3097
3098 ! ----------------------------------------------------------------------
3099 END SUBROUTINE NOPAC
3100 ! ----------------------------------------------------------------------
3101
3102 SUBROUTINE PENMAN (SFCTMP,SFCPRS,CH,T2V,TH2,PRCP,FDOWN,T24,SSOIL, &
3103 & Q2,Q2SAT,ETP,RCH,EPSCA,RR,SNOWNG,FRZGRA, &
3104 & DQSDT2,FLX2,EMISSI)
3105
3106 ! ----------------------------------------------------------------------
3107 ! SUBROUTINE PENMAN
3108 ! ----------------------------------------------------------------------
3109 ! CALCULATE POTENTIAL EVAPORATION FOR THE CURRENT POINT. VARIOUS
3110 ! PARTIAL SUMS/PRODUCTS ARE ALSO CALCULATED AND PASSED BACK TO THE
3111 ! CALLING ROUTINE FOR LATER USE.
3112 ! ----------------------------------------------------------------------
3113 IMPLICIT NONE
3114 LOGICAL, INTENT(IN) :: SNOWNG, FRZGRA
3115 REAL, INTENT(IN) :: CH, DQSDT2,FDOWN,PRCP, &
3116 Q2, Q2SAT,SSOIL, SFCPRS, SFCTMP, &
3117 T2V, TH2,EMISSI
3118 REAL, INTENT(OUT) :: EPSCA,ETP,FLX2,RCH,RR,T24
3119 REAL :: A, DELTA, FNET,RAD,RHO
3120
3121 REAL, PARAMETER :: ELCP = 2.4888E+3, LSUBC = 2.501000E+6,CP = 1004.6
3122
3123 ! ----------------------------------------------------------------------
3124 ! EXECUTABLE CODE BEGINS HERE:
3125 ! ----------------------------------------------------------------------
3126 ! ----------------------------------------------------------------------
3127 ! PREPARE PARTIAL QUANTITIES FOR PENMAN EQUATION.
3128 ! ----------------------------------------------------------------------
3129 FLX2 = 0.0
3130 DELTA = ELCP * DQSDT2
3131 T24 = SFCTMP * SFCTMP * SFCTMP * SFCTMP
3132 ! RR = T24 * 6.48E-8 / (SFCPRS * CH) + 1.0
3133 RR = EMISSI*T24 * 6.48E-8 / (SFCPRS * CH) + 1.0
3134 RHO = SFCPRS / (RD * T2V)
3135
3136 ! ----------------------------------------------------------------------
3137 ! ADJUST THE PARTIAL SUMS / PRODUCTS WITH THE LATENT HEAT
3138 ! EFFECTS CAUSED BY FALLING PRECIPITATION.
3139 ! ----------------------------------------------------------------------
3140 RCH = RHO * CP * CH
3141 IF (.NOT. SNOWNG) THEN
3142 IF (PRCP > 0.0) RR = RR + CPH2O * PRCP / RCH
3143 ELSE
3144 RR = RR + CPICE * PRCP / RCH
3145 END IF
3146
3147 ! ----------------------------------------------------------------------
3148 ! INCLUDE THE LATENT HEAT EFFECTS OF FRZNG RAIN CONVERTING TO ICE ON
3149 ! IMPACT IN THE CALCULATION OF FLX2 AND FNET.
3150 ! ----------------------------------------------------------------------
3151 ! FNET = FDOWN - SIGMA * T24- SSOIL
3152 FNET = FDOWN - EMISSI*SIGMA * T24- SSOIL
3153 IF (FRZGRA) THEN
3154 FLX2 = - LSUBF * PRCP
3155 FNET = FNET - FLX2
3156 ! ----------------------------------------------------------------------
3157 ! FINISH PENMAN EQUATION CALCULATIONS.
3158 ! ----------------------------------------------------------------------
3159 END IF
3160 RAD = FNET / RCH + TH2- SFCTMP
3161 A = ELCP * (Q2SAT - Q2)
3162 EPSCA = (A * RR + RAD * DELTA) / (DELTA + RR)
3163 ETP = EPSCA * RCH / LSUBC
3164
3165 ! ----------------------------------------------------------------------
3166 END SUBROUTINE PENMAN
3167 ! ----------------------------------------------------------------------
3168
3169 SUBROUTINE REDPRM (VEGTYP,SOILTYP,SLOPETYP,CFACTR,CMCMAX,RSMAX, &
3170 TOPT, &
3171 REFKDT,KDT,SBETA, SHDFAC,RSMIN,RGL,HS,ZBOT,FRZX, &
3172 PSISAT,SLOPE,SNUP,SALP,BEXP,DKSAT,DWSAT, &
3173 SMCMAX,SMCWLT,SMCREF,SMCDRY,F1,QUARTZ,FXEXP, &
3174 RTDIS,SLDPTH,ZSOIL, NROOT,NSOIL,Z0BRD,CZIL,LAI, &
3175 CSOIL,ALBBRD,PTU,LLANDUSE,LSOIL,LOCAL)
3176
3177 IMPLICIT NONE
3178 ! ----------------------------------------------------------------------
3179 ! Internally set (default valuess)
3180 ! all soil and vegetation parameters required for the execusion oF
3181 ! the Noah lsm are defined in VEGPARM.TBL, SOILPARM.TB, and GENPARM.TBL.
3182 ! ----------------------------------------------------------------------
3183 ! Vegetation parameters:
3184 ! ALBBRD: SFC background snow-free albedo
3185 ! CMXTBL: MAX CNPY Capacity
3186 ! Z0BRD: Background roughness length
3187 ! SHDFAC: Green vegetation fraction
3188 ! NROOT: Rooting depth
3189 ! RSMIN: Mimimum stomatal resistance
3190 ! RSMAX: Max. stomatal resistance
3191 ! RGL: Parameters used in radiation stress function
3192 ! HS: Parameter used in vapor pressure deficit functio
3193 ! TOPT: Optimum transpiration air temperature.
3194 ! CMCMAX: Maximum canopy water capacity
3195 ! CFACTR: Parameter used in the canopy inteception calculation
3196 ! SNUP: Threshold snow depth (in water equivalent m) that
3197 ! implies 100 percent snow cover
3198 ! LAI: Leaf area index
3199 !
3200 ! ----------------------------------------------------------------------
3201 ! Soil parameters:
3202 ! SMCMAX: MAX soil moisture content (porosity)
3203 ! SMCREF: Reference soil moisture (field capacity)
3204 ! SMCWLT: Wilting point soil moisture
3205 ! SMCWLT: Air dry soil moist content limits
3206 ! SSATPSI: SAT (saturation) soil potential
3207 ! DKSAT: SAT soil conductivity
3208 ! BEXP: B parameter
3209 ! SSATDW: SAT soil diffusivity
3210 ! F1: Soil thermal diffusivity/conductivity coef.
3211 ! QUARTZ: Soil quartz content
3212 ! Modified by F. Chen (12/22/97) to use the STATSGO soil map
3213 ! Modified By F. Chen (01/22/00) to include PLaya, Lava, and White San
3214 ! Modified By F. Chen (08/05/02) to include additional parameters for the Noah
3215 ! NOTE: SATDW = BB*SATDK*(SATPSI/MAXSMC)
3216 ! F11 = ALOG10(SATPSI) + BB*ALOG10(MAXSMC) + 2.0
3217 ! REFSMC1=MAXSMC*(5.79E-9/SATDK)**(1/(2*BB+3)) 5.79E-9 m/s= 0.5 mm
3218 ! REFSMC=REFSMC1+1./3.(MAXSMC-REFSMC1)
3219 ! WLTSMC1=MAXSMC*(200./SATPSI)**(-1./BB) (Wetzel and Chang, 198
3220 ! WLTSMC=WLTSMC1-0.5*WLTSMC1
3221 ! Note: the values for playa is set for it to have a thermal conductivit
3222 ! as sand and to have a hydrulic conductivity as clay
3223 !
3224 ! ----------------------------------------------------------------------
3225 ! Class parameter 'SLOPETYP' was included to estimate linear reservoir
3226 ! coefficient 'SLOPE' to the baseflow runoff out of the bottom layer.
3227 ! lowest class (slopetyp=0) means highest slope parameter = 1.
3228 ! definition of slopetyp from 'zobler' slope type:
3229 ! slope class percent slope
3230 ! 1 0-8
3231 ! 2 8-30
3232 ! 3 > 30
3233 ! 4 0-30
3234 ! 5 0-8 & > 30
3235 ! 6 8-30 & > 30
3236 ! 7 0-8, 8-30, > 30
3237 ! 9 GLACIAL ICE
3238 ! BLANK OCEAN/SEA
3239 ! SLOPE_DATA: linear reservoir coefficient
3240 ! SBETA_DATA: parameter used to caluculate vegetation effect on soil heat
3241 ! FXEXP_DAT: soil evaporation exponent used in DEVAP
3242 ! CSOIL_DATA: soil heat capacity [J M-3 K-1]
3243 ! SALP_DATA: shape parameter of distribution function of snow cover
3244 ! REFDK_DATA and REFKDT_DATA: parameters in the surface runoff parameteriz
3245 ! FRZK_DATA: frozen ground parameter
3246 ! ZBOT_DATA: depth[M] of lower boundary soil temperature
3247 ! CZIL_DATA: calculate roughness length of heat
3248 ! SMLOW_DATA and MHIGH_DATA: two soil moisture wilt, soil moisture referen
3249 ! parameters
3250 ! Set maximum number of soil-, veg-, and slopetyp in data statement.
3251 ! ----------------------------------------------------------------------
3252 INTEGER, PARAMETER :: MAX_SLOPETYP=30,MAX_SOILTYP=30,MAX_VEGTYP=30
3253 LOGICAL :: LOCAL
3254 CHARACTER (LEN=4), INTENT(IN):: LLANDUSE, LSOIL
3255
3256 ! Veg parameters
3257 INTEGER, INTENT(IN) :: VEGTYP
3258 INTEGER, INTENT(OUT) :: NROOT
3259 REAL, INTENT(OUT) :: HS,LAI,RSMIN,RGL,SHDFAC,SNUP,Z0BRD, &
3260 CMCMAX,RSMAX,TOPT,ALBBRD
3261 ! Soil parameters
3262 INTEGER, INTENT(IN) :: SOILTYP
3263 REAL, INTENT(OUT) :: BEXP,DKSAT,DWSAT,F1,QUARTZ,SMCDRY, &
3264 SMCMAX,SMCREF,SMCWLT,PSISAT
3265 ! General parameters
3266 INTEGER, INTENT(IN) :: SLOPETYP,NSOIL
3267 INTEGER :: I
3268
3269 REAL, INTENT(OUT) :: SLOPE,CZIL,SBETA,FXEXP, &
3270 CSOIL,SALP,FRZX,KDT,CFACTR, &
3271 ZBOT,REFKDT,PTU
3272 REAL,DIMENSION(1:NSOIL),INTENT(IN) :: SLDPTH,ZSOIL
3273 REAL,DIMENSION(1:NSOIL),INTENT(OUT):: RTDIS
3274 REAL :: FRZFACT,FRZK,REFDK
3275
3276 ! SAVE
3277 ! ----------------------------------------------------------------------
3278 !
3279 IF (SOILTYP .gt. SLCATS) THEN
3280 CALL wrf_error_fatal ( 'Warning: too many input soil types' )
3281 END IF
3282 IF (VEGTYP .gt. LUCATS) THEN
3283 CALL wrf_error_fatal ( 'Warning: too many input landuse types' )
3284 END IF
3285 IF (SLOPETYP .gt. SLPCATS) THEN
3286 CALL wrf_error_fatal ( 'Warning: too many input slope types' )
3287 END IF
3288
3289 ! ----------------------------------------------------------------------
3290 ! SET-UP SOIL PARAMETERS
3291 ! ----------------------------------------------------------------------
3292 CSOIL = CSOIL_DATA
3293 BEXP = BB (SOILTYP)
3294 DKSAT = SATDK (SOILTYP)
3295 DWSAT = SATDW (SOILTYP)
3296 F1 = F11 (SOILTYP)
3297 PSISAT = SATPSI (SOILTYP)
3298 QUARTZ = QTZ (SOILTYP)
3299 SMCDRY = DRYSMC (SOILTYP)
3300 SMCMAX = MAXSMC (SOILTYP)
3301 SMCREF = REFSMC (SOILTYP)
3302 SMCWLT = WLTSMC (SOILTYP)
3303 ! ----------------------------------------------------------------------
3304 ! Set-up universal parameters (not dependent on SOILTYP, VEGTYP or
3305 ! SLOPETYP)
3306 ! ----------------------------------------------------------------------
3307 ZBOT = ZBOT_DATA
3308 SALP = SALP_DATA
3309 SBETA = SBETA_DATA
3310 REFDK = REFDK_DATA
3311 FRZK = FRZK_DATA
3312 FXEXP = FXEXP_DATA
3313 REFKDT = REFKDT_DATA
3314 PTU = 0. ! (not used yet) to satisify intent(out)
3315 KDT = REFKDT * DKSAT / REFDK
3316 CZIL = CZIL_DATA
3317 SLOPE = SLOPE_DATA (SLOPETYP)
3318
3319 ! ----------------------------------------------------------------------
3320 ! TO ADJUST FRZK PARAMETER TO ACTUAL SOIL TYPE: FRZK * FRZFACT
3321 ! ----------------------------------------------------------------------
3322 FRZFACT = (SMCMAX / SMCREF) * (0.412 / 0.468)
3323 FRZX = FRZK * FRZFACT
3324
3325 ! ----------------------------------------------------------------------
3326 ! SET-UP VEGETATION PARAMETERS
3327 ! ----------------------------------------------------------------------
3328 TOPT = TOPT_DATA
3329 CMCMAX = CMCMAX_DATA
3330 CFACTR = CFACTR_DATA
3331 RSMAX = RSMAX_DATA
3332 NROOT = NROTBL (VEGTYP)
3333 SNUP = SNUPTBL (VEGTYP)
3334 RSMIN = RSTBL (VEGTYP)
3335 RGL = RGLTBL (VEGTYP)
3336 HS = HSTBL (VEGTYP)
3337 LAI = LAITBL (VEGTYP)
3338 IF(LOCAL) THEN
3339 ALBBRD = ALBTBL(VEGTYP)
3340 Z0BRD = Z0TBL(VEGTYP)
3341 SHDFAC = SHDTBL(VEGTYP)
3342 ENDIF
3343
3344 IF (VEGTYP .eq. BARE) SHDFAC = 0.0
3345 IF (NROOT .gt. NSOIL) THEN
3346 WRITE (err_message,*) 'Error: too many root layers ', &
3347 NSOIL,NROOT
3348 CALL wrf_error_fatal ( err_message )
3349 ! ----------------------------------------------------------------------
3350 ! CALCULATE ROOT DISTRIBUTION. PRESENT VERSION ASSUMES UNIFORM
3351 ! DISTRIBUTION BASED ON SOIL LAYER DEPTHS.
3352 ! ----------------------------------------------------------------------
3353 END IF
3354 DO I = 1,NROOT
3355 RTDIS (I) = - SLDPTH (I)/ ZSOIL (NROOT)
3356 ! ----------------------------------------------------------------------
3357 ! SET-UP SLOPE PARAMETER
3358 ! ----------------------------------------------------------------------
3359 END DO
3360
3361 ! print*,'end of PRMRED'
3362 ! print*,'VEGTYP',VEGTYP,'SOILTYP',SOILTYP,'SLOPETYP',SLOPETYP, &
3363 ! & 'CFACTR',CFACTR,'CMCMAX',CMCMAX,'RSMAX',RSMAX,'TOPT',TOPT, &
3364 ! & 'REFKDT',REFKDT,'KDT',KDT,'SBETA',SBETA, 'SHDFAC',SHDFAC, &
3365 ! & 'RSMIN',RSMIN,'RGL',RGL,'HS',HS,'ZBOT',ZBOT,'FRZX',FRZX, &
3366 ! & 'PSISAT',PSISAT,'SLOPE',SLOPE,'SNUP',SNUP,'SALP',SALP,'BEXP', &
3367 ! & BEXP, &
3368 ! & 'DKSAT',DKSAT,'DWSAT',DWSAT, &
3369 ! & 'SMCMAX',SMCMAX,'SMCWLT',SMCWLT,'SMCREF',SMCREF,'SMCDRY',SMCDRY, &
3370 ! & 'F1',F1,'QUARTZ',QUARTZ,'FXEXP',FXEXP, &
3371 ! & 'RTDIS',RTDIS,'SLDPTH',SLDPTH,'ZSOIL',ZSOIL, 'NROOT',NROOT, &
3372 ! & 'NSOIL',NSOIL,'Z0',Z0,'CZIL',CZIL,'LAI',LAI, &
3373 ! & 'CSOIL',CSOIL,'PTU',PTU, &
3374 ! & 'LOCAL', LOCAL
3375
3376 END SUBROUTINE REDPRM
3377
3378 SUBROUTINE ROSR12 (P,A,B,C,D,DELTA,NSOIL)
3379
3380 ! ----------------------------------------------------------------------
3381 ! SUBROUTINE ROSR12
3382 ! ----------------------------------------------------------------------
3383 ! INVERT (SOLVE) THE TRI-DIAGONAL MATRIX PROBLEM SHOWN BELOW:
3384 ! ### ### ### ### ### ###
3385 ! #B(1), C(1), 0 , 0 , 0 , . . . , 0 # # # # #
3386 ! #A(2), B(2), C(2), 0 , 0 , . . . , 0 # # # # #
3387 ! # 0 , A(3), B(3), C(3), 0 , . . . , 0 # # # # D(3) #
3388 ! # 0 , 0 , A(4), B(4), C(4), . . . , 0 # # P(4) # # D(4) #
3389 ! # 0 , 0 , 0 , A(5), B(5), . . . , 0 # # P(5) # # D(5) #
3390 ! # . . # # . # = # . #
3391 ! # . . # # . # # . #
3392 ! # . . # # . # # . #
3393 ! # 0 , . . . , 0 , A(M-2), B(M-2), C(M-2), 0 # #P(M-2)# #D(M-2)#
3394 ! # 0 , . . . , 0 , 0 , A(M-1), B(M-1), C(M-1)# #P(M-1)# #D(M-1)#
3395 ! # 0 , . . . , 0 , 0 , 0 , A(M) , B(M) # # P(M) # # D(M) #
3396 ! ### ### ### ### ### ###
3397 ! ----------------------------------------------------------------------
3398 IMPLICIT NONE
3399
3400 INTEGER, INTENT(IN) :: NSOIL
3401 INTEGER :: K, KK
3402
3403 REAL, DIMENSION(1:NSOIL), INTENT(IN):: A, B, D
3404 REAL, DIMENSION(1:NSOIL),INTENT(INOUT):: C,P,DELTA
3405
3406 ! ----------------------------------------------------------------------
3407 ! INITIALIZE EQN COEF C FOR THE LOWEST SOIL LAYER
3408 ! ----------------------------------------------------------------------
3409 C (NSOIL) = 0.0
3410 P (1) = - C (1) / B (1)
3411 ! ----------------------------------------------------------------------
3412 ! SOLVE THE COEFS FOR THE 1ST SOIL LAYER
3413 ! ----------------------------------------------------------------------
3414
3415 ! ----------------------------------------------------------------------
3416 ! SOLVE THE COEFS FOR SOIL LAYERS 2 THRU NSOIL
3417 ! ----------------------------------------------------------------------
3418 DELTA (1) = D (1) / B (1)
3419 DO K = 2,NSOIL
3420 P (K) = - C (K) * ( 1.0 / (B (K) + A (K) * P (K -1)) )
3421 DELTA (K) = (D (K) - A (K)* DELTA (K -1))* (1.0/ (B (K) + A (K)&
3422 * P (K -1)))
3423 END DO
3424 ! ----------------------------------------------------------------------
3425 ! SET P TO DELTA FOR LOWEST SOIL LAYER
3426 ! ----------------------------------------------------------------------
3427 P (NSOIL) = DELTA (NSOIL)
3428
3429 ! ----------------------------------------------------------------------
3430 ! ADJUST P FOR SOIL LAYERS 2 THRU NSOIL
3431 ! ----------------------------------------------------------------------
3432 DO K = 2,NSOIL
3433 KK = NSOIL - K + 1
3434 P (KK) = P (KK) * P (KK +1) + DELTA (KK)
3435 END DO
3436 ! ----------------------------------------------------------------------
3437 END SUBROUTINE ROSR12
3438 ! ----------------------------------------------------------------------
3439
3440
3441 SUBROUTINE SHFLX (SSOIL,STC,SMC,SMCMAX,NSOIL,T1,DT,YY,ZZ1,ZSOIL, &
3442 TBOT,ZBOT,SMCWLT,PSISAT,SH2O,BEXP,F1,DF1,ICE, &
3443 QUARTZ,CSOIL,VEGTYP)
3444
3445 ! ----------------------------------------------------------------------
3446 ! SUBROUTINE SHFLX
3447 ! ----------------------------------------------------------------------
3448 ! UPDATE THE TEMPERATURE STATE OF THE SOIL COLUMN BASED ON THE THERMAL
3449 ! DIFFUSION EQUATION AND UPDATE THE FROZEN SOIL MOISTURE CONTENT BASED
3450 ! ON THE TEMPERATURE.
3451 ! ----------------------------------------------------------------------
3452 IMPLICIT NONE
3453
3454 INTEGER, INTENT(IN) :: ICE, NSOIL, VEGTYP
3455 INTEGER :: I
3456
3457 REAL, INTENT(IN) :: BEXP,CSOIL,DF1,DT,F1,PSISAT,QUARTZ, &
3458 SMCMAX, SMCWLT, TBOT,YY, ZBOT,ZZ1
3459 REAL, INTENT(INOUT) :: T1
3460 REAL, INTENT(OUT) :: SSOIL
3461 REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SMC,ZSOIL
3462 REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: SH2O
3463 REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: STC
3464 REAL, DIMENSION(1:NSOIL) :: AI, BI, CI, STCF,RHSTS
3465 REAL, PARAMETER :: T0 = 273.15
3466
3467 ! ----------------------------------------------------------------------
3468 ! HRT ROUTINE CALCS THE RIGHT HAND SIDE OF THE SOIL TEMP DIF EQN
3469 ! ----------------------------------------------------------------------
3470
3471 ! ----------------------------------------------------------------------
3472 ! SEA-ICE CASE
3473 ! ----------------------------------------------------------------------
3474 IF (ICE == 1) THEN
3475
3476 CALL HRTICE (RHSTS,STC,NSOIL,ZSOIL,YY,ZZ1,DF1,AI,BI,CI)
3477
3478 CALL HSTEP (STCF,STC,RHSTS,DT,NSOIL,AI,BI,CI)
3479
3480 ! ----------------------------------------------------------------------
3481 ! LAND-MASS CASE
3482 ! ----------------------------------------------------------------------
3483 ELSE
3484 CALL HRT (RHSTS,STC,SMC,SMCMAX,NSOIL,ZSOIL,YY,ZZ1,TBOT, &
3485 ZBOT,PSISAT,SH2O,DT, &
3486 BEXP,F1,DF1,QUARTZ,CSOIL,AI,BI,CI,VEGTYP)
3487
3488 CALL HSTEP (STCF,STC,RHSTS,DT,NSOIL,AI,BI,CI)
3489 END IF
3490 DO I = 1,NSOIL
3491 STC (I) = STCF (I)
3492 END DO
3493
3494 ! ----------------------------------------------------------------------
3495 ! IN THE NO SNOWPACK CASE (VIA ROUTINE NOPAC BRANCH,) UPDATE THE GRND
3496 ! (SKIN) TEMPERATURE HERE IN RESPONSE TO THE UPDATED SOIL TEMPERATURE
3497 ! PROFILE ABOVE. (NOTE: INSPECTION OF ROUTINE SNOPAC SHOWS THAT T1
3498 ! BELOW IS A DUMMY VARIABLE ONLY, AS SKIN TEMPERATURE IS UPDATED
3499 ! DIFFERENTLY IN ROUTINE SNOPAC)
3500 ! ----------------------------------------------------------------------
3501 ! ----------------------------------------------------------------------
3502 ! CALCULATE SURFACE SOIL HEAT FLUX
3503 ! ----------------------------------------------------------------------
3504 T1 = (YY + (ZZ1- 1.0) * STC (1)) / ZZ1
3505 SSOIL = DF1 * (STC (1) - T1) / (0.5 * ZSOIL (1))
3506
3507 ! ----------------------------------------------------------------------
3508 END SUBROUTINE SHFLX
3509 ! ----------------------------------------------------------------------
3510
3511 SUBROUTINE SMFLX (SMC,NSOIL,CMC,DT,PRCP1,ZSOIL, &
3512 & SH2O,SLOPE,KDT,FRZFACT, &
3513 & SMCMAX,BEXP,SMCWLT,DKSAT,DWSAT, &
3514 & SHDFAC,CMCMAX, &
3515 & RUNOFF1,RUNOFF2,RUNOFF3, &
3516 & EDIR,EC,ET, &
3517 & DRIP)
3518
3519 ! ----------------------------------------------------------------------
3520 ! SUBROUTINE SMFLX
3521 ! ----------------------------------------------------------------------
3522 ! CALCULATE SOIL MOISTURE FLUX. THE SOIL MOISTURE CONTENT (SMC - A PER
3523 ! UNIT VOLUME MEASUREMENT) IS A DEPENDENT VARIABLE THAT IS UPDATED WITH
3524 ! PROGNOSTIC EQNS. THE CANOPY MOISTURE CONTENT (CMC) IS ALSO UPDATED.
3525 ! FROZEN GROUND VERSION: NEW STATES ADDED: SH2O, AND FROZEN GROUND
3526 ! CORRECTION FACTOR, FRZFACT AND PARAMETER SLOPE.
3527 ! ----------------------------------------------------------------------
3528 IMPLICIT NONE
3529
3530 INTEGER, INTENT(IN) :: NSOIL
3531 INTEGER :: I,K
3532
3533 REAL, INTENT(IN) :: BEXP, CMCMAX, DKSAT,DWSAT, DT, EC, EDIR, &
3534 KDT, PRCP1, SHDFAC, SLOPE, SMCMAX, SMCWLT
3535 REAL, INTENT(OUT) :: DRIP, RUNOFF1, RUNOFF2, RUNOFF3
3536 REAL, INTENT(INOUT) :: CMC
3537 REAL, DIMENSION(1:NSOIL), INTENT(IN) :: ET,ZSOIL
3538 REAL, DIMENSION(1:NSOIL), INTENT(INOUT):: SMC, SH2O
3539 REAL, DIMENSION(1:NSOIL) :: AI, BI, CI, STCF,RHSTS, RHSTT, &
3540 SICE, SH2OA, SH2OFG
3541 REAL :: DUMMY, EXCESS,FRZFACT,PCPDRP,RHSCT,TRHSCT
3542
3543 ! ----------------------------------------------------------------------
3544 ! EXECUTABLE CODE BEGINS HERE.
3545 ! ----------------------------------------------------------------------
3546 ! ----------------------------------------------------------------------
3547 ! COMPUTE THE RIGHT HAND SIDE OF THE CANOPY EQN TERM ( RHSCT )
3548 ! ----------------------------------------------------------------------
3549 DUMMY = 0.
3550
3551 ! ----------------------------------------------------------------------
3552 ! CONVERT RHSCT (A RATE) TO TRHSCT (AN AMOUNT) AND ADD IT TO EXISTING
3553 ! CMC. IF RESULTING AMT EXCEEDS MAX CAPACITY, IT BECOMES DRIP AND WILL
3554 ! FALL TO THE GRND.
3555 ! ----------------------------------------------------------------------
3556 RHSCT = SHDFAC * PRCP1- EC
3557 DRIP = 0.
3558 TRHSCT = DT * RHSCT
3559 EXCESS = CMC + TRHSCT
3560
3561 ! ----------------------------------------------------------------------
3562 ! PCPDRP IS THE COMBINED PRCP1 AND DRIP (FROM CMC) THAT GOES INTO THE
3563 ! SOIL
3564 ! ----------------------------------------------------------------------
3565 IF (EXCESS > CMCMAX) DRIP = EXCESS - CMCMAX
3566 PCPDRP = (1. - SHDFAC) * PRCP1+ DRIP / DT
3567
3568 ! ----------------------------------------------------------------------
3569 ! STORE ICE CONTENT AT EACH SOIL LAYER BEFORE CALLING SRT and SSTEP
3570 !
3571 DO I = 1,NSOIL
3572 SICE (I) = SMC (I) - SH2O (I)
3573 END DO
3574 ! ----------------------------------------------------------------------
3575 ! CALL SUBROUTINES SRT AND SSTEP TO SOLVE THE SOIL MOISTURE
3576 ! TENDENCY EQUATIONS.
3577 ! IF THE INFILTRATING PRECIP RATE IS NONTRIVIAL,
3578 ! (WE CONSIDER NONTRIVIAL TO BE A PRECIP TOTAL OVER THE TIME STEP
3579 ! EXCEEDING ONE ONE-THOUSANDTH OF THE WATER HOLDING CAPACITY OF
3580 ! THE FIRST SOIL LAYER)
3581 ! THEN CALL THE SRT/SSTEP SUBROUTINE PAIR TWICE IN THE MANNER OF
3582 ! TIME SCHEME "F" (IMPLICIT STATE, AVERAGED COEFFICIENT)
3583 ! OF SECTION 2 OF KALNAY AND KANAMITSU (1988, MWR, VOL 116,
3584 ! PAGES 1945-1958)TO MINIMIZE 2-DELTA-T OSCILLATIONS IN THE
3585 ! SOIL MOISTURE VALUE OF THE TOP SOIL LAYER THAT CAN ARISE BECAUSE
3586 ! OF THE EXTREME NONLINEAR DEPENDENCE OF THE SOIL HYDRAULIC
3587 ! DIFFUSIVITY COEFFICIENT AND THE HYDRAULIC CONDUCTIVITY ON THE
3588 ! SOIL MOISTURE STATE
3589 ! OTHERWISE CALL THE SRT/SSTEP SUBROUTINE PAIR ONCE IN THE MANNER OF
3590 ! TIME SCHEME "D" (IMPLICIT STATE, EXPLICIT COEFFICIENT)
3591 ! OF SECTION 2 OF KALNAY AND KANAMITSU
3592 ! PCPDRP IS UNITS OF KG/M**2/S OR MM/S, ZSOIL IS NEGATIVE DEPTH IN M
3593 ! ----------------------------------------------------------------------
3594 ! IF ( PCPDRP .GT. 0.0 ) THEN
3595
3596 ! ----------------------------------------------------------------------
3597 ! FROZEN GROUND VERSION:
3598 ! SMC STATES REPLACED BY SH2O STATES IN SRT SUBR. SH2O & SICE STATES
3599 ! INC&UDED IN SSTEP SUBR. FROZEN GROUND CORRECTION FACTOR, FRZFACT
3600 ! ADDED. ALL WATER BALANCE CALCULATIONS USING UNFROZEN WATER
3601 ! ----------------------------------------------------------------------
3602 IF ( (PCPDRP * DT) > (0.001*1000.0* (- ZSOIL (1))* SMCMAX) ) THEN
3603 CALL SRT (RHSTT,EDIR,ET,SH2O,SH2O,NSOIL,PCPDRP,ZSOIL, &
3604 DWSAT,DKSAT,SMCMAX,BEXP,RUNOFF1, &
3605 RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZFACT,SICE,AI,BI,CI)
3606 CALL SSTEP (SH2OFG,SH2O,DUMMY,RHSTT,RHSCT,DT,NSOIL,SMCMAX, &
3607 CMCMAX,RUNOFF3,ZSOIL,SMC,SICE,AI,BI,CI)
3608 DO K = 1,NSOIL
3609 SH2OA (K) = (SH2O (K) + SH2OFG (K)) * 0.5
3610 END DO
3611 CALL SRT (RHSTT,EDIR,ET,SH2O,SH2OA,NSOIL,PCPDRP,ZSOIL, &
3612 DWSAT,DKSAT,SMCMAX,BEXP,RUNOFF1, &
3613 RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZFACT,SICE,AI,BI,CI)
3614 CALL SSTEP (SH2O,SH2O,CMC,RHSTT,RHSCT,DT,NSOIL,SMCMAX, &
3615 CMCMAX,RUNOFF3,ZSOIL,SMC,SICE,AI,BI,CI)
3616
3617 ELSE
3618 CALL SRT (RHSTT,EDIR,ET,SH2O,SH2O,NSOIL,PCPDRP,ZSOIL, &
3619 DWSAT,DKSAT,SMCMAX,BEXP,RUNOFF1, &
3620 RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZFACT,SICE,AI,BI,CI)
3621 CALL SSTEP (SH2O,SH2O,CMC,RHSTT,RHSCT,DT,NSOIL,SMCMAX, &
3622 CMCMAX,RUNOFF3,ZSOIL,SMC,SICE,AI,BI,CI)
3623 ! RUNOF = RUNOFF
3624
3625 END IF
3626
3627 ! ----------------------------------------------------------------------
3628 END SUBROUTINE SMFLX
3629 ! ----------------------------------------------------------------------
3630
3631
3632 SUBROUTINE SNFRAC (SNEQV,SNUP,SALP,SNOWH,SNCOVR)
3633
3634 ! ----------------------------------------------------------------------
3635 ! SUBROUTINE SNFRAC
3636 ! ----------------------------------------------------------------------
3637 ! CALCULATE SNOW FRACTION (0 -> 1)
3638 ! SNEQV SNOW WATER EQUIVALENT (M)
3639 ! SNUP THRESHOLD SNEQV DEPTH ABOVE WHICH SNCOVR=1
3640 ! SALP TUNING PARAMETER
3641 ! SNCOVR FRACTIONAL SNOW COVER
3642 ! ----------------------------------------------------------------------
3643 IMPLICIT NONE
3644
3645 REAL, INTENT(IN) :: SNEQV,SNUP,SALP,SNOWH
3646 REAL, INTENT(OUT) :: SNCOVR
3647 REAL :: RSNOW, Z0N
3648
3649 ! ----------------------------------------------------------------------
3650 ! SNUP IS VEG-CLASS DEPENDENT SNOWDEPTH THRESHHOLD (SET IN ROUTINE
3651 ! REDPRM) ABOVE WHICH SNOCVR=1.
3652 ! ----------------------------------------------------------------------
3653 IF (SNEQV < SNUP) THEN
3654 RSNOW = SNEQV / SNUP
3655 SNCOVR = 1. - ( EXP ( - SALP * RSNOW) - RSNOW * EXP ( - SALP))
3656 ELSE
3657 SNCOVR = 1.0
3658 END IF
3659
3660 ! FORMULATION OF DICKINSON ET AL. 1986
3661 ! Z0N = 0.035
3662
3663 ! SNCOVR=SNOWH/(SNOWH + 5*Z0N)
3664
3665 ! FORMULATION OF MARSHALL ET AL. 1994
3666 ! SNCOVR=SNEQV/(SNEQV + 2*Z0N)
3667
3668 ! ----------------------------------------------------------------------
3669 END SUBROUTINE SNFRAC
3670 ! ----------------------------------------------------------------------
3671
3672 SUBROUTINE SNKSRC (TSNSR,TAVG,SMC,SH2O,ZSOIL,NSOIL, &
3673 & SMCMAX,PSISAT,BEXP,DT,K,QTOT)
3674
3675 ! ----------------------------------------------------------------------
3676 ! SUBROUTINE SNKSRC
3677 ! ----------------------------------------------------------------------
3678 ! CALCULATE SINK/SOURCE TERM OF THE TERMAL DIFFUSION EQUATION. (SH2O) IS
3679 ! AVAILABLE LIQUED WATER.
3680 ! ----------------------------------------------------------------------
3681 IMPLICIT NONE
3682
3683 INTEGER, INTENT(IN) :: K,NSOIL
3684
3685 REAL, INTENT(IN) :: BEXP, DT, PSISAT, QTOT, SMC, SMCMAX, &
3686 TAVG
3687 REAL, INTENT(INOUT) :: SH2O
3688
3689 REAL, DIMENSION(1:NSOIL), INTENT(IN):: ZSOIL
3690
3691 REAL :: DF, DZ, DZH, FREE, TSNSR, &
3692 TDN, TM, TUP, TZ, X0, XDN, XH2O, XUP
3693
3694 REAL, PARAMETER :: DH2O = 1.0000E3, HLICE = 3.3350E5, &
3695 T0 = 2.7315E2
3696
3697 IF (K == 1) THEN
3698 DZ = - ZSOIL (1)
3699 ELSE
3700 DZ = ZSOIL (K -1) - ZSOIL (K)
3701 END IF
3702 ! ----------------------------------------------------------------------
3703 ! VIA FUNCTION FRH2O, COMPUTE POTENTIAL OR 'EQUILIBRIUM' UNFROZEN
3704 ! SUPERCOOLED FREE WATER FOR GIVEN SOIL TYPE AND SOIL LAYER TEMPERATURE.
3705 ! FUNCTION FRH20 INVOKES EQN (17) FROM V. KOREN ET AL (1999, JGR, VOL.
3706 ! 104, PG 19573). (ASIDE: LATTER EQN IN JOURNAL IN CENTIGRADE UNITS.
3707 ! ROUTINE FRH2O USE FORM OF EQN IN KELVIN UNITS.)
3708 ! ----------------------------------------------------------------------
3709 ! FREE = FRH2O(TAVG,SMC,SH2O,SMCMAX,BEXP,PSISAT)
3710
3711 ! ----------------------------------------------------------------------
3712 ! IN NEXT BLOCK OF CODE, INVOKE EQN 18 OF V. KOREN ET AL (1999, JGR,
3713 ! VOL. 104, PG 19573.) THAT IS, FIRST ESTIMATE THE NEW AMOUNTOF LIQUID
3714 ! WATER, 'XH2O', IMPLIED BY THE SUM OF (1) THE LIQUID WATER AT THE BEGIN
3715 ! OF CURRENT TIME STEP, AND (2) THE FREEZE OF THAW CHANGE IN LIQUID
3716 ! WATER IMPLIED BY THE HEAT FLUX 'QTOT' PASSED IN FROM ROUTINE HRT.
3717 ! SECOND, DETERMINE IF XH2O NEEDS TO BE BOUNDED BY 'FREE' (EQUIL AMT) OR
3718 ! IF 'FREE' NEEDS TO BE BOUNDED BY XH2O.
3719 ! ----------------------------------------------------------------------
3720 CALL FRH2O (FREE,TAVG,SMC,SH2O,SMCMAX,BEXP,PSISAT)
3721
3722 ! ----------------------------------------------------------------------
3723 ! FIRST, IF FREEZING AND REMAINING LIQUID LESS THAN LOWER BOUND, THEN
3724 ! REDUCE EXTENT OF FREEZING, THEREBY LETTING SOME OR ALL OF HEAT FLUX
3725 ! QTOT COOL THE SOIL TEMP LATER IN ROUTINE HRT.
3726 ! ----------------------------------------------------------------------
3727 XH2O = SH2O + QTOT * DT / (DH2O * HLICE * DZ)
3728 IF ( XH2O < SH2O .AND. XH2O < FREE) THEN
3729 IF ( FREE > SH2O ) THEN
3730 XH2O = SH2O
3731 ELSE
3732 XH2O = FREE
3733 END IF
3734 END IF
3735 ! ----------------------------------------------------------------------
3736 ! SECOND, IF THAWING AND THE INCREASE IN LIQUID WATER GREATER THAN UPPER
3737 ! BOUND, THEN REDUCE EXTENT OF THAW, THEREBY LETTING SOME OR ALL OF HEAT
3738 ! FLUX QTOT WARM THE SOIL TEMP LATER IN ROUTINE HRT.
3739 ! ----------------------------------------------------------------------
3740 IF ( XH2O > SH2O .AND. XH2O > FREE ) THEN
3741 IF ( FREE < SH2O ) THEN
3742 XH2O = SH2O
3743 ELSE
3744 XH2O = FREE
3745 END IF
3746 END IF
3747
3748 ! ----------------------------------------------------------------------
3749 ! CALCULATE PHASE-CHANGE HEAT SOURCE/SINK TERM FOR USE IN ROUTINE HRT
3750 ! AND UPDATE LIQUID WATER TO REFLCET FINAL FREEZE/THAW INCREMENT.
3751 ! ----------------------------------------------------------------------
3752 ! SNKSRC = -DH2O*HLICE*DZ*(XH2O-SH2O)/DT
3753 IF (XH2O < 0.) XH2O = 0.
3754 IF (XH2O > SMC) XH2O = SMC
3755 TSNSR = - DH2O * HLICE * DZ * (XH2O - SH2O)/ DT
3756 SH2O = XH2O
3757
3758 ! ----------------------------------------------------------------------
3759 END SUBROUTINE SNKSRC
3760 ! ----------------------------------------------------------------------
3761
3762 SUBROUTINE SNOPAC (ETP,ETA,PRCP,PRCPF,SNOWNG,SMC,SMCMAX,SMCWLT, &
3763 SMCREF,SMCDRY,CMC,CMCMAX,NSOIL,DT, &
3764 SBETA,DF1, &
3765 Q2,T1,SFCTMP,T24,TH2,FDOWN,F1,SSOIL,STC,EPSCA,&
3766 SFCPRS,BEXP,PC,RCH,RR,CFACTR,SNCOVR,ESD,SNDENS,&
3767 SNOWH,SH2O,SLOPE,KDT,FRZFACT,PSISAT, &
3768 ZSOIL,DWSAT,DKSAT,TBOT,ZBOT,SHDFAC,RUNOFF1, &
3769 RUNOFF2,RUNOFF3,EDIR,EC,ET,ETT,NROOT,SNOMLT, &
3770 ICE,RTDIS,QUARTZ,FXEXP,CSOIL, &
3771 BETA,DRIP,DEW,FLX1,FLX2,FLX3,ESNOW,ETNS,EMISSI,&
3772 VEGTYP)
3773
3774 ! ----------------------------------------------------------------------
3775 ! SUBROUTINE SNOPAC
3776 ! ----------------------------------------------------------------------
3777 ! CALCULATE SOIL MOISTURE AND HEAT FLUX VALUES & UPDATE SOIL MOISTURE
3778 ! CONTENT AND SOIL HEAT CONTENT VALUES FOR THE CASE WHEN A SNOW PACK IS
3779 ! PRESENT.
3780 ! ----------------------------------------------------------------------
3781 IMPLICIT NONE
3782
3783 INTEGER, INTENT(IN) :: ICE, NROOT, NSOIL,VEGTYP
3784 INTEGER :: K
3785 LOGICAL, INTENT(IN) :: SNOWNG
3786 REAL, INTENT(IN) :: BEXP,CFACTR, CMCMAX,CSOIL,DF1,DKSAT, &
3787 DT,DWSAT, EPSCA,ETP,FDOWN,F1,FXEXP, &
3788 FRZFACT,KDT,PC, PRCP,PSISAT,Q2,QUARTZ, &
3789 RCH,RR,SBETA,SFCPRS, SFCTMP, SHDFAC, &
3790 SLOPE,SMCDRY,SMCMAX,SMCREF,SMCWLT, T24, &
3791 TBOT,TH2,ZBOT,EMISSI
3792 REAL, INTENT(INOUT) :: CMC, BETA, ESD,FLX2,PRCPF,SNOWH,SNCOVR, &
3793 SNDENS, T1
3794 REAL, INTENT(OUT) :: DEW,DRIP,EC,EDIR, ETNS, ESNOW,ETT, &
3795 FLX1,FLX3, RUNOFF1,RUNOFF2,RUNOFF3, &
3796 SSOIL,SNOMLT
3797 REAL, DIMENSION(1:NSOIL),INTENT(IN) :: RTDIS,ZSOIL
3798 REAL, DIMENSION(1:NSOIL),INTENT(OUT) :: ET
3799 REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: SMC,SH2O,STC
3800 REAL, DIMENSION(1:NSOIL) :: ET1
3801 REAL :: DENOM,DSOIL,DTOT,EC1,EDIR1,ESDFLX,ETA, &
3802 ETT1, ESNOW1, ESNOW2, ETA1,ETP1,ETP2, &
3803 ETP3, ETNS1, ETANRG, ETAX, EX, FLX3X, &
3804 FRCSNO,FRCSOI, PRCP1, QSAT,RSNOW, SEH, &
3805 SNCOND,SSOIL1, T11,T12, T12A, T12AX, &
3806 T12B, T14, YY, ZZ1, EMISSI_S
3807
3808 REAL, PARAMETER :: ESDMIN = 1.E-6, LSUBC = 2.501000E+6, &
3809 LSUBS = 2.83E+6, TFREEZ = 273.15, &
3810 SNOEXP = 1.0
3811
3812 ! ----------------------------------------------------------------------
3813 ! EXECUTABLE CODE BEGINS HERE:
3814 ! INITIALIZE EVAP TERMS.
3815 ! ----------------------------------------------------------------------
3816 ! conversions:
3817 ! ESNOW [KG M-2 S-1]
3818 ! ESDFLX [KG M-2 S-1] .le. ESNOW
3819 ! ESNOW1 [M S-1]
3820 ! ESNOW2 [M]
3821 ! ETP [KG M-2 S-1]
3822 ! ETP1 [M S-1]
3823 ! ETP2 [M]
3824 ! ----------------------------------------------------------------------
3825 EDIR = 0.
3826 EDIR1 = 0.
3827 EC1 = 0.
3828 EC = 0.
3829 EMISSI_S=0.9 ! For snow
3830 DO K = 1,NSOIL
3831 ET (K) = 0.
3832 ET1 (K) = 0.
3833 END DO
3834 ETT = 0.
3835 ETT1 = 0.
3836 ETNS = 0.
3837 ETNS1 = 0.
3838 ESNOW = 0.
3839 ESNOW1 = 0.
3840 ESNOW2 = 0.
3841
3842 ! CONVERT POTENTIAL EVAP (ETP) FROM KG M-2 S-1 TO M S-1 AND THEN TO AN
3843 ! AMOUNT (M) GIVEN TIMESTEP (DT) AND CALL IT AN EFFECTIVE SNOWPACK
3844 ! REDUCTION AMOUNT, ETP2 (M). THIS IS THE AMOUNT THE SNOWPACK WOULD BE
3845 ! REDUCED DUE TO EVAPORATION FROM THE SNOW SFC DURING THE TIMESTEP.
3846 ! EVAPORATION WILL PROCEED AT THE POTENTIAL RATE UNLESS THE SNOW DEPTH
3847 ! IS LESS THAN THE EXPECTED SNOWPACK REDUCTION.
3848 ! IF SEAICE (ICE=1), BETA REMAINS=1.
3849 ! ----------------------------------------------------------------------
3850 ! PRCP1 = PRCP1*0.001
3851 ! change name of PRCP1 to PRCPF
3852 DEW = 0.
3853
3854 PRCP1 = PRCPF *0.001
3855 ETP1 = ETP * 0.001
3856 ESNOW = ETP * SNCOVR
3857 ! write(*,*) 'ESNOW,ESDFLX=',ESNOW,ESDFLX
3858 ! ESDFLX = ESD *1000./ DT
3859 ! IF (ESDFLX .lt. ESNOW) THEN
3860 !ek ESD = 0.
3861 ! ESNOW = ESDFLX
3862 ! ESD = 0.
3863 ! ELSE
3864 ESNOW1 = ESNOW * 0.001
3865 !ek ESD = ESD - ESNOW2
3866 ESNOW2 = ESNOW1 * DT
3867 ! ESD = ESD- ESNOW2
3868 ! END IF
3869 ! ETP2 = ETP * 0.001 * DT
3870 ! IF (ICE .NE. 1) THEN
3871 ! IF (ESD .LT. ETP2) THEN
3872 ! BETA = ESD / ETP2
3873 ! ENDIF
3874 ! ENDIF
3875
3876 ! ----------------------------------------------------------------------
3877 ! IF ETP<0 (DOWNWARD) THEN DEWFALL (=FROSTFALL IN THIS CASE).
3878 ! ----------------------------------------------------------------------
3879 BETA = 1.0
3880 IF (ETP <= 0.0) THEN
3881 IF(ETP == 0.) BETA = 0.0
3882 DEW = - ETP * 0.001
3883 ! ETP1 = 0.0
3884 ELSE
3885 IF (ICE /= 1) THEN
3886 ! write(*,*) 'ETP1=',ETP1
3887 IF (SNCOVR < 1.) THEN
3888 CALL EVAPO (ETA1,SMC,NSOIL,CMC,ETP1,DT,ZSOIL, &
3889 SH2O, &
3890 SMCMAX,BEXP,PC,SMCWLT,DKSAT,DWSAT, &
3891 SMCREF,SHDFAC,CMCMAX, &
3892 SMCDRY,CFACTR, &
3893 EDIR1,EC1,ET1,ETT1,SFCTMP,Q2,NROOT,RTDIS, &
3894 FXEXP)
3895 EDIR1 = EDIR1* (1. - SNCOVR)
3896 EC1 = EC1* (1. - SNCOVR)
3897 DO K = 1,NSOIL
3898 ET1 (K) = ET1 (K)* (1. - SNCOVR)
3899 END DO
3900 ETT1 = ETT1*(1.-SNCOVR)
3901 ETNS1 = EDIR1+ EC1+ ETT1
3902 EDIR = EDIR1*1000.
3903 EC = EC1*1000.
3904 DO K = 1,NSOIL
3905 ET (K) = ET1 (K)*1000.
3906 END DO
3907 ETT = ETT1*1000.
3908 ETNS = ETNS1*1000.
3909 ! write(*,*) EDIR*2.5E+9,
3910 ! & EC*2.5E+9,
3911 ! & ET(1)*2.5E+9,
3912 ! & ET(2)*2.5E+9,
3913 ! & ET(3)*2.5E+9,
3914 ! & ET(4)*2.5E+9,
3915 ! & ETT*2.5E+9
3916 ! write(*,*) 'SNCOVR=',SNCOVR
3917 ETNS = ETNS1*1000.
3918 ! write(*,*) 'ESNOW[W/M2],ETNS[W/M2]=',
3919 ! & ESNOW*LSUBS,ETNS*LSUBC
3920 ! write(*,*) 'ETP[W/M2],ETANRG=',
3921 ! & ETP*LSUBS,ETANRG
3922 ! ----------------------------------------------------------------------
3923
3924 ! end IF (SNCOVR .lt. 1.)
3925 END IF
3926
3927 ! end IF (ICE .ne. 1)
3928 END IF
3929
3930 ! end IF (ETP .le. 0.0)
3931 END IF
3932 ! ----------------------------------------------------------------------
3933 ! COMPUTE TOTAL EVAPORATION, SNOW AND NON-SNOW
3934 ! also compute energy units (ETANRG) needed later below
3935 ! ----------------------------------------------------------------------
3936
3937 ! write(*,*)'ESNOW*LSUBS,ETNS*LSUBC,ETANRG=',
3938 ! & ESNOW*LSUBS,ETNS*LSUBC,ETANRG
3939 ETANRG = ESNOW * LSUBS + ETNS * LSUBC
3940
3941 ! ----------------------------------------------------------------------
3942 ! IF PRECIP IS FALLING, CALCULATE HEAT FLUX FROM SNOW SFC TO NEWLY
3943 ! ACCUMULATING PRECIP. NOTE THAT THIS REFLECTS THE FLUX APPROPRIATE FOR
3944 ! THE NOT-YET-UPDATED SKIN TEMPERATURE (T1). ASSUMES TEMPERATURE OF THE
3945 ! SNOWFALL STRIKING THE GROUND IS =SFCTMP (LOWEST MODEL LEVEL AIR TEMP).
3946 ! ----------------------------------------------------------------------
3947 FLX1 = 0.0
3948 IF (SNOWNG) THEN
3949 FLX1 = CPICE * PRCP * (T1- SFCTMP)
3950 ELSE
3951 IF (PRCP > 0.0) FLX1 = CPH2O * PRCP * (T1- SFCTMP)
3952 ! ----------------------------------------------------------------------
3953 ! CALCULATE AN 'EFFECTIVE SNOW-GRND SFC TEMP' (T12) BASED ON HEAT FLUXES
3954 ! BETWEEN THE SNOW PACK AND THE SOIL AND ON NET RADIATION.
3955 ! INCLUDE FLX1 (PRECIP-SNOW SFC) AND FLX2 (FREEZING RAIN LATENT HEAT)
3956 ! FLUXES. FLX1 FROM ABOVE, FLX2 BROUGHT IN VIA COMMOM BLOCK RITE.
3957 ! FLX2 REFLECTS FREEZING RAIN LATENT HEAT FLUX USING T1 CALCULATED IN
3958 ! PENMAN.
3959 ! ----------------------------------------------------------------------
3960 END IF
3961 DSOIL = - (0.5 * ZSOIL (1))
3962 DTOT = SNOWH + DSOIL
3963 !ek T12A = ( (FDOWN-FLX1-FLX2-SIGMA*T24)/RCH
3964 !ek & + TH2 - SFCTMP - BETA*EPSCA ) / RR
3965 !ek T12AX = ( (FDOWN-FLX1-FLX2-SIGMA*T24)/RCH
3966 !ek & + TH2 - SFCTMP - ETANRG/RCH ) / RR
3967 DENOM = 1.0+ DF1 / (DTOT * RR * RCH)
3968 ! write(*,*) 'T12A,T12AX=',T12A,T12AX
3969 !Place for emiss, snow emiss=0.90
3970 ! T12A = ( (FDOWN - FLX1- FLX2- SIGMA * T24)/ RCH &
3971 T12A = ( (FDOWN - FLX1- FLX2- EMISSI_S*SIGMA * T24)/ RCH &
3972 & + TH2- SFCTMP - ETANRG / RCH ) / RR
3973 T12B = DF1 * STC (1) / (DTOT * RR * RCH)
3974
3975 ! ----------------------------------------------------------------------
3976 ! IF THE 'EFFECTIVE SNOW-GRND SFC TEMP' IS AT OR BELOW FREEZING, NO SNOW
3977 ! MELT WILL OCCUR. SET THE SKIN TEMP TO THIS EFFECTIVE TEMP. REDUCE
3978 ! (BY SUBLIMINATION ) OR INCREASE (BY FROST) THE DEPTH OF THE SNOWPACK,
3979 ! DEPENDING ON SIGN OF ETP.
3980 ! UPDATE SOIL HEAT FLUX (SSOIL) USING NEW SKIN TEMPERATURE (T1)
3981 ! SINCE NO SNOWMELT, SET ACCUMULATED SNOWMELT TO ZERO, SET 'EFFECTIVE'
3982 ! PRECIP FROM SNOWMELT TO ZERO, SET PHASE-CHANGE HEAT FLUX FROM SNOWMELT
3983 ! TO ZERO.
3984 ! ----------------------------------------------------------------------
3985 ! SUB-FREEZING BLOCK
3986 ! ----------------------------------------------------------------------
3987 T12 = (SFCTMP + T12A + T12B) / DENOM
3988 IF (T12 <= TFREEZ) THEN
3989 T1 = T12
3990 SSOIL = DF1 * (T1- STC (1)) / DTOT
3991 ! ESD = MAX (0.0, ESD- ETP2)
3992 ESD = MAX(0.0, ESD-ESNOW2)
3993 FLX3 = 0.0
3994 EX = 0.0
3995
3996 SNOMLT = 0.0
3997 ! ----------------------------------------------------------------------
3998 ! IF THE 'EFFECTIVE SNOW-GRND SFC TEMP' IS ABOVE FREEZING, SNOW MELT
3999 ! WILL OCCUR. CALL THE SNOW MELT RATE,EX AND AMT, SNOMLT. REVISE THE
4000 ! EFFECTIVE SNOW DEPTH. REVISE THE SKIN TEMP BECAUSE IT WOULD HAVE CHGD
4001 ! DUE TO THE LATENT HEAT RELEASED BY THE MELTING. CALC THE LATENT HEAT
4002 ! RELEASED, FLX3. SET THE EFFECTIVE PRECIP, PRCP1 TO THE SNOW MELT RATE,
4003 ! EX FOR USE IN SMFLX. ADJUSTMENT TO T1 TO ACCOUNT FOR SNOW PATCHES.
4004 ! CALCULATE QSAT VALID AT FREEZING POINT. NOTE THAT ESAT (SATURATION
4005 ! VAPOR PRESSURE) VALUE OF 6.11E+2 USED HERE IS THAT VALID AT FRZZING
4006 ! POINT. NOTE THAT ETP FROM CALL PENMAN IN SFLX IS IGNORED HERE IN
4007 ! FAVOR OF BULK ETP OVER 'OPEN WATER' AT FREEZING TEMP.
4008 ! UPDATE SOIL HEAT FLUX (S) USING NEW SKIN TEMPERATURE (T1)
4009 ! ----------------------------------------------------------------------
4010 ! ABOVE FREEZING BLOCK
4011 ! ----------------------------------------------------------------------
4012 !ek2 T1 = TFREEZ * SNCOVR + T12 * (1.0 - SNCOVR)
4013 ELSE
4014 !ek QSAT = (0.622*6.11E2)/(SFCPRS-0.378*6.11E2)
4015 !ek ETP = RCH*(QSAT-Q2)/CP
4016 !ek ETP2 = ETP*0.001*DT
4017 T1 = TFREEZ * SNCOVR ** SNOEXP + T12 * (1.0- SNCOVR ** SNOEXP)
4018 BETA = 1.0
4019
4020 ! ----------------------------------------------------------------------
4021 ! IF POTENTIAL EVAP (SUBLIMATION) GREATER THAN DEPTH OF SNOWPACK.
4022 ! BETA<1
4023 ! SNOWPACK HAS SUBLIMATED AWAY, SET DEPTH TO ZERO.
4024 ! ----------------------------------------------------------------------
4025 SSOIL = DF1 * (T1- STC (1)) / DTOT
4026 ! IF (ESD .le. ETP2) THEN
4027 ! BETA = ESD / ETP2
4028 IF (ESD <= ESNOW2) THEN
4029 ESD = 0.0
4030 EX = 0.0
4031 SNOMLT = 0.0
4032 ! ----------------------------------------------------------------------
4033 ! POTENTIAL EVAP (SUBLIMATION) LESS THAN DEPTH OF SNOWPACK, RETAIN
4034 ! BETA=1.
4035 ! SNOWPACK (ESD) REDUCED BY POTENTIAL EVAP RATE
4036 ! ETP3 (CONVERT TO FLUX)
4037 ! ----------------------------------------------------------------------
4038 ELSE
4039 ! ESD = ESD- ETP2
4040 ESD = ESD-ESNOW2
4041 ETP3 = ETP * LSUBC
4042 ! WRITE (*,*) 'SNCOVR=',SNCOVR
4043 ! WRITE (*,*) 'ETP3,ETANRG=',ETP3,ETANRG,'ETP',ETP,'LSUBC',LSUBC
4044 SEH = RCH * (T1- TH2)
4045 T14 = T1* T1
4046 !ek FLX3 = FDOWN - FLX1 - FLX2 - SIGMA*T14 - SSOIL - SEH - ETP3
4047 !ek FLX3X = FDOWN - FLX1 - FLX2 - SIGMA*T14 - SSOIL - SEH - ETANRG
4048 T14 = T14* T14
4049 ! FLX3 = FDOWN - FLX1- FLX2- SIGMA * T14- SSOIL - SEH - ETANRG
4050 FLX3 = FDOWN - FLX1- FLX2- EMISSI_S*SIGMA * T14- SSOIL- SEH- &
4051 & ETANRG
4052 ! WRITE (*,*) 'FLX3,FLX3X=',FLX3,FLX3X
4053 IF (FLX3 <= 0.0) FLX3 = 0.0
4054 ! ----------------------------------------------------------------------
4055 ! SNOWMELT REDUCTION DEPENDING ON SNOW COVER
4056 ! IF SNOW COVER LESS THAN 5% NO SNOWMELT REDUCTION
4057 ! ***NOTE: DOES 'IF' BELOW FAIL TO MATCH THE MELT WATER WITH THE MELT
4058 ! ENERGY?
4059 ! ----------------------------------------------------------------------
4060 EX = FLX3*0.001/ LSUBF
4061 ! IF (SNCOVR .gt. 0.05) EX = EX * SNCOVR
4062
4063 ! ----------------------------------------------------------------------
4064 ! ESDMIN REPRESENTS A SNOWPACK DEPTH THRESHOLD VALUE BELOW WHICH WE
4065 ! CHOOSE NOT TO RETAIN ANY SNOWPACK, AND INSTEAD INCLUDE IT IN SNOWMELT.
4066 ! ----------------------------------------------------------------------
4067 SNOMLT = EX * DT
4068 IF (ESD- SNOMLT >= ESDMIN) THEN
4069 ESD = ESD- SNOMLT
4070 ! ----------------------------------------------------------------------
4071 ! SNOWMELT EXCEEDS SNOW DEPTH
4072 ! ----------------------------------------------------------------------
4073 ELSE
4074 EX = ESD / DT
4075 FLX3 = EX *1000.0* LSUBF
4076 SNOMLT = ESD
4077
4078 ESD = 0.0
4079 ! ----------------------------------------------------------------------
4080 ! END OF 'ESD .LE. ETP2' IF-BLOCK
4081 ! ----------------------------------------------------------------------
4082 END IF
4083 END IF
4084
4085 ! ----------------------------------------------------------------------
4086 ! END OF 'T12 .LE. TFREEZ' IF-BLOCK
4087 ! ----------------------------------------------------------------------
4088 PRCP1 = PRCP1+ EX
4089 ! ----------------------------------------------------------------------
4090 ! FINAL BETA NOW IN HAND, SO COMPUTE EVAPORATION. EVAP EQUALS ETP
4091 ! UNLESS BETA<1.
4092 ! ----------------------------------------------------------------------
4093 !ek ETA = BETA*ETP
4094 ! write(*,*) 'ETA,ETAX,ETAX/ETA=',ETA,ETAX,ETAX/ETA
4095 ! IF (ETAX/ETA .GE. 1.0) THEN
4096 ! write(*,*) '*********************'
4097 ! write(*,*) '*********************'
4098 ! write(*,*) '*********************'
4099 ! write(*,*) 'ETAX/ETA .GE. 1.0 !!!'
4100 ! write(*,*) '*********************'
4101 ! write(*,*) '*********************'
4102 ! write(*,*) '*********************'
4103 ! ENDIF
4104
4105 ! ----------------------------------------------------------------------
4106 ! SET THE EFFECTIVE POTNL EVAPOTRANSP (ETP1) TO ZERO SINCE THIS IS SNOW
4107 ! CASE, SO SURFACE EVAP NOT CALCULATED FROM EDIR, EC, OR ETT IN SMFLX
4108 ! (BELOW).
4109 ! IF SEAICE (ICE=1) SKIP CALL TO SMFLX.
4110 ! SMFLX RETURNS UPDATED SOIL MOISTURE VALUES. IN THIS, THE SNOW PACK
4111 ! CASE, ETA1 IS NOT USED IN CALCULATION OF EVAP.
4112 ! ----------------------------------------------------------------------
4113 ! ETP1 = 0.0
4114 END IF
4115 IF (ICE /= 1) THEN
4116 CALL SMFLX (SMC,NSOIL,CMC,DT,PRCP1,ZSOIL, &
4117 SH2O,SLOPE,KDT,FRZFACT, &
4118 SMCMAX,BEXP,SMCWLT,DKSAT,DWSAT, &
4119 SHDFAC,CMCMAX, &
4120 RUNOFF1,RUNOFF2,RUNOFF3, &
4121 EDIR1,EC1,ET1, &
4122 DRIP)
4123 ! ----------------------------------------------------------------------
4124 ! BEFORE CALL SHFLX IN THIS SNOWPACK CASE, SET ZZ1 AND YY ARGUMENTS TO
4125 ! SPECIAL VALUES THAT ENSURE THAT GROUND HEAT FLUX CALCULATED IN SHFLX
4126 ! MATCHES THAT ALREADY COMPUTER FOR BELOW THE SNOWPACK, THUS THE SFC
4127 ! HEAT FLUX TO BE COMPUTED IN SHFLX WILL EFFECTIVELY BE THE FLUX AT THE
4128 ! SNOW TOP SURFACE. T11 IS A DUMMY ARGUEMENT SO WE WILL NOT USE THE
4129 ! SKIN TEMP VALUE AS REVISED BY SHFLX.
4130 ! ----------------------------------------------------------------------
4131 END IF
4132 ZZ1 = 1.0
4133 YY = STC (1) -0.5* SSOIL * ZSOIL (1)* ZZ1/ DF1
4134
4135 ! ----------------------------------------------------------------------
4136 ! SHFLX WILL CALC/UPDATE THE SOIL TEMPS. NOTE: THE SUB-SFC HEAT FLUX
4137 ! (SSOIL1) AND THE SKIN TEMP (T11) OUTPUT FROM THIS SHFLX CALL ARE NOT
4138 ! USED IN ANY SUBSEQUENT CALCULATIONS. RATHER, THEY ARE DUMMY VARIABLES
4139 ! HERE IN THE SNOPAC CASE, SINCE THE SKIN TEMP AND SUB-SFC HEAT FLUX ARE
4140 ! UPDATED INSTEAD NEAR THE BEGINNING OF THE CALL TO SNOPAC.
4141 ! ----------------------------------------------------------------------
4142 T11 = T1
4143 CALL SHFLX (SSOIL1,STC,SMC,SMCMAX,NSOIL,T11,DT,YY,ZZ1,ZSOIL, &
4144 TBOT,ZBOT,SMCWLT,PSISAT,SH2O,BEXP,F1,DF1,ICE, &
4145 QUARTZ,CSOIL,VEGTYP)
4146
4147 ! ----------------------------------------------------------------------
4148 ! SNOW DEPTH AND DENSITY ADJUSTMENT BASED ON SNOW COMPACTION. YY IS
4149 ! ASSUMED TO BE THE SOIL TEMPERTURE AT THE TOP OF THE SOIL COLUMN.
4150 ! ----------------------------------------------------------------------
4151 IF (ESD > 0.) THEN
4152 CALL SNOWPACK (ESD,DT,SNOWH,SNDENS,T1,YY)
4153 ELSE
4154 ESD = 0.
4155 SNOWH = 0.
4156 SNDENS = 0.
4157 SNCOND = 1.
4158 END IF
4159 ! ----------------------------------------------------------------------
4160 END SUBROUTINE SNOPAC
4161 ! ----------------------------------------------------------------------
4162
4163
4164 SUBROUTINE SNOWPACK (ESD,DTSEC,SNOWH,SNDENS,TSNOW,TSOIL)
4165
4166 ! ----------------------------------------------------------------------
4167 ! SUBROUTINE SNOWPACK
4168 ! ----------------------------------------------------------------------
4169 ! CALCULATE COMPACTION OF SNOWPACK UNDER CONDITIONS OF INCREASING SNOW
4170 ! DENSITY, AS OBTAINED FROM AN APPROXIMATE SOLUTION OF E. ANDERSON'S
4171 ! DIFFERENTIAL EQUATION (3.29), NOAA TECHNICAL REPORT NWS 19, BY VICTOR
4172 ! KOREN, 03/25/95.
4173 ! ----------------------------------------------------------------------
4174 ! ESD WATER EQUIVALENT OF SNOW (M)
4175 ! DTSEC TIME STEP (SEC)
4176 ! SNOWH SNOW DEPTH (M)
4177 ! SNDENS SNOW DENSITY (G/CM3=DIMENSIONLESS FRACTION OF H2O DENSITY)
4178 ! TSNOW SNOW SURFACE TEMPERATURE (K)
4179 ! TSOIL SOIL SURFACE TEMPERATURE (K)
4180
4181 ! SUBROUTINE WILL RETURN NEW VALUES OF SNOWH AND SNDENS
4182 ! ----------------------------------------------------------------------
4183 IMPLICIT NONE
4184
4185 INTEGER :: IPOL, J
4186 REAL, INTENT(IN) :: ESD, DTSEC,TSNOW,TSOIL
4187 REAL, INTENT(INOUT) :: SNOWH, SNDENS
4188 REAL :: BFAC,DSX,DTHR,DW,SNOWHC,PEXP, &
4189 TAVGC,TSNOWC,TSOILC,ESDC,ESDCX
4190 REAL, PARAMETER :: C1 = 0.01, C2 = 21.0, G = 9.81, &
4191 KN = 4000.0
4192 ! ----------------------------------------------------------------------
4193 ! CONVERSION INTO SIMULATION UNITS
4194 ! ----------------------------------------------------------------------
4195 SNOWHC = SNOWH *100.
4196 ESDC = ESD *100.
4197 DTHR = DTSEC /3600.
4198 TSNOWC = TSNOW -273.15
4199 TSOILC = TSOIL -273.15
4200
4201 ! ----------------------------------------------------------------------
4202 ! CALCULATING OF AVERAGE TEMPERATURE OF SNOW PACK
4203 ! ----------------------------------------------------------------------
4204 ! ----------------------------------------------------------------------
4205 ! CALCULATING OF SNOW DEPTH AND DENSITY AS A RESULT OF COMPACTION
4206 ! SNDENS=DS0*(EXP(BFAC*ESD)-1.)/(BFAC*ESD)
4207 ! BFAC=DTHR*C1*EXP(0.08*TAVGC-C2*DS0)
4208 ! NOTE: BFAC*ESD IN SNDENS EQN ABOVE HAS TO BE CAREFULLY TREATED
4209 ! NUMERICALLY BELOW:
4210 ! C1 IS THE FRACTIONAL INCREASE IN DENSITY (1/(CM*HR))
4211 ! C2 IS A CONSTANT (CM3/G) KOJIMA ESTIMATED AS 21 CMS/G
4212 ! ----------------------------------------------------------------------
4213 TAVGC = 0.5* (TSNOWC + TSOILC)
4214 IF (ESDC > 1.E-2) THEN
4215 ESDCX = ESDC
4216 ELSE
4217 ESDCX = 1.E-2
4218 END IF
4219
4220 ! DSX = SNDENS*((DEXP(BFAC*ESDC)-1.)/(BFAC*ESDC))
4221 ! ----------------------------------------------------------------------
4222 ! THE FUNCTION OF THE FORM (e**x-1)/x IMBEDDED IN ABOVE EXPRESSION
4223 ! FOR DSX WAS CAUSING NUMERICAL DIFFICULTIES WHEN THE DENOMINATOR "x"
4224 ! (I.E. BFAC*ESDC) BECAME ZERO OR APPROACHED ZERO (DESPITE THE FACT THAT
4225 ! THE ANALYTICAL FUNCTION (e**x-1)/x HAS A WELL DEFINED LIMIT AS
4226 ! "x" APPROACHES ZERO), HENCE BELOW WE REPLACE THE (e**x-1)/x
4227 ! EXPRESSION WITH AN EQUIVALENT, NUMERICALLY WELL-BEHAVED
4228 ! POLYNOMIAL EXPANSION.
4229
4230 ! NUMBER OF TERMS OF POLYNOMIAL EXPANSION, AND HENCE ITS ACCURACY,
4231 ! IS GOVERNED BY ITERATION LIMIT "IPOL".
4232 ! IPOL GREATER THAN 9 ONLY MAKES A DIFFERENCE ON DOUBLE
4233 ! PRECISION (RELATIVE ERRORS GIVEN IN PERCENT %).
4234 ! IPOL=9, FOR REL.ERROR <~ 1.6 E-6 % (8 SIGNIFICANT DIGITS)
4235 ! IPOL=8, FOR REL.ERROR <~ 1.8 E-5 % (7 SIGNIFICANT DIGITS)
4236 ! IPOL=7, FOR REL.ERROR <~ 1.8 E-4 % ...
4237 ! ----------------------------------------------------------------------
4238 BFAC = DTHR * C1* EXP (0.08* TAVGC - C2* SNDENS)
4239 IPOL = 4
4240 PEXP = 0.
4241 ! PEXP = (1. + PEXP)*BFAC*ESDC/REAL(J+1)
4242 DO J = IPOL,1, -1
4243 PEXP = (1. + PEXP)* BFAC * ESDCX / REAL (J +1)
4244 END DO
4245
4246 PEXP = PEXP + 1.
4247 ! ----------------------------------------------------------------------
4248 ! ABOVE LINE ENDS POLYNOMIAL SUBSTITUTION
4249 ! ----------------------------------------------------------------------
4250 ! END OF KOREAN FORMULATION
4251
4252 ! BASE FORMULATION (COGLEY ET AL., 1990)
4253 ! CONVERT DENSITY FROM G/CM3 TO KG/M3
4254 ! DSM=SNDENS*1000.0
4255
4256 ! DSX=DSM+DTSEC*0.5*DSM*G*ESD/
4257 ! & (1E7*EXP(-0.02*DSM+KN/(TAVGC+273.16)-14.643))
4258
4259 ! & CONVERT DENSITY FROM KG/M3 TO G/CM3
4260 ! DSX=DSX/1000.0
4261
4262 ! END OF COGLEY ET AL. FORMULATION
4263
4264 ! ----------------------------------------------------------------------
4265 ! SET UPPER/LOWER LIMIT ON SNOW DENSITY
4266 ! ----------------------------------------------------------------------
4267 DSX = SNDENS * (PEXP)
4268 IF (DSX > 0.40) DSX = 0.40
4269 IF (DSX < 0.05) DSX = 0.05
4270 ! ----------------------------------------------------------------------
4271 ! UPDATE OF SNOW DEPTH AND DENSITY DEPENDING ON LIQUID WATER DURING
4272 ! SNOWMELT. ASSUMED THAT 13% OF LIQUID WATER CAN BE STORED IN SNOW PER
4273 ! DAY DURING SNOWMELT TILL SNOW DENSITY 0.40.
4274 ! ----------------------------------------------------------------------
4275 SNDENS = DSX
4276 IF (TSNOWC >= 0.) THEN
4277 DW = 0.13* DTHR /24.
4278 SNDENS = SNDENS * (1. - DW) + DW
4279 IF (SNDENS >= 0.40) SNDENS = 0.40
4280 ! ----------------------------------------------------------------------
4281 ! CALCULATE SNOW DEPTH (CM) FROM SNOW WATER EQUIVALENT AND SNOW DENSITY.
4282 ! CHANGE SNOW DEPTH UNITS TO METERS
4283 ! ----------------------------------------------------------------------
4284 END IF
4285 SNOWHC = ESDC / SNDENS
4286 SNOWH = SNOWHC *0.01
4287
4288 ! ----------------------------------------------------------------------
4289 END SUBROUTINE SNOWPACK
4290 ! ----------------------------------------------------------------------
4291
4292 SUBROUTINE SNOWZ0 (SNCOVR,Z0, Z0BRD)
4293
4294 ! ----------------------------------------------------------------------
4295 ! SUBROUTINE SNOWZ0
4296 ! ----------------------------------------------------------------------
4297 ! CALCULATE TOTAL ROUGHNESS LENGTH OVER SNOW
4298 ! SNCOVR FRACTIONAL SNOW COVER
4299 ! Z0 ROUGHNESS LENGTH (m)
4300 ! Z0S SNOW ROUGHNESS LENGTH:=0.001 (m)
4301 ! ----------------------------------------------------------------------
4302 IMPLICIT NONE
4303 REAL, INTENT(IN) :: SNCOVR, Z0BRD
4304 REAL, INTENT(OUT) :: Z0
4305 REAL, PARAMETER :: Z0S=0.001
4306
4307 ! CURRENT NOAH LSM CONDITION - MBEK, 09-OCT-2001
4308 ! Z0S = Z0BRD
4309 Z0 = (1.- SNCOVR)* Z0BRD + SNCOVR * Z0S
4310 ! ----------------------------------------------------------------------
4311 END SUBROUTINE SNOWZ0
4312 ! ----------------------------------------------------------------------
4313
4314
4315 SUBROUTINE SNOW_NEW (TEMP,NEWSN,SNOWH,SNDENS)
4316
4317 ! ----------------------------------------------------------------------
4318 ! SUBROUTINE SNOW_NEW
4319 ! ----------------------------------------------------------------------
4320 ! CALCULATE SNOW DEPTH AND DENSITITY TO ACCOUNT FOR THE NEW SNOWFALL.
4321 ! NEW VALUES OF SNOW DEPTH & DENSITY RETURNED.
4322
4323 ! TEMP AIR TEMPERATURE (K)
4324 ! NEWSN NEW SNOWFALL (M)
4325 ! SNOWH SNOW DEPTH (M)
4326 ! SNDENS SNOW DENSITY (G/CM3=DIMENSIONLESS FRACTION OF H2O DENSITY)
4327 ! ----------------------------------------------------------------------
4328 IMPLICIT NONE
4329 REAL, INTENT(IN) :: NEWSN, TEMP
4330 REAL, INTENT(INOUT) :: SNDENS, SNOWH
4331 REAL :: DSNEW, HNEWC, SNOWHC,NEWSNC,TEMPC
4332
4333 ! ----------------------------------------------------------------------
4334 ! CONVERSION INTO SIMULATION UNITS
4335 ! ----------------------------------------------------------------------
4336 SNOWHC = SNOWH *100.
4337 NEWSNC = NEWSN *100.
4338
4339 ! ----------------------------------------------------------------------
4340 ! CALCULATING NEW SNOWFALL DENSITY DEPENDING ON TEMPERATURE
4341 ! EQUATION FROM GOTTLIB L. 'A GENERAL RUNOFF MODEL FOR SNOWCOVERED
4342 ! AND GLACIERIZED BASIN', 6TH NORDIC HYDROLOGICAL CONFERENCE,
4343 ! VEMADOLEN, SWEDEN, 1980, 172-177PP.
4344 !-----------------------------------------------------------------------
4345 TEMPC = TEMP -273.15
4346 IF (TEMPC <= -15.) THEN
4347 DSNEW = 0.05
4348 ELSE
4349 DSNEW = 0.05+0.0017* (TEMPC +15.)**1.5
4350 END IF
4351 ! ----------------------------------------------------------------------
4352 ! ADJUSTMENT OF SNOW DENSITY DEPENDING ON NEW SNOWFALL
4353 ! ----------------------------------------------------------------------
4354 HNEWC = NEWSNC / DSNEW
4355 IF (SNOWHC + HNEWC .LT. 1.0E-3) THEN
4356 SNDENS = MAX(DSNEW,SNDENS)
4357 ELSE
4358 SNDENS = (SNOWHC * SNDENS + HNEWC * DSNEW)/ (SNOWHC + HNEWC)
4359 ENDIF
4360 SNOWHC = SNOWHC + HNEWC
4361 SNOWH = SNOWHC *0.01
4362
4363 ! ----------------------------------------------------------------------
4364 END SUBROUTINE SNOW_NEW
4365 ! ----------------------------------------------------------------------
4366
4367 SUBROUTINE SRT (RHSTT,EDIR,ET,SH2O,SH2OA,NSOIL,PCPDRP, &
4368 ZSOIL,DWSAT,DKSAT,SMCMAX,BEXP,RUNOFF1, &
4369 RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZX,SICE,AI,BI,CI)
4370
4371 ! ----------------------------------------------------------------------
4372 ! SUBROUTINE SRT
4373 ! ----------------------------------------------------------------------
4374 ! CALCULATE THE RIGHT HAND SIDE OF THE TIME TENDENCY TERM OF THE SOIL
4375 ! WATER DIFFUSION EQUATION. ALSO TO COMPUTE ( PREPARE ) THE MATRIX
4376 ! COEFFICIENTS FOR THE TRI-DIAGONAL MATRIX OF THE IMPLICIT TIME SCHEME.
4377 ! ----------------------------------------------------------------------
4378 IMPLICIT NONE
4379 INTEGER, INTENT(IN) :: NSOIL
4380 INTEGER :: IALP1, IOHINF, J, JJ, K, KS
4381 REAL, INTENT(IN) :: BEXP, DKSAT, DT, DWSAT, EDIR, FRZX, &
4382 KDT, PCPDRP, SLOPE, SMCMAX, SMCWLT
4383 REAL, INTENT(OUT) :: RUNOFF1, RUNOFF2
4384 REAL, DIMENSION(1:NSOIL), INTENT(IN) :: ET, SH2O, SH2OA, SICE, &
4385 ZSOIL
4386 REAL, DIMENSION(1:NSOIL), INTENT(OUT) :: RHSTT
4387 REAL, DIMENSION(1:NSOIL), INTENT(OUT) :: AI, BI, CI
4388 REAL, DIMENSION(1:NSOIL) :: DMAX
4389 REAL :: ACRT, DD, DDT, DDZ, DDZ2, DENOM, &
4390 DENOM2,DICE, DSMDZ, DSMDZ2, DT1, &
4391 FCR,INFMAX,MXSMC,MXSMC2,NUMER,PDDUM, &
4392 PX, SICEMAX,SLOPX, SMCAV, SSTT, &
4393 SUM, VAL, WCND, WCND2, WDF, WDF2
4394 INTEGER, PARAMETER :: CVFRZ = 3
4395
4396 ! ----------------------------------------------------------------------
4397 ! FROZEN GROUND VERSION:
4398 ! REFERENCE FROZEN GROUND PARAMETER, CVFRZ, IS A SHAPE PARAMETER OF
4399 ! AREAL DISTRIBUTION FUNCTION OF SOIL ICE CONTENT WHICH EQUALS 1/CV.
4400 ! CV IS A COEFFICIENT OF SPATIAL VARIATION OF SOIL ICE CONTENT. BASED
4401 ! ON FIELD DATA CV DEPENDS ON AREAL MEAN OF FROZEN DEPTH, AND IT CLOSE
4402 ! TO CONSTANT = 0.6 IF AREAL MEAN FROZEN DEPTH IS ABOVE 20 CM. THAT IS
4403 ! WHY PARAMETER CVFRZ = 3 (INT{1/0.6*0.6}).
4404 ! CURRENT LOGIC DOESN'T ALLOW CVFRZ BE BIGGER THAN 3
4405 ! ----------------------------------------------------------------------
4406
4407 ! ----------------------------------------------------------------------
4408 ! DETERMINE RAINFALL INFILTRATION RATE AND RUNOFF. INCLUDE THE
4409 ! INFILTRATION FORMULE FROM SCHAAKE AND KOREN MODEL.
4410 ! MODIFIED BY Q DUAN
4411 ! ----------------------------------------------------------------------
4412 ! ----------------------------------------------------------------------
4413 ! LET SICEMAX BE THE GREATEST, IF ANY, FROZEN WATER CONTENT WITHIN SOIL
4414 ! LAYERS.
4415 ! ----------------------------------------------------------------------
4416 IOHINF = 1
4417 SICEMAX = 0.0
4418 DO KS = 1,NSOIL
4419 IF (SICE (KS) > SICEMAX) SICEMAX = SICE (KS)
4420 ! ----------------------------------------------------------------------
4421 ! DETERMINE RAINFALL INFILTRATION RATE AND RUNOFF
4422 ! ----------------------------------------------------------------------
4423 END DO
4424 PDDUM = PCPDRP
4425 RUNOFF1 = 0.0
4426
4427 ! ----------------------------------------------------------------------
4428 ! MODIFIED BY Q. DUAN, 5/16/94
4429 ! ----------------------------------------------------------------------
4430 ! IF (IOHINF == 1) THEN
4431
4432 IF (PCPDRP /= 0.0) THEN
4433 DT1 = DT /86400.
4434 SMCAV = SMCMAX - SMCWLT
4435
4436 ! ----------------------------------------------------------------------
4437 ! FROZEN GROUND VERSION:
4438 ! ----------------------------------------------------------------------
4439 DMAX (1)= - ZSOIL (1)* SMCAV
4440
4441 DICE = - ZSOIL (1) * SICE (1)
4442 DMAX (1)= DMAX (1)* (1.0- (SH2OA (1) + SICE (1) - SMCWLT)/ &
4443 SMCAV)
4444
4445 DD = DMAX (1)
4446
4447 ! ----------------------------------------------------------------------
4448 ! FROZEN GROUND VERSION:
4449 ! ----------------------------------------------------------------------
4450 DO KS = 2,NSOIL
4451
4452 DICE = DICE+ ( ZSOIL (KS -1) - ZSOIL (KS) ) * SICE (KS)
4453 DMAX (KS) = (ZSOIL (KS -1) - ZSOIL (KS))* SMCAV
4454 DMAX (KS) = DMAX (KS)* (1.0- (SH2OA (KS) + SICE (KS) &
4455 - SMCWLT)/ SMCAV)
4456 DD = DD+ DMAX (KS)
4457 ! ----------------------------------------------------------------------
4458 ! VAL = (1.-EXP(-KDT*SQRT(DT1)))
4459 ! IN BELOW, REMOVE THE SQRT IN ABOVE
4460 ! ----------------------------------------------------------------------
4461 END DO
4462 VAL = (1. - EXP ( - KDT * DT1))
4463 DDT = DD * VAL
4464 PX = PCPDRP * DT
4465 IF (PX < 0.0) PX = 0.0
4466
4467 ! ----------------------------------------------------------------------
4468 ! FROZEN GROUND VERSION:
4469 ! REDUCTION OF INFILTRATION BASED ON FROZEN GROUND PARAMETERS
4470 ! ----------------------------------------------------------------------
4471 INFMAX = (PX * (DDT / (PX + DDT)))/ DT
4472 FCR = 1.
4473 IF (DICE > 1.E-2) THEN
4474 ACRT = CVFRZ * FRZX / DICE
4475 SUM = 1.
4476 IALP1 = CVFRZ - 1
4477 DO J = 1,IALP1
4478 K = 1
4479 DO JJ = J +1,IALP1
4480 K = K * JJ
4481 END DO
4482 SUM = SUM + (ACRT ** ( CVFRZ - J)) / FLOAT (K)
4483 END DO
4484 FCR = 1. - EXP ( - ACRT) * SUM
4485 END IF
4486
4487 ! ----------------------------------------------------------------------
4488 ! CORRECTION OF INFILTRATION LIMITATION:
4489 ! IF INFMAX .LE. HYDROLIC CONDUCTIVITY ASSIGN INFMAX THE VALUE OF
4490 ! HYDROLIC CONDUCTIVITY
4491 ! ----------------------------------------------------------------------
4492 ! MXSMC = MAX ( SH2OA(1), SH2OA(2) )
4493 INFMAX = INFMAX * FCR
4494
4495 MXSMC = SH2OA (1)
4496 CALL WDFCND (WDF,WCND,MXSMC,SMCMAX,BEXP,DKSAT,DWSAT, &
4497 SICEMAX)
4498 INFMAX = MAX (INFMAX,WCND)
4499
4500 INFMAX = MIN (INFMAX,PX)
4501 IF (PCPDRP > INFMAX) THEN
4502 RUNOFF1 = PCPDRP - INFMAX
4503 PDDUM = INFMAX
4504 END IF
4505 ! ----------------------------------------------------------------------
4506 ! TO AVOID SPURIOUS DRAINAGE BEHAVIOR, 'UPSTREAM DIFFERENCING' IN LINE
4507 ! BELOW REPLACED WITH NEW APPROACH IN 2ND LINE:
4508 ! 'MXSMC = MAX(SH2OA(1), SH2OA(2))'
4509 ! ----------------------------------------------------------------------
4510 END IF
4511
4512 MXSMC = SH2OA (1)
4513 CALL WDFCND (WDF,WCND,MXSMC,SMCMAX,BEXP,DKSAT,DWSAT, &
4514 SICEMAX)
4515 ! ----------------------------------------------------------------------
4516 ! CALC THE MATRIX COEFFICIENTS AI, BI, AND CI FOR THE TOP LAYER
4517 ! ----------------------------------------------------------------------
4518 DDZ = 1. / ( - .5 * ZSOIL (2) )
4519 AI (1) = 0.0
4520 BI (1) = WDF * DDZ / ( - ZSOIL (1) )
4521
4522 ! ----------------------------------------------------------------------
4523 ! CALC RHSTT FOR THE TOP LAYER AFTER CALC'NG THE VERTICAL SOIL MOISTURE
4524 ! GRADIENT BTWN THE TOP AND NEXT TO TOP LAYERS.
4525 ! ----------------------------------------------------------------------
4526 CI (1) = - BI (1)
4527 DSMDZ = ( SH2O (1) - SH2O (2) ) / ( - .5 * ZSOIL (2) )
4528 RHSTT (1) = (WDF * DSMDZ + WCND- PDDUM + EDIR + ET (1))/ ZSOIL (1)
4529
4530 ! ----------------------------------------------------------------------
4531 ! INITIALIZE DDZ2
4532 ! ----------------------------------------------------------------------
4533 SSTT = WDF * DSMDZ + WCND+ EDIR + ET (1)
4534
4535 ! ----------------------------------------------------------------------
4536 ! LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABV PROCESS
4537 ! ----------------------------------------------------------------------
4538 DDZ2 = 0.0
4539 DO K = 2,NSOIL
4540 DENOM2 = (ZSOIL (K -1) - ZSOIL (K))
4541 IF (K /= NSOIL) THEN
4542
4543 ! ----------------------------------------------------------------------
4544 ! AGAIN, TO AVOID SPURIOUS DRAINAGE BEHAVIOR, 'UPSTREAM DIFFERENCING' IN
4545 ! LINE BELOW REPLACED WITH NEW APPROACH IN 2ND LINE:
4546 ! 'MXSMC2 = MAX (SH2OA(K), SH2OA(K+1))'
4547 ! ----------------------------------------------------------------------
4548 SLOPX = 1.
4549
4550 MXSMC2 = SH2OA (K)
4551 CALL WDFCND (WDF2,WCND2,MXSMC2,SMCMAX,BEXP,DKSAT,DWSAT, &
4552 SICEMAX)
4553 ! -----------------------------------------------------------------------
4554 ! CALC SOME PARTIAL PRODUCTS FOR LATER USE IN CALC'NG RHSTT
4555 ! ----------------------------------------------------------------------
4556 DENOM = (ZSOIL (K -1) - ZSOIL (K +1))
4557
4558 ! ----------------------------------------------------------------------
4559 ! CALC THE MATRIX COEF, CI, AFTER CALC'NG ITS PARTIAL PRODUCT
4560 ! ----------------------------------------------------------------------
4561 DSMDZ2 = (SH2O (K) - SH2O (K +1)) / (DENOM * 0.5)
4562 DDZ2 = 2.0 / DENOM
4563 CI (K) = - WDF2 * DDZ2 / DENOM2
4564
4565 ELSE
4566 ! ----------------------------------------------------------------------
4567 ! SLOPE OF BOTTOM LAYER IS INTRODUCED
4568 ! ----------------------------------------------------------------------
4569
4570 ! ----------------------------------------------------------------------
4571 ! RETRIEVE THE SOIL WATER DIFFUSIVITY AND HYDRAULIC CONDUCTIVITY FOR
4572 ! THIS LAYER
4573 ! ----------------------------------------------------------------------
4574 SLOPX = SLOPE
4575 CALL WDFCND (WDF2,WCND2,SH2OA (NSOIL),SMCMAX,BEXP,DKSAT,DWSAT, &
4576 SICEMAX)
4577
4578 ! ----------------------------------------------------------------------
4579 ! CALC A PARTIAL PRODUCT FOR LATER USE IN CALC'NG RHSTT
4580 ! ----------------------------------------------------------------------
4581
4582 ! ----------------------------------------------------------------------
4583 ! SET MATRIX COEF CI TO ZERO
4584 ! ----------------------------------------------------------------------
4585 DSMDZ2 = 0.0
4586 CI (K) = 0.0
4587 ! ----------------------------------------------------------------------
4588 ! CALC RHSTT FOR THIS LAYER AFTER CALC'NG ITS NUMERATOR
4589 ! ----------------------------------------------------------------------
4590 END IF
4591 NUMER = (WDF2 * DSMDZ2) + SLOPX * WCND2- (WDF * DSMDZ) &
4592 - WCND+ ET (K)
4593
4594 ! ----------------------------------------------------------------------
4595 ! CALC MATRIX COEFS, AI, AND BI FOR THIS LAYER
4596 ! ----------------------------------------------------------------------
4597 RHSTT (K) = NUMER / ( - DENOM2)
4598 AI (K) = - WDF * DDZ / DENOM2
4599
4600 ! ----------------------------------------------------------------------
4601 ! RESET VALUES OF WDF, WCND, DSMDZ, AND DDZ FOR LOOP TO NEXT LYR
4602 ! RUNOFF2: SUB-SURFACE OR BASEFLOW RUNOFF
4603 ! ----------------------------------------------------------------------
4604 BI (K) = - ( AI (K) + CI (K) )
4605 IF (K .eq. NSOIL) THEN
4606 RUNOFF2 = SLOPX * WCND2
4607 END IF
4608 IF (K .ne. NSOIL) THEN
4609 WDF = WDF2
4610 WCND = WCND2
4611 DSMDZ = DSMDZ2
4612 DDZ = DDZ2
4613 END IF
4614 END DO
4615 ! ----------------------------------------------------------------------
4616 END SUBROUTINE SRT
4617 ! ----------------------------------------------------------------------
4618
4619 SUBROUTINE SSTEP (SH2OOUT,SH2OIN,CMC,RHSTT,RHSCT,DT, &
4620 NSOIL,SMCMAX,CMCMAX,RUNOFF3,ZSOIL,SMC,SICE, &
4621 AI,BI,CI)
4622
4623 ! ----------------------------------------------------------------------
4624 ! SUBROUTINE SSTEP
4625 ! ----------------------------------------------------------------------
4626 ! CALCULATE/UPDATE SOIL MOISTURE CONTENT VALUES AND CANOPY MOISTURE
4627 ! CONTENT VALUES.
4628 ! ----------------------------------------------------------------------
4629 IMPLICIT NONE
4630 INTEGER, INTENT(IN) :: NSOIL
4631 INTEGER :: I, K, KK11
4632
4633 REAL, INTENT(IN) :: CMCMAX, DT, SMCMAX
4634 REAL, INTENT(OUT) :: RUNOFF3
4635 REAL, INTENT(INOUT) :: CMC
4636 REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SH2OIN, SICE, ZSOIL
4637 REAL, DIMENSION(1:NSOIL), INTENT(OUT) :: SH2OOUT
4638 REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: RHSTT, SMC
4639 REAL, DIMENSION(1:NSOIL), INTENT(INOUT) :: AI, BI, CI
4640 REAL, DIMENSION(1:NSOIL) :: RHSTTin
4641 REAL, DIMENSION(1:NSOIL) :: CIin
4642 REAL :: DDZ, RHSCT, STOT, WPLUS
4643
4644 ! ----------------------------------------------------------------------
4645 ! CREATE 'AMOUNT' VALUES OF VARIABLES TO BE INPUT TO THE
4646 ! TRI-DIAGONAL MATRIX ROUTINE.
4647 ! ----------------------------------------------------------------------
4648 DO K = 1,NSOIL
4649 RHSTT (K) = RHSTT (K) * DT
4650 AI (K) = AI (K) * DT
4651 BI (K) = 1. + BI (K) * DT
4652 CI (K) = CI (K) * DT
4653 END DO
4654 ! ----------------------------------------------------------------------
4655 ! COPY VALUES FOR INPUT VARIABLES BEFORE CALL TO ROSR12
4656 ! ----------------------------------------------------------------------
4657 DO K = 1,NSOIL
4658 RHSTTin (K) = RHSTT (K)
4659 END DO
4660 DO K = 1,NSOIL
4661 CIin (K) = CI (K)
4662 END DO
4663 ! ----------------------------------------------------------------------
4664 ! CALL ROSR12 TO SOLVE THE TRI-DIAGONAL MATRIX
4665 ! ----------------------------------------------------------------------
4666 CALL ROSR12 (CI,AI,BI,CIin,RHSTTin,RHSTT,NSOIL)
4667 ! ----------------------------------------------------------------------
4668 ! SUM THE PREVIOUS SMC VALUE AND THE MATRIX SOLUTION TO GET A
4669 ! NEW VALUE. MIN ALLOWABLE VALUE OF SMC WILL BE 0.02.
4670 ! RUNOFF3: RUNOFF WITHIN SOIL LAYERS
4671 ! ----------------------------------------------------------------------
4672 WPLUS = 0.0
4673 RUNOFF3 = 0.
4674
4675 DDZ = - ZSOIL (1)
4676 DO K = 1,NSOIL
4677 IF (K /= 1) DDZ = ZSOIL (K - 1) - ZSOIL (K)
4678 SH2OOUT (K) = SH2OIN (K) + CI (K) + WPLUS / DDZ
4679 STOT = SH2OOUT (K) + SICE (K)
4680 IF (STOT > SMCMAX) THEN
4681 IF (K .eq. 1) THEN
4682 DDZ = - ZSOIL (1)
4683 ELSE
4684 KK11 = K - 1
4685 DDZ = - ZSOIL (K) + ZSOIL (KK11)
4686 END IF
4687 WPLUS = (STOT - SMCMAX) * DDZ
4688 ELSE
4689 WPLUS = 0.
4690 END IF
4691 SMC (K) = MAX ( MIN (STOT,SMCMAX),0.02 )
4692 SH2OOUT (K) = MAX ( (SMC (K) - SICE (K)),0.0)
4693 END DO
4694
4695 ! ----------------------------------------------------------------------
4696 ! UPDATE CANOPY WATER CONTENT/INTERCEPTION (CMC). CONVERT RHSCT TO
4697 ! AN 'AMOUNT' VALUE AND ADD TO PREVIOUS CMC VALUE TO GET NEW CMC.
4698 ! ----------------------------------------------------------------------
4699 RUNOFF3 = WPLUS
4700 CMC = CMC + DT * RHSCT
4701 IF (CMC < 1.E-20) CMC = 0.0
4702 CMC = MIN (CMC,CMCMAX)
4703
4704 ! ----------------------------------------------------------------------
4705 END SUBROUTINE SSTEP
4706 ! ----------------------------------------------------------------------
4707
4708 SUBROUTINE TBND (TU,TB,ZSOIL,ZBOT,K,NSOIL,TBND1)
4709
4710 ! ----------------------------------------------------------------------
4711 ! SUBROUTINE TBND
4712 ! ----------------------------------------------------------------------
4713 ! CALCULATE TEMPERATURE ON THE BOUNDARY OF THE LAYER BY INTERPOLATION OF
4714 ! THE MIDDLE LAYER TEMPERATURES
4715 ! ----------------------------------------------------------------------
4716 IMPLICIT NONE
4717 INTEGER, INTENT(IN) :: NSOIL
4718 INTEGER :: K
4719 REAL, INTENT(IN) :: TB, TU, ZBOT
4720 REAL, INTENT(OUT) :: TBND1
4721 REAL, DIMENSION(1:NSOIL), INTENT(IN) :: ZSOIL
4722 REAL :: ZB, ZUP
4723 REAL, PARAMETER :: T0 = 273.15
4724
4725 ! ----------------------------------------------------------------------
4726 ! USE SURFACE TEMPERATURE ON THE TOP OF THE FIRST LAYER
4727 ! ----------------------------------------------------------------------
4728 IF (K == 1) THEN
4729 ZUP = 0.
4730 ELSE
4731 ZUP = ZSOIL (K -1)
4732 END IF
4733 ! ----------------------------------------------------------------------
4734 ! USE DEPTH OF THE CONSTANT BOTTOM TEMPERATURE WHEN INTERPOLATE
4735 ! TEMPERATURE INTO THE LAST LAYER BOUNDARY
4736 ! ----------------------------------------------------------------------
4737 IF (K == NSOIL) THEN
4738 ZB = 2.* ZBOT - ZSOIL (K)
4739 ELSE
4740 ZB = ZSOIL (K +1)
4741 END IF
4742 ! ----------------------------------------------------------------------
4743 ! LINEAR INTERPOLATION BETWEEN THE AVERAGE LAYER TEMPERATURES
4744 ! ----------------------------------------------------------------------
4745
4746 TBND1 = TU + (TB - TU)* (ZUP - ZSOIL (K))/ (ZUP - ZB)
4747 ! ----------------------------------------------------------------------
4748 END SUBROUTINE TBND
4749 ! ----------------------------------------------------------------------
4750
4751
4752 SUBROUTINE TDFCND ( DF, SMC, QZ, SMCMAX, SH2O)
4753
4754 ! ----------------------------------------------------------------------
4755 ! SUBROUTINE TDFCND
4756 ! ----------------------------------------------------------------------
4757 ! CALCULATE THERMAL DIFFUSIVITY AND CONDUCTIVITY OF THE SOIL FOR A GIVEN
4758 ! POINT AND TIME.
4759 ! ----------------------------------------------------------------------
4760 ! PETERS-LIDARD APPROACH (PETERS-LIDARD et al., 1998)
4761 ! June 2001 CHANGES: FROZEN SOIL CONDITION.
4762 ! ----------------------------------------------------------------------
4763 IMPLICIT NONE
4764 REAL, INTENT(IN) :: QZ, SMC, SMCMAX, SH2O
4765 REAL, INTENT(OUT) :: DF
4766 REAL :: AKE, GAMMD, THKDRY, THKICE, THKO, &
4767 THKQTZ,THKSAT,THKS,THKW,SATRATIO,XU, &
4768 XUNFROZ
4769
4770 ! ----------------------------------------------------------------------
4771 ! WE NOW GET QUARTZ AS AN INPUT ARGUMENT (SET IN ROUTINE REDPRM):
4772 ! DATA QUARTZ /0.82, 0.10, 0.25, 0.60, 0.52,
4773 ! & 0.35, 0.60, 0.40, 0.82/
4774 ! ----------------------------------------------------------------------
4775 ! IF THE SOIL HAS ANY MOISTURE CONTENT COMPUTE A PARTIAL SUM/PRODUCT
4776 ! OTHERWISE USE A CONSTANT VALUE WHICH WORKS WELL WITH MOST SOILS
4777 ! ----------------------------------------------------------------------
4778 ! THKW ......WATER THERMAL CONDUCTIVITY
4779 ! THKQTZ ....THERMAL CONDUCTIVITY FOR QUARTZ
4780 ! THKO ......THERMAL CONDUCTIVITY FOR OTHER SOIL COMPONENTS
4781 ! THKS ......THERMAL CONDUCTIVITY FOR THE SOLIDS COMBINED(QUARTZ+OTHER)
4782 ! THKICE ....ICE THERMAL CONDUCTIVITY
4783 ! SMCMAX ....POROSITY (= SMCMAX)
4784 ! QZ .........QUARTZ CONTENT (SOIL TYPE DEPENDENT)
4785 ! ----------------------------------------------------------------------
4786 ! USE AS IN PETERS-LIDARD, 1998 (MODIF. FROM JOHANSEN, 1975).
4787
4788 ! PABLO GRUNMANN, 08/17/98
4789 ! REFS.:
4790 ! FAROUKI, O.T.,1986: THERMAL PROPERTIES OF SOILS. SERIES ON ROCK
4791 ! AND SOIL MECHANICS, VOL. 11, TRANS TECH, 136 PP.
4792 ! JOHANSEN, O., 1975: THERMAL CONDUCTIVITY OF SOILS. PH.D. THESIS,
4793 ! UNIVERSITY OF TRONDHEIM,
4794 ! PETERS-LIDARD, C. D., ET AL., 1998: THE EFFECT OF SOIL THERMAL
4795 ! CONDUCTIVITY PARAMETERIZATION ON SURFACE ENERGY FLUXES
4796 ! AND TEMPERATURES. JOURNAL OF THE ATMOSPHERIC SCIENCES,
4797 ! VOL. 55, PP. 1209-1224.
4798 ! ----------------------------------------------------------------------
4799 ! NEEDS PARAMETERS
4800 ! POROSITY(SOIL TYPE):
4801 ! POROS = SMCMAX
4802 ! SATURATION RATIO:
4803 ! PARAMETERS W/(M.K)
4804 SATRATIO = SMC / SMCMAX
4805 THKICE = 2.2
4806 THKW = 0.57
4807 ! IF (QZ .LE. 0.2) THKO = 3.0
4808 THKO = 2.0
4809 ! SOLIDS' CONDUCTIVITY
4810 ! QUARTZ' CONDUCTIVITY
4811 THKQTZ = 7.7
4812
4813 ! UNFROZEN FRACTION (FROM 1., i.e., 100%LIQUID, TO 0. (100% FROZEN))
4814 THKS = (THKQTZ ** QZ)* (THKO ** (1. - QZ))
4815
4816 ! UNFROZEN VOLUME FOR SATURATION (POROSITY*XUNFROZ)
4817 XUNFROZ = SH2O / SMC
4818 ! SATURATED THERMAL CONDUCTIVITY
4819 XU = XUNFROZ * SMCMAX
4820
4821 ! DRY DENSITY IN KG/M3
4822 THKSAT = THKS ** (1. - SMCMAX)* THKICE ** (SMCMAX - XU)* THKW ** &
4823 (XU)
4824
4825 ! DRY THERMAL CONDUCTIVITY IN W.M-1.K-1
4826 GAMMD = (1. - SMCMAX)*2700.
4827
4828 THKDRY = (0.135* GAMMD+ 64.7)/ (2700. - 0.947* GAMMD)
4829 ! FROZEN
4830 IF ( (SH2O + 0.0005) < SMC ) THEN
4831 AKE = SATRATIO
4832 ! UNFROZEN
4833 ! RANGE OF VALIDITY FOR THE KERSTEN NUMBER (AKE)
4834 ELSE
4835
4836 ! KERSTEN NUMBER (USING "FINE" FORMULA, VALID FOR SOILS CONTAINING AT
4837 ! LEAST 5% OF PARTICLES WITH DIAMETER LESS THAN 2.E-6 METERS.)
4838 ! (FOR "COARSE" FORMULA, SEE PETERS-LIDARD ET AL., 1998).
4839
4840 IF ( SATRATIO > 0.1 ) THEN
4841
4842 AKE = LOG10 (SATRATIO) + 1.0
4843
4844 ! USE K = KDRY
4845 ELSE
4846
4847 AKE = 0.0
4848 END IF
4849 ! THERMAL CONDUCTIVITY
4850
4851 END IF
4852
4853 DF = AKE * (THKSAT - THKDRY) + THKDRY
4854 ! ----------------------------------------------------------------------
4855 END SUBROUTINE TDFCND
4856 ! ----------------------------------------------------------------------
4857
4858 SUBROUTINE TMPAVG (TAVG,TUP,TM,TDN,ZSOIL,NSOIL,K)
4859
4860 ! ----------------------------------------------------------------------
4861 ! SUBROUTINE TMPAVG
4862 ! ----------------------------------------------------------------------
4863 ! CALCULATE SOIL LAYER AVERAGE TEMPERATURE (TAVG) IN FREEZING/THAWING
4864 ! LAYER USING UP, DOWN, AND MIDDLE LAYER TEMPERATURES (TUP, TDN, TM),
4865 ! WHERE TUP IS AT TOP BOUNDARY OF LAYER, TDN IS AT BOTTOM BOUNDARY OF
4866 ! LAYER. TM IS LAYER PROGNOSTIC STATE TEMPERATURE.
4867 ! ----------------------------------------------------------------------
4868 IMPLICIT NONE
4869 INTEGER K
4870
4871 INTEGER NSOIL
4872 REAL DZ
4873 REAL DZH
4874 REAL T0
4875 REAL TAVG
4876 REAL TDN
4877 REAL TM
4878 REAL TUP
4879 REAL X0
4880 REAL XDN
4881 REAL XUP
4882
4883 REAL ZSOIL (NSOIL)
4884
4885 ! ----------------------------------------------------------------------
4886 PARAMETER (T0 = 2.7315E2)
4887 IF (K .eq. 1) THEN
4888 DZ = - ZSOIL (1)
4889 ELSE
4890 DZ = ZSOIL (K -1) - ZSOIL (K)
4891 END IF
4892
4893 DZH = DZ *0.5
4894 IF (TUP .lt. T0) THEN
4895 IF (TM .lt. T0) THEN
4896 ! ----------------------------------------------------------------------
4897 ! TUP, TM, TDN < T0
4898 ! ----------------------------------------------------------------------
4899 IF (TDN .lt. T0) THEN
4900 TAVG = (TUP + 2.0* TM + TDN)/ 4.0
4901 ! ----------------------------------------------------------------------
4902 ! TUP & TM < T0, TDN .ge. T0
4903 ! ----------------------------------------------------------------------
4904 ELSE
4905 X0 = (T0- TM) * DZH / (TDN - TM)
4906 TAVG = 0.5 * (TUP * DZH + TM * (DZH + X0) + T0* ( &
4907 & 2.* DZH - X0)) / DZ
4908 END IF
4909 ELSE
4910 ! ----------------------------------------------------------------------
4911 ! TUP < T0, TM .ge. T0, TDN < T0
4912 ! ----------------------------------------------------------------------
4913 IF (TDN .lt. T0) THEN
4914 XUP = (T0- TUP) * DZH / (TM - TUP)
4915 XDN = DZH - (T0- TM) * DZH / (TDN - TM)
4916 TAVG = 0.5 * (TUP * XUP + T0* (2.* DZ - XUP - XDN) &
4917 & + TDN * XDN) / DZ
4918 ! ----------------------------------------------------------------------
4919 ! TUP < T0, TM .ge. T0, TDN .ge. T0
4920 ! ----------------------------------------------------------------------
4921 ELSE
4922 XUP = (T0- TUP) * DZH / (TM - TUP)
4923 TAVG = 0.5 * (TUP * XUP + T0* (2.* DZ - XUP)) / DZ
4924 END IF
4925 END IF
4926 ELSE
4927 IF (TM .lt. T0) THEN
4928 ! ----------------------------------------------------------------------
4929 ! TUP .ge. T0, TM < T0, TDN < T0
4930 ! ----------------------------------------------------------------------
4931 IF (TDN .lt. T0) THEN
4932 XUP = DZH - (T0- TUP) * DZH / (TM - TUP)
4933 TAVG = 0.5 * (T0* (DZ - XUP) + TM * (DZH + XUP) &
4934 & + TDN * DZH) / DZ
4935 ! ----------------------------------------------------------------------
4936 ! TUP .ge. T0, TM < T0, TDN .ge. T0
4937 ! ----------------------------------------------------------------------
4938 ELSE
4939 XUP = DZH - (T0- TUP) * DZH / (TM - TUP)
4940 XDN = (T0- TM) * DZH / (TDN - TM)
4941 TAVG = 0.5 * (T0* (2.* DZ - XUP - XDN) + TM * &
4942 & (XUP + XDN)) / DZ
4943 END IF
4944 ELSE
4945 ! ----------------------------------------------------------------------
4946 ! TUP .ge. T0, TM .ge. T0, TDN < T0
4947 ! ----------------------------------------------------------------------
4948 IF (TDN .lt. T0) THEN
4949 XDN = DZH - (T0- TM) * DZH / (TDN - TM)
4950 TAVG = (T0* (DZ - XDN) +0.5* (T0+ TDN)* XDN) / DZ
4951 ! ----------------------------------------------------------------------
4952 ! TUP .ge. T0, TM .ge. T0, TDN .ge. T0
4953 ! ----------------------------------------------------------------------
4954 ELSE
4955 TAVG = (TUP + 2.0* TM + TDN) / 4.0
4956 END IF
4957 END IF
4958 END IF
4959 ! ----------------------------------------------------------------------
4960 END SUBROUTINE TMPAVG
4961 ! ----------------------------------------------------------------------
4962
4963 SUBROUTINE TRANSP (ET,NSOIL,ETP1,SMC,CMC,ZSOIL,SHDFAC,SMCWLT, &
4964 & CMCMAX,PC,CFACTR,SMCREF,SFCTMP,Q2,NROOT, &
4965 & RTDIS)
4966
4967 ! ----------------------------------------------------------------------
4968 ! SUBROUTINE TRANSP
4969 ! ----------------------------------------------------------------------
4970 ! CALCULATE TRANSPIRATION FOR THE VEG CLASS.
4971 ! ----------------------------------------------------------------------
4972 IMPLICIT NONE
4973 INTEGER I
4974 INTEGER K
4975 INTEGER NSOIL
4976
4977 INTEGER NROOT
4978 REAL CFACTR
4979 REAL CMC
4980 REAL CMCMAX
4981 REAL DENOM
4982 REAL ET (NSOIL)
4983 REAL ETP1
4984 REAL ETP1A
4985 !.....REAL PART(NSOIL)
4986 REAL GX (7)
4987 REAL PC
4988 REAL Q2
4989 REAL RTDIS (NSOIL)
4990 REAL RTX
4991 REAL SFCTMP
4992 REAL SGX
4993 REAL SHDFAC
4994 REAL SMC (NSOIL)
4995 REAL SMCREF
4996 REAL SMCWLT
4997
4998 ! ----------------------------------------------------------------------
4999 ! INITIALIZE PLANT TRANSP TO ZERO FOR ALL SOIL LAYERS.
5000 ! ----------------------------------------------------------------------
5001 REAL ZSOIL (NSOIL)
5002 DO K = 1,NSOIL
5003 ET (K) = 0.
5004 ! ----------------------------------------------------------------------
5005 ! CALCULATE AN 'ADJUSTED' POTENTIAL TRANSPIRATION
5006 ! IF STATEMENT BELOW TO AVOID TANGENT LINEAR PROBLEMS NEAR ZERO
5007 ! NOTE: GX AND OTHER TERMS BELOW REDISTRIBUTE TRANSPIRATION BY LAYER,
5008 ! ET(K), AS A FUNCTION OF SOIL MOISTURE AVAILABILITY, WHILE PRESERVING
5009 ! TOTAL ETP1A.
5010 ! ----------------------------------------------------------------------
5011 END DO
5012 IF (CMC .ne. 0.0) THEN
5013 ETP1A = SHDFAC * PC * ETP1 * (1.0- (CMC / CMCMAX) ** CFACTR)
5014 ELSE
5015 ETP1A = SHDFAC * PC * ETP1
5016 END IF
5017 SGX = 0.0
5018 DO I = 1,NROOT
5019 GX (I) = ( SMC (I) - SMCWLT ) / ( SMCREF - SMCWLT )
5020 GX (I) = MAX ( MIN ( GX (I), 1. ), 0. )
5021 SGX = SGX + GX (I)
5022 END DO
5023
5024 SGX = SGX / NROOT
5025 DENOM = 0.
5026 DO I = 1,NROOT
5027 RTX = RTDIS (I) + GX (I) - SGX
5028 GX (I) = GX (I) * MAX ( RTX, 0. )
5029 DENOM = DENOM + GX (I)
5030 END DO
5031
5032 IF (DENOM .le. 0.0) DENOM = 1.
5033 DO I = 1,NROOT
5034 ET (I) = ETP1A * GX (I) / DENOM
5035 ! ----------------------------------------------------------------------
5036 ! ABOVE CODE ASSUMES A VERTICALLY UNIFORM ROOT DISTRIBUTION
5037 ! CODE BELOW TESTS A VARIABLE ROOT DISTRIBUTION
5038 ! ----------------------------------------------------------------------
5039 ! ET(1) = ( ZSOIL(1) / ZSOIL(NROOT) ) * GX * ETP1A
5040 ! ET(1) = ( ZSOIL(1) / ZSOIL(NROOT) ) * ETP1A
5041 ! ----------------------------------------------------------------------
5042 ! USING ROOT DISTRIBUTION AS WEIGHTING FACTOR
5043 ! ----------------------------------------------------------------------
5044 ! ET(1) = RTDIS(1) * ETP1A
5045 ! ET(1) = ETP1A * PART(1)
5046 ! ----------------------------------------------------------------------
5047 ! LOOP DOWN THRU THE SOIL LAYERS REPEATING THE OPERATION ABOVE,
5048 ! BUT USING THE THICKNESS OF THE SOIL LAYER (RATHER THAN THE
5049 ! ABSOLUTE DEPTH OF EACH LAYER) IN THE FINAL CALCULATION.
5050 ! ----------------------------------------------------------------------
5051 ! DO K = 2,NROOT
5052 ! GX = ( SMC(K) - SMCWLT ) / ( SMCREF - SMCWLT )
5053 ! GX = MAX ( MIN ( GX, 1. ), 0. )
5054 ! TEST CANOPY RESISTANCE
5055 ! GX = 1.0
5056 ! ET(K) = ((ZSOIL(K)-ZSOIL(K-1))/ZSOIL(NROOT))*GX*ETP1A
5057 ! ET(K) = ((ZSOIL(K)-ZSOIL(K-1))/ZSOIL(NROOT))*ETP1A
5058 ! ----------------------------------------------------------------------
5059 ! USING ROOT DISTRIBUTION AS WEIGHTING FACTOR
5060 ! ----------------------------------------------------------------------
5061 ! ET(K) = RTDIS(K) * ETP1A
5062 ! ET(K) = ETP1A*PART(K)
5063 ! END DO
5064 END DO
5065 ! ----------------------------------------------------------------------
5066 END SUBROUTINE TRANSP
5067 ! ----------------------------------------------------------------------
5068
5069 SUBROUTINE WDFCND (WDF,WCND,SMC,SMCMAX,BEXP,DKSAT,DWSAT, &
5070 & SICEMAX)
5071
5072 ! ----------------------------------------------------------------------
5073 ! SUBROUTINE WDFCND
5074 ! ----------------------------------------------------------------------
5075 ! CALCULATE SOIL WATER DIFFUSIVITY AND SOIL HYDRAULIC CONDUCTIVITY.
5076 ! ----------------------------------------------------------------------
5077 IMPLICIT NONE
5078 REAL BEXP
5079 REAL DKSAT
5080 REAL DWSAT
5081 REAL EXPON
5082 REAL FACTR1
5083 REAL FACTR2
5084 REAL SICEMAX
5085 REAL SMC
5086 REAL SMCMAX
5087 REAL VKwgt
5088 REAL WCND
5089
5090 ! ----------------------------------------------------------------------
5091 ! CALC THE RATIO OF THE ACTUAL TO THE MAX PSBL SOIL H2O CONTENT
5092 ! ----------------------------------------------------------------------
5093 REAL WDF
5094 FACTR1 = 0.2 / SMCMAX
5095
5096 ! ----------------------------------------------------------------------
5097 ! PREP AN EXPNTL COEF AND CALC THE SOIL WATER DIFFUSIVITY
5098 ! ----------------------------------------------------------------------
5099 FACTR2 = SMC / SMCMAX
5100 EXPON = BEXP + 2.0
5101
5102 ! ----------------------------------------------------------------------
5103 ! FROZEN SOIL HYDRAULIC DIFFUSIVITY. VERY SENSITIVE TO THE VERTICAL
5104 ! GRADIENT OF UNFROZEN WATER. THE LATTER GRADIENT CAN BECOME VERY
5105 ! EXTREME IN FREEZING/THAWING SITUATIONS, AND GIVEN THE RELATIVELY
5106 ! FEW AND THICK SOIL LAYERS, THIS GRADIENT SUFFERES SERIOUS
5107 ! TRUNCTION ERRORS YIELDING ERRONEOUSLY HIGH VERTICAL TRANSPORTS OF
5108 ! UNFROZEN WATER IN BOTH DIRECTIONS FROM HUGE HYDRAULIC DIFFUSIVITY.
5109 ! THEREFORE, WE FOUND WE HAD TO ARBITRARILY CONSTRAIN WDF
5110 ! --
5111 ! VERSION D_10CM: ........ FACTR1 = 0.2/SMCMAX
5112 ! WEIGHTED APPROACH...................... PABLO GRUNMANN, 28_SEP_1999.
5113 ! ----------------------------------------------------------------------
5114 WDF = DWSAT * FACTR2 ** EXPON
5115 IF (SICEMAX .gt. 0.0) THEN
5116 VKWGT = 1./ (1. + (500.* SICEMAX)**3.)
5117 WDF = VKWGT * WDF + (1. - VKWGT)* DWSAT * FACTR1** EXPON
5118 ! ----------------------------------------------------------------------
5119 ! RESET THE EXPNTL COEF AND CALC THE HYDRAULIC CONDUCTIVITY
5120 ! ----------------------------------------------------------------------
5121 END IF
5122 EXPON = (2.0 * BEXP) + 3.0
5123 WCND = DKSAT * FACTR2 ** EXPON
5124
5125 ! ----------------------------------------------------------------------
5126 END SUBROUTINE WDFCND
5127 ! ----------------------------------------------------------------------
5128
5129 SUBROUTINE SFCDIF_off (ZLM,Z0,THZ0,THLM,SFCSPD,CZIL,AKMS,AKHS)
5130
5131 ! ----------------------------------------------------------------------
5132 ! SUBROUTINE SFCDIF (renamed SFCDIF_off to avoid clash with Eta PBL)
5133 ! ----------------------------------------------------------------------
5134 ! CALCULATE SURFACE LAYER EXCHANGE COEFFICIENTS VIA ITERATIVE PROCESS.
5135 ! SEE CHEN ET AL (1997, BLM)
5136 ! ----------------------------------------------------------------------
5137
5138 IMPLICIT NONE
5139 REAL WWST, WWST2, G, VKRM, EXCM, BETA, BTG, ELFC, WOLD, WNEW
5140 REAL PIHF, EPSU2, EPSUST, EPSIT, EPSA, ZTMIN, ZTMAX, HPBL, &
5141 & SQVISC
5142 REAL RIC, RRIC, FHNEU, RFC, RFAC, ZZ, PSLMU, PSLMS, PSLHU, &
5143 & PSLHS
5144 REAL XX, PSPMU, YY, PSPMS, PSPHU, PSPHS, ZLM, Z0, THZ0, THLM
5145 REAL SFCSPD, CZIL, AKMS, AKHS, ZILFC, ZU, ZT, RDZ, CXCH
5146 REAL DTHV, DU2, BTGH, WSTAR2, USTAR, ZSLU, ZSLT, RLOGU, RLOGT
5147 REAL RLMO, ZETALT, ZETALU, ZETAU, ZETAT, XLU4, XLT4, XU4, XT4
5148 !CC ......REAL ZTFC
5149
5150 REAL XLU, XLT, XU, XT, PSMZ, SIMM, PSHZ, SIMH, USTARK, RLMN, &
5151 & RLMA
5152
5153 INTEGER ITRMX, ILECH, ITR
5154 PARAMETER &
5155 & (WWST = 1.2,WWST2 = WWST * WWST,G = 9.8,VKRM = 0.40, &
5156 & EXCM = 0.001 &
5157 & ,BETA = 1./270.,BTG = BETA * G,ELFC = VKRM * BTG &
5158 & ,WOLD =.15,WNEW = 1. - WOLD,ITRMX = 05, &
5159 & PIHF = 3.14159265/2.)
5160 PARAMETER &
5161 & (EPSU2 = 1.E-4,EPSUST = 0.07,EPSIT = 1.E-4,EPSA = 1.E-8 &
5162 & ,ZTMIN = -5.,ZTMAX = 1.,HPBL = 1000.0 &
5163 & ,SQVISC = 258.2)
5164 PARAMETER &
5165 & (RIC = 0.183,RRIC = 1.0/ RIC,FHNEU = 0.8,RFC = 0.191 &
5166 & ,RFAC = RIC / (FHNEU * RFC * RFC))
5167
5168 ! ----------------------------------------------------------------------
5169 ! NOTE: THE TWO CODE BLOCKS BELOW DEFINE FUNCTIONS
5170 ! ----------------------------------------------------------------------
5171 ! LECH'S SURFACE FUNCTIONS
5172 ! ----------------------------------------------------------------------
5173 PSLMU (ZZ)= -0.96* log (1.0-4.5* ZZ)
5174 PSLMS (ZZ)= ZZ * RRIC -2.076* (1. -1./ (ZZ +1.))
5175 PSLHU (ZZ)= -0.96* log (1.0-4.5* ZZ)
5176
5177 ! ----------------------------------------------------------------------
5178 ! PAULSON'S SURFACE FUNCTIONS
5179 ! ----------------------------------------------------------------------
5180 PSLHS (ZZ)= ZZ * RFAC -2.076* (1. -1./ (ZZ +1.))
5181 PSPMU (XX)= -2.* log ( (XX +1.)*0.5) - log ( (XX * XX +1.)*0.5) &
5182 & +2.* ATAN (XX) &
5183 &- PIHF
5184 PSPMS (YY)= 5.* YY
5185 PSPHU (XX)= -2.* log ( (XX * XX +1.)*0.5)
5186
5187 ! ----------------------------------------------------------------------
5188 ! THIS ROUTINE SFCDIF CAN HANDLE BOTH OVER OPEN WATER (SEA, OCEAN) AND
5189 ! OVER SOLID SURFACE (LAND, SEA-ICE).
5190 ! ----------------------------------------------------------------------
5191 PSPHS (YY)= 5.* YY
5192
5193 ! ----------------------------------------------------------------------
5194 ! ZTFC: RATIO OF ZOH/ZOM LESS OR EQUAL THAN 1
5195 ! C......ZTFC=0.1
5196 ! CZIL: CONSTANT C IN Zilitinkevich, S. S.1995,:NOTE ABOUT ZT
5197 ! ----------------------------------------------------------------------
5198 ILECH = 0
5199
5200 ! ----------------------------------------------------------------------
5201 ZILFC = - CZIL * VKRM * SQVISC
5202 ! C.......ZT=Z0*ZTFC
5203 ZU = Z0
5204 RDZ = 1./ ZLM
5205 CXCH = EXCM * RDZ
5206 DTHV = THLM - THZ0
5207
5208 ! ----------------------------------------------------------------------
5209 ! BELJARS CORRECTION OF USTAR
5210 ! ----------------------------------------------------------------------
5211 DU2 = MAX (SFCSPD * SFCSPD,EPSU2)
5212 !cc If statements to avoid TANGENT LINEAR problems near zero
5213 BTGH = BTG * HPBL
5214 IF (BTGH * AKHS * DTHV .ne. 0.0) THEN
5215 WSTAR2 = WWST2* ABS (BTGH * AKHS * DTHV)** (2./3.)
5216 ELSE
5217 WSTAR2 = 0.0
5218 END IF
5219
5220 ! ----------------------------------------------------------------------
5221 ! ZILITINKEVITCH APPROACH FOR ZT
5222 ! ----------------------------------------------------------------------
5223 USTAR = MAX (SQRT (AKMS * SQRT (DU2+ WSTAR2)),EPSUST)
5224
5225 ! ----------------------------------------------------------------------
5226 ZT = EXP (ZILFC * SQRT (USTAR * Z0))* Z0
5227 ZSLU = ZLM + ZU
5228 ! PRINT*,'ZSLT=',ZSLT
5229 ! PRINT*,'ZLM=',ZLM
5230 ! PRINT*,'ZT=',ZT
5231
5232 ZSLT = ZLM + ZT
5233 RLOGU = log (ZSLU / ZU)
5234
5235 RLOGT = log (ZSLT / ZT)
5236 ! PRINT*,'RLMO=',RLMO
5237 ! PRINT*,'ELFC=',ELFC
5238 ! PRINT*,'AKHS=',AKHS
5239 ! PRINT*,'DTHV=',DTHV
5240 ! PRINT*,'USTAR=',USTAR
5241
5242 RLMO = ELFC * AKHS * DTHV / USTAR **3
5243 ! ----------------------------------------------------------------------
5244 ! 1./MONIN-OBUKKHOV LENGTH-SCALE
5245 ! ----------------------------------------------------------------------
5246 DO ITR = 1,ITRMX
5247 ZETALT = MAX (ZSLT * RLMO,ZTMIN)
5248 RLMO = ZETALT / ZSLT
5249 ZETALU = ZSLU * RLMO
5250 ZETAU = ZU * RLMO
5251
5252 ZETAT = ZT * RLMO
5253 IF (ILECH .eq. 0) THEN
5254 IF (RLMO .lt. 0.)THEN
5255 XLU4 = 1. -16.* ZETALU
5256 XLT4 = 1. -16.* ZETALT
5257 XU4 = 1. -16.* ZETAU
5258
5259 XT4 = 1. -16.* ZETAT
5260 XLU = SQRT (SQRT (XLU4))
5261 XLT = SQRT (SQRT (XLT4))
5262 XU = SQRT (SQRT (XU4))
5263
5264 XT = SQRT (SQRT (XT4))
5265 ! PRINT*,'-----------1------------'
5266 ! PRINT*,'PSMZ=',PSMZ
5267 ! PRINT*,'PSPMU(ZETAU)=',PSPMU(ZETAU)
5268 ! PRINT*,'XU=',XU
5269 ! PRINT*,'------------------------'
5270 PSMZ = PSPMU (XU)
5271 SIMM = PSPMU (XLU) - PSMZ + RLOGU
5272 PSHZ = PSPHU (XT)
5273 SIMH = PSPHU (XLT) - PSHZ + RLOGT
5274 ELSE
5275 ZETALU = MIN (ZETALU,ZTMAX)
5276 ZETALT = MIN (ZETALT,ZTMAX)
5277 ! PRINT*,'-----------2------------'
5278 ! PRINT*,'PSMZ=',PSMZ
5279 ! PRINT*,'PSPMS(ZETAU)=',PSPMS(ZETAU)
5280 ! PRINT*,'ZETAU=',ZETAU
5281 ! PRINT*,'------------------------'
5282 PSMZ = PSPMS (ZETAU)
5283 SIMM = PSPMS (ZETALU) - PSMZ + RLOGU
5284 PSHZ = PSPHS (ZETAT)
5285 SIMH = PSPHS (ZETALT) - PSHZ + RLOGT
5286 END IF
5287 ! ----------------------------------------------------------------------
5288 ! LECH'S FUNCTIONS
5289 ! ----------------------------------------------------------------------
5290 ELSE
5291 IF (RLMO .lt. 0.)THEN
5292 ! PRINT*,'-----------3------------'
5293 ! PRINT*,'PSMZ=',PSMZ
5294 ! PRINT*,'PSLMU(ZETAU)=',PSLMU(ZETAU)
5295 ! PRINT*,'ZETAU=',ZETAU
5296 ! PRINT*,'------------------------'
5297 PSMZ = PSLMU (ZETAU)
5298 SIMM = PSLMU (ZETALU) - PSMZ + RLOGU
5299 PSHZ = PSLHU (ZETAT)
5300 SIMH = PSLHU (ZETALT) - PSHZ + RLOGT
5301 ELSE
5302 ZETALU = MIN (ZETALU,ZTMAX)
5303
5304 ZETALT = MIN (ZETALT,ZTMAX)
5305 ! PRINT*,'-----------4------------'
5306 ! PRINT*,'PSMZ=',PSMZ
5307 ! PRINT*,'PSLMS(ZETAU)=',PSLMS(ZETAU)
5308 ! PRINT*,'ZETAU=',ZETAU
5309 ! PRINT*,'------------------------'
5310 PSMZ = PSLMS (ZETAU)
5311 SIMM = PSLMS (ZETALU) - PSMZ + RLOGU
5312 PSHZ = PSLHS (ZETAT)
5313 SIMH = PSLHS (ZETALT) - PSHZ + RLOGT
5314 END IF
5315 ! ----------------------------------------------------------------------
5316 ! BELJAARS CORRECTION FOR USTAR
5317 ! ----------------------------------------------------------------------
5318 END IF
5319
5320 ! ----------------------------------------------------------------------
5321 ! ZILITINKEVITCH FIX FOR ZT
5322 ! ----------------------------------------------------------------------
5323 USTAR = MAX (SQRT (AKMS * SQRT (DU2+ WSTAR2)),EPSUST)
5324
5325 ZT = EXP (ZILFC * SQRT (USTAR * Z0))* Z0
5326 ZSLT = ZLM + ZT
5327 !-----------------------------------------------------------------------
5328 RLOGT = log (ZSLT / ZT)
5329 USTARK = USTAR * VKRM
5330 AKMS = MAX (USTARK / SIMM,CXCH)
5331 !-----------------------------------------------------------------------
5332 ! IF STATEMENTS TO AVOID TANGENT LINEAR PROBLEMS NEAR ZERO
5333 !-----------------------------------------------------------------------
5334 AKHS = MAX (USTARK / SIMH,CXCH)
5335 IF (BTGH * AKHS * DTHV .ne. 0.0) THEN
5336 WSTAR2 = WWST2* ABS (BTGH * AKHS * DTHV)** (2./3.)
5337 ELSE
5338 WSTAR2 = 0.0
5339 END IF
5340 !-----------------------------------------------------------------------
5341 RLMN = ELFC * AKHS * DTHV / USTAR **3
5342 !-----------------------------------------------------------------------
5343 ! IF(ABS((RLMN-RLMO)/RLMA).LT.EPSIT) GO TO 110
5344 !-----------------------------------------------------------------------
5345 RLMA = RLMO * WOLD+ RLMN * WNEW
5346 !-----------------------------------------------------------------------
5347 RLMO = RLMA
5348 ! PRINT*,'----------------------------'
5349 ! PRINT*,'SFCDIF OUTPUT ! ! ! ! ! ! ! ! ! ! ! !'
5350
5351 ! PRINT*,'ZLM=',ZLM
5352 ! PRINT*,'Z0=',Z0
5353 ! PRINT*,'THZ0=',THZ0
5354 ! PRINT*,'THLM=',THLM
5355 ! PRINT*,'SFCSPD=',SFCSPD
5356 ! PRINT*,'CZIL=',CZIL
5357 ! PRINT*,'AKMS=',AKMS
5358 ! PRINT*,'AKHS=',AKHS
5359 ! PRINT*,'----------------------------'
5360
5361 END DO
5362 ! ----------------------------------------------------------------------
5363 END SUBROUTINE SFCDIF_off
5364 ! ----------------------------------------------------------------------
5365
5366 END MODULE module_sf_noahlsm