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 BETA = ETA/ETP
1382 ! ----------------------------------------------------------------------
1383
1384 ! ----------------------------------------------------------------------
1385 ! CONVERT THE SIGN OF SOIL HEAT FLUX SO THAT:
1386 ! SSOIL>0: WARM THE SURFACE (NIGHT TIME)
1387 ! SSOIL<0: COOL THE SURFACE (DAY TIME)
1388 ! ----------------------------------------------------------------------
1389 SSOIL = -1.0*SSOIL
1390
1391 ! ----------------------------------------------------------------------
1392 ! CONVERT RUNOFF3 (INTERNAL LAYER RUNOFF FROM SUPERSAT) FROM M TO M S-1
1393 ! AND ADD TO SUBSURFACE RUNOFF/DRAINAGE/BASEFLOW
1394 ! ----------------------------------------------------------------------
1395 RUNOFF3 = RUNOFF3/DT
1396 RUNOFF2 = RUNOFF2+RUNOFF3
1397
1398 ! ----------------------------------------------------------------------
1399 ! TOTAL COLUMN SOIL MOISTURE IN METERS (SOILM) AND ROOT-ZONE
1400 ! SOIL MOISTURE AVAILABILITY (FRACTION) RELATIVE TO POROSITY/SATURATION
1401 ! ----------------------------------------------------------------------
1402 SOILM = -1.0*SMC(1)*ZSOIL(1)
1403 DO K = 2,NSOIL
1404 SOILM = SOILM+SMC(K)*(ZSOIL(K-1)-ZSOIL(K))
1405 END DO
1406 SOILWM = -1.0*(SMCMAX-SMCWLT)*ZSOIL(1)
1407 SOILWW = -1.0*(SMC(1)-SMCWLT)*ZSOIL(1)
1408 DO K = 2,NROOT
1409 SOILWM = SOILWM+(SMCMAX-SMCWLT)*(ZSOIL(K-1)-ZSOIL(K))
1410 SOILWW = SOILWW+(SMC(K)-SMCWLT)*(ZSOIL(K-1)-ZSOIL(K))
1411 END DO
1412 SOILW = SOILWW/SOILWM
1413
1414 ! ----------------------------------------------------------------------
1415 ! END SUBROUTINE SFLX
1416 ! ----------------------------------------------------------------------
1417 END SUBROUTINE SFLX
1418
1419 SUBROUTINE ALCALC (ALB,SNOALB,SHDFAC,SHDMIN,SNCOVR,TSNOW,ALBEDO)
1420
1421 IMPLICIT NONE
1422
1423 ! ----------------------------------------------------------------------
1424 ! CALCULATE ALBEDO INCLUDING SNOW EFFECT (0 -> 1)
1425 ! ALB SNOWFREE ALBEDO
1426 ! SNOALB MAXIMUM (DEEP) SNOW ALBEDO
1427 ! SHDFAC AREAL FRACTIONAL COVERAGE OF GREEN VEGETATION
1428 ! SHDMIN MINIMUM AREAL FRACTIONAL COVERAGE OF GREEN VEGETATION
1429 ! SNCOVR FRACTIONAL SNOW COVER
1430 ! ALBEDO SURFACE ALBEDO INCLUDING SNOW EFFECT
1431 ! TSNOW SNOW SURFACE TEMPERATURE (K)
1432 ! ----------------------------------------------------------------------
1433 REAL ALB, SNOALB, SHDFAC, SHDMIN, SNCOVR, ALBEDO, TSNOW
1434
1435 ! ----------------------------------------------------------------------
1436 ! SNOALB IS ARGUMENT REPRESENTING MAXIMUM ALBEDO OVER DEEP SNOW,
1437 ! AS PASSED INTO SFLX, AND ADAPTED FROM THE SATELLITE-BASED MAXIMUM
1438 ! SNOW ALBEDO FIELDS PROVIDED BY D. ROBINSON AND G. KUKLA
1439 ! (1985, JCAM, VOL 24, 402-411)
1440 ! ----------------------------------------------------------------------
1441 ! changed in version 2.6 on June 2nd 2003
1442 ! ALBEDO = ALB + (1.0-(SHDFAC-SHDMIN))*SNCOVR*(SNOALB-ALB)
1443 ALBEDO = ALB + SNCOVR*(SNOALB-ALB)
1444 IF (ALBEDO .GT. SNOALB) ALBEDO=SNOALB
1445
1446 ! BASE FORMULATION (DICKINSON ET AL., 1986, COGLEY ET AL., 1990)
1447 ! IF (TSNOW.LE.263.16) THEN
1448 ! ALBEDO=SNOALB
1449 ! ELSE
1450 ! IF (TSNOW.LT.273.16) THEN
1451 ! TM=0.1*(TSNOW-263.16)
1452 ! ALBEDO=0.5*((0.9-0.2*(TM**3))+(0.8-0.16*(TM**3)))
1453 ! ELSE
1454 ! ALBEDO=0.67
1455 ! ENDIF
1456 ! ENDIF
1457
1458 ! ISBA FORMULATION (VERSEGHY, 1991; BAKER ET AL., 1990)
1459 ! IF (TSNOW.LT.273.16) THEN
1460 ! ALBEDO=SNOALB-0.008*DT/86400
1461 ! ELSE
1462 ! ALBEDO=(SNOALB-0.5)*EXP(-0.24*DT/86400)+0.5
1463 ! ENDIF
1464
1465 ! ----------------------------------------------------------------------
1466 ! END SUBROUTINE ALCALC
1467 ! ----------------------------------------------------------------------
1468 END SUBROUTINE ALCALC
1469
1470 SUBROUTINE CANRES (SOLAR,CH,SFCTMP,Q2,SFCPRS,SMC,ZSOIL,NSOIL, &
1471 & SMCWLT,SMCREF,RSMIN,RC,PC,NROOT,Q2SAT,DQSDT2, &
1472 & TOPT,RSMAX,RGL,HS,XLAI, &
1473 & RCS,RCT,RCQ,RCSOIL)
1474
1475 IMPLICIT NONE
1476
1477 ! ----------------------------------------------------------------------
1478 ! SUBROUTINE CANRES
1479 ! ----------------------------------------------------------------------
1480 ! CALCULATE CANOPY RESISTANCE WHICH DEPENDS ON INCOMING SOLAR RADIATION,
1481 ! AIR TEMPERATURE, ATMOSPHERIC WATER VAPOR PRESSURE DEFICIT AT THE
1482 ! LOWEST MODEL LEVEL, AND SOIL MOISTURE (PREFERABLY UNFROZEN SOIL
1483 ! MOISTURE RATHER THAN TOTAL)
1484 ! ----------------------------------------------------------------------
1485 ! SOURCE: JARVIS (1976), NOILHAN AND PLANTON (1989, MWR), JACQUEMIN AND
1486 ! NOILHAN (1990, BLM)
1487 ! SEE ALSO: CHEN ET AL (1996, JGR, VOL 101(D3), 7251-7268), EQNS 12-14
1488 ! AND TABLE 2 OF SEC. 3.1.2
1489 ! ----------------------------------------------------------------------
1490 ! INPUT:
1491 ! SOLAR INCOMING SOLAR RADIATION
1492 ! CH SURFACE EXCHANGE COEFFICIENT FOR HEAT AND MOISTURE
1493 ! SFCTMP AIR TEMPERATURE AT 1ST LEVEL ABOVE GROUND
1494 ! Q2 AIR HUMIDITY AT 1ST LEVEL ABOVE GROUND
1495 ! Q2SAT SATURATION AIR HUMIDITY AT 1ST LEVEL ABOVE GROUND
1496 ! DQSDT2 SLOPE OF SATURATION HUMIDITY FUNCTION WRT TEMP
1497 ! SFCPRS SURFACE PRESSURE
1498 ! SMC VOLUMETRIC SOIL MOISTURE
1499 ! ZSOIL SOIL DEPTH (NEGATIVE SIGN, AS IT IS BELOW GROUND)
1500 ! NSOIL NO. OF SOIL LAYERS
1501 ! NROOT NO. OF SOIL LAYERS IN ROOT ZONE (1.LE.NROOT.LE.NSOIL)
1502 ! XLAI LEAF AREA INDEX
1503 ! SMCWLT WILTING POINT
1504 ! SMCREF REFERENCE SOIL MOISTURE (WHERE SOIL WATER DEFICIT STRESS
1505 ! SETS IN)
1506 ! RSMIN, RSMAX, TOPT, RGL, HS ARE CANOPY STRESS PARAMETERS SET IN
1507 ! SURBOUTINE REDPRM
1508 ! OUTPUT:
1509 ! PC PLANT COEFFICIENT
1510 ! RC CANOPY RESISTANCE
1511 ! ----------------------------------------------------------------------
1512 INTEGER NSOLD
1513 PARAMETER(NSOLD = 20)
1514
1515 INTEGER K
1516 INTEGER NROOT
1517 INTEGER NSOIL
1518
1519 REAL CH
1520 REAL CP
1521 REAL DELTA
1522 REAL DQSDT2
1523 REAL FF
1524 REAL GX
1525 REAL HS
1526 REAL P
1527 REAL PART(NSOLD)
1528 REAL PC
1529 REAL Q2
1530 REAL Q2SAT
1531 REAL RC
1532 REAL RSMIN
1533 REAL RCQ
1534 REAL RCS
1535 REAL RCSOIL
1536 REAL RCT
1537 REAL RD
1538 REAL RGL
1539 REAL RR
1540 REAL RSMAX
1541 REAL SFCPRS
1542 REAL SFCTMP
1543 REAL SIGMA
1544 REAL SLV
1545 REAL SMC(NSOIL)
1546 REAL SMCREF
1547 REAL SMCWLT
1548 REAL SOLAR
1549 REAL TOPT
1550 REAL SLVCP
1551 REAL ST1
1552 REAL TAIR4
1553 REAL XLAI
1554 REAL ZSOIL(NSOIL)
1555
1556 PARAMETER(CP = 1004.5)
1557 PARAMETER(RD = 287.04)
1558 PARAMETER(SIGMA = 5.67E-8)
1559 PARAMETER(SLV = 2.501000E6)
1560
1561 ! ----------------------------------------------------------------------
1562 ! INITIALIZE CANOPY RESISTANCE MULTIPLIER TERMS.
1563 ! ----------------------------------------------------------------------
1564 RCS = 0.0
1565 RCT = 0.0
1566 RCQ = 0.0
1567 RCSOIL = 0.0
1568 RC = 0.0
1569
1570 ! ----------------------------------------------------------------------
1571 ! CONTRIBUTION DUE TO INCOMING SOLAR RADIATION
1572 ! ----------------------------------------------------------------------
1573 FF = 0.55*2.0*SOLAR/(RGL*XLAI)
1574 RCS = (FF + RSMIN/RSMAX) / (1.0 + FF)
1575 RCS = MAX(RCS,0.0001)
1576
1577 ! ----------------------------------------------------------------------
1578 ! CONTRIBUTION DUE TO AIR TEMPERATURE AT FIRST MODEL LEVEL ABOVE GROUND
1579 ! RCT EXPRESSION FROM NOILHAN AND PLANTON (1989, MWR).
1580 ! ----------------------------------------------------------------------
1581 RCT = 1.0 - 0.0016*((TOPT-SFCTMP)**2.0)
1582 RCT = MAX(RCT,0.0001)
1583
1584 ! ----------------------------------------------------------------------
1585 ! CONTRIBUTION DUE TO VAPOR PRESSURE DEFICIT AT FIRST MODEL LEVEL.
1586 ! RCQ EXPRESSION FROM SSIB
1587 ! ----------------------------------------------------------------------
1588 RCQ = 1.0/(1.0+HS*(Q2SAT-Q2))
1589 RCQ = MAX(RCQ,0.01)
1590
1591 ! ----------------------------------------------------------------------
1592 ! CONTRIBUTION DUE TO SOIL MOISTURE AVAILABILITY.
1593 ! DETERMINE CONTRIBUTION FROM EACH SOIL LAYER, THEN ADD THEM UP.
1594 ! ----------------------------------------------------------------------
1595 GX = (SMC(1) - SMCWLT) / (SMCREF - SMCWLT)
1596 IF (GX .GT. 1.) GX = 1.
1597 IF (GX .LT. 0.) GX = 0.
1598
1599 ! ----------------------------------------------------------------------
1600 ! USE SOIL DEPTH AS WEIGHTING FACTOR
1601 ! ----------------------------------------------------------------------
1602 PART(1) = (ZSOIL(1)/ZSOIL(NROOT)) * GX
1603 ! ----------------------------------------------------------------------
1604 ! USE ROOT DISTRIBUTION AS WEIGHTING FACTOR
1605 ! PART(1) = RTDIS(1) * GX
1606 ! ----------------------------------------------------------------------
1607 IF (NROOT .GT. 1) THEN
1608 DO K = 2,NROOT
1609 GX = (SMC(K) - SMCWLT) / (SMCREF - SMCWLT)
1610 IF (GX .GT. 1.) GX = 1.
1611 IF (GX .LT. 0.) GX = 0.
1612 ! ----------------------------------------------------------------------
1613 ! USE SOIL DEPTH AS WEIGHTING FACTOR
1614 ! ----------------------------------------------------------------------
1615 PART(K) = ((ZSOIL(K)-ZSOIL(K-1))/ZSOIL(NROOT)) * GX
1616 ! ----------------------------------------------------------------------
1617 ! USE ROOT DISTRIBUTION AS WEIGHTING FACTOR
1618 ! PART(K) = RTDIS(K) * GX
1619 ! ----------------------------------------------------------------------
1620 END DO
1621 ENDIF
1622
1623 DO K = 1,NROOT
1624 RCSOIL = RCSOIL+PART(K)
1625 END DO
1626 RCSOIL = MAX(RCSOIL,0.0001)
1627
1628 ! ----------------------------------------------------------------------
1629 ! DETERMINE CANOPY RESISTANCE DUE TO ALL FACTORS. CONVERT CANOPY
1630 ! RESISTANCE (RC) TO PLANT COEFFICIENT (PC) TO BE USED WITH POTENTIAL
1631 ! EVAP IN DETERMINING ACTUAL EVAP. PC IS DETERMINED BY:
1632 ! PC * LINERIZED PENMAN POTENTIAL EVAP =
1633 ! PENMAN-MONTEITH ACTUAL EVAPORATION (CONTAINING RC TERM).
1634 ! ----------------------------------------------------------------------
1635 RC = RSMIN/(XLAI*RCS*RCT*RCQ*RCSOIL)
1636
1637 ! TAIR4 = SFCTMP**4.
1638 ! ST1 = (4.*SIGMA*RD)/CP
1639 ! SLVCP = SLV/CP
1640 ! RR = ST1*TAIR4/(SFCPRS*CH) + 1.0
1641 RR = (4.*SIGMA*RD/CP)*(SFCTMP**4.)/(SFCPRS*CH) + 1.0
1642 DELTA = (SLV/CP)*DQSDT2
1643
1644 PC = (RR+DELTA)/(RR*(1.+RC*CH)+DELTA)
1645
1646 ! ----------------------------------------------------------------------
1647 ! END SUBROUTINE CANRES
1648 ! ----------------------------------------------------------------------
1649 END SUBROUTINE CANRES
1650
1651 SUBROUTINE DEVAP (EDIR1,ETP1,SMC,ZSOIL,SHDFAC,SMCMAX,BEXP, &
1652 & DKSAT,DWSAT,SMCDRY,SMCREF,SMCWLT,FXEXP)
1653
1654 IMPLICIT NONE
1655
1656 ! ----------------------------------------------------------------------
1657 ! SUBROUTINE DEVAP
1658 ! ----------------------------------------------------------------------
1659 ! CALCULATE DIRECT SOIL EVAPORATION
1660 ! ----------------------------------------------------------------------
1661 REAL BEXP
1662 ! REAL DEVAP
1663 REAL EDIR1
1664 REAL DKSAT
1665 REAL DWSAT
1666 REAL ETP1
1667 REAL FX
1668 REAL FXEXP
1669 REAL SHDFAC
1670 REAL SMC
1671 REAL SMCDRY
1672 REAL SMCMAX
1673 REAL ZSOIL
1674 REAL SMCREF
1675 REAL SMCWLT
1676 REAL SRATIO
1677
1678 ! ----------------------------------------------------------------------
1679 ! DIRECT EVAP A FUNCTION OF RELATIVE SOIL MOISTURE AVAILABILITY, LINEAR
1680 ! WHEN FXEXP=1.
1681 ! FX > 1 REPRESENTS DEMAND CONTROL
1682 ! FX < 1 REPRESENTS FLUX CONTROL
1683 ! ----------------------------------------------------------------------
1684 SRATIO = (SMC - SMCDRY) / (SMCMAX - SMCDRY)
1685 IF (SRATIO .GT. 0.) THEN
1686 FX = SRATIO**FXEXP
1687 FX = MAX ( MIN ( FX, 1. ) ,0. )
1688 ELSE
1689 FX = 0.
1690 ENDIF
1691
1692 ! ----------------------------------------------------------------------
1693 ! ALLOW FOR THE DIRECT-EVAP-REDUCING EFFECT OF SHADE
1694 ! ----------------------------------------------------------------------
1695 ! DEVAP = FX * ( 1.0 - SHDFAC ) * ETP1
1696 EDIR1 = FX * ( 1.0 - SHDFAC ) * ETP1
1697
1698 ! ----------------------------------------------------------------------
1699 ! END SUBROUTINE DEVAP
1700 ! ----------------------------------------------------------------------
1701 END SUBROUTINE DEVAP
1702
1703 SUBROUTINE EVAPO (ETA1,SMC,NSOIL,CMC,ETP1,DT,ZSOIL, &
1704 & SH2O, &
1705 & SMCMAX,BEXP,PC,SMCWLT,DKSAT,DWSAT, &
1706 & SMCREF,SHDFAC,CMCMAX, &
1707 & SMCDRY,CFACTR, &
1708 & EDIR1,EC1,ET1,ETT1,SFCTMP,Q2,NROOT,RTDIS,FXEXP)
1709
1710 IMPLICIT NONE
1711
1712 ! ----------------------------------------------------------------------
1713 ! SUBROUTINE EVAPO
1714 ! ----------------------------------------------------------------------
1715 ! CALCULATE SOIL MOISTURE FLUX. THE SOIL MOISTURE CONTENT (SMC - A PER
1716 ! UNIT VOLUME MEASUREMENT) IS A DEPENDENT VARIABLE THAT IS UPDATED WITH
1717 ! PROGNOSTIC EQNS. THE CANOPY MOISTURE CONTENT (CMC) IS ALSO UPDATED.
1718 ! FROZEN GROUND VERSION: NEW STATES ADDED: SH2O, AND FROZEN GROUND
1719 ! CORRECTION FACTOR, FRZFACT AND PARAMETER SLOPE.
1720 ! ----------------------------------------------------------------------
1721 INTEGER NSOLD
1722 PARAMETER(NSOLD = 20)
1723
1724 INTEGER I
1725 INTEGER K
1726 INTEGER NSOIL
1727 INTEGER NROOT
1728
1729 REAL BEXP
1730 REAL CFACTR
1731 REAL CMC
1732 REAL CMC2MS
1733 REAL CMCMAX
1734 ! REAL DEVAP
1735 REAL DKSAT
1736 REAL DT
1737 REAL DWSAT
1738 REAL EC1
1739 REAL EDIR1
1740 REAL ET1(NSOIL)
1741 REAL ETA1
1742 REAL ETP1
1743 REAL ETT1
1744 REAL FXEXP
1745 REAL PC
1746 REAL Q2
1747 REAL RTDIS(NSOIL)
1748 REAL SFCTMP
1749 REAL SHDFAC
1750 REAL SMC(NSOIL)
1751 REAL SH2O(NSOIL)
1752 REAL SMCDRY
1753 REAL SMCMAX
1754 REAL SMCREF
1755 REAL SMCWLT
1756 REAL ZSOIL(NSOIL)
1757
1758 ! ----------------------------------------------------------------------
1759 ! EXECUTABLE CODE BEGINS HERE IF THE POTENTIAL EVAPOTRANSPIRATION IS
1760 ! GREATER THAN ZERO.
1761 ! ----------------------------------------------------------------------
1762 EDIR1 = 0.
1763 EC1 = 0.
1764 DO K = 1,NSOIL
1765 ET1(K) = 0.
1766 END DO
1767 ETT1 = 0.
1768
1769 IF (ETP1 .GT. 0.0) THEN
1770
1771 ! ----------------------------------------------------------------------
1772 ! RETRIEVE DIRECT EVAPORATION FROM SOIL SURFACE. CALL THIS FUNCTION
1773 ! ONLY IF VEG COVER NOT COMPLETE.
1774 ! FROZEN GROUND VERSION: SH2O STATES REPLACE SMC STATES.
1775 ! ----------------------------------------------------------------------
1776 IF (SHDFAC .LT. 1.) THEN
1777 CALL DEVAP (EDIR1,ETP1,SH2O(1),ZSOIL(1),SHDFAC,SMCMAX, &
1778 ! EDIR = DEVAP(ETP1,SH2O(1),ZSOIL(1),SHDFAC,SMCMAX, &
1779 & BEXP,DKSAT,DWSAT,SMCDRY,SMCREF,SMCWLT,FXEXP)
1780 ENDIF
1781
1782 ! ----------------------------------------------------------------------
1783 ! INITIALIZE PLANT TOTAL TRANSPIRATION, RETRIEVE PLANT TRANSPIRATION,
1784 ! AND ACCUMULATE IT FOR ALL SOIL LAYERS.
1785 ! ----------------------------------------------------------------------
1786 IF (SHDFAC.GT.0.0) THEN
1787
1788 CALL TRANSP (ET1,NSOIL,ETP1,SH2O,CMC,ZSOIL,SHDFAC,SMCWLT, &
1789 & CMCMAX,PC,CFACTR,SMCREF,SFCTMP,Q2,NROOT,RTDIS)
1790
1791 DO K = 1,NSOIL
1792 ETT1 = ETT1 + ET1(K)
1793 END DO
1794
1795 ! ----------------------------------------------------------------------
1796 ! CALCULATE CANOPY EVAPORATION.
1797 ! IF STATEMENTS TO AVOID TANGENT LINEAR PROBLEMS NEAR CMC=0.0.
1798 ! ----------------------------------------------------------------------
1799 IF (CMC .GT. 0.0) THEN
1800 EC1 = SHDFAC * ( ( CMC / CMCMAX ) ** CFACTR ) * ETP1
1801 ELSE
1802 EC1 = 0.0
1803 ENDIF
1804
1805 ! ----------------------------------------------------------------------
1806 ! EC SHOULD BE LIMITED BY THE TOTAL AMOUNT OF AVAILABLE WATER ON THE
1807 ! CANOPY. -F.CHEN, 18-OCT-1994
1808 ! ----------------------------------------------------------------------
1809 CMC2MS = CMC / DT
1810 EC1 = MIN ( CMC2MS, EC1 )
1811 ENDIF
1812 ENDIF
1813
1814 ! ----------------------------------------------------------------------
1815 ! TOTAL UP EVAP AND TRANSP TYPES TO OBTAIN ACTUAL EVAPOTRANSP
1816 ! ----------------------------------------------------------------------
1817 ETA1 = EDIR1 + ETT1 + EC1
1818
1819 ! ----------------------------------------------------------------------
1820 ! END SUBROUTINE EVAPO
1821 ! ----------------------------------------------------------------------
1822 END SUBROUTINE EVAPO
1823
1824 SUBROUTINE HRT (RHSTS,STC,SMC,SMCMAX,NSOIL,ZSOIL,YY,ZZ1, &
1825 & TBOT,ZBOT,PSISAT,SH2O,DT,BEXP, &
1826 & F1,DF1,QUARTZ,CSOIL,AI,BI,CI)
1827
1828 IMPLICIT NONE
1829
1830 ! ----------------------------------------------------------------------
1831 ! SUBROUTINE HRT
1832 ! ----------------------------------------------------------------------
1833 ! CALCULATE THE RIGHT HAND SIDE OF THE TIME TENDENCY TERM OF THE SOIL
1834 ! THERMAL DIFFUSION EQUATION. ALSO TO COMPUTE ( PREPARE ) THE MATRIX
1835 ! COEFFICIENTS FOR THE TRI-DIAGONAL MATRIX OF THE IMPLICIT TIME SCHEME.
1836 ! ----------------------------------------------------------------------
1837 INTEGER NSOLD
1838 PARAMETER(NSOLD = 20)
1839
1840 LOGICAL ITAVG
1841
1842 INTEGER I
1843 INTEGER K
1844 INTEGER NSOIL
1845
1846 ! ----------------------------------------------------------------------
1847 ! DECLARE WORK ARRAYS NEEDED IN TRI-DIAGONAL IMPLICIT SOLVER
1848 ! ----------------------------------------------------------------------
1849 REAL AI(NSOLD)
1850 REAL BI(NSOLD)
1851 REAL CI(NSOLD)
1852
1853 ! ----------------------------------------------------------------------
1854 ! DECLARATIONS
1855 ! ----------------------------------------------------------------------
1856 REAL BEXP
1857 REAL CAIR
1858 REAL CH2O
1859 REAL CICE
1860 REAL CSOIL
1861 REAL DDZ
1862 REAL DDZ2
1863 REAL DENOM
1864 REAL DF1
1865 REAL DF1N
1866 REAL DF1K
1867 REAL DT
1868 REAL DTSDZ
1869 REAL DTSDZ2
1870 REAL F1
1871 REAL HCPCT
1872 REAL PSISAT
1873 REAL QUARTZ
1874 REAL QTOT
1875 REAL RHSTS(NSOIL)
1876 REAL SSOIL
1877 REAL SICE
1878 REAL SMC(NSOIL)
1879 REAL SH2O(NSOIL)
1880 REAL SMCMAX
1881 ! REAL SNKSRC
1882 REAL STC(NSOIL)
1883 REAL T0
1884 REAL TAVG
1885 REAL TBK
1886 REAL TBK1
1887 REAL TBOT
1888 REAL ZBOT
1889 REAL TSNSR
1890 REAL TSURF
1891 REAL YY
1892 REAL ZSOIL(NSOIL)
1893 REAL ZZ1
1894
1895 PARAMETER(T0 = 273.15)
1896
1897 ! ----------------------------------------------------------------------
1898 ! SET SPECIFIC HEAT CAPACITIES OF AIR, WATER, ICE, SOIL MINERAL
1899 ! ----------------------------------------------------------------------
1900 PARAMETER(CAIR = 1004.0)
1901 PARAMETER(CH2O = 4.2E6)
1902 PARAMETER(CICE = 2.106E6)
1903 ! NOTE: CSOIL NOW SET IN ROUTINE REDPRM AND PASSED IN
1904 ! PARAMETER(CSOIL = 1.26E6)
1905
1906 ! ----------------------------------------------------------------------
1907 ! INITIALIZE LOGICAL FOR SOIL LAYER TEMPERATURE AVERAGING.
1908 ! ----------------------------------------------------------------------
1909 ITAVG = .TRUE.
1910 ! ITAVG = .FALSE.
1911
1912 ! ----------------------------------------------------------------------
1913 ! BEGIN SECTION FOR TOP SOIL LAYER
1914 ! ----------------------------------------------------------------------
1915 ! CALC THE HEAT CAPACITY OF THE TOP SOIL LAYER
1916 ! ----------------------------------------------------------------------
1917 HCPCT = SH2O(1)*CH2O + (1.0-SMCMAX)*CSOIL + (SMCMAX-SMC(1))*CAIR &
1918 & + ( SMC(1) - SH2O(1) )*CICE
1919
1920 ! ----------------------------------------------------------------------
1921 ! CALC THE MATRIX COEFFICIENTS AI, BI, AND CI FOR THE TOP LAYER
1922 ! ----------------------------------------------------------------------
1923 DDZ = 1.0 / ( -0.5 * ZSOIL(2) )
1924 AI(1) = 0.0
1925 CI(1) = (DF1 * DDZ) / (ZSOIL(1) * HCPCT)
1926 BI(1) = -CI(1) + DF1 / (0.5 * ZSOIL(1) * ZSOIL(1)*HCPCT*ZZ1)
1927
1928 ! ----------------------------------------------------------------------
1929 ! CALCULATE THE VERTICAL SOIL TEMP GRADIENT BTWN THE 1ST AND 2ND SOIL
1930 ! LAYERS. THEN CALCULATE THE SUBSURFACE HEAT FLUX. USE THE TEMP
1931 ! GRADIENT AND SUBSFC HEAT FLUX TO CALC "RIGHT-HAND SIDE TENDENCY
1932 ! TERMS", OR "RHSTS", FOR TOP SOIL LAYER.
1933 ! ----------------------------------------------------------------------
1934 DTSDZ = (STC(1) - STC(2)) / (-0.5 * ZSOIL(2))
1935 SSOIL = DF1 * (STC(1) - YY) / (0.5 * ZSOIL(1) * ZZ1)
1936 RHSTS(1) = (DF1 * DTSDZ - SSOIL) / (ZSOIL(1) * HCPCT)
1937
1938 ! ----------------------------------------------------------------------
1939 ! NEXT CAPTURE THE VERTICAL DIFFERENCE OF THE HEAT FLUX AT TOP AND
1940 ! BOTTOM OF FIRST SOIL LAYER FOR USE IN HEAT FLUX CONSTRAINT APPLIED TO
1941 ! POTENTIAL SOIL FREEZING/THAWING IN ROUTINE SNKSRC.
1942 ! ----------------------------------------------------------------------
1943 QTOT = SSOIL - DF1*DTSDZ
1944
1945 ! ----------------------------------------------------------------------
1946 ! IF TEMPERATURE AVERAGING INVOKED (ITAVG=TRUE; ELSE SKIP):
1947 ! SET TEMP "TSURF" AT TOP OF SOIL COLUMN (FOR USE IN FREEZING SOIL
1948 ! PHYSICS LATER IN FUNCTION SUBROUTINE SNKSRC). IF SNOWPACK CONTENT IS
1949 ! ZERO, THEN TSURF EXPRESSION BELOW GIVES TSURF = SKIN TEMP. IF
1950 ! SNOWPACK IS NONZERO (HENCE ARGUMENT ZZ1=1), THEN TSURF EXPRESSION
1951 ! BELOW YIELDS SOIL COLUMN TOP TEMPERATURE UNDER SNOWPACK. THEN
1952 ! CALCULATE TEMPERATURE AT BOTTOM INTERFACE OF 1ST SOIL LAYER FOR USE
1953 ! LATER IN FUNCTION SUBROUTINE SNKSRC
1954 ! ----------------------------------------------------------------------
1955 IF (ITAVG) THEN
1956 TSURF = (YY + (ZZ1-1) * STC(1)) / ZZ1
1957 CALL TBND (STC(1),STC(2),ZSOIL,ZBOT,1,NSOIL,TBK)
1958 ENDIF
1959
1960 ! ----------------------------------------------------------------------
1961 ! CALCULATE FROZEN WATER CONTENT IN 1ST SOIL LAYER.
1962 ! ----------------------------------------------------------------------
1963 SICE = SMC(1) - SH2O(1)
1964
1965 ! ----------------------------------------------------------------------
1966 ! IF FROZEN WATER PRESENT OR ANY OF LAYER-1 MID-POINT OR BOUNDING
1967 ! INTERFACE TEMPERATURES BELOW FREEZING, THEN CALL SNKSRC TO
1968 ! COMPUTE HEAT SOURCE/SINK (AND CHANGE IN FROZEN WATER CONTENT)
1969 ! DUE TO POSSIBLE SOIL WATER PHASE CHANGE
1970 ! ----------------------------------------------------------------------
1971 IF ( (SICE .GT. 0.) .OR. (TSURF .LT. T0) .OR. &
1972 & (STC(1) .LT. T0) .OR. (TBK .LT. T0) ) THEN
1973
1974 IF (ITAVG) THEN
1975 CALL TMPAVG(TAVG,TSURF,STC(1),TBK,ZSOIL,NSOIL,1)
1976 ELSE
1977 TAVG = STC(1)
1978 ENDIF
1979 TSNSR = SNKSRC (TAVG,SMC(1),SH2O(1), &
1980 & ZSOIL,NSOIL,SMCMAX,PSISAT,BEXP,DT,1,QTOT)
1981
1982 RHSTS(1) = RHSTS(1) - TSNSR / ( ZSOIL(1) * HCPCT )
1983 ENDIF
1984
1985 ! ----------------------------------------------------------------------
1986 ! THIS ENDS SECTION FOR TOP SOIL LAYER.
1987 ! ----------------------------------------------------------------------
1988 ! INITIALIZE DDZ2
1989 ! ----------------------------------------------------------------------
1990 DDZ2 = 0.0
1991
1992 ! ----------------------------------------------------------------------
1993 ! LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABOVE PROCESS
1994 ! (EXCEPT SUBSFC OR "GROUND" HEAT FLUX NOT REPEATED IN LOWER LAYERS)
1995 ! ----------------------------------------------------------------------
1996 DF1K = DF1
1997 DO K = 2,NSOIL
1998
1999 ! ----------------------------------------------------------------------
2000 ! CALCULATE HEAT CAPACITY FOR THIS SOIL LAYER.
2001 ! ----------------------------------------------------------------------
2002 HCPCT = SH2O(K)*CH2O +(1.0-SMCMAX)*CSOIL +(SMCMAX-SMC(K))*CAIR &
2003 & + ( SMC(K) - SH2O(K) )*CICE
2004
2005 IF (K .NE. NSOIL) THEN
2006 ! ----------------------------------------------------------------------
2007 ! THIS SECTION FOR LAYER 2 OR GREATER, BUT NOT LAST LAYER.
2008 ! ----------------------------------------------------------------------
2009 ! CALCULATE THERMAL DIFFUSIVITY FOR THIS LAYER.
2010 ! ----------------------------------------------------------------------
2011 CALL TDFCND (DF1N,SMC(K),QUARTZ,SMCMAX,SH2O(K))
2012
2013 ! ----------------------------------------------------------------------
2014 ! CALC THE VERTICAL SOIL TEMP GRADIENT THRU THIS LAYER
2015 ! ----------------------------------------------------------------------
2016 DENOM = 0.5 * ( ZSOIL(K-1) - ZSOIL(K+1) )
2017 DTSDZ2 = ( STC(K) - STC(K+1) ) / DENOM
2018
2019 ! ----------------------------------------------------------------------
2020 ! CALC THE MATRIX COEF, CI, AFTER CALC'NG ITS PARTIAL PRODUCT
2021 ! ----------------------------------------------------------------------
2022 DDZ2 = 2. / (ZSOIL(K-1) - ZSOIL(K+1))
2023 CI(K) = -DF1N * DDZ2 / ((ZSOIL(K-1) - ZSOIL(K)) * HCPCT)
2024
2025 ! ----------------------------------------------------------------------
2026 ! IF TEMPERATURE AVERAGING INVOKED (ITAVG=TRUE; ELSE SKIP): CALCULATE
2027 ! TEMP AT BOTTOM OF LAYER.
2028 ! ----------------------------------------------------------------------
2029 IF (ITAVG) THEN
2030 CALL TBND (STC(K),STC(K+1),ZSOIL,ZBOT,K,NSOIL,TBK1)
2031 ENDIF
2032 ELSE
2033
2034 ! ----------------------------------------------------------------------
2035 ! SPECIAL CASE OF BOTTOM SOIL LAYER: CALCULATE THERMAL DIFFUSIVITY FOR
2036 ! BOTTOM LAYER.
2037 ! ----------------------------------------------------------------------
2038 CALL TDFCND (DF1N,SMC(K),QUARTZ,SMCMAX,SH2O(K))
2039
2040 ! ----------------------------------------------------------------------
2041 ! CALC THE VERTICAL SOIL TEMP GRADIENT THRU BOTTOM LAYER.
2042 ! ----------------------------------------------------------------------
2043 DENOM = .5 * (ZSOIL(K-1) + ZSOIL(K)) - ZBOT
2044 DTSDZ2 = (STC(K)-TBOT) / DENOM
2045
2046 ! ----------------------------------------------------------------------
2047 ! SET MATRIX COEF, CI TO ZERO IF BOTTOM LAYER.
2048 ! ----------------------------------------------------------------------
2049 CI(K) = 0.
2050
2051 ! ----------------------------------------------------------------------
2052 ! IF TEMPERATURE AVERAGING INVOKED (ITAVG=TRUE; ELSE SKIP): CALCULATE
2053 ! TEMP AT BOTTOM OF LAST LAYER.
2054 ! ----------------------------------------------------------------------
2055 IF (ITAVG) THEN
2056 CALL TBND (STC(K),TBOT,ZSOIL,ZBOT,K,NSOIL,TBK1)
2057 ENDIF
2058
2059 ENDIF
2060 ! ----------------------------------------------------------------------
2061 ! THIS ENDS SPECIAL LOOP FOR BOTTOM LAYER.
2062 ! ----------------------------------------------------------------------
2063 ! CALCULATE RHSTS FOR THIS LAYER AFTER CALC'NG A PARTIAL PRODUCT.
2064 ! ----------------------------------------------------------------------
2065 DENOM = ( ZSOIL(K) - ZSOIL(K-1) ) * HCPCT
2066 RHSTS(K) = ( DF1N * DTSDZ2 - DF1K * DTSDZ ) / DENOM
2067 QTOT = -1.0*DENOM*RHSTS(K)
2068 SICE = SMC(K) - SH2O(K)
2069
2070 IF ( (SICE .GT. 0.) .OR. (TBK .LT. T0) .OR. &
2071 & (STC(K) .LT. T0) .OR. (TBK1 .LT. T0) ) THEN
2072
2073 IF (ITAVG) THEN
2074 CALL TMPAVG(TAVG,TBK,STC(K),TBK1,ZSOIL,NSOIL,K)
2075 ELSE
2076 TAVG = STC(K)
2077 ENDIF
2078 TSNSR = SNKSRC(TAVG,SMC(K),SH2O(K),ZSOIL,NSOIL, &
2079 & SMCMAX,PSISAT,BEXP,DT,K,QTOT)
2080 RHSTS(K) = RHSTS(K) - TSNSR / DENOM
2081 ENDIF
2082
2083 ! ----------------------------------------------------------------------
2084 ! CALC MATRIX COEFS, AI, AND BI FOR THIS LAYER.
2085 ! ----------------------------------------------------------------------
2086 AI(K) = - DF1 * DDZ / ((ZSOIL(K-1) - ZSOIL(K)) * HCPCT)
2087 BI(K) = -(AI(K) + CI(K))
2088
2089 ! ----------------------------------------------------------------------
2090 ! RESET VALUES OF DF1, DTSDZ, DDZ, AND TBK FOR LOOP TO NEXT SOIL LAYER.
2091 ! ----------------------------------------------------------------------
2092 TBK = TBK1
2093 DF1K = DF1N
2094 DTSDZ = DTSDZ2
2095 DDZ = DDZ2
2096 END DO
2097
2098 ! ----------------------------------------------------------------------
2099 ! END SUBROUTINE HRT
2100 ! ----------------------------------------------------------------------
2101 END SUBROUTINE HRT
2102
2103 SUBROUTINE HRTICE (RHSTS,STC,NSOIL,ZSOIL,YY,ZZ1,DF1,AI,BI,CI)
2104
2105 IMPLICIT NONE
2106
2107 ! ----------------------------------------------------------------------
2108 ! SUBROUTINE HRTICE
2109 ! ----------------------------------------------------------------------
2110 ! CALCULATE THE RIGHT HAND SIDE OF THE TIME TENDENCY TERM OF THE SOIL
2111 ! THERMAL DIFFUSION EQUATION IN THE CASE OF SEA-ICE PACK. ALSO TO
2112 ! COMPUTE (PREPARE) THE MATRIX COEFFICIENTS FOR THE TRI-DIAGONAL MATRIX
2113 ! OF THE IMPLICIT TIME SCHEME.
2114 ! ----------------------------------------------------------------------
2115 INTEGER NSOLD
2116 PARAMETER(NSOLD = 20)
2117
2118 INTEGER K
2119 INTEGER NSOIL
2120
2121 REAL AI(NSOLD)
2122 REAL BI(NSOLD)
2123 REAL CI(NSOLD)
2124
2125 REAL DDZ
2126 REAL DDZ2
2127 REAL DENOM
2128 REAL DF1
2129 REAL DTSDZ
2130 REAL DTSDZ2
2131 REAL HCPCT
2132 REAL RHSTS(NSOIL)
2133 REAL SSOIL
2134 REAL STC(NSOIL)
2135 REAL TBOT
2136 REAL YY
2137 REAL ZBOT
2138 REAL ZSOIL(NSOIL)
2139 REAL ZZ1
2140
2141 DATA TBOT /271.16/
2142
2143 ! ----------------------------------------------------------------------
2144 ! SET A NOMINAL UNIVERSAL VALUE OF THE SEA-ICE SPECIFIC HEAT CAPACITY,
2145 ! HCPCT = 1880.0*917.0.
2146 ! ----------------------------------------------------------------------
2147 PARAMETER(HCPCT = 1.72396E+6)
2148
2149 ! ----------------------------------------------------------------------
2150 ! THE INPUT ARGUMENT DF1 IS A UNIVERSALLY CONSTANT VALUE OF SEA-ICE
2151 ! THERMAL DIFFUSIVITY, SET IN ROUTINE SNOPAC AS DF1 = 2.2.
2152 ! ----------------------------------------------------------------------
2153 ! SET ICE PACK DEPTH. USE TBOT AS ICE PACK LOWER BOUNDARY TEMPERATURE
2154 ! (THAT OF UNFROZEN SEA WATER AT BOTTOM OF SEA ICE PACK). ASSUME ICE
2155 ! PACK IS OF N=NSOIL LAYERS SPANNING A UNIFORM CONSTANT ICE PACK
2156 ! THICKNESS AS DEFINED BY ZSOIL(NSOIL) IN ROUTINE SFLX.
2157 ! ----------------------------------------------------------------------
2158 ZBOT = ZSOIL(NSOIL)
2159
2160 ! ----------------------------------------------------------------------
2161 ! CALC THE MATRIX COEFFICIENTS AI, BI, AND CI FOR THE TOP LAYER
2162 ! ----------------------------------------------------------------------
2163 DDZ = 1.0 / ( -0.5 * ZSOIL(2) )
2164 AI(1) = 0.0
2165 CI(1) = (DF1 * DDZ) / (ZSOIL(1) * HCPCT)
2166 BI(1) = -CI(1) + DF1/(0.5 * ZSOIL(1) * ZSOIL(1) * HCPCT * ZZ1)
2167
2168 ! ----------------------------------------------------------------------
2169 ! CALC THE VERTICAL SOIL TEMP GRADIENT BTWN THE TOP AND 2ND SOIL LAYERS.
2170 ! RECALC/ADJUST THE SOIL HEAT FLUX. USE THE GRADIENT AND FLUX TO CALC
2171 ! RHSTS FOR THE TOP SOIL LAYER.
2172 ! ----------------------------------------------------------------------
2173 DTSDZ = ( STC(1) - STC(2) ) / ( -0.5 * ZSOIL(2) )
2174 SSOIL = DF1 * ( STC(1) - YY ) / ( 0.5 * ZSOIL(1) * ZZ1 )
2175 RHSTS(1) = ( DF1 * DTSDZ - SSOIL ) / ( ZSOIL(1) * HCPCT )
2176
2177 ! ----------------------------------------------------------------------
2178 ! INITIALIZE DDZ2
2179 ! ----------------------------------------------------------------------
2180 DDZ2 = 0.0
2181
2182 ! ----------------------------------------------------------------------
2183 ! LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABOVE PROCESS
2184 ! ----------------------------------------------------------------------
2185 DO K = 2,NSOIL
2186 IF (K .NE. NSOIL) THEN
2187
2188 ! ----------------------------------------------------------------------
2189 ! CALC THE VERTICAL SOIL TEMP GRADIENT THRU THIS LAYER.
2190 ! ----------------------------------------------------------------------
2191 DENOM = 0.5 * ( ZSOIL(K-1) - ZSOIL(K+1) )
2192 DTSDZ2 = ( STC(K) - STC(K+1) ) / DENOM
2193
2194 ! ----------------------------------------------------------------------
2195 ! CALC THE MATRIX COEF, CI, AFTER CALC'NG ITS PARTIAL PRODUCT.
2196 ! ----------------------------------------------------------------------
2197 DDZ2 = 2. / (ZSOIL(K-1) - ZSOIL(K+1))
2198 CI(K) = -DF1 * DDZ2 / ((ZSOIL(K-1) - ZSOIL(K)) * HCPCT)
2199 ELSE
2200
2201 ! ----------------------------------------------------------------------
2202 ! CALC THE VERTICAL SOIL TEMP GRADIENT THRU THE LOWEST LAYER.
2203 ! ----------------------------------------------------------------------
2204 DTSDZ2 = (STC(K)-TBOT)/(.5 * (ZSOIL(K-1) + ZSOIL(K))-ZBOT)
2205
2206 ! ----------------------------------------------------------------------
2207 ! SET MATRIX COEF, CI TO ZERO.
2208 ! ----------------------------------------------------------------------
2209 CI(K) = 0.
2210 ENDIF
2211
2212 ! ----------------------------------------------------------------------
2213 ! CALC RHSTS FOR THIS LAYER AFTER CALC'NG A PARTIAL PRODUCT.
2214 ! ----------------------------------------------------------------------
2215 DENOM = ( ZSOIL(K) - ZSOIL(K-1) ) * HCPCT
2216 RHSTS(K) = ( DF1 * DTSDZ2 - DF1 * DTSDZ ) / DENOM
2217
2218 ! ----------------------------------------------------------------------
2219 ! CALC MATRIX COEFS, AI, AND BI FOR THIS LAYER.
2220 ! ----------------------------------------------------------------------
2221 AI(K) = - DF1 * DDZ / ((ZSOIL(K-1) - ZSOIL(K)) * HCPCT)
2222 BI(K) = -(AI(K) + CI(K))
2223
2224 ! ----------------------------------------------------------------------
2225 ! RESET VALUES OF DTSDZ AND DDZ FOR LOOP TO NEXT SOIL LYR.
2226 ! ----------------------------------------------------------------------
2227 DTSDZ = DTSDZ2
2228 DDZ = DDZ2
2229
2230 END DO
2231 ! ----------------------------------------------------------------------
2232 ! END SUBROUTINE HRTICE
2233 ! ----------------------------------------------------------------------
2234 END SUBROUTINE HRTICE
2235
2236 SUBROUTINE HSTEP (STCOUT,STCIN,RHSTS,DT,NSOIL,AI,BI,CI)
2237
2238 IMPLICIT NONE
2239
2240 ! ----------------------------------------------------------------------
2241 ! SUBROUTINE HSTEP
2242 ! ----------------------------------------------------------------------
2243 ! CALCULATE/UPDATE THE SOIL TEMPERATURE FIELD.
2244 ! ----------------------------------------------------------------------
2245 INTEGER NSOLD
2246 PARAMETER(NSOLD = 20)
2247
2248 INTEGER K
2249 INTEGER NSOIL
2250
2251 REAL AI(NSOLD)
2252 REAL BI(NSOLD)
2253 REAL CI(NSOLD)
2254 REAL CIin(NSOLD)
2255 REAL DT
2256 REAL RHSTS(NSOIL)
2257 REAL RHSTSin(NSOIL)
2258 REAL STCIN(NSOIL)
2259 REAL STCOUT(NSOIL)
2260
2261 ! ----------------------------------------------------------------------
2262 ! CREATE FINITE DIFFERENCE VALUES FOR USE IN ROSR12 ROUTINE
2263 ! ----------------------------------------------------------------------
2264 DO K = 1,NSOIL
2265 RHSTS(K) = RHSTS(K) * DT
2266 AI(K) = AI(K) * DT
2267 BI(K) = 1. + BI(K) * DT
2268 CI(K) = CI(K) * DT
2269 END DO
2270
2271 ! ----------------------------------------------------------------------
2272 ! COPY VALUES FOR INPUT VARIABLES BEFORE CALL TO ROSR12
2273 ! ----------------------------------------------------------------------
2274 DO K = 1,NSOIL
2275 RHSTSin(K) = RHSTS(K)
2276 END DO
2277 DO K = 1,NSOIL
2278 CIin(K) = CI(K)
2279 END DO
2280
2281 ! ----------------------------------------------------------------------
2282 ! SOLVE THE TRI-DIAGONAL MATRIX EQUATION
2283 ! ----------------------------------------------------------------------
2284 CALL ROSR12(CI,AI,BI,CIin,RHSTSin,RHSTS,NSOIL)
2285
2286 ! ----------------------------------------------------------------------
2287 ! CALC/UPDATE THE SOIL TEMPS USING MATRIX SOLUTION
2288 ! ----------------------------------------------------------------------
2289 DO K = 1,NSOIL
2290 STCOUT(K) = STCIN(K) + CI(K)
2291 END DO
2292
2293 ! ----------------------------------------------------------------------
2294 ! END SUBROUTINE HSTEP
2295 ! ----------------------------------------------------------------------
2296 END SUBROUTINE HSTEP
2297
2298 SUBROUTINE NOPAC(ETP,ETA,PRCP,SMC,SMCMAX,SMCWLT, &
2299 & SMCREF,SMCDRY,CMC,CMCMAX,NSOIL,DT,SHDFAC, &
2300 & SBETA,Q2,T1,SFCTMP,T24,TH2,FDOWN,F1,SSOIL, &
2301 & STC,EPSCA,BEXP,PC,RCH,RR,CFACTR, &
2302 & SH2O,SLOPE,KDT,FRZFACT,PSISAT,ZSOIL, &
2303 & DKSAT,DWSAT,TBOT,ZBOT,RUNOFF1,RUNOFF2, &
2304 & RUNOFF3,EDIR,EC,ET,ETT,NROOT,ICE,RTDIS, &
2305 & QUARTZ,FXEXP,CSOIL, &
2306 & BETA,DRIP,DEW,FLX1,FLX2,FLX3)
2307
2308 IMPLICIT NONE
2309
2310 ! ----------------------------------------------------------------------
2311 ! SUBROUTINE NOPAC
2312 ! ----------------------------------------------------------------------
2313 ! CALCULATE SOIL MOISTURE AND HEAT FLUX VALUES AND UPDATE SOIL MOISTURE
2314 ! CONTENT AND SOIL HEAT CONTENT VALUES FOR THE CASE WHEN NO SNOW PACK IS
2315 ! PRESENT.
2316 ! ----------------------------------------------------------------------
2317 INTEGER ICE
2318 INTEGER NROOT
2319 INTEGER NSOIL
2320
2321 REAL BEXP
2322 REAL BETA
2323 REAL CFACTR
2324 REAL CMC
2325 REAL CMCMAX
2326 REAL CP
2327 REAL CSOIL
2328 REAL DEW
2329 REAL DF1
2330 REAL DKSAT
2331 REAL DRIP
2332 REAL DT
2333 REAL DWSAT
2334 REAL EC
2335 REAL EDIR
2336 REAL EPSCA
2337 REAL ETA
2338 REAL ETA1
2339 REAL ETP
2340 REAL ETP1
2341 REAL ET(NSOIL)
2342 REAL ETT
2343 REAL FDOWN
2344 REAL F1
2345 REAL FXEXP
2346 REAL FLX1
2347 REAL FLX2
2348 REAL FLX3
2349 REAL FRZFACT
2350 REAL KDT
2351 REAL PC
2352 REAL PRCP
2353 REAL PRCP1
2354 REAL PSISAT
2355 REAL Q2
2356 REAL QUARTZ
2357 REAL RCH
2358 REAL RR
2359 REAL RTDIS(NSOIL)
2360 REAL RUNOFF1
2361 REAL RUNOFF2
2362 REAL RUNOFF3
2363 REAL SSOIL
2364 REAL SBETA
2365 REAL SFCTMP
2366 REAL SHDFAC
2367 REAL SH2O(NSOIL)
2368 REAL SIGMA
2369 REAL SLOPE
2370 REAL SMC(NSOIL)
2371 REAL SMCDRY
2372 REAL SMCMAX
2373 REAL SMCREF
2374 REAL SMCWLT
2375 REAL STC(NSOIL)
2376 REAL T1
2377 REAL T24
2378 REAL TBOT
2379 REAL TH2
2380 REAL YY
2381 REAL YYNUM
2382 REAL ZBOT
2383 REAL ZSOIL(NSOIL)
2384 REAL ZZ1
2385
2386 REAL EC1
2387 REAL EDIR1
2388 REAL ET1(NSOIL)
2389 REAL ETT1
2390
2391 INTEGER K
2392
2393 PARAMETER(CP = 1004.5)
2394 PARAMETER(SIGMA = 5.67E-8)
2395
2396 ! ----------------------------------------------------------------------
2397 ! EXECUTABLE CODE BEGINS HERE:
2398 ! CONVERT ETP FROM KG M-2 S-1 TO MS-1 AND INITIALIZE DEW.
2399 ! ----------------------------------------------------------------------
2400 PRCP1 = PRCP * 0.001
2401 ETP1 = ETP * 0.001
2402 DEW = 0.0
2403
2404 EDIR = 0.
2405 EDIR1 = 0.
2406 EC = 0.
2407 EC1 = 0.
2408 DO K = 1,NSOIL
2409 ET(K) = 0.
2410 ET1(K) = 0.
2411 END DO
2412 ETT = 0.
2413 ETT1 = 0.
2414
2415 IF (ETP .GT. 0.0) THEN
2416
2417 ! ----------------------------------------------------------------------
2418 ! CONVERT PRCP FROM 'KG M-2 S-1' TO 'M S-1'.
2419 ! ----------------------------------------------------------------------
2420 CALL EVAPO (ETA1,SMC,NSOIL,CMC,ETP1,DT,ZSOIL, &
2421 & SH2O,SMCMAX,BEXP,PC,SMCWLT,DKSAT,DWSAT, &
2422 & SMCREF,SHDFAC,CMCMAX, &
2423 & SMCDRY,CFACTR, &
2424 & EDIR1,EC1,ET1,ETT1,SFCTMP,Q2,NROOT,RTDIS,FXEXP)
2425 CALL SMFLX (SMC,NSOIL,CMC,DT,PRCP1,ZSOIL, &
2426 & SH2O,SLOPE,KDT,FRZFACT, &
2427 & SMCMAX,BEXP,SMCWLT,DKSAT,DWSAT, &
2428 & SHDFAC,CMCMAX, &
2429 & RUNOFF1,RUNOFF2,RUNOFF3, &
2430 & EDIR1,EC1,ET1, &
2431 & DRIP)
2432
2433 ! ----------------------------------------------------------------------
2434 ! CONVERT MODELED EVAPOTRANSPIRATION FM M S-1 TO KG M-2 S-1
2435 ! ----------------------------------------------------------------------
2436 ! ETA = ETA1 * 1000.0
2437
2438 ! ----------------------------------------------------------------------
2439 ! EDIR = EDIR1 * 1000.0
2440 ! EC = EC1 * 1000.0
2441 ! ETT = ETT1 * 1000.0
2442 ! ET(1) = ET1(1) * 1000.0
2443 ! ET(2) = ET1(2) * 1000.0
2444 ! ET(3) = ET1(3) * 1000.0
2445 ! ET(4) = ET1(4) * 1000.0
2446 ! ----------------------------------------------------------------------
2447
2448 ELSE
2449
2450 ! ----------------------------------------------------------------------
2451 ! IF ETP < 0, ASSUME DEW FORMS (TRANSFORM ETP1 INTO DEW AND REINITIALIZE
2452 ! ETP1 TO ZERO).
2453 ! ----------------------------------------------------------------------
2454 DEW = -ETP1
2455 ! ETP1 = 0.0
2456
2457 ! ----------------------------------------------------------------------
2458 ! CONVERT PRCP FROM 'KG M-2 S-1' TO 'M S-1' AND ADD DEW AMOUNT.
2459 ! ----------------------------------------------------------------------
2460 PRCP1 = PRCP1 + DEW
2461 !
2462 ! CALL EVAPO (ETA1,SMC,NSOIL,CMC,ETP1,DT,ZSOIL,
2463 ! & SH2O,SMCMAX,BEXP,PC,SMCWLT,DKSAT,DWSAT,
2464 ! & SMCREF,SHDFAC,CMCMAX,
2465 ! & SMCDRY,CFACTR,
2466 ! & EDIR1,EC1,ET1,ETT,SFCTMP,Q2,NROOT,RTDIS,FXEXP)
2467 CALL SMFLX (SMC,NSOIL,CMC,DT,PRCP1,ZSOIL, &
2468 & SH2O,SLOPE,KDT,FRZFACT, &
2469 & SMCMAX,BEXP,SMCWLT,DKSAT,DWSAT, &
2470 & SHDFAC,CMCMAX, &
2471 & RUNOFF1,RUNOFF2,RUNOFF3, &
2472 & EDIR1,EC1,ET1, &
2473 & DRIP)
2474
2475 ! ----------------------------------------------------------------------
2476 ! CONVERT MODELED EVAPOTRANSPIRATION FROM 'M S-1' TO 'KG M-2 S-1'.
2477 ! ----------------------------------------------------------------------
2478 ! ETA = ETA1 * 1000.0
2479
2480 ! ----------------------------------------------------------------------
2481 ! EDIR = EDIR1 * 1000.0
2482 ! EC = EC1 * 1000.0
2483 ! ETT = ETT1 * 1000.0
2484 ! ET(1) = ET1(1) * 1000.0
2485 ! ET(2) = ET1(2) * 1000.0
2486 ! ET(3) = ET1(3) * 1000.0
2487 ! ET(4) = ET1(4) * 1000.0
2488 ! ----------------------------------------------------------------------
2489
2490 ENDIF
2491
2492 ! ----------------------------------------------------------------------
2493 ! CONVERT MODELED EVAPOTRANSPIRATION FM M S-1 TO KG M-2 S-1
2494 ! ----------------------------------------------------------------------
2495 ETA = ETA1 * 1000.0
2496
2497 ! ----------------------------------------------------------------------
2498 EDIR = EDIR1 * 1000.0
2499 EC = EC1 * 1000.0
2500 DO K = 1,NSOIL
2501 ET(K) = ET1(K) * 1000.0
2502 ! ET(1) = ET1(1) * 1000.0
2503 ! ET(2) = ET1(2) * 1000.0
2504 ! ET(3) = ET1(3) * 1000.0
2505 ! ET(4) = ET1(4) * 1000.0
2506 ENDDO
2507 ETT = ETT1 * 1000.0
2508 ! ----------------------------------------------------------------------
2509
2510 ! ----------------------------------------------------------------------
2511 ! BASED ON ETP AND E VALUES, DETERMINE BETA
2512 ! ----------------------------------------------------------------------
2513 IF ( ETP .LE. 0.0 ) THEN
2514 BETA = 0.0
2515 IF ( ETP .LT. 0.0 ) THEN
2516 BETA = 1.0
2517 ! ETA = ETP
2518 ENDIF
2519 ELSE
2520 BETA = ETA / ETP
2521 ENDIF
2522
2523 ! ----------------------------------------------------------------------
2524 ! GET SOIL THERMAL DIFFUXIVITY/CONDUCTIVITY FOR TOP SOIL LYR,
2525 ! CALC. ADJUSTED TOP LYR SOIL TEMP AND ADJUSTED SOIL FLUX, THEN
2526 ! CALL SHFLX TO COMPUTE/UPDATE SOIL HEAT FLUX AND SOIL TEMPS.
2527 ! ----------------------------------------------------------------------
2528 CALL TDFCND (DF1,SMC(1),QUARTZ,SMCMAX,SH2O(1))
2529
2530 ! ----------------------------------------------------------------------
2531 ! VEGETATION GREENNESS FRACTION REDUCTION IN SUBSURFACE HEAT FLUX
2532 ! VIA REDUCTION FACTOR, WHICH IS CONVENIENT TO APPLY HERE TO THERMAL
2533 ! DIFFUSIVITY THAT IS LATER USED IN HRT TO COMPUTE SUB SFC HEAT FLUX
2534 ! (SEE ADDITIONAL COMMENTS ON VEG EFFECT SUB-SFC HEAT FLX IN
2535 ! ROUTINE SFLX)
2536 ! ----------------------------------------------------------------------
2537 DF1 = DF1 * EXP(SBETA*SHDFAC)
2538
2539 ! ----------------------------------------------------------------------
2540 ! COMPUTE INTERMEDIATE TERMS PASSED TO ROUTINE HRT (VIA ROUTINE
2541 ! SHFLX BELOW) FOR USE IN COMPUTING SUBSURFACE HEAT FLUX IN HRT
2542 ! ----------------------------------------------------------------------
2543 YYNUM = FDOWN - SIGMA * T24
2544 YY = SFCTMP + (YYNUM/RCH+TH2-SFCTMP-BETA*EPSCA) / RR
2545 ZZ1 = DF1 / ( -0.5 * ZSOIL(1) * RCH * RR ) + 1.0
2546
2547 CALL SHFLX (SSOIL,STC,SMC,SMCMAX,NSOIL,T1,DT,YY,ZZ1,ZSOIL, &
2548 & TBOT,ZBOT,SMCWLT,PSISAT,SH2O,BEXP,F1,DF1,ICE, &
2549 & QUARTZ,CSOIL)
2550
2551 ! ----------------------------------------------------------------------
2552 ! SET FLX1 AND FLX3 (SNOPACK PHASE CHANGE HEAT FLUXES) TO ZERO SINCE
2553 ! THEY ARE NOT USED HERE IN SNOPAC. FLX2 (FREEZING RAIN HEAT FLUX) WAS
2554 ! SIMILARLY INITIALIZED IN THE PENMAN ROUTINE.
2555 ! ----------------------------------------------------------------------
2556 FLX1 = 0.0
2557 FLX3 = 0.0
2558
2559 ! ----------------------------------------------------------------------
2560 ! END SUBROUTINE NOPAC
2561 ! ----------------------------------------------------------------------
2562 END SUBROUTINE NOPAC
2563
2564 SUBROUTINE PENMAN (SFCTMP,SFCPRS,CH,T2V,TH2,PRCP,FDOWN,T24,SSOIL, &
2565 & Q2,Q2SAT,ETP,RCH,EPSCA,RR,SNOWNG,FRZGRA, &
2566 & DQSDT2,FLX2)
2567
2568 IMPLICIT NONE
2569
2570 ! ----------------------------------------------------------------------
2571 ! SUBROUTINE PENMAN
2572 ! ----------------------------------------------------------------------
2573 ! CALCULATE POTENTIAL EVAPORATION FOR THE CURRENT POINT. VARIOUS
2574 ! PARTIAL SUMS/PRODUCTS ARE ALSO CALCULATED AND PASSED BACK TO THE
2575 ! CALLING ROUTINE FOR LATER USE.
2576 ! ----------------------------------------------------------------------
2577 LOGICAL SNOWNG
2578 LOGICAL FRZGRA
2579
2580 REAL A
2581 REAL BETA
2582 REAL CH
2583 REAL CP
2584 REAL CPH2O
2585 REAL CPICE
2586 REAL DELTA
2587 REAL DQSDT2
2588 REAL ELCP
2589 REAL EPSCA
2590 REAL ETP
2591 REAL FDOWN
2592 REAL FLX2
2593 REAL FNET
2594 REAL LSUBC
2595 REAL LSUBF
2596 REAL PRCP
2597 REAL Q2
2598 REAL Q2SAT
2599 REAL R
2600 REAL RAD
2601 REAL RCH
2602 REAL RHO
2603 REAL RR
2604 REAL SSOIL
2605 REAL SFCPRS
2606 REAL SFCTMP
2607 REAL SIGMA
2608 REAL T24
2609 REAL T2V
2610 REAL TH2
2611
2612 PARAMETER(CP = 1004.6)
2613 PARAMETER(CPH2O = 4.218E+3)
2614 PARAMETER(CPICE = 2.106E+3)
2615 PARAMETER(R = 287.04)
2616 PARAMETER(ELCP = 2.4888E+3)
2617 PARAMETER(LSUBF = 3.335E+5)
2618 PARAMETER(LSUBC = 2.501000E+6)
2619 PARAMETER(SIGMA = 5.67E-8)
2620
2621 ! ----------------------------------------------------------------------
2622 ! EXECUTABLE CODE BEGINS HERE:
2623 ! ----------------------------------------------------------------------
2624 FLX2 = 0.0
2625
2626 ! ----------------------------------------------------------------------
2627 ! PREPARE PARTIAL QUANTITIES FOR PENMAN EQUATION.
2628 ! ----------------------------------------------------------------------
2629 DELTA = ELCP * DQSDT2
2630 T24 = SFCTMP * SFCTMP * SFCTMP * SFCTMP
2631 RR = T24 * 6.48E-8 /(SFCPRS * CH) + 1.0
2632 RHO = SFCPRS / (R * T2V)
2633 RCH = RHO * CP * CH
2634
2635 ! ----------------------------------------------------------------------
2636 ! ADJUST THE PARTIAL SUMS / PRODUCTS WITH THE LATENT HEAT
2637 ! EFFECTS CAUSED BY FALLING PRECIPITATION.
2638 ! ----------------------------------------------------------------------
2639 IF (.NOT. SNOWNG) THEN
2640 IF (PRCP .GT. 0.0) RR = RR + CPH2O*PRCP/RCH
2641 ELSE
2642 RR = RR + CPICE*PRCP/RCH
2643 ENDIF
2644
2645 FNET = FDOWN - SIGMA*T24 - SSOIL
2646
2647 ! ----------------------------------------------------------------------
2648 ! INCLUDE THE LATENT HEAT EFFECTS OF FRZNG RAIN CONVERTING TO ICE ON
2649 ! IMPACT IN THE CALCULATION OF FLX2 AND FNET.
2650 ! ----------------------------------------------------------------------
2651 IF (FRZGRA) THEN
2652 FLX2 = -LSUBF * PRCP
2653 FNET = FNET - FLX2
2654 ENDIF
2655
2656 ! ----------------------------------------------------------------------
2657 ! FINISH PENMAN EQUATION CALCULATIONS.
2658 ! ----------------------------------------------------------------------
2659 RAD = FNET/RCH + TH2 - SFCTMP
2660 A = ELCP * (Q2SAT - Q2)
2661 EPSCA = (A*RR + RAD*DELTA) / (DELTA + RR)
2662 ETP = EPSCA * RCH / LSUBC
2663
2664 ! ----------------------------------------------------------------------
2665 ! END SUBROUTINE PENMAN
2666 ! ----------------------------------------------------------------------
2667 END SUBROUTINE PENMAN
2668
2669 SUBROUTINE REDPRM ( &
2670 & VEGTYP,SOILTYP,SLOPETYP, &
2671 & CFACTR,CMCMAX,RSMAX,TOPT,REFKDT,KDT,SBETA, &
2672 & SHDFAC,RSMIN,RGL,HS,ZBOT,FRZX,PSISAT,SLOPE, &
2673 & SNUP,SALP,BEXP,DKSAT,DWSAT, &
2674 & SMCMAX,SMCWLT,SMCREF, &
2675 & SMCDRY,F1,QUARTZ,FXEXP,RTDIS,SLDPTH,ZSOIL, &
2676 & NROOT,NSOIL,Z0,CZIL,LAI,CSOIL,PTU)
2677
2678 IMPLICIT NONE
2679 ! ----------------------------------------------------------------------
2680 ! SUBROUTINE REDPRM
2681 ! ----------------------------------------------------------------------
2682 ! INTERNALLY SET (DEFAULT VALUESS), OR OPTIONALLY READ-IN VIA NAMELIST
2683 ! I/O, ALL SOIL AND VEGETATION PARAMETERS REQUIRED FOR THE EXECUSION OF
2684 ! THE NOAH LSM.
2685 !
2686 ! OPTIONAL NON-DEFAULT PARAMETERS CAN BE READ IN, ACCOMMODATING UP TO 30
2687 ! SOIL, VEG, OR SLOPE CLASSES, IF THE DEFAULT MAX NUMBER OF SOIL, VEG,
2688 ! AND/OR SLOPE TYPES IS RESET.
2689 !
2690 ! FUTURE UPGRADES OF ROUTINE REDPRM MUST EXPAND TO INCORPORATE SOME OF
2691 ! THE EMPIRICAL PARAMETERS OF THE FROZEN SOIL AND SNOWPACK PHYSICS (SUCH
2692 ! AS IN ROUTINES FRH2O, SNOWPACK, AND SNOW_NEW) NOT YET SET IN THIS
2693 ! REDPRM ROUTINE, BUT RATHER SET IN LOWER LEVEL SUBROUTINES.
2694 !
2695 ! SET MAXIMUM NUMBER OF SOIL-, VEG-, AND SLOPETYP IN DATA STATEMENT.
2696 ! ----------------------------------------------------------------------
2697 INTEGER MAX_SLOPETYP
2698 INTEGER MAX_SOILTYP
2699 INTEGER MAX_VEGTYP
2700
2701 PARAMETER(MAX_SLOPETYP = 30)
2702 PARAMETER(MAX_SOILTYP = 30)
2703 PARAMETER(MAX_VEGTYP = 30)
2704
2705 ! ----------------------------------------------------------------------
2706 ! NUMBER OF DEFINED SOIL-, VEG-, AND SLOPETYPS USED.
2707 ! ----------------------------------------------------------------------
2708 INTEGER DEFINED_VEG
2709 INTEGER DEFINED_SOIL
2710 INTEGER DEFINED_SLOPE
2711
2712 DATA DEFINED_VEG/27/
2713 DATA DEFINED_SOIL/19/
2714 DATA DEFINED_SLOPE/9/
2715
2716 ! ----------------------------------------------------------------------
2717 ! SET-UP SOIL PARAMETERS FOR GIVEN SOIL TYPE
2718 ! INPUT: SOLTYP: SOIL TYPE (INTEGER INDEX)
2719 ! OUTPUT: SOIL PARAMETERS:
2720 ! MAXSMC: MAX SOIL MOISTURE CONTENT (POROSITY)
2721 ! REFSMC: REFERENCE SOIL MOISTURE (ONSET OF SOIL MOISTURE
2722 ! STRESS IN TRANSPIRATION)
2723 ! WLTSMC: WILTING PT SOIL MOISTURE CONTENTS
2724 ! DRYSMC: AIR DRY SOIL MOIST CONTENT LIMITS
2725 ! SATPSI: SATURATED SOIL POTENTIAL
2726 ! SATDK: SATURATED SOIL HYDRAULIC CONDUCTIVITY
2727 ! BB: THE 'B' PARAMETER
2728 ! SATDW: SATURATED SOIL DIFFUSIVITY
2729 ! F11: USED TO COMPUTE SOIL DIFFUSIVITY/CONDUCTIVITY
2730 ! QUARTZ: SOIL QUARTZ CONTENT
2731 ! ----------------------------------------------------------------------
2732 ! SOIL STATSGO
2733 ! TYPE CLASS
2734 ! ---- -------
2735 ! 1 SAND
2736 ! 2 LOAMY SAND
2737 ! 3 SANDY LOAM
2738 ! 4 SILT LOAM
2739 ! 5 SILT
2740 ! 6 LOAM
2741 ! 7 SANDY CLAY LOAM
2742 ! 8 SILTY CLAY LOAM
2743 ! 9 CLAY LOAM
2744 ! 10 SANDY CLAY
2745 ! 11 SILTY CLAY
2746 ! 12 CLAY
2747 ! 13 ORGANIC MATERIAL
2748 ! 14 WATER
2749 ! 15 BEDROCK
2750 ! 16 OTHER(land-ice)
2751 ! 17 PLAYA
2752 ! 18 LAVA
2753 ! 19 WHITE SAND
2754 ! ----------------------------------------------------------------------
2755
2756 REAL BB(MAX_SOILTYP)
2757 REAL DRYSMC(MAX_SOILTYP)
2758 REAL F11(MAX_SOILTYP)
2759 REAL MAXSMC(MAX_SOILTYP)
2760 REAL REFSMC(MAX_SOILTYP)
2761 REAL SATPSI(MAX_SOILTYP)
2762 REAL SATDK(MAX_SOILTYP)
2763 REAL SATDW(MAX_SOILTYP)
2764 REAL WLTSMC(MAX_SOILTYP)
2765 REAL QTZ(MAX_SOILTYP)
2766
2767 REAL BEXP
2768 REAL DKSAT
2769 REAL DWSAT
2770 REAL F1
2771 REAL PTU
2772 REAL QUARTZ
2773 REAL REFSMC1
2774 REAL SMCDRY
2775 REAL SMCMAX
2776 REAL SMCREF
2777 REAL SMCWLT
2778 REAL WLTSMC1
2779
2780 ! ----------------------------------------------------------------------
2781 ! SOIL TEXTURE-RELATED ARRAYS.
2782 ! ----------------------------------------------------------------------
2783 DATA MAXSMC/0.395, 0.421, 0.434, 0.476, 0.476, 0.439, &
2784 & 0.404, 0.464, 0.465, 0.406, 0.468, 0.457, &
2785 & 0.464, 0.464, 0.200, 0.421, 0.457, 0.200, &
2786 & 0.395, 0.000, 0.000, 0.000, 0.000, 0.000, &
2787 & 0.000, 0.000, 0.000, 0.000, 0.000, 0.000/
2788 ! ----------------------------------------------------------------------
2789 DATA SATPSI/0.0350, 0.0363, 0.1413, 0.7586, 0.7586, 0.3548, &
2790 & 0.1349, 0.6166, 0.2630, 0.0977, 0.3236, 0.4677, &
2791 & 0.3548, 0.3548, 0.0350, 0.0363, 0.4677, 0.0350, &
2792 & 0.0350, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, &
2793 & 0.000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000/
2794 ! ----------------------------------------------------------------------
2795 DATA SATDK /1.7600E-4, 1.4078E-5, 5.2304E-6, 2.8089E-6, 2.8089E-6,&
2796 & 3.3770E-6, 4.4518E-6, 2.0348E-6, 2.4464E-6, 7.2199E-6,&
2797 & 1.3444E-6, 9.7394E-7, 3.3770E-6, 3.3770E-6, 1.4078E-5,&
2798 & 1.4078E-5, 9.7394E-7, 1.4078E-5, 1.7600E-4, 0.0,&
2799 & 0.0, 0.0, 0.0, 0.0, 0.0,&
2800 & 0.0, 0.0, 0.0, 0.0, 0.0/
2801 ! ----------------------------------------------------------------------
2802 DATA BB /4.05, 4.26, 4.74, 5.33, 5.33, 5.25, &
2803 & 6.77, 8.72, 8.17, 10.73, 10.39, 11.55, &
2804 & 5.25, 5.25, 4.05, 4.26, 11.55, 4.05, &
2805 & 4.05, 0.00, 0.00, 0.00, 0.00, 0.00, &
2806 & 0.00, 0.00, 0.00, 0.00, 0.00, 0.00/
2807 ! ----------------------------------------------------------------------
2808 DATA QTZ /0.92, 0.82, 0.60, 0.25, 0.10, 0.40, &
2809 & 0.60, 0.10, 0.35, 0.52, 0.10, 0.25, &
2810 & 0.05, 0.05, 0.07, 0.25, 0.60, 0.52, &
2811 & 0.92, 0.00, 0.00, 0.00, 0.00, 0.00, &
2812 & 0.00, 0.00, 0.00, 0.00, 0.00, 0.00/
2813
2814 ! ----------------------------------------------------------------------
2815 ! THE FOLLOWING 5 PARAMETERS ARE DERIVED LATER IN REDPRM.F FROM THE SOIL
2816 ! DATA, AND ARE JUST GIVEN HERE FOR REFERENCE AND TO FORCE STATIC
2817 ! STORAGE ALLOCATION. -DAG LOHMANN, FEB. 2001
2818 ! ----------------------------------------------------------------------
2819 ! !!!!!!!!!!!!!! The following values in the table are NOT used
2820 ! !!!!!!!!!!!!!! and are just given for reference
2821 DATA REFSMC/0.196, 0.248, 0.282, 0.332, 0.332, 0.301, &
2822 & 0.293, 0.368, 0.361, 0.320, 0.388, 0.389, &
2823 & 0.319, 0.000, 0.116, 0.248, 0.389, 0.116, &
2824 & 0.196, 0.000, 0.000, 0.000, 0.000, 0.000, &
2825 & 0.000, 0.000, 0.000, 0.000, 0.000, 0.000/
2826 ! !!!!!!!!!!!!!! The following values in the table are NOT used
2827 ! !!!!!!!!!!!!!! and are just given for reference
2828 DATA WLTSMC/0.023, 0.028, 0.047, 0.084, 0.084, 0.066, &
2829 & 0.069, 0.120, 0.103, 0.100, 0.126, 0.135, &
2830 & 0.069, 0.000, 0.012, 0.028, 0.135, 0.012, &
2831 & 0.023, 0.000, 0.000, 0.000, 0.000, 0.000, &
2832 & 0.000, 0.000, 0.000, 0.000, 0.000, 0.000/
2833 ! !!!!!!!!!!!!!! The following values in the table are NOT used
2834 ! !!!!!!!!!!!!!! and are just given for reference
2835 DATA DRYSMC/0.023, 0.028, 0.047, 0.084, 0.084, 0.066, &
2836 & 0.069, 0.120, 0.103, 0.100, 0.126, 0.135, &
2837 & 0.069, 0.000, 0.012, 0.028, 0.135, 0.012, &
2838 & 0.023, 0.000, 0.000, 0.000, 0.000, 0.000, &
2839 & 0.000, 0.000, 0.000, 0.000, 0.000, 0.000/
2840
2841 ! !!!!!!!!!!!!!! The following values in the table are NOT used
2842 ! !!!!!!!!!!!!!! and are just given for reference
2843 DATA SATDW /0.632E-4, 0.517E-5, 0.807E-5, 0.239E-4, 0.239E-4, &
2844 & 0.143E-4, 0.101E-4, 0.236E-4, 0.113E-4, 0.186E-4, &
2845 & 0.966E-5, 0.115E-4, 0.136E-4, 0.0, 0.998E-5, &
2846 & 0.517E-5, 0.115E-4, 0.998E-5, 0.632E-4, 0.0, &
2847 & 0.0, 0.0, 0.0, 0.0, 0.0, &
2848 & 0.0, 0.0, 0.0, 0.0, 0.0/
2849 ! !!!!!!!!!!!!!! The following values in the table are NOT used
2850 ! !!!!!!!!!!!!!! and are just given for reference
2851 DATA F11 /-1.090, -1.041, -0.568, 0.162, 0.162, -0.327, &
2852 & -1.535, -1.118, -1.297, -3.211, -1.916, -2.258, &
2853 & -0.201, 0.000, -2.287, -1.041, -2.258, -2.287, &
2854 & -1.090, 0.000, 0.000, 0.000, 0.000, 0.000, &
2855 & 0.000, 0.000, 0.000, 0.000, 0.000, 0.000/
2856
2857 ! ----------------------------------------------------------------------
2858 ! SET-UP VEGETATION PARAMETERS FOR A GIVEN VEGETAION TYPE:
2859 ! INPUT: VEGTYP = VEGETATION TYPE (INTEGER INDEX)
2860 ! OUPUT: VEGETATION PARAMETERS
2861 ! SHDFAC: VEGETATION GREENNESS FRACTION
2862 ! RSMIN: MIMIMUM STOMATAL RESISTANCE
2863 ! RGL: PARAMETER USED IN SOLAR RAD TERM OF
2864 ! CANOPY RESISTANCE FUNCTION
2865 ! HS: PARAMETER USED IN VAPOR PRESSURE DEFICIT TERM OF
2866 ! CANOPY RESISTANCE FUNCTION
2867 ! SNUP: THRESHOLD SNOW DEPTH (IN WATER EQUIVALENT M) THAT
2868 ! IMPLIES 100% SNOW COVER
2869 ! ----------------------------------------------------------------------
2870 ! CLASS USGS-WRF VEGETATION/SURFACE TYPE
2871 ! 1 Urban and Built-Up Land
2872 ! 2 Dryland Cropland and Pasture
2873 ! 3 Irrigated Cropland and Pasture
2874 ! 4 Mixed Dryland/Irrigated Cropland and Pasture
2875 ! 5 Cropland/Grassland Mosaic
2876 ! 6 Cropland/Woodland Mosaic
2877 ! 7 Grassland
2878 ! 8 Shrubland
2879 ! 9 Mixed Shrubland/Grassland
2880 ! 10 Savanna
2881 ! 11 Deciduous Broadleaf Forest
2882 ! 12 Deciduous Needleleaf Forest
2883 ! 13 Evergreen Broadleaf Forest
2884 ! 14 Evergreen Needleleaf Forest
2885 ! 15 Mixed Forest
2886 ! 16 Water Bodies
2887 ! 17 Herbaceous Wetland
2888 ! 18 Wooded Wetland
2889 ! 19 Barren or Sparsely Vegetated
2890 ! 20 Herbaceous Tundra
2891 ! 21 Wooded Tundra
2892 ! 22 Mixed Tundra
2893 ! 23 Bare Ground Tundra
2894 ! 24 Snow or Ice
2895 ! 25 Playa
2896 ! 26 Lava
2897 ! 27 White Sand
2898 ! ----------------------------------------------------------------------
2899
2900 INTEGER NROOT
2901 INTEGER NROOT_DATA(MAX_VEGTYP)
2902
2903 REAL FRZFACT
2904 REAL HS
2905 REAL HSTBL(MAX_VEGTYP)
2906 REAL LAI
2907 REAL LAI_DATA(MAX_VEGTYP)
2908 REAL PSISAT
2909 REAL RSMIN
2910 REAL RGL
2911 REAL RGLTBL(MAX_VEGTYP)
2912 REAL RSMTBL(MAX_VEGTYP)
2913 REAL SHDFAC
2914 REAL SNUP
2915 REAL SNUPX(MAX_VEGTYP)
2916 REAL Z0
2917 REAL Z0_DATA(MAX_VEGTYP)
2918
2919 ! ----------------------------------------------------------------------
2920 ! VEGETATION CLASS-RELATED ARRAYS
2921 ! ----------------------------------------------------------------------
2922 ! DATA NROOT_DATA /2,3,3,3,3,3,3,3,3,3,
2923 ! & 4,4,4,4,4,2,2,2,2,3,
2924 ! & 3,3,2,2,2,2,2,0,0,0/
2925 DATA NROOT_DATA /1,3,3,3,3,3,3,3,3,3, &
2926 & 4,4,4,4,4,0,2,2,1,3, &
2927 & 3,3,2,1,1,1,1,0,0,0/
2928 DATA RSMTBL /200.0, 70.0, 70.0, 70.0, 70.0, 70.0, &
2929 & 70.0, 300.0, 170.0, 70.0, 100.0, 150.0, &
2930 & 150.0, 125.0, 125.0, 100.0, 40.0, 100.0, &
2931 & 300.0, 150.0, 150.0, 150.0, 200.0, 200.0, &
2932 & 40.0, 100.0, 300.0, 0.0, 0.0, 0.0/
2933 DATA RGLTBL /100.0, 100.0, 100.0, 100.0, 100.0, 65.0, &
2934 & 100.0, 100.0, 100.0, 65.0, 30.0, 30.0, &
2935 & 30.0, 30.0, 30.0, 30.0, 100.0, 30.0, &
2936 & 100.0, 100.0, 100.0, 100.0, 100.0, 100.0, &
2937 & 100.0, 100.0, 100.0, 0.0, 0.0, 0.0/
2938 DATA HSTBL /42.00, 36.25, 36.25, 36.25, 36.25, 44.14, &
2939 & 36.35, 42.00, 39.18, 54.53, 54.53, 47.35, &
2940 & 41.69, 47.35, 51.93, 51.75, 60.00, 51.93, &
2941 & 42.00, 42.00, 42.00, 42.00, 42.00, 42.00, &
2942 & 36.25, 42.00, 42.00, 0.00, 0.00, 0.00/
2943 DATA SNUPX /0.020, 0.020, 0.020, 0.020, 0.020, 0.020, &
2944 & 0.020, 0.020, 0.020, 0.040, 0.040, 0.040, &
2945 & 0.040, 0.040, 0.040, 0.010, 0.013, 0.020, &
2946 & 0.013, 0.020, 0.020, 0.020, 0.020, 0.013, &
2947 & 0.013, 0.013, 0.013, 0.000, 0.000, 0.000/
2948 DATA Z0_DATA / 1.00, 0.07, 0.07, 0.07, 0.07, 0.15, &
2949 & 0.08, 0.03, 0.05, 0.86, 0.80, 0.85, &
2950 & 2.65, 1.09, 0.80, 0.001, 0.04, 0.05, &
2951 & 0.01, 0.04, 0.06, 0.05, 0.03, 0.001, &
2952 & 0.01, 0.15, 0.01, 0.00, 0.00, 0.00/
2953 DATA LAI_DATA /4.0, 4.0, 4.0, 4.0, 4.0, 4.0, &
2954 & 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, &
2955 & 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, &
2956 & 4.0, 4.0, 4.0, 4.0, 4.0, 4.0, &
2957 & 4.0, 4.0, 4.0, 0.0, 0.0, 0.0/
2958
2959 ! ----------------------------------------------------------------------
2960 ! CLASS PARAMETER 'SLOPETYP' WAS INCLUDED TO ESTIMATE LINEAR RESERVOIR
2961 ! COEFFICIENT 'SLOPE' TO THE BASEFLOW RUNOFF OUT OF THE BOTTOM LAYER.
2962 ! LOWEST CLASS (SLOPETYP=0) MEANS HIGHEST SLOPE PARAMETER = 1.
2963 ! DEFINITION OF SLOPETYP FROM 'ZOBLER' SLOPE TYPE:
2964 ! SLOPE CLASS PERCENT SLOPE
2965 ! 1 0-8
2966 ! 2 8-30
2967 ! 3 > 30
2968 ! 4 0-30
2969 ! 5 0-8 & > 30
2970 ! 6 8-30 & > 30
2971 ! 7 0-8, 8-30, > 30
2972 ! 9 GLACIAL ICE
2973 ! BLANK OCEAN/SEA
2974 ! ----------------------------------------------------------------------
2975 ! NOTE:
2976 ! CLASS 9 FROM 'ZOBLER' FILE SHOULD BE REPLACED BY 8 AND 'BLANK' 9
2977 ! ----------------------------------------------------------------------
2978 REAL SLOPE
2979 REAL SLOPE_DATA(MAX_SLOPETYP)
2980
2981 DATA SLOPE_DATA /0.1, 0.6, 1.0, 0.35, 0.55, 0.8, &
2982 & 0.63, 0.0, 0.0, 0.0, 0.0, 0.0, &
2983 & 0.0 , 0.0, 0.0, 0.0, 0.0, 0.0, &
2984 & 0.0 , 0.0, 0.0, 0.0, 0.0, 0.0, &
2985 & 0.0 , 0.0, 0.0, 0.0, 0.0, 0.0/
2986
2987 ! ----------------------------------------------------------------------
2988 ! SET NAMELIST FILE NAME
2989 ! ----------------------------------------------------------------------
2990 CHARACTER*50 NAMELIST_NAME
2991
2992 ! ----------------------------------------------------------------------
2993 ! SET UNIVERSAL PARAMETERS (NOT DEPENDENT ON SOIL, VEG, SLOPE TYPE)
2994 ! ----------------------------------------------------------------------
2995 INTEGER I
2996 INTEGER NSOIL
2997 INTEGER SLOPETYP
2998 INTEGER SOILTYP
2999 INTEGER VEGTYP
3000
3001 INTEGER BARE
3002 ! DATA BARE /11/
3003 DATA BARE /19/
3004
3005 LOGICAL LPARAM
3006 DATA LPARAM /.TRUE./
3007
3008 LOGICAL LFIRST
3009 DATA LFIRST /.TRUE./
3010
3011 ! ----------------------------------------------------------------------
3012 ! PARAMETER USED TO CALCULATE ROUGHNESS LENGTH OF HEAT.
3013 ! ----------------------------------------------------------------------
3014 REAL CZIL
3015 REAL CZIL_DATA
3016 ! changed in version 2.6 June 2nd 2003
3017 ! DATA CZIL_DATA /0.2/
3018 DATA CZIL_DATA /0.1/
3019
3020 ! ----------------------------------------------------------------------
3021 ! PARAMETER USED TO CALUCULATE VEGETATION EFFECT ON SOIL HEAT FLUX.
3022 ! ----------------------------------------------------------------------
3023 REAL SBETA
3024 REAL SBETA_DATA
3025 DATA SBETA_DATA /-2.0/
3026
3027 ! ----------------------------------------------------------------------
3028 ! BARE SOIL EVAPORATION EXPONENT USED IN DEVAP.
3029 ! ----------------------------------------------------------------------
3030 REAL FXEXP
3031 REAL FXEXP_DATA
3032 DATA FXEXP_DATA /2.0/
3033
3034 ! ----------------------------------------------------------------------
3035 ! SOIL HEAT CAPACITY [J M-3 K-1]
3036 ! ----------------------------------------------------------------------
3037 REAL CSOIL
3038 REAL CSOIL_DATA
3039 ! DATA CSOIL_DATA /1.26E+6/
3040 DATA CSOIL_DATA /2.00E+6/
3041
3042 ! ----------------------------------------------------------------------
3043 ! SPECIFY SNOW DISTRIBUTION SHAPE PARAMETER SALP - SHAPE PARAMETER OF
3044 ! DISTRIBUTION FUNCTION OF SNOW COVER. FROM ANDERSON'S DATA (HYDRO-17)
3045 ! BEST FIT IS WHEN SALP = 2.6
3046 ! ----------------------------------------------------------------------
3047 REAL SALP
3048 REAL SALP_DATA
3049 ! changed for version 2.6 June 2nd 2003
3050 ! DATA SALP_DATA /2.6/
3051 DATA SALP_DATA /4.0/
3052
3053 ! ----------------------------------------------------------------------
3054 ! KDT IS DEFINED BY REFERENCE REFKDT AND DKSAT; REFDK=2.E-6 IS THE SAT.
3055 ! DK. VALUE FOR THE SOIL TYPE 2
3056 ! ----------------------------------------------------------------------
3057 REAL REFDK
3058 REAL REFDK_DATA
3059 DATA REFDK_DATA /2.0E-6/
3060
3061 REAL REFKDT
3062 REAL REFKDT_DATA
3063 DATA REFKDT_DATA /3.0/
3064
3065 REAL FRZX
3066 REAL KDT
3067
3068 ! ----------------------------------------------------------------------
3069 ! FROZEN GROUND PARAMETER, FRZK, DEFINITION: ICE CONTENT THRESHOLD ABOVE
3070 ! WHICH FROZEN SOIL IS IMPERMEABLE REFERENCE VALUE OF THIS PARAMETER FOR
3071 ! THE LIGHT CLAY SOIL (TYPE=3) FRZK = 0.15 M.
3072 ! ----------------------------------------------------------------------
3073 REAL FRZK
3074 REAL FRZK_DATA
3075 DATA FRZK_DATA /0.15/
3076
3077 REAL RTDIS(NSOIL)
3078 REAL SLDPTH(NSOIL)
3079 REAL ZSOIL(NSOIL)
3080
3081 ! ----------------------------------------------------------------------
3082 ! SET TWO CANOPY WATER PARAMETERS.
3083 ! ----------------------------------------------------------------------
3084 REAL CFACTR
3085 REAL CFACTR_DATA
3086 DATA CFACTR_DATA /0.5/
3087
3088 REAL CMCMAX
3089 REAL CMCMAX_DATA
3090 DATA CMCMAX_DATA /0.5E-3/
3091
3092 ! ----------------------------------------------------------------------
3093 ! SET MAX. STOMATAL RESISTANCE.
3094 ! ----------------------------------------------------------------------
3095 REAL RSMAX
3096 REAL RSMAX_DATA
3097 DATA RSMAX_DATA /5000.0/
3098
3099 ! ----------------------------------------------------------------------
3100 ! SET OPTIMUM TRANSPIRATION AIR TEMPERATURE.
3101 ! ----------------------------------------------------------------------
3102 REAL TOPT
3103 REAL TOPT_DATA
3104 DATA TOPT_DATA /298.0/
3105
3106 ! ----------------------------------------------------------------------
3107 ! SPECIFY DEPTH[M] OF LOWER BOUNDARY SOIL TEMPERATURE.
3108 ! ----------------------------------------------------------------------
3109 REAL ZBOT
3110 REAL ZBOT_DATA
3111 ! changed for version 2.5.2
3112 ! DATA ZBOT_DATA /-3.0/
3113 DATA ZBOT_DATA /-8.0/
3114
3115 ! ----------------------------------------------------------------------
3116 ! SET TWO SOIL MOISTURE WILT, SOIL MOISTURE REFERENCE PARAMETERS
3117 ! ----------------------------------------------------------------------
3118 REAL SMLOW
3119 REAL SMLOW_DATA
3120 DATA SMLOW_DATA /0.5/
3121
3122 REAL SMHIGH
3123 REAL SMHIGH_DATA
3124 ! changed in 2.6 from 3 to 6 on June 2nd 2003
3125 DATA SMHIGH_DATA /3.0/
3126 ! DATA SMHIGH_DATA /6.0/
3127
3128 ! ----------------------------------------------------------------------
3129 ! NAMELIST DEFINITION:
3130 ! ----------------------------------------------------------------------
3131 NAMELIST /SOIL_VEG/ SLOPE_DATA, RSMTBL, RGLTBL, HSTBL, SNUPX, &
3132 & BB, DRYSMC, F11, MAXSMC, REFSMC, SATPSI, SATDK, SATDW, &
3133 & WLTSMC, QTZ, LPARAM, ZBOT_DATA, SALP_DATA, CFACTR_DATA, &
3134 & CMCMAX_DATA, SBETA_DATA, RSMAX_DATA, TOPT_DATA, &
3135 & REFDK_DATA, FRZK_DATA, BARE, DEFINED_VEG, DEFINED_SOIL, &
3136 & DEFINED_SLOPE, FXEXP_DATA, NROOT_DATA, REFKDT_DATA, Z0_DATA, &
3137 & CZIL_DATA, LAI_DATA, CSOIL_DATA
3138
3139 ! ----------------------------------------------------------------------
3140 ! READ NAMELIST FILE TO OVERRIDE DEFAULT PARAMETERS ONLY ONCE.
3141 ! NAMELIST_NAME must be 50 characters or less.
3142 ! ----------------------------------------------------------------------
3143 IF (LFIRST) THEN
3144 ! WRITE(*,*) 'READ NAMELIST'
3145 ! OPEN(58, FILE = 'namelist_filename.txt')
3146 ! READ(58,'(A)') NAMELIST_NAME
3147 ! CLOSE(58)
3148 ! WRITE(*,*) 'Namelist Filename is ', NAMELIST_NAME
3149 ! OPEN(59, FILE = NAMELIST_NAME)
3150 ! 50 CONTINUE
3151 ! READ(59, SOIL_VEG, END=100)
3152 ! IF (LPARAM) GOTO 50
3153 ! 100 CONTINUE
3154 ! CLOSE(59)
3155 ! WRITE(*,NML=SOIL_VEG)
3156 LFIRST = .FALSE.
3157 IF (DEFINED_SOIL .GT. MAX_SOILTYP) THEN
3158 WRITE(*,*) 'Warning: DEFINED_SOIL too large in namelist'
3159 STOP 222
3160 ENDIF
3161 IF (DEFINED_VEG .GT. MAX_VEGTYP) THEN
3162 WRITE(*,*) 'Warning: DEFINED_VEG too large in namelist'
3163 STOP 222
3164 ENDIF
3165 IF (DEFINED_SLOPE .GT. MAX_SLOPETYP) THEN
3166 WRITE(*,*) 'Warning: DEFINED_SLOPE too large in namelist'
3167 STOP 222
3168 ENDIF
3169
3170 SMLOW = SMLOW_DATA
3171 SMHIGH = SMHIGH_DATA
3172
3173 DO I = 1,DEFINED_SOIL
3174 SATDW(I) = BB(I)*SATDK(I)*(SATPSI(I)/MAXSMC(I))
3175 F11(I) = ALOG10(SATPSI(I)) + BB(I)*ALOG10(MAXSMC(I)) + 2.0
3176 REFSMC1 = MAXSMC(I)*(5.79E-9/SATDK(I)) &
3177 & **(1.0/(2.0*BB(I)+3.0))
3178 REFSMC(I) = REFSMC1 + (MAXSMC(I)-REFSMC1) / SMHIGH
3179 WLTSMC1 = MAXSMC(I) * (200.0/SATPSI(I))**(-1.0/BB(I))
3180 WLTSMC(I) = WLTSMC1 - SMLOW * WLTSMC1
3181
3182 ! ----------------------------------------------------------------------
3183 ! CURRENT VERSION DRYSMC VALUES THAT EQUATE TO WLTSMC.
3184 ! FUTURE VERSION COULD LET DRYSMC BE INDEPENDENTLY SET VIA NAMELIST.
3185 ! ----------------------------------------------------------------------
3186 DRYSMC(I) = WLTSMC(I)
3187 END DO
3188
3189 ! ----------------------------------------------------------------------
3190 ! END LFIRST BLOCK
3191 ! ----------------------------------------------------------------------
3192 ENDIF
3193
3194 IF (SOILTYP .GT. DEFINED_SOIL) THEN
3195 WRITE(*,*) 'Warning: too many soil types'
3196 STOP 333
3197 ENDIF
3198 IF (VEGTYP .GT. DEFINED_VEG) THEN
3199 WRITE(*,*) 'Warning: too many veg types'
3200 STOP 333
3201 ENDIF
3202 IF (SLOPETYP .GT. DEFINED_SLOPE) THEN
3203 WRITE(*,*) 'Warning: too many slope types'
3204 STOP 333
3205 ENDIF
3206
3207 ! ----------------------------------------------------------------------
3208 ! SET-UP UNIVERSAL PARAMETERS (NOT DEPENDENT ON SOILTYP, VEGTYP OR
3209 ! SLOPETYP)
3210 ! ----------------------------------------------------------------------
3211 ZBOT = ZBOT_DATA
3212 SALP = SALP_DATA
3213 CFACTR = CFACTR_DATA
3214 CMCMAX = CMCMAX_DATA
3215 SBETA = SBETA_DATA
3216 RSMAX = RSMAX_DATA
3217 TOPT = TOPT_DATA
3218 REFDK = REFDK_DATA
3219 FRZK = FRZK_DATA
3220 FXEXP = FXEXP_DATA
3221 REFKDT = REFKDT_DATA
3222 CZIL = CZIL_DATA
3223 CSOIL = CSOIL_DATA
3224
3225 ! ----------------------------------------------------------------------
3226 ! SET-UP SOIL PARAMETERS
3227 ! ----------------------------------------------------------------------
3228 BEXP = BB(SOILTYP)
3229 DKSAT = SATDK(SOILTYP)
3230 DWSAT = SATDW(SOILTYP)
3231 F1 = F11(SOILTYP)
3232 KDT = REFKDT * DKSAT/REFDK
3233 PSISAT = SATPSI(SOILTYP)
3234 QUARTZ = QTZ(SOILTYP)
3235 SMCDRY = DRYSMC(SOILTYP)
3236 SMCMAX = MAXSMC(SOILTYP)
3237 SMCREF = REFSMC(SOILTYP)
3238 SMCWLT = WLTSMC(SOILTYP)
3239 FRZFACT = (SMCMAX / SMCREF) * (0.412 / 0.468)
3240
3241 ! ----------------------------------------------------------------------
3242 ! TO ADJUST FRZK PARAMETER TO ACTUAL SOIL TYPE: FRZK * FRZFACT
3243 ! ----------------------------------------------------------------------
3244 FRZX = FRZK * FRZFACT
3245
3246 ! ----------------------------------------------------------------------
3247 ! SET-UP VEGETATION PARAMETERS
3248 ! ----------------------------------------------------------------------
3249 NROOT = NROOT_DATA(VEGTYP)
3250 SNUP = SNUPX(VEGTYP)
3251 RSMIN = RSMTBL(VEGTYP)
3252 RGL = RGLTBL(VEGTYP)
3253 HS = HSTBL(VEGTYP)
3254 Z0 = Z0_DATA(VEGTYP)
3255 LAI = LAI_DATA(VEGTYP)
3256 IF (VEGTYP .EQ. BARE) SHDFAC = 0.0
3257
3258 IF (NROOT .GT. NSOIL) THEN
3259 WRITE(*,*) 'Warning: too many root layers'
3260 STOP 333
3261 ENDIF
3262
3263 ! ----------------------------------------------------------------------
3264 ! CALCULATE ROOT DISTRIBUTION. PRESENT VERSION ASSUMES UNIFORM
3265 ! DISTRIBUTION BASED ON SOIL LAYER DEPTHS.
3266 ! ----------------------------------------------------------------------
3267 DO I = 1,NROOT
3268 RTDIS(I) = -SLDPTH(I)/ZSOIL(NROOT)
3269 END DO
3270
3271 ! ----------------------------------------------------------------------
3272 ! SET-UP SLOPE PARAMETER
3273 ! ----------------------------------------------------------------------
3274 SLOPE = SLOPE_DATA(SLOPETYP)
3275
3276 ! ----------------------------------------------------------------------
3277 ! END SUBROUTINE REDPRM
3278 ! ----------------------------------------------------------------------
3279 END SUBROUTINE REDPRM
3280
3281 SUBROUTINE ROSR12 (P,A,B,C,D,DELTA,NSOIL)
3282
3283 IMPLICIT NONE
3284
3285 ! ----------------------------------------------------------------------
3286 ! SUBROUTINE ROSR12
3287 ! ----------------------------------------------------------------------
3288 ! INVERT (SOLVE) THE TRI-DIAGONAL MATRIX PROBLEM SHOWN BELOW:
3289 ! ### ### ### ### ### ###
3290 ! #B(1), C(1), 0 , 0 , 0 , . . . , 0 # # # # #
3291 ! #A(2), B(2), C(2), 0 , 0 , . . . , 0 # # # # #
3292 ! # 0 , A(3), B(3), C(3), 0 , . . . , 0 # # # # D(3) #
3293 ! # 0 , 0 , A(4), B(4), C(4), . . . , 0 # # P(4) # # D(4) #
3294 ! # 0 , 0 , 0 , A(5), B(5), . . . , 0 # # P(5) # # D(5) #
3295 ! # . . # # . # = # . #
3296 ! # . . # # . # # . #
3297 ! # . . # # . # # . #
3298 ! # 0 , . . . , 0 , A(M-2), B(M-2), C(M-2), 0 # #P(M-2)# #D(M-2)#
3299 ! # 0 , . . . , 0 , 0 , A(M-1), B(M-1), C(M-1)# #P(M-1)# #D(M-1)#
3300 ! # 0 , . . . , 0 , 0 , 0 , A(M) , B(M) # # P(M) # # D(M) #
3301 ! ### ### ### ### ### ###
3302 ! ----------------------------------------------------------------------
3303 INTEGER K
3304 INTEGER KK
3305 INTEGER NSOIL
3306
3307 REAL A(NSOIL)
3308 REAL B(NSOIL)
3309 REAL C(NSOIL)
3310 REAL D(NSOIL)
3311 REAL DELTA(NSOIL)
3312 REAL P(NSOIL)
3313
3314 ! ----------------------------------------------------------------------
3315 ! INITIALIZE EQN COEF C FOR THE LOWEST SOIL LAYER
3316 ! ----------------------------------------------------------------------
3317 C(NSOIL) = 0.0
3318
3319 ! ----------------------------------------------------------------------
3320 ! SOLVE THE COEFS FOR THE 1ST SOIL LAYER
3321 ! ----------------------------------------------------------------------
3322 P(1) = -C(1) / B(1)
3323 DELTA(1) = D(1) / B(1)
3324
3325 ! ----------------------------------------------------------------------
3326 ! SOLVE THE COEFS FOR SOIL LAYERS 2 THRU NSOIL
3327 ! ----------------------------------------------------------------------
3328 DO K = 2,NSOIL
3329 P(K) = -C(K) * ( 1.0 / (B(K) + A (K) * P(K-1)) )
3330 DELTA(K) = (D(K)-A(K)*DELTA(K-1))*(1.0/(B(K)+A(K)*P(K-1)))
3331 END DO
3332
3333 ! ----------------------------------------------------------------------
3334 ! SET P TO DELTA FOR LOWEST SOIL LAYER
3335 ! ----------------------------------------------------------------------
3336 P(NSOIL) = DELTA(NSOIL)
3337
3338 ! ----------------------------------------------------------------------
3339 ! ADJUST P FOR SOIL LAYERS 2 THRU NSOIL
3340 ! ----------------------------------------------------------------------
3341 DO K = 2,NSOIL
3342 KK = NSOIL - K + 1
3343 P(KK) = P(KK) * P(KK+1) + DELTA(KK)
3344 END DO
3345
3346 ! ----------------------------------------------------------------------
3347 ! END SUBROUTINE ROSR12
3348 ! ----------------------------------------------------------------------
3349 END SUBROUTINE ROSR12
3350
3351 SUBROUTINE SHFLX (SSOIL,STC,SMC,SMCMAX,NSOIL,T1,DT,YY,ZZ1,ZSOIL, &
3352 & TBOT,ZBOT,SMCWLT,PSISAT,SH2O,BEXP,F1,DF1,ICE, &
3353 & QUARTZ,CSOIL)
3354
3355 IMPLICIT NONE
3356
3357 ! ----------------------------------------------------------------------
3358 ! SUBROUTINE SHFLX
3359 ! ----------------------------------------------------------------------
3360 ! UPDATE THE TEMPERATURE STATE OF THE SOIL COLUMN BASED ON THE THERMAL
3361 ! DIFFUSION EQUATION AND UPDATE THE FROZEN SOIL MOISTURE CONTENT BASED
3362 ! ON THE TEMPERATURE.
3363 ! ----------------------------------------------------------------------
3364 INTEGER NSOLD
3365 PARAMETER(NSOLD = 20)
3366
3367 INTEGER I
3368 INTEGER ICE
3369 INTEGER IFRZ
3370 INTEGER NSOIL
3371
3372 REAL AI(NSOLD)
3373 REAL BI(NSOLD)
3374 REAL CI(NSOLD)
3375
3376 REAL BEXP
3377 REAL CSOIL
3378 REAL DF1
3379 REAL DT
3380 REAL F1
3381 REAL PSISAT
3382 REAL QUARTZ
3383 REAL RHSTS(NSOLD)
3384 REAL SSOIL
3385 REAL SH2O(NSOIL)
3386 REAL SMC(NSOIL)
3387 REAL SMCMAX
3388 REAL SMCWLT
3389 REAL STC(NSOIL)
3390 REAL STCF(NSOLD)
3391 REAL T0
3392 REAL T1
3393 REAL TBOT
3394 REAL YY
3395 REAL ZBOT
3396 REAL ZSOIL(NSOIL)
3397 REAL ZZ1
3398
3399 PARAMETER(T0 = 273.15)
3400
3401 ! ----------------------------------------------------------------------
3402 ! HRT ROUTINE CALCS THE RIGHT HAND SIDE OF THE SOIL TEMP DIF EQN
3403 ! ----------------------------------------------------------------------
3404 IF (ICE.EQ.1) THEN
3405
3406 ! ----------------------------------------------------------------------
3407 ! SEA-ICE CASE
3408 ! ----------------------------------------------------------------------
3409 CALL HRTICE (RHSTS,STC,NSOIL,ZSOIL,YY,ZZ1,DF1,AI,BI,CI)
3410
3411 CALL HSTEP (STCF,STC,RHSTS,DT,NSOIL,AI,BI,CI)
3412
3413 ELSE
3414
3415 ! ----------------------------------------------------------------------
3416 ! LAND-MASS CASE
3417 ! ----------------------------------------------------------------------
3418 CALL HRT (RHSTS,STC,SMC,SMCMAX,NSOIL,ZSOIL,YY,ZZ1,TBOT, &
3419 & ZBOT,PSISAT,SH2O,DT, &
3420 & BEXP,F1,DF1,QUARTZ,CSOIL,AI,BI,CI)
3421
3422 CALL HSTEP (STCF,STC,RHSTS,DT,NSOIL,AI,BI,CI)
3423
3424 ENDIF
3425
3426 DO I = 1,NSOIL
3427 STC(I) = STCF(I)
3428 END DO
3429
3430 ! ----------------------------------------------------------------------
3431 ! IN THE NO SNOWPACK CASE (VIA ROUTINE NOPAC BRANCH,) UPDATE THE GRND
3432 ! (SKIN) TEMPERATURE HERE IN RESPONSE TO THE UPDATED SOIL TEMPERATURE
3433 ! PROFILE ABOVE. (NOTE: INSPECTION OF ROUTINE SNOPAC SHOWS THAT T1
3434 ! BELOW IS A DUMMY VARIABLE ONLY, AS SKIN TEMPERATURE IS UPDATED
3435 ! DIFFERENTLY IN ROUTINE SNOPAC)
3436 ! ----------------------------------------------------------------------
3437 T1 = (YY + (ZZ1 - 1.0) * STC(1)) / ZZ1
3438
3439 ! ----------------------------------------------------------------------
3440 ! CALCULATE SURFACE SOIL HEAT FLUX
3441 ! ----------------------------------------------------------------------
3442 SSOIL = DF1 * (STC(1) - T1) / (0.5 * ZSOIL(1))
3443
3444 ! ----------------------------------------------------------------------
3445 ! END SUBROUTINE SHFLX
3446 ! ----------------------------------------------------------------------
3447 END SUBROUTINE SHFLX
3448
3449 SUBROUTINE SMFLX (SMC,NSOIL,CMC,DT,PRCP1,ZSOIL, &
3450 & SH2O,SLOPE,KDT,FRZFACT, &
3451 & SMCMAX,BEXP,SMCWLT,DKSAT,DWSAT, &
3452 & SHDFAC,CMCMAX, &
3453 & RUNOFF1,RUNOFF2,RUNOFF3, &
3454 & EDIR1,EC1,ET1, &
3455 & DRIP)
3456
3457 IMPLICIT NONE
3458
3459 ! ----------------------------------------------------------------------
3460 ! SUBROUTINE SMFLX
3461 ! ----------------------------------------------------------------------
3462 ! CALCULATE SOIL MOISTURE FLUX. THE SOIL MOISTURE CONTENT (SMC - A PER
3463 ! UNIT VOLUME MEASUREMENT) IS A DEPENDENT VARIABLE THAT IS UPDATED WITH
3464 ! PROGNOSTIC EQNS. THE CANOPY MOISTURE CONTENT (CMC) IS ALSO UPDATED.
3465 ! FROZEN GROUND VERSION: NEW STATES ADDED: SH2O, AND FROZEN GROUND
3466 ! CORRECTION FACTOR, FRZFACT AND PARAMETER SLOPE.
3467 ! ----------------------------------------------------------------------
3468 INTEGER NSOLD
3469 PARAMETER(NSOLD = 20)
3470
3471 INTEGER I
3472 INTEGER K
3473 INTEGER NSOIL
3474
3475 REAL AI(NSOLD)
3476 REAL BI(NSOLD)
3477 REAL CI(NSOLD)
3478
3479 REAL BEXP
3480 REAL CMC
3481 REAL CMCMAX
3482 REAL DKSAT
3483 REAL DRIP
3484 REAL DT
3485 REAL DUMMY
3486 REAL DWSAT
3487 REAL EC1
3488 REAL EDIR1
3489 REAL ET1(NSOIL)
3490 REAL EXCESS
3491 REAL FRZFACT
3492 REAL KDT
3493 REAL PCPDRP
3494 REAL PRCP1
3495 REAL RHSCT
3496 REAL RHSTT(NSOLD)
3497 REAL RUNOFF1
3498 REAL RUNOFF2
3499 REAL RUNOFF3
3500 REAL SHDFAC
3501 REAL SMC(NSOIL)
3502 REAL SH2O(NSOIL)
3503 REAL SICE(NSOLD)
3504 REAL SH2OA(NSOLD)
3505 REAL SH2OFG(NSOLD)
3506 REAL SLOPE
3507 REAL SMCMAX
3508 REAL SMCWLT
3509 REAL TRHSCT
3510 REAL ZSOIL(NSOIL)
3511
3512 ! ----------------------------------------------------------------------
3513 ! EXECUTABLE CODE BEGINS HERE.
3514 ! ----------------------------------------------------------------------
3515 DUMMY = 0.
3516
3517 ! ----------------------------------------------------------------------
3518 ! COMPUTE THE RIGHT HAND SIDE OF THE CANOPY EQN TERM ( RHSCT )
3519 ! ----------------------------------------------------------------------
3520 RHSCT = SHDFAC * PRCP1 - EC1
3521
3522 ! ----------------------------------------------------------------------
3523 ! CONVERT RHSCT (A RATE) TO TRHSCT (AN AMOUNT) AND ADD IT TO EXISTING
3524 ! CMC. IF RESULTING AMT EXCEEDS MAX CAPACITY, IT BECOMES DRIP AND WILL
3525 ! FALL TO THE GRND.
3526 ! ----------------------------------------------------------------------
3527 DRIP = 0.
3528 TRHSCT = DT * RHSCT
3529 EXCESS = CMC + TRHSCT
3530 IF (EXCESS .GT. CMCMAX) DRIP = EXCESS - CMCMAX
3531
3532 ! ----------------------------------------------------------------------
3533 ! PCPDRP IS THE COMBINED PRCP1 AND DRIP (FROM CMC) THAT GOES INTO THE
3534 ! SOIL
3535 ! ----------------------------------------------------------------------
3536 PCPDRP = (1. - SHDFAC) * PRCP1 + DRIP / DT
3537
3538 ! ----------------------------------------------------------------------
3539 ! STORE ICE CONTENT AT EACH SOIL LAYER BEFORE CALLING SRT & SSTEP
3540 ! ----------------------------------------------------------------------
3541 DO I = 1,NSOIL
3542 SICE(I) = SMC(I) - SH2O(I)
3543 END DO
3544
3545 ! ----------------------------------------------------------------------
3546 ! CALL SUBROUTINES SRT AND SSTEP TO SOLVE THE SOIL MOISTURE
3547 ! TENDENCY EQUATIONS.
3548 !
3549 ! IF THE INFILTRATING PRECIP RATE IS NONTRIVIAL,
3550 ! (WE CONSIDER NONTRIVIAL TO BE A PRECIP TOTAL OVER THE TIME STEP
3551 ! EXCEEDING ONE ONE-THOUSANDTH OF THE WATER HOLDING CAPACITY OF
3552 ! THE FIRST SOIL LAYER)
3553 ! THEN CALL THE SRT/SSTEP SUBROUTINE PAIR TWICE IN THE MANNER OF
3554 ! TIME SCHEME "F" (IMPLICIT STATE, AVERAGED COEFFICIENT)
3555 ! OF SECTION 2 OF KALNAY AND KANAMITSU (1988, MWR, VOL 116,
3556 ! PAGES 1945-1958)TO MINIMIZE 2-DELTA-T OSCILLATIONS IN THE
3557 ! SOIL MOISTURE VALUE OF THE TOP SOIL LAYER THAT CAN ARISE BECAUSE
3558 ! OF THE EXTREME NONLINEAR DEPENDENCE OF THE SOIL HYDRAULIC
3559 ! DIFFUSIVITY COEFFICIENT AND THE HYDRAULIC CONDUCTIVITY ON THE
3560 ! SOIL MOISTURE STATE
3561 ! OTHERWISE CALL THE SRT/SSTEP SUBROUTINE PAIR ONCE IN THE MANNER OF
3562 ! TIME SCHEME "D" (IMPLICIT STATE, EXPLICIT COEFFICIENT)
3563 ! OF SECTION 2 OF KALNAY AND KANAMITSU
3564 ! PCPDRP IS UNITS OF KG/M**2/S OR MM/S, ZSOIL IS NEGATIVE DEPTH IN M
3565 ! ----------------------------------------------------------------------
3566 ! IF ( PCPDRP .GT. 0.0 ) THEN
3567 IF ( (PCPDRP*DT) .GT. (0.001*1000.0*(-ZSOIL(1))*SMCMAX) ) THEN
3568
3569 ! ----------------------------------------------------------------------
3570 ! FROZEN GROUND VERSION:
3571 ! SMC STATES REPLACED BY SH2O STATES IN SRT SUBR. SH2O & SICE STATES
3572 ! INCLUDED IN SSTEP SUBR. FROZEN GROUND CORRECTION FACTOR, FRZFACT
3573 ! ADDED. ALL WATER BALANCE CALCULATIONS USING UNFROZEN WATER
3574 ! ----------------------------------------------------------------------
3575 CALL SRT (RHSTT,EDIR1,ET1,SH2O,SH2O,NSOIL,PCPDRP,ZSOIL, &
3576 & DWSAT,DKSAT,SMCMAX,BEXP,RUNOFF1, &
3577 & RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZFACT,SICE,AI,BI,CI)
3578
3579 CALL SSTEP (SH2OFG,SH2O,DUMMY,RHSTT,RHSCT,DT,NSOIL,SMCMAX, &
3580 & CMCMAX,RUNOFF3,ZSOIL,SMC,SICE,AI,BI,CI)
3581
3582 DO K = 1,NSOIL
3583 SH2OA(K) = (SH2O(K) + SH2OFG(K)) * 0.5
3584 END DO
3585
3586 CALL SRT (RHSTT,EDIR1,ET1,SH2O,SH2OA,NSOIL,PCPDRP,ZSOIL, &
3587 & DWSAT,DKSAT,SMCMAX,BEXP,RUNOFF1, &
3588 & RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZFACT,SICE,AI,BI,CI)
3589
3590 CALL SSTEP (SH2O,SH2O,CMC,RHSTT,RHSCT,DT,NSOIL,SMCMAX, &
3591 & CMCMAX,RUNOFF3,ZSOIL,SMC,SICE,AI,BI,CI)
3592
3593 ELSE
3594
3595 CALL SRT (RHSTT,EDIR1,ET1,SH2O,SH2O,NSOIL,PCPDRP,ZSOIL, &
3596 & DWSAT,DKSAT,SMCMAX,BEXP,RUNOFF1, &
3597 & RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZFACT,SICE,AI,BI,CI)
3598
3599 CALL SSTEP (SH2O,SH2O,CMC,RHSTT,RHSCT,DT,NSOIL,SMCMAX, &
3600 & CMCMAX,RUNOFF3,ZSOIL,SMC,SICE,AI,BI,CI)
3601
3602 ENDIF
3603
3604 ! RUNOF = RUNOFF
3605
3606 ! ----------------------------------------------------------------------
3607 ! END SUBROUTINE SMFLX
3608 ! ----------------------------------------------------------------------
3609 END SUBROUTINE SMFLX
3610
3611 SUBROUTINE SNFRAC (SNEQV,SNUP,SALP,SNOWH,SNCOVR)
3612
3613 IMPLICIT NONE
3614
3615 ! ----------------------------------------------------------------------
3616 ! SUBROUTINE SNFRAC
3617 ! ----------------------------------------------------------------------
3618 ! CALCULATE SNOW FRACTION (0 -> 1)
3619 ! SNEQV SNOW WATER EQUIVALENT (M)
3620 ! SNUP THRESHOLD SNEQV DEPTH ABOVE WHICH SNCOVR=1
3621 ! SALP TUNING PARAMETER
3622 ! SNCOVR FRACTIONAL SNOW COVER
3623 ! ----------------------------------------------------------------------
3624 REAL SNEQV, SNUP, SALP, SNCOVR, RSNOW, Z0N, SNOWH
3625
3626 ! ----------------------------------------------------------------------
3627 ! SNUP IS VEG-CLASS DEPENDENT SNOWDEPTH THRESHHOLD (SET IN ROUTINE
3628 ! REDPRM) ABOVE WHICH SNOCVR=1.
3629 ! ----------------------------------------------------------------------
3630 IF (SNEQV .LT. SNUP) THEN
3631 RSNOW = SNEQV/SNUP
3632 SNCOVR = 1. - ( EXP(-SALP*RSNOW) - RSNOW*EXP(-SALP))
3633 ELSE
3634 SNCOVR = 1.0
3635 ENDIF
3636
3637 Z0N=0.035
3638 ! FORMULATION OF DICKINSON ET AL. 1986
3639
3640 ! SNCOVR=SNOWH/(SNOWH + 5*Z0N)
3641
3642 ! FORMULATION OF MARSHALL ET AL. 1994
3643 ! SNCOVR=SNEQV/(SNEQV + 2*Z0N)
3644
3645 ! ----------------------------------------------------------------------
3646 ! END SUBROUTINE SNFRAC
3647 ! ----------------------------------------------------------------------
3648 END SUBROUTINE SNFRAC
3649
3650 SUBROUTINE SNOPAC (ETP,ETA,PRCP,PRCP1,SNOWNG,SMC,SMCMAX,SMCWLT, &
3651 & SMCREF,SMCDRY,CMC,CMCMAX,NSOIL,DT, &
3652 & SBETA,DF1, &
3653 & Q2,T1,SFCTMP,T24,TH2,FDOWN,F1,SSOIL,STC,EPSCA, &
3654 & SFCPRS,BEXP,PC,RCH,RR,CFACTR,SNCOVR,ESD,SNDENS,&
3655 & SNOWH,SH2O,SLOPE,KDT,FRZFACT,PSISAT,SNUP, &
3656 & ZSOIL,DWSAT,DKSAT,TBOT,ZBOT,SHDFAC,RUNOFF1, &
3657 & RUNOFF2,RUNOFF3,EDIR,EC,ET,ETT,NROOT,SNOMLT, &
3658 & ICE,RTDIS,QUARTZ,FXEXP,CSOIL, &
3659 & BETA,DRIP,DEW,FLX1,FLX2,FLX3,ESNOW)
3660
3661 IMPLICIT NONE
3662
3663 ! ----------------------------------------------------------------------
3664 ! SUBROUTINE SNOPAC
3665 ! ----------------------------------------------------------------------
3666 ! CALCULATE SOIL MOISTURE AND HEAT FLUX VALUES & UPDATE SOIL MOISTURE
3667 ! CONTENT AND SOIL HEAT CONTENT VALUES FOR THE CASE WHEN A SNOW PACK IS
3668 ! PRESENT.
3669 ! ----------------------------------------------------------------------
3670 INTEGER ICE
3671 INTEGER NROOT
3672 INTEGER NSOIL
3673
3674 LOGICAL SNOWNG
3675
3676 REAL BEXP
3677 REAL BETA
3678 REAL CFACTR
3679 REAL CMC
3680 REAL CMCMAX
3681 REAL CP
3682 REAL CPH2O
3683 REAL CPICE
3684 REAL CSOIL
3685 REAL DENOM
3686 REAL DEW
3687 REAL DF1
3688 REAL DKSAT
3689 REAL DRIP
3690 REAL DSOIL
3691 REAL DTOT
3692 REAL DT
3693 REAL DWSAT
3694 REAL EC
3695 REAL EDIR
3696 REAL EPSCA
3697 REAL ESD
3698 REAL ESDMIN
3699 REAL EXPSNO
3700 REAL EXPSOI
3701 REAL ETA
3702 REAL ETA1
3703 REAL ETP
3704 REAL ETP1
3705 REAL ETP2
3706 REAL ET(NSOIL)
3707 REAL ETT
3708 REAL EX
3709 REAL EXPFAC
3710 REAL FDOWN
3711 REAL FXEXP
3712 REAL FLX1
3713 REAL FLX2
3714 REAL FLX3
3715 REAL F1
3716 REAL KDT
3717 REAL LSUBF
3718 REAL LSUBC
3719 REAL LSUBS
3720 REAL PC
3721 REAL PRCP
3722 REAL PRCP1
3723 REAL Q2
3724 REAL RCH
3725 REAL RR
3726 REAL RTDIS(NSOIL)
3727 REAL SSOIL
3728 REAL SBETA
3729 REAL SSOIL1
3730 REAL SFCTMP
3731 REAL SHDFAC
3732 REAL SIGMA
3733 REAL SMC(NSOIL)
3734 REAL SH2O(NSOIL)
3735 REAL SMCDRY
3736 REAL SMCMAX
3737 REAL SMCREF
3738 REAL SMCWLT
3739 REAL SNOMLT
3740 REAL SNOWH
3741 REAL STC(NSOIL)
3742 REAL T1
3743 REAL T11
3744 REAL T12
3745 REAL T12A
3746 REAL T12B
3747 REAL T24
3748 REAL TBOT
3749 REAL ZBOT
3750 REAL TH2
3751 REAL YY
3752 REAL ZSOIL(NSOIL)
3753 REAL ZZ1
3754 REAL TFREEZ
3755 REAL SALP
3756 REAL SFCPRS
3757 REAL SLOPE
3758 REAL FRZFACT
3759 REAL PSISAT
3760 REAL SNUP
3761 REAL RUNOFF1
3762 REAL RUNOFF2
3763 REAL RUNOFF3
3764 REAL QUARTZ
3765 REAL SNDENS
3766 REAL SNCOND
3767 REAL RSNOW
3768 REAL SNCOVR
3769 REAL QSAT
3770 REAL ETP3
3771 REAL SEH
3772 REAL T14
3773 ! REAL CSNOW
3774
3775 REAL EC1
3776 REAL EDIR1
3777 REAL ET1(NSOIL)
3778 REAL ETT1
3779
3780 REAL ETNS
3781 REAL ETNS1
3782 REAL ESNOW
3783 REAL ESNOW1
3784 REAL ESNOW2
3785 REAL ETANRG
3786
3787 INTEGER K
3788
3789 REAL SNOEXP
3790
3791 PARAMETER(CP = 1004.5)
3792 PARAMETER(CPH2O = 4.218E+3)
3793 PARAMETER(CPICE = 2.106E+3)
3794 PARAMETER(ESDMIN = 1.E-6)
3795 PARAMETER(LSUBF = 3.335E+5)
3796 PARAMETER(LSUBC = 2.501000E+6)
3797 PARAMETER(LSUBS = 2.83E+6)
3798 PARAMETER(SIGMA = 5.67E-8)
3799 PARAMETER(TFREEZ = 273.15)
3800
3801 ! DATA SNOEXP /1.0/
3802 DATA SNOEXP /2.0/
3803
3804 ! ----------------------------------------------------------------------
3805 ! EXECUTABLE CODE BEGINS HERE:
3806 ! CONVERT POTENTIAL EVAP (ETP) FROM KG M-2 S-1 TO M S-1 AND THEN TO AN
3807 ! AMOUNT (M) GIVEN TIMESTEP (DT) AND CALL IT AN EFFECTIVE SNOWPACK
3808 ! REDUCTION AMOUNT, ETP2 (M). THIS IS THE AMOUNT THE SNOWPACK WOULD BE
3809 ! REDUCED DUE TO EVAPORATION FROM THE SNOW SFC DURING THE TIMESTEP.
3810 ! EVAPORATION WILL PROCEED AT THE POTENTIAL RATE UNLESS THE SNOW DEPTH
3811 ! IS LESS THAN THE EXPECTED SNOWPACK REDUCTION.
3812 ! IF SEAICE (ICE=1), BETA REMAINS=1.
3813 ! ----------------------------------------------------------------------
3814 PRCP1 = PRCP1*0.001
3815
3816 ! ETP2 = ETP * 0.001 * DT
3817 BETA = 1.0
3818 IF (ICE .NE. 1) THEN
3819 IF (ESD .LT. ETP2) THEN
3820 ! BETA = ESD / ETP2
3821 ENDIF
3822 ENDIF
3823
3824 ! ----------------------------------------------------------------------
3825 EDIR = 0.0
3826 EDIR1 = 0.0
3827 EC = 0.0
3828 EC1 = 0.0
3829 DO K = 1,NSOIL
3830 ET(K) = 0.0
3831 ET1(K) = 0.0
3832 ENDDO
3833 ETT = 0.0
3834 ETT1 = 0.0
3835 ETNS = 0.0
3836 ETNS1 = 0.0
3837 ESNOW = 0.0
3838 ESNOW1 = 0.0
3839 ESNOW2 = 0.0
3840 ! ----------------------------------------------------------------------
3841
3842 ! ----------------------------------------------------------------------
3843 ! IF ETP<0 (DOWNWARD) THEN DEWFALL (=FROSTFALL IN THIS CASE).
3844 ! ----------------------------------------------------------------------
3845 DEW = 0.0
3846 ETP1 = ETP*0.001
3847 IF (ETP .LT. 0.0) THEN
3848 ! DEW = -ETP * 0.001
3849 DEW = -ETP1
3850 ! ESNOW2 = ETP * 0.001 * DT
3851 ESNOW2 = ETP1 * DT
3852 ETANRG = ETP*((1.-SNCOVR)*LSUBC + SNCOVR*LSUBS)
3853 ! ENDIF
3854 ELSE
3855 ! ----------------------------------------------------------------------
3856 ! ETP1 = 0.0
3857 ! ETP1 = ETP*0.001
3858 IF (ICE .NE. 1) THEN
3859 IF (SNCOVR .LT. 1.) THEN
3860 ! CALL EVAPO (ETA1,SMC,NSOIL,CMC,ETP1,DT,ZSOIL,
3861 CALL EVAPO (ETNS1,SMC,NSOIL,CMC,ETP1,DT,ZSOIL, &
3862 & SH2O,SMCMAX,BEXP,PC,SMCWLT,DKSAT,DWSAT, &
3863 & SMCREF,SHDFAC,CMCMAX, &
3864 & SMCDRY,CFACTR, &
3865 & EDIR1,EC1,ET1,ETT1,SFCTMP,Q2,NROOT,RTDIS,FXEXP)
3866 ! ENDIF
3867 ! ----------------------------------------------------------------------
3868 EDIR1 = EDIR1*(1.-SNCOVR)
3869 EC1 = EC1*(1.-SNCOVR)
3870 DO K = 1,NSOIL
3871 ET1(K) = ET1(K)*(1.-SNCOVR)
3872 END DO
3873 ETT1 = ETT1*(1.-SNCOVR)
3874 ETNS1 = ETNS1*(1.-SNCOVR)
3875 ! ----------------------------------------------------------------------
3876 EDIR = EDIR1 * 1000.0
3877 EC = EC1 * 1000.0
3878 DO K = 1,NSOIL
3879 ET(K) = ET1(K) * 1000.0
3880 END DO
3881 ETT = ETT1 * 1000.0
3882 ETNS = ETNS1 * 1000.0
3883 ! ----------------------------------------------------------------------
3884 ENDIF
3885 ESNOW = ETP*SNCOVR
3886 ! ESNOW1 = ETP*0.001
3887 ESNOW1 = ESNOW*0.001
3888 ESNOW2 = ESNOW1*DT
3889 ETANRG = ESNOW*LSUBS + ETNS*LSUBC
3890 ENDIF
3891 ENDIF
3892
3893 ! ----------------------------------------------------------------------
3894 ! IF PRECIP IS FALLING, CALCULATE HEAT FLUX FROM SNOW SFC TO NEWLY
3895 ! ACCUMULATING PRECIP. NOTE THAT THIS REFLECTS THE FLUX APPROPRIATE FOR
3896 ! THE NOT-YET-UPDATED SKIN TEMPERATURE (T1). ASSUMES TEMPERATURE OF THE
3897 ! SNOWFALL STRIKING THE GOUND IS =SFCTMP (LOWEST MODEL LEVEL AIR TEMP).
3898 ! ----------------------------------------------------------------------
3899 FLX1 = 0.0
3900 IF (SNOWNG) THEN
3901 FLX1 = CPICE * PRCP * (T1 - SFCTMP)
3902 ELSE
3903 IF (PRCP .GT. 0.0) FLX1 = CPH2O * PRCP * (T1 - SFCTMP)
3904 ENDIF
3905
3906 ! ----------------------------------------------------------------------
3907 ! CALCULATE AN 'EFFECTIVE SNOW-GRND SFC TEMP' (T12) BASED ON HEAT FLUXES
3908 ! BETWEEN THE SNOW PACK AND THE SOIL AND ON NET RADIATION.
3909 ! INCLUDE FLX1 (PRECIP-SNOW SFC) AND FLX2 (FREEZING RAIN LATENT HEAT)
3910 ! FLUXES.
3911 ! FLX2 REFLECTS FREEZING RAIN LATENT HEAT FLUX USING T1 CALCULATED IN
3912 ! PENMAN.
3913 ! ----------------------------------------------------------------------
3914 DSOIL = -(0.5 * ZSOIL(1))
3915 DTOT = SNOWH + DSOIL
3916 DENOM = 1.0 + DF1 / (DTOT * RR * RCH)
3917 ! T12A = ( (FDOWN-FLX1-FLX2-SIGMA*T24)/RCH
3918 ! & + TH2 - SFCTMP - BETA*EPSCA ) / RR
3919 ! T12A = ( (FDOWN-FLX1-FLX2-SIGMA*T24)/RCH
3920 ! M.Ek, 24Nov04, add snow emissivity
3921 T12A = ((FDOWN-FLX1-FLX2 &
3922 & -(0.95*SNCOVR+(1.0-SNCOVR))*SIGMA*T24)/RCH &
3923 & + TH2 - SFCTMP - ETANRG/RCH ) / RR
3924 T12B = DF1 * STC(1) / (DTOT * RR * RCH)
3925 T12 = (SFCTMP + T12A + T12B) / DENOM
3926
3927 ! ----------------------------------------------------------------------
3928 ! IF THE 'EFFECTIVE SNOW-GRND SFC TEMP' IS AT OR BELOW FREEZING, NO SNOW
3929 ! MELT WILL OCCUR. SET THE SKIN TEMP TO THIS EFFECTIVE TEMP. REDUCE
3930 ! (BY SUBLIMINATION ) OR INCREASE (BY FROST) THE DEPTH OF THE SNOWPACK,
3931 ! DEPENDING ON SIGN OF ETP.
3932 ! UPDATE SOIL HEAT FLUX (SSOIL) USING NEW SKIN TEMPERATURE (T1)
3933 ! SINCE NO SNOWMELT, SET ACCUMULATED SNOWMELT TO ZERO, SET 'EFFECTIVE'
3934 ! PRECIP FROM SNOWMELT TO ZERO, SET PHASE-CHANGE HEAT FLUX FROM SNOWMELT
3935 ! TO ZERO.
3936 ! ----------------------------------------------------------------------
3937 IF (T12 .LE. TFREEZ) THEN
3938 T1 = T12
3939 SSOIL = DF1 * (T1 - STC(1)) / DTOT
3940 ! ESD = MAX(0.0, ESD-ETP2)
3941 ESD = MAX(0.0, ESD-ESNOW2)
3942 FLX3 = 0.0
3943 EX = 0.0
3944 SNOMLT = 0.0
3945
3946 ELSE
3947 ! ----------------------------------------------------------------------
3948 ! IF THE 'EFFECTIVE SNOW-GRND SFC TEMP' IS ABOVE FREEZING, SNOW MELT
3949 ! WILL OCCUR. CALL THE SNOW MELT RATE,EX AND AMT, SNOMLT. REVISE THE
3950 ! EFFECTIVE SNOW DEPTH. REVISE THE SKIN TEMP BECAUSE IT WOULD HAVE CHGD
3951 ! DUE TO THE LATENT HEAT RELEASED BY THE MELTING. CALC THE LATENT HEAT
3952 ! RELEASED, FLX3. SET THE EFFECTIVE PRECIP, PRCP1 TO THE SNOW MELT RATE,
3953 ! EX FOR USE IN SMFLX. ADJUSTMENT TO T1 TO ACCOUNT FOR SNOW PATCHES.
3954 ! CALCULATE QSAT VALID AT FREEZING POINT. NOTE THAT ESAT (SATURATION
3955 ! VAPOR PRESSURE) VALUE OF 6.11E+2 USED HERE IS THAT VALID AT FRZZING
3956 ! POINT. NOTE THAT ETP FROM CALL PENMAN IN SFLX IS IGNORED HERE IN
3957 ! FAVOR OF BULK ETP OVER 'OPEN WATER' AT FREEZING TEMP.
3958 ! UPDATE SOIL HEAT FLUX (S) USING NEW SKIN TEMPERATURE (T1)
3959 ! ----------------------------------------------------------------------
3960 ! T1 = TFREEZ * SNCOVR + T12 * (1.0 - SNCOVR)
3961 ! mek Feb2004
3962 ! non-linear weighting of snow vs non-snow covered portions of gridbox
3963 ! so with SNOEXP = 2.0 (>1), surface skin temperature is higher than for
3964 ! the linear case (SNOEXP = 1).
3965 T1 = TFREEZ * SNCOVR**SNOEXP + T12 * (1.0 - SNCOVR**SNOEXP)
3966 ! QSAT = (0.622*6.11E2)/(SFCPRS-0.378*6.11E2)
3967 ! ETP = RCH*(QSAT-Q2)/CP
3968 ! ETP2 = ETP*0.001*DT
3969 BETA = 1.0
3970 SSOIL = DF1 * (T1 - STC(1)) / DTOT
3971
3972 ! ----------------------------------------------------------------------
3973 ! IF POTENTIAL EVAP (SUBLIMATION) GREATER THAN DEPTH OF SNOWPACK.
3974 ! BETA<1
3975 ! SNOWPACK HAS SUBLIMATED AWAY, SET DEPTH TO ZERO.
3976 ! ----------------------------------------------------------------------
3977 ! IF (ESD .LE. ETP2) THEN
3978 ! IF (ESD .LE. ESNOW2) THEN
3979 IF (ESD-ESNOW2 .LE. ESDMIN) THEN
3980 ! BETA = ESD / ETP2
3981 ESD = 0.0
3982 EX = 0.0
3983 SNOMLT = 0.0
3984
3985 ELSE
3986 ! ----------------------------------------------------------------------
3987 ! POTENTIAL EVAP (SUBLIMATION) LESS THAN DEPTH OF SNOWPACK, RETAIN
3988 ! BETA=1.
3989 ! SNOWPACK (ESD) REDUCED BY POTENTIAL EVAP RATE
3990 ! ETP3 (CONVERT TO FLUX)
3991 ! ----------------------------------------------------------------------
3992 ! ESD = ESD-ETP2
3993 ESD = ESD-ESNOW2
3994 ! ETP3 = ETP*LSUBC
3995 SEH = RCH*(T1-TH2)
3996 T14 = T1*T1
3997 T14 = T14*T14
3998 ! FLX3 = FDOWN - FLX1 - FLX2 - SIGMA*T14 - SSOIL - SEH - ETP3
3999 ! FLX3 = FDOWN - FLX1 - FLX2 - SIGMA*T14 - SSOIL - SEH - ETANRG
4000 ! M.Ek, 24Nov04, add snow emissivity
4001 FLX3 = FDOWN - FLX1 - FLX2 - &
4002 & (0.95*SNCOVR+(1.0-SNCOVR))*SIGMA*T14 - SSOIL - SEH - ETANRG
4003 IF (FLX3 .LE. 0.0) FLX3 = 0.0
4004 EX = FLX3*0.001/LSUBF
4005
4006 ! ----------------------------------------------------------------------
4007 ! SNOWMELT REDUCTION DEPENDING ON SNOW COVER
4008 ! IF SNOW COVER LESS THAN 5% NO SNOWMELT REDUCTION
4009 ! ***NOTE: DOES 'IF' BELOW FAIL TO MATCH THE MELT WATER WITH THE MELT
4010 ! ENERGY?
4011 ! ----------------------------------------------------------------------
4012 ! IF (SNCOVR .GT. 0.05) EX = EX * SNCOVR
4013 SNOMLT = EX * DT
4014
4015 ! ----------------------------------------------------------------------
4016 ! ESDMIN REPRESENTS A SNOWPACK DEPTH THRESHOLD VALUE BELOW WHICH WE
4017 ! CHOOSE NOT TO RETAIN ANY SNOWPACK, AND INSTEAD INCLUDE IT IN SNOWMELT.
4018 ! ----------------------------------------------------------------------
4019 IF (ESD-SNOMLT .GE. ESDMIN) THEN
4020 ESD = ESD - SNOMLT
4021
4022 ELSE
4023 ! ----------------------------------------------------------------------
4024 ! SNOWMELT EXCEEDS SNOW DEPTH
4025 ! ----------------------------------------------------------------------
4026 EX = ESD/DT
4027 FLX3 = EX*1000.0*LSUBF
4028 SNOMLT = ESD
4029 ESD = 0.0
4030
4031 ENDIF
4032 ! ----------------------------------------------------------------------
4033 ! END OF 'ESD .LE. ETP2' IF-BLOCK
4034 ! ----------------------------------------------------------------------
4035 ENDIF
4036
4037 PRCP1 = PRCP1 + EX
4038
4039 ! ----------------------------------------------------------------------
4040 ! END OF 'T12 .LE. TFREEZ' IF-BLOCK
4041 ! ----------------------------------------------------------------------
4042 ENDIF
4043
4044 ! ----------------------------------------------------------------------
4045 ! FINAL BETA NOW IN HAND, SO COMPUTE EVAPORATION. EVAP EQUALS ETP
4046 ! UNLESS BETA<1.
4047 ! ----------------------------------------------------------------------
4048 ! ETA = BETA*ETP
4049
4050 ! ----------------------------------------------------------------------
4051 ! SET THE EFFECTIVE POTNL EVAPOTRANSP (ETP1) TO ZERO SINCE THIS IS SNOW
4052 ! CASE, SO SURFACE EVAP NOT CALCULATED FROM EDIR, EC, OR ETT IN SMFLX
4053 ! (BELOW).
4054 ! IF SEAICE (ICE=1) SKIP CALL TO SMFLX.
4055 ! SMFLX RETURNS UPDATED SOIL MOISTURE VALUES. IN THIS, THE SNOW PACK
4056 ! CASE, ETA1 IS NOT USED IN CALCULATION OF EVAP.
4057 ! ----------------------------------------------------------------------
4058 ! ETP1 = 0.0
4059 IF (ICE .NE. 1) THEN
4060 ! CALL EVAPO (ETA1,SMC,NSOIL,CMC,ETP1,DT,ZSOIL,
4061 ! & SH2O,SMCMAX,BEXP,PC,SMCWLT,DKSAT,DWSAT,
4062 ! & SMCREF,SHDFAC,CMCMAX,
4063 ! & SMCDRY,CFACTR,
4064 ! & EDIR1,EC1,ET1,ETT1,SFCTMP,Q2,NROOT,RTDIS,FXEXP)
4065 CALL SMFLX (SMC,NSOIL,CMC,DT,PRCP1,ZSOIL, &
4066 & SH2O,SLOPE,KDT,FRZFACT, &
4067 & SMCMAX,BEXP,SMCWLT,DKSAT,DWSAT, &
4068 & SHDFAC,CMCMAX, &
4069 & RUNOFF1,RUNOFF2,RUNOFF3, &
4070 & EDIR1,EC1,ET1, &
4071 & DRIP)
4072
4073 ENDIF
4074
4075 ! ----------------------------------------------------------------------
4076 ! EDIR = EDIR1 * 1000.0
4077 ! EC = EC1 * 1000.0
4078 ! ETT = ETT1 * 1000.0
4079 ! ET(1) = ET1(1) * 1000.0
4080 ! ET(2) = ET1(2) * 1000.0
4081 ! ET(3) = ET1(3) * 1000.0
4082 ! ET(4) = ET1(4) * 1000.0
4083 ! ----------------------------------------------------------------------
4084
4085 ! ----------------------------------------------------------------------
4086 ! BEFORE CALL SHFLX IN THIS SNOWPACK CASE, SET ZZ1 AND YY ARGUMENTS TO
4087 ! SPECIAL VALUES THAT ENSURE THAT GROUND HEAT FLUX CALCULATED IN SHFLX
4088 ! MATCHES THAT ALREADY COMPUTER FOR BELOW THE SNOWPACK, THUS THE SFC
4089 ! HEAT FLUX TO BE COMPUTED IN SHFLX WILL EFFECTIVELY BE THE FLUX AT THE
4090 ! SNOW TOP SURFACE. T11 IS A DUMMY ARGUEMENT SO WE WILL NOT USE THE
4091 ! SKIN TEMP VALUE AS REVISED BY SHFLX.
4092 ! ----------------------------------------------------------------------
4093 ZZ1 = 1.0
4094 YY = STC(1)-0.5*SSOIL*ZSOIL(1)*ZZ1/DF1
4095 T11 = T1
4096
4097 ! ----------------------------------------------------------------------
4098 ! SHFLX WILL CALC/UPDATE THE SOIL TEMPS. NOTE: THE SUB-SFC HEAT FLUX
4099 ! (SSOIL1) AND THE SKIN TEMP (T11) OUTPUT FROM THIS SHFLX CALL ARE NOT
4100 ! USED IN ANY SUBSEQUENT CALCULATIONS. RATHER, THEY ARE DUMMY VARIABLES
4101 ! HERE IN THE SNOPAC CASE, SINCE THE SKIN TEMP AND SUB-SFC HEAT FLUX ARE
4102 ! UPDATED INSTEAD NEAR THE BEGINNING OF THE CALL TO SNOPAC.
4103 ! ----------------------------------------------------------------------
4104 CALL SHFLX (SSOIL1,STC,SMC,SMCMAX,NSOIL,T11,DT,YY,ZZ1,ZSOIL, &
4105 & TBOT,ZBOT,SMCWLT,PSISAT,SH2O,BEXP,F1,DF1,ICE, &
4106 & QUARTZ,CSOIL)
4107
4108 ! ----------------------------------------------------------------------
4109 ! SNOW DEPTH AND DENSITY ADJUSTMENT BASED ON SNOW COMPACTION. YY IS
4110 ! ASSUMED TO BE THE SOIL TEMPERTURE AT THE TOP OF THE SOIL COLUMN.
4111 ! ----------------------------------------------------------------------
4112 IF (ESD .GT. 0.) THEN
4113 CALL SNOWPACK (ESD,DT,SNOWH,SNDENS,T1,YY)
4114 ELSE
4115 ESD = 0.
4116 SNOWH = 0.
4117 SNDENS = 0.
4118 SNCOND = 1.
4119 SNCOVR = 0.
4120 ENDIF
4121
4122 ! ----------------------------------------------------------------------
4123 ! END SUBROUTINE SNOPAC
4124 ! ----------------------------------------------------------------------
4125 END SUBROUTINE SNOPAC
4126
4127 SUBROUTINE SNOWPACK (ESD,DTSEC,SNOWH,SNDENS,TSNOW,TSOIL)
4128
4129 IMPLICIT NONE
4130
4131 ! ----------------------------------------------------------------------
4132 ! SUBROUTINE SNOWPACK
4133 ! ----------------------------------------------------------------------
4134 ! CALCULATE COMPACTION OF SNOWPACK UNDER CONDITIONS OF INCREASING SNOW
4135 ! DENSITY, AS OBTAINED FROM AN APPROXIMATE SOLUTION OF E. ANDERSON'S
4136 ! DIFFERENTIAL EQUATION (3.29), NOAA TECHNICAL REPORT NWS 19, BY VICTOR
4137 ! KOREN, 03/25/95.
4138 ! ----------------------------------------------------------------------
4139 ! ESD WATER EQUIVALENT OF SNOW (M)
4140 ! DTSEC TIME STEP (SEC)
4141 ! SNOWH SNOW DEPTH (M)
4142 ! SNDENS SNOW DENSITY (G/CM3=DIMENSIONLESS FRACTION OF H2O DENSITY)
4143 ! TSNOW SNOW SURFACE TEMPERATURE (K)
4144 ! TSOIL SOIL SURFACE TEMPERATURE (K)
4145 !
4146 ! SUBROUTINE WILL RETURN NEW VALUES OF SNOWH AND SNDENS
4147 ! ----------------------------------------------------------------------
4148 INTEGER IPOL, J
4149
4150 REAL BFAC,C1,C2,SNDENS,DSX,DTHR,DTSEC,DW,SNOWHC,SNOWH,PEXP,TAVGC, &
4151 & TSNOW,TSNOWC,TSOIL,TSOILC,ESD,ESDC,ESDCX,G,KN
4152
4153 PARAMETER(C1 = 0.01, C2=21.0, G=9.81, KN=4000.0)
4154
4155 ! ----------------------------------------------------------------------
4156 ! CONVERSION INTO SIMULATION UNITS
4157 ! ----------------------------------------------------------------------
4158 SNOWHC = SNOWH*100.
4159 ESDC = ESD*100.
4160 DTHR = DTSEC/3600.
4161 TSNOWC = TSNOW-273.15
4162 TSOILC = TSOIL-273.15
4163
4164 ! ----------------------------------------------------------------------
4165 ! CALCULATING OF AVERAGE TEMPERATURE OF SNOW PACK
4166 ! ----------------------------------------------------------------------
4167 TAVGC = 0.5*(TSNOWC+TSOILC)
4168
4169 ! ----------------------------------------------------------------------
4170 ! CALCULATING OF SNOW DEPTH AND DENSITY AS A RESULT OF COMPACTION
4171 ! SNDENS=DS0*(EXP(BFAC*ESD)-1.)/(BFAC*ESD)
4172 ! BFAC=DTHR*C1*EXP(0.08*TAVGC-C2*DS0)
4173 ! NOTE: BFAC*ESD IN SNDENS EQN ABOVE HAS TO BE CAREFULLY TREATED
4174 ! NUMERICALLY BELOW:
4175 ! C1 IS THE FRACTIONAL INCREASE IN DENSITY (1/(CM*HR))
4176 ! C2 IS A CONSTANT (CM3/G) KOJIMA ESTIMATED AS 21 CMS/G
4177 ! ----------------------------------------------------------------------
4178 IF (ESDC .GT. 1.E-2) THEN
4179 ESDCX = ESDC
4180 ELSE
4181 ESDCX = 1.E-2
4182 ENDIF
4183 BFAC = DTHR*C1*EXP(0.08*TAVGC-C2*SNDENS)
4184
4185 ! DSX = SNDENS*((DEXP(BFAC*ESDC)-1.)/(BFAC*ESDC))
4186 ! ----------------------------------------------------------------------
4187 ! THE FUNCTION OF THE FORM (e**x-1)/x IMBEDDED IN ABOVE EXPRESSION
4188 ! FOR DSX WAS CAUSING NUMERICAL DIFFICULTIES WHEN THE DENOMINATOR "x"
4189 ! (I.E. BFAC*ESDC) BECAME ZERO OR APPROACHED ZERO (DESPITE THE FACT THAT
4190 ! THE ANALYTICAL FUNCTION (e**x-1)/x HAS A WELL DEFINED LIMIT AS
4191 ! "x" APPROACHES ZERO), HENCE BELOW WE REPLACE THE (e**x-1)/x
4192 ! EXPRESSION WITH AN EQUIVALENT, NUMERICALLY WELL-BEHAVED
4193 ! POLYNOMIAL EXPANSION.
4194 !
4195 ! NUMBER OF TERMS OF POLYNOMIAL EXPANSION, AND HENCE ITS ACCURACY,
4196 ! IS GOVERNED BY ITERATION LIMIT "IPOL".
4197 ! IPOL GREATER THAN 9 ONLY MAKES A DIFFERENCE ON DOUBLE
4198 ! PRECISION (RELATIVE ERRORS GIVEN IN PERCENT %).
4199 ! IPOL=9, FOR REL.ERROR <~ 1.6 E-6 % (8 SIGNIFICANT DIGITS)
4200 ! IPOL=8, FOR REL.ERROR <~ 1.8 E-5 % (7 SIGNIFICANT DIGITS)
4201 ! IPOL=7, FOR REL.ERROR <~ 1.8 E-4 % ...
4202 ! ----------------------------------------------------------------------
4203 IPOL = 4
4204 PEXP = 0.
4205 DO J = IPOL,1,-1
4206 ! PEXP = (1. + PEXP)*BFAC*ESDC/REAL(J+1)
4207 PEXP = (1. + PEXP)*BFAC*ESDCX/REAL(J+1)
4208 END DO
4209 PEXP = PEXP + 1.
4210
4211 DSX = SNDENS*(PEXP)
4212 ! ----------------------------------------------------------------------
4213 ! ABOVE LINE ENDS POLYNOMIAL SUBSTITUTION
4214 ! ----------------------------------------------------------------------
4215 ! END OF KOREAN FORMULATION
4216
4217 ! BASE FORMULATION (COGLEY ET AL., 1990)
4218 ! CONVERT DENSITY FROM G/CM3 TO KG/M3
4219 ! DSM=SNDENS*1000.0
4220
4221 ! DSX=DSM+DTSEC*0.5*DSM*G*ESD/
4222 ! & (1E7*EXP(-0.02*DSM+KN/(TAVGC+273.16)-14.643))
4223
4224 ! CONVERT DENSITY FROM KG/M3 TO G/CM3
4225 ! DSX=DSX/1000.0
4226
4227 ! END OF COGLEY ET AL. FORMULATION
4228
4229 ! ----------------------------------------------------------------------
4230 ! SET UPPER/LOWER LIMIT ON SNOW DENSITY
4231 ! ----------------------------------------------------------------------
4232 IF (DSX .GT. 0.40) DSX = 0.40
4233 IF (DSX .LT. 0.05) DSX = 0.05
4234 SNDENS = DSX
4235 ! ----------------------------------------------------------------------
4236 ! UPDATE OF SNOW DEPTH AND DENSITY DEPENDING ON LIQUID WATER DURING
4237 ! SNOWMELT. ASSUMED THAT 13% OF LIQUID WATER CAN BE STORED IN SNOW PER
4238 ! DAY DURING SNOWMELT TILL SNOW DENSITY 0.40.
4239 ! ----------------------------------------------------------------------
4240 IF (TSNOWC .GE. 0.) THEN
4241 DW = 0.13*DTHR/24.
4242 SNDENS = SNDENS*(1.-DW)+DW
4243 IF (SNDENS .GT. 0.40) SNDENS = 0.40
4244 ENDIF
4245
4246 ! ----------------------------------------------------------------------
4247 ! CALCULATE SNOW DEPTH (CM) FROM SNOW WATER EQUIVALENT AND SNOW DENSITY.
4248 ! CHANGE SNOW DEPTH UNITS TO METERS
4249 ! ----------------------------------------------------------------------
4250 SNOWHC = ESDC/SNDENS
4251 SNOWH = SNOWHC*0.01
4252
4253 ! ----------------------------------------------------------------------
4254 ! END SUBROUTINE SNOWPACK
4255 ! ----------------------------------------------------------------------
4256 END SUBROUTINE SNOWPACK
4257
4258 SUBROUTINE SNOWZ0 (SNCOVR,Z0)
4259
4260 IMPLICIT NONE
4261
4262 ! ----------------------------------------------------------------------
4263 ! SUBROUTINE SNOWZ0
4264 ! ----------------------------------------------------------------------
4265 ! CALCULATE TOTAL ROUGHNESS LENGTH OVER SNOW
4266 ! SNCOVR FRACTIONAL SNOW COVER
4267 ! Z0 ROUGHNESS LENGTH (m)
4268 ! Z0S SNOW ROUGHNESS LENGTH:=0.001 (m)
4269 ! ----------------------------------------------------------------------
4270 REAL SNCOVR, Z0, Z0S
4271 ! PARAMETER (Z0S=0.001)
4272
4273 ! CURRENT NOAH LSM CONDITION - MBEK, 09-OCT-2001
4274 Z0S = Z0
4275 !
4276 Z0 = (1-SNCOVR)*Z0 + SNCOVR*Z0S
4277 ! ----------------------------------------------------------------------
4278 ! END SUBROUTINE SNOWZ0
4279 ! ----------------------------------------------------------------------
4280 END SUBROUTINE SNOWZ0
4281
4282 SUBROUTINE SNOW_NEW (TEMP,NEWSN,SNOWH,SNDENS)
4283
4284 IMPLICIT NONE
4285
4286 ! ----------------------------------------------------------------------
4287 ! SUBROUTINE SNOW_NEW
4288 ! ----------------------------------------------------------------------
4289 ! CALCULATE SNOW DEPTH AND DENSITITY TO ACCOUNT FOR THE NEW SNOWFALL.
4290 ! NEW VALUES OF SNOW DEPTH & DENSITY RETURNED.
4291 !
4292 ! TEMP AIR TEMPERATURE (K)
4293 ! NEWSN NEW SNOWFALL (M)
4294 ! SNOWH SNOW DEPTH (M)
4295 ! SNDENS SNOW DENSITY (G/CM3=DIMENSIONLESS FRACTION OF H2O DENSITY)
4296 ! ----------------------------------------------------------------------
4297 REAL SNDENS
4298 REAL DSNEW
4299 REAL SNOWHC
4300 REAL HNEWC
4301 REAL SNOWH
4302 REAL NEWSN
4303 REAL NEWSNC
4304 REAL TEMP
4305 REAL TEMPC
4306
4307 ! ----------------------------------------------------------------------
4308 ! CONVERSION INTO SIMULATION UNITS
4309 ! ----------------------------------------------------------------------
4310 SNOWHC = SNOWH*100.
4311 NEWSNC = NEWSN*100.
4312 TEMPC = TEMP-273.15
4313
4314 ! ----------------------------------------------------------------------
4315 ! CALCULATING NEW SNOWFALL DENSITY DEPENDING ON TEMPERATURE
4316 ! EQUATION FROM GOTTLIB L. 'A GENERAL RUNOFF MODEL FOR SNOWCOVERED
4317 ! AND GLACIERIZED BASIN', 6TH NORDIC HYDROLOGICAL CONFERENCE,
4318 ! VEMADOLEN, SWEDEN, 1980, 172-177PP.
4319 !-----------------------------------------------------------------------
4320 IF (TEMPC .LE. -15.) THEN
4321 DSNEW = 0.05
4322 ELSE
4323 DSNEW = 0.05+0.0017*(TEMPC+15.)**1.5
4324 ENDIF
4325
4326 ! ----------------------------------------------------------------------
4327 ! ADJUSTMENT OF SNOW DENSITY DEPENDING ON NEW SNOWFALL
4328 ! ----------------------------------------------------------------------
4329 HNEWC = NEWSNC/DSNEW
4330 SNDENS = (SNOWHC*SNDENS+HNEWC*DSNEW)/(SNOWHC+HNEWC)
4331 SNOWHC = SNOWHC+HNEWC
4332 SNOWH = SNOWHC*0.01
4333
4334 ! ----------------------------------------------------------------------
4335 ! END SUBROUTINE SNOW_NEW
4336 ! ----------------------------------------------------------------------
4337 END SUBROUTINE SNOW_NEW
4338
4339 SUBROUTINE SRT (RHSTT,EDIR,ET,SH2O,SH2OA,NSOIL,PCPDRP, &
4340 & ZSOIL,DWSAT,DKSAT,SMCMAX,BEXP,RUNOFF1, &
4341 & RUNOFF2,DT,SMCWLT,SLOPE,KDT,FRZX,SICE,AI,BI,CI)
4342
4343 IMPLICIT NONE
4344
4345 ! ----------------------------------------------------------------------
4346 ! SUBROUTINE SRT
4347 ! ----------------------------------------------------------------------
4348 ! CALCULATE THE RIGHT HAND SIDE OF THE TIME TENDENCY TERM OF THE SOIL
4349 ! WATER DIFFUSION EQUATION. ALSO TO COMPUTE ( PREPARE ) THE MATRIX
4350 ! COEFFICIENTS FOR THE TRI-DIAGONAL MATRIX OF THE IMPLICIT TIME SCHEME.
4351 ! ----------------------------------------------------------------------
4352 INTEGER NSOLD
4353 PARAMETER(NSOLD = 20)
4354
4355 INTEGER CVFRZ
4356 INTEGER IALP1
4357 INTEGER IOHINF
4358 INTEGER J
4359 INTEGER JJ
4360 INTEGER K
4361 INTEGER KS
4362 INTEGER NSOIL
4363
4364 REAL ACRT
4365 REAL AI(NSOLD)
4366 REAL BEXP
4367 REAL BI(NSOLD)
4368 REAL CI(NSOLD)
4369 REAL DD
4370 REAL DDT
4371 REAL DDZ
4372 REAL DDZ2
4373 REAL DENOM
4374 REAL DENOM2
4375 REAL DICE
4376 REAL DKSAT
4377 REAL DMAX(NSOLD)
4378 REAL DSMDZ
4379 REAL DSMDZ2
4380 REAL DT
4381 REAL DT1
4382 REAL DWSAT
4383 REAL EDIR
4384 REAL ET(NSOIL)
4385 REAL FCR
4386 REAL FRZX
4387 REAL INFMAX
4388 REAL KDT
4389 REAL MXSMC
4390 REAL MXSMC2
4391 REAL NUMER
4392 REAL PCPDRP
4393 REAL PDDUM
4394 REAL PX
4395 REAL RHSTT(NSOIL)
4396 REAL RUNOFF1
4397 REAL RUNOFF2
4398 REAL SH2O(NSOIL)
4399 REAL SH2OA(NSOIL)
4400 REAL SICE(NSOIL)
4401 REAL SICEMAX
4402 REAL SLOPE
4403 REAL SLOPX
4404 REAL SMCAV
4405 REAL SMCMAX
4406 REAL SMCWLT
4407 REAL SSTT
4408 REAL SUM
4409 REAL VAL
4410 REAL WCND
4411 REAL WCND2
4412 REAL WDF
4413 REAL WDF2
4414 REAL ZSOIL(NSOIL)
4415
4416 ! ----------------------------------------------------------------------
4417 ! FROZEN GROUND VERSION:
4418 ! REFERENCE FROZEN GROUND PARAMETER, CVFRZ, IS A SHAPE PARAMETER OF
4419 ! AREAL DISTRIBUTION FUNCTION OF SOIL ICE CONTENT WHICH EQUALS 1/CV.
4420 ! CV IS A COEFFICIENT OF SPATIAL VARIATION OF SOIL ICE CONTENT. BASED
4421 ! ON FIELD DATA CV DEPENDS ON AREAL MEAN OF FROZEN DEPTH, AND IT CLOSE
4422 ! TO CONSTANT = 0.6 IF AREAL MEAN FROZEN DEPTH IS ABOVE 20 CM. THAT IS
4423 ! WHY PARAMETER CVFRZ = 3 (INT{1/0.6*0.6}).
4424 ! CURRENT LOGIC DOESN'T ALLOW CVFRZ BE BIGGER THAN 3
4425 ! ----------------------------------------------------------------------
4426 PARAMETER(CVFRZ = 3)
4427
4428 ! ----------------------------------------------------------------------
4429 ! DETERMINE RAINFALL INFILTRATION RATE AND RUNOFF. INCLUDE THE
4430 ! INFILTRATION FORMULE FROM SCHAAKE AND KOREN MODEL.
4431 ! MODIFIED BY Q DUAN
4432 ! ----------------------------------------------------------------------
4433 IOHINF=1
4434
4435 ! ----------------------------------------------------------------------
4436 ! LET SICEMAX BE THE GREATEST, IF ANY, FROZEN WATER CONTENT WITHIN SOIL
4437 ! LAYERS.
4438 ! ----------------------------------------------------------------------
4439 SICEMAX = 0.0
4440 DO KS=1,NSOIL
4441 IF (SICE(KS) .GT. SICEMAX) SICEMAX = SICE(KS)
4442 END DO
4443
4444 ! ----------------------------------------------------------------------
4445 ! DETERMINE RAINFALL INFILTRATION RATE AND RUNOFF
4446 ! ----------------------------------------------------------------------
4447 PDDUM = PCPDRP
4448 RUNOFF1 = 0.0
4449 IF (PCPDRP .NE. 0.0) THEN
4450
4451 ! ----------------------------------------------------------------------
4452 ! MODIFIED BY Q. DUAN, 5/16/94
4453 ! ----------------------------------------------------------------------
4454 ! IF (IOHINF .EQ. 1) THEN
4455
4456 DT1 = DT/86400.
4457 SMCAV = SMCMAX - SMCWLT
4458 DMAX(1)=-ZSOIL(1)*SMCAV
4459
4460 ! ----------------------------------------------------------------------
4461 ! FROZEN GROUND VERSION:
4462 ! ----------------------------------------------------------------------
4463 DICE = -ZSOIL(1) * SICE(1)
4464
4465 DMAX(1)=DMAX(1)*(1.0 - (SH2OA(1)+SICE(1)-SMCWLT)/SMCAV)
4466 DD=DMAX(1)
4467
4468 DO KS=2,NSOIL
4469
4470 ! ----------------------------------------------------------------------
4471 ! FROZEN GROUND VERSION:
4472 ! ----------------------------------------------------------------------
4473 DICE = DICE + ( ZSOIL(KS-1) - ZSOIL(KS) ) * SICE(KS)
4474
4475 DMAX(KS) = (ZSOIL(KS-1)-ZSOIL(KS))*SMCAV
4476 DMAX(KS) = DMAX(KS)*(1.0 - (SH2OA(KS)+SICE(KS)-SMCWLT)/SMCAV)
4477 DD = DD+DMAX(KS)
4478 END DO
4479
4480 ! ----------------------------------------------------------------------
4481 ! VAL = (1.-EXP(-KDT*SQRT(DT1)))
4482 ! IN BELOW, REMOVE THE SQRT IN ABOVE
4483 ! ----------------------------------------------------------------------
4484 VAL = (1.-EXP(-KDT*DT1))
4485 DDT = DD*VAL
4486 PX = PCPDRP*DT
4487 IF (PX .LT. 0.0) PX = 0.0
4488 INFMAX = (PX*(DDT/(PX+DDT)))/DT
4489
4490 ! ----------------------------------------------------------------------
4491 ! FROZEN GROUND VERSION:
4492 ! REDUCTION OF INFILTRATION BASED ON FROZEN GROUND PARAMETERS
4493 ! ----------------------------------------------------------------------
4494 FCR = 1.
4495 IF (DICE .GT. 1.E-2) THEN
4496 ACRT = CVFRZ * FRZX / DICE
4497 SUM = 1.
4498 IALP1 = CVFRZ - 1
4499 DO J = 1,IALP1
4500 K = 1
4501 DO JJ = J+1,IALP1
4502 K = K * JJ
4503 END DO
4504 SUM = SUM + (ACRT ** ( CVFRZ-J)) / FLOAT (K)
4505 END DO
4506 FCR = 1. - EXP(-ACRT) * SUM
4507 ENDIF
4508 INFMAX = INFMAX * FCR
4509
4510 ! ----------------------------------------------------------------------
4511 ! CORRECTION OF INFILTRATION LIMITATION:
4512 ! IF INFMAX .LE. HYDROLIC CONDUCTIVITY ASSIGN INFMAX THE VALUE OF
4513 ! HYDROLIC CONDUCTIVITY
4514 ! ----------------------------------------------------------------------
4515 ! MXSMC = MAX ( SH2OA(1), SH2OA(2) )
4516 MXSMC = SH2OA(1)
4517
4518 CALL WDFCND (WDF,WCND,MXSMC,SMCMAX,BEXP,DKSAT,DWSAT, &
4519 & SICEMAX)
4520
4521 INFMAX = MAX(INFMAX,WCND)
4522 INFMAX = MIN(INFMAX,PX)
4523
4524 IF (PCPDRP .GT. INFMAX) THEN
4525 RUNOFF1 = PCPDRP - INFMAX
4526 PDDUM = INFMAX
4527 ENDIF
4528
4529 ENDIF
4530
4531 ! ----------------------------------------------------------------------
4532 ! TO AVOID SPURIOUS DRAINAGE BEHAVIOR, 'UPSTREAM DIFFERENCING' IN LINE
4533 ! BELOW REPLACED WITH NEW APPROACH IN 2ND LINE:
4534 ! 'MXSMC = MAX(SH2OA(1), SH2OA(2))'
4535 ! ----------------------------------------------------------------------
4536 MXSMC = SH2OA(1)
4537
4538 CALL WDFCND (WDF,WCND,MXSMC,SMCMAX,BEXP,DKSAT,DWSAT, &
4539 & SICEMAX)
4540
4541 ! ----------------------------------------------------------------------
4542 ! CALC THE MATRIX COEFFICIENTS AI, BI, AND CI FOR THE TOP LAYER
4543 ! ----------------------------------------------------------------------
4544 DDZ = 1. / ( -.5 * ZSOIL(2) )
4545 AI(1) = 0.0
4546 BI(1) = WDF * DDZ / ( -ZSOIL(1) )
4547 CI(1) = -BI(1)
4548
4549 ! ----------------------------------------------------------------------
4550 ! CALC RHSTT FOR THE TOP LAYER AFTER CALC'NG THE VERTICAL SOIL MOISTURE
4551 ! GRADIENT BTWN THE TOP AND NEXT TO TOP LAYERS.
4552 ! ----------------------------------------------------------------------
4553 DSMDZ = ( SH2O(1) - SH2O(2) ) / ( -.5 * ZSOIL(2) )
4554 RHSTT(1) = (WDF * DSMDZ + WCND - PDDUM + EDIR + ET(1))/ZSOIL(1)
4555 SSTT = WDF * DSMDZ + WCND + EDIR + ET(1)
4556
4557 ! ----------------------------------------------------------------------
4558 ! INITIALIZE DDZ2
4559 ! ----------------------------------------------------------------------
4560 DDZ2 = 0.0
4561
4562 ! ----------------------------------------------------------------------
4563 ! LOOP THRU THE REMAINING SOIL LAYERS, REPEATING THE ABV PROCESS
4564 ! ----------------------------------------------------------------------
4565 DO K = 2,NSOIL
4566 DENOM2 = (ZSOIL(K-1) - ZSOIL(K))
4567 IF (K .NE. NSOIL) THEN
4568 SLOPX = 1.
4569
4570 ! ----------------------------------------------------------------------
4571 ! AGAIN, TO AVOID SPURIOUS DRAINAGE BEHAVIOR, 'UPSTREAM DIFFERENCING' IN
4572 ! LINE BELOW REPLACED WITH NEW APPROACH IN 2ND LINE:
4573 ! 'MXSMC2 = MAX (SH2OA(K), SH2OA(K+1))'
4574 ! ----------------------------------------------------------------------
4575 MXSMC2 = SH2OA(K)
4576
4577 CALL WDFCND (WDF2,WCND2,MXSMC2,SMCMAX,BEXP,DKSAT,DWSAT, &
4578 & SICEMAX)
4579
4580 ! ----------------------------------------------------------------------
4581 ! CALC SOME PARTIAL PRODUCTS FOR LATER USE IN CALC'NG RHSTT
4582 ! ----------------------------------------------------------------------
4583 DENOM = (ZSOIL(K-1) - ZSOIL(K+1))
4584 DSMDZ2 = (SH2O(K) - SH2O(K+1)) / (DENOM * 0.5)
4585
4586 ! ----------------------------------------------------------------------
4587 ! CALC THE MATRIX COEF, CI, AFTER CALC'NG ITS PARTIAL PRODUCT
4588 ! ----------------------------------------------------------------------
4589 DDZ2 = 2.0 / DENOM
4590 CI(K) = -WDF2 * DDZ2 / DENOM2
4591 ELSE
4592
4593 ! ----------------------------------------------------------------------
4594 ! SLOPE OF BOTTOM LAYER IS INTRODUCED
4595 ! ----------------------------------------------------------------------
4596 SLOPX = SLOPE
4597
4598 ! ----------------------------------------------------------------------
4599 ! RETRIEVE THE SOIL WATER DIFFUSIVITY AND HYDRAULIC CONDUCTIVITY FOR
4600 ! THIS LAYER
4601 ! ----------------------------------------------------------------------
4602 CALL WDFCND (WDF2,WCND2,SH2OA(NSOIL),SMCMAX,BEXP,DKSAT,DWSAT, &
4603 & SICEMAX)
4604
4605 ! ----------------------------------------------------------------------
4606 ! CALC A PARTIAL PRODUCT FOR LATER USE IN CALC'NG RHSTT
4607 ! ----------------------------------------------------------------------
4608 DSMDZ2 = 0.0
4609
4610 ! ----------------------------------------------------------------------
4611 ! SET MATRIX COEF CI TO ZERO
4612 ! ----------------------------------------------------------------------
4613 CI(K) = 0.0
4614 ENDIF
4615
4616 ! ----------------------------------------------------------------------
4617 ! CALC RHSTT FOR THIS LAYER AFTER CALC'NG ITS NUMERATOR
4618 ! ----------------------------------------------------------------------
4619 NUMER = (WDF2 * DSMDZ2) + SLOPX * WCND2 - (WDF * DSMDZ) &
4620 & - WCND + ET(K)
4621 RHSTT(K) = NUMER / (-DENOM2)
4622
4623 ! ----------------------------------------------------------------------
4624 ! CALC MATRIX COEFS, AI, AND BI FOR THIS LAYER
4625 ! ----------------------------------------------------------------------
4626 AI(K) = -WDF * DDZ / DENOM2
4627 BI(K) = -( AI(K) + CI(K) )
4628
4629 ! ----------------------------------------------------------------------
4630 ! RESET VALUES OF WDF, WCND, DSMDZ, AND DDZ FOR LOOP TO NEXT LYR
4631 ! RUNOFF2: SUB-SURFACE OR BASEFLOW RUNOFF
4632 ! ----------------------------------------------------------------------
4633 IF (K .EQ. NSOIL) THEN
4634 RUNOFF2 = SLOPX * WCND2
4635 ENDIF
4636
4637 IF (K .NE. NSOIL) THEN
4638 WDF = WDF2
4639 WCND = WCND2
4640 DSMDZ = DSMDZ2
4641 DDZ = DDZ2
4642 ENDIF
4643 END DO
4644
4645 ! ----------------------------------------------------------------------
4646 ! END SUBROUTINE SRT
4647 ! ----------------------------------------------------------------------
4648 END SUBROUTINE SRT
4649
4650 SUBROUTINE SSTEP (SH2OOUT,SH2OIN,CMC,RHSTT,RHSCT,DT, &
4651 & NSOIL,SMCMAX,CMCMAX,RUNOFF3,ZSOIL,SMC,SICE, &
4652 & AI,BI,CI)
4653
4654 IMPLICIT NONE
4655
4656 ! ----------------------------------------------------------------------
4657 ! SUBROUTINE SSTEP
4658 ! ----------------------------------------------------------------------
4659 ! CALCULATE/UPDATE SOIL MOISTURE CONTENT VALUES AND CANOPY MOISTURE
4660 ! CONTENT VALUES.
4661 ! ----------------------------------------------------------------------
4662 INTEGER NSOLD
4663 PARAMETER(NSOLD = 20)
4664
4665 INTEGER I
4666 INTEGER K
4667 INTEGER KK11
4668 INTEGER NSOIL
4669
4670 REAL AI(NSOLD)
4671 REAL BI(NSOLD)
4672 REAL CI(NSOLD)
4673 REAL CIin(NSOLD)
4674 REAL CMC
4675 REAL CMCMAX
4676 REAL DDZ
4677 REAL DT
4678 REAL RHSCT
4679 REAL RHSTT(NSOIL)
4680 REAL RHSTTin(NSOIL)
4681 REAL RUNOFF3
4682 REAL SH2OIN(NSOIL)
4683 REAL SH2OOUT(NSOIL)
4684 REAL SICE(NSOIL)
4685 REAL SMC(NSOIL)
4686 REAL SMCMAX
4687 REAL STOT
4688 REAL WPLUS
4689 REAL ZSOIL(NSOIL)
4690
4691 ! ----------------------------------------------------------------------
4692 ! CREATE 'AMOUNT' VALUES OF VARIABLES TO BE INPUT TO THE
4693 ! TRI-DIAGONAL MATRIX ROUTINE.
4694 ! ----------------------------------------------------------------------
4695 DO K = 1,NSOIL
4696 RHSTT(K) = RHSTT(K) * DT
4697 AI(K) = AI(K) * DT
4698 BI(K) = 1. + BI(K) * DT
4699 CI(K) = CI(K) * DT
4700 END DO
4701
4702 ! ----------------------------------------------------------------------
4703 ! COPY VALUES FOR INPUT VARIABLES BEFORE CALL TO ROSR12
4704 ! ----------------------------------------------------------------------
4705 DO K = 1,NSOIL
4706 RHSTTin(K) = RHSTT(K)
4707 END DO
4708 DO K = 1,NSOIL
4709 CIin(K) = CI(K)
4710 END DO
4711
4712 ! ----------------------------------------------------------------------
4713 ! CALL ROSR12 TO SOLVE THE TRI-DIAGONAL MATRIX
4714 ! ----------------------------------------------------------------------
4715 CALL ROSR12 (CI,AI,BI,CIin,RHSTTin,RHSTT,NSOIL)
4716
4717 ! ----------------------------------------------------------------------
4718 ! SUM THE PREVIOUS SMC VALUE AND THE MATRIX SOLUTION TO GET A
4719 ! NEW VALUE. MIN ALLOWABLE VALUE OF SMC WILL BE 0.02.
4720 ! RUNOFF3: RUNOFF WITHIN SOIL LAYERS
4721 ! ----------------------------------------------------------------------
4722 WPLUS = 0.0
4723 RUNOFF3 = 0.
4724 DDZ = -ZSOIL(1)
4725
4726 DO K = 1,NSOIL
4727 IF (K .NE. 1) DDZ = ZSOIL(K - 1) - ZSOIL(K)
4728 SH2OOUT(K) = SH2OIN(K) + CI(K) + WPLUS / DDZ
4729
4730 STOT = SH2OOUT(K) + SICE(K)
4731 IF (STOT .GT. SMCMAX) THEN
4732 IF (K .EQ. 1) THEN
4733 DDZ = -ZSOIL(1)
4734 ELSE
4735 KK11 = K - 1
4736 DDZ = -ZSOIL(K) + ZSOIL(KK11)
4737 ENDIF
4738 WPLUS = (STOT-SMCMAX) * DDZ
4739 ELSE
4740 WPLUS = 0.
4741 ENDIF
4742 SMC(K) = MAX ( MIN(STOT,SMCMAX),0.02 )
4743 SH2OOUT(K) = MAX((SMC(K)-SICE(K)),0.0)
4744 END DO
4745
4746 RUNOFF3 = WPLUS
4747
4748 ! ----------------------------------------------------------------------
4749 ! UPDATE CANOPY WATER CONTENT/INTERCEPTION (CMC). CONVERT RHSCT TO
4750 ! AN 'AMOUNT' VALUE AND ADD TO PREVIOUS CMC VALUE TO GET NEW CMC.
4751 ! ----------------------------------------------------------------------
4752 CMC = CMC + DT * RHSCT
4753 IF (CMC .LT. 1.E-20) CMC=0.0
4754 CMC = MIN(CMC,CMCMAX)
4755
4756 ! ----------------------------------------------------------------------
4757 ! END SUBROUTINE SSTEP
4758 ! ----------------------------------------------------------------------
4759 END SUBROUTINE SSTEP
4760
4761 SUBROUTINE TBND (TU,TB,ZSOIL,ZBOT,K,NSOIL,TBND1)
4762
4763 IMPLICIT NONE
4764
4765 ! ----------------------------------------------------------------------
4766 ! SUBROUTINE TBND
4767 ! ----------------------------------------------------------------------
4768 ! CALCULATE TEMPERATURE ON THE BOUNDARY OF THE LAYER BY INTERPOLATION OF
4769 ! THE MIDDLE LAYER TEMPERATURES
4770 ! ----------------------------------------------------------------------
4771 INTEGER NSOIL
4772 INTEGER K
4773
4774 REAL TBND1
4775 REAL T0
4776 REAL TU
4777 REAL TB
4778 REAL ZB
4779 REAL ZBOT
4780 REAL ZUP
4781 REAL ZSOIL (NSOIL)
4782
4783 PARAMETER(T0 = 273.15)
4784
4785 ! ----------------------------------------------------------------------
4786 ! USE SURFACE TEMPERATURE ON THE TOP OF THE FIRST LAYER
4787 ! ----------------------------------------------------------------------
4788 IF (K .EQ. 1) THEN
4789 ZUP = 0.
4790 ELSE
4791 ZUP = ZSOIL(K-1)
4792 ENDIF
4793
4794 ! ----------------------------------------------------------------------
4795 ! USE DEPTH OF THE CONSTANT BOTTOM TEMPERATURE WHEN INTERPOLATE
4796 ! TEMPERATURE INTO THE LAST LAYER BOUNDARY
4797 ! ----------------------------------------------------------------------
4798 IF (K .EQ. NSOIL) THEN
4799 ZB = 2.*ZBOT-ZSOIL(K)
4800 ELSE
4801 ZB = ZSOIL(K+1)
4802 ENDIF
4803
4804 ! ----------------------------------------------------------------------
4805 ! LINEAR INTERPOLATION BETWEEN THE AVERAGE LAYER TEMPERATURES
4806 ! ----------------------------------------------------------------------
4807 TBND1 = TU+(TB-TU)*(ZUP-ZSOIL(K))/(ZUP-ZB)
4808
4809 ! ----------------------------------------------------------------------
4810 ! END SUBROUTINE TBND
4811 ! ----------------------------------------------------------------------
4812 END SUBROUTINE TBND
4813
4814 SUBROUTINE TDFCND ( DF, SMC, QZ, SMCMAX, SH2O)
4815
4816 IMPLICIT NONE
4817
4818 ! ----------------------------------------------------------------------
4819 ! SUBROUTINE TDFCND
4820 ! ----------------------------------------------------------------------
4821 ! CALCULATE THERMAL DIFFUSIVITY AND CONDUCTIVITY OF THE SOIL FOR A GIVEN
4822 ! POINT AND TIME.
4823 ! ----------------------------------------------------------------------
4824 ! PETERS-LIDARD APPROACH (PETERS-LIDARD et al., 1998)
4825 ! June 2001 CHANGES: FROZEN SOIL CONDITION.
4826 ! ----------------------------------------------------------------------
4827 REAL DF
4828 REAL GAMMD
4829 REAL THKDRY
4830 REAL AKE
4831 REAL THKICE
4832 REAL THKO
4833 REAL THKQTZ
4834 REAL THKSAT
4835 REAL THKS
4836 REAL THKW
4837 REAL QZ
4838 REAL SATRATIO
4839 REAL SH2O
4840 REAL SMC
4841 REAL SMCMAX
4842 REAL XU
4843 REAL XUNFROZ
4844
4845 ! ----------------------------------------------------------------------
4846 ! WE NOW GET QUARTZ AS AN INPUT ARGUMENT (SET IN ROUTINE REDPRM):
4847 ! DATA QUARTZ /0.82, 0.10, 0.25, 0.60, 0.52,
4848 ! & 0.35, 0.60, 0.40, 0.82/
4849 ! ----------------------------------------------------------------------
4850 ! IF THE SOIL HAS ANY MOISTURE CONTENT COMPUTE A PARTIAL SUM/PRODUCT
4851 ! OTHERWISE USE A CONSTANT VALUE WHICH WORKS WELL WITH MOST SOILS
4852 ! ----------------------------------------------------------------------
4853 ! THKW ......WATER THERMAL CONDUCTIVITY
4854 ! THKQTZ ....THERMAL CONDUCTIVITY FOR QUARTZ
4855 ! THKO ......THERMAL CONDUCTIVITY FOR OTHER SOIL COMPONENTS
4856 ! THKS ......THERMAL CONDUCTIVITY FOR THE SOLIDS COMBINED(QUARTZ+OTHER)
4857 ! THKICE ....ICE THERMAL CONDUCTIVITY
4858 ! SMCMAX ....POROSITY (= SMCMAX)
4859 ! QZ .........QUARTZ CONTENT (SOIL TYPE DEPENDENT)
4860 ! ----------------------------------------------------------------------
4861 ! USE AS IN PETERS-LIDARD, 1998 (MODIF. FROM JOHANSEN, 1975).
4862 !
4863 ! PABLO GRUNMANN, 08/17/98
4864 ! REFS.:
4865 ! FAROUKI, O.T.,1986: THERMAL PROPERTIES OF SOILS. SERIES ON ROCK
4866 ! AND SOIL MECHANICS, VOL. 11, TRANS TECH, 136 PP.
4867 ! JOHANSEN, O., 1975: THERMAL CONDUCTIVITY OF SOILS. PH.D. THESIS,
4868 ! UNIVERSITY OF TRONDHEIM,
4869 ! PETERS-LIDARD, C. D., ET AL., 1998: THE EFFECT OF SOIL THERMAL
4870 ! CONDUCTIVITY PARAMETERIZATION ON SURFACE ENERGY FLUXES
4871 ! AND TEMPERATURES. JOURNAL OF THE ATMOSPHERIC SCIENCES,
4872 ! VOL. 55, PP. 1209-1224.
4873 ! ----------------------------------------------------------------------
4874 ! NEEDS PARAMETERS
4875 ! POROSITY(SOIL TYPE):
4876 ! POROS = SMCMAX
4877 ! SATURATION RATIO:
4878 SATRATIO = SMC/SMCMAX
4879
4880 ! PARAMETERS W/(M.K)
4881 THKICE = 2.2
4882 THKW = 0.57
4883 THKO = 2.0
4884 ! IF (QZ .LE. 0.2) THKO = 3.0
4885 THKQTZ = 7.7
4886 ! SOLIDS' CONDUCTIVITY
4887 THKS = (THKQTZ**QZ)*(THKO**(1.- QZ))
4888
4889 ! UNFROZEN FRACTION (FROM 1., i.e., 100%LIQUID, TO 0. (100% FROZEN))
4890 XUNFROZ = (SH2O + 1.E-9) / (SMC + 1.E-9)
4891
4892 ! UNFROZEN VOLUME FOR SATURATION (POROSITY*XUNFROZ)
4893 XU=XUNFROZ*SMCMAX
4894 ! SATURATED THERMAL CONDUCTIVITY
4895 THKSAT = THKS**(1.-SMCMAX)*THKICE**(SMCMAX-XU)*THKW**(XU)
4896
4897 ! DRY DENSITY IN KG/M3
4898 GAMMD = (1. - SMCMAX)*2700.
4899
4900 ! DRY THERMAL CONDUCTIVITY IN W.M-1.K-1
4901 THKDRY = (0.135*GAMMD + 64.7)/(2700. - 0.947*GAMMD)
4902
4903 IF ( (SH2O + 0.0005) .LT. SMC ) THEN
4904 ! FROZEN
4905 AKE = SATRATIO
4906 ELSE
4907 ! UNFROZEN
4908 ! RANGE OF VALIDITY FOR THE KERSTEN NUMBER (AKE)
4909 IF ( SATRATIO .GT. 0.1 ) THEN
4910
4911 ! KERSTEN NUMBER (USING "FINE" FORMULA, VALID FOR SOILS CONTAINING AT
4912 ! LEAST 5% OF PARTICLES WITH DIAMETER LESS THAN 2.E-6 METERS.)
4913 ! (FOR "COARSE" FORMULA, SEE PETERS-LIDARD ET AL., 1998).
4914
4915 AKE = LOG10(SATRATIO) + 1.0
4916
4917 ELSE
4918
4919 ! USE K = KDRY
4920 AKE = 0.0
4921
4922 ENDIF
4923 ENDIF
4924
4925 ! THERMAL CONDUCTIVITY
4926
4927 DF = AKE*(THKSAT - THKDRY) + THKDRY
4928
4929 ! ----------------------------------------------------------------------
4930 ! END SUBROUTINE TDFCND
4931 ! ----------------------------------------------------------------------
4932 END SUBROUTINE TDFCND
4933
4934 SUBROUTINE TMPAVG (TAVG,TUP,TM,TDN,ZSOIL,NSOIL,K)
4935
4936 IMPLICIT NONE
4937
4938 ! ----------------------------------------------------------------------
4939 ! SUBROUTINE TMPAVG
4940 ! ----------------------------------------------------------------------
4941 ! CALCULATE SOIL LAYER AVERAGE TEMPERATURE (TAVG) IN FREEZING/THAWING
4942 ! LAYER USING UP, DOWN, AND MIDDLE LAYER TEMPERATURES (TUP, TDN, TM),
4943 ! WHERE TUP IS AT TOP BOUNDARY OF LAYER, TDN IS AT BOTTOM BOUNDARY OF
4944 ! LAYER. TM IS LAYER PROGNOSTIC STATE TEMPERATURE.
4945 ! ----------------------------------------------------------------------
4946 INTEGER K
4947 INTEGER NSOIL
4948
4949 REAL DZ
4950 REAL DZH
4951 REAL T0
4952 REAL TAVG
4953 REAL TDN
4954 REAL TM
4955 REAL TUP
4956 REAL X0
4957 REAL XDN
4958 REAL XUP
4959 REAL ZSOIL (NSOIL)
4960
4961 PARAMETER(T0 = 2.7315E2)
4962
4963 ! ----------------------------------------------------------------------
4964 IF (K .EQ. 1) THEN
4965 DZ = -ZSOIL(1)
4966 ELSE
4967 DZ = ZSOIL(K-1)-ZSOIL(K)
4968 ENDIF
4969
4970 DZH=DZ*0.5
4971
4972 IF (TUP .LT. T0) THEN
4973 IF (TM .LT. T0) THEN
4974 IF (TDN .LT. T0) THEN
4975 ! ----------------------------------------------------------------------
4976 ! TUP, TM, TDN < T0
4977 ! ----------------------------------------------------------------------
4978 TAVG = (TUP + 2.0*TM + TDN)/ 4.0
4979 ELSE
4980 ! ----------------------------------------------------------------------
4981 ! TUP & TM < T0, TDN >= T0
4982 ! ----------------------------------------------------------------------
4983 X0 = (T0 - TM) * DZH / (TDN - TM)
4984 TAVG = 0.5 * (TUP*DZH+TM*(DZH+X0)+T0*(2.*DZH-X0)) / DZ
4985 ENDIF
4986 ELSE
4987 IF (TDN .LT. T0) THEN
4988 ! ----------------------------------------------------------------------
4989 ! TUP < T0, TM >= T0, TDN < T0
4990 ! ----------------------------------------------------------------------
4991 XUP = (T0-TUP) * DZH / (TM-TUP)
4992 XDN = DZH - (T0-TM) * DZH / (TDN-TM)
4993 TAVG = 0.5 * (TUP*XUP+T0*(2.*DZ-XUP-XDN)+TDN*XDN) / DZ
4994 ELSE
4995 ! ----------------------------------------------------------------------
4996 ! TUP < T0, TM >= T0, TDN >= T0
4997 ! ----------------------------------------------------------------------
4998 XUP = (T0-TUP) * DZH / (TM-TUP)
4999 TAVG = 0.5 * (TUP*XUP+T0*(2.*DZ-XUP)) / DZ
5000 ENDIF
5001 ENDIF
5002 ELSE
5003 IF (TM .LT. T0) THEN
5004 IF (TDN .LT. T0) THEN
5005 ! ----------------------------------------------------------------------
5006 ! TUP >= T0, TM < T0, TDN < T0
5007 ! ----------------------------------------------------------------------
5008 XUP = DZH - (T0-TUP) * DZH / (TM-TUP)
5009 TAVG = 0.5 * (T0*(DZ-XUP)+TM*(DZH+XUP)+TDN*DZH) / DZ
5010 ELSE
5011 ! ----------------------------------------------------------------------
5012 ! TUP >= T0, TM < T0, TDN >= T0
5013 ! ----------------------------------------------------------------------
5014 XUP = DZH - (T0-TUP) * DZH / (TM-TUP)
5015 XDN = (T0-TM) * DZH / (TDN-TM)
5016 TAVG = 0.5 * (T0*(2.*DZ-XUP-XDN)+TM*(XUP+XDN)) / DZ
5017 ENDIF
5018 ELSE
5019 IF (TDN .LT. T0) THEN
5020 ! ----------------------------------------------------------------------
5021 ! TUP >= T0, TM >= T0, TDN < T0
5022 ! ----------------------------------------------------------------------
5023 XDN = DZH - (T0-TM) * DZH / (TDN-TM)
5024 TAVG = (T0*(DZ-XDN)+0.5*(T0+TDN)*XDN) / DZ
5025 ELSE
5026 ! ----------------------------------------------------------------------
5027 ! TUP >= T0, TM >= T0, TDN >= T0
5028 ! ----------------------------------------------------------------------
5029 TAVG = (TUP + 2.0*TM + TDN) / 4.0
5030 ENDIF
5031 ENDIF
5032 ENDIF
5033 ! ----------------------------------------------------------------------
5034 ! END SUBROUTINE TMPAVG
5035 ! ----------------------------------------------------------------------
5036 END SUBROUTINE TMPAVG
5037
5038 SUBROUTINE TRANSP (ET1,NSOIL,ETP1,SMC,CMC,ZSOIL,SHDFAC,SMCWLT, &
5039 & CMCMAX,PC,CFACTR,SMCREF,SFCTMP,Q2,NROOT,RTDIS)
5040
5041 IMPLICIT NONE
5042
5043 ! ----------------------------------------------------------------------
5044 ! SUBROUTINE TRANSP
5045 ! ----------------------------------------------------------------------
5046 ! CALCULATE TRANSPIRATION FOR THE VEG CLASS.
5047 ! ----------------------------------------------------------------------
5048 INTEGER I
5049 INTEGER K
5050 INTEGER NSOIL
5051 INTEGER NROOT
5052
5053 REAL CFACTR
5054 REAL CMC
5055 REAL CMCMAX
5056 REAL DENOM
5057 REAL ET1(NSOIL)
5058 REAL ETP1
5059 REAL ETP1A
5060 REAL GX (7)
5061 !.....REAL PART(NSOIL)
5062 REAL PC
5063 REAL Q2
5064 REAL RTDIS(NSOIL)
5065 REAL RTX
5066 REAL SFCTMP
5067 REAL SGX
5068 REAL SHDFAC
5069 REAL SMC(NSOIL)
5070 REAL SMCREF
5071 REAL SMCWLT
5072 REAL ZSOIL(NSOIL)
5073
5074 ! ----------------------------------------------------------------------
5075 ! INITIALIZE PLANT TRANSP TO ZERO FOR ALL SOIL LAYERS.
5076 ! ----------------------------------------------------------------------
5077 DO K = 1,NSOIL
5078 ET1(K) = 0.
5079 END DO
5080
5081 ! ----------------------------------------------------------------------
5082 ! CALCULATE AN 'ADJUSTED' POTENTIAL TRANSPIRATION
5083 ! IF STATEMENT BELOW TO AVOID TANGENT LINEAR PROBLEMS NEAR ZERO
5084 ! NOTE: GX AND OTHER TERMS BELOW REDISTRIBUTE TRANSPIRATION BY LAYER,
5085 ! ET(K), AS A FUNCTION OF SOIL MOISTURE AVAILABILITY, WHILE PRESERVING
5086 ! TOTAL ETP1A.
5087 ! ----------------------------------------------------------------------
5088 IF (CMC .NE. 0.0) THEN
5089 ETP1A = SHDFAC * PC * ETP1 * (1.0 - (CMC /CMCMAX) ** CFACTR)
5090 ELSE
5091 ETP1A = SHDFAC * PC * ETP1
5092 ENDIF
5093
5094 SGX = 0.0
5095 DO I = 1,NROOT
5096 GX(I) = ( SMC(I) - SMCWLT ) / ( SMCREF - SMCWLT )
5097 GX(I) = MAX ( MIN ( GX(I), 1. ), 0. )
5098 SGX = SGX + GX (I)
5099 END DO
5100 SGX = SGX / NROOT
5101
5102 DENOM = 0.
5103 DO I = 1,NROOT
5104 RTX = RTDIS(I) + GX(I) - SGX
5105 GX(I) = GX(I) * MAX ( RTX, 0. )
5106 DENOM = DENOM + GX(I)
5107 END DO
5108 IF (DENOM .LE. 0.0) DENOM = 1.
5109
5110 DO I = 1,NROOT
5111 ET1(I) = ETP1A * GX(I) / DENOM
5112 END DO
5113
5114 ! ----------------------------------------------------------------------
5115 ! ABOVE CODE ASSUMES A VERTICALLY UNIFORM ROOT DISTRIBUTION
5116 ! CODE BELOW TESTS A VARIABLE ROOT DISTRIBUTION
5117 ! ----------------------------------------------------------------------
5118 ! ET(1) = ( ZSOIL(1) / ZSOIL(NROOT) ) * GX * ETP1A
5119 ! ET(1) = ( ZSOIL(1) / ZSOIL(NROOT) ) * ETP1A
5120 ! ----------------------------------------------------------------------
5121 ! USING ROOT DISTRIBUTION AS WEIGHTING FACTOR
5122 ! ----------------------------------------------------------------------
5123 ! ET(1) = RTDIS(1) * ETP1A
5124 ! ET(1) = ETP1A * PART(1)
5125 ! ----------------------------------------------------------------------
5126 ! LOOP DOWN THRU THE SOIL LAYERS REPEATING THE OPERATION ABOVE,
5127 ! BUT USING THE THICKNESS OF THE SOIL LAYER (RATHER THAN THE
5128 ! ABSOLUTE DEPTH OF EACH LAYER) IN THE FINAL CALCULATION.
5129 ! ----------------------------------------------------------------------
5130 ! DO K = 2,NROOT
5131 ! GX = ( SMC(K) - SMCWLT ) / ( SMCREF - SMCWLT )
5132 ! GX = MAX ( MIN ( GX, 1. ), 0. )
5133 ! TEST CANOPY RESISTANCE
5134 ! GX = 1.0
5135 ! ET(K) = ((ZSOIL(K)-ZSOIL(K-1))/ZSOIL(NROOT))*GX*ETP1A
5136 ! ET(K) = ((ZSOIL(K)-ZSOIL(K-1))/ZSOIL(NROOT))*ETP1A
5137 ! ----------------------------------------------------------------------
5138 ! USING ROOT DISTRIBUTION AS WEIGHTING FACTOR
5139 ! ----------------------------------------------------------------------
5140 ! ET(K) = RTDIS(K) * ETP1A
5141 ! ET(K) = ETP1A*PART(K)
5142 ! END DO
5143 ! ----------------------------------------------------------------------
5144 ! END SUBROUTINE TRANSP
5145 ! ----------------------------------------------------------------------
5146 END SUBROUTINE TRANSP
5147
5148 SUBROUTINE WDFCND (WDF,WCND,SMC,SMCMAX,BEXP,DKSAT,DWSAT, &
5149 & SICEMAX)
5150
5151 IMPLICIT NONE
5152
5153 ! ----------------------------------------------------------------------
5154 ! SUBROUTINE WDFCND
5155 ! ----------------------------------------------------------------------
5156 ! CALCULATE SOIL WATER DIFFUSIVITY AND SOIL HYDRAULIC CONDUCTIVITY.
5157 ! ----------------------------------------------------------------------
5158 REAL BEXP
5159 REAL DKSAT
5160 REAL DWSAT
5161 REAL EXPON
5162 REAL FACTR1
5163 REAL FACTR2
5164 REAL SICEMAX
5165 REAL SMC
5166 REAL SMCMAX
5167 REAL VKwgt
5168 REAL WCND
5169 REAL WDF
5170
5171 ! ----------------------------------------------------------------------
5172 ! CALC THE RATIO OF THE ACTUAL TO THE MAX PSBL SOIL H2O CONTENT
5173 ! ----------------------------------------------------------------------
5174 SMC = SMC
5175 SMCMAX = SMCMAX
5176 FACTR1 = 0.2 / SMCMAX
5177 FACTR2 = SMC / SMCMAX
5178
5179 ! ----------------------------------------------------------------------
5180 ! PREP AN EXPNTL COEF AND CALC THE SOIL WATER DIFFUSIVITY
5181 ! ----------------------------------------------------------------------
5182 EXPON = BEXP + 2.0
5183 WDF = DWSAT * FACTR2 ** EXPON
5184
5185 ! ----------------------------------------------------------------------
5186 ! FROZEN SOIL HYDRAULIC DIFFUSIVITY. VERY SENSITIVE TO THE VERTICAL
5187 ! GRADIENT OF UNFROZEN WATER. THE LATTER GRADIENT CAN BECOME VERY
5188 ! EXTREME IN FREEZING/THAWING SITUATIONS, AND GIVEN THE RELATIVELY
5189 ! FEW AND THICK SOIL LAYERS, THIS GRADIENT SUFFERES SERIOUS
5190 ! TRUNCTION ERRORS YIELDING ERRONEOUSLY HIGH VERTICAL TRANSPORTS OF
5191 ! UNFROZEN WATER IN BOTH DIRECTIONS FROM HUGE HYDRAULIC DIFFUSIVITY.
5192 ! THEREFORE, WE FOUND WE HAD TO ARBITRARILY CONSTRAIN WDF
5193 ! --
5194 ! VERSION D_10CM: ........ FACTR1 = 0.2/SMCMAX
5195 ! WEIGHTED APPROACH...................... PABLO GRUNMANN, 28_SEP_1999.
5196 ! ----------------------------------------------------------------------
5197 IF (SICEMAX .GT. 0.0) THEN
5198 VKWGT = 1./(1.+(500.*SICEMAX)**3.)
5199 WDF = VKWGT*WDF + (1.- VKWGT)*DWSAT*FACTR1**EXPON
5200 ENDIF
5201
5202 ! ----------------------------------------------------------------------
5203 ! RESET THE EXPNTL COEF AND CALC THE HYDRAULIC CONDUCTIVITY
5204 ! ----------------------------------------------------------------------
5205 EXPON = (2.0 * BEXP) + 3.0
5206 WCND = DKSAT * FACTR2 ** EXPON
5207
5208 ! ----------------------------------------------------------------------
5209 ! END SUBROUTINE WDFCND
5210 ! ----------------------------------------------------------------------
5211 END SUBROUTINE WDFCND
5212
5213 SUBROUTINE nmmlsminit(isn,XICE,VEGFRA,SNOW,SNOWC,CANWAT,SMSTAV, &
5214 SMSTOT, SFCRUNOFF,UDRUNOFF,GRDFLX,ACSNOW, &
5215 ACSNOM,IVGTYP,ISLTYP,TSLB,SMOIS,DZS,SFCEVP, & ! STEMP
5216 TMN, &
5217 num_soil_layers, &
5218 allowed_to_read, &
5219 ids,ide, jds,jde, kds,kde, &
5220 ims,ime, jms,jme, kms,kme, &
5221 its,ite, jts,jte, kts,kte )
5222
5223 IMPLICIT NONE
5224
5225 ! Arguments
5226 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
5227 ims,ime, jms,jme, kms,kme, &
5228 its,ite, jts,jte, kts,kte
5229
5230 INTEGER, INTENT(IN) :: num_soil_layers
5231
5232 REAL, DIMENSION( num_soil_layers), INTENT(IN) :: DZS
5233
5234 REAL, DIMENSION( ims:ime, num_soil_layers, jms:jme ) , &
5235 INTENT(INOUT) :: SMOIS, &
5236 TSLB !STEMP
5237
5238 REAL, DIMENSION( ims:ime, jms:jme ) , &
5239 INTENT(INOUT) :: SNOW, &
5240 SNOWC, &
5241 CANWAT, &
5242 SMSTAV, &
5243 SMSTOT, &
5244 SFCRUNOFF, &
5245 UDRUNOFF, &
5246 SFCEVP, &
5247 GRDFLX, &
5248 ACSNOW, &
5249 XICE, &
5250 VEGFRA, &
5251 TMN, &
5252 ACSNOM
5253
5254 INTEGER, DIMENSION( ims:ime, jms:jme ) , &
5255 INTENT(INOUT) :: IVGTYP, &
5256 ISLTYP
5257
5258 !
5259
5260 INTEGER, INTENT(IN) :: isn
5261 LOGICAL, INTENT(IN) :: allowed_to_read
5262 ! Local
5263 INTEGER :: iseason
5264 INTEGER :: icm,jcm,itf,jtf
5265 INTEGER :: I,J,L
5266
5267
5268 itf=min0(ite,ide-1)
5269 jtf=min0(jte,jde-1)
5270
5271 icm = ide/2
5272 jcm = jde/2
5273
5274 iseason=isn
5275
5276 DO J=jts,jtf
5277 DO I=its,itf
5278 ! SNOW(i,j)=0.
5279 SNOWC(i,j)=0.
5280 ! SMSTAV(i,j)=
5281 ! SMSTOT(i,j)=
5282 ! SFCRUNOFF(i,j)=
5283 ! UDRUNOFF(i,j)=
5284 ! GRDFLX(i,j)=
5285 ! ACSNOW(i,j)=
5286 ! ACSNOM(i,j)=
5287 ENDDO
5288 ENDDO
5289
5290 END SUBROUTINE nmmlsminit
5291
5292 FUNCTION CSNOW (DSNOW)
5293
5294 IMPLICIT NONE
5295
5296 ! ----------------------------------------------------------------------
5297 ! FUNCTION CSNOW
5298 ! ----------------------------------------------------------------------
5299 ! CALCULATE SNOW TERMAL CONDUCTIVITY
5300 ! ----------------------------------------------------------------------
5301 REAL C
5302 REAL DSNOW
5303 REAL CSNOW
5304 REAL UNIT
5305
5306 PARAMETER(UNIT = 0.11631)
5307
5308 ! ----------------------------------------------------------------------
5309 ! CSNOW IN UNITS OF CAL/(CM*HR*C), RETURNED IN W/(M*C)
5310 ! BASIC VERSION IS DYACHKOVA EQUATION (1960), FOR RANGE 0.1-0.4
5311 ! ----------------------------------------------------------------------
5312 C=0.328*10**(2.25*DSNOW)
5313 ! CSNOW=UNIT*C
5314 ! MEK JAN 2006, DOUBLE SNOW THERMAL CONDUCTIVITY
5315 CSNOW=2.0*UNIT*C
5316
5317 ! ----------------------------------------------------------------------
5318 ! DE VAUX EQUATION (1933), IN RANGE 0.1-0.6
5319 ! ----------------------------------------------------------------------
5320 ! CSNOW=0.0293*(1.+100.*DSNOW**2)
5321
5322 ! ----------------------------------------------------------------------
5323 ! E. ANDERSEN FROM FLERCHINGER
5324 ! ----------------------------------------------------------------------
5325 ! CSNOW=0.021+2.51*DSNOW**2
5326
5327 ! ----------------------------------------------------------------------
5328 ! END FUNCTION CSNOW
5329 ! ----------------------------------------------------------------------
5330 END FUNCTION CSNOW
5331
5332 FUNCTION FRH2O (TKELV,SMC,SH2O,SMCMAX,BEXP,PSIS)
5333
5334 IMPLICIT NONE
5335
5336 ! ----------------------------------------------------------------------
5337 ! FUNCTION FRH2O
5338 ! ----------------------------------------------------------------------
5339 ! CALCULATE AMOUNT OF SUPERCOOLED LIQUID SOIL WATER CONTENT IF
5340 ! TEMPERATURE IS BELOW 273.15K (T0). REQUIRES NEWTON-TYPE ITERATION TO
5341 ! SOLVE THE NONLINEAR IMPLICIT EQUATION GIVEN IN EQN 17 OF KOREN ET AL
5342 ! (1999, JGR, VOL 104(D16), 19569-19585).
5343 ! ----------------------------------------------------------------------
5344 ! NEW VERSION (JUNE 2001): MUCH FASTER AND MORE ACCURATE NEWTON
5345 ! ITERATION ACHIEVED BY FIRST TAKING LOG OF EQN CITED ABOVE -- LESS THAN
5346 ! 4 (TYPICALLY 1 OR 2) ITERATIONS ACHIEVES CONVERGENCE. ALSO, EXPLICIT
5347 ! 1-STEP SOLUTION OPTION FOR SPECIAL CASE OF PARAMETER CK=0, WHICH
5348 ! REDUCES THE ORIGINAL IMPLICIT EQUATION TO A SIMPLER EXPLICIT FORM,
5349 ! KNOWN AS THE "FLERCHINGER EQN". IMPROVED HANDLING OF SOLUTION IN THE
5350 ! LIMIT OF FREEZING POINT TEMPERATURE T0.
5351 ! ----------------------------------------------------------------------
5352 ! INPUT:
5353 !
5354 ! TKELV.........TEMPERATURE (Kelvin)
5355 ! SMC...........TOTAL SOIL MOISTURE CONTENT (VOLUMETRIC)
5356 ! SH2O..........LIQUID SOIL MOISTURE CONTENT (VOLUMETRIC)
5357 ! SMCMAX........SATURATION SOIL MOISTURE CONTENT (FROM REDPRM)
5358 ! B.............SOIL TYPE "B" PARAMETER (FROM REDPRM)
5359 ! PSIS..........SATURATED SOIL MATRIC POTENTIAL (FROM REDPRM)
5360 !
5361 ! OUTPUT:
5362 ! FRH2O.........SUPERCOOLED LIQUID WATER CONTENT
5363 ! ----------------------------------------------------------------------
5364 REAL BEXP
5365 REAL BLIM
5366 REAL BX
5367 REAL CK
5368 REAL DENOM
5369 REAL DF
5370 REAL DH2O
5371 REAL DICE
5372 REAL DSWL
5373 REAL ERROR
5374 REAL FK
5375 REAL FRH2O
5376 REAL GS
5377 REAL HLICE
5378 REAL PSIS
5379 REAL SH2O
5380 REAL SMC
5381 REAL SMCMAX
5382 REAL SWL
5383 REAL SWLK
5384 REAL TKELV
5385 REAL T0
5386
5387 INTEGER NLOG
5388 INTEGER KCOUNT
5389
5390 PARAMETER(CK = 8.0)
5391 ! PARAMETER(CK = 0.0)
5392 PARAMETER(BLIM = 5.5)
5393 PARAMETER(ERROR = 0.005)
5394
5395 PARAMETER(HLICE = 3.335E5)
5396 PARAMETER(GS = 9.81)
5397 PARAMETER(DICE = 920.0)
5398 PARAMETER(DH2O = 1000.0)
5399 PARAMETER(T0 = 273.15)
5400
5401 ! ----------------------------------------------------------------------
5402 ! LIMITS ON PARAMETER B: B < 5.5 (use parameter BLIM)
5403 ! SIMULATIONS SHOWED IF B > 5.5 UNFROZEN WATER CONTENT IS
5404 ! NON-REALISTICALLY HIGH AT VERY LOW TEMPERATURES.
5405 ! ----------------------------------------------------------------------
5406 BX = BEXP
5407 IF (BEXP .GT. BLIM) BX = BLIM
5408
5409 ! ----------------------------------------------------------------------
5410 ! INITIALIZING ITERATIONS COUNTER AND ITERATIVE SOLUTION FLAG.
5411 ! ----------------------------------------------------------------------
5412 NLOG=0
5413 KCOUNT=0
5414
5415 ! ----------------------------------------------------------------------
5416 ! IF TEMPERATURE NOT SIGNIFICANTLY BELOW FREEZING (T0), SH2O = SMC
5417 ! ----------------------------------------------------------------------
5418 IF (TKELV .GT. (T0 - 1.E-3)) THEN
5419 FRH2O = SMC
5420 ELSE
5421 IF (CK .NE. 0.0) THEN
5422
5423 ! ----------------------------------------------------------------------
5424 ! OPTION 1: ITERATED SOLUTION FOR NONZERO CK
5425 ! IN KOREN ET AL, JGR, 1999, EQN 17
5426 ! ----------------------------------------------------------------------
5427 ! INITIAL GUESS FOR SWL (frozen content)
5428 ! ----------------------------------------------------------------------
5429 SWL = SMC-SH2O
5430
5431 ! ----------------------------------------------------------------------
5432 ! KEEP WITHIN BOUNDS.
5433 ! ----------------------------------------------------------------------
5434 IF (SWL .GT. (SMC-0.02)) SWL = SMC-0.02
5435 IF (SWL .LT. 0.) SWL = 0.
5436
5437 ! ----------------------------------------------------------------------
5438 ! START OF ITERATIONS
5439 ! ----------------------------------------------------------------------
5440 DO WHILE ( (NLOG .LT. 10) .AND. (KCOUNT .EQ. 0) )
5441 NLOG = NLOG+1
5442 DF = ALOG(( PSIS*GS/HLICE ) * ( ( 1.+CK*SWL )**2. ) * &
5443 & ( SMCMAX/(SMC-SWL) )**BX) - ALOG(-(TKELV-T0)/TKELV)
5444 DENOM = 2. * CK / ( 1.+CK*SWL ) + BX / ( SMC - SWL )
5445 SWLK = SWL - DF/DENOM
5446 ! ----------------------------------------------------------------------
5447 ! BOUNDS USEFUL FOR MATHEMATICAL SOLUTION.
5448 ! ----------------------------------------------------------------------
5449 IF (SWLK .GT. (SMC-0.02)) SWLK = SMC - 0.02
5450 IF (SWLK .LT. 0.) SWLK = 0.
5451
5452 ! ----------------------------------------------------------------------
5453 ! MATHEMATICAL SOLUTION BOUNDS APPLIED.
5454 ! ----------------------------------------------------------------------
5455 DSWL = ABS(SWLK-SWL)
5456 SWL = SWLK
5457
5458 ! ----------------------------------------------------------------------
5459 ! IF MORE THAN 10 ITERATIONS, USE EXPLICIT METHOD (CK=0 APPROX.)
5460 ! WHEN DSWL LESS OR EQ. ERROR, NO MORE ITERATIONS REQUIRED.
5461 ! ----------------------------------------------------------------------
5462 IF ( DSWL .LE. ERROR ) THEN
5463 KCOUNT = KCOUNT+1
5464 ENDIF
5465 END DO
5466
5467 ! ----------------------------------------------------------------------
5468 ! END OF ITERATIONS
5469 ! ----------------------------------------------------------------------
5470 ! BOUNDS APPLIED WITHIN DO-BLOCK ARE VALID FOR PHYSICAL SOLUTION.
5471 ! ----------------------------------------------------------------------
5472 FRH2O = SMC - SWL
5473
5474 ! ----------------------------------------------------------------------
5475 ! END OPTION 1
5476 ! ----------------------------------------------------------------------
5477 ENDIF
5478
5479 ! ----------------------------------------------------------------------
5480 ! OPTION 2: EXPLICIT SOLUTION FOR FLERCHINGER EQ. i.e. CK=0
5481 ! IN KOREN ET AL., JGR, 1999, EQN 17
5482 ! APPLY PHYSICAL BOUNDS TO FLERCHINGER SOLUTION
5483 ! ----------------------------------------------------------------------
5484 IF (KCOUNT .EQ. 0) THEN
5485 Print*,'Flerchinger used in NEW version. Iterations=',NLOG
5486 FK = (((HLICE/(GS*(-PSIS)))* &
5487 & ((TKELV-T0)/TKELV))**(-1/BX))*SMCMAX
5488 IF (FK .LT. 0.02) FK = 0.02
5489 FRH2O = MIN (FK, SMC)
5490 ! ----------------------------------------------------------------------
5491 ! END OPTION 2
5492 ! ----------------------------------------------------------------------
5493 ENDIF
5494
5495 ENDIF
5496
5497 ! ----------------------------------------------------------------------
5498 ! END FUNCTION FRH2O
5499 ! ----------------------------------------------------------------------
5500 END FUNCTION FRH2O
5501
5502 FUNCTION SNKSRC (TAVG,SMC,SH2O,ZSOIL,NSOIL, &
5503 & SMCMAX,PSISAT,BEXP,DT,K,QTOT)
5504
5505 IMPLICIT NONE
5506
5507 ! ----------------------------------------------------------------------
5508 ! FUNCTION SNKSRC
5509 ! ----------------------------------------------------------------------
5510 ! CALCULATE SINK/SOURCE TERM OF THE TERMAL DIFFUSION EQUATION. (SH2O) IS
5511 ! AVAILABLE LIQUED WATER.
5512 ! ----------------------------------------------------------------------
5513 INTEGER K
5514 INTEGER NSOIL
5515
5516 REAL BEXP
5517 REAL DF
5518 REAL DH2O
5519 REAL DT
5520 REAL DZ
5521 REAL DZH
5522 REAL FREE
5523 ! REAL FRH2O
5524 REAL HLICE
5525 REAL PSISAT
5526 REAL QTOT
5527 REAL SH2O
5528 REAL SMC
5529 REAL SMCMAX
5530 REAL SNKSRC
5531 REAL T0
5532 REAL TAVG
5533 REAL TDN
5534 REAL TM
5535 REAL TUP
5536 REAL TZ
5537 REAL X0
5538 REAL XDN
5539 REAL XH2O
5540 REAL XUP
5541 REAL ZSOIL (NSOIL)
5542
5543 PARAMETER(DH2O = 1.0000E3)
5544 PARAMETER(HLICE = 3.3350E5)
5545 PARAMETER(T0 = 2.7315E2)
5546
5547 IF (K .EQ. 1) THEN
5548 DZ = -ZSOIL(1)
5549 ELSE
5550 DZ = ZSOIL(K-1)-ZSOIL(K)
5551 ENDIF
5552
5553 ! ----------------------------------------------------------------------
5554 ! VIA FUNCTION FRH2O, COMPUTE POTENTIAL OR 'EQUILIBRIUM' UNFROZEN
5555 ! SUPERCOOLED FREE WATER FOR GIVEN SOIL TYPE AND SOIL LAYER TEMPERATURE.
5556 ! FUNCTION FRH20 INVOKES EQN (17) FROM V. KOREN ET AL (1999, JGR, VOL.
5557 ! 104, PG 19573). (ASIDE: LATTER EQN IN JOURNAL IN CENTIGRADE UNITS.
5558 ! ROUTINE FRH2O USE FORM OF EQN IN KELVIN UNITS.)
5559 ! ----------------------------------------------------------------------
5560 FREE = FRH2O(TAVG,SMC,SH2O,SMCMAX,BEXP,PSISAT)
5561
5562 ! ----------------------------------------------------------------------
5563 ! IN NEXT BLOCK OF CODE, INVOKE EQN 18 OF V. KOREN ET AL (1999, JGR,
5564 ! VOL. 104, PG 19573.) THAT IS, FIRST ESTIMATE THE NEW AMOUNTOF LIQUID
5565 ! WATER, 'XH2O', IMPLIED BY THE SUM OF (1) THE LIQUID WATER AT THE BEGIN
5566 ! OF CURRENT TIME STEP, AND (2) THE FREEZE OF THAW CHANGE IN LIQUID
5567 ! WATER IMPLIED BY THE HEAT FLUX 'QTOT' PASSED IN FROM ROUTINE HRT.
5568 ! SECOND, DETERMINE IF XH2O NEEDS TO BE BOUNDED BY 'FREE' (EQUIL AMT) OR
5569 ! IF 'FREE' NEEDS TO BE BOUNDED BY XH2O.
5570 ! ----------------------------------------------------------------------
5571 XH2O = SH2O + QTOT*DT/(DH2O*HLICE*DZ)
5572
5573 ! ----------------------------------------------------------------------
5574 ! FIRST, IF FREEZING AND REMAINING LIQUID LESS THAN LOWER BOUND, THEN
5575 ! REDUCE EXTENT OF FREEZING, THEREBY LETTING SOME OR ALL OF HEAT FLUX
5576 ! QTOT COOL THE SOIL TEMP LATER IN ROUTINE HRT.
5577 ! ----------------------------------------------------------------------
5578 IF ( XH2O .LT. SH2O .AND. XH2O .LT. FREE) THEN
5579 IF ( FREE .GT. SH2O ) THEN
5580 XH2O = SH2O
5581 ELSE
5582 XH2O = FREE
5583 ENDIF
5584 ENDIF
5585
5586 ! ----------------------------------------------------------------------
5587 ! SECOND, IF THAWING AND THE INCREASE IN LIQUID WATER GREATER THAN UPPER
5588 ! BOUND, THEN REDUCE EXTENT OF THAW, THEREBY LETTING SOME OR ALL OF HEAT
5589 ! FLUX QTOT WARM THE SOIL TEMP LATER IN ROUTINE HRT.
5590 ! ----------------------------------------------------------------------
5591 IF ( XH2O .GT. SH2O .AND. XH2O .GT. FREE ) THEN
5592 IF ( FREE .LT. SH2O ) THEN
5593 XH2O = SH2O
5594 ELSE
5595 XH2O = FREE
5596 ENDIF
5597 ENDIF
5598
5599 IF (XH2O .LT. 0.) XH2O = 0.
5600 IF (XH2O .GT. SMC) XH2O = SMC
5601
5602 ! ----------------------------------------------------------------------
5603 ! CALCULATE PHASE-CHANGE HEAT SOURCE/SINK TERM FOR USE IN ROUTINE HRT
5604 ! AND UPDATE LIQUID WATER TO REFLCET FINAL FREEZE/THAW INCREMENT.
5605 ! ----------------------------------------------------------------------
5606 SNKSRC = -DH2O*HLICE*DZ*(XH2O-SH2O)/DT
5607 SH2O = XH2O
5608
5609 ! ----------------------------------------------------------------------
5610 ! END FUNCTION SNKSRC
5611 ! ----------------------------------------------------------------------
5612 END FUNCTION SNKSRC
5613
5614 END MODULE module_sf_lsm_nmm
5615