| blob | 84e445c44e48f6ab6576df1f81a326dbe8336c43 |
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 do k=1,nzs
740 smfr3d(i,k,j) = smfrkeep(k)
741 keepfr3dflag(i,k,j) = keepfr (k)
742 enddo
744 !tgs add together dew and cloud at the ground surface
745 qcg(i,j)=qcg(i,j)+dew(i,j)
747 Z0 (I,J) = ZNT (I,J)
748 SFCEXC (I,J) = TKMS
749 patmb=P8w(i,1,j)*1.e-2
750 Q2SAT=QSN(TABS,TBQ)/PATMB
751 QSFC(I,J) = QVG(I,J)/(1.+QVG(I,J))
752 ! for MYJ surface and PBL scheme
753 ! if (myj) then
754 ! MYJSFC expects QSFC as actual specific humidity at the surface
755 IF((QVATM.GE.Q2SAT*0.95).AND.QVATM.LT.qvg(I,J))THEN
756 CHKLOWQ(I,J)=0.
757 ELSE
758 CHKLOWQ(I,J)=1.
759 ENDIF
760 ! else
761 ! CHKLOWQ(I,J)=1.
762 ! endif
764 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
765 if(CHKLOWQ(I,J).eq.0.) then
766 print *,'i,j,CHKLOWQ', &
767 i,j,CHKLOWQ(I,J)
768 endif
769 ENDIF
771 ! SNOW is in [mm], SNWE is in [m]; CANWAT is in mm, CANWATR is in m
772 SNOW (i,j) = SNWE*1000.
773 SNOWH (I,J) = SNHEI
774 CANWAT (I,J) = CANWATR*1000.
775 INFILTR(I,J) = INFILTRP
777 MAVAIL (i,j) = LMAVAIL(I,J)
778 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
779 print *,' LAND, I=,J=, QFX, HFX after SFCTMP', i,j,lh(i,j),hfx(i,j)
780 ENDIF
781 !!! QFX (I,J) = LH(I,J)/LV
782 SFCEVP (I,J) = SFCEVP (I,J) + QFX (I,J) * DT
783 GRDFLX (I,J) = -1. * sflx(I,J)
785 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
786 print *,' QFX after change, LH ', i,j, QFX(i,j),LH(I,J)
787 ENDIF
788 !--- SNOWC snow cover flag
789 SNOWC(I,J)=SNOWFRAC
792 !--- get 3d soil fields
793 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
794 print *,'LAND, i,j,tso1d,soilm1d - end of time step', &
795 i,j,tso1d,soilm1d
796 ENDIF
798 !--- end of a land or sea ice point
799 ENDIF
801 ENDDO
803 ENDDO
805 !-----------------------------------------------------------------
806 END SUBROUTINE LSMRUC
807 !-----------------------------------------------------------------
811 SUBROUTINE SFCTMP (delt,ktau,conflx,i,j, &
812 !--- input variables
813 nzs,nddzs,nroot,meltfactor, &
814 ILAND,ISOIL,XLAND,IVGTYP,PRCPMS, &
815 NEWSNMS,SNWE,SNHEI,SNOWFRAC,RHOSN,RHONEWSN, &
816 PATM,TABS,QVATM,QCATM,rho, &
817 GLW,GSW,EMISS,QKMS,TKMS,PC, &
818 MAVAIL,CST,VEGFRA,ALB,ZNT, &
819 ALB_SNOW,ALB_SNOW_FREE, &
820 MYJ,SEAICE, &
821 !--- soil fixed fields
822 QWRTZ,rhocs,dqm,qmin,ref,wilt,psis,bclh,ksat, &
823 sat,cn,zsmain,zshalf,DTDZS,DTDZS2,tbq, &
824 !--- constants
825 cp,rovcp,g0,lv,stbolt,cw,c1sn,c2sn, &
826 KQWRTZ,KICE,KWT, &
827 !--- output variables
828 snweprint,snheiprint,rsm, &
829 soilm1d,ts1d,smfrkeep,keepfr,soilt,soilt1, &
830 tsnav,dew,qvg,qsg,qcg, &
831 SMELT,SNOH,SNFLX,SNOM,ACSNOW, &
832 edir1,ec1,ett1,eeta,qfx,hfx,s,sublim, &
833 evapl,prcpl,runoff1,runoff2,soilice, &
834 soiliqw,infiltr)
835 !-----------------------------------------------------------------
836 IMPLICIT NONE
837 !-----------------------------------------------------------------
839 !--- input variables
841 INTEGER, INTENT(IN ) :: i,j,nroot,ktau,nzs , &
842 nddzs !nddzs=2*(nzs-2)
844 REAL, INTENT(IN ) :: DELT,CONFLX,meltfactor
845 REAL, INTENT(IN ) :: C1SN,C2SN
846 LOGICAL, INTENT(IN ) :: myj
847 !--- 3-D Atmospheric variables
848 REAL , &
849 INTENT(IN ) :: PATM, &
850 TABS, &
851 QVATM, &
852 QCATM
853 REAL , &
854 INTENT(IN ) :: GLW, &
855 GSW, &
856 PC, &
857 VEGFRA, &
858 ALB_SNOW_FREE, &
859 SEAICE, &
860 XLAND, &
861 RHO, &
862 QKMS, &
863 TKMS
865 INTEGER, INTENT(IN ) :: IVGTYP
866 !--- 2-D variables
867 REAL , &
868 INTENT(INOUT) :: EMISS, &
869 MAVAIL, &
870 SNOWFRAC, &
871 ALB_SNOW, &
872 ALB, &
873 CST
875 !--- soil properties
876 REAL :: &
877 RHOCS, &
878 BCLH, &
879 DQM, &
880 KSAT, &
881 PSIS, &
882 QMIN, &
883 QWRTZ, &
884 REF, &
885 SAT, &
886 WILT
888 REAL, INTENT(IN ) :: CN, &
889 CW, &
890 CP, &
891 ROVCP, &
892 G0, &
893 LV, &
894 STBOLT, &
895 KQWRTZ, &
896 KICE, &
897 KWT
899 REAL, DIMENSION(1:NZS), INTENT(IN) :: ZSMAIN, &
900 ZSHALF, &
901 DTDZS2
904 REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS
906 REAL, DIMENSION(1:4001), INTENT(IN) :: TBQ
909 !--- input/output variables
910 !-------- 3-d soil moisture and temperature
911 REAL, DIMENSION( 1:nzs ) , &
912 INTENT(INOUT) :: TS1D, &
913 SOILM1D, &
914 SMFRKEEP
915 REAL, DIMENSION( 1:nzs ) , &
916 INTENT(INOUT) :: KEEPFR
919 INTEGER, INTENT(INOUT) :: ILAND,ISOIL
921 !-------- 2-d variables
922 REAL , &
923 INTENT(INOUT) :: DEW, &
924 EDIR1, &
925 EC1, &
926 ETT1, &
927 EETA, &
928 EVAPL, &
929 INFILTR, &
930 RHOSN, &
931 RHONEWSN, &
932 SUBLIM, &
933 PRCPL, &
934 QVG, &
935 QSG, &
936 QCG, &
937 QFX, &
938 HFX, &
939 S, &
940 RUNOFF1, &
941 RUNOFF2, &
942 ACSNOW, &
943 SNWE, &
944 SNHEI, &
945 SMELT, &
946 SNOM, &
947 SNOH, &
948 SNFLX, &
949 SOILT, &
950 SOILT1, &
951 TSNAV, &
952 ZNT
954 REAL, DIMENSION(1:NZS) :: &
955 tice, &
956 rhosice, &
957 capice, &
958 thdifice
959 !-------- 1-d variables
960 REAL, DIMENSION(1:NZS), INTENT(OUT) :: SOILICE, &
961 SOILIQW
965 REAL, INTENT(OUT) :: RSM, &
966 SNWEPRINT, &
967 SNHEIPRINT
968 !--- Local variables
970 INTEGER :: K,ILNB
972 REAL :: BSN, XSN , &
973 RAINF, SNTH, NEWSN, PRCPMS, NEWSNMS , &
974 T3, UPFLUX, XINET
975 REAL :: snhei_crit, keep_snow_albedo
977 REAL :: RNET,GSWNEW,EMISSN,ZNTSN
978 REAL :: VEGFRAC
979 real :: cice, albice, albsn
981 !-----------------------------------------------------------------
982 integer, parameter :: ilsnow=99
984 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
985 print *,' in SFCTMP',i,j,nzs,nddzs,nroot, &
986 SNWE,RHOSN,SNOM,SMELT,TS1D
987 ENDIF
989 NEWSN=0.
990 RAINF = 0.
991 RSM=0.
992 INFILTR=0.
993 VEGFRAC=0.01*VEGFRA
994 if(VEGFRAC.le.0.01) VEGFRAC=0.
995 !---initialize local arrays for sea ice
996 do k=1,nzs
997 tice(k) = 0.
998 rhosice(k) = 0.
999 cice = 0.
1000 capice(k) = 0.
1001 thdifice(k) = 0.
1002 enddo
1004 GSWnew=GSW
1005 ALBice=ALB_SNOW_FREE
1006 !--- sea ice properties
1007 !--- N.N Zubov "Arctic Ice"
1008 !--- no salinity dependence because we consider the ice pack
1009 !--- to be old and to have low salinity (0.0002)
1010 if(SEAICE.ge.0.5) then
1011 do k=1,nzs
1012 tice(k) = ts1d(k) - 273.15
1013 rhosice(k) = 917.6/(1-0.000165*tice(k))
1014 cice = 2115.85 +7.7948*tice(k)
1015 capice(k) = cice*rhosice(k)
1016 thdifice(k) = 2.260872/capice(k)
1017 enddo
1018 !-- SEA ICE ALB dependence on ice temperature. When ice temperature is
1019 !-- below critical value of -10C - no change to albedo.
1020 !-- If temperature is higher that -10C then albedo is decreasing.
1021 !-- The minimum albedo at t=0C for ice is 0.1 less.
1022 GSWNEW=GSW/(1.-ALB)
1023 ALBice = MIN(ALB_SNOW_FREE,MAX(ALB_SNOW_FREE - 0.1, &
1024 ALB_SNOW_FREE - 0.1*(tice(1)+10.)/10. ))
1025 GSWNEW=GSW*(1.-ALBice)
1026 endif
1028 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
1029 print *,'I,J,KTAU,QKMS,TKMS', i,j,ktau,qkms,tkms
1030 print *,'GSW,GSWnew,GLW,SOILT,EMISS,ALB,ALBice,SNWE',&
1031 GSW,GSWnew,GLW,SOILT,EMISS,ALB,ALBice,SNWE
1032 ENDIF
1034 SNHEI = SNWE * 1000. / RHOSN
1035 !--------------
1036 T3 = STBOLT*SOILT*SOILT*SOILT
1037 UPFLUX = T3 *SOILT
1038 XINET = EMISS*(GLW-UPFLUX)
1039 RNET = GSWnew + XINET
1041 !Calculate the amount (m) of fresh snow
1043 if(snhei.gt.0.0081*1.e3/rhosn) then
1044 !*** Correct snow density for current temperature (Koren et al. 1999)
1045 BSN=delt/3600.*c1sn*exp(0.08*tsnav-c2sn*rhosn*1.e-3)
1046 if(bsn*snwe*100..lt.1.e-4) goto 777
1047 XSN=rhosn*(exp(bsn*snwe*100.)-1.)/(bsn*snwe*100.)
1048 rhosn=MIN(MAX(100.,XSN),400.)
1049 ! rhosn=MIN(MAX(50.,XSN),400.)
1050 777 continue
1052 else
1053 rhosn =200.
1054 rhonewsn =200.
1055 endif
1057 newsn=newsnms*delt
1058 !---- ACSNOW - accumulation of snow water [m]
1059 acsnow=acsnow+newsn
1061 IF(NEWSN.GE.1.E-8) THEN
1062 !*** Calculate fresh snow density (t > -15C, else MIN value)
1063 !*** Eq. 10 from Koren et al. (1999)
1065 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
1066 print *, 'THERE IS NEW SNOW, newsn', newsn
1067 ENDIF
1068 if(tabs.lt.258.15) then
1069 ! rhonewsn=50.
1070 rhonewsn=100.
1072 else
1073 rhonewsn=1.e3*max((0.10+0.0017*(Tabs-273.15+15.)**1.5) &
1074 , 0.10)
1075 ! rhonewsn=1.e3*max((0.05+0.0017*(Tabs-273.15+15.)**1.5) &
1076 ! , 0.05)
1077 rhonewsn=MIN(rhonewsn,400.)
1078 ! rhonewsn=100.
1079 endif
1082 !*** Define average snow density of the snow pack considering
1083 !*** the amount of fresh snow (eq. 9 in Koren et al.(1999)
1084 !*** without snow melt )
1085 xsn=(rhosn*snwe+rhonewsn*newsn)/ &
1086 (snwe+newsn)
1087 rhosn=MIN(MAX(100.,XSN),400.)
1088 ! rhosn=MIN(MAX(50.,XSN),400.)
1090 snwe=snwe+newsn
1091 snhei=snwe*1.E3/rhosn
1092 NEWSN=NEWSN*1.E3/rhosn
1093 endif
1095 IF(PRCPMS.NE.0.) THEN
1097 ! PRCPMS is liquid precipitation rate
1098 ! RAINF is a flag used for calculation of rain water
1099 ! heat content contribution into heat budget equation. Rain's temperature
1100 ! is set equal to air temperature at the first atmospheric
1101 ! level.
1103 RAINF=1.
1104 ENDIF
1106 IF(SNHEI.GT.0.0) THEN
1107 !--- Set of surface parameters should be changed to snow values for grid
1108 !--- points where the snow cover exceeds snow threshold of 2 cm
1110 SNHEI_CRIT=0.01601*1.e3/rhosn
1111 SNOWFRAC=MIN(1.,SNHEI/(2.*SNHEI_CRIT))
1112 !tgs NEW - low limit on snow fraction
1113 if(SNOWFRAC.lt.0.01) snowfrac=0.
1114 !--- EMISS = 0.91 for snow
1115 EMISS = EMISS*(1.-snowfrac)+0.91*snowfrac
1117 KEEP_SNOW_ALBEDO = 0.
1118 IF (NEWSN.GT.0.) KEEP_SNOW_ALBEDO = 1.
1120 !--- GSW in-coming solar
1121 GSWNEW=GSW/(1.-ALB)
1123 IF(SEAICE .LT. 0.5) THEN
1124 !----- SNOW on soil
1125 !-- ALB dependence on snow depth
1126 ALBsn = MAX(keep_snow_albedo*alb_snow, &
1127 MIN((alb_snow_free + &
1128 (alb_snow - alb_snow_free) * &
1129 (snhei/(2.*SNHEI_CRIT))), alb_snow))
1131 !-- ALB dependence on snow temperature. When snow temperature is
1132 !-- below critical value of -10C - no change to albedo.
1133 !-- If temperature is higher that -10C then albedo is decreasing.
1134 !-- The minimum albedo at t=0C for snow on land is 15% less than
1135 !-- albedo of temperatures below -10C.
1136 if(albsn.lt.alb_snow)then
1137 ALB=ALBsn
1138 else
1139 ALB = MIN(ALBSN,MAX(ALBSN - 0.15*(soilt - 263.15)/ &
1140 (273.15-263.15), ALBSN - 0.15))
1141 endif
1142 ELSE
1143 !----- SNOW on ice
1144 ALBsn = MAX(keep_snow_albedo*alb_snow, &
1145 MIN((albice + (alb_snow - albice) * &
1146 (snhei/(2.*SNHEI_CRIT))), alb_snow))
1148 !-- ALB dependence on snow temperature. When snow temperature is
1149 !-- below critical value of -10C - no change to albedo.
1150 !-- If temperature is higher that -10C then albedo is decreasing.
1151 !-- The minimum albedo at t=0C for snow on ice is 0.15 less.
1152 if(albsn.lt.alb_snow)then
1153 ALB=ALBsn
1154 else
1155 ALB = MIN(ALBSN,MAX(ALBSN - 0.15*(soilt - 263.15)/ &
1156 (273.15-263.15), ALBSN - 0.15))
1157 endif
1159 ENDIF
1161 !--- recompute absorbed solar radiation and net radiation
1162 !--- for new value of albedo
1163 gswnew=gswnew*(1.-alb)
1165 XINET = EMISS*(GLW-UPFLUX)
1166 RNET = GSWnew + XINET
1168 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
1169 print *,'I,J,GSW,GSWnew,GLW,UPFLUX,ALB',&
1170 i,j,GSW,GSWnew,GLW,UPFLUX,ALB
1171 ENDIF
1174 if (SEAICE .LT. 0.5) then
1175 ! LAND
1176 CALL SNOWSOIL ( & !--- input variables
1177 i,j,isoil,delt,ktau,conflx,nzs,nddzs,nroot, &
1178 meltfactor,rhonewsn, & ! new
1179 ILAND,PRCPMS,RAINF,NEWSN,snhei,SNWE,snowfrac, &
1180 RHOSN,PATM,QVATM,QCATM, &
1181 GLW,GSWnew,EMISS,RNET,IVGTYP, &
1182 QKMS,TKMS,PC,CST, &
1183 RHO,VEGFRAC,ALB,ZNT, &
1184 MYJ, &
1185 !--- soil fixed fields
1186 QWRTZ,rhocs,dqm,qmin,ref,wilt,psis,bclh,ksat, &
1187 sat,cn,zsmain,zshalf,DTDZS,DTDZS2,tbq, &
1188 !--- constants
1189 lv,CP,rovcp,G0,cw,stbolt,tabs, &
1190 KQWRTZ,KICE,KWT, &
1191 !--- output variables
1192 ilnb,snweprint,snheiprint,rsm, &
1193 soilm1d,ts1d,smfrkeep,keepfr, &
1194 dew,soilt,soilt1,tsnav,qvg,qsg,qcg, &
1195 SMELT,SNOH,SNFLX,SNOM,edir1,ec1,ett1,eeta, &
1196 qfx,hfx,s,sublim,prcpl,runoff1,runoff2, &
1197 mavail,soilice,soiliqw,infiltr )
1198 else
1199 ! SEA ICE
1200 CALL SNOWSEAICE ( &
1201 i,j,isoil,delt,ktau,conflx,nzs,nddzs, &
1202 meltfactor,rhonewsn, & ! new
1203 ILAND,PRCPMS,RAINF,NEWSN,snhei,SNWE,snowfrac, &
1204 RHOSN,PATM,QVATM,QCATM, &
1205 GLW,GSWnew,EMISS,RNET, &
1206 QKMS,TKMS,RHO, &
1207 !--- sea ice parameters
1208 ALB,ZNT, &
1209 tice,rhosice,capice,thdifice, &
1210 zsmain,zshalf,DTDZS,DTDZS2,tbq, &
1211 !--- constants
1212 lv,CP,rovcp,cw,stbolt,tabs, &
1213 !--- output variables
1214 ilnb,snweprint,snheiprint,rsm,ts1d, &
1215 dew,soilt,soilt1,tsnav,qvg,qsg,qcg, &
1216 SMELT,SNOH,SNFLX,SNOM,eeta, &
1217 qfx,hfx,s,sublim,prcpl &
1218 )
1219 edir1 = eeta
1220 ec1 = 0.
1221 ett1 = 0.
1222 runoff1 = smelt
1223 runoff2 = 0.
1224 mavail = 1.
1225 infiltr=0.
1226 cst=0.
1227 do k=1,nzs
1228 soilm1d(k)=1.
1229 soiliqw(k)=0.
1230 soilice(k)=1.
1231 smfrkeep(k)=1.
1232 keepfr(k)=0.
1233 enddo
1234 endif
1236 if(snhei.eq.0.) then
1237 !--- all snow is melted
1238 alb=alb_snow_free
1239 endif
1241 ELSE
1242 !--- no snow
1243 snheiprint=0.
1244 snweprint=0.
1246 if(SEAICE .LT. 0.5) then
1247 ! LAND
1248 CALL SOIL( &
1249 !--- input variables
1250 i,j,iland,isoil,delt,ktau,conflx,nzs,nddzs,nroot, &
1251 PRCPMS,RAINF,PATM,QVATM,QCATM,GLW,GSWnew, &
1252 EMISS,RNET,QKMS,TKMS,PC,cst,rho,vegfrac, &
1253 !--- soil fixed fields
1254 QWRTZ,rhocs,dqm,qmin,ref,wilt, &
1255 psis,bclh,ksat,sat,cn, &
1256 zsmain,zshalf,DTDZS,DTDZS2,tbq, &
1257 !--- constants
1258 lv,CP,rovcp,G0,cw,stbolt,tabs, &
1259 KQWRTZ,KICE,KWT, &
1260 !--- output variables
1261 soilm1d,ts1d,smfrkeep,keepfr, &
1262 dew,soilt,qvg,qsg,qcg,edir1,ec1, &
1263 ett1,eeta,qfx,hfx,s,evapl,prcpl,runoff1, &
1264 runoff2,mavail,soilice,soiliqw, &
1265 infiltr)
1266 else
1267 ! SEA ICE
1268 CALL SICE( &
1269 !--- input variables
1270 i,j,iland,isoil,delt,ktau,conflx,nzs,nddzs,nroot, &
1271 PRCPMS,RAINF,PATM,QVATM,QCATM,GLW,GSWnew, &
1272 EMISS,RNET,QKMS,TKMS,rho, &
1273 !--- sea ice parameters
1274 tice,rhosice,capice,thdifice, &
1275 zsmain,zshalf,DTDZS,DTDZS2,tbq, &
1276 !--- constants
1277 lv,CP,rovcp,cw,stbolt,tabs, &
1278 !--- output variables
1279 ts1d,dew,soilt,qvg,qsg,qcg, &
1280 eeta,qfx,hfx,s,evapl,prcpl &
1281 )
1282 edir1 = eeta
1283 ec1 = 0.
1284 ett1 = 0.
1285 runoff1 = prcpms
1286 runoff2 = 0.
1287 mavail = 1.
1288 infiltr=0.
1289 cst=0.
1290 do k=1,nzs
1291 soilm1d(k)=1.
1292 soiliqw(k)=0.
1293 soilice(k)=1.
1294 smfrkeep(k)=1.
1295 keepfr(k)=0.
1296 enddo
1297 endif
1299 ENDIF
1300 ! ENDIF
1302 !
1304 ! RETURN
1305 ! END
1306 !---------------------------------------------------------------
1307 END SUBROUTINE SFCTMP
1308 !---------------------------------------------------------------
1311 FUNCTION QSN(TN,T)
1312 !****************************************************************
1313 REAL, DIMENSION(1:4001), INTENT(IN ) :: T
1314 REAL, INTENT(IN ) :: TN
1316 REAL QSN, R,R1,R2
1317 INTEGER I
1319 R=(TN-173.15)/.05+1.
1320 I=INT(R)
1321 IF(I.GE.1) goto 10
1322 I=1
1323 R=1.
1324 10 IF(I.LE.4000) GOTO 20
1325 I=4000
1326 R=4001.
1327 20 R1=T(I)
1328 R2=R-I
1329 QSN=(T(I+1)-R1)*R2 + R1
1330 ! print *,' in QSN, I,R,R1,R2,T(I+1),TN, QSN', I,R,r1,r2,t(i+1),tn,QSN
1331 ! RETURN
1332 ! END
1333 !-----------------------------------------------------------------------
1334 END FUNCTION QSN
1335 !------------------------------------------------------------------------
1338 SUBROUTINE SOIL ( &
1339 !--- input variables
1340 i,j,iland,isoil,delt,ktau,conflx,nzs,nddzs,nroot,&
1341 PRCPMS,RAINF,PATM,QVATM,QCATM, &
1342 GLW,GSW,EMISS,RNET, &
1343 QKMS,TKMS,PC,cst,rho,vegfrac, &
1344 !--- soil fixed fields
1345 QWRTZ,rhocs,dqm,qmin,ref,wilt,psis,bclh,ksat, &
1346 sat,cn,zsmain,zshalf,DTDZS,DTDZS2,tbq, &
1347 !--- constants
1348 xlv,CP,rovcp,G0_P,cw,stbolt,TABS, &
1349 KQWRTZ,KICE,KWT, &
1350 !--- output variables
1351 soilmois,tso,smfrkeep,keepfr, &
1352 dew,soilt,qvg,qsg,qcg, &
1353 edir1,ec1,ett1,eeta,qfx,hfx,s,evapl, &
1354 prcpl,runoff1,runoff2,mavail,soilice, &
1355 soiliqw,infiltrp)
1357 !*************************************************************
1358 ! Energy and moisture budget for vegetated surfaces
1359 ! without snow, heat diffusion and Richards eqns. in
1360 ! soil
1361 !
1362 ! DELT - time step (s)
1363 ! ktau - numver of time step
1364 ! CONFLX - depth of constant flux layer (m)
1365 ! J,I - the location of grid point
1366 ! IME, JME, KME, NZS - dimensions of the domain
1367 ! NROOT - number of levels within the root zone
1368 ! PRCPMS - precipitation rate in m/s
1369 ! PATM - pressure [bar]
1370 ! QVATM,QCATM - cloud and water vapor mixing ratio (kg/kg)
1371 ! at the first atm. level
1372 ! GLW, GSW - incoming longwave and absorbed shortwave
1373 ! radiation at the surface (W/m^2)
1374 ! EMISS,RNET - emissivity of the ground surface (0-1) and net
1375 ! radiation at the surface (W/m^2)
1376 ! QKMS - exchange coefficient for water vapor in the
1377 ! surface layer (m/s)
1378 ! TKMS - exchange coefficient for heat in the surface
1379 ! layer (m/s)
1380 ! PC - plant coefficient (resistance) (0-1)
1381 ! RHO - density of atmosphere near sueface (kg/m^3)
1382 ! VEGFRAC - greeness fraction
1383 ! RHOCS - volumetric heat capacity of dry soil
1384 ! DQM, QMIN - porosity minus residual soil moisture QMIN (m^3/m^3)
1385 ! REF, WILT - field capacity soil moisture and the
1386 ! wilting point (m^3/m^3)
1387 ! PSIS - matrix potential at saturation (m)
1388 ! BCLH - exponent for Clapp-Hornberger parameterization
1389 ! KSAT - saturated hydraulic conductivity (m/s)
1390 ! SAT - maximum value of water intercepted by canopy (m)
1391 ! CN - exponent for calculation of canopy water
1392 ! ZSMAIN - main levels in soil (m)
1393 ! ZSHALF - middle of the soil layers (m)
1394 ! DTDZS,DTDZS2 - dt/(2.*dzshalf*dzmain) and dt/dzshalf in soil
1395 ! TBQ - table to define saturated mixing ration
1396 ! of water vapor for given temperature and pressure
1397 ! SOILMOIS,TSO - soil moisture (m^3/m^3) and temperature (K)
1398 ! DEW - dew in kg/m^2s
1399 ! SOILT - skin temperature (K)
1400 ! QSG,QVG,QCG - saturated mixing ratio, mixing ratio of
1401 ! water vapor and cloud at the ground
1402 ! surface, respectively (kg/kg)
1403 ! EDIR1, EC1, ETT1, EETA - direct evaporation, evaporation of
1404 ! canopy water, transpiration in kg m-2 s-1 and total
1405 ! evaporation in m s-1.
1406 ! QFX, HFX - latent and sensible heat fluxes (W/m^2)
1407 ! S - soil heat flux in the top layer (W/m^2)
1408 ! RUNOFF - surface runoff (m/s)
1409 ! RUNOFF2 - underground runoff (m)
1410 ! MAVAIL - moisture availability in the top soil layer (0-1)
1411 ! INFILTRP - infiltration flux from the top of soil domain (m/s)
1412 !
1413 !*****************************************************************
1414 IMPLICIT NONE
1415 !-----------------------------------------------------------------
1417 !--- input variables
1419 INTEGER, INTENT(IN ) :: nroot,ktau,nzs , &
1420 nddzs !nddzs=2*(nzs-2)
1421 INTEGER, INTENT(IN ) :: i,j,iland,isoil
1422 REAL, INTENT(IN ) :: DELT,CONFLX
1423 !--- 3-D Atmospheric variables
1424 REAL, &
1425 INTENT(IN ) :: PATM, &
1426 QVATM, &
1427 QCATM
1428 !--- 2-D variables
1429 REAL, &
1430 INTENT(IN ) :: GLW, &
1431 GSW, &
1432 EMISS, &
1433 RHO, &
1434 PC, &
1435 VEGFRAC, &
1436 QKMS, &
1437 TKMS
1439 !--- soil properties
1440 REAL, &
1441 INTENT(IN ) :: RHOCS, &
1442 BCLH, &
1443 DQM, &
1444 KSAT, &
1445 PSIS, &
1446 QMIN, &
1447 QWRTZ, &
1448 REF, &
1449 WILT
1451 REAL, INTENT(IN ) :: CN, &
1452 CW, &
1453 KQWRTZ, &
1454 KICE, &
1455 KWT, &
1456 XLV, &
1457 g0_p
1460 REAL, DIMENSION(1:NZS), INTENT(IN) :: ZSMAIN, &
1461 ZSHALF, &
1462 DTDZS2
1464 REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS
1466 REAL, DIMENSION(1:4001), INTENT(IN) :: TBQ
1469 !--- input/output variables
1470 !-------- 3-d soil moisture and temperature
1471 REAL, DIMENSION( 1:nzs ) , &
1472 INTENT(INOUT) :: TSO, &
1473 SOILMOIS, &
1474 SMFRKEEP
1476 REAL, DIMENSION( 1:nzs ) , &
1477 INTENT(INOUT) :: KEEPFR
1479 !-------- 2-d variables
1480 REAL, &
1481 INTENT(INOUT) :: DEW, &
1482 CST, &
1483 EDIR1, &
1484 EC1, &
1485 ETT1, &
1486 EETA, &
1487 EVAPL, &
1488 PRCPL, &
1489 MAVAIL, &
1490 QVG, &
1491 QSG, &
1492 QCG, &
1493 RNET, &
1494 QFX, &
1495 HFX, &
1496 S, &
1497 SAT, &
1498 RUNOFF1, &
1499 RUNOFF2, &
1500 SOILT
1502 !-------- 1-d variables
1503 REAL, DIMENSION(1:NZS), INTENT(OUT) :: SOILICE, &
1504 SOILIQW
1506 !--- Local variables
1508 REAL :: INFILTRP, transum , &
1509 RAINF, PRCPMS , &
1510 TABS, T3, UPFLUX, XINET
1511 REAL :: CP,rovcp,G0,LV,STBOLT,xlmelt,dzstop , &
1512 can,epot,fac,fltot,ft,fq,hft , &
1513 q1,ras,rhoice,sph , &
1514 trans,zn,ci,cvw,tln,tavln,pi , &
1515 DD1,CMC2MS,DRYCAN,WETCAN , &
1516 INFMAX,RIW
1517 REAL, DIMENSION(1:NZS) :: transp,cap,diffu,hydro , &
1518 thdif,tranf,tav,soilmoism , &
1519 soilicem,soiliqwm,detal , &
1520 fwsat,lwsat,told,smold
1522 REAL :: drip
1524 INTEGER :: nzs1,nzs2,k
1526 !-----------------------------------------------------------------
1528 !-- define constants
1529 ! STBOLT=5.670151E-8
1530 RHOICE=900.
1531 CI=RHOICE*2100.
1532 XLMELT=3.35E+5
1533 cvw=cw
1535 SAT=0.0004
1536 prcpl=prcpms
1538 !--- Initializing local arrays
1539 DO K=1,NZS
1540 TRANSP (K)=0.
1541 soilmoism(k)=0.
1542 soilice (k)=0.
1543 soiliqw (k)=0.
1544 soilicem (k)=0.
1545 soiliqwm (k)=0.
1546 lwsat (k)=0.
1547 fwsat (k)=0.
1548 tav (k)=0.
1549 cap (k)=0.
1550 thdif (k)=0.
1551 diffu (k)=0.
1552 hydro (k)=0.
1553 tranf (k)=0.
1554 detal (k)=0.
1555 told (k)=0.
1556 smold (k)=0.
1557 ENDDO
1559 NZS1=NZS-1
1560 NZS2=NZS-2
1561 dzstop=1./(zsmain(2)-zsmain(1))
1562 RAS=RHO*1.E-3
1563 RIW=rhoice*1.e-3
1565 !--- Computation of volumetric content of ice in soil
1567 DO K=1,NZS
1568 !- main levels
1569 tln=log(tso(k)/273.15)
1570 if(tln.lt.0.) then
1571 soiliqw(k)=(dqm+qmin)*(XLMELT* &
1572 (tso(k)-273.15)/tso(k)/9.81/psis) &
1573 **(-1./bclh)-qmin
1574 soiliqw(k)=max(0.,soiliqw(k))
1575 soiliqw(k)=min(soiliqw(k),soilmois(k))
1576 soilice(k)=(soilmois(k)-soiliqw(k))/RIW
1578 !---- melting and freezing is balanced, soil ice cannot increase
1579 if(keepfr(k).eq.1.) then
1580 soilice(k)=min(soilice(k),smfrkeep(k))
1581 soiliqw(k)=max(0.,soilmois(k)-soilice(k)*riw)
1582 endif
1584 else
1585 soilice(k)=0.
1586 soiliqw(k)=soilmois(k)
1587 endif
1589 ENDDO
1591 DO K=1,NZS1
1592 !- middle of soil layers
1593 tav(k)=0.5*(tso(k)+tso(k+1))
1594 soilmoism(k)=0.5*(soilmois(k)+soilmois(k+1))
1595 tavln=log(tav(k)/273.15)
1597 if(tavln.lt.0.) then
1598 soiliqwm(k)=(dqm+qmin)*(XLMELT* &
1599 (tav(k)-273.15)/tav(k)/9.81/psis) &
1600 **(-1./bclh)-qmin
1601 fwsat(k)=dqm-soiliqwm(k)
1602 lwsat(k)=soiliqwm(k)+qmin
1603 soiliqwm(k)=max(0.,soiliqwm(k))
1604 soiliqwm(k)=min(soiliqwm(k), soilmoism(k))
1605 soilicem(k)=(soilmoism(k)-soiliqwm(k))/riw
1606 !---- melting and freezing is balanced, soil ice cannot increase
1607 if(keepfr(k).eq.1.) then
1608 soilicem(k)=min(soilicem(k), &
1609 0.5*(smfrkeep(k)+smfrkeep(k+1)))
1610 soiliqwm(k)=max(0.,soilmoism(k)-soilicem(k)*riw)
1611 fwsat(k)=dqm-soiliqwm(k)
1612 lwsat(k)=soiliqwm(k)+qmin
1613 endif
1615 else
1616 soilicem(k)=0.
1617 soiliqwm(k)=soilmoism(k)
1618 lwsat(k)=dqm+qmin
1619 fwsat(k)=0.
1620 endif
1622 ENDDO
1624 do k=1,nzs
1625 if(soilice(k).gt.0.) then
1626 smfrkeep(k)=soilice(k)
1627 else
1628 smfrkeep(k)=soilmois(k)/riw
1629 endif
1630 enddo
1632 !******************************************************************
1633 ! SOILPROP computes thermal diffusivity, and diffusional and
1634 ! hydraulic condeuctivities
1635 !******************************************************************
1636 CALL SOILPROP( &
1637 !--- input variables
1638 nzs,fwsat,lwsat,tav,keepfr, &
1639 soilmois,soiliqw,soilice, &
1640 soilmoism,soiliqwm,soilicem, &
1641 !--- soil fixed fields
1642 QWRTZ,rhocs,dqm,qmin,psis,bclh,ksat, &
1643 !--- constants
1644 riw,xlmelt,CP,G0_P,cvw,ci, &
1645 kqwrtz,kice,kwt, &
1646 !--- output variables
1647 thdif,diffu,hydro,cap)
1649 !********************************************************************
1650 !--- CALCULATION OF CANOPY WATER (Smirnova et al., 1996, EQ.16) AND DEW
1652 DRIP=0.
1653 DD1=0.
1655 FQ=QKMS
1657 Q1=-QKMS*RAS*(QVATM - QSG)
1659 DEW=0.
1660 IF(QVATM.GE.QSG)THEN
1661 DEW=FQ*(QVATM-QSG)
1662 ENDIF
1663 IF(DEW.NE.0.)THEN
1664 DD1=CST+DELT*(PRCPMS +DEW*RAS)*vegfrac
1665 ELSE
1666 DD1=CST+ &
1667 DELT*(PRCPMS+RAS*FQ*(QVATM-QSG) &
1668 *(CST/SAT)**CN)*vegfrac
1669 ENDIF
1671 IF(DD1.LT.0.) DD1=0.
1672 if(vegfrac.eq.0.)then
1673 cst=0.
1674 drip=0.
1675 endif
1676 IF (vegfrac.GT.0.) THEN
1677 CST=DD1
1678 IF(CST.GT.SAT) THEN
1679 CST=SAT
1680 DRIP=DD1-SAT
1681 ENDIF
1682 ENDIF
1684 !--- WETCAN is the fraction of vegetated area covered by canopy
1685 !--- water, and DRYCAN is the fraction of vegetated area where
1686 !--- transpiration may take place.
1688 WETCAN=(CST/SAT)**CN
1689 DRYCAN=1.-WETCAN
1691 !**************************************************************
1692 ! TRANSF computes transpiration function
1693 !**************************************************************
1694 CALL TRANSF( &
1695 !--- input variables
1696 nzs,nroot,soiliqw,tabs, &
1697 !--- soil fixed fields
1698 dqm,qmin,ref,wilt,zshalf, &
1699 !--- output variables
1700 tranf,transum)
1703 !--- Save soil temp and moisture from the beginning of time step
1704 do k=1,nzs
1705 told(k)=tso(k)
1706 smold(k)=soilmois(k)
1707 enddo
1709 !**************************************************************
1710 ! SOILTEMP soilves heat budget and diffusion eqn. in soil
1711 !**************************************************************
1713 CALL SOILTEMP( &
1714 !--- input variables
1715 i,j,iland,isoil, &
1716 delt,ktau,conflx,nzs,nddzs,nroot, &
1717 PRCPMS,RAINF, &
1718 PATM,TABS,QVATM,QCATM,EMISS,RNET, &
1719 QKMS,TKMS,PC,rho,vegfrac, &
1720 thdif,cap,drycan,wetcan, &
1721 transum,dew,mavail, &
1722 !--- soil fixed fields
1723 dqm,qmin,bclh,zsmain,zshalf,DTDZS,tbq, &
1724 !--- constants
1725 xlv,CP,G0_P,cvw,stbolt, &
1726 !--- output variables
1727 tso,soilt,qvg,qsg,qcg)
1729 !************************************************************************
1731 !--- CALCULATION OF DEW USING NEW VALUE OF QSG OR TRANSP IF NO DEW
1732 ETT1=0.
1733 DEW=0.
1735 IF(QVATM.GE.QSG)THEN
1736 DEW=QKMS*(QVATM-QSG)
1737 DO K=1,NZS
1738 TRANSP(K)=0.
1739 ENDDO
1740 ELSE
1741 DO K=1,NROOT
1742 TRANSP(K)=VEGFRAC*RAS*QKMS* &
1743 (QVATM-QSG)* &
1744 PC*TRANF(K)*DRYCAN/ZSHALF(NROOT+1)
1745 IF(TRANSP(K).GT.0.) TRANSP(K)=0.
1746 ETT1=ETT1-TRANSP(K)
1747 ENDDO
1748 DO k=nroot+1,nzs
1749 transp(k)=0.
1750 enddo
1751 ENDIF
1753 !-- Recalculating of volumetric content of frozen water in soil
1754 DO K=1,NZS
1755 !- main levels
1756 tln=log(tso(k)/273.15)
1757 if(tln.lt.0.) then
1758 soiliqw(k)=(dqm+qmin)*(XLMELT* &
1759 (tso(k)-273.15)/tso(k)/9.81/psis) &
1760 **(-1./bclh)-qmin
1761 soiliqw(k)=max(0.,soiliqw(k))
1762 soiliqw(k)=min(soiliqw(k),soilmois(k))
1763 soilice(k)=(soilmois(k)-soiliqw(k))/riw
1764 !---- melting and freezing is balanced, soil ice cannot increase
1765 if(keepfr(k).eq.1.) then
1766 soilice(k)=min(soilice(k),smfrkeep(k))
1767 soiliqw(k)=max(0.,soilmois(k)-soilice(k)*riw)
1768 endif
1770 else
1771 soilice(k)=0.
1772 soiliqw(k)=soilmois(k)
1773 endif
1774 ENDDO
1776 !*************************************************************************
1777 ! SOILMOIST solves moisture budget (Smirnova et al., 1996, EQ.22,28)
1778 ! and Richards eqn.
1779 !*************************************************************************
1780 CALL SOILMOIST ( &
1781 !-- input
1782 delt,nzs,nddzs,DTDZS,DTDZS2,RIW, &
1783 zsmain,zshalf,diffu,hydro, &
1784 QSG,QVG,QCG,QCATM,QVATM,-PRCPMS, &
1785 QKMS,TRANSP,DRIP,DEW,0.,SOILICE,VEGFRAC, &
1786 ! QKMS,TRANSP,DRIP,DEW,0.,SOILICE,VEGFRAC, &
1787 !-- soil properties
1788 DQM,QMIN,REF,KSAT,RAS,INFMAX, &
1789 !-- output
1790 SOILMOIS,SOILIQW,MAVAIL,RUNOFF1, &
1791 RUNOFF2,INFILTRP)
1793 !--- KEEPFR is 1 when the temperature and moisture in soil
1794 !--- are both increasing. In this case soil ice should not
1795 !--- be increasing according to the freezing curve.
1796 !--- Some part of ice is melted, but additional water is
1797 !--- getting frozen. Thus, only structure of frozen soil is
1798 !--- changed, and phase changes are not affecting the heat
1799 !--- transfer. This situation may happen when it rains on the
1800 !--- frozen soil.
1802 do k=1,nzs
1803 if (soilice(k).gt.0.) then
1804 if(tso(k).gt.told(k).and.soilmois(k).gt.smold(k)) then
1805 keepfr(k)=1.
1806 else
1807 keepfr(k)=0.
1808 endif
1809 endif
1810 enddo
1811 !--- THE DIAGNOSTICS OF SURFACE FLUXES
1813 T3 = STBOLT*SOILT*SOILT*SOILT
1814 UPFLUX = T3 *SOILT
1815 XINET = EMISS*(GLW-UPFLUX)
1816 RNET = GSW + XINET
1817 HFT=-TKMS*CP*RHO*(TABS-SOILT) &
1818 *(P1000mb*0.00001/Patm)**ROVCP
1819 Q1=-QKMS*RAS*(QVATM - QSG)
1820 IF (Q1.LE.0.) THEN
1821 ! --- condensation
1822 EC1=0.
1823 EDIR1=0.
1824 ETT1=0.
1825 !-- moisture flux for coupling with PBL
1826 EETA=-QKMS*RAS*(QVATM/(1.+QVATM) - QSG/(1.+QSG))*1.E3
1827 QFX= XLV*EETA
1828 !-- actual moisture flux from RUC LSM
1829 EETA= RHO*DEW
1830 ELSE
1831 ! --- evaporation
1832 EDIR1 =-(1.-vegfrac)*QKMS*RAS* &
1833 (QVATM-QVG)
1834 EC1 = Q1 * WETCAN
1835 CMC2MS=CST/DELT
1836 if(EC1.gt.CMC2MS) cst=0.
1837 EC1=MIN(CMC2MS,EC1)*vegfrac
1838 !-- moisture flux for coupling with PBL
1839 EETA=-QKMS*RAS*(QVATM/(1.+QVATM) - QVG/(1.+QVG))*1.E3
1840 ! to convert from kg m-2 s-1 to m s-1: 1/rho water=1.e-3************
1841 QFX= XLV * EETA
1842 !-- actual moisture flux from RUC LSM
1843 EETA = (EDIR1 + EC1 + ETT1)*1.E3
1844 ENDIF
1845 EVAPL=EETA
1846 S=THDIF(1)*CAP(1)*DZSTOP*(TSO(1)-TSO(2))
1847 HFX=HFT
1848 FLTOT=RNET-HFT-QFX-S
1850 222 CONTINUE
1852 1123 FORMAT(I5,8F12.3)
1853 1133 FORMAT(I7,8E12.4)
1854 123 format(i6,f6.2,7f8.1)
1855 122 FORMAT(1X,2I3,6F8.1,F8.3,F8.2)
1858 ! RETURN
1859 ! END
1860 !-------------------------------------------------------------------
1861 END SUBROUTINE SOIL
1862 !-------------------------------------------------------------------
1864 SUBROUTINE SICE ( &
1865 !--- input variables
1866 i,j,iland,isoil,delt,ktau,conflx,nzs,nddzs,nroot, &
1867 PRCPMS,RAINF,PATM,QVATM,QCATM,GLW,GSW, &
1868 EMISS,RNET,QKMS,TKMS,rho, &
1869 !--- sea ice parameters
1870 tice,rhosice,capice,thdifice, &
1871 zsmain,zshalf,DTDZS,DTDZS2,tbq, &
1872 !--- constants
1873 xlv,CP,rovcp,cw,stbolt,tabs, &
1874 !--- output variables
1875 tso,dew,soilt,qvg,qsg,qcg, &
1876 eeta,qfx,hfx,s,evapl,prcpl &
1877 )
1879 !*****************************************************************
1880 ! Energy budget and heat diffusion eqns. for
1881 ! sea ice
1882 !*************************************************************
1884 IMPLICIT NONE
1885 !-----------------------------------------------------------------
1887 !--- input variables
1889 INTEGER, INTENT(IN ) :: nroot,ktau,nzs , &
1890 nddzs !nddzs=2*(nzs-2)
1891 INTEGER, INTENT(IN ) :: i,j,iland,isoil
1892 REAL, INTENT(IN ) :: DELT,CONFLX
1893 !--- 3-D Atmospheric variables
1894 REAL, &
1895 INTENT(IN ) :: PATM, &
1896 QVATM, &
1897 QCATM
1898 !--- 2-D variables
1899 REAL, &
1900 INTENT(IN ) :: GLW, &
1901 GSW, &
1902 EMISS, &
1903 RHO, &
1904 QKMS, &
1905 TKMS
1906 !--- sea ice properties
1907 REAL, DIMENSION(1:NZS) , &
1908 INTENT(IN ) :: &
1909 tice, &
1910 rhosice, &
1911 capice, &
1912 thdifice
1915 REAL, INTENT(IN ) :: &
1916 CW, &
1917 XLV
1920 REAL, DIMENSION(1:NZS), INTENT(IN) :: ZSMAIN, &
1921 ZSHALF, &
1922 DTDZS2
1924 REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS
1926 REAL, DIMENSION(1:4001), INTENT(IN) :: TBQ
1929 !--- input/output variables
1930 !----soil temperature
1931 REAL, DIMENSION( 1:nzs ), INTENT(INOUT) :: TSO
1932 !-------- 2-d variables
1933 REAL, &
1934 INTENT(INOUT) :: DEW, &
1935 EETA, &
1936 EVAPL, &
1937 PRCPL, &
1938 QVG, &
1939 QSG, &
1940 QCG, &
1941 RNET, &
1942 QFX, &
1943 HFX, &
1944 S, &
1945 SOILT
1947 !--- Local variables
1948 REAL :: x,x1,x2,x4,tn,denom
1949 REAL :: RAINF, PRCPMS , &
1950 TABS, T3, UPFLUX, XINET
1952 REAL :: CP,rovcp,G0,LV,STBOLT,xlmelt,dzstop , &
1953 epot,fltot,ft,fq,hft,ras,cvw
1955 REAL :: FKT,D1,D2,D9,D10,DID,R211,R21,R22,R6,R7,D11 , &
1956 PI,H,FKQ,R210,AA,BB,PP,Q1,QS1,TS1,TQ2,TX2 , &
1957 TDENOM
1959 REAL :: AA1,RHCS
1962 REAL, DIMENSION(1:NZS) :: cotso,rhtso
1964 INTEGER :: nzs1,nzs2,k,k1,kn,kk
1966 !-----------------------------------------------------------------
1968 !-- define constants
1969 ! STBOLT=5.670151E-8
1970 XLMELT=3.35E+5
1971 cvw=cw
1973 prcpl=prcpms
1975 NZS1=NZS-1
1976 NZS2=NZS-2
1977 dzstop=1./(zsmain(2)-zsmain(1))
1978 RAS=RHO*1.E-3
1980 do k=1,nzs
1981 cotso(k)=0.
1982 rhtso(k)=0.
1983 enddo
1985 cotso(1)=0.
1986 rhtso(1)=TSO(NZS)
1988 DO 33 K=1,NZS2
1989 KN=NZS-K
1990 K1=2*KN-3
1991 X1=DTDZS(K1)*THDIFICE(KN-1)
1992 X2=DTDZS(K1+1)*THDIFICE(KN)
1993 FT=TSO(KN)+X1*(TSO(KN-1)-TSO(KN)) &
1994 -X2*(TSO(KN)-TSO(KN+1))
1995 DENOM=1.+X1+X2-X2*cotso(K)
1996 cotso(K+1)=X1/DENOM
1997 rhtso(K+1)=(FT+X2*rhtso(K))/DENOM
1998 33 CONTINUE
2000 !************************************************************************
2001 !--- THE HEAT BALANCE EQUATION (Smirnova et al., 1996, EQ. 21,26)
2002 RHCS=CAPICE(1)
2003 H=1.
2004 FKT=TKMS
2005 D1=cotso(NZS1)
2006 D2=rhtso(NZS1)
2007 TN=SOILT
2008 D9=THDIFICE(1)*RHCS*dzstop
2009 D10=TKMS*CP*RHO
2010 R211=.5*CONFLX/DELT
2011 R21=R211*CP*RHO
2012 R22=.5/(THDIFICE(1)*DELT*dzstop**2)
2013 R6=EMISS *STBOLT*.5*TN**4
2014 R7=R6/TN
2015 D11=RNET+R6
2016 TDENOM=D9*(1.-D1+R22)+D10+R21+R7 &
2017 +RAINF*CVW*PRCPMS
2018 FKQ=QKMS*RHO
2019 R210=R211*RHO
2020 AA=XLS*(FKQ+R210)/TDENOM
2021 BB=(D10*TABS+R21*TN+XLS*(QVATM*FKQ &
2022 +R210*QVG)+D11+D9*(D2+R22*TN) &
2023 +RAINF*CVW*PRCPMS*max(273.15,TABS) &
2024 )/TDENOM
2025 AA1=AA
2026 PP=PATM*1.E3
2027 AA1=AA1/PP
2028 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
2029 PRINT *,' VILKA-SEAICE1'
2030 print *,'D10,TABS,R21,TN,QVATM,FKQ', &
2031 D10,TABS,R21,TN,QVATM,FKQ
2032 print *,'RNET, EMISS, STBOLT, SOILT',RNET, EMISS, STBOLT, SOILT
2033 print *,'R210,QVG,D11,D9,D2,R22,RAINF,CVW,PRCPMS,TDENOM', &
2034 R210,QVG,D11,D9,D2,R22,RAINF,CVW,PRCPMS,TDENOM
2035 print *,'tn,aa1,bb,pp,fkq,r210', &
2036 tn,aa1,bb,pp,fkq,r210
2037 ENDIF
2038 CALL VILKA(TN,AA1,BB,PP,QS1,TS1,TBQ,KTAU,i,j,iland,isoil)
2039 !--- it is saturation over sea ice
2040 QVG=QS1
2041 QSG=QS1
2042 TSO(1)=min(271.4,TS1)
2043 QCG=0.
2044 !--- sea ice melting is not included in this simple approach
2045 !--- SOILT - skin temperature
2046 SOILT=TSO(1)
2047 !---- Final solution for soil temperature - TSO
2048 DO K=2,NZS
2049 KK=NZS-K+1
2050 TSO(K)=rhtso(KK)+cotso(KK)*TSO(K-1)
2051 END DO
2052 !--- CALCULATION OF DEW USING NEW VALUE OF QSG OR TRANSP IF NO DEW
2053 DEW=0.
2055 !--- THE DIAGNOSTICS OF SURFACE FLUXES
2056 T3 = STBOLT*SOILT*SOILT*SOILT
2057 UPFLUX = T3 *SOILT
2058 XINET = EMISS*(GLW-UPFLUX)
2059 RNET = GSW + XINET
2060 HFT=-TKMS*CP*RHO*(TABS-SOILT) &
2061 *(P1000mb*0.00001/Patm)**ROVCP
2062 Q1=-QKMS*RAS*(QVATM - QSG)
2063 IF (Q1.LE.0.) THEN
2064 ! --- condensation
2065 !-- moisture flux for coupling with PBL
2066 EETA=-QKMS*RAS*(QVATM/(1.+QVATM) - QSG/(1.+QSG))*1.E3
2067 QFX= XLS*EETA
2068 !-- actual moisture flux from RUC LSM
2069 DEW=QKMS*(QVATM-QSG)
2070 EETA= RHO*DEW
2071 ELSE
2072 ! --- evaporation
2073 !-- moisture flux for coupling with PBL
2074 EETA=-QKMS*RAS*(QVATM/(1.+QVATM) - QVG/(1.+QVG))*1.E3
2075 ! EETA=-QKMS*RAS*(QVATM/(1.+QVATM) - QSG/(1.+QSG))*1.E3
2076 ! to convert from kg m-2 s-1 to m s-1: 1/rho water=1.e-3************
2077 QFX= XLS * EETA
2078 !-- actual moisture flux from RUC LSM
2079 EETA = Q1*1.E3
2080 ENDIF
2081 EVAPL=EETA
2082 S=THDIFICE(1)*CAPICE(1)*DZSTOP*(TSO(1)-TSO(2))
2083 HFX=HFT
2084 FLTOT=RNET-HFT-QFX-S
2086 !-------------------------------------------------------------------
2087 END SUBROUTINE SICE
2088 !-------------------------------------------------------------------
2092 SUBROUTINE SNOWSOIL ( &
2093 !--- input variables
2094 i,j,isoil,delt,ktau,conflx,nzs,nddzs,nroot, &
2095 meltfactor,rhonewsn, & ! new
2096 ILAND,PRCPMS,RAINF,NEWSNOW,snhei,SNWE,SNOWFRAC, &
2097 RHOSN, &
2098 PATM,QVATM,QCATM, &
2099 GLW,GSW,EMISS,RNET,IVGTYP, &
2100 QKMS,TKMS,PC,cst,rho,vegfrac,alb,znt, &
2101 MYJ, &
2102 !--- soil fixed fields
2103 QWRTZ,rhocs,dqm,qmin,ref,wilt,psis,bclh,ksat, &
2104 sat,cn,zsmain,zshalf,DTDZS,DTDZS2,tbq, &
2105 !--- constants
2106 xlv,CP,rovcp,G0_P,cw,stbolt,TABS, &
2107 KQWRTZ,KICE,KWT, &
2108 !--- output variables
2109 ilnb,snweprint,snheiprint,rsm, &
2110 soilmois,tso,smfrkeep,keepfr, &
2111 dew,soilt,soilt1,tsnav, &
2112 qvg,qsg,qcg,SMELT,SNOH,SNFLX,SNOM, &
2113 edir1,ec1,ett1,eeta,qfx,hfx,s,sublim, &
2114 prcpl,runoff1,runoff2,mavail,soilice, &
2115 soiliqw,infiltrp )
2117 !***************************************************************
2118 ! Energy and moisture budget for snow, heat diffusion eqns.
2119 ! in snow and soil, Richards eqn. for soil covered with snow
2120 !
2121 ! DELT - time step (s)
2122 ! ktau - numver of time step
2123 ! CONFLX - depth of constant flux layer (m)
2124 ! J,I - the location of grid point
2125 ! IME, JME, NZS - dimensions of the domain
2126 ! NROOT - number of levels within the root zone
2127 ! PRCPMS - precipitation rate in m/s
2128 ! NEWSNOW - pcpn in soilid form (m)
2129 ! SNHEI, SNWE - snow height and snow water equivalent (m)
2130 ! RHOSN - snow density (kg/m-3)
2131 ! PATM - pressure (bar)
2132 ! QVATM,QCATM - cloud and water vapor mixing ratio
2133 ! at the first atm. level (kg/kg)
2134 ! GLW, GSW - incoming longwave and absorbed shortwave
2135 ! radiation at the surface (W/m^2)
2136 ! EMISS,RNET - emissivity (0-1) of the ground surface and net
2137 ! radiation at the surface (W/m^2)
2138 ! QKMS - exchange coefficient for water vapor in the
2139 ! surface layer (m/s)
2140 ! TKMS - exchange coefficient for heat in the surface
2141 ! layer (m/s)
2142 ! PC - plant coefficient (resistance) (0-1)
2143 ! RHO - density of atmosphere near surface (kg/m^3)
2144 ! VEGFRAC - greeness fraction (0-1)
2145 ! RHOCS - volumetric heat capacity of dry soil (J/m^3/K)
2146 ! DQM, QMIN - porosity minus residual soil moisture QMIN (m^3/m^3)
2147 ! REF, WILT - field capacity soil moisture and the
2148 ! wilting point (m^3/m^3)
2149 ! PSIS - matrix potential at saturation (m)
2150 ! BCLH - exponent for Clapp-Hornberger parameterization
2151 ! KSAT - saturated hydraulic conductivity (m/s)
2152 ! SAT - maximum value of water intercepted by canopy (m)
2153 ! CN - exponent for calculation of canopy water
2154 ! ZSMAIN - main levels in soil (m)
2155 ! ZSHALF - middle of the soil layers (m)
2156 ! DTDZS,DTDZS2 - dt/(2.*dzshalf*dzmain) and dt/dzshalf in soil
2157 ! TBQ - table to define saturated mixing ration
2158 ! of water vapor for given temperature and pressure
2159 ! ilnb - number of layers in snow
2160 ! rsm - liquid water inside snow pack (m)
2161 ! SOILMOIS,TSO - soil moisture (m^3/m^3) and temperature (K)
2162 ! DEW - dew in (kg/m^2 s)
2163 ! SOILT - skin temperature (K)
2164 ! SOILT1 - snow temperature at 7.5 cm depth (K)
2165 ! TSNAV - average temperature of snow pack (C)
2166 ! QSG,QVG,QCG - saturated mixing ratio, mixing ratio of
2167 ! water vapor and cloud at the ground
2168 ! surface, respectively (kg/kg)
2169 ! EDIR1, EC1, ETT1, EETA - direct evaporation, evaporation of
2170 ! canopy water, transpiration (kg m-2 s-1) and total
2171 ! evaporation in (m s-1).
2172 ! QFX, HFX - latent and sensible heat fluxes (W/m^2)
2173 ! S - soil heat flux in the top layer (W/m^2)
2174 ! SUBLIM - snow sublimation (kg/m^2/s)
2175 ! RUNOFF1 - surface runoff (m/s)
2176 ! RUNOFF2 - underground runoff (m)
2177 ! MAVAIL - moisture availability in the top soil layer (0-1)
2178 ! SOILICE - content of soil ice in soil layers (m^3/m^3)
2179 ! SOILIQW - lliquid water in soil layers (m^3/m^3)
2180 ! INFILTRP - infiltration flux from the top of soil domain (m/s)
2181 ! XINET - net long-wave radiation (W/m^2)
2182 !
2183 !*******************************************************************
2185 IMPLICIT NONE
2186 !-------------------------------------------------------------------
2187 !--- input variables
2189 INTEGER, INTENT(IN ) :: nroot,ktau,nzs , &
2190 nddzs !nddzs=2*(nzs-2)
2191 INTEGER, INTENT(IN ) :: i,j,isoil
2193 REAL, INTENT(IN ) :: DELT,CONFLX,PRCPMS , &
2194 RAINF,NEWSNOW,RHONEWSN,meltfactor
2196 LOGICAL, INTENT(IN ) :: myj
2198 !--- 3-D Atmospheric variables
2199 REAL, &
2200 INTENT(IN ) :: PATM, &
2201 QVATM, &
2202 QCATM
2203 !--- 2-D variables
2204 REAL , &
2205 INTENT(IN ) :: GLW, &
2206 GSW, &
2207 RHO, &
2208 PC, &
2209 VEGFRAC, &
2210 QKMS, &
2211 TKMS
2213 INTEGER, INTENT(IN ) :: IVGTYP
2214 !--- soil properties
2215 REAL , &
2216 INTENT(IN ) :: RHOCS, &
2217 BCLH, &
2218 DQM, &
2219 KSAT, &
2220 PSIS, &
2221 QMIN, &
2222 QWRTZ, &
2223 REF, &
2224 SAT, &
2225 WILT
2227 REAL, INTENT(IN ) :: CN, &
2228 CW, &
2229 XLV, &
2230 G0_P, &
2231 KQWRTZ, &
2232 KICE, &
2233 KWT
2236 REAL, DIMENSION(1:NZS), INTENT(IN) :: ZSMAIN, &
2237 ZSHALF, &
2238 DTDZS2
2240 REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS
2242 REAL, DIMENSION(1:4001), INTENT(IN) :: TBQ
2245 !--- input/output variables
2246 !-------- 3-d soil moisture and temperature
2247 REAL, DIMENSION( 1:nzs ) , &
2248 INTENT(INOUT) :: TSO, &
2249 SOILMOIS, &
2250 SMFRKEEP
2252 REAL, DIMENSION( 1:nzs ) , &
2253 INTENT(INOUT) :: KEEPFR
2256 INTEGER, INTENT(INOUT) :: ILAND
2259 !-------- 2-d variables
2260 REAL , &
2261 INTENT(INOUT) :: DEW, &
2262 CST, &
2263 EDIR1, &
2264 EC1, &
2265 ETT1, &
2266 EETA, &
2267 RHOSN, &
2268 SUBLIM, &
2269 PRCPL, &
2270 ALB, &
2271 EMISS, &
2272 ZNT, &
2273 MAVAIL, &
2274 QVG, &
2275 QSG, &
2276 QCG, &
2277 QFX, &
2278 HFX, &
2279 S, &
2280 RUNOFF1, &
2281 RUNOFF2, &
2282 SNWE, &
2283 SNHEI, &
2284 SMELT, &
2285 SNOM, &
2286 SNOH, &
2287 SNFLX, &
2288 SOILT, &
2289 SOILT1, &
2290 SNOWFRAC, &
2291 TSNAV
2293 INTEGER, INTENT(INOUT) :: ILNB
2295 !-------- 1-d variables
2296 REAL, DIMENSION(1:NZS), INTENT(OUT) :: SOILICE, &
2297 SOILIQW
2299 REAL, INTENT(OUT) :: RSM, &
2300 SNWEPRINT, &
2301 SNHEIPRINT
2302 !--- Local variables
2305 INTEGER :: nzs1,nzs2,k
2307 REAL :: INFILTRP, TRANSUM , &
2308 SNTH, NEWSN , &
2309 TABS, T3, UPFLUX, XINET , &
2310 BETA, SNWEPR,EPDT,PP
2311 REAL :: CP,rovcp,G0,LV,xlvm,STBOLT,xlmelt,dzstop , &
2312 can,epot,fac,fltot,ft,fq,hft , &
2313 q1,ras,rhoice,sph , &
2314 trans,zn,ci,cvw,tln,tavln,pi , &
2315 DD1,CMC2MS,DRYCAN,WETCAN , &
2316 INFMAX,RIW,DELTSN,H,UMVEG
2318 REAL, DIMENSION(1:NZS) :: transp,cap,diffu,hydro , &
2319 thdif,tranf,tav,soilmoism , &
2320 soilicem,soiliqwm,detal , &
2321 fwsat,lwsat,told,smold
2322 REAL :: drip
2324 REAL :: RNET
2326 !-----------------------------------------------------------------
2328 cvw=cw
2329 XLMELT=3.35E+5
2330 !-- the next line calculates heat of sublimation of water vapor
2331 XLVm=XLV+XLMELT
2332 ! STBOLT=5.670151E-8
2334 !--- SNOW flag -- 99
2335 ILAND=99
2337 !--- DELTSN - is the threshold for splitting the snow layer into 2 layers.
2338 !--- With snow density 400 kg/m^3, this threshold is equal to 7.5 cm,
2339 !--- equivalent to 0.03 m SNWE. For other snow densities the threshold is
2340 !--- computed using SNWE=0.03 m and current snow density.
2341 !--- SNTH - the threshold below which the snow layer is combined with
2342 !--- the top soil layer. SNTH is computed using snwe=0.016 m, and
2343 !--- equals 4 cm for snow density 400 kg/m^3.
2345 DELTSN=0.0301*1.e3/rhosn
2346 snth=0.01601*1.e3/rhosn
2348 !tgs when the snow depth is marginlly higher than DELTSN,
2349 ! reset DELTSN to half of snow depth.
2350 IF(SNHEI.GT.DELTSN+SNTH) THEN
2351 ! 2-layer model
2352 if(snhei-deltsn-snth.lt.snth) deltsn=0.5*(snhei-snth)
2353 ENDIF
2355 RHOICE=900.
2356 CI=RHOICE*2100.
2357 RAS=RHO*1.E-3
2358 RIW=rhoice*1.e-3
2359 MAVAIL=1.
2360 RSM=0.
2362 DO K=1,NZS
2363 TRANSP (K)=0.
2364 soilmoism (k)=0.
2365 soiliqwm (k)=0.
2366 soilice (k)=0.
2367 soilicem (k)=0.
2368 lwsat (k)=0.
2369 fwsat (k)=0.
2370 tav (k)=0.
2371 cap (k)=0.
2372 diffu (k)=0.
2373 hydro (k)=0.
2374 thdif (k)=0.
2375 tranf (k)=0.
2376 detal (k)=0.
2377 told (k)=0.
2378 smold (k)=0.
2379 ENDDO
2381 snweprint=0.
2382 snheiprint=0.
2383 prcpl=prcpms
2385 !*** DELTSN is the depth of the top layer of snow where
2386 !*** there is a temperature gradient, the rest of the snow layer
2387 !*** is considered to have constant temperature
2390 NZS1=NZS-1
2391 NZS2=NZS-2
2392 DZSTOP=1./(zsmain(2)-zsmain(1))
2394 !----- THE CALCULATION OF THERMAL DIFFUSIVITY, DIFFUSIONAL AND ---
2395 !----- HYDRAULIC CONDUCTIVITY (SMIRNOVA ET AL. 1996, EQ.2,5,6) ---
2396 !tgs - the following loop is added to define the amount of frozen
2397 !tgs - water in soil if there is any
2398 DO K=1,NZS
2400 tln=log(tso(k)/273.15)
2401 if(tln.lt.0.) then
2402 soiliqw(k)=(dqm+qmin)*(XLMELT* &
2403 (tso(k)-273.15)/tso(k)/9.81/psis) &
2404 **(-1./bclh)-qmin
2405 soiliqw(k)=max(0.,soiliqw(k))
2406 soiliqw(k)=min(soiliqw(k),soilmois(k))
2407 soilice(k)=(soilmois(k)-soiliqw(k))/riw
2409 !---- melting and freezing is balanced, soil ice cannot increase
2410 if(keepfr(k).eq.1.) then
2411 soilice(k)=min(soilice(k),smfrkeep(k))
2412 soiliqw(k)=max(0.,soilmois(k)-soilice(k)*rhoice*1.e-3)
2413 endif
2415 else
2416 soilice(k)=0.
2417 soiliqw(k)=soilmois(k)
2418 endif
2420 ENDDO
2422 DO K=1,NZS1
2424 tav(k)=0.5*(tso(k)+tso(k+1))
2425 soilmoism(k)=0.5*(soilmois(k)+soilmois(k+1))
2426 tavln=log(tav(k)/273.15)
2428 if(tavln.lt.0.) then
2429 soiliqwm(k)=(dqm+qmin)*(XLMELT* &
2430 (tav(k)-273.15)/tav(k)/9.81/psis) &
2431 **(-1./bclh)-qmin
2432 fwsat(k)=dqm-soiliqwm(k)
2433 lwsat(k)=soiliqwm(k)+qmin
2434 soiliqwm(k)=max(0.,soiliqwm(k))
2435 soiliqwm(k)=min(soiliqwm(k), soilmoism(k))
2436 soilicem(k)=(soilmoism(k)-soiliqwm(k))/riw
2437 !---- melting and freezing is balanced, soil ice cannot increase
2438 if(keepfr(k).eq.1.) then
2439 soilicem(k)=min(soilicem(k), &
2440 0.5*(smfrkeep(k)+smfrkeep(k+1)))
2441 soiliqwm(k)=max(0.,soilmoism(k)-soilicem(k)*riw)
2442 fwsat(k)=dqm-soiliqwm(k)
2443 lwsat(k)=soiliqwm(k)+qmin
2444 endif
2446 else
2447 soilicem(k)=0.
2448 soiliqwm(k)=soilmoism(k)
2449 lwsat(k)=dqm+qmin
2450 fwsat(k)=0.
2452 endif
2453 ENDDO
2455 do k=1,nzs
2456 if(soilice(k).gt.0.) then
2457 smfrkeep(k)=soilice(k)
2458 else
2459 smfrkeep(k)=soilmois(k)/riw
2460 endif
2461 enddo
2463 !******************************************************************
2464 ! SOILPROP computes thermal diffusivity, and diffusional and
2465 ! hydraulic condeuctivities
2466 !******************************************************************
2467 CALL SOILPROP( &
2468 !--- input variables
2469 nzs,fwsat,lwsat,tav,keepfr, &
2470 soilmois,soiliqw,soilice, &
2471 soilmoism,soiliqwm,soilicem, &
2472 !--- soil fixed fields
2473 QWRTZ,rhocs,dqm,qmin,psis,bclh,ksat, &
2474 !--- constants
2475 riw,xlmelt,CP,G0_P,cvw,ci, &
2476 kqwrtz,kice,kwt, &
2477 !--- output variables
2478 thdif,diffu,hydro,cap)
2480 !********************************************************************
2481 !--- CALCULATION OF CANOPY WATER (Smirnova et al., 1996, EQ.16) AND DEW
2483 DRIP=0.
2484 SMELT=0.
2485 DD1=0.
2486 H=1.
2488 FQ=QKMS
2491 !--- If vegfrac.ne.0. then part of falling snow can be
2492 !--- intercepted by the canopy.
2494 DEW=0.
2495 UMVEG=1.-vegfrac
2496 EPOT = -FQ*(QVATM-QSG)
2498 IF(vegfrac.EQ.0.) then
2499 cst=0.
2500 drip=0.
2501 ELSE
2502 IF(EPOT.GE.0.) THEN
2503 ! Evaporation
2504 ! DD1=CST+(NEWSNOW*RHOSN*1.E-3 &
2505 DD1=CST+(NEWSNOW*RHOnewSN*1.E-3 &
2506 !-- this change will not let liquid waer be intercepted by the canopy....
2507 -DELT*(RAS*EPOT &
2508 ! -DELT*(-PRCPMS+RAS*EPOT &
2509 *(CST/SAT)**CN)) *vegfrac
2510 ELSE
2511 ! Sublimation
2512 DEW = - EPOT
2513 ! DD1=CST+(NEWSNOW*RHOSN*1.E-3+delt*( &
2514 DD1=CST+(NEWSNOW*RHOnewSN*1.E-3+delt*( &
2515 ! DD1=CST+(NEWSNOW*RHOSN*1.E-3+delt*(PRCPMS &
2516 +DEW*RAS)) *vegfrac
2517 ENDIF
2519 IF(DD1.LT.0.) DD1=0.
2520 IF (vegfrac.GT.0.) THEN
2521 CST=DD1
2522 IF(CST.GT.SAT*vegfrac) THEN
2523 CST=SAT*vegfrac
2524 DRIP=DD1-SAT*vegfrac
2525 ENDIF
2526 ENDIF
2529 !--- In SFCTMP NEWSNOW is added to SNHEI as if there is no vegetation
2530 !--- With vegetation part of NEWSNOW can be intercepted by canopy until
2531 !--- the saturation is reached. After the canopy saturation is reached
2532 !--- DRIP in the solid form will be added to SNOW cover.
2534 SNWE=(SNHEI-vegfrac*NEWSNOW)*RHOSN*1.E-3 &
2535 + DRIP
2537 ENDIF
2539 DRIP=0.
2540 SNHEI=SNWE*1.e3/RHOSN
2541 SNWEPR=SNWE
2543 ! check if all snow can evaporate during DT
2544 BETA=1.
2545 EPDT = EPOT * RAS *DELT*UMVEG
2546 IF(SNWEPR.LE.EPDT) THEN
2547 BETA=SNWEPR/max(1.e-8,EPDT)
2548 SNWE=0.
2549 SNHEI=0.
2550 ENDIF
2552 WETCAN=(CST/SAT)**CN
2553 DRYCAN=1.-WETCAN
2555 !**************************************************************
2556 ! TRANSF computes transpiration function
2557 !**************************************************************
2558 CALL TRANSF( &
2559 !--- input variables
2560 nzs,nroot,soiliqw,tabs, &
2561 !--- soil fixed fields
2562 dqm,qmin,ref,wilt,zshalf, &
2563 !--- output variables
2564 tranf,transum)
2566 !--- Save soil temp and moisture from the beginning of time step
2567 do k=1,nzs
2568 told(k)=tso(k)
2569 smold(k)=soilmois(k)
2570 enddo
2572 !**************************************************************
2573 ! SNOWTEMP solves heat budget and diffusion eqn. in soil
2574 !**************************************************************
2576 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
2577 print *, 'TSO before calling SNOWTEMP: ', tso
2578 ENDIF
2579 CALL SNOWTEMP( &
2580 !--- input variables
2581 i,j,iland,isoil, &
2582 delt,ktau,conflx,nzs,nddzs,nroot, &
2583 snwe,snwepr,snhei,newsnow,snowfrac, &
2584 beta,deltsn,snth,rhosn,rhonewsn,meltfactor, & ! add meltfactor
2585 PRCPMS,RAINF, &
2586 PATM,TABS,QVATM,QCATM, &
2587 GLW,GSW,EMISS,RNET, &
2588 QKMS,TKMS,PC,rho,vegfrac, &
2589 thdif,cap,drycan,wetcan,cst, &
2590 tranf,transum,dew,mavail, &
2591 !--- soil fixed fields
2592 dqm,qmin,psis,bclh, &
2593 zsmain,zshalf,DTDZS,tbq, &
2594 !--- constants
2595 xlvm,CP,rovcp,G0_P,cvw,stbolt, &
2596 !--- output variables
2597 snweprint,snheiprint,rsm, &
2598 tso,soilt,soilt1,tsnav,qvg,qsg,qcg, &
2599 smelt,snoh,snflx,ilnb)
2601 !************************************************************************
2602 !--- RECALCULATION OF DEW USING NEW VALUE OF QSG OR TRANSP IF NO DEW
2603 DEW=0.
2604 ETT1=0.
2605 PP=PATM*1.E3
2606 QSG= QSN(SOILT,TBQ)/PP
2607 EPOT = -FQ*(QVATM-QSG)
2608 IF(EPOT.GE.0.) THEN
2609 ! Evaporation
2610 DO K=1,NROOT
2611 TRANSP(K)=vegfrac*RAS*FQ*(QVATM-QSG) &
2612 *PC*tranf(K)*DRYCAN/zshalf(NROOT+1)
2613 IF(TRANSP(K).GT.0.) TRANSP(K)=0.
2614 ETT1=ETT1-TRANSP(K)
2615 ENDDO
2616 DO k=nroot+1,nzs
2617 transp(k)=0.
2618 enddo
2620 ELSE
2621 ! Sublimation
2622 DEW=-EPOT
2623 DO K=1,NZS
2624 TRANSP(K)=0.
2625 ENDDO
2626 ETT1=0.
2627 ENDIF
2629 !-- recalculating of frozen water in soil
2630 DO K=1,NZS
2631 tln=log(tso(k)/273.15)
2632 if(tln.lt.0.) then
2633 soiliqw(k)=(dqm+qmin)*(XLMELT* &
2634 (tso(k)-273.15)/tso(k)/9.81/psis) &
2635 **(-1./bclh)-qmin
2636 soiliqw(k)=max(0.,soiliqw(k))
2637 soiliqw(k)=min(soiliqw(k),soilmois(k))
2638 soilice(k)=(soilmois(k)-soiliqw(k))/riw
2639 !---- melting and freezing is balanced, soil ice cannot increase
2640 if(keepfr(k).eq.1.) then
2641 soilice(k)=min(soilice(k),smfrkeep(k))
2642 soiliqw(k)=max(0.,soilmois(k)-soilice(k)*riw)
2643 endif
2645 else
2646 soilice(k)=0.
2647 soiliqw(k)=soilmois(k)
2648 endif
2649 ENDDO
2651 !*************************************************************************
2652 !--- TQCAN FOR SOLUTION OF MOISTURE BALANCE (Smirnova et al. 1996, EQ.22,28)
2653 ! AND TSO,ETA PROFILES
2654 !*************************************************************************
2655 CALL SOILMOIST ( &
2656 !-- input
2657 delt,nzs,nddzs,DTDZS,DTDZS2,RIW, &
2658 zsmain,zshalf,diffu,hydro, &
2659 QSG,QVG,QCG,QCATM,QVATM,-PRCPMS, &
2660 0.,TRANSP,0., &
2661 0.,SMELT,soilice,vegfrac, &
2662 !-- soil properties
2663 DQM,QMIN,REF,KSAT,RAS,INFMAX, &
2664 !-- output
2665 SOILMOIS,SOILIQW,MAVAIL,RUNOFF1, &
2666 RUNOFF2,infiltrp)
2668 ! 4 Nov 07 - update CST for snow melt
2669 if(snwe.ne.0.) then
2670 CST=(1.-min(1.,smelt/snwe))*CST
2671 else
2672 CST=0.
2673 endif
2675 !-- Restore land-use parameters if snow is less than threshold
2676 IF(SNHEI.EQ.0.) then
2677 tsnav=soilt-273.15
2678 CALL SNOWFREE(ivgtyp,myj,emiss, &
2679 znt,iland)
2680 smelt=smelt+snwe/delt
2681 rsm=0.
2682 ! snwe=0.
2683 ENDIF
2685 ! 21apr2009
2686 ! SNOM goes into the passed-in ACSNOM variable in the grid derived type
2687 SNOM=SNOM+SMELT*DELT*1.e3
2689 !--- KEEPFR is 1 when the temperature and moisture in soil
2690 !--- are both increasing. In this case soil ice should not
2691 !--- be increasing according to the freezing curve.
2692 !--- Some part of ice is melted, but additional water is
2693 !--- getting frozen. Thus, only structure of frozen soil is
2694 !--- changed, and phase changes are not affecting the heat
2695 !--- transfer. This situation may happen when it rains on the
2696 !--- frozen soil.
2698 do k=1,nzs
2699 if (soilice(k).gt.0.) then
2700 if(tso(k).gt.told(k).and.soilmois(k).gt.smold(k)) then
2701 keepfr(k)=1.
2702 else
2703 keepfr(k)=0.
2704 endif
2705 endif
2706 enddo
2707 !--- THE DIAGNOSTICS OF SURFACE FLUXES
2709 T3 = STBOLT*SOILT*SOILT*SOILT
2710 UPFLUX = T3 *SOILT
2711 XINET = EMISS*(GLW-UPFLUX)
2712 RNET = GSW + XINET
2713 HFT=-TKMS*CP*RHO*(TABS-SOILT) &
2714 *(P1000mb*0.00001/Patm)**ROVCP
2715 Q1 = - FQ*RAS* (QVATM - QSG)
2717 IF (Q1.LT.0.) THEN
2718 ! --- condensation
2719 EDIR1=0.
2720 EC1=0.
2721 ETT1=0.
2722 ! --- condensation
2723 !-- moisture flux for coupling with PBL
2724 EETA=-QKMS*RAS*(QVATM/(1.+QVATM) - QSG/(1.+QSG))*1.E3
2725 QFX= XLVm*EETA
2726 !-- actual moisture flux from RUC LSM
2727 DEW=QKMS*(QVATM-QSG)
2728 EETA= RHO*DEW
2729 ELSE
2730 ! --- evaporation
2731 EDIR1 = Q1*UMVEG *BETA
2732 EC1 = Q1 * WETCAN
2733 CMC2MS=CST/DELT
2734 if(EC1.gt.CMC2MS) cst=0.
2735 EC1=MIN(CMC2MS,EC1)*vegfrac
2736 !-- moisture flux for coupling with PBL
2737 EETA=-QKMS*RAS*(QVATM/(1.+QVATM) - QVG/(1.+QVG))*1.E3
2738 ! EETA=-QKMS*RAS*(QVATM/(1.+QVATM) - QSG/(1.+QSG))*1.E3
2739 ! to convert from kg m-2 s-1 to m s-1: 1/rho water=1.e-3************
2740 QFX= XLVm * EETA
2741 !-- actual moisture flux from RUC LSM
2742 EETA = (EDIR1 + EC1 + ETT1)*1.E3
2743 ENDIF
2744 s=THDIF(1)*CAP(1)*dzstop*(tso(1)-tso(2))
2745 HFX=HFT
2746 FLTOT=RNET-HFT-QFX-S
2748 222 CONTINUE
2750 1123 FORMAT(I5,8F12.3)
2751 1133 FORMAT(I7,8E12.4)
2752 123 format(i6,f6.2,7f8.1)
2753 122 FORMAT(1X,2I3,6F8.1,F8.3,F8.2)
2756 ! RETURN
2757 ! END
2758 !-------------------------------------------------------------------
2759 END SUBROUTINE SNOWSOIL
2760 !-------------------------------------------------------------------
2762 SUBROUTINE SNOWSEAICE( &
2763 i,j,isoil,delt,ktau,conflx,nzs,nddzs, &
2764 meltfactor,rhonewsn, & ! new
2765 ILAND,PRCPMS,RAINF,NEWSNOW,snhei,SNWE,snowfrac, &
2766 RHOSN,PATM,QVATM,QCATM, &
2767 GLW,GSW,EMISS,RNET, &
2768 QKMS,TKMS,RHO, &
2769 !--- sea ice parameters
2770 ALB,ZNT, &
2771 tice,rhosice,capice,thdifice, &
2772 zsmain,zshalf,DTDZS,DTDZS2,tbq, &
2773 !--- constants
2774 xlv,CP,rovcp,cw,stbolt,tabs, &
2775 !--- output variables
2776 ilnb,snweprint,snheiprint,rsm,tso, &
2777 dew,soilt,soilt1,tsnav,qvg,qsg,qcg, &
2778 SMELT,SNOH,SNFLX,SNOM,eeta, &
2779 qfx,hfx,s,sublim,prcpl &
2780 )
2781 !***************************************************************
2782 ! Solving energy budget for snow on sea ice and heat diffusion
2783 ! eqns. in snow and sea ice
2784 !***************************************************************
2787 IMPLICIT NONE
2788 !-------------------------------------------------------------------
2789 !--- input variables
2791 INTEGER, INTENT(IN ) :: ktau,nzs , &
2792 nddzs !nddzs=2*(nzs-2)
2793 INTEGER, INTENT(IN ) :: i,j,isoil
2795 REAL, INTENT(IN ) :: DELT,CONFLX,PRCPMS , &
2796 RAINF,NEWSNOW,RHONEWSN,meltfactor
2797 real :: rhonewcsn
2799 !--- 3-D Atmospheric variables
2800 REAL, &
2801 INTENT(IN ) :: PATM, &
2802 QVATM, &
2803 QCATM
2804 !--- 2-D variables
2805 REAL , &
2806 INTENT(IN ) :: GLW, &
2807 GSW, &
2808 RHO, &
2809 QKMS, &
2810 TKMS
2812 !--- sea ice properties
2813 REAL, DIMENSION(1:NZS) , &
2814 INTENT(IN ) :: &
2815 tice, &
2816 rhosice, &
2817 capice, &
2818 thdifice
2820 REAL, INTENT(IN ) :: &
2821 CW, &
2822 XLV
2824 REAL, DIMENSION(1:NZS), INTENT(IN) :: ZSMAIN, &
2825 ZSHALF, &
2826 DTDZS2
2828 REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS
2830 REAL, DIMENSION(1:4001), INTENT(IN) :: TBQ
2832 !--- input/output variables
2833 !-------- 3-d soil moisture and temperature
2834 REAL, DIMENSION( 1:nzs ) , &
2835 INTENT(INOUT) :: TSO
2837 INTEGER, INTENT(INOUT) :: ILAND
2840 !-------- 2-d variables
2841 REAL , &
2842 INTENT(INOUT) :: DEW, &
2843 EETA, &
2844 RHOSN, &
2845 SUBLIM, &
2846 PRCPL, &
2847 ALB, &
2848 EMISS, &
2849 ZNT, &
2850 QVG, &
2851 QSG, &
2852 QCG, &
2853 QFX, &
2854 HFX, &
2855 S, &
2856 SNWE, &
2857 SNHEI, &
2858 SMELT, &
2859 SNOM, &
2860 SNOH, &
2861 SNFLX, &
2862 SOILT, &
2863 SOILT1, &
2864 SNOWFRAC, &
2865 TSNAV
2867 INTEGER, INTENT(INOUT) :: ILNB
2869 REAL, INTENT(OUT) :: RSM, &
2870 SNWEPRINT, &
2871 SNHEIPRINT
2872 !--- Local variables
2875 INTEGER :: nzs1,nzs2,k,k1,kn,kk
2876 REAL :: x,x1,x2,dzstop,ft,tn,denom
2878 REAL :: SNTH, NEWSN , &
2879 TABS, T3, UPFLUX, XINET , &
2880 BETA, SNWEPR,EPDT,PP
2881 REAL :: CP,rovcp,G0,LV,xlvm,STBOLT,xlmelt , &
2882 epot,fltot,fq,hft,q1,ras,rhoice,ci,cvw , &
2883 RIW,DELTSN,H
2885 REAL :: rhocsn,thdifsn, &
2886 xsn,ddzsn,x1sn,d1sn,d2sn,d9sn,r22sn
2888 REAL :: cotsn,rhtsn,xsn1,ddzsn1,x1sn1,ftsnow,denomsn
2889 REAL :: fso,fsn, &
2890 FKT,D1,D2,D9,D10,DID,R211,R21,R22,R6,R7,D11, &
2891 FKQ,R210,AA,BB,QS1,TS1,TQ2,TX2, &
2892 TDENOM,AA1,RHCS,H1,TSOB, SNPRIM, &
2893 SNODIF,SOH,TNOLD,QGOLD,SNOHGNEW
2894 REAL, DIMENSION(1:NZS) :: cotso,rhtso
2896 REAL :: RNET,rsmfrac,soiltfrac,hsn
2897 integer :: nmelt
2900 !-----------------------------------------------------------------
2901 XLMELT=3.35E+5
2902 !-- the next line calculates heat of sublimation of water vapor
2903 XLVm=XLV+XLMELT
2904 ! STBOLT=5.670151E-8
2906 !--- SNOW flag -- 99
2907 ILAND=99
2909 !--- DELTSN - is the threshold for splitting the snow layer into 2 layers.
2910 !--- With snow density 400 kg/m^3, this threshold is equal to 7.5 cm,
2911 !--- equivalent to 0.03 m SNWE. For other snow densities the threshold is
2912 !--- computed using SNWE=0.03 m and current snow density.
2913 !--- SNTH - the threshold below which the snow layer is combined with
2914 !--- the top sea ice layer. SNTH is computed using snwe=0.016 m, and
2915 !--- equals 4 cm for snow density 400 kg/m^3.
2917 DELTSN=0.0301*1.e3/rhosn
2918 snth=0.01601*1.e3/rhosn
2920 !tgs when the snow depth is marginlly higher than DELTSN,
2921 ! reset DELTSN to half of snow depth.
2922 IF(SNHEI.GT.DELTSN+SNTH) THEN
2923 ! 2-layer model
2924 if(snhei-deltsn-snth.lt.snth) deltsn=0.5*(snhei-snth)
2925 ENDIF
2928 RHOICE=900.
2929 CI=RHOICE*2100.
2930 RAS=RHO*1.E-3
2931 RIW=rhoice*1.e-3
2932 RSM=0.
2934 XLMELT=3.35E+5
2935 RHOCSN=2090.* RHOSN
2936 !18apr08 - add rhonewcsn
2937 RHOnewCSN=2090.* RHOnewSN
2938 THDIFSN = 0.265/RHOCSN
2939 RAS=RHO*1.E-3
2941 SOILTFRAC=SOILT
2943 SMELT=0.
2944 SOH=0.
2945 SNODIF=0.
2946 SNOH=0.
2947 SNOHGNEW=0.
2948 RSM = 0.
2949 RSMFRAC = 0.
2950 fsn=1.
2951 fso=0.
2952 hsn=snhei
2954 NZS1=NZS-1
2955 NZS2=NZS-2
2957 QGOLD=QVG
2958 TNOLD=SOILT
2959 DZSTOP=1./(ZSMAIN(2)-ZSMAIN(1))
2961 snweprint=0.
2962 snheiprint=0.
2963 prcpl=prcpms
2965 !*** DELTSN is the depth of the top layer of snow where
2966 !*** there is a temperature gradient, the rest of the snow layer
2967 !*** is considered to have constant temperature
2970 H=1.
2971 SMELT=0.
2973 FQ=QKMS
2974 SNHEI=SNWE*1.e3/RHOSN
2975 SNWEPR=SNWE
2977 ! check if all snow can evaporate during DT
2978 BETA=1.
2979 EPDT = EPOT * RAS *DELT
2980 IF(SNWEPR.LE.EPDT) THEN
2981 BETA=SNWEPR/max(1.e-8,EPDT)
2982 SNWE=0.
2983 SNHEI=0.
2984 ENDIF
2986 !******************************************************************************
2987 ! COEFFICIENTS FOR THOMAS ALGORITHM FOR TSO
2988 !******************************************************************************
2990 cotso(1)=0.
2991 rhtso(1)=TSO(NZS)
2992 DO 33 K=1,NZS2
2993 KN=NZS-K
2994 K1=2*KN-3
2995 X1=DTDZS(K1)*THDIFICE(KN-1)
2996 X2=DTDZS(K1+1)*THDIFICE(KN)
2997 FT=TSO(KN)+X1*(TSO(KN-1)-TSO(KN)) &
2998 -X2*(TSO(KN)-TSO(KN+1))
2999 DENOM=1.+X1+X2-X2*cotso(K)
3000 cotso(K+1)=X1/DENOM
3001 rhtso(K+1)=(FT+X2*rhtso(K))/DENOM
3002 33 CONTINUE
3003 !--- THE NZS element in COTSO and RHTSO will be for snow
3004 !--- There will be 2 layers in snow if it is deeper than DELTSN+SNTH
3005 IF(SNHEI.GT.SNTH) then
3006 if(snhei.le.DELTSN+SNTH) then
3007 !-- 1-layer snow model
3008 ilnb=1
3009 snprim=snhei
3010 soilt1=tso(1)
3011 tsob=tso(1)
3012 XSN = DELT/2./(zshalf(2)+0.5*SNPRIM)
3013 DDZSN = XSN / SNPRIM
3014 X1SN = DDZSN * thdifsn
3015 X2 = DTDZS(1)*THDIFICE(1)
3016 FT = TSO(1)+X1SN*(SOILT-TSO(1)) &
3017 -X2*(TSO(1)-TSO(2))
3018 DENOM = 1. + X1SN + X2 -X2*cotso(NZS1)
3019 cotso(NZS)=X1SN/DENOM
3020 rhtso(NZS)=(FT+X2*rhtso(NZS1))/DENOM
3021 cotsn=cotso(NZS)
3022 rhtsn=rhtso(NZS)
3023 !*** Average temperature of snow pack (C)
3024 tsnav=0.5*(soilt+tso(1)) &
3025 -273.15
3027 else
3028 !-- 2 layers in snow, SOILT1 is temperasture at DELTSN depth
3029 ilnb=2
3030 snprim=deltsn
3031 tsob=soilt1
3032 XSN = DELT/2./(0.5*SNHEI)
3033 XSN1= DELT/2./(zshalf(2)+0.5*(SNHEI-DELTSN))
3034 DDZSN = XSN / DELTSN
3035 DDZSN1 = XSN1 / (SNHEI-DELTSN)
3036 X1SN = DDZSN * thdifsn
3037 X1SN1 = DDZSN1 * thdifsn
3038 X2 = DTDZS(1)*THDIFICE(1)
3039 FT = TSO(1)+X1SN1*(SOILT1-TSO(1)) &
3040 -X2*(TSO(1)-TSO(2))
3041 DENOM = 1. + X1SN1 + X2 - X2*cotso(NZS1)
3042 cotso(nzs)=x1sn1/denom
3043 rhtso(nzs)=(ft+x2*rhtso(nzs1))/denom
3044 ftsnow = soilt1+x1sn*(soilt-soilt1) &
3045 -x1sn1*(soilt1-tso(1))
3046 denomsn = 1. + X1SN + X1SN1 - X1SN1*cotso(NZS)
3047 cotsn=x1sn/denomsn
3048 rhtsn=(ftsnow+X1SN1*rhtso(NZS))/denomsn
3049 !*** Average temperature of snow pack (C)
3050 tsnav=0.5/snhei*((soilt+soilt1)*deltsn &
3051 +(soilt1+tso(1))*(SNHEI-DELTSN)) &
3052 -273.15
3053 endif
3054 ENDIF
3056 IF(SNHEI.LE.SNTH.AND.SNHEI.GT.0.) then
3057 !--- snow is too thin to be treated separately, therefore it
3058 !--- is combined with the first sea ice layer.
3059 fsn=SNHEI/(SNHEI+zsmain(2))
3060 fso=1.-fsn
3061 soilt1=tso(1)
3062 tsob=tso(2)
3063 snprim=SNHEI+zsmain(2)
3064 XSN = DELT/2./((zshalf(3)-zsmain(2))+0.5*snprim)
3065 DDZSN = XSN /snprim
3066 X1SN = DDZSN * (fsn*thdifsn+fso*thdifice(1))
3067 X2=DTDZS(2)*THDIFICE(2)
3068 FT=TSO(2)+X1SN*(SOILT-TSO(2))- &
3069 X2*(TSO(2)-TSO(3))
3070 denom = 1. + x1sn + x2 - x2*cotso(nzs-2)
3071 cotso(nzs1) = x1sn/denom
3072 rhtso(nzs1)=(FT+X2*rhtso(NZS-2))/denom
3073 tsnav=0.5*(soilt+tso(1)) &
3074 -273.15
3075 ENDIF
3077 !************************************************************************
3078 !--- THE HEAT BALANCE EQUATION
3079 !18apr08 nmelt is the flag for melting, and SNOH is heat of snow phase changes
3080 nmelt=0
3081 SNOH=0.
3083 EPOT=-QKMS*(QVATM-QSG)
3084 RHCS=CAPICE(1)
3085 H=1.
3086 FKT=TKMS
3087 D1=cotso(NZS1)
3088 D2=rhtso(NZS1)
3089 TN=SOILT
3090 D9=THDIFICE(1)*RHCS*dzstop
3091 D10=TKMS*CP*RHO
3092 R211=.5*CONFLX/DELT
3093 R21=R211*CP*RHO
3094 R22=.5/(THDIFICE(1)*DELT*dzstop**2)
3095 R6=EMISS *STBOLT*.5*TN**4
3096 R7=R6/TN
3097 D11=RNET+R6
3099 IF(SNHEI.GT.SNTH) THEN
3100 D1SN = cotsn
3101 D2SN = rhtsn
3102 D9SN= THDIFSN*RHOCSN / SNPRIM
3103 R22SN = SNPRIM*SNPRIM*0.5/(THDIFSN*DELT)
3104 ENDIF
3106 IF(SNHEI.LE.SNTH.AND.SNHEI.GT.0.) then
3107 !--- thin snow is combined with sea ice
3108 D1SN = D1
3109 D2SN = D2
3110 D9SN = (fsn*THDIFSN*RHOCSN+fso*THDIFICE(1)*RHCS)/ &
3111 snprim
3112 R22SN = snprim*snprim*0.5 &
3113 /((fsn*THDIFSN+fso*THDIFICE(1))*delt)
3114 ENDIF
3116 IF(SNHEI.eq.0.)then
3117 !--- all snow is sublimated
3118 D9SN = D9
3119 R22SN = R22
3120 D1SN = D1
3121 D2SN = D2
3122 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
3123 print *,' SNHEI = 0, D9SN,R22SN,D1SN,D2SN: ',D9SN,R22SN,D1SN,D2SN
3124 ENDIF
3125 ENDIF
3127 !---- TDENOM for snow
3128 !18apr08 - the iteration start point
3129 212 continue
3130 TDENOM = D9SN*(1.-D1SN +R22SN)+D10+R21+R7 &
3131 +RAINF*CVW*PRCPMS &
3132 +RHOnewCSN*NEWSNOW/DELT
3134 FKQ=QKMS*RHO
3135 R210=R211*RHO
3136 AA=XLVM*(BETA*FKQ+R210)/TDENOM
3137 BB=(D10*TABS+R21*TN+XLVM*(QVATM* &
3138 (BETA*FKQ) &
3139 +R210*QVG)+D11+D9SN*(D2SN+R22SN*TN) &
3140 +RAINF*CVW*PRCPMS*max(273.15,TABS) &
3141 + RHOnewCSN*NEWSNOW/DELT*min(273.15,TABS) &
3142 !18apr08 - add heat of snow phase change
3143 -SNOH &
3144 )/TDENOM
3145 AA1=AA
3146 PP=PATM*1.E3
3147 AA1=AA1/PP
3148 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
3149 print *,'VILKA-SNOW on SEAICE'
3150 print *,'tn,aa1,bb,pp,fkq,r210', &
3151 tn,aa1,bb,pp,fkq,r210
3152 ENDIF
3154 CALL VILKA(TN,AA1,BB,PP,QS1,TS1,TBQ,KTAU,i,j,iland,isoil)
3155 !--- it is saturation over snow
3156 QVG=QS1
3157 QSG=QS1
3158 QCG=0.
3160 !--- SOILT - skin temperature
3161 SOILT=TS1
3163 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
3164 print *,' AFTER VILKA-SNOW on SEAICE'
3165 print *,' TS1,QS1: ', ts1,qs1
3166 ENDIF
3167 ! Solution for temperature at 7.5 cm depth and snow-seaice interface
3168 IF(SNHEI.GT.SNTH) THEN
3169 if(snhei.gt.DELTSN+SNTH) then
3170 !-- 2-layer snow model
3171 SOILT1=rhtsn+cotsn*SOILT
3172 TSO(1)=min(271.4,(rhtso(NZS)+cotso(NZS)*SOILT1))
3173 tsob=soilt1
3174 else
3175 !-- 1 layer in snow
3176 TSO(1)=min(271.4,(rhtso(NZS)+cotso(NZS)*SOILT))
3177 SOILT1=TSO(1)
3178 tsob=tso(1)
3179 endif
3180 ELSE
3181 tso(1)=min(271.4,(tso(2)+(soilt-tso(2))*fso))
3182 SOILT1=TSO(1)
3183 tsob=tso(2)
3184 ENDIF
3185 IF(SNHEI.EQ.0.) THEN
3186 !-- all snow is evaporated
3187 SOILT=min(271.4,SOILT)
3188 TSO(1)=SOILT
3189 SOILT1=SOILT
3190 tsob=SOILT
3191 ENDIF
3192 !---- Final solution for TSO in sea ice
3193 DO K=2,NZS
3194 KK=NZS-K+1
3195 TSO(K)=rhtso(KK)+cotso(KK)*TSO(K-1)
3196 END DO
3197 !--- For thin snow layer combined with the top sea ice layer
3198 !--- TSO is computed by linear inmterpolation between SOILT
3199 !--- and TSO(2)
3201 if(nmelt.eq.1) go to 220
3203 !--- IF SOILT > 273.15 F then melting of snow can happen
3204 IF(SOILT.GT.273.15.AND.SNHEI.GT.0.) THEN
3205 nmelt = 1
3206 soiltfrac=snowfrac*273.15+(1.-snowfrac)*SOILT
3207 QSG= QSN(soiltfrac,TBQ)/PP
3208 QVG=QSG
3209 T3 = STBOLT*SOILTfrac*SOILTfrac*SOILTfrac
3210 UPFLUX = T3 * SOILTfrac
3211 XINET = EMISS*(GLW-UPFLUX)
3212 RNET = GSW + XINET
3213 EPOT = -QKMS*(QVATM-QSG)
3214 Q1=EPOT*RAS
3216 IF (Q1.LE.0.) THEN
3217 ! --- condensation
3218 DEW=-EPOT
3220 QFX= XLVM*RHO*DEW
3221 EETA=QFX/XLVM
3222 ELSE
3223 ! --- evaporation
3224 EETA = Q1 * BETA *1.E3
3225 ! to convert from kg m-2 s-1 to m s-1: 1/rho water=1.e-3************
3226 QFX= - XLVM * EETA
3227 ENDIF
3229 HFX=D10*(TABS-soiltfrac)
3231 IF(SNHEI.GT.SNTH)then
3232 SOH=thdifsn*RHOCSN*(soiltfrac-TSOB)/SNPRIM
3233 SNFLX=SOH
3234 ELSE
3235 SOH=(fsn*thdifsn*rhocsn+fso*thdifice(1)*rhcs)* &
3236 (soiltfrac-TSOB)/snprim
3237 SNFLX=SOH
3238 ENDIF
3239 X= (R21+D9SN*R22SN)*(soiltfrac-TNOLD) + &
3240 XLVM*R210*(QSG-QGOLD)
3241 !-- SNOH is energy flux of snow phase change
3242 SNOH=RNET+QFX +HFX &
3243 +RHOnewCSN*NEWSNOW/DELT*(min(273.15,TABS)-TN) &
3244 -SOH-X+RAINF*CVW*PRCPMS* &
3245 (max(273.15,TABS)-TN)
3246 SNOH=AMAX1(0.,SNOH)
3247 !-- SMELT is speed of melting in M/S
3248 SMELT= SNOH /XLMELT*1.E-3
3249 SMELT=AMIN1(SMELT,SNWEPR/DELT-BETA*EPOT*RAS)
3250 !18apr08 - Egglston limit
3251 ! SMELT= amin1 (smelt, 5.6E-8*meltfactor*max(1.,(soilt-273.15)))
3252 !-- pre 17 nov09 SMELT= amin1 (smelt, 5.6E-7*meltfactor*max(1.,(soilt-273.15)))
3253 SMELT= amin1 (smelt, 9.6E-8*meltfactor*max(1.,(soilt-273.15)))
3254 !*** From Koren et al. (1999) 13% of snow melt stays in the snow pack
3255 !!! rsm=0.13*smelt*delt
3256 if(snwepr.gt.0.) then
3257 rsmfrac=min(0.18,(max(0.08,0.10/snwe*0.13)))
3258 ! else
3259 ! rsmfrac=0.13
3260 endif
3262 rsm=rsmfrac*smelt*delt
3264 SNOHGNEW=SMELT*XLMELT*1.E3
3265 SNODIF=AMAX1(0.,(SNOH-SNOHGNEW))
3267 SNOH=SNOHGNEW
3269 !18apr08 rsm is part of melted water that stays in snow as liquid
3270 SMELT=AMAX1(0.,SMELT-rsm/delt)
3272 !18apr08 - if snow melt occurred then go into iteration for energy budget
3273 ! solution
3274 !-- correction of liquid equivalent of snow depth
3275 !-- due to evaporation and snow melt
3276 SNWE = AMAX1(0.,(SNWEPR- &
3277 (SMELT+BETA*EPOT*RAS)*DELT &
3278 ) )
3280 !--- If all snow melts, then 13% of snow melt we kept in the
3281 !--- snow pack should be added back to snow melt and infiltrate
3282 !--- into soil.
3283 if(snwe.le.rsm) then
3284 smelt=smelt+rsm/delt
3285 snwe=0.
3286 rsm=0.
3287 else
3288 !*** Correct snow density on effect of snow melt, melted
3289 !*** from the top of the snow. 13% of melted water
3290 !*** remains in the pack and changes its density.
3291 !*** Eq. 9 (with my correction) in Koren et al. (1999)
3293 if(snwe.gt.0.) then
3294 xsn=(rhosn*(snwe-rsm)+1.e3*rsm)/ &
3295 snwe
3296 rhosn=MIN(XSN,400.)
3298 RHOCSN=2090.* RHOSN
3299 thdifsn = 0.265/RHOCSN
3300 endif
3302 endif
3304 !--- If there is no snow melting then just evaporation
3305 !--- or condensation cxhanges SNWE
3306 ELSE
3307 EPOT=-QKMS*(QVATM-QSG)
3308 SNWE = AMAX1(0.,(SNWEPR- &
3309 BETA*EPOT*RAS*DELT))
3311 ENDIF
3312 !*** Correct snow density on effect of snow melt, melted
3313 !*** from the top of the snow. 13% of melted water
3314 !*** remains in the pack and changes its density.
3315 !*** Eq. 9 (with my correction) in Koren et al. (1999)
3317 SNHEI=SNWE *1.E3 / RHOSN
3319 snweprint=snwe
3320 ! &
3321 !--- if VEGFRAC.ne.0. then some snow stays on the canopy
3322 !--- and should be added to SNWE for water conservation
3323 ! 4 Nov 07 +VEGFRAC*cst
3324 snheiprint=snweprint*1.E3 / RHOSN
3326 if(nmelt.eq.1) goto 212
3327 220 continue
3329 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
3330 print *, 'snweprint : ',snweprint
3331 print *, 'D9SN,SOILT,TSOB : ', D9SN,SOILT,TSOB
3332 ENDIF
3333 !--- Compute flux in the top snow layer
3334 SNFLX=D9SN*(SOILT-TSOB)
3335 IF(SNHEI.GT.0.) THEN
3336 if(ilnb.gt.1) then
3337 tsnav=0.5/snhei*((soilt+soilt1)*deltsn &
3338 +(soilt1+tso(1))*(SNHEI-DELTSN)) &
3339 -273.15
3340 else
3341 tsnav=0.5*(soilt+tso(1)) - 273.15
3342 endif
3343 ENDIF
3344 !--- RECALCULATION OF DEW USING NEW VALUE OF QSG
3345 DEW=0.
3346 PP=PATM*1.E3
3347 QSG= QSN(SOILT,TBQ)/PP
3348 EPOT = -FQ*(QVATM-QSG)
3349 IF(EPOT.LT.0.) THEN
3350 ! Sublimation
3351 DEW=-EPOT
3352 ENDIF
3353 !-- Restore sea-ice parameters if snow is less than threshold
3354 IF(SNHEI.EQ.0.) then
3355 tsnav=soilt-273.15
3356 smelt=smelt+snwe/delt
3357 rsm=0.
3358 emiss=1.
3359 znt=0.011
3360 alb=0.55
3361 ENDIF
3363 SNOM=SNOM+SMELT*DELT*1.e3
3365 !--- THE DIAGNOSTICS OF SURFACE FLUXES
3367 T3 = STBOLT*SOILT*SOILT*SOILT
3368 UPFLUX = T3 *SOILT
3369 XINET = EMISS*(GLW-UPFLUX)
3370 RNET = GSW + XINET
3371 HFT=-TKMS*CP*RHO*(TABS-SOILT) &
3372 *(P1000mb*0.00001/Patm)**ROVCP
3373 Q1 = - FQ*RAS* (QVATM - QSG)
3374 IF (Q1.LE.0.) THEN
3375 ! --- condensation
3376 !-- moisture flux for coupling with PBL
3377 EETA=-QKMS*RAS*(QVATM/(1.+QVATM) - QSG/(1.+QSG))*1.E3
3378 QFX= XLVm*EETA
3379 !-- actual moisture flux from RUC LSM
3380 DEW=QKMS*(QVATM-QSG)
3381 EETA= RHO*DEW
3382 sublim = EETA
3383 ELSE
3384 ! --- evaporation
3385 !-- moisture flux for coupling with PBL
3386 EETA=-QKMS*RAS*(QVATM/(1.+QVATM) - QVG/(1.+QVG))*1.E3
3387 ! EETA=-QKMS*RAS*(QVATM/(1.+QVATM) - QSG/(1.+QSG))*1.E3
3388 ! to convert from kg m-2 s-1 to m s-1: 1/rho water=1.e-3************
3389 QFX= XLVm * EETA
3390 !-- actual moisture flux from RUC LSM
3391 EETA = Q1*1.E3
3392 sublim = EETA
3393 ENDIF
3395 s=THDIFICE(1)*CAPICE(1)*dzstop*(tso(1)-tso(2))
3396 ! s=D9SN*(SOILT-TSOB)
3397 HFX=HFT
3398 FLTOT=RNET-HFT-QFX-S
3399 !------------------------------------------------------------------------
3400 !------------------------------------------------------------------------
3401 END SUBROUTINE SNOWSEAICE
3402 !------------------------------------------------------------------------
3405 SUBROUTINE SOILTEMP( &
3406 !--- input variables
3407 i,j,iland,isoil, &
3408 delt,ktau,conflx,nzs,nddzs,nroot, &
3409 PRCPMS,RAINF,PATM,TABS,QVATM,QCATM, &
3410 EMISS,RNET, &
3411 QKMS,TKMS,PC,RHO,VEGFRAC, &
3412 THDIF,CAP,DRYCAN,WETCAN, &
3413 TRANSUM,DEW,MAVAIL, &
3414 !--- soil fixed fields
3415 DQM,QMIN,BCLH, &
3416 ZSMAIN,ZSHALF,DTDZS,TBQ, &
3417 !--- constants
3418 XLV,CP,G0_P,CVW,STBOLT, &
3419 !--- output variables
3420 TSO,SOILT,QVG,QSG,QCG)
3422 !*************************************************************
3423 ! Energy budget equation and heat diffusion eqn are
3424 ! solved here and
3425 !
3426 ! DELT - time step (s)
3427 ! ktau - numver of time step
3428 ! CONFLX - depth of constant flux layer (m)
3429 ! IME, JME, KME, NZS - dimensions of the domain
3430 ! NROOT - number of levels within the root zone
3431 ! PRCPMS - precipitation rate in m/s
3432 ! COTSO, RHTSO - coefficients for implicit solution of
3433 ! heat diffusion equation
3434 ! THDIF - thermal diffusivity (m^2/s)
3435 ! QSG,QVG,QCG - saturated mixing ratio, mixing ratio of
3436 ! water vapor and cloud at the ground
3437 ! surface, respectively (kg/kg)
3438 ! PATM - pressure [bar]
3439 ! QC3D,QV3D - cloud and water vapor mixing ratio
3440 ! at the first atm. level (kg/kg)
3441 ! EMISS,RNET - emissivity (0-1) of the ground surface and net
3442 ! radiation at the surface (W/m^2)
3443 ! QKMS - exchange coefficient for water vapor in the
3444 ! surface layer (m/s)
3445 ! TKMS - exchange coefficient for heat in the surface
3446 ! layer (m/s)
3447 ! PC - plant coefficient (resistance)
3448 ! RHO - density of atmosphere near surface (kg/m^3)
3449 ! VEGFRAC - greeness fraction (0-1)
3450 ! CAP - volumetric heat capacity (J/m^3/K)
3451 ! DRYCAN - dry fraction of vegetated area where
3452 ! transpiration may take place (0-1)
3453 ! WETCAN - fraction of vegetated area covered by canopy
3454 ! water (0-1)
3455 ! TRANSUM - transpiration function integrated over the
3456 ! rooting zone (m)
3457 ! DEW - dew in kg/m^2s
3458 ! MAVAIL - fraction of maximum soil moisture in the top
3459 ! layer (0-1)
3460 ! ZSMAIN - main levels in soil (m)
3461 ! ZSHALF - middle of the soil layers (m)
3462 ! DTDZS - dt/(2.*dzshalf*dzmain)
3463 ! TBQ - table to define saturated mixing ration
3464 ! of water vapor for given temperature and pressure
3465 ! TSO - soil temperature (K)
3466 ! SOILT - skin temperature (K)
3467 !
3468 !****************************************************************
3470 IMPLICIT NONE
3471 !-----------------------------------------------------------------
3473 !--- input variables
3475 INTEGER, INTENT(IN ) :: nroot,ktau,nzs , &
3476 nddzs !nddzs=2*(nzs-2)
3477 INTEGER, INTENT(IN ) :: i,j,iland,isoil
3478 REAL, INTENT(IN ) :: DELT,CONFLX,PRCPMS, RAINF
3479 REAL, INTENT(INOUT) :: DRYCAN,WETCAN,TRANSUM
3480 !--- 3-D Atmospheric variables
3481 REAL, &
3482 INTENT(IN ) :: PATM, &
3483 QVATM, &
3484 QCATM
3485 !--- 2-D variables
3486 REAL , &
3487 INTENT(IN ) :: &
3488 EMISS, &
3489 RHO, &
3490 RNET, &
3491 PC, &
3492 VEGFRAC, &
3493 DEW, &
3494 QKMS, &
3495 TKMS
3497 !--- soil properties
3498 REAL , &
3499 INTENT(IN ) :: &
3500 BCLH, &
3501 DQM, &
3502 QMIN
3504 REAL, INTENT(IN ) :: CP, &
3505 CVW, &
3506 XLV, &
3507 STBOLT, &
3508 TABS, &
3509 G0_P
3512 REAL, DIMENSION(1:NZS), INTENT(IN) :: ZSMAIN, &
3513 ZSHALF, &
3514 THDIF, &
3515 CAP
3517 REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS
3519 REAL, DIMENSION(1:4001), INTENT(IN) :: TBQ
3522 !--- input/output variables
3523 !-------- 3-d soil moisture and temperature
3524 REAL, DIMENSION( 1:nzs ) , &
3525 INTENT(INOUT) :: TSO
3527 !-------- 2-d variables
3528 REAL , &
3529 INTENT(INOUT) :: &
3530 MAVAIL, &
3531 QVG, &
3532 QSG, &
3533 QCG, &
3534 SOILT
3537 !--- Local variables
3539 REAL :: x,x1,x2,x4,dzstop,can,ft,sph , &
3540 tn,trans,umveg,denom
3542 REAL :: FKT,D1,D2,D9,D10,DID,R211,R21,R22,R6,R7,D11 , &
3543 PI,H,FKQ,R210,AA,BB,PP,Q1,QS1,TS1,TQ2,TX2 , &
3544 TDENOM
3546 REAL :: C,CC,AA1,RHCS,H1
3548 REAL, DIMENSION(1:NZS) :: cotso,rhtso
3550 INTEGER :: nzs1,nzs2,k,k1,kn,kk
3552 !-----------------------------------------------------------------
3555 NZS1=NZS-1
3556 NZS2=NZS-2
3557 dzstop=1./(ZSMAIN(2)-ZSMAIN(1))
3559 do k=1,nzs
3560 cotso(k)=0.
3561 rhtso(k)=0.
3562 enddo
3563 !******************************************************************************
3564 ! COEFFICIENTS FOR THOMAS ALGORITHM FOR TSO
3565 !******************************************************************************
3566 ! did=2.*(ZSMAIN(nzs)-ZSHALF(nzs))
3567 ! h1=DTDZS(8)*THDIF(nzs-1)*(ZSHALF(nzs)-ZSHALF(nzs-1))/did
3568 ! cotso(1)=h1/(1.+h1)
3569 ! rhtso(1)=(tso(nzs)+h1*(tso(nzs-1)-tso(nzs)))/
3570 ! 1 (1.+h1)
3571 cotso(1)=0.
3572 rhtso(1)=TSO(NZS)
3573 DO 33 K=1,NZS2
3574 KN=NZS-K
3575 K1=2*KN-3
3576 X1=DTDZS(K1)*THDIF(KN-1)
3577 X2=DTDZS(K1+1)*THDIF(KN)
3578 FT=TSO(KN)+X1*(TSO(KN-1)-TSO(KN)) &
3579 -X2*(TSO(KN)-TSO(KN+1))
3580 DENOM=1.+X1+X2-X2*cotso(K)
3581 cotso(K+1)=X1/DENOM
3582 rhtso(K+1)=(FT+X2*rhtso(K))/DENOM
3583 33 CONTINUE
3585 !************************************************************************
3586 !--- THE HEAT BALANCE EQUATION (Smirnova et al., 1996, EQ. 21,26)
3588 RHCS=CAP(1)
3589 H=MAVAIL
3590 IF(DEW.NE.0.)THEN
3591 DRYCAN=0.
3592 WETCAN=1.
3593 ENDIf
3594 TRANS=PC*TRANSUM*DRYCAN/ZSHALF(NROOT+1)
3595 CAN=WETCAN+TRANS
3596 UMVEG=1.-VEGFRAC
3597 FKT=TKMS
3598 D1=cotso(NZS1)
3599 D2=rhtso(NZS1)
3600 TN=SOILT
3601 D9=THDIF(1)*RHCS*dzstop
3602 D10=TKMS*CP*RHO
3603 R211=.5*CONFLX/DELT
3604 R21=R211*CP*RHO
3605 R22=.5/(THDIF(1)*DELT*dzstop**2)
3606 R6=EMISS *STBOLT*.5*TN**4
3607 R7=R6/TN
3608 D11=RNET+R6
3609 TDENOM=D9*(1.-D1+R22)+D10+R21+R7 &
3610 +RAINF*CVW*PRCPMS
3611 FKQ=QKMS*RHO
3612 R210=R211*RHO
3613 C=VEGFRAC*FKQ*CAN
3614 CC=C*XLV/TDENOM
3615 AA=XLV*(FKQ*UMVEG+R210)/TDENOM
3616 BB=(D10*TABS+R21*TN+XLV*(QVATM* &
3617 (FKQ*UMVEG+C) &
3618 +R210*QVG)+D11+D9*(D2+R22*TN) &
3619 +RAINF*CVW*PRCPMS*max(273.15,TABS) &
3620 )/TDENOM
3621 AA1=AA+CC
3622 PP=PATM*1.E3
3623 AA1=AA1/PP
3624 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
3625 PRINT *,' VILKA-1'
3626 print *,'D10,TABS,R21,TN,QVATM,FKQ,UMVEG,VEGFRAC,CAN', &
3627 D10,TABS,R21,TN,QVATM,FKQ,UMVEG,VEGFRAC,CAN
3628 print *,'RNET, EMISS, STBOLT, SOILT',RNET, EMISS, STBOLT, SOILT
3629 print *,'R210,QVG,D11,D9,D2,R22,RAINF,CVW,PRCPMS,TDENOM', &
3630 R210,QVG,D11,D9,D2,R22,RAINF,CVW,PRCPMS,TDENOM
3631 print *,'tn,aa1,bb,pp,umveg,fkq,r210,vegfrac', &
3632 tn,aa1,bb,pp,umveg,fkq,r210,vegfrac
3633 ENDIF
3634 CALL VILKA(TN,AA1,BB,PP,QS1,TS1,TBQ,KTAU,i,j,iland,isoil)
3635 TQ2=QVATM
3636 TX2=TQ2*(1.-H)
3637 Q1=TX2+H*QS1
3638 IF(Q1.LT.QS1) GOTO 100
3639 !--- if no saturation - goto 100
3640 !--- if saturation - goto 90
3641 90 QVG=QS1
3642 QSG=QS1
3643 TSO(1)=TS1
3644 QCG=max(0.,Q1-QS1)
3645 GOTO 200
3646 100 BB=BB-AA*TX2
3647 AA=(AA*H+CC)/PP
3648 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
3649 PRINT *,' VILKA-2'
3650 print *,'D10,TABS,R21,TN,QVATM,FKQ,UMVEG,VEGFRAC,CAN', &
3651 D10,TABS,R21,TN,QVATM,FKQ,UMVEG,VEGFRAC,CAN
3652 print *,'R210,QVG,D11,D9,D2,R22,RAINF,CVW,PRCPMS,TDENOM', &
3653 R210,QVG,D11,D9,D2,R22,RAINF,CVW,PRCPMS,TDENOM
3655 print *,'tn,aa1,bb,pp,umveg,fkq,r210,vegfrac', &
3656 tn,aa1,bb,pp,umveg,fkq,r210,vegfrac
3657 ENDIF
3659 CALL VILKA(TN,AA,BB,PP,QS1,TS1,TBQ,KTAU,i,j,iland,isoil)
3660 Q1=TX2+H*QS1
3661 IF(Q1.GT.QS1) GOTO 90
3662 QSG=QS1
3663 QVG=Q1
3664 TSO(1)=TS1
3665 QCG=0.
3666 200 CONTINUE
3668 !--- SOILT - skin temperature
3669 SOILT=TS1
3671 !---- Final solution for soil temperature - TSO
3672 DO K=2,NZS
3673 KK=NZS-K+1
3674 TSO(K)=rhtso(KK)+cotso(KK)*TSO(K-1)
3675 END DO
3677 !--------------------------------------------------------------------
3678 END SUBROUTINE SOILTEMP
3679 !--------------------------------------------------------------------
3682 SUBROUTINE SNOWTEMP( &
3683 !--- input variables
3684 i,j,iland,isoil, &
3685 delt,ktau,conflx,nzs,nddzs,nroot, &
3686 snwe,snwepr,snhei,newsnow,snowfrac, &
3687 beta,deltsn,snth,rhosn,rhonewsn,meltfactor, & ! add meltfactor
3688 PRCPMS,RAINF, &
3689 PATM,TABS,QVATM,QCATM, &
3690 GLW,GSW,EMISS,RNET, &
3691 QKMS,TKMS,PC,RHO,VEGFRAC, &
3692 THDIF,CAP,DRYCAN,WETCAN,CST, &
3693 TRANF,TRANSUM,DEW,MAVAIL, &
3694 !--- soil fixed fields
3695 DQM,QMIN,PSIS,BCLH, &
3696 ZSMAIN,ZSHALF,DTDZS,TBQ, &
3697 !--- constants
3698 XLVM,CP,rovcp,G0_P,CVW,STBOLT, &
3699 !--- output variables
3700 SNWEPRINT,SNHEIPRINT,RSM, &
3701 TSO,SOILT,SOILT1,TSNAV,QVG,QSG,QCG, &
3702 SMELT,SNOH,SNFLX,ILNB)
3704 !********************************************************************
3705 ! Energy budget equation and heat diffusion eqn are
3706 ! solved here to obtain snow and soil temperatures
3707 !
3708 ! DELT - time step (s)
3709 ! ktau - numver of time step
3710 ! CONFLX - depth of constant flux layer (m)
3711 ! IME, JME, KME, NZS - dimensions of the domain
3712 ! NROOT - number of levels within the root zone
3713 ! PRCPMS - precipitation rate in m/s
3714 ! COTSO, RHTSO - coefficients for implicit solution of
3715 ! heat diffusion equation
3716 ! THDIF - thermal diffusivity (W/m/K)
3717 ! QSG,QVG,QCG - saturated mixing ratio, mixing ratio of
3718 ! water vapor and cloud at the ground
3719 ! surface, respectively (kg/kg)
3720 ! PATM - pressure [bar]
3721 ! QCATM,QVATM - cloud and water vapor mixing ratio
3722 ! at the first atm. level (kg/kg)
3723 ! EMISS,RNET - emissivity (0-1) of the ground surface and net
3724 ! radiation at the surface (W/m^2)
3725 ! QKMS - exchange coefficient for water vapor in the
3726 ! surface layer (m/s)
3727 ! TKMS - exchange coefficient for heat in the surface
3728 ! layer (m/s)
3729 ! PC - plant coefficient (resistance)
3730 ! RHO - density of atmosphere near surface (kg/m^3)
3731 ! VEGFRAC - greeness fraction (0-1)
3732 ! CAP - volumetric heat capacity (J/m^3/K)
3733 ! DRYCAN - dry fraction of vegetated area where
3734 ! transpiration may take place (0-1)
3735 ! WETCAN - fraction of vegetated area covered by canopy
3736 ! water (0-1)
3737 ! TRANSUM - transpiration function integrated over the
3738 ! rooting zone (m)
3739 ! DEW - dew in kg/m^2/s
3740 ! MAVAIL - fraction of maximum soil moisture in the top
3741 ! layer (0-1)
3742 ! ZSMAIN - main levels in soil (m)
3743 ! ZSHALF - middle of the soil layers (m)
3744 ! DTDZS - dt/(2.*dzshalf*dzmain)
3745 ! TBQ - table to define saturated mixing ration
3746 ! of water vapor for given temperature and pressure
3747 ! TSO - soil temperature (K)
3748 ! SOILT - skin temperature (K)
3749 !
3750 !*********************************************************************
3752 IMPLICIT NONE
3753 !---------------------------------------------------------------------
3754 !--- input variables
3756 INTEGER, INTENT(IN ) :: nroot,ktau,nzs , &
3757 nddzs !nddzs=2*(nzs-2)
3759 INTEGER, INTENT(IN ) :: i,j,iland,isoil
3760 REAL, INTENT(IN ) :: DELT,CONFLX,PRCPMS , &
3761 RAINF,NEWSNOW,DELTSN,SNTH , &
3762 TABS,TRANSUM,SNWEPR , &
3763 rhonewsn,meltfactor
3764 real :: rhonewcsn
3766 !--- 3-D Atmospheric variables
3767 REAL, &
3768 INTENT(IN ) :: PATM, &
3769 QVATM, &
3770 QCATM
3771 !--- 2-D variables
3772 REAL , &
3773 INTENT(IN ) :: GLW, &
3774 GSW, &
3775 RHO, &
3776 PC, &
3777 VEGFRAC, &
3778 QKMS, &
3779 TKMS
3781 !--- soil properties
3782 REAL , &
3783 INTENT(IN ) :: &
3784 BCLH, &
3785 DQM, &
3786 PSIS, &
3787 QMIN
3789 REAL, INTENT(IN ) :: CP, &
3790 ROVCP, &
3791 CVW, &
3792 STBOLT, &
3793 XLVM, &
3794 G0_P
3797 REAL, DIMENSION(1:NZS), INTENT(IN) :: ZSMAIN, &
3798 ZSHALF, &
3799 THDIF, &
3800 CAP, &
3801 TRANF
3803 REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS
3805 REAL, DIMENSION(1:4001), INTENT(IN) :: TBQ
3808 !--- input/output variables
3809 !-------- 3-d soil moisture and temperature
3810 REAL, DIMENSION( 1:nzs ) , &
3811 INTENT(INOUT) :: TSO
3814 !-------- 2-d variables
3815 REAL , &
3816 INTENT(INOUT) :: DEW, &
3817 CST, &
3818 RHOSN, &
3819 EMISS, &
3820 MAVAIL, &
3821 QVG, &
3822 QSG, &
3823 QCG, &
3824 SNWE, &
3825 SNHEI, &
3826 SNOWFRAC, &
3827 SMELT, &
3828 SNOH, &
3829 SNFLX, &
3830 SOILT, &
3831 SOILT1, &
3832 TSNAV
3834 REAL, INTENT(INOUT) :: DRYCAN, WETCAN
3836 REAL, INTENT(OUT) :: RSM, &
3837 SNWEPRINT, &
3838 SNHEIPRINT
3839 INTEGER, INTENT(OUT) :: ilnb
3840 !--- Local variables
3843 INTEGER :: nzs1,nzs2,k,k1,kn,kk
3845 REAL :: x,x1,x2,x4,dzstop,can,ft,sph, &
3846 tn,trans,umveg,denom
3848 REAL :: cotsn,rhtsn,xsn1,ddzsn1,x1sn1,ftsnow,denomsn
3850 REAL :: t3,upflux,xinet,ras, &
3851 xlmelt,rhocsn,thdifsn, &
3852 beta,epot,xsn,ddzsn,x1sn,d1sn,d2sn,d9sn,r22sn
3854 REAL :: fso,fsn, &
3855 FKT,D1,D2,D9,D10,DID,R211,R21,R22,R6,R7,D11, &
3856 PI,H,FKQ,R210,AA,BB,PP,Q1,QS1,TS1,TQ2,TX2, &
3857 TDENOM,C,CC,AA1,RHCS,H1, &
3858 tsob, snprim, sh1, sh2, &
3859 smeltg,snohg,snodif,soh, &
3860 CMC2MS,TNOLD,QGOLD,SNOHGNEW
3862 REAL, DIMENSION(1:NZS) :: transp,cotso,rhtso
3863 REAL :: edir1, &
3864 ec1, &
3865 ett1, &
3866 eeta, &
3867 s, &
3868 qfx, &
3869 hfx
3871 REAL :: RNET,rsmfrac,soiltfrac,hsn
3872 integer :: nmelt
3874 !-----------------------------------------------------------------
3876 do k=1,nzs
3877 transp (k)=0.
3878 cotso (k)=0.
3879 rhtso (k)=0.
3880 enddo
3882 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
3883 print *, 'SNOWTEMP: SNHEI,SNTH,SOILT1: ',SNHEI,SNTH,SOILT1,soilt
3884 ENDIF
3885 XLMELT=3.35E+5
3886 RHOCSN=2090.* RHOSN
3887 !18apr08 - add rhonewcsn
3888 RHOnewCSN=2090.* RHOnewSN
3889 THDIFSN = 0.265/RHOCSN
3890 RAS=RHO*1.E-3
3892 SOILTFRAC=SOILT
3894 SMELT=0.
3895 SOH=0.
3896 SMELTG=0.
3897 SNOHG=0.
3898 SNODIF=0.
3899 RSM = 0.
3900 RSMFRAC = 0.
3901 fsn=1.
3902 fso=0.
3903 hsn=snhei
3905 NZS1=NZS-1
3906 NZS2=NZS-2
3908 QGOLD=QVG
3909 TNOLD=SOILT
3910 DZSTOP=1./(ZSMAIN(2)-ZSMAIN(1))
3912 !******************************************************************************
3913 ! COEFFICIENTS FOR THOMAS ALGORITHM FOR TSO
3914 !******************************************************************************
3915 ! did=2.*(ZSMAIN(nzs)-ZSHALF(nzs))
3916 ! h1=DTDZS(8)*THDIF(nzs-1)*(ZSHALF(nzs)-ZSHALF(nzs-1))/did
3917 ! cotso(1)=h1/(1.+h1)
3918 ! rhtso(1)=(tso(nzs)+h1*(tso(nzs-1)-tso(nzs)))/
3919 ! 1 (1.+h1)
3921 cotso(1)=0.
3922 rhtso(1)=TSO(NZS)
3923 DO 33 K=1,NZS2
3924 KN=NZS-K
3925 K1=2*KN-3
3926 X1=DTDZS(K1)*THDIF(KN-1)
3927 X2=DTDZS(K1+1)*THDIF(KN)
3928 FT=TSO(KN)+X1*(TSO(KN-1)-TSO(KN)) &
3929 -X2*(TSO(KN)-TSO(KN+1))
3930 DENOM=1.+X1+X2-X2*cotso(K)
3931 cotso(K+1)=X1/DENOM
3932 rhtso(K+1)=(FT+X2*rhtso(K))/DENOM
3933 33 CONTINUE
3934 !--- THE NZS element in COTSO and RHTSO will be for snow
3935 !--- There will be 2 layers in snow if it is deeper than DELTSN+SNTH
3936 IF(SNHEI.GT.SNTH) then
3937 if(snhei.le.DELTSN+SNTH) then
3938 !-- 1-layer snow model
3939 ilnb=1
3940 snprim=snhei
3941 tsob=tso(1)
3942 XSN = DELT/2./(zshalf(2)+0.5*SNPRIM)
3943 DDZSN = XSN / SNPRIM
3944 X1SN = DDZSN * thdifsn
3945 X2 = DTDZS(1)*THDIF(1)
3946 FT = TSO(1)+X1SN*(SOILT-TSO(1)) &
3947 -X2*(TSO(1)-TSO(2))
3948 DENOM = 1. + X1SN + X2 -X2*cotso(NZS1)
3949 cotso(NZS)=X1SN/DENOM
3950 rhtso(NZS)=(FT+X2*rhtso(NZS1))/DENOM
3951 cotsn=cotso(NZS)
3952 rhtsn=rhtso(NZS)
3953 !*** Average temperature of snow pack (C)
3954 tsnav=0.5*(soilt+tso(1)) &
3955 -273.15
3957 else
3958 !-- 2 layers in snow, SOILT1 is temperasture at DELTSN depth
3959 ilnb=2
3960 snprim=deltsn
3961 tsob=soilt1
3962 XSN = DELT/2./(0.5*SNPRIM)
3963 XSN1= DELT/2./(zshalf(2)+0.5*(SNHEI-DELTSN))
3964 DDZSN = XSN / DELTSN
3965 DDZSN1 = XSN1 / (SNHEI-DELTSN)
3966 X1SN = DDZSN * thdifsn
3967 X1SN1 = DDZSN1 * thdifsn
3968 X2 = DTDZS(1)*THDIF(1)
3969 FT = TSO(1)+X1SN1*(SOILT1-TSO(1)) &
3970 -X2*(TSO(1)-TSO(2))
3971 DENOM = 1. + X1SN1 + X2 - X2*cotso(NZS1)
3972 cotso(nzs)=x1sn1/denom
3973 rhtso(nzs)=(ft+x2*rhtso(nzs1))/denom
3974 ftsnow = soilt1+x1sn*(soilt-soilt1) &
3975 -x1sn1*(soilt1-tso(1))
3976 denomsn = 1. + X1SN + X1SN1 - X1SN1*cotso(NZS)
3977 cotsn=x1sn/denomsn
3978 rhtsn=(ftsnow+X1SN1*rhtso(NZS))/denomsn
3979 !*** Average temperature of snow pack (C)
3980 tsnav=0.5/snhei*((soilt+soilt1)*deltsn &
3981 +(soilt1+tso(1))*(SNHEI-DELTSN)) &
3982 -273.15
3983 endif
3984 ENDIF
3986 IF(SNHEI.LE.SNTH.AND.SNHEI.GT.0.) then
3987 !--- snow is too thin to be treated separately, therefore it
3988 !--- is combined with the first soil layer.
3989 fsn=SNHEI/(SNHEI+zsmain(2))
3990 fso=1.-fsn
3991 tsob=tso(2)
3992 snprim=SNHEI+zsmain(2)
3993 XSN = DELT/2./((zshalf(3)-zsmain(2))+0.5*snprim)
3994 DDZSN = XSN /snprim
3995 X1SN = DDZSN * (fsn*thdifsn+fso*thdif(1))
3996 X2=DTDZS(2)*THDIF(2)
3997 FT=TSO(2)+X1SN*(SOILT-TSO(2))- &
3998 X2*(TSO(2)-TSO(3))
3999 denom = 1. + x1sn + x2 - x2*cotso(nzs-2)
4000 cotso(nzs1) = x1sn/denom
4001 rhtso(nzs1)=(FT+X2*rhtso(NZS-2))/denom
4002 tsnav=0.5*(soilt+tso(1)) &
4003 -273.15
4004 ENDIF
4006 !************************************************************************
4007 !--- THE HEAT BALANCE EQUATION (Smirnova et al. 1996, EQ. 21,26)
4008 !18apr08 nmelt is the flag for melting, and SNOH is heat of snow phase changes
4009 nmelt=0
4010 SNOH=0.
4013 ETT1=0.
4014 EPOT=-QKMS*(QVATM-QSG)
4015 RHCS=CAP(1)
4016 H=MAVAIL
4017 IF(DEW.NE.0.)THEN
4018 DRYCAN=0.
4019 WETCAN=1.
4020 ENDIF
4021 TRANS=PC*TRANSUM*DRYCAN/ZSHALF(NROOT+1)
4022 CAN=WETCAN+TRANS
4023 UMVEG=1.-VEGFRAC
4024 FKT=TKMS
4025 D1=cotso(NZS1)
4026 D2=rhtso(NZS1)
4027 TN=SOILT
4028 D9=THDIF(1)*RHCS*dzstop
4029 D10=TKMS*CP*RHO
4030 R211=.5*CONFLX/DELT
4031 R21=R211*CP*RHO
4032 R22=.5/(THDIF(1)*DELT*dzstop**2)
4033 R6=EMISS *STBOLT*.5*TN**4
4034 R7=R6/TN
4035 D11=RNET+R6
4037 IF(SNHEI.GT.SNTH) THEN
4038 D1SN = cotsn
4039 D2SN = rhtsn
4040 D9SN= THDIFSN*RHOCSN / SNPRIM
4041 R22SN = SNPRIM*SNPRIM*0.5/(THDIFSN*DELT)
4042 ENDIF
4044 IF(SNHEI.LE.SNTH.AND.SNHEI.GT.0.) then
4045 !--- thin snow is combined with soil
4046 D1SN = D1
4047 D2SN = D2
4048 D9SN = (fsn*THDIFSN*RHOCSN+fso*THDIF(1)*RHCS)/ &
4049 snprim
4050 R22SN = snprim*snprim*0.5 &
4051 /((fsn*THDIFSN+fso*THDIF(1))*delt)
4052 ENDIF
4054 IF(SNHEI.eq.0.)then
4055 !--- all snow is sublimated
4056 D9SN = D9
4057 R22SN = R22
4058 D1SN = D1
4059 D2SN = D2
4060 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
4061 print *,' SNHEI = 0, D9SN,R22SN,D1SN,D2SN: ',D9SN,R22SN,D1SN,D2SN
4062 ENDIF
4063 ENDIF
4065 !---- TDENOM for snow
4066 !18apr08 - the iteration start point
4067 212 continue
4068 TDENOM = D9SN*(1.-D1SN +R22SN)+D10+R21+R7 &
4069 +RAINF*CVW*PRCPMS &
4070 +RHOnewCSN*NEWSNOW/DELT
4072 FKQ=QKMS*RHO
4073 R210=R211*RHO
4074 C=VEGFRAC*FKQ*CAN
4075 CC=C*XLVM/TDENOM
4076 AA=XLVM*(BETA*FKQ*UMVEG+R210)/TDENOM
4077 BB=(D10*TABS+R21*TN+XLVM*(QVATM* &
4078 (BETA*FKQ*UMVEG+C) &
4079 +R210*QVG)+D11+D9SN*(D2SN+R22SN*TN) &
4080 +RAINF*CVW*PRCPMS*max(273.15,TABS) &
4081 + RHOnewCSN*NEWSNOW/DELT*min(273.15,TABS) &
4082 !18apr08 - added heat of snow phase change computed in the first iteration
4083 -SNOH &
4084 )/TDENOM
4085 AA1=AA+CC
4086 PP=PATM*1.E3
4087 AA1=AA1/PP
4088 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
4089 print *,'VILKA-SNOW'
4090 print *,'tn,aa1,bb,pp,umveg,fkq,r210,vegfrac', &
4091 tn,aa1,bb,pp,umveg,fkq,r210,vegfrac
4092 ENDIF
4094 CALL VILKA(TN,AA1,BB,PP,QS1,TS1,TBQ,KTAU,i,j,iland,isoil)
4095 !--- it is saturation over snow
4096 QVG=QS1
4097 QSG=QS1
4098 QCG=0.
4099 !--- SOILT - skin temperature
4100 SOILT=TS1
4102 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
4103 print *,' AFTER VILKA-SNOW'
4104 print *,' TS1,QS1: ', ts1,qs1
4105 ENDIF
4107 ! Solution for temperature at 7.5 cm depth and snow-soil interface
4108 IF(SNHEI.GT.SNTH) THEN
4109 if(snhei.gt.DELTSN+SNTH) then
4110 !-- 2-layer snow model
4111 SOILT1=rhtsn+cotsn*SOILT
4112 TSO(1)=rhtso(NZS)+cotso(NZS)*SOILT1
4113 tsob=soilt1
4114 else
4115 !-- 1 layer in snow
4116 TSO(1)=rhtso(NZS)+cotso(NZS)*SOILT
4117 SOILT1=TSO(1)
4118 tsob=tso(1)
4119 endif
4120 ELSE
4121 tso(1)=tso(2)+(soilt-tso(2))*fso
4122 SOILT1=TSO(1)
4123 tsob=tso(2)
4124 ENDIF
4126 IF(snhei.eq.0.) THEN
4127 !-- all snow is evaporated
4128 TSO(1)=SOILT
4129 SOILT1=SOILT
4130 tsob=SOILT
4131 ENDIF
4133 !---- Final solution for TSO
4134 DO K=2,NZS
4135 KK=NZS-K+1
4136 TSO(K)=rhtso(KK)+cotso(KK)*TSO(K-1)
4137 END DO
4138 !--- For thin snow layer combined with the top soil layer
4139 !--- TSO is computed by linear inmterpolation between SOILT
4140 !--- and TSO(2)
4141 if(nmelt.eq.1) go to 220
4143 !--- IF SOILT > 273.15 F then melting of snow can happen
4144 IF(SOILT.GT.273.15.AND.SNHEI.GT.0.) THEN
4145 nmelt = 1
4146 soiltfrac=snowfrac*273.15+(1.-snowfrac)*SOILT
4147 QSG= QSN(soiltfrac,TBQ)/PP
4148 QVG=QSG
4149 T3 = STBOLT*SOILTfrac*SOILTfrac*SOILTfrac
4150 UPFLUX = T3 * SOILTfrac
4151 XINET = EMISS*(GLW-UPFLUX)
4152 RNET = GSW + XINET
4153 EPOT = -QKMS*(QVATM-QSG)
4154 Q1=EPOT*RAS
4156 IF (Q1.LE.0.) THEN
4157 ! --- condensation
4158 DEW=-EPOT
4159 DO K=1,NZS
4160 TRANSP(K)=0.
4161 ENDDO
4163 QFX= XLVM*RHO*DEW
4164 EETA=QFX/XLVM
4165 ELSE
4166 ! --- evaporation
4167 DO K=1,NROOT
4168 TRANSP(K)=-VEGFRAC*q1 &
4169 *PC*TRANF(K)*DRYCAN/zshalf(NROOT+1)
4170 IF(TRANSP(K).GT.0.) TRANSP(K)=0.
4171 ETT1=ETT1-TRANSP(K)
4172 ENDDO
4173 DO k=nroot+1,nzs
4174 transp(k)=0.
4175 enddo
4177 EDIR1 = Q1*UMVEG * BETA
4178 EC1 = Q1 * WETCAN *VEGFRAC
4179 CMC2MS=CST/DELT
4180 EC1=MIN(CMC2MS,EC1)
4181 EETA = (EDIR1 + EC1 + ETT1)*1.E3
4182 ! to convert from kg m-2 s-1 to m s-1: 1/rho water=1.e-3************
4183 QFX= - XLVM * EETA
4184 ENDIF
4186 HFX=D10*(TABS-soiltfrac)
4188 IF(SNHEI.GT.SNTH)then
4189 SOH=thdifsn*RHOCSN*(soiltfrac-TSOB)/SNPRIM
4190 SNFLX=SOH
4191 ELSE
4192 SOH=(fsn*thdifsn*rhocsn+fso*thdif(1)*rhcs)* &
4193 (soiltfrac-TSOB)/snprim
4194 SNFLX=SOH
4195 ENDIF
4197 X= (R21+D9SN*R22SN)*(soiltfrac-TNOLD) + &
4198 XLVM*R210*(QSG-QGOLD)
4199 !-- SNOH is energy flux of snow phase change
4200 SNOH=RNET+QFX +HFX &
4201 +RHOnewCSN*NEWSNOW/DELT*(min(273.15,TABS)-TN) &
4202 -SOH-X+RAINF*CVW*PRCPMS* &
4203 (max(273.15,TABS)-TN)
4204 SNOH=AMAX1(0.,SNOH)
4205 !-- SMELT is speed of melting in M/S
4206 SMELT= SNOH /XLMELT*1.E-3
4207 ! SMELT=AMIN1(SMELT,SNWEPR/DELT-BETA*EPOT*RAS*UMVEG)
4208 SMELT=AMIN1(SMELT,SNWEPR/DELT-BETA*EPOT*RAS)
4209 SMELT=AMAX1(0.,SMELT)
4211 !*** From Koren et al. (1999) 13% of snow melt stays in the snow pack
4212 !!! rsm=0.13*smelt*delt
4213 if(snwepr.gt.0.) then
4214 rsmfrac=min(0.18,(max(0.08,0.10/snwepr*0.13)))
4215 ! else
4216 ! rsmfrac=0.13
4217 endif
4219 rsm=rsmfrac*smelt*delt
4220 !18apr08 - Egglston limit
4221 SMELT= amin1 (smelt, 5.6E-7*meltfactor*max(1.,(soilt-273.15)))
4222 ! SMELT= amin1 (smelt, 5.6E-8*meltfactor*max(1.,(soilt-273.15)))
4224 SNOHGNEW=SMELT*XLMELT*1.E3
4225 SNODIF=AMAX1(0.,(SNOH-SNOHGNEW))
4227 SNOH=SNOHGNEW
4229 !18apr08 rsm is part of melted water that stays in snow as liquid
4230 SMELT=AMAX1(0.,SMELT-rsm/delt)
4232 !-- correction of liquid equivalent of snow depth
4233 !-- due to evaporation and snow melt
4234 SNWE = AMAX1(0.,(SNWEPR- &
4235 (SMELT+BETA*EPOT*RAS)*DELT &
4236 ! (SMELT+BETA*EPOT*RAS*UMVEG)*DELT &
4237 ) )
4239 !--- If all snow melts, then 13% of snow melt we kept in the
4240 !--- snow pack should be added back to snow melt and infiltrate
4241 !--- into soil.
4242 if(snwe.le.rsm) then
4243 smelt=smelt+rsm/delt
4244 snwe=0.
4245 rsm=0.
4246 else
4247 !*** Correct snow density on effect of snow melt, melted
4248 !*** from the top of the snow. 13% of melted water
4249 !*** remains in the pack and changes its density.
4250 !*** Eq. 9 (with my correction) in Koren et al. (1999)
4251 if(snwe.gt.0.) then
4252 xsn=(rhosn*(snwe-rsm)+1.e3*rsm)/ &
4253 snwe
4254 rhosn=MIN(XSN,400.)
4256 RHOCSN=2090.* RHOSN
4257 thdifsn = 0.265/RHOCSN
4258 endif
4260 endif
4262 !--- If there is no snow melting then just evaporation
4263 !--- or condensation cxhanges SNWE
4264 ELSE
4265 EPOT=-QKMS*(QVATM-QSG)
4266 SNWE = AMAX1(0.,(SNWEPR- &
4267 BETA*EPOT*RAS*DELT))
4268 ! BETA*EPOT*RAS*UMVEG*DELT))
4270 ENDIF
4271 !*** Correct snow density on effect of snow melt, melted
4272 !*** from the top of the snow. 13% of melted water
4273 !*** remains in the pack and changes its density.
4274 !*** Eq. 9 (with my correction) in Koren et al. (1999)
4276 SNHEI=SNWE *1.E3 / RHOSN
4278 !18apr08 - if snow melt occurred then go into iteration for energy budget
4279 ! solution
4280 if(nmelt.eq.1) goto 212
4281 220 continue
4282 !-- Snow melt from the top is done. But if ground surface temperature
4283 !-- is above freezing snow can melt from the bottom. The following
4284 !-- piece of code will check if bottom melting is possible.
4286 IF(TSO(1).GT.273.15.AND.SNHEI.GT.0.) THEN
4287 if (snhei.GE.deltsn+snth) then
4288 hsn = snhei - deltsn
4289 else
4290 hsn = snhei
4291 endif
4293 soiltfrac=snowfrac*273.15+(1.-snowfrac)*TSO(1)
4295 SNOHG=(TSO(1)-soiltfrac)*(RHCS*zshalf(2)+ &
4296 RHOCSN*0.5*hsn) / DELT
4297 SNOHG=AMAX1(0.,SNOHG)
4298 SNODIF=0.
4299 SMELTG=SNOHG/XLMELT*1.E-3
4300 ! Egglston - empirical limit on snow melt from the bottom of snow pack
4301 SMELTG=AMIN1(SMELTG, 5.8e-9)
4303 if(SNWE-SMELTG*DELT.ge.rsm) then
4304 SNWE = AMAX1(0.,SNWE-SMELTG*DELT)
4305 else
4306 smeltg=snwe/delt
4307 snwe=0.
4308 rsm=0.
4309 hsn=0.
4310 endif
4312 SNOHGNEW=SMELTG*XLMELT*1.e3
4313 SNODIF=AMAX1(0.,(SNOHG-SNOHGNEW))
4314 TSO(1)=soiltfrac &
4315 + SNODIF/(RHCS*zshalf(2)+ RHOCSN*0.5*hsn)* DELT
4316 SMELT=SMELT+SMELTG
4317 SNOH=SNOH+SNOHGNEW
4319 ENDIF
4321 SNHEI=SNWE *1.E3 / RHOSN
4323 snweprint=snwe
4324 ! &
4325 !--- if VEGFRAC.ne.0. then some snow stays on the canopy
4326 !--- and should be added to SNWE for water conservation
4327 ! 4 Nov 07 +VEGFRAC*cst
4328 snheiprint=snweprint*1.E3 / RHOSN
4330 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
4331 print *, 'snweprint : ',snweprint
4332 print *, 'D9SN,SOILT,TSOB : ', D9SN,SOILT,TSOB
4333 ENDIF
4334 !--- Compute flux in the top snow layer
4335 SNFLX=D9SN*(SOILT-TSOB)
4336 IF(SNHEI.GT.0.) THEN
4337 if(ilnb.gt.1) then
4338 tsnav=0.5/snhei*((soilt+soilt1)*deltsn &
4339 +(soilt1+tso(1))*(SNHEI-DELTSN)) &
4340 -273.15
4341 else
4342 tsnav=0.5*(soilt+tso(1)) - 273.15
4343 endif
4344 ENDIF
4346 ! return
4347 ! end
4348 !------------------------------------------------------------------------
4349 END SUBROUTINE SNOWTEMP
4350 !------------------------------------------------------------------------
4353 SUBROUTINE SOILMOIST ( &
4354 !--input parameters
4355 DELT,NZS,NDDZS,DTDZS,DTDZS2,RIW, &
4356 ZSMAIN,ZSHALF,DIFFU,HYDRO, &
4357 QSG,QVG,QCG,QCATM,QVATM,PRCP, &
4358 QKMS,TRANSP,DRIP, &
4359 DEW,SMELT,SOILICE,VEGFRAC, &
4360 !--soil properties
4361 DQM,QMIN,REF,KSAT,RAS,INFMAX, &
4362 !--output
4363 SOILMOIS,SOILIQW,MAVAIL,RUNOFF,RUNOFF2,INFILTRP)
4364 !*************************************************************************
4365 ! moisture balance equation and Richards eqn.
4366 ! are solved here
4367 !
4368 ! DELT - time step (s)
4369 ! IME,JME,NZS - dimensions of soil domain
4370 ! ZSMAIN - main levels in soil (m)
4371 ! ZSHALF - middle of the soil layers (m)
4372 ! DTDZS - dt/(2.*dzshalf*dzmain)
4373 ! DTDZS2 - dt/(2.*dzshalf)
4374 ! DIFFU - diffusional conductivity (m^2/s)
4375 ! HYDRO - hydraulic conductivity (m/s)
4376 ! QSG,QVG,QCG - saturated mixing ratio, mixing ratio of
4377 ! water vapor and cloud at the ground
4378 ! surface, respectively (kg/kg)
4379 ! QCATM,QVATM - cloud and water vapor mixing ratio
4380 ! at the first atm. level (kg/kg)
4381 ! PRCP - precipitation rate in m/s
4382 ! QKMS - exchange coefficient for water vapor in the
4383 ! surface layer (m/s)
4384 ! TRANSP - transpiration from the soil layers (m/s)
4385 ! DRIP - liquid water dripping from the canopy to soil (m)
4386 ! DEW - dew in kg/m^2s
4387 ! SMELT - melting rate in m/s
4388 ! SOILICE - volumetric content of ice in soil (m^3/m^3)
4389 ! SOILIQW - volumetric content of liquid water in soil (m^3/m^3)
4390 ! VEGFRAC - greeness fraction (0-1)
4391 ! RAS - ration of air density to soil density
4392 ! INFMAX - maximum infiltration rate (kg/m^2/s)
4393 !
4394 ! SOILMOIS - volumetric soil moisture, 6 levels (m^3/m^3)
4395 ! MAVAIL - fraction of maximum soil moisture in the top
4396 ! layer (0-1)
4397 ! RUNOFF - surface runoff (m/s)
4398 ! RUNOFF2 - underground runoff (m)
4399 ! INFILTRP - point infiltration flux into soil (m/s)
4400 ! /(snow bottom runoff) (mm/s)
4401 !
4402 ! COSMC, RHSMC - coefficients for implicit solution of
4403 ! Richards equation
4404 !******************************************************************
4405 IMPLICIT NONE
4406 !------------------------------------------------------------------
4407 !--- input variables
4408 REAL, INTENT(IN ) :: DELT
4409 INTEGER, INTENT(IN ) :: NZS,NDDZS
4411 ! input variables
4413 REAL, DIMENSION(1:NZS), INTENT(IN ) :: ZSMAIN, &
4414 ZSHALF, &
4415 DIFFU, &
4416 HYDRO, &
4417 TRANSP, &
4418 SOILICE, &
4419 DTDZS2
4421 REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS
4423 REAL, INTENT(IN ) :: QSG,QVG,QCG,QCATM,QVATM , &
4424 QKMS,VEGFRAC,DRIP,PRCP , &
4425 DEW,SMELT , &
4426 DQM,QMIN,REF,KSAT,RAS,RIW
4428 ! output
4430 REAL, DIMENSION( 1:nzs ) , &
4432 INTENT(INOUT) :: SOILMOIS,SOILIQW
4434 REAL, INTENT(INOUT) :: MAVAIL,RUNOFF,RUNOFF2,INFILTRP, &
4435 INFMAX
4437 ! local variables
4439 REAL, DIMENSION( 1:nzs ) :: COSMC,RHSMC
4441 REAL :: DZS,R1,R2,R3,R4,R5,R6,R7,R8,R9,R10
4442 REAL :: REFKDT,REFDK,DELT1,F1MAX,F2MAX
4443 REAL :: F1,F2,FD,KDT,VAL,DDT,PX,FK,FKMAX
4444 REAL :: QQ,UMVEG,INFMAX1,TRANS
4445 REAL :: TOTLIQ,FLX,FLXSAT,QTOT
4446 REAL :: DID,X1,X2,X4,DENOM,Q2,Q4
4447 REAL :: dice,fcr,acrt,frzx,sum,cvfrz
4449 INTEGER :: NZS1,NZS2,K,KK,K1,KN,ialp1,jj,jk
4451 !******************************************************************************
4452 ! COEFFICIENTS FOR THOMAS ALGORITHM FOR SOILMOIS
4453 !******************************************************************************
4454 NZS1=NZS-1
4455 NZS2=NZS-2
4457 118 format(6(10Pf23.19))
4459 do k=1,nzs
4460 cosmc(k)=0.
4461 rhsmc(k)=0.
4462 enddo
4464 DID=(ZSMAIN(NZS)-ZSHALF(NZS))
4465 X1=ZSMAIN(NZS)-ZSMAIN(NZS1)
4467 !7may09 DID=(ZSMAIN(NZS)-ZSHALF(NZS))*2.
4468 ! DENOM=DID/DELT+DIFFU(NZS1)/X1
4469 ! COSMC(1)=DIFFU(NZS1)/X1/DENOM
4470 ! RHSMC(1)=(SOILMOIS(NZS)*DID/DELT
4471 ! 1 +TRANSP(NZS)-(HYDRO(NZS)*SOILMOIS(NZS)
4472 ! 1 -HYDRO(NZS1)*SOILMOIS(NZS1))*DID
4473 ! 1 /X1) /DENOM
4475 DENOM=(1.+DIFFU(nzs1)/X1/DID*DELT+HYDRO(NZS)/(2.*DID)*DELT)
4476 COSMC(1)=DELT*(DIFFU(nzs1)/DID/X1 &
4477 +HYDRO(NZS1)/2./DID)/DENOM
4478 RHSMC(1)=(SOILIQW(NZS)+TRANSP(NZS)*DELT/ &
4479 DID)/DENOM
4481 DO 330 K=1,NZS2
4482 KN=NZS-K
4483 K1=2*KN-3
4484 X4=2.*DTDZS(K1)*DIFFU(KN-1)
4485 X2=2.*DTDZS(K1+1)*DIFFU(KN)
4486 Q4=X4+HYDRO(KN-1)*DTDZS2(KN-1)
4487 Q2=X2-HYDRO(KN+1)*DTDZS2(KN-1)
4488 DENOM=1.+X2+X4-Q2*COSMC(K)
4489 COSMC(K+1)=Q4/DENOM
4490 330 RHSMC(K+1)=(SOILIQW(KN)+Q2*RHSMC(K) &
4491 +TRANSP(KN) &
4492 /(ZSHALF(KN+1)-ZSHALF(KN)) &
4493 *DELT)/DENOM
4495 ! --- MOISTURE BALANCE BEGINS HERE
4497 TRANS=TRANSP(1)
4498 UMVEG=1.-VEGFRAC
4500 RUNOFF=0.
4501 RUNOFF2=0.
4502 DZS=ZSMAIN(2)
4503 R1=COSMC(NZS1)
4504 R2= RHSMC(NZS1)
4505 R3=DIFFU(1)/DZS
4506 R4=R3+HYDRO(1)*.5
4507 R5=R3-HYDRO(2)*.5
4508 R6=QKMS*RAS
4509 !-- Total liquid water available on the top of soil domain
4510 !-- Without snow - 3 sources of water: precipitation,
4511 !-- water dripping from the canopy and dew
4512 !-- With snow - only one source of water - snow melt
4514 191 format (f23.19)
4516 TOTLIQ=UMVEG*PRCP-DRIP/DELT-UMVEG*DEW*RAS-SMELT
4519 FLX=TOTLIQ
4520 INFILTRP=TOTLIQ
4522 ! ----------- FROZEN GROUND VERSION -------------------------
4523 ! REFERENCE FROZEN GROUND PARAMETER, CVFRZ, IS A SHAPE PARAMETER OF
4524 ! AREAL DISTRIBUTION FUNCTION OF SOIL ICE CONTENT WHICH EQUALS 1/CV.
4525 ! CV IS A COEFFICIENT OF SPATIAL VARIATION OF SOIL ICE CONTENT.
4526 ! BASED ON FIELD DATA CV DEPENDS ON AREAL MEAN OF FROZEN DEPTH, AND IT
4527 ! CLOSE TO CONSTANT = 0.6 IF AREAL MEAN FROZEN DEPTH IS ABOVE 20 CM.
4528 ! THAT IS WHY PARAMETER CVFRZ = 3 (INT{1/0.6*0.6})
4529 !
4530 ! Current logic doesn't allow CVFRZ be bigger than 3
4531 CVFRZ = 3.
4533 !-- SCHAAKE/KOREN EXPRESSION for calculation of max infiltration
4534 REFKDT=3.
4535 REFDK=3.4341E-6
4536 DELT1=DELT/86400.
4537 F1MAX=DQM*ZSHALF(2)
4538 F2MAX=DQM*(ZSHALF(3)-ZSHALF(2))
4539 F1=F1MAX*(1.-SOILMOIS(1)/DQM)
4540 DICE=SOILICE(1)*ZSHALF(2)
4541 FD=F1
4542 do k=2,nzs1
4543 DICE=DICE+(ZSHALF(k+1)-ZSHALF(k))*SOILICE(K)
4544 FKMAX=DQM*(ZSHALF(k+1)-ZSHALF(k))
4545 FK=FKMAX*(1.-SOILMOIS(k)/DQM)
4546 FD=FD+FK
4547 enddo
4548 KDT=REFKDT*KSAT/REFDK
4549 VAL=(1.-EXP(-KDT*DELT1))
4550 DDT = FD*VAL
4551 PX= - TOTLIQ * DELT
4552 IF(PX.LT.0.0) PX = 0.0
4553 INFMAX1 = (PX*(DDT/(PX+DDT)))/DELT
4554 ! print *,'INFMAX1=,ksat',infmax1,ksat,f1,f2,kdt,val,ddt,px
4555 ! ----------- FROZEN GROUND VERSION --------------------------
4556 ! REDUCTION OF INFILTRATION BASED ON FROZEN GROUND PARAMETERS
4557 !
4558 ! ------------------------------------------------------------------
4560 FRZX= 0.15*((dqm+qmin)/ref) * (0.412 / 0.468)
4561 FCR = 1.
4562 IF ( DICE .GT. 1.E-2) THEN
4563 ACRT = CVFRZ * FRZX / DICE
4564 SUM = 1.
4565 IALP1 = CVFRZ - 1
4566 DO JK = 1,IALP1
4567 K = 1
4568 DO JJ = JK+1, IALP1
4569 K = K * JJ
4570 END DO
4571 SUM = SUM + (ACRT ** ( CVFRZ-JK)) / FLOAT (K)
4572 END DO
4573 FCR = 1. - EXP(-ACRT) * SUM
4574 END IF
4575 ! print *,'FCR--------',fcr
4576 INFMAX1 = INFMAX1* FCR
4577 ! -------------------------------------------------------------------
4579 INFMAX = MAX(INFMAX1,HYDRO(1)*SOILIQW(1))
4580 INFMAX = MIN(INFMAX, -TOTLIQ)
4582 !----
4583 IF (-TOTLIQ.GT.INFMAX)THEN
4584 RUNOFF=-TOTLIQ-INFMAX
4585 FLX=-INFMAX
4586 ENDIF
4587 ! INFILTRP is total infiltration flux in M/S
4588 INFILTRP=FLX
4589 ! Solution of moisture budget
4590 R7=.5*DZS/DELT
4591 R4=R4+R7
4592 FLX=FLX-SOILIQW(1)*R7
4593 R8=UMVEG*R6
4594 QTOT=QVATM+QCATM
4595 R9=TRANS
4596 R10=QTOT-QSG
4597 !-- evaporation regime
4598 IF(R10.LE.0.) THEN
4599 QQ=(R5*R2-FLX+R9)/(R4-R5*R1-R10*R8/(REF-QMIN))
4600 FLXSAT=-DQM*(R4-R5*R1-R10*R8/(REF-QMIN)) &
4601 +R5*R2+R9
4602 ELSE
4603 !-- dew formation regime
4604 QQ=(R2*R5-FLX+R8*(QTOT-QCG-QVG)+R9)/(R4-R1*R5)
4605 FLXSAT=-DQM*(R4-R1*R5)+R2*R5+R8*(QTOT-QVG-QCG)+R9
4606 END IF
4608 IF(QQ.LT.0.) THEN
4609 SOILIQW (1)=1.e-8
4610 SOILMOIS(1)=1.e-8+soilice(1)*riw
4612 ELSE IF(QQ.GT.DQM) THEN
4613 !-- saturation
4614 SOILIQW (1)=DQM
4615 SOILMOIS(1)=DQM
4616 RUNOFF2=(FLXSAT-FLX)*DELT
4617 RUNOFF=RUNOFF+(FLXSAT-FLX)
4618 ELSE
4619 SOILIQW (1)=min(dqm,max(1.e-8,QQ))
4620 SOILMOIS(1)=max(1.e-8,QQ)+soilice(1)*riw
4621 END IF
4623 !--- FINAL SOLUTION FOR SOILMOIS
4624 DO K=2,NZS
4625 KK=NZS-K+1
4626 QQ=COSMC(KK)*SOILIQW(K-1)+RHSMC(KK)
4628 IF (QQ.LT.0.) THEN
4629 SOILIQW(K) =1.e-8
4630 SOILMOIS(K)=1.e-8 + soilice(k)*riw
4632 ELSE IF(QQ.GT.DQM) THEN
4633 !-- saturation
4634 SOILIQW (K)=DQM
4635 SOILMOIS(K)=DQM
4636 IF(K.EQ.NZS)THEN
4637 RUNOFF2=RUNOFF2+(QQ-DQM)*(ZSMAIN(K)-ZSHALF(K))
4638 ELSE
4639 RUNOFF2=RUNOFF2+(QQ-DQM)*(ZSHALF(K+1)-ZSHALF(K))
4640 ENDIF
4641 ELSE
4642 SOILIQW (K)=min(dqm,max(1.e-8,QQ))
4643 SOILMOIS(K)=max(1.e-8,QQ)+soilice(k)*riw
4644 END IF
4645 END DO
4646 ! MAVAIL=min(1.,SOILMOIS(1)/(REF-QMIN))
4647 MAVAIL=min(1.,SOILMOIS(1)/DQM)
4648 if (MAVAIL.EQ.0.) MAVAIL=.00001
4650 ! RETURN
4651 ! END
4652 !-------------------------------------------------------------------
4653 END SUBROUTINE SOILMOIST
4654 !-------------------------------------------------------------------
4657 SUBROUTINE SOILPROP( &
4658 !--- input variables
4659 nzs,fwsat,lwsat,tav,keepfr, &
4660 soilmois,soiliqw,soilice, &
4661 soilmoism,soiliqwm,soilicem, &
4662 !--- soil fixed fields
4663 QWRTZ,rhocs,dqm,qmin,psis,bclh,ksat, &
4664 !--- constants
4665 riw,xlmelt,CP,G0_P,cvw,ci, &
4666 kqwrtz,kice,kwt, &
4667 !--- output variables
4668 thdif,diffu,hydro,cap)
4670 !******************************************************************
4671 ! SOILPROP computes thermal diffusivity, and diffusional and
4672 ! hydraulic condeuctivities
4673 !******************************************************************
4674 ! NX,NY,NZS - dimensions of soil domain
4675 ! FWSAT, LWSAT - volumetric content of frozen and liquid water
4676 ! for saturated condition at given temperatures (m^3/m^3)
4677 ! TAV - temperature averaged for soil layers (K)
4678 ! SOILMOIS - volumetric soil moisture at the main soil levels (m^3/m^3)
4679 ! SOILMOISM - volumetric soil moisture averaged for layers (m^3/m^3)
4680 ! SOILIQWM - volumetric liquid soil moisture averaged for layers (m^3/m^3)
4681 ! SOILICEM - volumetric content of soil ice averaged for layers (m^3/m^3)
4682 ! THDIF - thermal diffusivity for soil layers (W/m/K)
4683 ! DIFFU - diffusional conductivity (m^2/s)
4684 ! HYDRO - hydraulic conductivity (m/s)
4685 ! CAP - volumetric heat capacity (J/m^3/K)
4686 !
4687 !******************************************************************
4689 IMPLICIT NONE
4690 !-----------------------------------------------------------------
4692 !--- soil properties
4693 INTEGER, INTENT(IN ) :: NZS
4694 REAL , &
4695 INTENT(IN ) :: RHOCS, &
4696 BCLH, &
4697 DQM, &
4698 KSAT, &
4699 PSIS, &
4700 QWRTZ, &
4701 QMIN
4703 REAL, DIMENSION( 1:nzs ) , &
4704 INTENT(IN ) :: SOILMOIS, &
4705 keepfr
4708 REAL, INTENT(IN ) :: CP, &
4709 CVW, &
4710 RIW, &
4711 kqwrtz, &
4712 kice, &
4713 kwt, &
4714 XLMELT, &
4715 G0_P
4719 !--- output variables
4720 REAL, DIMENSION(1:NZS) , &
4721 INTENT(INOUT) :: cap,diffu,hydro , &
4722 thdif,tav , &
4723 soilmoism , &
4724 soiliqw,soilice , &
4725 soilicem,soiliqwm , &
4726 fwsat,lwsat
4728 !--- local variables
4729 REAL, DIMENSION(1:NZS) :: hk,detal,kasat,kjpl
4731 REAL :: x,x1,x2,x4,ws,wd,fact,fach,facd,psif,ci
4732 REAL :: tln,tavln,tn,pf,a,am,ame,h
4733 INTEGER :: nzs1,k
4735 !-- for Johansen thermal conductivity
4736 REAL :: kzero,gamd,kdry,kas,x5,sr,ke
4739 nzs1=nzs-1
4741 !-- Constants for Johansen (1975) thermal conductivity
4742 kzero =2. ! if qwrtz > 0.2
4745 do k=1,nzs
4746 detal (k)=0.
4747 kasat (k)=0.
4748 kjpl (k)=0.
4749 hk (k)=0.
4750 enddo
4752 ws=dqm+qmin
4753 x1=xlmelt/(g0_p*psis)
4754 x2=x1/bclh*ws
4755 x4=(bclh+1.)/bclh
4756 !--- Next 3 lines are for Johansen thermal conduct.
4757 gamd=(1.-ws)*2700.
4758 kdry=(0.135*gamd+64.7)/(2700.-0.947*gamd)
4759 kas=kqwrtz**qwrtz*kzero**(1.-qwrtz)
4761 DO K=1,NZS1
4762 tn=tav(k) - 273.15
4763 wd=ws - riw*soilicem(k)
4764 psif=psis*100.*(wd/(soiliqwm(k)+qmin))**bclh &
4765 * (ws/wd)**3.
4766 !--- PSIF should be in [CM] to compute PF
4767 pf=log10(abs(psif))
4768 fact=1.+riw*soilicem(k)
4769 !--- HK is for McCumber thermal conductivity
4770 IF(PF.LE.5.2) THEN
4771 HK(K)=420.*EXP(-(PF+2.7))*fact
4772 ELSE
4773 HK(K)=.1744*fact
4774 END IF
4776 IF(soilicem(k).NE.0.AND.TN.LT.0.) then
4777 !--- DETAL is taking care of energy spent on freezing or released from
4778 ! melting of soil water
4780 DETAL(K)=273.15*X2/(TAV(K)*TAV(K))* &
4781 (TAV(K)/(X1*TN))**X4
4783 if(keepfr(k).eq.1.) then
4784 detal(k)=0.
4785 endif
4787 ENDIF
4789 !--- Next 10 lines calculate Johansen thermal conductivity KJPL
4790 kasat(k)=kas**(1.-ws)*kice**fwsat(k) &
4791 *kwt**lwsat(k)
4793 X5=(soilmoism(k)+qmin)/ws
4794 if(soilicem(k).eq.0.) then
4795 sr=max(0.101,x5)
4796 ke=log10(sr)+1.
4797 !--- next 2 lines - for coarse soils
4798 ! sr=max(0.0501,x5)
4799 ! ke=0.7*log10(sr)+1.
4800 else
4801 ke=x5
4802 endif
4804 kjpl(k)=ke*(kasat(k)-kdry)+kdry
4806 !--- CAP -volumetric heat capacity
4807 CAP(K)=(1.-WS)*RHOCS &
4808 + (soiliqwm(K)+qmin)*CVW &
4809 + soilicem(K)*CI &
4810 + (dqm-soilmoism(k))*CP*1.2 &
4811 - DETAL(K)*1.e3*xlmelt
4813 a=RIW*soilicem(K)
4815 if((ws-a).lt.0.12)then
4816 diffu(K)=0.
4817 else
4818 H=max(0.,(soilmoism(K)-a)/(max(1.e-8,(dqm-a))))
4819 facd=1.
4820 if(a.ne.0.)facd=1.-a/max(1.e-8,soilmoism(K))
4821 ame=max(1.e-8,dqm-riw*soilicem(K))
4822 !--- DIFFU is diffusional conductivity of soil water
4823 diffu(K)=-BCLH*KSAT*PSIS/ame* &
4824 (dqm/ame)**3. &
4825 *H**(BCLH+2.)*facd
4826 endif
4828 ! diffu(K)=-BCLH*KSAT*PSIS/dqm &
4829 ! *H**(BCLH+2.)
4832 !--- thdif - thermal diffusivity
4833 ! thdif(K)=HK(K)/CAP(K)
4834 !--- Use thermal conductivity from Johansen (1975)
4835 thdif(K)=KJPL(K)/CAP(K)
4837 END DO
4839 DO K=1,NZS
4841 if((ws-riw*soilice(k)).lt.0.12)then
4842 hydro(k)=0.
4843 else
4844 fach=1.
4845 if(soilice(k).ne.0.) &
4846 fach=1.-riw*soilice(k)/max(1.e-8,soilmois(k))
4847 am=max(1.e-8,dqm-riw*soilice(k))
4848 !--- HYDRO is hydraulic conductivity of soil water
4849 hydro(K)=KSAT/am* &
4850 (soiliqw(K)/am) &
4851 **(2.*BCLH+2.) &
4852 * fach
4853 endif
4855 ENDDO
4857 ! RETURN
4858 ! END
4860 !-----------------------------------------------------------------------
4861 END SUBROUTINE SOILPROP
4862 !-----------------------------------------------------------------------
4865 SUBROUTINE TRANSF( &
4866 !--- input variables
4867 nzs,nroot,soiliqw,tabs, &
4868 !--- soil fixed fields
4869 dqm,qmin,ref,wilt,zshalf, &
4870 !--- output variables
4871 tranf,transum)
4873 !-------------------------------------------------------------------
4874 !--- TRANF(K) - THE TRANSPIRATION FUNCTION (Smirnova et al. 1996, EQ. 18,19)
4875 !*******************************************************************
4876 ! NX,NY,NZS - dimensions of soil domain
4877 ! SOILIQW - volumetric liquid soil moisture at the main levels (m^3/m^3)
4878 ! TRANF - the transpiration function at levels (m)
4879 ! TRANSUM - transpiration function integrated over the rooting zone (m)
4880 !
4881 !*******************************************************************
4882 IMPLICIT NONE
4883 !-------------------------------------------------------------------
4885 !--- input variables
4887 INTEGER, INTENT(IN ) :: nroot,nzs
4889 REAL , &
4890 INTENT(IN ) :: TABS
4891 !--- soil properties
4892 REAL , &
4893 INTENT(IN ) :: DQM, &
4894 QMIN, &
4895 REF, &
4896 WILT
4898 REAL, DIMENSION(1:NZS), INTENT(IN) :: soiliqw, &
4899 ZSHALF
4901 !-- output
4902 REAL, DIMENSION(1:NZS), INTENT(OUT) :: TRANF
4903 REAL, INTENT(OUT) :: TRANSUM
4905 !-- local variables
4906 REAL :: totliq, did
4907 INTEGER :: k
4909 !-- for non-linear root distribution
4910 REAL :: gx,sm1,sm2,sm3,sm4,ap0,ap1,ap2,ap3,ap4
4911 REAL :: FTEM
4912 REAL, DIMENSION(1:NZS) :: PART
4913 !--------------------------------------------------------------------
4915 do k=1,nzs
4916 part(k)=0.
4917 enddo
4919 transum=0.
4920 totliq=soiliqw(1)+qmin
4921 sm1=totliq
4922 sm2=sm1*sm1
4923 sm3=sm2*sm1
4924 sm4=sm3*sm1
4925 ap0=0.299
4926 ap1=-8.152
4927 ap2=61.653
4928 ap3=-115.876
4929 ap4=59.656
4930 gx=ap0+ap1*sm1+ap2*sm2+ap3*sm3+ap4*sm4
4931 if(totliq.ge.ref) gx=1.
4932 if(totliq.le.0.) gx=0.
4933 if(gx.gt.1.) gx=1.
4934 if(gx.lt.0.) gx=0.
4935 DID=zshalf(2)
4936 part(1)=DID*gx
4937 !--- air temperature function
4938 ! Avissar (1985) and AX 7/95
4939 IF (TABS .LE. 302.15) THEN
4940 FTEM = 1.0 / (1.0 + EXP(-0.41 * (TABS - 282.05)))
4941 ELSE
4942 FTEM = 1.0 / (1.0 + EXP(0.5 * (TABS - 314.0)))
4943 ENDIF
4944 !---
4945 IF(TOTLIQ.GT.REF) THEN
4946 TRANF(1)=DID
4947 ELSE IF(TOTLIQ.LE.WILT) THEN
4948 TRANF(1)=0.
4949 ELSE
4950 TRANF(1)=(TOTLIQ-WILT)/(REF-WILT)*DID
4951 ENDIF
4952 !-- uncomment next line for non-linear root distribution
4953 !cc TRANF(1)=part(1)
4954 TRANF(1)=TRANF(1)*FTEM
4956 DO K=2,NROOT
4957 totliq=soiliqw(k)+qmin
4958 sm1=totliq
4959 sm2=sm1*sm1
4960 sm3=sm2*sm1
4961 sm4=sm3*sm1
4962 gx=ap0+ap1*sm1+ap2*sm2+ap3*sm3+ap4*sm4
4963 if(totliq.ge.ref) gx=1.
4964 if(totliq.le.0.) gx=0.
4965 if(gx.gt.1.) gx=1.
4966 if(gx.lt.0.) gx=0.
4967 DID=zshalf(K+1)-zshalf(K)
4968 part(k)=did*gx
4969 IF(totliq.GE.REF) THEN
4970 TRANF(K)=DID
4971 ELSE IF(totliq.LE.WILT) THEN
4972 TRANF(K)=0.
4973 ELSE
4974 TRANF(K)=(totliq-WILT) &
4975 /(REF-WILT)*DID
4976 ENDIF
4977 !-- uncomment next line for non-linear root distribution
4978 !cc TRANF(k)=part(k)
4979 TRANF(k)=TRANF(k)*FTEM
4980 END DO
4982 !-- TRANSUM - total for the rooting zone
4983 transum=0.
4984 DO K=1,NROOT
4985 transum=transum+tranf(k)
4986 END DO
4988 ! RETURN
4989 ! END
4990 !-----------------------------------------------------------------
4991 END SUBROUTINE TRANSF
4992 !-----------------------------------------------------------------
4995 SUBROUTINE VILKA(TN,D1,D2,PP,QS,TS,TT,NSTEP,ii,j,iland,isoil)
4996 !--------------------------------------------------------------
4997 !--- VILKA finds the solution of energy budget at the surface
4998 !--- using table T,QS computed from Clausius-Klapeiron
4999 !--------------------------------------------------------------
5000 REAL, DIMENSION(1:4001), INTENT(IN ) :: TT
5001 REAL, INTENT(IN ) :: TN,D1,D2,PP
5002 INTEGER, INTENT(IN ) :: NSTEP,ii,j,iland,isoil
5004 REAL, INTENT(OUT ) :: QS, TS
5006 REAL :: F1,T1,T2,RN
5007 INTEGER :: I,I1
5009 I=(TN-1.7315E2)/.05+1
5010 T1=173.1+FLOAT(I)*.05
5011 F1=T1+D1*TT(I)-D2
5012 I1=I-F1/(.05+D1*(TT(I+1)-TT(I)))
5013 I=I1
5014 IF(I.GT.4000.OR.I.LT.1) GOTO 1
5015 10 I1=I
5016 T1=173.1+FLOAT(I)*.05
5017 F1=T1+D1*TT(I)-D2
5018 RN=F1/(.05+D1*(TT(I+1)-TT(I)))
5019 I=I-INT(RN)
5020 IF(I.GT.4000.OR.I.LT.1) GOTO 1
5021 IF(I1.NE.I) GOTO 10
5022 TS=T1-.05*RN
5023 QS=(TT(I)+(TT(I)-TT(I+1))*RN)/PP
5024 GOTO 20
5025 1 PRINT *,' AVOST IN VILKA '
5026 ! WRITE(12,*)'AVOST',TN,D1,D2,PP,NSTEP
5027 PRINT *,TN,D1,D2,PP,NSTEP,I,TT(i),ii,j,iland,isoil
5028 CALL wrf_error_fatal (' AVOST IN VILKA ' )
5029 20 CONTINUE
5030 ! RETURN
5031 ! END
5032 !-----------------------------------------------------------------------
5033 END SUBROUTINE VILKA
5034 !-----------------------------------------------------------------------
5037 SUBROUTINE SOILVEGIN ( IVGTYP,ISLTYP,MYJ, &
5038 IFOREST,EMISS,PC,ZNT,QWRTZ, &
5039 RHOCS,BCLH,DQM,KSAT,PSIS,QMIN,REF,WILT )
5041 !************************************************************************
5042 ! Set-up soil and vegetation Parameters in the case when
5043 ! snow disappears during the forecast and snow parameters
5044 ! shold be replaced by surface parameters according to
5045 ! soil and vegetation types in this point.
5046 !
5047 ! Output:
5048 !
5049 !
5050 ! Soil parameters:
5051 ! DQM: MAX soil moisture content - MIN (m^3/m^3)
5052 ! REF: Reference soil moisture (m^3/m^3)
5053 ! WILT: Wilting PT soil moisture contents (m^3/m^3)
5054 ! QMIN: Air dry soil moist content limits (m^3/m^3)
5055 ! PSIS: SAT soil potential coefs. (m)
5056 ! KSAT: SAT soil diffusivity/conductivity coefs. (m/s)
5057 ! BCLH: Soil diffusivity/conductivity exponent.
5058 !
5059 ! ************************************************************************
5060 IMPLICIT NONE
5061 !---------------------------------------------------------------------------
5062 integer, parameter :: nsoilclas=19
5063 integer, parameter :: nvegclas=24+3
5064 integer, parameter :: iwater=16
5065 integer, parameter :: ilsnow=99
5068 !--- soiltyp classification according to STATSGO(nclasses=16)
5069 !
5070 ! 1 SAND SAND
5071 ! 2 LOAMY SAND LOAMY SAND
5072 ! 3 SANDY LOAM SANDY LOAM
5073 ! 4 SILT LOAM SILTY LOAM
5074 ! 5 SILT SILTY LOAM
5075 ! 6 LOAM LOAM
5076 ! 7 SANDY CLAY LOAM SANDY CLAY LOAM
5077 ! 8 SILTY CLAY LOAM SILTY CLAY LOAM
5078 ! 9 CLAY LOAM CLAY LOAM
5079 ! 10 SANDY CLAY SANDY CLAY
5080 ! 11 SILTY CLAY SILTY CLAY
5081 ! 12 CLAY LIGHT CLAY
5082 ! 13 ORGANIC MATERIALS LOAM
5083 ! 14 WATER
5084 ! 15 BEDROCK
5085 ! Bedrock is reclassified as class 14
5086 ! 16 OTHER (land-ice)
5087 ! 17 Playa
5088 ! 18 Lava
5089 ! 19 White Sand
5090 !
5091 !----------------------------------------------------------------------
5092 REAL LQMA(nsoilclas),LRHC(nsoilclas), &
5093 LPSI(nsoilclas),LQMI(nsoilclas), &
5094 LBCL(nsoilclas),LKAS(nsoilclas), &
5095 LWIL(nsoilclas),LREF(nsoilclas), &
5096 DATQTZ(nsoilclas)
5097 !-- LQMA Rawls et al.[1982]
5098 ! DATA LQMA /0.417, 0.437, 0.453, 0.501, 0.486, 0.463, 0.398,
5099 ! & 0.471, 0.464, 0.430, 0.479, 0.475, 0.439, 1.0, 0.20, 0.401/
5100 !---
5101 !-- Clapp, R. and G. Hornberger, 1978: Empirical equations for some soil
5102 ! hydraulic properties, Water Resour. Res., 14, 601-604.
5104 !-- Clapp et al. [1978]
5105 DATA LQMA /0.395, 0.410, 0.435, 0.485, 0.485, 0.451, 0.420, &
5106 0.477, 0.476, 0.426, 0.492, 0.482, 0.451, 1.0, &
5107 0.20, 0.435, 0.468, 0.200, 0.339/
5109 !-- LREF Rawls et al.[1982]
5110 ! DATA LREF /0.091, 0.125, 0.207, 0.330, 0.360, 0.270, 0.255,
5111 ! & 0.366, 0.318, 0.339, 0.387, 0.396, 0.329, 1.0, 0.108, 0.283/
5113 !-- Clapp et al. [1978]
5114 DATA LREF /0.174, 0.179, 0.249, 0.369, 0.369, 0.314, 0.299, &
5115 0.357, 0.391, 0.316, 0.409, 0.400, 0.314, 1., &
5116 0.1, 0.249, 0.454, 0.17, 0.236/
5118 !-- LWIL Rawls et al.[1982]
5119 ! DATA LWIL/0.033, 0.055, 0.095, 0.133, 0.133, 0.117, 0.148,
5120 ! & 0.208, 0.197, 0.239, 0.250, 0.272, 0.066, 0.0, 0.006, 0.029/
5122 !-- Clapp et al. [1978]
5123 DATA LWIL/0.068, 0.075, 0.114, 0.179, 0.179, 0.155, 0.175, &
5124 0.218, 0.250, 0.219, 0.283, 0.286, 0.155, 0.0, &
5125 0.006, 0.114, 0.030, 0.006, 0.01/
5127 ! DATA LQMI/0.010, 0.028, 0.047, 0.084, 0.084, 0.066, 0.067,
5128 ! & 0.120, 0.103, 0.100, 0.126, 0.138, 0.066, 0.0, 0.006, 0.028/
5130 !-- Carsel and Parrish [1988]
5131 DATA LQMI/0.045, 0.057, 0.065, 0.067, 0.034, 0.078, 0.10, &
5132 0.089, 0.095, 0.10, 0.070, 0.068, 0.078, 0.0, &
5133 0.004, 0.065, 0.020, 0.004, 0.008/
5135 !-- LPSI Cosby et al[1984]
5136 ! DATA LPSI/0.060, 0.036, 0.141, 0.759, 0.759, 0.355, 0.135,
5137 ! & 0.617, 0.263, 0.098, 0.324, 0.468, 0.355, 0.0, 0.069, 0.036/
5138 ! & 0.617, 0.263, 0.098, 0.324, 0.468, 0.355, 0.0, 0.069, 0.036/
5140 !-- Clapp et al. [1978]
5141 DATA LPSI/0.121, 0.090, 0.218, 0.786, 0.786, 0.478, 0.299, &
5142 0.356, 0.630, 0.153, 0.490, 0.405, 0.478, 0.0, &
5143 0.121, 0.218, 0.468, 0.069, 0.069/
5145 !-- LKAS Rawls et al.[1982]
5146 ! DATA LKAS/5.83E-5, 1.70E-5, 7.19E-6, 1.89E-6, 1.89E-6,
5147 ! & 3.67E-6, 1.19E-6, 4.17E-7, 6.39E-7, 3.33E-7, 2.50E-7,
5148 ! & 1.67E-7, 3.38E-6, 0.0, 1.41E-4, 1.41E-5/
5150 !-- Clapp et al. [1978]
5151 DATA LKAS/1.76E-4, 1.56E-4, 3.47E-5, 7.20E-6, 7.20E-6, &
5152 6.95E-6, 6.30E-6, 1.70E-6, 2.45E-6, 2.17E-6, &
5153 1.03E-6, 1.28E-6, 6.95E-6, 0.0, 1.41E-4, &
5154 3.47E-5, 1.28E-6, 1.41E-4, 1.76E-4/
5156 !-- LBCL Cosby et al [1984]
5157 ! DATA LBCL/2.79, 4.26, 4.74, 5.33, 5.33, 5.25, 6.66,
5158 ! & 8.72, 8.17, 10.73, 10.39, 11.55, 5.25, 0.0, 2.79, 4.26/
5160 !-- Clapp et al. [1978]
5161 DATA LBCL/4.05, 4.38, 4.90, 5.30, 5.30, 5.39, 7.12, &
5162 7.75, 8.52, 10.40, 10.40, 11.40, 5.39, 0.0, &
5163 4.05, 4.90, 11.55, 2.79, 2.79/
5165 DATA LRHC /1.47,1.41,1.34,1.27,1.27,1.21,1.18,1.32,1.23, &
5166 1.18,1.15,1.09,1.21,4.18,2.03,2.10,1.09,2.03,1.47/
5168 DATA DATQTZ/0.92,0.82,0.60,0.25,0.10,0.40,0.60,0.10,0.35, &
5169 0.52,0.10,0.25,0.00,0.,0.60,0.0,0.25,0.60,0.92/
5171 !--------------------------------------------------------------------------
5172 !
5173 ! USGS Vegetation Types
5174 !
5175 ! 1: Urban and Built-Up Land
5176 ! 2: Dryland Cropland and Pasture
5177 ! 3: Irrigated Cropland and Pasture
5178 ! 4: Mixed Dryland/Irrigated Cropland and Pasture
5179 ! 5: Cropland/Grassland Mosaic
5180 ! 6: Cropland/Woodland Mosaic
5181 ! 7: Grassland
5182 ! 8: Shrubland
5183 ! 9: Mixed Shrubland/Grassland
5184 ! 10: Savanna
5185 ! 11: Deciduous Broadleaf Forest
5186 ! 12: Deciduous Needleleaf Forest
5187 ! 13: Evergreen Broadleaf Forest
5188 ! 14: Evergreen Needleleaf Fores
5189 ! 15: Mixed Forest
5190 ! 16: Water Bodies
5191 ! 17: Herbaceous Wetland
5192 ! 18: Wooded Wetland
5193 ! 19: Barren or Sparsely Vegetated
5194 ! 20: Herbaceous Tundra
5195 ! 21: Wooded Tundra
5196 ! 22: Mixed Tundra
5197 ! 23: Bare Ground Tundra
5198 ! 24: Snow or Ice
5199 !
5200 ! 25: Playa
5201 ! 26: Lava
5202 ! 27: White Sand
5205 !---- Below are the arrays for the vegetation parameters
5206 REAL LALB(nvegclas),LMOI(nvegclas),LEMI(nvegclas), &
5207 LROU(nvegclas),LTHI(nvegclas),LSIG(nvegclas), &
5208 LPC(nvegclas), NROTBL(nvegclas)
5210 !************************************************************************
5211 !---- vegetation parameters
5212 !
5213 !-- USGS model
5214 !
5215 DATA LALB/.18,.17,.18,.18,.18,.16,.19,.22,.20,.20,.16,.14, &
5216 .12,.12,.13,.08,.14,.14,.25,.15,.15,.15,.25,.55, &
5217 .30,.16,.60 /
5218 DATA LEMI/.88,4*.92,.93,.92,.88,.9,.92,.93,.94, &
5219 .95,.95,.94,.98,.95,.95,.85,.92,.93,.92,.85,.95, &
5220 .85,.85,.90 /
5221 !-- Roughness length is changed for forests and some others
5222 ! DATA LROU/.5,.06,.075,.065,.05,.2,.075,.1,.11,.15,.8,.85, &
5223 ! 2.0,1.0,.563,.0001,.2,.4,.05,.1,.15,.1,.065,.05/
5224 DATA LROU/.5,.06,.075,.065,.05,.2,.075,.1,.11,.15,.5,.5, &
5225 .5,.5,.5,.0001,.2,.4,.05,.1,.15,.1,.065,.05, &
5226 .01,.15,.01 /
5228 DATA LMOI/.1,.3,.5,.25,.25,.35,.15,.1,.15,.15,.3,.3, &
5229 .5,.3,.3,1.,.6,.35,.02,.5,.5,.5,.02,.95,.40,.50,.40/
5230 !
5231 !---- still needs to be corrected
5232 !
5233 ! DATA LPC/ 15*.8,0.,.8,.8,.5,.5,.5,.5,.5,.0/
5234 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, &
5235 0.3,0.3,0.4,0.4,0.3,0.,.3,0.,0./
5237 ! used in RUC DATA LPC /0.6,6*0.8,0.7,0.75,6*0.8,0.,0.8,0.8, &
5238 ! 0.5,0.7,0.6,0.7,0.5,0./
5241 !***************************************************************************
5244 INTEGER :: &
5245 IVGTYP, &
5246 ISLTYP
5248 LOGICAL, INTENT(IN ) :: myj
5250 REAL , &
5251 INTENT ( OUT) :: pc
5253 REAL , &
5254 INTENT (INOUT ) :: emiss, &
5255 znt
5256 !--- soil properties
5257 REAL , &
5258 INTENT( OUT) :: RHOCS, &
5259 BCLH, &
5260 DQM, &
5261 KSAT, &
5262 PSIS, &
5263 QMIN, &
5264 QWRTZ, &
5265 REF, &
5266 WILT
5268 INTEGER, DIMENSION( 1:(nvegclas) ) , &
5269 INTENT ( OUT) :: iforest
5273 INTEGER, DIMENSION( 1:(nvegclas) ) :: if1
5274 INTEGER :: kstart, kfin, lstart, lfin
5275 INTEGER :: i,j,k
5277 !***********************************************************************
5278 ! DATA ZS1/0.0,0.05,0.20,0.40,1.6,3.0/ ! o - levels in soil
5279 ! DATA ZS2/0.0,0.025,0.125,0.30,1.,2.3/ ! x - levels in soil
5280 DATA IF1/12*0,1,1,1,12*0/
5282 do k=1,nvegclas
5283 iforest(k)=if1(k)
5284 enddo
5287 ! EMISS = LEMI(IVGTYP)
5288 EMISS = LEMITBL(IVGTYP)
5289 ! When MYJ sfc scheme is used - better use recommended in MYJSFCINIT
5290 ! values of roughness length, and not redefine it here.
5291 ! The table in this routine is the one we use in RUC with RUC LSM.
5293 !tgs 3 Oct 2007 - MYJSFCINIT roughness produced unrealistic fluxes,
5294 !tgs - will not use it any more!
5295 !tgs if (.not. myj) then
5296 ! ZNT = LROU(IVGTYP)
5297 ZNT = Z0TBL(IVGTYP)
5298 !tgs endif
5300 ! PC = LPC (IVGTYP)
5301 PC = PCTBL(IVGTYP)
5302 ! RHOCS = LRHC(ISLTYP)*1.E6
5303 RHOCS = HC(ISLTYP)*1.E6
5305 ! parameters from SOILPARM.TBL
5306 BCLH = BB(ISLTYP)
5307 DQM = MAXSMC(ISLTYP)- &
5308 DRYSMC(ISLTYP)
5309 KSAT = SATDK(ISLTYP)
5310 PSIS = - SATPSI(ISLTYP)
5311 QMIN = DRYSMC(ISLTYP)
5312 REF = REFSMC(ISLTYP)
5313 WILT = WLTSMC(ISLTYP)
5314 QWRTZ = QTZ(ISLTYP)
5316 ! parameters from the look-up tables
5317 ! BCLH = LBCL(ISLTYP)
5318 ! DQM = LQMA(ISLTYP)- &
5319 ! LQMI(ISLTYP)
5320 ! KSAT = LKAS(ISLTYP)
5321 ! PSIS = - LPSI(ISLTYP)
5322 ! QMIN = LQMI(ISLTYP)
5323 ! REF = LREF(ISLTYP)
5324 ! WILT = LWIL(ISLTYP)
5325 ! QWRTZ = DATQTZ(ISLTYP)
5327 !--------------------------------------------------------------------------
5328 END SUBROUTINE SOILVEGIN
5329 !--------------------------------------------------------------------------
5332 SUBROUTINE SNOWFREE (ivgtyp,myj,emiss,znt,iland)
5333 !************************************************************************
5334 ! Set-up soil and vegetation Parameters in the case when
5335 ! snow disappears during the forecast and snow parameters
5336 ! shold be replaced by surface parameters according to
5337 ! soil and vegetation types in this point.
5338 !
5339 !***************************************************************************
5340 IMPLICIT NONE
5341 !---------------------------------------------------------------------------
5342 integer, parameter :: nvegclas=24+3
5345 INTEGER :: IVGTYP
5347 LOGICAL, INTENT(IN ) :: myj
5349 REAL, INTENT(INOUT) :: &
5350 emiss, &
5351 znt
5352 INTEGER, INTENT(INOUT) :: ILAND
5354 !---- Below are the arrays for the vegetation parameters
5355 REAL, DIMENSION( 1:nvegclas ) :: LALB, &
5356 LEMI, &
5357 LROU_MYJ,&
5358 LROU
5360 !************************************************************************
5361 !-- USGS model
5362 !
5363 DATA LALB/.18,.17,.18,.18,.18,.16,.19,.22,.20,.20,.16,.14, &
5364 .12,.12,.13,.08,.14,.14,.25,.15,.15,.15,.25,.55, &
5365 .30,.16,.60/
5366 DATA LEMI/.88,4*.92,.93,.92,.88,.9,.92,.93,.94, &
5367 .95,.95,.94,.98,.95,.95,.85,.92,.93,.92,.85,.95, &
5368 .85,.85,.90 /
5369 !-- Roughness length is changed for forests and some others
5370 ! next 2 lines - table from RUC
5371 ! DATA LROU/.5,.06,.075,.065,.05,.2,.075,.1,.11,.15,.8,.85, &
5372 ! 2.0,1.0,.563,.0001,.2,.4,.05,.1,.15,.1,.065,.05/
5374 DATA LROU/.5,.06,.075,.065,.05,.2,.075,.1,.11,.15,.5,.5, &
5375 .5,.5,.5,.0001,.2,.4,.05,.1,.15,.1,.065,.05, &
5376 .01,.15,.01 /
5378 ! With MYJSFC better use the table from MYJSFCINIT
5379 DATA LROU_MYJ/1.0,.07,.07,.07,.07,.15,.08,.03,.05,.86,.8,.85, &
5380 2.65,1.09,.8,.001,.04,.05,.01,.04,.06,.05,.03,.001, &
5381 .01,.15,.01 /
5385 !--------------------------------------------------------------------------
5387 EMISS = LEMITBL(IVGTYP)
5388 !tgs 3 Oct 2007 - LROU_MYJ gives unrealistic surface fluxes with RUC LSM with RUC LSM
5389 ! if(myj) then
5390 ! ZNT = LROU_MYJ(IVGTYP)
5391 ! else
5392 ZNT = Z0TBL(IVGTYP)
5393 !!! ZNT = LROU(IVGTYP)
5394 ! endif
5395 ILAND = IVGTYP
5396 ! ---
5398 ! RETURN
5399 ! END
5400 !--------------------------------------------------------------------------
5401 END SUBROUTINE SNOWFREE
5403 SUBROUTINE RUCLSMINIT( SH2O,SMFR3D,TSLB,SMOIS,ISLTYP,IVGTYP, &
5404 XICE,mavail,nzs, iswater, isice, restart, &
5405 allowed_to_read , &
5406 ids,ide, jds,jde, kds,kde, &
5407 ims,ime, jms,jme, kms,kme, &
5408 its,ite, jts,jte, kts,kte )
5409 #ifdef WRF_CHEM
5410 USE module_data_gocart_dust
5411 #endif
5412 IMPLICIT NONE
5415 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
5416 ims,ime, jms,jme, kms,kme, &
5417 its,ite, jts,jte, kts,kte, &
5418 nzs, iswater, isice
5420 REAL, DIMENSION( ims:ime, 1:nzs, jms:jme ) , &
5421 INTENT(IN) :: TSLB, &
5422 SMOIS
5424 INTEGER, DIMENSION( ims:ime, jms:jme ) , &
5425 INTENT(INOUT) :: ISLTYP,IVGTYP
5427 REAL, DIMENSION( ims:ime, 1:nzs, jms:jme ) , &
5428 INTENT(INOUT) :: SMFR3D, &
5429 SH2O
5431 REAL, DIMENSION( ims:ime, jms:jme ) , &
5432 INTENT(INOUT) :: XICE,MAVAIL
5434 REAL, DIMENSION ( 1:nzs ) :: SOILIQW
5436 LOGICAL , INTENT(IN) :: restart, allowed_to_read
5438 !
5439 INTEGER :: I,J,L,itf,jtf
5440 REAL :: RIW,XLMELT,TLN,DQM,REF,PSIS,QMIN,BCLH
5442 INTEGER :: errflag
5444 ! itf=min0(ite,ide-1)
5445 ! jtf=min0(jte,jde-1)
5448 RIW=900.*1.e-3
5449 XLMELT=3.35E+5
5451 ! initialize three LSM related tables
5452 IF ( allowed_to_read ) THEN
5453 CALL wrf_message( 'INITIALIZE THREE LSM RELATED TABLES' )
5454 CALL RUCLSM_PARM_INIT
5455 ENDIF
5457 #ifdef WRF_CHEM
5458 !
5459 ! need this parameter for dust parameterization in wrf/chem
5460 !
5461 do I=1,NSLTYPE
5462 porosity(i)=maxsmc(i)
5463 enddo
5464 #endif
5466 IF(.not.restart)THEN
5468 itf=min0(ite,ide-1)
5469 jtf=min0(jte,jde-1)
5471 errflag = 0
5472 DO j = jts,jtf
5473 DO i = its,itf
5474 IF ( ISLTYP( i,j ) .LT. 1 ) THEN
5475 errflag = 1
5476 WRITE(err_message,*)"module_sf_ruclsm.F: lsminit: out of range ISLTYP ",i,j,ISLTYP( i,j )
5477 CALL wrf_message(err_message)
5478 ENDIF
5479 ENDDO
5480 ENDDO
5481 IF ( errflag .EQ. 1 ) THEN
5482 CALL wrf_error_fatal( "module_sf_ruclsm.F: lsminit: out of range value "// &
5483 "of ISLTYP. Is this field in the input?" )
5484 ENDIF
5486 DO J=jts,jtf
5487 DO I=its,itf
5489 ! CALL SOILIN ( ISLTYP(I,J), DQM, REF, PSIS, QMIN, BCLH )
5492 !--- Computation of volumetric content of ice in soil
5493 !--- and initialize MAVAIL
5494 DQM = MAXSMC (ISLTYP(I,J)) - &
5495 DRYSMC (ISLTYP(I,J))
5496 REF = REFSMC (ISLTYP(I,J))
5497 PSIS = - SATPSI (ISLTYP(I,J))
5498 QMIN = DRYSMC (ISLTYP(I,J))
5499 BCLH = BB (ISLTYP(I,J))
5502 !!! IF (.not.restart) THEN
5504 IF(xice(i,j).gt.0.) THEN
5505 !-- for ice
5506 DO L=1,NZS
5507 smfr3d(i,l,j)=1.
5508 sh2o(i,l,j)=0.
5509 mavail(i,j) = 1.
5510 ENDDO
5511 ELSE
5512 if(isltyp(i,j).ne.14 ) then
5513 !-- land
5514 mavail(i,j) = max(0.00001,min(1.,smois(i,1,j)/dqm))
5515 ! mavail(i,j) = max(0.00001,min(1.,smois(i,1,j)/(ref-qmin)))
5516 DO L=1,NZS
5517 !-- for land points initialize soil ice
5518 tln=log(TSLB(i,l,j)/273.15)
5520 if(tln.lt.0.) then
5521 soiliqw(l)=(dqm+qmin)*(XLMELT* &
5522 (tslb(i,l,j)-273.15)/tslb(i,l,j)/9.81/psis) &
5523 **(-1./bclh)-qmin
5524 soiliqw(l)=max(0.,soiliqw(l))
5525 soiliqw(l)=min(soiliqw(l),smois(i,l,j))
5526 sh2o(i,l,j)=soiliqw(l)
5527 smfr3d(i,l,j)=(smois(i,l,j)-soiliqw(l))/RIW
5529 else
5530 smfr3d(i,l,j)=0.
5531 sh2o(i,l,j)=smois(i,l,j)
5532 endif
5533 ENDDO
5535 else
5536 !-- for water ISLTYP=14
5537 DO L=1,NZS
5538 smfr3d(i,l,j)=0.
5539 sh2o(i,l,j)=1.
5540 mavail(i,j) = 1.
5541 ENDDO
5542 endif
5543 ENDIF
5545 ENDDO
5546 ENDDO
5548 ENDIF
5550 END SUBROUTINE ruclsminit
5551 !
5552 !-----------------------------------------------------------------
5553 SUBROUTINE RUCLSM_PARM_INIT
5554 !-----------------------------------------------------------------
5556 character*8 :: MMINLU, MMINSL
5558 MMINLU='USGS-RUC'
5559 MMINSL='STAS-RUC'
5560 call RUCLSM_SOILVEGPARM( MMINLU, MMINSL)
5562 !-----------------------------------------------------------------
5563 END SUBROUTINE RUCLSM_PARM_INIT
5564 !-----------------------------------------------------------------
5566 !-----------------------------------------------------------------
5567 SUBROUTINE RUCLSM_SOILVEGPARM( MMINLU, MMINSL)
5568 !-----------------------------------------------------------------
5570 IMPLICIT NONE
5572 integer :: LUMATCH, IINDEX, LC, NUM_SLOPE
5573 integer :: ierr
5574 INTEGER , PARAMETER :: OPEN_OK = 0
5576 character*8 :: MMINLU, MMINSL
5577 character*128 :: mess , message, vege_parm_string
5578 logical, external :: wrf_dm_on_monitor
5581 !-----SPECIFY VEGETATION RELATED CHARACTERISTICS :
5582 ! ALBBCK: SFC albedo (in percentage)
5583 ! Z0: Roughness length (m)
5584 ! LEMI: Emissivity
5585 ! PC: Plant coefficient for transpiration function
5586 ! -- the rest of the parameters are read in but not used currently
5587 ! SHDFAC: Green vegetation fraction (in percentage)
5588 ! Note: The ALBEDO, Z0, and SHDFAC values read from the following table
5589 ! ALBEDO, amd Z0 are specified in LAND-USE TABLE; and SHDFAC is
5590 ! the monthly green vegetation data
5591 ! CMXTBL: MAX CNPY Capacity (m)
5592 ! NROTBL: Rooting depth (layer)
5593 ! RSMIN: Mimimum stomatal resistance (s m-1)
5594 ! RSMAX: Max. stomatal resistance (s m-1)
5595 ! RGL: Parameters used in radiation stress function
5596 ! HS: Parameter used in vapor pressure deficit functio
5597 ! TOPT: Optimum transpiration air temperature. (K)
5598 ! CMCMAX: Maximum canopy water capacity
5599 ! CFACTR: Parameter used in the canopy inteception calculati
5600 ! SNUP: Threshold snow depth (in water equivalent m) that
5601 ! implies 100% snow cover
5602 ! LAI: Leaf area index (dimensionless)
5603 ! MAXALB: Upper bound on maximum albedo over deep snow
5604 !
5605 !-----READ IN VEGETAION PROPERTIES FROM VEGPARM.TBL
5606 !
5608 IF ( wrf_dm_on_monitor() ) THEN
5610 OPEN(19, FILE='VEGPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
5611 IF(ierr .NE. OPEN_OK ) THEN
5612 WRITE(message,FMT='(A)') &
5613 'module_sf_ruclsm.F: soil_veg_gen_parm: failure opening VEGPARM.TBL'
5614 CALL wrf_error_fatal ( message )
5615 END IF
5617 WRITE ( mess, * ) 'INPUT VEGPARM FOR ',MMINLU
5618 CALL wrf_message( mess )
5620 LUMATCH=0
5622 2000 FORMAT (A8)
5623 READ (19,'(A)') vege_parm_string
5624 outer : DO
5625 READ (19,2000,END=2002)LUTYPE
5626 READ (19,*)LUCATS,IINDEX
5628 WRITE( mess , * ) 'VEGPARM FOR ',LUTYPE,' FOUND', LUCATS,' CATEGORIES'
5629 CALL wrf_message( mess )
5631 IF(LUTYPE.NE.MMINLU)THEN ! Skip over the undesired table
5632 write ( mess , * ) 'Skipping ', LUTYPE, ' table'
5633 CALL wrf_message( mess )
5634 DO LC=1,LUCATS
5635 READ (19,*)
5636 ENDDO
5637 inner : DO ! Find the next "Vegetation Parameters"
5638 READ (19,'(A)',END=2002) vege_parm_string
5639 IF (TRIM(vege_parm_string) .EQ. "Vegetation Parameters") THEN
5640 EXIT inner
5641 END IF
5642 ENDDO inner
5643 ELSE
5644 LUMATCH=1
5645 write ( mess , * ) 'Found ', LUTYPE, ' table'
5646 CALL wrf_message( mess )
5647 EXIT outer ! Found the table, read the data
5648 END IF
5650 ENDDO outer
5652 IF (LUMATCH == 1) then
5653 write ( mess , * ) 'Reading ',LUTYPE,' table'
5654 CALL wrf_message( mess )
5655 DO LC=1,LUCATS
5656 READ (19,*)IINDEX,ALBTBL(LC),Z0TBL(LC),LEMITBL(LC),PCTBL(LC), &
5657 SHDTBL(LC),NROTBL(LC),RSTBL(LC),RGLTBL(LC), &
5658 HSTBL(LC),SNUPTBL(LC),LAITBL(LC),MAXALB(LC)
5659 ENDDO
5660 !
5661 READ (19,*)
5662 READ (19,*)TOPT_DATA
5663 READ (19,*)
5664 READ (19,*)CMCMAX_DATA
5665 READ (19,*)
5666 READ (19,*)CFACTR_DATA
5667 READ (19,*)
5668 READ (19,*)RSMAX_DATA
5669 READ (19,*)
5670 READ (19,*)BARE
5671 READ (19,*)
5672 READ (19,*)NATURAL
5673 ENDIF
5675 2002 CONTINUE
5676 CLOSE (19)
5678 IF (LUMATCH == 0) then
5679 CALL wrf_error_fatal ("Land Use Dataset '"//MMINLU//"' not found in VEGPARM.TBL.")
5680 ENDIF
5682 END IF
5684 CALL wrf_dm_bcast_string ( LUTYPE , 8 )
5685 CALL wrf_dm_bcast_integer ( LUCATS , 1 )
5686 CALL wrf_dm_bcast_integer ( IINDEX , 1 )
5687 CALL wrf_dm_bcast_integer ( LUMATCH , 1 )
5688 CALL wrf_dm_bcast_real ( ALBTBL , NLUS )
5689 CALL wrf_dm_bcast_real ( Z0TBL , NLUS )
5690 CALL wrf_dm_bcast_real ( LEMITBL , NLUS )
5691 CALL wrf_dm_bcast_real ( PCTBL , NLUS )
5692 CALL wrf_dm_bcast_real ( SHDTBL , NLUS )
5693 CALL wrf_dm_bcast_real ( NROTBL , NLUS )
5694 CALL wrf_dm_bcast_real ( RSTBL , NLUS )
5695 CALL wrf_dm_bcast_real ( RGLTBL , NLUS )
5696 CALL wrf_dm_bcast_real ( HSTBL , NLUS )
5697 CALL wrf_dm_bcast_real ( SNUPTBL , NLUS )
5698 CALL wrf_dm_bcast_real ( LAITBL , NLUS )
5699 CALL wrf_dm_bcast_real ( MAXALB , NLUS )
5700 CALL wrf_dm_bcast_real ( TOPT_DATA , 1 )
5701 CALL wrf_dm_bcast_real ( CMCMAX_DATA , 1 )
5702 CALL wrf_dm_bcast_real ( CFACTR_DATA , 1 )
5703 CALL wrf_dm_bcast_real ( RSMAX_DATA , 1 )
5704 CALL wrf_dm_bcast_integer ( BARE , 1 )
5706 !
5707 !-----READ IN SOIL PROPERTIES FROM SOILPARM.TBL
5708 !
5709 IF ( wrf_dm_on_monitor() ) THEN
5710 OPEN(19, FILE='SOILPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
5711 IF(ierr .NE. OPEN_OK ) THEN
5712 WRITE(message,FMT='(A)') &
5713 'module_sf_ruclsm.F: soil_veg_gen_parm: failure opening SOILPARM.TBL'
5714 CALL wrf_error_fatal ( message )
5715 END IF
5717 WRITE(mess,*) 'INPUT SOIL TEXTURE CLASSIFICAION = ',MMINSL
5718 CALL wrf_message( mess )
5720 LUMATCH=0
5722 READ (19,*)
5723 READ (19,2000,END=2003)SLTYPE
5724 READ (19,*)SLCATS,IINDEX
5725 IF(SLTYPE.NE.MMINSL)THEN
5726 DO LC=1,SLCATS
5727 READ (19,*) IINDEX,BB(LC),DRYSMC(LC),HC(LC),MAXSMC(LC),&
5728 REFSMC(LC),SATPSI(LC),SATDK(LC), SATDW(LC), &
5729 WLTSMC(LC), QTZ(LC)
5730 ENDDO
5731 ENDIF
5732 READ (19,*)
5733 READ (19,2000,END=2003)SLTYPE
5734 READ (19,*)SLCATS,IINDEX
5736 IF(SLTYPE.EQ.MMINSL)THEN
5737 WRITE( mess , * ) 'SOIL TEXTURE CLASSIFICATION = ',SLTYPE,' FOUND', &
5738 SLCATS,' CATEGORIES'
5739 CALL wrf_message ( mess )
5740 LUMATCH=1
5741 ENDIF
5742 IF(SLTYPE.EQ.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
5750 2003 CONTINUE
5752 CLOSE (19)
5753 ENDIF
5755 CALL wrf_dm_bcast_integer ( LUMATCH , 1 )
5756 CALL wrf_dm_bcast_string ( SLTYPE , 8 )
5757 CALL wrf_dm_bcast_string ( MMINSL , 8 ) ! since this is reset above, see oct2 ^
5758 CALL wrf_dm_bcast_integer ( SLCATS , 1 )
5759 CALL wrf_dm_bcast_integer ( IINDEX , 1 )
5760 CALL wrf_dm_bcast_real ( BB , NSLTYPE )
5761 CALL wrf_dm_bcast_real ( DRYSMC , NSLTYPE )
5762 CALL wrf_dm_bcast_real ( HC , NSLTYPE )
5763 CALL wrf_dm_bcast_real ( MAXSMC , NSLTYPE )
5764 CALL wrf_dm_bcast_real ( REFSMC , NSLTYPE )
5765 CALL wrf_dm_bcast_real ( SATPSI , NSLTYPE )
5766 CALL wrf_dm_bcast_real ( SATDK , NSLTYPE )
5767 CALL wrf_dm_bcast_real ( SATDW , NSLTYPE )
5768 CALL wrf_dm_bcast_real ( WLTSMC , NSLTYPE )
5769 CALL wrf_dm_bcast_real ( QTZ , NSLTYPE )
5771 IF(LUMATCH.EQ.0)THEN
5772 CALL wrf_message( 'SOIl TEXTURE IN INPUT FILE DOES NOT ' )
5773 CALL wrf_message( 'MATCH SOILPARM TABLE' )
5774 CALL wrf_error_fatal ( 'INCONSISTENT OR MISSING SOILPARM FILE' )
5775 ENDIF
5777 !
5778 !-----READ IN GENERAL PARAMETERS FROM GENPARM.TBL
5779 !
5780 IF ( wrf_dm_on_monitor() ) THEN
5781 OPEN(19, FILE='GENPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
5782 IF(ierr .NE. OPEN_OK ) THEN
5783 WRITE(message,FMT='(A)') &
5784 'module_sf_ruclsm.F: soil_veg_gen_parm: failure opening GENPARM.TBL'
5785 CALL wrf_error_fatal ( message )
5786 END IF
5788 READ (19,*)
5789 READ (19,*)
5790 READ (19,*) NUM_SLOPE
5792 SLPCATS=NUM_SLOPE
5794 DO LC=1,SLPCATS
5795 READ (19,*)SLOPE_DATA(LC)
5796 ENDDO
5798 READ (19,*)
5799 READ (19,*)SBETA_DATA
5800 READ (19,*)
5801 READ (19,*)FXEXP_DATA
5802 READ (19,*)
5803 READ (19,*)CSOIL_DATA
5804 READ (19,*)
5805 READ (19,*)SALP_DATA
5806 READ (19,*)
5807 READ (19,*)REFDK_DATA
5808 READ (19,*)
5809 READ (19,*)REFKDT_DATA
5810 READ (19,*)
5811 READ (19,*)FRZK_DATA
5812 READ (19,*)
5813 READ (19,*)ZBOT_DATA
5814 READ (19,*)
5815 READ (19,*)CZIL_DATA
5816 READ (19,*)
5817 READ (19,*)SMLOW_DATA
5818 READ (19,*)
5819 READ (19,*)SMHIGH_DATA
5820 CLOSE (19)
5821 ENDIF
5823 CALL wrf_dm_bcast_integer ( NUM_SLOPE , 1 )
5824 CALL wrf_dm_bcast_integer ( SLPCATS , 1 )
5825 CALL wrf_dm_bcast_real ( SLOPE_DATA , NSLOPE )
5826 CALL wrf_dm_bcast_real ( SBETA_DATA , 1 )
5827 CALL wrf_dm_bcast_real ( FXEXP_DATA , 1 )
5828 CALL wrf_dm_bcast_real ( CSOIL_DATA , 1 )
5829 CALL wrf_dm_bcast_real ( SALP_DATA , 1 )
5830 CALL wrf_dm_bcast_real ( REFDK_DATA , 1 )
5831 CALL wrf_dm_bcast_real ( REFKDT_DATA , 1 )
5832 CALL wrf_dm_bcast_real ( FRZK_DATA , 1 )
5833 CALL wrf_dm_bcast_real ( ZBOT_DATA , 1 )
5834 CALL wrf_dm_bcast_real ( CZIL_DATA , 1 )
5835 CALL wrf_dm_bcast_real ( SMLOW_DATA , 1 )
5836 CALL wrf_dm_bcast_real ( SMHIGH_DATA , 1 )
5839 !-----------------------------------------------------------------
5840 END SUBROUTINE RUCLSM_SOILVEGPARM
5841 !-----------------------------------------------------------------
5844 SUBROUTINE SOILIN (ISLTYP, DQM, REF, PSIS, QMIN, BCLH )
5846 !--- soiltyp classification according to STATSGO(nclasses=16)
5847 !
5848 ! 1 SAND SAND
5849 ! 2 LOAMY SAND LOAMY SAND
5850 ! 3 SANDY LOAM SANDY LOAM
5851 ! 4 SILT LOAM SILTY LOAM
5852 ! 5 SILT SILTY LOAM
5853 ! 6 LOAM LOAM
5854 ! 7 SANDY CLAY LOAM SANDY CLAY LOAM
5855 ! 8 SILTY CLAY LOAM SILTY CLAY LOAM
5856 ! 9 CLAY LOAM CLAY LOAM
5857 ! 10 SANDY CLAY SANDY CLAY
5858 ! 11 SILTY CLAY SILTY CLAY
5859 ! 12 CLAY LIGHT CLAY
5860 ! 13 ORGANIC MATERIALS LOAM
5861 ! 14 WATER
5862 ! 15 BEDROCK
5863 ! Bedrock is reclassified as class 14
5864 ! 16 OTHER (land-ice)
5865 ! extra classes from Fei Chen
5866 ! 17 Playa
5867 ! 18 Lava
5868 ! 19 White Sand
5869 !
5870 !----------------------------------------------------------------------
5871 integer, parameter :: nsoilclas=19
5873 integer, intent ( in) :: isltyp
5874 real, intent ( out) :: dqm,ref,qmin,psis
5876 REAL LQMA(nsoilclas),LREF(nsoilclas),LBCL(nsoilclas), &
5877 LPSI(nsoilclas),LQMI(nsoilclas)
5879 !-- LQMA Rawls et al.[1982]
5880 ! DATA LQMA /0.417, 0.437, 0.453, 0.501, 0.486, 0.463, 0.398,
5881 ! & 0.471, 0.464, 0.430, 0.479, 0.475, 0.439, 1.0, 0.20, 0.401/
5882 !---
5883 !-- Clapp, R. and G. Hornberger, Empirical equations for some soil
5884 ! hydraulic properties, Water Resour. Res., 14,601-604,1978.
5885 !-- Clapp et al. [1978]
5886 DATA LQMA /0.395, 0.410, 0.435, 0.485, 0.485, 0.451, 0.420, &
5887 0.477, 0.476, 0.426, 0.492, 0.482, 0.451, 1.0, &
5888 0.20, 0.435, 0.468, 0.200, 0.339/
5890 !-- Clapp et al. [1978]
5891 DATA LREF /0.174, 0.179, 0.249, 0.369, 0.369, 0.314, 0.299, &
5892 0.357, 0.391, 0.316, 0.409, 0.400, 0.314, 1., &
5893 0.1, 0.249, 0.454, 0.17, 0.236/
5895 !-- Carsel and Parrish [1988]
5896 DATA LQMI/0.045, 0.057, 0.065, 0.067, 0.034, 0.078, 0.10, &
5897 0.089, 0.095, 0.10, 0.070, 0.068, 0.078, 0.0, &
5898 0.004, 0.065, 0.020, 0.004, 0.008/
5900 !-- Clapp et al. [1978]
5901 DATA LPSI/0.121, 0.090, 0.218, 0.786, 0.786, 0.478, 0.299, &
5902 0.356, 0.630, 0.153, 0.490, 0.405, 0.478, 0.0, &
5903 0.121, 0.218, 0.468, 0.069, 0.069/
5905 !-- Clapp et al. [1978]
5906 DATA LBCL/4.05, 4.38, 4.90, 5.30, 5.30, 5.39, 7.12, &
5907 7.75, 8.52, 10.40, 10.40, 11.40, 5.39, 0.0, &
5908 4.05, 4.90, 11.55, 2.79, 2.79/
5911 DQM = LQMA(ISLTYP)- &
5912 LQMI(ISLTYP)
5913 REF = LREF(ISLTYP)
5914 PSIS = - LPSI(ISLTYP)
5915 QMIN = LQMI(ISLTYP)
5916 BCLH = LBCL(ISLTYP)
5918 END SUBROUTINE SOILIN
5920 END MODULE module_sf_ruclsm