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