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