module_sf_lsm_nmm.F
References to this file elsewhere.
1 !WRF:MODEL_LAYER:PHYSICS
2 !
3 MODULE MODULE_SF_LSM_NMM
4
5 USE MODULE_MPP
6 USE MODULE_MODEL_CONSTANTS
7
8 REAL, SAVE :: SCFX(30)
9
10 INTEGER, SAVE :: ISEASON
11 CHARACTER*256 :: errmess
12
13 CONTAINS
14
15 !-----------------------------------------------------------------------
16 SUBROUTINE NMMLSM(DZ8W,QV3D,P8W3D,RHO3D, &
17 & T3D,TH3D,TSK,CHS, &
18 & HFX,QFX,QGH,GSW,GLW,ELFLX,RMOL, & ! added for WRF CHEM
19 & SMSTAV,SMSTOT,SFCRUNOFF, &
20 & UDRUNOFF,IVGTYP,ISLTYP,VEGFRA,SFCEVP,POTEVP, &
21 & GRDFLX,SFCEXC,ACSNOW,ACSNOM,SNOPCX, &
22 & ALBSF,TMN,XLAND,XICE,QZ0, &
23 & TH2,Q2,SNOWC,CHS2,QSFC,TBOT,CHKLOWQ,RAINBL, &
24 & NUM_SOIL_LAYERS,DT,DZS,ITIMESTEP, &
25 & SMOIS,TSLB,SNOW,CANWAT,CPM,ROVCP,SR, &
26 & ALB,SNOALB,SMLIQ,SNOWH, &
27 & IDS,IDE, JDS,JDE, KDS,KDE, &
28 & IMS,IME, JMS,JME, KMS,KME, &
29 & ITS,ITE, JTS,JTE, KTS,KTE )
30 !-----------------------------------------------------------------------
31 !-----------------------------------------------------------------------
32 IMPLICIT NONE
33 !-----------------------------------------------------------------------
34 !-----------------------------------------------------------------------
35 !-- DZ8W thickness of layers (m)
36 !-- T3D temperature (K)
37 !-- QV3D 3D water vapor mixing ratio (Kg/Kg)
38 !-- P8W3D 3D pressure on layer interfaces (Pa)
39 !-- FLHC exchange coefficient for heat (m/s)
40 !-- FLQC exchange coefficient for moisture (m/s)
41 !-- PSFC surface pressure (Pa)
42 !-- XLAND land mask (1 for land, 2 for water)
43 !-- TMN soil temperature at lower boundary (K)
44 !-- HFX upward heat flux at the surface (W/m^2)
45 !-- QFX upward moisture flux at the surface (kg/m^2/s)
46 !-- TSK surface temperature (K)
47 !-- GSW NET downward short wave flux at ground surface (W/m^2)
48 !-- GLW downward long wave flux at ground surface (W/m^2)
49 !-- ELFLX actual latent heat flux (w m-2: positive, if up from surface)
50 !-- SFCEVP accumulated surface evaporation (W/m^2)
51 !-- POTEVP accumulated potential evaporation (W/m^2)
52 !-- CAPG heat capacity for soil (J/K/m^3)
53 !-- THC thermal inertia (Cal/cm/K/s^0.5)
54 !-- TBOT bottom soil temperature (local yearly-mean sfc air temperature)
55 !-- SNOWC flag indicating snow coverage (1 for snow cover)
56 !-- EMISS surface emissivity (between 0 and 1)
57 !-- DELTSM time step (second)
58 !-- ROVCP R/CP
59 !-- SR fraction of frozen precip (0.0 to 1.0)
60 !-- XLV latent heat of melting (J/kg)
61 !-- DTMIN time step (minute)
62 !-- IFSNOW ifsnow=1 for snow-cover effects
63 !-- SVP1 constant for saturation vapor pressure (kPa)
64 !-- SVP2 constant for saturation vapor pressure (dimensionless)
65 !-- SVP3 constant for saturation vapor pressure (K)
66 !-- SVPT0 constant for saturation vapor pressure (K)
67 !-- EP1 constant for virtual temperature (R_v/R_d - 1) (dimensionless)
68 !-- EP2 constant for specific humidity calculation
69 ! (R_d/R_v) (dimensionless)
70 !-- KARMAN Von Karman constant
71 !-- EOMEG angular velocity of earth's rotation (rad/s)
72 !-- STBOLT Stefan-Boltzmann constant (W/m^2/K^4)
73 !-- STEM soil temperature in 5-layer model
74 !-- ZS depths of centers of soil layers
75 !-- DZS thicknesses of soil layers
76 !-- num_soil_layers the number of soil layers
77 !-- ACSNOW accumulated snowfall (water equivalent) (mm)
78 !-- ACSNOM accumulated snowmelt (water equivalent) (mm)
79 !-- SNOPCX snow phase change heat flux (W/m^2)
80 !-- ids start index for i in domain
81 !-- ide end index for i in domain
82 !-- jds start index for j in domain
83 !-- jde end index for j in domain
84 !-- kds start index for k in domain
85 !-- kde end index for k in domain
86 !-- ims start index for i in memory
87 !-- ime end index for i in memory
88 !-- jms start index for j in memory
89 !-- jme end index for j in memory
90 !-- kms start index for k in memory
91 !-- kme end index for k in memory
92 !-- its start index for i in tile
93 !-- ite end index for i in tile
94 !-- jts start index for j in tile
95 !-- jte end index for j in tile
96 !-- kts start index for k in tile
97 !-- kte end index for k in tile
98 !-----------------------------------------------------------------------
99 INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE, &
100 & IMS,IME,JMS,JME,KMS,KME, &
101 & ITS,ITE,JTS,JTE,KTS,KTE
102 !
103 INTEGER,INTENT(IN) :: NUM_SOIL_LAYERS,ITIMESTEP
104 !
105 REAL,INTENT(IN) :: DT,ROVCP
106 !
107 REAL,DIMENSION(IMS:IME,1:NUM_SOIL_LAYERS,JMS:JME), &
108 & INTENT(INOUT) :: SMOIS, & ! new
109 SMLIQ, & ! new
110 TSLB !
111
112 REAL,DIMENSION(1:NUM_SOIL_LAYERS),INTENT(IN) :: DZS
113 !
114 REAL,DIMENSION(ims:ime,jms:jme),INTENT(INOUT) :: &
115 & TSK, & !was TGB (temperature)
116 & HFX, &
117 & QFX, &
118 & QSFC,&
119 & SNOW, & !new
120 & SNOWH, & !new
121 & ALB, &
122 & SNOALB, &
123 & ALBSF, &
124 & SNOWC, &
125 & CANWAT, & ! new
126 & SMSTAV, &
127 & SMSTOT, &
128 & SFCRUNOFF, &
129 & UDRUNOFF, &
130 & SFCEVP, &
131 & POTEVP, &
132 & GRDFLX, &
133 & ACSNOW, &
134 & ACSNOM, &
135 & SNOPCX, &
136 & Q2, &
137 & TH2, &
138 & SFCEXC
139
140 INTEGER,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: IVGTYP, &
141 ISLTYP
142
143 REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: TMN, &
144 XLAND, &
145 XICE, &
146 VEGFRA, &
147 GSW, &
148 GLW, &
149 QZ0, &
150 SR
151
152 REAL,DIMENSION(IMS:IME,KMS:KME,JMS:JME),INTENT(IN) :: QV3D, &
153 P8W3D, &
154 RHO3D, &
155 TH3D, &
156 T3D, &
157 DZ8W
158
159 !
160 REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: RAINBL
161 !
162 REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: CHS2, &
163 CHS, &
164 QGH, &
165 CPM
166 !
167 REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(OUT) :: TBOT
168 !
169 REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(OUT) :: CHKLOWQ, &
170 ELFLX
171 ! added for WRF-CHEM, 20041205, JM -- not used in this routine as yet
172 REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(INOUT) :: RMOL
173
174 ! LOCAL VARS
175
176 REAL,DIMENSION(ITS:ITE) :: QV1D, &
177 & T1D, &
178 & TH1D, &
179 & ZA1D, &
180 & P8W1D, &
181 & PSFC1D, &
182 & RHO1D, &
183 & PREC1D
184
185 INTEGER :: I,J
186 REAL :: RATIOMX
187 !-----------------------------------------------------------------------
188 !-----------------------------------------------------------------------
189
190 DO J=JTS,JTE
191
192 DO I=ITS,ITE
193 T1D(I) = T3D(I,1,J)
194 TH1D(I) = TH3D(I,1,J)
195 !!! QV1D(I) = QV3D(I,1,J)
196 RATIOMX = QV3D(I,1,J)
197 QV1D(I) = RATIOMX/(1.+RATIOMX)
198 P8W1D(I) = (P8W3D(I,KTS+1,j)+P8W3D(i,KTS,j))*0.5
199 PSFC1D(I) = P8W3D(I,1,J)
200 ZA1D(I) = 0.5*DZ8W(I,1,J)
201 RHO1D(I) = RHO3D(I,1,J)
202 PREC1D(I) = RAINBL(I,J)/DT
203 ENDDO
204
205 !FLHC = SFCEXC
206
207 !-----------------------------------------------------------------------
208 CALL SURFCE(J,ZA1D,QV1D,P8W1D,PSFC1D,RHO1D,T1D,TH1D,TSK, &
209 CHS(IMS,J),PREC1D,HFX,QFX,QGH(IMS,J),GSW,GLW, &
210 SMSTAV,SMSTOT,SFCRUNOFF, &
211 UDRUNOFF,IVGTYP,ISLTYP,VEGFRA,SFCEVP,POTEVP,GRDFLX, &
212 ELFLX,SFCEXC,ACSNOW,ACSNOM,SNOPCX, &
213 ALBSF,TMN,XLAND,XICE,QZ0, &
214 TH2,Q2,SNOWC,CHS2(IMS,J),QSFC,TBOT,CHKLOWQ, &
215 NUM_SOIL_LAYERS,DT,DZS,ITIMESTEP, &
216 SMOIS,TSLB,SNOW,CANWAT,CPM(IMS,J),ROVCP,SR, &
217 ALB,SNOALB,SMLIQ,SNOWH, &
218 IMS,IME,JMS,JME,KMS,KME, &
219 ITS,ITE,JTS,JTE,KTS,KTE )
220 !
221 ENDDO
222
223 END SUBROUTINE NMMLSM
224
225 !-----------------------------------------------------------------------
226 !-----------------------------------------------------------------------
227 SUBROUTINE SURFCE(J,ZA,QV,P8W,PSFC,RHO,T,TH,TSK,CHS,PREC,HFX,QFX, &
228 QGH,GSW,GLW,SMSTAV,SMSTOT,SFCRUNOFF,UDRUNOFF, &
229 IVGTYP,ISLTYP,VEGFRA,SFCEVP,POTEVP,GRDFLX, &
230 ELFLX,SFCEXC,ACSNOW,ACSNOM,SNOPCX, &
231 ALBSF,TMN,XLAND,XICE,QZ0, &
232 TH2,Q2,SNOWC,CHS2,QSFC,TBOT,CHKLOWQ, &
233 NUM_SOIL_LAYERS,DT,DZS,ITIMESTEP, &
234 SMOIS,TSLB,SNOW,CANWAT,CPM,ROVCP,SR, &
235 ALB,SNOALB,SMLIQ,SNOWH, &
236 IMS,IME,JMS,JME,KMS,KME, &
237 ITS,ITE,JTS,JTE,KTS,KTE )
238 !------------------------------------------------------------------------
239 IMPLICIT NONE
240 !------------------------------------------------------------------------
241 !$$$ SUBPROGRAM DOCUMENTATION BLOCK
242 ! . . .
243 ! SUBPROGRAM: SURFCE CALCULATE SURFACE CONDITIONS
244 ! PRGRMMR: F. CHEN DATE: 97-12-06
245 !
246 ! ABSTRACT:
247 ! THIS ROUTINE IS THE DRIVER FOR COMPUTATION OF GROUND CONDITIONS
248 ! BY USING A LAND SURFACE MODEL (LSM).
249 !
250 ! PROGRAM HISTORY LOG:
251 ! 97-12-06 CHEN - ORIGINATOR
252 !
253 ! REFERENCES:
254 ! PAN AND MAHRT (1987) BOUN. LAYER METEOR.
255 ! CHEN ET AL. (1996) J. GEOPHYS. RES.
256 ! CHEN ET AL. (1997) BOUN. LAYER METEOR.
257 ! CHEN and Dudhia (2000) Mon. Wea. Rev.
258 !
259 ! SUBPROGRAMS CALLED:
260 ! SFLX
261 !
262 ! SET LOCAL PARAMETERS.
263 !----------------------------------------------------------------------
264 INTEGER, INTENT(IN ) :: IMS,IME, JMS,JME, KMS,KME, &
265 ITS,ITE, JTS,JTE, KTS,KTE, &
266 J,ITIMESTEP
267
268 INTEGER , INTENT(IN) :: NUM_SOIL_LAYERS
269
270 REAL, INTENT(IN ) :: DT,ROVCP
271
272 REAL, DIMENSION(1:num_soil_layers), INTENT(IN)::DZS
273
274
275 REAL, PARAMETER :: PQ0=379.90516
276 REAL, PARAMETER :: TRESH=.95E0,A2=17.2693882,A3=273.16,A4=35.86, &
277 T0=273.16E0,T1=274.16E0,ROW=1.E3, &
278 ELWV=2.50E6,ELIV=XLS,ELIW=XLF, &
279 A23M4=A2*(A3-A4), RLIVWV=ELIV/ELWV, &
280 ROWLIW=ROW*ELIW,ROWLIV=ROW*ELIV,CAPA=R_D/CP
281
282 INTEGER, PARAMETER :: NROOT=3
283 !
284 REAL, DIMENSION( ims:ime , 1:num_soil_layers, jms:jme ), &
285 INTENT(INOUT) :: SMOIS, & ! new
286 SMLIQ, & ! new
287 TSLB ! new !STEMP
288
289
290 REAL, DIMENSION( ims:ime, jms:jme ) , &
291 INTENT(INOUT) :: TSK, & !was TGB (temperature)
292 HFX, & !new
293 QFX, & !new
294 QSFC,& !new
295 SNOW, & !new
296 SNOWH, & !new
297 ALB, &
298 SNOALB, &
299 ALBSF, &
300 SNOWC, &
301 CANWAT, & ! new
302 SMSTAV, &
303 SMSTOT, &
304 SFCRUNOFF, &
305 UDRUNOFF, &
306 SFCEVP, &
307 POTEVP, &
308 GRDFLX, &
309 ACSNOW, &
310 ACSNOM, &
311 SNOPCX
312
313 INTEGER, DIMENSION( ims:ime, jms:jme ) , &
314 INTENT(IN ) :: IVGTYP, &
315 ISLTYP
316
317 REAL, DIMENSION( ims:ime, jms:jme ) , &
318 INTENT(IN ) :: TMN, &
319 XLAND, &
320 XICE, &
321 VEGFRA, &
322 GSW, &
323 GLW, &
324 QZ0, &
325 SR
326
327 REAL, DIMENSION( ims:ime, jms:jme ) , &
328 INTENT(INOUT) :: Q2, &
329 TH2, &
330 SFCEXC
331
332 REAL, DIMENSION( ims:ime, jms:jme ) , &
333 INTENT(OUT) :: TBOT
334
335
336 REAL, DIMENSION( ims:ime, jms:jme ) , &
337 INTENT(OUT) :: CHKLOWQ, &
338 ELFLX
339
340 REAL, DIMENSION( ims:ime ) , &
341 INTENT(IN ) :: QGH, &
342 CHS, &
343 CPM, &
344 CHS2
345
346 ! MODULE-LOCAL VARIABLES, DEFINED IN SUBROUTINE LSM
347 REAL, DIMENSION( its:ite ) , &
348 INTENT(IN ) :: ZA, &
349 TH, &
350 QV, &
351 T, &
352 p8w, &
353 PSFC, &
354 rho, &
355 PREC ! one time step in mm
356
357 REAL, DIMENSION( its:ite ) :: TGDSA
358
359 ! LOCAL VARS
360
361 REAL, DIMENSION(1:num_soil_layers) :: SMLIQ1D,SMOIS1D,STEMP1D
362
363 !----------------------------------------------------------------------
364 !*** DECLARATIONS FOR IMPLICIT NONE
365
366 REAL :: APELM,APES,FDTLIW,FDTW,Q2SAT,Z,FK,SOLDN,SFCTMP,SFCTH2, &
367 SFCPRS,PRCP,Q2K,DQSDTK,SATFLG,TBOTK,CHK,VGFRCK,T1K,LWDN, &
368 CMCK,Q2M,SNODPK,PLFLX,HFLX,GFLX,RNOF1K, &
369 RNOF2K,Q1K,SMELTK,SOILQW,SOILQM,T2K,PRESK,CHFF,STIMESTEP, &
370 ALB1D,SNOALB1D,SNOWH1D,ALBSF1D,SOLNET,FFROZP, &
371 DUM1,DUM2,DUM3,DUM4,DUM5,DUM6,DUM7
372
373 INTEGER :: I,K,NS,ICE,IVGTPK,ISLTPK,ISPTPK,NOOUT,NSOIL,LZ
374
375 !----------------------------------------------------------------------
376 !***********************************************************************
377 ! START SURFCE HERE
378 !***
379 !*** SET CONSTANTS CALCULATED HERE FOR CLARITY.
380 !***
381 FDTLIW=DT/ROWLIW
382 ! FDTLIV=DT/ROWLIV
383 FDTW=DT/(XLV*RHOWATER)
384 !***
385 !*** SET LSM CONSTANTS AND TIME INDEPENDENT VARIABLES
386 !*** INITIALIZE LSM HISTORICAL VARIABLES
387 !***
388 !-----------------------------------------------------------------------
389
390 NSOIL=num_soil_layers
391
392 IF(ITIMESTEP.EQ.1)THEN
393 DO 50 I=its,ite
394 !*** SET ZERO-VALUE FOR SOME OUTPUT DIAGNOSTIC ARRAYS
395 IF((XLAND(I,J)-1.5).GE.0.)THEN
396 ! check sea-ice point
397 IF(XICE(I,J).EQ.1.)PRINT*,' sea-ice at water point, I=',I, &
398 'J=',J
399 !*** Open Water Case
400 SMSTAV(I,J)=1.0
401 SMSTOT(I,J)=1.0
402 DO NS=1,NSOIL
403 SMOIS(I,NS,J)=1.0
404 TSLB(I,NS,J)=273.16 !STEMP
405 ENDDO
406 ELSE
407 IF(XICE(I,J).EQ.1.)THEN
408 !*** SEA-ICE CASE
409 SMSTAV(I,J)=1.0
410 SMSTOT(I,J)=1.0
411 DO NS=1,NSOIL
412 SMOIS(I,NS,J)=1.0
413 ENDDO
414 ENDIF
415 ENDIF
416 !
417 50 CONTINUE
418 ENDIF
419 !-----------------------------------------------------------------------
420 DO 100 I=its,ite
421 ! SFCPRS=(A(KL)*PSB(I,J)+PTOP+PP3D(I,J,KL)*0.001)*1.E3
422 SFCPRS=p8w(I) !Pressure in middle of lowest layer
423 Q2SAT=QGH(I)
424 ! CHKLOWQ(I,J)=1.
425 CHFF=CHS(I)*RHO(I)*CPM(I)
426 !CHK*RHO*CP
427 ! TGDSA: potential T
428 TGDSA(I)=TSK(I,J)*(1.E5/SFCPRS)**ROVCP
429 !
430 !*** CHECK FOR SATURATION AT THE LOWEST MODEL LEVEL
431 !
432 Q2K=QV(I)
433 APES=(1.E5/PSFC(I))**CAPA
434 !
435 IF((Q2K.GE.Q2SAT*TRESH).AND.Q2K.LT.QZ0(I,J))THEN
436 SATFLG=0.
437 CHKLOWQ(I,J)=0.
438 ELSE
439 SATFLG=1.0
440 CHKLOWQ(I,J)=1.
441 ENDIF
442 !
443 TBOT(I,J)=273.16
444 !***
445 !*** LOADING AND UNLOADING MM5/LSM LAND SOIL VARIABLES
446 !***
447 IF((XLAND(I,J)-1.5).GE.0.)THEN
448 !*** Water
449 HFX(I,J)=HFX(I,J)/APES
450 QFX(I,J)=QFX(I,J)*SATFLG
451 SFCEVP(I,J)=SFCEVP(I,J)+QFX(I,J)*DT
452 ELSE
453 !*** LAND OR SEA-ICE
454 !ATEC ICE=INT(XICE(I,J)+0.3)
455 IF (XICE(I,J) .GT. 0.5) THEN
456 ICE=1
457 ELSE
458 ICE=0
459 ENDIF
460 !
461 Q2K=MIN(QV(I),Q2SAT)
462 Z=ZA(I)
463 ! FK=GSW(I,J)+GLW(I,J)
464 LWDN=GLW(I,J)
465 !
466 !*** GSW is net downward shortwave
467 !
468 ! SOLNET=GSW(I,J)
469 !
470 !*** GSW is total downward shortwave
471 !
472 SOLDN=GSW(I,J)
473 !
474 !*** Simple use of albedo to determine total incoming solar shortwave SOLDN
475 !*** (no solar zenith angle correction)
476 !
477 ! SOLDN=SOLNET/(1.-ALB(I,J))
478 SOLNET=SOLDN*(1.-ALB(I,J))
479 !
480 ALBSF1D=ALBSF(I,J)
481 SNOALB1D=SNOALB(I,J)
482 SFCTMP=T(I)
483 !!! SFCTH2=SFCTMP+(0.0097545*Z)
484 APELM=(1.E5/SFCPRS)**CAPA
485 SFCTH2=SFCTMP*APELM
486 SFCTH2=SFCTH2/APES
487 PRCP=PREC(I)
488 !!! Q2K=QV(I)
489 !!! Q2SAT=PQ0/SFCPRS*EXP(A2*(SFCTMP-A3)/(SFCTMP-A4))
490 DQSDTK=Q2SAT*A23M4/(SFCTMP-A4)**2
491 IF(ICE.EQ.0)THEN
492 TBOTK=TMN(I,J)
493 ELSE
494 TBOTK=271.16
495 ENDIF
496 CHK=CHS(I)
497 IVGTPK=IVGTYP(I,J)
498 IF(IVGTPK.EQ.0)IVGTPK=13
499 ISLTPK=ISLTYP(I,J)
500 IF(ISLTPK.EQ.0)ISLTPK=9
501 ! hardwire slope type (ISPTPK)=1
502 ISPTPK=1
503 VGFRCK=VEGFRA(I,J)/100.
504 IF(IVGTPK.EQ.25) VGFRCK=0.0001
505 IF(ISLTPK.EQ.14.AND.XICE(I,J).EQ.0.)THEN
506 PRINT*,' SOIL TYPE FOUND TO BE WATER AT A LAND-POINT'
507 PRINT*,i,j,'RESET SOIL in surfce.F'
508 ! ISLTYP(I,J)=7
509 ISLTPK=7
510 ENDIF
511 T1K=TSK(I,J)
512 CMCK=CANWAT(I,J)
513 !*** convert snow depth from mm to meter
514 SNODPK=SNOW(I,J)*0.001
515 SNOWH1D=SNOWH(I,J)*0.001
516 !
517 !*** fraction of frozen precip
518 !
519 FFROZP=SR(I,J)
520 !
521 DO 70 NS=1,NSOIL
522 SMOIS1D(NS)=SMOIS(I,NS,J)
523 SMLIQ1D(NS)=SMLIQ(I,NS,J)
524 STEMP1D(NS)=TSLB(I,NS,J) !STEMP
525 70 CONTINUE
526
527 !
528 ! print*,'BF SFLX','ISLTPK',ISLTPK,'IVGTPK=',IVGTPK,'SMOIS1D',&
529 ! SMOIS1D,'STEMP1',STEMP1D,'VGFRCK',VGFRCK
530 !-----------------------------------------------------------------------
531 ! old WRF call to SFLX
532 ! CALL SFLX(ICE,SATFLG,DT,Z,NSOIL,NROOT,DZS,FK,SOLDN,SFCPRS, &
533 ! PRCP,SFCTMP,SFCTH2,Q2K,Q2SAT,DQSDTK,TBOTK,CHK,CHFF, &
534 ! IVGTPK,ISLTPK,VGFRCK,PLFLX,ELFLX,HFLX,GFLX,RNOF1K,RNOF2K,&
535 ! Q1K,SMELTK,T1K,CMCK,SMOIS1D,STEMP1D,SNODPK,SOILQW,SOILQM)
536 !-----------------------------------------------------------------------
537 ! ----------------------------------------------------------------------
538 ! Ek 12 June 2002 - NEW CALL SFLX
539 ! ops Eta call to SFLX ...'tailor' this to WRF
540 ! CALL SFLX
541 ! I (ICE,DTK,Z,NSOIL,SLDPTH,
542 ! I LWDN,SOLDN,SFCPRS,PRCP,SFCTMP,SFCTH2,Q2K,SFCSPD,Q2SAT,DQSDTK,
543 ! I IVGTPK,ISLTPK,ISPTPK,
544 ! I VGFRCK,PTU,TBOT,ALB,SNOALB,
545 ! 2 CMCK,T1K,STCK,SMCK,SH2OK,SNOWH,SNODPK,ALB2D,CHK,CMK,
546 ! O PLFLX,ELFLX,HFLX,GFLX,RNOF1K,RNOF2K,Q1K,SMELTK,
547 ! O SOILQW,SOILQM,DUM1,DUM2,DUM3,DUM4)
548 !-----------------------------------------------------------------------
549 CALL SFLX &
550 (FFROZP,ICE,DT,Z,NSOIL,DZS, &
551 LWDN,SOLDN,SOLNET,SFCPRS,PRCP,SFCTMP,SFCTH2,Q2K,DUM5,Q2SAT, &
552 DQSDTK,IVGTPK,ISLTPK,ISPTPK, &
553 VGFRCK,DUM6,TBOTK,ALBSF1D,SNOALB1D, &
554 CMCK,T1K,STEMP1D,SMOIS1D,SMLIQ1D,SNOWH1D,SNODPK,ALB1D,CHK,DUM7, &
555 PLFLX,ELFLX(I,J),HFLX,GFLX,RNOF1K,RNOF2K,Q1K,SMELTK, &
556 SOILQW,SOILQM,DUM1,DUM2,DUM3,DUM4,I,J)
557 !-----------------------------------------------------------------------
558 !*** DIAGNOSTICS
559 ! Convert the water unit into mm
560 SFCRUNOFF(I,J)=SFCRUNOFF(I,J)+RNOF1K*DT*1000.0
561 UDRUNOFF(I,J)=UDRUNOFF(I,J)+RNOF2K*DT*1000.0
562 SMSTAV(I,J)=SOILQW
563
564 !mp
565 if (abs(SMSTAV(I,J)) .lt. 3.5) then
566 else
567 write(errmess,*) 'bad SMSTAV: ', I,J,SMSTAV(I,J)
568 CALL wrf_message( errmess )
569 endif
570 !mp
571
572 SMSTOT(I,J)=SOILQM*1000.
573 SFCEXC(I,J)=CHK
574 ! IF(SNOB(I,J).GT.0..OR.SICE(I,J).GT.0.)THEN
575 ! QFC1(I,J)=QFC1(I,J)*RLIVWV
576 ! ENDIF
577 IF(FFROZP.GT.0.5)THEN
578 ACSNOW(I,J)=ACSNOW(I,J)+PREC(I)*DT
579 ENDIF
580 IF(SNOW(I,J).GT.0.)THEN
581 ACSNOM(I,J)=ACSNOM(I,J)+SMELTK*1000.
582 SNOPCX(I,J)=SNOPCX(I,J)-SMELTK/FDTLIW
583 ENDIF
584 POTEVP(I,J)=POTEVP(I,J)+PLFLX*FDTW
585 ! POTFLX(I,J)=POTFLX(I,J)-PLFLX
586 !*** WRF LOWER BOUNDARY CONDITIONS
587 GRDFLX(I,J)=GFLX
588 HFX(I,J)=HFLX
589 QFX(I,J)=ELFLX(I,J)/ELWV
590 SFCEVP(I,J)=SFCEVP(I,J)+QFX(I,J)*DT
591 TSK(I,J)=T1K
592 T2K=T1K-HFX(I,J)/(RHO(I)*CPM(I)*CHS2(I))
593 TH2(I,J)=T2K*(1.E5/SFCPRS)**ROVCP
594 Q2M=Q1K-QFX(I,J)/(RHO(I)*CHS2(I))
595 !!!!!! Q2(I,J)=Q2M
596 !!!!!! Q2(I,J)=Q2K
597 ! t2k=th2k/(1.E5/SFCPRS)**ROVCP
598 ! QS(I,J)=Q1K
599 !!! QSFC(I,J)=Q1K
600 !*** UPDATE STATE VARIABLES
601 SNOW(I,J)=SNODPK*1000.0
602 SNOWH(I,J)=SNOWH1D*1000.0
603 CANWAT(I,J)=CMCK
604 IF(SNOW(I,J).GT.1.0)THEN
605 ! ALB(I,J)=0.01*ALBD(IVGTPK,ISEASON)*(1.+SCFX(IVGTPK))
606 SNOWC(I,J)=1.0
607 ELSE
608 ! ALB(I,J)=0.01*ALBD(IVGTPK,ISEASON)
609 SNOWC(I,J)=0.0
610 ENDIF
611 ! update albedo
612 ALB(I,J)=ALB1D
613 ! update bottom soil temperature
614 TBOT(I,J)=TBOTK
615
616 DO 80 NS=1,NSOIL
617 SMOIS(I,NS,J)=SMOIS1D(NS)
618 SMLIQ(I,NS,J)=SMLIQ1D(NS)
619 TSLB(I,NS,J)=STEMP1D(NS) ! STEMP
620 80 CONTINUE
621 ENDIF
622 #if 0
623 NOOUT=0
624
625 IF((ITIMESTEP.EQ.1.OR.MOD(ITIMESTEP,1).EQ.0) &
626 .AND. I .EQ.29.AND.J.EQ.23) THEN
627 ! print*, 'GLW',GLW(I,J),'GSW',GSW(I,J)
628 print*, 'T2K',T2K,'T1K',T1K,'HFX',HFX(I,J),'RHO',RHO(I),'CPM',CPM(I), &
629 'CHS2',CHS2(I),'soil T',STEMP1D,'soil m', SMOIS1D
630 ! print*,'Q2M',Q2M,'Q1K',Q1K,'QFX',QFX(I,J),'RHO',RHO(I),'CHS2',CHS2(I),'latent',ELFLX
631 ENDIF
632
633 IF(NOOUT.EQ.1)GOTO 100
634 ! write output to 29
635 IF(ITIMESTEP.EQ.1.AND.I.EQ.1.AND.J.EQ.1) &
636 WRITE (29,*)&
637 'itimestep ',' FK ',' SOLDN ',' SFCPR ', &
638 ' SFCTMP ',' Q2K ',' TBOTK ', &
639 ' CHK ',' ELFLX ',' HFLX ',' GFLX ', &
640 ' RNOF1K ',' RNOF2K ',' T1K ',' CMCK ', &
641 ' SMCK1 ',' SMCK2 ',' SMCK3 ',' SMCK4 ', &
642 ' STCK1 ',' STCK2 ',' STCK3 ',' STCK4 ', &
643 ' SNODPK ',' T2 ', &
644 ' Q2 ',' SMSTOT ',' SFCEVP ', ' RAIN'
645 IF((ITIMESTEP.EQ.1.OR.MOD(ITIMESTEP,1).EQ.0) &
646 .AND. I .EQ.29.AND.J.EQ.23) THEN
647 print *,'outputting at itimestep =', itimestep
648 STIMESTEP=FLOAT(itimestep)
649 WRITE (29,1029)STIMESTEP,FK,SOLDN,SFCPRS/100.,SFCTMP,1000.* &
650 Q2K,TBOTK,1000.*CHK,ELFLX(i,j),HFLX,GFLX,SFCRUNOFF(I,J)&
651 ,UDRUNOFF(I,J),T1K,CMCK,SMOIS1D,STEMP1D,SNODPK, &
652 ! ,UDRUNOFF(I,J),T1K,CMCK,SMOIS1D(3),SMOIS1D(7),SMOIS1D(11),&
653 ! SMOIS1D(14),STEMP1D(3), STEMP1D(7),STEMP1D(11), &
654 ! STEMP1D(14), SNODPK, &
655 T2K,1000.*Q2M,SMSTOT(I,J),SFCEVP(I,J),PRCP
656 1029 FORMAT (29F10.4)
657 ! IF(ITIMESTEP.EQ.0)WRITE (39,*)' P ',' T ', &
658 ! ' TH ',' Q ',' U ',' V ', &
659 ! ' QC '
660 ! WRITE (39,1039)itimestep
661 ! DO K=kts,kte
662 ! WRITE (39,1039)PRESK,TX(I,K),THX(I,K),1000.*QX(I,K),UX(I,K),&
663 ! VX(I,K),1000.*QCX(I,K)
664 1039 FORMAT (7F10.5)
665 ! ENDDO
666 ENDIF
667 !
668 #endif
669 100 CONTINUE
670 !
671 !-----------------------------------------------------------------------
672 END SUBROUTINE SURFCE
673 !-----------------------------------------------------------------------
674
675 SUBROUTINE SFLX ( &
676 FFROZP,ICE,DT,ZLVL,NSOIL,SLDPTH, &
677 LWDN,SOLDN,SOLNET,SFCPRS,PRCP,SFCTMP,TH2,Q2,SFCSPD,Q2SAT, &
678 DQSDT2,VEGTYP,SOILTYP,SLOPETYP, &
679 SHDFAC,PTU,TBOT,ALB,SNOALB, &
680 CMC,T1,STC,SMC,SH2O,SNOWH,SNEQV,ALBEDO,CH,CM, &
681 ETP,ETA,SHEAT,SSOIL,RUNOFF1,RUNOFF2,Q1,SNOMLT, &
682 SOILW,SOILM,SMCWLT,SMCDRY,SMCREF,SMCMAX,I,J)
683 ! ----------------------------------------------------------------------
684 ! & ETA,SHEAT, &
685 ! ----------------------------------------------------------------------
686 ! OUTPUTS, DIAGNOSTICS, PARAMETERS BELOW GENERALLY NOT NECESSARY WHEN
687 ! COUPLED WITH E.G. A NWP MODEL (SUCH AS THE NOAA/NWS/NCEP MESOSCALE ETA
688 ! MODEL). OTHER APPLICATIONS MAY REQUIRE DIFFERENT OUTPUT VARIABLES.
689 ! ----------------------------------------------------------------------
690 ! & EC,EDIR,ET,ETT,ESNOW,DRIP,DEW, &
691 ! & BETA,ETP,SSOIL, &
692 ! & FLX1,FLX2,FLX3, &
693 ! & SNOMLT,SNCOVR, &
694 ! & RUNOFF1,RUNOFF2,RUNOFF3, &
695 ! & RC,PC,RSMIN,XLAI,RCS,RCT,RCQ,RCSOIL, &
696 ! & SOILW,SOILM, &
697 ! & SMCWLT,SMCDRY,SMCREF,SMCMAX,NROOT,I,J)
698
699 IMPLICIT NONE
700
701 ! ----------------------------------------------------------------------
702 ! SUBROUTINE SFLX - VERSION 2.7 - June 2nd 2003
703 ! ----------------------------------------------------------------------
704 ! SUB-DRIVER FOR "NOAH/OSU LSM" FAMILY OF PHYSICS SUBROUTINES FOR A
705 ! SOIL/VEG/SNOWPACK LAND-SURFACE MODEL TO UPDATE SOIL MOISTURE, SOIL
706 ! ICE, SOIL TEMPERATURE, SKIN TEMPERATURE, SNOWPACK WATER CONTENT,
707 ! SNOWDEPTH, AND ALL TERMS OF THE SURFACE ENERGY BALANCE AND SURFACE
708 ! WATER BALANCE (EXCLUDING INPUT ATMOSPHERIC FORCINGS OF DOWNWARD
709 ! RADIATION AND PRECIP)
710 ! ----------------------------------------------------------------------
711 ! SFLX ARGUMENT LIST KEY:
712 ! ----------------------------------------------------------------------
713 ! C CONFIGURATION INFORMATION
714 ! F FORCING DATA
715 ! I OTHER (INPUT) FORCING DATA
716 ! S SURFACE CHARACTERISTICS
717 ! H HISTORY (STATE) VARIABLES
718 ! O OUTPUT VARIABLES
719 ! D DIAGNOSTIC OUTPUT
720 ! ----------------------------------------------------------------------
721 ! 1. CONFIGURATION INFORMATION (C):
722 ! ----------------------------------------------------------------------
723 ! ICE SEA-ICE FLAG (=1: SEA-ICE, =0: LAND)
724 ! DT TIMESTEP (SEC) (DT SHOULD NOT EXCEED 3600 SECS, RECOMMEND
725 ! 1800 SECS OR LESS)
726 ! ZLVL HEIGHT (M) ABOVE GROUND OF ATMOSPHERIC FORCING VARIABLES
727 ! NSOIL NUMBER OF SOIL LAYERS (AT LEAST 2, AND NOT GREATER THAN
728 ! PARAMETER NSOLD SET BELOW)
729 ! SLDPTH THE THICKNESS OF EACH SOIL LAYER (M)
730 ! ----------------------------------------------------------------------
731 ! 2. FORCING DATA (F):
732 ! ----------------------------------------------------------------------
733 ! LWDN LW DOWNWARD RADIATION (W M-2; POSITIVE, NOT NET LONGWAVE)
734 ! SOLDN SOLAR DOWNWARD RADIATION (W M-2; POSITIVE, NOT NET SOLAR)
735 ! SFCPRS PRESSURE AT HEIGHT ZLVL ABOVE GROUND (PASCALS)
736 ! PRCP PRECIP RATE (KG M-2 S-1) (NOTE, THIS IS A RATE)
737 ! SFCTMP AIR TEMPERATURE (K) AT HEIGHT ZLVL ABOVE GROUND
738 ! TH2 AIR POTENTIAL TEMPERATURE (K) AT HEIGHT ZLVL ABOVE GROUND
739 ! Q2 MIXING RATIO AT HEIGHT ZLVL ABOVE GROUND (KG KG-1)
740 ! ----------------------------------------------------------------------
741 ! 3. OTHER FORCING (INPUT) DATA (I):
742 ! ----------------------------------------------------------------------
743 ! SFCSPD WIND SPEED (M S-1) AT HEIGHT ZLVL ABOVE GROUND
744 ! Q2SAT SAT MIXING RATIO AT HEIGHT ZLVL ABOVE GROUND (KG KG-1)
745 ! DQSDT2 SLOPE OF SAT SPECIFIC HUMIDITY CURVE AT T=SFCTMP
746 ! (KG KG-1 K-1)
747 ! ----------------------------------------------------------------------
748 ! 4. CANOPY/SOIL CHARACTERISTICS (S):
749 ! ----------------------------------------------------------------------
750 ! VEGTYP VEGETATION TYPE (INTEGER INDEX)
751 ! SOILTYP SOIL TYPE (INTEGER INDEX)
752 ! SLOPETYP CLASS OF SFC SLOPE (INTEGER INDEX)
753 ! SHDFAC AREAL FRACTIONAL COVERAGE OF GREEN VEGETATION
754 ! (FRACTION= 0.0-1.0)
755 ! SHDMIN MINIMUM AREAL FRACTIONAL COVERAGE OF GREEN VEGETATION
756 ! (FRACTION= 0.0-1.0) <= SHDFAC
757 ! PTU PHOTO THERMAL UNIT (PLANT PHENOLOGY FOR ANNUALS/CROPS)
758 ! (NOT YET USED, BUT PASSED TO REDPRM FOR FUTURE USE IN
759 ! VEG PARMS)
760 ! ALB BACKROUND SNOW-FREE SURFACE ALBEDO (FRACTION), FOR JULIAN
761 ! DAY OF YEAR (USUALLY FROM TEMPORAL INTERPOLATION OF
762 ! MONTHLY MEAN VALUES' CALLING PROG MAY OR MAY NOT
763 ! INCLUDE DIURNAL SUN ANGLE EFFECT)
764 ! SNOALB UPPER BOUND ON MAXIMUM ALBEDO OVER DEEP SNOW (E.G. FROM
765 ! ROBINSON AND KUKLA, 1985, J. CLIM. & APPL. METEOR.)
766 ! TBOT BOTTOM SOIL TEMPERATURE (LOCAL YEARLY-MEAN SFC AIR
767 ! TEMPERATURE)
768 ! ----------------------------------------------------------------------
769 ! 5. HISTORY (STATE) VARIABLES (H):
770 ! ----------------------------------------------------------------------
771 ! CMC CANOPY MOISTURE CONTENT (M)
772 ! T1 GROUND/CANOPY/SNOWPACK) EFFECTIVE SKIN TEMPERATURE (K)
773 ! STC(NSOIL) SOIL TEMP (K)
774 ! SMC(NSOIL) TOTAL SOIL MOISTURE CONTENT (VOLUMETRIC FRACTION)
775 ! SH2O(NSOIL) UNFROZEN SOIL MOISTURE CONTENT (VOLUMETRIC FRACTION)
776 ! NOTE: FROZEN SOIL MOISTURE = SMC - SH2O
777 ! SNOWH ACTUAL SNOW DEPTH (M)
778 ! SNEQV LIQUID WATER-EQUIVALENT SNOW DEPTH (M)
779 ! NOTE: SNOW DENSITY = SNEQV/SNOWH
780 ! ALBEDO SURFACE ALBEDO INCLUDING SNOW EFFECT (UNITLESS FRACTION)
781 ! =SNOW-FREE ALBEDO (ALB) WHEN SNEQV=0, OR
782 ! =FCT(MSNOALB,ALB,VEGTYP,SHDFAC,SHDMIN) WHEN SNEQV>0
783 ! CH SURFACE EXCHANGE COEFFICIENT FOR HEAT AND MOISTURE
784 ! (M S-1); NOTE: CH IS TECHNICALLY A CONDUCTANCE SINCE
785 ! IT HAS BEEN MULTIPLIED BY WIND SPEED.
786 ! CM SURFACE EXCHANGE COEFFICIENT FOR MOMENTUM (M S-1); NOTE:
787 ! CM IS TECHNICALLY A CONDUCTANCE SINCE IT HAS BEEN
788 ! MULTIPLIED BY WIND SPEED. CM IS NOT NEEDED IN SFLX
789 ! ----------------------------------------------------------------------
790 ! 6. OUTPUT (O):
791 ! ----------------------------------------------------------------------
792 ! OUTPUT VARIABLES NECESSARY FOR A COUPLED NUMERICAL WEATHER PREDICTION
793 ! MODEL, E.G. NOAA/NWS/NCEP MESOSCALE ETA MODEL. FOR THIS APPLICATION,
794 ! THE REMAINING OUTPUT/DIAGNOSTIC/PARAMETER BLOCKS BELOW ARE NOT
795 ! NECESSARY. OTHER APPLICATIONS MAY REQUIRE DIFFERENT OUTPUT VARIABLES.
796 ! ETA ACTUAL LATENT HEAT FLUX (W M-2: NEGATIVE, IF UP FROM
797 ! SURFACE)
798 ! SHEAT SENSIBLE HEAT FLUX (W M-2: NEGATIVE, IF UPWARD FROM
799 ! SURFACE)
800 ! ----------------------------------------------------------------------
801 ! EC CANOPY WATER EVAPORATION (W M-2)
802 ! EDIR DIRECT SOIL EVAPORATION (W M-2)
803 ! ET(NSOIL) PLANT TRANSPIRATION FROM A PARTICULAR ROOT (SOIL) LAYER
804 ! (W M-2)
805 ! ETT TOTAL PLANT TRANSPIRATION (W M-2)
806 ! ESNOW SUBLIMATION FROM SNOWPACK (W M-2)
807 ! DRIP THROUGH-FALL OF PRECIP AND/OR DEW IN EXCESS OF CANOPY
808 ! WATER-HOLDING CAPACITY (M)
809 ! DEW DEWFALL (OR FROSTFALL FOR T<273.15) (M)
810 ! ----------------------------------------------------------------------
811 ! BETA RATIO OF ACTUAL/POTENTIAL EVAP (DIMENSIONLESS)
812 ! ETP POTENTIAL EVAPORATION (W M-2)
813 ! SSOIL SOIL HEAT FLUX (W M-2: NEGATIVE IF DOWNWARD FROM SURFACE)
814 ! ----------------------------------------------------------------------
815 ! FLX1 PRECIP-SNOW SFC (W M-2)
816 ! FLX2 FREEZING RAIN LATENT HEAT FLUX (W M-2)
817 ! FLX3 PHASE-CHANGE HEAT FLUX FROM SNOWMELT (W M-2)
818 ! ----------------------------------------------------------------------
819 ! SNOMLT SNOW MELT (M) (WATER EQUIVALENT)
820 ! SNCOVR FRACTIONAL SNOW COVER (UNITLESS FRACTION, 0-1)
821 ! ----------------------------------------------------------------------
822 ! RUNOFF1 SURFACE RUNOFF (M S-1), NOT INFILTRATING THE SURFACE
823 ! RUNOFF2 SUBSURFACE RUNOFF (M S-1), DRAINAGE OUT BOTTOM OF LAST
824 ! SOIL LAYER
825 ! RUNOFF3 NUMERICAL TRUNCTATION IN EXCESS OF POROSITY (SMCMAX)
826 ! FOR A GIVEN SOIL LAYER AT THE END OF A TIME STEP
827 ! ----------------------------------------------------------------------
828 ! RC CANOPY RESISTANCE (S M-1)
829 ! PC PLANT COEFFICIENT (UNITLESS FRACTION, 0-1) WHERE PC*ETP
830 ! = ACTUAL TRANSPIRATION
831 ! XLAI LEAF AREA INDEX (DIMENSIONLESS)
832 ! RSMIN MINIMUM CANOPY RESISTANCE (S M-1)
833 ! RCS INCOMING SOLAR RC FACTOR (DIMENSIONLESS)
834 ! RCT AIR TEMPERATURE RC FACTOR (DIMENSIONLESS)
835 ! RCQ ATMOS VAPOR PRESSURE DEFICIT RC FACTOR (DIMENSIONLESS)
836 ! RCSOIL SOIL MOISTURE RC FACTOR (DIMENSIONLESS)
837 ! ----------------------------------------------------------------------
838 ! 7. DIAGNOSTIC OUTPUT (D):
839 ! ----------------------------------------------------------------------
840 ! SOILW AVAILABLE SOIL MOISTURE IN ROOT ZONE (UNITLESS FRACTION
841 ! BETWEEN SMCWLT AND SMCMAX)
842 ! SOILM TOTAL SOIL COLUMN MOISTURE CONTENT (FROZEN+UNFROZEN) (M)
843 ! ----------------------------------------------------------------------
844 ! 8. PARAMETERS (P):
845 ! ----------------------------------------------------------------------
846 ! SMCWLT WILTING POINT (VOLUMETRIC)
847 ! SMCDRY DRY SOIL MOISTURE THRESHOLD WHERE DIRECT EVAP FRM TOP
848 ! LAYER ENDS (VOLUMETRIC)
849 ! SMCREF SOIL MOISTURE THRESHOLD WHERE TRANSPIRATION BEGINS TO
850 ! STRESS (VOLUMETRIC)
851 ! SMCMAX POROSITY, I.E. SATURATED VALUE OF SOIL MOISTURE
852 ! (VOLUMETRIC)
853 ! NROOT NUMBER OF ROOT LAYERS, A FUNCTION OF VEG TYPE, DETERMINED
854 ! IN SUBROUTINE REDPRM.
855 ! ----------------------------------------------------------------------
856 INTEGER NSOLD
857 PARAMETER(NSOLD = 20)
858
859 ! ----------------------------------------------------------------------
860 ! DECLARATIONS - LOGICAL
861 ! ----------------------------------------------------------------------
862 LOGICAL FRZGRA
863 LOGICAL SATURATED
864 LOGICAL SNOWNG
865
866 ! ----------------------------------------------------------------------
867 ! DECLARATIONS - INTEGER
868 ! ----------------------------------------------------------------------
869 INTEGER ICE
870 INTEGER K
871 INTEGER KZ
872 INTEGER NSOIL
873 INTEGER NROOT
874 INTEGER SLOPETYP
875 INTEGER SOILTYP
876 INTEGER VEGTYP
877 INTEGER I
878 INTEGER J
879
880 ! ----------------------------------------------------------------------
881 ! DECLARATIONS - REAL
882 ! ----------------------------------------------------------------------
883 REAL ALBEDO
884 REAL ALB
885 REAL BEXP
886 REAL BETA
887 REAL CFACTR
888 REAL CH
889 REAL CM
890 REAL CMC
891 REAL CMCMAX
892 REAL CP
893 ! REAL CSNOW
894 REAL CSOIL
895 REAL CZIL
896 REAL DEW
897 REAL DF1
898 REAL DF1H
899 REAL DF1A
900 REAL DKSAT
901 REAL DT
902 REAL DWSAT
903 REAL DQSDT2
904 REAL DSOIL
905 REAL DTOT
906 REAL DRIP
907 REAL EC
908 REAL EDIR
909 REAL ESNOW
910 REAL ET(NSOIL)
911 REAL ETT
912 REAL FRCSNO
913 REAL FRCSOI
914 REAL EPSCA
915 REAL ETA
916 REAL ETP
917 REAL FDOWN
918 REAL F1
919 REAL FLX1
920 REAL FLX2
921 REAL FLX3
922 REAL FXEXP
923 REAL FRZX
924 REAL SHEAT
925 REAL HS
926 REAL KDT
927 REAL LWDN
928 REAL LVH2O
929 REAL PC
930 REAL PRCP
931 REAL PTU
932 REAL PRCP1
933 REAL PSISAT
934 REAL Q2
935 REAL Q2SAT
936 REAL QUARTZ
937 REAL R
938 REAL RCH
939 REAL REFKDT
940 REAL RR
941 REAL RTDIS(NSOLD)
942 REAL RUNOFF1
943 REAL RUNOFF2
944 REAL RGL
945 REAL RUNOFF3
946 REAL RSMAX
947 REAL RC
948 REAL RSMIN
949 REAL RCQ
950 REAL RCS
951 REAL RCSOIL
952 REAL RCT
953 REAL RSNOW
954 REAL SNDENS
955 REAL SNCOND
956 REAL SSOIL
957 REAL SBETA
958 REAL SFCPRS
959 REAL SFCSPD
960 REAL SFCTMP
961 REAL SHDFAC
962 REAL SHDMIN
963 REAL SH2O(NSOIL)
964 REAL SLDPTH(NSOIL)
965 REAL SMCDRY
966 REAL SMCMAX
967 REAL SMCREF
968 REAL SMCWLT
969 REAL SMC(NSOIL)
970 REAL SNEQV
971 REAL SNCOVR
972 REAL SNOWH
973 REAL SN_NEW
974 REAL SLOPE
975 REAL SNUP
976 REAL SALP
977 REAL SNOALB
978 REAL STC(NSOIL)
979 REAL SNOMLT
980 REAL SOLDN
981 REAL SOILM
982 REAL SOILW
983 REAL SOILWM
984 REAL SOILWW
985 REAL T1
986 REAL T1V
987 REAL T24
988 REAL T2V
989 REAL TBOT
990 REAL TH2
991 REAL TH2V
992 REAL TOPT
993 REAL TFREEZ
994 REAL TSNOW
995 REAL XLAI
996 REAL ZLVL
997 REAL ZBOT
998 REAL Z0
999 REAL ZSOIL(NSOLD)
1000
1001 REAL FFROZP
1002 REAL SOLNET
1003 REAL LSUBS
1004
1005 REAL Q1
1006
1007 ! ----------------------------------------------------------------------
1008 ! DECLARATIONS - PARAMETERS
1009 ! ----------------------------------------------------------------------
1010 PARAMETER(TFREEZ = 273.15)
1011 PARAMETER(LVH2O = 2.501E+6)
1012 PARAMETER(LSUBS = 2.83E+6)
1013 PARAMETER(R = 287.04)
1014 PARAMETER(CP = 1004.5)
1015
1016 ! ----------------------------------------------------------------------
1017 ! INITIALIZATION
1018 ! ----------------------------------------------------------------------
1019 RUNOFF1 = 0.0
1020 RUNOFF2 = 0.0
1021 RUNOFF3 = 0.0
1022 SNOMLT = 0.0
1023
1024 ! ----------------------------------------------------------------------
1025 ! THE VARIABLE "ICE" IS A FLAG DENOTING SEA-ICE CASE
1026 ! ----------------------------------------------------------------------
1027 IF (ICE .EQ. 1) THEN
1028
1029 ! ----------------------------------------------------------------------
1030 ! SEA-ICE LAYERS ARE EQUAL THICKNESS AND SUM TO 3 METERS
1031 ! ----------------------------------------------------------------------
1032 DO KZ = 1,NSOIL
1033 ZSOIL(KZ) = -3.*FLOAT(KZ)/FLOAT(NSOIL)
1034 END DO
1035
1036 ELSE
1037
1038 ! ----------------------------------------------------------------------
1039 ! CALCULATE DEPTH (NEGATIVE) BELOW GROUND FROM TOP SKIN SFC TO BOTTOM OF
1040 ! EACH SOIL LAYER. NOTE: SIGN OF ZSOIL IS NEGATIVE (DENOTING BELOW
1041 ! GROUND)
1042 ! ----------------------------------------------------------------------
1043 ZSOIL(1) = -SLDPTH(1)
1044 DO KZ = 2,NSOIL
1045 ZSOIL(KZ) = -SLDPTH(KZ)+ZSOIL(KZ-1)
1046 END DO
1047
1048 ENDIF
1049
1050 ! ----------------------------------------------------------------------
1051 ! NEXT IS CRUCIAL CALL TO SET THE LAND-SURFACE PARAMETERS, INCLUDING
1052 ! SOIL-TYPE AND VEG-TYPE DEPENDENT PARAMETERS.
1053 ! ----------------------------------------------------------------------
1054 CALL REDPRM (VEGTYP,SOILTYP,SLOPETYP, &
1055 & CFACTR,CMCMAX,RSMAX,TOPT,REFKDT,KDT,SBETA, &
1056 & SHDFAC,RSMIN,RGL,HS,ZBOT,FRZX,PSISAT,SLOPE, &
1057 & SNUP,SALP,BEXP,DKSAT,DWSAT,SMCMAX,SMCWLT,SMCREF, &
1058 & SMCDRY,F1,QUARTZ,FXEXP,RTDIS,SLDPTH,ZSOIL, &
1059 & NROOT,NSOIL,Z0,CZIL,XLAI,CSOIL,PTU)
1060
1061 ! ----------------------------------------------------------------------
1062 ! INITIALIZE PRECIPITATION LOGICALS.
1063 ! ----------------------------------------------------------------------
1064 SNOWNG = .FALSE.
1065 FRZGRA = .FALSE.
1066
1067 ! ----------------------------------------------------------------------
1068 ! IF SEA-ICE CASE, ASSIGN DEFAULT WATER-EQUIV SNOW ON TOP
1069 ! ----------------------------------------------------------------------
1070 IF (ICE .EQ. 1) THEN
1071 SNEQV = 0.01
1072 SNOWH = 0.05
1073 SNDENS = SNEQV/SNOWH
1074 ENDIF
1075
1076 ! ----------------------------------------------------------------------
1077 ! IF INPUT SNOWPACK IS NONZERO, THEN COMPUTE SNOW DENSITY "SNDENS" AND
1078 ! SNOW THERMAL CONDUCTIVITY "SNCOND" (NOTE THAT CSNOW IS A FUNCTION
1079 ! SUBROUTINE)
1080 ! ----------------------------------------------------------------------
1081 IF (SNEQV .EQ. 0.0) THEN
1082 SNDENS = 0.0
1083 SNOWH = 0.0
1084 SNCOND = 1.0
1085 ELSE
1086 SNDENS = SNEQV/SNOWH
1087 SNCOND = CSNOW(SNDENS)
1088 ENDIF
1089
1090 ! ----------------------------------------------------------------------
1091 ! DETERMINE IF IT'S PRECIPITATING AND WHAT KIND OF PRECIP IT IS.
1092 ! IF IT'S PRCPING AND THE AIR TEMP IS COLDER THAN 0 C, IT'S SNOWING!
1093 ! IF IT'S PRCPING AND THE AIR TEMP IS WARMER THAN 0 C, BUT THE GRND
1094 ! TEMP IS COLDER THAN 0 C, FREEZING RAIN IS PRESUMED TO BE FALLING.
1095 ! ----------------------------------------------------------------------
1096 IF (PRCP .GT. 0.0) THEN
1097 ! IF (SFCTMP .LE. TFREEZ) THEN
1098 IF (FFROZP .GT. 0.5) THEN
1099 SNOWNG = .TRUE.
1100 ELSE
1101 IF (T1 .LE. TFREEZ) FRZGRA = .TRUE.
1102 ENDIF
1103 ENDIF
1104
1105 ! ----------------------------------------------------------------------
1106 ! IF EITHER PRCP FLAG IS SET, DETERMINE NEW SNOWFALL (CONVERTING PRCP
1107 ! RATE FROM KG M-2 S-1 TO A LIQUID EQUIV SNOW DEPTH IN METERS) AND ADD
1108 ! IT TO THE EXISTING SNOWPACK.
1109 ! NOTE THAT SINCE ALL PRECIP IS ADDED TO SNOWPACK, NO PRECIP INFILTRATES
1110 ! INTO THE SOIL SO THAT PRCP1 IS SET TO ZERO.
1111 ! ----------------------------------------------------------------------
1112 IF ( (SNOWNG) .OR. (FRZGRA) ) THEN
1113 SN_NEW = PRCP * DT * 0.001
1114 SNEQV = SNEQV + SN_NEW
1115 PRCP1 = 0.0
1116
1117 ! ----------------------------------------------------------------------
1118 ! UPDATE SNOW DENSITY BASED ON NEW SNOWFALL, USING OLD AND NEW SNOW.
1119 ! UPDATE SNOW THERMAL CONDUCTIVITY
1120 ! ----------------------------------------------------------------------
1121 CALL SNOW_NEW (SFCTMP,SN_NEW,SNOWH,SNDENS)
1122 SNCOND = CSNOW (SNDENS)
1123 ELSE
1124
1125 ! ----------------------------------------------------------------------
1126 ! PRECIP IS LIQUID (RAIN), HENCE SAVE IN THE PRECIP VARIABLE THAT
1127 ! LATER CAN WHOLELY OR PARTIALLY INFILTRATE THE SOIL (ALONG WITH
1128 ! ANY CANOPY "DRIP" ADDED TO THIS LATER)
1129 ! ----------------------------------------------------------------------
1130 PRCP1 = PRCP
1131
1132 ENDIF
1133
1134 ! ----------------------------------------------------------------------
1135 ! DETERMINE SNOWCOVER AND ALBEDO OVER LAND.
1136 ! ----------------------------------------------------------------------
1137 IF (ICE .EQ. 0) THEN
1138
1139 ! ----------------------------------------------------------------------
1140 ! IF SNOW DEPTH=0, SET SNOW FRACTION=0, ALBEDO=SNOW FREE ALBEDO.
1141 ! ----------------------------------------------------------------------
1142 IF (SNEQV .EQ. 0.0) THEN
1143 SNCOVR = 0.0
1144 ALBEDO = ALB
1145
1146 ELSE
1147 ! ----------------------------------------------------------------------
1148 ! DETERMINE SNOW FRACTIONAL COVERAGE.
1149 ! DETERMINE SURFACE ALBEDO MODIFICATION DUE TO SNOWDEPTH STATE.
1150 ! ----------------------------------------------------------------------
1151 CALL SNFRAC (SNEQV,SNUP,SALP,SNOWH,SNCOVR)
1152 ! MEK JAN 2006, LIMIT SNOW COVER TO A MAXIMUM FRACTION OF 0.98
1153 SNCOVR = MIN(SNCOVR,0.98)
1154 CALL ALCALC (ALB,SNOALB,SHDFAC,SHDMIN,SNCOVR,TSNOW,ALBEDO)
1155 ENDIF
1156
1157 ELSE
1158 ! ----------------------------------------------------------------------
1159 ! SNOW COVER, ALBEDO OVER SEA-ICE
1160 ! ----------------------------------------------------------------------
1161 SNCOVR = 1.0
1162 ! changed in version 2.6 on June 2nd 2003
1163 ! ALBEDO = 0.60
1164 ALBEDO = 0.65
1165 ENDIF
1166
1167 ! ----------------------------------------------------------------------
1168 ! THERMAL CONDUCTIVITY FOR SEA-ICE CASE
1169 ! ----------------------------------------------------------------------
1170 IF (ICE .EQ. 1) THEN
1171 DF1 = 2.2
1172
1173 ELSE
1174
1175 ! ----------------------------------------------------------------------
1176 ! NEXT CALCULATE THE SUBSURFACE HEAT FLUX, WHICH FIRST REQUIRES
1177 ! CALCULATION OF THE THERMAL DIFFUSIVITY. TREATMENT OF THE
1178 ! LATTER FOLLOWS THAT ON PAGES 148-149 FROM "HEAT TRANSFER IN
1179 ! COLD CLIMATES", BY V. J. LUNARDINI (PUBLISHED IN 1981
1180 ! BY VAN NOSTRAND REINHOLD CO.) I.E. TREATMENT OF TWO CONTIGUOUS
1181 ! "PLANE PARALLEL" MEDIUMS (NAMELY HERE THE FIRST SOIL LAYER
1182 ! AND THE SNOWPACK LAYER, IF ANY). THIS DIFFUSIVITY TREATMENT
1183 ! BEHAVES WELL FOR BOTH ZERO AND NONZERO SNOWPACK, INCLUDING THE
1184 ! LIMIT OF VERY THIN SNOWPACK. THIS TREATMENT ALSO ELIMINATES
1185 ! THE NEED TO IMPOSE AN ARBITRARY UPPER BOUND ON SUBSURFACE
1186 ! HEAT FLUX WHEN THE SNOWPACK BECOMES EXTREMELY THIN.
1187 ! ----------------------------------------------------------------------
1188 ! FIRST CALCULATE THERMAL DIFFUSIVITY OF TOP SOIL LAYER, USING
1189 ! BOTH THE FROZEN AND LIQUID SOIL MOISTURE, FOLLOWING THE
1190 ! SOIL THERMAL DIFFUSIVITY FUNCTION OF PETERS-LIDARD ET AL.
1191 ! (1998,JAS, VOL 55, 1209-1224), WHICH REQUIRES THE SPECIFYING
1192 ! THE QUARTZ CONTENT OF THE GIVEN SOIL CLASS (SEE ROUTINE REDPRM)
1193 ! ----------------------------------------------------------------------
1194 CALL TDFCND (DF1,SMC(1),QUARTZ,SMCMAX,SH2O(1))
1195
1196 ! ----------------------------------------------------------------------
1197 ! NEXT ADD SUBSURFACE HEAT FLUX REDUCTION EFFECT FROM THE
1198 ! OVERLYING GREEN CANOPY, ADAPTED FROM SECTION 2.1.2 OF
1199 ! PETERS-LIDARD ET AL. (1997, JGR, VOL 102(D4))
1200 ! ----------------------------------------------------------------------
1201 DF1 = DF1 * EXP(SBETA*SHDFAC)
1202 ENDIF
1203
1204 ! ----------------------------------------------------------------------
1205 ! FINALLY "PLANE PARALLEL" SNOWPACK EFFECT FOLLOWING
1206 ! V.J. LINARDINI REFERENCE CITED ABOVE. NOTE THAT DTOT IS
1207 ! COMBINED DEPTH OF SNOWDEPTH AND THICKNESS OF FIRST SOIL LAYER
1208 ! ----------------------------------------------------------------------
1209 DSOIL = -(0.5 * ZSOIL(1))
1210
1211 IF (SNEQV .EQ. 0.) THEN
1212 SSOIL = DF1 * (T1 - STC(1) ) / DSOIL
1213 ELSE
1214 DTOT = SNOWH + DSOIL
1215 FRCSNO = SNOWH/DTOT
1216 FRCSOI = DSOIL/DTOT
1217 !
1218 ! 1. HARMONIC MEAN (SERIES FLOW)
1219 ! DF1 = (SNCOND*DF1)/(FRCSOI*SNCOND+FRCSNO*DF1)
1220 DF1H = (SNCOND*DF1)/(FRCSOI*SNCOND+FRCSNO*DF1)
1221 ! 2. ARITHMETIC MEAN (PARALLEL FLOW)
1222 ! DF1 = FRCSNO*SNCOND + FRCSOI*DF1
1223 DF1A = FRCSNO*SNCOND + FRCSOI*DF1
1224 !
1225 ! 3. GEOMETRIC MEAN (INTERMEDIATE BETWEEN HARMONIC AND ARITHMETIC MEAN)
1226 ! DF1 = (SNCOND**FRCSNO)*(DF1**FRCSOI)
1227 ! TEST - MBEK, 10 Jan 2002
1228 ! weigh DF by snow fraction
1229 ! DF1 = DF1H*SNCOVR + DF1A*(1.0-SNCOVR)
1230 ! DF1 = DF1H*SNCOVR + DF1*(1.0-SNCOVR)
1231 DF1 = DF1A*SNCOVR + DF1*(1.0-SNCOVR)
1232
1233 ! ----------------------------------------------------------------------
1234 ! CALCULATE SUBSURFACE HEAT FLUX, SSOIL, FROM FINAL THERMAL DIFFUSIVITY
1235 ! OF SURFACE MEDIUMS, DF1 ABOVE, AND SKIN TEMPERATURE AND TOP
1236 ! MID-LAYER SOIL TEMPERATURE
1237 ! ----------------------------------------------------------------------
1238 SSOIL = DF1 * (T1 - STC(1) ) / DTOT
1239 ENDIF
1240
1241 ! MEK -- DEBUG -- AUG 2005
1242 ! WRITE(*,*) 'T1,STC(1),DSOIL=',T1,STC(1),DSOIL
1243 ! WRITE(*,*) 'DF1,SBETA,SHDFAC=',DF1,SBETA,SHDFAC
1244 ! WRITE(*,*) 'SSOIL=',SSOIL
1245
1246 ! ----------------------------------------------------------------------
1247 ! DETERMINE SURFACE ROUGHNESS OVER SNOWPACK USING SNOW CONDITION FROM
1248 ! THE PREVIOUS TIMESTEP.
1249 ! ----------------------------------------------------------------------
1250 IF (SNCOVR .GT. 0.) THEN
1251 CALL SNOWZ0 (SNCOVR,Z0)
1252 ENDIF
1253
1254 ! ----------------------------------------------------------------------
1255 ! NEXT CALL ROUTINE SFCDIF TO CALCULATE THE SFC EXCHANGE COEF (CH) FOR
1256 ! HEAT AND MOISTURE.
1257 !
1258 ! NOTE !!!
1259 ! COMMENT OUT CALL SFCDIF, IF SFCDIF ALREADY CALLED IN CALLING PROGRAM
1260 ! (SUCH AS IN COUPLED ATMOSPHERIC MODEL).
1261 !
1262 ! NOTE !!!
1263 ! DO NOT CALL SFCDIF UNTIL AFTER ABOVE CALL TO REDPRM, IN CASE
1264 ! ALTERNATIVE VALUES OF ROUGHNESS LENGTH (Z0) AND ZILINTINKEVICH COEF
1265 ! (CZIL) ARE SET THERE VIA NAMELIST I/O.
1266 !
1267 ! NOTE !!!
1268 ! ROUTINE SFCDIF RETURNS A CH THAT REPRESENTS THE WIND SPD TIMES THE
1269 ! "ORIGINAL" NONDIMENSIONAL "Ch" TYPICAL IN LITERATURE. HENCE THE CH
1270 ! RETURNED FROM SFCDIF HAS UNITS OF M/S. THE IMPORTANT COMPANION
1271 ! COEFFICIENT OF CH, CARRIED HERE AS "RCH", IS THE CH FROM SFCDIF TIMES
1272 ! AIR DENSITY AND PARAMETER "CP". "RCH" IS COMPUTED IN "CALL PENMAN".
1273 ! RCH RATHER THAN CH IS THE COEFF USUALLY INVOKED LATER IN EQNS.
1274 !
1275 ! NOTE !!!
1276 ! SFCDIF ALSO RETURNS THE SURFACE EXCHANGE COEFFICIENT FOR MOMENTUM, CM,
1277 ! ALSO KNOWN AS THE SURFACE DRAGE COEFFICIENT, BUT CM IS NOT USED HERE.
1278 ! ----------------------------------------------------------------------
1279 ! CALC VIRTUAL TEMPS AND VIRTUAL POTENTIAL TEMPS NEEDED BY SUBROUTINES
1280 ! SFCDIF AND PENMAN.
1281 ! ----------------------------------------------------------------------
1282 T2V = SFCTMP * (1.0 + 0.61 * Q2 )
1283 ! ----------------------------------------------------------------------
1284 ! COMMENT OUT BELOW 2 LINES IF CALL SFCDIF IS COMMENTED OUT, I.E. IN THE
1285 ! COUPLED MODEL.
1286 ! ----------------------------------------------------------------------
1287 ! T1V = T1 * (1.0 + 0.61 * Q2)
1288 ! TH2V = TH2 * (1.0 + 0.61 * Q2)
1289 !
1290 ! CALL SFCDIF (ZLVL,Z0,T1V,TH2V,SFCSPD,CZIL,CM,CH)
1291
1292 ! ----------------------------------------------------------------------
1293 ! CALCULATE TOTAL DOWNWARD RADIATION (SOLAR PLUS LONGWAVE) NEEDED IN
1294 ! PENMAN EP SUBROUTINE THAT FOLLOWS
1295 ! ----------------------------------------------------------------------
1296 ! FDOWN = SOLDN*(1.0-ALBEDO) + LWDN
1297 FDOWN = SOLNET + LWDN
1298
1299 ! ----------------------------------------------------------------------
1300 ! CALL PENMAN SUBROUTINE TO CALCULATE POTENTIAL EVAPORATION (ETP), AND
1301 ! OTHER PARTIAL PRODUCTS AND SUMS SAVE IN COMMON/RITE FOR LATER
1302 ! CALCULATIONS.
1303 ! ----------------------------------------------------------------------
1304 CALL PENMAN (SFCTMP,SFCPRS,CH,T2V,TH2,PRCP,FDOWN,T24,SSOIL, &
1305 & Q2,Q2SAT,ETP,RCH,EPSCA,RR,SNOWNG,FRZGRA, &
1306 & DQSDT2,FLX2)
1307
1308 ! ----------------------------------------------------------------------
1309 ! CALL CANRES TO CALCULATE THE CANOPY RESISTANCE AND CONVERT IT INTO PC
1310 ! IF NONZERO GREENNESS FRACTION
1311 ! ----------------------------------------------------------------------
1312 IF (SHDFAC .GT. 0.) THEN
1313
1314 ! ----------------------------------------------------------------------
1315 ! FROZEN GROUND EXTENSION: TOTAL SOIL WATER "SMC" WAS REPLACED
1316 ! BY UNFROZEN SOIL WATER "SH2O" IN CALL TO CANRES BELOW
1317 ! ----------------------------------------------------------------------
1318 CALL CANRES (SOLDN,CH,SFCTMP,Q2,SFCPRS,SH2O,ZSOIL,NSOIL, &
1319 & SMCWLT,SMCREF,RSMIN,RC,PC,NROOT,Q2SAT,DQSDT2, &
1320 & TOPT,RSMAX,RGL,HS,XLAI, &
1321 & RCS,RCT,RCQ,RCSOIL)
1322
1323 ENDIF
1324
1325 ! ----------------------------------------------------------------------
1326 ! NOW DECIDE MAJOR PATHWAY BRANCH TO TAKE DEPENDING ON WHETHER SNOWPACK
1327 ! EXISTS OR NOT:
1328 ! ----------------------------------------------------------------------
1329 ESNOW = 0.0
1330 IF (SNEQV .EQ. 0.0) THEN
1331 CALL NOPAC (ETP,ETA,PRCP,SMC,SMCMAX,SMCWLT, &
1332 & SMCREF,SMCDRY,CMC,CMCMAX,NSOIL,DT,SHDFAC, &
1333 & SBETA,Q2,T1,SFCTMP,T24,TH2,FDOWN,F1,SSOIL, &
1334 & STC,EPSCA,BEXP,PC,RCH,RR,CFACTR, &
1335 & SH2O,SLOPE,KDT,FRZX,PSISAT,ZSOIL, &
1336 & DKSAT,DWSAT,TBOT,ZBOT,RUNOFF1,RUNOFF2, &
1337 & RUNOFF3,EDIR,EC,ET,ETT,NROOT,ICE,RTDIS, &
1338 & QUARTZ,FXEXP,CSOIL, &
1339 & BETA,DRIP,DEW,FLX1,FLX2,FLX3)
1340 ELSE
1341 CALL SNOPAC (ETP,ETA,PRCP,PRCP1,SNOWNG,SMC,SMCMAX,SMCWLT, &
1342 & SMCREF,SMCDRY,CMC,CMCMAX,NSOIL,DT, &
1343 & SBETA,DF1, &
1344 & Q2,T1,SFCTMP,T24,TH2,FDOWN,F1,SSOIL,STC,EPSCA, &
1345 & SFCPRS,BEXP,PC,RCH,RR,CFACTR,SNCOVR,SNEQV,SNDENS, &
1346 & SNOWH,SH2O,SLOPE,KDT,FRZX,PSISAT,SNUP, &
1347 & ZSOIL,DWSAT,DKSAT,TBOT,ZBOT,SHDFAC,RUNOFF1, &
1348 & RUNOFF2,RUNOFF3,EDIR,EC,ET,ETT,NROOT,SNOMLT, &
1349 & ICE,RTDIS,QUARTZ,FXEXP,CSOIL, &
1350 & BETA,DRIP,DEW,FLX1,FLX2,FLX3,ESNOW)
1351 ! ESNOW = ETA
1352 ENDIF
1353
1354 ! ----------------------------------------------------------------------
1355 ! PREPARE SENSIBLE HEAT (H) FOR RETURN TO PARENT MODEL
1356 ! ----------------------------------------------------------------------
1357 SHEAT = -(CH * CP * SFCPRS)/(R * T2V) * ( TH2 - T1 )
1358
1359 ! ----------------------------------------------------------------------
1360 ! CONVERT UNITS AND/OR SIGN OF TOTAL EVAP (ETA), POTENTIAL EVAP (ETP),
1361 ! SUBSURFACE HEAT FLUX (S), AND RUNOFFS FOR WHAT PARENT MODEL EXPECTS
1362 ! CONVERT ETA FROM KG M-2 S-1 TO W M-2
1363 ! ----------------------------------------------------------------------
1364 ! ETA = ETA*LVH2O
1365 ! ETP = ETP*LVH2O
1366
1367 ! ----------------------------------------------------------------------
1368 EDIR = EDIR * LVH2O
1369 EC = EC * LVH2O
1370 DO K=1,4
1371 ET(K) = ET(K) * LVH2O
1372 ENDDO
1373 ETT = ETT * LVH2O
1374 ESNOW = ESNOW * LSUBS
1375 ETP = ETP*((1.-SNCOVR)*LVH2O + SNCOVR*LSUBS)
1376 IF (ETP .GT. 0.) THEN
1377 ETA = EDIR + EC + ETT + ESNOW
1378 ELSE
1379 ETA = ETP
1380 ENDIF
1381 ! ----------------------------------------------------------------------
1382 ! DETERMINE BETA (RATIO OF ACTUAL TO POTENTIAL EVAP)
1383 ! ----------------------------------------------------------------------
1384 IF (ETP == 0.0) THEN
1385 BETA = 0.0
1386 ELSE
1387 BETA = ETA/ETP
1388 ENDIF
1389
1390 ! ----------------------------------------------------------------------
1391
1392 ! ----------------------------------------------------------------------
1393 ! CONVERT THE SIGN OF SOIL HEAT FLUX SO THAT:
1394 ! SSOIL>0: WARM THE SURFACE (NIGHT TIME)
1395 ! SSOIL<0: COOL THE SURFACE (DAY TIME)
1396 ! ----------------------------------------------------------------------
1397 SSOIL = -1.0*SSOIL
1398
1399 ! ----------------------------------------------------------------------
1400 ! CONVERT RUNOFF3 (INTERNAL LAYER RUNOFF FROM SUPERSAT) FROM M TO M S-1
1401 ! AND ADD TO SUBSURFACE RUNOFF/DRAINAGE/BASEFLOW
1402 ! ----------------------------------------------------------------------
1403 RUNOFF3 = RUNOFF3/DT
1404 RUNOFF2 = RUNOFF2+RUNOFF3
1405
1406 ! ----------------------------------------------------------------------
1407 ! TOTAL COLUMN SOIL MOISTURE IN METERS (SOILM) AND ROOT-ZONE
1408 ! SOIL MOISTURE AVAILABILITY (FRACTION) RELATIVE TO POROSITY/SATURATION
1409 ! ----------------------------------------------------------------------
1410 SOILM = -1.0*SMC(1)*ZSOIL(1)
1411 DO K = 2,NSOIL
1412 SOILM = SOILM+SMC(K)*(ZSOIL(K-1)-ZSOIL(K))
1413 END DO
1414 SOILWM = -1.0*(SMCMAX-SMCWLT)*ZSOIL(1)
1415 SOILWW = -1.0*(SMC(1)-SMCWLT)*ZSOIL(1)
1416 DO K = 2,NROOT
1417 SOILWM = SOILWM+(SMCMAX-SMCWLT)*(ZSOIL(K-1)-ZSOIL(K))
1418 SOILWW = SOILWW+(SMC(K)-SMCWLT)*(ZSOIL(K-1)-ZSOIL(K))
1419 END DO
1420 SOILW = SOILWW/SOILWM
1421
1422 ! ----------------------------------------------------------------------
1423 ! END SUBROUTINE SFLX
1424 ! ----------------------------------------------------------------------
1425 END SUBROUTINE SFLX
1426
1427 SUBROUTINE ALCALC (ALB,SNOALB,SHDFAC,SHDMIN,SNCOVR,TSNOW,ALBEDO)
1428
1429 IMPLICIT NONE
1430
1431 ! ----------------------------------------------------------------------
1432 ! CALCULATE ALBEDO INCLUDING SNOW EFFECT (0 -> 1)
1433 ! ALB SNOWFREE ALBEDO
1434 ! SNOALB MAXIMUM (DEEP) SNOW ALBEDO
1435 ! SHDFAC AREAL FRACTIONAL COVERAGE OF GREEN VEGETATION
1436 ! SHDMIN MINIMUM AREAL FRACTIONAL COVERAGE OF GREEN VEGETATION
1437 ! SNCOVR FRACTIONAL SNOW COVER
1438 ! ALBEDO SURFACE ALBEDO INCLUDING SNOW EFFECT
1439 ! TSNOW SNOW SURFACE TEMPERATURE (K)
1440 ! ----------------------------------------------------------------------
1441 REAL ALB, SNOALB, SHDFAC, SHDMIN, SNCOVR, ALBEDO, TSNOW
1442
1443 ! ----------------------------------------------------------------------
1444 ! SNOALB IS ARGUMENT REPRESENTING MAXIMUM ALBEDO OVER DEEP SNOW,
1445 ! AS PASSED INTO SFLX, AND ADAPTED FROM THE SATELLITE-BASED MAXIMUM
1446 ! SNOW ALBEDO FIELDS PROVIDED BY D. ROBINSON AND G. KUKLA
1447 ! (1985, JCAM, VOL 24, 402-411)
1448 ! ----------------------------------------------------------------------
1449 ! changed in version 2.6 on June 2nd 2003
1450 ! ALBEDO = ALB + (1.0-(SHDFAC-SHDMIN))*SNCOVR*(SNOALB-ALB)
1451 ALBEDO = ALB + SNCOVR*(SNOALB-ALB)
1452 IF (ALBEDO .GT. SNOALB) ALBEDO=SNOALB
1453
1454 ! BASE FORMULATION (DICKINSON ET AL., 1986, COGLEY ET AL., 1990)
1455 ! IF (TSNOW.LE.263.16) THEN
1456 ! ALBEDO=SNOALB
1457 ! ELSE
1458 ! IF (TSNOW.LT.273.16) THEN
1459 ! TM=0.1*(TSNOW-263.16)
1460 ! ALBEDO=0.5*((0.9-0.2*(TM**3))+(0.8-0.16*(TM**3)))
1461 ! ELSE
1462 ! ALBEDO=0.67
1463 ! ENDIF
1464 ! ENDIF
1465
1466 ! ISBA FORMULATION (VERSEGHY, 1991; BAKER ET AL., 1990)
1467 ! IF (TSNOW.LT.273.16) THEN
1468 ! ALBEDO=SNOALB-0.008*DT/86400
1469 ! ELSE
1470 ! ALBEDO=(SNOALB-0.5)*EXP(-0.24*DT/86400)+0.5
1471 ! ENDIF
1472
1473 ! ----------------------------------------------------------------------
1474 ! END SUBROUTINE ALCALC
1475 ! ----------------------------------------------------------------------
1476 END SUBROUTINE ALCALC
1477
1478 SUBROUTINE CANRES (SOLAR,CH,SFCTMP,Q2,SFCPRS,SMC,ZSOIL,NSOIL, &
1479 & SMCWLT,SMCREF,RSMIN,RC,PC,NROOT,Q2SAT,DQSDT2, &
1480 & TOPT,RSMAX,RGL,HS,XLAI, &
1481 & RCS,RCT,RCQ,RCSOIL)
1482
1483 IMPLICIT NONE
1484
1485 ! ----------------------------------------------------------------------
1486 ! SUBROUTINE CANRES
1487 ! ----------------------------------------------------------------------
1488 ! CALCULATE CANOPY RESISTANCE WHICH DEPENDS ON INCOMING SOLAR RADIATION,
1489 ! AIR TEMPERATURE, ATMOSPHERIC WATER VAPOR PRESSURE DEFICIT AT THE
1490 ! LOWEST MODEL LEVEL, AND SOIL MOISTURE (PREFERABLY UNFROZEN SOIL
1491 ! MOISTURE RATHER THAN TOTAL)
1492 ! ----------------------------------------------------------------------
1493 ! SOURCE: JARVIS (1976), NOILHAN AND PLANTON (1989, MWR), JACQUEMIN AND
1494 ! NOILHAN (1990, BLM)
1495 ! SEE ALSO: CHEN ET AL (1996, JGR, VOL 101(D3), 7251-7268), EQNS 12-14
1496 ! AND TABLE 2 OF SEC. 3.1.2
1497 ! ----------------------------------------------------------------------
1498 ! INPUT:
1499 ! SOLAR INCOMING SOLAR RADIATION
1500 ! CH SURFACE EXCHANGE COEFFICIENT FOR HEAT AND MOISTURE
1501 ! SFCTMP AIR TEMPERATURE AT 1ST LEVEL ABOVE GROUND
1502 ! Q2 AIR HUMIDITY AT 1ST LEVEL ABOVE GROUND
1503 ! Q2SAT SATURATION AIR HUMIDITY AT 1ST LEVEL ABOVE GROUND
1504 ! DQSDT2 SLOPE OF SATURATION HUMIDITY FUNCTION WRT TEMP
1505 ! SFCPRS SURFACE PRESSURE
1506 ! SMC VOLUMETRIC SOIL MOISTURE
1507 ! ZSOIL SOIL DEPTH (NEGATIVE SIGN, AS IT IS BELOW GROUND)
1508 ! NSOIL NO. OF SOIL LAYERS
1509 ! NROOT NO. OF SOIL LAYERS IN ROOT ZONE (1.LE.NROOT.LE.NSOIL)
1510 ! XLAI LEAF AREA INDEX
1511 ! SMCWLT WILTING POINT
1512 ! SMCREF REFERENCE SOIL MOISTURE (WHERE SOIL WATER DEFICIT STRESS
1513 ! SETS IN)
1514 ! RSMIN, RSMAX, TOPT, RGL, HS ARE CANOPY STRESS PARAMETERS SET IN
1515 ! SURBOUTINE REDPRM
1516 ! OUTPUT:
1517 ! PC PLANT COEFFICIENT
1518 ! RC CANOPY RESISTANCE
1519 ! ----------------------------------------------------------------------
1520 INTEGER NSOLD
1521 PARAMETER(NSOLD = 20)
1522
1523 INTEGER K
1524 INTEGER NROOT
1525 INTEGER NSOIL
1526
1527 REAL CH
1528 REAL CP
1529 REAL DELTA
1530 REAL DQSDT2
1531 REAL FF
1532 REAL GX
1533 REAL HS
1534 REAL P
1535 REAL PART(NSOLD)
1536 REAL PC
1537 REAL Q2
1538 REAL Q2SAT
1539 REAL RC
1540 REAL RSMIN
1541 REAL RCQ
1542 REAL RCS
1543 REAL RCSOIL
1544 REAL RCT
1545 REAL RD
1546 REAL RGL
1547 REAL RR
1548 REAL RSMAX
1549 REAL SFCPRS
1550 REAL SFCTMP
1551 REAL SIGMA
1552 REAL SLV
1553 REAL SMC(NSOIL)
1554 REAL SMCREF
1555 REAL SMCWLT
1556 REAL SOLAR
1557 REAL TOPT
1558 REAL SLVCP
1559 REAL ST1
1560 REAL TAIR4
1561 REAL XLAI
1562 REAL ZSOIL(NSOIL)
1563
1564 PARAMETER(CP = 1004.5)
1565 PARAMETER(RD = 287.04)
1566 PARAMETER(SIGMA = 5.67E-8)
1567 PARAMETER(SLV = 2.501000E6)
1568
1569 ! ----------------------------------------------------------------------
1570 ! INITIALIZE CANOPY RESISTANCE MULTIPLIER TERMS.
1571 ! ----------------------------------------------------------------------
1572 RCS = 0.0
1573 RCT = 0.0
1574 RCQ = 0.0
1575 RCSOIL = 0.0
1576 RC = 0.0
1577
1578 ! ----------------------------------------------------------------------
1579 ! CONTRIBUTION DUE TO INCOMING SOLAR RADIATION
1580 ! ----------------------------------------------------------------------
1581 FF = 0.55*2.0*SOLAR/(RGL*XLAI)
1582 RCS = (FF + RSMIN/RSMAX) / (1.0 + FF)
1583 RCS = MAX(RCS,0.0001)
1584
1585 ! ----------------------------------------------------------------------
1586 ! CONTRIBUTION DUE TO AIR TEMPERATURE AT FIRST MODEL LEVEL ABOVE GROUND
1587 ! RCT EXPRESSION FROM NOILHAN AND PLANTON (1989, MWR).
1588 ! ----------------------------------------------------------------------
1589 RCT = 1.0 - 0.0016*((TOPT-SFCTMP)**2.0)
1590 RCT = MAX(RCT,0.0001)
1591
1592 ! ----------------------------------------------------------------------
1593 ! CONTRIBUTION DUE TO VAPOR PRESSURE DEFICIT AT FIRST MODEL LEVEL.
1594 ! RCQ EXPRESSION FROM SSIB
1595 ! ----------------------------------------------------------------------
1596 RCQ = 1.0/(1.0+HS*(Q2SAT-Q2))
1597 RCQ = MAX(RCQ,0.01)
1598
1599 ! ----------------------------------------------------------------------
1600 ! CONTRIBUTION DUE TO SOIL MOISTURE AVAILABILITY.
1601 ! DETERMINE CONTRIBUTION FROM EACH SOIL LAYER, THEN ADD THEM UP.
1602 ! ----------------------------------------------------------------------
1603 GX = (SMC(1) - SMCWLT) / (SMCREF - SMCWLT)
1604 IF (GX .GT. 1.) GX = 1.
1605 IF (GX .LT. 0.) GX = 0.
1606
1607 ! ----------------------------------------------------------------------
1608 ! USE SOIL DEPTH AS WEIGHTING FACTOR
1609 ! ----------------------------------------------------------------------
1610 PART(1) = (ZSOIL(1)/ZSOIL(NROOT)) * GX
1611 ! ----------------------------------------------------------------------
1612 ! USE ROOT DISTRIBUTION AS WEIGHTING FACTOR
1613 ! PART(1) = RTDIS(1) * GX
1614 ! ----------------------------------------------------------------------
1615 IF (NROOT .GT. 1) THEN
1616 DO K = 2,NROOT
1617 GX = (SMC(K) - SMCWLT) / (SMCREF - SMCWLT)
1618 IF (GX .GT. 1.) GX = 1.
1619 IF (GX .LT. 0.) GX = 0.
1620 ! ----------------------------------------------------------------------
1621 ! USE SOIL DEPTH AS WEIGHTING FACTOR
1622 ! ----------------------------------------------------------------------
1623 PART(K) = ((ZSOIL(K)-ZSOIL(K-1))/ZSOIL(NROOT)) * GX
1624 ! ----------------------------------------------------------------------
1625 ! USE ROOT DISTRIBUTION AS WEIGHTING FACTOR
1626 ! PART(K) = RTDIS(K) * GX
1627 ! ----------------------------------------------------------------------
1628 END DO
1629 ENDIF
1630
1631 DO K = 1,NROOT
1632 RCSOIL = RCSOIL+PART(K)
1633 END DO
1634 RCSOIL = MAX(RCSOIL,0.0001)
1635
1636 ! ----------------------------------------------------------------------
1637 ! DETERMINE CANOPY RESISTANCE DUE TO ALL FACTORS. CONVERT CANOPY
1638 ! RESISTANCE (RC) TO PLANT COEFFICIENT (PC) TO BE USED WITH POTENTIAL
1639 ! EVAP IN DETERMINING ACTUAL EVAP. PC IS DETERMINED BY:
1640 ! PC * LINERIZED PENMAN POTENTIAL EVAP =
1641 ! PENMAN-MONTEITH ACTUAL EVAPORATION (CONTAINING RC TERM).
1642 ! ----------------------------------------------------------------------
1643 RC = RSMIN/(XLAI*RCS*RCT*RCQ*RCSOIL)
1644
1645 ! TAIR4 = SFCTMP**4.
1646 ! ST1 = (4.*SIGMA*RD)/CP
1647 ! SLVCP = SLV/CP
1648 ! RR = ST1*TAIR4/(SFCPRS*CH) + 1.0
1649 RR = (4.*SIGMA*RD/CP)*(SFCTMP**4.)/(SFCPRS*CH) + 1.0
1650 DELTA = (SLV/CP)*DQSDT2
1651
1652 PC = (RR+DELTA)/(RR*(1.+RC*CH)+DELTA)
1653
1654 ! ----------------------------------------------------------------------
1655 ! END SUBROUTINE CANRES
1656 ! ----------------------------------------------------------------------
1657 END SUBROUTINE CANRES
1658
1659 SUBROUTINE DEVAP (EDIR1,ETP1,SMC,ZSOIL,SHDFAC,SMCMAX,BEXP, &
1660 & DKSAT,DWSAT,SMCDRY,SMCREF,SMCWLT,FXEXP)
1661
1662 IMPLICIT NONE
1663
1664 ! ----------------------------------------------------------------------
1665 ! SUBROUTINE DEVAP
1666 ! ----------------------------------------------------------------------
1667 ! CALCULATE DIRECT SOIL EVAPORATION
1668 ! ----------------------------------------------------------------------
1669 REAL BEXP
1670 ! REAL DEVAP
1671 REAL EDIR1
1672 REAL DKSAT
1673 REAL DWSAT
1674 REAL ETP1
1675 REAL FX
1676 REAL FXEXP
1677 REAL SHDFAC
1678 REAL SMC
1679 REAL SMCDRY
1680 REAL SMCMAX
1681 REAL ZSOIL
1682 REAL SMCREF
1683 REAL SMCWLT
1684 REAL SRATIO
1685
1686 ! ----------------------------------------------------------------------
1687 ! DIRECT EVAP A FUNCTION OF RELATIVE SOIL MOISTURE AVAILABILITY, LINEAR
1688 ! WHEN FXEXP=1.
1689 ! FX > 1 REPRESENTS DEMAND CONTROL
1690 ! FX < 1 REPRESENTS FLUX CONTROL
1691 ! ----------------------------------------------------------------------
1692 SRATIO = (SMC - SMCDRY) / (SMCMAX - SMCDRY)
1693 IF (SRATIO .GT. 0.) THEN
1694 FX = SRATIO**FXEXP
1695 FX = MAX ( MIN ( FX, 1. ) ,0. )
1696 ELSE
1697 FX = 0.
1698 ENDIF
1699
1700 ! ----------------------------------------------------------------------
1701 ! ALLOW FOR THE DIRECT-EVAP-REDUCING EFFECT OF SHADE
1702 ! ----------------------------------------------------------------------
1703 ! DEVAP = FX * ( 1.0 - SHDFAC ) * ETP1
1704 EDIR1 = FX * ( 1.0 - SHDFAC ) * ETP1
1705
1706 ! ----------------------------------------------------------------------
1707 ! END SUBROUTINE DEVAP
1708 ! ----------------------------------------------------------------------
1709 END SUBROUTINE DEVAP
1710
1711 SUBROUTINE EVAPO (ETA1,SMC,NSOIL,CMC,ETP1,DT,ZSOIL, &
1712 & SH2O, &
1713 & SMCMAX,BEXP,PC,SMCWLT,DKSAT,DWSAT, &
1714 & SMCREF,SHDFAC,CMCMAX, &
1715 & SMCDRY,CFACTR, &
1716 & EDIR1,EC1,ET1,ETT1,SFCTMP,Q2,NROOT,RTDIS,FXEXP)
1717
1718 IMPLICIT NONE
1719
1720 ! ----------------------------------------------------------------------
1721 ! SUBROUTINE EVAPO
1722 ! ----------------------------------------------------------------------
1723 ! CALCULATE SOIL MOISTURE FLUX. THE SOIL MOISTURE CONTENT (SMC - A PER
1724 ! UNIT VOLUME MEASUREMENT) IS A DEPENDENT VARIABLE THAT IS UPDATED WITH
1725 ! PROGNOSTIC EQNS. THE CANOPY MOISTURE CONTENT (CMC) IS ALSO UPDATED.
1726 ! FROZEN GROUND VERSION: NEW STATES ADDED: SH2O, AND FROZEN GROUND
1727 ! CORRECTION FACTOR, FRZFACT AND PARAMETER SLOPE.
1728 ! ----------------------------------------------------------------------
1729 INTEGER NSOLD
1730 PARAMETER(NSOLD = 20)
1731
1732 INTEGER I
1733 INTEGER K
1734 INTEGER NSOIL
1735 INTEGER NROOT
1736
1737 REAL BEXP
1738 REAL CFACTR
1739 REAL CMC
1740 REAL CMC2MS
1741 REAL CMCMAX
1742 ! REAL DEVAP
1743 REAL DKSAT
1744 REAL DT
1745 REAL DWSAT
1746 REAL EC1
1747 REAL EDIR1
1748 REAL ET1(NSOIL)
1749 REAL ETA1
1750 REAL ETP1
1751 REAL ETT1
1752 REAL FXEXP
1753 REAL PC
1754 REAL Q2
1755 REAL RTDIS(NSOIL)
1756 REAL SFCTMP
1757 REAL SHDFAC
1758 REAL SMC(NSOIL)
1759 REAL SH2O(NSOIL)
1760 REAL SMCDRY
1761 REAL SMCMAX
1762 REAL SMCREF
1763 REAL SMCWLT
1764 REAL ZSOIL(NSOIL)
1765
1766 ! ----------------------------------------------------------------------
1767 ! EXECUTABLE CODE BEGINS HERE IF THE POTENTIAL EVAPOTRANSPIRATION IS
1768 ! GREATER THAN ZERO.
1769 ! ----------------------------------------------------------------------
1770 EDIR1 = 0.
1771 EC1 = 0.
1772 DO K = 1,NSOIL
1773 ET1(K) = 0.
1774 END DO
1775 ETT1 = 0.
1776
1777 IF (ETP1 .GT. 0.0) THEN
1778
1779 ! ----------------------------------------------------------------------
1780 ! RETRIEVE DIRECT EVAPORATION FROM SOIL SURFACE. CALL THIS FUNCTION
1781 ! ONLY IF VEG COVER NOT COMPLETE.
1782 ! FROZEN GROUND VERSION: SH2O STATES REPLACE SMC STATES.
1783 ! ----------------------------------------------------------------------
1784 IF (SHDFAC .LT. 1.) THEN
1785 CALL DEVAP (EDIR1,ETP1,SH2O(1),ZSOIL(1),SHDFAC,SMCMAX, &
1786 ! EDIR = DEVAP(ETP1,SH2O(1),ZSOIL(1),SHDFAC,SMCMAX, &
1787 & BEXP,DKSAT,DWSAT,SMCDRY,SMCREF,SMCWLT,FXEXP)
1788 ENDIF
1789
1790 ! ----------------------------------------------------------------------
1791 ! INITIALIZE PLANT TOTAL TRANSPIRATION, RETRIEVE PLANT TRANSPIRATION,
1792 ! AND ACCUMULATE IT FOR ALL SOIL LAYERS.
1793 ! ----------------------------------------------------------------------
1794 IF (SHDFAC.GT.0.0) THEN
1795
1796 CALL TRANSP (ET1,NSOIL,ETP1,SH2O,CMC,ZSOIL,SHDFAC,SMCWLT, &
1797 & CMCMAX,PC,CFACTR,SMCREF,SFCTMP,Q2,NROOT,RTDIS)
1798
1799 DO K = 1,NSOIL
1800 ETT1 = ETT1 + ET1(K)
1801 END DO
1802
1803 ! ----------------------------------------------------------------------
1804 ! CALCULATE CANOPY EVAPORATION.
1805 ! IF STATEMENTS TO AVOID TANGENT LINEAR PROBLEMS NEAR CMC=0.0.
1806 ! ----------------------------------------------------------------------
1807 IF (CMC .GT. 0.0) THEN
1808 EC1 = SHDFAC * ( ( CMC / CMCMAX ) ** CFACTR ) * ETP1
1809 ELSE
1810 EC1 = 0.0
1811 ENDIF
1812
1813 ! ----------------------------------------------------------------------
1814 ! EC SHOULD BE LIMITED BY THE TOTAL AMOUNT OF AVAILABLE WATER ON THE
1815 ! CANOPY. -F.CHEN, 18-OCT-1994
1816 ! ----------------------------------------------------------------------
1817 CMC2MS = CMC / DT
1818 EC1 = MIN ( CMC2MS, EC1 )
1819 ENDIF
1820 ENDIF
1821
1822 ! ----------------------------------------------------------------------
1823 ! TOTAL UP EVAP AND TRANSP TYPES TO OBTAIN ACTUAL EVAPOTRANSP
1824 ! ----------------------------------------------------------------------
1825 ETA1 = EDIR1 + ETT1 + EC1
1826
1827 ! ----------------------------------------------------------------------
1828 ! END SUBROUTINE EVAPO
1829 ! ----------------------------------------------------------------------
1830 END SUBROUTINE EVAPO
1831
1832 SUBROUTINE HRT (RHSTS,STC,SMC,SMCMAX,NSOIL,ZSOIL,YY,ZZ1, &
1833 & TBOT,ZBOT,PSISAT,SH2O,DT,BEXP, &
1834 & F1,DF1,QUARTZ,CSOIL,AI,BI,CI)
1835
1836 IMPLICIT NONE
1837
1838 ! ----------------------------------------------------------------------
1839 ! SUBROUTINE HRT
1840 ! ----------------------------------------------------------------------
1841 ! CALCULATE THE RIGHT HAND SIDE OF THE TIME TENDENCY TERM OF THE SOIL
1842 ! THERMAL DIFFUSION EQUATION. ALSO TO COMPUTE ( PREPARE ) THE MATRIX
1843 ! COEFFICIENTS FOR THE TRI-DIAGONAL MATRIX OF THE IMPLICIT TIME SCHEME.
1844 ! ----------------------------------------------------------------------
1845 INTEGER NSOLD
1846 PARAMETER(NSOLD = 20)
1847
1848 LOGICAL ITAVG
1849
1850 INTEGER I
1851 INTEGER K
1852 INTEGER NSOIL
1853
1854 ! ----------------------------------------------------------------------
1855 ! DECLARE WORK ARRAYS NEEDED IN TRI-DIAGONAL IMPLICIT SOLVER
1856 ! ----------------------------------------------------------------------
1857 REAL AI(NSOLD)
1858 REAL BI(NSOLD)
1859 REAL CI(NSOLD)
1860
1861 ! ----------------------------------------------------------------------
1862 ! DECLARATIONS
1863 ! ----------------------------------------------------------------------
1864 REAL BEXP
1865 REAL CAIR
1866 REAL CH2O
1867 REAL CICE
1868 REAL CSOIL
1869 REAL DDZ
1870 REAL DDZ2
1871 REAL DENOM
1872 REAL DF1
1873 REAL DF1N
1874 REAL DF1K
1875 REAL DT
1876 REAL DTSDZ
1877 REAL DTSDZ2
1878 REAL F1
1879 REAL HCPCT
1880 REAL PSISAT
1881 REAL QUARTZ
1882 REAL QTOT
1883 REAL RHSTS(NSOIL)
1884 REAL SSOIL
1885 REAL SICE
1886 REAL SMC(NSOIL)
1887 REAL SH2O(NSOIL)
1888 REAL SMCMAX
1889 ! REAL SNKSRC
1890 REAL STC(NSOIL)
1891 REAL T0
1892 REAL TAVG
1893 REAL TBK
1894 REAL TBK1
1895 REAL TBOT
1896 REAL ZBOT
1897 REAL TSNSR
1898 REAL TSURF
1899 REAL YY
1900 REAL ZSOIL(NSOIL)
1901 REAL ZZ1
1902
1903 PARAMETER(T0 = 273.15)
1904
1905 ! ----------------------------------------------------------------------
1906 ! SET SPECIFIC HEAT CAPACITIES OF AIR, WATER, ICE, SOIL MINERAL
1907 ! ----------------------------------------------------------------------
1908 PARAMETER(CAIR = 1004.0)
1909 PARAMETER(CH2O = 4.2E6)
1910 PARAMETER(CICE = 2.106E6)
1911 ! NOTE: CSOIL NOW SET IN ROUTINE REDPRM AND PASSED IN
1912 ! PARAMETER(CSOIL = 1.26E6)
1913
1914 ! ----------------------------------------------------------------------
1915 ! INITIALIZE LOGICAL FOR SOIL LAYER TEMPERATURE AVERAGING.
1916 ! ----------------------------------------------------------------------
1917 ITAVG = .TRUE.
1918 ! ITAVG = .FALSE.
1919
1920 ! ----------------------------------------------------------------------
1921 ! BEGIN SECTION FOR TOP SOIL LAYER
1922 ! ----------------------------------------------------------------------
1923 ! CALC THE HEAT CAPACITY OF THE TOP SOIL LAYER
1924 ! ----------------------------------------------------------------------
1925 HCPCT = SH2O(1)*CH2O + (1.0-SMCMAX)*CSOIL + (SMCMAX-SMC(1))*CAIR &
1926 & + ( SMC(1) - SH2O(1) )*CICE
1927
1928 ! ----------------------------------------------------------------------
1929 ! CALC THE MATRIX COEFFICIENTS AI, BI, AND CI FOR THE TOP LAYER
1930 ! ----------------------------------------------------------------------
1931 DDZ = 1.0 / ( -0.5 * ZSOIL(2) )
1932 AI(1) = 0.0
1933 CI(1) = (DF1 * DDZ) / (ZSOIL(1) * HCPCT)
1934 BI(1) = -CI(1) + DF1 / (0.5 * ZSOIL(1) * ZSOIL(1)*HCPCT*ZZ1)
1935
1936 ! ----------------------------------------------------------------------
1937 ! CALCULATE THE VERTICAL SOIL TEMP GRADIENT BTWN THE 1ST AND 2ND SOIL
1938 ! LAYERS. THEN CALCULATE THE SUBSURFACE HEAT FLUX. USE THE TEMP
1939 ! GRADIENT AND SUBSFC HEAT FLUX TO CALC "RIGHT-HAND SIDE TENDENCY
1940 ! TERMS", OR "RHSTS", FOR TOP SOIL LAYER.
1941 ! ----------------------------------------------------------------------
1942 DTSDZ = (STC(1) - STC(2)) / (-0.5 * ZSOIL(2))
1943 SSOIL = DF1 * (STC(1) - YY) / (0.5 * ZSOIL(1) * ZZ1)
1944 RHSTS(1) = (DF1 * DTSDZ - SSOIL) / (ZSOIL(1) * HCPCT)
1945
1946 ! ----------------------------------------------------------------------
1947 ! NEXT CAPTURE THE VERTICAL DIFFERENCE OF THE HEAT FLUX AT TOP AND
1948 ! BOTTOM OF FIRST SOIL LAYER FOR USE IN HEAT FLUX CONSTRAINT APPLIED TO
1949 ! POTENTIAL SOIL FREEZING/THAWING IN ROUTINE SNKSRC.
1950 ! ----------------------------------------------------------------------
1951 QTOT = SSOIL - DF1*DTSDZ
1952
1953 ! ----------------------------------------------------------------------
1954 ! IF TEMPERATURE AVERAGING INVOKED (ITAVG=TRUE; ELSE SKIP):
1955 ! SET TEMP "TSURF" AT TOP OF SOIL COLUMN (FOR USE IN FREEZING SOIL
1956 ! PHYSICS LATER IN FUNCTION SUBROUTINE SNKSRC). IF SNOWPACK CONTENT IS
1957 ! ZERO, THEN TSURF EXPRESSION BELOW GIVES TSURF = SKIN TEMP. IF
1958 ! SNOWPACK IS NONZERO (HENCE ARGUMENT ZZ1=1), THEN TSURF EXPRESSION
1959 ! BELOW YIELDS SOIL COLUMN TOP TEMPERATURE UNDER SNOWPACK. THEN
1960 ! CALCULATE TEMPERATURE AT BOTTOM INTERFACE OF 1ST SOIL LAYER FOR USE
1961 ! LATER IN FUNCTION SUBROUTINE SNKSRC
1962 ! ----------------------------------------------------------------------
1963 IF (ITAVG) THEN
1964 TSURF = (YY + (ZZ1-1) * STC(1)) / ZZ1
1965 CALL TBND (STC(1),STC(2),ZSOIL,ZBOT,1,NSOIL,TBK)
1966 ENDIF
1967
1968 ! ----------------------------------------------------------------------
1969 ! CALCULATE FROZEN WATER CONTENT IN 1ST SOIL LAYER.
1970 ! ----------------------------------------------------------------------
1971 SICE = SMC(1) - SH2O(1)
1972
1973 ! ----------------------------------------------------------------------
1974 ! IF FROZEN WATER PRESENT OR ANY OF LAYER-1 MID-POINT OR BOUNDING
1975 ! INTERFACE TEMPERATURES BELOW FREEZING, THEN CALL SNKSRC TO
1976 ! COMPUTE HEAT SOURCE/SINK (AND CHANGE IN FROZEN WATER CONTENT)
1977 ! DUE TO POSSIBLE SOIL WATER PHASE CHANGE
1978 ! ----------------------------------------------------------------------
1979 IF ( (SICE .GT. 0.) .OR. (TSURF .LT. T0) .OR. &
1980 & (STC(1) .LT. T0) .OR. (TBK .LT. T0) ) THEN
1981
1982 IF (ITAVG) THEN
1983 CALL TMPAVG(TAVG,TSURF,STC(1),TBK,ZSOIL,NSOIL,1)
1984 ELSE
1985 TAVG = STC(1)
1986 ENDIF
1987 TSNSR = SNKSRC (TAVG,SMC(1),SH2O(1), &
1988 & ZSOIL,NSOIL,SMCMAX,PSISAT,BEXP,DT,1,QTOT)
1989
1990 RHSTS(1) = RHSTS(1) - TSNSR / ( ZSOIL(1) * HCPCT )
1991 ENDIF
1992
1993 ! ----------------------------------------------------------------------
1994 ! THIS ENDS SECTION FOR TOP SOIL LAYER.
1995 ! ----------------------------------------------------------------------
1996 ! INITIALIZE DDZ2
1997 ! ----------------------------------------------------------------------
1998 DDZ2 = 0.0
1999
2000 ! ----------------------------------------------------------------------
2001 ! LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABOVE PROCESS
2002 ! (EXCEPT SUBSFC OR "GROUND" HEAT FLUX NOT REPEATED IN LOWER LAYERS)
2003 ! ----------------------------------------------------------------------
2004 DF1K = DF1
2005 DO K = 2,NSOIL
2006
2007 ! ----------------------------------------------------------------------
2008 ! CALCULATE HEAT CAPACITY FOR THIS SOIL LAYER.
2009 ! ----------------------------------------------------------------------
2010 HCPCT = SH2O(K)*CH2O +(1.0-SMCMAX)*CSOIL +(SMCMAX-SMC(K))*CAIR &
2011 & + ( SMC(K) - SH2O(K) )*CICE
2012
2013 IF (K .NE. NSOIL) THEN
2014 ! ----------------------------------------------------------------------
2015 ! THIS SECTION FOR LAYER 2 OR GREATER, BUT NOT LAST LAYER.
2016 ! ----------------------------------------------------------------------
2017 ! CALCULATE THERMAL DIFFUSIVITY FOR THIS LAYER.
2018 ! ----------------------------------------------------------------------
2019 CALL TDFCND (DF1N,SMC(K),QUARTZ,SMCMAX,SH2O(K))
2020
2021 ! ----------------------------------------------------------------------
2022 ! CALC THE VERTICAL SOIL TEMP GRADIENT THRU THIS LAYER
2023 ! ----------------------------------------------------------------------
2024 DENOM = 0.5 * ( ZSOIL(K-1) - ZSOIL(K+1) )
2025 DTSDZ2 = ( STC(K) - STC(K+1) ) / DENOM
2026
2027 ! ----------------------------------------------------------------------
2028 ! CALC THE MATRIX COEF, CI, AFTER CALC'NG ITS PARTIAL PRODUCT
2029 ! ----------------------------------------------------------------------
2030 DDZ2 = 2. / (ZSOIL(K-1) - ZSOIL(K+1))
2031 CI(K) = -DF1N * DDZ2 / ((ZSOIL(K-1) - ZSOIL(K)) * HCPCT)
2032
2033 ! ----------------------------------------------------------------------
2034 ! IF TEMPERATURE AVERAGING INVOKED (ITAVG=TRUE; ELSE SKIP): CALCULATE
2035 ! TEMP AT BOTTOM OF LAYER.
2036 ! ----------------------------------------------------------------------
2037 IF (ITAVG) THEN
2038 CALL TBND (STC(K),STC(K+1),ZSOIL,ZBOT,K,NSOIL,TBK1)
2039 ENDIF
2040 ELSE
2041
2042 ! ----------------------------------------------------------------------
2043 ! SPECIAL CASE OF BOTTOM SOIL LAYER: CALCULATE THERMAL DIFFUSIVITY FOR
2044 ! BOTTOM LAYER.
2045 ! ----------------------------------------------------------------------
2046 CALL TDFCND (DF1N,SMC(K),QUARTZ,SMCMAX,SH2O(K))
2047
2048 ! ----------------------------------------------------------------------
2049 ! CALC THE VERTICAL SOIL TEMP GRADIENT THRU BOTTOM LAYER.
2050 ! ----------------------------------------------------------------------
2051 DENOM = .5 * (ZSOIL(K-1) + ZSOIL(K)) - ZBOT
2052 DTSDZ2 = (STC(K)-TBOT) / DENOM
2053
2054 ! ----------------------------------------------------------------------
2055 ! SET MATRIX COEF, CI TO ZERO IF BOTTOM LAYER.
2056 ! ----------------------------------------------------------------------
2057 CI(K) = 0.
2058
2059 ! ----------------------------------------------------------------------
2060 ! IF TEMPERATURE AVERAGING INVOKED (ITAVG=TRUE; ELSE SKIP): CALCULATE
2061 ! TEMP AT BOTTOM OF LAST LAYER.
2062 ! ----------------------------------------------------------------------
2063 IF (ITAVG) THEN
2064 CALL TBND (STC(K),TBOT,ZSOIL,ZBOT,K,NSOIL,TBK1)
2065 ENDIF
2066
2067 ENDIF
2068 ! ----------------------------------------------------------------------
2069 ! THIS ENDS SPECIAL LOOP FOR BOTTOM LAYER.
2070 ! ----------------------------------------------------------------------
2071 ! CALCULATE RHSTS FOR THIS LAYER AFTER CALC'NG A PARTIAL PRODUCT.
2072 ! ----------------------------------------------------------------------
2073 DENOM = ( ZSOIL(K) - ZSOIL(K-1) ) * HCPCT
2074 RHSTS(K) = ( DF1N * DTSDZ2 - DF1K * DTSDZ ) / DENOM
2075 QTOT = -1.0*DENOM*RHSTS(K)
2076 SICE = SMC(K) - SH2O(K)
2077
2078 IF ( (SICE .GT. 0.) .OR. (TBK .LT. T0) .OR. &
2079 & (STC(K) .LT. T0) .OR. (TBK1 .LT. T0) ) THEN
2080
2081 IF (ITAVG) THEN
2082 CALL TMPAVG(TAVG,TBK,STC(K),TBK1,ZSOIL,NSOIL,K)
2083 ELSE
2084 TAVG = STC(K)
2085 ENDIF
2086 TSNSR = SNKSRC(TAVG,SMC(K),SH2O(K),ZSOIL,NSOIL, &
2087 & SMCMAX,PSISAT,BEXP,DT,K,QTOT)
2088 RHSTS(K) = RHSTS(K) - TSNSR / DENOM
2089 ENDIF
2090
2091 ! ----------------------------------------------------------------------
2092 ! CALC MATRIX COEFS, AI, AND BI FOR THIS LAYER.
2093 ! ----------------------------------------------------------------------
2094 AI(K) = - DF1 * DDZ / ((ZSOIL(K-1) - ZSOIL(K)) * HCPCT)
2095 BI(K) = -(AI(K) + CI(K))
2096
2097 ! ----------------------------------------------------------------------
2098 ! RESET VALUES OF DF1, DTSDZ, DDZ, AND TBK FOR LOOP TO NEXT SOIL LAYER.
2099 ! ----------------------------------------------------------------------
2100 TBK = TBK1
2101 DF1K = DF1N
2102 DTSDZ = DTSDZ2
2103 DDZ = DDZ2
2104 END DO
2105
2106 ! ----------------------------------------------------------------------
2107 ! END SUBROUTINE HRT
2108 ! ----------------------------------------------------------------------
2109 END SUBROUTINE HRT
2110
2111 SUBROUTINE HRTICE (RHSTS,STC,NSOIL,ZSOIL,YY,ZZ1,DF1,AI,BI,CI)
2112
2113 IMPLICIT NONE
2114
2115 ! ----------------------------------------------------------------------
2116 ! SUBROUTINE HRTICE
2117 ! ----------------------------------------------------------------------
2118 ! CALCULATE THE RIGHT HAND SIDE OF THE TIME TENDENCY TERM OF THE SOIL
2119 ! THERMAL DIFFUSION EQUATION IN THE CASE OF SEA-ICE PACK. ALSO TO
2120 ! COMPUTE (PREPARE) THE MATRIX COEFFICIENTS FOR THE TRI-DIAGONAL MATRIX
2121 ! OF THE IMPLICIT TIME SCHEME.
2122 ! ----------------------------------------------------------------------
2123 INTEGER NSOLD
2124 PARAMETER(NSOLD = 20)
2125
2126 INTEGER K
2127 INTEGER NSOIL
2128
2129 REAL AI(NSOLD)
2130 REAL BI(NSOLD)
2131 REAL CI(NSOLD)
2132
2133 REAL DDZ
2134 REAL DDZ2
2135 REAL DENOM
2136 REAL DF1
2137 REAL DTSDZ
2138 REAL DTSDZ2
2139 REAL HCPCT
2140 REAL RHSTS(NSOIL)
2141 REAL SSOIL
2142 REAL STC(NSOIL)
2143 REAL TBOT
2144 REAL YY
2145 REAL ZBOT
2146 REAL ZSOIL(NSOIL)
2147 REAL ZZ1
2148
2149 DATA TBOT /271.16/
2150
2151 ! ----------------------------------------------------------------------
2152 ! SET A NOMINAL UNIVERSAL VALUE OF THE SEA-ICE SPECIFIC HEAT CAPACITY,
2153 ! HCPCT = 1880.0*917.0.
2154 ! ----------------------------------------------------------------------
2155 PARAMETER(HCPCT = 1.72396E+6)
2156
2157 ! ----------------------------------------------------------------------
2158 ! THE INPUT ARGUMENT DF1 IS A UNIVERSALLY CONSTANT VALUE OF SEA-ICE
2159 ! THERMAL DIFFUSIVITY, SET IN ROUTINE SNOPAC AS DF1 = 2.2.
2160 ! ----------------------------------------------------------------------
2161 ! SET ICE PACK DEPTH. USE TBOT AS ICE PACK LOWER BOUNDARY TEMPERATURE
2162 ! (THAT OF UNFROZEN SEA WATER AT BOTTOM OF SEA ICE PACK). ASSUME ICE
2163 ! PACK IS OF N=NSOIL LAYERS SPANNING A UNIFORM CONSTANT ICE PACK
2164 ! THICKNESS AS DEFINED BY ZSOIL(NSOIL) IN ROUTINE SFLX.
2165 ! ----------------------------------------------------------------------
2166 ZBOT = ZSOIL(NSOIL)
2167
2168 ! ----------------------------------------------------------------------
2169 ! CALC THE MATRIX COEFFICIENTS AI, BI, AND CI FOR THE TOP LAYER
2170 ! ----------------------------------------------------------------------
2171 DDZ = 1.0 / ( -0.5 * ZSOIL(2) )
2172 AI(1) = 0.0
2173 CI(1) = (DF1 * DDZ) / (ZSOIL(1) * HCPCT)
2174 BI(1) = -CI(1) + DF1/(0.5 * ZSOIL(1) * ZSOIL(1) * HCPCT * ZZ1)
2175
2176 ! ----------------------------------------------------------------------
2177 ! CALC THE VERTICAL SOIL TEMP GRADIENT BTWN THE TOP AND 2ND SOIL LAYERS.
2178 ! RECALC/ADJUST THE SOIL HEAT FLUX. USE THE GRADIENT AND FLUX TO CALC
2179 ! RHSTS FOR THE TOP SOIL LAYER.
2180 ! ----------------------------------------------------------------------
2181 DTSDZ = ( STC(1) - STC(2) ) / ( -0.5 * ZSOIL(2) )
2182 SSOIL = DF1 * ( STC(1) - YY ) / ( 0.5 * ZSOIL(1) * ZZ1 )
2183 RHSTS(1) = ( DF1 * DTSDZ - SSOIL ) / ( ZSOIL(1) * HCPCT )
2184
2185 ! ----------------------------------------------------------------------
2186 ! INITIALIZE DDZ2
2187 ! ----------------------------------------------------------------------
2188 DDZ2 = 0.0
2189
2190 ! ----------------------------------------------------------------------
2191 ! LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABOVE PROCESS
2192 ! ----------------------------------------------------------------------
2193 DO K = 2,NSOIL
2194 IF (K .NE. NSOIL) THEN
2195
2196 ! ----------------------------------------------------------------------
2197 ! CALC THE VERTICAL SOIL TEMP GRADIENT THRU THIS LAYER.
2198 ! ----------------------------------------------------------------------
2199 DENOM = 0.5 * ( ZSOIL(K-1) - ZSOIL(K+1) )
2200 DTSDZ2 = ( STC(K) - STC(K+1) ) / DENOM
2201
2202 ! ----------------------------------------------------------------------
2203 ! CALC THE MATRIX COEF, CI, AFTER CALC'NG ITS PARTIAL PRODUCT.
2204 ! ----------------------------------------------------------------------
2205 DDZ2 = 2. / (ZSOIL(K-1) - ZSOIL(K+1))
2206 CI(K) = -DF1 * DDZ2 / ((ZSOIL(K-1) - ZSOIL(K)) * HCPCT)
2207 ELSE
2208
2209 ! ----------------------------------------------------------------------
2210 ! CALC THE VERTICAL SOIL TEMP GRADIENT THRU THE LOWEST LAYER.
2211 ! ----------------------------------------------------------------------
2212 DTSDZ2 = (STC(K)-TBOT)/(.5 * (ZSOIL(K-1) + ZSOIL(K))-ZBOT)
2213
2214 ! ----------------------------------------------------------------------
2215 ! SET MATRIX COEF, CI TO ZERO.
2216 ! ----------------------------------------------------------------------
2217 CI(K) = 0.
2218 ENDIF
2219
2220 ! ----------------------------------------------------------------------
2221 ! CALC RHSTS FOR THIS LAYER AFTER CALC'NG A PARTIAL PRODUCT.
2222 ! ----------------------------------------------------------------------
2223 DENOM = ( ZSOIL(K) - ZSOIL(K-1) ) * HCPCT
2224 RHSTS(K) = ( DF1 * DTSDZ2 - DF1 * DTSDZ ) / DENOM
2225
2226 ! ----------------------------------------------------------------------
2227 ! CALC MATRIX COEFS, AI, AND BI FOR THIS LAYER.
2228 ! ----------------------------------------------------------------------
2229 AI(K) = - DF1 * DDZ / ((ZSOIL(K-1) - ZSOIL(K)) * HCPCT)
2230 BI(K) = -(AI(K) + CI(K))
2231
2232 ! ----------------------------------------------------------------------
2233 ! RESET VALUES OF DTSDZ AND DDZ FOR LOOP TO NEXT SOIL LYR.
2234 ! ----------------------------------------------------------------------
2235 DTSDZ = DTSDZ2
2236 DDZ = DDZ2
2237
2238 END DO
2239 ! ----------------------------------------------------------------------
2240 ! END SUBROUTINE HRTICE
2241 ! ----------------------------------------------------------------------
2242 END SUBROUTINE HRTICE
2243
2244 SUBROUTINE HSTEP (STCOUT,STCIN,RHSTS,DT,NSOIL,AI,BI,CI)
2245
2246 IMPLICIT NONE
2247
2248 ! ----------------------------------------------------------------------
2249 ! SUBROUTINE HSTEP
2250 ! ----------------------------------------------------------------------
2251 ! CALCULATE/UPDATE THE SOIL TEMPERATURE FIELD.
2252 ! ----------------------------------------------------------------------
2253 INTEGER NSOLD
2254 PARAMETER(NSOLD = 20)
2255
2256 INTEGER K
2257 INTEGER NSOIL
2258
2259 REAL AI(NSOLD)
2260 REAL BI(NSOLD)
2261 REAL CI(NSOLD)
2262 REAL CIin(NSOLD)
2263 REAL DT
2264 REAL RHSTS(NSOIL)
2265 REAL RHSTSin(NSOIL)
2266 REAL STCIN(NSOIL)
2267 REAL STCOUT(NSOIL)
2268
2269 ! ----------------------------------------------------------------------
2270 ! CREATE FINITE DIFFERENCE VALUES FOR USE IN ROSR12 ROUTINE
2271 ! ----------------------------------------------------------------------
2272 DO K = 1,NSOIL
2273 RHSTS(K) = RHSTS(K) * DT
2274 AI(K) = AI(K) * DT
2275 BI(K) = 1. + BI(K) * DT
2276 CI(K) = CI(K) * DT
2277 END DO
2278
2279 ! ----------------------------------------------------------------------
2280 ! COPY VALUES FOR INPUT VARIABLES BEFORE CALL TO ROSR12
2281 ! ----------------------------------------------------------------------
2282 DO K = 1,NSOIL
2283 RHSTSin(K) = RHSTS(K)
2284 END DO
2285 DO K = 1,NSOIL
2286 CIin(K) = CI(K)
2287 END DO
2288
2289 ! ----------------------------------------------------------------------
2290 ! SOLVE THE TRI-DIAGONAL MATRIX EQUATION
2291 ! ----------------------------------------------------------------------
2292 CALL ROSR12(CI,AI,BI,CIin,RHSTSin,RHSTS,NSOIL)
2293
2294 ! ----------------------------------------------------------------------
2295 ! CALC/UPDATE THE SOIL TEMPS USING MATRIX SOLUTION
2296 ! ----------------------------------------------------------------------
2297 DO K = 1,NSOIL
2298 STCOUT(K) = STCIN(K) + CI(K)
2299 END DO
2300
2301 ! ----------------------------------------------------------------------
2302 ! END SUBROUTINE HSTEP
2303 ! ----------------------------------------------------------------------
2304 END SUBROUTINE HSTEP
2305
2306 SUBROUTINE NOPAC(ETP,ETA,PRCP,SMC,SMCMAX,SMCWLT, &
2307 & SMCREF,SMCDRY,CMC,CMCMAX,NSOIL,DT,SHDFAC, &
2308 & SBETA,Q2,T1,SFCTMP,T24,TH2,FDOWN,F1,SSOIL, &
2309 & STC,EPSCA,BEXP,PC,RCH,RR,CFACTR, &
2310 & SH2O,SLOPE,KDT,FRZFACT,PSISAT,ZSOIL, &
2311 & DKSAT,DWSAT,TBOT,ZBOT,RUNOFF1,RUNOFF2, &
2312 & RUNOFF3,EDIR,EC,ET,ETT,NROOT,ICE,RTDIS, &
2313 & QUARTZ,FXEXP,CSOIL, &
2314 & BETA,DRIP,DEW,FLX1,FLX2,FLX3)
2315
2316 IMPLICIT NONE
2317
2318 ! ----------------------------------------------------------------------
2319 ! SUBROUTINE NOPAC
2320 ! ----------------------------------------------------------------------
2321 ! CALCULATE SOIL MOISTURE AND HEAT FLUX VALUES AND UPDATE SOIL MOISTURE
2322 ! CONTENT AND SOIL HEAT CONTENT VALUES FOR THE CASE WHEN NO SNOW PACK IS
2323 ! PRESENT.
2324 ! ----------------------------------------------------------------------
2325 INTEGER ICE
2326 INTEGER NROOT
2327 INTEGER NSOIL
2328
2329 REAL BEXP
2330 REAL BETA
2331 REAL CFACTR
2332 REAL CMC
2333 REAL CMCMAX
2334 REAL CP
2335 REAL CSOIL
2336 REAL DEW
2337 REAL DF1
2338 REAL DKSAT
2339 REAL DRIP
2340 REAL DT
2341 REAL DWSAT
2342 REAL EC
2343 REAL EDIR
2344 REAL EPSCA
2345 REAL ETA
2346 REAL ETA1
2347 REAL ETP
2348 REAL ETP1
2349 REAL ET(NSOIL)
2350 REAL ETT
2351 REAL FDOWN
2352 REAL F1
2353 REAL FXEXP
2354 REAL FLX1
2355 REAL FLX2
2356 REAL FLX3
2357 REAL FRZFACT
2358 REAL KDT
2359 REAL PC
2360 REAL PRCP
2361 REAL PRCP1
2362 REAL PSISAT
2363 REAL Q2
2364 REAL QUARTZ
2365 REAL RCH
2366 REAL RR
2367 REAL RTDIS(NSOIL)
2368 REAL RUNOFF1
2369 REAL RUNOFF2
2370 REAL RUNOFF3
2371 REAL SSOIL
2372 REAL SBETA
2373 REAL SFCTMP
2374 REAL SHDFAC
2375 REAL SH2O(NSOIL)
2376 REAL SIGMA
2377 REAL SLOPE
2378 REAL SMC(NSOIL)
2379 REAL SMCDRY
2380 REAL SMCMAX
2381 REAL SMCREF
2382 REAL SMCWLT
2383 REAL STC(NSOIL)
2384 REAL T1
2385 REAL T24
2386 REAL TBOT
2387 REAL TH2
2388 REAL YY
2389 REAL YYNUM
2390 REAL ZBOT
2391 REAL ZSOIL(NSOIL)
2392 REAL ZZ1
2393
2394 REAL EC1
2395 REAL EDIR1
2396 REAL ET1(NSOIL)
2397 REAL ETT1
2398
2399 INTEGER K
2400
2401 PARAMETER(CP = 1004.5)
2402 PARAMETER(SIGMA = 5.67E-8)
2403
2404 ! ----------------------------------------------------------------------
2405 ! EXECUTABLE CODE BEGINS HERE:
2406 ! CONVERT ETP FROM KG M-2 S-1 TO MS-1 AND INITIALIZE DEW.
2407 ! ----------------------------------------------------------------------
2408 PRCP1 = PRCP * 0.001
2409 ETP1 = ETP * 0.001
2410 DEW = 0.0
2411
2412 EDIR = 0.
2413 EDIR1 = 0.
2414 EC = 0.
2415 EC1 = 0.
2416 DO K = 1,NSOIL
2417 ET(K) = 0.
2418 ET1(K) = 0.
2419 END DO
2420 ETT = 0.
2421 ETT1 = 0.
2422
2423 IF (ETP .GT. 0.0) THEN
2424
2425 ! ----------------------------------------------------------------------
2426 ! CONVERT PRCP FROM 'KG M-2 S-1' TO 'M S-1'.
2427 ! ----------------------------------------------------------------------
2428 CALL EVAPO (ETA1,SMC,NSOIL,CMC,ETP1,DT,ZSOIL, &
2429 & SH2O,SMCMAX,BEXP,PC,SMCWLT,DKSAT,DWSAT, &
2430 & SMCREF,SHDFAC,CMCMAX, &
2431 & SMCDRY,CFACTR, &
2432 & EDIR1,EC1,ET1,ETT1,SFCTMP,Q2,NROOT,RTDIS,FXEXP)
2433 CALL SMFLX (SMC,NSOIL,CMC,DT,PRCP1,ZSOIL, &
2434 & SH2O,SLOPE,KDT,FRZFACT, &
2435 & SMCMAX,BEXP,SMCWLT,DKSAT,DWSAT, &
2436 & SHDFAC,CMCMAX, &
2437 & RUNOFF1,RUNOFF2,RUNOFF3, &
2438 & EDIR1,EC1,ET1, &
2439 & DRIP)
2440
2441 ! ----------------------------------------------------------------------
2442 ! CONVERT MODELED EVAPOTRANSPIRATION FM M S-1 TO KG M-2 S-1
2443 ! ----------------------------------------------------------------------
2444 ETA = ETA1 * 1000.0
2445
2446 ! ----------------------------------------------------------------------
2447 ! EDIR = EDIR1 * 1000.0
2448 ! EC = EC1 * 1000.0
2449 ! ETT = ETT1 * 1000.0
2450 ! ET(1) = ET1(1) * 1000.0
2451 ! ET(2) = ET1(2) * 1000.0
2452 ! ET(3) = ET1(3) * 1000.0
2453 ! ET(4) = ET1(4) * 1000.0
2454 ! ----------------------------------------------------------------------
2455
2456 ELSE
2457
2458 ! ----------------------------------------------------------------------
2459 ! IF ETP < 0, ASSUME DEW FORMS (TRANSFORM ETP1 INTO DEW AND REINITIALIZE
2460 ! ETP1 TO ZERO).
2461 ! ----------------------------------------------------------------------
2462 DEW = -ETP1
2463 ! ETP1 = 0.0
2464
2465 ! ----------------------------------------------------------------------
2466 ! CONVERT PRCP FROM 'KG M-2 S-1' TO 'M S-1' AND ADD DEW AMOUNT.
2467 ! ----------------------------------------------------------------------
2468 PRCP1 = PRCP1 + DEW
2469 !
2470 ! CALL EVAPO (ETA1,SMC,NSOIL,CMC,ETP1,DT,ZSOIL,
2471 ! & SH2O,SMCMAX,BEXP,PC,SMCWLT,DKSAT,DWSAT,
2472 ! & SMCREF,SHDFAC,CMCMAX,
2473 ! & SMCDRY,CFACTR,
2474 ! & EDIR1,EC1,ET1,ETT,SFCTMP,Q2,NROOT,RTDIS,FXEXP)
2475 CALL SMFLX (SMC,NSOIL,CMC,DT,PRCP1,ZSOIL, &
2476 & SH2O,SLOPE,KDT,FRZFACT, &
2477 & SMCMAX,BEXP,SMCWLT,DKSAT,DWSAT, &
2478 & SHDFAC,CMCMAX, &
2479 & RUNOFF1,RUNOFF2,RUNOFF3, &
2480 & EDIR1,EC1,ET1, &
2481 & DRIP)
2482
2483 ! ----------------------------------------------------------------------
2484 ! CONVERT MODELED EVAPOTRANSPIRATION FROM 'M S-1' TO 'KG M-2 S-1'.
2485 ! ----------------------------------------------------------------------
2486 ! ETA = ETA1 * 1000.0
2487
2488 ! ----------------------------------------------------------------------
2489 ! EDIR = EDIR1 * 1000.0
2490 ! EC = EC1 * 1000.0
2491 ! ETT = ETT1 * 1000.0
2492 ! ET(1) = ET1(1) * 1000.0
2493 ! ET(2) = ET1(2) * 1000.0
2494 ! ET(3) = ET1(3) * 1000.0
2495 ! ET(4) = ET1(4) * 1000.0
2496 ! ----------------------------------------------------------------------
2497
2498 ENDIF
2499
2500 ! ----------------------------------------------------------------------
2501 ! CONVERT MODELED EVAPOTRANSPIRATION FM M S-1 TO KG M-2 S-1
2502 ! ----------------------------------------------------------------------
2503 ! ETA = ETA1 * 1000.0
2504
2505 ! ----------------------------------------------------------------------
2506 EDIR = EDIR1 * 1000.0
2507 EC = EC1 * 1000.0
2508 DO K = 1,NSOIL
2509 ET(K) = ET1(K) * 1000.0
2510 ! ET(1) = ET1(1) * 1000.0
2511 ! ET(2) = ET1(2) * 1000.0
2512 ! ET(3) = ET1(3) * 1000.0
2513 ! ET(4) = ET1(4) * 1000.0
2514 ENDDO
2515 ETT = ETT1 * 1000.0
2516 ! ----------------------------------------------------------------------
2517
2518 ! ----------------------------------------------------------------------
2519 ! BASED ON ETP AND E VALUES, DETERMINE BETA
2520 ! ----------------------------------------------------------------------
2521 IF ( ETP .LE. 0.0 ) THEN
2522 BETA = 0.0
2523 IF ( ETP .LT. 0.0 ) THEN
2524 BETA = 1.0
2525 ETA = ETP
2526 ENDIF
2527 ELSE
2528 BETA = ETA / ETP
2529 ENDIF
2530
2531 ! ----------------------------------------------------------------------
2532 ! GET SOIL THERMAL DIFFUXIVITY/CONDUCTIVITY FOR TOP SOIL LYR,
2533 ! CALC. ADJUSTED TOP LYR SOIL TEMP AND ADJUSTED SOIL FLUX, THEN
2534 ! CALL SHFLX TO COMPUTE/UPDATE SOIL HEAT FLUX AND SOIL TEMPS.
2535 ! ----------------------------------------------------------------------
2536 CALL TDFCND (DF1,SMC(1),QUARTZ,SMCMAX,SH2O(1))
2537
2538 ! ----------------------------------------------------------------------
2539 ! VEGETATION GREENNESS FRACTION REDUCTION IN SUBSURFACE HEAT FLUX
2540 ! VIA REDUCTION FACTOR, WHICH IS CONVENIENT TO APPLY HERE TO THERMAL
2541 ! DIFFUSIVITY THAT IS LATER USED IN HRT TO COMPUTE SUB SFC HEAT FLUX
2542 ! (SEE ADDITIONAL COMMENTS ON VEG EFFECT SUB-SFC HEAT FLX IN
2543 ! ROUTINE SFLX)
2544 ! ----------------------------------------------------------------------
2545 DF1 = DF1 * EXP(SBETA*SHDFAC)
2546
2547 ! ----------------------------------------------------------------------
2548 ! COMPUTE INTERMEDIATE TERMS PASSED TO ROUTINE HRT (VIA ROUTINE
2549 ! SHFLX BELOW) FOR USE IN COMPUTING SUBSURFACE HEAT FLUX IN HRT
2550 ! ----------------------------------------------------------------------
2551 YYNUM = FDOWN - SIGMA * T24
2552 YY = SFCTMP + (YYNUM/RCH+TH2-SFCTMP-BETA*EPSCA) / RR
2553 ZZ1 = DF1 / ( -0.5 * ZSOIL(1) * RCH * RR ) + 1.0
2554
2555 CALL SHFLX (SSOIL,STC,SMC,SMCMAX,NSOIL,T1,DT,YY,ZZ1,ZSOIL, &
2556 & TBOT,ZBOT,SMCWLT,PSISAT,SH2O,BEXP,F1,DF1,ICE, &
2557 & QUARTZ,CSOIL)
2558
2559 ! ----------------------------------------------------------------------
2560 ! SET FLX1 AND FLX3 (SNOPACK PHASE CHANGE HEAT FLUXES) TO ZERO SINCE
2561 ! THEY ARE NOT USED HERE IN SNOPAC. FLX2 (FREEZING RAIN HEAT FLUX) WAS
2562 ! SIMILARLY INITIALIZED IN THE PENMAN ROUTINE.
2563 ! ----------------------------------------------------------------------
2564 FLX1 = 0.0
2565 FLX3 = 0.0
2566
2567 ! ----------------------------------------------------------------------
2568 ! END SUBROUTINE NOPAC
2569 ! ----------------------------------------------------------------------
2570 END SUBROUTINE NOPAC
2571
2572 SUBROUTINE PENMAN (SFCTMP,SFCPRS,CH,T2V,TH2,PRCP,FDOWN,T24,SSOIL, &
2573 & Q2,Q2SAT,ETP,RCH,EPSCA,RR,SNOWNG,FRZGRA, &
2574 & DQSDT2,FLX2)
2575
2576 IMPLICIT NONE
2577
2578 ! ----------------------------------------------------------------------
2579 ! SUBROUTINE PENMAN
2580 ! ----------------------------------------------------------------------
2581 ! CALCULATE POTENTIAL EVAPORATION FOR THE CURRENT POINT. VARIOUS
2582 ! PARTIAL SUMS/PRODUCTS ARE ALSO CALCULATED AND PASSED BACK TO THE
2583 ! CALLING ROUTINE FOR LATER USE.
2584 ! ----------------------------------------------------------------------
2585 LOGICAL SNOWNG
2586 LOGICAL FRZGRA
2587
2588 REAL A
2589 REAL BETA
2590 REAL CH
2591 REAL CP
2592 REAL CPH2O
2593 REAL CPICE
2594 REAL DELTA
2595 REAL DQSDT2
2596 REAL ELCP
2597 REAL EPSCA
2598 REAL ETP
2599 REAL FDOWN
2600 REAL FLX2
2601 REAL FNET
2602 REAL LSUBC
2603 REAL LSUBF
2604 REAL PRCP
2605 REAL Q2
2606 REAL Q2SAT
2607 REAL R
2608 REAL RAD
2609 REAL RCH
2610 REAL RHO
2611 REAL RR
2612 REAL SSOIL
2613 REAL SFCPRS
2614 REAL SFCTMP
2615 REAL SIGMA
2616 REAL T24
2617 REAL T2V
2618 REAL TH2
2619
2620 PARAMETER(CP = 1004.6)
2621 PARAMETER(CPH2O = 4.218E+3)
2622 PARAMETER(CPICE = 2.106E+3)
2623 PARAMETER(R = 287.04)
2624 PARAMETER(ELCP = 2.4888E+3)
2625 PARAMETER(LSUBF = 3.335E+5)
2626 PARAMETER(LSUBC = 2.501000E+6)
2627 PARAMETER(SIGMA = 5.67E-8)
2628
2629 ! ----------------------------------------------------------------------
2630 ! EXECUTABLE CODE BEGINS HERE:
2631 ! ----------------------------------------------------------------------
2632 FLX2 = 0.0
2633
2634 ! ----------------------------------------------------------------------
2635 ! PREPARE PARTIAL QUANTITIES FOR PENMAN EQUATION.
2636 ! ----------------------------------------------------------------------
2637 DELTA = ELCP * DQSDT2
2638 T24 = SFCTMP * SFCTMP * SFCTMP * SFCTMP
2639 RR = T24 * 6.48E-8 /(SFCPRS * CH) + 1.0
2640 RHO = SFCPRS / (R * T2V)
2641 RCH = RHO * CP * CH
2642
2643 ! ----------------------------------------------------------------------
2644 ! ADJUST THE PARTIAL SUMS / PRODUCTS WITH THE LATENT HEAT
2645 ! EFFECTS CAUSED BY FALLING PRECIPITATION.
2646 ! ----------------------------------------------------------------------
2647 IF (.NOT. SNOWNG) THEN
2648 IF (PRCP .GT. 0.0) RR = RR + CPH2O*PRCP/RCH
2649 ELSE
2650 RR = RR + CPICE*PRCP/RCH
2651 ENDIF
2652
2653 FNET = FDOWN - SIGMA*T24 - SSOIL
2654
2655 ! ----------------------------------------------------------------------
2656 ! INCLUDE THE LATENT HEAT EFFECTS OF FRZNG RAIN CONVERTING TO ICE ON
2657 ! IMPACT IN THE CALCULATION OF FLX2 AND FNET.
2658 ! ----------------------------------------------------------------------
2659 IF (FRZGRA) THEN
2660 FLX2 = -LSUBF * PRCP
2661 FNET = FNET - FLX2
2662 ENDIF
2663
2664 ! ----------------------------------------------------------------------
2665 ! FINISH PENMAN EQUATION CALCULATIONS.
2666 ! ----------------------------------------------------------------------
2667 RAD = FNET/RCH + TH2 - SFCTMP
2668 A = ELCP * (Q2SAT - Q2)
2669 EPSCA = (A*RR + RAD*DELTA) / (DELTA + RR)
2670 ETP = EPSCA * RCH / LSUBC
2671
2672 ! ----------------------------------------------------------------------
2673 ! END SUBROUTINE PENMAN
2674 ! ----------------------------------------------------------------------
2675 END SUBROUTINE PENMAN
2676
2677 SUBROUTINE REDPRM ( &
2678 & VEGTYP,SOILTYP,SLOPETYP, &
2679 & CFACTR,CMCMAX,RSMAX,TOPT,REFKDT,KDT,SBETA, &
2680 & SHDFAC,RSMIN,RGL,HS,ZBOT,FRZX,PSISAT,SLOPE, &
2681 & SNUP,SALP,BEXP,DKSAT,DWSAT, &
2682 & SMCMAX,SMCWLT,SMCREF, &
2683 & SMCDRY,F1,QUARTZ,FXEXP,RTDIS,SLDPTH,ZSOIL, &
2684 & NROOT,NSOIL,Z0,CZIL,LAI,CSOIL,PTU)
2685
2686 IMPLICIT NONE
2687 ! ----------------------------------------------------------------------
2688 ! SUBROUTINE REDPRM
2689 ! ----------------------------------------------------------------------
2690 ! INTERNALLY SET (DEFAULT VALUESS), OR OPTIONALLY READ-IN VIA NAMELIST
2691 ! I/O, ALL SOIL AND VEGETATION PARAMETERS REQUIRED FOR THE EXECUSION OF
2692 ! THE NOAH LSM.
2693 !
2694 ! OPTIONAL NON-DEFAULT PARAMETERS CAN BE READ IN, ACCOMMODATING UP TO 30
2695 ! SOIL, VEG, OR SLOPE CLASSES, IF THE DEFAULT MAX NUMBER OF SOIL, VEG,
2696 ! AND/OR SLOPE TYPES IS RESET.
2697 !
2698 ! FUTURE UPGRADES OF ROUTINE REDPRM MUST EXPAND TO INCORPORATE SOME OF
2699 ! THE EMPIRICAL PARAMETERS OF THE FROZEN SOIL AND SNOWPACK PHYSICS (SUCH
2700 ! AS IN ROUTINES FRH2O, SNOWPACK, AND SNOW_NEW) NOT YET SET IN THIS
2701 ! REDPRM ROUTINE, BUT RATHER SET IN LOWER LEVEL SUBROUTINES.
2702 !
2703 ! SET MAXIMUM NUMBER OF SOIL-, VEG-, AND SLOPETYP IN DATA STATEMENT.
2704 ! ----------------------------------------------------------------------
2705 INTEGER MAX_SLOPETYP
2706 INTEGER MAX_SOILTYP
2707 INTEGER MAX_VEGTYP
2708
2709 PARAMETER(MAX_SLOPETYP = 30)
2710 PARAMETER(MAX_SOILTYP = 30)
2711 PARAMETER(MAX_VEGTYP = 30)
2712
2713 ! ----------------------------------------------------------------------
2714 ! NUMBER OF DEFINED SOIL-, VEG-, AND SLOPETYPS USED.
2715 ! ----------------------------------------------------------------------
2716 INTEGER DEFINED_VEG
2717 INTEGER DEFINED_SOIL
2718 INTEGER DEFINED_SLOPE
2719
2720 DATA DEFINED_VEG/27/
2721 DATA DEFINED_SOIL/19/
2722 DATA DEFINED_SLOPE/9/
2723
2724 ! ----------------------------------------------------------------------
2725 ! SET-UP SOIL PARAMETERS FOR GIVEN SOIL TYPE
2726 ! INPUT: SOLTYP: SOIL TYPE (INTEGER INDEX)
2727 ! OUTPUT: SOIL PARAMETERS:
2728 ! MAXSMC: MAX SOIL MOISTURE CONTENT (POROSITY)
2729 ! REFSMC: REFERENCE SOIL MOISTURE (ONSET OF SOIL MOISTURE
2730 ! STRESS IN TRANSPIRATION)
2731 ! WLTSMC: WILTING PT SOIL MOISTURE CONTENTS
2732 ! DRYSMC: AIR DRY SOIL MOIST CONTENT LIMITS
2733 ! SATPSI: SATURATED SOIL POTENTIAL
2734 ! SATDK: SATURATED SOIL HYDRAULIC CONDUCTIVITY
2735 ! BB: THE 'B' PARAMETER
2736 ! SATDW: SATURATED SOIL DIFFUSIVITY
2737 ! F11: USED TO COMPUTE SOIL DIFFUSIVITY/CONDUCTIVITY
2738 ! QUARTZ: SOIL QUARTZ CONTENT
2739 ! ----------------------------------------------------------------------
2740 ! SOIL STATSGO
2741 ! TYPE CLASS
2742 ! ---- -------
2743 ! 1 SAND
2744 ! 2 LOAMY SAND
2745 ! 3 SANDY LOAM
2746 ! 4 SILT LOAM
2747 ! 5 SILT
2748 ! 6 LOAM
2749 ! 7 SANDY CLAY LOAM
2750 ! 8 SILTY CLAY LOAM
2751 ! 9 CLAY LOAM
2752 ! 10 SANDY CLAY
2753 ! 11 SILTY CLAY
2754 ! 12 CLAY
2755 ! 13 ORGANIC MATERIAL
2756 ! 14 WATER
2757 ! 15 BEDROCK
2758 ! 16 OTHER(land-ice)
2759 ! 17 PLAYA
2760 ! 18 LAVA
2761 ! 19 WHITE SAND
2762 ! ----------------------------------------------------------------------
2763
2764 REAL BB(MAX_SOILTYP)
2765 REAL DRYSMC(MAX_SOILTYP)
2766 REAL F11(MAX_SOILTYP)
2767 REAL MAXSMC(MAX_SOILTYP)
2768 REAL REFSMC(MAX_SOILTYP)
2769 REAL SATPSI(MAX_SOILTYP)
2770 REAL SATDK(MAX_SOILTYP)
2771 REAL SATDW(MAX_SOILTYP)
2772 REAL WLTSMC(MAX_SOILTYP)
2773 REAL QTZ(MAX_SOILTYP)
2774
2775 REAL BEXP
2776 REAL DKSAT
2777 REAL DWSAT
2778 REAL F1
2779 REAL PTU
2780 REAL QUARTZ
2781 REAL REFSMC1
2782 REAL SMCDRY
2783 REAL SMCMAX
2784 REAL SMCREF
2785 REAL SMCWLT
2786 REAL WLTSMC1
2787
2788 ! ----------------------------------------------------------------------
2789 ! SOIL TEXTURE-RELATED ARRAYS.
2790 ! ----------------------------------------------------------------------
2791 DATA MAXSMC/0.395, 0.421, 0.434, 0.476, 0.476, 0.439, &
2792 & 0.404, 0.464, 0.465, 0.406, 0.468, 0.457, &
2793 & 0.464, 0.464, 0.200, 0.421, 0.457, 0.200, &
2794 & 0.395, 0.000, 0.000, 0.000, 0.000, 0.000, &
2795 & 0.000, 0.000, 0.000, 0.000, 0.000, 0.000/
2796 ! ----------------------------------------------------------------------
2797 DATA SATPSI/0.0350, 0.0363, 0.1413, 0.7586, 0.7586, 0.3548, &
2798 & 0.1349, 0.6166, 0.2630, 0.0977, 0.3236, 0.4677, &
2799 & 0.3548, 0.3548, 0.0350, 0.0363, 0.4677, 0.0350, &
2800 & 0.0350, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, &
2801 & 0.000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000/
2802 ! ----------------------------------------------------------------------
2803 DATA SATDK /1.7600E-4, 1.4078E-5, 5.2304E-6, 2.8089E-6, 2.8089E-6,&
2804 & 3.3770E-6, 4.4518E-6, 2.0348E-6, 2.4464E-6, 7.2199E-6,&
2805 & 1.3444E-6, 9.7394E-7, 3.3770E-6, 3.3770E-6, 1.4078E-5,&
2806 & 1.4078E-5, 9.7394E-7, 1.4078E-5, 1.7600E-4, 0.0,&
2807 & 0.0, 0.0, 0.0, 0.0, 0.0,&
2808 & 0.0, 0.0, 0.0, 0.0, 0.0/
2809 ! ----------------------------------------------------------------------
2810 DATA BB /4.05, 4.26, 4.74, 5.33, 5.33, 5.25, &
2811 & 6.77, 8.72, 8.17, 10.73, 10.39, 11.55, &
2812 & 5.25, 5.25, 4.05, 4.26, 11.55, 4.05, &
2813 & 4.05, 0.00, 0.00, 0.00, 0.00, 0.00, &
2814 & 0.00, 0.00, 0.00, 0.00, 0.00, 0.00/
2815 ! ----------------------------------------------------------------------
2816 DATA QTZ /0.92, 0.82, 0.60, 0.25, 0.10, 0.40, &
2817 & 0.60, 0.10, 0.35, 0.52, 0.10, 0.25, &
2818 & 0.05, 0.05, 0.07, 0.25, 0.60, 0.52, &
2819 & 0.92, 0.00, 0.00, 0.00, 0.00, 0.00, &
2820 & 0.00, 0.00, 0.00, 0.00, 0.00, 0.00/
2821
2822 ! ----------------------------------------------------------------------
2823 ! THE FOLLOWING 5 PARAMETERS ARE DERIVED LATER IN REDPRM.F FROM THE SOIL
2824 ! DATA, AND ARE JUST GIVEN HERE FOR REFERENCE AND TO FORCE STATIC
2825 ! STORAGE ALLOCATION. -DAG LOHMANN, FEB. 2001
2826 ! ----------------------------------------------------------------------
2827 ! !!!!!!!!!!!!!! The following values in the table are NOT used
2828 ! !!!!!!!!!!!!!! and are just given for reference
2829 DATA REFSMC/0.196, 0.248, 0.282, 0.332, 0.332, 0.301, &
2830 & 0.293, 0.368, 0.361, 0.320, 0.388, 0.389, &
2831 & 0.319, 0.000, 0.116, 0.248, 0.389, 0.116, &
2832 & 0.196, 0.000, 0.000, 0.000, 0.000, 0.000, &
2833 & 0.000, 0.000, 0.000, 0.000, 0.000, 0.000/
2834 ! !!!!!!!!!!!!!! The following values in the table are NOT used
2835 ! !!!!!!!!!!!!!! and are just given for reference
2836 DATA WLTSMC/0.023, 0.028, 0.047, 0.084, 0.084, 0.066, &
2837 & 0.069, 0.120, 0.103, 0.100, 0.126, 0.135, &
2838 & 0.069, 0.000, 0.012, 0.028, 0.135, 0.012, &
2839 & 0.023, 0.000, 0.000, 0.000, 0.000, 0.000, &
2840 & 0.000, 0.000, 0.000, 0.000, 0.000, 0.000/
2841 ! !!!!!!!!!!!!!! The following values in the table are NOT used
2842 ! !!!!!!!!!!!!!! and are just given for reference
2843 DATA DRYSMC/0.023, 0.028, 0.047, 0.084, 0.084, 0.066, &
2844 & 0.069, 0.120, 0.103, 0.100, 0.126, 0.135, &
2845 & 0.069, 0.000, 0.012, 0.028, 0.135, 0.012, &
2846 & 0.023, 0.000, 0.000, 0.000, 0.000, 0.000, &
2847 & 0.000, 0.000, 0.000, 0.000, 0.000, 0.000/
2848
2849 ! !!!!!!!!!!!!!! The following values in the table are NOT used
2850 ! !!!!!!!!!!!!!! and are just given for reference
2851 DATA SATDW /0.632E-4, 0.517E-5, 0.807E-5, 0.239E-4, 0.239E-4, &
2852 & 0.143E-4, 0.101E-4, 0.236E-4, 0.113E-4, 0.186E-4, &
2853 & 0.966E-5, 0.115E-4, 0.136E-4, 0.0, 0.998E-5, &
2854 & 0.517E-5, 0.115E-4, 0.998E-5, 0.632E-4, 0.0, &
2855 & 0.0, 0.0, 0.0, 0.0, 0.0, &
2856 & 0.0, 0.0, 0.0, 0.0, 0.0/
2857 ! !!!!!!!!!!!!!! The following values in the table are NOT used
2858 ! !!!!!!!!!!!!!! and are just given for reference
2859 DATA F11 /-1.090, -1.041, -0.568, 0.162, 0.162, -0.327, &
2860 & -1.535, -1.118, -1.297, -3.211, -1.916, -2.258, &
2861 & -0.201, 0.000, -2.287, -1.041, -2.258, -2.287, &
2862 & -1.090, 0.000, 0.000, 0.000, 0.000, 0.000, &
2863 & 0.000, 0.000, 0.000, 0.000, 0.000, 0.000/
2864
2865 ! ----------------------------------------------------------------------
2866 ! SET-UP VEGETATION PARAMETERS FOR A GIVEN VEGETAION TYPE:
2867 ! INPUT: VEGTYP = VEGETATION TYPE (INTEGER INDEX)
2868 ! OUPUT: VEGETATION PARAMETERS
2869 ! SHDFAC: VEGETATION GREENNESS FRACTION
2870 ! RSMIN: MIMIMUM STOMATAL RESISTANCE
2871 ! RGL: PARAMETER USED IN SOLAR RAD TERM OF
2872 ! CANOPY RESISTANCE FUNCTION
2873 ! HS: PARAMETER USED IN VAPOR PRESSURE DEFICIT TERM OF
2874 ! CANOPY RESISTANCE FUNCTION
2875 ! SNUP: THRESHOLD SNOW DEPTH (IN WATER EQUIVALENT M) THAT
2876 ! IMPLIES 100% SNOW COVER
2877 ! ----------------------------------------------------------------------
2878 ! CLASS USGS-WRF VEGETATION/SURFACE TYPE
2879 ! 1 Urban and Built-Up Land
2880 ! 2 Dryland Cropland and Pasture
2881 ! 3 Irrigated Cropland and Pasture
2882 ! 4 Mixed Dryland/Irrigated Cropland and Pasture
2883 ! 5 Cropland/Grassland Mosaic
2884 ! 6 Cropland/Woodland Mosaic
2885 ! 7 Grassland
2886 ! 8 Shrubland
2887 ! 9 Mixed Shrubland/Grassland
2888 ! 10 Savanna
2889 ! 11 Deciduous Broadleaf Forest
2890 ! 12 Deciduous Needleleaf Forest
2891 ! 13 Evergreen Broadleaf Forest
2892 ! 14 Evergreen Needleleaf Forest
2893 ! 15 Mixed Forest
2894 ! 16 Water Bodies
2895 ! 17 Herbaceous Wetland
2896 ! 18 Wooded Wetland
2897 ! 19 Barren or Sparsely Vegetated
2898 ! 20 Herbaceous Tundra
2899 ! 21 Wooded Tundra
2900 ! 22 Mixed Tundra
2901 ! 23 Bare Ground Tundra
2902 ! 24 Snow or Ice
2903 ! 25 Playa
2904 ! 26 Lava
2905 ! 27 White Sand
2906 ! ----------------------------------------------------------------------
2907
2908 INTEGER NROOT
2909 INTEGER NROOT_DATA(MAX_VEGTYP)
2910
2911 REAL FRZFACT
2912 REAL HS
2913 REAL HSTBL(MAX_VEGTYP)
2914 REAL LAI
2915 REAL LAI_DATA(MAX_VEGTYP)
2916 REAL PSISAT
2917 REAL RSMIN
2918 REAL RGL
2919 REAL RGLTBL(MAX_VEGTYP)
2920 REAL RSMTBL(MAX_VEGTYP)
2921 REAL SHDFAC
2922 REAL SNUP
2923 REAL SNUPX(MAX_VEGTYP)
2924 REAL Z0
2925 REAL Z0_DATA(MAX_VEGTYP)
2926
2927 ! ----------------------------------------------------------------------
2928 ! VEGETATION CLASS-RELATED ARRAYS
2929 ! ----------------------------------------------------------------------
2930 ! DATA NROOT_DATA /2,3,3,3,3,3,3,3,3,3,
2931 ! & 4,4,4,4,4,2,2,2,2,3,
2932 ! & 3,3,2,2,2,2,2,0,0,0/
2933 DATA NROOT_DATA /1,3,3,3,3,3,3,3,3,3, &
2934 & 4,4,4,4,4,0,2,2,1,3, &
2935 & 3,3,2,1,1,1,1,0,0,0/
2936 DATA RSMTBL /200.0, 70.0, 70.0, 70.0, 70.0, 70.0, &
2937 & 70.0, 300.0, 170.0, 70.0, 100.0, 150.0, &
2938 & 150.0, 125.0, 125.0, 100.0, 40.0, 100.0, &
2939 & 300.0, 150.0, 150.0, 150.0, 200.0, 200.0, &
2940 & 40.0, 100.0, 300.0, 0.0, 0.0, 0.0/
2941 DATA RGLTBL /100.0, 100.0, 100.0, 100.0, 100.0, 65.0, &
2942 & 100.0, 100.0, 100.0, 65.0, 30.0, 30.0, &
2943 & 30.0, 30.0, 30.0, 30.0, 100.0, 30.0, &
2944 & 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, &
2945 & 100.0, 100.0, 100.0, 0.0, 0.0, 0.0/
2946 DATA HSTBL /42.00, 36.25, 36.25, 36.25, 36.25, 44.14, &
2947 & 36.35, 42.00, 39.18, 54.53, 54.53, 47.35, &
2948 & 41.69, 47.35, 51.93, 51.75, 60.00, 51.93, &
2949 & 42.00, 42.00, 42.00, 42.00, 42.00, 42.00, &
2950 & 36.25, 42.00, 42.00, 0.00, 0.00, 0.00/
2951 DATA SNUPX /0.020, 0.020, 0.020, 0.020, 0.020, 0.020, &
2952 & 0.020, 0.020, 0.020, 0.040, 0.040, 0.040, &
2953 & 0.040, 0.040, 0.040, 0.010, 0.013, 0.020, &
2954 & 0.013, 0.020, 0.020, 0.020, 0.020, 0.013, &
2955 & 0.013, 0.013, 0.013, 0.000, 0.000, 0.000/
2956 DATA Z0_DATA / 1.00, 0.07, 0.07, 0.07, 0.07, 0.15, &
2957 & 0.08, 0.03, 0.05, 0.86, 0.80, 0.85, &
2958 & 2.65, 1.09, 0.80, 0.001, 0.04, 0.05, &
2959 & 0.01, 0.04, 0.06, 0.05, 0.03, 0.001, &
2960 & 0.01, 0.15, 0.01, 0.00, 0.00, 0.00/
2961 DATA LAI_DATA /4.0, 4.0, 4.0, 4.0, 4.0, 4.0, &
2962 & 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, &
2963 & 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, &
2964 & 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, &
2965 & 4.0, 4.0, 4.0, 0.0, 0.0, 0.0/
2966
2967 ! ----------------------------------------------------------------------
2968 ! CLASS PARAMETER 'SLOPETYP' WAS INCLUDED TO ESTIMATE LINEAR RESERVOIR
2969 ! COEFFICIENT 'SLOPE' TO THE BASEFLOW RUNOFF OUT OF THE BOTTOM LAYER.
2970 ! LOWEST CLASS (SLOPETYP=0) MEANS HIGHEST SLOPE PARAMETER = 1.
2971 ! DEFINITION OF SLOPETYP FROM 'ZOBLER' SLOPE TYPE:
2972 ! SLOPE CLASS PERCENT SLOPE
2973 ! 1 0-8
2974 ! 2 8-30
2975 ! 3 > 30
2976 ! 4 0-30
2977 ! 5 0-8 & > 30
2978 ! 6 8-30 & > 30
2979 ! 7 0-8, 8-30, > 30
2980 ! 9 GLACIAL ICE
2981 ! BLANK OCEAN/SEA
2982 ! ----------------------------------------------------------------------
2983 ! NOTE:
2984 ! CLASS 9 FROM 'ZOBLER' FILE SHOULD BE REPLACED BY 8 AND 'BLANK' 9
2985 ! ----------------------------------------------------------------------
2986 REAL SLOPE
2987 REAL SLOPE_DATA(MAX_SLOPETYP)
2988
2989 DATA SLOPE_DATA /0.1, 0.6, 1.0, 0.35, 0.55, 0.8, &
2990 & 0.63, 0.0, 0.0, 0.0, 0.0, 0.0, &
2991 & 0.0 , 0.0, 0.0, 0.0, 0.0, 0.0, &
2992 & 0.0 , 0.0, 0.0, 0.0, 0.0, 0.0, &
2993 & 0.0 , 0.0, 0.0, 0.0, 0.0, 0.0/
2994
2995 ! ----------------------------------------------------------------------
2996 ! SET NAMELIST FILE NAME
2997 ! ----------------------------------------------------------------------
2998 CHARACTER*50 NAMELIST_NAME
2999
3000 ! ----------------------------------------------------------------------
3001 ! SET UNIVERSAL PARAMETERS (NOT DEPENDENT ON SOIL, VEG, SLOPE TYPE)
3002 ! ----------------------------------------------------------------------
3003 INTEGER I
3004 INTEGER NSOIL
3005 INTEGER SLOPETYP
3006 INTEGER SOILTYP
3007 INTEGER VEGTYP
3008
3009 INTEGER BARE
3010 ! DATA BARE /11/
3011 DATA BARE /19/
3012
3013 LOGICAL LPARAM
3014 DATA LPARAM /.TRUE./
3015
3016 LOGICAL LFIRST
3017 DATA LFIRST /.TRUE./
3018
3019 ! ----------------------------------------------------------------------
3020 ! PARAMETER USED TO CALCULATE ROUGHNESS LENGTH OF HEAT.
3021 ! ----------------------------------------------------------------------
3022 REAL CZIL
3023 REAL CZIL_DATA
3024 ! changed in version 2.6 June 2nd 2003
3025 ! DATA CZIL_DATA /0.2/
3026 DATA CZIL_DATA /0.1/
3027
3028 ! ----------------------------------------------------------------------
3029 ! PARAMETER USED TO CALUCULATE VEGETATION EFFECT ON SOIL HEAT FLUX.
3030 ! ----------------------------------------------------------------------
3031 REAL SBETA
3032 REAL SBETA_DATA
3033 DATA SBETA_DATA /-2.0/
3034
3035 ! ----------------------------------------------------------------------
3036 ! BARE SOIL EVAPORATION EXPONENT USED IN DEVAP.
3037 ! ----------------------------------------------------------------------
3038 REAL FXEXP
3039 REAL FXEXP_DATA
3040 DATA FXEXP_DATA /2.0/
3041
3042 ! ----------------------------------------------------------------------
3043 ! SOIL HEAT CAPACITY [J M-3 K-1]
3044 ! ----------------------------------------------------------------------
3045 REAL CSOIL
3046 REAL CSOIL_DATA
3047 ! DATA CSOIL_DATA /1.26E+6/
3048 DATA CSOIL_DATA /2.00E+6/
3049
3050 ! ----------------------------------------------------------------------
3051 ! SPECIFY SNOW DISTRIBUTION SHAPE PARAMETER SALP - SHAPE PARAMETER OF
3052 ! DISTRIBUTION FUNCTION OF SNOW COVER. FROM ANDERSON'S DATA (HYDRO-17)
3053 ! BEST FIT IS WHEN SALP = 2.6
3054 ! ----------------------------------------------------------------------
3055 REAL SALP
3056 REAL SALP_DATA
3057 ! changed for version 2.6 June 2nd 2003
3058 ! DATA SALP_DATA /2.6/
3059 DATA SALP_DATA /4.0/
3060
3061 ! ----------------------------------------------------------------------
3062 ! KDT IS DEFINED BY REFERENCE REFKDT AND DKSAT; REFDK=2.E-6 IS THE SAT.
3063 ! DK. VALUE FOR THE SOIL TYPE 2
3064 ! ----------------------------------------------------------------------
3065 REAL REFDK
3066 REAL REFDK_DATA
3067 DATA REFDK_DATA /2.0E-6/
3068
3069 REAL REFKDT
3070 REAL REFKDT_DATA
3071 DATA REFKDT_DATA /3.0/
3072
3073 REAL FRZX
3074 REAL KDT
3075
3076 ! ----------------------------------------------------------------------
3077 ! FROZEN GROUND PARAMETER, FRZK, DEFINITION: ICE CONTENT THRESHOLD ABOVE
3078 ! WHICH FROZEN SOIL IS IMPERMEABLE REFERENCE VALUE OF THIS PARAMETER FOR
3079 ! THE LIGHT CLAY SOIL (TYPE=3) FRZK = 0.15 M.
3080 ! ----------------------------------------------------------------------
3081 REAL FRZK
3082 REAL FRZK_DATA
3083 DATA FRZK_DATA /0.15/
3084
3085 REAL RTDIS(NSOIL)
3086 REAL SLDPTH(NSOIL)
3087 REAL ZSOIL(NSOIL)
3088
3089 ! ----------------------------------------------------------------------
3090 ! SET TWO CANOPY WATER PARAMETERS.
3091 ! ----------------------------------------------------------------------
3092 REAL CFACTR
3093 REAL CFACTR_DATA
3094 DATA CFACTR_DATA /0.5/
3095
3096 REAL CMCMAX
3097 REAL CMCMAX_DATA
3098 DATA CMCMAX_DATA /0.5E-3/
3099
3100 ! ----------------------------------------------------------------------
3101 ! SET MAX. STOMATAL RESISTANCE.
3102 ! ----------------------------------------------------------------------
3103 REAL RSMAX
3104 REAL RSMAX_DATA
3105 DATA RSMAX_DATA /5000.0/
3106
3107 ! ----------------------------------------------------------------------
3108 ! SET OPTIMUM TRANSPIRATION AIR TEMPERATURE.
3109 ! ----------------------------------------------------------------------
3110 REAL TOPT
3111 REAL TOPT_DATA
3112 DATA TOPT_DATA /298.0/
3113
3114 ! ----------------------------------------------------------------------
3115 ! SPECIFY DEPTH[M] OF LOWER BOUNDARY SOIL TEMPERATURE.
3116 ! ----------------------------------------------------------------------
3117 REAL ZBOT
3118 REAL ZBOT_DATA
3119 ! changed for version 2.5.2
3120 ! DATA ZBOT_DATA /-3.0/
3121 DATA ZBOT_DATA /-8.0/
3122
3123 ! ----------------------------------------------------------------------
3124 ! SET TWO SOIL MOISTURE WILT, SOIL MOISTURE REFERENCE PARAMETERS
3125 ! ----------------------------------------------------------------------
3126 REAL SMLOW
3127 REAL SMLOW_DATA
3128 DATA SMLOW_DATA /0.5/
3129
3130 REAL SMHIGH
3131 REAL SMHIGH_DATA
3132 ! changed in 2.6 from 3 to 6 on June 2nd 2003
3133 DATA SMHIGH_DATA /3.0/
3134 ! DATA SMHIGH_DATA /6.0/
3135
3136 ! ----------------------------------------------------------------------
3137 ! NAMELIST DEFINITION:
3138 ! ----------------------------------------------------------------------
3139 NAMELIST /SOIL_VEG/ SLOPE_DATA, RSMTBL, RGLTBL, HSTBL, SNUPX, &
3140 & BB, DRYSMC, F11, MAXSMC, REFSMC, SATPSI, SATDK, SATDW, &
3141 & WLTSMC, QTZ, LPARAM, ZBOT_DATA, SALP_DATA, CFACTR_DATA, &
3142 & CMCMAX_DATA, SBETA_DATA, RSMAX_DATA, TOPT_DATA, &
3143 & REFDK_DATA, FRZK_DATA, BARE, DEFINED_VEG, DEFINED_SOIL, &
3144 & DEFINED_SLOPE, FXEXP_DATA, NROOT_DATA, REFKDT_DATA, Z0_DATA, &
3145 & CZIL_DATA, LAI_DATA, CSOIL_DATA
3146
3147 ! ----------------------------------------------------------------------
3148 ! READ NAMELIST FILE TO OVERRIDE DEFAULT PARAMETERS ONLY ONCE.
3149 ! NAMELIST_NAME must be 50 characters or less.
3150 ! ----------------------------------------------------------------------
3151 IF (LFIRST) THEN
3152 ! WRITE(*,*) 'READ NAMELIST'
3153 ! OPEN(58, FILE = 'namelist_filename.txt')
3154 ! READ(58,'(A)') NAMELIST_NAME
3155 ! CLOSE(58)
3156 ! WRITE(*,*) 'Namelist Filename is ', NAMELIST_NAME
3157 ! OPEN(59, FILE = NAMELIST_NAME)
3158 ! 50 CONTINUE
3159 ! READ(59, SOIL_VEG, END=100)
3160 ! IF (LPARAM) GOTO 50
3161 ! 100 CONTINUE
3162 ! CLOSE(59)
3163 ! WRITE(*,NML=SOIL_VEG)
3164 LFIRST = .FALSE.
3165 IF (DEFINED_SOIL .GT. MAX_SOILTYP) THEN
3166 WRITE(*,*) 'Warning: DEFINED_SOIL too large in namelist'
3167 STOP 222
3168 ENDIF
3169 IF (DEFINED_VEG .GT. MAX_VEGTYP) THEN
3170 WRITE(*,*) 'Warning: DEFINED_VEG too large in namelist'
3171 STOP 222
3172 ENDIF
3173 IF (DEFINED_SLOPE .GT. MAX_SLOPETYP) THEN
3174 WRITE(*,*) 'Warning: DEFINED_SLOPE too large in namelist'
3175 STOP 222
3176 ENDIF
3177
3178 SMLOW = SMLOW_DATA
3179 SMHIGH = SMHIGH_DATA
3180
3181 DO I = 1,DEFINED_SOIL
3182 SATDW(I) = BB(I)*SATDK(I)*(SATPSI(I)/MAXSMC(I))
3183 F11(I) = ALOG10(SATPSI(I)) + BB(I)*ALOG10(MAXSMC(I)) + 2.0
3184 REFSMC1 = MAXSMC(I)*(5.79E-9/SATDK(I)) &
3185 & **(1.0/(2.0*BB(I)+3.0))
3186 REFSMC(I) = REFSMC1 + (MAXSMC(I)-REFSMC1) / SMHIGH
3187 WLTSMC1 = MAXSMC(I) * (200.0/SATPSI(I))**(-1.0/BB(I))
3188 WLTSMC(I) = WLTSMC1 - SMLOW * WLTSMC1
3189
3190 ! ----------------------------------------------------------------------
3191 ! CURRENT VERSION DRYSMC VALUES THAT EQUATE TO WLTSMC.
3192 ! FUTURE VERSION COULD LET DRYSMC BE INDEPENDENTLY SET VIA NAMELIST.
3193 ! ----------------------------------------------------------------------
3194 DRYSMC(I) = WLTSMC(I)
3195 END DO
3196
3197 ! ----------------------------------------------------------------------
3198 ! END LFIRST BLOCK
3199 ! ----------------------------------------------------------------------
3200 ENDIF
3201
3202 IF (SOILTYP .GT. DEFINED_SOIL) THEN
3203 WRITE(*,*) 'Warning: too many soil types'
3204 STOP 333
3205 ENDIF
3206 IF (VEGTYP .GT. DEFINED_VEG) THEN
3207 WRITE(*,*) 'Warning: too many veg types'
3208 STOP 333
3209 ENDIF
3210 IF (SLOPETYP .GT. DEFINED_SLOPE) THEN
3211 WRITE(*,*) 'Warning: too many slope types'
3212 STOP 333
3213 ENDIF
3214
3215 ! ----------------------------------------------------------------------
3216 ! SET-UP UNIVERSAL PARAMETERS (NOT DEPENDENT ON SOILTYP, VEGTYP OR
3217 ! SLOPETYP)
3218 ! ----------------------------------------------------------------------
3219 ZBOT = ZBOT_DATA
3220 SALP = SALP_DATA
3221 CFACTR = CFACTR_DATA
3222 CMCMAX = CMCMAX_DATA
3223 SBETA = SBETA_DATA
3224 RSMAX = RSMAX_DATA
3225 TOPT = TOPT_DATA
3226 REFDK = REFDK_DATA
3227 FRZK = FRZK_DATA
3228 FXEXP = FXEXP_DATA
3229 REFKDT = REFKDT_DATA
3230 CZIL = CZIL_DATA
3231 CSOIL = CSOIL_DATA
3232
3233 ! ----------------------------------------------------------------------
3234 ! SET-UP SOIL PARAMETERS
3235 ! ----------------------------------------------------------------------
3236 BEXP = BB(SOILTYP)
3237 DKSAT = SATDK(SOILTYP)
3238 DWSAT = SATDW(SOILTYP)
3239 F1 = F11(SOILTYP)
3240 KDT = REFKDT * DKSAT/REFDK
3241 PSISAT = SATPSI(SOILTYP)
3242 QUARTZ = QTZ(SOILTYP)
3243 SMCDRY = DRYSMC(SOILTYP)
3244 SMCMAX = MAXSMC(SOILTYP)
3245 SMCREF = REFSMC(SOILTYP)
3246 SMCWLT = WLTSMC(SOILTYP)
3247 FRZFACT = (SMCMAX / SMCREF) * (0.412 / 0.468)
3248
3249 ! ----------------------------------------------------------------------
3250 ! TO ADJUST FRZK PARAMETER TO ACTUAL SOIL TYPE: FRZK * FRZFACT
3251 ! ----------------------------------------------------------------------
3252 FRZX = FRZK * FRZFACT
3253
3254 ! ----------------------------------------------------------------------
3255 ! SET-UP VEGETATION PARAMETERS
3256 ! ----------------------------------------------------------------------
3257 NROOT = NROOT_DATA(VEGTYP)
3258 SNUP = SNUPX(VEGTYP)
3259 RSMIN = RSMTBL(VEGTYP)
3260 RGL = RGLTBL(VEGTYP)
3261 HS = HSTBL(VEGTYP)
3262 Z0 = Z0_DATA(VEGTYP)
3263 LAI = LAI_DATA(VEGTYP)
3264 IF (VEGTYP .EQ. BARE) SHDFAC = 0.0
3265
3266 IF (NROOT .GT. NSOIL) THEN
3267 WRITE(*,*) 'Warning: too many root layers'
3268 STOP 333
3269 ENDIF
3270
3271 ! ----------------------------------------------------------------------
3272 ! CALCULATE ROOT DISTRIBUTION. PRESENT VERSION ASSUMES UNIFORM
3273 ! DISTRIBUTION BASED ON SOIL LAYER DEPTHS.
3274 ! ----------------------------------------------------------------------
3275 DO I = 1,NROOT
3276 RTDIS(I) = -SLDPTH(I)/ZSOIL(NROOT)
3277 END DO
3278
3279 ! ----------------------------------------------------------------------
3280 ! SET-UP SLOPE PARAMETER
3281 ! ----------------------------------------------------------------------
3282 SLOPE = SLOPE_DATA(SLOPETYP)
3283
3284 ! ----------------------------------------------------------------------
3285 ! END SUBROUTINE REDPRM
3286 ! ----------------------------------------------------------------------
3287 END SUBROUTINE REDPRM
3288
3289 SUBROUTINE ROSR12 (P,A,B,C,D,DELTA,NSOIL)
3290
3291 IMPLICIT NONE
3292
3293 ! ----------------------------------------------------------------------
3294 ! SUBROUTINE ROSR12
3295 ! ----------------------------------------------------------------------
3296 ! INVERT (SOLVE) THE TRI-DIAGONAL MATRIX PROBLEM SHOWN BELOW:
3297 ! ### ### ### ### ### ###
3298 ! #B(1), C(1), 0 , 0 , 0 , . . . , 0 # # # # #
3299 ! #A(2), B(2), C(2), 0 , 0 , . . . , 0 # # # # #
3300 ! # 0 , A(3), B(3), C(3), 0 , . . . , 0 # # # # D(3) #
3301 ! # 0 , 0 , A(4), B(4), C(4), . . . , 0 # # P(4) # # D(4) #
3302 ! # 0 , 0 , 0 , A(5), B(5), . . . , 0 # # P(5) # # D(5) #
3303 ! # . . # # . # = # . #
3304 ! # . . # # . # # . #
3305 ! # . . # # . # # . #
3306 ! # 0 , . . . , 0 , A(M-2), B(M-2), C(M-2), 0 # #P(M-2)# #D(M-2)#
3307 ! # 0 , . . . , 0 , 0 , A(M-1), B(M-1), C(M-1)# #P(M-1)# #D(M-1)#
3308 ! # 0 , . . . , 0 , 0 , 0 , A(M) , B(M) # # P(M) # # D(M) #
3309 ! ### ### ### ### ### ###
3310 ! ----------------------------------------------------------------------
3311 INTEGER K
3312 INTEGER KK
3313 INTEGER NSOIL
3314
3315 REAL A(NSOIL)
3316 REAL B(NSOIL)
3317 REAL C(NSOIL)
3318 REAL D(NSOIL)
3319 REAL DELTA(NSOIL)
3320 REAL P(NSOIL)
3321
3322 ! ----------------------------------------------------------------------
3323 ! INITIALIZE EQN COEF C FOR THE LOWEST SOIL LAYER
3324 ! ----------------------------------------------------------------------
3325 C(NSOIL) = 0.0
3326
3327 ! ----------------------------------------------------------------------
3328 ! SOLVE THE COEFS FOR THE 1ST SOIL LAYER
3329 ! ----------------------------------------------------------------------
3330 P(1) = -C(1) / B(1)
3331 DELTA(1) = D(1) / B(1)
3332
3333 ! ----------------------------------------------------------------------
3334 ! SOLVE THE COEFS FOR SOIL LAYERS 2 THRU NSOIL
3335 ! ----------------------------------------------------------------------
3336 DO K = 2,NSOIL
3337 P(K) = -C(K) * ( 1.0 / (B(K) + A (K) * P(K-1)) )
3338 DELTA(K) = (D(K)-A(K)*DELTA(K-1))*(1.0/(B(K)+A(K)*P(K-1)))
3339 END DO
3340
3341 ! ----------------------------------------------------------------------
3342 ! SET P TO DELTA FOR LOWEST SOIL LAYER
3343 ! ----------------------------------------------------------------------
3344 P(NSOIL) = DELTA(NSOIL)
3345
3346 ! ----------------------------------------------------------------------
3347 ! ADJUST P FOR SOIL LAYERS 2 THRU NSOIL
3348 ! ----------------------------------------------------------------------
3349 DO K = 2,NSOIL
3350 KK = NSOIL - K + 1
3351 P(KK) = P(KK) * P(KK+1) + DELTA(KK)
3352 END DO
3353
3354 ! ----------------------------------------------------------------------
3355 ! END SUBROUTINE ROSR12
3356 ! ----------------------------------------------------------------------
3357 END SUBROUTINE ROSR12
3358
3359 SUBROUTINE SHFLX (SSOIL,STC,SMC,SMCMAX,NSOIL,T1,DT,YY,ZZ1,ZSOIL, &
3360 & TBOT,ZBOT,SMCWLT,PSISAT,SH2O,BEXP,F1,DF1,ICE, &
3361 & QUARTZ,CSOIL)
3362
3363 IMPLICIT NONE
3364
3365 ! ----------------------------------------------------------------------
3366 ! SUBROUTINE SHFLX
3367 ! ----------------------------------------------------------------------
3368 ! UPDATE THE TEMPERATURE STATE OF THE SOIL COLUMN BASED ON THE THERMAL
3369 ! DIFFUSION EQUATION AND UPDATE THE FROZEN SOIL MOISTURE CONTENT BASED
3370 ! ON THE TEMPERATURE.
3371 ! ----------------------------------------------------------------------
3372 INTEGER NSOLD
3373 PARAMETER(NSOLD = 20)
3374
3375 INTEGER I
3376 INTEGER ICE
3377 INTEGER IFRZ
3378 INTEGER NSOIL
3379
3380 REAL AI(NSOLD)
3381 REAL BI(NSOLD)
3382 REAL CI(NSOLD)
3383
3384 REAL BEXP
3385 REAL CSOIL
3386 REAL DF1
3387 REAL DT
3388 REAL F1
3389 REAL PSISAT
3390 REAL QUARTZ
3391 REAL RHSTS(NSOLD)
3392 REAL SSOIL
3393 REAL SH2O(NSOIL)
3394 REAL SMC(NSOIL)
3395 REAL SMCMAX
3396 REAL SMCWLT
3397 REAL STC(NSOIL)
3398 REAL STCF(NSOLD)
3399 REAL T0
3400 REAL T1
3401 REAL TBOT
3402 REAL YY
3403 REAL ZBOT
3404 REAL ZSOIL(NSOIL)
3405 REAL ZZ1
3406
3407 PARAMETER(T0 = 273.15)
3408
3409 ! ----------------------------------------------------------------------
3410 ! HRT ROUTINE CALCS THE RIGHT HAND SIDE OF THE SOIL TEMP DIF EQN
3411 ! ----------------------------------------------------------------------
3412 IF (ICE.EQ.1) THEN
3413
3414 ! ----------------------------------------------------------------------
3415 ! SEA-ICE CASE
3416 ! ----------------------------------------------------------------------
3417 CALL HRTICE (RHSTS,STC,NSOIL,ZSOIL,YY,ZZ1,DF1,AI,BI,CI)
3418
3419 CALL HSTEP (STCF,STC,RHSTS,DT,NSOIL,AI,BI,CI)
3420
3421 ELSE
3422
3423 ! ----------------------------------------------------------------------
3424 ! LAND-MASS CASE
3425 ! ----------------------------------------------------------------------
3426 CALL HRT (RHSTS,STC,SMC,SMCMAX,NSOIL,ZSOIL,YY,ZZ1,TBOT, &
3427 & ZBOT,PSISAT,SH2O,DT, &
3428 & BEXP,F1,DF1,QUARTZ,CSOIL,AI,BI,CI)
3429
3430 CALL HSTEP (STCF,STC,RHSTS,DT,NSOIL,AI,BI,CI)
3431
3432 ENDIF
3433
3434 DO I = 1,NSOIL
3435 STC(I) = STCF(I)
3436 END DO
3437
3438 ! ----------------------------------------------------------------------
3439 ! IN THE NO SNOWPACK CASE (VIA ROUTINE NOPAC BRANCH,) UPDATE THE GRND
3440 ! (SKIN) TEMPERATURE HERE IN RESPONSE TO THE UPDATED SOIL TEMPERATURE
3441 ! PROFILE ABOVE. (NOTE: INSPECTION OF ROUTINE SNOPAC SHOWS THAT T1
3442 ! BELOW IS A DUMMY VARIABLE ONLY, AS SKIN TEMPERATURE IS UPDATED
3443 ! DIFFERENTLY IN ROUTINE SNOPAC)
3444 ! ----------------------------------------------------------------------
3445 T1 = (YY + (ZZ1 - 1.0) * STC(1)) / ZZ1
3446
3447 ! ----------------------------------------------------------------------
3448 ! CALCULATE SURFACE SOIL HEAT FLUX
3449 ! ----------------------------------------------------------------------
3450 SSOIL = DF1 * (STC(1) - T1) / (0.5 * ZSOIL(1))
3451
3452 ! ----------------------------------------------------------------------
3453 ! END SUBROUTINE SHFLX
3454 ! ----------------------------------------------------------------------
3455 END SUBROUTINE SHFLX
3456
3457 SUBROUTINE SMFLX (SMC,NSOIL,CMC,DT,PRCP1,ZSOIL, &
3458 & SH2O,SLOPE,KDT,FRZFACT, &
3459 & SMCMAX,BEXP,SMCWLT,DKSAT,DWSAT, &
3460 & SHDFAC,CMCMAX, &
3461 & RUNOFF1,RUNOFF2,RUNOFF3, &
3462 & EDIR1,EC1,ET1, &
3463 & DRIP)
3464
3465 IMPLICIT NONE
3466
3467 ! ----------------------------------------------------------------------
3468 ! SUBROUTINE SMFLX
3469 ! ----------------------------------------------------------------------
3470 ! CALCULATE SOIL MOISTURE FLUX. THE SOIL MOISTURE CONTENT (SMC - A PER
3471 ! UNIT VOLUME MEASUREMENT) IS A DEPENDENT VARIABLE THAT IS UPDATED WITH
3472 ! PROGNOSTIC EQNS. THE CANOPY MOISTURE CONTENT (CMC) IS ALSO UPDATED.
3473 ! FROZEN GROUND VERSION: NEW STATES ADDED: SH2O, AND FROZEN GROUND
3474 ! CORRECTION FACTOR, FRZFACT AND PARAMETER SLOPE.
3475 ! ----------------------------------------------------------------------
3476 INTEGER NSOLD
3477 PARAMETER(NSOLD = 20)
3478
3479 INTEGER I
3480 INTEGER K
3481 INTEGER NSOIL
3482
3483 REAL AI(NSOLD)
3484 REAL BI(NSOLD)
3485 REAL CI(NSOLD)
3486
3487 REAL BEXP
3488 REAL CMC
3489 REAL CMCMAX
3490 REAL DKSAT
3491 REAL DRIP
3492 REAL DT
3493 REAL DUMMY
3494 REAL DWSAT
3495 REAL EC1
3496 REAL EDIR1
3497 REAL ET1(NSOIL)
3498 REAL EXCESS
3499 REAL FRZFACT
3500 REAL KDT
3501 REAL PCPDRP
3502 REAL PRCP1
3503 REAL RHSCT
3504 REAL RHSTT(NSOLD)
3505 REAL RUNOFF1
3506 REAL RUNOFF2
3507 REAL RUNOFF3
3508 REAL SHDFAC
3509 REAL SMC(NSOIL)
3510 REAL SH2O(NSOIL)
3511 REAL SICE(NSOLD)
3512 REAL SH2OA(NSOLD)
3513 REAL SH2OFG(NSOLD)
3514 REAL SLOPE
3515 REAL SMCMAX
3516 REAL SMCWLT
3517 REAL TRHSCT
3518 REAL ZSOIL(NSOIL)
3519
3520 ! ----------------------------------------------------------------------
3521 ! EXECUTABLE CODE BEGINS HERE.
3522 ! ----------------------------------------------------------------------
3523 DUMMY = 0.
3524
3525 ! ----------------------------------------------------------------------
3526 ! COMPUTE THE RIGHT HAND SIDE OF THE CANOPY EQN TERM ( RHSCT )
3527 ! ----------------------------------------------------------------------
3528 RHSCT = SHDFAC * PRCP1 - EC1
3529
3530 ! ----------------------------------------------------------------------
3531 ! CONVERT RHSCT (A RATE) TO TRHSCT (AN AMOUNT) AND ADD IT TO EXISTING
3532 ! CMC. IF RESULTING AMT EXCEEDS MAX CAPACITY, IT BECOMES DRIP AND WILL
3533 ! FALL TO THE GRND.
3534 ! ----------------------------------------------------------------------
3535 DRIP = 0.
3536 TRHSCT = DT * RHSCT
3537 EXCESS = CMC + TRHSCT
3538 IF (EXCESS .GT. CMCMAX) DRIP = EXCESS - CMCMAX
3539
3540 ! ----------------------------------------------------------------------
3541 ! PCPDRP IS THE COMBINED PRCP1 AND DRIP (FROM CMC) THAT GOES INTO THE
3542 ! SOIL
3543 ! ----------------------------------------------------------------------
3544 PCPDRP = (1. - SHDFAC) * PRCP1 + DRIP / DT
3545
3546 ! ----------------------------------------------------------------------
3547 ! STORE ICE CONTENT AT EACH SOIL LAYER BEFORE CALLING SRT & SSTEP
3548 ! ----------------------------------------------------------------------
3549 DO I = 1,NSOIL
3550 SICE(I) = SMC(I) - SH2O(I)
3551 END DO
3552
3553 ! ----------------------------------------------------------------------
3554 ! CALL SUBROUTINES SRT AND SSTEP TO SOLVE THE SOIL MOISTURE
3555 ! TENDENCY EQUATIONS.
3556 !
3557 ! IF THE INFILTRATING PRECIP RATE IS NONTRIVIAL,
3558 ! (WE CONSIDER NONTRIVIAL TO BE A PRECIP TOTAL OVER THE TIME STEP
3559 ! EXCEEDING ONE ONE-THOUSANDTH OF THE WATER HOLDING CAPACITY OF
3560 ! THE FIRST SOIL LAYER)
3561 ! THEN CALL THE SRT/SSTEP SUBROUTINE PAIR TWICE IN THE MANNER OF
3562 ! TIME SCHEME "F" (IMPLICIT STATE, AVERAGED COEFFICIENT)
3563 ! OF SECTION 2 OF KALNAY AND KANAMITSU (1988, MWR, VOL 116,
3564 ! PAGES 1945-1958)TO MINIMIZE 2-DELTA-T OSCILLATIONS IN THE
3565 ! SOIL MOISTURE VALUE OF THE TOP SOIL LAYER THAT CAN ARISE BECAUSE
3566 ! OF THE EXTREME NONLINEAR DEPENDENCE OF THE SOIL HYDRAULIC
3567 ! DIFFUSIVITY COEFFICIENT AND THE HYDRAULIC CONDUCTIVITY ON THE
3568 ! SOIL MOISTURE STATE
3569 ! OTHERWISE CALL THE SRT/SSTEP SUBROUTINE PAIR ONCE IN THE MANNER OF
3570 ! TIME SCHEME "D" (IMPLICIT STATE, EXPLICIT COEFFICIENT)
3571 ! OF SECTION 2 OF KALNAY AND KANAMITSU
3572 ! PCPDRP IS UNITS OF KG/M**2/S OR MM/S, ZSOIL IS NEGATIVE DEPTH IN M
3573 ! ----------------------------------------------------------------------
3574 ! IF ( PCPDRP .GT. 0.0 ) THEN
3575 IF ( (PCPDRP*DT) .GT. (0.001*1000.0*(-ZSOIL(1))*SMCMAX) ) THEN
3576
3577 ! ----------------------------------------------------------------------
3578 ! FROZEN GROUND VERSION:
3579 ! SMC STATES REPLACED BY SH2O STATES IN SRT SUBR. SH2O & SICE STATES
3580 ! INCLUDED IN SSTEP SUBR. FROZEN GROUND CORRECTION FACTOR, FRZFACT
3581 ! ADDED. ALL WATER BALANCE CALCULATIONS USING UNFROZEN WATER
3582 ! ----------------------------------------------------------------------
3583 CALL SRT (RHSTT,EDIR1,ET1,SH2O,SH2O,NSOIL,PCPDRP,ZSOIL, &
3584 & DWSAT,DKSAT,SMCMAX,BEXP,RUNOFF1, &
3585 & RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZFACT,SICE,AI,BI,CI)
3586
3587 CALL SSTEP (SH2OFG,SH2O,DUMMY,RHSTT,RHSCT,DT,NSOIL,SMCMAX, &
3588 & CMCMAX,RUNOFF3,ZSOIL,SMC,SICE,AI,BI,CI)
3589
3590 DO K = 1,NSOIL
3591 SH2OA(K) = (SH2O(K) + SH2OFG(K)) * 0.5
3592 END DO
3593
3594 CALL SRT (RHSTT,EDIR1,ET1,SH2O,SH2OA,NSOIL,PCPDRP,ZSOIL, &
3595 & DWSAT,DKSAT,SMCMAX,BEXP,RUNOFF1, &
3596 & RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZFACT,SICE,AI,BI,CI)
3597
3598 CALL SSTEP (SH2O,SH2O,CMC,RHSTT,RHSCT,DT,NSOIL,SMCMAX, &
3599 & CMCMAX,RUNOFF3,ZSOIL,SMC,SICE,AI,BI,CI)
3600
3601 ELSE
3602
3603 CALL SRT (RHSTT,EDIR1,ET1,SH2O,SH2O,NSOIL,PCPDRP,ZSOIL, &
3604 & DWSAT,DKSAT,SMCMAX,BEXP,RUNOFF1, &
3605 & RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZFACT,SICE,AI,BI,CI)
3606
3607 CALL SSTEP (SH2O,SH2O,CMC,RHSTT,RHSCT,DT,NSOIL,SMCMAX, &
3608 & CMCMAX,RUNOFF3,ZSOIL,SMC,SICE,AI,BI,CI)
3609
3610 ENDIF
3611
3612 ! RUNOF = RUNOFF
3613
3614 ! ----------------------------------------------------------------------
3615 ! END SUBROUTINE SMFLX
3616 ! ----------------------------------------------------------------------
3617 END SUBROUTINE SMFLX
3618
3619 SUBROUTINE SNFRAC (SNEQV,SNUP,SALP,SNOWH,SNCOVR)
3620
3621 IMPLICIT NONE
3622
3623 ! ----------------------------------------------------------------------
3624 ! SUBROUTINE SNFRAC
3625 ! ----------------------------------------------------------------------
3626 ! CALCULATE SNOW FRACTION (0 -> 1)
3627 ! SNEQV SNOW WATER EQUIVALENT (M)
3628 ! SNUP THRESHOLD SNEQV DEPTH ABOVE WHICH SNCOVR=1
3629 ! SALP TUNING PARAMETER
3630 ! SNCOVR FRACTIONAL SNOW COVER
3631 ! ----------------------------------------------------------------------
3632 REAL SNEQV, SNUP, SALP, SNCOVR, RSNOW, Z0N, SNOWH
3633
3634 ! ----------------------------------------------------------------------
3635 ! SNUP IS VEG-CLASS DEPENDENT SNOWDEPTH THRESHHOLD (SET IN ROUTINE
3636 ! REDPRM) ABOVE WHICH SNOCVR=1.
3637 ! ----------------------------------------------------------------------
3638 IF (SNEQV .LT. SNUP) THEN
3639 RSNOW = SNEQV/SNUP
3640 SNCOVR = 1. - ( EXP(-SALP*RSNOW) - RSNOW*EXP(-SALP))
3641 ELSE
3642 SNCOVR = 1.0
3643 ENDIF
3644
3645 Z0N=0.035
3646 ! FORMULATION OF DICKINSON ET AL. 1986
3647
3648 ! SNCOVR=SNOWH/(SNOWH + 5*Z0N)
3649
3650 ! FORMULATION OF MARSHALL ET AL. 1994
3651 ! SNCOVR=SNEQV/(SNEQV + 2*Z0N)
3652
3653 ! ----------------------------------------------------------------------
3654 ! END SUBROUTINE SNFRAC
3655 ! ----------------------------------------------------------------------
3656 END SUBROUTINE SNFRAC
3657
3658 SUBROUTINE SNOPAC (ETP,ETA,PRCP,PRCP1,SNOWNG,SMC,SMCMAX,SMCWLT, &
3659 & SMCREF,SMCDRY,CMC,CMCMAX,NSOIL,DT, &
3660 & SBETA,DF1, &
3661 & Q2,T1,SFCTMP,T24,TH2,FDOWN,F1,SSOIL,STC,EPSCA, &
3662 & SFCPRS,BEXP,PC,RCH,RR,CFACTR,SNCOVR,ESD,SNDENS,&
3663 & SNOWH,SH2O,SLOPE,KDT,FRZFACT,PSISAT,SNUP, &
3664 & ZSOIL,DWSAT,DKSAT,TBOT,ZBOT,SHDFAC,RUNOFF1, &
3665 & RUNOFF2,RUNOFF3,EDIR,EC,ET,ETT,NROOT,SNOMLT, &
3666 & ICE,RTDIS,QUARTZ,FXEXP,CSOIL, &
3667 & BETA,DRIP,DEW,FLX1,FLX2,FLX3,ESNOW)
3668
3669 IMPLICIT NONE
3670
3671 ! ----------------------------------------------------------------------
3672 ! SUBROUTINE SNOPAC
3673 ! ----------------------------------------------------------------------
3674 ! CALCULATE SOIL MOISTURE AND HEAT FLUX VALUES & UPDATE SOIL MOISTURE
3675 ! CONTENT AND SOIL HEAT CONTENT VALUES FOR THE CASE WHEN A SNOW PACK IS
3676 ! PRESENT.
3677 ! ----------------------------------------------------------------------
3678 INTEGER ICE
3679 INTEGER NROOT
3680 INTEGER NSOIL
3681
3682 LOGICAL SNOWNG
3683
3684 REAL BEXP
3685 REAL BETA
3686 REAL CFACTR
3687 REAL CMC
3688 REAL CMCMAX
3689 REAL CP
3690 REAL CPH2O
3691 REAL CPICE
3692 REAL CSOIL
3693 REAL DENOM
3694 REAL DEW
3695 REAL DF1
3696 REAL DKSAT
3697 REAL DRIP
3698 REAL DSOIL
3699 REAL DTOT
3700 REAL DT
3701 REAL DWSAT
3702 REAL EC
3703 REAL EDIR
3704 REAL EPSCA
3705 REAL ESD
3706 REAL ESDMIN
3707 REAL EXPSNO
3708 REAL EXPSOI
3709 REAL ETA
3710 REAL ETA1
3711 REAL ETP
3712 REAL ETP1
3713 REAL ETP2
3714 REAL ET(NSOIL)
3715 REAL ETT
3716 REAL EX
3717 REAL EXPFAC
3718 REAL FDOWN
3719 REAL FXEXP
3720 REAL FLX1
3721 REAL FLX2
3722 REAL FLX3
3723 REAL F1
3724 REAL KDT
3725 REAL LSUBF
3726 REAL LSUBC
3727 REAL LSUBS
3728 REAL PC
3729 REAL PRCP
3730 REAL PRCP1
3731 REAL Q2
3732 REAL RCH
3733 REAL RR
3734 REAL RTDIS(NSOIL)
3735 REAL SSOIL
3736 REAL SBETA
3737 REAL SSOIL1
3738 REAL SFCTMP
3739 REAL SHDFAC
3740 REAL SIGMA
3741 REAL SMC(NSOIL)
3742 REAL SH2O(NSOIL)
3743 REAL SMCDRY
3744 REAL SMCMAX
3745 REAL SMCREF
3746 REAL SMCWLT
3747 REAL SNOMLT
3748 REAL SNOWH
3749 REAL STC(NSOIL)
3750 REAL T1
3751 REAL T11
3752 REAL T12
3753 REAL T12A
3754 REAL T12B
3755 REAL T24
3756 REAL TBOT
3757 REAL ZBOT
3758 REAL TH2
3759 REAL YY
3760 REAL ZSOIL(NSOIL)
3761 REAL ZZ1
3762 REAL TFREEZ
3763 REAL SALP
3764 REAL SFCPRS
3765 REAL SLOPE
3766 REAL FRZFACT
3767 REAL PSISAT
3768 REAL SNUP
3769 REAL RUNOFF1
3770 REAL RUNOFF2
3771 REAL RUNOFF3
3772 REAL QUARTZ
3773 REAL SNDENS
3774 REAL SNCOND
3775 REAL RSNOW
3776 REAL SNCOVR
3777 REAL QSAT
3778 REAL ETP3
3779 REAL SEH
3780 REAL T14
3781 ! REAL CSNOW
3782
3783 REAL EC1
3784 REAL EDIR1
3785 REAL ET1(NSOIL)
3786 REAL ETT1
3787
3788 REAL ETNS
3789 REAL ETNS1
3790 REAL ESNOW
3791 REAL ESNOW1
3792 REAL ESNOW2
3793 REAL ETANRG
3794
3795 INTEGER K
3796
3797 REAL SNOEXP
3798
3799 PARAMETER(CP = 1004.5)
3800 PARAMETER(CPH2O = 4.218E+3)
3801 PARAMETER(CPICE = 2.106E+3)
3802 PARAMETER(ESDMIN = 1.E-6)
3803 PARAMETER(LSUBF = 3.335E+5)
3804 PARAMETER(LSUBC = 2.501000E+6)
3805 PARAMETER(LSUBS = 2.83E+6)
3806 PARAMETER(SIGMA = 5.67E-8)
3807 PARAMETER(TFREEZ = 273.15)
3808
3809 ! DATA SNOEXP /1.0/
3810 DATA SNOEXP /2.0/
3811
3812 ! ----------------------------------------------------------------------
3813 ! EXECUTABLE CODE BEGINS HERE:
3814 ! CONVERT POTENTIAL EVAP (ETP) FROM KG M-2 S-1 TO M S-1 AND THEN TO AN
3815 ! AMOUNT (M) GIVEN TIMESTEP (DT) AND CALL IT AN EFFECTIVE SNOWPACK
3816 ! REDUCTION AMOUNT, ETP2 (M). THIS IS THE AMOUNT THE SNOWPACK WOULD BE
3817 ! REDUCED DUE TO EVAPORATION FROM THE SNOW SFC DURING THE TIMESTEP.
3818 ! EVAPORATION WILL PROCEED AT THE POTENTIAL RATE UNLESS THE SNOW DEPTH
3819 ! IS LESS THAN THE EXPECTED SNOWPACK REDUCTION.
3820 ! IF SEAICE (ICE=1), BETA REMAINS=1.
3821 ! ----------------------------------------------------------------------
3822 PRCP1 = PRCP1*0.001
3823
3824 ! ETP2 = ETP * 0.001 * DT
3825 ! BETA = 1.0
3826 ! IF (ICE .NE. 1) THEN
3827 ! IF (ESD .LT. ETP2) THEN
3828 ! BETA = ESD / ETP2
3829 ! ENDIF
3830 ! ENDIF
3831
3832 ! ----------------------------------------------------------------------
3833 EDIR = 0.0
3834 EDIR1 = 0.0
3835 EC = 0.0
3836 EC1 = 0.0
3837 DO K = 1,NSOIL
3838 ET(K) = 0.0
3839 ET1(K) = 0.0
3840 ENDDO
3841 ETT = 0.0
3842 ETT1 = 0.0
3843 ETNS = 0.0
3844 ETNS1 = 0.0
3845 ESNOW = 0.0
3846 ESNOW1 = 0.0
3847 ESNOW2 = 0.0
3848 ! ----------------------------------------------------------------------
3849
3850 ! ----------------------------------------------------------------------
3851 ! IF ETP<0 (DOWNWARD) THEN DEWFALL (=FROSTFALL IN THIS CASE).
3852 ! ----------------------------------------------------------------------
3853 DEW = 0.0
3854 ETP1 = ETP*0.001
3855 IF (ETP .LT. 0.0) THEN
3856 ! DEW = -ETP * 0.001
3857 DEW = -ETP1
3858 ! ESNOW2 = ETP * 0.001 * DT
3859 ESNOW2 = ETP1 * DT
3860 ETANRG = ETP*((1.-SNCOVR)*LSUBC + SNCOVR*LSUBS)
3861 ! ENDIF
3862 ELSE
3863 ! ----------------------------------------------------------------------
3864 ! ETP1 = 0.0
3865 ! ETP1 = ETP*0.001
3866 IF (ICE .NE. 1) THEN
3867 IF (SNCOVR .LT. 1.) THEN
3868 ! CALL EVAPO (ETA1,SMC,NSOIL,CMC,ETP1,DT,ZSOIL,
3869 CALL EVAPO (ETNS1,SMC,NSOIL,CMC,ETP1,DT,ZSOIL, &
3870 & SH2O,SMCMAX,BEXP,PC,SMCWLT,DKSAT,DWSAT, &
3871 & SMCREF,SHDFAC,CMCMAX, &
3872 & SMCDRY,CFACTR, &
3873 & EDIR1,EC1,ET1,ETT1,SFCTMP,Q2,NROOT,RTDIS,FXEXP)
3874 ! ENDIF
3875 ! ----------------------------------------------------------------------
3876 EDIR1 = EDIR1*(1.-SNCOVR)
3877 EC1 = EC1*(1.-SNCOVR)
3878 DO K = 1,NSOIL
3879 ET1(K) = ET1(K)*(1.-SNCOVR)
3880 END DO
3881 ETT1 = ETT1*(1.-SNCOVR)
3882 ETNS1 = ETNS1*(1.-SNCOVR)
3883 ! ----------------------------------------------------------------------
3884 EDIR = EDIR1 * 1000.0
3885 EC = EC1 * 1000.0
3886 DO K = 1,NSOIL
3887 ET(K) = ET1(K) * 1000.0
3888 END DO
3889 ETT = ETT1 * 1000.0
3890 ETNS = ETNS1 * 1000.0
3891 ! ----------------------------------------------------------------------
3892 ENDIF
3893 ! ESNOW = ETP*SNCOVR
3894 ! ESNOW1 = ETP*0.001
3895 ! ESNOW1 = ESNOW*0.001
3896 ! ESNOW2 = ESNOW1*DT
3897 ! ETANRG = ESNOW*LSUBS + ETNS*LSUBC
3898 ENDIF
3899 ESNOW = ETP*SNCOVR
3900 ESNOW1 = ESNOW*0.001
3901 ESNOW2 = ESNOW1*DT
3902 ETANRG = ESNOW*LSUBS + ETNS*LSUBC
3903 ENDIF
3904
3905 ! ----------------------------------------------------------------------
3906 ! IF PRECIP IS FALLING, CALCULATE HEAT FLUX FROM SNOW SFC TO NEWLY
3907 ! ACCUMULATING PRECIP. NOTE THAT THIS REFLECTS THE FLUX APPROPRIATE FOR
3908 ! THE NOT-YET-UPDATED SKIN TEMPERATURE (T1). ASSUMES TEMPERATURE OF THE
3909 ! SNOWFALL STRIKING THE GOUND IS =SFCTMP (LOWEST MODEL LEVEL AIR TEMP).
3910 ! ----------------------------------------------------------------------
3911 FLX1 = 0.0
3912 IF (SNOWNG) THEN
3913 FLX1 = CPICE * PRCP * (T1 - SFCTMP)
3914 ELSE
3915 IF (PRCP .GT. 0.0) FLX1 = CPH2O * PRCP * (T1 - SFCTMP)
3916 ENDIF
3917
3918 ! ----------------------------------------------------------------------
3919 ! CALCULATE AN 'EFFECTIVE SNOW-GRND SFC TEMP' (T12) BASED ON HEAT FLUXES
3920 ! BETWEEN THE SNOW PACK AND THE SOIL AND ON NET RADIATION.
3921 ! INCLUDE FLX1 (PRECIP-SNOW SFC) AND FLX2 (FREEZING RAIN LATENT HEAT)
3922 ! FLUXES.
3923 ! FLX2 REFLECTS FREEZING RAIN LATENT HEAT FLUX USING T1 CALCULATED IN
3924 ! PENMAN.
3925 ! ----------------------------------------------------------------------
3926 DSOIL = -(0.5 * ZSOIL(1))
3927 DTOT = SNOWH + DSOIL
3928 DENOM = 1.0 + DF1 / (DTOT * RR * RCH)
3929 ! T12A = ( (FDOWN-FLX1-FLX2-SIGMA*T24)/RCH
3930 ! & + TH2 - SFCTMP - BETA*EPSCA ) / RR
3931 ! T12A = ( (FDOWN-FLX1-FLX2-SIGMA*T24)/RCH
3932 ! M.Ek, 24Nov04, add snow emissivity
3933 T12A = ((FDOWN-FLX1-FLX2 &
3934 & -(0.95*SNCOVR+(1.0-SNCOVR))*SIGMA*T24)/RCH &
3935 & + TH2 - SFCTMP - ETANRG/RCH ) / RR
3936 T12B = DF1 * STC(1) / (DTOT * RR * RCH)
3937 T12 = (SFCTMP + T12A + T12B) / DENOM
3938
3939 ! ----------------------------------------------------------------------
3940 ! IF THE 'EFFECTIVE SNOW-GRND SFC TEMP' IS AT OR BELOW FREEZING, NO SNOW
3941 ! MELT WILL OCCUR. SET THE SKIN TEMP TO THIS EFFECTIVE TEMP. REDUCE
3942 ! (BY SUBLIMINATION ) OR INCREASE (BY FROST) THE DEPTH OF THE SNOWPACK,
3943 ! DEPENDING ON SIGN OF ETP.
3944 ! UPDATE SOIL HEAT FLUX (SSOIL) USING NEW SKIN TEMPERATURE (T1)
3945 ! SINCE NO SNOWMELT, SET ACCUMULATED SNOWMELT TO ZERO, SET 'EFFECTIVE'
3946 ! PRECIP FROM SNOWMELT TO ZERO, SET PHASE-CHANGE HEAT FLUX FROM SNOWMELT
3947 ! TO ZERO.
3948 ! ----------------------------------------------------------------------
3949 IF (T12 .LE. TFREEZ) THEN
3950 T1 = T12
3951 SSOIL = DF1 * (T1 - STC(1)) / DTOT
3952 ! ESD = MAX(0.0, ESD-ETP2)
3953 ESD = MAX(0.0, ESD-ESNOW2)
3954 FLX3 = 0.0
3955 EX = 0.0
3956 SNOMLT = 0.0
3957
3958 ELSE
3959 ! ----------------------------------------------------------------------
3960 ! IF THE 'EFFECTIVE SNOW-GRND SFC TEMP' IS ABOVE FREEZING, SNOW MELT
3961 ! WILL OCCUR. CALL THE SNOW MELT RATE,EX AND AMT, SNOMLT. REVISE THE
3962 ! EFFECTIVE SNOW DEPTH. REVISE THE SKIN TEMP BECAUSE IT WOULD HAVE CHGD
3963 ! DUE TO THE LATENT HEAT RELEASED BY THE MELTING. CALC THE LATENT HEAT
3964 ! RELEASED, FLX3. SET THE EFFECTIVE PRECIP, PRCP1 TO THE SNOW MELT RATE,
3965 ! EX FOR USE IN SMFLX. ADJUSTMENT TO T1 TO ACCOUNT FOR SNOW PATCHES.
3966 ! CALCULATE QSAT VALID AT FREEZING POINT. NOTE THAT ESAT (SATURATION
3967 ! VAPOR PRESSURE) VALUE OF 6.11E+2 USED HERE IS THAT VALID AT FRZZING
3968 ! POINT. NOTE THAT ETP FROM CALL PENMAN IN SFLX IS IGNORED HERE IN
3969 ! FAVOR OF BULK ETP OVER 'OPEN WATER' AT FREEZING TEMP.
3970 ! UPDATE SOIL HEAT FLUX (S) USING NEW SKIN TEMPERATURE (T1)
3971 ! ----------------------------------------------------------------------
3972 ! T1 = TFREEZ * SNCOVR + T12 * (1.0 - SNCOVR)
3973 ! mek Feb2004
3974 ! non-linear weighting of snow vs non-snow covered portions of gridbox
3975 ! so with SNOEXP = 2.0 (>1), surface skin temperature is higher than for
3976 ! the linear case (SNOEXP = 1).
3977 T1 = TFREEZ * SNCOVR**SNOEXP + T12 * (1.0 - SNCOVR**SNOEXP)
3978 ! QSAT = (0.622*6.11E2)/(SFCPRS-0.378*6.11E2)
3979 ! ETP = RCH*(QSAT-Q2)/CP
3980 ! ETP2 = ETP*0.001*DT
3981 ! BETA = 1.0
3982 SSOIL = DF1 * (T1 - STC(1)) / DTOT
3983
3984 ! ----------------------------------------------------------------------
3985 ! IF POTENTIAL EVAP (SUBLIMATION) GREATER THAN DEPTH OF SNOWPACK.
3986 ! BETA<1
3987 ! SNOWPACK HAS SUBLIMATED AWAY, SET DEPTH TO ZERO.
3988 ! ----------------------------------------------------------------------
3989 ! IF (ESD .LE. ETP2) THEN
3990 ! IF (ESD .LE. ESNOW2) THEN
3991 IF (ESD-ESNOW2 .LE. ESDMIN) THEN
3992 ! BETA = ESD / ETP2
3993 ESD = 0.0
3994 EX = 0.0
3995 SNOMLT = 0.0
3996
3997 ELSE
3998 ! ----------------------------------------------------------------------
3999 ! POTENTIAL EVAP (SUBLIMATION) LESS THAN DEPTH OF SNOWPACK, RETAIN
4000 ! BETA=1.
4001 ! SNOWPACK (ESD) REDUCED BY POTENTIAL EVAP RATE
4002 ! ETP3 (CONVERT TO FLUX)
4003 ! ----------------------------------------------------------------------
4004 ! ESD = ESD-ETP2
4005 ESD = ESD-ESNOW2
4006 ! ETP3 = ETP*LSUBC
4007 SEH = RCH*(T1-TH2)
4008 T14 = T1*T1
4009 T14 = T14*T14
4010 ! FLX3 = FDOWN - FLX1 - FLX2 - SIGMA*T14 - SSOIL - SEH - ETP3
4011 ! FLX3 = FDOWN - FLX1 - FLX2 - SIGMA*T14 - SSOIL - SEH - ETANRG
4012 ! M.Ek, 24Nov04, add snow emissivity
4013 FLX3 = FDOWN - FLX1 - FLX2 - &
4014 & (0.95*SNCOVR+(1.0-SNCOVR))*SIGMA*T14 - SSOIL - SEH - ETANRG
4015 IF (FLX3 .LE. 0.0) FLX3 = 0.0
4016 EX = FLX3*0.001/LSUBF
4017
4018 ! ----------------------------------------------------------------------
4019 ! SNOWMELT REDUCTION DEPENDING ON SNOW COVER
4020 ! IF SNOW COVER LESS THAN 5% NO SNOWMELT REDUCTION
4021 ! ***NOTE: DOES 'IF' BELOW FAIL TO MATCH THE MELT WATER WITH THE MELT
4022 ! ENERGY?
4023 ! ----------------------------------------------------------------------
4024 ! IF (SNCOVR .GT. 0.05) EX = EX * SNCOVR
4025 SNOMLT = EX * DT
4026
4027 ! ----------------------------------------------------------------------
4028 ! ESDMIN REPRESENTS A SNOWPACK DEPTH THRESHOLD VALUE BELOW WHICH WE
4029 ! CHOOSE NOT TO RETAIN ANY SNOWPACK, AND INSTEAD INCLUDE IT IN SNOWMELT.
4030 ! ----------------------------------------------------------------------
4031 IF (ESD-SNOMLT .GE. ESDMIN) THEN
4032 ESD = ESD - SNOMLT
4033
4034 ELSE
4035 ! ----------------------------------------------------------------------
4036 ! SNOWMELT EXCEEDS SNOW DEPTH
4037 ! ----------------------------------------------------------------------
4038 EX = ESD/DT
4039 FLX3 = EX*1000.0*LSUBF
4040 SNOMLT = ESD
4041 ESD = 0.0
4042
4043 ENDIF
4044 ! ----------------------------------------------------------------------
4045 ! END OF 'ESD .LE. ETP2' IF-BLOCK
4046 ! ----------------------------------------------------------------------
4047 ENDIF
4048
4049 PRCP1 = PRCP1 + EX
4050
4051 ! ----------------------------------------------------------------------
4052 ! END OF 'T12 .LE. TFREEZ' IF-BLOCK
4053 ! ----------------------------------------------------------------------
4054 ENDIF
4055
4056 ! ----------------------------------------------------------------------
4057 ! FINAL BETA NOW IN HAND, SO COMPUTE EVAPORATION. EVAP EQUALS ETP
4058 ! UNLESS BETA<1.
4059 ! ----------------------------------------------------------------------
4060 ! ETA = BETA*ETP
4061
4062 ! ----------------------------------------------------------------------
4063 ! SET THE EFFECTIVE POTNL EVAPOTRANSP (ETP1) TO ZERO SINCE THIS IS SNOW
4064 ! CASE, SO SURFACE EVAP NOT CALCULATED FROM EDIR, EC, OR ETT IN SMFLX
4065 ! (BELOW).
4066 ! IF SEAICE (ICE=1) SKIP CALL TO SMFLX.
4067 ! SMFLX RETURNS UPDATED SOIL MOISTURE VALUES. IN THIS, THE SNOW PACK
4068 ! CASE, ETA1 IS NOT USED IN CALCULATION OF EVAP.
4069 ! ----------------------------------------------------------------------
4070 ! ETP1 = 0.0
4071 IF (ICE .NE. 1) THEN
4072 ! CALL EVAPO (ETA1,SMC,NSOIL,CMC,ETP1,DT,ZSOIL,
4073 ! & SH2O,SMCMAX,BEXP,PC,SMCWLT,DKSAT,DWSAT,
4074 ! & SMCREF,SHDFAC,CMCMAX,
4075 ! & SMCDRY,CFACTR,
4076 ! & EDIR1,EC1,ET1,ETT1,SFCTMP,Q2,NROOT,RTDIS,FXEXP)
4077 CALL SMFLX (SMC,NSOIL,CMC,DT,PRCP1,ZSOIL, &
4078 & SH2O,SLOPE,KDT,FRZFACT, &
4079 & SMCMAX,BEXP,SMCWLT,DKSAT,DWSAT, &
4080 & SHDFAC,CMCMAX, &
4081 & RUNOFF1,RUNOFF2,RUNOFF3, &
4082 & EDIR1,EC1,ET1, &
4083 & DRIP)
4084
4085 ENDIF
4086
4087 ! ----------------------------------------------------------------------
4088 ! EDIR = EDIR1 * 1000.0
4089 ! EC = EC1 * 1000.0
4090 ! ETT = ETT1 * 1000.0
4091 ! ET(1) = ET1(1) * 1000.0
4092 ! ET(2) = ET1(2) * 1000.0
4093 ! ET(3) = ET1(3) * 1000.0
4094 ! ET(4) = ET1(4) * 1000.0
4095 ! ----------------------------------------------------------------------
4096
4097 ! ----------------------------------------------------------------------
4098 ! BEFORE CALL SHFLX IN THIS SNOWPACK CASE, SET ZZ1 AND YY ARGUMENTS TO
4099 ! SPECIAL VALUES THAT ENSURE THAT GROUND HEAT FLUX CALCULATED IN SHFLX
4100 ! MATCHES THAT ALREADY COMPUTER FOR BELOW THE SNOWPACK, THUS THE SFC
4101 ! HEAT FLUX TO BE COMPUTED IN SHFLX WILL EFFECTIVELY BE THE FLUX AT THE
4102 ! SNOW TOP SURFACE. T11 IS A DUMMY ARGUEMENT SO WE WILL NOT USE THE
4103 ! SKIN TEMP VALUE AS REVISED BY SHFLX.
4104 ! ----------------------------------------------------------------------
4105 ZZ1 = 1.0
4106 YY = STC(1)-0.5*SSOIL*ZSOIL(1)*ZZ1/DF1
4107 T11 = T1
4108
4109 ! ----------------------------------------------------------------------
4110 ! SHFLX WILL CALC/UPDATE THE SOIL TEMPS. NOTE: THE SUB-SFC HEAT FLUX
4111 ! (SSOIL1) AND THE SKIN TEMP (T11) OUTPUT FROM THIS SHFLX CALL ARE NOT
4112 ! USED IN ANY SUBSEQUENT CALCULATIONS. RATHER, THEY ARE DUMMY VARIABLES
4113 ! HERE IN THE SNOPAC CASE, SINCE THE SKIN TEMP AND SUB-SFC HEAT FLUX ARE
4114 ! UPDATED INSTEAD NEAR THE BEGINNING OF THE CALL TO SNOPAC.
4115 ! ----------------------------------------------------------------------
4116 CALL SHFLX (SSOIL1,STC,SMC,SMCMAX,NSOIL,T11,DT,YY,ZZ1,ZSOIL, &
4117 & TBOT,ZBOT,SMCWLT,PSISAT,SH2O,BEXP,F1,DF1,ICE, &
4118 & QUARTZ,CSOIL)
4119
4120 ! ----------------------------------------------------------------------
4121 ! SNOW DEPTH AND DENSITY ADJUSTMENT BASED ON SNOW COMPACTION. YY IS
4122 ! ASSUMED TO BE THE SOIL TEMPERTURE AT THE TOP OF THE SOIL COLUMN.
4123 ! ----------------------------------------------------------------------
4124 IF (ESD .GT. 0.) THEN
4125 CALL SNOWPACK (ESD,DT,SNOWH,SNDENS,T1,YY)
4126 ELSE
4127 ESD = 0.
4128 SNOWH = 0.
4129 SNDENS = 0.
4130 SNCOND = 1.
4131 SNCOVR = 0.
4132 ENDIF
4133
4134 ! ----------------------------------------------------------------------
4135 ! END SUBROUTINE SNOPAC
4136 ! ----------------------------------------------------------------------
4137 END SUBROUTINE SNOPAC
4138
4139 SUBROUTINE SNOWPACK (ESD,DTSEC,SNOWH,SNDENS,TSNOW,TSOIL)
4140
4141 IMPLICIT NONE
4142
4143 ! ----------------------------------------------------------------------
4144 ! SUBROUTINE SNOWPACK
4145 ! ----------------------------------------------------------------------
4146 ! CALCULATE COMPACTION OF SNOWPACK UNDER CONDITIONS OF INCREASING SNOW
4147 ! DENSITY, AS OBTAINED FROM AN APPROXIMATE SOLUTION OF E. ANDERSON'S
4148 ! DIFFERENTIAL EQUATION (3.29), NOAA TECHNICAL REPORT NWS 19, BY VICTOR
4149 ! KOREN, 03/25/95.
4150 ! ----------------------------------------------------------------------
4151 ! ESD WATER EQUIVALENT OF SNOW (M)
4152 ! DTSEC TIME STEP (SEC)
4153 ! SNOWH SNOW DEPTH (M)
4154 ! SNDENS SNOW DENSITY (G/CM3=DIMENSIONLESS FRACTION OF H2O DENSITY)
4155 ! TSNOW SNOW SURFACE TEMPERATURE (K)
4156 ! TSOIL SOIL SURFACE TEMPERATURE (K)
4157 !
4158 ! SUBROUTINE WILL RETURN NEW VALUES OF SNOWH AND SNDENS
4159 ! ----------------------------------------------------------------------
4160 INTEGER IPOL, J
4161
4162 REAL BFAC,C1,C2,SNDENS,DSX,DTHR,DTSEC,DW,SNOWHC,SNOWH,PEXP,TAVGC, &
4163 & TSNOW,TSNOWC,TSOIL,TSOILC,ESD,ESDC,ESDCX,G,KN
4164
4165 PARAMETER(C1 = 0.01, C2=21.0, G=9.81, KN=4000.0)
4166
4167 ! ----------------------------------------------------------------------
4168 ! CONVERSION INTO SIMULATION UNITS
4169 ! ----------------------------------------------------------------------
4170 SNOWHC = SNOWH*100.
4171 ESDC = ESD*100.
4172 DTHR = DTSEC/3600.
4173 TSNOWC = TSNOW-273.15
4174 TSOILC = TSOIL-273.15
4175
4176 ! ----------------------------------------------------------------------
4177 ! CALCULATING OF AVERAGE TEMPERATURE OF SNOW PACK
4178 ! ----------------------------------------------------------------------
4179 TAVGC = 0.5*(TSNOWC+TSOILC)
4180
4181 ! ----------------------------------------------------------------------
4182 ! CALCULATING OF SNOW DEPTH AND DENSITY AS A RESULT OF COMPACTION
4183 ! SNDENS=DS0*(EXP(BFAC*ESD)-1.)/(BFAC*ESD)
4184 ! BFAC=DTHR*C1*EXP(0.08*TAVGC-C2*DS0)
4185 ! NOTE: BFAC*ESD IN SNDENS EQN ABOVE HAS TO BE CAREFULLY TREATED
4186 ! NUMERICALLY BELOW:
4187 ! C1 IS THE FRACTIONAL INCREASE IN DENSITY (1/(CM*HR))
4188 ! C2 IS A CONSTANT (CM3/G) KOJIMA ESTIMATED AS 21 CMS/G
4189 ! ----------------------------------------------------------------------
4190 IF (ESDC .GT. 1.E-2) THEN
4191 ESDCX = ESDC
4192 ELSE
4193 ESDCX = 1.E-2
4194 ENDIF
4195 BFAC = DTHR*C1*EXP(0.08*TAVGC-C2*SNDENS)
4196
4197 ! DSX = SNDENS*((DEXP(BFAC*ESDC)-1.)/(BFAC*ESDC))
4198 ! ----------------------------------------------------------------------
4199 ! THE FUNCTION OF THE FORM (e**x-1)/x IMBEDDED IN ABOVE EXPRESSION
4200 ! FOR DSX WAS CAUSING NUMERICAL DIFFICULTIES WHEN THE DENOMINATOR "x"
4201 ! (I.E. BFAC*ESDC) BECAME ZERO OR APPROACHED ZERO (DESPITE THE FACT THAT
4202 ! THE ANALYTICAL FUNCTION (e**x-1)/x HAS A WELL DEFINED LIMIT AS
4203 ! "x" APPROACHES ZERO), HENCE BELOW WE REPLACE THE (e**x-1)/x
4204 ! EXPRESSION WITH AN EQUIVALENT, NUMERICALLY WELL-BEHAVED
4205 ! POLYNOMIAL EXPANSION.
4206 !
4207 ! NUMBER OF TERMS OF POLYNOMIAL EXPANSION, AND HENCE ITS ACCURACY,
4208 ! IS GOVERNED BY ITERATION LIMIT "IPOL".
4209 ! IPOL GREATER THAN 9 ONLY MAKES A DIFFERENCE ON DOUBLE
4210 ! PRECISION (RELATIVE ERRORS GIVEN IN PERCENT %).
4211 ! IPOL=9, FOR REL.ERROR <~ 1.6 E-6 % (8 SIGNIFICANT DIGITS)
4212 ! IPOL=8, FOR REL.ERROR <~ 1.8 E-5 % (7 SIGNIFICANT DIGITS)
4213 ! IPOL=7, FOR REL.ERROR <~ 1.8 E-4 % ...
4214 ! ----------------------------------------------------------------------
4215 IPOL = 4
4216 PEXP = 0.
4217 DO J = IPOL,1,-1
4218 ! PEXP = (1. + PEXP)*BFAC*ESDC/REAL(J+1)
4219 PEXP = (1. + PEXP)*BFAC*ESDCX/REAL(J+1)
4220 END DO
4221 PEXP = PEXP + 1.
4222
4223 DSX = SNDENS*(PEXP)
4224 ! ----------------------------------------------------------------------
4225 ! ABOVE LINE ENDS POLYNOMIAL SUBSTITUTION
4226 ! ----------------------------------------------------------------------
4227 ! END OF KOREAN FORMULATION
4228
4229 ! BASE FORMULATION (COGLEY ET AL., 1990)
4230 ! CONVERT DENSITY FROM G/CM3 TO KG/M3
4231 ! DSM=SNDENS*1000.0
4232
4233 ! DSX=DSM+DTSEC*0.5*DSM*G*ESD/
4234 ! & (1E7*EXP(-0.02*DSM+KN/(TAVGC+273.16)-14.643))
4235
4236 ! CONVERT DENSITY FROM KG/M3 TO G/CM3
4237 ! DSX=DSX/1000.0
4238
4239 ! END OF COGLEY ET AL. FORMULATION
4240
4241 ! ----------------------------------------------------------------------
4242 ! SET UPPER/LOWER LIMIT ON SNOW DENSITY
4243 ! ----------------------------------------------------------------------
4244 IF (DSX .GT. 0.40) DSX = 0.40
4245 IF (DSX .LT. 0.05) DSX = 0.05
4246 SNDENS = DSX
4247 ! ----------------------------------------------------------------------
4248 ! UPDATE OF SNOW DEPTH AND DENSITY DEPENDING ON LIQUID WATER DURING
4249 ! SNOWMELT. ASSUMED THAT 13% OF LIQUID WATER CAN BE STORED IN SNOW PER
4250 ! DAY DURING SNOWMELT TILL SNOW DENSITY 0.40.
4251 ! ----------------------------------------------------------------------
4252 IF (TSNOWC .GE. 0.) THEN
4253 DW = 0.13*DTHR/24.
4254 SNDENS = SNDENS*(1.-DW)+DW
4255 IF (SNDENS .GT. 0.40) SNDENS = 0.40
4256 ENDIF
4257
4258 ! ----------------------------------------------------------------------
4259 ! CALCULATE SNOW DEPTH (CM) FROM SNOW WATER EQUIVALENT AND SNOW DENSITY.
4260 ! CHANGE SNOW DEPTH UNITS TO METERS
4261 ! ----------------------------------------------------------------------
4262 SNOWHC = ESDC/SNDENS
4263 SNOWH = SNOWHC*0.01
4264
4265 ! ----------------------------------------------------------------------
4266 ! END SUBROUTINE SNOWPACK
4267 ! ----------------------------------------------------------------------
4268 END SUBROUTINE SNOWPACK
4269
4270 SUBROUTINE SNOWZ0 (SNCOVR,Z0)
4271
4272 IMPLICIT NONE
4273
4274 ! ----------------------------------------------------------------------
4275 ! SUBROUTINE SNOWZ0
4276 ! ----------------------------------------------------------------------
4277 ! CALCULATE TOTAL ROUGHNESS LENGTH OVER SNOW
4278 ! SNCOVR FRACTIONAL SNOW COVER
4279 ! Z0 ROUGHNESS LENGTH (m)
4280 ! Z0S SNOW ROUGHNESS LENGTH:=0.001 (m)
4281 ! ----------------------------------------------------------------------
4282 REAL SNCOVR, Z0, Z0S
4283 ! PARAMETER (Z0S=0.001)
4284
4285 ! CURRENT NOAH LSM CONDITION - MBEK, 09-OCT-2001
4286 Z0S = Z0
4287 !
4288 Z0 = (1-SNCOVR)*Z0 + SNCOVR*Z0S
4289 ! ----------------------------------------------------------------------
4290 ! END SUBROUTINE SNOWZ0
4291 ! ----------------------------------------------------------------------
4292 END SUBROUTINE SNOWZ0
4293
4294 SUBROUTINE SNOW_NEW (TEMP,NEWSN,SNOWH,SNDENS)
4295
4296 IMPLICIT NONE
4297
4298 ! ----------------------------------------------------------------------
4299 ! SUBROUTINE SNOW_NEW
4300 ! ----------------------------------------------------------------------
4301 ! CALCULATE SNOW DEPTH AND DENSITITY TO ACCOUNT FOR THE NEW SNOWFALL.
4302 ! NEW VALUES OF SNOW DEPTH & DENSITY RETURNED.
4303 !
4304 ! TEMP AIR TEMPERATURE (K)
4305 ! NEWSN NEW SNOWFALL (M)
4306 ! SNOWH SNOW DEPTH (M)
4307 ! SNDENS SNOW DENSITY (G/CM3=DIMENSIONLESS FRACTION OF H2O DENSITY)
4308 ! ----------------------------------------------------------------------
4309 REAL SNDENS
4310 REAL DSNEW
4311 REAL SNOWHC
4312 REAL HNEWC
4313 REAL SNOWH
4314 REAL NEWSN
4315 REAL NEWSNC
4316 REAL TEMP
4317 REAL TEMPC
4318
4319 ! ----------------------------------------------------------------------
4320 ! CONVERSION INTO SIMULATION UNITS
4321 ! ----------------------------------------------------------------------
4322 SNOWHC = SNOWH*100.
4323 NEWSNC = NEWSN*100.
4324 TEMPC = TEMP-273.15
4325
4326 ! ----------------------------------------------------------------------
4327 ! CALCULATING NEW SNOWFALL DENSITY DEPENDING ON TEMPERATURE
4328 ! EQUATION FROM GOTTLIB L. 'A GENERAL RUNOFF MODEL FOR SNOWCOVERED
4329 ! AND GLACIERIZED BASIN', 6TH NORDIC HYDROLOGICAL CONFERENCE,
4330 ! VEMADOLEN, SWEDEN, 1980, 172-177PP.
4331 !-----------------------------------------------------------------------
4332 IF (TEMPC .LE. -15.) THEN
4333 DSNEW = 0.05
4334 ELSE
4335 DSNEW = 0.05+0.0017*(TEMPC+15.)**1.5
4336 ENDIF
4337
4338 ! ----------------------------------------------------------------------
4339 ! ADJUSTMENT OF SNOW DENSITY DEPENDING ON NEW SNOWFALL
4340 ! ----------------------------------------------------------------------
4341 HNEWC = NEWSNC/DSNEW
4342 SNDENS = (SNOWHC*SNDENS+HNEWC*DSNEW)/(SNOWHC+HNEWC)
4343 SNOWHC = SNOWHC+HNEWC
4344 SNOWH = SNOWHC*0.01
4345
4346 ! ----------------------------------------------------------------------
4347 ! END SUBROUTINE SNOW_NEW
4348 ! ----------------------------------------------------------------------
4349 END SUBROUTINE SNOW_NEW
4350
4351 SUBROUTINE SRT (RHSTT,EDIR,ET,SH2O,SH2OA,NSOIL,PCPDRP, &
4352 & ZSOIL,DWSAT,DKSAT,SMCMAX,BEXP,RUNOFF1, &
4353 & RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZX,SICE,AI,BI,CI)
4354
4355 IMPLICIT NONE
4356
4357 ! ----------------------------------------------------------------------
4358 ! SUBROUTINE SRT
4359 ! ----------------------------------------------------------------------
4360 ! CALCULATE THE RIGHT HAND SIDE OF THE TIME TENDENCY TERM OF THE SOIL
4361 ! WATER DIFFUSION EQUATION. ALSO TO COMPUTE ( PREPARE ) THE MATRIX
4362 ! COEFFICIENTS FOR THE TRI-DIAGONAL MATRIX OF THE IMPLICIT TIME SCHEME.
4363 ! ----------------------------------------------------------------------
4364 INTEGER NSOLD
4365 PARAMETER(NSOLD = 20)
4366
4367 INTEGER CVFRZ
4368 INTEGER IALP1
4369 INTEGER IOHINF
4370 INTEGER J
4371 INTEGER JJ
4372 INTEGER K
4373 INTEGER KS
4374 INTEGER NSOIL
4375
4376 REAL ACRT
4377 REAL AI(NSOLD)
4378 REAL BEXP
4379 REAL BI(NSOLD)
4380 REAL CI(NSOLD)
4381 REAL DD
4382 REAL DDT
4383 REAL DDZ
4384 REAL DDZ2
4385 REAL DENOM
4386 REAL DENOM2
4387 REAL DICE
4388 REAL DKSAT
4389 REAL DMAX(NSOLD)
4390 REAL DSMDZ
4391 REAL DSMDZ2
4392 REAL DT
4393 REAL DT1
4394 REAL DWSAT
4395 REAL EDIR
4396 REAL ET(NSOIL)
4397 REAL FCR
4398 REAL FRZX
4399 REAL INFMAX
4400 REAL KDT
4401 REAL MXSMC
4402 REAL MXSMC2
4403 REAL NUMER
4404 REAL PCPDRP
4405 REAL PDDUM
4406 REAL PX
4407 REAL RHSTT(NSOIL)
4408 REAL RUNOFF1
4409 REAL RUNOFF2
4410 REAL SH2O(NSOIL)
4411 REAL SH2OA(NSOIL)
4412 REAL SICE(NSOIL)
4413 REAL SICEMAX
4414 REAL SLOPE
4415 REAL SLOPX
4416 REAL SMCAV
4417 REAL SMCMAX
4418 REAL SMCWLT
4419 REAL SSTT
4420 REAL SUM
4421 REAL VAL
4422 REAL WCND
4423 REAL WCND2
4424 REAL WDF
4425 REAL WDF2
4426 REAL ZSOIL(NSOIL)
4427
4428 ! ----------------------------------------------------------------------
4429 ! FROZEN GROUND VERSION:
4430 ! REFERENCE FROZEN GROUND PARAMETER, CVFRZ, IS A SHAPE PARAMETER OF
4431 ! AREAL DISTRIBUTION FUNCTION OF SOIL ICE CONTENT WHICH EQUALS 1/CV.
4432 ! CV IS A COEFFICIENT OF SPATIAL VARIATION OF SOIL ICE CONTENT. BASED
4433 ! ON FIELD DATA CV DEPENDS ON AREAL MEAN OF FROZEN DEPTH, AND IT CLOSE
4434 ! TO CONSTANT = 0.6 IF AREAL MEAN FROZEN DEPTH IS ABOVE 20 CM. THAT IS
4435 ! WHY PARAMETER CVFRZ = 3 (INT{1/0.6*0.6}).
4436 ! CURRENT LOGIC DOESN'T ALLOW CVFRZ BE BIGGER THAN 3
4437 ! ----------------------------------------------------------------------
4438 PARAMETER(CVFRZ = 3)
4439
4440 ! ----------------------------------------------------------------------
4441 ! DETERMINE RAINFALL INFILTRATION RATE AND RUNOFF. INCLUDE THE
4442 ! INFILTRATION FORMULE FROM SCHAAKE AND KOREN MODEL.
4443 ! MODIFIED BY Q DUAN
4444 ! ----------------------------------------------------------------------
4445 IOHINF=1
4446
4447 ! ----------------------------------------------------------------------
4448 ! LET SICEMAX BE THE GREATEST, IF ANY, FROZEN WATER CONTENT WITHIN SOIL
4449 ! LAYERS.
4450 ! ----------------------------------------------------------------------
4451 SICEMAX = 0.0
4452 DO KS=1,NSOIL
4453 IF (SICE(KS) .GT. SICEMAX) SICEMAX = SICE(KS)
4454 END DO
4455
4456 ! ----------------------------------------------------------------------
4457 ! DETERMINE RAINFALL INFILTRATION RATE AND RUNOFF
4458 ! ----------------------------------------------------------------------
4459 PDDUM = PCPDRP
4460 RUNOFF1 = 0.0
4461 IF (PCPDRP .NE. 0.0) THEN
4462
4463 ! ----------------------------------------------------------------------
4464 ! MODIFIED BY Q. DUAN, 5/16/94
4465 ! ----------------------------------------------------------------------
4466 ! IF (IOHINF .EQ. 1) THEN
4467
4468 DT1 = DT/86400.
4469 SMCAV = SMCMAX - SMCWLT
4470 DMAX(1)=-ZSOIL(1)*SMCAV
4471
4472 ! ----------------------------------------------------------------------
4473 ! FROZEN GROUND VERSION:
4474 ! ----------------------------------------------------------------------
4475 DICE = -ZSOIL(1) * SICE(1)
4476
4477 DMAX(1)=DMAX(1)*(1.0 - (SH2OA(1)+SICE(1)-SMCWLT)/SMCAV)
4478 DD=DMAX(1)
4479
4480 DO KS=2,NSOIL
4481
4482 ! ----------------------------------------------------------------------
4483 ! FROZEN GROUND VERSION:
4484 ! ----------------------------------------------------------------------
4485 DICE = DICE + ( ZSOIL(KS-1) - ZSOIL(KS) ) * SICE(KS)
4486
4487 DMAX(KS) = (ZSOIL(KS-1)-ZSOIL(KS))*SMCAV
4488 DMAX(KS) = DMAX(KS)*(1.0 - (SH2OA(KS)+SICE(KS)-SMCWLT)/SMCAV)
4489 DD = DD+DMAX(KS)
4490 END DO
4491
4492 ! ----------------------------------------------------------------------
4493 ! VAL = (1.-EXP(-KDT*SQRT(DT1)))
4494 ! IN BELOW, REMOVE THE SQRT IN ABOVE
4495 ! ----------------------------------------------------------------------
4496 VAL = (1.-EXP(-KDT*DT1))
4497 DDT = DD*VAL
4498 PX = PCPDRP*DT
4499 IF (PX .LT. 0.0) PX = 0.0
4500 INFMAX = (PX*(DDT/(PX+DDT)))/DT
4501
4502 ! ----------------------------------------------------------------------
4503 ! FROZEN GROUND VERSION:
4504 ! REDUCTION OF INFILTRATION BASED ON FROZEN GROUND PARAMETERS
4505 ! ----------------------------------------------------------------------
4506 FCR = 1.
4507 IF (DICE .GT. 1.E-2) THEN
4508 ACRT = CVFRZ * FRZX / DICE
4509 SUM = 1.
4510 IALP1 = CVFRZ - 1
4511 DO J = 1,IALP1
4512 K = 1
4513 DO JJ = J+1,IALP1
4514 K = K * JJ
4515 END DO
4516 SUM = SUM + (ACRT ** ( CVFRZ-J)) / FLOAT (K)
4517 END DO
4518 FCR = 1. - EXP(-ACRT) * SUM
4519 ENDIF
4520 INFMAX = INFMAX * FCR
4521
4522 ! ----------------------------------------------------------------------
4523 ! CORRECTION OF INFILTRATION LIMITATION:
4524 ! IF INFMAX .LE. HYDROLIC CONDUCTIVITY ASSIGN INFMAX THE VALUE OF
4525 ! HYDROLIC CONDUCTIVITY
4526 ! ----------------------------------------------------------------------
4527 ! MXSMC = MAX ( SH2OA(1), SH2OA(2) )
4528 MXSMC = SH2OA(1)
4529
4530 CALL WDFCND (WDF,WCND,MXSMC,SMCMAX,BEXP,DKSAT,DWSAT, &
4531 & SICEMAX)
4532
4533 INFMAX = MAX(INFMAX,WCND)
4534 INFMAX = MIN(INFMAX,PX)
4535
4536 IF (PCPDRP .GT. INFMAX) THEN
4537 RUNOFF1 = PCPDRP - INFMAX
4538 PDDUM = INFMAX
4539 ENDIF
4540
4541 ENDIF
4542
4543 ! ----------------------------------------------------------------------
4544 ! TO AVOID SPURIOUS DRAINAGE BEHAVIOR, 'UPSTREAM DIFFERENCING' IN LINE
4545 ! BELOW REPLACED WITH NEW APPROACH IN 2ND LINE:
4546 ! 'MXSMC = MAX(SH2OA(1), SH2OA(2))'
4547 ! ----------------------------------------------------------------------
4548 MXSMC = SH2OA(1)
4549
4550 CALL WDFCND (WDF,WCND,MXSMC,SMCMAX,BEXP,DKSAT,DWSAT, &
4551 & SICEMAX)
4552
4553 ! ----------------------------------------------------------------------
4554 ! CALC THE MATRIX COEFFICIENTS AI, BI, AND CI FOR THE TOP LAYER
4555 ! ----------------------------------------------------------------------
4556 DDZ = 1. / ( -.5 * ZSOIL(2) )
4557 AI(1) = 0.0
4558 BI(1) = WDF * DDZ / ( -ZSOIL(1) )
4559 CI(1) = -BI(1)
4560
4561 ! ----------------------------------------------------------------------
4562 ! CALC RHSTT FOR THE TOP LAYER AFTER CALC'NG THE VERTICAL SOIL MOISTURE
4563 ! GRADIENT BTWN THE TOP AND NEXT TO TOP LAYERS.
4564 ! ----------------------------------------------------------------------
4565 DSMDZ = ( SH2O(1) - SH2O(2) ) / ( -.5 * ZSOIL(2) )
4566 RHSTT(1) = (WDF * DSMDZ + WCND - PDDUM + EDIR + ET(1))/ZSOIL(1)
4567 SSTT = WDF * DSMDZ + WCND + EDIR + ET(1)
4568
4569 ! ----------------------------------------------------------------------
4570 ! INITIALIZE DDZ2
4571 ! ----------------------------------------------------------------------
4572 DDZ2 = 0.0
4573
4574 ! ----------------------------------------------------------------------
4575 ! LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABV PROCESS
4576 ! ----------------------------------------------------------------------
4577 DO K = 2,NSOIL
4578 DENOM2 = (ZSOIL(K-1) - ZSOIL(K))
4579 IF (K .NE. NSOIL) THEN
4580 SLOPX = 1.
4581
4582 ! ----------------------------------------------------------------------
4583 ! AGAIN, TO AVOID SPURIOUS DRAINAGE BEHAVIOR, 'UPSTREAM DIFFERENCING' IN
4584 ! LINE BELOW REPLACED WITH NEW APPROACH IN 2ND LINE:
4585 ! 'MXSMC2 = MAX (SH2OA(K), SH2OA(K+1))'
4586 ! ----------------------------------------------------------------------
4587 MXSMC2 = SH2OA(K)
4588
4589 CALL WDFCND (WDF2,WCND2,MXSMC2,SMCMAX,BEXP,DKSAT,DWSAT, &
4590 & SICEMAX)
4591
4592 ! ----------------------------------------------------------------------
4593 ! CALC SOME PARTIAL PRODUCTS FOR LATER USE IN CALC'NG RHSTT
4594 ! ----------------------------------------------------------------------
4595 DENOM = (ZSOIL(K-1) - ZSOIL(K+1))
4596 DSMDZ2 = (SH2O(K) - SH2O(K+1)) / (DENOM * 0.5)
4597
4598 ! ----------------------------------------------------------------------
4599 ! CALC THE MATRIX COEF, CI, AFTER CALC'NG ITS PARTIAL PRODUCT
4600 ! ----------------------------------------------------------------------
4601 DDZ2 = 2.0 / DENOM
4602 CI(K) = -WDF2 * DDZ2 / DENOM2
4603 ELSE
4604
4605 ! ----------------------------------------------------------------------
4606 ! SLOPE OF BOTTOM LAYER IS INTRODUCED
4607 ! ----------------------------------------------------------------------
4608 SLOPX = SLOPE
4609
4610 ! ----------------------------------------------------------------------
4611 ! RETRIEVE THE SOIL WATER DIFFUSIVITY AND HYDRAULIC CONDUCTIVITY FOR
4612 ! THIS LAYER
4613 ! ----------------------------------------------------------------------
4614 CALL WDFCND (WDF2,WCND2,SH2OA(NSOIL),SMCMAX,BEXP,DKSAT,DWSAT, &
4615 & SICEMAX)
4616
4617 ! ----------------------------------------------------------------------
4618 ! CALC A PARTIAL PRODUCT FOR LATER USE IN CALC'NG RHSTT
4619 ! ----------------------------------------------------------------------
4620 DSMDZ2 = 0.0
4621
4622 ! ----------------------------------------------------------------------
4623 ! SET MATRIX COEF CI TO ZERO
4624 ! ----------------------------------------------------------------------
4625 CI(K) = 0.0
4626 ENDIF
4627
4628 ! ----------------------------------------------------------------------
4629 ! CALC RHSTT FOR THIS LAYER AFTER CALC'NG ITS NUMERATOR
4630 ! ----------------------------------------------------------------------
4631 NUMER = (WDF2 * DSMDZ2) + SLOPX * WCND2 - (WDF * DSMDZ) &
4632 & - WCND + ET(K)
4633 RHSTT(K) = NUMER / (-DENOM2)
4634
4635 ! ----------------------------------------------------------------------
4636 ! CALC MATRIX COEFS, AI, AND BI FOR THIS LAYER
4637 ! ----------------------------------------------------------------------
4638 AI(K) = -WDF * DDZ / DENOM2
4639 BI(K) = -( AI(K) + CI(K) )
4640
4641 ! ----------------------------------------------------------------------
4642 ! RESET VALUES OF WDF, WCND, DSMDZ, AND DDZ FOR LOOP TO NEXT LYR
4643 ! RUNOFF2: SUB-SURFACE OR BASEFLOW RUNOFF
4644 ! ----------------------------------------------------------------------
4645 IF (K .EQ. NSOIL) THEN
4646 RUNOFF2 = SLOPX * WCND2
4647 ENDIF
4648
4649 IF (K .NE. NSOIL) THEN
4650 WDF = WDF2
4651 WCND = WCND2
4652 DSMDZ = DSMDZ2
4653 DDZ = DDZ2
4654 ENDIF
4655 END DO
4656
4657 ! ----------------------------------------------------------------------
4658 ! END SUBROUTINE SRT
4659 ! ----------------------------------------------------------------------
4660 END SUBROUTINE SRT
4661
4662 SUBROUTINE SSTEP (SH2OOUT,SH2OIN,CMC,RHSTT,RHSCT,DT, &
4663 & NSOIL,SMCMAX,CMCMAX,RUNOFF3,ZSOIL,SMC,SICE, &
4664 & AI,BI,CI)
4665
4666 IMPLICIT NONE
4667
4668 ! ----------------------------------------------------------------------
4669 ! SUBROUTINE SSTEP
4670 ! ----------------------------------------------------------------------
4671 ! CALCULATE/UPDATE SOIL MOISTURE CONTENT VALUES AND CANOPY MOISTURE
4672 ! CONTENT VALUES.
4673 ! ----------------------------------------------------------------------
4674 INTEGER NSOLD
4675 PARAMETER(NSOLD = 20)
4676
4677 INTEGER I
4678 INTEGER K
4679 INTEGER KK11
4680 INTEGER NSOIL
4681
4682 REAL AI(NSOLD)
4683 REAL BI(NSOLD)
4684 REAL CI(NSOLD)
4685 REAL CIin(NSOLD)
4686 REAL CMC
4687 REAL CMCMAX
4688 REAL DDZ
4689 REAL DT
4690 REAL RHSCT
4691 REAL RHSTT(NSOIL)
4692 REAL RHSTTin(NSOIL)
4693 REAL RUNOFF3
4694 REAL SH2OIN(NSOIL)
4695 REAL SH2OOUT(NSOIL)
4696 REAL SICE(NSOIL)
4697 REAL SMC(NSOIL)
4698 REAL SMCMAX
4699 REAL STOT
4700 REAL WPLUS
4701 REAL ZSOIL(NSOIL)
4702
4703 ! ----------------------------------------------------------------------
4704 ! CREATE 'AMOUNT' VALUES OF VARIABLES TO BE INPUT TO THE
4705 ! TRI-DIAGONAL MATRIX ROUTINE.
4706 ! ----------------------------------------------------------------------
4707 DO K = 1,NSOIL
4708 RHSTT(K) = RHSTT(K) * DT
4709 AI(K) = AI(K) * DT
4710 BI(K) = 1. + BI(K) * DT
4711 CI(K) = CI(K) * DT
4712 END DO
4713
4714 ! ----------------------------------------------------------------------
4715 ! COPY VALUES FOR INPUT VARIABLES BEFORE CALL TO ROSR12
4716 ! ----------------------------------------------------------------------
4717 DO K = 1,NSOIL
4718 RHSTTin(K) = RHSTT(K)
4719 END DO
4720 DO K = 1,NSOIL
4721 CIin(K) = CI(K)
4722 END DO
4723
4724 ! ----------------------------------------------------------------------
4725 ! CALL ROSR12 TO SOLVE THE TRI-DIAGONAL MATRIX
4726 ! ----------------------------------------------------------------------
4727 CALL ROSR12 (CI,AI,BI,CIin,RHSTTin,RHSTT,NSOIL)
4728
4729 ! ----------------------------------------------------------------------
4730 ! SUM THE PREVIOUS SMC VALUE AND THE MATRIX SOLUTION TO GET A
4731 ! NEW VALUE. MIN ALLOWABLE VALUE OF SMC WILL BE 0.02.
4732 ! RUNOFF3: RUNOFF WITHIN SOIL LAYERS
4733 ! ----------------------------------------------------------------------
4734 WPLUS = 0.0
4735 RUNOFF3 = 0.
4736 DDZ = -ZSOIL(1)
4737
4738 DO K = 1,NSOIL
4739 IF (K .NE. 1) DDZ = ZSOIL(K - 1) - ZSOIL(K)
4740 SH2OOUT(K) = SH2OIN(K) + CI(K) + WPLUS / DDZ
4741
4742 STOT = SH2OOUT(K) + SICE(K)
4743 IF (STOT .GT. SMCMAX) THEN
4744 IF (K .EQ. 1) THEN
4745 DDZ = -ZSOIL(1)
4746 ELSE
4747 KK11 = K - 1
4748 DDZ = -ZSOIL(K) + ZSOIL(KK11)
4749 ENDIF
4750 WPLUS = (STOT-SMCMAX) * DDZ
4751 ELSE
4752 WPLUS = 0.
4753 ENDIF
4754 SMC(K) = MAX ( MIN(STOT,SMCMAX),0.02 )
4755 SH2OOUT(K) = MAX((SMC(K)-SICE(K)),0.0)
4756 END DO
4757
4758 RUNOFF3 = WPLUS
4759
4760 ! ----------------------------------------------------------------------
4761 ! UPDATE CANOPY WATER CONTENT/INTERCEPTION (CMC). CONVERT RHSCT TO
4762 ! AN 'AMOUNT' VALUE AND ADD TO PREVIOUS CMC VALUE TO GET NEW CMC.
4763 ! ----------------------------------------------------------------------
4764 CMC = CMC + DT * RHSCT
4765 IF (CMC .LT. 1.E-20) CMC=0.0
4766 CMC = MIN(CMC,CMCMAX)
4767
4768 ! ----------------------------------------------------------------------
4769 ! END SUBROUTINE SSTEP
4770 ! ----------------------------------------------------------------------
4771 END SUBROUTINE SSTEP
4772
4773 SUBROUTINE TBND (TU,TB,ZSOIL,ZBOT,K,NSOIL,TBND1)
4774
4775 IMPLICIT NONE
4776
4777 ! ----------------------------------------------------------------------
4778 ! SUBROUTINE TBND
4779 ! ----------------------------------------------------------------------
4780 ! CALCULATE TEMPERATURE ON THE BOUNDARY OF THE LAYER BY INTERPOLATION OF
4781 ! THE MIDDLE LAYER TEMPERATURES
4782 ! ----------------------------------------------------------------------
4783 INTEGER NSOIL
4784 INTEGER K
4785
4786 REAL TBND1
4787 REAL T0
4788 REAL TU
4789 REAL TB
4790 REAL ZB
4791 REAL ZBOT
4792 REAL ZUP
4793 REAL ZSOIL (NSOIL)
4794
4795 PARAMETER(T0 = 273.15)
4796
4797 ! ----------------------------------------------------------------------
4798 ! USE SURFACE TEMPERATURE ON THE TOP OF THE FIRST LAYER
4799 ! ----------------------------------------------------------------------
4800 IF (K .EQ. 1) THEN
4801 ZUP = 0.
4802 ELSE
4803 ZUP = ZSOIL(K-1)
4804 ENDIF
4805
4806 ! ----------------------------------------------------------------------
4807 ! USE DEPTH OF THE CONSTANT BOTTOM TEMPERATURE WHEN INTERPOLATE
4808 ! TEMPERATURE INTO THE LAST LAYER BOUNDARY
4809 ! ----------------------------------------------------------------------
4810 IF (K .EQ. NSOIL) THEN
4811 ZB = 2.*ZBOT-ZSOIL(K)
4812 ELSE
4813 ZB = ZSOIL(K+1)
4814 ENDIF
4815
4816 ! ----------------------------------------------------------------------
4817 ! LINEAR INTERPOLATION BETWEEN THE AVERAGE LAYER TEMPERATURES
4818 ! ----------------------------------------------------------------------
4819 TBND1 = TU+(TB-TU)*(ZUP-ZSOIL(K))/(ZUP-ZB)
4820
4821 ! ----------------------------------------------------------------------
4822 ! END SUBROUTINE TBND
4823 ! ----------------------------------------------------------------------
4824 END SUBROUTINE TBND
4825
4826 SUBROUTINE TDFCND ( DF, SMC, QZ, SMCMAX, SH2O)
4827
4828 IMPLICIT NONE
4829
4830 ! ----------------------------------------------------------------------
4831 ! SUBROUTINE TDFCND
4832 ! ----------------------------------------------------------------------
4833 ! CALCULATE THERMAL DIFFUSIVITY AND CONDUCTIVITY OF THE SOIL FOR A GIVEN
4834 ! POINT AND TIME.
4835 ! ----------------------------------------------------------------------
4836 ! PETERS-LIDARD APPROACH (PETERS-LIDARD et al., 1998)
4837 ! June 2001 CHANGES: FROZEN SOIL CONDITION.
4838 ! ----------------------------------------------------------------------
4839 REAL DF
4840 REAL GAMMD
4841 REAL THKDRY
4842 REAL AKE
4843 REAL THKICE
4844 REAL THKO
4845 REAL THKQTZ
4846 REAL THKSAT
4847 REAL THKS
4848 REAL THKW
4849 REAL QZ
4850 REAL SATRATIO
4851 REAL SH2O
4852 REAL SMC
4853 REAL SMCMAX
4854 REAL XU
4855 REAL XUNFROZ
4856
4857 ! ----------------------------------------------------------------------
4858 ! WE NOW GET QUARTZ AS AN INPUT ARGUMENT (SET IN ROUTINE REDPRM):
4859 ! DATA QUARTZ /0.82, 0.10, 0.25, 0.60, 0.52,
4860 ! & 0.35, 0.60, 0.40, 0.82/
4861 ! ----------------------------------------------------------------------
4862 ! IF THE SOIL HAS ANY MOISTURE CONTENT COMPUTE A PARTIAL SUM/PRODUCT
4863 ! OTHERWISE USE A CONSTANT VALUE WHICH WORKS WELL WITH MOST SOILS
4864 ! ----------------------------------------------------------------------
4865 ! THKW ......WATER THERMAL CONDUCTIVITY
4866 ! THKQTZ ....THERMAL CONDUCTIVITY FOR QUARTZ
4867 ! THKO ......THERMAL CONDUCTIVITY FOR OTHER SOIL COMPONENTS
4868 ! THKS ......THERMAL CONDUCTIVITY FOR THE SOLIDS COMBINED(QUARTZ+OTHER)
4869 ! THKICE ....ICE THERMAL CONDUCTIVITY
4870 ! SMCMAX ....POROSITY (= SMCMAX)
4871 ! QZ .........QUARTZ CONTENT (SOIL TYPE DEPENDENT)
4872 ! ----------------------------------------------------------------------
4873 ! USE AS IN PETERS-LIDARD, 1998 (MODIF. FROM JOHANSEN, 1975).
4874 !
4875 ! PABLO GRUNMANN, 08/17/98
4876 ! REFS.:
4877 ! FAROUKI, O.T.,1986: THERMAL PROPERTIES OF SOILS. SERIES ON ROCK
4878 ! AND SOIL MECHANICS, VOL. 11, TRANS TECH, 136 PP.
4879 ! JOHANSEN, O., 1975: THERMAL CONDUCTIVITY OF SOILS. PH.D. THESIS,
4880 ! UNIVERSITY OF TRONDHEIM,
4881 ! PETERS-LIDARD, C. D., ET AL., 1998: THE EFFECT OF SOIL THERMAL
4882 ! CONDUCTIVITY PARAMETERIZATION ON SURFACE ENERGY FLUXES
4883 ! AND TEMPERATURES. JOURNAL OF THE ATMOSPHERIC SCIENCES,
4884 ! VOL. 55, PP. 1209-1224.
4885 ! ----------------------------------------------------------------------
4886 ! NEEDS PARAMETERS
4887 ! POROSITY(SOIL TYPE):
4888 ! POROS = SMCMAX
4889 ! SATURATION RATIO:
4890 SATRATIO = SMC/SMCMAX
4891
4892 ! PARAMETERS W/(M.K)
4893 THKICE = 2.2
4894 THKW = 0.57
4895 THKO = 2.0
4896 ! IF (QZ .LE. 0.2) THKO = 3.0
4897 THKQTZ = 7.7
4898 ! SOLIDS' CONDUCTIVITY
4899 THKS = (THKQTZ**QZ)*(THKO**(1.- QZ))
4900
4901 ! UNFROZEN FRACTION (FROM 1., i.e., 100%LIQUID, TO 0. (100% FROZEN))
4902 XUNFROZ = (SH2O + 1.E-9) / (SMC + 1.E-9)
4903
4904 ! UNFROZEN VOLUME FOR SATURATION (POROSITY*XUNFROZ)
4905 XU=XUNFROZ*SMCMAX
4906 ! SATURATED THERMAL CONDUCTIVITY
4907 THKSAT = THKS**(1.-SMCMAX)*THKICE**(SMCMAX-XU)*THKW**(XU)
4908
4909 ! DRY DENSITY IN KG/M3
4910 GAMMD = (1. - SMCMAX)*2700.
4911
4912 ! DRY THERMAL CONDUCTIVITY IN W.M-1.K-1
4913 THKDRY = (0.135*GAMMD + 64.7)/(2700. - 0.947*GAMMD)
4914
4915 IF ( (SH2O + 0.0005) .LT. SMC ) THEN
4916 ! FROZEN
4917 AKE = SATRATIO
4918 ELSE
4919 ! UNFROZEN
4920 ! RANGE OF VALIDITY FOR THE KERSTEN NUMBER (AKE)
4921 IF ( SATRATIO .GT. 0.1 ) THEN
4922
4923 ! KERSTEN NUMBER (USING "FINE" FORMULA, VALID FOR SOILS CONTAINING AT
4924 ! LEAST 5% OF PARTICLES WITH DIAMETER LESS THAN 2.E-6 METERS.)
4925 ! (FOR "COARSE" FORMULA, SEE PETERS-LIDARD ET AL., 1998).
4926
4927 AKE = LOG10(SATRATIO) + 1.0
4928
4929 ELSE
4930
4931 ! USE K = KDRY
4932 AKE = 0.0
4933
4934 ENDIF
4935 ENDIF
4936
4937 ! THERMAL CONDUCTIVITY
4938
4939 DF = AKE*(THKSAT - THKDRY) + THKDRY
4940
4941 ! ----------------------------------------------------------------------
4942 ! END SUBROUTINE TDFCND
4943 ! ----------------------------------------------------------------------
4944 END SUBROUTINE TDFCND
4945
4946 SUBROUTINE TMPAVG (TAVG,TUP,TM,TDN,ZSOIL,NSOIL,K)
4947
4948 IMPLICIT NONE
4949
4950 ! ----------------------------------------------------------------------
4951 ! SUBROUTINE TMPAVG
4952 ! ----------------------------------------------------------------------
4953 ! CALCULATE SOIL LAYER AVERAGE TEMPERATURE (TAVG) IN FREEZING/THAWING
4954 ! LAYER USING UP, DOWN, AND MIDDLE LAYER TEMPERATURES (TUP, TDN, TM),
4955 ! WHERE TUP IS AT TOP BOUNDARY OF LAYER, TDN IS AT BOTTOM BOUNDARY OF
4956 ! LAYER. TM IS LAYER PROGNOSTIC STATE TEMPERATURE.
4957 ! ----------------------------------------------------------------------
4958 INTEGER K
4959 INTEGER NSOIL
4960
4961 REAL DZ
4962 REAL DZH
4963 REAL T0
4964 REAL TAVG
4965 REAL TDN
4966 REAL TM
4967 REAL TUP
4968 REAL X0
4969 REAL XDN
4970 REAL XUP
4971 REAL ZSOIL (NSOIL)
4972
4973 PARAMETER(T0 = 2.7315E2)
4974
4975 ! ----------------------------------------------------------------------
4976 IF (K .EQ. 1) THEN
4977 DZ = -ZSOIL(1)
4978 ELSE
4979 DZ = ZSOIL(K-1)-ZSOIL(K)
4980 ENDIF
4981
4982 DZH=DZ*0.5
4983
4984 IF (TUP .LT. T0) THEN
4985 IF (TM .LT. T0) THEN
4986 IF (TDN .LT. T0) THEN
4987 ! ----------------------------------------------------------------------
4988 ! TUP, TM, TDN < T0
4989 ! ----------------------------------------------------------------------
4990 TAVG = (TUP + 2.0*TM + TDN)/ 4.0
4991 ELSE
4992 ! ----------------------------------------------------------------------
4993 ! TUP & TM < T0, TDN >= T0
4994 ! ----------------------------------------------------------------------
4995 X0 = (T0 - TM) * DZH / (TDN - TM)
4996 TAVG = 0.5 * (TUP*DZH+TM*(DZH+X0)+T0*(2.*DZH-X0)) / DZ
4997 ENDIF
4998 ELSE
4999 IF (TDN .LT. T0) THEN
5000 ! ----------------------------------------------------------------------
5001 ! TUP < T0, TM >= T0, TDN < T0
5002 ! ----------------------------------------------------------------------
5003 XUP = (T0-TUP) * DZH / (TM-TUP)
5004 XDN = DZH - (T0-TM) * DZH / (TDN-TM)
5005 TAVG = 0.5 * (TUP*XUP+T0*(2.*DZ-XUP-XDN)+TDN*XDN) / DZ
5006 ELSE
5007 ! ----------------------------------------------------------------------
5008 ! TUP < T0, TM >= T0, TDN >= T0
5009 ! ----------------------------------------------------------------------
5010 XUP = (T0-TUP) * DZH / (TM-TUP)
5011 TAVG = 0.5 * (TUP*XUP+T0*(2.*DZ-XUP)) / DZ
5012 ENDIF
5013 ENDIF
5014 ELSE
5015 IF (TM .LT. T0) THEN
5016 IF (TDN .LT. T0) THEN
5017 ! ----------------------------------------------------------------------
5018 ! TUP >= T0, TM < T0, TDN < T0
5019 ! ----------------------------------------------------------------------
5020 XUP = DZH - (T0-TUP) * DZH / (TM-TUP)
5021 TAVG = 0.5 * (T0*(DZ-XUP)+TM*(DZH+XUP)+TDN*DZH) / DZ
5022 ELSE
5023 ! ----------------------------------------------------------------------
5024 ! TUP >= T0, TM < T0, TDN >= T0
5025 ! ----------------------------------------------------------------------
5026 XUP = DZH - (T0-TUP) * DZH / (TM-TUP)
5027 XDN = (T0-TM) * DZH / (TDN-TM)
5028 TAVG = 0.5 * (T0*(2.*DZ-XUP-XDN)+TM*(XUP+XDN)) / DZ
5029 ENDIF
5030 ELSE
5031 IF (TDN .LT. T0) THEN
5032 ! ----------------------------------------------------------------------
5033 ! TUP >= T0, TM >= T0, TDN < T0
5034 ! ----------------------------------------------------------------------
5035 XDN = DZH - (T0-TM) * DZH / (TDN-TM)
5036 TAVG = (T0*(DZ-XDN)+0.5*(T0+TDN)*XDN) / DZ
5037 ELSE
5038 ! ----------------------------------------------------------------------
5039 ! TUP >= T0, TM >= T0, TDN >= T0
5040 ! ----------------------------------------------------------------------
5041 TAVG = (TUP + 2.0*TM + TDN) / 4.0
5042 ENDIF
5043 ENDIF
5044 ENDIF
5045 ! ----------------------------------------------------------------------
5046 ! END SUBROUTINE TMPAVG
5047 ! ----------------------------------------------------------------------
5048 END SUBROUTINE TMPAVG
5049
5050 SUBROUTINE TRANSP (ET1,NSOIL,ETP1,SMC,CMC,ZSOIL,SHDFAC,SMCWLT, &
5051 & CMCMAX,PC,CFACTR,SMCREF,SFCTMP,Q2,NROOT,RTDIS)
5052
5053 IMPLICIT NONE
5054
5055 ! ----------------------------------------------------------------------
5056 ! SUBROUTINE TRANSP
5057 ! ----------------------------------------------------------------------
5058 ! CALCULATE TRANSPIRATION FOR THE VEG CLASS.
5059 ! ----------------------------------------------------------------------
5060 INTEGER I
5061 INTEGER K
5062 INTEGER NSOIL
5063 INTEGER NROOT
5064
5065 REAL CFACTR
5066 REAL CMC
5067 REAL CMCMAX
5068 REAL DENOM
5069 REAL ET1(NSOIL)
5070 REAL ETP1
5071 REAL ETP1A
5072 REAL GX (7)
5073 !.....REAL PART(NSOIL)
5074 REAL PC
5075 REAL Q2
5076 REAL RTDIS(NSOIL)
5077 REAL RTX
5078 REAL SFCTMP
5079 REAL SGX
5080 REAL SHDFAC
5081 REAL SMC(NSOIL)
5082 REAL SMCREF
5083 REAL SMCWLT
5084 REAL ZSOIL(NSOIL)
5085
5086 ! ----------------------------------------------------------------------
5087 ! INITIALIZE PLANT TRANSP TO ZERO FOR ALL SOIL LAYERS.
5088 ! ----------------------------------------------------------------------
5089 DO K = 1,NSOIL
5090 ET1(K) = 0.
5091 END DO
5092
5093 ! ----------------------------------------------------------------------
5094 ! CALCULATE AN 'ADJUSTED' POTENTIAL TRANSPIRATION
5095 ! IF STATEMENT BELOW TO AVOID TANGENT LINEAR PROBLEMS NEAR ZERO
5096 ! NOTE: GX AND OTHER TERMS BELOW REDISTRIBUTE TRANSPIRATION BY LAYER,
5097 ! ET(K), AS A FUNCTION OF SOIL MOISTURE AVAILABILITY, WHILE PRESERVING
5098 ! TOTAL ETP1A.
5099 ! ----------------------------------------------------------------------
5100 IF (CMC .NE. 0.0) THEN
5101 ETP1A = SHDFAC * PC * ETP1 * (1.0 - (CMC /CMCMAX) ** CFACTR)
5102 ELSE
5103 ETP1A = SHDFAC * PC * ETP1
5104 ENDIF
5105
5106 SGX = 0.0
5107 DO I = 1,NROOT
5108 GX(I) = ( SMC(I) - SMCWLT ) / ( SMCREF - SMCWLT )
5109 GX(I) = MAX ( MIN ( GX(I), 1. ), 0. )
5110 SGX = SGX + GX (I)
5111 END DO
5112 SGX = SGX / NROOT
5113
5114 DENOM = 0.
5115 DO I = 1,NROOT
5116 RTX = RTDIS(I) + GX(I) - SGX
5117 GX(I) = GX(I) * MAX ( RTX, 0. )
5118 DENOM = DENOM + GX(I)
5119 END DO
5120 IF (DENOM .LE. 0.0) DENOM = 1.
5121
5122 DO I = 1,NROOT
5123 ET1(I) = ETP1A * GX(I) / DENOM
5124 END DO
5125
5126 ! ----------------------------------------------------------------------
5127 ! ABOVE CODE ASSUMES A VERTICALLY UNIFORM ROOT DISTRIBUTION
5128 ! CODE BELOW TESTS A VARIABLE ROOT DISTRIBUTION
5129 ! ----------------------------------------------------------------------
5130 ! ET(1) = ( ZSOIL(1) / ZSOIL(NROOT) ) * GX * ETP1A
5131 ! ET(1) = ( ZSOIL(1) / ZSOIL(NROOT) ) * ETP1A
5132 ! ----------------------------------------------------------------------
5133 ! USING ROOT DISTRIBUTION AS WEIGHTING FACTOR
5134 ! ----------------------------------------------------------------------
5135 ! ET(1) = RTDIS(1) * ETP1A
5136 ! ET(1) = ETP1A * PART(1)
5137 ! ----------------------------------------------------------------------
5138 ! LOOP DOWN THRU THE SOIL LAYERS REPEATING THE OPERATION ABOVE,
5139 ! BUT USING THE THICKNESS OF THE SOIL LAYER (RATHER THAN THE
5140 ! ABSOLUTE DEPTH OF EACH LAYER) IN THE FINAL CALCULATION.
5141 ! ----------------------------------------------------------------------
5142 ! DO K = 2,NROOT
5143 ! GX = ( SMC(K) - SMCWLT ) / ( SMCREF - SMCWLT )
5144 ! GX = MAX ( MIN ( GX, 1. ), 0. )
5145 ! TEST CANOPY RESISTANCE
5146 ! GX = 1.0
5147 ! ET(K) = ((ZSOIL(K)-ZSOIL(K-1))/ZSOIL(NROOT))*GX*ETP1A
5148 ! ET(K) = ((ZSOIL(K)-ZSOIL(K-1))/ZSOIL(NROOT))*ETP1A
5149 ! ----------------------------------------------------------------------
5150 ! USING ROOT DISTRIBUTION AS WEIGHTING FACTOR
5151 ! ----------------------------------------------------------------------
5152 ! ET(K) = RTDIS(K) * ETP1A
5153 ! ET(K) = ETP1A*PART(K)
5154 ! END DO
5155 ! ----------------------------------------------------------------------
5156 ! END SUBROUTINE TRANSP
5157 ! ----------------------------------------------------------------------
5158 END SUBROUTINE TRANSP
5159
5160 SUBROUTINE WDFCND (WDF,WCND,SMC,SMCMAX,BEXP,DKSAT,DWSAT, &
5161 & SICEMAX)
5162
5163 IMPLICIT NONE
5164
5165 ! ----------------------------------------------------------------------
5166 ! SUBROUTINE WDFCND
5167 ! ----------------------------------------------------------------------
5168 ! CALCULATE SOIL WATER DIFFUSIVITY AND SOIL HYDRAULIC CONDUCTIVITY.
5169 ! ----------------------------------------------------------------------
5170 REAL BEXP
5171 REAL DKSAT
5172 REAL DWSAT
5173 REAL EXPON
5174 REAL FACTR1
5175 REAL FACTR2
5176 REAL SICEMAX
5177 REAL SMC
5178 REAL SMCMAX
5179 REAL VKwgt
5180 REAL WCND
5181 REAL WDF
5182
5183 ! ----------------------------------------------------------------------
5184 ! CALC THE RATIO OF THE ACTUAL TO THE MAX PSBL SOIL H2O CONTENT
5185 ! ----------------------------------------------------------------------
5186 SMC = SMC
5187 SMCMAX = SMCMAX
5188 FACTR1 = 0.2 / SMCMAX
5189 FACTR2 = SMC / SMCMAX
5190
5191 ! ----------------------------------------------------------------------
5192 ! PREP AN EXPNTL COEF AND CALC THE SOIL WATER DIFFUSIVITY
5193 ! ----------------------------------------------------------------------
5194 EXPON = BEXP + 2.0
5195 WDF = DWSAT * FACTR2 ** EXPON
5196
5197 ! ----------------------------------------------------------------------
5198 ! FROZEN SOIL HYDRAULIC DIFFUSIVITY. VERY SENSITIVE TO THE VERTICAL
5199 ! GRADIENT OF UNFROZEN WATER. THE LATTER GRADIENT CAN BECOME VERY
5200 ! EXTREME IN FREEZING/THAWING SITUATIONS, AND GIVEN THE RELATIVELY
5201 ! FEW AND THICK SOIL LAYERS, THIS GRADIENT SUFFERES SERIOUS
5202 ! TRUNCTION ERRORS YIELDING ERRONEOUSLY HIGH VERTICAL TRANSPORTS OF
5203 ! UNFROZEN WATER IN BOTH DIRECTIONS FROM HUGE HYDRAULIC DIFFUSIVITY.
5204 ! THEREFORE, WE FOUND WE HAD TO ARBITRARILY CONSTRAIN WDF
5205 ! --
5206 ! VERSION D_10CM: ........ FACTR1 = 0.2/SMCMAX
5207 ! WEIGHTED APPROACH...................... PABLO GRUNMANN, 28_SEP_1999.
5208 ! ----------------------------------------------------------------------
5209 IF (SICEMAX .GT. 0.0) THEN
5210 VKWGT = 1./(1.+(500.*SICEMAX)**3.)
5211 WDF = VKWGT*WDF + (1.- VKWGT)*DWSAT*FACTR1**EXPON
5212 ENDIF
5213
5214 ! ----------------------------------------------------------------------
5215 ! RESET THE EXPNTL COEF AND CALC THE HYDRAULIC CONDUCTIVITY
5216 ! ----------------------------------------------------------------------
5217 EXPON = (2.0 * BEXP) + 3.0
5218 WCND = DKSAT * FACTR2 ** EXPON
5219
5220 ! ----------------------------------------------------------------------
5221 ! END SUBROUTINE WDFCND
5222 ! ----------------------------------------------------------------------
5223 END SUBROUTINE WDFCND
5224
5225 SUBROUTINE nmmlsminit(isn,XICE,VEGFRA,SNOW,SNOWC,CANWAT,SMSTAV, &
5226 SMSTOT, SFCRUNOFF,UDRUNOFF,GRDFLX,ACSNOW, &
5227 ACSNOM,IVGTYP,ISLTYP,TSLB,SMOIS,DZS,SFCEVP, & ! STEMP
5228 TMN, &
5229 num_soil_layers, &
5230 allowed_to_read, &
5231 ids,ide, jds,jde, kds,kde, &
5232 ims,ime, jms,jme, kms,kme, &
5233 its,ite, jts,jte, kts,kte )
5234
5235 IMPLICIT NONE
5236
5237 ! Arguments
5238 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
5239 ims,ime, jms,jme, kms,kme, &
5240 its,ite, jts,jte, kts,kte
5241
5242 INTEGER, INTENT(IN) :: num_soil_layers
5243
5244 REAL, DIMENSION( num_soil_layers), INTENT(IN) :: DZS
5245
5246 REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , &
5247 INTENT(INOUT) :: SMOIS, &
5248 TSLB !STEMP
5249
5250 REAL, DIMENSION( ims:ime, jms:jme ) , &
5251 INTENT(INOUT) :: SNOW, &
5252 SNOWC, &
5253 CANWAT, &
5254 SMSTAV, &
5255 SMSTOT, &
5256 SFCRUNOFF, &
5257 UDRUNOFF, &
5258 SFCEVP, &
5259 GRDFLX, &
5260 ACSNOW, &
5261 XICE, &
5262 VEGFRA, &
5263 TMN, &
5264 ACSNOM
5265
5266 INTEGER, DIMENSION( ims:ime, jms:jme ) , &
5267 INTENT(INOUT) :: IVGTYP, &
5268 ISLTYP
5269
5270 !
5271
5272 INTEGER, INTENT(IN) :: isn
5273 LOGICAL, INTENT(IN) :: allowed_to_read
5274 ! Local
5275 INTEGER :: iseason
5276 INTEGER :: icm,jcm,itf,jtf
5277 INTEGER :: I,J,L
5278
5279
5280 itf=min0(ite,ide-1)
5281 jtf=min0(jte,jde-1)
5282
5283 icm = ide/2
5284 jcm = jde/2
5285
5286 iseason=isn
5287
5288 DO J=jts,jtf
5289 DO I=its,itf
5290 ! SNOW(i,j)=0.
5291 SNOWC(i,j)=0.
5292 ! SMSTAV(i,j)=
5293 ! SMSTOT(i,j)=
5294 ! SFCRUNOFF(i,j)=
5295 ! UDRUNOFF(i,j)=
5296 ! GRDFLX(i,j)=
5297 ! ACSNOW(i,j)=
5298 ! ACSNOM(i,j)=
5299 ENDDO
5300 ENDDO
5301
5302 END SUBROUTINE nmmlsminit
5303
5304 FUNCTION CSNOW (DSNOW)
5305
5306 IMPLICIT NONE
5307
5308 ! ----------------------------------------------------------------------
5309 ! FUNCTION CSNOW
5310 ! ----------------------------------------------------------------------
5311 ! CALCULATE SNOW TERMAL CONDUCTIVITY
5312 ! ----------------------------------------------------------------------
5313 REAL C
5314 REAL DSNOW
5315 REAL CSNOW
5316 REAL UNIT
5317
5318 PARAMETER(UNIT = 0.11631)
5319
5320 ! ----------------------------------------------------------------------
5321 ! CSNOW IN UNITS OF CAL/(CM*HR*C), RETURNED IN W/(M*C)
5322 ! BASIC VERSION IS DYACHKOVA EQUATION (1960), FOR RANGE 0.1-0.4
5323 ! ----------------------------------------------------------------------
5324 C=0.328*10**(2.25*DSNOW)
5325 ! CSNOW=UNIT*C
5326 ! MEK JAN 2006, DOUBLE SNOW THERMAL CONDUCTIVITY
5327 CSNOW=2.0*UNIT*C
5328
5329 ! ----------------------------------------------------------------------
5330 ! DE VAUX EQUATION (1933), IN RANGE 0.1-0.6
5331 ! ----------------------------------------------------------------------
5332 ! CSNOW=0.0293*(1.+100.*DSNOW**2)
5333
5334 ! ----------------------------------------------------------------------
5335 ! E. ANDERSEN FROM FLERCHINGER
5336 ! ----------------------------------------------------------------------
5337 ! CSNOW=0.021+2.51*DSNOW**2
5338
5339 ! ----------------------------------------------------------------------
5340 ! END FUNCTION CSNOW
5341 ! ----------------------------------------------------------------------
5342 END FUNCTION CSNOW
5343
5344 FUNCTION FRH2O (TKELV,SMC,SH2O,SMCMAX,BEXP,PSIS)
5345
5346 IMPLICIT NONE
5347
5348 ! ----------------------------------------------------------------------
5349 ! FUNCTION FRH2O
5350 ! ----------------------------------------------------------------------
5351 ! CALCULATE AMOUNT OF SUPERCOOLED LIQUID SOIL WATER CONTENT IF
5352 ! TEMPERATURE IS BELOW 273.15K (T0). REQUIRES NEWTON-TYPE ITERATION TO
5353 ! SOLVE THE NONLINEAR IMPLICIT EQUATION GIVEN IN EQN 17 OF KOREN ET AL
5354 ! (1999, JGR, VOL 104(D16), 19569-19585).
5355 ! ----------------------------------------------------------------------
5356 ! NEW VERSION (JUNE 2001): MUCH FASTER AND MORE ACCURATE NEWTON
5357 ! ITERATION ACHIEVED BY FIRST TAKING LOG OF EQN CITED ABOVE -- LESS THAN
5358 ! 4 (TYPICALLY 1 OR 2) ITERATIONS ACHIEVES CONVERGENCE. ALSO, EXPLICIT
5359 ! 1-STEP SOLUTION OPTION FOR SPECIAL CASE OF PARAMETER CK=0, WHICH
5360 ! REDUCES THE ORIGINAL IMPLICIT EQUATION TO A SIMPLER EXPLICIT FORM,
5361 ! KNOWN AS THE "FLERCHINGER EQN". IMPROVED HANDLING OF SOLUTION IN THE
5362 ! LIMIT OF FREEZING POINT TEMPERATURE T0.
5363 ! ----------------------------------------------------------------------
5364 ! INPUT:
5365 !
5366 ! TKELV.........TEMPERATURE (Kelvin)
5367 ! SMC...........TOTAL SOIL MOISTURE CONTENT (VOLUMETRIC)
5368 ! SH2O..........LIQUID SOIL MOISTURE CONTENT (VOLUMETRIC)
5369 ! SMCMAX........SATURATION SOIL MOISTURE CONTENT (FROM REDPRM)
5370 ! B.............SOIL TYPE "B" PARAMETER (FROM REDPRM)
5371 ! PSIS..........SATURATED SOIL MATRIC POTENTIAL (FROM REDPRM)
5372 !
5373 ! OUTPUT:
5374 ! FRH2O.........SUPERCOOLED LIQUID WATER CONTENT
5375 ! ----------------------------------------------------------------------
5376 REAL BEXP
5377 REAL BLIM
5378 REAL BX
5379 REAL CK
5380 REAL DENOM
5381 REAL DF
5382 REAL DH2O
5383 REAL DICE
5384 REAL DSWL
5385 REAL ERROR
5386 REAL FK
5387 REAL FRH2O
5388 REAL GS
5389 REAL HLICE
5390 REAL PSIS
5391 REAL SH2O
5392 REAL SMC
5393 REAL SMCMAX
5394 REAL SWL
5395 REAL SWLK
5396 REAL TKELV
5397 REAL T0
5398
5399 INTEGER NLOG
5400 INTEGER KCOUNT
5401
5402 PARAMETER(CK = 8.0)
5403 ! PARAMETER(CK = 0.0)
5404 PARAMETER(BLIM = 5.5)
5405 PARAMETER(ERROR = 0.005)
5406
5407 PARAMETER(HLICE = 3.335E5)
5408 PARAMETER(GS = 9.81)
5409 PARAMETER(DICE = 920.0)
5410 PARAMETER(DH2O = 1000.0)
5411 PARAMETER(T0 = 273.15)
5412
5413 ! ----------------------------------------------------------------------
5414 ! LIMITS ON PARAMETER B: B < 5.5 (use parameter BLIM)
5415 ! SIMULATIONS SHOWED IF B > 5.5 UNFROZEN WATER CONTENT IS
5416 ! NON-REALISTICALLY HIGH AT VERY LOW TEMPERATURES.
5417 ! ----------------------------------------------------------------------
5418 BX = BEXP
5419 IF (BEXP .GT. BLIM) BX = BLIM
5420
5421 ! ----------------------------------------------------------------------
5422 ! INITIALIZING ITERATIONS COUNTER AND ITERATIVE SOLUTION FLAG.
5423 ! ----------------------------------------------------------------------
5424 NLOG=0
5425 KCOUNT=0
5426
5427 ! ----------------------------------------------------------------------
5428 ! IF TEMPERATURE NOT SIGNIFICANTLY BELOW FREEZING (T0), SH2O = SMC
5429 ! ----------------------------------------------------------------------
5430 IF (TKELV .GT. (T0 - 1.E-3)) THEN
5431 FRH2O = SMC
5432 ELSE
5433 IF (CK .NE. 0.0) THEN
5434
5435 ! ----------------------------------------------------------------------
5436 ! OPTION 1: ITERATED SOLUTION FOR NONZERO CK
5437 ! IN KOREN ET AL, JGR, 1999, EQN 17
5438 ! ----------------------------------------------------------------------
5439 ! INITIAL GUESS FOR SWL (frozen content)
5440 ! ----------------------------------------------------------------------
5441 SWL = SMC-SH2O
5442
5443 ! ----------------------------------------------------------------------
5444 ! KEEP WITHIN BOUNDS.
5445 ! ----------------------------------------------------------------------
5446 IF (SWL .GT. (SMC-0.02)) SWL = SMC-0.02
5447 IF (SWL .LT. 0.) SWL = 0.
5448
5449 ! ----------------------------------------------------------------------
5450 ! START OF ITERATIONS
5451 ! ----------------------------------------------------------------------
5452 DO WHILE ( (NLOG .LT. 10) .AND. (KCOUNT .EQ. 0) )
5453 NLOG = NLOG+1
5454 DF = ALOG(( PSIS*GS/HLICE ) * ( ( 1.+CK*SWL )**2. ) * &
5455 & ( SMCMAX/(SMC-SWL) )**BX) - ALOG(-(TKELV-T0)/TKELV)
5456 DENOM = 2. * CK / ( 1.+CK*SWL ) + BX / ( SMC - SWL )
5457 SWLK = SWL - DF/DENOM
5458 ! ----------------------------------------------------------------------
5459 ! BOUNDS USEFUL FOR MATHEMATICAL SOLUTION.
5460 ! ----------------------------------------------------------------------
5461 IF (SWLK .GT. (SMC-0.02)) SWLK = SMC - 0.02
5462 IF (SWLK .LT. 0.) SWLK = 0.
5463
5464 ! ----------------------------------------------------------------------
5465 ! MATHEMATICAL SOLUTION BOUNDS APPLIED.
5466 ! ----------------------------------------------------------------------
5467 DSWL = ABS(SWLK-SWL)
5468 SWL = SWLK
5469
5470 ! ----------------------------------------------------------------------
5471 ! IF MORE THAN 10 ITERATIONS, USE EXPLICIT METHOD (CK=0 APPROX.)
5472 ! WHEN DSWL LESS OR EQ. ERROR, NO MORE ITERATIONS REQUIRED.
5473 ! ----------------------------------------------------------------------
5474 IF ( DSWL .LE. ERROR ) THEN
5475 KCOUNT = KCOUNT+1
5476 ENDIF
5477 END DO
5478
5479 ! ----------------------------------------------------------------------
5480 ! END OF ITERATIONS
5481 ! ----------------------------------------------------------------------
5482 ! BOUNDS APPLIED WITHIN DO-BLOCK ARE VALID FOR PHYSICAL SOLUTION.
5483 ! ----------------------------------------------------------------------
5484 FRH2O = SMC - SWL
5485
5486 ! ----------------------------------------------------------------------
5487 ! END OPTION 1
5488 ! ----------------------------------------------------------------------
5489 ENDIF
5490
5491 ! ----------------------------------------------------------------------
5492 ! OPTION 2: EXPLICIT SOLUTION FOR FLERCHINGER EQ. i.e. CK=0
5493 ! IN KOREN ET AL., JGR, 1999, EQN 17
5494 ! APPLY PHYSICAL BOUNDS TO FLERCHINGER SOLUTION
5495 ! ----------------------------------------------------------------------
5496 IF (KCOUNT .EQ. 0) THEN
5497 Print*,'Flerchinger used in NEW version. Iterations=',NLOG
5498 FK = (((HLICE/(GS*(-PSIS)))* &
5499 & ((TKELV-T0)/TKELV))**(-1/BX))*SMCMAX
5500 IF (FK .LT. 0.02) FK = 0.02
5501 FRH2O = MIN (FK, SMC)
5502 ! ----------------------------------------------------------------------
5503 ! END OPTION 2
5504 ! ----------------------------------------------------------------------
5505 ENDIF
5506
5507 ENDIF
5508
5509 ! ----------------------------------------------------------------------
5510 ! END FUNCTION FRH2O
5511 ! ----------------------------------------------------------------------
5512 END FUNCTION FRH2O
5513
5514 FUNCTION SNKSRC (TAVG,SMC,SH2O,ZSOIL,NSOIL, &
5515 & SMCMAX,PSISAT,BEXP,DT,K,QTOT)
5516
5517 IMPLICIT NONE
5518
5519 ! ----------------------------------------------------------------------
5520 ! FUNCTION SNKSRC
5521 ! ----------------------------------------------------------------------
5522 ! CALCULATE SINK/SOURCE TERM OF THE TERMAL DIFFUSION EQUATION. (SH2O) IS
5523 ! AVAILABLE LIQUED WATER.
5524 ! ----------------------------------------------------------------------
5525 INTEGER K
5526 INTEGER NSOIL
5527
5528 REAL BEXP
5529 REAL DF
5530 REAL DH2O
5531 REAL DT
5532 REAL DZ
5533 REAL DZH
5534 REAL FREE
5535 ! REAL FRH2O
5536 REAL HLICE
5537 REAL PSISAT
5538 REAL QTOT
5539 REAL SH2O
5540 REAL SMC
5541 REAL SMCMAX
5542 REAL SNKSRC
5543 REAL T0
5544 REAL TAVG
5545 REAL TDN
5546 REAL TM
5547 REAL TUP
5548 REAL TZ
5549 REAL X0
5550 REAL XDN
5551 REAL XH2O
5552 REAL XUP
5553 REAL ZSOIL (NSOIL)
5554
5555 PARAMETER(DH2O = 1.0000E3)
5556 PARAMETER(HLICE = 3.3350E5)
5557 PARAMETER(T0 = 2.7315E2)
5558
5559 IF (K .EQ. 1) THEN
5560 DZ = -ZSOIL(1)
5561 ELSE
5562 DZ = ZSOIL(K-1)-ZSOIL(K)
5563 ENDIF
5564
5565 ! ----------------------------------------------------------------------
5566 ! VIA FUNCTION FRH2O, COMPUTE POTENTIAL OR 'EQUILIBRIUM' UNFROZEN
5567 ! SUPERCOOLED FREE WATER FOR GIVEN SOIL TYPE AND SOIL LAYER TEMPERATURE.
5568 ! FUNCTION FRH20 INVOKES EQN (17) FROM V. KOREN ET AL (1999, JGR, VOL.
5569 ! 104, PG 19573). (ASIDE: LATTER EQN IN JOURNAL IN CENTIGRADE UNITS.
5570 ! ROUTINE FRH2O USE FORM OF EQN IN KELVIN UNITS.)
5571 ! ----------------------------------------------------------------------
5572 FREE = FRH2O(TAVG,SMC,SH2O,SMCMAX,BEXP,PSISAT)
5573
5574 ! ----------------------------------------------------------------------
5575 ! IN NEXT BLOCK OF CODE, INVOKE EQN 18 OF V. KOREN ET AL (1999, JGR,
5576 ! VOL. 104, PG 19573.) THAT IS, FIRST ESTIMATE THE NEW AMOUNTOF LIQUID
5577 ! WATER, 'XH2O', IMPLIED BY THE SUM OF (1) THE LIQUID WATER AT THE BEGIN
5578 ! OF CURRENT TIME STEP, AND (2) THE FREEZE OF THAW CHANGE IN LIQUID
5579 ! WATER IMPLIED BY THE HEAT FLUX 'QTOT' PASSED IN FROM ROUTINE HRT.
5580 ! SECOND, DETERMINE IF XH2O NEEDS TO BE BOUNDED BY 'FREE' (EQUIL AMT) OR
5581 ! IF 'FREE' NEEDS TO BE BOUNDED BY XH2O.
5582 ! ----------------------------------------------------------------------
5583 XH2O = SH2O + QTOT*DT/(DH2O*HLICE*DZ)
5584
5585 ! ----------------------------------------------------------------------
5586 ! FIRST, IF FREEZING AND REMAINING LIQUID LESS THAN LOWER BOUND, THEN
5587 ! REDUCE EXTENT OF FREEZING, THEREBY LETTING SOME OR ALL OF HEAT FLUX
5588 ! QTOT COOL THE SOIL TEMP LATER IN ROUTINE HRT.
5589 ! ----------------------------------------------------------------------
5590 IF ( XH2O .LT. SH2O .AND. XH2O .LT. FREE) THEN
5591 IF ( FREE .GT. SH2O ) THEN
5592 XH2O = SH2O
5593 ELSE
5594 XH2O = FREE
5595 ENDIF
5596 ENDIF
5597
5598 ! ----------------------------------------------------------------------
5599 ! SECOND, IF THAWING AND THE INCREASE IN LIQUID WATER GREATER THAN UPPER
5600 ! BOUND, THEN REDUCE EXTENT OF THAW, THEREBY LETTING SOME OR ALL OF HEAT
5601 ! FLUX QTOT WARM THE SOIL TEMP LATER IN ROUTINE HRT.
5602 ! ----------------------------------------------------------------------
5603 IF ( XH2O .GT. SH2O .AND. XH2O .GT. FREE ) THEN
5604 IF ( FREE .LT. SH2O ) THEN
5605 XH2O = SH2O
5606 ELSE
5607 XH2O = FREE
5608 ENDIF
5609 ENDIF
5610
5611 IF (XH2O .LT. 0.) XH2O = 0.
5612 IF (XH2O .GT. SMC) XH2O = SMC
5613
5614 ! ----------------------------------------------------------------------
5615 ! CALCULATE PHASE-CHANGE HEAT SOURCE/SINK TERM FOR USE IN ROUTINE HRT
5616 ! AND UPDATE LIQUID WATER TO REFLCET FINAL FREEZE/THAW INCREMENT.
5617 ! ----------------------------------------------------------------------
5618 SNKSRC = -DH2O*HLICE*DZ*(XH2O-SH2O)/DT
5619 SH2O = XH2O
5620
5621 ! ----------------------------------------------------------------------
5622 ! END FUNCTION SNKSRC
5623 ! ----------------------------------------------------------------------
5624 END FUNCTION SNKSRC
5625
5626 END MODULE module_sf_lsm_nmm
5627