| blob | e7e44dec839ee6ae676fdacafe1d7d1b07361529 |
1 #define LSMRUC_DBG_LVL 3000
2 !WRF:MODEL_LAYER:PHYSICS
3 !
4 MODULE module_sf_ruclsm
6 USE module_model_constants
7 USE module_wrf_error
9 ! VEGETATION PARAMETERS
10 INTEGER :: LUCATS , BARE, NATURAL
11 integer, PARAMETER :: NLUS=50
12 CHARACTER*8 LUTYPE
13 INTEGER, DIMENSION(1:NLUS) :: NROTBL
14 real, dimension(1:NLUS) :: SNUPTBL, RSTBL, RGLTBL, HSTBL, LAITBL, &
15 ALBTBL, Z0TBL, LEMITBL, PCTBL, SHDTBL, MAXALB
16 REAL :: TOPT_DATA,CMCMAX_DATA,CFACTR_DATA,RSMAX_DATA
18 ! SOIL PARAMETERS
19 INTEGER :: SLCATS
20 INTEGER, PARAMETER :: NSLTYPE=30
21 CHARACTER*8 SLTYPE
22 REAL, DIMENSION (1:NSLTYPE) :: BB,DRYSMC,HC, &
23 MAXSMC, REFSMC,SATPSI,SATDK,SATDW, WLTSMC,QTZ
25 ! LSM GENERAL PARAMETERS
26 INTEGER :: SLPCATS
27 INTEGER, PARAMETER :: NSLOPE=30
28 REAL, DIMENSION (1:NSLOPE) :: SLOPE_DATA
29 REAL :: SBETA_DATA,FXEXP_DATA,CSOIL_DATA,SALP_DATA,REFDK_DATA, &
30 REFKDT_DATA,FRZK_DATA,ZBOT_DATA, SMLOW_DATA,SMHIGH_DATA, &
31 CZIL_DATA
33 CHARACTER*256 :: err_message
36 CONTAINS
38 !-----------------------------------------------------------------
39 SUBROUTINE LSMRUC( &
40 DT,KTAU,NSL,ZS, &
41 RAINBL,SNOW,SNOWH,SNOWC,FRZFRAC,frpcpn, &
42 Z3D,P8W,T3D,QV3D,QC3D,RHO3D, & !p8W in [PA]
43 GLW,GSW,EMISS,CHKLOWQ, CHS, &
44 FLQC,FLHC,MAVAIL,CANWAT,VEGFRA,ALB,ZNT, &
45 Z0,SNOALB,ALBBCK, & !new
46 QSFC,QSG,QVG,QCG,DEW,SOILT1,TSNAV, &
47 TBOT,IVGTYP,ISLTYP,XLAND, &
48 ISICE,XICE,XICE_THRESHOLD, &
49 CP,ROVCP,G0,LV,STBOLT, &
50 SOILMOIS,SH2O,SMAVAIL,SMMAX, &
51 TSO,SOILT,HFX,QFX,LH, &
52 SFCRUNOFF,UDRUNOFF,SFCEXC, &
53 SFCEVP,GRDFLX,ACSNOW,SNOM, &
54 SMFR3D,KEEPFR3DFLAG, &
55 myj, &
56 ids,ide, jds,jde, kds,kde, &
57 ims,ime, jms,jme, kms,kme, &
58 its,ite, jts,jte, kts,kte )
59 !-----------------------------------------------------------------
60 IMPLICIT NONE
61 !-----------------------------------------------------------------
62 !
63 ! The RUC LSM model is described in:
64 ! Smirnova, T.G., J.M. Brown, and S.G. Benjamin, 1997:
65 ! Performance of different soil model configurations in simulating
66 ! ground surface temperature and surface fluxes.
67 ! Mon. Wea. Rev. 125, 1870-1884.
68 ! Smirnova, T.G., J.M. Brown, and D. Kim, 2000: Parameterization of
69 ! cold-season processes in the MAPS land-surface scheme.
70 ! J. Geophys. Res. 105, 4077-4086.
71 !-----------------------------------------------------------------
72 !-- DT time step (second)
73 ! ktau - number of time step
74 ! NSL - number of soil layers
75 ! NZS - number of levels in soil
76 ! ZS - depth of soil levels (m)
77 !-- RAINBL - accumulated rain in [mm] between the PBL calls
78 !-- RAINNCV one time step grid scale precipitation (mm/step)
79 ! SNOW - snow water equivalent [mm]
80 ! FRAZFRAC - fraction of frozen precipitation
81 !-- SNOWC flag indicating snow coverage (1 for snow cover)
82 !-- Z3D heights (m)
83 !-- P8W 3D pressure (Pa)
84 !-- T3D temperature (K)
85 !-- QV3D 3D water vapor mixing ratio (Kg/Kg)
86 ! QC3D - 3D cloud water mixing ratio (Kg/Kg)
87 ! RHO3D - 3D air density (kg/m^3)
88 !-- GLW downward long wave flux at ground surface (W/m^2)
89 !-- GSW absorbed short wave flux at ground surface (W/m^2)
90 !-- EMISS surface emissivity (between 0 and 1)
91 ! FLQC - surface exchange coefficient for moisture (kg/m^2/s)
92 ! FLHC - surface exchange coefficient for heat [W/m^2/s/degreeK]
93 ! SFCEXC - surface exchange coefficient for heat [m/s]
94 ! CANWAT - CANOPY MOISTURE CONTENT (mm)
95 ! VEGFRA - vegetation fraction (between 0 and 1)
96 ! ALB - surface albedo (between 0 and 1)
97 ! SNOALB - maximum snow albedo (between 0 and 1)
98 ! ALBBCK - snow-free albedo (between 0 and 1)
99 ! ZNT - roughness length [m]
100 !-- TBOT soil temperature at lower boundary (K)
101 ! IVGTYP - USGS vegetation type (24 classes)
102 ! ISLTYP - STASGO soil type (16 classes)
103 !-- XLAND land mask (1 for land, 2 for water)
104 !-- CP heat capacity at constant pressure for dry air (J/kg/K)
105 !-- G0 acceleration due to gravity (m/s^2)
106 !-- LV latent heat of melting (J/kg)
107 !-- STBOLT Stefan-Boltzmann constant (W/m^2/K^4)
108 ! SOILMOIS - soil moisture content (volumetric fraction)
109 ! TSO - soil temp (K)
110 !-- SOILT surface temperature (K)
111 !-- HFX upward heat flux at the surface (W/m^2)
112 !-- QFX upward moisture flux at the surface (kg/m^2/s)
113 !-- LH upward latent heat flux (W/m^2)
114 ! SFCRUNOFF - ground surface runoff [mm]
115 ! UDRUNOFF - underground runoff [mm]
116 ! SFCEVP - total evaporation in [kg/m^2]
117 ! GRDFLX - soil heat flux (W/m^2: negative, if downward from surface)
118 ! ACSNOW - accumulation of snow water [m]
119 !-- CHKLOWQ - is either 0 or 1 (so far set equal to 1).
120 !-- used only in MYJPBL.
121 !-- tice - sea ice temperture (C)
122 !-- rhosice - sea ice density (kg m^-3)
123 !-- capice - sea ice volumetric heat capacity (J/m^3/K)
124 !-- thdifice - sea ice thermal diffusivity (m^2/s)
125 !--
126 !-- ims start index for i in memory
127 !-- ime end index for i in memory
128 !-- jms start index for j in memory
129 !-- jme end index for j in memory
130 !-- kms start index for k in memory
131 !-- kme end index for k in memory
132 !-------------------------------------------------------------------------
133 ! INTEGER, PARAMETER :: nzss=5
134 ! INTEGER, PARAMETER :: nddzs=2*(nzss-2)
136 INTEGER, PARAMETER :: nvegclas=24+3
138 REAL, INTENT(IN ) :: DT
139 LOGICAL, INTENT(IN ) :: myj,frpcpn
140 INTEGER, INTENT(IN ) :: ktau, nsl, isice, &
141 ims,ime, jms,jme, kms,kme, &
142 ids,ide, jds,jde, kds,kde, &
143 its,ite, jts,jte, kts,kte
145 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
146 INTENT(IN ) :: QV3D, &
147 QC3D, &
148 p8w, &
149 rho3D, &
150 T3D, &
151 z3D
153 REAL, DIMENSION( ims:ime , jms:jme ), &
154 INTENT(IN ) :: RAINBL, &
155 GLW, &
156 GSW, &
157 FLHC, &
158 FLQC, &
159 CHS , &
160 EMISS, &
161 XICE, &
162 XLAND, &
163 ALBBCK, &
164 VEGFRA, &
165 TBOT
167 REAL, DIMENSION( 1:nsl), INTENT(IN ) :: ZS
169 REAL, DIMENSION( ims:ime , jms:jme ), &
170 INTENT(INOUT) :: &
171 SNOW, & !new
172 SNOWH, &
173 SNOWC, &
174 CANWAT, & ! new
175 SNOALB, &
176 ALB, &
177 MAVAIL, &
178 SFCEXC, &
179 Z0 , &
180 ZNT
182 REAL, DIMENSION( ims:ime , jms:jme ), &
183 INTENT(IN ) :: &
184 FRZFRAC
186 INTEGER, DIMENSION( ims:ime , jms:jme ), &
187 INTENT(IN ) :: IVGTYP, &
188 ISLTYP
190 REAL, INTENT(IN ) :: CP,ROVCP,G0,LV,STBOLT,XICE_threshold
192 REAL, DIMENSION( ims:ime , 1:nsl, jms:jme ) , &
193 INTENT(INOUT) :: SOILMOIS,SH2O,TSO
195 REAL, DIMENSION( ims:ime, jms:jme ) , &
196 INTENT(INOUT) :: SOILT, &
197 HFX, &
198 QFX, &
199 LH, &
200 SFCEVP, &
201 SFCRUNOFF, &
202 UDRUNOFF, &
203 GRDFLX, &
204 ACSNOW, &
205 SNOM, &
206 QVG, &
207 QCG, &
208 DEW, &
209 QSFC, &
210 QSG, &
211 CHKLOWQ, &
212 SOILT1, &
213 TSNAV
215 REAL, DIMENSION( ims:ime, jms:jme ) , &
216 INTENT(INOUT) :: SMAVAIL, &
217 SMMAX
219 REAL, DIMENSION( its:ite, jts:jte ) :: &
220 PC, &
221 RUNOFF1, &
222 RUNOFF2, &
223 EMISSL, &
224 ZNTL, &
225 LMAVAIL, &
226 SMELT, &
227 SNOH, &
228 SNFLX, &
229 EDIR, &
230 EC, &
231 ETT, &
232 SUBLIM, &
233 sflx, &
234 EVAPL, &
235 PRCPL, &
236 SEAICE, &
237 INFILTR
239 !--- soil/snow properties
240 REAL, DIMENSION( ims:ime, 1:nsl, jms:jme) &
241 :: KEEPFR3DFLAG, &
242 SMFR3D
244 REAL &
245 :: RHOCS, &
246 RHOSN, &
247 RHONEWSN, &
248 BCLH, &
249 DQM, &
250 KSAT, &
251 PSIS, &
252 QMIN, &
253 QWRTZ, &
254 REF, &
255 WILT, &
256 CANWATR, &
257 SNOWFRAC, &
258 SNHEI, &
259 SNWE
261 REAL :: CN, &
262 SAT,CW, &
263 C1SN, &
264 C2SN, &
265 KQWRTZ, &
266 KICE, &
267 KWT
270 REAL, DIMENSION(1:NSL) :: ZSMAIN, &
271 ZSHALF, &
272 DTDZS2
274 REAL, DIMENSION(1:2*(nsl-2)) :: DTDZS
276 REAL, DIMENSION(1:4001) :: TBQ
279 REAL, DIMENSION( 1:nsl ) :: SOILM1D, &
280 TSO1D, &
281 SOILICE, &
282 SOILIQW, &
283 SMFRKEEP
285 REAL, DIMENSION( 1:nsl ) :: KEEPFR
288 REAL :: RSM, &
289 SNWEPRINT, &
290 SNHEIPRINT
292 REAL :: PRCPMS, &
293 NEWSNMS, &
294 PATM, &
295 PATMB, &
296 TABS, &
297 QVATM, &
298 QCATM, &
299 Q2SAT, &
300 CONFLX, &
301 RHO, &
302 QKMS, &
303 TKMS, &
304 INFILTRP
305 REAL :: cq,r61,r273,arp,brp,x,evs,eis
307 REAL :: meltfactor
308 INTEGER :: NROOT
309 INTEGER :: ILAND,ISOIL
311 INTEGER, DIMENSION( 1:(nvegclas) ) :: IFOREST
313 INTEGER :: I,J,K,NZS,NZS1,NDDZS
314 INTEGER :: k1,l,k2,kp,km
315 CHARACTER (LEN=132) :: message
317 !-----------------------------------------------------------------
319 NZS=NSL
320 NDDZS=2*(nzs-2)
322 !---- table TBQ is for resolution of balance equation in VILKA
323 CQ=173.15-.05
324 R273=1./273.15
325 R61=6.1153*0.62198
326 ARP=77455.*41.9/461.525
327 BRP=64.*41.9/461.525
329 DO K=1,4001
330 CQ=CQ+.05
331 ! TBQ(K)=R61*EXP(ARP*(R273-1./CQ)-BRP*LOG(CQ*R273))
332 EVS=EXP(17.67*(CQ-273.15)/(CQ-29.65))
333 EIS=EXP(22.514-6.15E3/CQ)
334 if(CQ.ge.273.15) then
335 ! tbq is in mb
336 tbq(k) = R61*evs
337 else
338 tbq(k) = R61*eis
339 endif
341 END DO
343 !--- Initialize soil/vegetation parameters
344 !--- This is temporary until SI is added to mass coordinate ---!!!!!
346 #if ( NMM_CORE == 1 )
347 if(ktau+1.eq.1) then
348 #else
349 if(ktau.eq.1) then
350 #endif
351 DO J=jts,jte
352 DO i=its,ite
353 do k=1,nsl
354 ! smfr3d (i,k,j)=soilmois(i,k,j)/900.*1.e3
355 ! sh2o (i,k,j)=soilmois(i,k,j)-smfr3d(i,k,j)/1.e3*900.
356 keepfr3dflag(i,k,j)=0.
357 enddo
358 !--- initializing to zero snow fraction
359 snowc(i,j) = min(1.,snowh(i,j)/0.05)
360 !--- initializing inside snow temp if it is not defined
361 IF((soilt1(i,j) .LT. 170.) .or. (soilt1(i,j) .GT.400.)) THEN
362 IF(snowc(i,j).gt.0.1) THEN
363 soilt1(i,j)=0.5*(soilt(i,j)+tso(i,1,j))
364 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
365 WRITE ( message , FMT='(A,F8.3,2I6)' ) &
366 'Temperature inside snow is initialized in RUCLSM ', soilt1(i,j),i,j
367 CALL wrf_debug ( 0 , message )
368 ENDIF
369 ELSE
370 soilt1(i,j) = soilt(i,j)
371 ENDIF
372 ENDIF
373 tsnav(i,j) =0.5*(soilt(i,j)+tso(i,1,j))-273.
374 qcg (i,j) =0.
375 patmb=P8w(i,kms,j)*1.e-2
376 QSG (i,j) = QSN(SOILT(i,j),TBQ)/PATMB
377 qvg (i,j) = QSG(i,j)*mavail(i,j)
378 ! qvg (i,j) =qv3d(i,1,j)
379 ! qsfc(i,j) = qsg(i,j)/(1.+qsg(i,j))
380 qsfc(i,j) = qvg(i,j)/(1.+qvg(i,j))
381 SMELT(i,j) = 0.
382 SNOM (i,j) = 0.
383 SNFLX(i,j) = 0.
384 DEW (i,j) = 0.
385 PC (i,j) = 0.
386 zntl (i,j) = 0.
387 RUNOFF1(i,j) = 0.
388 RUNOFF2(i,j) = 0.
389 SFCRUNOFF(i,j) = 0.
390 UDRUNOFF(i,j) = 0.
391 emissl (i,j) = 0.
392 ! Temporarily!!!
393 canwat(i,j)=0.
395 ! For RUC LSM CHKLOWQ needed for MYJPBL should
396 ! 1 because is actual specific humidity at the surface, and
397 ! not the saturation value
398 chklowq(i,j) = 1.
399 infiltr(i,j) = 0.
400 snoh (i,j) = 0.
401 edir (i,j) = 0.
402 ec (i,j) = 0.
403 ett (i,j) = 0.
404 sublim(i,j) = 0.
405 sflx (i,j) = 0.
406 evapl (i,j) = 0.
407 prcpl (i,j) = 0.
408 ENDDO
409 ENDDO
411 do k=1,nsl
412 soilice(k)=0.
413 soiliqw(k)=0.
414 enddo
415 endif
417 !-----------------------------------------------------------------
419 PRCPMS = 0.
422 DO J=jts,jte
424 DO i=its,ite
426 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
427 print *,' IN LSMRUC ','ims,ime,jms,jme,its,ite,jts,jte,nzs', &
428 ims,ime,jms,jme,its,ite,jts,jte,nzs
429 print *,' IVGTYP, ISLTYP ', ivgtyp(i,j),isltyp(i,j)
430 print *,' MAVAIL ', mavail(i,j)
431 print *,' SOILT,QVG,P8w',soilt(i,j),qvg(i,j),p8w(i,1,j)
432 print *, 'LSMRUC, I,J,xland, QFX,HFX from SFCLAY',i,j,xland(i,j), &
433 qfx(i,j),hfx(i,j)
434 print *, ' GSW, GLW =',gsw(i,j),glw(i,j)
435 print *, 'SOILT, TSO start of time step =',soilt(i,j),(tso(i,k,j),k=1,nsl)
436 print *, 'SOILMOIS start of time step =',(soilmois(i,k,j),k=1,nsl)
437 print *, 'SMFROZEN start of time step =',(smfr3d(i,k,j),k=1,nsl)
438 print *, ' I,J=, after SFCLAY CHS,FLHC ',i,j,chs(i,j),flhc(i,j)
439 print *, 'LSMRUC, IVGTYP,ISLTYP,ZNT,ALB = ', ivgtyp(i,j),isltyp(i,j),znt(i,j),alb(i,j),i,j
440 print *, 'LSMRUC I,J,DT,RAINBL =',I,J,dt,RAINBL(i,j)
441 print *, 'XLAND ---->, ivgtype,isoiltyp,i,j',xland(i,j),ivgtyp(i,j),isltyp(i,j),i,j
442 ENDIF
445 ILAND = IVGTYP(i,j)
446 ISOIL = ISLTYP(I,J)
447 TABS = T3D(i,kms,j)
448 QVATM = QV3D(i,kms,j)
449 QCATM = QC3D(i,kms,j)
450 PATM = P8w(i,kms,j)*1.e-5
451 !-- Z3D(1) is thickness between first full sigma level and the surface,
452 !-- but first mass level is at the half of the first sigma level
453 !-- (u and v are also at the half of first sigma level)
454 CONFLX = Z3D(i,kms,j)*0.5
455 RHO = RHO3D(I,kms,J)
456 !--- 1*e-3 is to convert from mm/s to m/s
457 IF(FRPCPN) THEN
458 PRCPMS = (RAINBL(i,j)/DT*1.e-3)*(1-FRZFRAC(I,J))
459 NEWSNMS = (RAINBL(i,j)/DT*1.e-3)*FRZFRAC(I,J)
460 ELSE
461 if (tabs.le.273.15) then
462 PRCPMS = 0.
463 NEWSNMS = RAINBL(i,j)/DT*1.e-3
464 else
465 PRCPMS = RAINBL(i,j)/DT*1.e-3
466 NEWSNMS = 0.
467 endif
468 ENDIF
469 ! if (myj) then
470 QKMS=CHS(i,j)
471 TKMS=CHS(i,j)
472 ! else
473 !--- convert exchange coeff QKMS to [m/s]
474 ! QKMS=FLQC(I,J)/RHO/MAVAIL(I,J)
475 ! TKMS=FLHC(I,J)/RHO/CP
476 ! endif
477 !--- convert incoming snow and canwat from mm to m
478 SNWE=SNOW(I,J)*1.E-3
479 SNHEI=SNOWH(I,J)
480 CANWATR=CANWAT(I,J)*1.E-3
481 SNOWFRAC=SNOWC(I,J)
483 !-----
484 zsmain(1)=0.
485 zshalf(1)=0.
486 do k=2,nzs
487 zsmain(k)= zs(k)
488 zshalf(k)=0.5*(zsmain(k-1) + zsmain(k))
489 enddo
492 !------------------------------------------------------------
493 !----- DDZS and DSDZ1 are for implicit solution of soil eqns.
494 !-------------------------------------------------------------
495 NZS1=NZS-1
496 !-----
497 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
498 print *,' DT,NZS1, ZSMAIN, ZSHALF --->', dt,nzs1,zsmain,zshalf
499 ENDIF
501 DO K=2,NZS1
502 K1=2*K-3
503 K2=K1+1
504 X=DT/2./(ZSHALF(K+1)-ZSHALF(K))
505 DTDZS(K1)=X/(ZSMAIN(K)-ZSMAIN(K-1))
506 DTDZS2(K-1)=X
507 DTDZS(K2)=X/(ZSMAIN(K+1)-ZSMAIN(K))
508 END DO
511 CN=0.5 ! exponent
512 SAT=0.0004 ! canopy water saturated
514 CW =4.183E6
517 !--- Constants used in Johansen soil thermal
518 !--- conductivity method
520 KQWRTZ=7.7
521 KICE=2.2
522 KWT=0.57
524 !***********************************************************************
525 !--- Constants for snow density calculations C1SN and C2SN
527 c1sn=0.026
528 ! c1sn=0.01
529 c2sn=21.
531 !***********************************************************************
533 NROOT= 4
534 ! ! rooting depth
536 RHONEWSN = 200.
537 if(SNOWH(i,j).gt.0.) then
538 RHOSN = SNOW(i,j)/SNOWH(i,j)
539 else
540 RHOSN = 300.
541 endif
543 !--- initializing soil and surface properties
544 CALL SOILVEGIN ( ILAND,ISOIL,MYJ,IFOREST, &
545 EMISSL(I,J),PC(i,j),ZNT(I,J),QWRTZ, &
546 ! EMISSL(I,J),PC(i,j),ZNTL(I,J),QWRTZ, &
547 RHOCS,BCLH,DQM,KSAT,PSIS,QMIN,REF,WILT )
549 !-- definition of number of soil levels in the rooting zone
550 IF(iforest(ivgtyp(i,j)).ne.1) THEN
551 !---- all vegetation types except evergreen and mixed forests
552 !18apr08 - define meltfactor for Egglston melting limit:
553 ! for open areas factor is 2, and for forests - factor is 1.5
554 ! This will make limit on snow melting smaller and let snow stay
555 ! longer in the forests.
556 meltfactor = 2.0
558 do k=2,nzs
559 if(zsmain(k).ge.0.4) then
560 NROOT=K
561 goto 111
562 endif
563 enddo
564 ELSE
565 !---- evergreen and mixed forests
566 !18apr08 - define meltfactor
567 meltfactor = 1.5
569 do k=2,nzs
570 if(zsmain(k).ge.1.1) then
571 NROOT=K
572 goto 111
573 endif
574 enddo
575 ENDIF
576 111 continue
578 !-----
579 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
580 print *,' ZS, ZSMAIN, ZSHALF, CONFLX --->', zs,zsmain,zshalf,conflx
581 print *,'NROOT, iforest, ivgtyp, i,j ', nroot,iforest(ivgtyp(i,j)),ivgtyp(I,J),I,J
582 ENDIF
584 !*** SET ZERO-VALUE FOR SOME OUTPUT DIAGNOSTIC ARRAYS
586 IF((XLAND(I,J)-1.5).GE.0.)THEN
587 !-- Water
588 SMAVAIL(I,J)=1.0
589 SMMAX(I,J)=1.0
590 SNOW(I,J)=0.0
591 LMAVAIL(I,J)=1.0
593 ILAND=16
594 ISOIL=14
596 patmb=P8w(i,1,j)*1.e-2
597 qvg (i,j) = QSN(SOILT(i,j),TBQ)/PATMB
598 qsfc(i,j) = qvg(i,j)/(1.+qvg(i,j))
599 CHKLOWQ(I,J)=1.
600 Q2SAT=QSN(TABS,TBQ)/PATMB
602 DO K=1,NZS
603 SOILMOIS(I,K,J)=1.0
604 SH2O (I,K,J)=1.0
605 TSO(I,K,J)= SOILT(I,J)
606 ENDDO
608 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
609 PRINT*,' water point, I=',I, &
610 'J=',J, 'SOILT=', SOILT(i,j)
611 ENDIF
613 ELSE
615 ! LAND POINT OR SEA ICE
616 if(xice(i,j).ge.xice_threshold) then
617 ! if(IVGTYP(i,j).eq.isice) then
618 SEAICE(i,j)=1.
619 else
620 SEAICE(i,j)=0.
621 endif
623 IF(SEAICE(I,J).GT.0.5)THEN
624 !-- Sea-ice case
625 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
626 PRINT*,' sea-ice at water point, I=',I, &
627 'J=',J
628 ENDIF
629 ILAND = 24
630 ISOIL = 16
631 ZNT(I,J) = 0.011
632 snoalb(i,j) = 0.75
633 dqm = 1.
634 ref = 1.
635 qmin = 0.
636 wilt = 0.
637 emissl(i,j) = 1.0
639 DO K=1,NZS
640 soilmois(i,k,j) = 1.
641 smfr3d(i,k,j) = 1.
642 sh2o(i,k,j) = 0.
643 keepfr3dflag(i,k,j) = 0.
644 ENDDO
645 ENDIF
647 ! Attention!!!! RUC LSM uses soil moisture content minus residual (minimum
648 ! or dry soil moisture content for a given soil type) as a state variable.
650 DO k=1,nzs
651 ! soilm1d - soil moisture content minus residual [m**3/m**3]
652 soilm1d (k) = min(max(0.,soilmois(i,k,j)),dqm)
653 tso1d (k) = tso(i,k,j)
654 soiliqw (k) = sh2o(i,k,j)
655 ENDDO
657 do k=1,nzs
658 smfrkeep(k) = smfr3d(i,k,j)
659 keepfr (k) = keepfr3dflag(i,k,j)
660 enddo
662 ! LMAVAIL(I,J)=max(0.00001,min(1.,soilmois(i,1,j)/(REF-QMIN)))
663 LMAVAIL(I,J)=max(0.00001,min(1.,soilmois(i,1,j)/dqm))
665 ! extract dew from the cloud water at the surface
666 QCG(I,J)=QCG(I,J)-DEW(I,J)
668 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
669 print *,'LAND, i,j,tso1d,soilm1d,PATM,TABS,QVATM,QCATM,RHO', &
670 i,j,tso1d,soilm1d,PATM,TABS,QVATM,QCATM,RHO
671 print *,'CONFLX =',CONFLX
672 print *,'SMFRKEEP,KEEPFR ',SMFRKEEP,KEEPFR
673 ENDIF
675 !-----------------------------------------------------------------
676 CALL SFCTMP (dt,ktau,conflx,i,j, &
677 !--- input variables
678 nzs,nddzs,nroot,meltfactor, & !added meltfactor
679 iland,isoil,xland(i,j),ivgtyp(i,j),PRCPMS, &
680 NEWSNMS,SNWE,SNHEI,SNOWFRAC,RHOSN,RHONEWSN, &
681 PATM,TABS,QVATM,QCATM,RHO, &
682 GLW(I,J),GSW(I,J),EMISSL(I,J), &
683 QKMS,TKMS,PC(I,J),LMAVAIL(I,J), &
684 canwatr,vegfra(I,J),alb(I,J),znt(I,J), &
685 snoalb(i,j),albbck(i,j), & !new
686 myj,seaice(i,j), &
687 !--- soil fixed fields
688 QWRTZ, &
689 rhocs,dqm,qmin,ref, &
690 wilt,psis,bclh,ksat, &
691 sat,cn,zsmain,zshalf,DTDZS,DTDZS2,tbq, &
692 !--- constants
693 cp,rovcp,g0,lv,stbolt,cw,c1sn,c2sn, &
694 KQWRTZ,KICE,KWT, &
695 !--- output variables
696 snweprint,snheiprint,rsm, &
697 soilm1d,tso1d,smfrkeep,keepfr, &
698 soilt(I,J),soilt1(i,j),tsnav(i,j),dew(I,J), &
699 qvg(I,J),qsg(I,J),qcg(I,J),SMELT(I,J), &
700 SNOH(I,J),SNFLX(I,J),SNOM(I,J),ACSNOW(I,J), &
701 edir(I,J),ec(I,J),ett(I,J),qfx(I,J), &
702 lh(I,J),hfx(I,J),sflx(I,J),sublim(I,J), &
703 evapl(I,J),prcpl(I,J),runoff1(I,J), &
704 runoff2(I,J),soilice,soiliqw,infiltrp)
705 !-----------------------------------------------------------------
707 !*** DIAGNOSTICS
708 !--- available and maximum soil moisture content in the soil
709 !--- domain
711 smavail(i,j) = 0.
712 smmax (i,j) = 0.
714 do k=1,nzs-1
715 smavail(i,j)=smavail(i,j)+(qmin+soilm1d(k))* &
716 (zshalf(k+1)-zshalf(k))
717 smmax (i,j) =smmax (i,j)+(qmin+dqm)* &
718 (zshalf(k+1)-zshalf(k))
719 enddo
721 smavail(i,j)=smavail(i,j)+(qmin+soilm1d(nzs))* &
722 (zsmain(nzs)-zshalf(nzs))
723 smmax (i,j) =smmax (i,j)+(qmin+dqm)* &
724 (zsmain(nzs)-zshalf(nzs))
726 !--- Convert the water unit into mm
727 SFCRUNOFF(I,J) = SFCRUNOFF(I,J)+RUNOFF1(I,J)*DT*1000.0
728 UDRUNOFF (I,J) = UDRUNOFF(I,J)+RUNOFF2(I,J)*1000.0
729 SMAVAIL (I,J) = SMAVAIL(I,J) * 1000.
730 SMMAX (I,J) = SMMAX(I,J) * 1000.
732 do k=1,nzs
734 soilmois(i,k,j) = soilm1d(k)
735 sh2o (i,k,j) = soiliqw(k)
736 tso(i,k,j) = tso1d(k)
737 enddo
739 tso(i,nzs,j) = tbot(i,j)
741 do k=1,nzs
742 smfr3d(i,k,j) = smfrkeep(k)
743 keepfr3dflag(i,k,j) = keepfr (k)
744 enddo
746 !tgs add together dew and cloud at the ground surface
747 qcg(i,j)=qcg(i,j)+dew(i,j)
749 Z0 (I,J) = ZNT (I,J)
750 SFCEXC (I,J) = TKMS
751 patmb=P8w(i,1,j)*1.e-2
752 Q2SAT=QSN(TABS,TBQ)/PATMB
753 QSFC(I,J) = QVG(I,J)/(1.+QVG(I,J))
754 ! for MYJ surface and PBL scheme
755 ! if (myj) then
756 ! MYJSFC expects QSFC as actual specific humidity at the surface
757 IF((QVATM.GE.Q2SAT*0.95).AND.QVATM.LT.qvg(I,J))THEN
758 CHKLOWQ(I,J)=0.
759 ELSE
760 CHKLOWQ(I,J)=1.
761 ENDIF
762 ! else
763 ! CHKLOWQ(I,J)=1.
764 ! endif
766 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
767 if(CHKLOWQ(I,J).eq.0.) then
768 print *,'i,j,CHKLOWQ', &
769 i,j,CHKLOWQ(I,J)
770 endif
771 ENDIF
773 ! SNOW is in [mm], SNWE is in [m]; CANWAT is in mm, CANWATR is in m
774 SNOW (i,j) = SNWE*1000.
775 SNOWH (I,J) = SNHEI
776 CANWAT (I,J) = CANWATR*1000.
777 INFILTR(I,J) = INFILTRP
779 MAVAIL (i,j) = LMAVAIL(I,J)
780 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
781 print *,' LAND, I=,J=, QFX, HFX after SFCTMP', i,j,lh(i,j),hfx(i,j)
782 ENDIF
783 !!! QFX (I,J) = LH(I,J)/LV
784 SFCEVP (I,J) = SFCEVP (I,J) + QFX (I,J) * DT
785 GRDFLX (I,J) = -1. * sflx(I,J)
787 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
788 print *,' QFX after change, LH ', i,j, QFX(i,j),LH(I,J)
789 ENDIF
790 !--- SNOWC snow cover flag
791 if(snowfrac > 0. .and. (xice(i,j).ge.xice_threshold .and. xice(i,j) < 1.)) then
792 SNOWFRAC = SNOWFRAC*XICE(I,J)
793 endif
795 SNOWC(I,J)=SNOWFRAC
798 !--- get 3d soil fields
799 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
800 print *,'LAND, i,j,tso1d,soilm1d - end of time step', &
801 i,j,tso1d,soilm1d
802 ENDIF
804 !--- end of a land or sea ice point
805 ENDIF
807 ENDDO
809 ENDDO
811 !-----------------------------------------------------------------
812 END SUBROUTINE LSMRUC
813 !-----------------------------------------------------------------
817 SUBROUTINE SFCTMP (delt,ktau,conflx,i,j, &
818 !--- input variables
819 nzs,nddzs,nroot,meltfactor, &
820 ILAND,ISOIL,XLAND,IVGTYP,PRCPMS, &
821 NEWSNMS,SNWE,SNHEI,SNOWFRAC,RHOSN,RHONEWSN, &
822 PATM,TABS,QVATM,QCATM,rho, &
823 GLW,GSW,EMISS,QKMS,TKMS,PC, &
824 MAVAIL,CST,VEGFRA,ALB,ZNT, &
825 ALB_SNOW,ALB_SNOW_FREE, &
826 MYJ,SEAICE, &
827 !--- soil fixed fields
828 QWRTZ,rhocs,dqm,qmin,ref,wilt,psis,bclh,ksat, &
829 sat,cn,zsmain,zshalf,DTDZS,DTDZS2,tbq, &
830 !--- constants
831 cp,rovcp,g0,lv,stbolt,cw,c1sn,c2sn, &
832 KQWRTZ,KICE,KWT, &
833 !--- output variables
834 snweprint,snheiprint,rsm, &
835 soilm1d,ts1d,smfrkeep,keepfr,soilt,soilt1, &
836 tsnav,dew,qvg,qsg,qcg, &
837 SMELT,SNOH,SNFLX,SNOM,ACSNOW, &
838 edir1,ec1,ett1,eeta,qfx,hfx,s,sublim, &
839 evapl,prcpl,runoff1,runoff2,soilice, &
840 soiliqw,infiltr)
841 !-----------------------------------------------------------------
842 IMPLICIT NONE
843 !-----------------------------------------------------------------
845 !--- input variables
847 INTEGER, INTENT(IN ) :: i,j,nroot,ktau,nzs , &
848 nddzs !nddzs=2*(nzs-2)
850 REAL, INTENT(IN ) :: DELT,CONFLX,meltfactor
851 REAL, INTENT(IN ) :: C1SN,C2SN
852 LOGICAL, INTENT(IN ) :: myj
853 !--- 3-D Atmospheric variables
854 REAL , &
855 INTENT(IN ) :: PATM, &
856 TABS, &
857 QVATM, &
858 QCATM
859 REAL , &
860 INTENT(IN ) :: GLW, &
861 GSW, &
862 PC, &
863 VEGFRA, &
864 ALB_SNOW_FREE, &
865 SEAICE, &
866 XLAND, &
867 RHO, &
868 QKMS, &
869 TKMS
871 INTEGER, INTENT(IN ) :: IVGTYP
872 !--- 2-D variables
873 REAL , &
874 INTENT(INOUT) :: EMISS, &
875 MAVAIL, &
876 SNOWFRAC, &
877 ALB_SNOW, &
878 ALB, &
879 CST
881 !--- soil properties
882 REAL :: &
883 RHOCS, &
884 BCLH, &
885 DQM, &
886 KSAT, &
887 PSIS, &
888 QMIN, &
889 QWRTZ, &
890 REF, &
891 SAT, &
892 WILT
894 REAL, INTENT(IN ) :: CN, &
895 CW, &
896 CP, &
897 ROVCP, &
898 G0, &
899 LV, &
900 STBOLT, &
901 KQWRTZ, &
902 KICE, &
903 KWT
905 REAL, DIMENSION(1:NZS), INTENT(IN) :: ZSMAIN, &
906 ZSHALF, &
907 DTDZS2
910 REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS
912 REAL, DIMENSION(1:4001), INTENT(IN) :: TBQ
915 !--- input/output variables
916 !-------- 3-d soil moisture and temperature
917 REAL, DIMENSION( 1:nzs ) , &
918 INTENT(INOUT) :: TS1D, &
919 SOILM1D, &
920 SMFRKEEP
921 REAL, DIMENSION( 1:nzs ) , &
922 INTENT(INOUT) :: KEEPFR
925 INTEGER, INTENT(INOUT) :: ILAND,ISOIL
927 !-------- 2-d variables
928 REAL , &
929 INTENT(INOUT) :: DEW, &
930 EDIR1, &
931 EC1, &
932 ETT1, &
933 EETA, &
934 EVAPL, &
935 INFILTR, &
936 RHOSN, &
937 RHONEWSN, &
938 SUBLIM, &
939 PRCPL, &
940 QVG, &
941 QSG, &
942 QCG, &
943 QFX, &
944 HFX, &
945 S, &
946 RUNOFF1, &
947 RUNOFF2, &
948 ACSNOW, &
949 SNWE, &
950 SNHEI, &
951 SMELT, &
952 SNOM, &
953 SNOH, &
954 SNFLX, &
955 SOILT, &
956 SOILT1, &
957 TSNAV, &
958 ZNT
960 REAL, DIMENSION(1:NZS) :: &
961 tice, &
962 rhosice, &
963 capice, &
964 thdifice
965 !-------- 1-d variables
966 REAL, DIMENSION(1:NZS), INTENT(OUT) :: SOILICE, &
967 SOILIQW
971 REAL, INTENT(OUT) :: RSM, &
972 SNWEPRINT, &
973 SNHEIPRINT
974 !--- Local variables
976 INTEGER :: K,ILNB
978 REAL :: BSN, XSN , &
979 RAINF, SNTH, NEWSN, PRCPMS, NEWSNMS , &
980 T3, UPFLUX, XINET
981 REAL :: snhei_crit, keep_snow_albedo
983 REAL :: RNET,GSWNEW,EMISSN,ZNTSN
984 REAL :: VEGFRAC
985 real :: cice, albice, albsn
987 !-----------------------------------------------------------------
988 integer, parameter :: ilsnow=99
990 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
991 print *,' in SFCTMP',i,j,nzs,nddzs,nroot, &
992 SNWE,RHOSN,SNOM,SMELT,TS1D
993 ENDIF
995 NEWSN=0.
996 RAINF = 0.
997 RSM=0.
998 INFILTR=0.
999 VEGFRAC=0.01*VEGFRA
1000 if(VEGFRAC.le.0.01) VEGFRAC=0.
1001 !---initialize local arrays for sea ice
1002 do k=1,nzs
1003 tice(k) = 0.
1004 rhosice(k) = 0.
1005 cice = 0.
1006 capice(k) = 0.
1007 thdifice(k) = 0.
1008 enddo
1010 GSWnew=GSW
1011 ALBice=ALB_SNOW_FREE
1012 !--- sea ice properties
1013 !--- N.N Zubov "Arctic Ice"
1014 !--- no salinity dependence because we consider the ice pack
1015 !--- to be old and to have low salinity (0.0002)
1016 if(SEAICE.ge.0.5) then
1017 do k=1,nzs
1018 tice(k) = ts1d(k) - 273.15
1019 rhosice(k) = 917.6/(1-0.000165*tice(k))
1020 cice = 2115.85 +7.7948*tice(k)
1021 capice(k) = cice*rhosice(k)
1022 thdifice(k) = 2.260872/capice(k)
1023 enddo
1024 !-- SEA ICE ALB dependence on ice temperature. When ice temperature is
1025 !-- below critical value of -10C - no change to albedo.
1026 !-- If temperature is higher that -10C then albedo is decreasing.
1027 !-- The minimum albedo at t=0C for ice is 0.1 less.
1028 GSWNEW=GSW/(1.-ALB)
1029 ALBice = MIN(ALB_SNOW_FREE,MAX(ALB_SNOW_FREE - 0.1, &
1030 ALB_SNOW_FREE - 0.1*(tice(1)+10.)/10. ))
1031 GSWNEW=GSW*(1.-ALBice)
1032 endif
1034 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
1035 print *,'I,J,KTAU,QKMS,TKMS', i,j,ktau,qkms,tkms
1036 print *,'GSW,GSWnew,GLW,SOILT,EMISS,ALB,ALBice,SNWE',&
1037 GSW,GSWnew,GLW,SOILT,EMISS,ALB,ALBice,SNWE
1038 ENDIF
1040 SNHEI = SNWE * 1000. / RHOSN
1041 !--------------
1042 T3 = STBOLT*SOILT*SOILT*SOILT
1043 UPFLUX = T3 *SOILT
1044 XINET = EMISS*(GLW-UPFLUX)
1045 RNET = GSWnew + XINET
1047 !Calculate the amount (m) of fresh snow
1049 if(snhei.gt.0.0081*1.e3/rhosn) then
1050 !*** Correct snow density for current temperature (Koren et al. 1999)
1051 BSN=delt/3600.*c1sn*exp(0.08*tsnav-c2sn*rhosn*1.e-3)
1052 if(bsn*snwe*100..lt.1.e-4) goto 777
1053 XSN=rhosn*(exp(bsn*snwe*100.)-1.)/(bsn*snwe*100.)
1054 rhosn=MIN(MAX(100.,XSN),400.)
1055 ! rhosn=MIN(MAX(50.,XSN),400.)
1056 777 continue
1058 else
1059 rhosn =200.
1060 rhonewsn =200.
1061 endif
1063 newsn=newsnms*delt
1064 !---- ACSNOW - accumulation of snow water [m]
1065 acsnow=acsnow+newsn
1067 IF(NEWSN.GE.1.E-8) THEN
1068 !*** Calculate fresh snow density (t > -15C, else MIN value)
1069 !*** Eq. 10 from Koren et al. (1999)
1071 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
1072 print *, 'THERE IS NEW SNOW, newsn', newsn
1073 ENDIF
1074 if(tabs.lt.258.15) then
1075 ! rhonewsn=50.
1076 rhonewsn=100.
1078 else
1079 rhonewsn=1.e3*max((0.10+0.0017*(Tabs-273.15+15.)**1.5) &
1080 , 0.10)
1081 ! rhonewsn=1.e3*max((0.05+0.0017*(Tabs-273.15+15.)**1.5) &
1082 ! , 0.05)
1083 rhonewsn=MIN(rhonewsn,400.)
1084 ! rhonewsn=100.
1085 endif
1088 !*** Define average snow density of the snow pack considering
1089 !*** the amount of fresh snow (eq. 9 in Koren et al.(1999)
1090 !*** without snow melt )
1091 xsn=(rhosn*snwe+rhonewsn*newsn)/ &
1092 (snwe+newsn)
1093 rhosn=MIN(MAX(100.,XSN),400.)
1094 ! rhosn=MIN(MAX(50.,XSN),400.)
1096 snwe=snwe+newsn
1097 snhei=snwe*1.E3/rhosn
1098 NEWSN=NEWSN*1.E3/rhosn
1099 endif
1101 IF(PRCPMS.NE.0.) THEN
1103 ! PRCPMS is liquid precipitation rate
1104 ! RAINF is a flag used for calculation of rain water
1105 ! heat content contribution into heat budget equation. Rain's temperature
1106 ! is set equal to air temperature at the first atmospheric
1107 ! level.
1109 RAINF=1.
1110 ENDIF
1112 IF(SNHEI.GT.0.0) THEN
1113 !--- Set of surface parameters should be changed to snow values for grid
1114 !--- points where the snow cover exceeds snow threshold of 2 cm
1116 SNHEI_CRIT=0.01601*1.e3/rhosn
1117 SNOWFRAC=MIN(1.,SNHEI/(2.*SNHEI_CRIT))
1118 !tgs NEW - low limit on snow fraction
1119 if(SNOWFRAC.lt.0.01) snowfrac=0.
1120 !--- EMISS = 0.91 for snow
1121 EMISS = EMISS*(1.-snowfrac)+0.91*snowfrac
1123 KEEP_SNOW_ALBEDO = 0.
1124 IF (NEWSN.GT.0.) KEEP_SNOW_ALBEDO = 1.
1126 !--- GSW in-coming solar
1127 GSWNEW=GSW/(1.-ALB)
1129 IF(SEAICE .LT. 0.5) THEN
1130 !----- SNOW on soil
1131 !-- ALB dependence on snow depth
1132 ALBsn = MAX(keep_snow_albedo*alb_snow, &
1133 MIN((alb_snow_free + &
1134 (alb_snow - alb_snow_free) * &
1135 (snhei/(2.*SNHEI_CRIT))), alb_snow))
1137 !-- ALB dependence on snow temperature. When snow temperature is
1138 !-- below critical value of -10C - no change to albedo.
1139 !-- If temperature is higher that -10C then albedo is decreasing.
1140 !-- The minimum albedo at t=0C for snow on land is 15% less than
1141 !-- albedo of temperatures below -10C.
1142 if(albsn.lt.alb_snow)then
1143 ALB=ALBsn
1144 else
1145 ALB = MIN(ALBSN,MAX(ALBSN - 0.15*(soilt - 263.15)/ &
1146 (273.15-263.15), ALBSN - 0.15))
1147 endif
1148 ELSE
1149 !----- SNOW on ice
1150 ALBsn = MAX(keep_snow_albedo*alb_snow, &
1151 MIN((albice + (alb_snow - albice) * &
1152 (snhei/(2.*SNHEI_CRIT))), alb_snow))
1154 !-- ALB dependence on snow temperature. When snow temperature is
1155 !-- below critical value of -10C - no change to albedo.
1156 !-- If temperature is higher that -10C then albedo is decreasing.
1157 !-- The minimum albedo at t=0C for snow on ice is 0.15 less.
1158 if(albsn.lt.alb_snow)then
1159 ALB=ALBsn
1160 else
1161 ALB = MIN(ALBSN,MAX(ALBSN - 0.15*(soilt - 263.15)/ &
1162 (273.15-263.15), ALBSN - 0.15))
1163 endif
1165 ENDIF
1167 !--- recompute absorbed solar radiation and net radiation
1168 !--- for new value of albedo
1169 gswnew=gswnew*(1.-alb)
1171 XINET = EMISS*(GLW-UPFLUX)
1172 RNET = GSWnew + XINET
1174 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
1175 print *,'I,J,GSW,GSWnew,GLW,UPFLUX,ALB',&
1176 i,j,GSW,GSWnew,GLW,UPFLUX,ALB
1177 ENDIF
1180 if (SEAICE .LT. 0.5) then
1181 ! LAND
1182 CALL SNOWSOIL ( & !--- input variables
1183 i,j,isoil,delt,ktau,conflx,nzs,nddzs,nroot, &
1184 meltfactor,rhonewsn, & ! new
1185 ILAND,PRCPMS,RAINF,NEWSN,snhei,SNWE,snowfrac, &
1186 RHOSN,PATM,QVATM,QCATM, &
1187 GLW,GSWnew,EMISS,RNET,IVGTYP, &
1188 QKMS,TKMS,PC,CST, &
1189 RHO,VEGFRAC,ALB,ZNT, &
1190 MYJ, &
1191 !--- soil fixed fields
1192 QWRTZ,rhocs,dqm,qmin,ref,wilt,psis,bclh,ksat, &
1193 sat,cn,zsmain,zshalf,DTDZS,DTDZS2,tbq, &
1194 !--- constants
1195 lv,CP,rovcp,G0,cw,stbolt,tabs, &
1196 KQWRTZ,KICE,KWT, &
1197 !--- output variables
1198 ilnb,snweprint,snheiprint,rsm, &
1199 soilm1d,ts1d,smfrkeep,keepfr, &
1200 dew,soilt,soilt1,tsnav,qvg,qsg,qcg, &
1201 SMELT,SNOH,SNFLX,SNOM,edir1,ec1,ett1,eeta, &
1202 qfx,hfx,s,sublim,prcpl,runoff1,runoff2, &
1203 mavail,soilice,soiliqw,infiltr )
1204 else
1205 ! SEA ICE
1206 CALL SNOWSEAICE ( &
1207 i,j,isoil,delt,ktau,conflx,nzs,nddzs, &
1208 meltfactor,rhonewsn, & ! new
1209 ILAND,PRCPMS,RAINF,NEWSN,snhei,SNWE,snowfrac, &
1210 RHOSN,PATM,QVATM,QCATM, &
1211 GLW,GSWnew,EMISS,RNET, &
1212 QKMS,TKMS,RHO, &
1213 !--- sea ice parameters
1214 ALB,ZNT, &
1215 tice,rhosice,capice,thdifice, &
1216 zsmain,zshalf,DTDZS,DTDZS2,tbq, &
1217 !--- constants
1218 lv,CP,rovcp,cw,stbolt,tabs, &
1219 !--- output variables
1220 ilnb,snweprint,snheiprint,rsm,ts1d, &
1221 dew,soilt,soilt1,tsnav,qvg,qsg,qcg, &
1222 SMELT,SNOH,SNFLX,SNOM,eeta, &
1223 qfx,hfx,s,sublim,prcpl &
1224 )
1225 edir1 = eeta
1226 ec1 = 0.
1227 ett1 = 0.
1228 runoff1 = smelt
1229 runoff2 = 0.
1230 mavail = 1.
1231 infiltr=0.
1232 cst=0.
1233 do k=1,nzs
1234 soilm1d(k)=1.
1235 soiliqw(k)=0.
1236 soilice(k)=1.
1237 smfrkeep(k)=1.
1238 keepfr(k)=0.
1239 enddo
1240 endif
1242 if(snhei.eq.0.) then
1243 !--- all snow is melted
1244 alb=alb_snow_free
1245 endif
1247 ELSE
1248 !--- no snow
1249 snheiprint=0.
1250 snweprint=0.
1252 if(SEAICE .LT. 0.5) then
1253 ! LAND
1254 CALL SOIL( &
1255 !--- input variables
1256 i,j,iland,isoil,delt,ktau,conflx,nzs,nddzs,nroot, &
1257 PRCPMS,RAINF,PATM,QVATM,QCATM,GLW,GSWnew, &
1258 EMISS,RNET,QKMS,TKMS,PC,cst,rho,vegfrac, &
1259 !--- soil fixed fields
1260 QWRTZ,rhocs,dqm,qmin,ref,wilt, &
1261 psis,bclh,ksat,sat,cn, &
1262 zsmain,zshalf,DTDZS,DTDZS2,tbq, &
1263 !--- constants
1264 lv,CP,rovcp,G0,cw,stbolt,tabs, &
1265 KQWRTZ,KICE,KWT, &
1266 !--- output variables
1267 soilm1d,ts1d,smfrkeep,keepfr, &
1268 dew,soilt,qvg,qsg,qcg,edir1,ec1, &
1269 ett1,eeta,qfx,hfx,s,evapl,prcpl,runoff1, &
1270 runoff2,mavail,soilice,soiliqw, &
1271 infiltr)
1272 else
1273 ! SEA ICE
1274 CALL SICE( &
1275 !--- input variables
1276 i,j,iland,isoil,delt,ktau,conflx,nzs,nddzs,nroot, &
1277 PRCPMS,RAINF,PATM,QVATM,QCATM,GLW,GSWnew, &
1278 EMISS,RNET,QKMS,TKMS,rho, &
1279 !--- sea ice parameters
1280 tice,rhosice,capice,thdifice, &
1281 zsmain,zshalf,DTDZS,DTDZS2,tbq, &
1282 !--- constants
1283 lv,CP,rovcp,cw,stbolt,tabs, &
1284 !--- output variables
1285 ts1d,dew,soilt,qvg,qsg,qcg, &
1286 eeta,qfx,hfx,s,evapl,prcpl &
1287 )
1288 edir1 = eeta
1289 ec1 = 0.
1290 ett1 = 0.
1291 runoff1 = prcpms
1292 runoff2 = 0.
1293 mavail = 1.
1294 infiltr=0.
1295 cst=0.
1296 do k=1,nzs
1297 soilm1d(k)=1.
1298 soiliqw(k)=0.
1299 soilice(k)=1.
1300 smfrkeep(k)=1.
1301 keepfr(k)=0.
1302 enddo
1303 endif
1305 ENDIF
1306 ! ENDIF
1308 !
1310 ! RETURN
1311 ! END
1312 !---------------------------------------------------------------
1313 END SUBROUTINE SFCTMP
1314 !---------------------------------------------------------------
1317 FUNCTION QSN(TN,T)
1318 !****************************************************************
1319 REAL, DIMENSION(1:4001), INTENT(IN ) :: T
1320 REAL, INTENT(IN ) :: TN
1322 REAL QSN, R,R1,R2
1323 INTEGER I
1325 R=(TN-173.15)/.05+1.
1326 I=INT(R)
1327 IF(I.GE.1) goto 10
1328 I=1
1329 R=1.
1330 10 IF(I.LE.4000) GOTO 20
1331 I=4000
1332 R=4001.
1333 20 R1=T(I)
1334 R2=R-I
1335 QSN=(T(I+1)-R1)*R2 + R1
1336 ! print *,' in QSN, I,R,R1,R2,T(I+1),TN, QSN', I,R,r1,r2,t(i+1),tn,QSN
1337 ! RETURN
1338 ! END
1339 !-----------------------------------------------------------------------
1340 END FUNCTION QSN
1341 !------------------------------------------------------------------------
1344 SUBROUTINE SOIL ( &
1345 !--- input variables
1346 i,j,iland,isoil,delt,ktau,conflx,nzs,nddzs,nroot,&
1347 PRCPMS,RAINF,PATM,QVATM,QCATM, &
1348 GLW,GSW,EMISS,RNET, &
1349 QKMS,TKMS,PC,cst,rho,vegfrac, &
1350 !--- soil fixed fields
1351 QWRTZ,rhocs,dqm,qmin,ref,wilt,psis,bclh,ksat, &
1352 sat,cn,zsmain,zshalf,DTDZS,DTDZS2,tbq, &
1353 !--- constants
1354 xlv,CP,rovcp,G0_P,cw,stbolt,TABS, &
1355 KQWRTZ,KICE,KWT, &
1356 !--- output variables
1357 soilmois,tso,smfrkeep,keepfr, &
1358 dew,soilt,qvg,qsg,qcg, &
1359 edir1,ec1,ett1,eeta,qfx,hfx,s,evapl, &
1360 prcpl,runoff1,runoff2,mavail,soilice, &
1361 soiliqw,infiltrp)
1363 !*************************************************************
1364 ! Energy and moisture budget for vegetated surfaces
1365 ! without snow, heat diffusion and Richards eqns. in
1366 ! soil
1367 !
1368 ! DELT - time step (s)
1369 ! ktau - numver of time step
1370 ! CONFLX - depth of constant flux layer (m)
1371 ! J,I - the location of grid point
1372 ! IME, JME, KME, NZS - dimensions of the domain
1373 ! NROOT - number of levels within the root zone
1374 ! PRCPMS - precipitation rate in m/s
1375 ! PATM - pressure [bar]
1376 ! QVATM,QCATM - cloud and water vapor mixing ratio (kg/kg)
1377 ! at the first atm. level
1378 ! GLW, GSW - incoming longwave and absorbed shortwave
1379 ! radiation at the surface (W/m^2)
1380 ! EMISS,RNET - emissivity of the ground surface (0-1) and net
1381 ! radiation at the surface (W/m^2)
1382 ! QKMS - exchange coefficient for water vapor in the
1383 ! surface layer (m/s)
1384 ! TKMS - exchange coefficient for heat in the surface
1385 ! layer (m/s)
1386 ! PC - plant coefficient (resistance) (0-1)
1387 ! RHO - density of atmosphere near sueface (kg/m^3)
1388 ! VEGFRAC - greeness fraction
1389 ! RHOCS - volumetric heat capacity of dry soil
1390 ! DQM, QMIN - porosity minus residual soil moisture QMIN (m^3/m^3)
1391 ! REF, WILT - field capacity soil moisture and the
1392 ! wilting point (m^3/m^3)
1393 ! PSIS - matrix potential at saturation (m)
1394 ! BCLH - exponent for Clapp-Hornberger parameterization
1395 ! KSAT - saturated hydraulic conductivity (m/s)
1396 ! SAT - maximum value of water intercepted by canopy (m)
1397 ! CN - exponent for calculation of canopy water
1398 ! ZSMAIN - main levels in soil (m)
1399 ! ZSHALF - middle of the soil layers (m)
1400 ! DTDZS,DTDZS2 - dt/(2.*dzshalf*dzmain) and dt/dzshalf in soil
1401 ! TBQ - table to define saturated mixing ration
1402 ! of water vapor for given temperature and pressure
1403 ! SOILMOIS,TSO - soil moisture (m^3/m^3) and temperature (K)
1404 ! DEW - dew in kg/m^2s
1405 ! SOILT - skin temperature (K)
1406 ! QSG,QVG,QCG - saturated mixing ratio, mixing ratio of
1407 ! water vapor and cloud at the ground
1408 ! surface, respectively (kg/kg)
1409 ! EDIR1, EC1, ETT1, EETA - direct evaporation, evaporation of
1410 ! canopy water, transpiration in kg m-2 s-1 and total
1411 ! evaporation in m s-1.
1412 ! QFX, HFX - latent and sensible heat fluxes (W/m^2)
1413 ! S - soil heat flux in the top layer (W/m^2)
1414 ! RUNOFF - surface runoff (m/s)
1415 ! RUNOFF2 - underground runoff (m)
1416 ! MAVAIL - moisture availability in the top soil layer (0-1)
1417 ! INFILTRP - infiltration flux from the top of soil domain (m/s)
1418 !
1419 !*****************************************************************
1420 IMPLICIT NONE
1421 !-----------------------------------------------------------------
1423 !--- input variables
1425 INTEGER, INTENT(IN ) :: nroot,ktau,nzs , &
1426 nddzs !nddzs=2*(nzs-2)
1427 INTEGER, INTENT(IN ) :: i,j,iland,isoil
1428 REAL, INTENT(IN ) :: DELT,CONFLX
1429 !--- 3-D Atmospheric variables
1430 REAL, &
1431 INTENT(IN ) :: PATM, &
1432 QVATM, &
1433 QCATM
1434 !--- 2-D variables
1435 REAL, &
1436 INTENT(IN ) :: GLW, &
1437 GSW, &
1438 EMISS, &
1439 RHO, &
1440 PC, &
1441 VEGFRAC, &
1442 QKMS, &
1443 TKMS
1445 !--- soil properties
1446 REAL, &
1447 INTENT(IN ) :: RHOCS, &
1448 BCLH, &
1449 DQM, &
1450 KSAT, &
1451 PSIS, &
1452 QMIN, &
1453 QWRTZ, &
1454 REF, &
1455 WILT
1457 REAL, INTENT(IN ) :: CN, &
1458 CW, &
1459 KQWRTZ, &
1460 KICE, &
1461 KWT, &
1462 XLV, &
1463 g0_p
1466 REAL, DIMENSION(1:NZS), INTENT(IN) :: ZSMAIN, &
1467 ZSHALF, &
1468 DTDZS2
1470 REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS
1472 REAL, DIMENSION(1:4001), INTENT(IN) :: TBQ
1475 !--- input/output variables
1476 !-------- 3-d soil moisture and temperature
1477 REAL, DIMENSION( 1:nzs ) , &
1478 INTENT(INOUT) :: TSO, &
1479 SOILMOIS, &
1480 SMFRKEEP
1482 REAL, DIMENSION( 1:nzs ) , &
1483 INTENT(INOUT) :: KEEPFR
1485 !-------- 2-d variables
1486 REAL, &
1487 INTENT(INOUT) :: DEW, &
1488 CST, &
1489 EDIR1, &
1490 EC1, &
1491 ETT1, &
1492 EETA, &
1493 EVAPL, &
1494 PRCPL, &
1495 MAVAIL, &
1496 QVG, &
1497 QSG, &
1498 QCG, &
1499 RNET, &
1500 QFX, &
1501 HFX, &
1502 S, &
1503 SAT, &
1504 RUNOFF1, &
1505 RUNOFF2, &
1506 SOILT
1508 !-------- 1-d variables
1509 REAL, DIMENSION(1:NZS), INTENT(OUT) :: SOILICE, &
1510 SOILIQW
1512 !--- Local variables
1514 REAL :: INFILTRP, transum , &
1515 RAINF, PRCPMS , &
1516 TABS, T3, UPFLUX, XINET
1517 REAL :: CP,rovcp,G0,LV,STBOLT,xlmelt,dzstop , &
1518 can,epot,fac,fltot,ft,fq,hft , &
1519 q1,ras,rhoice,sph , &
1520 trans,zn,ci,cvw,tln,tavln,pi , &
1521 DD1,CMC2MS,DRYCAN,WETCAN , &
1522 INFMAX,RIW
1523 REAL, DIMENSION(1:NZS) :: transp,cap,diffu,hydro , &
1524 thdif,tranf,tav,soilmoism , &
1525 soilicem,soiliqwm,detal , &
1526 fwsat,lwsat,told,smold
1528 REAL :: drip
1530 INTEGER :: nzs1,nzs2,k
1532 !-----------------------------------------------------------------
1534 !-- define constants
1535 ! STBOLT=5.670151E-8
1536 RHOICE=900.
1537 CI=RHOICE*2100.
1538 XLMELT=3.35E+5
1539 cvw=cw
1541 SAT=0.0004
1542 prcpl=prcpms
1544 !--- Initializing local arrays
1545 DO K=1,NZS
1546 TRANSP (K)=0.
1547 soilmoism(k)=0.
1548 soilice (k)=0.
1549 soiliqw (k)=0.
1550 soilicem (k)=0.
1551 soiliqwm (k)=0.
1552 lwsat (k)=0.
1553 fwsat (k)=0.
1554 tav (k)=0.
1555 cap (k)=0.
1556 thdif (k)=0.
1557 diffu (k)=0.
1558 hydro (k)=0.
1559 tranf (k)=0.
1560 detal (k)=0.
1561 told (k)=0.
1562 smold (k)=0.
1563 ENDDO
1565 NZS1=NZS-1
1566 NZS2=NZS-2
1567 dzstop=1./(zsmain(2)-zsmain(1))
1568 RAS=RHO*1.E-3
1569 RIW=rhoice*1.e-3
1571 !--- Computation of volumetric content of ice in soil
1573 DO K=1,NZS
1574 !- main levels
1575 tln=log(tso(k)/273.15)
1576 if(tln.lt.0.) then
1577 soiliqw(k)=(dqm+qmin)*(XLMELT* &
1578 (tso(k)-273.15)/tso(k)/9.81/psis) &
1579 **(-1./bclh)-qmin
1580 soiliqw(k)=max(0.,soiliqw(k))
1581 soiliqw(k)=min(soiliqw(k),soilmois(k))
1582 soilice(k)=(soilmois(k)-soiliqw(k))/RIW
1584 !---- melting and freezing is balanced, soil ice cannot increase
1585 if(keepfr(k).eq.1.) then
1586 soilice(k)=min(soilice(k),smfrkeep(k))
1587 soiliqw(k)=max(0.,soilmois(k)-soilice(k)*riw)
1588 endif
1590 else
1591 soilice(k)=0.
1592 soiliqw(k)=soilmois(k)
1593 endif
1595 ENDDO
1597 DO K=1,NZS1
1598 !- middle of soil layers
1599 tav(k)=0.5*(tso(k)+tso(k+1))
1600 soilmoism(k)=0.5*(soilmois(k)+soilmois(k+1))
1601 tavln=log(tav(k)/273.15)
1603 if(tavln.lt.0.) then
1604 soiliqwm(k)=(dqm+qmin)*(XLMELT* &
1605 (tav(k)-273.15)/tav(k)/9.81/psis) &
1606 **(-1./bclh)-qmin
1607 fwsat(k)=dqm-soiliqwm(k)
1608 lwsat(k)=soiliqwm(k)+qmin
1609 soiliqwm(k)=max(0.,soiliqwm(k))
1610 soiliqwm(k)=min(soiliqwm(k), soilmoism(k))
1611 soilicem(k)=(soilmoism(k)-soiliqwm(k))/riw
1612 !---- melting and freezing is balanced, soil ice cannot increase
1613 if(keepfr(k).eq.1.) then
1614 soilicem(k)=min(soilicem(k), &
1615 0.5*(smfrkeep(k)+smfrkeep(k+1)))
1616 soiliqwm(k)=max(0.,soilmoism(k)-soilicem(k)*riw)
1617 fwsat(k)=dqm-soiliqwm(k)
1618 lwsat(k)=soiliqwm(k)+qmin
1619 endif
1621 else
1622 soilicem(k)=0.
1623 soiliqwm(k)=soilmoism(k)
1624 lwsat(k)=dqm+qmin
1625 fwsat(k)=0.
1626 endif
1628 ENDDO
1630 do k=1,nzs
1631 if(soilice(k).gt.0.) then
1632 smfrkeep(k)=soilice(k)
1633 else
1634 smfrkeep(k)=soilmois(k)/riw
1635 endif
1636 enddo
1638 !******************************************************************
1639 ! SOILPROP computes thermal diffusivity, and diffusional and
1640 ! hydraulic condeuctivities
1641 !******************************************************************
1642 CALL SOILPROP( &
1643 !--- input variables
1644 nzs,fwsat,lwsat,tav,keepfr, &
1645 soilmois,soiliqw,soilice, &
1646 soilmoism,soiliqwm,soilicem, &
1647 !--- soil fixed fields
1648 QWRTZ,rhocs,dqm,qmin,psis,bclh,ksat, &
1649 !--- constants
1650 riw,xlmelt,CP,G0_P,cvw,ci, &
1651 kqwrtz,kice,kwt, &
1652 !--- output variables
1653 thdif,diffu,hydro,cap)
1655 !********************************************************************
1656 !--- CALCULATION OF CANOPY WATER (Smirnova et al., 1996, EQ.16) AND DEW
1658 DRIP=0.
1659 DD1=0.
1661 FQ=QKMS
1663 Q1=-QKMS*RAS*(QVATM - QSG)
1665 DEW=0.
1666 IF(QVATM.GE.QSG)THEN
1667 DEW=FQ*(QVATM-QSG)
1668 ENDIF
1669 IF(DEW.NE.0.)THEN
1670 DD1=CST+DELT*(PRCPMS +DEW*RAS)*vegfrac
1671 ELSE
1672 DD1=CST+ &
1673 DELT*(PRCPMS+RAS*FQ*(QVATM-QSG) &
1674 *(CST/SAT)**CN)*vegfrac
1675 ENDIF
1677 IF(DD1.LT.0.) DD1=0.
1678 if(vegfrac.eq.0.)then
1679 cst=0.
1680 drip=0.
1681 endif
1682 IF (vegfrac.GT.0.) THEN
1683 CST=DD1
1684 IF(CST.GT.SAT) THEN
1685 CST=SAT
1686 DRIP=DD1-SAT
1687 ENDIF
1688 ENDIF
1690 !--- WETCAN is the fraction of vegetated area covered by canopy
1691 !--- water, and DRYCAN is the fraction of vegetated area where
1692 !--- transpiration may take place.
1694 WETCAN=(CST/SAT)**CN
1695 DRYCAN=1.-WETCAN
1697 !**************************************************************
1698 ! TRANSF computes transpiration function
1699 !**************************************************************
1700 CALL TRANSF( &
1701 !--- input variables
1702 nzs,nroot,soiliqw,tabs, &
1703 !--- soil fixed fields
1704 dqm,qmin,ref,wilt,zshalf, &
1705 !--- output variables
1706 tranf,transum)
1709 !--- Save soil temp and moisture from the beginning of time step
1710 do k=1,nzs
1711 told(k)=tso(k)
1712 smold(k)=soilmois(k)
1713 enddo
1715 !**************************************************************
1716 ! SOILTEMP soilves heat budget and diffusion eqn. in soil
1717 !**************************************************************
1719 CALL SOILTEMP( &
1720 !--- input variables
1721 i,j,iland,isoil, &
1722 delt,ktau,conflx,nzs,nddzs,nroot, &
1723 PRCPMS,RAINF, &
1724 PATM,TABS,QVATM,QCATM,EMISS,RNET, &
1725 QKMS,TKMS,PC,rho,vegfrac, &
1726 thdif,cap,drycan,wetcan, &
1727 transum,dew,mavail, &
1728 !--- soil fixed fields
1729 dqm,qmin,bclh,zsmain,zshalf,DTDZS,tbq, &
1730 !--- constants
1731 xlv,CP,G0_P,cvw,stbolt, &
1732 !--- output variables
1733 tso,soilt,qvg,qsg,qcg)
1735 !************************************************************************
1737 !--- CALCULATION OF DEW USING NEW VALUE OF QSG OR TRANSP IF NO DEW
1738 ETT1=0.
1739 DEW=0.
1741 IF(QVATM.GE.QSG)THEN
1742 DEW=QKMS*(QVATM-QSG)
1743 DO K=1,NZS
1744 TRANSP(K)=0.
1745 ENDDO
1746 ELSE
1747 DO K=1,NROOT
1748 TRANSP(K)=VEGFRAC*RAS*QKMS* &
1749 (QVATM-QSG)* &
1750 PC*TRANF(K)*DRYCAN/ZSHALF(NROOT+1)
1751 IF(TRANSP(K).GT.0.) TRANSP(K)=0.
1752 ETT1=ETT1-TRANSP(K)
1753 ENDDO
1754 DO k=nroot+1,nzs
1755 transp(k)=0.
1756 enddo
1757 ENDIF
1759 !-- Recalculating of volumetric content of frozen water in soil
1760 DO K=1,NZS
1761 !- main levels
1762 tln=log(tso(k)/273.15)
1763 if(tln.lt.0.) then
1764 soiliqw(k)=(dqm+qmin)*(XLMELT* &
1765 (tso(k)-273.15)/tso(k)/9.81/psis) &
1766 **(-1./bclh)-qmin
1767 soiliqw(k)=max(0.,soiliqw(k))
1768 soiliqw(k)=min(soiliqw(k),soilmois(k))
1769 soilice(k)=(soilmois(k)-soiliqw(k))/riw
1770 !---- melting and freezing is balanced, soil ice cannot increase
1771 if(keepfr(k).eq.1.) then
1772 soilice(k)=min(soilice(k),smfrkeep(k))
1773 soiliqw(k)=max(0.,soilmois(k)-soilice(k)*riw)
1774 endif
1776 else
1777 soilice(k)=0.
1778 soiliqw(k)=soilmois(k)
1779 endif
1780 ENDDO
1782 !*************************************************************************
1783 ! SOILMOIST solves moisture budget (Smirnova et al., 1996, EQ.22,28)
1784 ! and Richards eqn.
1785 !*************************************************************************
1786 CALL SOILMOIST ( &
1787 !-- input
1788 delt,nzs,nddzs,DTDZS,DTDZS2,RIW, &
1789 zsmain,zshalf,diffu,hydro, &
1790 QSG,QVG,QCG,QCATM,QVATM,-PRCPMS, &
1791 QKMS,TRANSP,DRIP,DEW,0.,SOILICE,VEGFRAC, &
1792 ! QKMS,TRANSP,DRIP,DEW,0.,SOILICE,VEGFRAC, &
1793 !-- soil properties
1794 DQM,QMIN,REF,KSAT,RAS,INFMAX, &
1795 !-- output
1796 SOILMOIS,SOILIQW,MAVAIL,RUNOFF1, &
1797 RUNOFF2,INFILTRP)
1799 !--- KEEPFR is 1 when the temperature and moisture in soil
1800 !--- are both increasing. In this case soil ice should not
1801 !--- be increasing according to the freezing curve.
1802 !--- Some part of ice is melted, but additional water is
1803 !--- getting frozen. Thus, only structure of frozen soil is
1804 !--- changed, and phase changes are not affecting the heat
1805 !--- transfer. This situation may happen when it rains on the
1806 !--- frozen soil.
1808 do k=1,nzs
1809 if (soilice(k).gt.0.) then
1810 if(tso(k).gt.told(k).and.soilmois(k).gt.smold(k)) then
1811 keepfr(k)=1.
1812 else
1813 keepfr(k)=0.
1814 endif
1815 endif
1816 enddo
1817 !--- THE DIAGNOSTICS OF SURFACE FLUXES
1819 T3 = STBOLT*SOILT*SOILT*SOILT
1820 UPFLUX = T3 *SOILT
1821 XINET = EMISS*(GLW-UPFLUX)
1822 RNET = GSW + XINET
1823 HFT=-TKMS*CP*RHO*(TABS-SOILT) &
1824 *(P1000mb*0.00001/Patm)**ROVCP
1825 Q1=-QKMS*RAS*(QVATM - QSG)
1826 IF (Q1.LE.0.) THEN
1827 ! --- condensation
1828 EC1=0.
1829 EDIR1=0.
1830 ETT1=0.
1831 !-- moisture flux for coupling with PBL
1832 EETA=-QKMS*RAS*(QVATM/(1.+QVATM) - QSG/(1.+QSG))*1.E3
1833 QFX= XLV*EETA
1834 !-- actual moisture flux from RUC LSM
1835 EETA= RHO*DEW
1836 ELSE
1837 ! --- evaporation
1838 EDIR1 =-(1.-vegfrac)*QKMS*RAS* &
1839 (QVATM-QVG)
1840 EC1 = Q1 * WETCAN
1841 CMC2MS=CST/DELT
1842 if(EC1.gt.CMC2MS) cst=0.
1843 EC1=MIN(CMC2MS,EC1)*vegfrac
1844 !-- moisture flux for coupling with PBL
1845 EETA=-QKMS*RAS*(QVATM/(1.+QVATM) - QVG/(1.+QVG))*1.E3
1846 ! to convert from kg m-2 s-1 to m s-1: 1/rho water=1.e-3************
1847 QFX= XLV * EETA
1848 !-- actual moisture flux from RUC LSM
1849 EETA = (EDIR1 + EC1 + ETT1)*1.E3
1850 ENDIF
1851 EVAPL=EETA
1852 S=THDIF(1)*CAP(1)*DZSTOP*(TSO(1)-TSO(2))
1853 HFX=HFT
1854 FLTOT=RNET-HFT-QFX-S
1856 222 CONTINUE
1858 1123 FORMAT(I5,8F12.3)
1859 1133 FORMAT(I7,8E12.4)
1860 123 format(i6,f6.2,7f8.1)
1861 122 FORMAT(1X,2I3,6F8.1,F8.3,F8.2)
1864 ! RETURN
1865 ! END
1866 !-------------------------------------------------------------------
1867 END SUBROUTINE SOIL
1868 !-------------------------------------------------------------------
1870 SUBROUTINE SICE ( &
1871 !--- input variables
1872 i,j,iland,isoil,delt,ktau,conflx,nzs,nddzs,nroot, &
1873 PRCPMS,RAINF,PATM,QVATM,QCATM,GLW,GSW, &
1874 EMISS,RNET,QKMS,TKMS,rho, &
1875 !--- sea ice parameters
1876 tice,rhosice,capice,thdifice, &
1877 zsmain,zshalf,DTDZS,DTDZS2,tbq, &
1878 !--- constants
1879 xlv,CP,rovcp,cw,stbolt,tabs, &
1880 !--- output variables
1881 tso,dew,soilt,qvg,qsg,qcg, &
1882 eeta,qfx,hfx,s,evapl,prcpl &
1883 )
1885 !*****************************************************************
1886 ! Energy budget and heat diffusion eqns. for
1887 ! sea ice
1888 !*************************************************************
1890 IMPLICIT NONE
1891 !-----------------------------------------------------------------
1893 !--- input variables
1895 INTEGER, INTENT(IN ) :: nroot,ktau,nzs , &
1896 nddzs !nddzs=2*(nzs-2)
1897 INTEGER, INTENT(IN ) :: i,j,iland,isoil
1898 REAL, INTENT(IN ) :: DELT,CONFLX
1899 !--- 3-D Atmospheric variables
1900 REAL, &
1901 INTENT(IN ) :: PATM, &
1902 QVATM, &
1903 QCATM
1904 !--- 2-D variables
1905 REAL, &
1906 INTENT(IN ) :: GLW, &
1907 GSW, &
1908 EMISS, &
1909 RHO, &
1910 QKMS, &
1911 TKMS
1912 !--- sea ice properties
1913 REAL, DIMENSION(1:NZS) , &
1914 INTENT(IN ) :: &
1915 tice, &
1916 rhosice, &
1917 capice, &
1918 thdifice
1921 REAL, INTENT(IN ) :: &
1922 CW, &
1923 XLV
1926 REAL, DIMENSION(1:NZS), INTENT(IN) :: ZSMAIN, &
1927 ZSHALF, &
1928 DTDZS2
1930 REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS
1932 REAL, DIMENSION(1:4001), INTENT(IN) :: TBQ
1935 !--- input/output variables
1936 !----soil temperature
1937 REAL, DIMENSION( 1:nzs ), INTENT(INOUT) :: TSO
1938 !-------- 2-d variables
1939 REAL, &
1940 INTENT(INOUT) :: DEW, &
1941 EETA, &
1942 EVAPL, &
1943 PRCPL, &
1944 QVG, &
1945 QSG, &
1946 QCG, &
1947 RNET, &
1948 QFX, &
1949 HFX, &
1950 S, &
1951 SOILT
1953 !--- Local variables
1954 REAL :: x,x1,x2,x4,tn,denom
1955 REAL :: RAINF, PRCPMS , &
1956 TABS, T3, UPFLUX, XINET
1958 REAL :: CP,rovcp,G0,LV,STBOLT,xlmelt,dzstop , &
1959 epot,fltot,ft,fq,hft,ras,cvw
1961 REAL :: FKT,D1,D2,D9,D10,DID,R211,R21,R22,R6,R7,D11 , &
1962 PI,H,FKQ,R210,AA,BB,PP,Q1,QS1,TS1,TQ2,TX2 , &
1963 TDENOM
1965 REAL :: AA1,RHCS
1968 REAL, DIMENSION(1:NZS) :: cotso,rhtso
1970 INTEGER :: nzs1,nzs2,k,k1,kn,kk
1972 !-----------------------------------------------------------------
1974 !-- define constants
1975 ! STBOLT=5.670151E-8
1976 XLMELT=3.35E+5
1977 cvw=cw
1979 prcpl=prcpms
1981 NZS1=NZS-1
1982 NZS2=NZS-2
1983 dzstop=1./(zsmain(2)-zsmain(1))
1984 RAS=RHO*1.E-3
1986 do k=1,nzs
1987 cotso(k)=0.
1988 rhtso(k)=0.
1989 enddo
1991 cotso(1)=0.
1992 rhtso(1)=TSO(NZS)
1994 DO 33 K=1,NZS2
1995 KN=NZS-K
1996 K1=2*KN-3
1997 X1=DTDZS(K1)*THDIFICE(KN-1)
1998 X2=DTDZS(K1+1)*THDIFICE(KN)
1999 FT=TSO(KN)+X1*(TSO(KN-1)-TSO(KN)) &
2000 -X2*(TSO(KN)-TSO(KN+1))
2001 DENOM=1.+X1+X2-X2*cotso(K)
2002 cotso(K+1)=X1/DENOM
2003 rhtso(K+1)=(FT+X2*rhtso(K))/DENOM
2004 33 CONTINUE
2006 !************************************************************************
2007 !--- THE HEAT BALANCE EQUATION (Smirnova et al., 1996, EQ. 21,26)
2008 RHCS=CAPICE(1)
2009 H=1.
2010 FKT=TKMS
2011 D1=cotso(NZS1)
2012 D2=rhtso(NZS1)
2013 TN=SOILT
2014 D9=THDIFICE(1)*RHCS*dzstop
2015 D10=TKMS*CP*RHO
2016 R211=.5*CONFLX/DELT
2017 R21=R211*CP*RHO
2018 R22=.5/(THDIFICE(1)*DELT*dzstop**2)
2019 R6=EMISS *STBOLT*.5*TN**4
2020 R7=R6/TN
2021 D11=RNET+R6
2022 TDENOM=D9*(1.-D1+R22)+D10+R21+R7 &
2023 +RAINF*CVW*PRCPMS
2024 FKQ=QKMS*RHO
2025 R210=R211*RHO
2026 AA=XLS*(FKQ+R210)/TDENOM
2027 BB=(D10*TABS+R21*TN+XLS*(QVATM*FKQ &
2028 +R210*QVG)+D11+D9*(D2+R22*TN) &
2029 +RAINF*CVW*PRCPMS*max(273.15,TABS) &
2030 )/TDENOM
2031 AA1=AA
2032 PP=PATM*1.E3
2033 AA1=AA1/PP
2034 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
2035 PRINT *,' VILKA-SEAICE1'
2036 print *,'D10,TABS,R21,TN,QVATM,FKQ', &
2037 D10,TABS,R21,TN,QVATM,FKQ
2038 print *,'RNET, EMISS, STBOLT, SOILT',RNET, EMISS, STBOLT, SOILT
2039 print *,'R210,QVG,D11,D9,D2,R22,RAINF,CVW,PRCPMS,TDENOM', &
2040 R210,QVG,D11,D9,D2,R22,RAINF,CVW,PRCPMS,TDENOM
2041 print *,'tn,aa1,bb,pp,fkq,r210', &
2042 tn,aa1,bb,pp,fkq,r210
2043 ENDIF
2044 CALL VILKA(TN,AA1,BB,PP,QS1,TS1,TBQ,KTAU,i,j,iland,isoil)
2045 !--- it is saturation over sea ice
2046 QVG=QS1
2047 QSG=QS1
2048 TSO(1)=min(271.4,TS1)
2049 QCG=0.
2050 !--- sea ice melting is not included in this simple approach
2051 !--- SOILT - skin temperature
2052 SOILT=TSO(1)
2053 !---- Final solution for soil temperature - TSO
2054 DO K=2,NZS
2055 KK=NZS-K+1
2056 TSO(K)=min(271.4,rhtso(KK)+cotso(KK)*TSO(K-1))
2057 END DO
2058 !--- CALCULATION OF DEW USING NEW VALUE OF QSG OR TRANSP IF NO DEW
2059 DEW=0.
2061 !--- THE DIAGNOSTICS OF SURFACE FLUXES
2062 T3 = STBOLT*SOILT*SOILT*SOILT
2063 UPFLUX = T3 *SOILT
2064 XINET = EMISS*(GLW-UPFLUX)
2065 RNET = GSW + XINET
2066 HFT=-TKMS*CP*RHO*(TABS-SOILT) &
2067 *(P1000mb*0.00001/Patm)**ROVCP
2068 Q1=-QKMS*RAS*(QVATM - QSG)
2069 IF (Q1.LE.0.) THEN
2070 ! --- condensation
2071 !-- moisture flux for coupling with PBL
2072 EETA=-QKMS*RAS*(QVATM/(1.+QVATM) - QSG/(1.+QSG))*1.E3
2073 QFX= XLS*EETA
2074 !-- actual moisture flux from RUC LSM
2075 DEW=QKMS*(QVATM-QSG)
2076 EETA= RHO*DEW
2077 ELSE
2078 ! --- evaporation
2079 !-- moisture flux for coupling with PBL
2080 EETA=-QKMS*RAS*(QVATM/(1.+QVATM) - QVG/(1.+QVG))*1.E3
2081 ! EETA=-QKMS*RAS*(QVATM/(1.+QVATM) - QSG/(1.+QSG))*1.E3
2082 ! to convert from kg m-2 s-1 to m s-1: 1/rho water=1.e-3************
2083 QFX= XLS * EETA
2084 !-- actual moisture flux from RUC LSM
2085 EETA = Q1*1.E3
2086 ENDIF
2087 EVAPL=EETA
2088 S=THDIFICE(1)*CAPICE(1)*DZSTOP*(TSO(1)-TSO(2))
2089 HFX=HFT
2090 FLTOT=RNET-HFT-QFX-S
2092 !-------------------------------------------------------------------
2093 END SUBROUTINE SICE
2094 !-------------------------------------------------------------------
2098 SUBROUTINE SNOWSOIL ( &
2099 !--- input variables
2100 i,j,isoil,delt,ktau,conflx,nzs,nddzs,nroot, &
2101 meltfactor,rhonewsn, & ! new
2102 ILAND,PRCPMS,RAINF,NEWSNOW,snhei,SNWE,SNOWFRAC, &
2103 RHOSN, &
2104 PATM,QVATM,QCATM, &
2105 GLW,GSW,EMISS,RNET,IVGTYP, &
2106 QKMS,TKMS,PC,cst,rho,vegfrac,alb,znt, &
2107 MYJ, &
2108 !--- soil fixed fields
2109 QWRTZ,rhocs,dqm,qmin,ref,wilt,psis,bclh,ksat, &
2110 sat,cn,zsmain,zshalf,DTDZS,DTDZS2,tbq, &
2111 !--- constants
2112 xlv,CP,rovcp,G0_P,cw,stbolt,TABS, &
2113 KQWRTZ,KICE,KWT, &
2114 !--- output variables
2115 ilnb,snweprint,snheiprint,rsm, &
2116 soilmois,tso,smfrkeep,keepfr, &
2117 dew,soilt,soilt1,tsnav, &
2118 qvg,qsg,qcg,SMELT,SNOH,SNFLX,SNOM, &
2119 edir1,ec1,ett1,eeta,qfx,hfx,s,sublim, &
2120 prcpl,runoff1,runoff2,mavail,soilice, &
2121 soiliqw,infiltrp )
2123 !***************************************************************
2124 ! Energy and moisture budget for snow, heat diffusion eqns.
2125 ! in snow and soil, Richards eqn. for soil covered with snow
2126 !
2127 ! DELT - time step (s)
2128 ! ktau - numver of time step
2129 ! CONFLX - depth of constant flux layer (m)
2130 ! J,I - the location of grid point
2131 ! IME, JME, NZS - dimensions of the domain
2132 ! NROOT - number of levels within the root zone
2133 ! PRCPMS - precipitation rate in m/s
2134 ! NEWSNOW - pcpn in soilid form (m)
2135 ! SNHEI, SNWE - snow height and snow water equivalent (m)
2136 ! RHOSN - snow density (kg/m-3)
2137 ! PATM - pressure (bar)
2138 ! QVATM,QCATM - cloud and water vapor mixing ratio
2139 ! at the first atm. level (kg/kg)
2140 ! GLW, GSW - incoming longwave and absorbed shortwave
2141 ! radiation at the surface (W/m^2)
2142 ! EMISS,RNET - emissivity (0-1) of the ground surface and net
2143 ! radiation at the surface (W/m^2)
2144 ! QKMS - exchange coefficient for water vapor in the
2145 ! surface layer (m/s)
2146 ! TKMS - exchange coefficient for heat in the surface
2147 ! layer (m/s)
2148 ! PC - plant coefficient (resistance) (0-1)
2149 ! RHO - density of atmosphere near surface (kg/m^3)
2150 ! VEGFRAC - greeness fraction (0-1)
2151 ! RHOCS - volumetric heat capacity of dry soil (J/m^3/K)
2152 ! DQM, QMIN - porosity minus residual soil moisture QMIN (m^3/m^3)
2153 ! REF, WILT - field capacity soil moisture and the
2154 ! wilting point (m^3/m^3)
2155 ! PSIS - matrix potential at saturation (m)
2156 ! BCLH - exponent for Clapp-Hornberger parameterization
2157 ! KSAT - saturated hydraulic conductivity (m/s)
2158 ! SAT - maximum value of water intercepted by canopy (m)
2159 ! CN - exponent for calculation of canopy water
2160 ! ZSMAIN - main levels in soil (m)
2161 ! ZSHALF - middle of the soil layers (m)
2162 ! DTDZS,DTDZS2 - dt/(2.*dzshalf*dzmain) and dt/dzshalf in soil
2163 ! TBQ - table to define saturated mixing ration
2164 ! of water vapor for given temperature and pressure
2165 ! ilnb - number of layers in snow
2166 ! rsm - liquid water inside snow pack (m)
2167 ! SOILMOIS,TSO - soil moisture (m^3/m^3) and temperature (K)
2168 ! DEW - dew in (kg/m^2 s)
2169 ! SOILT - skin temperature (K)
2170 ! SOILT1 - snow temperature at 7.5 cm depth (K)
2171 ! TSNAV - average temperature of snow pack (C)
2172 ! QSG,QVG,QCG - saturated mixing ratio, mixing ratio of
2173 ! water vapor and cloud at the ground
2174 ! surface, respectively (kg/kg)
2175 ! EDIR1, EC1, ETT1, EETA - direct evaporation, evaporation of
2176 ! canopy water, transpiration (kg m-2 s-1) and total
2177 ! evaporation in (m s-1).
2178 ! QFX, HFX - latent and sensible heat fluxes (W/m^2)
2179 ! S - soil heat flux in the top layer (W/m^2)
2180 ! SUBLIM - snow sublimation (kg/m^2/s)
2181 ! RUNOFF1 - surface runoff (m/s)
2182 ! RUNOFF2 - underground runoff (m)
2183 ! MAVAIL - moisture availability in the top soil layer (0-1)
2184 ! SOILICE - content of soil ice in soil layers (m^3/m^3)
2185 ! SOILIQW - lliquid water in soil layers (m^3/m^3)
2186 ! INFILTRP - infiltration flux from the top of soil domain (m/s)
2187 ! XINET - net long-wave radiation (W/m^2)
2188 !
2189 !*******************************************************************
2191 IMPLICIT NONE
2192 !-------------------------------------------------------------------
2193 !--- input variables
2195 INTEGER, INTENT(IN ) :: nroot,ktau,nzs , &
2196 nddzs !nddzs=2*(nzs-2)
2197 INTEGER, INTENT(IN ) :: i,j,isoil
2199 REAL, INTENT(IN ) :: DELT,CONFLX,PRCPMS , &
2200 RAINF,NEWSNOW,RHONEWSN,meltfactor
2202 LOGICAL, INTENT(IN ) :: myj
2204 !--- 3-D Atmospheric variables
2205 REAL, &
2206 INTENT(IN ) :: PATM, &
2207 QVATM, &
2208 QCATM
2209 !--- 2-D variables
2210 REAL , &
2211 INTENT(IN ) :: GLW, &
2212 GSW, &
2213 RHO, &
2214 PC, &
2215 VEGFRAC, &
2216 QKMS, &
2217 TKMS
2219 INTEGER, INTENT(IN ) :: IVGTYP
2220 !--- soil properties
2221 REAL , &
2222 INTENT(IN ) :: RHOCS, &
2223 BCLH, &
2224 DQM, &
2225 KSAT, &
2226 PSIS, &
2227 QMIN, &
2228 QWRTZ, &
2229 REF, &
2230 SAT, &
2231 WILT
2233 REAL, INTENT(IN ) :: CN, &
2234 CW, &
2235 XLV, &
2236 G0_P, &
2237 KQWRTZ, &
2238 KICE, &
2239 KWT
2242 REAL, DIMENSION(1:NZS), INTENT(IN) :: ZSMAIN, &
2243 ZSHALF, &
2244 DTDZS2
2246 REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS
2248 REAL, DIMENSION(1:4001), INTENT(IN) :: TBQ
2251 !--- input/output variables
2252 !-------- 3-d soil moisture and temperature
2253 REAL, DIMENSION( 1:nzs ) , &
2254 INTENT(INOUT) :: TSO, &
2255 SOILMOIS, &
2256 SMFRKEEP
2258 REAL, DIMENSION( 1:nzs ) , &
2259 INTENT(INOUT) :: KEEPFR
2262 INTEGER, INTENT(INOUT) :: ILAND
2265 !-------- 2-d variables
2266 REAL , &
2267 INTENT(INOUT) :: DEW, &
2268 CST, &
2269 EDIR1, &
2270 EC1, &
2271 ETT1, &
2272 EETA, &
2273 RHOSN, &
2274 SUBLIM, &
2275 PRCPL, &
2276 ALB, &
2277 EMISS, &
2278 ZNT, &
2279 MAVAIL, &
2280 QVG, &
2281 QSG, &
2282 QCG, &
2283 QFX, &
2284 HFX, &
2285 S, &
2286 RUNOFF1, &
2287 RUNOFF2, &
2288 SNWE, &
2289 SNHEI, &
2290 SMELT, &
2291 SNOM, &
2292 SNOH, &
2293 SNFLX, &
2294 SOILT, &
2295 SOILT1, &
2296 SNOWFRAC, &
2297 TSNAV
2299 INTEGER, INTENT(INOUT) :: ILNB
2301 !-------- 1-d variables
2302 REAL, DIMENSION(1:NZS), INTENT(OUT) :: SOILICE, &
2303 SOILIQW
2305 REAL, INTENT(OUT) :: RSM, &
2306 SNWEPRINT, &
2307 SNHEIPRINT
2308 !--- Local variables
2311 INTEGER :: nzs1,nzs2,k
2313 REAL :: INFILTRP, TRANSUM , &
2314 SNTH, NEWSN , &
2315 TABS, T3, UPFLUX, XINET , &
2316 BETA, SNWEPR,EPDT,PP
2317 REAL :: CP,rovcp,G0,LV,xlvm,STBOLT,xlmelt,dzstop , &
2318 can,epot,fac,fltot,ft,fq,hft , &
2319 q1,ras,rhoice,sph , &
2320 trans,zn,ci,cvw,tln,tavln,pi , &
2321 DD1,CMC2MS,DRYCAN,WETCAN , &
2322 INFMAX,RIW,DELTSN,H,UMVEG
2324 REAL, DIMENSION(1:NZS) :: transp,cap,diffu,hydro , &
2325 thdif,tranf,tav,soilmoism , &
2326 soilicem,soiliqwm,detal , &
2327 fwsat,lwsat,told,smold
2328 REAL :: drip
2330 REAL :: RNET
2332 !-----------------------------------------------------------------
2334 cvw=cw
2335 XLMELT=3.35E+5
2336 !-- the next line calculates heat of sublimation of water vapor
2337 XLVm=XLV+XLMELT
2338 ! STBOLT=5.670151E-8
2340 !--- SNOW flag -- 99
2341 ILAND=99
2343 !--- DELTSN - is the threshold for splitting the snow layer into 2 layers.
2344 !--- With snow density 400 kg/m^3, this threshold is equal to 7.5 cm,
2345 !--- equivalent to 0.03 m SNWE. For other snow densities the threshold is
2346 !--- computed using SNWE=0.03 m and current snow density.
2347 !--- SNTH - the threshold below which the snow layer is combined with
2348 !--- the top soil layer. SNTH is computed using snwe=0.016 m, and
2349 !--- equals 4 cm for snow density 400 kg/m^3.
2351 DELTSN=0.0301*1.e3/rhosn
2352 snth=0.01601*1.e3/rhosn
2354 !tgs when the snow depth is marginlly higher than DELTSN,
2355 ! reset DELTSN to half of snow depth.
2356 IF(SNHEI.GE.DELTSN+SNTH) THEN
2357 ! 2-layer model
2358 if(snhei-deltsn-snth.lt.snth) deltsn=0.5*(snhei-snth)
2359 ENDIF
2361 RHOICE=900.
2362 CI=RHOICE*2100.
2363 RAS=RHO*1.E-3
2364 RIW=rhoice*1.e-3
2365 MAVAIL=1.
2366 RSM=0.
2368 DO K=1,NZS
2369 TRANSP (K)=0.
2370 soilmoism (k)=0.
2371 soiliqwm (k)=0.
2372 soilice (k)=0.
2373 soilicem (k)=0.
2374 lwsat (k)=0.
2375 fwsat (k)=0.
2376 tav (k)=0.
2377 cap (k)=0.
2378 diffu (k)=0.
2379 hydro (k)=0.
2380 thdif (k)=0.
2381 tranf (k)=0.
2382 detal (k)=0.
2383 told (k)=0.
2384 smold (k)=0.
2385 ENDDO
2387 snweprint=0.
2388 snheiprint=0.
2389 prcpl=prcpms
2391 !*** DELTSN is the depth of the top layer of snow where
2392 !*** there is a temperature gradient, the rest of the snow layer
2393 !*** is considered to have constant temperature
2396 NZS1=NZS-1
2397 NZS2=NZS-2
2398 DZSTOP=1./(zsmain(2)-zsmain(1))
2400 !----- THE CALCULATION OF THERMAL DIFFUSIVITY, DIFFUSIONAL AND ---
2401 !----- HYDRAULIC CONDUCTIVITY (SMIRNOVA ET AL. 1996, EQ.2,5,6) ---
2402 !tgs - the following loop is added to define the amount of frozen
2403 !tgs - water in soil if there is any
2404 DO K=1,NZS
2406 tln=log(tso(k)/273.15)
2407 if(tln.lt.0.) then
2408 soiliqw(k)=(dqm+qmin)*(XLMELT* &
2409 (tso(k)-273.15)/tso(k)/9.81/psis) &
2410 **(-1./bclh)-qmin
2411 soiliqw(k)=max(0.,soiliqw(k))
2412 soiliqw(k)=min(soiliqw(k),soilmois(k))
2413 soilice(k)=(soilmois(k)-soiliqw(k))/riw
2415 !---- melting and freezing is balanced, soil ice cannot increase
2416 if(keepfr(k).eq.1.) then
2417 soilice(k)=min(soilice(k),smfrkeep(k))
2418 soiliqw(k)=max(0.,soilmois(k)-soilice(k)*rhoice*1.e-3)
2419 endif
2421 else
2422 soilice(k)=0.
2423 soiliqw(k)=soilmois(k)
2424 endif
2426 ENDDO
2428 DO K=1,NZS1
2430 tav(k)=0.5*(tso(k)+tso(k+1))
2431 soilmoism(k)=0.5*(soilmois(k)+soilmois(k+1))
2432 tavln=log(tav(k)/273.15)
2434 if(tavln.lt.0.) then
2435 soiliqwm(k)=(dqm+qmin)*(XLMELT* &
2436 (tav(k)-273.15)/tav(k)/9.81/psis) &
2437 **(-1./bclh)-qmin
2438 fwsat(k)=dqm-soiliqwm(k)
2439 lwsat(k)=soiliqwm(k)+qmin
2440 soiliqwm(k)=max(0.,soiliqwm(k))
2441 soiliqwm(k)=min(soiliqwm(k), soilmoism(k))
2442 soilicem(k)=(soilmoism(k)-soiliqwm(k))/riw
2443 !---- melting and freezing is balanced, soil ice cannot increase
2444 if(keepfr(k).eq.1.) then
2445 soilicem(k)=min(soilicem(k), &
2446 0.5*(smfrkeep(k)+smfrkeep(k+1)))
2447 soiliqwm(k)=max(0.,soilmoism(k)-soilicem(k)*riw)
2448 fwsat(k)=dqm-soiliqwm(k)
2449 lwsat(k)=soiliqwm(k)+qmin
2450 endif
2452 else
2453 soilicem(k)=0.
2454 soiliqwm(k)=soilmoism(k)
2455 lwsat(k)=dqm+qmin
2456 fwsat(k)=0.
2458 endif
2459 ENDDO
2461 do k=1,nzs
2462 if(soilice(k).gt.0.) then
2463 smfrkeep(k)=soilice(k)
2464 else
2465 smfrkeep(k)=soilmois(k)/riw
2466 endif
2467 enddo
2469 !******************************************************************
2470 ! SOILPROP computes thermal diffusivity, and diffusional and
2471 ! hydraulic condeuctivities
2472 !******************************************************************
2473 CALL SOILPROP( &
2474 !--- input variables
2475 nzs,fwsat,lwsat,tav,keepfr, &
2476 soilmois,soiliqw,soilice, &
2477 soilmoism,soiliqwm,soilicem, &
2478 !--- soil fixed fields
2479 QWRTZ,rhocs,dqm,qmin,psis,bclh,ksat, &
2480 !--- constants
2481 riw,xlmelt,CP,G0_P,cvw,ci, &
2482 kqwrtz,kice,kwt, &
2483 !--- output variables
2484 thdif,diffu,hydro,cap)
2486 !********************************************************************
2487 !--- CALCULATION OF CANOPY WATER (Smirnova et al., 1996, EQ.16) AND DEW
2489 DRIP=0.
2490 SMELT=0.
2491 DD1=0.
2492 H=1.
2494 FQ=QKMS
2497 !--- If vegfrac.ne.0. then part of falling snow can be
2498 !--- intercepted by the canopy.
2500 DEW=0.
2501 UMVEG=1.-vegfrac
2502 EPOT = -FQ*(QVATM-QSG)
2504 IF(vegfrac.EQ.0.) then
2505 cst=0.
2506 drip=0.
2507 ELSE
2508 IF(EPOT.GE.0.) THEN
2509 ! Evaporation
2510 ! DD1=CST+(NEWSNOW*RHOSN*1.E-3 &
2511 DD1=CST+(NEWSNOW*RHOnewSN*1.E-3 &
2512 !-- this change will not let liquid waer be intercepted by the canopy....
2513 -DELT*(RAS*EPOT &
2514 ! -DELT*(-PRCPMS+RAS*EPOT &
2515 *(CST/SAT)**CN)) *vegfrac
2516 ELSE
2517 ! Sublimation
2518 DEW = - EPOT
2519 ! DD1=CST+(NEWSNOW*RHOSN*1.E-3+delt*( &
2520 DD1=CST+(NEWSNOW*RHOnewSN*1.E-3+delt*( &
2521 ! DD1=CST+(NEWSNOW*RHOSN*1.E-3+delt*(PRCPMS &
2522 +DEW*RAS)) *vegfrac
2523 ENDIF
2525 IF(DD1.LT.0.) DD1=0.
2526 IF (vegfrac.GT.0.) THEN
2527 CST=DD1
2528 IF(CST.GT.SAT*vegfrac) THEN
2529 CST=SAT*vegfrac
2530 DRIP=DD1-SAT*vegfrac
2531 ENDIF
2532 ENDIF
2535 !--- In SFCTMP NEWSNOW is added to SNHEI as if there is no vegetation
2536 !--- With vegetation part of NEWSNOW can be intercepted by canopy until
2537 !--- the saturation is reached. After the canopy saturation is reached
2538 !--- DRIP in the solid form will be added to SNOW cover.
2540 SNWE=(SNHEI-vegfrac*NEWSNOW)*RHOSN*1.E-3 &
2541 + DRIP
2543 ENDIF
2545 DRIP=0.
2546 SNHEI=SNWE*1.e3/RHOSN
2547 SNWEPR=SNWE
2549 ! check if all snow can evaporate during DT
2550 BETA=1.
2551 EPDT = EPOT * RAS *DELT*UMVEG
2552 IF(SNWEPR.LE.EPDT) THEN
2553 BETA=SNWEPR/max(1.e-8,EPDT)
2554 SNWE=0.
2555 SNHEI=0.
2556 ENDIF
2558 WETCAN=(CST/SAT)**CN
2559 DRYCAN=1.-WETCAN
2561 !**************************************************************
2562 ! TRANSF computes transpiration function
2563 !**************************************************************
2564 CALL TRANSF( &
2565 !--- input variables
2566 nzs,nroot,soiliqw,tabs, &
2567 !--- soil fixed fields
2568 dqm,qmin,ref,wilt,zshalf, &
2569 !--- output variables
2570 tranf,transum)
2572 !--- Save soil temp and moisture from the beginning of time step
2573 do k=1,nzs
2574 told(k)=tso(k)
2575 smold(k)=soilmois(k)
2576 enddo
2578 !**************************************************************
2579 ! SNOWTEMP solves heat budget and diffusion eqn. in soil
2580 !**************************************************************
2582 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
2583 print *, 'TSO before calling SNOWTEMP: ', tso
2584 ENDIF
2585 CALL SNOWTEMP( &
2586 !--- input variables
2587 i,j,iland,isoil, &
2588 delt,ktau,conflx,nzs,nddzs,nroot, &
2589 snwe,snwepr,snhei,newsnow,snowfrac, &
2590 beta,deltsn,snth,rhosn,rhonewsn,meltfactor, & ! add meltfactor
2591 PRCPMS,RAINF, &
2592 PATM,TABS,QVATM,QCATM, &
2593 GLW,GSW,EMISS,RNET, &
2594 QKMS,TKMS,PC,rho,vegfrac, &
2595 thdif,cap,drycan,wetcan,cst, &
2596 tranf,transum,dew,mavail, &
2597 !--- soil fixed fields
2598 dqm,qmin,psis,bclh, &
2599 zsmain,zshalf,DTDZS,tbq, &
2600 !--- constants
2601 xlvm,CP,rovcp,G0_P,cvw,stbolt, &
2602 !--- output variables
2603 snweprint,snheiprint,rsm, &
2604 tso,soilt,soilt1,tsnav,qvg,qsg,qcg, &
2605 smelt,snoh,snflx,ilnb)
2607 !************************************************************************
2608 !--- RECALCULATION OF DEW USING NEW VALUE OF QSG OR TRANSP IF NO DEW
2609 DEW=0.
2610 ETT1=0.
2611 PP=PATM*1.E3
2612 QSG= QSN(SOILT,TBQ)/PP
2613 EPOT = -FQ*(QVATM-QSG)
2614 IF(EPOT.GE.0.) THEN
2615 ! Evaporation
2616 DO K=1,NROOT
2617 TRANSP(K)=vegfrac*RAS*FQ*(QVATM-QSG) &
2618 *PC*tranf(K)*DRYCAN/zshalf(NROOT+1)
2619 IF(TRANSP(K).GT.0.) TRANSP(K)=0.
2620 ETT1=ETT1-TRANSP(K)
2621 ENDDO
2622 DO k=nroot+1,nzs
2623 transp(k)=0.
2624 enddo
2626 ELSE
2627 ! Sublimation
2628 DEW=-EPOT
2629 DO K=1,NZS
2630 TRANSP(K)=0.
2631 ENDDO
2632 ETT1=0.
2633 ENDIF
2635 !-- recalculating of frozen water in soil
2636 DO K=1,NZS
2637 tln=log(tso(k)/273.15)
2638 if(tln.lt.0.) then
2639 soiliqw(k)=(dqm+qmin)*(XLMELT* &
2640 (tso(k)-273.15)/tso(k)/9.81/psis) &
2641 **(-1./bclh)-qmin
2642 soiliqw(k)=max(0.,soiliqw(k))
2643 soiliqw(k)=min(soiliqw(k),soilmois(k))
2644 soilice(k)=(soilmois(k)-soiliqw(k))/riw
2645 !---- melting and freezing is balanced, soil ice cannot increase
2646 if(keepfr(k).eq.1.) then
2647 soilice(k)=min(soilice(k),smfrkeep(k))
2648 soiliqw(k)=max(0.,soilmois(k)-soilice(k)*riw)
2649 endif
2651 else
2652 soilice(k)=0.
2653 soiliqw(k)=soilmois(k)
2654 endif
2655 ENDDO
2657 !*************************************************************************
2658 !--- TQCAN FOR SOLUTION OF MOISTURE BALANCE (Smirnova et al. 1996, EQ.22,28)
2659 ! AND TSO,ETA PROFILES
2660 !*************************************************************************
2661 CALL SOILMOIST ( &
2662 !-- input
2663 delt,nzs,nddzs,DTDZS,DTDZS2,RIW, &
2664 zsmain,zshalf,diffu,hydro, &
2665 QSG,QVG,QCG,QCATM,QVATM,-PRCPMS, &
2666 0.,TRANSP,0., &
2667 0.,SMELT,soilice,vegfrac, &
2668 !-- soil properties
2669 DQM,QMIN,REF,KSAT,RAS,INFMAX, &
2670 !-- output
2671 SOILMOIS,SOILIQW,MAVAIL,RUNOFF1, &
2672 RUNOFF2,infiltrp)
2674 ! 4 Nov 07 - update CST for snow melt
2675 if(snwe.ne.0.) then
2676 CST=(1.-min(1.,smelt/snwe))*CST
2677 else
2678 CST=0.
2679 endif
2681 !-- Restore land-use parameters if snow is less than threshold
2682 IF(SNHEI.EQ.0.) then
2683 tsnav=soilt-273.15
2684 CALL SNOWFREE(ivgtyp,myj,emiss, &
2685 znt,iland)
2686 smelt=smelt+snwe/delt
2687 rsm=0.
2688 ! snwe=0.
2689 ENDIF
2691 ! 21apr2009
2692 ! SNOM goes into the passed-in ACSNOM variable in the grid derived type
2693 SNOM=SNOM+SMELT*DELT*1.e3
2695 !--- KEEPFR is 1 when the temperature and moisture in soil
2696 !--- are both increasing. In this case soil ice should not
2697 !--- be increasing according to the freezing curve.
2698 !--- Some part of ice is melted, but additional water is
2699 !--- getting frozen. Thus, only structure of frozen soil is
2700 !--- changed, and phase changes are not affecting the heat
2701 !--- transfer. This situation may happen when it rains on the
2702 !--- frozen soil.
2704 do k=1,nzs
2705 if (soilice(k).gt.0.) then
2706 if(tso(k).gt.told(k).and.soilmois(k).gt.smold(k)) then
2707 keepfr(k)=1.
2708 else
2709 keepfr(k)=0.
2710 endif
2711 endif
2712 enddo
2713 !--- THE DIAGNOSTICS OF SURFACE FLUXES
2715 T3 = STBOLT*SOILT*SOILT*SOILT
2716 UPFLUX = T3 *SOILT
2717 XINET = EMISS*(GLW-UPFLUX)
2718 RNET = GSW + XINET
2719 HFT=-TKMS*CP*RHO*(TABS-SOILT) &
2720 *(P1000mb*0.00001/Patm)**ROVCP
2721 Q1 = - FQ*RAS* (QVATM - QSG)
2723 IF (Q1.LT.0.) THEN
2724 ! --- condensation
2725 EDIR1=0.
2726 EC1=0.
2727 ETT1=0.
2728 ! --- condensation
2729 !-- moisture flux for coupling with PBL
2730 EETA=-QKMS*RAS*(QVATM/(1.+QVATM) - QSG/(1.+QSG))*1.E3
2731 QFX= XLVm*EETA
2732 !-- actual moisture flux from RUC LSM
2733 DEW=QKMS*(QVATM-QSG)
2734 EETA= RHO*DEW
2735 ELSE
2736 ! --- evaporation
2737 EDIR1 = Q1*UMVEG *BETA
2738 EC1 = Q1 * WETCAN
2739 CMC2MS=CST/DELT
2740 if(EC1.gt.CMC2MS) cst=0.
2741 EC1=MIN(CMC2MS,EC1)*vegfrac
2742 !-- moisture flux for coupling with PBL
2743 EETA=-QKMS*RAS*(QVATM/(1.+QVATM) - QVG/(1.+QVG))*1.E3
2744 ! EETA=-QKMS*RAS*(QVATM/(1.+QVATM) - QSG/(1.+QSG))*1.E3
2745 ! to convert from kg m-2 s-1 to m s-1: 1/rho water=1.e-3************
2746 QFX= XLVm * EETA
2747 !-- actual moisture flux from RUC LSM
2748 EETA = (EDIR1 + EC1 + ETT1)*1.E3
2749 ENDIF
2750 s=THDIF(1)*CAP(1)*dzstop*(tso(1)-tso(2))
2751 HFX=HFT
2752 FLTOT=RNET-HFT-QFX-S
2754 222 CONTINUE
2756 1123 FORMAT(I5,8F12.3)
2757 1133 FORMAT(I7,8E12.4)
2758 123 format(i6,f6.2,7f8.1)
2759 122 FORMAT(1X,2I3,6F8.1,F8.3,F8.2)
2762 ! RETURN
2763 ! END
2764 !-------------------------------------------------------------------
2765 END SUBROUTINE SNOWSOIL
2766 !-------------------------------------------------------------------
2768 SUBROUTINE SNOWSEAICE( &
2769 i,j,isoil,delt,ktau,conflx,nzs,nddzs, &
2770 meltfactor,rhonewsn, & ! new
2771 ILAND,PRCPMS,RAINF,NEWSNOW,snhei,SNWE,snowfrac, &
2772 RHOSN,PATM,QVATM,QCATM, &
2773 GLW,GSW,EMISS,RNET, &
2774 QKMS,TKMS,RHO, &
2775 !--- sea ice parameters
2776 ALB,ZNT, &
2777 tice,rhosice,capice,thdifice, &
2778 zsmain,zshalf,DTDZS,DTDZS2,tbq, &
2779 !--- constants
2780 xlv,CP,rovcp,cw,stbolt,tabs, &
2781 !--- output variables
2782 ilnb,snweprint,snheiprint,rsm,tso, &
2783 dew,soilt,soilt1,tsnav,qvg,qsg,qcg, &
2784 SMELT,SNOH,SNFLX,SNOM,eeta, &
2785 qfx,hfx,s,sublim,prcpl &
2786 )
2787 !***************************************************************
2788 ! Solving energy budget for snow on sea ice and heat diffusion
2789 ! eqns. in snow and sea ice
2790 !***************************************************************
2793 IMPLICIT NONE
2794 !-------------------------------------------------------------------
2795 !--- input variables
2797 INTEGER, INTENT(IN ) :: ktau,nzs , &
2798 nddzs !nddzs=2*(nzs-2)
2799 INTEGER, INTENT(IN ) :: i,j,isoil
2801 REAL, INTENT(IN ) :: DELT,CONFLX,PRCPMS , &
2802 RAINF,NEWSNOW,RHONEWSN,meltfactor
2803 real :: rhonewcsn
2805 !--- 3-D Atmospheric variables
2806 REAL, &
2807 INTENT(IN ) :: PATM, &
2808 QVATM, &
2809 QCATM
2810 !--- 2-D variables
2811 REAL , &
2812 INTENT(IN ) :: GLW, &
2813 GSW, &
2814 RHO, &
2815 QKMS, &
2816 TKMS
2818 !--- sea ice properties
2819 REAL, DIMENSION(1:NZS) , &
2820 INTENT(IN ) :: &
2821 tice, &
2822 rhosice, &
2823 capice, &
2824 thdifice
2826 REAL, INTENT(IN ) :: &
2827 CW, &
2828 XLV
2830 REAL, DIMENSION(1:NZS), INTENT(IN) :: ZSMAIN, &
2831 ZSHALF, &
2832 DTDZS2
2834 REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS
2836 REAL, DIMENSION(1:4001), INTENT(IN) :: TBQ
2838 !--- input/output variables
2839 !-------- 3-d soil moisture and temperature
2840 REAL, DIMENSION( 1:nzs ) , &
2841 INTENT(INOUT) :: TSO
2843 INTEGER, INTENT(INOUT) :: ILAND
2846 !-------- 2-d variables
2847 REAL , &
2848 INTENT(INOUT) :: DEW, &
2849 EETA, &
2850 RHOSN, &
2851 SUBLIM, &
2852 PRCPL, &
2853 ALB, &
2854 EMISS, &
2855 ZNT, &
2856 QVG, &
2857 QSG, &
2858 QCG, &
2859 QFX, &
2860 HFX, &
2861 S, &
2862 SNWE, &
2863 SNHEI, &
2864 SMELT, &
2865 SNOM, &
2866 SNOH, &
2867 SNFLX, &
2868 SOILT, &
2869 SOILT1, &
2870 SNOWFRAC, &
2871 TSNAV
2873 INTEGER, INTENT(INOUT) :: ILNB
2875 REAL, INTENT(OUT) :: RSM, &
2876 SNWEPRINT, &
2877 SNHEIPRINT
2878 !--- Local variables
2881 INTEGER :: nzs1,nzs2,k,k1,kn,kk
2882 REAL :: x,x1,x2,dzstop,ft,tn,denom
2884 REAL :: SNTH, NEWSN , &
2885 TABS, T3, UPFLUX, XINET , &
2886 BETA, SNWEPR,EPDT,PP
2887 REAL :: CP,rovcp,G0,LV,xlvm,STBOLT,xlmelt , &
2888 epot,fltot,fq,hft,q1,ras,rhoice,ci,cvw , &
2889 RIW,DELTSN,H
2891 REAL :: rhocsn,thdifsn, &
2892 xsn,ddzsn,x1sn,d1sn,d2sn,d9sn,r22sn
2894 REAL :: cotsn,rhtsn,xsn1,ddzsn1,x1sn1,ftsnow,denomsn
2895 REAL :: fso,fsn, &
2896 FKT,D1,D2,D9,D10,DID,R211,R21,R22,R6,R7,D11, &
2897 FKQ,R210,AA,BB,QS1,TS1,TQ2,TX2, &
2898 TDENOM,AA1,RHCS,H1,TSOB, SNPRIM, &
2899 SNODIF,SOH,TNOLD,QGOLD,SNOHGNEW
2900 REAL, DIMENSION(1:NZS) :: cotso,rhtso
2902 REAL :: RNET,rsmfrac,soiltfrac,hsn
2903 integer :: nmelt
2906 !-----------------------------------------------------------------
2907 XLMELT=3.35E+5
2908 !-- the next line calculates heat of sublimation of water vapor
2909 XLVm=XLV+XLMELT
2910 ! STBOLT=5.670151E-8
2912 !--- SNOW flag -- 99
2913 ILAND=99
2915 !--- DELTSN - is the threshold for splitting the snow layer into 2 layers.
2916 !--- With snow density 400 kg/m^3, this threshold is equal to 7.5 cm,
2917 !--- equivalent to 0.03 m SNWE. For other snow densities the threshold is
2918 !--- computed using SNWE=0.03 m and current snow density.
2919 !--- SNTH - the threshold below which the snow layer is combined with
2920 !--- the top sea ice layer. SNTH is computed using snwe=0.016 m, and
2921 !--- equals 4 cm for snow density 400 kg/m^3.
2923 DELTSN=0.0301*1.e3/rhosn
2924 snth=0.01601*1.e3/rhosn
2926 !tgs when the snow depth is marginlly higher than DELTSN,
2927 ! reset DELTSN to half of snow depth.
2928 IF(SNHEI.GE.DELTSN+SNTH) THEN
2929 ! 2-layer model
2930 if(snhei-deltsn-snth.lt.snth) deltsn=0.5*(snhei-snth)
2931 ENDIF
2934 RHOICE=900.
2935 CI=RHOICE*2100.
2936 RAS=RHO*1.E-3
2937 RIW=rhoice*1.e-3
2938 RSM=0.
2940 XLMELT=3.35E+5
2941 RHOCSN=2090.* RHOSN
2942 !18apr08 - add rhonewcsn
2943 RHOnewCSN=2090.* RHOnewSN
2944 THDIFSN = 0.265/RHOCSN
2945 RAS=RHO*1.E-3
2947 SOILTFRAC=SOILT
2949 SMELT=0.
2950 SOH=0.
2951 SNODIF=0.
2952 SNOH=0.
2953 SNOHGNEW=0.
2954 RSM = 0.
2955 RSMFRAC = 0.
2956 fsn=1.
2957 fso=0.
2958 hsn=snhei
2959 cvw=cw
2961 NZS1=NZS-1
2962 NZS2=NZS-2
2964 QGOLD=QVG
2965 TNOLD=SOILT
2966 DZSTOP=1./(ZSMAIN(2)-ZSMAIN(1))
2968 snweprint=0.
2969 snheiprint=0.
2970 prcpl=prcpms
2972 !*** DELTSN is the depth of the top layer of snow where
2973 !*** there is a temperature gradient, the rest of the snow layer
2974 !*** is considered to have constant temperature
2977 H=1.
2978 SMELT=0.
2980 FQ=QKMS
2981 SNHEI=SNWE*1.e3/RHOSN
2982 SNWEPR=SNWE
2984 ! check if all snow can evaporate during DT
2985 BETA=1.
2986 EPDT = EPOT * RAS *DELT
2987 IF(SNWEPR.LE.EPDT) THEN
2988 BETA=SNWEPR/max(1.e-8,EPDT)
2989 SNWE=0.
2990 SNHEI=0.
2991 ENDIF
2993 !******************************************************************************
2994 ! COEFFICIENTS FOR THOMAS ALGORITHM FOR TSO
2995 !******************************************************************************
2997 cotso(1)=0.
2998 rhtso(1)=TSO(NZS)
2999 DO 33 K=1,NZS2
3000 KN=NZS-K
3001 K1=2*KN-3
3002 X1=DTDZS(K1)*THDIFICE(KN-1)
3003 X2=DTDZS(K1+1)*THDIFICE(KN)
3004 FT=TSO(KN)+X1*(TSO(KN-1)-TSO(KN)) &
3005 -X2*(TSO(KN)-TSO(KN+1))
3006 DENOM=1.+X1+X2-X2*cotso(K)
3007 cotso(K+1)=X1/DENOM
3008 rhtso(K+1)=(FT+X2*rhtso(K))/DENOM
3009 33 CONTINUE
3010 !--- THE NZS element in COTSO and RHTSO will be for snow
3011 !--- There will be 2 layers in snow if it is deeper than DELTSN+SNTH
3012 IF(SNHEI.GE.SNTH) then
3013 if(snhei.le.DELTSN+SNTH) then
3014 !-- 1-layer snow model
3015 ilnb=1
3016 snprim=snhei
3017 soilt1=tso(1)
3018 tsob=tso(1)
3019 XSN = DELT/2./(zshalf(2)+0.5*SNPRIM)
3020 DDZSN = XSN / SNPRIM
3021 X1SN = DDZSN * thdifsn
3022 X2 = DTDZS(1)*THDIFICE(1)
3023 FT = TSO(1)+X1SN*(SOILT-TSO(1)) &
3024 -X2*(TSO(1)-TSO(2))
3025 DENOM = 1. + X1SN + X2 -X2*cotso(NZS1)
3026 cotso(NZS)=X1SN/DENOM
3027 rhtso(NZS)=(FT+X2*rhtso(NZS1))/DENOM
3028 cotsn=cotso(NZS)
3029 rhtsn=rhtso(NZS)
3030 !*** Average temperature of snow pack (C)
3031 tsnav=0.5*(soilt+tso(1)) &
3032 -273.15
3034 else
3035 !-- 2 layers in snow, SOILT1 is temperasture at DELTSN depth
3036 ilnb=2
3037 snprim=deltsn
3038 tsob=soilt1
3039 XSN = DELT/2./(0.5*SNHEI)
3040 XSN1= DELT/2./(zshalf(2)+0.5*(SNHEI-DELTSN))
3041 DDZSN = XSN / DELTSN
3042 DDZSN1 = XSN1 / (SNHEI-DELTSN)
3043 X1SN = DDZSN * thdifsn
3044 X1SN1 = DDZSN1 * thdifsn
3045 X2 = DTDZS(1)*THDIFICE(1)
3046 FT = TSO(1)+X1SN1*(SOILT1-TSO(1)) &
3047 -X2*(TSO(1)-TSO(2))
3048 DENOM = 1. + X1SN1 + X2 - X2*cotso(NZS1)
3049 cotso(nzs)=x1sn1/denom
3050 rhtso(nzs)=(ft+x2*rhtso(nzs1))/denom
3051 ftsnow = soilt1+x1sn*(soilt-soilt1) &
3052 -x1sn1*(soilt1-tso(1))
3053 denomsn = 1. + X1SN + X1SN1 - X1SN1*cotso(NZS)
3054 cotsn=x1sn/denomsn
3055 rhtsn=(ftsnow+X1SN1*rhtso(NZS))/denomsn
3056 !*** Average temperature of snow pack (C)
3057 tsnav=0.5/snhei*((soilt+soilt1)*deltsn &
3058 +(soilt1+tso(1))*(SNHEI-DELTSN)) &
3059 -273.15
3060 endif
3061 ENDIF
3063 IF(SNHEI.LT.SNTH.AND.SNHEI.GT.0.) then
3064 !--- snow is too thin to be treated separately, therefore it
3065 !--- is combined with the first sea ice layer.
3066 fsn=SNHEI/(SNHEI+zsmain(2))
3067 fso=1.-fsn
3068 soilt1=tso(1)
3069 tsob=tso(2)
3070 snprim=SNHEI+zsmain(2)
3071 XSN = DELT/2./((zshalf(3)-zsmain(2))+0.5*snprim)
3072 DDZSN = XSN /snprim
3073 X1SN = DDZSN * (fsn*thdifsn+fso*thdifice(1))
3074 X2=DTDZS(2)*THDIFICE(2)
3075 FT=TSO(2)+X1SN*(SOILT-TSO(2))- &
3076 X2*(TSO(2)-TSO(3))
3077 denom = 1. + x1sn + x2 - x2*cotso(nzs-2)
3078 cotso(nzs1) = x1sn/denom
3079 rhtso(nzs1)=(FT+X2*rhtso(NZS-2))/denom
3080 tsnav=0.5*(soilt+tso(1)) &
3081 -273.15
3082 ENDIF
3084 !************************************************************************
3085 !--- THE HEAT BALANCE EQUATION
3086 !18apr08 nmelt is the flag for melting, and SNOH is heat of snow phase changes
3087 nmelt=0
3088 SNOH=0.
3090 EPOT=-QKMS*(QVATM-QSG)
3091 RHCS=CAPICE(1)
3092 H=1.
3093 FKT=TKMS
3094 D1=cotso(NZS1)
3095 D2=rhtso(NZS1)
3096 TN=SOILT
3097 D9=THDIFICE(1)*RHCS*dzstop
3098 D10=TKMS*CP*RHO
3099 R211=.5*CONFLX/DELT
3100 R21=R211*CP*RHO
3101 R22=.5/(THDIFICE(1)*DELT*dzstop**2)
3102 R6=EMISS *STBOLT*.5*TN**4
3103 R7=R6/TN
3104 D11=RNET+R6
3106 IF(SNHEI.GE.SNTH) THEN
3107 if(snhei.le.DELTSN+SNTH) then
3108 !--- 1-layer snow
3109 D1SN = cotso(NZS)
3110 D2SN = rhtso(NZS)
3111 else
3112 !--- 2-layer snow
3113 D1SN = cotsn
3114 D2SN = rhtsn
3115 endif
3116 D9SN= THDIFSN*RHOCSN / SNPRIM
3117 R22SN = SNPRIM*SNPRIM*0.5/(THDIFSN*DELT)
3118 ENDIF
3120 IF(SNHEI.LT.SNTH.AND.SNHEI.GT.0.) then
3121 !--- thin snow is combined with sea ice
3122 D1SN = D1
3123 D2SN = D2
3124 D9SN = (fsn*THDIFSN*RHOCSN+fso*THDIFICE(1)*RHCS)/ &
3125 snprim
3126 R22SN = snprim*snprim*0.5 &
3127 /((fsn*THDIFSN+fso*THDIFICE(1))*delt)
3128 ENDIF
3130 IF(SNHEI.eq.0.)then
3131 !--- all snow is sublimated
3132 D9SN = D9
3133 R22SN = R22
3134 D1SN = D1
3135 D2SN = D2
3136 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
3137 print *,' SNHEI = 0, D9SN,R22SN,D1SN,D2SN: ',D9SN,R22SN,D1SN,D2SN
3138 ENDIF
3139 ENDIF
3141 !---- TDENOM for snow
3142 !18apr08 - the iteration start point
3143 212 continue
3144 TDENOM = D9SN*(1.-D1SN +R22SN)+D10+R21+R7 &
3145 +RAINF*CVW*PRCPMS &
3146 +RHOnewCSN*NEWSNOW/DELT
3148 FKQ=QKMS*RHO
3149 R210=R211*RHO
3150 AA=XLVM*(BETA*FKQ+R210)/TDENOM
3151 BB=(D10*TABS+R21*TN+XLVM*(QVATM* &
3152 (BETA*FKQ) &
3153 +R210*QVG)+D11+D9SN*(D2SN+R22SN*TN) &
3154 +RAINF*CVW*PRCPMS*max(273.15,TABS) &
3155 + RHOnewCSN*NEWSNOW/DELT*min(273.15,TABS) &
3156 !18apr08 - add heat of snow phase change
3157 -SNOH &
3158 )/TDENOM
3159 AA1=AA
3160 PP=PATM*1.E3
3161 AA1=AA1/PP
3162 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
3163 print *,'VILKA-SNOW on SEAICE'
3164 print *,'tn,aa1,bb,pp,fkq,r210', &
3165 tn,aa1,bb,pp,fkq,r210
3166 ENDIF
3168 CALL VILKA(TN,AA1,BB,PP,QS1,TS1,TBQ,KTAU,i,j,iland,isoil)
3169 !--- it is saturation over snow
3170 QVG=QS1
3171 QSG=QS1
3172 QCG=0.
3174 !--- SOILT - skin temperature
3175 SOILT=TS1
3177 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
3178 print *,' AFTER VILKA-SNOW on SEAICE'
3179 print *,' TS1,QS1: ', ts1,qs1
3180 ENDIF
3181 ! Solution for temperature at 7.5 cm depth and snow-seaice interface
3182 IF(SNHEI.GE.SNTH) THEN
3183 if(snhei.gt.DELTSN+SNTH) then
3184 !-- 2-layer snow model
3185 SOILT1=rhtsn+cotsn*SOILT
3186 TSO(1)=min(271.4,(rhtso(NZS)+cotso(NZS)*SOILT1))
3187 tsob=soilt1
3188 else
3189 !-- 1 layer in snow
3190 TSO(1)=min(271.4,(rhtso(NZS)+cotso(NZS)*SOILT))
3191 SOILT1=TSO(1)
3192 tsob=tso(1)
3193 endif
3194 ELSE
3195 TSO(1)=SOILT
3196 SOILT1=SOILT
3197 tsob=SOILT
3198 ENDIF
3199 !---- Final solution for TSO in sea ice
3200 DO K=2,NZS
3201 KK=NZS-K+1
3202 TSO(K)=min(271.4,(rhtso(KK)+cotso(KK)*TSO(K-1)))
3203 END DO
3204 !--- For thin snow layer combined with the top sea ice layer
3205 !--- TSO is computed by linear inmterpolation between SOILT
3206 !--- and TSO(2)
3208 if(nmelt.eq.1) go to 220
3210 !--- IF SOILT > 273.15 F then melting of snow can happen
3211 IF(SOILT.GT.273.15.AND.SNHEI.GT.0.) THEN
3212 nmelt = 1
3213 soiltfrac=snowfrac*273.15+(1.-snowfrac)*SOILT
3214 QSG= QSN(soiltfrac,TBQ)/PP
3215 QVG=QSG
3216 T3 = STBOLT*SOILTfrac*SOILTfrac*SOILTfrac
3217 UPFLUX = T3 * SOILTfrac
3218 XINET = EMISS*(GLW-UPFLUX)
3219 RNET = GSW + XINET
3220 EPOT = -QKMS*(QVATM-QSG)
3221 Q1=EPOT*RAS
3223 IF (Q1.LE.0.) THEN
3224 ! --- condensation
3225 DEW=-EPOT
3227 QFX= XLVM*RHO*DEW
3228 EETA=QFX/XLVM
3229 ELSE
3230 ! --- evaporation
3231 EETA = Q1 * BETA *1.E3
3232 ! to convert from kg m-2 s-1 to m s-1: 1/rho water=1.e-3************
3233 QFX= - XLVM * EETA
3234 ENDIF
3236 HFX=D10*(TABS-soiltfrac)
3238 IF(SNHEI.GE.SNTH)then
3239 SOH=thdifsn*RHOCSN*(soiltfrac-TSOB)/SNPRIM
3240 SNFLX=SOH
3241 ELSE
3242 SOH=(fsn*thdifsn*rhocsn+fso*thdifice(1)*rhcs)* &
3243 (soiltfrac-TSOB)/snprim
3244 SNFLX=SOH
3245 ENDIF
3246 X= (R21+D9SN*R22SN)*(soiltfrac-TNOLD) + &
3247 XLVM*R210*(QSG-QGOLD)
3248 !-- SNOH is energy flux of snow phase change
3249 SNOH=RNET+QFX +HFX &
3250 +RHOnewCSN*NEWSNOW/DELT*(min(273.15,TABS)-TN) &
3251 -SOH-X+RAINF*CVW*PRCPMS* &
3252 (max(273.15,TABS)-TN)
3253 SNOH=AMAX1(0.,SNOH)
3254 !-- SMELT is speed of melting in M/S
3255 SMELT= SNOH /XLMELT*1.E-3
3256 SMELT=AMIN1(SMELT,SNWEPR/DELT-BETA*EPOT*RAS)
3257 SMELT=AMAX1(0.,SMELT)
3259 !18apr08 - Egglston limit
3260 SMELT= amin1 (smelt, 5.6E-7*meltfactor*max(1.,(soilt-273.15)))
3261 ! SMELT= amin1 (smelt, 5.6E-8*meltfactor*max(1.,(soilt-273.15)))
3262 !*** From Koren et al. (1999) 13% of snow melt stays in the snow pack
3263 !!! rsm=0.13*smelt*delt
3264 if(snwepr.gt.0.) then
3265 rsmfrac=min(0.18,(max(0.08,0.10/snwepr*0.13)))
3266 ! else
3267 ! rsmfrac=0.13
3268 endif
3270 rsm=rsmfrac*smelt*delt
3271 !18apr08 rsm is part of melted water that stays in snow as liquid
3272 SMELT=AMAX1(0.,SMELT-rsm/delt)
3274 SNOHGNEW=SMELT*XLMELT*1.E3
3275 SNODIF=AMAX1(0.,(SNOH-SNOHGNEW))
3277 SNOH=SNOHGNEW
3280 !18apr08 - if snow melt occurred then go into iteration for energy budget
3281 ! solution
3282 !-- correction of liquid equivalent of snow depth
3283 !-- due to evaporation and snow melt
3284 SNWE = AMAX1(0.,(SNWEPR- &
3285 (SMELT+BETA*EPOT*RAS)*DELT &
3286 ) )
3288 !--- If all snow melts, then 13% of snow melt we kept in the
3289 !--- snow pack should be added back to snow melt and infiltrate
3290 !--- into soil.
3291 if(snwe.le.rsm) then
3292 smelt=smelt+rsm/delt
3293 snwe=0.
3294 rsm=0.
3295 else
3296 !*** Correct snow density on effect of snow melt, melted
3297 !*** from the top of the snow. 13% of melted water
3298 !*** remains in the pack and changes its density.
3299 !*** Eq. 9 (with my correction) in Koren et al. (1999)
3301 if(snwe.gt.0.) then
3302 xsn=(rhosn*(snwe-rsm)+1.e3*rsm)/ &
3303 snwe
3304 rhosn=MIN(XSN,400.)
3306 RHOCSN=2090.* RHOSN
3307 thdifsn = 0.265/RHOCSN
3308 endif
3310 endif
3312 !--- If there is no snow melting then just evaporation
3313 !--- or condensation cxhanges SNWE
3314 ELSE
3315 EPOT=-QKMS*(QVATM-QSG)
3316 SNWE = AMAX1(0.,(SNWEPR- &
3317 BETA*EPOT*RAS*DELT))
3319 ENDIF
3320 !*** Correct snow density on effect of snow melt, melted
3321 !*** from the top of the snow. 13% of melted water
3322 !*** remains in the pack and changes its density.
3323 !*** Eq. 9 (with my correction) in Koren et al. (1999)
3325 SNHEI=SNWE *1.E3 / RHOSN
3327 snweprint=snwe
3328 ! &
3329 !--- if VEGFRAC.ne.0. then some snow stays on the canopy
3330 !--- and should be added to SNWE for water conservation
3331 ! 4 Nov 07 +VEGFRAC*cst
3332 snheiprint=snweprint*1.E3 / RHOSN
3334 if(nmelt.eq.1) goto 212
3335 220 continue
3337 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
3338 print *, 'snweprint : ',snweprint
3339 print *, 'D9SN,SOILT,TSOB : ', D9SN,SOILT,TSOB
3340 ENDIF
3341 !--- Compute flux in the top snow layer
3342 SNFLX=D9SN*(SOILT-TSOB)
3343 IF(SNHEI.GT.0.) THEN
3344 if(ilnb.gt.1) then
3345 tsnav=0.5/snhei*((soilt+soilt1)*deltsn &
3346 +(soilt1+tso(1))*(SNHEI-DELTSN)) &
3347 -273.15
3348 else
3349 tsnav=0.5*(soilt+tso(1)) - 273.15
3350 endif
3351 ENDIF
3352 !--- RECALCULATION OF DEW USING NEW VALUE OF QSG
3353 DEW=0.
3354 PP=PATM*1.E3
3355 QSG= QSN(SOILT,TBQ)/PP
3356 EPOT = -FQ*(QVATM-QSG)
3357 IF(EPOT.LT.0.) THEN
3358 ! Sublimation
3359 DEW=-EPOT
3360 ENDIF
3361 !-- Restore sea-ice parameters if snow is less than threshold
3362 IF(SNHEI.EQ.0.) then
3363 tsnav=soilt-273.15
3364 smelt=smelt+snwe/delt
3365 rsm=0.
3366 emiss=1.
3367 znt=0.011
3368 alb=0.55
3369 ENDIF
3371 SNOM=SNOM+SMELT*DELT*1.e3
3373 !--- THE DIAGNOSTICS OF SURFACE FLUXES
3375 T3 = STBOLT*SOILT*SOILT*SOILT
3376 UPFLUX = T3 *SOILT
3377 XINET = EMISS*(GLW-UPFLUX)
3378 RNET = GSW + XINET
3379 HFT=-TKMS*CP*RHO*(TABS-SOILT) &
3380 *(P1000mb*0.00001/Patm)**ROVCP
3381 Q1 = - FQ*RAS* (QVATM - QSG)
3382 IF (Q1.LT.0.) THEN
3383 ! --- condensation
3384 !-- moisture flux for coupling with PBL
3385 EETA=-QKMS*RAS*(QVATM/(1.+QVATM) - QSG/(1.+QSG))*1.E3
3386 QFX= XLVm*EETA
3387 !-- actual moisture flux from RUC LSM
3388 DEW=QKMS*(QVATM-QSG)
3389 EETA= RHO*DEW
3390 sublim = EETA
3391 ELSE
3392 ! --- evaporation
3393 !-- moisture flux for coupling with PBL
3394 EETA=-QKMS*RAS*(QVATM/(1.+QVATM) - QVG/(1.+QVG))*1.E3
3395 ! EETA=-QKMS*RAS*(QVATM/(1.+QVATM) - QSG/(1.+QSG))*1.E3
3396 ! to convert from kg m-2 s-1 to m s-1: 1/rho water=1.e-3************
3397 QFX= XLVm * EETA
3398 !-- actual moisture flux from RUC LSM
3399 EETA = Q1*1.E3
3400 sublim = EETA
3401 ENDIF
3403 s=THDIFICE(1)*CAPICE(1)*dzstop*(tso(1)-tso(2))
3404 ! s=D9SN*(SOILT-TSOB)
3405 HFX=HFT
3406 FLTOT=RNET-HFT-QFX-S
3407 !------------------------------------------------------------------------
3408 !------------------------------------------------------------------------
3409 END SUBROUTINE SNOWSEAICE
3410 !------------------------------------------------------------------------
3413 SUBROUTINE SOILTEMP( &
3414 !--- input variables
3415 i,j,iland,isoil, &
3416 delt,ktau,conflx,nzs,nddzs,nroot, &
3417 PRCPMS,RAINF,PATM,TABS,QVATM,QCATM, &
3418 EMISS,RNET, &
3419 QKMS,TKMS,PC,RHO,VEGFRAC, &
3420 THDIF,CAP,DRYCAN,WETCAN, &
3421 TRANSUM,DEW,MAVAIL, &
3422 !--- soil fixed fields
3423 DQM,QMIN,BCLH, &
3424 ZSMAIN,ZSHALF,DTDZS,TBQ, &
3425 !--- constants
3426 XLV,CP,G0_P,CVW,STBOLT, &
3427 !--- output variables
3428 TSO,SOILT,QVG,QSG,QCG)
3430 !*************************************************************
3431 ! Energy budget equation and heat diffusion eqn are
3432 ! solved here and
3433 !
3434 ! DELT - time step (s)
3435 ! ktau - numver of time step
3436 ! CONFLX - depth of constant flux layer (m)
3437 ! IME, JME, KME, NZS - dimensions of the domain
3438 ! NROOT - number of levels within the root zone
3439 ! PRCPMS - precipitation rate in m/s
3440 ! COTSO, RHTSO - coefficients for implicit solution of
3441 ! heat diffusion equation
3442 ! THDIF - thermal diffusivity (m^2/s)
3443 ! QSG,QVG,QCG - saturated mixing ratio, mixing ratio of
3444 ! water vapor and cloud at the ground
3445 ! surface, respectively (kg/kg)
3446 ! PATM - pressure [bar]
3447 ! QC3D,QV3D - cloud and water vapor mixing ratio
3448 ! at the first atm. level (kg/kg)
3449 ! EMISS,RNET - emissivity (0-1) of the ground surface and net
3450 ! radiation at the surface (W/m^2)
3451 ! QKMS - exchange coefficient for water vapor in the
3452 ! surface layer (m/s)
3453 ! TKMS - exchange coefficient for heat in the surface
3454 ! layer (m/s)
3455 ! PC - plant coefficient (resistance)
3456 ! RHO - density of atmosphere near surface (kg/m^3)
3457 ! VEGFRAC - greeness fraction (0-1)
3458 ! CAP - volumetric heat capacity (J/m^3/K)
3459 ! DRYCAN - dry fraction of vegetated area where
3460 ! transpiration may take place (0-1)
3461 ! WETCAN - fraction of vegetated area covered by canopy
3462 ! water (0-1)
3463 ! TRANSUM - transpiration function integrated over the
3464 ! rooting zone (m)
3465 ! DEW - dew in kg/m^2s
3466 ! MAVAIL - fraction of maximum soil moisture in the top
3467 ! layer (0-1)
3468 ! ZSMAIN - main levels in soil (m)
3469 ! ZSHALF - middle of the soil layers (m)
3470 ! DTDZS - dt/(2.*dzshalf*dzmain)
3471 ! TBQ - table to define saturated mixing ration
3472 ! of water vapor for given temperature and pressure
3473 ! TSO - soil temperature (K)
3474 ! SOILT - skin temperature (K)
3475 !
3476 !****************************************************************
3478 IMPLICIT NONE
3479 !-----------------------------------------------------------------
3481 !--- input variables
3483 INTEGER, INTENT(IN ) :: nroot,ktau,nzs , &
3484 nddzs !nddzs=2*(nzs-2)
3485 INTEGER, INTENT(IN ) :: i,j,iland,isoil
3486 REAL, INTENT(IN ) :: DELT,CONFLX,PRCPMS, RAINF
3487 REAL, INTENT(INOUT) :: DRYCAN,WETCAN,TRANSUM
3488 !--- 3-D Atmospheric variables
3489 REAL, &
3490 INTENT(IN ) :: PATM, &
3491 QVATM, &
3492 QCATM
3493 !--- 2-D variables
3494 REAL , &
3495 INTENT(IN ) :: &
3496 EMISS, &
3497 RHO, &
3498 RNET, &
3499 PC, &
3500 VEGFRAC, &
3501 DEW, &
3502 QKMS, &
3503 TKMS
3505 !--- soil properties
3506 REAL , &
3507 INTENT(IN ) :: &
3508 BCLH, &
3509 DQM, &
3510 QMIN
3512 REAL, INTENT(IN ) :: CP, &
3513 CVW, &
3514 XLV, &
3515 STBOLT, &
3516 TABS, &
3517 G0_P
3520 REAL, DIMENSION(1:NZS), INTENT(IN) :: ZSMAIN, &
3521 ZSHALF, &
3522 THDIF, &
3523 CAP
3525 REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS
3527 REAL, DIMENSION(1:4001), INTENT(IN) :: TBQ
3530 !--- input/output variables
3531 !-------- 3-d soil moisture and temperature
3532 REAL, DIMENSION( 1:nzs ) , &
3533 INTENT(INOUT) :: TSO
3535 !-------- 2-d variables
3536 REAL , &
3537 INTENT(INOUT) :: &
3538 MAVAIL, &
3539 QVG, &
3540 QSG, &
3541 QCG, &
3542 SOILT
3545 !--- Local variables
3547 REAL :: x,x1,x2,x4,dzstop,can,ft,sph , &
3548 tn,trans,umveg,denom
3550 REAL :: FKT,D1,D2,D9,D10,DID,R211,R21,R22,R6,R7,D11 , &
3551 PI,H,FKQ,R210,AA,BB,PP,Q1,QS1,TS1,TQ2,TX2 , &
3552 TDENOM
3554 REAL :: C,CC,AA1,RHCS,H1
3556 REAL, DIMENSION(1:NZS) :: cotso,rhtso
3558 INTEGER :: nzs1,nzs2,k,k1,kn,kk
3560 !-----------------------------------------------------------------
3563 NZS1=NZS-1
3564 NZS2=NZS-2
3565 dzstop=1./(ZSMAIN(2)-ZSMAIN(1))
3567 do k=1,nzs
3568 cotso(k)=0.
3569 rhtso(k)=0.
3570 enddo
3571 !******************************************************************************
3572 ! COEFFICIENTS FOR THOMAS ALGORITHM FOR TSO
3573 !******************************************************************************
3574 ! did=2.*(ZSMAIN(nzs)-ZSHALF(nzs))
3575 ! h1=DTDZS(8)*THDIF(nzs-1)*(ZSHALF(nzs)-ZSHALF(nzs-1))/did
3576 ! cotso(1)=h1/(1.+h1)
3577 ! rhtso(1)=(tso(nzs)+h1*(tso(nzs-1)-tso(nzs)))/
3578 ! 1 (1.+h1)
3579 cotso(1)=0.
3580 rhtso(1)=TSO(NZS)
3581 DO 33 K=1,NZS2
3582 KN=NZS-K
3583 K1=2*KN-3
3584 X1=DTDZS(K1)*THDIF(KN-1)
3585 X2=DTDZS(K1+1)*THDIF(KN)
3586 FT=TSO(KN)+X1*(TSO(KN-1)-TSO(KN)) &
3587 -X2*(TSO(KN)-TSO(KN+1))
3588 DENOM=1.+X1+X2-X2*cotso(K)
3589 cotso(K+1)=X1/DENOM
3590 rhtso(K+1)=(FT+X2*rhtso(K))/DENOM
3591 33 CONTINUE
3593 !************************************************************************
3594 !--- THE HEAT BALANCE EQUATION (Smirnova et al., 1996, EQ. 21,26)
3596 RHCS=CAP(1)
3597 H=MAVAIL
3598 IF(DEW.NE.0.)THEN
3599 DRYCAN=0.
3600 WETCAN=1.
3601 ENDIf
3602 TRANS=PC*TRANSUM*DRYCAN/ZSHALF(NROOT+1)
3603 CAN=WETCAN+TRANS
3604 UMVEG=1.-VEGFRAC
3605 FKT=TKMS
3606 D1=cotso(NZS1)
3607 D2=rhtso(NZS1)
3608 TN=SOILT
3609 D9=THDIF(1)*RHCS*dzstop
3610 D10=TKMS*CP*RHO
3611 R211=.5*CONFLX/DELT
3612 R21=R211*CP*RHO
3613 R22=.5/(THDIF(1)*DELT*dzstop**2)
3614 R6=EMISS *STBOLT*.5*TN**4
3615 R7=R6/TN
3616 D11=RNET+R6
3617 TDENOM=D9*(1.-D1+R22)+D10+R21+R7 &
3618 +RAINF*CVW*PRCPMS
3619 FKQ=QKMS*RHO
3620 R210=R211*RHO
3621 C=VEGFRAC*FKQ*CAN
3622 CC=C*XLV/TDENOM
3623 AA=XLV*(FKQ*UMVEG+R210)/TDENOM
3624 BB=(D10*TABS+R21*TN+XLV*(QVATM* &
3625 (FKQ*UMVEG+C) &
3626 +R210*QVG)+D11+D9*(D2+R22*TN) &
3627 +RAINF*CVW*PRCPMS*max(273.15,TABS) &
3628 )/TDENOM
3629 AA1=AA+CC
3630 PP=PATM*1.E3
3631 AA1=AA1/PP
3632 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
3633 PRINT *,' VILKA-1'
3634 print *,'D10,TABS,R21,TN,QVATM,FKQ,UMVEG,VEGFRAC,CAN', &
3635 D10,TABS,R21,TN,QVATM,FKQ,UMVEG,VEGFRAC,CAN
3636 print *,'RNET, EMISS, STBOLT, SOILT',RNET, EMISS, STBOLT, SOILT
3637 print *,'R210,QVG,D11,D9,D2,R22,RAINF,CVW,PRCPMS,TDENOM', &
3638 R210,QVG,D11,D9,D2,R22,RAINF,CVW,PRCPMS,TDENOM
3639 print *,'tn,aa1,bb,pp,umveg,fkq,r210,vegfrac', &
3640 tn,aa1,bb,pp,umveg,fkq,r210,vegfrac
3641 ENDIF
3642 CALL VILKA(TN,AA1,BB,PP,QS1,TS1,TBQ,KTAU,i,j,iland,isoil)
3643 TQ2=QVATM
3644 TX2=TQ2*(1.-H)
3645 Q1=TX2+H*QS1
3646 IF(Q1.LT.QS1) GOTO 100
3647 !--- if no saturation - goto 100
3648 !--- if saturation - goto 90
3649 90 QVG=QS1
3650 QSG=QS1
3651 TSO(1)=TS1
3652 QCG=max(0.,Q1-QS1)
3653 GOTO 200
3654 100 BB=BB-AA*TX2
3655 AA=(AA*H+CC)/PP
3656 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
3657 PRINT *,' VILKA-2'
3658 print *,'D10,TABS,R21,TN,QVATM,FKQ,UMVEG,VEGFRAC,CAN', &
3659 D10,TABS,R21,TN,QVATM,FKQ,UMVEG,VEGFRAC,CAN
3660 print *,'R210,QVG,D11,D9,D2,R22,RAINF,CVW,PRCPMS,TDENOM', &
3661 R210,QVG,D11,D9,D2,R22,RAINF,CVW,PRCPMS,TDENOM
3663 print *,'tn,aa1,bb,pp,umveg,fkq,r210,vegfrac', &
3664 tn,aa1,bb,pp,umveg,fkq,r210,vegfrac
3665 ENDIF
3667 CALL VILKA(TN,AA,BB,PP,QS1,TS1,TBQ,KTAU,i,j,iland,isoil)
3668 Q1=TX2+H*QS1
3669 IF(Q1.GT.QS1) GOTO 90
3670 QSG=QS1
3671 QVG=Q1
3672 TSO(1)=TS1
3673 QCG=0.
3674 200 CONTINUE
3676 !--- SOILT - skin temperature
3677 SOILT=TS1
3679 !---- Final solution for soil temperature - TSO
3680 DO K=2,NZS
3681 KK=NZS-K+1
3682 TSO(K)=rhtso(KK)+cotso(KK)*TSO(K-1)
3683 END DO
3685 !--------------------------------------------------------------------
3686 END SUBROUTINE SOILTEMP
3687 !--------------------------------------------------------------------
3690 SUBROUTINE SNOWTEMP( &
3691 !--- input variables
3692 i,j,iland,isoil, &
3693 delt,ktau,conflx,nzs,nddzs,nroot, &
3694 snwe,snwepr,snhei,newsnow,snowfrac, &
3695 beta,deltsn,snth,rhosn,rhonewsn,meltfactor, & ! add meltfactor
3696 PRCPMS,RAINF, &
3697 PATM,TABS,QVATM,QCATM, &
3698 GLW,GSW,EMISS,RNET, &
3699 QKMS,TKMS,PC,RHO,VEGFRAC, &
3700 THDIF,CAP,DRYCAN,WETCAN,CST, &
3701 TRANF,TRANSUM,DEW,MAVAIL, &
3702 !--- soil fixed fields
3703 DQM,QMIN,PSIS,BCLH, &
3704 ZSMAIN,ZSHALF,DTDZS,TBQ, &
3705 !--- constants
3706 XLVM,CP,rovcp,G0_P,CVW,STBOLT, &
3707 !--- output variables
3708 SNWEPRINT,SNHEIPRINT,RSM, &
3709 TSO,SOILT,SOILT1,TSNAV,QVG,QSG,QCG, &
3710 SMELT,SNOH,SNFLX,ILNB)
3712 !********************************************************************
3713 ! Energy budget equation and heat diffusion eqn are
3714 ! solved here to obtain snow and soil temperatures
3715 !
3716 ! DELT - time step (s)
3717 ! ktau - numver of time step
3718 ! CONFLX - depth of constant flux layer (m)
3719 ! IME, JME, KME, NZS - dimensions of the domain
3720 ! NROOT - number of levels within the root zone
3721 ! PRCPMS - precipitation rate in m/s
3722 ! COTSO, RHTSO - coefficients for implicit solution of
3723 ! heat diffusion equation
3724 ! THDIF - thermal diffusivity (W/m/K)
3725 ! QSG,QVG,QCG - saturated mixing ratio, mixing ratio of
3726 ! water vapor and cloud at the ground
3727 ! surface, respectively (kg/kg)
3728 ! PATM - pressure [bar]
3729 ! QCATM,QVATM - cloud and water vapor mixing ratio
3730 ! at the first atm. level (kg/kg)
3731 ! EMISS,RNET - emissivity (0-1) of the ground surface and net
3732 ! radiation at the surface (W/m^2)
3733 ! QKMS - exchange coefficient for water vapor in the
3734 ! surface layer (m/s)
3735 ! TKMS - exchange coefficient for heat in the surface
3736 ! layer (m/s)
3737 ! PC - plant coefficient (resistance)
3738 ! RHO - density of atmosphere near surface (kg/m^3)
3739 ! VEGFRAC - greeness fraction (0-1)
3740 ! CAP - volumetric heat capacity (J/m^3/K)
3741 ! DRYCAN - dry fraction of vegetated area where
3742 ! transpiration may take place (0-1)
3743 ! WETCAN - fraction of vegetated area covered by canopy
3744 ! water (0-1)
3745 ! TRANSUM - transpiration function integrated over the
3746 ! rooting zone (m)
3747 ! DEW - dew in kg/m^2/s
3748 ! MAVAIL - fraction of maximum soil moisture in the top
3749 ! layer (0-1)
3750 ! ZSMAIN - main levels in soil (m)
3751 ! ZSHALF - middle of the soil layers (m)
3752 ! DTDZS - dt/(2.*dzshalf*dzmain)
3753 ! TBQ - table to define saturated mixing ration
3754 ! of water vapor for given temperature and pressure
3755 ! TSO - soil temperature (K)
3756 ! SOILT - skin temperature (K)
3757 !
3758 !*********************************************************************
3760 IMPLICIT NONE
3761 !---------------------------------------------------------------------
3762 !--- input variables
3764 INTEGER, INTENT(IN ) :: nroot,ktau,nzs , &
3765 nddzs !nddzs=2*(nzs-2)
3767 INTEGER, INTENT(IN ) :: i,j,iland,isoil
3768 REAL, INTENT(IN ) :: DELT,CONFLX,PRCPMS , &
3769 RAINF,NEWSNOW,DELTSN,SNTH , &
3770 TABS,TRANSUM,SNWEPR , &
3771 rhonewsn,meltfactor
3772 real :: rhonewcsn
3774 !--- 3-D Atmospheric variables
3775 REAL, &
3776 INTENT(IN ) :: PATM, &
3777 QVATM, &
3778 QCATM
3779 !--- 2-D variables
3780 REAL , &
3781 INTENT(IN ) :: GLW, &
3782 GSW, &
3783 RHO, &
3784 PC, &
3785 VEGFRAC, &
3786 QKMS, &
3787 TKMS
3789 !--- soil properties
3790 REAL , &
3791 INTENT(IN ) :: &
3792 BCLH, &
3793 DQM, &
3794 PSIS, &
3795 QMIN
3797 REAL, INTENT(IN ) :: CP, &
3798 ROVCP, &
3799 CVW, &
3800 STBOLT, &
3801 XLVM, &
3802 G0_P
3805 REAL, DIMENSION(1:NZS), INTENT(IN) :: ZSMAIN, &
3806 ZSHALF, &
3807 THDIF, &
3808 CAP, &
3809 TRANF
3811 REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS
3813 REAL, DIMENSION(1:4001), INTENT(IN) :: TBQ
3816 !--- input/output variables
3817 !-------- 3-d soil moisture and temperature
3818 REAL, DIMENSION( 1:nzs ) , &
3819 INTENT(INOUT) :: TSO
3822 !-------- 2-d variables
3823 REAL , &
3824 INTENT(INOUT) :: DEW, &
3825 CST, &
3826 RHOSN, &
3827 EMISS, &
3828 MAVAIL, &
3829 QVG, &
3830 QSG, &
3831 QCG, &
3832 SNWE, &
3833 SNHEI, &
3834 SNOWFRAC, &
3835 SMELT, &
3836 SNOH, &
3837 SNFLX, &
3838 SOILT, &
3839 SOILT1, &
3840 TSNAV
3842 REAL, INTENT(INOUT) :: DRYCAN, WETCAN
3844 REAL, INTENT(OUT) :: RSM, &
3845 SNWEPRINT, &
3846 SNHEIPRINT
3847 INTEGER, INTENT(OUT) :: ilnb
3848 !--- Local variables
3851 INTEGER :: nzs1,nzs2,k,k1,kn,kk
3853 REAL :: x,x1,x2,x4,dzstop,can,ft,sph, &
3854 tn,trans,umveg,denom
3856 REAL :: cotsn,rhtsn,xsn1,ddzsn1,x1sn1,ftsnow,denomsn
3858 REAL :: t3,upflux,xinet,ras, &
3859 xlmelt,rhocsn,thdifsn, &
3860 beta,epot,xsn,ddzsn,x1sn,d1sn,d2sn,d9sn,r22sn
3862 REAL :: fso,fsn, &
3863 FKT,D1,D2,D9,D10,DID,R211,R21,R22,R6,R7,D11, &
3864 PI,H,FKQ,R210,AA,BB,PP,Q1,QS1,TS1,TQ2,TX2, &
3865 TDENOM,C,CC,AA1,RHCS,H1, &
3866 tsob, snprim, sh1, sh2, &
3867 smeltg,snohg,snodif,soh, &
3868 CMC2MS,TNOLD,QGOLD,SNOHGNEW
3870 REAL, DIMENSION(1:NZS) :: transp,cotso,rhtso
3871 REAL :: edir1, &
3872 ec1, &
3873 ett1, &
3874 eeta, &
3875 s, &
3876 qfx, &
3877 hfx
3879 REAL :: RNET,rsmfrac,soiltfrac,hsn
3880 integer :: nmelt
3882 !-----------------------------------------------------------------
3884 do k=1,nzs
3885 transp (k)=0.
3886 cotso (k)=0.
3887 rhtso (k)=0.
3888 enddo
3890 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
3891 print *, 'SNOWTEMP: SNHEI,SNTH,SOILT1: ',SNHEI,SNTH,SOILT1,soilt
3892 ENDIF
3893 XLMELT=3.35E+5
3894 RHOCSN=2090.* RHOSN
3895 !18apr08 - add rhonewcsn
3896 RHOnewCSN=2090.* RHOnewSN
3897 THDIFSN = 0.265/RHOCSN
3898 RAS=RHO*1.E-3
3900 SOILTFRAC=SOILT
3902 SMELT=0.
3903 SOH=0.
3904 SMELTG=0.
3905 SNOHG=0.
3906 SNODIF=0.
3907 RSM = 0.
3908 RSMFRAC = 0.
3909 fsn=1.
3910 fso=0.
3911 hsn=snhei
3913 NZS1=NZS-1
3914 NZS2=NZS-2
3916 QGOLD=QVG
3917 TNOLD=SOILT
3918 DZSTOP=1./(ZSMAIN(2)-ZSMAIN(1))
3920 !******************************************************************************
3921 ! COEFFICIENTS FOR THOMAS ALGORITHM FOR TSO
3922 !******************************************************************************
3923 ! did=2.*(ZSMAIN(nzs)-ZSHALF(nzs))
3924 ! h1=DTDZS(8)*THDIF(nzs-1)*(ZSHALF(nzs)-ZSHALF(nzs-1))/did
3925 ! cotso(1)=h1/(1.+h1)
3926 ! rhtso(1)=(tso(nzs)+h1*(tso(nzs-1)-tso(nzs)))/
3927 ! 1 (1.+h1)
3929 cotso(1)=0.
3930 rhtso(1)=TSO(NZS)
3931 DO 33 K=1,NZS2
3932 KN=NZS-K
3933 K1=2*KN-3
3934 X1=DTDZS(K1)*THDIF(KN-1)
3935 X2=DTDZS(K1+1)*THDIF(KN)
3936 FT=TSO(KN)+X1*(TSO(KN-1)-TSO(KN)) &
3937 -X2*(TSO(KN)-TSO(KN+1))
3938 DENOM=1.+X1+X2-X2*cotso(K)
3939 cotso(K+1)=X1/DENOM
3940 rhtso(K+1)=(FT+X2*rhtso(K))/DENOM
3941 33 CONTINUE
3942 !--- THE NZS element in COTSO and RHTSO will be for snow
3943 !--- There will be 2 layers in snow if it is deeper than DELTSN+SNTH
3944 IF(SNHEI.GE.SNTH) then
3945 if(snhei.le.DELTSN+SNTH) then
3946 !-- 1-layer snow model
3947 ilnb=1
3948 snprim=snhei
3949 tsob=tso(1)
3950 XSN = DELT/2./(zshalf(2)+0.5*SNPRIM)
3951 DDZSN = XSN / SNPRIM
3952 X1SN = DDZSN * thdifsn
3953 X2 = DTDZS(1)*THDIF(1)
3954 FT = TSO(1)+X1SN*(SOILT-TSO(1)) &
3955 -X2*(TSO(1)-TSO(2))
3956 DENOM = 1. + X1SN + X2 -X2*cotso(NZS1)
3957 cotso(NZS)=X1SN/DENOM
3958 rhtso(NZS)=(FT+X2*rhtso(NZS1))/DENOM
3959 cotsn=cotso(NZS)
3960 rhtsn=rhtso(NZS)
3961 !*** Average temperature of snow pack (C)
3962 tsnav=0.5*(soilt+tso(1)) &
3963 -273.15
3965 else
3966 !-- 2 layers in snow, SOILT1 is temperasture at DELTSN depth
3967 ilnb=2
3968 snprim=deltsn
3969 tsob=soilt1
3970 XSN = DELT/2./(0.5*SNPRIM)
3971 XSN1= DELT/2./(zshalf(2)+0.5*(SNHEI-DELTSN))
3972 DDZSN = XSN / DELTSN
3973 DDZSN1 = XSN1 / (SNHEI-DELTSN)
3974 X1SN = DDZSN * thdifsn
3975 X1SN1 = DDZSN1 * thdifsn
3976 X2 = DTDZS(1)*THDIF(1)
3977 FT = TSO(1)+X1SN1*(SOILT1-TSO(1)) &
3978 -X2*(TSO(1)-TSO(2))
3979 DENOM = 1. + X1SN1 + X2 - X2*cotso(NZS1)
3980 cotso(nzs)=x1sn1/denom
3981 rhtso(nzs)=(ft+x2*rhtso(nzs1))/denom
3982 ftsnow = soilt1+x1sn*(soilt-soilt1) &
3983 -x1sn1*(soilt1-tso(1))
3984 denomsn = 1. + X1SN + X1SN1 - X1SN1*cotso(NZS)
3985 cotsn=x1sn/denomsn
3986 rhtsn=(ftsnow+X1SN1*rhtso(NZS))/denomsn
3987 !*** Average temperature of snow pack (C)
3988 tsnav=0.5/snhei*((soilt+soilt1)*deltsn &
3989 +(soilt1+tso(1))*(SNHEI-DELTSN)) &
3990 -273.15
3991 endif
3992 ENDIF
3994 IF(SNHEI.LT.SNTH.AND.SNHEI.GT.0.) then
3995 !--- snow is too thin to be treated separately, therefore it
3996 !--- is combined with the first soil layer.
3997 fsn=SNHEI/(SNHEI+zsmain(2))
3998 fso=1.-fsn
3999 soilt1=tso(1)
4000 tsob=tso(2)
4001 snprim=SNHEI+zsmain(2)
4002 XSN = DELT/2./((zshalf(3)-zsmain(2))+0.5*snprim)
4003 DDZSN = XSN /snprim
4004 X1SN = DDZSN * (fsn*thdifsn+fso*thdif(1))
4005 X2=DTDZS(2)*THDIF(2)
4006 FT=TSO(2)+X1SN*(SOILT-TSO(2))- &
4007 X2*(TSO(2)-TSO(3))
4008 denom = 1. + x1sn + x2 - x2*cotso(nzs-2)
4009 cotso(nzs1) = x1sn/denom
4010 rhtso(nzs1)=(FT+X2*rhtso(NZS-2))/denom
4011 tsnav=0.5*(soilt+tso(1)) &
4012 -273.15
4013 ENDIF
4015 !************************************************************************
4016 !--- THE HEAT BALANCE EQUATION (Smirnova et al. 1996, EQ. 21,26)
4017 !18apr08 nmelt is the flag for melting, and SNOH is heat of snow phase changes
4018 nmelt=0
4019 SNOH=0.
4022 ETT1=0.
4023 EPOT=-QKMS*(QVATM-QSG)
4024 RHCS=CAP(1)
4025 H=MAVAIL
4026 IF(DEW.NE.0.)THEN
4027 DRYCAN=0.
4028 WETCAN=1.
4029 ENDIF
4030 TRANS=PC*TRANSUM*DRYCAN/ZSHALF(NROOT+1)
4031 CAN=WETCAN+TRANS
4032 UMVEG=1.-VEGFRAC
4033 FKT=TKMS
4034 D1=cotso(NZS1)
4035 D2=rhtso(NZS1)
4036 TN=SOILT
4037 D9=THDIF(1)*RHCS*dzstop
4038 D10=TKMS*CP*RHO
4039 R211=.5*CONFLX/DELT
4040 R21=R211*CP*RHO
4041 R22=.5/(THDIF(1)*DELT*dzstop**2)
4042 R6=EMISS *STBOLT*.5*TN**4
4043 R7=R6/TN
4044 D11=RNET+R6
4046 IF(SNHEI.GE.SNTH) THEN
4047 if(snhei.le.DELTSN+SNTH) then
4048 !--- 1-layer snow
4049 D1SN = cotso(NZS)
4050 D2SN = rhtso(NZS)
4051 else
4052 !--- 2-layer snow
4053 D1SN = cotsn
4054 D2SN = rhtsn
4055 endif
4056 D9SN= THDIFSN*RHOCSN / SNPRIM
4057 R22SN = SNPRIM*SNPRIM*0.5/(THDIFSN*DELT)
4058 ENDIF
4060 IF(SNHEI.LT.SNTH.AND.SNHEI.GT.0.) then
4061 !--- thin snow is combined with soil
4062 D1SN = D1
4063 D2SN = D2
4064 D9SN = (fsn*THDIFSN*RHOCSN+fso*THDIF(1)*RHCS)/ &
4065 snprim
4066 R22SN = snprim*snprim*0.5 &
4067 /((fsn*THDIFSN+fso*THDIF(1))*delt)
4068 ENDIF
4070 IF(SNHEI.eq.0.)then
4071 !--- all snow is sublimated
4072 D9SN = D9
4073 R22SN = R22
4074 D1SN = D1
4075 D2SN = D2
4076 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
4077 print *,' SNHEI = 0, D9SN,R22SN,D1SN,D2SN: ',D9SN,R22SN,D1SN,D2SN
4078 ENDIF
4079 ENDIF
4081 !---- TDENOM for snow
4082 !18apr08 - the iteration start point
4083 212 continue
4084 TDENOM = D9SN*(1.-D1SN +R22SN)+D10+R21+R7 &
4085 +RAINF*CVW*PRCPMS &
4086 +RHOnewCSN*NEWSNOW/DELT
4088 FKQ=QKMS*RHO
4089 R210=R211*RHO
4090 C=VEGFRAC*FKQ*CAN
4091 CC=C*XLVM/TDENOM
4092 AA=XLVM*(BETA*FKQ*UMVEG+R210)/TDENOM
4093 BB=(D10*TABS+R21*TN+XLVM*(QVATM* &
4094 (BETA*FKQ*UMVEG+C) &
4095 +R210*QVG)+D11+D9SN*(D2SN+R22SN*TN) &
4096 +RAINF*CVW*PRCPMS*max(273.15,TABS) &
4097 + RHOnewCSN*NEWSNOW/DELT*min(273.15,TABS) &
4098 !18apr08 - added heat of snow phase change computed in the first iteration
4099 -SNOH &
4100 )/TDENOM
4101 AA1=AA+CC
4102 PP=PATM*1.E3
4103 AA1=AA1/PP
4104 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
4105 print *,'VILKA-SNOW'
4106 print *,'tn,aa1,bb,pp,umveg,fkq,r210,vegfrac', &
4107 tn,aa1,bb,pp,umveg,fkq,r210,vegfrac
4108 ENDIF
4110 CALL VILKA(TN,AA1,BB,PP,QS1,TS1,TBQ,KTAU,i,j,iland,isoil)
4111 !--- it is saturation over snow
4112 QVG=QS1
4113 QSG=QS1
4114 QCG=0.
4115 !--- SOILT - skin temperature
4116 SOILT=TS1
4118 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
4119 print *,' AFTER VILKA-SNOW'
4120 print *,' TS1,QS1: ', ts1,qs1
4121 ENDIF
4123 ! Solution for temperature at 7.5 cm depth and snow-soil interface
4124 IF(SNHEI.GE.SNTH) THEN
4125 if(snhei.gt.DELTSN+SNTH) then
4126 !-- 2-layer snow model
4127 SOILT1=rhtsn+cotsn*SOILT
4128 TSO(1)=rhtso(NZS)+cotso(NZS)*SOILT1
4129 tsob=soilt1
4130 else
4131 !-- 1 layer in snow
4132 TSO(1)=rhtso(NZS)+cotso(NZS)*SOILT
4133 SOILT1=TSO(1)
4134 tsob=tso(1)
4135 endif
4136 ELSE
4137 TSO(1)=SOILT
4138 SOILT1=SOILT
4139 tsob=SOILT
4140 ENDIF
4142 !---- Final solution for TSO
4143 DO K=2,NZS
4144 KK=NZS-K+1
4145 TSO(K)=rhtso(KK)+cotso(KK)*TSO(K-1)
4146 END DO
4147 !--- For thin snow layer combined with the top soil layer
4148 !--- TSO is computed by linear inmterpolation between SOILT
4149 !--- and TSO(2)
4151 if(SNHEI.LT.SNTH.AND.SNHEI.GT.0.)then
4152 tso(1)=tso(2)+(soilt-tso(2))*fso
4153 SOILT1=TSO(1)
4154 tsob=tso(2)
4155 !!! tsob=tso(1)
4156 endif
4158 if(nmelt.eq.1) go to 220
4160 !--- IF SOILT > 273.15 F then melting of snow can happen
4161 IF(SOILT.GT.273.15.AND.SNHEI.GT.0.) THEN
4162 nmelt = 1
4163 soiltfrac=snowfrac*273.15+(1.-snowfrac)*SOILT
4164 QSG= QSN(soiltfrac,TBQ)/PP
4165 QVG=QSG
4166 T3 = STBOLT*SOILTfrac*SOILTfrac*SOILTfrac
4167 UPFLUX = T3 * SOILTfrac
4168 XINET = EMISS*(GLW-UPFLUX)
4169 RNET = GSW + XINET
4170 EPOT = -QKMS*(QVATM-QSG)
4171 Q1=EPOT*RAS
4173 IF (Q1.LE.0.) THEN
4174 ! --- condensation
4175 DEW=-EPOT
4176 DO K=1,NZS
4177 TRANSP(K)=0.
4178 ENDDO
4180 QFX= XLVM*RHO*DEW
4181 EETA=QFX/XLVM
4182 ELSE
4183 ! --- evaporation
4184 DO K=1,NROOT
4185 TRANSP(K)=-VEGFRAC*q1 &
4186 *PC*TRANF(K)*DRYCAN/zshalf(NROOT+1)
4187 IF(TRANSP(K).GT.0.) TRANSP(K)=0.
4188 ETT1=ETT1-TRANSP(K)
4189 ENDDO
4190 DO k=nroot+1,nzs
4191 transp(k)=0.
4192 enddo
4194 EDIR1 = Q1*UMVEG * BETA
4195 EC1 = Q1 * WETCAN *VEGFRAC
4196 CMC2MS=CST/DELT
4197 EC1=MIN(CMC2MS,EC1)
4198 EETA = (EDIR1 + EC1 + ETT1)*1.E3
4199 ! to convert from kg m-2 s-1 to m s-1: 1/rho water=1.e-3************
4200 QFX= - XLVM * EETA
4201 ENDIF
4203 HFX=D10*(TABS-soiltfrac)
4205 IF(SNHEI.GE.SNTH)then
4206 SOH=thdifsn*RHOCSN*(soiltfrac-TSOB)/SNPRIM
4207 SNFLX=SOH
4208 ELSE
4209 SOH=(fsn*thdifsn*rhocsn+fso*thdif(1)*rhcs)* &
4210 (soiltfrac-TSOB)/snprim
4211 SNFLX=SOH
4212 ENDIF
4214 X= (R21+D9SN*R22SN)*(soiltfrac-TNOLD) + &
4215 XLVM*R210*(QSG-QGOLD)
4216 !-- SNOH is energy flux of snow phase change
4217 SNOH=RNET+QFX +HFX &
4218 +RHOnewCSN*NEWSNOW/DELT*(min(273.15,TABS)-TN) &
4219 -SOH-X+RAINF*CVW*PRCPMS* &
4220 (max(273.15,TABS)-TN)
4221 SNOH=AMAX1(0.,SNOH)
4222 !-- SMELT is speed of melting in M/S
4223 SMELT= SNOH /XLMELT*1.E-3
4224 ! SMELT=AMIN1(SMELT,SNWEPR/DELT-BETA*EPOT*RAS*UMVEG)
4225 SMELT=AMIN1(SMELT,SNWEPR/DELT-BETA*EPOT*RAS)
4226 SMELT=AMAX1(0.,SMELT)
4228 !18apr08 - Egglston limit
4229 SMELT= amin1 (smelt, 5.6E-7*meltfactor*max(1.,(soilt-273.15)))
4230 ! SMELT= amin1 (smelt, 5.6E-8*meltfactor*max(1.,(soilt-273.15)))
4232 !*** From Koren et al. (1999) 13% of snow melt stays in the snow pack
4233 !!! rsm=0.13*smelt*delt
4234 if(snwepr.gt.0.) then
4235 rsmfrac=min(0.18,(max(0.08,0.10/snwepr*0.13)))
4236 ! else
4237 ! rsmfrac=0.13
4238 endif
4240 rsm=rsmfrac*smelt*delt
4241 !18apr08 rsm is part of melted water that stays in snow as liquid
4242 SMELT=AMAX1(0.,SMELT-rsm/delt)
4244 SNOHGNEW=SMELT*XLMELT*1.E3
4245 SNODIF=AMAX1(0.,(SNOH-SNOHGNEW))
4247 SNOH=SNOHGNEW
4249 !-- correction of liquid equivalent of snow depth
4250 !-- due to evaporation and snow melt
4251 SNWE = AMAX1(0.,(SNWEPR- &
4252 (SMELT+BETA*EPOT*RAS)*DELT &
4253 ! (SMELT+BETA*EPOT*RAS*UMVEG)*DELT &
4254 ) )
4256 !--- If all snow melts, then 13% of snow melt we kept in the
4257 !--- snow pack should be added back to snow melt and infiltrate
4258 !--- into soil.
4259 if(snwe.le.rsm) then
4260 smelt=smelt+rsm/delt
4261 snwe=0.
4262 rsm=0.
4263 else
4264 !*** Correct snow density on effect of snow melt, melted
4265 !*** from the top of the snow. 13% of melted water
4266 !*** remains in the pack and changes its density.
4267 !*** Eq. 9 (with my correction) in Koren et al. (1999)
4268 if(snwe.gt.0.) then
4269 xsn=(rhosn*(snwe-rsm)+1.e3*rsm)/ &
4270 snwe
4271 rhosn=MIN(XSN,400.)
4273 RHOCSN=2090.* RHOSN
4274 thdifsn = 0.265/RHOCSN
4275 endif
4277 endif
4279 !--- If there is no snow melting then just evaporation
4280 !--- or condensation cxhanges SNWE
4281 ELSE
4282 EPOT=-QKMS*(QVATM-QSG)
4283 SNWE = AMAX1(0.,(SNWEPR- &
4284 BETA*EPOT*RAS*DELT))
4285 ! BETA*EPOT*RAS*UMVEG*DELT))
4287 ENDIF
4288 !*** Correct snow density on effect of snow melt, melted
4289 !*** from the top of the snow. 13% of melted water
4290 !*** remains in the pack and changes its density.
4291 !*** Eq. 9 (with my correction) in Koren et al. (1999)
4293 SNHEI=SNWE *1.E3 / RHOSN
4295 !18apr08 - if snow melt occurred then go into iteration for energy budget
4296 ! solution
4297 if(nmelt.eq.1) goto 212
4298 220 continue
4299 !-- Snow melt from the top is done. But if ground surface temperature
4300 !-- is above freezing snow can melt from the bottom. The following
4301 !-- piece of code will check if bottom melting is possible.
4303 IF(TSO(1).GT.273.15.AND.SNHEI.GT.0.) THEN
4304 if (snhei.GE.deltsn+snth) then
4305 hsn = snhei - deltsn
4306 else
4307 hsn = snhei
4308 endif
4310 soiltfrac=snowfrac*273.15+(1.-snowfrac)*TSO(1)
4312 SNOHG=(TSO(1)-soiltfrac)*(RHCS*zshalf(2)+ &
4313 RHOCSN*0.5*hsn) / DELT
4314 SNOHG=AMAX1(0.,SNOHG)
4315 SNODIF=0.
4316 SMELTG=SNOHG/XLMELT*1.E-3
4317 ! Egglston - empirical limit on snow melt from the bottom of snow pack
4318 SMELTG=AMIN1(SMELTG, 5.8e-9)
4320 if(SNWE-SMELTG*DELT.ge.rsm) then
4321 SNWE = AMAX1(0.,SNWE-SMELTG*DELT)
4322 else
4323 smeltg=snwe/delt
4324 snwe=0.
4325 rsm=0.
4326 hsn=0.
4327 endif
4329 SNOHGNEW=SMELTG*XLMELT*1.e3
4330 SNODIF=AMAX1(0.,(SNOHG-SNOHGNEW))
4331 TSO(1)=soiltfrac &
4332 + SNODIF/(RHCS*zshalf(2)+ RHOCSN*0.5*hsn)* DELT
4333 SMELT=SMELT+SMELTG
4334 SNOH=SNOH+SNOHGNEW
4336 ENDIF
4338 SNHEI=SNWE *1.E3 / RHOSN
4340 snweprint=snwe
4341 ! &
4342 !--- if VEGFRAC.ne.0. then some snow stays on the canopy
4343 !--- and should be added to SNWE for water conservation
4344 ! 4 Nov 07 +VEGFRAC*cst
4345 snheiprint=snweprint*1.E3 / RHOSN
4347 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
4348 print *, 'snweprint : ',snweprint
4349 print *, 'D9SN,SOILT,TSOB : ', D9SN,SOILT,TSOB
4350 ENDIF
4351 !--- Compute flux in the top snow layer
4352 SNFLX=D9SN*(SOILT-TSOB)
4353 IF(SNHEI.GT.0.) THEN
4354 if(ilnb.gt.1) then
4355 tsnav=0.5/snhei*((soilt+soilt1)*deltsn &
4356 +(soilt1+tso(1))*(SNHEI-DELTSN)) &
4357 -273.15
4358 else
4359 tsnav=0.5*(soilt+tso(1)) - 273.15
4360 endif
4361 ENDIF
4363 ! return
4364 ! end
4365 !------------------------------------------------------------------------
4366 END SUBROUTINE SNOWTEMP
4367 !------------------------------------------------------------------------
4370 SUBROUTINE SOILMOIST ( &
4371 !--input parameters
4372 DELT,NZS,NDDZS,DTDZS,DTDZS2,RIW, &
4373 ZSMAIN,ZSHALF,DIFFU,HYDRO, &
4374 QSG,QVG,QCG,QCATM,QVATM,PRCP, &
4375 QKMS,TRANSP,DRIP, &
4376 DEW,SMELT,SOILICE,VEGFRAC, &
4377 !--soil properties
4378 DQM,QMIN,REF,KSAT,RAS,INFMAX, &
4379 !--output
4380 SOILMOIS,SOILIQW,MAVAIL,RUNOFF,RUNOFF2,INFILTRP)
4381 !*************************************************************************
4382 ! moisture balance equation and Richards eqn.
4383 ! are solved here
4384 !
4385 ! DELT - time step (s)
4386 ! IME,JME,NZS - dimensions of soil domain
4387 ! ZSMAIN - main levels in soil (m)
4388 ! ZSHALF - middle of the soil layers (m)
4389 ! DTDZS - dt/(2.*dzshalf*dzmain)
4390 ! DTDZS2 - dt/(2.*dzshalf)
4391 ! DIFFU - diffusional conductivity (m^2/s)
4392 ! HYDRO - hydraulic conductivity (m/s)
4393 ! QSG,QVG,QCG - saturated mixing ratio, mixing ratio of
4394 ! water vapor and cloud at the ground
4395 ! surface, respectively (kg/kg)
4396 ! QCATM,QVATM - cloud and water vapor mixing ratio
4397 ! at the first atm. level (kg/kg)
4398 ! PRCP - precipitation rate in m/s
4399 ! QKMS - exchange coefficient for water vapor in the
4400 ! surface layer (m/s)
4401 ! TRANSP - transpiration from the soil layers (m/s)
4402 ! DRIP - liquid water dripping from the canopy to soil (m)
4403 ! DEW - dew in kg/m^2s
4404 ! SMELT - melting rate in m/s
4405 ! SOILICE - volumetric content of ice in soil (m^3/m^3)
4406 ! SOILIQW - volumetric content of liquid water in soil (m^3/m^3)
4407 ! VEGFRAC - greeness fraction (0-1)
4408 ! RAS - ration of air density to soil density
4409 ! INFMAX - maximum infiltration rate (kg/m^2/s)
4410 !
4411 ! SOILMOIS - volumetric soil moisture, 6 levels (m^3/m^3)
4412 ! MAVAIL - fraction of maximum soil moisture in the top
4413 ! layer (0-1)
4414 ! RUNOFF - surface runoff (m/s)
4415 ! RUNOFF2 - underground runoff (m)
4416 ! INFILTRP - point infiltration flux into soil (m/s)
4417 ! /(snow bottom runoff) (mm/s)
4418 !
4419 ! COSMC, RHSMC - coefficients for implicit solution of
4420 ! Richards equation
4421 !******************************************************************
4422 IMPLICIT NONE
4423 !------------------------------------------------------------------
4424 !--- input variables
4425 REAL, INTENT(IN ) :: DELT
4426 INTEGER, INTENT(IN ) :: NZS,NDDZS
4428 ! input variables
4430 REAL, DIMENSION(1:NZS), INTENT(IN ) :: ZSMAIN, &
4431 ZSHALF, &
4432 DIFFU, &
4433 HYDRO, &
4434 TRANSP, &
4435 SOILICE, &
4436 DTDZS2
4438 REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS
4440 REAL, INTENT(IN ) :: QSG,QVG,QCG,QCATM,QVATM , &
4441 QKMS,VEGFRAC,DRIP,PRCP , &
4442 DEW,SMELT , &
4443 DQM,QMIN,REF,KSAT,RAS,RIW
4445 ! output
4447 REAL, DIMENSION( 1:nzs ) , &
4449 INTENT(INOUT) :: SOILMOIS,SOILIQW
4451 REAL, INTENT(INOUT) :: MAVAIL,RUNOFF,RUNOFF2,INFILTRP, &
4452 INFMAX
4454 ! local variables
4456 REAL, DIMENSION( 1:nzs ) :: COSMC,RHSMC
4458 REAL :: DZS,R1,R2,R3,R4,R5,R6,R7,R8,R9,R10
4459 REAL :: REFKDT,REFDK,DELT1,F1MAX,F2MAX
4460 REAL :: F1,F2,FD,KDT,VAL,DDT,PX,FK,FKMAX
4461 REAL :: QQ,UMVEG,INFMAX1,TRANS
4462 REAL :: TOTLIQ,FLX,FLXSAT,QTOT
4463 REAL :: DID,X1,X2,X4,DENOM,Q2,Q4
4464 REAL :: dice,fcr,acrt,frzx,sum,cvfrz
4466 INTEGER :: NZS1,NZS2,K,KK,K1,KN,ialp1,jj,jk
4468 !******************************************************************************
4469 ! COEFFICIENTS FOR THOMAS ALGORITHM FOR SOILMOIS
4470 !******************************************************************************
4471 NZS1=NZS-1
4472 NZS2=NZS-2
4474 118 format(6(10Pf23.19))
4476 do k=1,nzs
4477 cosmc(k)=0.
4478 rhsmc(k)=0.
4479 enddo
4481 DID=(ZSMAIN(NZS)-ZSHALF(NZS))
4482 X1=ZSMAIN(NZS)-ZSMAIN(NZS1)
4484 !7may09 DID=(ZSMAIN(NZS)-ZSHALF(NZS))*2.
4485 ! DENOM=DID/DELT+DIFFU(NZS1)/X1
4486 ! COSMC(1)=DIFFU(NZS1)/X1/DENOM
4487 ! RHSMC(1)=(SOILMOIS(NZS)*DID/DELT
4488 ! 1 +TRANSP(NZS)-(HYDRO(NZS)*SOILMOIS(NZS)
4489 ! 1 -HYDRO(NZS1)*SOILMOIS(NZS1))*DID
4490 ! 1 /X1) /DENOM
4492 DENOM=(1.+DIFFU(nzs1)/X1/DID*DELT+HYDRO(NZS)/(2.*DID)*DELT)
4493 COSMC(1)=DELT*(DIFFU(nzs1)/DID/X1 &
4494 +HYDRO(NZS1)/2./DID)/DENOM
4495 RHSMC(1)=(SOILIQW(NZS)+TRANSP(NZS)*DELT/ &
4496 DID)/DENOM
4498 DO 330 K=1,NZS2
4499 KN=NZS-K
4500 K1=2*KN-3
4501 X4=2.*DTDZS(K1)*DIFFU(KN-1)
4502 X2=2.*DTDZS(K1+1)*DIFFU(KN)
4503 Q4=X4+HYDRO(KN-1)*DTDZS2(KN-1)
4504 Q2=X2-HYDRO(KN+1)*DTDZS2(KN-1)
4505 DENOM=1.+X2+X4-Q2*COSMC(K)
4506 COSMC(K+1)=Q4/DENOM
4507 330 RHSMC(K+1)=(SOILIQW(KN)+Q2*RHSMC(K) &
4508 +TRANSP(KN) &
4509 /(ZSHALF(KN+1)-ZSHALF(KN)) &
4510 *DELT)/DENOM
4512 ! --- MOISTURE BALANCE BEGINS HERE
4514 TRANS=TRANSP(1)
4515 UMVEG=1.-VEGFRAC
4517 RUNOFF=0.
4518 RUNOFF2=0.
4519 DZS=ZSMAIN(2)
4520 R1=COSMC(NZS1)
4521 R2= RHSMC(NZS1)
4522 R3=DIFFU(1)/DZS
4523 R4=R3+HYDRO(1)*.5
4524 R5=R3-HYDRO(2)*.5
4525 R6=QKMS*RAS
4526 !-- Total liquid water available on the top of soil domain
4527 !-- Without snow - 3 sources of water: precipitation,
4528 !-- water dripping from the canopy and dew
4529 !-- With snow - only one source of water - snow melt
4531 191 format (f23.19)
4533 TOTLIQ=UMVEG*PRCP-DRIP/DELT-UMVEG*DEW*RAS-SMELT
4536 FLX=TOTLIQ
4537 INFILTRP=TOTLIQ
4539 ! ----------- FROZEN GROUND VERSION -------------------------
4540 ! REFERENCE FROZEN GROUND PARAMETER, CVFRZ, IS A SHAPE PARAMETER OF
4541 ! AREAL DISTRIBUTION FUNCTION OF SOIL ICE CONTENT WHICH EQUALS 1/CV.
4542 ! CV IS A COEFFICIENT OF SPATIAL VARIATION OF SOIL ICE CONTENT.
4543 ! BASED ON FIELD DATA CV DEPENDS ON AREAL MEAN OF FROZEN DEPTH, AND IT
4544 ! CLOSE TO CONSTANT = 0.6 IF AREAL MEAN FROZEN DEPTH IS ABOVE 20 CM.
4545 ! THAT IS WHY PARAMETER CVFRZ = 3 (INT{1/0.6*0.6})
4546 !
4547 ! Current logic doesn't allow CVFRZ be bigger than 3
4548 CVFRZ = 3.
4550 !-- SCHAAKE/KOREN EXPRESSION for calculation of max infiltration
4551 REFKDT=3.
4552 REFDK=3.4341E-6
4553 DELT1=DELT/86400.
4554 F1MAX=DQM*ZSHALF(2)
4555 F2MAX=DQM*(ZSHALF(3)-ZSHALF(2))
4556 F1=F1MAX*(1.-SOILMOIS(1)/DQM)
4557 DICE=SOILICE(1)*ZSHALF(2)
4558 FD=F1
4559 do k=2,nzs1
4560 DICE=DICE+(ZSHALF(k+1)-ZSHALF(k))*SOILICE(K)
4561 FKMAX=DQM*(ZSHALF(k+1)-ZSHALF(k))
4562 FK=FKMAX*(1.-SOILMOIS(k)/DQM)
4563 FD=FD+FK
4564 enddo
4565 KDT=REFKDT*KSAT/REFDK
4566 VAL=(1.-EXP(-KDT*DELT1))
4567 DDT = FD*VAL
4568 PX= - TOTLIQ * DELT
4569 IF(PX.LT.0.0) PX = 0.0
4570 INFMAX1 = (PX*(DDT/(PX+DDT)))/DELT
4571 ! print *,'INFMAX1=,ksat',infmax1,ksat,f1,f2,kdt,val,ddt,px
4572 ! ----------- FROZEN GROUND VERSION --------------------------
4573 ! REDUCTION OF INFILTRATION BASED ON FROZEN GROUND PARAMETERS
4574 !
4575 ! ------------------------------------------------------------------
4577 FRZX= 0.15*((dqm+qmin)/ref) * (0.412 / 0.468)
4578 FCR = 1.
4579 IF ( DICE .GT. 1.E-2) THEN
4580 ACRT = CVFRZ * FRZX / DICE
4581 SUM = 1.
4582 IALP1 = CVFRZ - 1
4583 DO JK = 1,IALP1
4584 K = 1
4585 DO JJ = JK+1, IALP1
4586 K = K * JJ
4587 END DO
4588 SUM = SUM + (ACRT ** ( CVFRZ-JK)) / FLOAT (K)
4589 END DO
4590 FCR = 1. - EXP(-ACRT) * SUM
4591 END IF
4592 ! print *,'FCR--------',fcr
4593 INFMAX1 = INFMAX1* FCR
4594 ! -------------------------------------------------------------------
4596 INFMAX = MAX(INFMAX1,HYDRO(1)*SOILIQW(1))
4597 INFMAX = MIN(INFMAX, -TOTLIQ)
4599 !----
4600 IF (-TOTLIQ.GT.INFMAX)THEN
4601 RUNOFF=-TOTLIQ-INFMAX
4602 FLX=-INFMAX
4603 ENDIF
4604 ! INFILTRP is total infiltration flux in M/S
4605 INFILTRP=FLX
4606 ! Solution of moisture budget
4607 R7=.5*DZS/DELT
4608 R4=R4+R7
4609 FLX=FLX-SOILIQW(1)*R7
4610 R8=UMVEG*R6
4611 QTOT=QVATM+QCATM
4612 R9=TRANS
4613 R10=QTOT-QSG
4614 !-- evaporation regime
4615 IF(R10.LE.0.) THEN
4616 QQ=(R5*R2-FLX+R9)/(R4-R5*R1-R10*R8/(REF-QMIN))
4617 FLXSAT=-DQM*(R4-R5*R1-R10*R8/(REF-QMIN)) &
4618 +R5*R2+R9
4619 ELSE
4620 !-- dew formation regime
4621 QQ=(R2*R5-FLX+R8*(QTOT-QCG-QVG)+R9)/(R4-R1*R5)
4622 FLXSAT=-DQM*(R4-R1*R5)+R2*R5+R8*(QTOT-QVG-QCG)+R9
4623 END IF
4625 IF(QQ.LT.0.) THEN
4626 SOILIQW (1)=1.e-8
4627 SOILMOIS(1)=1.e-8+soilice(1)*riw
4629 ELSE IF(QQ.GT.DQM) THEN
4630 !-- saturation
4631 SOILIQW (1)=DQM
4632 SOILMOIS(1)=DQM
4633 RUNOFF2=(FLXSAT-FLX)*DELT
4634 RUNOFF=RUNOFF+(FLXSAT-FLX)
4635 ELSE
4636 SOILIQW (1)=min(dqm,max(1.e-8,QQ))
4637 SOILMOIS(1)=max(1.e-8,QQ)+soilice(1)*riw
4638 END IF
4640 !--- FINAL SOLUTION FOR SOILMOIS
4641 DO K=2,NZS
4642 KK=NZS-K+1
4643 QQ=COSMC(KK)*SOILIQW(K-1)+RHSMC(KK)
4645 IF (QQ.LT.0.) THEN
4646 SOILIQW(K) =1.e-8
4647 SOILMOIS(K)=1.e-8 + soilice(k)*riw
4649 ELSE IF(QQ.GT.DQM) THEN
4650 !-- saturation
4651 SOILIQW (K)=DQM
4652 SOILMOIS(K)=DQM
4653 IF(K.EQ.NZS)THEN
4654 RUNOFF2=RUNOFF2+(QQ-DQM)*(ZSMAIN(K)-ZSHALF(K))
4655 ELSE
4656 RUNOFF2=RUNOFF2+(QQ-DQM)*(ZSHALF(K+1)-ZSHALF(K))
4657 ENDIF
4658 ELSE
4659 SOILIQW (K)=min(dqm,max(1.e-8,QQ))
4660 SOILMOIS(K)=max(1.e-8,QQ)+soilice(k)*riw
4661 END IF
4662 END DO
4663 ! MAVAIL=min(1.,SOILMOIS(1)/(REF-QMIN))
4664 MAVAIL=min(1.,SOILMOIS(1)/DQM)
4665 if (MAVAIL.EQ.0.) MAVAIL=.00001
4667 ! RETURN
4668 ! END
4669 !-------------------------------------------------------------------
4670 END SUBROUTINE SOILMOIST
4671 !-------------------------------------------------------------------
4674 SUBROUTINE SOILPROP( &
4675 !--- input variables
4676 nzs,fwsat,lwsat,tav,keepfr, &
4677 soilmois,soiliqw,soilice, &
4678 soilmoism,soiliqwm,soilicem, &
4679 !--- soil fixed fields
4680 QWRTZ,rhocs,dqm,qmin,psis,bclh,ksat, &
4681 !--- constants
4682 riw,xlmelt,CP,G0_P,cvw,ci, &
4683 kqwrtz,kice,kwt, &
4684 !--- output variables
4685 thdif,diffu,hydro,cap)
4687 !******************************************************************
4688 ! SOILPROP computes thermal diffusivity, and diffusional and
4689 ! hydraulic condeuctivities
4690 !******************************************************************
4691 ! NX,NY,NZS - dimensions of soil domain
4692 ! FWSAT, LWSAT - volumetric content of frozen and liquid water
4693 ! for saturated condition at given temperatures (m^3/m^3)
4694 ! TAV - temperature averaged for soil layers (K)
4695 ! SOILMOIS - volumetric soil moisture at the main soil levels (m^3/m^3)
4696 ! SOILMOISM - volumetric soil moisture averaged for layers (m^3/m^3)
4697 ! SOILIQWM - volumetric liquid soil moisture averaged for layers (m^3/m^3)
4698 ! SOILICEM - volumetric content of soil ice averaged for layers (m^3/m^3)
4699 ! THDIF - thermal diffusivity for soil layers (W/m/K)
4700 ! DIFFU - diffusional conductivity (m^2/s)
4701 ! HYDRO - hydraulic conductivity (m/s)
4702 ! CAP - volumetric heat capacity (J/m^3/K)
4703 !
4704 !******************************************************************
4706 IMPLICIT NONE
4707 !-----------------------------------------------------------------
4709 !--- soil properties
4710 INTEGER, INTENT(IN ) :: NZS
4711 REAL , &
4712 INTENT(IN ) :: RHOCS, &
4713 BCLH, &
4714 DQM, &
4715 KSAT, &
4716 PSIS, &
4717 QWRTZ, &
4718 QMIN
4720 REAL, DIMENSION( 1:nzs ) , &
4721 INTENT(IN ) :: SOILMOIS, &
4722 keepfr
4725 REAL, INTENT(IN ) :: CP, &
4726 CVW, &
4727 RIW, &
4728 kqwrtz, &
4729 kice, &
4730 kwt, &
4731 XLMELT, &
4732 G0_P
4736 !--- output variables
4737 REAL, DIMENSION(1:NZS) , &
4738 INTENT(INOUT) :: cap,diffu,hydro , &
4739 thdif,tav , &
4740 soilmoism , &
4741 soiliqw,soilice , &
4742 soilicem,soiliqwm , &
4743 fwsat,lwsat
4745 !--- local variables
4746 REAL, DIMENSION(1:NZS) :: hk,detal,kasat,kjpl
4748 REAL :: x,x1,x2,x4,ws,wd,fact,fach,facd,psif,ci
4749 REAL :: tln,tavln,tn,pf,a,am,ame,h
4750 INTEGER :: nzs1,k
4752 !-- for Johansen thermal conductivity
4753 REAL :: kzero,gamd,kdry,kas,x5,sr,ke
4756 nzs1=nzs-1
4758 !-- Constants for Johansen (1975) thermal conductivity
4759 kzero =2. ! if qwrtz > 0.2
4762 do k=1,nzs
4763 detal (k)=0.
4764 kasat (k)=0.
4765 kjpl (k)=0.
4766 hk (k)=0.
4767 enddo
4769 ws=dqm+qmin
4770 x1=xlmelt/(g0_p*psis)
4771 x2=x1/bclh*ws
4772 x4=(bclh+1.)/bclh
4773 !--- Next 3 lines are for Johansen thermal conduct.
4774 gamd=(1.-ws)*2700.
4775 kdry=(0.135*gamd+64.7)/(2700.-0.947*gamd)
4776 kas=kqwrtz**qwrtz*kzero**(1.-qwrtz)
4778 DO K=1,NZS1
4779 tn=tav(k) - 273.15
4780 wd=ws - riw*soilicem(k)
4781 psif=psis*100.*(wd/(soiliqwm(k)+qmin))**bclh &
4782 * (ws/wd)**3.
4783 !--- PSIF should be in [CM] to compute PF
4784 pf=log10(abs(psif))
4785 fact=1.+riw*soilicem(k)
4786 !--- HK is for McCumber thermal conductivity
4787 IF(PF.LE.5.2) THEN
4788 HK(K)=420.*EXP(-(PF+2.7))*fact
4789 ELSE
4790 HK(K)=.1744*fact
4791 END IF
4793 IF(soilicem(k).NE.0.AND.TN.LT.0.) then
4794 !--- DETAL is taking care of energy spent on freezing or released from
4795 ! melting of soil water
4797 DETAL(K)=273.15*X2/(TAV(K)*TAV(K))* &
4798 (TAV(K)/(X1*TN))**X4
4800 if(keepfr(k).eq.1.) then
4801 detal(k)=0.
4802 endif
4804 ENDIF
4806 !--- Next 10 lines calculate Johansen thermal conductivity KJPL
4807 kasat(k)=kas**(1.-ws)*kice**fwsat(k) &
4808 *kwt**lwsat(k)
4810 X5=(soilmoism(k)+qmin)/ws
4811 if(soilicem(k).eq.0.) then
4812 sr=max(0.101,x5)
4813 ke=log10(sr)+1.
4814 !--- next 2 lines - for coarse soils
4815 ! sr=max(0.0501,x5)
4816 ! ke=0.7*log10(sr)+1.
4817 else
4818 ke=x5
4819 endif
4821 kjpl(k)=ke*(kasat(k)-kdry)+kdry
4823 !--- CAP -volumetric heat capacity
4824 CAP(K)=(1.-WS)*RHOCS &
4825 + (soiliqwm(K)+qmin)*CVW &
4826 + soilicem(K)*CI &
4827 + (dqm-soilmoism(k))*CP*1.2 &
4828 - DETAL(K)*1.e3*xlmelt
4830 a=RIW*soilicem(K)
4832 if((ws-a).lt.0.12)then
4833 diffu(K)=0.
4834 else
4835 H=max(0.,(soilmoism(K)-a)/(max(1.e-8,(dqm-a))))
4836 facd=1.
4837 if(a.ne.0.)facd=1.-a/max(1.e-8,soilmoism(K))
4838 ame=max(1.e-8,dqm-riw*soilicem(K))
4839 !--- DIFFU is diffusional conductivity of soil water
4840 diffu(K)=-BCLH*KSAT*PSIS/ame* &
4841 (dqm/ame)**3. &
4842 *H**(BCLH+2.)*facd
4843 endif
4845 ! diffu(K)=-BCLH*KSAT*PSIS/dqm &
4846 ! *H**(BCLH+2.)
4849 !--- thdif - thermal diffusivity
4850 ! thdif(K)=HK(K)/CAP(K)
4851 !--- Use thermal conductivity from Johansen (1975)
4852 thdif(K)=KJPL(K)/CAP(K)
4854 END DO
4856 DO K=1,NZS
4858 if((ws-riw*soilice(k)).lt.0.12)then
4859 hydro(k)=0.
4860 else
4861 fach=1.
4862 if(soilice(k).ne.0.) &
4863 fach=1.-riw*soilice(k)/max(1.e-8,soilmois(k))
4864 am=max(1.e-8,dqm-riw*soilice(k))
4865 !--- HYDRO is hydraulic conductivity of soil water
4866 hydro(K)=KSAT/am* &
4867 (soiliqw(K)/am) &
4868 **(2.*BCLH+2.) &
4869 * fach
4870 endif
4872 ENDDO
4874 ! RETURN
4875 ! END
4877 !-----------------------------------------------------------------------
4878 END SUBROUTINE SOILPROP
4879 !-----------------------------------------------------------------------
4882 SUBROUTINE TRANSF( &
4883 !--- input variables
4884 nzs,nroot,soiliqw,tabs, &
4885 !--- soil fixed fields
4886 dqm,qmin,ref,wilt,zshalf, &
4887 !--- output variables
4888 tranf,transum)
4890 !-------------------------------------------------------------------
4891 !--- TRANF(K) - THE TRANSPIRATION FUNCTION (Smirnova et al. 1996, EQ. 18,19)
4892 !*******************************************************************
4893 ! NX,NY,NZS - dimensions of soil domain
4894 ! SOILIQW - volumetric liquid soil moisture at the main levels (m^3/m^3)
4895 ! TRANF - the transpiration function at levels (m)
4896 ! TRANSUM - transpiration function integrated over the rooting zone (m)
4897 !
4898 !*******************************************************************
4899 IMPLICIT NONE
4900 !-------------------------------------------------------------------
4902 !--- input variables
4904 INTEGER, INTENT(IN ) :: nroot,nzs
4906 REAL , &
4907 INTENT(IN ) :: TABS
4908 !--- soil properties
4909 REAL , &
4910 INTENT(IN ) :: DQM, &
4911 QMIN, &
4912 REF, &
4913 WILT
4915 REAL, DIMENSION(1:NZS), INTENT(IN) :: soiliqw, &
4916 ZSHALF
4918 !-- output
4919 REAL, DIMENSION(1:NZS), INTENT(OUT) :: TRANF
4920 REAL, INTENT(OUT) :: TRANSUM
4922 !-- local variables
4923 REAL :: totliq, did
4924 INTEGER :: k
4926 !-- for non-linear root distribution
4927 REAL :: gx,sm1,sm2,sm3,sm4,ap0,ap1,ap2,ap3,ap4
4928 REAL :: FTEM
4929 REAL, DIMENSION(1:NZS) :: PART
4930 !--------------------------------------------------------------------
4932 do k=1,nzs
4933 part(k)=0.
4934 enddo
4936 transum=0.
4937 totliq=soiliqw(1)+qmin
4938 sm1=totliq
4939 sm2=sm1*sm1
4940 sm3=sm2*sm1
4941 sm4=sm3*sm1
4942 ap0=0.299
4943 ap1=-8.152
4944 ap2=61.653
4945 ap3=-115.876
4946 ap4=59.656
4947 gx=ap0+ap1*sm1+ap2*sm2+ap3*sm3+ap4*sm4
4948 if(totliq.ge.ref) gx=1.
4949 if(totliq.le.0.) gx=0.
4950 if(gx.gt.1.) gx=1.
4951 if(gx.lt.0.) gx=0.
4952 DID=zshalf(2)
4953 part(1)=DID*gx
4954 !--- air temperature function
4955 ! Avissar (1985) and AX 7/95
4956 IF (TABS .LE. 302.15) THEN
4957 FTEM = 1.0 / (1.0 + EXP(-0.41 * (TABS - 282.05)))
4958 ELSE
4959 FTEM = 1.0 / (1.0 + EXP(0.5 * (TABS - 314.0)))
4960 ENDIF
4961 !---
4962 IF(TOTLIQ.GT.REF) THEN
4963 TRANF(1)=DID
4964 ELSE IF(TOTLIQ.LE.WILT) THEN
4965 TRANF(1)=0.
4966 ELSE
4967 TRANF(1)=(TOTLIQ-WILT)/(REF-WILT)*DID
4968 ENDIF
4969 !-- uncomment next line for non-linear root distribution
4970 !cc TRANF(1)=part(1)
4971 TRANF(1)=TRANF(1)*FTEM
4973 DO K=2,NROOT
4974 totliq=soiliqw(k)+qmin
4975 sm1=totliq
4976 sm2=sm1*sm1
4977 sm3=sm2*sm1
4978 sm4=sm3*sm1
4979 gx=ap0+ap1*sm1+ap2*sm2+ap3*sm3+ap4*sm4
4980 if(totliq.ge.ref) gx=1.
4981 if(totliq.le.0.) gx=0.
4982 if(gx.gt.1.) gx=1.
4983 if(gx.lt.0.) gx=0.
4984 DID=zshalf(K+1)-zshalf(K)
4985 part(k)=did*gx
4986 IF(totliq.GE.REF) THEN
4987 TRANF(K)=DID
4988 ELSE IF(totliq.LE.WILT) THEN
4989 TRANF(K)=0.
4990 ELSE
4991 TRANF(K)=(totliq-WILT) &
4992 /(REF-WILT)*DID
4993 ENDIF
4994 !-- uncomment next line for non-linear root distribution
4995 !cc TRANF(k)=part(k)
4996 TRANF(k)=TRANF(k)*FTEM
4997 END DO
4999 !-- TRANSUM - total for the rooting zone
5000 transum=0.
5001 DO K=1,NROOT
5002 transum=transum+tranf(k)
5003 END DO
5005 ! RETURN
5006 ! END
5007 !-----------------------------------------------------------------
5008 END SUBROUTINE TRANSF
5009 !-----------------------------------------------------------------
5012 SUBROUTINE VILKA(TN,D1,D2,PP,QS,TS,TT,NSTEP,ii,j,iland,isoil)
5013 !--------------------------------------------------------------
5014 !--- VILKA finds the solution of energy budget at the surface
5015 !--- using table T,QS computed from Clausius-Klapeiron
5016 !--------------------------------------------------------------
5017 REAL, DIMENSION(1:4001), INTENT(IN ) :: TT
5018 REAL, INTENT(IN ) :: TN,D1,D2,PP
5019 INTEGER, INTENT(IN ) :: NSTEP,ii,j,iland,isoil
5021 REAL, INTENT(OUT ) :: QS, TS
5023 REAL :: F1,T1,T2,RN
5024 INTEGER :: I,I1
5026 I=(TN-1.7315E2)/.05+1
5027 T1=173.1+FLOAT(I)*.05
5028 F1=T1+D1*TT(I)-D2
5029 I1=I-F1/(.05+D1*(TT(I+1)-TT(I)))
5030 I=I1
5031 IF(I.GT.4000.OR.I.LT.1) GOTO 1
5032 10 I1=I
5033 T1=173.1+FLOAT(I)*.05
5034 F1=T1+D1*TT(I)-D2
5035 RN=F1/(.05+D1*(TT(I+1)-TT(I)))
5036 I=I-INT(RN)
5037 IF(I.GT.4000.OR.I.LT.1) GOTO 1
5038 IF(I1.NE.I) GOTO 10
5039 TS=T1-.05*RN
5040 QS=(TT(I)+(TT(I)-TT(I+1))*RN)/PP
5041 GOTO 20
5042 1 PRINT *,' AVOST IN VILKA '
5043 ! WRITE(12,*)'AVOST',TN,D1,D2,PP,NSTEP
5044 PRINT *,TN,D1,D2,PP,NSTEP,I,TT(i),ii,j,iland,isoil
5045 CALL wrf_error_fatal (' AVOST IN VILKA ' )
5046 20 CONTINUE
5047 ! RETURN
5048 ! END
5049 !-----------------------------------------------------------------------
5050 END SUBROUTINE VILKA
5051 !-----------------------------------------------------------------------
5054 SUBROUTINE SOILVEGIN ( IVGTYP,ISLTYP,MYJ, &
5055 IFOREST,EMISS,PC,ZNT,QWRTZ, &
5056 RHOCS,BCLH,DQM,KSAT,PSIS,QMIN,REF,WILT )
5058 !************************************************************************
5059 ! Set-up soil and vegetation Parameters in the case when
5060 ! snow disappears during the forecast and snow parameters
5061 ! shold be replaced by surface parameters according to
5062 ! soil and vegetation types in this point.
5063 !
5064 ! Output:
5065 !
5066 !
5067 ! Soil parameters:
5068 ! DQM: MAX soil moisture content - MIN (m^3/m^3)
5069 ! REF: Reference soil moisture (m^3/m^3)
5070 ! WILT: Wilting PT soil moisture contents (m^3/m^3)
5071 ! QMIN: Air dry soil moist content limits (m^3/m^3)
5072 ! PSIS: SAT soil potential coefs. (m)
5073 ! KSAT: SAT soil diffusivity/conductivity coefs. (m/s)
5074 ! BCLH: Soil diffusivity/conductivity exponent.
5075 !
5076 ! ************************************************************************
5077 IMPLICIT NONE
5078 !---------------------------------------------------------------------------
5079 integer, parameter :: nsoilclas=19
5080 integer, parameter :: nvegclas=24+3
5081 integer, parameter :: iwater=16
5082 integer, parameter :: ilsnow=99
5085 !--- soiltyp classification according to STATSGO(nclasses=16)
5086 !
5087 ! 1 SAND SAND
5088 ! 2 LOAMY SAND LOAMY SAND
5089 ! 3 SANDY LOAM SANDY LOAM
5090 ! 4 SILT LOAM SILTY LOAM
5091 ! 5 SILT SILTY LOAM
5092 ! 6 LOAM LOAM
5093 ! 7 SANDY CLAY LOAM SANDY CLAY LOAM
5094 ! 8 SILTY CLAY LOAM SILTY CLAY LOAM
5095 ! 9 CLAY LOAM CLAY LOAM
5096 ! 10 SANDY CLAY SANDY CLAY
5097 ! 11 SILTY CLAY SILTY CLAY
5098 ! 12 CLAY LIGHT CLAY
5099 ! 13 ORGANIC MATERIALS LOAM
5100 ! 14 WATER
5101 ! 15 BEDROCK
5102 ! Bedrock is reclassified as class 14
5103 ! 16 OTHER (land-ice)
5104 ! 17 Playa
5105 ! 18 Lava
5106 ! 19 White Sand
5107 !
5108 !----------------------------------------------------------------------
5109 REAL LQMA(nsoilclas),LRHC(nsoilclas), &
5110 LPSI(nsoilclas),LQMI(nsoilclas), &
5111 LBCL(nsoilclas),LKAS(nsoilclas), &
5112 LWIL(nsoilclas),LREF(nsoilclas), &
5113 DATQTZ(nsoilclas)
5114 !-- LQMA Rawls et al.[1982]
5115 ! DATA LQMA /0.417, 0.437, 0.453, 0.501, 0.486, 0.463, 0.398,
5116 ! & 0.471, 0.464, 0.430, 0.479, 0.475, 0.439, 1.0, 0.20, 0.401/
5117 !---
5118 !-- Clapp, R. and G. Hornberger, 1978: Empirical equations for some soil
5119 ! hydraulic properties, Water Resour. Res., 14, 601-604.
5121 !-- Clapp et al. [1978]
5122 DATA LQMA /0.395, 0.410, 0.435, 0.485, 0.485, 0.451, 0.420, &
5123 0.477, 0.476, 0.426, 0.492, 0.482, 0.451, 1.0, &
5124 0.20, 0.435, 0.468, 0.200, 0.339/
5126 !-- LREF Rawls et al.[1982]
5127 ! DATA LREF /0.091, 0.125, 0.207, 0.330, 0.360, 0.270, 0.255,
5128 ! & 0.366, 0.318, 0.339, 0.387, 0.396, 0.329, 1.0, 0.108, 0.283/
5130 !-- Clapp et al. [1978]
5131 DATA LREF /0.174, 0.179, 0.249, 0.369, 0.369, 0.314, 0.299, &
5132 0.357, 0.391, 0.316, 0.409, 0.400, 0.314, 1., &
5133 0.1, 0.249, 0.454, 0.17, 0.236/
5135 !-- LWIL Rawls et al.[1982]
5136 ! DATA LWIL/0.033, 0.055, 0.095, 0.133, 0.133, 0.117, 0.148,
5137 ! & 0.208, 0.197, 0.239, 0.250, 0.272, 0.066, 0.0, 0.006, 0.029/
5139 !-- Clapp et al. [1978]
5140 DATA LWIL/0.068, 0.075, 0.114, 0.179, 0.179, 0.155, 0.175, &
5141 0.218, 0.250, 0.219, 0.283, 0.286, 0.155, 0.0, &
5142 0.006, 0.114, 0.030, 0.006, 0.01/
5144 ! DATA LQMI/0.010, 0.028, 0.047, 0.084, 0.084, 0.066, 0.067,
5145 ! & 0.120, 0.103, 0.100, 0.126, 0.138, 0.066, 0.0, 0.006, 0.028/
5147 !-- Carsel and Parrish [1988]
5148 DATA LQMI/0.045, 0.057, 0.065, 0.067, 0.034, 0.078, 0.10, &
5149 0.089, 0.095, 0.10, 0.070, 0.068, 0.078, 0.0, &
5150 0.004, 0.065, 0.020, 0.004, 0.008/
5152 !-- LPSI Cosby et al[1984]
5153 ! DATA LPSI/0.060, 0.036, 0.141, 0.759, 0.759, 0.355, 0.135,
5154 ! & 0.617, 0.263, 0.098, 0.324, 0.468, 0.355, 0.0, 0.069, 0.036/
5155 ! & 0.617, 0.263, 0.098, 0.324, 0.468, 0.355, 0.0, 0.069, 0.036/
5157 !-- Clapp et al. [1978]
5158 DATA LPSI/0.121, 0.090, 0.218, 0.786, 0.786, 0.478, 0.299, &
5159 0.356, 0.630, 0.153, 0.490, 0.405, 0.478, 0.0, &
5160 0.121, 0.218, 0.468, 0.069, 0.069/
5162 !-- LKAS Rawls et al.[1982]
5163 ! DATA LKAS/5.83E-5, 1.70E-5, 7.19E-6, 1.89E-6, 1.89E-6,
5164 ! & 3.67E-6, 1.19E-6, 4.17E-7, 6.39E-7, 3.33E-7, 2.50E-7,
5165 ! & 1.67E-7, 3.38E-6, 0.0, 1.41E-4, 1.41E-5/
5167 !-- Clapp et al. [1978]
5168 DATA LKAS/1.76E-4, 1.56E-4, 3.47E-5, 7.20E-6, 7.20E-6, &
5169 6.95E-6, 6.30E-6, 1.70E-6, 2.45E-6, 2.17E-6, &
5170 1.03E-6, 1.28E-6, 6.95E-6, 0.0, 1.41E-4, &
5171 3.47E-5, 1.28E-6, 1.41E-4, 1.76E-4/
5173 !-- LBCL Cosby et al [1984]
5174 ! DATA LBCL/2.79, 4.26, 4.74, 5.33, 5.33, 5.25, 6.66,
5175 ! & 8.72, 8.17, 10.73, 10.39, 11.55, 5.25, 0.0, 2.79, 4.26/
5177 !-- Clapp et al. [1978]
5178 DATA LBCL/4.05, 4.38, 4.90, 5.30, 5.30, 5.39, 7.12, &
5179 7.75, 8.52, 10.40, 10.40, 11.40, 5.39, 0.0, &
5180 4.05, 4.90, 11.55, 2.79, 2.79/
5182 DATA LRHC /1.47,1.41,1.34,1.27,1.27,1.21,1.18,1.32,1.23, &
5183 1.18,1.15,1.09,1.21,4.18,2.03,2.10,1.09,2.03,1.47/
5185 DATA DATQTZ/0.92,0.82,0.60,0.25,0.10,0.40,0.60,0.10,0.35, &
5186 0.52,0.10,0.25,0.00,0.,0.60,0.0,0.25,0.60,0.92/
5188 !--------------------------------------------------------------------------
5189 !
5190 ! USGS Vegetation Types
5191 !
5192 ! 1: Urban and Built-Up Land
5193 ! 2: Dryland Cropland and Pasture
5194 ! 3: Irrigated Cropland and Pasture
5195 ! 4: Mixed Dryland/Irrigated Cropland and Pasture
5196 ! 5: Cropland/Grassland Mosaic
5197 ! 6: Cropland/Woodland Mosaic
5198 ! 7: Grassland
5199 ! 8: Shrubland
5200 ! 9: Mixed Shrubland/Grassland
5201 ! 10: Savanna
5202 ! 11: Deciduous Broadleaf Forest
5203 ! 12: Deciduous Needleleaf Forest
5204 ! 13: Evergreen Broadleaf Forest
5205 ! 14: Evergreen Needleleaf Fores
5206 ! 15: Mixed Forest
5207 ! 16: Water Bodies
5208 ! 17: Herbaceous Wetland
5209 ! 18: Wooded Wetland
5210 ! 19: Barren or Sparsely Vegetated
5211 ! 20: Herbaceous Tundra
5212 ! 21: Wooded Tundra
5213 ! 22: Mixed Tundra
5214 ! 23: Bare Ground Tundra
5215 ! 24: Snow or Ice
5216 !
5217 ! 25: Playa
5218 ! 26: Lava
5219 ! 27: White Sand
5222 !---- Below are the arrays for the vegetation parameters
5223 REAL LALB(nvegclas),LMOI(nvegclas),LEMI(nvegclas), &
5224 LROU(nvegclas),LTHI(nvegclas),LSIG(nvegclas), &
5225 LPC(nvegclas), NROTBL(nvegclas)
5227 !************************************************************************
5228 !---- vegetation parameters
5229 !
5230 !-- USGS model
5231 !
5232 DATA LALB/.18,.17,.18,.18,.18,.16,.19,.22,.20,.20,.16,.14, &
5233 .12,.12,.13,.08,.14,.14,.25,.15,.15,.15,.25,.55, &
5234 .30,.16,.60 /
5235 DATA LEMI/.88,4*.92,.93,.92,.88,.9,.92,.93,.94, &
5236 .95,.95,.94,.98,.95,.95,.85,.92,.93,.92,.85,.95, &
5237 .85,.85,.90 /
5238 !-- Roughness length is changed for forests and some others
5239 ! DATA LROU/.5,.06,.075,.065,.05,.2,.075,.1,.11,.15,.8,.85, &
5240 ! 2.0,1.0,.563,.0001,.2,.4,.05,.1,.15,.1,.065,.05/
5241 DATA LROU/.5,.06,.075,.065,.05,.2,.075,.1,.11,.15,.5,.5, &
5242 .5,.5,.5,.0001,.2,.4,.05,.1,.15,.1,.065,.05, &
5243 .01,.15,.01 /
5245 DATA LMOI/.1,.3,.5,.25,.25,.35,.15,.1,.15,.15,.3,.3, &
5246 .5,.3,.3,1.,.6,.35,.02,.5,.5,.5,.02,.95,.40,.50,.40/
5247 !
5248 !---- still needs to be corrected
5249 !
5250 ! DATA LPC/ 15*.8,0.,.8,.8,.5,.5,.5,.5,.5,.0/
5251 DATA LPC /0.4,0.3,0.4,0.4,0.4,0.4,0.4,0.4,0.4,0.4,5*0.55,0.,0.55,0.55, &
5252 0.3,0.3,0.4,0.4,0.3,0.,.3,0.,0./
5254 ! used in RUC DATA LPC /0.6,6*0.8,0.7,0.75,6*0.8,0.,0.8,0.8, &
5255 ! 0.5,0.7,0.6,0.7,0.5,0./
5258 !***************************************************************************
5261 INTEGER :: &
5262 IVGTYP, &
5263 ISLTYP
5265 LOGICAL, INTENT(IN ) :: myj
5267 REAL , &
5268 INTENT ( OUT) :: pc
5270 REAL , &
5271 INTENT (INOUT ) :: emiss, &
5272 znt
5273 !--- soil properties
5274 REAL , &
5275 INTENT( OUT) :: RHOCS, &
5276 BCLH, &
5277 DQM, &
5278 KSAT, &
5279 PSIS, &
5280 QMIN, &
5281 QWRTZ, &
5282 REF, &
5283 WILT
5285 INTEGER, DIMENSION( 1:(nvegclas) ) , &
5286 INTENT ( OUT) :: iforest
5290 INTEGER, DIMENSION( 1:(nvegclas) ) :: if1
5291 INTEGER :: kstart, kfin, lstart, lfin
5292 INTEGER :: i,j,k
5294 !***********************************************************************
5295 ! DATA ZS1/0.0,0.05,0.20,0.40,1.6,3.0/ ! o - levels in soil
5296 ! DATA ZS2/0.0,0.025,0.125,0.30,1.,2.3/ ! x - levels in soil
5297 DATA IF1/12*0,1,1,1,12*0/
5299 do k=1,nvegclas
5300 iforest(k)=if1(k)
5301 enddo
5304 ! EMISS = LEMI(IVGTYP)
5305 EMISS = LEMITBL(IVGTYP)
5306 ! When MYJ sfc scheme is used - better use recommended in MYJSFCINIT
5307 ! values of roughness length, and not redefine it here.
5308 ! The table in this routine is the one we use in RUC with RUC LSM.
5310 !tgs 3 Oct 2007 - MYJSFCINIT roughness produced unrealistic fluxes,
5311 !tgs - will not use it any more!
5312 !tgs if (.not. myj) then
5313 ! ZNT = LROU(IVGTYP)
5314 ZNT = Z0TBL(IVGTYP)
5315 !tgs endif
5317 ! PC = LPC (IVGTYP)
5318 PC = PCTBL(IVGTYP)
5319 ! RHOCS = LRHC(ISLTYP)*1.E6
5320 RHOCS = HC(ISLTYP)*1.E6
5322 ! parameters from SOILPARM.TBL
5323 BCLH = BB(ISLTYP)
5324 DQM = MAXSMC(ISLTYP)- &
5325 DRYSMC(ISLTYP)
5326 KSAT = SATDK(ISLTYP)
5327 PSIS = - SATPSI(ISLTYP)
5328 QMIN = DRYSMC(ISLTYP)
5329 REF = REFSMC(ISLTYP)
5330 WILT = WLTSMC(ISLTYP)
5331 QWRTZ = QTZ(ISLTYP)
5333 ! parameters from the look-up tables
5334 ! BCLH = LBCL(ISLTYP)
5335 ! DQM = LQMA(ISLTYP)- &
5336 ! LQMI(ISLTYP)
5337 ! KSAT = LKAS(ISLTYP)
5338 ! PSIS = - LPSI(ISLTYP)
5339 ! QMIN = LQMI(ISLTYP)
5340 ! REF = LREF(ISLTYP)
5341 ! WILT = LWIL(ISLTYP)
5342 ! QWRTZ = DATQTZ(ISLTYP)
5344 !--------------------------------------------------------------------------
5345 END SUBROUTINE SOILVEGIN
5346 !--------------------------------------------------------------------------
5349 SUBROUTINE SNOWFREE (ivgtyp,myj,emiss,znt,iland)
5350 !************************************************************************
5351 ! Set-up soil and vegetation Parameters in the case when
5352 ! snow disappears during the forecast and snow parameters
5353 ! shold be replaced by surface parameters according to
5354 ! soil and vegetation types in this point.
5355 !
5356 !***************************************************************************
5357 IMPLICIT NONE
5358 !---------------------------------------------------------------------------
5359 integer, parameter :: nvegclas=24+3
5362 INTEGER :: IVGTYP
5364 LOGICAL, INTENT(IN ) :: myj
5366 REAL, INTENT(INOUT) :: &
5367 emiss, &
5368 znt
5369 INTEGER, INTENT(INOUT) :: ILAND
5371 !---- Below are the arrays for the vegetation parameters
5372 REAL, DIMENSION( 1:nvegclas ) :: LALB, &
5373 LEMI, &
5374 LROU_MYJ,&
5375 LROU
5377 !************************************************************************
5378 !-- USGS model
5379 !
5380 DATA LALB/.18,.17,.18,.18,.18,.16,.19,.22,.20,.20,.16,.14, &
5381 .12,.12,.13,.08,.14,.14,.25,.15,.15,.15,.25,.55, &
5382 .30,.16,.60/
5383 DATA LEMI/.88,4*.92,.93,.92,.88,.9,.92,.93,.94, &
5384 .95,.95,.94,.98,.95,.95,.85,.92,.93,.92,.85,.95, &
5385 .85,.85,.90 /
5386 !-- Roughness length is changed for forests and some others
5387 ! next 2 lines - table from RUC
5388 ! DATA LROU/.5,.06,.075,.065,.05,.2,.075,.1,.11,.15,.8,.85, &
5389 ! 2.0,1.0,.563,.0001,.2,.4,.05,.1,.15,.1,.065,.05/
5391 DATA LROU/.5,.06,.075,.065,.05,.2,.075,.1,.11,.15,.5,.5, &
5392 .5,.5,.5,.0001,.2,.4,.05,.1,.15,.1,.065,.05, &
5393 .01,.15,.01 /
5395 ! With MYJSFC better use the table from MYJSFCINIT
5396 DATA LROU_MYJ/1.0,.07,.07,.07,.07,.15,.08,.03,.05,.86,.8,.85, &
5397 2.65,1.09,.8,.001,.04,.05,.01,.04,.06,.05,.03,.001, &
5398 .01,.15,.01 /
5402 !--------------------------------------------------------------------------
5404 EMISS = LEMITBL(IVGTYP)
5405 !tgs 3 Oct 2007 - LROU_MYJ gives unrealistic surface fluxes with RUC LSM with RUC LSM
5406 ! if(myj) then
5407 ! ZNT = LROU_MYJ(IVGTYP)
5408 ! else
5409 ZNT = Z0TBL(IVGTYP)
5410 !!! ZNT = LROU(IVGTYP)
5411 ! endif
5412 ILAND = IVGTYP
5413 ! ---
5415 ! RETURN
5416 ! END
5417 !--------------------------------------------------------------------------
5418 END SUBROUTINE SNOWFREE
5420 SUBROUTINE RUCLSMINIT( SH2O,SMFR3D,TSLB,SMOIS,ISLTYP,IVGTYP, &
5421 XICE,mavail,nzs, iswater, isice, restart, &
5422 allowed_to_read , &
5423 ids,ide, jds,jde, kds,kde, &
5424 ims,ime, jms,jme, kms,kme, &
5425 its,ite, jts,jte, kts,kte )
5426 #ifdef WRF_CHEM
5427 USE module_data_gocart_dust
5428 #endif
5429 IMPLICIT NONE
5432 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
5433 ims,ime, jms,jme, kms,kme, &
5434 its,ite, jts,jte, kts,kte, &
5435 nzs, iswater, isice
5437 REAL, DIMENSION( ims:ime, 1:nzs, jms:jme ) , &
5438 INTENT(IN) :: TSLB, &
5439 SMOIS
5441 INTEGER, DIMENSION( ims:ime, jms:jme ) , &
5442 INTENT(INOUT) :: ISLTYP,IVGTYP
5444 REAL, DIMENSION( ims:ime, 1:nzs, jms:jme ) , &
5445 INTENT(INOUT) :: SMFR3D, &
5446 SH2O
5448 REAL, DIMENSION( ims:ime, jms:jme ) , &
5449 INTENT(INOUT) :: XICE,MAVAIL
5451 REAL, DIMENSION ( 1:nzs ) :: SOILIQW
5453 LOGICAL , INTENT(IN) :: restart, allowed_to_read
5455 !
5456 INTEGER :: I,J,L,itf,jtf
5457 REAL :: RIW,XLMELT,TLN,DQM,REF,PSIS,QMIN,BCLH
5459 INTEGER :: errflag
5461 ! itf=min0(ite,ide-1)
5462 ! jtf=min0(jte,jde-1)
5465 RIW=900.*1.e-3
5466 XLMELT=3.35E+5
5468 ! initialize three LSM related tables
5469 IF ( allowed_to_read ) THEN
5470 CALL wrf_message( 'INITIALIZE THREE LSM RELATED TABLES' )
5471 CALL RUCLSM_PARM_INIT
5472 ENDIF
5474 #ifdef WRF_CHEM
5475 !
5476 ! need this parameter for dust parameterization in wrf/chem
5477 !
5478 do I=1,NSLTYPE
5479 porosity(i)=maxsmc(i)
5480 enddo
5481 #endif
5483 IF(.not.restart)THEN
5485 itf=min0(ite,ide-1)
5486 jtf=min0(jte,jde-1)
5488 errflag = 0
5489 DO j = jts,jtf
5490 DO i = its,itf
5491 IF ( ISLTYP( i,j ) .LT. 1 ) THEN
5492 errflag = 1
5493 WRITE(err_message,*)"module_sf_ruclsm.F: lsminit: out of range ISLTYP ",i,j,ISLTYP( i,j )
5494 CALL wrf_message(err_message)
5495 ENDIF
5496 ENDDO
5497 ENDDO
5498 IF ( errflag .EQ. 1 ) THEN
5499 CALL wrf_error_fatal( "module_sf_ruclsm.F: lsminit: out of range value "// &
5500 "of ISLTYP. Is this field in the input?" )
5501 ENDIF
5503 DO J=jts,jtf
5504 DO I=its,itf
5506 ! CALL SOILIN ( ISLTYP(I,J), DQM, REF, PSIS, QMIN, BCLH )
5509 !--- Computation of volumetric content of ice in soil
5510 !--- and initialize MAVAIL
5511 DQM = MAXSMC (ISLTYP(I,J)) - &
5512 DRYSMC (ISLTYP(I,J))
5513 REF = REFSMC (ISLTYP(I,J))
5514 PSIS = - SATPSI (ISLTYP(I,J))
5515 QMIN = DRYSMC (ISLTYP(I,J))
5516 BCLH = BB (ISLTYP(I,J))
5519 !!! IF (.not.restart) THEN
5521 IF(xice(i,j).gt.0.) THEN
5522 !-- for ice
5523 DO L=1,NZS
5524 smfr3d(i,l,j)=1.
5525 sh2o(i,l,j)=0.
5526 mavail(i,j) = 1.
5527 ENDDO
5528 ELSE
5529 if(isltyp(i,j).ne.14 ) then
5530 !-- land
5531 mavail(i,j) = max(0.00001,min(1.,smois(i,1,j)/dqm))
5532 ! mavail(i,j) = max(0.00001,min(1.,smois(i,1,j)/(ref-qmin)))
5533 DO L=1,NZS
5534 !-- for land points initialize soil ice
5535 tln=log(TSLB(i,l,j)/273.15)
5537 if(tln.lt.0.) then
5538 soiliqw(l)=(dqm+qmin)*(XLMELT* &
5539 (tslb(i,l,j)-273.15)/tslb(i,l,j)/9.81/psis) &
5540 **(-1./bclh)-qmin
5541 soiliqw(l)=max(0.,soiliqw(l))
5542 soiliqw(l)=min(soiliqw(l),smois(i,l,j))
5543 sh2o(i,l,j)=soiliqw(l)
5544 smfr3d(i,l,j)=(smois(i,l,j)-soiliqw(l))/RIW
5546 else
5547 smfr3d(i,l,j)=0.
5548 sh2o(i,l,j)=smois(i,l,j)
5549 endif
5550 ENDDO
5552 else
5553 !-- for water ISLTYP=14
5554 DO L=1,NZS
5555 smfr3d(i,l,j)=0.
5556 sh2o(i,l,j)=1.
5557 mavail(i,j) = 1.
5558 ENDDO
5559 endif
5560 ENDIF
5562 ENDDO
5563 ENDDO
5565 ENDIF
5567 END SUBROUTINE ruclsminit
5568 !
5569 !-----------------------------------------------------------------
5570 SUBROUTINE RUCLSM_PARM_INIT
5571 !-----------------------------------------------------------------
5573 character*8 :: MMINLU, MMINSL
5575 MMINLU='USGS-RUC'
5576 MMINSL='STAS-RUC'
5577 call RUCLSM_SOILVEGPARM( MMINLU, MMINSL)
5579 !-----------------------------------------------------------------
5580 END SUBROUTINE RUCLSM_PARM_INIT
5581 !-----------------------------------------------------------------
5583 !-----------------------------------------------------------------
5584 SUBROUTINE RUCLSM_SOILVEGPARM( MMINLU, MMINSL)
5585 !-----------------------------------------------------------------
5587 IMPLICIT NONE
5589 integer :: LUMATCH, IINDEX, LC, NUM_SLOPE
5590 integer :: ierr
5591 INTEGER , PARAMETER :: OPEN_OK = 0
5593 character*8 :: MMINLU, MMINSL
5594 character*128 :: mess , message, vege_parm_string
5595 logical, external :: wrf_dm_on_monitor
5598 !-----SPECIFY VEGETATION RELATED CHARACTERISTICS :
5599 ! ALBBCK: SFC albedo (in percentage)
5600 ! Z0: Roughness length (m)
5601 ! LEMI: Emissivity
5602 ! PC: Plant coefficient for transpiration function
5603 ! -- the rest of the parameters are read in but not used currently
5604 ! SHDFAC: Green vegetation fraction (in percentage)
5605 ! Note: The ALBEDO, Z0, and SHDFAC values read from the following table
5606 ! ALBEDO, amd Z0 are specified in LAND-USE TABLE; and SHDFAC is
5607 ! the monthly green vegetation data
5608 ! CMXTBL: MAX CNPY Capacity (m)
5609 ! NROTBL: Rooting depth (layer)
5610 ! RSMIN: Mimimum stomatal resistance (s m-1)
5611 ! RSMAX: Max. stomatal resistance (s m-1)
5612 ! RGL: Parameters used in radiation stress function
5613 ! HS: Parameter used in vapor pressure deficit functio
5614 ! TOPT: Optimum transpiration air temperature. (K)
5615 ! CMCMAX: Maximum canopy water capacity
5616 ! CFACTR: Parameter used in the canopy inteception calculati
5617 ! SNUP: Threshold snow depth (in water equivalent m) that
5618 ! implies 100% snow cover
5619 ! LAI: Leaf area index (dimensionless)
5620 ! MAXALB: Upper bound on maximum albedo over deep snow
5621 !
5622 !-----READ IN VEGETAION PROPERTIES FROM VEGPARM.TBL
5623 !
5625 IF ( wrf_dm_on_monitor() ) THEN
5627 OPEN(19, FILE='VEGPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
5628 IF(ierr .NE. OPEN_OK ) THEN
5629 WRITE(message,FMT='(A)') &
5630 'module_sf_ruclsm.F: soil_veg_gen_parm: failure opening VEGPARM.TBL'
5631 CALL wrf_error_fatal ( message )
5632 END IF
5634 WRITE ( mess, * ) 'INPUT VEGPARM FOR ',MMINLU
5635 CALL wrf_message( mess )
5637 LUMATCH=0
5639 2000 FORMAT (A8)
5640 READ (19,'(A)') vege_parm_string
5641 outer : DO
5642 READ (19,2000,END=2002)LUTYPE
5643 READ (19,*)LUCATS,IINDEX
5645 WRITE( mess , * ) 'VEGPARM FOR ',LUTYPE,' FOUND', LUCATS,' CATEGORIES'
5646 CALL wrf_message( mess )
5648 IF(LUTYPE.NE.MMINLU)THEN ! Skip over the undesired table
5649 write ( mess , * ) 'Skipping ', LUTYPE, ' table'
5650 CALL wrf_message( mess )
5651 DO LC=1,LUCATS
5652 READ (19,*)
5653 ENDDO
5654 inner : DO ! Find the next "Vegetation Parameters"
5655 READ (19,'(A)',END=2002) vege_parm_string
5656 IF (TRIM(vege_parm_string) .EQ. "Vegetation Parameters") THEN
5657 EXIT inner
5658 END IF
5659 ENDDO inner
5660 ELSE
5661 LUMATCH=1
5662 write ( mess , * ) 'Found ', LUTYPE, ' table'
5663 CALL wrf_message( mess )
5664 EXIT outer ! Found the table, read the data
5665 END IF
5667 ENDDO outer
5669 IF (LUMATCH == 1) then
5670 write ( mess , * ) 'Reading ',LUTYPE,' table'
5671 CALL wrf_message( mess )
5672 DO LC=1,LUCATS
5673 READ (19,*)IINDEX,ALBTBL(LC),Z0TBL(LC),LEMITBL(LC),PCTBL(LC), &
5674 SHDTBL(LC),NROTBL(LC),RSTBL(LC),RGLTBL(LC), &
5675 HSTBL(LC),SNUPTBL(LC),LAITBL(LC),MAXALB(LC)
5676 ENDDO
5677 !
5678 READ (19,*)
5679 READ (19,*)TOPT_DATA
5680 READ (19,*)
5681 READ (19,*)CMCMAX_DATA
5682 READ (19,*)
5683 READ (19,*)CFACTR_DATA
5684 READ (19,*)
5685 READ (19,*)RSMAX_DATA
5686 READ (19,*)
5687 READ (19,*)BARE
5688 READ (19,*)
5689 READ (19,*)NATURAL
5690 ENDIF
5692 2002 CONTINUE
5693 CLOSE (19)
5695 IF (LUMATCH == 0) then
5696 CALL wrf_error_fatal ("Land Use Dataset '"//MMINLU//"' not found in VEGPARM.TBL.")
5697 ENDIF
5699 END IF
5701 CALL wrf_dm_bcast_string ( LUTYPE , 8 )
5702 CALL wrf_dm_bcast_integer ( LUCATS , 1 )
5703 CALL wrf_dm_bcast_integer ( IINDEX , 1 )
5704 CALL wrf_dm_bcast_integer ( LUMATCH , 1 )
5705 CALL wrf_dm_bcast_real ( ALBTBL , NLUS )
5706 CALL wrf_dm_bcast_real ( Z0TBL , NLUS )
5707 CALL wrf_dm_bcast_real ( LEMITBL , NLUS )
5708 CALL wrf_dm_bcast_real ( PCTBL , NLUS )
5709 CALL wrf_dm_bcast_real ( SHDTBL , NLUS )
5710 CALL wrf_dm_bcast_real ( NROTBL , NLUS )
5711 CALL wrf_dm_bcast_real ( RSTBL , NLUS )
5712 CALL wrf_dm_bcast_real ( RGLTBL , NLUS )
5713 CALL wrf_dm_bcast_real ( HSTBL , NLUS )
5714 CALL wrf_dm_bcast_real ( SNUPTBL , NLUS )
5715 CALL wrf_dm_bcast_real ( LAITBL , NLUS )
5716 CALL wrf_dm_bcast_real ( MAXALB , NLUS )
5717 CALL wrf_dm_bcast_real ( TOPT_DATA , 1 )
5718 CALL wrf_dm_bcast_real ( CMCMAX_DATA , 1 )
5719 CALL wrf_dm_bcast_real ( CFACTR_DATA , 1 )
5720 CALL wrf_dm_bcast_real ( RSMAX_DATA , 1 )
5721 CALL wrf_dm_bcast_integer ( BARE , 1 )
5723 !
5724 !-----READ IN SOIL PROPERTIES FROM SOILPARM.TBL
5725 !
5726 IF ( wrf_dm_on_monitor() ) THEN
5727 OPEN(19, FILE='SOILPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
5728 IF(ierr .NE. OPEN_OK ) THEN
5729 WRITE(message,FMT='(A)') &
5730 'module_sf_ruclsm.F: soil_veg_gen_parm: failure opening SOILPARM.TBL'
5731 CALL wrf_error_fatal ( message )
5732 END IF
5734 WRITE(mess,*) 'INPUT SOIL TEXTURE CLASSIFICATION = ',MMINSL
5735 CALL wrf_message( mess )
5737 LUMATCH=0
5739 READ (19,*)
5740 READ (19,2000,END=2003)SLTYPE
5741 READ (19,*)SLCATS,IINDEX
5742 IF(SLTYPE.NE.MMINSL)THEN
5743 DO LC=1,SLCATS
5744 READ (19,*) IINDEX,BB(LC),DRYSMC(LC),HC(LC),MAXSMC(LC),&
5745 REFSMC(LC),SATPSI(LC),SATDK(LC), SATDW(LC), &
5746 WLTSMC(LC), QTZ(LC)
5747 ENDDO
5748 ENDIF
5749 READ (19,*)
5750 READ (19,2000,END=2003)SLTYPE
5751 READ (19,*)SLCATS,IINDEX
5753 IF(SLTYPE.EQ.MMINSL)THEN
5754 WRITE( mess , * ) 'SOIL TEXTURE CLASSIFICATION = ',SLTYPE,' FOUND', &
5755 SLCATS,' CATEGORIES'
5756 CALL wrf_message ( mess )
5757 LUMATCH=1
5758 ENDIF
5759 IF(SLTYPE.EQ.MMINSL)THEN
5760 DO LC=1,SLCATS
5761 READ (19,*) IINDEX,BB(LC),DRYSMC(LC),HC(LC),MAXSMC(LC),&
5762 REFSMC(LC),SATPSI(LC),SATDK(LC), SATDW(LC), &
5763 WLTSMC(LC), QTZ(LC)
5764 ENDDO
5765 ENDIF
5767 2003 CONTINUE
5769 CLOSE (19)
5770 ENDIF
5772 CALL wrf_dm_bcast_integer ( LUMATCH , 1 )
5773 CALL wrf_dm_bcast_string ( SLTYPE , 8 )
5774 CALL wrf_dm_bcast_string ( MMINSL , 8 ) ! since this is reset above, see oct2 ^
5775 CALL wrf_dm_bcast_integer ( SLCATS , 1 )
5776 CALL wrf_dm_bcast_integer ( IINDEX , 1 )
5777 CALL wrf_dm_bcast_real ( BB , NSLTYPE )
5778 CALL wrf_dm_bcast_real ( DRYSMC , NSLTYPE )
5779 CALL wrf_dm_bcast_real ( HC , NSLTYPE )
5780 CALL wrf_dm_bcast_real ( MAXSMC , NSLTYPE )
5781 CALL wrf_dm_bcast_real ( REFSMC , NSLTYPE )
5782 CALL wrf_dm_bcast_real ( SATPSI , NSLTYPE )
5783 CALL wrf_dm_bcast_real ( SATDK , NSLTYPE )
5784 CALL wrf_dm_bcast_real ( SATDW , NSLTYPE )
5785 CALL wrf_dm_bcast_real ( WLTSMC , NSLTYPE )
5786 CALL wrf_dm_bcast_real ( QTZ , NSLTYPE )
5788 IF(LUMATCH.EQ.0)THEN
5789 CALL wrf_message( 'SOIl TEXTURE IN INPUT FILE DOES NOT ' )
5790 CALL wrf_message( 'MATCH SOILPARM TABLE' )
5791 CALL wrf_error_fatal ( 'INCONSISTENT OR MISSING SOILPARM FILE' )
5792 ENDIF
5794 !
5795 !-----READ IN GENERAL PARAMETERS FROM GENPARM.TBL
5796 !
5797 IF ( wrf_dm_on_monitor() ) THEN
5798 OPEN(19, FILE='GENPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
5799 IF(ierr .NE. OPEN_OK ) THEN
5800 WRITE(message,FMT='(A)') &
5801 'module_sf_ruclsm.F: soil_veg_gen_parm: failure opening GENPARM.TBL'
5802 CALL wrf_error_fatal ( message )
5803 END IF
5805 READ (19,*)
5806 READ (19,*)
5807 READ (19,*) NUM_SLOPE
5809 SLPCATS=NUM_SLOPE
5811 DO LC=1,SLPCATS
5812 READ (19,*)SLOPE_DATA(LC)
5813 ENDDO
5815 READ (19,*)
5816 READ (19,*)SBETA_DATA
5817 READ (19,*)
5818 READ (19,*)FXEXP_DATA
5819 READ (19,*)
5820 READ (19,*)CSOIL_DATA
5821 READ (19,*)
5822 READ (19,*)SALP_DATA
5823 READ (19,*)
5824 READ (19,*)REFDK_DATA
5825 READ (19,*)
5826 READ (19,*)REFKDT_DATA
5827 READ (19,*)
5828 READ (19,*)FRZK_DATA
5829 READ (19,*)
5830 READ (19,*)ZBOT_DATA
5831 READ (19,*)
5832 READ (19,*)CZIL_DATA
5833 READ (19,*)
5834 READ (19,*)SMLOW_DATA
5835 READ (19,*)
5836 READ (19,*)SMHIGH_DATA
5837 CLOSE (19)
5838 ENDIF
5840 CALL wrf_dm_bcast_integer ( NUM_SLOPE , 1 )
5841 CALL wrf_dm_bcast_integer ( SLPCATS , 1 )
5842 CALL wrf_dm_bcast_real ( SLOPE_DATA , NSLOPE )
5843 CALL wrf_dm_bcast_real ( SBETA_DATA , 1 )
5844 CALL wrf_dm_bcast_real ( FXEXP_DATA , 1 )
5845 CALL wrf_dm_bcast_real ( CSOIL_DATA , 1 )
5846 CALL wrf_dm_bcast_real ( SALP_DATA , 1 )
5847 CALL wrf_dm_bcast_real ( REFDK_DATA , 1 )
5848 CALL wrf_dm_bcast_real ( REFKDT_DATA , 1 )
5849 CALL wrf_dm_bcast_real ( FRZK_DATA , 1 )
5850 CALL wrf_dm_bcast_real ( ZBOT_DATA , 1 )
5851 CALL wrf_dm_bcast_real ( CZIL_DATA , 1 )
5852 CALL wrf_dm_bcast_real ( SMLOW_DATA , 1 )
5853 CALL wrf_dm_bcast_real ( SMHIGH_DATA , 1 )
5856 !-----------------------------------------------------------------
5857 END SUBROUTINE RUCLSM_SOILVEGPARM
5858 !-----------------------------------------------------------------
5861 SUBROUTINE SOILIN (ISLTYP, DQM, REF, PSIS, QMIN, BCLH )
5863 !--- soiltyp classification according to STATSGO(nclasses=16)
5864 !
5865 ! 1 SAND SAND
5866 ! 2 LOAMY SAND LOAMY SAND
5867 ! 3 SANDY LOAM SANDY LOAM
5868 ! 4 SILT LOAM SILTY LOAM
5869 ! 5 SILT SILTY LOAM
5870 ! 6 LOAM LOAM
5871 ! 7 SANDY CLAY LOAM SANDY CLAY LOAM
5872 ! 8 SILTY CLAY LOAM SILTY CLAY LOAM
5873 ! 9 CLAY LOAM CLAY LOAM
5874 ! 10 SANDY CLAY SANDY CLAY
5875 ! 11 SILTY CLAY SILTY CLAY
5876 ! 12 CLAY LIGHT CLAY
5877 ! 13 ORGANIC MATERIALS LOAM
5878 ! 14 WATER
5879 ! 15 BEDROCK
5880 ! Bedrock is reclassified as class 14
5881 ! 16 OTHER (land-ice)
5882 ! extra classes from Fei Chen
5883 ! 17 Playa
5884 ! 18 Lava
5885 ! 19 White Sand
5886 !
5887 !----------------------------------------------------------------------
5888 integer, parameter :: nsoilclas=19
5890 integer, intent ( in) :: isltyp
5891 real, intent ( out) :: dqm,ref,qmin,psis
5893 REAL LQMA(nsoilclas),LREF(nsoilclas),LBCL(nsoilclas), &
5894 LPSI(nsoilclas),LQMI(nsoilclas)
5896 !-- LQMA Rawls et al.[1982]
5897 ! DATA LQMA /0.417, 0.437, 0.453, 0.501, 0.486, 0.463, 0.398,
5898 ! & 0.471, 0.464, 0.430, 0.479, 0.475, 0.439, 1.0, 0.20, 0.401/
5899 !---
5900 !-- Clapp, R. and G. Hornberger, Empirical equations for some soil
5901 ! hydraulic properties, Water Resour. Res., 14,601-604,1978.
5902 !-- Clapp et al. [1978]
5903 DATA LQMA /0.395, 0.410, 0.435, 0.485, 0.485, 0.451, 0.420, &
5904 0.477, 0.476, 0.426, 0.492, 0.482, 0.451, 1.0, &
5905 0.20, 0.435, 0.468, 0.200, 0.339/
5907 !-- Clapp et al. [1978]
5908 DATA LREF /0.174, 0.179, 0.249, 0.369, 0.369, 0.314, 0.299, &
5909 0.357, 0.391, 0.316, 0.409, 0.400, 0.314, 1., &
5910 0.1, 0.249, 0.454, 0.17, 0.236/
5912 !-- Carsel and Parrish [1988]
5913 DATA LQMI/0.045, 0.057, 0.065, 0.067, 0.034, 0.078, 0.10, &
5914 0.089, 0.095, 0.10, 0.070, 0.068, 0.078, 0.0, &
5915 0.004, 0.065, 0.020, 0.004, 0.008/
5917 !-- Clapp et al. [1978]
5918 DATA LPSI/0.121, 0.090, 0.218, 0.786, 0.786, 0.478, 0.299, &
5919 0.356, 0.630, 0.153, 0.490, 0.405, 0.478, 0.0, &
5920 0.121, 0.218, 0.468, 0.069, 0.069/
5922 !-- Clapp et al. [1978]
5923 DATA LBCL/4.05, 4.38, 4.90, 5.30, 5.30, 5.39, 7.12, &
5924 7.75, 8.52, 10.40, 10.40, 11.40, 5.39, 0.0, &
5925 4.05, 4.90, 11.55, 2.79, 2.79/
5928 DQM = LQMA(ISLTYP)- &
5929 LQMI(ISLTYP)
5930 REF = LREF(ISLTYP)
5931 PSIS = - LPSI(ISLTYP)
5932 QMIN = LQMI(ISLTYP)
5933 BCLH = LBCL(ISLTYP)
5935 END SUBROUTINE SOILIN
5937 END MODULE module_sf_ruclsm