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