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