module_sf_ruclsm.F
References to this file elsewhere.
1 #define LSMRUC_DBG_LVL 3000
2 !WRF:MODEL_LAYER:PHYSICS
3 !
4 MODULE module_sf_ruclsm
5
6 USE module_model_constants
7 USE module_wrf_error
8
9 ! VEGETATION PARAMETERS
10 INTEGER :: LUCATS , BARE
11 integer, PARAMETER :: NLUS=50
12 CHARACTER*8 LUTYPE
13 INTEGER, DIMENSION(1:NLUS) :: NROTBL
14 real, dimension(1:NLUS) :: SNUPTBL, RSTBL, RGLTBL, HSTBL, LAITBL, &
15 ALBTBL, Z0TBL, LEMITBL, PCTBL, SHDTBL, MAXALB
16 REAL :: TOPT_DATA,CMCMAX_DATA,CFACTR_DATA,RSMAX_DATA
17
18 ! SOIL PARAMETERS
19 INTEGER :: SLCATS
20 INTEGER, PARAMETER :: NSLTYPE=30
21 CHARACTER*8 SLTYPE
22 REAL, DIMENSION (1:NSLTYPE) :: BB,DRYSMC,HC, &
23 MAXSMC, REFSMC,SATPSI,SATDK,SATDW, WLTSMC,QTZ
24
25 ! LSM GENERAL PARAMETERS
26 INTEGER :: SLPCATS
27 INTEGER, PARAMETER :: NSLOPE=30
28 REAL, DIMENSION (1:NSLOPE) :: SLOPE_DATA
29 REAL :: SBETA_DATA,FXEXP_DATA,CSOIL_DATA,SALP_DATA,REFDK_DATA, &
30 REFKDT_DATA,FRZK_DATA,ZBOT_DATA, SMLOW_DATA,SMHIGH_DATA, &
31 CZIL_DATA
32
33 CHARACTER*256 :: err_message
34
35
36 CONTAINS
37
38 !-----------------------------------------------------------------
39 SUBROUTINE LSMRUC( &
40 DT,KTAU,NSL,ZS, &
41 RAINBL,SNOW,SNOWH,SNOWC,FRZFRAC,frpcpn, &
42 Z3D,P8W,T3D,QV3D,QC3D,RHO3D, & !p8W in [PA]
43 GLW,GSW,EMISS,CHKLOWQ, CHS, &
44 FLQC,FLHC,MAVAIL,CANWAT,VEGFRA,ALB,ZNT, &
45 SNOALB,ALBBCK, & !new
46 QSFC,QSG,QVG,QCG,SOILT1,TSNAV, &
47 TBOT,IVGTYP,ISLTYP,XLAND,XICE, &
48 CP,G0,LV,STBOLT, &
49 SOILMOIS,SH2O,SMAVAIL,SMMAX, &
50 TSO,SOILT,HFX,QFX,LH, &
51 SFCRUNOFF,UDRUNOFF,SFCEXC, &
52 SFCEVP,GRDFLX,ACSNOW, &
53 SMFR3D,KEEPFR3DFLAG, &
54 myj, &
55 ids,ide, jds,jde, kds,kde, &
56 ims,ime, jms,jme, kms,kme, &
57 its,ite, jts,jte, kts,kte )
58 !-----------------------------------------------------------------
59 IMPLICIT NONE
60 !-----------------------------------------------------------------
61 !
62 ! The RUC LSM model is described in:
63 ! Smirnova, T.G., J.M. Brown, and S.G. Benjamin, 1997:
64 ! Performance of different soil model configurations in simulating
65 ! ground surface temperature and surface fluxes.
66 ! Mon. Wea. Rev. 125, 1870-1884.
67 ! Smirnova, T.G., J.M. Brown, and D. Kim, 2000: Parameterization of
68 ! cold-season processes in the MAPS land-surface scheme.
69 ! J. Geophys. Res. 105, 4077-4086.
70 !-----------------------------------------------------------------
71 !-- DT time step (second)
72 ! ktau - number of time step
73 ! NSL - number of soil layers
74 ! NZS - number of levels in soil
75 ! ZS - depth of soil levels (m)
76 !-- RAINBL - accumulated rain in [mm] between the PBL calls
77 !-- RAINNCV one time step grid scale precipitation (mm/step)
78 ! SNOW - snow water equivalent [mm]
79 ! FRAZFRAC - fraction of frozen precipitation
80 !-- SNOWC flag indicating snow coverage (1 for snow cover)
81 !-- Z3D heights (m)
82 !-- P8W 3D pressure (Pa)
83 !-- T3D temperature (K)
84 !-- QV3D 3D water vapor mixing ratio (Kg/Kg)
85 ! QC3D - 3D cloud water mixing ratio (Kg/Kg)
86 ! RHO3D - 3D air density (kg/m^3)
87 !-- GLW downward long wave flux at ground surface (W/m^2)
88 !-- GSW absorbed short wave flux at ground surface (W/m^2)
89 !-- EMISS surface emissivity (between 0 and 1)
90 ! FLQC - surface exchange coefficient for moisture (kg/m^2/s)
91 ! FLHC - surface exchange coefficient for heat [W/m^2/s/degreeK]
92 ! SFCEXC - surface exchange coefficient for heat [m/s]
93 ! CANWAT - CANOPY MOISTURE CONTENT (mm)
94 ! VEGFRA - vegetation fraction (between 0 and 1)
95 ! ALB - surface albedo (between 0 and 1)
96 ! SNOALB - maximum snow albedo (between 0 and 1)
97 ! ALBBCK - snow-free albedo (between 0 and 1)
98 ! ZNT - roughness length [m]
99 !-- TBOT soil temperature at lower boundary (K)
100 ! IVGTYP - USGS vegetation type (24 classes)
101 ! ISLTYP - STASGO soil type (16 classes)
102 !-- XLAND land mask (1 for land, 2 for water)
103 !-- CP heat capacity at constant pressure for dry air (J/kg/K)
104 !-- G0 acceleration due to gravity (m/s^2)
105 !-- LV latent heat of melting (J/kg)
106 !-- STBOLT Stefan-Boltzmann constant (W/m^2/K^4)
107 ! SOILMOIS - soil moisture content (volumetric fraction)
108 ! TSO - soil temp (K)
109 !-- SOILT surface temperature (K)
110 !-- HFX upward heat flux at the surface (W/m^2)
111 !-- QFX upward moisture flux at the surface (kg/m^2/s)
112 !-- LH upward latent heat flux (W/m^2)
113 ! SFCRUNOFF - ground surface runoff [mm]
114 ! UDRUNOFF - underground runoff [mm]
115 ! SFCEVP - total evaporation in [kg/m^2]
116 ! GRDFLX - soil heat flux (W/m^2: negative, if downward from surface)
117 ! ACSNOW - accumulation of snow water [m]
118 !-- CHKLOWQ - is either 0 or 1 (so far set equal to 1).
119 !-- used only in MYJPBL.
120 !-- ims start index for i in memory
121 !-- ime end index for i in memory
122 !-- jms start index for j in memory
123 !-- jme end index for j in memory
124 !-- kms start index for k in memory
125 !-- kme end index for k in memory
126 !-------------------------------------------------------------------------
127 ! INTEGER, PARAMETER :: nzss=5
128 ! INTEGER, PARAMETER :: nddzs=2*(nzss-2)
129
130 INTEGER, PARAMETER :: nvegclas=24+3
131
132 REAL, INTENT(IN ) :: DT
133 LOGICAL, INTENT(IN ) :: myj,frpcpn
134 INTEGER, INTENT(IN ) :: ktau, nsl, &
135 ims,ime, jms,jme, kms,kme, &
136 ids,ide, jds,jde, kds,kde, &
137 its,ite, jts,jte, kts,kte
138
139 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
140 INTENT(IN ) :: QV3D, &
141 QC3D, &
142 p8w, &
143 rho3D, &
144 T3D, &
145 z3D
146
147 REAL, DIMENSION( ims:ime , jms:jme ), &
148 INTENT(IN ) :: RAINBL, &
149 GLW, &
150 GSW, &
151 SNOALB, &
152 ALBBCK, &
153 FLHC, &
154 FLQC, &
155 CHS , &
156 EMISS, &
157 XICE, &
158 XLAND, &
159 VEGFRA, &
160 TBOT
161
162 REAL, DIMENSION( 1:nsl), INTENT(IN ) :: ZS
163
164 REAL, DIMENSION( ims:ime , jms:jme ), &
165 INTENT(INOUT) :: &
166 SNOW, & !new
167 SNOWH, &
168 SNOWC, &
169 CANWAT, & ! new
170 ALB, &
171 MAVAIL, &
172 SFCEXC, &
173 ZNT
174
175 REAL, DIMENSION( ims:ime , jms:jme ), &
176 INTENT(IN ) :: &
177 FRZFRAC
178
179 INTEGER, DIMENSION( ims:ime , jms:jme ), &
180 INTENT(IN ) :: IVGTYP, &
181 ISLTYP
182
183 REAL, INTENT(IN ) :: CP,G0,LV,STBOLT
184
185 REAL, DIMENSION( ims:ime , 1:nsl, jms:jme ) , &
186 INTENT(INOUT) :: SOILMOIS,SH2O,TSO
187
188 REAL, DIMENSION( ims:ime, jms:jme ) , &
189 INTENT(INOUT) :: SOILT, &
190 HFX, &
191 QFX, &
192 LH, &
193 SFCEVP, &
194 SFCRUNOFF, &
195 UDRUNOFF, &
196 GRDFLX, &
197 ACSNOW, &
198 QVG, &
199 QCG, &
200 QSFC, &
201 QSG, &
202 CHKLOWQ, &
203 SOILT1, &
204 TSNAV
205
206 REAL, DIMENSION( ims:ime, jms:jme ) , &
207 INTENT(INOUT) :: SMAVAIL, &
208 SMMAX
209
210 REAL, DIMENSION( its:ite, jts:jte ) :: DEW, &
211 PC, &
212 RUNOFF1, &
213 RUNOFF2, &
214 EMISSL, &
215 ZNTL, &
216 LMAVAIL, &
217 SMELT, &
218 SNOH, &
219 SNFLX, &
220 SNOM, &
221 EDIR, &
222 EC, &
223 ETT, &
224 SUBLIM, &
225 sflx, &
226 EVAPL, &
227 PRCPL, &
228 XICED, &
229 INFILTR
230
231 !--- soil/snow properties
232 REAL, DIMENSION( ims:ime, 1:nsl, jms:jme) &
233 :: KEEPFR3DFLAG, &
234 SMFR3D
235
236 REAL &
237 :: RHOCS, &
238 RHOSN, &
239 BCLH, &
240 DQM, &
241 KSAT, &
242 PSIS, &
243 QMIN, &
244 QWRTZ, &
245 REF, &
246 WILT, &
247 CANWATR, &
248 SNOWFRAC, &
249 SNHEI, &
250 SNWE
251
252 REAL :: CN, &
253 SAT,CW, &
254 C1SN, &
255 C2SN, &
256 KQWRTZ, &
257 KICE, &
258 KWT
259
260
261 REAL, DIMENSION(1:NSL) :: ZSMAIN, &
262 ZSHALF, &
263 DTDZS2
264
265 REAL, DIMENSION(1:2*(nsl-2)) :: DTDZS
266
267 REAL, DIMENSION(1:4001) :: TBQ
268
269
270 REAL, DIMENSION( 1:nsl ) :: SOILM1D, &
271 TSO1D, &
272 SOILICE, &
273 SOILIQW, &
274 SMFRKEEP
275
276 REAL, DIMENSION( 1:nsl ) :: KEEPFR
277
278
279 REAL :: RSM, &
280 SNWEPRINT, &
281 SNHEIPRINT
282
283 REAL :: PRCPMS, &
284 NEWSNMS, &
285 PATM, &
286 TABS, &
287 QVATM, &
288 QCATM, &
289 Q2SAT, &
290 SATFLG, &
291 CONFLX, &
292 RHO, &
293 QKMS, &
294 TKMS, &
295 INFILTRP
296 REAL :: cq,r61,r273,arp,brp,x,evs,eis
297
298 INTEGER :: NROOT
299 INTEGER :: ILAND,ISOIL
300
301 INTEGER, DIMENSION ( 1:nvegclas ) :: IFOREST
302
303 INTEGER :: I,J,K,NZS,NZS1,NDDZS
304 INTEGER :: k1,l,k2,kp,km
305
306
307 !-----------------------------------------------------------------
308
309 NZS=NSL
310 NDDZS=2*(nzs-2)
311
312 !---- table TBQ is for resolution of balance equation in VILKA
313 CQ=173.15-.05
314 R273=1./273.15
315 R61=6.1153*0.62198
316 ARP=77455.*41.9/461.525
317 BRP=64.*41.9/461.525
318
319 DO K=1,4001
320 CQ=CQ+.05
321 ! TBQ(K)=R61*EXP(ARP*(R273-1./CQ)-BRP*LOG(CQ*R273))
322 EVS=EXP(17.67*(CQ-273.15)/(CQ-29.65))
323 EIS=EXP(22.514-6.15E3/CQ)
324 if(CQ.ge.273.15) then
325 ! tbq is in mb
326 tbq(k) = R61*evs
327 else
328 tbq(k) = R61*eis
329 endif
330
331 END DO
332
333 !--- Initialize soil/vegetation parameters
334 !--- This is temporary until SI is added to mass coordinate ---!!!!!
335
336 #if ( NMM_CORE == 1 )
337 if(ktau+1.eq.1) then
338 #else
339 if(ktau.eq.1) then
340 #endif
341 DO J=jts,jte
342 DO i=its,ite
343 do k=1,nsl
344 ! smfr3d (i,k,j)=soilmois(i,k,j)/900.*1.e3
345 keepfr3dflag(i,k,j)=0.
346 sh2o (i,k,j)=0.
347 enddo
348 !--- initializing to zero snow fraction
349 snowc(i,j) = min(1.,snowh(i,j)/0.1)
350 !--- initializing of snow temp
351 soilt1(i,j)=soilt(i,j)
352 tsnav(i,j) =0.5*(soilt(i,j)+tso(i,1,j))-273.
353 qcg (i,j) =0.
354 patm=P8w(i,kms,j)*1.e-2
355 QSG (i,j) = QSN(SOILT(i,j),TBQ)/PATM
356 qvg (i,j) = QSG(i,j)*mavail(i,j)
357 ! qvg (i,j) =qv3d(i,kms,j)
358 qsfc(i,j) = qsg(i,j)/(1.+qsg(i,j))
359 SMELT(i,j) = 0.
360 SNOM (i,j) = 0.
361 SNFLX(i,j) = 0.
362 DEW (i,j) = 0.
363 PC (i,j) = 0.
364 zntl (i,j) = 0.
365 RUNOFF1(i,j) = 0.
366 RUNOFF2(i,j) = 0.
367 emissl (i,j) = 0.
368 ! Temporarily!!!
369 canwat(i,j)=0.
370
371 ! For RUC LSM CHKLOWQ needed for MYJPBL should
372 ! 1 because is actual specific humidity at the surface, and
373 ! not the saturation value
374 chklowq(i,j) = 1.
375 infiltr(i,j) = 0.
376 snoh (i,j) = 0.
377 edir (i,j) = 0.
378 ec (i,j) = 0.
379 ett (i,j) = 0.
380 sublim(i,j) = 0.
381 sflx (i,j) = 0.
382 evapl (i,j) = 0.
383 prcpl (i,j) = 0.
384 ENDDO
385 ENDDO
386
387 do k=1,nsl
388 soilice(k)=0.
389 soiliqw(k)=0.
390 enddo
391 endif
392
393 !-----------------------------------------------------------------
394
395 PRCPMS = 0.
396 ! NROOT = 4
397
398
399 DO J=jts,jte
400
401 DO i=its,ite
402
403 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
404 print *,' IN LSMRUC ','ims,ime,jms,jme,its,ite,jts,jte,nzs', &
405 ims,ime,jms,jme,its,ite,jts,jte,nzs
406 print *,' IVGTYP, ISLTYP ', ivgtyp(i,j),isltyp(i,j)
407 print *,' MAVAIL ', mavail(i,j)
408 print *,' SOILT,QVG,P8w',soilt(i,j),qvg(i,j),p8w(i,1,j)
409 print *, 'LSMRUC, I,J,xland, QFX,HFX from SFCLAY',i,j,xland(i,j), &
410 qfx(i,j),hfx(i,j)
411 print *, ' GSW, GLW =',gsw(i,j),glw(i,j)
412 print *, 'SOILT, TSO start of time step =',soilt(i,j),(tso(i,k,j),k=1,nsl)
413 print *, 'SOILMOIS start of time step =',(soilmois(i,k,j),k=1,nsl)
414 print *, 'SMFROZEN start of time step =',(smfr3d(i,k,j),k=1,nsl)
415 print *, ' I,J=, after SFCLAY CHS,FLHC ',i,j,chs(i,j),flhc(i,j)
416 print *, 'LSMRUC, IVGTYP,ISLTYP,ZNT,ALB = ', ivgtyp(i,j),isltyp(i,j),znt(i,j),alb(i,j),i,j
417 print *, 'LSMRUC I,J,DT,RAINBL =',I,J,dt,RAINBL(i,j)
418 print *, 'XLAND ---->, ivgtype,isoiltyp,i,j',xland(i,j),ivgtyp(i,j),isltyp(i,j),i,j
419 ENDIF
420
421
422 ILAND = IVGTYP(i,j)
423 ISOIL = ISLTYP(I,J)
424 TABS = T3D(i,kms,j)
425 QVATM = QV3D(i,kms,j)
426 QCATM = QC3D(i,kms,j)
427 PATM = P8w(i,kms,j)*1.e-5
428 !---- what height is the first level?---- check!!!!!
429 !-- need to de-stagger from w levels to P levels
430 CONFLX = Z3D(i,kms,j)
431 ! CONFLX = 0.5*Z3D(i,kms,j)
432 ! CONFLX = 5.
433 RHO = RHO3D(I,kms,J)
434 !--- 1*e-3 is to convert from mm/s to m/s
435 IF(FRPCPN) THEN
436 PRCPMS = (RAINBL(i,j)/DT*1.e-3)*(1-FRZFRAC(I,J))
437 NEWSNMS = (RAINBL(i,j)/DT*1.e-3)*FRZFRAC(I,J)
438 ELSE
439 if (tabs.le.273.15) then
440 PRCPMS = 0.
441 NEWSNMS = RAINBL(i,j)/DT*1.e-3
442 else
443 PRCPMS = RAINBL(i,j)/DT*1.e-3
444 NEWSNMS = 0.
445 endif
446 ENDIF
447 !--- rooting depth is 5 levels for forests
448 ! if(iforest(ivgtyp(i,j)).eq.1) nroot=5
449 !--- convert exchange coeff to [m/s]
450 ! QKMS=FLQC(I,J)/RHO/MAVAIL(I,J)
451 if (myj) then
452 QKMS=CHS(i,j)
453 TKMS=CHS(i,j)
454 else
455 QKMS=FLQC(I,J)/RHO/MAVAIL(I,J)
456 TKMS=FLHC(I,J)/RHO/CP
457 endif
458 !--- convert incoming snow and canwat from mm to m
459 SNWE=SNOW(I,J)*1.E-3
460 SNHEI=SNOWH(I,J)
461 CANWATR=CANWAT(I,J)*1.E-3
462 SNOWFRAC=SNOWC(I,J)
463
464 !-----
465 zsmain(1)=0.
466 zshalf(1)=0.
467 do k=2,nzs
468 zsmain(k)= zs(k)
469 zshalf(k)=0.5*(zsmain(k-1) + zsmain(k))
470 enddo
471
472 !-- definition of number of soil levels in the rooting zone
473 IF(iforest(ivgtyp(i,j)).ne.1) THEN
474 !---- non forests and deciduous forests
475 do k=2,nzs
476 if(zsmain(k).ge.0.4) then
477 NROOT=K
478 goto 111
479 endif
480 enddo
481
482 ELSE
483 !---- evergreen and mixed forests
484 do k=2,nzs
485 if(zsmain(k).ge.1.1) then
486 NROOT=K
487 goto 111
488 endif
489 enddo
490 ENDIF
491 111 continue
492
493 !-----
494 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
495 print *,' ZS, ZSMAIN, ZSHALF, CONFLX --->', zs,zsmain,zshalf,conflx
496 print *,'NROOT, iforest, ivgtyp, i,j ', nroot,iforest(ivgtyp(i,j)),ivgtyp(I,J),I,J
497 ENDIF
498
499 !------------------------------------------------------------
500 !----- DDZS and DSDZ1 are for implicit soilution of soil eqns.
501 !-------------------------------------------------------------
502 NZS1=NZS-1
503 !-----
504 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
505 print *,' DT,NZS1, ZSMAIN, ZSHALF --->', dt,nzs1,zsmain,zshalf
506 ENDIF
507
508 DO K=2,NZS1
509 K1=2*K-3
510 K2=K1+1
511 X=DT/2./(ZSHALF(K+1)-ZSHALF(K))
512 DTDZS(K1)=X/(ZSMAIN(K)-ZSMAIN(K-1))
513 DTDZS2(K-1)=X
514 DTDZS(K2)=X/(ZSMAIN(K+1)-ZSMAIN(K))
515 END DO
516
517
518 CN=0.5 ! exponent
519 SAT=0.0004 ! canopy water saturated
520
521 CW =4.183E6
522
523
524 !--- Constants used in Johansen soil thermal
525 !--- conductivity method
526
527 KQWRTZ=7.7
528 KICE=2.2
529 KWT=0.57
530
531 !***********************************************************************
532 !--- Constants for snow density calculations C1SN and C2SN
533
534 c1sn=0.026
535 ! c1sn=0.01
536 c2sn=21.
537
538 !***********************************************************************
539
540 NROOT= 4
541 ! ! rooting depth
542
543 if(SNOWH(i,j).gt.0.) then
544 RHOSN = SNOW(i,j)/SNOWH(i,j)
545 else
546 RHOSN = 300.
547 endif
548
549 !--- initializing soil and surface properties
550 CALL SOILVEGIN ( ILAND,ISOIL,MYJ,IFOREST, &
551 EMISSL(I,J),PC(i,j),ZNT(I,J),QWRTZ, &
552 ! EMISSL(I,J),PC(i,j),ZNTL(I,J),QWRTZ, &
553 RHOCS,BCLH,DQM,KSAT,PSIS,QMIN,REF,WILT )
554
555 !*** SET ZERO-VALUE FOR SOME OUTPUT DIAGNOSTIC ARRAYS
556
557
558 IF((XLAND(I,J)-1.5).GE.0.)THEN
559 !-- Water point
560 SMAVAIL(I,J)=1.0
561 SMMAX(I,J)=1.0
562 ! SNOW(I,J)=0.0
563 LMAVAIL(I,J)=1.0
564
565 ILAND=16
566 ISOIL=14
567
568 patm=P8w(i,kms,j)*1.e-2
569 qvg (i,j) = QSN(SOILT(i,j),TBQ)/PATM
570 qsfc(i,j) = qvg(i,j)/(1.+qvg(i,j))
571 CHKLOWQ(I,J)=1.
572 Q2SAT=QSN(TABS,TBQ)/PATM
573
574 DO K=1,NZS
575 SOILMOIS(I,K,J)=1.0
576 TSO(I,K,J)= SOILT(I,J)
577 ENDDO
578
579 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
580 PRINT*,' water point, I=',I, &
581 'J=',J, 'SOILT=', SOILT(i,j)
582 ENDIF
583 !--- decide if this water point is ice:
584 ! if(tabs.le.271.) then
585 if(xice(i,j).gt.0.5) then
586 ! if(soilt(i,j).le.271.or.xice(i,j).eq.1.) then
587 ! if(tabs.le.271.or.xice(i,j).eq.1.) then
588 XICED(i,j)=1.
589 else
590 XICED(i,j)=0.
591 endif
592
593 IF(XICED(I,J).NE.1.) SNOW(I,J)=0.
594 IF(XICED(I,J).GT.0.5)THEN
595 !-- Sea-ice case
596 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
597 PRINT*,' sea-ice at water point, I=',I, &
598 'J=',J
599 ENDIF
600 ILAND = 24
601 ISOIL = 16
602
603 SMAVAIL(I,J)=1.0
604 SMMAX(I,J)=1.0
605 LMAVAIL(I,J)=1.0
606 ! SOILT(I,J) = MIN(273.15,SOILT(I,J))
607
608 DO K=1,NZS
609 SOILMOIS(I,K,J)=1.0
610 TSO(I,K,J)= MIN(273.15,SOILT(I,J))
611 ENDDO
612 ENDIF
613
614 ! for MYJ surface and PBL scheme
615 if (myj) then
616 IF((QVATM.GE.Q2SAT*0.95).AND.QVATM.LT.qvg(I,J))THEN
617 ! IF((QVATM.GE.Q2SAT*0.95).AND.QVATM.LT.qsg(I,J))THEN
618 SATFLG=0.
619 ELSE
620 SATFLG=1.0
621 ENDIF
622 else
623 SATFLG=1.0
624 endif
625 QFX(I,J)=QFX(I,J)*SATFLG
626
627
628 ELSE
629
630 !-- Land point
631 ! Attention!!!! RUC LSM uses soil moisture content minus residual (minimum
632 ! soil moisture content for a given soil type) as a state variable.
633 ! If the WRF model is initialized from the RUC background model, then the
634 ! soil moisture variable is consistent with the RUC LSM.
635 ! If the WRF model is initialized from another background model (ETA, GFS...)
636 ! then the residual value should be subtracted when the 1-d array of soil
637 ! moisture is initialized before the call to SFCTMP, and after SFCTMP qmin
638 ! should be added back in.
639 !
640 ! soilm1d (k) = min(max(0.,soilmois(i,k,j)-qmin(i,j)),dqm(i,j))
641
642
643 DO k=1,nzs
644 ! soilm1d - soil moisture content minus residual [m**3/m**3]
645 soilm1d (k) = min(max(0.,soilmois(i,k,j)),dqm)
646 tso1d (k) = tso(i,k,j)
647 ENDDO
648
649 do k=1,nzs
650 smfrkeep(k) = smfr3d(i,k,j)
651 keepfr (k) = keepfr3dflag(i,k,j)
652 enddo
653
654 ! LMAVAIL(I,J)=max(0.00001,min(1.,soilmois(i,1,j)/(REF-QMIN)))
655 LMAVAIL(I,J)=max(0.00001,min(1.,soilmois(i,1,j)/dqm))
656
657 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
658 print *,'LAND, i,j,tso1d,soilm1d,PATM,TABS,QVATM,QCATM,RHO', &
659 i,j,tso1d,soilm1d,PATM,TABS,QVATM,QCATM,RHO
660 print *,'CONFLX =',CONFLX
661 print *,'SMFRKEEP,KEEPFR ',SMFRKEEP,KEEPFR
662 ENDIF
663
664 !-----------------------------------------------------------------
665 CALL SFCTMP (dt,ktau,conflx,i,j, &
666 !--- input variables
667 nzs,nddzs,nroot, &
668 iland,isoil,xland(i,j),ivgtyp(i,j), &
669 PRCPMS,NEWSNMS,SNWE,SNHEI,SNOWFRAC,RHOSN, &
670 PATM,TABS,QVATM,QCATM,RHO, &
671 GLW(I,J),GSW(I,J),EMISSL(I,J), &
672 QKMS,TKMS,PC(I,J),LMAVAIL(I,J), &
673 canwatr,vegfra(I,J),alb(I,J),znt(I,J), &
674 snoalb(i,j),albbck(i,j), & !new
675 myj, &
676 !--- soil fixed fields
677 QWRTZ, &
678 rhocs,dqm,qmin,ref, &
679 wilt,psis,bclh,ksat, &
680 sat,cn,zsmain,zshalf,DTDZS,DTDZS2,tbq, &
681 !--- constants
682 cp,g0,lv,stbolt,cw,c1sn,c2sn, &
683 KQWRTZ,KICE,KWT, &
684 !--- output variables
685 snweprint,snheiprint,rsm, &
686 soilm1d,tso1d,smfrkeep,keepfr, &
687 soilt(I,J),soilt1(i,j),tsnav(i,j),dew(I,J), &
688 qvg(I,J),qsg(I,J),qcg(I,J),SMELT(I,J), &
689 SNOH(I,J),SNFLX(I,J),SNOM(I,J),ACSNOW(I,J), &
690 edir(I,J),ec(I,J),ett(I,J),sfcevp(I,J), &
691 lh(I,J),hfx(I,J),sflx(I,J),sublim(I,J), &
692 evapl(I,J),prcpl(I,J),runoff1(I,J), &
693 runoff2(I,J),soilice,soiliqw,infiltrp)
694 !-----------------------------------------------------------------
695
696 !*** DIAGNOSTICS
697 !--- available and maximum soil moisture content in the soil
698 !--- domain
699 smavail(i,j) = 0.
700 smmax (i,j) = 0.
701
702 do k=1,nzs-1
703 smavail(i,j)=smavail(i,j)+(qmin+soilm1d(k))* &
704 (zshalf(k+1)-zshalf(k))
705 smmax (i,j) =smmax (i,j)+(qmin+dqm)* &
706 (zshalf(k+1)-zshalf(k))
707 enddo
708
709 smavail(i,j)=smavail(i,j)+(qmin+soilm1d(nzs))* &
710 (zsmain(nzs)-zshalf(nzs))
711 smmax (i,j) =smmax (i,j)+(qmin+dqm)* &
712 (zsmain(nzs)-zshalf(nzs))
713
714 !--- Convert the water unit into mm
715 SFCRUNOFF(I,J) = SFCRUNOFF(I,J)+RUNOFF1(I,J)*DT*1000.0
716 UDRUNOFF (I,J) = UDRUNOFF(I,J)+RUNOFF2(I,J)*1000.0
717 SMAVAIL (I,J) = SMAVAIL(I,J) * 1000.
718 SMMAX (I,J) = SMMAX(I,J) * 1000.
719 SFCEXC (I,J) = TKMS
720 ! SFCLAY expects QSFC as saturation mixing ration at surface
721 ! QSFC(I,J) = QSG(I,J)
722 ! MYJSFC expects QSFC as actual specific humidity at the surface
723 QSFC(I,J) = QVG(I,J)/(1.+QVG(I,J))
724 patm=P8w(i,kms,j)*1.e-2
725 Q2SAT=QSN(TABS,TBQ)/PATM
726 ! for MYJ surface and PBL scheme
727 if (myj) then
728 IF((QVATM.GE.Q2SAT*0.95).AND.QVATM.LT.qvg(I,J))THEN
729 CHKLOWQ(I,J)=0.
730 ELSE
731 CHKLOWQ(I,J)=1.
732 ENDIF
733 else
734 CHKLOWQ(I,J)=1.
735 endif
736
737 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
738 if(CHKLOWQ(I,J).eq.0.) then
739 print *,'i,j,CHKLOWQ', &
740 i,j,CHKLOWQ(I,J)
741 endif
742 ENDIF
743
744 MAVAIL (i,j) = LMAVAIL(I,J)
745 ! SNOW is in [mm], SNWE is in [m]; CANWAT is in mm, CANWATR is in m
746 SNOW (i,j) = SNWE*1000.
747 SNOWH (I,J) = SNHEI
748 CANWAT (I,J) = CANWATR*1000.
749
750 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
751 print *,' LAND, I=,J=, QFX, HFX after SFCTMP', i,j,lh(i,j),hfx(i,j)
752 ENDIF
753 QFX (I,J) = LH(I,J)/LV
754 GRDFLX (I,J) = -1. * sflx(I,J)
755
756 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
757 print *,' QFX after change, LH ', i,j, QFX(i,j),LH(I,J)
758 ENDIF
759 !--- SNOWC snow cover flag
760 SNOWC(I,J)=SNOWFRAC
761
762 ! IF(SNOWH(I,J).GT.0.02)THEN
763 ! SNOWC(I,J)=1.0
764 ! ELSE
765 ! SNOWC(I,J)=0.0
766 ! ENDIF
767
768 INFILTR(I,J) = INFILTRP
769
770 !--- get 3d soil fields
771 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
772 print *,'LAND, i,j,tso1d,soilm1d - end of time step', &
773 i,j,tso1d,soilm1d
774 ENDIF
775
776 do k=1,nzs
777
778 soilmois(i,k,j) = soilm1d(k)
779 sh2o (i,k,j) = soiliqw(k)
780 tso(i,k,j) = tso1d(k)
781 enddo
782
783 do k=1,nzs
784 smfr3d(i,k,j) = smfrkeep(k)
785 keepfr3dflag(i,k,j) = keepfr (k)
786 enddo
787
788 ENDIF
789
790 ENDDO
791
792 ENDDO
793
794 !-----------------------------------------------------------------
795 END SUBROUTINE LSMRUC
796 !-----------------------------------------------------------------
797
798
799
800 SUBROUTINE SFCTMP (delt,ktau,conflx,i,j, &
801 !--- input variables
802 nzs,nddzs,nroot, &
803 ILAND,ISOIL,XLAND,IVGTYP, &
804 PRCPMS,NEWSNMS,SNWE,SNHEI,SNOWFRAC,RHOSN, &
805 PATM,TABS,QVATM,QCATM,rho, &
806 GLW,GSW,EMISS,QKMS,TKMS,PC, &
807 MAVAIL,CST,VEGFRA,ALB,ZNT, &
808 ALB_SNOW,ALB_SNOW_FREE, &
809 MYJ, &
810 !--- soil fixed fields
811 QWRTZ,rhocs,dqm,qmin,ref,wilt,psis,bclh,ksat, &
812 sat,cn,zsmain,zshalf,DTDZS,DTDZS2,tbq, &
813 !--- constants
814 cp,g0,lv,stbolt,cw,c1sn,c2sn, &
815 KQWRTZ,KICE,KWT, &
816 !--- output variables
817 snweprint,snheiprint,rsm, &
818 soilm1d,ts1d,smfrkeep,keepfr,soilt,soilt1, &
819 tsnav,dew,qvg,qsg,qcg, &
820 SMELT,SNOH,SNFLX,SNOM,ACSNOW, &
821 edir1,ec1,ett1,eeta,qfx,hfx,s,sublim, &
822 evapl,prcpl,runoff1,runoff2,soilice, &
823 soiliqw,infiltr)
824 !-----------------------------------------------------------------
825 IMPLICIT NONE
826 !-----------------------------------------------------------------
827
828 !--- input variables
829
830 INTEGER, INTENT(IN ) :: i,j,nroot,ktau,nzs , &
831 nddzs !nddzs=2*(nzs-2)
832
833 REAL, INTENT(IN ) :: DELT,CONFLX
834 REAL, INTENT(IN ) :: C1SN,C2SN
835 LOGICAL, INTENT(IN ) :: myj
836 !--- 3-D Atmospheric variables
837 REAL , &
838 INTENT(IN ) :: PATM, &
839 TABS, &
840 QVATM, &
841 QCATM
842 REAL , &
843 INTENT(IN ) :: GLW, &
844 GSW, &
845 PC, &
846 ALB_SNOW, &
847 ALB_SNOW_FREE, &
848 VEGFRA, &
849 XLAND, &
850 RHO, &
851 QKMS, &
852 TKMS
853
854 INTEGER, INTENT(IN ) :: IVGTYP
855 !--- 2-D variables
856 REAL , &
857 INTENT(INOUT) :: EMISS, &
858 MAVAIL, &
859 SNOWFRAC, &
860 ALB, &
861 CST
862
863 !--- soil properties
864 REAL :: &
865 RHOCS, &
866 BCLH, &
867 DQM, &
868 KSAT, &
869 PSIS, &
870 QMIN, &
871 QWRTZ, &
872 REF, &
873 SAT, &
874 WILT
875
876 REAL, INTENT(IN ) :: CN, &
877 CW, &
878 CP, &
879 G0, &
880 LV, &
881 STBOLT, &
882 KQWRTZ, &
883 KICE, &
884 KWT
885
886 REAL, DIMENSION(1:NZS), INTENT(IN) :: ZSMAIN, &
887 ZSHALF, &
888 DTDZS2
889
890 REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS
891
892 REAL, DIMENSION(1:4001), INTENT(IN) :: TBQ
893
894
895 !--- input/output variables
896 !-------- 3-d soil moisture and temperature
897 REAL, DIMENSION( 1:nzs ) , &
898 INTENT(INOUT) :: TS1D, &
899 SOILM1D, &
900 SMFRKEEP
901 REAL, DIMENSION( 1:nzs ) , &
902 INTENT(INOUT) :: KEEPFR
903
904
905 INTEGER, INTENT(INOUT) :: ILAND,ISOIL
906
907 !-------- 2-d variables
908 REAL , &
909 INTENT(INOUT) :: DEW, &
910 EDIR1, &
911 EC1, &
912 ETT1, &
913 EETA, &
914 EVAPL, &
915 INFILTR, &
916 RHOSN, &
917 SUBLIM, &
918 PRCPL, &
919 QVG, &
920 QSG, &
921 QCG, &
922 QFX, &
923 HFX, &
924 S, &
925 RUNOFF1, &
926 RUNOFF2, &
927 ACSNOW, &
928 SNWE, &
929 SNHEI, &
930 SMELT, &
931 SNOM, &
932 SNOH, &
933 SNFLX, &
934 SOILT, &
935 SOILT1, &
936 TSNAV, &
937 ZNT
938
939 !-------- 1-d variables
940 REAL, DIMENSION(1:NZS), INTENT(OUT) :: SOILICE, &
941 SOILIQW
942
943 REAL, INTENT(OUT) :: RSM, &
944 SNWEPRINT, &
945 SNHEIPRINT
946 !--- Local variables
947
948 INTEGER :: K,ILNB
949
950 REAL :: BSN, XSN, RHONEWSN , &
951 RAINF, SNTH, NEWSN, PRCPMS, NEWSNMS , &
952 T3, UPFLUX, XINET
953 REAL :: snhei_crit, keep_snow_albedo
954
955 REAL :: RNET,GSWNEW,EMISSN,ALBSN,ZNTSN
956 REAL :: VEGFRAC
957
958 !-----------------------------------------------------------------
959 integer, parameter :: ilsnow=99
960
961 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
962 print *,' in SFCTMP',i,j,nzs,nddzs,nroot, &
963 SNWE,RHOSN,SNOM,SMELT,TS1D
964 ENDIF
965 ! print *,' in SFCTMP',i,j,nzs,nddzs,nroot, &
966 ! IVGTYP,ISOIL,ILAND, &
967 ! PRCPMS,SNWE,RHOSN, &
968 ! PATM,TABS,QVATM,QCATM,rho
969 ! GLW,GSW,EMISS,QKMS,TKMS,PC, &
970 ! cst,vegfrac,alb,znt, &
971 !--- soil fixed fields
972 ! QWRTZ,rhocs,dqm,qmin,ref,wilt,psis,bclh,ksat, &
973 ! sat,cn,zsmain,zshalf,DTDZS,DTDZS2,tbq, &
974 !--- constants
975 ! cp,g0,lv,stbolt,cw,c1sn,c2sn, &
976 ! KQWRTZ,KICE,KWT
977
978 NEWSN=0.
979 RAINF = 0.
980 RSM=0.
981 INFILTR=0.
982 VEGFRAC=0.01*VEGFRA
983
984 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
985 print *,'I,J,KTAU,QKMS,TKMS', i,j,ktau,qkms,tkms
986 print *,'GSW, GLW, SOILT, STBOLT, EMISS', &
987 GSW, GLW, SOILT, STBOLT, EMISS
988 ENDIF
989
990
991 SNHEI = SNWE * 1000. / RHOSN
992 !--------------
993 T3 = STBOLT*SOILT*SOILT*SOILT
994 UPFLUX = T3 *SOILT
995 XINET = EMISS*(GLW-UPFLUX)
996 RNET = GSW + XINET
997
998 !Calculate the amount (m) of fresh snow
999
1000 if(snhei.gt.0.0081*1.e3/rhosn) then
1001 !*** Correct snow density for current temperature (Koren et al. 1999)
1002 BSN=delt/3600.*c1sn*exp(0.08*tsnav-c2sn*rhosn*1.e-3)
1003 if(bsn*snwe*100..lt.1.e-4) goto 777
1004 XSN=rhosn*(exp(bsn*snwe*100.)-1.)/(bsn*snwe*100.)
1005 rhosn=MIN(MAX(100.,XSN),400.)
1006 ! rhosn=MIN(MAX(50.,XSN),400.)
1007 777 continue
1008
1009 else
1010 rhosn =200.
1011 rhonewsn =100.
1012 endif
1013
1014 ! IF(TABS.LE.273.15)THEN
1015
1016 newsn=newsnms*delt
1017 !--- consider for now that all PRCPMS went into snow
1018 ! prcpms = 0.
1019 !---- ACSNOW - accumulation of snow water [m]
1020 acsnow=acsnow+newsn
1021
1022 IF(NEWSN.GE.1.E-8) THEN
1023 !*** Calculate fresh snow density (t > -15C, else MIN value)
1024 !*** Eq. 10 from Koren et al. (1999)
1025
1026 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
1027 print *, 'THERE IS NEW SNOW, newsn', newsn
1028 ENDIF
1029 if(tabs.lt.258.15) then
1030 ! rhonewsn=50.
1031 rhonewsn=100.
1032
1033 else
1034 rhonewsn=1.e3*max((0.10+0.0017*(Tabs-273.15+15.)**1.5) &
1035 , 0.10)
1036 ! rhonewsn=1.e3*max((0.05+0.0017*(Tabs-273.15+15.)**1.5) &
1037 ! , 0.05)
1038 rhonewsn=MIN(rhonewsn,400.)
1039 ! rhonewsn=100.
1040 endif
1041
1042
1043 !*** Define average snow density of the snow pack considering
1044 !*** the amount of fresh snow (eq. 9 in Koren et al.(1999)
1045 !*** without snow melt )
1046 xsn=(rhosn*snwe+rhonewsn*newsn)/ &
1047 (snwe+newsn)
1048 rhosn=MIN(MAX(100.,XSN),400.)
1049 ! rhosn=MIN(MAX(50.,XSN),400.)
1050
1051 snwe=snwe+newsn
1052 snhei=snwe*1.E3/rhosn
1053 NEWSN=NEWSN*1.E3/rhosn
1054 endif
1055
1056 ! ELSE
1057 !--- TABS is above freezing. Needed precip rates from microphysics
1058 !--- to do a better job with mixed phase precip.
1059
1060 ! NEWSN = 0.
1061 !
1062 ! ENDIF
1063
1064 IF(PRCPMS.NE.0.) THEN
1065
1066 ! PRCPMS is liquid precipitation rate
1067 ! RAINF is a flag used for calculation of rain water
1068 ! heat content contribution into heat budget equation. Rain's temperature
1069 ! is set equal to air temperature at the first atmospheric
1070 ! level.
1071
1072 RAINF=1.
1073 ENDIF
1074
1075 ! IF((XLAND-1.5).GE.0.)THEN
1076 ! IF(ILAND.EQ.16) THEN
1077 ! SNHEI=0.
1078 ! SNWE=0.
1079 ! ELSE
1080
1081 IF(SNHEI.GT.0.0) THEN
1082 !--- Set of surface parameters should be changed to snow values for grid
1083 !--- points where the snow cover exceeds snow threshold of 2 cm
1084 EMISS = 0.91
1085
1086 ! GSWNEW = GSW
1087 ! The following lines compute albedo depending on snow
1088 ! depth. For now commented out.
1089 ! alb_snow_free=0.2
1090 ! alb_snow=0.70
1091 ! SNHEI_CRIT=0.05
1092
1093 SNHEI_CRIT=0.01601*1.e3/rhosn
1094 SNOWFRAC=MIN(1.,SNHEI/(2.*SNHEI_CRIT))
1095
1096 KEEP_SNOW_ALBEDO = 0.
1097 IF (NEWSN.GT.0.) KEEP_SNOW_ALBEDO = 1.
1098
1099 !--- GSW in-coming solar
1100 GSWNEW=GSW/(1.-ALB)
1101
1102 ALB = MAX(keep_snow_albedo*alb_snow, &
1103 MIN((alb_snow_free + &
1104 (alb_snow - alb_snow_free) * &
1105 (snhei/(2.*SNHEI_CRIT))), alb_snow))
1106 !--- recompute absorbed solar radiation and net radiation
1107 !--- for new value of albedo
1108 gswnew=gswnew*(1.-alb)
1109
1110 XINET = EMISS*(GLW-UPFLUX)
1111 RNET = GSWnew + XINET
1112
1113 CALL SNOWSOIL ( & !--- input variables
1114 i,j,isoil,delt,ktau,conflx,nzs,nddzs,nroot, &
1115 ILAND,PRCPMS,RAINF,NEWSN,snhei,SNWE,snowfrac, &
1116 RHOSN,PATM,QVATM,QCATM, &
1117 GLW,GSWnew,EMISS,RNET,IVGTYP, &
1118 QKMS,TKMS,PC,CST, &
1119 RHO,VEGFRAC,ALB,ZNT, &
1120 MYJ, &
1121 !--- soil fixed fields
1122 QWRTZ,rhocs,dqm,qmin,ref,wilt,psis,bclh,ksat, &
1123 sat,cn,zsmain,zshalf,DTDZS,DTDZS2,tbq, &
1124 !--- constants
1125 lv,CP,G0,cw,stbolt,tabs, &
1126 KQWRTZ,KICE,KWT, &
1127 !--- output variables
1128 ilnb,snweprint,snheiprint,rsm, &
1129 soilm1d,ts1d,smfrkeep,keepfr, &
1130 dew,soilt,soilt1,tsnav,qvg,qsg,qcg, &
1131 SMELT,SNOH,SNFLX,SNOM,edir1,ec1,ett1,eeta, &
1132 qfx,hfx,s,sublim,prcpl,runoff1,runoff2, &
1133 mavail,soilice,soiliqw,infiltr )
1134
1135 if(snhei.eq.0.) then
1136
1137 ! if(snhei.le.2.e-2) then
1138 !--- all snow is melted
1139 ! gswnew=gswnew/(1.-alb)
1140 alb=alb_snow_free
1141 ! gswnew=gswnew*(1.-alb)
1142 endif
1143
1144 ELSE
1145
1146 snheiprint=0.
1147 snweprint=0.
1148
1149 CALL SOIL( &
1150 !--- input variables
1151 i,j,iland,isoil,delt,ktau,conflx,nzs,nddzs,nroot, &
1152 PRCPMS,RAINF,PATM,QVATM,QCATM,GLW,GSW, &
1153 EMISS,RNET,QKMS,TKMS,PC,cst,rho,vegfrac, &
1154 !--- soil fixed fields
1155 QWRTZ,rhocs,dqm,qmin,ref,wilt, &
1156 psis,bclh,ksat,sat,cn, &
1157 zsmain,zshalf,DTDZS,DTDZS2,tbq, &
1158 !--- constants
1159 lv,CP,G0,cw,stbolt,tabs, &
1160 KQWRTZ,KICE,KWT, &
1161 !--- output variables
1162 soilm1d,ts1d,smfrkeep,keepfr, &
1163 dew,soilt,qvg,qsg,qcg,edir1,ec1, &
1164 ett1,eeta,qfx,hfx,s,evapl,prcpl,runoff1, &
1165 runoff2,mavail,soilice,soiliqw, &
1166 infiltr)
1167
1168 ENDIF
1169 ! ENDIF
1170
1171 !
1172
1173 ! RETURN
1174 ! END
1175 !---------------------------------------------------------------
1176 END SUBROUTINE SFCTMP
1177 !---------------------------------------------------------------
1178
1179
1180 FUNCTION QSN(TN,T)
1181 !****************************************************************
1182 REAL, DIMENSION(1:4001), INTENT(IN ) :: T
1183 REAL, INTENT(IN ) :: TN
1184
1185 REAL QSN, R,R1,R2
1186 INTEGER I
1187
1188 R=(TN-173.15)/.05+1.
1189 I=INT(R)
1190 IF(I.GE.1) goto 10
1191 I=1
1192 R=1.
1193 10 IF(I.LE.4000) GOTO 20
1194 I=4000
1195 R=4001.
1196 20 R1=T(I)
1197 R2=R-I
1198 QSN=(T(I+1)-R1)*R2 + R1
1199 ! print *,' in QSN, I,R,R1,R2,T(I+1),TN, QSN', I,R,r1,r2,t(i+1),tn,QSN
1200 ! RETURN
1201 ! END
1202 !-----------------------------------------------------------------------
1203 END FUNCTION QSN
1204 !------------------------------------------------------------------------
1205
1206
1207 SUBROUTINE SOIL ( &
1208 !--- input variables
1209 i,j,iland,isoil,delt,ktau,conflx,nzs,nddzs,nroot,&
1210 PRCPMS,RAINF,PATM,QVATM,QCATM, &
1211 GLW,GSW,EMISS,RNET, &
1212 QKMS,TKMS,PC,cst,rho,vegfrac, &
1213 !--- soil fixed fields
1214 QWRTZ,rhocs,dqm,qmin,ref,wilt,psis,bclh,ksat, &
1215 sat,cn,zsmain,zshalf,DTDZS,DTDZS2,tbq, &
1216 !--- constants
1217 xlv,CP,G0_P,cw,stbolt,TABS, &
1218 KQWRTZ,KICE,KWT, &
1219 !--- output variables
1220 soilmois,tso,smfrkeep,keepfr, &
1221 dew,soilt,qvg,qsg,qcg, &
1222 edir1,ec1,ett1,eeta,qfx,hfx,s,evapl, &
1223 prcpl,runoff1,runoff2,mavail,soilice, &
1224 soiliqw,infiltrp)
1225
1226 !*************************************************************
1227 ! Energy and moisture budget for vegetated surfaces
1228 ! without snow, heat diffusion amf Richards eqns. in
1229 ! soil
1230 !
1231 ! DELT - time step (s)
1232 ! ktau - numver of time step
1233 ! CONFLX - depth of constant flux layer (m)
1234 ! J,I - the location of grid point
1235 ! IME, JME, KME, NZS - dimensions of the domain
1236 ! NROOT - number of levels within the root zone
1237 ! PRCPMS - precipitation rate in m/s
1238 ! PATM - pressure [bar]
1239 ! QVATM,QCATM - cloud and water vapor mixing ratio (kg/kg)
1240 ! at the first atm. level
1241 ! GLW, GSW - incoming longwave and absorbed shortwave
1242 ! radiation at the surface (W/m^2)
1243 ! EMISS,RNET - emissivity of the ground surface (0-1) and net
1244 ! radiation at the surface (W/m^2)
1245 ! QKMS - exchange coefficient for water vapor in the
1246 ! surface layer (m/s)
1247 ! TKMS - exchange coefficient for heat in the surface
1248 ! layer (m/s)
1249 ! PC - plant coefficient (resistance) (0-1)
1250 ! RHO - density of atmosphere near sueface (kg/m^3)
1251 ! VEGFRAC - greeness fraction
1252 ! RHOCS - volumetric heat capacity of dry soil
1253 ! DQM, QMIN - porosity minus residual soil moisture QMIN (m^3/m^3)
1254 ! REF, WILT - field capacity soil moisture and the
1255 ! wilting point (m^3/m^3)
1256 ! PSIS - matrix potential at saturation (m)
1257 ! BCLH - exponent for Clapp-Hornberger parameterization
1258 ! KSAT - saturated hydraulic conductivity (m/s)
1259 ! SAT - maximum value of water intercepted by canopy (m)
1260 ! CN - exponent for calculation of canopy water
1261 ! ZSMAIN - main levels in soil (m)
1262 ! ZSHALF - middle of the soil layers (m)
1263 ! DTDZS,DTDZS2 - dt/(2.*dzshalf*dzmain) and dt/dzshalf in soil
1264 ! TBQ - table to define saturated mixing ration
1265 ! of water vapor for given temperature and pressure
1266 ! SOILMOIS,TSO - soil moisture (m^3/m^3) and temperature (K)
1267 ! DEW - dew in kg/m^2s
1268 ! SOILT - skin temperature (K)
1269 ! QSG,QVG,QCG - saturated mixing ratio, mixing ratio of
1270 ! water vapor and cloud at the ground
1271 ! surface, respectively (kg/kg)
1272 ! EDIR1, EC1, ETT1, EETA - direct evaporation, evaporation of
1273 ! canopy water, transpiration in kg m-2 s-1 and total
1274 ! evaporation in m s-1.
1275 ! QFX, HFX - latent and sensible heat fluxes (W/m^2)
1276 ! S - soil heat flux in the top layer (W/m^2)
1277 ! RUNOFF - surface runoff (m/s)
1278 ! RUNOFF2 - underground runoff (m)
1279 ! MAVAIL - moisture availability in the top soil layer (0-1)
1280 ! INFILTRP - infiltration flux from the top of soil domain (m/s)
1281 !
1282 !*****************************************************************
1283 IMPLICIT NONE
1284 !-----------------------------------------------------------------
1285
1286 !--- input variables
1287
1288 INTEGER, INTENT(IN ) :: nroot,ktau,nzs , &
1289 nddzs !nddzs=2*(nzs-2)
1290 INTEGER, INTENT(IN ) :: i,j,iland,isoil
1291 REAL, INTENT(IN ) :: DELT,CONFLX
1292 !--- 3-D Atmospheric variables
1293 REAL, &
1294 INTENT(IN ) :: PATM, &
1295 QVATM, &
1296 QCATM
1297 !--- 2-D variables
1298 REAL, &
1299 INTENT(IN ) :: GLW, &
1300 GSW, &
1301 EMISS, &
1302 RHO, &
1303 PC, &
1304 VEGFRAC, &
1305 QKMS, &
1306 TKMS
1307
1308 !--- soil properties
1309 REAL, &
1310 INTENT(IN ) :: RHOCS, &
1311 BCLH, &
1312 DQM, &
1313 KSAT, &
1314 PSIS, &
1315 QMIN, &
1316 QWRTZ, &
1317 REF, &
1318 WILT
1319
1320 REAL, INTENT(IN ) :: CN, &
1321 CW, &
1322 KQWRTZ, &
1323 KICE, &
1324 KWT, &
1325 XLV, &
1326 g0_p
1327
1328
1329 REAL, DIMENSION(1:NZS), INTENT(IN) :: ZSMAIN, &
1330 ZSHALF, &
1331 DTDZS2
1332
1333 REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS
1334
1335 REAL, DIMENSION(1:4001), INTENT(IN) :: TBQ
1336
1337
1338 !--- input/output variables
1339 !-------- 3-d soil moisture and temperature
1340 REAL, DIMENSION( 1:nzs ) , &
1341 INTENT(INOUT) :: TSO, &
1342 SOILMOIS, &
1343 SMFRKEEP
1344
1345 REAL, DIMENSION( 1:nzs ) , &
1346 INTENT(INOUT) :: KEEPFR
1347
1348 !-------- 2-d variables
1349 REAL, &
1350 INTENT(INOUT) :: DEW, &
1351 CST, &
1352 EDIR1, &
1353 EC1, &
1354 ETT1, &
1355 EETA, &
1356 EVAPL, &
1357 PRCPL, &
1358 MAVAIL, &
1359 QVG, &
1360 QSG, &
1361 QCG, &
1362 RNET, &
1363 QFX, &
1364 HFX, &
1365 S, &
1366 SAT, &
1367 RUNOFF1, &
1368 RUNOFF2, &
1369 SOILT
1370
1371 !-------- 1-d variables
1372 REAL, DIMENSION(1:NZS), INTENT(OUT) :: SOILICE, &
1373 SOILIQW
1374
1375 !--- Local variables
1376
1377 REAL :: INFILTRP, transum , &
1378 RAINF, PRCPMS , &
1379 TABS, T3, UPFLUX, XINET
1380 REAL :: CP,G0,LV,STBOLT,xlmelt,dzstop , &
1381 can,epot,fac,fltot,ft,fq,hft , &
1382 q1,ras,rhoice,sph , &
1383 trans,zn,ci,cvw,tln,tavln,pi , &
1384 DD1,CMC2MS,DRYCAN,WETCAN , &
1385 INFMAX,RIW
1386 REAL, DIMENSION(1:NZS) :: transp,cap,diffu,hydro , &
1387 thdif,tranf,tav,soilmoism , &
1388 soilicem,soiliqwm,detal , &
1389 fwsat,lwsat,told,smold
1390
1391 REAL :: drip
1392
1393 INTEGER :: nzs1,nzs2,k
1394
1395 !-----------------------------------------------------------------
1396
1397 !-- define constants
1398 ! STBOLT=5.670151E-8
1399 RHOICE=900.
1400 CI=RHOICE*2100.
1401 XLMELT=3.335E+5
1402 cvw=cw
1403
1404 SAT=0.0004
1405 prcpl=prcpms
1406
1407 !--- Initializing local arrays
1408 DO K=1,NZS
1409 TRANSP (K)=0.
1410 soilmoism(k)=0.
1411 soilice (k)=0.
1412 soiliqw (k)=0.
1413 soilicem (k)=0.
1414 soiliqwm (k)=0.
1415 lwsat (k)=0.
1416 fwsat (k)=0.
1417 tav (k)=0.
1418 cap (k)=0.
1419 thdif (k)=0.
1420 diffu (k)=0.
1421 hydro (k)=0.
1422 tranf (k)=0.
1423 detal (k)=0.
1424 told (k)=0.
1425 smold (k)=0.
1426 ENDDO
1427
1428 NZS1=NZS-1
1429 NZS2=NZS-2
1430 dzstop=1./(zsmain(2)-zsmain(1))
1431 RAS=RHO*1.E-3
1432 RIW=rhoice*1.e-3
1433
1434 !--- Computation of volumetric content of ice in soil
1435
1436 DO K=1,NZS
1437 !- main levels
1438 tln=log(tso(k)/273.15)
1439 if(tln.lt.0.) then
1440 soiliqw(k)=(dqm+qmin)*(XLMELT* &
1441 (tso(k)-273.15)/tso(k)/9.81/psis) &
1442 **(-1./bclh)-qmin
1443 soiliqw(k)=max(0.,soiliqw(k))
1444 soiliqw(k)=min(soiliqw(k),soilmois(k))
1445 soilice(k)=(soilmois(k)-soiliqw(k))/RIW
1446
1447 !---- melting and freezing is balanced, soil ice cannot increase
1448 if(keepfr(k).eq.1.) then
1449 soilice(k)=min(soilice(k),smfrkeep(k))
1450 soiliqw(k)=max(0.,soilmois(k)-soilice(k)*riw)
1451 endif
1452
1453 else
1454 soilice(k)=0.
1455 soiliqw(k)=soilmois(k)
1456 endif
1457
1458 ENDDO
1459
1460 DO K=1,NZS1
1461 !- middle of soil layers
1462 tav(k)=0.5*(tso(k)+tso(k+1))
1463 soilmoism(k)=0.5*(soilmois(k)+soilmois(k+1))
1464 tavln=log(tav(k)/273.15)
1465
1466 if(tavln.lt.0.) then
1467 soiliqwm(k)=(dqm+qmin)*(XLMELT* &
1468 (tav(k)-273.15)/tav(k)/9.81/psis) &
1469 **(-1./bclh)-qmin
1470 fwsat(k)=dqm-soiliqwm(k)
1471 lwsat(k)=soiliqwm(k)+qmin
1472 soiliqwm(k)=max(0.,soiliqwm(k))
1473 soiliqwm(k)=min(soiliqwm(k), soilmoism(k))
1474 soilicem(k)=(soilmoism(k)-soiliqwm(k))/riw
1475 !---- melting and freezing is balanced, soil ice cannot increase
1476 if(keepfr(k).eq.1.) then
1477 soilicem(k)=min(soilicem(k), &
1478 0.5*(smfrkeep(k)+smfrkeep(k+1)))
1479 soiliqwm(k)=max(0.,soilmoism(k)-soilicem(k)*riw)
1480 fwsat(k)=dqm-soiliqwm(k)
1481 lwsat(k)=soiliqwm(k)+qmin
1482 endif
1483
1484 else
1485 soilicem(k)=0.
1486 soiliqwm(k)=soilmoism(k)
1487 lwsat(k)=dqm+qmin
1488 fwsat(k)=0.
1489 endif
1490
1491 ENDDO
1492
1493 do k=1,nzs
1494 if(soilice(k).gt.0.) then
1495 smfrkeep(k)=soilice(k)
1496 else
1497 smfrkeep(k)=soilmois(k)/riw
1498 endif
1499 enddo
1500
1501 !******************************************************************
1502 ! SOILPROP computes thermal diffusivity, and diffusional and
1503 ! hydraulic condeuctivities
1504 !******************************************************************
1505 CALL SOILPROP( &
1506 !--- input variables
1507 nzs,fwsat,lwsat,tav,keepfr, &
1508 soilmois,soiliqw,soilice, &
1509 soilmoism,soiliqwm,soilicem, &
1510 !--- soil fixed fields
1511 QWRTZ,rhocs,dqm,qmin,psis,bclh,ksat, &
1512 !--- constants
1513 riw,xlmelt,CP,G0_P,cvw,ci, &
1514 kqwrtz,kice,kwt, &
1515 !--- output variables
1516 thdif,diffu,hydro,cap)
1517
1518 !********************************************************************
1519 !--- CALCULATION OF CANOPY WATER (Smirnova et al., 1996, EQ.16) AND DEW
1520
1521 DRIP=0.
1522 DD1=0.
1523
1524 FQ=QKMS
1525
1526 Q1=-QKMS*RAS*(QVATM - QSG)
1527
1528 DEW=0.
1529 IF(QVATM.GE.QSG)THEN
1530 DEW=FQ*(QVATM-QSG)
1531 ENDIF
1532 IF(DEW.NE.0.)THEN
1533 DD1=CST+DELT*(PRCPMS +DEW*RAS)*vegfrac
1534 ELSE
1535 DD1=CST+ &
1536 DELT*(PRCPMS+RAS*FQ*(QVATM-QSG) &
1537 *(CST/SAT)**CN)*vegfrac
1538 ENDIF
1539
1540 IF(DD1.LT.0.) DD1=0.
1541 if(vegfrac.eq.0.)then
1542 cst=0.
1543 drip=0.
1544 endif
1545 IF (vegfrac.GT.0.) THEN
1546 CST=DD1
1547 IF(CST.GT.SAT) THEN
1548 CST=SAT
1549 DRIP=DD1-SAT
1550 ENDIF
1551 ENDIF
1552
1553 !--- WETCAN is the fraction of vegetated area covered by canopy
1554 !--- water, and DRYCAN is the fraction of vegetated area where
1555 !--- transpiration may take place.
1556
1557 WETCAN=(CST/SAT)**CN
1558 DRYCAN=1.-WETCAN
1559
1560 ! print *,'CST,DRIP',cst,drip
1561
1562 !**************************************************************
1563 ! TRANSF computes transpiration function
1564 !**************************************************************
1565 CALL TRANSF( &
1566 !--- input variables
1567 nzs,nroot,soiliqw, &
1568 !--- soil fixed fields
1569 dqm,qmin,ref,wilt,zshalf, &
1570 !--- output variables
1571 tranf,transum)
1572
1573
1574 !--- Save soil temp and moisture from the beginning of time step
1575 do k=1,nzs
1576 told(k)=tso(k)
1577 smold(k)=soilmois(k)
1578 enddo
1579
1580 !**************************************************************
1581 ! SOILTEMP soilves heat budget and diffusion eqn. in soil
1582 !**************************************************************
1583
1584 CALL SOILTEMP( &
1585 !--- input variables
1586 i,j,iland,isoil, &
1587 delt,ktau,conflx,nzs,nddzs,nroot, &
1588 PRCPMS,RAINF, &
1589 PATM,TABS,QVATM,QCATM,EMISS,RNET, &
1590 QKMS,TKMS,PC,rho,vegfrac, &
1591 thdif,cap,drycan,wetcan, &
1592 transum,dew,mavail, &
1593 !--- soil fixed fields
1594 dqm,qmin,bclh,zsmain,zshalf,DTDZS,tbq, &
1595 !--- constants
1596 xlv,CP,G0_P,cvw,stbolt, &
1597 !--- output variables
1598 tso,soilt,qvg,qsg,qcg)
1599
1600 !************************************************************************
1601
1602 !--- CALCULATION OF DEW USING NEW VALUE OF QSG OR TRANSP IF NO DEW
1603 ETT1=0.
1604 DEW=0.
1605
1606 IF(QVATM.GE.QSG)THEN
1607 DEW=QKMS*(QVATM-QSG)
1608 DO K=1,NZS
1609 TRANSP(K)=0.
1610 ENDDO
1611 ELSE
1612 DO K=1,NROOT
1613 TRANSP(K)=VEGFRAC*RAS*QKMS* &
1614 (QVATM-QSG)* &
1615 PC*TRANF(K)*DRYCAN/ZSHALF(NROOT+1)
1616 IF(TRANSP(K).GT.0.) TRANSP(K)=0.
1617 ETT1=ETT1-TRANSP(K)
1618 ENDDO
1619 DO k=nroot+1,nzs
1620 transp(k)=0.
1621 enddo
1622 ENDIF
1623
1624 !-- Recalculating of volumetric content of frozen water in soil
1625 DO K=1,NZS
1626 !- main levels
1627 tln=log(tso(k)/273.15)
1628 if(tln.lt.0.) then
1629 soiliqw(k)=(dqm+qmin)*(XLMELT* &
1630 (tso(k)-273.15)/tso(k)/9.81/psis) &
1631 **(-1./bclh)-qmin
1632 soiliqw(k)=max(0.,soiliqw(k))
1633 soiliqw(k)=min(soiliqw(k),soilmois(k))
1634 soilice(k)=(soilmois(k)-soiliqw(k))/riw
1635 !---- melting and freezing is balanced, soil ice cannot increase
1636 if(keepfr(k).eq.1.) then
1637 soilice(k)=min(soilice(k),smfrkeep(k))
1638 soiliqw(k)=max(0.,soilmois(k)-soilice(k)*riw)
1639 endif
1640
1641 else
1642 soilice(k)=0.
1643 soiliqw(k)=soilmois(k)
1644 endif
1645 ENDDO
1646
1647 INFMAX=999.
1648 !--- The threshold when the infiltration stops is:
1649 !--- volumetric content of unfrozen pores < 0.12
1650 if((dqm+qmin-riw*soilicem(1)).lt.0.12) &
1651 INFMAX=0.
1652
1653 !*************************************************************************
1654 ! SOILMOIST solves moisture budget (Smirnova et al., 1996, EQ.22,28)
1655 ! and Richards eqn.
1656 !*************************************************************************
1657 CALL SOILMOIST ( &
1658 !-- input
1659 delt,nzs,nddzs,DTDZS,DTDZS2, &
1660 zsmain,zshalf,diffu,hydro, &
1661 QSG,QVG,QCG,QCATM,QVATM,-PRCPMS, &
1662 QKMS,TRANSP,DRIP,DEW,0.,SOILICE,VEGFRAC, &
1663 !-- soil properties
1664 DQM,QMIN,REF,KSAT,RAS,INFMAX, &
1665 !-- output
1666 SOILMOIS,MAVAIL,RUNOFF1, &
1667 RUNOFF2,INFILTRP)
1668
1669 !--- KEEPFR is 1 when the temperature and moisture in soil
1670 !--- are both increasing. In this case soil ice should not
1671 !--- be increasing according to the freezing curve.
1672 !--- Some part of ice is melted, but additional water is
1673 !--- getting frozen. Thus, only structure of frozen soil is
1674 !--- changed, and phase changes are not affecting the heat
1675 !--- transfer. This situation may happen when it rains on the
1676 !--- frozen soil.
1677
1678 do k=1,nzs
1679 if (soilice(k).gt.0.) then
1680 if(tso(k).gt.told(k).and.soilmois(k).gt.smold(k)) then
1681 keepfr(k)=1.
1682 else
1683 keepfr(k)=0.
1684 endif
1685 endif
1686 enddo
1687 !--- THE DIAGNOSTICS OF SURFACE FLUXES
1688
1689 T3 = STBOLT*SOILT*SOILT*SOILT
1690 UPFLUX = T3 *SOILT
1691 XINET = EMISS*(GLW-UPFLUX)
1692 RNET = GSW + XINET
1693 HFT=-TKMS*CP*RHO*(TABS-SOILT)
1694 Q1=-QKMS*RAS*(QVATM - QSG)
1695 EDIR1 =-(1.-vegfrac)*QKMS*RAS* &
1696 (QVATM-QVG)
1697 IF (Q1.LE.0.) THEN
1698 ! --- condensation
1699 EC1=0.
1700 EDIR1=0.
1701 ETT1=0.
1702 EETA=0.
1703 QFX=- XLV*RHO*DEW
1704 ELSE
1705 ! --- evaporation
1706 EC1 = Q1 * WETCAN
1707 CMC2MS=CST/DELT
1708 if(EC1.gt.CMC2MS) cst=0.
1709 EC1=MIN(CMC2MS,EC1)*vegfrac
1710 EETA = (EDIR1 + EC1 + ETT1)*1.E3
1711 ! to convert from kg m-2 s-1 to m s-1: 1/rho water=1.e-3************
1712 QFX= XLV * EETA
1713 ENDIF
1714 EVAPL=QFX/XLV
1715 S=THDIF(1)*CAP(1)*DZSTOP*(TSO(1)-TSO(2))
1716 HFX=HFT
1717 FLTOT=RNET-HFT-QFX-S
1718
1719 222 CONTINUE
1720
1721 1123 FORMAT(I5,8F12.3)
1722 1133 FORMAT(I7,8E12.4)
1723 123 format(i6,f6.2,7f8.1)
1724 122 FORMAT(1X,2I3,6F8.1,F8.3,F8.2)
1725
1726
1727 ! RETURN
1728 ! END
1729 !-------------------------------------------------------------------
1730 END SUBROUTINE SOIL
1731 !-------------------------------------------------------------------
1732
1733
1734 SUBROUTINE SNOWSOIL ( &
1735 !--- input variables
1736 i,j,isoil,delt,ktau,conflx,nzs,nddzs,nroot, &
1737 ILAND,PRCPMS,RAINF,NEWSNOW,snhei,SNWE,SNOWFRAC, &
1738 RHOSN, &
1739 PATM,QVATM,QCATM, &
1740 GLW,GSW,EMISS,RNET,IVGTYP, &
1741 QKMS,TKMS,PC,cst,rho,vegfrac,alb,znt, &
1742 MYJ, &
1743 !--- soil fixed fields
1744 QWRTZ,rhocs,dqm,qmin,ref,wilt,psis,bclh,ksat, &
1745 sat,cn,zsmain,zshalf,DTDZS,DTDZS2,tbq, &
1746 !--- constants
1747 xlv,CP,G0_P,cw,stbolt,TABS, &
1748 KQWRTZ,KICE,KWT, &
1749 !--- output variables
1750 ilnb,snweprint,snheiprint,rsm, &
1751 soilmois,tso,smfrkeep,keepfr, &
1752 dew,soilt,soilt1,tsnav, &
1753 qvg,qsg,qcg,SMELT,SNOH,SNFLX,SNOM, &
1754 edir1,ec1,ett1,eeta,qfx,hfx,s,sublim, &
1755 prcpl,runoff1,runoff2,mavail,soilice, &
1756 soiliqw,infiltrp )
1757
1758 !***************************************************************
1759 ! Energy and moisture budget for snow, heat diffusion eqns.
1760 ! in snow and soil, Richards eqn. for soil covered with snow
1761 !
1762 ! DELT - time step (s)
1763 ! ktau - numver of time step
1764 ! CONFLX - depth of constant flux layer (m)
1765 ! J,I - the location of grid point
1766 ! IME, JME, NZS - dimensions of the domain
1767 ! NROOT - number of levels within the root zone
1768 ! PRCPMS - precipitation rate in m/s
1769 ! NEWSNOW - pcpn in soilid form (m)
1770 ! SNHEI, SNWE - snow height and snow water equivalent (m)
1771 ! RHOSN - snow density (kg/m-3)
1772 ! PATM - pressure (bar)
1773 ! QVATM,QCATM - cloud and water vapor mixing ratio
1774 ! at the first atm. level (kg/kg)
1775 ! GLW, GSW - incoming longwave and absorbed shortwave
1776 ! radiation at the surface (W/m^2)
1777 ! EMISS,RNET - emissivity (0-1) of the ground surface and net
1778 ! radiation at the surface (W/m^2)
1779 ! QKMS - exchange coefficient for water vapor in the
1780 ! surface layer (m/s)
1781 ! TKMS - exchange coefficient for heat in the surface
1782 ! layer (m/s)
1783 ! PC - plant coefficient (resistance) (0-1)
1784 ! RHO - density of atmosphere near surface (kg/m^3)
1785 ! VEGFRAC - greeness fraction (0-1)
1786 ! RHOCS - volumetric heat capacity of dry soil (J/m^3/K)
1787 ! DQM, QMIN - porosity minus residual soil moisture QMIN (m^3/m^3)
1788 ! REF, WILT - field capacity soil moisture and the
1789 ! wilting point (m^3/m^3)
1790 ! PSIS - matrix potential at saturation (m)
1791 ! BCLH - exponent for Clapp-Hornberger parameterization
1792 ! KSAT - saturated hydraulic conductivity (m/s)
1793 ! SAT - maximum value of water intercepted by canopy (m)
1794 ! CN - exponent for calculation of canopy water
1795 ! ZSMAIN - main levels in soil (m)
1796 ! ZSHALF - middle of the soil layers (m)
1797 ! DTDZS,DTDZS2 - dt/(2.*dzshalf*dzmain) and dt/dzshalf in soil
1798 ! TBQ - table to define saturated mixing ration
1799 ! of water vapor for given temperature and pressure
1800 ! ilnb - number of layers in snow
1801 ! rsm - liquid water inside snow pack (m)
1802 ! SOILMOIS,TSO - soil moisture (m^3/m^3) and temperature (K)
1803 ! DEW - dew in (kg/m^2 s)
1804 ! SOILT - skin temperature (K)
1805 ! SOILT1 - snow temperature at 7.5 cm depth (K)
1806 ! TSNAV - average temperature of snow pack (C)
1807 ! QSG,QVG,QCG - saturated mixing ratio, mixing ratio of
1808 ! water vapor and cloud at the ground
1809 ! surface, respectively (kg/kg)
1810 ! EDIR1, EC1, ETT1, EETA - direct evaporation, evaporation of
1811 ! canopy water, transpiration (kg m-2 s-1) and total
1812 ! evaporation in (m s-1).
1813 ! QFX, HFX - latent and sensible heat fluxes (W/m^2)
1814 ! S - soil heat flux in the top layer (W/m^2)
1815 ! SUBLIM - snow sublimation (kg/m^2/s)
1816 ! RUNOFF1 - surface runoff (m/s)
1817 ! RUNOFF2 - underground runoff (m)
1818 ! MAVAIL - moisture availability in the top soil layer (0-1)
1819 ! SOILICE - content of soil ice in soil layers (m^3/m^3)
1820 ! SOILIQW - lliquid water in soil layers (m^3/m^3)
1821 ! INFILTRP - infiltration flux from the top of soil domain (m/s)
1822 ! XINET - net long-wave radiation (W/m^2)
1823 !
1824 !*******************************************************************
1825
1826 IMPLICIT NONE
1827 !-------------------------------------------------------------------
1828 !--- input variables
1829
1830 INTEGER, INTENT(IN ) :: nroot,ktau,nzs , &
1831 nddzs !nddzs=2*(nzs-2)
1832 INTEGER, INTENT(IN ) :: i,j,isoil
1833
1834 REAL, INTENT(IN ) :: DELT,CONFLX,PRCPMS , &
1835 RAINF,NEWSNOW
1836
1837 LOGICAL, INTENT(IN ) :: myj
1838
1839 !--- 3-D Atmospheric variables
1840 REAL, &
1841 INTENT(IN ) :: PATM, &
1842 QVATM, &
1843 QCATM
1844 !--- 2-D variables
1845 REAL , &
1846 INTENT(IN ) :: GLW, &
1847 GSW, &
1848 RHO, &
1849 PC, &
1850 VEGFRAC, &
1851 QKMS, &
1852 TKMS
1853
1854 INTEGER, INTENT(IN ) :: IVGTYP
1855 !--- soil properties
1856 REAL , &
1857 INTENT(IN ) :: RHOCS, &
1858 BCLH, &
1859 DQM, &
1860 KSAT, &
1861 PSIS, &
1862 QMIN, &
1863 QWRTZ, &
1864 REF, &
1865 SAT, &
1866 WILT
1867
1868 REAL, INTENT(IN ) :: CN, &
1869 CW, &
1870 XLV, &
1871 G0_P, &
1872 KQWRTZ, &
1873 KICE, &
1874 KWT
1875
1876
1877 REAL, DIMENSION(1:NZS), INTENT(IN) :: ZSMAIN, &
1878 ZSHALF, &
1879 DTDZS2
1880
1881 REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS
1882
1883 REAL, DIMENSION(1:4001), INTENT(IN) :: TBQ
1884
1885
1886 !--- input/output variables
1887 !-------- 3-d soil moisture and temperature
1888 REAL, DIMENSION( 1:nzs ) , &
1889 INTENT(INOUT) :: TSO, &
1890 SOILMOIS, &
1891 SMFRKEEP
1892
1893 REAL, DIMENSION( 1:nzs ) , &
1894 INTENT(INOUT) :: KEEPFR
1895
1896
1897 INTEGER, INTENT(INOUT) :: ILAND
1898
1899
1900 !-------- 2-d variables
1901 REAL , &
1902 INTENT(INOUT) :: DEW, &
1903 CST, &
1904 EDIR1, &
1905 EC1, &
1906 ETT1, &
1907 EETA, &
1908 RHOSN, &
1909 SUBLIM, &
1910 PRCPL, &
1911 ALB, &
1912 EMISS, &
1913 ZNT, &
1914 MAVAIL, &
1915 QVG, &
1916 QSG, &
1917 QCG, &
1918 QFX, &
1919 HFX, &
1920 S, &
1921 RUNOFF1, &
1922 RUNOFF2, &
1923 SNWE, &
1924 SNHEI, &
1925 SMELT, &
1926 SNOM, &
1927 SNOH, &
1928 SNFLX, &
1929 SOILT, &
1930 SOILT1, &
1931 SNOWFRAC, &
1932 TSNAV
1933
1934 INTEGER, INTENT(INOUT) :: ILNB
1935
1936 !-------- 1-d variables
1937 REAL, DIMENSION(1:NZS), INTENT(OUT) :: SOILICE, &
1938 SOILIQW
1939
1940 REAL, INTENT(OUT) :: RSM, &
1941 SNWEPRINT, &
1942 SNHEIPRINT
1943 !--- Local variables
1944
1945
1946 INTEGER :: nzs1,nzs2,k
1947
1948 REAL :: INFILTRP, RHONEWSN,TRANSUM , &
1949 SNTH, NEWSN , &
1950 TABS, T3, UPFLUX, XINET , &
1951 BETA, SNWEPR,EPDT,PP
1952 REAL :: CP,G0,LV,xlvm,STBOLT,xlmelt,dzstop , &
1953 can,epot,fac,fltot,ft,fq,hft , &
1954 q1,ras,rhoice,sph , &
1955 trans,zn,ci,cvw,tln,tavln,pi , &
1956 DD1,CMC2MS,DRYCAN,WETCAN , &
1957 INFMAX,RIW,DELTSN,H,UMVEG
1958
1959 REAL, DIMENSION(1:NZS) :: transp,cap,diffu,hydro , &
1960 thdif,tranf,tav,soilmoism , &
1961 soilicem,soiliqwm,detal , &
1962 fwsat,lwsat,told,smold
1963 REAL :: drip
1964
1965 REAL :: RNET
1966
1967 !-----------------------------------------------------------------
1968
1969 cvw=cw
1970 XLMELT=3.335E+5
1971 !-- the next line calculates heat of sublimation of water vapor
1972 XLVm=XLV+XLMELT
1973 ! STBOLT=5.670151E-8
1974
1975 !--- SNOW flag -- 99
1976 ILAND=99
1977
1978 !--- DELTSN - is the threshold for splitting the snow layer into 2 layers.
1979 !--- With snow density 400 kg/m^3, this threshold is equal to 7.5 cm,
1980 !--- equivalent to 0.03 m SNWE. For other snow densities the threshold is
1981 !--- computed using SNWE=0.03 m and current snow density.
1982 !--- SNTH - the threshold below which the snow layer is combined with
1983 !--- the top soil layer. SNTH is computed using snwe=0.016 m, and
1984 !--- equals 4 cm for snow density 400 kg/m^3.
1985
1986 DELTSN=0.0301*1.e3/rhosn
1987 snth=0.01601*1.e3/rhosn
1988
1989 RHOICE=900.
1990 CI=RHOICE*2100.
1991 RAS=RHO*1.E-3
1992 RIW=rhoice*1.e-3
1993 MAVAIL=1.
1994 RSM=0.
1995
1996 DO K=1,NZS
1997 TRANSP (K)=0.
1998 soilmoism (k)=0.
1999 soiliqwm (k)=0.
2000 soilice (k)=0.
2001 soilicem (k)=0.
2002 lwsat (k)=0.
2003 fwsat (k)=0.
2004 tav (k)=0.
2005 cap (k)=0.
2006 diffu (k)=0.
2007 hydro (k)=0.
2008 thdif (k)=0.
2009 tranf (k)=0.
2010 detal (k)=0.
2011 told (k)=0.
2012 smold (k)=0.
2013 ENDDO
2014
2015 snweprint=0.
2016 snheiprint=0.
2017 prcpl=prcpms
2018
2019 !*** DELTSN is the depth of the top layer of snow where
2020 !*** there is a temperature gradient, the rest of the snow layer
2021 !*** is considered to have constant temperature
2022
2023
2024 NZS1=NZS-1
2025 NZS2=NZS-2
2026 DZSTOP=1./(zsmain(2)-zsmain(1))
2027
2028 !----- THE CALCULATION OF THERMAL DIFFUSIVITY, DIFFUSIONAL AND ---
2029 !----- HYDRAULIC CONDUCTIVITY (SMIRNOVA ET AL. 1996, EQ.2,5,6) ---
2030 !tgs - the following loop is added to define the amount of frozen
2031 !tgs - water in soil if there is any
2032 DO K=1,NZS
2033
2034 tln=log(tso(k)/273.15)
2035 if(tln.lt.0.) then
2036 soiliqw(k)=(dqm+qmin)*(XLMELT* &
2037 (tso(k)-273.15)/tso(k)/9.81/psis) &
2038 **(-1./bclh)-qmin
2039 soiliqw(k)=max(0.,soiliqw(k))
2040 soiliqw(k)=min(soiliqw(k),soilmois(k))
2041 soilice(k)=(soilmois(k)-soiliqw(k))/riw
2042
2043 !---- melting and freezing is balanced, soil ice cannot increase
2044 if(keepfr(k).eq.1.) then
2045 soilice(k)=min(soilice(k),smfrkeep(k))
2046 soiliqw(k)=max(0.,soilmois(k)-soilice(k)*rhoice*1.e-3)
2047 endif
2048
2049 else
2050 soilice(k)=0.
2051 soiliqw(k)=soilmois(k)
2052 endif
2053
2054 ENDDO
2055
2056 DO K=1,NZS1
2057
2058 tav(k)=0.5*(tso(k)+tso(k+1))
2059 soilmoism(k)=0.5*(soilmois(k)+soilmois(k+1))
2060 tavln=log(tav(k)/273.15)
2061
2062 if(tavln.lt.0.) then
2063 soiliqwm(k)=(dqm+qmin)*(XLMELT* &
2064 (tav(k)-273.15)/tav(k)/9.81/psis) &
2065 **(-1./bclh)-qmin
2066 fwsat(k)=dqm-soiliqwm(k)
2067 lwsat(k)=soiliqwm(k)+qmin
2068 soiliqwm(k)=max(0.,soiliqwm(k))
2069 soiliqwm(k)=min(soiliqwm(k), soilmoism(k))
2070 soilicem(k)=(soilmoism(k)-soiliqwm(k))/riw
2071 !---- melting and freezing is balanced, soil ice cannot increase
2072 if(keepfr(k).eq.1.) then
2073 soilicem(k)=min(soilicem(k), &
2074 0.5*(smfrkeep(k)+smfrkeep(k+1)))
2075 soiliqwm(k)=max(0.,soilmoism(k)-soilicem(k)*riw)
2076 fwsat(k)=dqm-soiliqwm(k)
2077 lwsat(k)=soiliqwm(k)+qmin
2078 endif
2079
2080 else
2081 soilicem(k)=0.
2082 soiliqwm(k)=soilmoism(k)
2083 lwsat(k)=dqm+qmin
2084 fwsat(k)=0.
2085
2086 endif
2087 ENDDO
2088
2089 do k=1,nzs
2090 if(soilice(k).gt.0.) then
2091 smfrkeep(k)=soilice(k)
2092 else
2093 smfrkeep(k)=soilmois(k)/riw
2094 endif
2095 enddo
2096
2097
2098 ! print *,'etaf,etal,etamf,etaml,lwsat,fwsat',
2099 ! 1 soilice,soiliqw,soilicem,soiliqwm,lwsat,fwsat
2100
2101 !******************************************************************
2102 ! SOILPROP computes thermal diffusivity, and diffusional and
2103 ! hydraulic condeuctivities
2104 !******************************************************************
2105 CALL SOILPROP( &
2106 !--- input variables
2107 nzs,fwsat,lwsat,tav,keepfr, &
2108 soilmois,soiliqw,soilice, &
2109 soilmoism,soiliqwm,soilicem, &
2110 !--- soil fixed fields
2111 QWRTZ,rhocs,dqm,qmin,psis,bclh,ksat, &
2112 !--- constants
2113 riw,xlmelt,CP,G0_P,cvw,ci, &
2114 kqwrtz,kice,kwt, &
2115 !--- output variables
2116 thdif,diffu,hydro,cap)
2117
2118 !********************************************************************
2119 !--- CALCULATION OF CANOPY WATER (Smirnova et al., 1996, EQ.16) AND DEW
2120
2121 DRIP=0.
2122 SMELT=0.
2123 DD1=0.
2124 H=1.
2125
2126 FQ=QKMS
2127
2128
2129 !--- If vegfrac.ne.0. then part of falling snow can be
2130 !--- intercepted by the canopy.
2131
2132 DEW=0.
2133 UMVEG=1.-vegfrac
2134 EPOT = -FQ*(QVATM-QSG)
2135
2136 IF(vegfrac.EQ.0.) then
2137 cst=0.
2138 drip=0.
2139 ELSE
2140 IF(EPOT.GE.0.) THEN
2141 ! Evaporation
2142 DD1=CST+(NEWSNOW*RHOSN*1.E-3 &
2143 !-- need to think more if we want this change....
2144 -DELT*(RAS*EPOT &
2145 ! -DELT*(-PRCPMS+RAS*EPOT &
2146 *(CST/SAT)**CN)) *vegfrac
2147 ELSE
2148 ! Sublimation
2149 DEW = - EPOT
2150 DD1=CST+(NEWSNOW*RHOSN*1.E-3+delt*( &
2151 ! DD1=CST+(NEWSNOW*RHOSN*1.E-3+delt*(PRCPMS &
2152 +DEW*RAS)) *vegfrac
2153 ENDIF
2154
2155 IF(DD1.LT.0.) DD1=0.
2156 IF (vegfrac.GT.0.) THEN
2157 CST=DD1
2158 IF(CST.GT.SAT*vegfrac) THEN
2159 CST=SAT*vegfrac
2160 DRIP=DD1-SAT*vegfrac
2161 ENDIF
2162 ENDIF
2163
2164
2165 !--- In SFCTMP NEWSNOW is added to SNHEI as if there is no vegetation
2166 !--- With vegetation part of NEWSNOW can be intercepted by canopy until
2167 !--- the saturation is reached. After the canopy saturation is reached
2168 !--- DRIP in the solid form will be added to SNOW cover.
2169
2170 !!! 4 Nov 07 SNWE=(SNHEI-vegfrac*NEWSNOW)*RHOSN*1.E-3 &
2171 !!! 4 Nov 07 + DRIP &
2172 ! - 10% of liquid precip could be added to snow water
2173 ! - this is based on SnowMIP2.
2174 ! - something more intelligent should be done to liquid water
2175 SNWE = SNWE &
2176 +0.10*prcpms*delt
2177
2178
2179 ENDIF
2180
2181 DRIP=0.
2182 SNHEI=SNWE*1.e3/RHOSN
2183 SNWEPR=SNWE
2184
2185 ! check if all snow can evaporate during DT
2186 BETA=1.
2187 EPDT = EPOT * RAS *DELT*UMVEG
2188 IF(SNWEPR.LE.EPDT) THEN
2189 BETA=SNWEPR/max(1.e-8,EPDT)
2190 SNWE=0.
2191 SNHEI=0.
2192 ENDIF
2193
2194 WETCAN=(CST/SAT)**CN
2195 DRYCAN=1.-WETCAN
2196
2197 !**************************************************************
2198 ! TRANSF computes transpiration function
2199 !**************************************************************
2200 CALL TRANSF( &
2201 !--- input variables
2202 nzs,nroot,soiliqw, &
2203 !--- soil fixed fields
2204 dqm,qmin,ref,wilt,zshalf, &
2205 !--- output variables
2206 tranf,transum)
2207
2208 !--- Save soil temp and moisture from the beginning of time step
2209 do k=1,nzs
2210 told(k)=tso(k)
2211 smold(k)=soilmois(k)
2212 enddo
2213
2214 !**************************************************************
2215 ! SOILTEMP soilves heat budget and diffusion eqn. in soil
2216 !**************************************************************
2217
2218 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
2219 print *, 'TSO before calling SNOWTEMP: ', tso
2220 ENDIF
2221 CALL SNOWTEMP( &
2222 !--- input variables
2223 i,j,iland,isoil, &
2224 delt,ktau,conflx,nzs,nddzs,nroot, &
2225 snwe,snwepr,snhei,newsnow,snowfrac, &
2226 beta,deltsn,snth,rhosn, &
2227 PRCPMS,RAINF, &
2228 PATM,TABS,QVATM,QCATM, &
2229 GLW,GSW,EMISS,RNET, &
2230 QKMS,TKMS,PC,rho,vegfrac, &
2231 thdif,cap,drycan,wetcan,cst, &
2232 tranf,transum,dew,mavail, &
2233 !--- soil fixed fields
2234 dqm,qmin,psis,bclh, &
2235 zsmain,zshalf,DTDZS,tbq, &
2236 !--- constants
2237 xlvm,CP,G0_P,cvw,stbolt, &
2238 !--- output variables
2239 snweprint,snheiprint,rsm, &
2240 tso,soilt,soilt1,tsnav,qvg,qsg,qcg, &
2241 smelt,snoh,snflx,ilnb)
2242
2243 !************************************************************************
2244 !--- RECALCULATION OF DEW USING NEW VALUE OF QSG OR TRANSP IF NO DEW
2245 DEW=0.
2246 ETT1=0.
2247 PP=PATM*1.E3
2248 QSG= QSN(SOILT,TBQ)/PP
2249 EPOT = -FQ*(QVATM-QSG)
2250 IF(EPOT.GE.0.) THEN
2251 ! Evaporation
2252 DO K=1,NROOT
2253 TRANSP(K)=vegfrac*RAS*FQ*(QVATM-QSG) &
2254 *PC*tranf(K)*DRYCAN/zshalf(NROOT+1)
2255 IF(TRANSP(K).GT.0.) TRANSP(K)=0.
2256 ETT1=ETT1-TRANSP(K)
2257 ENDDO
2258 DO k=nroot+1,nzs
2259 transp(k)=0.
2260 enddo
2261
2262 ELSE
2263 ! Sublimation
2264 DEW=-EPOT
2265 DO K=1,NZS
2266 TRANSP(K)=0.
2267 ENDDO
2268 ETT1=0.
2269 ENDIF
2270
2271 !-- recalculating of frozen water in soil
2272 DO K=1,NZS
2273 tln=log(tso(k)/273.15)
2274 if(tln.lt.0.) then
2275 soiliqw(k)=(dqm+qmin)*(XLMELT* &
2276 (tso(k)-273.15)/tso(k)/9.81/psis) &
2277 **(-1./bclh)-qmin
2278 soiliqw(k)=max(0.,soiliqw(k))
2279 soiliqw(k)=min(soiliqw(k),soilmois(k))
2280 soilice(k)=(soilmois(k)-soiliqw(k))/riw
2281 !---- melting and freezing is balanced, soil ice cannot increase
2282 if(keepfr(k).eq.1.) then
2283 soilice(k)=min(soilice(k),smfrkeep(k))
2284 soiliqw(k)=max(0.,soilmois(k)-soilice(k)*riw)
2285 endif
2286
2287 else
2288 soilice(k)=0.
2289 soiliqw(k)=soilmois(k)
2290 endif
2291 ENDDO
2292
2293 INFMAX=999.
2294 !--- The threshold when the infiltration stops is:
2295 !--- volumetric content of unfrozen pores < 0.12
2296 soilicem(1)=0.5*(soilice(1)+soilice(2))
2297 if((dqm+qmin-riw*soilicem(1)).lt.0.12) &
2298 INFMAX=0.
2299
2300 !*************************************************************************
2301 !--- TQCAN FOR SOLUTION OF MOISTURE BALANCE (Smirnova et al. 1996, EQ.22,28)
2302 ! AND TSO,ETA PROFILES
2303 !*************************************************************************
2304 CALL SOILMOIST ( &
2305 !-- input
2306 delt,nzs,nddzs,DTDZS,DTDZS2, &
2307 zsmain,zshalf,diffu,hydro, &
2308 ! 4 Nov 07 QSG,QVG,QCG,QCATM,QVATM,-0.9*PRCPMS/(1.-vegfrac), &
2309 QSG,QVG,QCG,QCATM,QVATM,-0.9*PRCPMS, &
2310 ! QSG,QVG,QCG,QCATM,QVATM,-PRCPMS, &
2311 0.,TRANSP,0., &
2312 0.,SMELT,soilice,vegfrac, &
2313 !-- soil properties
2314 DQM,QMIN,REF,KSAT,RAS,INFMAX, &
2315 !-- output
2316 soilmois,MAVAIL,RUNOFF1, &
2317 RUNOFF2,infiltrp)
2318
2319 ! 4 Nov 07 - update CST for snow melt
2320 if(snwe.ne.0.) then
2321 CST=(1.-min(1.,smelt/snwe))*CST
2322 else
2323 CST=0.
2324 endif
2325
2326 !-- Restore land-use parameters if snow is less than threshold
2327 IF(SNHEI.EQ.0.) then
2328 tsnav=soilt-273.15
2329 CALL SNOWFREE(ivgtyp,myj,emiss, &
2330 znt,iland)
2331 smelt=smelt+snwe/delt
2332 rsm=0.
2333 ! snwe=0.
2334 ENDIF
2335
2336 SNOM=SNOM+SMELT*DELT
2337
2338 !--- KEEPFR is 1 when the temperature and moisture in soil
2339 !--- are both increasing. In this case soil ice should not
2340 !--- be increasing according to the freezing curve.
2341 !--- Some part of ice is melted, but additional water is
2342 !--- getting frozen. Thus, only structure of frozen soil is
2343 !--- changed, and phase changes are not affecting the heat
2344 !--- transfer. This situation may happen when it rains on the
2345 !--- frozen soil.
2346
2347 do k=1,nzs
2348 if (soilice(k).gt.0.) then
2349 if(tso(k).gt.told(k).and.soilmois(k).gt.smold(k)) then
2350 keepfr(k)=1.
2351 else
2352 keepfr(k)=0.
2353 endif
2354 endif
2355 enddo
2356 !--- THE DIAGNOSTICS OF SURFACE FLUXES
2357
2358 T3 = STBOLT*SOILT*SOILT*SOILT
2359 UPFLUX = T3 *SOILT
2360 XINET = EMISS*(GLW-UPFLUX)
2361 RNET = GSW + XINET
2362 HFT=- TKMS*CP*RHO*(TABS-SOILT)
2363 Q1 = - FQ*RAS* (QVATM - QSG)
2364 EDIR1 = Q1*UMVEG *BETA
2365
2366 IF (Q1.LT.0.) THEN
2367 ! --- condensation
2368 EC1=0.
2369 EDIR1=0.
2370 ETT1=0.
2371 EETA=0.
2372 DEW=FQ*(QVATM-QSG)
2373 QFX= -XLVm*RHO*DEW
2374 sublim=QFX/XLVm
2375 ELSE
2376 ! --- evaporation
2377 EC1 = Q1 * WETCAN
2378 CMC2MS=CST/DELT
2379 if(EC1.gt.CMC2MS) cst=0.
2380 EC1=MIN(CMC2MS,EC1)*vegfrac
2381 EETA = (EDIR1 + EC1 + ETT1)*1.E3
2382 ! to convert from kg m-2 s-1 to m s-1: 1/rho water=1.e-3************
2383 QFX= XLVm * EETA
2384 sublim=(EDIR1 + EC1)*1.E3
2385 ENDIF
2386 s=THDIF(1)*CAP(1)*dzstop*(tso(1)-tso(2))
2387 HFX=HFT
2388 FLTOT=RNET-HFT-QFX-S
2389
2390 222 CONTINUE
2391
2392 1123 FORMAT(I5,8F12.3)
2393 1133 FORMAT(I7,8E12.4)
2394 123 format(i6,f6.2,7f8.1)
2395 122 FORMAT(1X,2I3,6F8.1,F8.3,F8.2)
2396
2397
2398 ! RETURN
2399 ! END
2400 !-------------------------------------------------------------------
2401 END SUBROUTINE SNOWSOIL
2402 !-------------------------------------------------------------------
2403
2404
2405 SUBROUTINE SOILTEMP( &
2406 !--- input variables
2407 i,j,iland,isoil, &
2408 delt,ktau,conflx,nzs,nddzs,nroot, &
2409 PRCPMS,RAINF,PATM,TABS,QVATM,QCATM, &
2410 EMISS,RNET, &
2411 QKMS,TKMS,PC,RHO,VEGFRAC, &
2412 THDIF,CAP,DRYCAN,WETCAN, &
2413 TRANSUM,DEW,MAVAIL, &
2414 !--- soil fixed fields
2415 DQM,QMIN,BCLH, &
2416 ZSMAIN,ZSHALF,DTDZS,TBQ, &
2417 !--- constants
2418 XLV,CP,G0_P,CVW,STBOLT, &
2419 !--- output variables
2420 TSO,SOILT,QVG,QSG,QCG)
2421
2422 !*************************************************************
2423 ! Energy budget equation and heat diffusion eqn are
2424 ! solved here and
2425 !
2426 ! DELT - time step (s)
2427 ! ktau - numver of time step
2428 ! CONFLX - depth of constant flux layer (m)
2429 ! IME, JME, KME, NZS - dimensions of the domain
2430 ! NROOT - number of levels within the root zone
2431 ! PRCPMS - precipitation rate in m/s
2432 ! COTSO, RHTSO - coefficients for implicit solution of
2433 ! heat diffusion equation
2434 ! THDIF - thermal diffusivity (m^2/s)
2435 ! QSG,QVG,QCG - saturated mixing ratio, mixing ratio of
2436 ! water vapor and cloud at the ground
2437 ! surface, respectively (kg/kg)
2438 ! PATM - pressure [bar]
2439 ! QC3D,QV3D - cloud and water vapor mixing ratio
2440 ! at the first atm. level (kg/kg)
2441 ! EMISS,RNET - emissivity (0-1) of the ground surface and net
2442 ! radiation at the surface (W/m^2)
2443 ! QKMS - exchange coefficient for water vapor in the
2444 ! surface layer (m/s)
2445 ! TKMS - exchange coefficient for heat in the surface
2446 ! layer (m/s)
2447 ! PC - plant coefficient (resistance)
2448 ! RHO - density of atmosphere near surface (kg/m^3)
2449 ! VEGFRAC - greeness fraction (0-1)
2450 ! CAP - volumetric heat capacity (J/m^3/K)
2451 ! DRYCAN - dry fraction of vegetated area where
2452 ! transpiration may take place (0-1)
2453 ! WETCAN - fraction of vegetated area covered by canopy
2454 ! water (0-1)
2455 ! TRANSUM - transpiration function integrated over the
2456 ! rooting zone (m)
2457 ! DEW - dew in kg/m^2s
2458 ! MAVAIL - fraction of maximum soil moisture in the top
2459 ! layer (0-1)
2460 ! ZSMAIN - main levels in soil (m)
2461 ! ZSHALF - middle of the soil layers (m)
2462 ! DTDZS - dt/(2.*dzshalf*dzmain)
2463 ! TBQ - table to define saturated mixing ration
2464 ! of water vapor for given temperature and pressure
2465 ! TSO - soil temperature (K)
2466 ! SOILT - skin temperature (K)
2467 !
2468 !****************************************************************
2469
2470 IMPLICIT NONE
2471 !-----------------------------------------------------------------
2472
2473 !--- input variables
2474
2475 INTEGER, INTENT(IN ) :: nroot,ktau,nzs , &
2476 nddzs !nddzs=2*(nzs-2)
2477 INTEGER, INTENT(IN ) :: i,j,iland,isoil
2478 REAL, INTENT(IN ) :: DELT,CONFLX,PRCPMS, RAINF
2479 REAL, INTENT(INOUT) :: DRYCAN,WETCAN,TRANSUM
2480 !--- 3-D Atmospheric variables
2481 REAL, &
2482 INTENT(IN ) :: PATM, &
2483 QVATM, &
2484 QCATM
2485 !--- 2-D variables
2486 REAL , &
2487 INTENT(IN ) :: &
2488 EMISS, &
2489 RHO, &
2490 RNET, &
2491 PC, &
2492 VEGFRAC, &
2493 DEW, &
2494 QKMS, &
2495 TKMS
2496
2497 !--- soil properties
2498 REAL , &
2499 INTENT(IN ) :: &
2500 BCLH, &
2501 DQM, &
2502 QMIN
2503
2504 REAL, INTENT(IN ) :: CP, &
2505 CVW, &
2506 XLV, &
2507 STBOLT, &
2508 TABS, &
2509 G0_P
2510
2511
2512 REAL, DIMENSION(1:NZS), INTENT(IN) :: ZSMAIN, &
2513 ZSHALF, &
2514 THDIF, &
2515 CAP
2516
2517 REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS
2518
2519 REAL, DIMENSION(1:4001), INTENT(IN) :: TBQ
2520
2521
2522 !--- input/output variables
2523 !-------- 3-d soil moisture and temperature
2524 REAL, DIMENSION( 1:nzs ) , &
2525 INTENT(INOUT) :: TSO
2526
2527 !-------- 2-d variables
2528 REAL , &
2529 INTENT(INOUT) :: &
2530 MAVAIL, &
2531 QVG, &
2532 QSG, &
2533 QCG, &
2534 SOILT
2535
2536
2537 !--- Local variables
2538
2539 REAL :: x,x1,x2,x4,dzstop,can,ft,sph , &
2540 tn,trans,umveg,denom
2541
2542 REAL :: FKT,D1,D2,D9,D10,DID,R211,R21,R22,R6,R7,D11 , &
2543 PI,H,FKQ,R210,AA,BB,PP,Q1,QS1,TS1,TQ2,TX2 , &
2544 TDENOM
2545
2546 REAL :: C,CC,AA1,RHCS,H1
2547
2548 REAL, DIMENSION(1:NZS) :: cotso,rhtso
2549
2550 INTEGER :: nzs1,nzs2,k,k1,kn,kk
2551
2552 !-----------------------------------------------------------------
2553
2554
2555 NZS1=NZS-1
2556 NZS2=NZS-2
2557 dzstop=1./(ZSMAIN(2)-ZSMAIN(1))
2558
2559 do k=1,nzs
2560 cotso(k)=0.
2561 rhtso(k)=0.
2562 enddo
2563 !******************************************************************************
2564 ! COEFFICIENTS FOR THOMAS ALGORITHM FOR TSO
2565 !******************************************************************************
2566 ! did=2.*(ZSMAIN(nzs)-ZSHALF(nzs))
2567 ! h1=DTDZS(8)*THDIF(nzs-1)*(ZSHALF(nzs)-ZSHALF(nzs-1))/did
2568 ! cotso(1)=h1/(1.+h1)
2569 ! rhtso(1)=(tso(nzs)+h1*(tso(nzs-1)-tso(nzs)))/
2570 ! 1 (1.+h1)
2571 cotso(1)=0.
2572 rhtso(1)=TSO(NZS)
2573 DO 33 K=1,NZS2
2574 KN=NZS-K
2575 K1=2*KN-3
2576 X1=DTDZS(K1)*THDIF(KN-1)
2577 X2=DTDZS(K1+1)*THDIF(KN)
2578 FT=TSO(KN)+X1*(TSO(KN-1)-TSO(KN)) &
2579 -X2*(TSO(KN)-TSO(KN+1))
2580 DENOM=1.+X1+X2-X2*cotso(K)
2581 cotso(K+1)=X1/DENOM
2582 rhtso(K+1)=(FT+X2*rhtso(K))/DENOM
2583 33 CONTINUE
2584
2585 !************************************************************************
2586 !--- THE HEAT BALANCE EQUATION (Smirnova et al., 1996, EQ. 21,26)
2587
2588 RHCS=CAP(1)
2589 H=MAVAIL
2590 IF(DEW.NE.0.)THEN
2591 DRYCAN=0.
2592 WETCAN=1.
2593 ENDIf
2594 TRANS=PC*TRANSUM*DRYCAN/ZSHALF(NROOT+1)
2595 CAN=WETCAN+TRANS
2596 UMVEG=1.-VEGFRAC
2597 FKT=TKMS
2598 D1=cotso(NZS1)
2599 D2=rhtso(NZS1)
2600 TN=SOILT
2601 D9=THDIF(1)*RHCS*dzstop
2602 D10=TKMS*CP*RHO
2603 R211=.5*CONFLX/DELT
2604 R21=R211*CP*RHO
2605 R22=.5/(THDIF(1)*DELT*dzstop**2)
2606 R6=EMISS *STBOLT*.5*TN**4
2607 R7=R6/TN
2608 D11=RNET+R6
2609 TDENOM=D9*(1.-D1+R22)+D10+R21+R7 &
2610 +RAINF*CVW*PRCPMS
2611 FKQ=QKMS*RHO
2612 R210=R211*RHO
2613 C=VEGFRAC*FKQ*CAN
2614 CC=C*XLV/TDENOM
2615 AA=XLV*(FKQ*UMVEG+R210)/TDENOM
2616 BB=(D10*TABS+R21*TN+XLV*(QVATM* &
2617 (FKQ*UMVEG+C) &
2618 +R210*QVG)+D11+D9*(D2+R22*TN) &
2619 +RAINF*CVW*PRCPMS*max(273.15,TABS) &
2620 )/TDENOM
2621 AA1=AA+CC
2622 PP=PATM*1.E3
2623 AA1=AA1/PP
2624 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
2625 PRINT *,' VILKA-1'
2626 print *,'D10,TABS,R21,TN,QVATM,FKQ,UMVEG,VEGFRAC,CAN', &
2627 D10,TABS,R21,TN,QVATM,FKQ,UMVEG,VEGFRAC,CAN
2628 print *,'RNET, EMISS, STBOLT, SOILT',RNET, EMISS, STBOLT, SOILT
2629 print *,'R210,QVG,D11,D9,D2,R22,RAINF,CVW,PRCPMS,TDENOM', &
2630 R210,QVG,D11,D9,D2,R22,RAINF,CVW,PRCPMS,TDENOM
2631 print *,'tn,aa1,bb,pp,umveg,fkq,r210,vegfrac', &
2632 tn,aa1,bb,pp,umveg,fkq,r210,vegfrac
2633 ENDIF
2634 CALL VILKA(TN,AA1,BB,PP,QS1,TS1,TBQ,KTAU,i,j,iland,isoil)
2635 TQ2=QVATM+QCATM
2636 TX2=TQ2*(1.-H)
2637 Q1=TX2+H*QS1
2638 IF(Q1.LT.QS1) GOTO 100
2639 !--- if no saturation - goto 100
2640 !--- if saturation - goto 90
2641 90 QVG=QS1
2642 QSG=QS1
2643 TSO(1)=TS1
2644 QCG=Q1-QS1
2645 GOTO 200
2646 100 BB=BB-AA*TX2
2647 AA=(AA*H+CC)/PP
2648 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
2649 PRINT *,' VILKA-2'
2650 print *,'D10,TABS,R21,TN,QVATM,FKQ,UMVEG,VEGFRAC,CAN', &
2651 D10,TABS,R21,TN,QVATM,FKQ,UMVEG,VEGFRAC,CAN
2652 print *,'R210,QVG,D11,D9,D2,R22,RAINF,CVW,PRCPMS,TDENOM', &
2653 R210,QVG,D11,D9,D2,R22,RAINF,CVW,PRCPMS,TDENOM
2654
2655 print *,'tn,aa1,bb,pp,umveg,fkq,r210,vegfrac', &
2656 tn,aa1,bb,pp,umveg,fkq,r210,vegfrac
2657 ENDIF
2658
2659 CALL VILKA(TN,AA,BB,PP,QS1,TS1,TBQ,KTAU,i,j,iland,isoil)
2660 Q1=TX2+H*QS1
2661 IF(Q1.GT.QS1) GOTO 90
2662 QSG=QS1
2663 QVG=Q1
2664 TSO(1)=TS1
2665 QCG=0.
2666 200 CONTINUE
2667
2668 !--- SOILT - skin temperature
2669 SOILT=TS1
2670
2671 !---- Final solution for soil temperature - TSO
2672 DO K=2,NZS
2673 KK=NZS-K+1
2674 TSO(K)=rhtso(KK)+cotso(KK)*TSO(K-1)
2675 END DO
2676
2677 ! return
2678 ! end
2679 !--------------------------------------------------------------------
2680 END SUBROUTINE SOILTEMP
2681 !--------------------------------------------------------------------
2682
2683
2684 SUBROUTINE SNOWTEMP( &
2685 !--- input variables
2686 i,j,iland,isoil, &
2687 delt,ktau,conflx,nzs,nddzs,nroot, &
2688 snwe,snwepr,snhei,newsnow,snowfrac, &
2689 beta,deltsn,snth,rhosn, &
2690 PRCPMS,RAINF, &
2691 PATM,TABS,QVATM,QCATM, &
2692 GLW,GSW,EMISS,RNET, &
2693 QKMS,TKMS,PC,RHO,VEGFRAC, &
2694 THDIF,CAP,DRYCAN,WETCAN,CST, &
2695 TRANF,TRANSUM,DEW,MAVAIL, &
2696 !--- soil fixed fields
2697 DQM,QMIN,PSIS,BCLH, &
2698 ZSMAIN,ZSHALF,DTDZS,TBQ, &
2699 !--- constants
2700 XLVM,CP,G0_P,CVW,STBOLT, &
2701 !--- output variables
2702 SNWEPRINT,SNHEIPRINT,RSM, &
2703 TSO,SOILT,SOILT1,TSNAV,QVG,QSG,QCG, &
2704 SMELT,SNOH,SNFLX,ILNB)
2705
2706 !********************************************************************
2707 ! Energy budget equation and heat diffusion eqn are
2708 ! solved here to obtain snow and soil temperatures
2709 !
2710 ! DELT - time step (s)
2711 ! ktau - numver of time step
2712 ! CONFLX - depth of constant flux layer (m)
2713 ! IME, JME, KME, NZS - dimensions of the domain
2714 ! NROOT - number of levels within the root zone
2715 ! PRCPMS - precipitation rate in m/s
2716 ! COTSO, RHTSO - coefficients for implicit solution of
2717 ! heat diffusion equation
2718 ! THDIF - thermal diffusivity (W/m/K)
2719 ! QSG,QVG,QCG - saturated mixing ratio, mixing ratio of
2720 ! water vapor and cloud at the ground
2721 ! surface, respectively (kg/kg)
2722 ! PATM - pressure [bar]
2723 ! QCATM,QVATM - cloud and water vapor mixing ratio
2724 ! at the first atm. level (kg/kg)
2725 ! EMISS,RNET - emissivity (0-1) of the ground surface and net
2726 ! radiation at the surface (W/m^2)
2727 ! QKMS - exchange coefficient for water vapor in the
2728 ! surface layer (m/s)
2729 ! TKMS - exchange coefficient for heat in the surface
2730 ! layer (m/s)
2731 ! PC - plant coefficient (resistance)
2732 ! RHO - density of atmosphere near surface (kg/m^3)
2733 ! VEGFRAC - greeness fraction (0-1)
2734 ! CAP - volumetric heat capacity (J/m^3/K)
2735 ! DRYCAN - dry fraction of vegetated area where
2736 ! transpiration may take place (0-1)
2737 ! WETCAN - fraction of vegetated area covered by canopy
2738 ! water (0-1)
2739 ! TRANSUM - transpiration function integrated over the
2740 ! rooting zone (m)
2741 ! DEW - dew in kg/m^2/s
2742 ! MAVAIL - fraction of maximum soil moisture in the top
2743 ! layer (0-1)
2744 ! ZSMAIN - main levels in soil (m)
2745 ! ZSHALF - middle of the soil layers (m)
2746 ! DTDZS - dt/(2.*dzshalf*dzmain)
2747 ! TBQ - table to define saturated mixing ration
2748 ! of water vapor for given temperature and pressure
2749 ! TSO - soil temperature (K)
2750 ! SOILT - skin temperature (K)
2751 !
2752 !*********************************************************************
2753
2754 IMPLICIT NONE
2755 !---------------------------------------------------------------------
2756 !--- input variables
2757
2758 INTEGER, INTENT(IN ) :: nroot,ktau,nzs , &
2759 nddzs !nddzs=2*(nzs-2)
2760
2761 INTEGER, INTENT(IN ) :: i,j,iland,isoil
2762 REAL, INTENT(IN ) :: DELT,CONFLX,PRCPMS , &
2763 RAINF,NEWSNOW,DELTSN,SNTH , &
2764 TABS,TRANSUM,SNWEPR
2765
2766 !--- 3-D Atmospheric variables
2767 REAL, &
2768 INTENT(IN ) :: PATM, &
2769 QVATM, &
2770 QCATM
2771 !--- 2-D variables
2772 REAL , &
2773 INTENT(IN ) :: GLW, &
2774 GSW, &
2775 RHO, &
2776 PC, &
2777 VEGFRAC, &
2778 QKMS, &
2779 TKMS
2780
2781 !--- soil properties
2782 REAL , &
2783 INTENT(IN ) :: &
2784 BCLH, &
2785 DQM, &
2786 PSIS, &
2787 QMIN
2788
2789 REAL, INTENT(IN ) :: CP, &
2790 CVW, &
2791 STBOLT, &
2792 XLVM, &
2793 G0_P
2794
2795
2796 REAL, DIMENSION(1:NZS), INTENT(IN) :: ZSMAIN, &
2797 ZSHALF, &
2798 THDIF, &
2799 CAP, &
2800 TRANF
2801
2802 REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS
2803
2804 REAL, DIMENSION(1:4001), INTENT(IN) :: TBQ
2805
2806
2807 !--- input/output variables
2808 !-------- 3-d soil moisture and temperature
2809 REAL, DIMENSION( 1:nzs ) , &
2810 INTENT(INOUT) :: TSO
2811
2812
2813 !-------- 2-d variables
2814 REAL , &
2815 INTENT(INOUT) :: DEW, &
2816 CST, &
2817 RHOSN, &
2818 EMISS, &
2819 MAVAIL, &
2820 QVG, &
2821 QSG, &
2822 QCG, &
2823 SNWE, &
2824 SNHEI, &
2825 SNOWFRAC, &
2826 SMELT, &
2827 SNOH, &
2828 SNFLX, &
2829 SOILT, &
2830 SOILT1, &
2831 TSNAV
2832
2833 REAL, INTENT(INOUT) :: DRYCAN, WETCAN
2834
2835 REAL, INTENT(OUT) :: RSM, &
2836 SNWEPRINT, &
2837 SNHEIPRINT
2838 INTEGER, INTENT(OUT) :: ilnb
2839 !--- Local variables
2840
2841
2842 INTEGER :: nzs1,nzs2,k,k1,kn,kk
2843
2844 REAL :: x,x1,x2,x4,dzstop,can,ft,sph, &
2845 tn,trans,umveg,denom
2846
2847 REAL :: cotsn,rhtsn,xsn1,ddzsn1,x1sn1,ftsnow,denomsn
2848
2849 REAL :: t3,upflux,xinet,ras, &
2850 xlmelt,rhocsn,thdifsn, &
2851 beta,epot,xsn,ddzsn,x1sn,d1sn,d2sn,d9sn,r22sn
2852
2853 REAL :: fso,fsn, &
2854 FKT,D1,D2,D9,D10,DID,R211,R21,R22,R6,R7,D11, &
2855 PI,H,FKQ,R210,AA,BB,PP,Q1,QS1,TS1,TQ2,TX2, &
2856 TDENOM,C,CC,AA1,RHCS,H1, &
2857 tsob, snprim, sh1, sh2, &
2858 smeltg,snohg,snodif,soh, &
2859 CMC2MS,TNOLD,QGOLD,SNOHGNEW
2860
2861 REAL, DIMENSION(1:NZS) :: transp,cotso,rhtso
2862 REAL :: edir1, &
2863 ec1, &
2864 ett1, &
2865 eeta, &
2866 s, &
2867 qfx, &
2868 hfx
2869
2870 REAL :: RNET,rsmfrac,soiltfrac,hsn
2871
2872 !-----------------------------------------------------------------
2873
2874 do k=1,nzs
2875 transp (k)=0.
2876 cotso (k)=0.
2877 rhtso (k)=0.
2878 enddo
2879
2880 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
2881 print *, 'SNOWTEMP: SNHEI,SNTH,SOILT1: ',SNHEI,SNTH,SOILT1,soilt
2882 ENDIF
2883 XLMELT=3.335E+5
2884 RHOCSN=2090.* RHOSN
2885 THDIFSN = 0.265/RHOCSN
2886 RAS=RHO*1.E-3
2887
2888 SOILTFRAC=SOILT
2889
2890 SMELT=0.
2891 SOH=0.
2892 SMELTG=0.
2893 SNOHG=0.
2894 SNODIF=0.
2895 RSM = 0.
2896 RSMFRAC = 0.
2897 fsn=1.
2898 fso=0.
2899 hsn=snhei
2900
2901 NZS1=NZS-1
2902 NZS2=NZS-2
2903
2904 QGOLD=QVG
2905 TNOLD=SOILT
2906 DZSTOP=1./(ZSMAIN(2)-ZSMAIN(1))
2907
2908 !******************************************************************************
2909 ! COEFFICIENTS FOR THOMAS ALGORITHM FOR TSO
2910 !******************************************************************************
2911 ! did=2.*(ZSMAIN(nzs)-ZSHALF(nzs))
2912 ! h1=DTDZS(8)*THDIF(nzs-1)*(ZSHALF(nzs)-ZSHALF(nzs-1))/did
2913 ! cotso(1)=h1/(1.+h1)
2914 ! rhtso(1)=(tso(nzs)+h1*(tso(nzs-1)-tso(nzs)))/
2915 ! 1 (1.+h1)
2916
2917 cotso(1)=0.
2918 rhtso(1)=TSO(NZS)
2919 DO 33 K=1,NZS2
2920 KN=NZS-K
2921 K1=2*KN-3
2922 X1=DTDZS(K1)*THDIF(KN-1)
2923 X2=DTDZS(K1+1)*THDIF(KN)
2924 FT=TSO(KN)+X1*(TSO(KN-1)-TSO(KN)) &
2925 -X2*(TSO(KN)-TSO(KN+1))
2926 DENOM=1.+X1+X2-X2*cotso(K)
2927 cotso(K+1)=X1/DENOM
2928 rhtso(K+1)=(FT+X2*rhtso(K))/DENOM
2929 33 CONTINUE
2930 !--- THE NZS element in COTSO and RHTSO will be for snow
2931 !--- There will be 2 layers in snow if it is deeper than DELTSN+SNTH
2932 IF(SNHEI.GE.SNTH) then
2933 ! if(snhei.le.DELTSN+DELTSN) then
2934 if(snhei.le.DELTSN+SNTH) then
2935 !-- 1-layer snow model
2936 ilnb=1
2937 snprim=snhei
2938 soilt1=tso(1)
2939 tsob=tso(1)
2940 XSN = DELT/2./(zshalf(2)+0.5*SNPRIM)
2941 DDZSN = XSN / SNPRIM
2942 X1SN = DDZSN * thdifsn
2943 X2 = DTDZS(1)*THDIF(1)
2944 FT = TSO(1)+X1SN*(SOILT-TSO(1)) &
2945 -X2*(TSO(1)-TSO(2))
2946 DENOM = 1. + X1SN + X2 -X2*cotso(NZS1)
2947 cotso(NZS)=X1SN/DENOM
2948 rhtso(NZS)=(FT+X2*rhtso(NZS1))/DENOM
2949 cotsn=cotso(NZS)
2950 rhtsn=rhtso(NZS)
2951 !*** Average temperature of snow pack (C)
2952 tsnav=0.5*(soilt+tso(1)) &
2953 -273.15
2954
2955 else
2956 !-- 2 layers in snow, SOILT1 is temperasture at DELTSN depth
2957 ilnb=2
2958 snprim=deltsn
2959 tsob=soilt1
2960 XSN = DELT/2./(0.5*SNHEI)
2961 XSN1= DELT/2./(zshalf(2)+0.5*(SNHEI-DELTSN))
2962 DDZSN = XSN / DELTSN
2963 DDZSN1 = XSN1 / (SNHEI-DELTSN)
2964 X1SN = DDZSN * thdifsn
2965 X1SN1 = DDZSN1 * thdifsn
2966 X2 = DTDZS(1)*THDIF(1)
2967 FT = TSO(1)+X1SN1*(SOILT1-TSO(1)) &
2968 -X2*(TSO(1)-TSO(2))
2969 DENOM = 1. + X1SN1 + X2 - X2*cotso(NZS1)
2970 cotso(nzs)=x1sn1/denom
2971 rhtso(nzs)=(ft+x2*rhtso(nzs1))/denom
2972 ftsnow = soilt1+x1sn*(soilt-soilt1) &
2973 -x1sn1*(soilt1-tso(1))
2974 denomsn = 1. + X1SN + X1SN1 - X1SN1*cotso(NZS)
2975 cotsn=x1sn/denomsn
2976 rhtsn=(ftsnow+X1SN1*rhtso(NZS))/denomsn
2977 !*** Average temperature of snow pack (C)
2978 tsnav=0.5/snhei*((soilt+soilt1)*deltsn &
2979 +(soilt1+tso(1))*(SNHEI-DELTSN)) &
2980 -273.15
2981 endif
2982 ENDIF
2983
2984 IF(SNHEI.LT.SNTH.AND.SNHEI.GT.0.) then
2985 !--- snow is too thin to be treated separately, therefore it
2986 !--- is combined with the first soil layer.
2987 fsn=SNHEI/(SNHEI+zsmain(2))
2988 fso=1.-fsn
2989 soilt1=tso(1)
2990 tsob=tso(2)
2991 snprim=SNHEI+zsmain(2)
2992 XSN = DELT/2./((zshalf(3)-zsmain(2))+0.5*snprim)
2993 DDZSN = XSN /snprim
2994 X1SN = DDZSN * (fsn*thdifsn+fso*thdif(1))
2995 X2=DTDZS(2)*THDIF(2)
2996 FT=TSO(2)+X1SN*(SOILT-TSO(2))- &
2997 X2*(TSO(2)-TSO(3))
2998 denom = 1. + x1sn + x2 - x2*cotso(nzs-2)
2999 cotso(nzs1) = x1sn/denom
3000 rhtso(nzs1)=(FT+X2*rhtso(NZS-2))/denom
3001 tsnav=0.5*(soilt+tso(1)) &
3002 -273.15
3003 ENDIF
3004
3005 !************************************************************************
3006 !--- THE HEAT BALANCE EQUATION (Smirnova et al. 1996, EQ. 21,26)
3007
3008 ETT1=0.
3009 EPOT=-QKMS*(QVATM-QSG)
3010 RHCS=CAP(1)
3011 H=MAVAIL
3012 IF(DEW.NE.0.)THEN
3013 DRYCAN=0.
3014 WETCAN=1.
3015 ENDIF
3016 TRANS=PC*TRANSUM*DRYCAN/ZSHALF(NROOT+1)
3017 CAN=WETCAN+TRANS
3018 UMVEG=1.-VEGFRAC
3019 FKT=TKMS
3020 D1=cotso(NZS1)
3021 D2=rhtso(NZS1)
3022 TN=SOILT
3023 D9=THDIF(1)*RHCS*dzstop
3024 D10=TKMS*CP*RHO
3025 R211=.5*CONFLX/DELT
3026 R21=R211*CP*RHO
3027 R22=.5/(THDIF(1)*DELT*dzstop**2)
3028 R6=EMISS *STBOLT*.5*TN**4
3029 R7=R6/TN
3030 D11=RNET+R6
3031
3032 IF(SNHEI.GE.SNTH) THEN
3033 ! if(snhei.le.DELTSN+DELTSN) then
3034 if(snhei.le.DELTSN+SNTH) then
3035 !--- 1-layer snow
3036 D1SN = cotso(NZS)
3037 D2SN = rhtso(NZS)
3038 else
3039 !--- 2-layer snow
3040 D1SN = cotsn
3041 D2SN = rhtsn
3042 endif
3043 D9SN= THDIFSN*RHOCSN / SNPRIM
3044 R22SN = SNPRIM*SNPRIM*0.5/(THDIFSN*DELT)
3045 ENDIF
3046
3047 IF(SNHEI.LT.SNTH.AND.SNHEI.GT.0.) then
3048 !--- thin snow is combined with soil
3049 D1SN = D1
3050 D2SN = D2
3051 D9SN = (fsn*THDIFSN*RHOCSN+fso*THDIF(1)*RHCS)/ &
3052 snprim
3053 R22SN = snprim*snprim*0.5 &
3054 /((fsn*THDIFSN+fso*THDIF(1))*delt)
3055 ENDIF
3056
3057 IF(SNHEI.eq.0.)then
3058 !--- all snow is sublimated
3059 D9SN = D9
3060 R22SN = R22
3061 D1SN = D1
3062 D2SN = D2
3063 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
3064 print *,' SNHEI = 0, D9SN,R22SN,D1SN,D2SN: ',D9SN,R22SN,D1SN,D2SN
3065 ENDIF
3066 ENDIF
3067
3068 !---- TDENOM for snow
3069
3070 TDENOM = D9SN*(1.-D1SN +R22SN)+D10+R21+R7 &
3071 +RAINF*CVW*PRCPMS &
3072 +RHOCSN*NEWSNOW/DELT
3073
3074 FKQ=QKMS*RHO
3075 R210=R211*RHO
3076 C=VEGFRAC*FKQ*CAN
3077 CC=C*XLVM/TDENOM
3078 AA=XLVM*(BETA*FKQ*UMVEG+R210)/TDENOM
3079 BB=(D10*TABS+R21*TN+XLVM*(QVATM* &
3080 (BETA*FKQ*UMVEG+C) &
3081 +R210*QVG)+D11+D9SN*(D2SN+R22SN*TN) &
3082 +RAINF*CVW*PRCPMS*max(273.15,TABS) &
3083 + RHOCSN*NEWSNOW/DELT*min(273.15,TABS) &
3084 )/TDENOM
3085 AA1=AA+CC
3086 PP=PATM*1.E3
3087 AA1=AA1/PP
3088 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
3089 print *,'VILKA-SNOW'
3090 print *,'tn,aa1,bb,pp,umveg,fkq,r210,vegfrac', &
3091 tn,aa1,bb,pp,umveg,fkq,r210,vegfrac
3092 ENDIF
3093
3094 CALL VILKA(TN,AA1,BB,PP,QS1,TS1,TBQ,KTAU,i,j,iland,isoil)
3095 TQ2=QVATM+QCATM
3096 TX2=TQ2*(1.-H)
3097 Q1=TX2+H*QS1
3098 !--- it is saturation over snow
3099 90 QVG=QS1
3100 QSG=QS1
3101 QCG=Q1-QS1
3102
3103 !--- SOILT - skin temperature
3104 SOILT=TS1
3105
3106 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
3107 print *,' AFTER VILKA-SNOW'
3108 print *,' TS1,QS1: ', ts1,qs1
3109 ENDIF
3110
3111 ! Solution for temperature at 7.5 cm depth and snow-soil interface
3112 IF(SNHEI.GE.SNTH) THEN
3113 ! if(snhei.gt.DELTSN+DELTSN) then
3114 if(snhei.gt.DELTSN+SNTH) then
3115 !-- 2-layer snow model
3116 SOILT1=rhtsn+cotsn*SOILT
3117 TSO(1)=rhtso(NZS)+cotso(NZS)*SOILT1
3118 tsob=soilt1
3119 else
3120 !-- 1 layer in snow
3121 TSO(1)=rhtso(NZS)+cotso(NZS)*SOILT
3122 SOILT1=TSO(1)
3123 tsob=tso(1)
3124 endif
3125 ELSE
3126 !-- all snow is evaporated
3127 TSO(1)=SOILT
3128 SOILT1=SOILT
3129 tsob=SOILT
3130 ENDIF
3131
3132 !---- Final solution for TSO
3133 DO K=2,NZS
3134 KK=NZS-K+1
3135 TSO(K)=rhtso(KK)+cotso(KK)*TSO(K-1)
3136 END DO
3137 !--- For thin snow layer combined with the top soil layer
3138 !--- TSO is computed by linear inmterpolation between SOILT
3139 !--- and TSO(2)
3140
3141 if(SNHEI.LT.SNTH.AND.SNHEI.GT.0.)then
3142 tso(1)=tso(2)+(soilt-tso(2))*fso
3143 SOILT1=TSO(1)
3144 tsob=tso(2)
3145 !!! tsob=tso(1)
3146 endif
3147
3148 !--- IF SOILT > 273.15 F then melting of snow can happen
3149 IF(SOILT.GE.273.15.AND.SNHEI.GT.0.) THEN
3150 !!! SOILT=273.15
3151 soiltfrac=snowfrac*273.15+(1.-snowfrac)*SOILT
3152 soilt=soiltfrac
3153 QSG= QSN(soilt,TBQ)/PP
3154 !!! QSG= QSN(273.15,TBQ)/PP
3155 QVG=QSG
3156 T3 = STBOLT*SOILT*SOILT*SOILT
3157 UPFLUX = T3 * SOILT
3158 XINET = EMISS*(GLW-UPFLUX)
3159 RNET = GSW + XINET
3160 EPOT = -QKMS*(QVATM-QSG)
3161 Q1=EPOT*RAS
3162
3163 IF (Q1.LE.0.) THEN
3164 ! --- condensation
3165 DEW=-EPOT
3166 DO K=1,NZS
3167 TRANSP(K)=0.
3168 ENDDO
3169
3170 QFX= XLVM*RHO*DEW
3171 ELSE
3172 ! --- evaporation
3173 DO K=1,NROOT
3174 TRANSP(K)=-VEGFRAC*q1 &
3175 *PC*TRANF(K)*DRYCAN/zshalf(NROOT+1)
3176 IF(TRANSP(K).GT.0.) TRANSP(K)=0.
3177 ETT1=ETT1-TRANSP(K)
3178 ENDDO
3179 DO k=nroot+1,nzs
3180 transp(k)=0.
3181 enddo
3182
3183 EDIR1 = Q1*UMVEG * BETA
3184 EC1 = Q1 * WETCAN *VEGFRAC
3185 CMC2MS=CST/DELT
3186 EC1=MIN(CMC2MS,EC1)
3187 EETA = (EDIR1 + EC1 + ETT1)*1.E3
3188 ! to convert from kg m-2 s-1 to m s-1: 1/rho water=1.e-3************
3189 QFX= - XLVM * EETA
3190 ENDIF
3191
3192 HFX=D10*(TABS-soilt)
3193 !!! HFX=D10*(TABS-273.15)
3194
3195 IF(SNHEI.GE.SNTH)then
3196 SOH=thdifsn*RHOCSN*(soilt-TSOB)/SNPRIM
3197 ! SOH=thdifsn*RHOCSN*(273.15-TSOB)/SNPRIM
3198 SNFLX=SOH
3199 ELSE
3200 SOH=(fsn*thdifsn*rhocsn+fso*thdif(1)*rhcs)* &
3201 (soilt-TSOB)/snprim
3202 !!! (273.15-TSOB)/snprim
3203 SNFLX=SOH
3204 ENDIF
3205
3206 X= (R21+D9SN*R22SN)*(soilt-TNOLD) + &
3207 !!! X= (R21+D9SN*R22SN)*(273.15-TNOLD) + &
3208 XLVM*R210*(QSG-QGOLD)
3209 !-- SNOH is energy flux of snow phase change
3210 SNOH=RNET+QFX +HFX &
3211 +RHOCSN*NEWSNOW/DELT*(min(273.15,TABS)-TN) &
3212 -SOH-X+RAINF*CVW*PRCPMS* &
3213 (max(273.15,TABS)-TN)
3214 SNOH=AMAX1(0.,SNOH)
3215 !-- SMELT is speed of melting in M/S
3216 SMELT= SNOH /XLMELT*1.E-3
3217 SMELT=AMIN1(SMELT,SNWEPR/DELT-BETA*EPOT*RAS)
3218 ! SMELT=AMIN1(SMELT,SNWEPR/DELT-BETA*EPOT*RAS*UMVEG)
3219
3220 SNOHGNEW=SMELT*XLMELT*1.E3
3221 SNODIF=AMAX1(0.,(SNOH-SNOHGNEW))
3222
3223 SNOH=SNOHGNEW
3224 ! SNOHSMELT*XLMELT*1.E3
3225
3226 !*** From Koren et al. (1999) 13% of snow melt stays in the snow pack
3227 !!! rsm=0.13*smelt*delt
3228 if(snwe.gt.0.) then
3229 rsmfrac=min(0.18,(max(0.08,0.10/snwe*0.13)))
3230 else
3231 rsmfrac=0.13
3232 endif
3233
3234 rsm=rsmfrac*smelt*delt
3235 SMELT=SMELT-rsm/delt
3236
3237 !-- correction of liquid equivalent of snow depth
3238 !-- due to evaporation and snow melt
3239 SNWE = AMAX1(0.,(SNWEPR- &
3240 (SMELT+BETA*EPOT*RAS)*DELT &
3241 ! (SMELT+BETA*EPOT*RAS*UMVEG)*DELT &
3242 ) )
3243
3244 !--- If all snow melts, then 13% of snow melt we kept in the
3245 !--- snow pack should be added back to snow melt and infiltrate
3246 !--- into soil.
3247 if(snwe.le.rsm) then
3248 smelt=smelt+rsm/delt
3249 snwe=0.
3250 rsm=0.
3251 SOILT=SNODIF*DELT/RHCS*ZSHALF(2) &
3252 +soiltfrac
3253 !!! +273.15
3254 else
3255 !*** Correct snow density on effect of snow melt, melted
3256 !*** from the top of the snow. 13% of melted water
3257 !*** remains in the pack and changes its density.
3258 !*** Eq. 9 (with my correction) in Koren et al. (1999)
3259
3260 if(snwe.gt.snth*rhosn*1.e-3) then
3261 xsn=(rhosn*(snwe-rsm)+1.e3*rsm)/ &
3262 snwe
3263 rhosn=MIN(XSN,400.)
3264
3265 RHOCSN=2090.* RHOSN
3266 thdifsn = 0.265/RHOCSN
3267 endif
3268
3269 SOILT=SOILTFRAC
3270
3271 endif
3272
3273 !--- If there is no snow melting then just evaporation
3274 !--- or condensation cxhanges SNWE
3275 ELSE
3276 EPOT=-QKMS*(QVATM-QSG)
3277 SNWE = AMAX1(0.,(SNWEPR- &
3278 BETA*EPOT*RAS*DELT))
3279 ! BETA*EPOT*RAS*UMVEG*DELT))
3280
3281 ENDIF
3282 !*** Correct snow density on effect of snow melt, melted
3283 !*** from the top of the snow. 13% of melted water
3284 !*** remains in the pack and changes its density.
3285 !*** Eq. 9 (with my correction) in Koren et al. (1999)
3286
3287 SNHEI=SNWE *1.E3 / RHOSN
3288
3289 !-- Snow melt from the top is done. But if ground surface temperature
3290 !-- is above freezing snow can melt from the bottom. The following
3291 !-- piece of code will check if bottom melting is possible.
3292
3293 IF(TSO(1).GE.273.15.AND.SNHEI.GT.0.) THEN
3294 soiltfrac=snowfrac*273.15+(1.-snowfrac)*TSO(1)
3295
3296 SNOHG=(TSO(1)-soiltfrac)*(RHCS*zshalf(2)+ &
3297 RHOCSN*0.5*SNHEI) / DELT
3298 SNOHG=AMAX1(0.,SNOHG)
3299 SNODIF=0.
3300 ! TSO(1)=273.15
3301 SMELTG=SNOHG/XLMELT*1.E-3
3302 ! SMELTG=AMIN1(SMELTG,SNWE/DELT)
3303 if(SNWE-SMELTG*DELT.ge.rsm) then
3304 ! SNWE = SNWE-SMELTG*DELT
3305 SNWE = AMAX1(0.,SNWE-SMELTG*DELT)
3306 else
3307 smeltg=snwe/delt
3308 snwe=0.
3309 rsm=0.
3310 endif
3311
3312 SNOHGNEW=SMELTG*XLMELT*1.e3
3313 SNODIF=AMAX1(0.,(SNOHG-SNOHGNEW))
3314 TSO(1)=soiltfrac
3315 if(snwe.eq.0.)then
3316 TSO(1)=SNODIF*DELT/RHCS*zshalf(2) + soiltfrac
3317 !!! TSO(1)=SNODIF*DELT/RHCS*zshalf(2) + 273.15
3318 endif
3319
3320 SMELT=SMELT+SMELTG
3321 SNOH=SNOH+SNOHGNEW
3322
3323 ENDIF
3324
3325 SNHEI=SNWE *1.E3 / RHOSN
3326
3327 snweprint=snwe
3328 ! &
3329 !--- if VEGFRAC.ne.0. then some snow stays on the canopy
3330 !--- and should be added to SNWE for water conservation
3331 ! 4 Nov 07 +VEGFRAC*cst
3332 snheiprint=snweprint*1.E3 / RHOSN
3333
3334 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
3335 print *, 'snweprint : ',snweprint
3336 print *, 'D9SN,SOILT,TSOB : ', D9SN,SOILT,TSOB
3337 ENDIF
3338 !--- Compute flux in the top snow layer
3339 SNFLX=D9SN*(SOILT-TSOB)
3340 IF(SNHEI.GT.0.) THEN
3341 if(ilnb.gt.1) then
3342 tsnav=0.5/snhei*((soilt+soilt1)*deltsn &
3343 +(soilt1+tso(1))*(SNHEI-DELTSN)) &
3344 -273.15
3345 else
3346 tsnav=0.5*(soilt+tso(1)) - 273.15
3347 endif
3348 ENDIF
3349
3350 ! return
3351 ! end
3352 !------------------------------------------------------------------------
3353 END SUBROUTINE SNOWTEMP
3354 !------------------------------------------------------------------------
3355
3356
3357 SUBROUTINE SOILMOIST ( &
3358 !--input parameters
3359 DELT,NZS,NDDZS,DTDZS,DTDZS2, &
3360 ZSMAIN,ZSHALF,DIFFU,HYDRO, &
3361 QSG,QVG,QCG,QCATM,QVATM,PRCP, &
3362 QKMS,TRANSP,DRIP, &
3363 DEW,SMELT,SOILICE,VEGFRAC, &
3364 !--soil properties
3365 DQM,QMIN,REF,KSAT,RAS,INFMAX, &
3366 !--output
3367 SOILMOIS,MAVAIL,RUNOFF,RUNOFF2,INFILTRP)
3368 !*************************************************************************
3369 ! moisture balance equation and Richards eqn.
3370 ! are solved here
3371 !
3372 ! DELT - time step (s)
3373 ! IME,JME,NZS - dimensions of soil domain
3374 ! ZSMAIN - main levels in soil (m)
3375 ! ZSHALF - middle of the soil layers (m)
3376 ! DTDZS - dt/(2.*dzshalf*dzmain)
3377 ! DTDZS2 - dt/(2.*dzshalf)
3378 ! DIFFU - diffusional conductivity (m^2/s)
3379 ! HYDRO - hydraulic conductivity (m/s)
3380 ! QSG,QVG,QCG - saturated mixing ratio, mixing ratio of
3381 ! water vapor and cloud at the ground
3382 ! surface, respectively (kg/kg)
3383 ! QCATM,QVATM - cloud and water vapor mixing ratio
3384 ! at the first atm. level (kg/kg)
3385 ! PRCP - precipitation rate in m/s
3386 ! QKMS - exchange coefficient for water vapor in the
3387 ! surface layer (m/s)
3388 ! TRANSP - transpiration from the soil layers (m/s)
3389 ! DRIP - liquid water dripping from the canopy to soil (m)
3390 ! DEW - dew in kg/m^2s
3391 ! SMELT - melting rate in m/s
3392 ! SOILICE - volumetric content of ice in soil (m^3/m^3)
3393 ! VEGFRAC - greeness fraction (0-1)
3394 ! RAS - ration of air density to soil density
3395 ! INFMAX - maximum infiltration rate (kg/m^2/s)
3396 !
3397 ! SOILMOIS - volumetric soil moisture, 6 levels (m^3/m^3)
3398 ! MAVAIL - fraction of maximum soil moisture in the top
3399 ! layer (0-1)
3400 ! RUNOFF - surface runoff (m/s)
3401 ! RUNOFF2 - underground runoff (m)
3402 ! INFILTRP - point infiltration flux into soil (m/s)
3403 ! /(snow bottom runoff) (mm/s)
3404 !
3405 ! COSMC, RHSMC - coefficients for implicit solution of
3406 ! Richards equation
3407 !******************************************************************
3408 IMPLICIT NONE
3409 !------------------------------------------------------------------
3410 !--- input variables
3411 REAL, INTENT(IN ) :: DELT
3412 INTEGER, INTENT(IN ) :: NZS,NDDZS
3413
3414 ! input variables
3415
3416 REAL, DIMENSION(1:NZS), INTENT(IN ) :: ZSMAIN, &
3417 ZSHALF, &
3418 DIFFU, &
3419 HYDRO, &
3420 TRANSP, &
3421 SOILICE, &
3422 DTDZS2
3423
3424 REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS
3425
3426 REAL, INTENT(IN ) :: QSG,QVG,QCG,QCATM,QVATM , &
3427 QKMS,VEGFRAC,DRIP,PRCP , &
3428 DEW,SMELT , &
3429 DQM,QMIN,REF,KSAT,RAS
3430
3431 ! output
3432
3433 REAL, DIMENSION( 1:nzs ) , &
3434
3435 INTENT(INOUT) :: SOILMOIS
3436
3437 REAL, INTENT(INOUT) :: MAVAIL,RUNOFF,RUNOFF2,INFILTRP, &
3438 INFMAX
3439
3440 ! local variables
3441
3442 REAL, DIMENSION( 1:nzs ) :: COSMC,RHSMC
3443
3444 REAL :: DZS,R1,R2,R3,R4,R5,R6,R7,R8,R9,R10
3445 REAL :: REFKDT,REFDK,DELT1,F1MAX,F2MAX
3446 REAL :: F1,F2,FD,KDT,VAL,DDT,PX
3447 REAL :: QQ,UMVEG,INFMAX1,TRANS
3448 REAL :: TOTLIQ,FLX,FLXSAT,QTOT
3449 REAL :: DID,X1,X2,X4,DENOM,Q2,Q4
3450 REAL :: dice,fcr,acrt,frzx,sum,cvfrz
3451
3452 INTEGER :: NZS1,NZS2,K,KK,K1,KN,ialp1,jj,jk
3453
3454 !******************************************************************************
3455 ! COEFFICIENTS FOR THOMAS ALGORITHM FOR SOILMOIS
3456 !******************************************************************************
3457 NZS1=NZS-1
3458 NZS2=NZS-2
3459
3460 118 format(6(10Pf23.19))
3461
3462 do k=1,nzs
3463 cosmc(k)=0.
3464 rhsmc(k)=0.
3465 enddo
3466
3467 DID=(ZSMAIN(NZS)-ZSHALF(NZS))*2.
3468 X1=ZSMAIN(NZS)-ZSMAIN(NZS1)
3469 ! DENOM=DID/DELT+DIFFU(NZS1)/X1
3470 ! COSMC(1)=DIFFU(NZS1)/X1/DENOM
3471 ! RHSMC(1)=(SOILMOIS(NZS)*DID/DELT
3472 ! 1 +TRANSP(NZS)-(HYDRO(NZS)*SOILMOIS(NZS)
3473 ! 1 -HYDRO(NZS1)*SOILMOIS(NZS1))*DID
3474 ! 1 /X1) /DENOM
3475
3476 DENOM=(1.+DIFFU(nzs1)/X1/DID*DELT+HYDRO(NZS)/(2.*DID)*DELT)
3477 COSMC(1)=DELT*(DIFFU(nzs1)/DID/X1 &
3478 +HYDRO(NZS1)/2./DID)/DENOM
3479 RHSMC(1)=(SOILMOIS(NZS)+TRANSP(NZS)*DELT/ &
3480 DID)/DENOM
3481
3482 DO 330 K=1,NZS2
3483 KN=NZS-K
3484 K1=2*KN-3
3485 X4=2.*DTDZS(K1)*DIFFU(KN-1)
3486 X2=2.*DTDZS(K1+1)*DIFFU(KN)
3487 Q4=X4+HYDRO(KN-1)*DTDZS2(KN-1)
3488 Q2=X2-HYDRO(KN+1)*DTDZS2(KN-1)
3489 DENOM=1.+X2+X4-Q2*COSMC(K)
3490 COSMC(K+1)=Q4/DENOM
3491 330 RHSMC(K+1)=(SOILMOIS(KN)+Q2*RHSMC(K) &
3492 +TRANSP(KN) &
3493 /(ZSHALF(KN+1)-ZSHALF(KN)) &
3494 *DELT)/DENOM
3495
3496 ! --- MOISTURE BALANCE BEGINS HERE
3497
3498 TRANS=TRANSP(1)
3499 UMVEG=1.-VEGFRAC
3500
3501 RUNOFF=0.
3502 RUNOFF2=0.
3503 DZS=ZSMAIN(2)
3504 R1=COSMC(NZS1)
3505 R2= RHSMC(NZS1)
3506 R3=DIFFU(1)/DZS
3507 R4=R3+HYDRO(1)*.5
3508 R5=R3-HYDRO(2)*.5
3509 R6=QKMS*RAS
3510 !-- Total liquid water available on the top of soil domain
3511 !-- Without snow - 3 sources of water: precipitation,
3512 !-- water dripping from the canopy and dew
3513 !-- With snow - only one source of water - snow melt
3514
3515 ! print *,'PRCP,DRIP,DEW,umveg,ras,smelt',
3516 ! 1 PRCP,DRIP,DEW,umveg,ras,smelt
3517 ! if (drip.ne.0.) then
3518 ! print *,'DRIP non-zero'
3519 ! write(6,191) drip
3520 ! write (6,191)soilmois(1)
3521 ! write (6,191)soilmois(2)
3522 ! endif
3523 191 format (f23.19)
3524
3525 TOTLIQ=UMVEG*PRCP-DRIP/DELT-UMVEG*DEW*RAS-SMELT
3526
3527
3528 FLX=TOTLIQ
3529 INFILTRP=TOTLIQ
3530
3531 ! ----------- FROZEN GROUND VERSION -------------------------
3532 ! REFERENCE FROZEN GROUND PARAMETER, CVFRZ, IS A SHAPE PARAMETER OF
3533 ! AREAL DISTRIBUTION FUNCTION OF SOIL ICE CONTENT WHICH EQUALS 1/CV.
3534 ! CV IS A COEFFICIENT OF SPATIAL VARIATION OF SOIL ICE CONTENT.
3535 ! BASED ON FIELD DATA CV DEPENDS ON AREAL MEAN OF FROZEN DEPTH, AND IT
3536 ! CLOSE TO CONSTANT = 0.6 IF AREAL MEAN FROZEN DEPTH IS ABOVE 20 CM.
3537 ! THAT IS WHY PARAMETER CVFRZ = 3 (INT{1/0.6*0.6})
3538 !
3539 ! Current logic doesn't allow CVFRZ be bigger than 3
3540 CVFRZ = 3.
3541
3542 !-- SCHAAKE/KOREN EXPRESSION for calculation of max infiltration
3543 REFKDT=3.
3544 REFDK=3.4341E-6
3545 DELT1=DELT/86400.
3546 F1MAX=DQM*ZSHALF(2)
3547 F2MAX=DQM*(ZSHALF(3)-ZSHALF(2))
3548 F1=F1MAX*(1.-SOILMOIS(1)/DQM)
3549 F2=F2MAX*(1.-SOILMOIS(2)/DQM)
3550 FD=F1+F2
3551 KDT=REFKDT*KSAT/REFDK
3552 VAL=(1.-EXP(-KDT*DELT1))
3553 DDT = FD*VAL
3554 PX= - TOTLIQ * DELT
3555 IF(PX.LT.0.0) PX = 0.0
3556 if(ddt.eq.0.) then
3557 infmax1=ksat
3558 else
3559 INFMAX1 = (PX*(DDT/(PX+DDT)))/DELT
3560 INFMAX1 = MIN(INFMAX1, KSAT)
3561 endif
3562 ! print *,'INFMAX1=,ksat',infmax1,ksat,f1,f2,kdt,val,ddt,px
3563 ! ----------- FROZEN GROUND VERSION --------------------------
3564 ! REDUCTION OF INFILTRATION BASED ON FROZEN GROUND PARAMETERS
3565 !
3566 ! ------------------------------------------------------------------
3567
3568 DICE = soilice(1)*zshalf(2)
3569 DO K=2,NZS1
3570 DICE = DICE + ( ZSHALF(K+1) - ZSHALF(K) ) * soilice(k)
3571 ENDDO
3572 FRZX= 0.28*((dqm+qmin)/ref) * (0.400 / 0.482)
3573 FCR = 1.
3574 IF ( DICE .GT. 1.E-2) THEN
3575 ACRT = CVFRZ * FRZX / DICE
3576 SUM = 1.
3577 IALP1 = CVFRZ - 1
3578 DO JK = 1,IALP1
3579 K = 1
3580 DO JJ = JK+1, IALP1
3581 K = K * JJ
3582 END DO
3583 SUM = SUM + (ACRT ** ( CVFRZ-JK)) / FLOAT (K)
3584 END DO
3585 FCR = 1. - EXP(-ACRT) * SUM
3586 END IF
3587 ! print *,'FCR--------',fcr
3588 INFMAX1 = INFMAX1* FCR
3589 INFMAX1 = MIN(INFMAX1, KSAT)
3590 ! -------------------------------------------------------------------
3591
3592 INFMAX = MIN(INFMAX,INFMAX1)
3593 !----
3594 IF (-TOTLIQ.GE.INFMAX)THEN
3595 RUNOFF=-TOTLIQ-INFMAX
3596 FLX=-INFMAX
3597 ENDIF
3598 ! INFILTRP is total infiltration flux in M/S
3599 INFILTRP=FLX
3600 ! print *,'PRCIP',infiltrp,flx,infmax
3601 ! Solution of moisture budget
3602 R7=.5*DZS/DELT
3603 R4=R4+R7
3604 FLX=FLX-SOILMOIS(1)*R7
3605 R8=UMVEG*R6
3606 QTOT=QVATM+QCATM
3607 R9=TRANS
3608 R10=QTOT-QSG
3609 !-- evaporation regime
3610 IF(R10.LE.0.) THEN
3611 QQ=(R5*R2-FLX+R9)/(R4-R5*R1-R10*R8/(REF-QMIN))
3612 FLXSAT=-DQM*(R4-R5*R1-R10*R8/(REF-QMIN)) &
3613 +R5*R2+R9
3614 ELSE
3615 !-- dew formation regime
3616 QQ=(R2*R5-FLX+R8*(QTOT-QCG-QVG)+R9)/(R4-R1*R5)
3617 FLXSAT=-DQM*(R4-R1*R5)+R2*R5+R8*(QTOT-QVG-QCG)+R9
3618 END IF
3619
3620 IF(QQ.LT.0.) THEN
3621 SOILMOIS(1)=0.
3622
3623 ELSE IF(QQ.GT.DQM) THEN
3624 !-- saturation
3625 SOILMOIS(1)=DQM
3626 RUNOFF2=runoff2+(FLXSAT-FLX)*DELT
3627 RUNOFF=RUNOFF+(FLXSAT-FLX)
3628 ELSE
3629 SOILMOIS(1)=max(1.e-8,QQ)
3630 END IF
3631
3632 !--- FINAL SOLUTION FOR SOILMOIS
3633 ! DO K=2,NZS
3634 DO K=2,NZS-1
3635 KK=NZS-K+1
3636 QQ=COSMC(KK)*SOILMOIS(K-1)+RHSMC(KK)
3637
3638 IF (QQ.LT.0.) THEN
3639 SOILMOIS(K)=0.
3640
3641 ELSE IF(QQ.GT.DQM) THEN
3642 !-- saturation
3643 SOILMOIS(K)=DQM
3644 IF(K.EQ.NZS)THEN
3645 RUNOFF2=RUNOFF2+(QQ-DQM)*(ZSMAIN(K)-ZSHALF(K))
3646 ELSE
3647 RUNOFF2=RUNOFF2+(QQ-DQM)*(ZSMAIN(K+1)-ZSHALF(K))
3648 ENDIF
3649 ELSE
3650 SOILMOIS(K)=max(1.e-8,QQ)
3651 END IF
3652 END DO
3653
3654 ! MAVAIL=min(1.,SOILMOIS(1)/(REF-QMIN))
3655 MAVAIL=min(1.,SOILMOIS(1)/DQM)
3656 if (MAVAIL.EQ.0.) MAVAIL=.00001
3657
3658 ! RETURN
3659 ! END
3660 !-------------------------------------------------------------------
3661 END SUBROUTINE SOILMOIST
3662 !-------------------------------------------------------------------
3663
3664
3665 SUBROUTINE SOILPROP( &
3666 !--- input variables
3667 nzs,fwsat,lwsat,tav,keepfr, &
3668 soilmois,soiliqw,soilice, &
3669 soilmoism,soiliqwm,soilicem, &
3670 !--- soil fixed fields
3671 QWRTZ,rhocs,dqm,qmin,psis,bclh,ksat, &
3672 !--- constants
3673 riw,xlmelt,CP,G0_P,cvw,ci, &
3674 kqwrtz,kice,kwt, &
3675 !--- output variables
3676 thdif,diffu,hydro,cap)
3677
3678 !******************************************************************
3679 ! SOILPROP computes thermal diffusivity, and diffusional and
3680 ! hydraulic condeuctivities
3681 !******************************************************************
3682 ! NX,NY,NZS - dimensions of soil domain
3683 ! FWSAT, LWSAT - volumetric content of frozen and liquid water
3684 ! for saturated condition at given temperatures (m^3/m^3)
3685 ! TAV - temperature averaged for soil layers (K)
3686 ! SOILMOIS - volumetric soil moisture at the main soil levels (m^3/m^3)
3687 ! SOILMOISM - volumetric soil moisture averaged for layers (m^3/m^3)
3688 ! SOILIQWM - volumetric liquid soil moisture averaged for layers (m^3/m^3)
3689 ! SOILICEM - volumetric content of soil ice averaged for layers (m^3/m^3)
3690 ! THDIF - thermal diffusivity for soil layers (W/m/K)
3691 ! DIFFU - diffusional conductivity (m^2/s)
3692 ! HYDRO - hydraulic conductivity (m/s)
3693 ! CAP - volumetric heat capacity (J/m^3/K)
3694 !
3695 !******************************************************************
3696
3697 IMPLICIT NONE
3698 !-----------------------------------------------------------------
3699
3700 !--- soil properties
3701 INTEGER, INTENT(IN ) :: NZS
3702 REAL , &
3703 INTENT(IN ) :: RHOCS, &
3704 BCLH, &
3705 DQM, &
3706 KSAT, &
3707 PSIS, &
3708 QWRTZ, &
3709 QMIN
3710
3711 REAL, DIMENSION( 1:nzs ) , &
3712 INTENT(IN ) :: SOILMOIS, &
3713 keepfr
3714
3715
3716 REAL, INTENT(IN ) :: CP, &
3717 CVW, &
3718 RIW, &
3719 kqwrtz, &
3720 kice, &
3721 kwt, &
3722 XLMELT, &
3723 G0_P
3724
3725
3726
3727 !--- output variables
3728 REAL, DIMENSION(1:NZS) , &
3729 INTENT(INOUT) :: cap,diffu,hydro , &
3730 thdif,tav , &
3731 soilmoism , &
3732 soiliqw,soilice , &
3733 soilicem,soiliqwm , &
3734 fwsat,lwsat
3735
3736 !--- local variables
3737 REAL, DIMENSION(1:NZS) :: hk,detal,kasat,kjpl
3738
3739 REAL :: x,x1,x2,x4,ws,wd,fact,fach,facd,psif,ci
3740 REAL :: tln,tavln,tn,pf,a,am,ame,h
3741 INTEGER :: nzs1,k
3742
3743 !-- for Johansen thermal conductivity
3744 REAL :: kzero,gamd,kdry,kas,x5,sr,ke
3745
3746
3747 nzs1=nzs-1
3748
3749 !-- Constants for Johansen (1975) thermal conductivity
3750 kzero =2. ! if qwrtz > 0.2
3751
3752
3753 do k=1,nzs
3754 detal (k)=0.
3755 kasat (k)=0.
3756 kjpl (k)=0.
3757 hk (k)=0.
3758 enddo
3759
3760 ws=dqm+qmin
3761 x1=xlmelt/(g0_p*psis)
3762 x2=x1/bclh*ws
3763 x4=(bclh+1.)/bclh
3764 !--- Next 3 lines are for Johansen thermal conduct.
3765 gamd=(1.-ws)*2700.
3766 kdry=(0.135*gamd+64.7)/(2700.-0.947*gamd)
3767 kas=kqwrtz**qwrtz*kzero**(1.-qwrtz)
3768
3769 DO K=1,NZS1
3770 tn=tav(k) - 273.15
3771 wd=ws - riw*soilicem(k)
3772 psif=psis*100.*(wd/(soiliqwm(k)+qmin))**bclh &
3773 * (ws/wd)**3.
3774 !--- PSIF should be in [CM] to compute PF
3775 pf=log10(abs(psif))
3776 fact=1.+riw*soilicem(k)
3777 !--- HK is for McCumber thermal conductivity
3778 IF(PF.LE.5.2) THEN
3779 HK(K)=420.*EXP(-(PF+2.7))*fact
3780 ELSE
3781 HK(K)=.1744*fact
3782 END IF
3783
3784 IF(soilicem(k).NE.0.AND.TN.LT.0.) then
3785 !--- DETAL is taking care of energy spent on freezing or released from
3786 ! melting of soil water
3787
3788 DETAL(K)=273.15*X2/(TAV(K)*TAV(K))* &
3789 (TAV(K)/(X1*TN))**X4
3790
3791 if(keepfr(k).eq.1.) then
3792 detal(k)=0.
3793 endif
3794
3795 ENDIF
3796
3797 !--- Next 10 lines calculate Johansen thermal conductivity KJPL
3798 kasat(k)=kas**(1.-ws)*kice**fwsat(k) &
3799 *kwt**lwsat(k)
3800
3801 X5=(soilmoism(k)+qmin)/ws
3802 if(soilicem(k).eq.0.) then
3803 sr=max(0.101,x5)
3804 ke=log10(sr)+1.
3805 !--- next 2 lines - for coarse soils
3806 ! sr=max(0.0501,x5)
3807 ! ke=0.7*log10(sr)+1.
3808 else
3809 ke=x5
3810 endif
3811
3812 kjpl(k)=ke*(kasat(k)-kdry)+kdry
3813
3814 !--- CAP -volumetric heat capacity
3815 CAP(K)=(1.-WS)*RHOCS &
3816 + (soiliqwm(K)+qmin)*CVW &
3817 + soilicem(K)*CI &
3818 + (dqm-soilmoism(k))*CP*1.2 &
3819 - DETAL(K)*1.e3*xlmelt
3820
3821 a=RIW*soilicem(K)
3822
3823 if((ws-a).lt.0.12)then
3824 diffu(K)=0.
3825 else
3826 H=max(0.,(soilmoism(K)-a)/(max(1.e-8,(dqm-a))))
3827 facd=1.
3828 if(a.ne.0.)facd=1.-a/max(1.e-8,soilmoism(K))
3829 ame=max(1.e-8,dqm-riw*soilicem(K))
3830 !--- DIFFU is diffusional conductivity of soil water
3831 diffu(K)=-BCLH*KSAT*PSIS/ame* &
3832 (dqm/ame)**3. &
3833 *H**(BCLH+2.)*facd
3834 endif
3835
3836 ! diffu(K)=-BCLH*KSAT*PSIS/dqm &
3837 ! *H**(BCLH+2.)
3838
3839
3840 !--- thdif - thermal diffusivity
3841 ! thdif(K)=HK(K)/CAP(K)
3842 !--- Use thermal conductivity from Johansen (1975)
3843 thdif(K)=KJPL(K)/CAP(K)
3844
3845 END DO
3846
3847 DO K=1,NZS
3848
3849 if((ws-riw*soilice(k)).lt.0.12)then
3850 hydro(k)=0.
3851 else
3852 fach=1.
3853 if(soilice(k).ne.0.) &
3854 fach=1.-riw*soilice(k)/max(1.e-8,soilmois(k))
3855 am=max(1.e-8,dqm-riw*soilice(k))
3856 !--- HYDRO is hydraulic conductivity of soil water
3857 hydro(K)=KSAT/am* &
3858 (soiliqw(K)/am) &
3859 **(2.*BCLH+2.) &
3860 * fach
3861 endif
3862
3863 ENDDO
3864
3865 ! RETURN
3866 ! END
3867
3868 !-----------------------------------------------------------------------
3869 END SUBROUTINE SOILPROP
3870 !-----------------------------------------------------------------------
3871
3872
3873 SUBROUTINE TRANSF( &
3874 !--- input variables
3875 nzs,nroot,soiliqw, &
3876 !--- soil fixed fields
3877 dqm,qmin,ref,wilt,zshalf, &
3878 !--- output variables
3879 tranf,transum)
3880
3881 !-------------------------------------------------------------------
3882 !--- TRANF(K) - THE TRANSPIRATION FUNCTION (Smirnova et al. 1996, EQ. 18,19)
3883 !*******************************************************************
3884 ! NX,NY,NZS - dimensions of soil domain
3885 ! SOILIQW - volumetric liquid soil moisture at the main levels (m^3/m^3)
3886 ! TRANF - the transpiration function at levels (m)
3887 ! TRANSUM - transpiration function integrated over the rooting zone (m)
3888 !
3889 !*******************************************************************
3890 IMPLICIT NONE
3891 !-------------------------------------------------------------------
3892
3893 !--- input variables
3894
3895 INTEGER, INTENT(IN ) :: nroot,nzs
3896
3897 !--- soil properties
3898 REAL , &
3899 INTENT(IN ) :: DQM, &
3900 QMIN, &
3901 REF, &
3902 WILT
3903
3904 REAL, DIMENSION(1:NZS), INTENT(IN) :: soiliqw, &
3905 ZSHALF
3906
3907 !-- output
3908 REAL, DIMENSION(1:NZS), INTENT(OUT) :: TRANF
3909 REAL, INTENT(OUT) :: TRANSUM
3910
3911 !-- local variables
3912 REAL :: totliq, did
3913 INTEGER :: k
3914
3915 !-- for non-linear root distribution
3916 REAL :: gx,sm1,sm2,sm3,sm4,ap0,ap1,ap2,ap3,ap4
3917 REAL, DIMENSION(1:NZS) :: PART
3918 !--------------------------------------------------------------------
3919
3920 do k=1,nzs
3921 part(k)=0.
3922 enddo
3923
3924 transum=0.
3925 totliq=soiliqw(1)+qmin
3926 sm1=totliq
3927 sm2=sm1*sm1
3928 sm3=sm2*sm1
3929 sm4=sm3*sm1
3930 ap0=0.299
3931 ap1=-8.152
3932 ap2=61.653
3933 ap3=-115.876
3934 ap4=59.656
3935 gx=ap0+ap1*sm1+ap2*sm2+ap3*sm3+ap4*sm4
3936 if(totliq.ge.ref) gx=1.
3937 if(totliq.le.0.) gx=0.
3938 if(gx.gt.1.) gx=1.
3939 if(gx.lt.0.) gx=0.
3940 DID=zshalf(2)
3941 part(1)=DID*gx
3942 IF(TOTLIQ.GT.REF) THEN
3943 TRANF(1)=DID
3944 ELSE IF(TOTLIQ.LE.WILT) THEN
3945 TRANF(1)=0.
3946 ELSE
3947 TRANF(1)=(TOTLIQ-WILT)/(REF-WILT)*DID
3948 ENDIF
3949 !-- uncomment next line for non-linear root distribution
3950 !cc TRANF(1)=part(1)
3951 DO K=2,NROOT
3952 totliq=soiliqw(k)+qmin
3953 sm1=totliq
3954 sm2=sm1*sm1
3955 sm3=sm2*sm1
3956 sm4=sm3*sm1
3957 gx=ap0+ap1*sm1+ap2*sm2+ap3*sm3+ap4*sm4
3958 if(totliq.ge.ref) gx=1.
3959 if(totliq.le.0.) gx=0.
3960 if(gx.gt.1.) gx=1.
3961 if(gx.lt.0.) gx=0.
3962 DID=zshalf(K+1)-zshalf(K)
3963 part(k)=did*gx
3964 IF(totliq.GE.REF) THEN
3965 TRANF(K)=DID
3966 ELSE IF(totliq.LE.WILT) THEN
3967 TRANF(K)=0.
3968 ELSE
3969 TRANF(K)=(totliq-WILT) &
3970 /(REF-WILT)*DID
3971 ENDIF
3972 !-- uncomment next line for non-linear root distribution
3973 !cc TRANF(k)=part(k)
3974 END DO
3975
3976 !-- TRANSUM - total for the rooting zone
3977 transum=0.
3978 DO K=1,NROOT
3979 transum=transum+tranf(k)
3980 END DO
3981
3982 ! RETURN
3983 ! END
3984 !-----------------------------------------------------------------
3985 END SUBROUTINE TRANSF
3986 !-----------------------------------------------------------------
3987
3988
3989 SUBROUTINE VILKA(TN,D1,D2,PP,QS,TS,TT,NSTEP,ii,j,iland,isoil)
3990 !--------------------------------------------------------------
3991 !--- VILKA finds the solution of energy budget at the surface
3992 !--- using table T,QS computed from Clausius-Klapeiron
3993 !--------------------------------------------------------------
3994 REAL, DIMENSION(1:4001), INTENT(IN ) :: TT
3995 REAL, INTENT(IN ) :: TN,D1,D2,PP
3996 INTEGER, INTENT(IN ) :: NSTEP,ii,j,iland,isoil
3997
3998 REAL, INTENT(OUT ) :: QS, TS
3999
4000 REAL :: F1,T1,T2,RN
4001 INTEGER :: I,I1
4002
4003 I=(TN-1.7315E2)/.05+1
4004 T1=173.1+FLOAT(I)*.05
4005 F1=T1+D1*TT(I)-D2
4006 I1=I-F1/(.05+D1*(TT(I+1)-TT(I)))
4007 I=I1
4008 IF(I.GT.4000.OR.I.LT.1) GOTO 1
4009 10 I1=I
4010 T1=173.1+FLOAT(I)*.05
4011 F1=T1+D1*TT(I)-D2
4012 RN=F1/(.05+D1*(TT(I+1)-TT(I)))
4013 I=I-INT(RN)
4014 IF(I.GT.4000.OR.I.LT.1) GOTO 1
4015 IF(I1.NE.I) GOTO 10
4016 TS=T1-.05*RN
4017 QS=(TT(I)+(TT(I)-TT(I+1))*RN)/PP
4018 GOTO 20
4019 1 PRINT *,' AVOST IN VILKA '
4020 ! WRITE(12,*)'AVOST',TN,D1,D2,PP,NSTEP
4021 PRINT *,TN,D1,D2,PP,NSTEP,I,TT(i),ii,j,iland,isoil
4022 CALL wrf_error_fatal (' AVOST IN VILKA ' )
4023 20 CONTINUE
4024 ! RETURN
4025 ! END
4026 !-----------------------------------------------------------------------
4027 END SUBROUTINE VILKA
4028 !-----------------------------------------------------------------------
4029
4030
4031 SUBROUTINE SOILVEGIN ( IVGTYP,ISLTYP,MYJ, &
4032 IFOREST,EMISS,PC,ZNT,QWRTZ, &
4033 RHOCS,BCLH,DQM,KSAT,PSIS,QMIN,REF,WILT )
4034
4035 !************************************************************************
4036 ! Set-up soil and vegetation Parameters in the case when
4037 ! snow disappears during the forecast and snow parameters
4038 ! shold be replaced by surface parameters according to
4039 ! soil and vegetation types in this point.
4040 !
4041 ! Output:
4042 !
4043 !
4044 ! Soil parameters:
4045 ! DQM: MAX soil moisture content - MIN (m^3/m^3)
4046 ! REF: Reference soil moisture (m^3/m^3)
4047 ! WILT: Wilting PT soil moisture contents (m^3/m^3)
4048 ! QMIN: Air dry soil moist content limits (m^3/m^3)
4049 ! PSIS: SAT soil potential coefs. (m)
4050 ! KSAT: SAT soil diffusivity/conductivity coefs. (m/s)
4051 ! BCLH: Soil diffusivity/conductivity exponent.
4052 !
4053 ! ************************************************************************
4054 IMPLICIT NONE
4055 !---------------------------------------------------------------------------
4056 integer, parameter :: nsoilclas=19
4057 integer, parameter :: nvegclas=24
4058 integer, parameter :: iwater=16
4059 integer, parameter :: ilsnow=99
4060
4061
4062 !--- soiltyp classification according to STATSGO(nclasses=16)
4063 !
4064 ! 1 SAND SAND
4065 ! 2 LOAMY SAND LOAMY SAND
4066 ! 3 SANDY LOAM SANDY LOAM
4067 ! 4 SILT LOAM SILTY LOAM
4068 ! 5 SILT SILTY LOAM
4069 ! 6 LOAM LOAM
4070 ! 7 SANDY CLAY LOAM SANDY CLAY LOAM
4071 ! 8 SILTY CLAY LOAM SILTY CLAY LOAM
4072 ! 9 CLAY LOAM CLAY LOAM
4073 ! 10 SANDY CLAY SANDY CLAY
4074 ! 11 SILTY CLAY SILTY CLAY
4075 ! 12 CLAY LIGHT CLAY
4076 ! 13 ORGANIC MATERIALS LOAM
4077 ! 14 WATER
4078 ! 15 BEDROCK
4079 ! Bedrock is reclassified as class 14
4080 ! 16 OTHER (land-ice)
4081 ! 17 Playa
4082 ! 18 Lava
4083 ! 19 White Sand
4084 !
4085 !----------------------------------------------------------------------
4086 REAL LQMA(nsoilclas),LRHC(nsoilclas), &
4087 LPSI(nsoilclas),LQMI(nsoilclas), &
4088 LBCL(nsoilclas),LKAS(nsoilclas), &
4089 LWIL(nsoilclas),LREF(nsoilclas), &
4090 DATQTZ(nsoilclas)
4091 !-- LQMA Rawls et al.[1982]
4092 ! DATA LQMA /0.417, 0.437, 0.453, 0.501, 0.486, 0.463, 0.398,
4093 ! & 0.471, 0.464, 0.430, 0.479, 0.475, 0.439, 1.0, 0.20, 0.401/
4094 !---
4095 !-- Clapp, R. and G. Hornberger, 1978: Empirical equations for some soil
4096 ! hydraulic properties, Water Resour. Res., 14, 601-604.
4097
4098 !-- Clapp et al. [1978]
4099 DATA LQMA /0.395, 0.410, 0.435, 0.485, 0.485, 0.451, 0.420, &
4100 0.477, 0.476, 0.426, 0.492, 0.482, 0.451, 1.0, &
4101 0.20, 0.435, 0.468, 0.200, 0.339/
4102
4103 !-- LREF Rawls et al.[1982]
4104 ! DATA LREF /0.091, 0.125, 0.207, 0.330, 0.360, 0.270, 0.255,
4105 ! & 0.366, 0.318, 0.339, 0.387, 0.396, 0.329, 1.0, 0.108, 0.283/
4106
4107 !-- Clapp et al. [1978]
4108 DATA LREF /0.174, 0.179, 0.249, 0.369, 0.369, 0.314, 0.299, &
4109 0.357, 0.391, 0.316, 0.409, 0.400, 0.314, 1., &
4110 0.1, 0.249, 0.454, 0.17, 0.236/
4111
4112 !-- LWIL Rawls et al.[1982]
4113 ! DATA LWIL/0.033, 0.055, 0.095, 0.133, 0.133, 0.117, 0.148,
4114 ! & 0.208, 0.197, 0.239, 0.250, 0.272, 0.066, 0.0, 0.006, 0.029/
4115
4116 !-- Clapp et al. [1978]
4117 DATA LWIL/0.068, 0.075, 0.114, 0.179, 0.179, 0.155, 0.175, &
4118 0.218, 0.250, 0.219, 0.283, 0.286, 0.155, 0.0, &
4119 0.006, 0.114, 0.030, 0.006, 0.01/
4120
4121 ! DATA LQMI/0.010, 0.028, 0.047, 0.084, 0.084, 0.066, 0.067,
4122 ! & 0.120, 0.103, 0.100, 0.126, 0.138, 0.066, 0.0, 0.006, 0.028/
4123
4124 !-- Carsel and Parrish [1988]
4125 DATA LQMI/0.045, 0.057, 0.065, 0.067, 0.034, 0.078, 0.10, &
4126 0.089, 0.095, 0.10, 0.070, 0.068, 0.078, 0.0, &
4127 0.004, 0.065, 0.020, 0.004, 0.008/
4128
4129 !-- LPSI Cosby et al[1984]
4130 ! DATA LPSI/0.060, 0.036, 0.141, 0.759, 0.759, 0.355, 0.135,
4131 ! & 0.617, 0.263, 0.098, 0.324, 0.468, 0.355, 0.0, 0.069, 0.036/
4132 ! & 0.617, 0.263, 0.098, 0.324, 0.468, 0.355, 0.0, 0.069, 0.036/
4133
4134 !-- Clapp et al. [1978]
4135 DATA LPSI/0.121, 0.090, 0.218, 0.786, 0.786, 0.478, 0.299, &
4136 0.356, 0.630, 0.153, 0.490, 0.405, 0.478, 0.0, &
4137 0.121, 0.218, 0.468, 0.069, 0.069/
4138
4139 !-- LKAS Rawls et al.[1982]
4140 ! DATA LKAS/5.83E-5, 1.70E-5, 7.19E-6, 1.89E-6, 1.89E-6,
4141 ! & 3.67E-6, 1.19E-6, 4.17E-7, 6.39E-7, 3.33E-7, 2.50E-7,
4142 ! & 1.67E-7, 3.38E-6, 0.0, 1.41E-4, 1.41E-5/
4143
4144 !-- Clapp et al. [1978]
4145 DATA LKAS/1.76E-4, 1.56E-4, 3.47E-5, 7.20E-6, 7.20E-6, &
4146 6.95E-6, 6.30E-6, 1.70E-6, 2.45E-6, 2.17E-6, &
4147 1.03E-6, 1.28E-6, 6.95E-6, 0.0, 1.41E-4, &
4148 3.47E-5, 1.28E-6, 1.41E-4, 1.76E-4/
4149
4150 !-- LBCL Cosby et al [1984]
4151 ! DATA LBCL/2.79, 4.26, 4.74, 5.33, 5.33, 5.25, 6.66,
4152 ! & 8.72, 8.17, 10.73, 10.39, 11.55, 5.25, 0.0, 2.79, 4.26/
4153
4154 !-- Clapp et al. [1978]
4155 DATA LBCL/4.05, 4.38, 4.90, 5.30, 5.30, 5.39, 7.12, &
4156 7.75, 8.52, 10.40, 10.40, 11.40, 5.39, 0.0, &
4157 4.05, 4.90, 11.55, 2.79, 2.79/
4158
4159 DATA LRHC /1.47,1.41,1.34,1.27,1.27,1.21,1.18,1.32,1.23, &
4160 1.18,1.15,1.09,1.21,4.18,2.03,2.10,1.09,2.03,1.47/
4161
4162 DATA DATQTZ/0.92,0.82,0.60,0.25,0.10,0.40,0.60,0.10,0.35, &
4163 0.52,0.10,0.25,0.00,0.,0.60,0.0,0.25,0.60,0.92/
4164
4165 !--------------------------------------------------------------------------
4166 !
4167 ! USGS Vegetation Types
4168 !
4169 ! 1: Urban and Built-Up Land
4170 ! 2: Dryland Cropland and Pasture
4171 ! 3: Irrigated Cropland and Pasture
4172 ! 4: Mixed Dryland/Irrigated Cropland and Pasture
4173 ! 5: Cropland/Grassland Mosaic
4174 ! 6: Cropland/Woodland Mosaic
4175 ! 7: Grassland
4176 ! 8: Shrubland
4177 ! 9: Mixed Shrubland/Grassland
4178 ! 10: Savanna
4179 ! 11: Deciduous Broadleaf Forest
4180 ! 12: Deciduous Needleleaf Forest
4181 ! 13: Evergreen Broadleaf Forest
4182 ! 14: Evergreen Needleleaf Fores
4183 ! 15: Mixed Forest
4184 ! 16: Water Bodies
4185 ! 17: Herbaceous Wetland
4186 ! 18: Wooded Wetland
4187 ! 19: Barren or Sparsely Vegetated
4188 ! 20: Herbaceous Tundra
4189 ! 21: Wooded Tundra
4190 ! 22: Mixed Tundra
4191 ! 23: Bare Ground Tundra
4192 ! 24: Snow or Ice
4193 !
4194 ! 25: Playa
4195 ! 26: Lava
4196 ! 27: White Sand
4197
4198
4199 !---- Below are the arrays for the vegetation parameters
4200 REAL LALB(nvegclas),LMOI(nvegclas),LEMI(nvegclas), &
4201 LROU(nvegclas),LTHI(nvegclas),LSIG(nvegclas), &
4202 LPC(nvegclas), NROTBL(nvegclas)
4203
4204 !************************************************************************
4205 !---- vegetation parameters
4206 !
4207 !-- USGS model
4208 !
4209 DATA LALB/.18,.17,.18,.18,.18,.16,.19,.22,.20,.20,.16,.14, &
4210 .12,.12,.13,.08,.14,.14,.25,.15,.15,.15,.25,.55/
4211 DATA LEMI/.88,4*.92,.93,.92,.88,.9,.92,.93,.94, &
4212 .95,.95,.94,.98,.95,.95,.85,.92,.93,.92,.85,.95/
4213 !-- Roughness length is changed for forests and some others
4214 ! DATA LROU/.5,.06,.075,.065,.05,.2,.075,.1,.11,.15,.8,.85, &
4215 ! 2.0,1.0,.563,.0001,.2,.4,.05,.1,.15,.1,.065,.05/
4216 DATA LROU/.5,.06,.075,.065,.05,.2,.075,.1,.11,.15,.5,.5, &
4217 .5,.5,.5,.0001,.2,.4,.05,.1,.15,.1,.065,.05/
4218
4219 DATA LMOI/.1,.3,.5,.25,.25,.35,.15,.1,.15,.15,.3,.3, &
4220 .5,.3,.3,1.,.6,.35,.02,.5,.5,.5,.02,.95/
4221 !
4222 !---- still needs to be corrected
4223 !
4224 ! DATA LPC/ 15*.8,0.,.8,.8,.5,.5,.5,.5,.5,.0/
4225 DATA LPC /0.4,0.3,0.4,0.4,0.4,0.4,0.4,0.4,0.4,0.4,5*0.55,0.,0.55,0.55, &
4226 0.3,0.3,0.4,0.4,0.3,0./
4227
4228 ! DATA LPC /0.6,6*0.8,0.7,0.75,6*0.8,0.,0.8,0.8, &
4229 ! 0.5,0.7,0.6,0.7,0.5,0./
4230
4231
4232 !***************************************************************************
4233
4234
4235 INTEGER :: &
4236 IVGTYP, &
4237 ISLTYP
4238
4239 LOGICAL, INTENT(IN ) :: myj
4240
4241 REAL , &
4242 INTENT ( OUT) :: pc
4243
4244 REAL , &
4245 INTENT (INOUT ) :: emiss, &
4246 znt
4247 !--- soil properties
4248 REAL , &
4249 INTENT( OUT) :: RHOCS, &
4250 BCLH, &
4251 DQM, &
4252 KSAT, &
4253 PSIS, &
4254 QMIN, &
4255 QWRTZ, &
4256 REF, &
4257 WILT
4258
4259 INTEGER, DIMENSION( 1:(nvegclas+3) ) , &
4260 INTENT ( OUT) :: iforest
4261
4262
4263
4264 INTEGER, DIMENSION( 1:(nvegclas+3) ) :: if1
4265 INTEGER :: kstart, kfin, lstart, lfin
4266 INTEGER :: i,j,k
4267
4268 !***********************************************************************
4269 ! DATA ZS1/0.0,0.05,0.20,0.40,1.6,3.0/ ! o - levels in soil
4270 ! DATA ZS2/0.0,0.025,0.125,0.30,1.,2.3/ ! x - levels in soil
4271 DATA IF1/12*0,1,1,1,12*0/
4272
4273 do k=1,nvegclas+3
4274 iforest(k)=if1(k)
4275 enddo
4276
4277
4278 ! EMISS = LEMI(IVGTYP)
4279 EMISS = LEMITBL(IVGTYP)
4280 ! When MYJ sfc scheme is used - better use recommended in MYJSFCINIT
4281 ! values of roughness length, and not redefine it here.
4282 ! The table in this routine is the one we use in RUC with RUC LSM.
4283
4284 !tgs 3 Oct 2007 - MYJSFCINIT roughness produced unrealistic fluxes,
4285 !tgs - will not use it any more!
4286 !tgs if (.not. myj) then
4287 ! ZNT = LROU(IVGTYP)
4288 ZNT = Z0TBL(IVGTYP)
4289 !tgs endif
4290
4291 ! PC = LPC (IVGTYP)
4292 PC = PCTBL(IVGTYP)
4293 ! RHOCS = LRHC(ISLTYP)*1.E6
4294 RHOCS = HC(ISLTYP)*1.E6
4295
4296 !tgs - parameters from SOILPARM.TBL
4297 BCLH = BB(ISLTYP)
4298 DQM = MAXSMC(ISLTYP)- &
4299 DRYSMC(ISLTYP)
4300 KSAT = SATDK(ISLTYP)
4301 PSIS = - SATPSI(ISLTYP)
4302 QMIN = DRYSMC(ISLTYP)
4303 REF = REFSMC(ISLTYP)
4304 WILT = WLTSMC(ISLTYP)
4305 QWRTZ = QTZ(ISLTYP)
4306
4307 ! BCLH = LBCL(ISLTYP)
4308 ! DQM = LQMA(ISLTYP)- &
4309 ! LQMI(ISLTYP)
4310 ! KSAT = LKAS(ISLTYP)
4311 ! PSIS = - LPSI(ISLTYP)
4312 ! QMIN = LQMI(ISLTYP)
4313 ! REF = LREF(ISLTYP)
4314 ! WILT = LWIL(ISLTYP)
4315 ! QWRTZ = DATQTZ(ISLTYP)
4316
4317 !--------------------------------------------------------------------------
4318 END SUBROUTINE SOILVEGIN
4319 !--------------------------------------------------------------------------
4320
4321
4322 SUBROUTINE SNOWFREE (ivgtyp,myj,emiss,znt,iland)
4323 !************************************************************************
4324 ! Set-up soil and vegetation Parameters in the case when
4325 ! snow disappears during the forecast and snow parameters
4326 ! shold be replaced by surface parameters according to
4327 ! soil and vegetation types in this point.
4328 !
4329 !***************************************************************************
4330 IMPLICIT NONE
4331 !---------------------------------------------------------------------------
4332 integer, parameter :: nvegclas=24
4333
4334
4335 INTEGER :: IVGTYP
4336
4337 LOGICAL, INTENT(IN ) :: myj
4338
4339 REAL, INTENT(INOUT) :: &
4340 emiss, &
4341 znt
4342 INTEGER, INTENT(INOUT) :: ILAND
4343
4344 !---- Below are the arrays for the vegetation parameters
4345 REAL, DIMENSION( 1:nvegclas ) :: LALB, &
4346 LEMI, &
4347 LROU_MYJ,&
4348 LROU
4349
4350 !************************************************************************
4351 !-- USGS model
4352 !
4353 DATA LALB/.18,.17,.18,.18,.18,.16,.19,.22,.20,.20,.16,.14, &
4354 .12,.12,.13,.08,.14,.14,.25,.15,.15,.15,.25,.55/
4355 DATA LEMI/.88,4*.92,.93,.92,.88,.9,.92,.93,.94, &
4356 .95,.95,.94,.98,.95,.95,.85,.92,.93,.92,.85,.95/
4357 !-- Roughness length is changed for forests and some others
4358 ! next 2 lines - table from RUC
4359 ! DATA LROU/.5,.06,.075,.065,.05,.2,.075,.1,.11,.15,.8,.85, &
4360 ! 2.0,1.0,.563,.0001,.2,.4,.05,.1,.15,.1,.065,.05/
4361
4362 DATA LROU/.5,.06,.075,.065,.05,.2,.075,.1,.11,.15,.5,.5, &
4363 .5,.5,.5,.0001,.2,.4,.05,.1,.15,.1,.065,.05/
4364
4365 ! With MYJSFC better use the table from MYJSFCINIT
4366 DATA LROU_MYJ/1.0,.07,.07,.07,.07,.15,.08,.03,.05,.86,.8,.85, &
4367 2.65,1.09,.8,.001,.04,.05,.01,.04,.06,.05,.03,.001/
4368
4369
4370
4371 !--------------------------------------------------------------------------
4372
4373 EMISS = LEMITBL(IVGTYP)
4374 !tgs 3 Oct 2007 - LROU_MYJ gives unrealistic surface fluxes with RUC LSM with RUC LSM
4375 ! if(myj) then
4376 ! ZNT = LROU_MYJ(IVGTYP)
4377 ! else
4378 ZNT = Z0TBL(IVGTYP)
4379 !!! ZNT = LROU(IVGTYP)
4380 ! endif
4381 ILAND = IVGTYP
4382 ! ---
4383
4384 ! RETURN
4385 ! END
4386 !--------------------------------------------------------------------------
4387 END SUBROUTINE SNOWFREE
4388
4389 SUBROUTINE LSMRUCINIT( SMFR3D,TSLB,SMOIS,ISLTYP,mavail, &
4390 nzs, restart, &
4391 allowed_to_read , &
4392 ids,ide, jds,jde, kds,kde, &
4393 ims,ime, jms,jme, kms,kme, &
4394 its,ite, jts,jte, kts,kte )
4395
4396 IMPLICIT NONE
4397
4398
4399 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
4400 ims,ime, jms,jme, kms,kme, &
4401 its,ite, jts,jte, kts,kte, &
4402 nzs
4403
4404 REAL, DIMENSION( ims:ime, 1:nzs, jms:jme ) , &
4405 INTENT(IN) :: TSLB, &
4406 SMOIS
4407
4408 INTEGER, DIMENSION( ims:ime, jms:jme ) , &
4409 INTENT(INOUT) :: ISLTYP
4410
4411 REAL, DIMENSION( ims:ime, 1:nzs, jms:jme ) , &
4412 INTENT(INOUT) :: SMFR3D
4413
4414 REAL, DIMENSION( ims:ime, jms:jme ) , &
4415 INTENT(INOUT) :: MAVAIL
4416
4417 REAL, DIMENSION ( 1:nzs ) :: SOILIQW
4418
4419 LOGICAL , INTENT(IN) :: restart, allowed_to_read
4420
4421 !
4422 INTEGER :: I,J,L,itf,jtf
4423 REAL :: RIW,XLMELT,TLN,DQM,REF,PSIS,QMIN,BCLH
4424
4425 INTEGER :: errflag
4426
4427 ! itf=min0(ite,ide-1)
4428 ! jtf=min0(jte,jde-1)
4429
4430
4431 RIW=900.*1.e-3
4432 XLMELT=3.335E+5
4433
4434 ! initialize three LSM related tables
4435 IF ( allowed_to_read ) THEN
4436 CALL wrf_message( 'INITIALIZE THREE LSM RELATED TABLES' )
4437 CALL RUCLSM_PARM_INIT
4438 ENDIF
4439
4440 IF(.not.restart)THEN
4441
4442 itf=min0(ite,ide-1)
4443 jtf=min0(jte,jde-1)
4444
4445 errflag = 0
4446 DO j = jts,jtf
4447 DO i = its,itf
4448 IF ( ISLTYP( i,j ) .LT. 1 ) THEN
4449 errflag = 1
4450 WRITE(err_message,*)"module_sf_ruclsm.F: lsminit: out of range ISLTYP ",i,j,ISLTYP( i,j )
4451 CALL wrf_message(err_message)
4452 ENDIF
4453 ENDDO
4454 ENDDO
4455 IF ( errflag .EQ. 1 ) THEN
4456 CALL wrf_error_fatal( "module_sf_ruclsm.F: lsminit: out of range value "// &
4457 "of ISLTYP. Is this field in the input?" )
4458 ENDIF
4459
4460 DO J=jts,jtf
4461 DO I=its,itf
4462
4463 ! CALL SOILIN ( ISLTYP(I,J), DQM, REF, PSIS, QMIN, BCLH )
4464
4465
4466 !--- Computation of volumetric content of ice in soil
4467 !--- and initialize MAVAIL
4468 DQM = MAXSMC (ISLTYP(I,J)) - &
4469 DRYSMC (ISLTYP(I,J))
4470 REF = REFSMC (ISLTYP(I,J))
4471 PSIS = - SATPSI (ISLTYP(I,J))
4472 QMIN = DRYSMC (ISLTYP(I,J))
4473 BCLH = BB (ISLTYP(I,J))
4474
4475
4476 !!! IF (.not.restart) THEN
4477
4478 if(isltyp(i,j).ne.14) then
4479 mavail(i,j) = max(0.00001,min(1.,smois(i,1,j)/dqm))
4480 ! mavail(i,j) = max(0.00001,min(1.,smois(i,1,j)/(ref-qmin)))
4481 else
4482 mavail(i,j) = 1.
4483 endif
4484 DO L=1,NZS
4485 if(isltyp(i,j).ne.14) then
4486 !-- for land points initialize soil ice
4487 tln=log(TSLB(i,l,j)/273.15)
4488
4489 if(tln.lt.0.) then
4490 soiliqw(l)=(dqm+qmin)*(XLMELT* &
4491 (tslb(i,l,j)-273.15)/tslb(i,l,j)/9.81/psis) &
4492 **(-1./bclh)-qmin
4493 soiliqw(l)=max(0.,soiliqw(l))
4494 soiliqw(l)=min(soiliqw(l),smois(i,l,j))
4495 smfr3d(i,l,j)=(smois(i,l,j)-soiliqw(l))/RIW
4496
4497 else
4498 smfr3d(i,l,j)=0.
4499 endif
4500 else
4501 !-- for water points
4502 smfr3d(i,l,j)=0.
4503 endif
4504
4505 ENDDO
4506 ! ENDIF
4507
4508 ENDDO
4509 ENDDO
4510
4511 ENDIF
4512
4513 END SUBROUTINE lsmrucinit
4514 !
4515 !-----------------------------------------------------------------
4516 SUBROUTINE RUCLSM_PARM_INIT
4517 !-----------------------------------------------------------------
4518
4519 character*8 :: MMINLU, MMINSL
4520
4521 MMINLU='USGS-RUC'
4522 MMINSL='STAS-RUC'
4523 call RUCLSM_SOILVEGPARM( MMINLU, MMINSL)
4524
4525 !-----------------------------------------------------------------
4526 END SUBROUTINE RUCLSM_PARM_INIT
4527 !-----------------------------------------------------------------
4528
4529 !-----------------------------------------------------------------
4530 SUBROUTINE RUCLSM_SOILVEGPARM( MMINLU, MMINSL)
4531 !-----------------------------------------------------------------
4532
4533 IMPLICIT NONE
4534
4535 integer :: LUMATCH, IINDEX, LC, NUM_SLOPE
4536 integer :: ierr
4537 INTEGER , PARAMETER :: OPEN_OK = 0
4538
4539 character*8 :: MMINLU, MMINSL
4540 character*128 :: mess , message
4541 logical, external :: wrf_dm_on_monitor
4542
4543
4544 !-----SPECIFY VEGETATION RELATED CHARACTERISTICS :
4545 ! ALBBCK: SFC albedo (in percentage)
4546 ! Z0: Roughness length (m)
4547 ! LEMI: Emissivity
4548 ! PC: Plant coefficient for transpiration function
4549 ! -- the rest of the parameters are read in but not used currently
4550 ! SHDFAC: Green vegetation fraction (in percentage)
4551 ! Note: The ALBEDO, Z0, and SHDFAC values read from the following table
4552 ! ALBEDO, amd Z0 are specified in LAND-USE TABLE; and SHDFAC is
4553 ! the monthly green vegetation data
4554 ! CMXTBL: MAX CNPY Capacity (m)
4555 ! NROTBL: Rooting depth (layer)
4556 ! RSMIN: Mimimum stomatal resistance (s m-1)
4557 ! RSMAX: Max. stomatal resistance (s m-1)
4558 ! RGL: Parameters used in radiation stress function
4559 ! HS: Parameter used in vapor pressure deficit functio
4560 ! TOPT: Optimum transpiration air temperature. (K)
4561 ! CMCMAX: Maximum canopy water capacity
4562 ! CFACTR: Parameter used in the canopy inteception calculati
4563 ! SNUP: Threshold snow depth (in water equivalent m) that
4564 ! implies 100% snow cover
4565 ! LAI: Leaf area index (dimensionless)
4566 ! MAXALB: Upper bound on maximum albedo over deep snow
4567 !
4568 !-----READ IN VEGETAION PROPERTIES FROM VEGPARM.TBL
4569 !
4570
4571 IF ( wrf_dm_on_monitor() ) THEN
4572
4573 OPEN(19, FILE='VEGPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
4574 IF(ierr .NE. OPEN_OK ) THEN
4575 WRITE(message,FMT='(A)') &
4576 'module_sf_ruclsm.F: soil_veg_gen_parm: failure opening VEGPARM.TBL'
4577 CALL wrf_error_fatal ( message )
4578 END IF
4579
4580 WRITE ( mess, * ) 'INPUT VEGPARM FOR ',MMINLU
4581 CALL wrf_message( mess )
4582
4583 LUMATCH=0
4584
4585 READ (19,*)
4586 READ (19,2000,END=2002)LUTYPE
4587 READ (19,*)LUCATS,IINDEX
4588 2000 FORMAT (A8)
4589
4590 IF(LUTYPE.NE.MMINLU)THEN
4591 DO LC=1,LUCATS
4592 READ (19,*)IINDEX,ALBTBL(LC),Z0TBL(LC), &
4593 SHDTBL(LC),NROTBL(LC),RSTBL(LC),RGLTBL(LC), &
4594 HSTBL(LC),SNUPTBL(LC),LAITBL(LC),MAXALB(LC)
4595 ENDDO
4596 READ (19,*)
4597 READ (19,*)TOPT_DATA
4598 READ (19,*)
4599 READ (19,*)CMCMAX_DATA
4600 READ (19,*)
4601 READ (19,*)CFACTR_DATA
4602 READ (19,*)
4603 READ (19,*)RSMAX_DATA
4604 READ (19,*)
4605 READ (19,*)BARE
4606 ENDIF
4607
4608 READ (19,*)
4609 READ (19,2000,END=2002)LUTYPE
4610 READ (19,*)LUCATS,IINDEX
4611
4612
4613 IF(LUTYPE.EQ.MMINLU)THEN
4614 WRITE( mess , * ) 'VEGPARM FOR ',LUTYPE,' FOUND', &
4615 LUCATS,' CATEGORIES'
4616 CALL wrf_message( mess )
4617 LUMATCH=1
4618 ENDIF
4619
4620 IF(LUTYPE.EQ.MMINLU)THEN
4621 DO LC=1,LUCATS
4622 READ (19,*)IINDEX,ALBTBL(LC),Z0TBL(LC),LEMITBL(LC),PCTBL(LC), &
4623 SHDTBL(LC),NROTBL(LC),RSTBL(LC),RGLTBL(LC), &
4624 HSTBL(LC),SNUPTBL(LC),LAITBL(LC),MAXALB(LC)
4625 ENDDO
4626 !
4627 READ (19,*)
4628 READ (19,*)TOPT_DATA
4629 READ (19,*)
4630 READ (19,*)CMCMAX_DATA
4631 READ (19,*)
4632 READ (19,*)CFACTR_DATA
4633 READ (19,*)
4634 READ (19,*)RSMAX_DATA
4635 READ (19,*)
4636 READ (19,*)BARE
4637 ENDIF
4638 !
4639 2002 CONTINUE
4640
4641 CLOSE (19)
4642 ENDIF
4643
4644 CALL wrf_dm_bcast_string ( LUTYPE , 8 )
4645 CALL wrf_dm_bcast_integer ( LUCATS , 1 )
4646 CALL wrf_dm_bcast_integer ( IINDEX , 1 )
4647 CALL wrf_dm_bcast_integer ( LUMATCH , 1 )
4648 CALL wrf_dm_bcast_real ( ALBTBL , NLUS )
4649 CALL wrf_dm_bcast_real ( Z0TBL , NLUS )
4650 CALL wrf_dm_bcast_real ( LEMITBL , NLUS )
4651 CALL wrf_dm_bcast_real ( PCTBL , NLUS )
4652 CALL wrf_dm_bcast_real ( SHDTBL , NLUS )
4653 CALL wrf_dm_bcast_real ( NROTBL , NLUS )
4654 CALL wrf_dm_bcast_real ( RSTBL , NLUS )
4655 CALL wrf_dm_bcast_real ( RGLTBL , NLUS )
4656 CALL wrf_dm_bcast_real ( HSTBL , NLUS )
4657 CALL wrf_dm_bcast_real ( SNUPTBL , NLUS )
4658 CALL wrf_dm_bcast_real ( LAITBL , NLUS )
4659 CALL wrf_dm_bcast_real ( MAXALB , NLUS )
4660 CALL wrf_dm_bcast_real ( TOPT_DATA , 1 )
4661 CALL wrf_dm_bcast_real ( CMCMAX_DATA , 1 )
4662 CALL wrf_dm_bcast_real ( CFACTR_DATA , 1 )
4663 CALL wrf_dm_bcast_real ( RSMAX_DATA , 1 )
4664 CALL wrf_dm_bcast_integer ( BARE , 1 )
4665
4666 !
4667 !-----READ IN SOIL PROPERTIES FROM SOILPARM.TBL
4668 !
4669 IF ( wrf_dm_on_monitor() ) THEN
4670 OPEN(19, FILE='SOILPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
4671 IF(ierr .NE. OPEN_OK ) THEN
4672 WRITE(message,FMT='(A)') &
4673 'module_sf_ruclsm.F: soil_veg_gen_parm: failure opening SOILPARM.TBL'
4674 CALL wrf_error_fatal ( message )
4675 END IF
4676
4677 WRITE(mess,*) 'INPUT SOIL TEXTURE CLASSIFICAION = ',MMINSL
4678 CALL wrf_message( mess )
4679
4680 LUMATCH=0
4681
4682 READ (19,*)
4683 READ (19,2000,END=2003)SLTYPE
4684 READ (19,*)SLCATS,IINDEX
4685 IF(SLTYPE.NE.MMINSL)THEN
4686 DO LC=1,SLCATS
4687 READ (19,*) IINDEX,BB(LC),DRYSMC(LC),HC(LC),MAXSMC(LC),&
4688 REFSMC(LC),SATPSI(LC),SATDK(LC), SATDW(LC), &
4689 WLTSMC(LC), QTZ(LC)
4690 ENDDO
4691 ENDIF
4692 READ (19,*)
4693 READ (19,2000,END=2003)SLTYPE
4694 READ (19,*)SLCATS,IINDEX
4695
4696 IF(SLTYPE.EQ.MMINSL)THEN
4697 WRITE( mess , * ) 'SOIL TEXTURE CLASSIFICATION = ',SLTYPE,' FOUND', &
4698 SLCATS,' CATEGORIES'
4699 CALL wrf_message ( mess )
4700 LUMATCH=1
4701 ENDIF
4702 IF(SLTYPE.EQ.MMINSL)THEN
4703 DO LC=1,SLCATS
4704 READ (19,*) IINDEX,BB(LC),DRYSMC(LC),HC(LC),MAXSMC(LC),&
4705 REFSMC(LC),SATPSI(LC),SATDK(LC), SATDW(LC), &
4706 WLTSMC(LC), QTZ(LC)
4707 ENDDO
4708 ENDIF
4709
4710 2003 CONTINUE
4711
4712 CLOSE (19)
4713 ENDIF
4714
4715 CALL wrf_dm_bcast_integer ( LUMATCH , 1 )
4716 CALL wrf_dm_bcast_string ( SLTYPE , 8 )
4717 CALL wrf_dm_bcast_string ( MMINSL , 8 ) ! since this is reset above, see oct2 ^
4718 CALL wrf_dm_bcast_integer ( SLCATS , 1 )
4719 CALL wrf_dm_bcast_integer ( IINDEX , 1 )
4720 CALL wrf_dm_bcast_real ( BB , NSLTYPE )
4721 CALL wrf_dm_bcast_real ( DRYSMC , NSLTYPE )
4722 CALL wrf_dm_bcast_real ( HC , NSLTYPE )
4723 CALL wrf_dm_bcast_real ( MAXSMC , NSLTYPE )
4724 CALL wrf_dm_bcast_real ( REFSMC , NSLTYPE )
4725 CALL wrf_dm_bcast_real ( SATPSI , NSLTYPE )
4726 CALL wrf_dm_bcast_real ( SATDK , NSLTYPE )
4727 CALL wrf_dm_bcast_real ( SATDW , NSLTYPE )
4728 CALL wrf_dm_bcast_real ( WLTSMC , NSLTYPE )
4729 CALL wrf_dm_bcast_real ( QTZ , NSLTYPE )
4730
4731 IF(LUMATCH.EQ.0)THEN
4732 CALL wrf_message( 'SOIl TEXTURE IN INPUT FILE DOES NOT ' )
4733 CALL wrf_message( 'MATCH SOILPARM TABLE' )
4734 CALL wrf_error_fatal ( 'INCONSISTENT OR MISSING SOILPARM FILE' )
4735 ENDIF
4736
4737 !
4738 !-----READ IN GENERAL PARAMETERS FROM GENPARM.TBL
4739 !
4740 IF ( wrf_dm_on_monitor() ) THEN
4741 OPEN(19, FILE='GENPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
4742 IF(ierr .NE. OPEN_OK ) THEN
4743 WRITE(message,FMT='(A)') &
4744 'module_sf_ruclsm.F: soil_veg_gen_parm: failure opening GENPARM.TBL'
4745 CALL wrf_error_fatal ( message )
4746 END IF
4747
4748 READ (19,*)
4749 READ (19,*)
4750 READ (19,*) NUM_SLOPE
4751
4752 SLPCATS=NUM_SLOPE
4753
4754 DO LC=1,SLPCATS
4755 READ (19,*)SLOPE_DATA(LC)
4756 ENDDO
4757
4758 READ (19,*)
4759 READ (19,*)SBETA_DATA
4760 READ (19,*)
4761 READ (19,*)FXEXP_DATA
4762 READ (19,*)
4763 READ (19,*)CSOIL_DATA
4764 READ (19,*)
4765 READ (19,*)SALP_DATA
4766 READ (19,*)
4767 READ (19,*)REFDK_DATA
4768 READ (19,*)
4769 READ (19,*)REFKDT_DATA
4770 READ (19,*)
4771 READ (19,*)FRZK_DATA
4772 READ (19,*)
4773 READ (19,*)ZBOT_DATA
4774 READ (19,*)
4775 READ (19,*)CZIL_DATA
4776 READ (19,*)
4777 READ (19,*)SMLOW_DATA
4778 READ (19,*)
4779 READ (19,*)SMHIGH_DATA
4780 CLOSE (19)
4781 ENDIF
4782
4783 CALL wrf_dm_bcast_integer ( NUM_SLOPE , 1 )
4784 CALL wrf_dm_bcast_integer ( SLPCATS , 1 )
4785 CALL wrf_dm_bcast_real ( SLOPE_DATA , NSLOPE )
4786 CALL wrf_dm_bcast_real ( SBETA_DATA , 1 )
4787 CALL wrf_dm_bcast_real ( FXEXP_DATA , 1 )
4788 CALL wrf_dm_bcast_real ( CSOIL_DATA , 1 )
4789 CALL wrf_dm_bcast_real ( SALP_DATA , 1 )
4790 CALL wrf_dm_bcast_real ( REFDK_DATA , 1 )
4791 CALL wrf_dm_bcast_real ( REFKDT_DATA , 1 )
4792 CALL wrf_dm_bcast_real ( FRZK_DATA , 1 )
4793 CALL wrf_dm_bcast_real ( ZBOT_DATA , 1 )
4794 CALL wrf_dm_bcast_real ( CZIL_DATA , 1 )
4795 CALL wrf_dm_bcast_real ( SMLOW_DATA , 1 )
4796 CALL wrf_dm_bcast_real ( SMHIGH_DATA , 1 )
4797
4798
4799 !-----------------------------------------------------------------
4800 END SUBROUTINE RUCLSM_SOILVEGPARM
4801 !-----------------------------------------------------------------
4802
4803
4804 SUBROUTINE SOILIN (ISLTYP, DQM, REF, PSIS, QMIN, BCLH )
4805
4806 !--- soiltyp classification according to STATSGO(nclasses=16)
4807 !
4808 ! 1 SAND SAND
4809 ! 2 LOAMY SAND LOAMY SAND
4810 ! 3 SANDY LOAM SANDY LOAM
4811 ! 4 SILT LOAM SILTY LOAM
4812 ! 5 SILT SILTY LOAM
4813 ! 6 LOAM LOAM
4814 ! 7 SANDY CLAY LOAM SANDY CLAY LOAM
4815 ! 8 SILTY CLAY LOAM SILTY CLAY LOAM
4816 ! 9 CLAY LOAM CLAY LOAM
4817 ! 10 SANDY CLAY SANDY CLAY
4818 ! 11 SILTY CLAY SILTY CLAY
4819 ! 12 CLAY LIGHT CLAY
4820 ! 13 ORGANIC MATERIALS LOAM
4821 ! 14 WATER
4822 ! 15 BEDROCK
4823 ! Bedrock is reclassified as class 14
4824 ! 16 OTHER (land-ice)
4825 ! extra classes from Fei Chen
4826 ! 17 Playa
4827 ! 18 Lava
4828 ! 19 White Sand
4829 !
4830 !----------------------------------------------------------------------
4831 integer, parameter :: nsoilclas=19
4832
4833 integer, intent ( in) :: isltyp
4834 real, intent ( out) :: dqm,ref,qmin,psis
4835
4836 REAL LQMA(nsoilclas),LREF(nsoilclas),LBCL(nsoilclas), &
4837 LPSI(nsoilclas),LQMI(nsoilclas)
4838
4839 !-- LQMA Rawls et al.[1982]
4840 ! DATA LQMA /0.417, 0.437, 0.453, 0.501, 0.486, 0.463, 0.398,
4841 ! & 0.471, 0.464, 0.430, 0.479, 0.475, 0.439, 1.0, 0.20, 0.401/
4842 !---
4843 !-- Clapp, R. and G. Hornberger, Empirical equations for some soil
4844 ! hydraulic properties, Water Resour. Res., 14,601-604,1978.
4845 !-- Clapp et al. [1978]
4846 DATA LQMA /0.395, 0.410, 0.435, 0.485, 0.485, 0.451, 0.420, &
4847 0.477, 0.476, 0.426, 0.492, 0.482, 0.451, 1.0, &
4848 0.20, 0.435, 0.468, 0.200, 0.339/
4849
4850 !-- Clapp et al. [1978]
4851 DATA LREF /0.174, 0.179, 0.249, 0.369, 0.369, 0.314, 0.299, &
4852 0.357, 0.391, 0.316, 0.409, 0.400, 0.314, 1., &
4853 0.1, 0.249, 0.454, 0.17, 0.236/
4854
4855 !-- Carsel and Parrish [1988]
4856 DATA LQMI/0.045, 0.057, 0.065, 0.067, 0.034, 0.078, 0.10, &
4857 0.089, 0.095, 0.10, 0.070, 0.068, 0.078, 0.0, &
4858 0.004, 0.065, 0.020, 0.004, 0.008/
4859
4860 !-- Clapp et al. [1978]
4861 DATA LPSI/0.121, 0.090, 0.218, 0.786, 0.786, 0.478, 0.299, &
4862 0.356, 0.630, 0.153, 0.490, 0.405, 0.478, 0.0, &
4863 0.121, 0.218, 0.468, 0.069, 0.069/
4864
4865 !-- Clapp et al. [1978]
4866 DATA LBCL/4.05, 4.38, 4.90, 5.30, 5.30, 5.39, 7.12, &
4867 7.75, 8.52, 10.40, 10.40, 11.40, 5.39, 0.0, &
4868 4.05, 4.90, 11.55, 2.79, 2.79/
4869
4870
4871 DQM = LQMA(ISLTYP)- &
4872 LQMI(ISLTYP)
4873 REF = LREF(ISLTYP)
4874 PSIS = - LPSI(ISLTYP)
4875 QMIN = LQMI(ISLTYP)
4876 BCLH = LBCL(ISLTYP)
4877
4878 END SUBROUTINE SOILIN
4879
4880 END MODULE module_sf_ruclsm