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