1 #define LSMRUC_DBG_LVL 3000
2 !WRF:MODEL_LAYER:PHYSICS
4 MODULE module_sf_ruclsm
6 ! Notes for perturbations of soil properties (Judith Berner)
7 ! Perturbations are applied in subroutine soilprob to array hydro;
8 ! soilprop is called from subroutine SFCTMP which is called from subroutine LSMRUC;
9 ! subroutine LSMRUC had two new 3D fields: pattern_spp_lsm (in) and field_sf(inout);
10 ! their vertical dimension is number of atmospheric levels (kms:kme) - (suboptimal, but easiest hack)
11 ! field_sf is used to pass perturbed fields of hydrop up to model (and output) driver;
12 ! in argument list to SFCTMP the arrays are passed as pattern_spp_lsm(i,1:nzs,j), and exist henceforth as
14 ! in the subroutines below SFCTMP (SNOW and SNOWSOIL) the fields are called rstochcol,fieldcol_sf
15 ! to reflect their dimension rstochcol (1:nzs)
18 USE module_model_constants
21 ! VEGETATION PARAMETERS
22 INTEGER :: LUCATS , BARE, NATURAL, CROP, URBAN
23 integer, PARAMETER :: NLUS=50
25 INTEGER, DIMENSION(1:NLUS) :: IFORTBL
26 real, dimension(1:NLUS) :: SNUPTBL, RSTBL, RGLTBL, HSTBL, LAITBL, &
27 ALBTBL, Z0TBL, LEMITBL, PCTBL, SHDTBL, MAXALB
28 REAL :: TOPT_DATA,CMCMAX_DATA,CFACTR_DATA,RSMAX_DATA
31 INTEGER, PARAMETER :: NSLTYPE=30
33 REAL, DIMENSION (1:NSLTYPE) :: BB,DRYSMC,HC, &
34 MAXSMC, REFSMC,SATPSI,SATDK,SATDW, WLTSMC,QTZ
36 ! LSM GENERAL PARAMETERS
38 INTEGER, PARAMETER :: NSLOPE=30
39 REAL, DIMENSION (1:NSLOPE) :: SLOPE_DATA
40 REAL :: SBETA_DATA,FXEXP_DATA,CSOIL_DATA,SALP_DATA,REFDK_DATA, &
41 REFKDT_DATA,FRZK_DATA,ZBOT_DATA, SMLOW_DATA,SMHIGH_DATA, &
44 CHARACTER*256 :: err_message
49 !-----------------------------------------------------------------
50 SUBROUTINE LSMRUC(spp_lsm, &
52 pattern_spp_lsm,field_sf, &
57 graupelncv,snowncv,rainncv, &
59 ZS,RAINBL,SNOW,SNOWH,SNOWC,FRZFRAC,frpcpn, &
60 rhosnf,precipfr, & ! pass it out to module_diagnostics
61 Z3D,P8W,T3D,QV3D,QC3D,RHO3D, & !p8W in [PA]
62 GLW,GSW,EMISS,CHKLOWQ, CHS, &
63 FLQC,FLHC,MAVAIL,CANWAT,VEGFRA,ALB,ZNT, &
64 Z0,SNOALB,ALBBCK,LAI, & !new
65 mminlu, landusef, nlcat, mosaic_lu, &
66 mosaic_soil, soilctop, nscat, & !new
67 QSFC,QSG,QVG,QCG,DEW,SOILT1,TSNAV, &
68 TBOT,IVGTYP,ISLTYP,XLAND, &
69 ISWATER,ISICE,XICE,XICE_THRESHOLD, &
70 CP,ROVCP,G0,LV,STBOLT, &
71 SOILMOIS,SH2O,SMAVAIL,SMMAX, &
72 TSO,SOILT,HFX,QFX,LH, &
73 SFCRUNOFF,UDRUNOFF,ACRUNOFF,SFCEXC, &
74 SFCEVP,GRDFLX,SNOWFALLAC,ACSNOW,SNOM, &
75 SMFR3D,KEEPFR3DFLAG, &
76 myjpbl,shdmin,shdmax,rdlai2d, &
77 ids,ide, jds,jde, kds,kde, &
78 ims,ime, jms,jme, kms,kme, &
79 its,ite, jts,jte, kts,kte )
80 !-----------------------------------------------------------------
82 !-----------------------------------------------------------------
84 ! The RUC LSM model is described in:
85 ! Smirnova, T.G., J.M. Brown, and S.G. Benjamin, 1997:
86 ! Performance of different soil model configurations in simulating
87 ! ground surface temperature and surface fluxes.
88 ! Mon. Wea. Rev. 125, 1870-1884.
89 ! Smirnova, T.G., J.M. Brown, and D. Kim, 2000: Parameterization of
90 ! cold-season processes in the MAPS land-surface scheme.
91 ! J. Geophys. Res. 105, 4077-4086.
92 !-----------------------------------------------------------------
93 !-- DT time step (second)
94 ! ktau - number of time step
95 ! NSL - number of soil layers
96 ! NZS - number of levels in soil
97 ! ZS - depth of soil levels (m)
98 !-- RAINBL - accumulated rain in [mm] between the PBL calls
99 !-- RAINNCV one time step grid scale precipitation (mm/step)
100 ! SNOW - snow water equivalent [mm]
101 ! FRAZFRAC - fraction of frozen precipitation
102 !-- PRECIPFR (mm) - time step frozen precipitation
103 !-- SNOWC flag indicating snow coverage (1 for snow cover)
105 !-- P8W 3D pressure (Pa)
106 !-- T3D temperature (K)
107 !-- QV3D 3D water vapor mixing ratio (Kg/Kg)
108 ! QC3D - 3D cloud water mixing ratio (Kg/Kg)
109 ! RHO3D - 3D air density (kg/m^3)
110 !-- GLW downward long wave flux at ground surface (W/m^2)
111 !-- GSW absorbed short wave flux at ground surface (W/m^2)
112 !-- EMISS surface emissivity (between 0 and 1)
113 ! FLQC - surface exchange coefficient for moisture (kg/m^2/s)
114 ! FLHC - surface exchange coefficient for heat [W/m^2/s/degreeK]
115 ! SFCEXC - surface exchange coefficient for heat [m/s]
116 ! CANWAT - CANOPY MOISTURE CONTENT (mm)
117 ! VEGFRA - vegetation fraction (between 0 and 100)
118 ! ALB - surface albedo (between 0 and 1)
119 ! SNOALB - maximum snow albedo (between 0 and 1)
120 ! ALBBCK - snow-free albedo (between 0 and 1)
121 ! ZNT - roughness length [m]
122 !-- TBOT soil temperature at lower boundary (K)
123 ! IVGTYP - USGS vegetation type (24 classes)
124 ! ISLTYP - STASGO soil type (16 classes)
125 !-- XLAND land mask (1 for land, 2 for water)
126 !-- CP heat capacity at constant pressure for dry air (J/kg/K)
127 !-- G0 acceleration due to gravity (m/s^2)
128 !-- LV latent heat of melting (J/kg)
129 !-- STBOLT Stefan-Boltzmann constant (W/m^2/K^4)
130 ! SOILMOIS - soil moisture content (volumetric fraction)
131 ! TSO - soil temp (K)
132 !-- SOILT surface temperature (K)
133 !-- HFX upward heat flux at the surface (W/m^2)
134 !-- QFX upward moisture flux at the surface (kg/m^2/s)
135 !-- LH upward latent heat flux (W/m^2)
136 ! SFCRUNOFF - ground surface runoff [mm]
137 ! UDRUNOFF - underground runoff [mm]
138 ! ACRUNOFF - run-total surface runoff [mm]
139 ! SFCEVP - total evaporation in [kg/m^2]
140 ! GRDFLX - soil heat flux (W/m^2: negative, if downward from surface)
141 ! SNOWFALLAC - run-total snowfall accumulation [m]
142 ! ACSNOW - run-toral SWE of snowfall [mm]
143 !-- CHKLOWQ - is either 0 or 1 (so far set equal to 1).
144 !-- used only in MYJPBL.
145 !-- tice - sea ice temperture (C)
146 !-- rhosice - sea ice density (kg m^-3)
147 !-- capice - sea ice volumetric heat capacity (J/m^3/K)
148 !-- thdifice - sea ice thermal diffusivity (m^2/s)
150 !-- ims start index for i in memory
151 !-- ime end index for i in memory
152 !-- jms start index for j in memory
153 !-- jme end index for j in memory
154 !-- kms start index for k in memory
155 !-- kme end index for k in memory
156 !-------------------------------------------------------------------------
157 ! INTEGER, PARAMETER :: nzss=5
158 ! INTEGER, PARAMETER :: nddzs=2*(nzss-2)
160 INTEGER, PARAMETER :: nvegclas=24+3
162 REAL, INTENT(IN ) :: DT
163 LOGICAL, INTENT(IN ) :: myjpbl,frpcpn
164 INTEGER, INTENT(IN ) :: spp_lsm
165 INTEGER, INTENT(IN ) :: NLCAT, NSCAT, mosaic_lu, mosaic_soil
166 INTEGER, INTENT(IN ) :: ktau, nsl, isice, iswater, &
167 ims,ime, jms,jme, kms,kme, &
168 ids,ide, jds,jde, kds,kde, &
169 its,ite, jts,jte, kts,kte
172 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),OPTIONAL:: pattern_spp_lsm
173 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),OPTIONAL:: field_sf
175 REAL, DIMENSION( ims:ime, 1 :nsl, jms:jme ) :: field_sf_loc
177 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
178 INTENT(IN ) :: QV3D, &
185 REAL, DIMENSION( ims:ime , jms:jme ), &
186 INTENT(IN ) :: RAINBL, &
200 REAL, DIMENSION( ims:ime , jms:jme ), &
201 INTENT(INOUT ) :: VEGFRA
205 REAL, OPTIONAL, DIMENSION( ims:ime , jms:jme ), &
206 INTENT(IN ) :: GRAUPELNCV, &
209 REAL, DIMENSION( ims:ime , jms:jme ), &
210 INTENT(IN ) :: lakemask
211 INTEGER, INTENT(IN ) :: LakeModel
214 REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN ):: SHDMAX
215 REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN ):: SHDMIN
216 LOGICAL, intent(in) :: rdlai2d
218 REAL, DIMENSION( 1:nsl), INTENT(IN ) :: ZS
220 REAL, DIMENSION( ims:ime , jms:jme ), &
235 REAL, DIMENSION( ims:ime , jms:jme ), &
239 INTEGER, DIMENSION( ims:ime , jms:jme ), &
240 INTENT(IN ) :: IVGTYP, &
242 CHARACTER(LEN=*), INTENT(IN ) :: MMINLU
243 REAL, DIMENSION( ims:ime , 1:nlcat, jms:jme ), INTENT(IN):: LANDUSEF
244 REAL, DIMENSION( ims:ime , 1:nscat, jms:jme ), INTENT(IN):: SOILCTOP
246 REAL, INTENT(IN ) :: CP,ROVCP,G0,LV,STBOLT,XICE_threshold
248 REAL, DIMENSION( ims:ime , 1:nsl, jms:jme ) , &
249 INTENT(INOUT) :: SOILMOIS,SH2O,TSO
251 REAL, DIMENSION( ims:ime, jms:jme ) , &
252 INTENT(INOUT) :: SOILT, &
272 REAL, DIMENSION( ims:ime, jms:jme ) , &
273 INTENT(INOUT) :: SMAVAIL, &
276 REAL, DIMENSION( its:ite, jts:jte ) :: &
296 ! Energy and water budget variables:
297 REAL, DIMENSION( its:ite, jts:jte ) :: &
307 REAL, DIMENSION( ims:ime, 1:nsl, jms:jme) &
311 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: &
312 RHOSNF, & !RHO of snowfall
313 PRECIPFR, & ! time-step frozen precip
315 !--- soil/snow properties
343 REAL, DIMENSION(1:NSL) :: ZSMAIN, &
347 REAL, DIMENSION(1:2*(nsl-2)) :: DTDZS
349 REAL, DIMENSION(1:5001) :: TBQ
352 REAL, DIMENSION( 1:nsl ) :: SOILM1D, &
358 REAL, DIMENSION( 1:nsl ) :: KEEPFR
360 REAL, DIMENSION( 1:nlcat ) :: lufrac
361 REAL, DIMENSION( 1:nscat ) :: soilfrac
389 REAL :: cq,r61,r273,arp,brp,x,evs,eis
392 REAL :: meltfactor, ac,as, wb
394 INTEGER :: ILAND,ISOIL,IFOREST
396 INTEGER :: I,J,K,NZS,NZS1,NDDZS
397 INTEGER :: k1,l,k2,kp,km
398 CHARACTER (LEN=132) :: message
400 REAL,DIMENSION(ims:ime,1:nsl,jms:jme) :: rstoch
402 REAL,DIMENSION(ims:ime,jms:jme)::EMISSO,VEGFRAO,ALBO,SNOALBO
403 REAL,DIMENSION(its:ite,jts:jte)::EMISSLO
405 !-----------------------------------------------------------------
417 rstoch(i,k,j) = pattern_spp_lsm(i,k,j)
418 field_sf_loc(i,k,j)=field_sf(i,k,j)
424 !---- table TBQ is for resolution of balance equation in VILKA
428 ARP=77455.*41.9/461.525
433 ! TBQ(K)=R61*EXP(ARP*(R273-1./CQ)-BRP*LOG(CQ*R273))
434 EVS=EXP(17.67*(CQ-273.15)/(CQ-29.65))
435 EIS=EXP(22.514-6.15E3/CQ)
436 if(CQ.ge.273.15) then
445 !--- Initialize soil/vegetation parameters
446 !--- This is temporary until SI is added to mass coordinate ---!!!!!
448 #if ( NMM_CORE == 1 )
456 keepfr3dflag(i,k,j)=0.
458 !--- initializing snow fraction, thereshold = 32 mm of snow water
459 ! or ~100 mm of snow height
461 snowc(i,j) = min(1.,snow(i,j)/32.)
462 if(snow(i,j).le.32.) soilt1(i,j)=tso(i,1,j)
463 !--- initializing inside snow temp if it is not defined
464 IF((soilt1(i,j) .LT. 170.) .or. (soilt1(i,j) .GT.400.)) THEN
465 IF(snow(i,j).gt.32.) THEN
466 soilt1(i,j)=0.5*(soilt(i,j)+tso(i,1,j))
467 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
468 WRITE ( message , FMT='(A,F8.3,2I6)' ) &
469 'Temperature inside snow is initialized in RUCLSM ', soilt1(i,j),i,j
470 CALL wrf_debug ( 0 , message )
473 soilt1(i,j) = tso(i,1,j)
476 tsnav(i,j) =0.5*(soilt(i,j)+tso(i,1,j))-273.15
478 patmb=P8w(i,kms,j)*1.e-2
479 QSG (i,j) = QSN(SOILT(i,j),TBQ)/PATMB
480 IF((qvg(i,j) .LE. 0.) .or. (qvg(i,j) .GT.0.1)) THEN
481 qvg (i,j) = QSG(i,j)*mavail(i,j)
482 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
483 WRITE ( message , FMT='(A,3F8.3,2I6)' ) &
484 'QVG is initialized in RUCLSM ', qvg(i,j),mavail(i,j),qsg(i,j),i,j
485 CALL wrf_debug ( 0 , message )
488 qsfc(i,j) = qvg(i,j)/(1.+qvg(i,j))
493 RHOSNF(i,j) = -1.e3 ! non-zero flag
506 waterbudget(i,j) = 0.
507 acwaterbudget(i,j) = 0.
513 ! For RUC LSM CHKLOWQ needed for MYJPBL should
514 ! 1 because is actual specific humidity at the surface, and
515 ! not the saturation value
536 !-----------------------------------------------------------------
550 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
551 print *,' IN LSMRUC ','ims,ime,jms,jme,its,ite,jts,jte,nzs', &
552 ims,ime,jms,jme,its,ite,jts,jte,nzs
553 print *,' IVGTYP, ISLTYP ', ivgtyp(i,j),isltyp(i,j)
554 print *,' MAVAIL ', mavail(i,j)
555 print *,' SOILT,QVG,P8w',soilt(i,j),qvg(i,j),p8w(i,1,j)
556 print *, 'LSMRUC, I,J,xland, QFX,HFX from SFCLAY',i,j,xland(i,j), &
558 print *, ' GSW, GLW =',gsw(i,j),glw(i,j)
559 print *, 'SOILT, TSO start of time step =',soilt(i,j),(tso(i,k,j),k=1,nsl)
560 print *, 'SOILMOIS start of time step =',(soilmois(i,k,j),k=1,nsl)
561 print *, 'SMFROZEN start of time step =',(smfr3d(i,k,j),k=1,nsl)
562 print *, ' I,J=, after SFCLAY CHS,FLHC ',i,j,chs(i,j),flhc(i,j)
563 print *, 'LSMRUC, IVGTYP,ISLTYP,ALB = ', ivgtyp(i,j),isltyp(i,j),alb(i,j),i,j
564 print *, 'LSMRUC I,J,DT,RAINBL =',I,J,dt,RAINBL(i,j)
565 print *, 'XLAND ---->, ivgtype,isoiltyp,i,j',xland(i,j),ivgtyp(i,j),isltyp(i,j),i,j
572 QVATM = QV3D(i,kms,j)
573 QCATM = QC3D(i,kms,j)
574 PATM = P8w(i,kms,j)*1.e-5
575 !-- Z3D(1) is thickness between first full sigma level and the surface,
576 !-- but first mass level is at the half of the first sigma level
577 !-- (u and v are also at the half of first sigma level)
578 CONFLX = Z3D(i,kms,j)*0.5
580 ! -- initialize snow, graupel and ice fractions in frozen precip
587 prcpncliq = rainncv(i,j)*(1.-frzfrac(i,j))
588 prcpncfr = rainncv(i,j)*frzfrac(i,j)
589 !- apply the same frozen precipitation fraction to convective precip
590 !tgs - 31 mar17 - add safety temperature check in case Thompson MP produces
591 ! frozen precip at T > 273.
592 if(frzfrac(i,j) > 0..and. tabs < 273.) then
593 prcpculiq = max(0.,(rainbl(i,j)-rainncv(i,j))*(1.-frzfrac(i,j)))
594 prcpcufr = max(0.,(rainbl(i,j)-rainncv(i,j))*frzfrac(i,j))
597 prcpcufr = max(0.,(rainbl(i,j)-rainncv(i,j)))
601 prcpculiq = max(0.,(rainbl(i,j)-rainncv(i,j)))
604 !--- 1*e-3 is to convert from mm/s to m/s
605 PRCPMS = (prcpncliq + prcpculiq)/DT*1.e-3
606 NEWSNMS = (prcpncfr + prcpcufr)/DT*1.e-3
608 IF ( PRESENT( graupelncv ) ) THEN
609 graupamt = graupelncv(i,j)
614 if((prcpncfr + prcpcufr) > 0.) then
615 ! -- calculate snow, graupel and ice fractions in falling frozen precip
616 snowrat=min(1.,max(0.,snowncv(i,j)/(prcpncfr + prcpcufr)))
617 grauprat=min(1.,max(0.,graupamt/(prcpncfr + prcpcufr)))
618 icerat=min(1.,max(0.,(prcpncfr-snowncv(i,j)-graupamt) &
619 /(prcpncfr + prcpcufr)))
620 curat=min(1.,max(0.,(prcpcufr/(prcpncfr + prcpcufr))))
623 PRCPMS = (RAINBL(i,j)/DT*1.e-3)*(1-FRZFRAC(I,J))
624 NEWSNMS = (RAINBL(i,j)/DT*1.e-3)*FRZFRAC(I,J)
625 if(newsnms == 0.) then
628 snowrat = min(1.,newsnms/(newsnms+prcpms))
633 if (tabs.le.273.15) then
635 NEWSNMS = RAINBL(i,j)/DT*1.e-3
636 !-- here no info about constituents of frozen precipitation,
637 !-- suppose it is all snow
640 PRCPMS = RAINBL(i,j)/DT*1.e-3
645 ! -- save time-step water equivalent of frozen precipitation in PRECIPFR array to be used in
647 precipfr(i,j) = NEWSNMS * DT *1.e3
649 !--- convert exchange coeff QKMS to [m/s]
650 QKMS=FLQC(I,J)/RHO/MAVAIL(I,J)
651 ! TKMS=FLHC(I,J)/RHO/CP
652 TKMS=FLHC(I,J)/RHO/(CP*(1.+0.84*QVATM)) ! mynnsfc uses CPM
653 !--- convert incoming snow and canwat from mm to m
656 CANWATR=CANWAT(I,J)*1.E-3
659 RHOSNFALL=RHOSNF(I,J)
667 zshalf(k)=0.5*(zsmain(k-1) + zsmain(k))
671 lufrac(k) = landusef(i,k,j)
674 soilfrac(k) = soilctop(i,k,j)
677 !------------------------------------------------------------
678 !----- DDZS and DSDZ1 are for implicit solution of soil eqns.
679 !-------------------------------------------------------------
682 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
683 print *,' DT,NZS1, ZSMAIN, ZSHALF --->', dt,nzs1,zsmain,zshalf
689 X=DT/2./(ZSHALF(K+1)-ZSHALF(K))
690 DTDZS(K1)=X/(ZSMAIN(K)-ZSMAIN(K-1))
692 DTDZS(K2)=X/(ZSMAIN(K+1)-ZSMAIN(K))
695 !27jul2011 - CN and SAT are defined in VEGPARM.TBL
697 ! SAT=0.0004 ! canopy water saturated
702 !--- Constants used in Johansen soil thermal
703 !--- conductivity method
709 !***********************************************************************
710 !--- Constants for snow density calculations C1SN and C2SN
716 !***********************************************************************
722 if(SNOW(i,j).gt.0. .and. SNOWH(i,j).gt.0.) then
723 RHOSN = SNOW(i,j)/SNOWH(i,j)
728 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
729 if(ktau.eq.1 .and.(i.eq.358.and.j.eq.260)) &
730 print *,'before SOILVEGIN - z0,znt(195,254)',z0(i,j),znt(i,j)
732 !--- initializing soil and surface properties
733 CALL SOILVEGIN ( mosaic_lu, mosaic_soil,soilfrac,nscat,shdmin(i,j),shdmax(i,j),&
734 NLCAT,ILAND,ISOIL,iswater,IFOREST,lufrac,VEGFRA(I,J), &
735 EMISSL(I,J),PC(I,J),ZNT(I,J),LAI(I,J),RDLAI2D, &
736 QWRTZ,RHOCS,BCLH,DQM,KSAT,PSIS,QMIN,REF,WILT,i,j )
737 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
738 if(ktau.eq.1 .and.(i.eq.358.and.j.eq.260)) &
739 print *,'after SOILVEGIN - z0,znt(375,254),lai(375,254)',z0(i,j),znt(i,j),lai(i,j)
741 if(ktau.eq.1 .and. (i.eq.358.and.j.eq.260)) then
742 print *,'NLCAT,iland,lufrac,EMISSL(I,J),PC(I,J),ZNT(I,J),LAI(I,J)', &
743 NLCAT,iland,lufrac,EMISSL(I,J),PC(I,J),ZNT(I,J),LAI(I,J),i,j
744 print *,'NSCAT,soilfrac,QWRTZ,RHOCS,BCLH,DQM,KSAT,PSIS,QMIN,REF,WILT',&
745 NSCAT,soilfrac,QWRTZ,RHOCS,BCLH,DQM,KSAT,PSIS,QMIN,REF,WILT,i,j
749 CN=CFACTR_DATA ! exponent
750 ! SAT=max(1.e-5,(min(5.e-4,(CMCMAX_DATA * (1.-exp(-0.5*lai(i,j))) * 0.01*VEGFRA(I,J))))) ! canopy water saturated
751 SAT = 5.e-4 ! units [m]
752 ! if(i==666.and.j==282) print *,'second 666,282 - sat',sat
754 !-- definition of number of soil levels in the rooting zone
755 ! IF(iforest(ivgtyp(i,j)).ne.1) THEN
756 IF(iforest.gt.2) THEN
757 !---- all vegetation types except evergreen and mixed forests
758 !18apr08 - define meltfactor for Egglston melting limit:
759 ! for open areas factor is 2, and for forests - factor is 0.85
760 ! This will make limit on snow melting smaller and let snow stay
761 ! longer in the forests.
765 if(zsmain(k).ge.0.4) then
771 !---- evergreen and mixed forests
772 !18apr08 - define meltfactor
774 ! 28 March 11 - Previously used value of metfactor= 1.5 needs to be further reduced
775 ! to compensate for low snow albedos in the forested areas.
776 ! Melting rate in forests will reduce.
780 if(zsmain(k).ge.1.1) then
789 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
790 print *,' ZNT, LAI, VEGFRA, SAT, EMIS, PC --->', &
791 ZNT(I,J),LAI(I,J),VEGFRA(I,J),SAT,EMISSL(I,J),PC(I,J)
792 print *,' ZS, ZSMAIN, ZSHALF, CONFLX, CN, SAT, --->', zs,zsmain,zshalf,conflx,cn,sat
793 print *,'NROOT, meltfactor, iforest, ivgtyp, i,j ', nroot,meltfactor,iforest,ivgtyp(I,J),I,J
794 ! print *,'NROOT, iforest, ivgtyp, i,j ', nroot,iforest(ivgtyp(i,j)),ivgtyp(I,J),I,J
797 !!*** SET ZERO-VALUE FOR SOME OUTPUT DIAGNOSTIC ARRAYS
798 ! if(i.eq.397.and.j.eq.562) then
799 ! print *,'RUC LSM - xland(i,j),xice(i,j),snow(i,j)',i,j,xland(i,j),xice(i,j),snow(i,j)
803 if(lakemodel==1. .and. lakemask(i,j)==1.) goto 2999
807 IF((XLAND(I,J)-1.5).GE.0.)THEN
819 patmb=P8w(i,1,j)*1.e-2
820 qvg (i,j) = QSN(SOILT(i,j),TBQ)/PATMB
821 qsfc(i,j) = qvg(i,j)/(1.+qvg(i,j))
823 Q2SAT=QSN(TABS,TBQ)/PATMB
828 TSO(I,K,J)= SOILT(I,J)
831 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
832 PRINT*,' water point, I=',I, &
833 'J=',J, 'SOILT=', SOILT(i,j)
838 ! LAND POINT OR SEA ICE
839 if(xice(i,j).ge.xice_threshold) then
840 ! if(IVGTYP(i,j).eq.isice) then
846 IF(SEAICE(I,J).GT.0.5)THEN
848 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
849 PRINT*,' sea-ice at water point, I=',I, &
863 patmb=P8w(i,1,j)*1.e-2
864 qvg (i,j) = QSN(SOILT(i,j),TBQ)/PATMB
866 qsfc(i,j) = qvg(i,j)/(1.+qvg(i,j))
872 keepfr3dflag(i,k,j) = 0.
873 tso(i,k,j) = min(271.4,tso(i,k,j))
877 ! Attention!!!! RUC LSM uses soil moisture content minus residual (minimum
878 ! or dry soil moisture content for a given soil type) as a state variable.
881 ! soilm1d - soil moisture content minus residual [m**3/m**3]
882 soilm1d (k) = min(max(0.,soilmois(i,k,j)-qmin),dqm)
883 ! soilm1d (k) = min(max(0.,soilmois(i,k,j)),dqm)
884 tso1d (k) = tso(i,k,j)
885 soiliqw (k) = min(max(0.,sh2o(i,k,j)-qmin),soilm1d(k))
886 soilice (k) =(soilm1d (k) - soiliqw (k))/0.9
890 smfrkeep(k) = smfr3d(i,k,j)
891 keepfr (k) = keepfr3dflag(i,k,j)
894 LMAVAIL(I,J)=max(0.00001,min(1.,soilm1d(1)/(REF-QMIN)))
895 ! LMAVAIL(I,J)=max(0.00001,min(1.,soilm1d(1)/dqm))
897 #if ( NMM_CORE == 1 )
903 ! extract dew from the cloud water at the surface
904 !30july13 QCG(I,J)=QCG(I,J)-DEW(I,J)/QKMS
907 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
908 print *,'LAND, i,j,tso1d,soilm1d,PATM,TABS,QVATM,QCATM,RHO', &
909 i,j,tso1d,soilm1d,PATM,TABS,QVATM,QCATM,RHO
910 print *,'CONFLX =',CONFLX
911 print *,'SMFRKEEP,KEEPFR ',SMFRKEEP,KEEPFR
916 smtotold(i,j)=smtotold(i,j)+(qmin+soilm1d(k))* &
917 (zshalf(k+1)-zshalf(k))
920 smtotold(i,j)=smtotold(i,j)+(qmin+soilm1d(nzs))* &
921 (zsmain(nzs)-zshalf(nzs))
923 canwatold(i,j) = canwatr
924 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
925 print *,'before SFCTMP, spp_lsm, rstoch, field_sf_loc', &
926 i,j,spp_lsm,(rstoch(i,k,j),k=1,nzs),(field_sf_loc(i,k,j),k=1,nzs)
928 !-----------------------------------------------------------------
929 CALL SFCTMP (spp_lsm,rstoch(i,:,j),field_sf_loc(i,:,j), &
930 dt,ktau,conflx,i,j, &
932 nzs,nddzs,nroot,meltfactor, & !added meltfactor
933 iland,isoil,xland(i,j),ivgtyp(i,j),isltyp(i,j), &
934 PRCPMS, NEWSNMS,SNWE,SNHEI,SNOWFRAC, &
935 RHOSN,RHONEWSN,RHOSNFALL, &
936 snowrat,grauprat,icerat,curat, &
937 PATM,TABS,QVATM,QCATM,RHO, &
938 GLW(I,J),GSW(I,J),EMISSL(I,J), &
939 QKMS,TKMS,PC(I,J),LMAVAIL(I,J), &
940 canwatr,vegfra(I,J),alb(I,J),znt(I,J), &
941 snoalb(i,j),albbck(i,j),lai(i,j), & !new
942 myjpbl,seaice(i,j),isice, &
943 !--- soil fixed fields
945 rhocs,dqm,qmin,ref, &
946 wilt,psis,bclh,ksat, &
947 sat,cn,zsmain,zshalf,DTDZS,DTDZS2,tbq, &
949 cp,rovcp,g0,lv,stbolt,cw,c1sn,c2sn, &
951 !--- output variables
952 snweprint,snheiprint,rsm, &
953 soilm1d,tso1d,smfrkeep,keepfr, &
954 soilt(I,J),soilt1(i,j),tsnav(i,j),dew(I,J), &
955 qvg(I,J),qsg(I,J),qcg(I,J),SMELT(I,J), &
956 SNOH(I,J),SNFLX(I,J),SNOM(I,J),SNOWFALLAC(I,J), &
957 ACSNOW(I,J),edir(I,J),ec(I,J),ett(I,J),qfx(I,J), &
958 lh(I,J),hfx(I,J),sflx(I,J),sublim(I,J), &
959 evapl(I,J),prcpl(I,J),budget(i,j),runoff1(i,j), &
960 runoff2(I,J),soilice,soiliqw,infiltrp,smf(i,j))
961 !-----------------------------------------------------------------
963 ! Fraction of cropland category in the grid box should not have soil moisture below
964 ! wilting point during the growing season.
965 ! Let's keep soil moisture 20% above wilting point for the fraction of grid box under
967 ! This change violates LSM moisture budget, but
968 ! can be considered as a compensation for irrigation not included into LSM.
970 IF (lufrac(crop) > 0 .and. lai(i,j) > 1.1) THEN
971 ! IF (ivgtyp(i,j) == crop .and. lai(i,j) > 1.1) THEN
974 cropsm=1.1*wilt - qmin
975 if(soilm1d(k) < cropsm*lufrac(crop)) then
976 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
977 print * ,'Soil moisture is below wilting in cropland category at time step',ktau &
978 ,'i,j,lufrac(crop),k,soilm1d(k),wilt,cropsm', &
979 i,j,lufrac(crop),k,soilm1d(k),wilt,cropsm
981 soilm1d(k) = cropsm*lufrac(crop)
982 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
983 print * ,'Added soil water to cropland category, i,j,k,soilm1d(k)',i,j,k,soilm1d(k)
988 ELSEIF (ivgtyp(i,j) == natural .and. lai(i,j) > 0.7) THEN
989 ! grassland: assume that 40% of grassland is irrigated cropland
991 cropsm=1.2*wilt - qmin
992 if(soilm1d(k) < cropsm*lufrac(natural)*0.4) then
993 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
994 print * ,'Soil moisture is below wilting in mixed grassland/cropland category at time step',ktau &
995 ,'i,j,lufrac(natural),k,soilm1d(k),wilt', &
996 i,j,lufrac(natural),k,soilm1d(k),wilt
998 soilm1d(k) = cropsm * lufrac(natural)*0.4
999 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
1000 print * ,'Added soil water to grassland category, i,j,k,soilm1d(k)',i,j,k,soilm1d(k)
1006 ! Fill in field_sf to pass perturbed field of hydraulic cond. up to model driver and output
1008 if (spp_lsm==1) then
1010 field_sf(i,k,j)=field_sf_loc(i,k,j)
1016 !--- available and maximum soil moisture content in the soil
1023 smavail(i,j)=smavail(i,j)+(qmin+soilm1d(k))* &
1024 (zshalf(k+1)-zshalf(k))
1025 smmax (i,j) =smmax (i,j)+(qmin+dqm)* &
1026 (zshalf(k+1)-zshalf(k))
1029 smavail(i,j)=smavail(i,j)+(qmin+soilm1d(nzs))* &
1030 (zsmain(nzs)-zshalf(nzs))
1031 smmax (i,j) =smmax (i,j)+(qmin+dqm)* &
1032 (zsmain(nzs)-zshalf(nzs))
1034 !--- Convert the water unit into mm
1035 SFCRUNOFF(I,J) = SFCRUNOFF(I,J)+RUNOFF1(I,J)*DT*1000.0
1036 UDRUNOFF (I,J) = UDRUNOFF(I,J)+RUNOFF2(I,J)*DT*1000.0
1037 ACRUNOFF(I,J) = ACRUNOFF(I,J)+RUNOFF1(I,J)*DT*1000.0
1038 SMAVAIL (I,J) = SMAVAIL(I,J) * 1000.
1039 SMMAX (I,J) = SMMAX(I,J) * 1000.
1040 smtotold (I,J) = smtotold(I,J) * 1000.
1044 ! soilmois(i,k,j) = soilm1d(k)
1045 soilmois(i,k,j) = soilm1d(k) + qmin
1046 sh2o (i,k,j) = min(soiliqw(k) + qmin,soilmois(i,k,j))
1047 tso(i,k,j) = tso1d(k)
1050 tso(i,nzs,j) = tbot(i,j)
1053 smfr3d(i,k,j) = smfrkeep(k)
1054 keepfr3dflag(i,k,j) = keepfr (k)
1057 !tgs add together dew and cloud at the ground surface
1058 !30july13 qcg(i,j)=qcg(i,j)+dew(i,j)/qkms
1060 Z0 (I,J) = ZNT (I,J)
1062 patmb=P8w(i,1,j)*1.e-2
1063 Q2SAT=QSN(TABS,TBQ)/PATMB
1064 QSFC(I,J) = QVG(I,J)/(1.+QVG(I,J))
1065 ! for MYJ PBL scheme
1066 IF((myjpbl).AND.(QVATM.GE.Q2SAT*0.95).AND.QVATM.LT.qvg(I,J))THEN
1072 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
1073 if(CHKLOWQ(I,J).eq.0.) then
1074 print *,'i,j,CHKLOWQ', &
1079 if(snow(i,j)==0.) EMISSL(i,j) = LEMITBL(IVGTYP(i,j))
1080 EMISS (I,J) = EMISSL(I,J)
1081 ! SNOW is in [mm], SNWE is in [m]; CANWAT is in mm, CANWATR is in m
1082 SNOW (i,j) = SNWE*1000.
1084 CANWAT (I,J) = CANWATR*1000.
1086 INFILTR(I,J) = INFILTRP
1088 MAVAIL (i,j) = LMAVAIL(I,J)
1089 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
1090 print *,' LAND, I=,J=, QFX, HFX after SFCTMP', i,j,lh(i,j),hfx(i,j)
1092 !!! QFX (I,J) = LH(I,J)/LV
1093 SFCEVP (I,J) = SFCEVP (I,J) + QFX (I,J) * DT
1094 GRDFLX (I,J) = -1. * sflx(I,J)
1096 ! if(smf(i,j) .ne.0.) then
1097 !tgs - SMF.NE.0. when there is phase change in the top soil layer
1098 ! The heat of soil water freezing/thawing is not computed explicitly
1099 ! and is responsible for the residual in the energy budget.
1100 ! print *,'Budget',budget(i,j),i,j,smf(i,j)
1103 !--- SNOWC snow cover flag
1104 if(snowfrac > 0. .and. xice(i,j).ge.xice_threshold ) then
1105 SNOWFRAC = SNOWFRAC*XICE(I,J)
1110 !--- RHOSNF - density of snowfall
1111 RHOSNF(I,J)=RHOSNFALL
1113 ! Accumulated moisture flux [kg/m^2]
1114 SFCEVP (I,J) = SFCEVP (I,J) + QFX (I,J) * DT
1116 !TEST!!!! for test put heat budget term in GRDFLX
1118 ! acbudget(i,j)=acbudget(i,j)+budget(i,j)-smf(i,j)
1119 ! GRDFLX (I,J) = acbudget(i,j)
1121 ! if(smf(i,j) .ne.0.) then
1122 !tgs - SMF.NE.0. when there is phase change in the top soil layer
1123 ! The heat of freezing/thawing of soil water is not computed explicitly
1124 ! and is responsible for the residual in the energy budget.
1126 ! budget(i,j)=budget(i,j)-smf(i,j)
1131 ac=max(0.,canwat(i,j)-canwatold(i,j))
1132 as=max(0.,snwe-snowold(i,j))
1133 wb =rainbl(i,j)+smelt(i,j)*dt*1.e3 & ! source
1135 -runoff1(i,j)*dt*1.e3-runoff2(i,j)*dt*1.e3 &
1136 -ac-as - (smavail(i,j)-smtotold(i,j))
1138 waterbudget(i,j)=rainbl(i,j)+smelt(i,j)*dt*1.e3 & ! source
1140 -runoff1(i,j)*dt*1.e3-runoff2(i,j)*dt*1.e3 &
1141 -ac-as - (smavail(i,j)-smtotold(i,j))
1144 ! waterbudget(i,j)=rainbl(i,j)-qfx(i,j)*dt-(smavail(i,j)-smtotold(i,j)) &
1145 acwaterbudget(i,j)=acwaterbudget(i,j)+waterbudget(i,j)
1147 !!!!TEST use LH to check water budget
1148 ! GRDFLX (I,J) = waterbudget(i,j)
1150 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
1151 print *,'Smf=',smf(i,j),i,j
1152 print *,'Budget',budget(i,j),i,j
1153 print *,'RUNOFF2= ', i,j,runoff2(i,j)
1154 print *,'Water budget ', i,j,waterbudget(i,j)
1155 print *,'rainbl,qfx*dt,runoff1,smelt*dt*1.e3,smchange', &
1156 i,j,rainbl(i,j),qfx(i,j)*dt,runoff1(i,j)*dt*1.e3, &
1157 smelt(i,j)*dt*1.e3, &
1158 (smavail(i,j)-smtotold(i,j))
1160 print *,'SNOW,SNOWold',i,j,snwe,snowold(i,j)
1161 print *,'SNOW-SNOWold',i,j,max(0.,snwe-snowold(i,j))
1162 print *,'CANWATold, canwat ',i,j,canwatold(i,j),canwat(i,j)
1163 print *,'canwat(i,j)-canwatold(i,j)',max(0.,canwat(i,j)-canwatold(i,j))
1167 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
1168 print *,'LAND, i,j,tso1d,soilm1d,soilt - end of time step', &
1169 i,j,tso1d,soilm1d,soilt(i,j)
1170 print *,'LAND, QFX, HFX after SFCTMP', i,j,lh(i,j),hfx(i,j)
1173 !--- end of a land or sea ice point
1175 2999 continue ! lakes
1181 !-----------------------------------------------------------------
1182 END SUBROUTINE LSMRUC
1183 !-----------------------------------------------------------------
1187 SUBROUTINE SFCTMP (spp_lsm,rstochcol,fieldcol_sf, &
1188 delt,ktau,conflx,i,j, &
1189 !--- input variables
1190 nzs,nddzs,nroot,meltfactor, &
1191 ILAND,ISOIL,XLAND,IVGTYP,ISLTYP,PRCPMS, &
1192 NEWSNMS,SNWE,SNHEI,SNOWFRAC, &
1193 RHOSN,RHONEWSN,RHOSNFALL, &
1194 snowrat,grauprat,icerat,curat, &
1195 PATM,TABS,QVATM,QCATM,rho, &
1196 GLW,GSW,EMISS,QKMS,TKMS,PC, &
1197 MAVAIL,CST,VEGFRA,ALB,ZNT, &
1198 ALB_SNOW,ALB_SNOW_FREE,lai, &
1200 !--- soil fixed fields
1201 QWRTZ,rhocs,dqm,qmin,ref,wilt,psis,bclh,ksat, &
1202 sat,cn,zsmain,zshalf,DTDZS,DTDZS2,tbq, &
1204 cp,rovcp,g0,lv,stbolt,cw,c1sn,c2sn, &
1206 !--- output variables
1207 snweprint,snheiprint,rsm, &
1208 soilm1d,ts1d,smfrkeep,keepfr,soilt,soilt1, &
1209 tsnav,dew,qvg,qsg,qcg, &
1210 SMELT,SNOH,SNFLX,SNOM,SNOWFALLAC,ACSNOW, &
1211 edir1,ec1,ett1,eeta,qfx,hfx,s,sublim, &
1212 evapl,prcpl,fltot,runoff1,runoff2,soilice, &
1213 soiliqw,infiltr,smf)
1214 !-----------------------------------------------------------------
1216 !-----------------------------------------------------------------
1218 !--- input variables
1220 INTEGER, INTENT(IN ) :: isice,i,j,nroot,ktau,nzs , &
1221 nddzs !nddzs=2*(nzs-2)
1223 REAL, INTENT(IN ) :: DELT,CONFLX,meltfactor
1224 REAL, INTENT(IN ) :: C1SN,C2SN
1225 LOGICAL, INTENT(IN ) :: myj
1226 !--- 3-D Atmospheric variables
1228 INTENT(IN ) :: PATM, &
1233 INTENT(IN ) :: GLW, &
1245 INTEGER, INTENT(IN ) :: IVGTYP, ISLTYP
1248 INTENT(INOUT) :: EMISS, &
1255 !--- soil properties
1268 REAL, INTENT(IN ) :: CN, &
1279 REAL, DIMENSION(1:NZS), INTENT(IN) :: ZSMAIN, &
1283 REAL, DIMENSION(1:NZS), INTENT(IN) :: rstochcol
1284 REAL, DIMENSION(1:NZS), INTENT(INOUT) :: fieldcol_sf
1287 REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS
1289 REAL, DIMENSION(1:5001), INTENT(IN) :: TBQ
1292 !--- input/output variables
1293 !-------- 3-d soil moisture and temperature
1294 REAL, DIMENSION( 1:nzs ) , &
1295 INTENT(INOUT) :: TS1D, &
1298 REAL, DIMENSION( 1:nzs ) , &
1299 INTENT(INOUT) :: KEEPFR
1301 REAL, DIMENSION(1:NZS), INTENT(INOUT) :: SOILICE, &
1305 INTEGER, INTENT(INOUT) :: ILAND,ISOIL
1308 !-------- 2-d variables
1310 INTENT(INOUT) :: DEW, &
1349 REAL, DIMENSION(1:NZS) :: &
1360 !-------- 1-d variables
1386 REAL, INTENT(INOUT) :: RSM, &
1389 INTEGER, INTENT(IN) :: spp_lsm
1390 !--- Local variables
1394 REAL :: BSN, XSN , &
1395 RAINF, SNTH, NEWSN, PRCPMS, NEWSNMS , &
1397 REAL :: snhei_crit, snhei_crit_newsn, keep_snow_albedo, SNOWFRACnewsn
1398 REAL :: newsnowratio, dd1
1400 REAL :: rhonewgr,rhonewice
1402 REAL :: RNET,GSWNEW,GSWIN,EMISSN,ZNTSN,EMISS_snowfree
1403 REAL :: VEGFRAC, snow_mosaic, snfr, vgfr
1404 real :: cice, albice, albsn, drip, dripsn, dripliq
1405 real :: interw, intersn, infwater, intwratio
1407 !-----------------------------------------------------------------
1408 integer, parameter :: ilsnow=99
1410 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
1411 print *,' in SFCTMP',i,j,nzs,nddzs,nroot, &
1412 SNWE,RHOSN,SNOM,SMELT,TS1D
1420 if(snhei == 0.) snowfrac=0.
1426 ! Jul 2016 - Avissar and Pielke (1989)
1427 ! This formulation depending on LAI defines relative contribution of the vegetation to
1428 ! the total heat fluxes between surface and atmosphere.
1429 ! With VEGFRA=100% and LAI=3, VEGFRAC=0.86 meaning that vegetation contributes
1430 ! only 86% of the total surface fluxes.
1431 ! VGFR=0.01*VEGFRA ! % --> fraction
1432 ! VEGFRAC=2.*lai*vgfr/(1.+2.*lai*vgfr)
1442 !---initialize local arrays for sea ice
1453 ALBice=ALB_SNOW_FREE
1456 EMISS_snowfree = LEMITBL(IVGTYP)
1458 !--- sea ice properties
1459 !--- N.N Zubov "Arctic Ice"
1460 !--- no salinity dependence because we consider the ice pack
1461 !--- to be old and to have low salinity (0.0002)
1462 if(SEAICE.ge.0.5) then
1464 tice(k) = ts1d(k) - 273.15
1465 rhosice(k) = 917.6/(1-0.000165*tice(k))
1466 cice = 2115.85 +7.7948*tice(k)
1467 capice(k) = cice*rhosice(k)
1468 thdifice(k) = 2.260872/capice(k)
1470 !-- SEA ICE ALB dependence on ice temperature. When ice temperature is
1471 !-- below critical value of -10C - no change to albedo.
1472 !-- If temperature is higher that -10C then albedo is decreasing.
1473 !-- The minimum albedo at t=0C for ice is 0.1 less.
1474 ALBice = MIN(ALB_SNOW_FREE,MAX(ALB_SNOW_FREE - 0.05, &
1475 ALB_SNOW_FREE - 0.1*(tice(1)+10.)/10. ))
1478 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
1479 ! print *,'I,J,KTAU,QKMS,TKMS', i,j,ktau,qkms,tkms
1480 print *,'alb_snow_free',ALB_SNOW_FREE
1481 print *,'GSW,GSWnew,GLW,SOILT,EMISS,ALB,ALBice,SNWE',&
1482 GSW,GSWnew,GLW,SOILT,EMISS,ALB,ALBice,SNWE
1485 if(snhei.gt.0.0081*1.e3/rhosn) then
1486 !*** Update snow density for current temperature (Koren et al. 1999)
1487 BSN=delt/3600.*c1sn*exp(0.08*min(0.,tsnav)-c2sn*rhosn*1.e-3)
1488 if(bsn*snwe*100..lt.1.e-4) goto 777
1489 XSN=rhosn*(exp(bsn*snwe*100.)-1.)/(bsn*snwe*100.)
1490 rhosn=MIN(MAX(58.8,XSN),500.) ! 13mar18 - switch from 76.9 to 58.8
1497 IF(NEWSN.GT.0.) THEN
1498 ! IF(NEWSN.GE.1.E-8) THEN
1500 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
1501 print *, 'THERE IS NEW SNOW, newsn', newsn
1504 newsnowratio = min(1.,newsn/(snwe+newsn))
1506 !--- 27 Feb 2014 - empirical formulations from John M. Brown
1507 ! rhonewsn=min(250.,rhowater/max(4.179,(13.*tanh((274.15-Tabs)*0.3333))))
1508 !--- 13 Mar 2018 - formulation from Trevor Elcott
1509 rhonewsn=min(125.,1000.0/max(8.,(17.*tanh((276.65-Tabs)*0.15))))
1510 rhonewgr=min(500.,rhowater/max(2.,(3.5*tanh((274.15-Tabs)*0.3333))))
1513 !--- compute density of "snowfall" from weighted contribution
1514 ! of snow, graupel and ice fractions
1516 rhosnfall = min(500.,max(58.8,(rhonewsn*snowrat + & ! 13mar18-switch from 76.9 to 58.8
1517 rhonewgr*grauprat + rhonewice*icerat + rhonewgr*curat)))
1519 ! from now on rhonewsn is the density of falling frozen precipitation
1522 !*** Define average snow density of the snow pack considering
1523 !*** the amount of fresh snow (eq. 9 in Koren et al.(1999)
1524 !*** without snow melt )
1525 xsn=(rhosn*snwe+rhonewsn*newsn)/ &
1527 rhosn=MIN(MAX(58.8,XSN),500.) ! 13mar18 - switch from 76.9 to 58.8
1529 ENDIF ! end NEWSN > 0.
1531 IF(PRCPMS.NE.0.) THEN
1533 ! PRCPMS is liquid precipitation rate
1534 ! RAINF is a flag used for calculation of rain water
1535 ! heat content contribution into heat budget equation. Rain's temperature
1536 ! is set equal to air temperature at the first atmospheric
1544 if(vegfrac > 0.01) then
1545 ! compute intercepted precipitation - Eq. 1 Lawrence et al.,
1546 ! J. of Hydrometeorology, 2006, CLM.
1547 interw=0.25*DELT*PRCPMS*(1.-exp(-0.5*lai))*vegfrac
1548 intersn=0.25*NEWSN*(1.-exp(-0.5*lai))*vegfrac
1549 infwater=PRCPMS - interw/delt
1550 if((interw+intersn) > 0.) then
1551 intwratio=interw/(interw+intersn)
1554 ! Update water/snow intercepted by the canopy
1555 dd1=CST + interw + intersn
1567 endif ! vegfrac > 0.01
1569 ! SNHEI_CRIT is a threshold for fractional snow
1570 SNHEI_CRIT=0.01601*1.e3/rhosn
1571 SNHEI_CRIT_newsn=0.0005*1.e3/rhosn
1572 ! snowfrac from the previous time step
1573 SNOWFRAC=MIN(1.,SNHEI/(2.*SNHEI_CRIT))
1574 if(snowfrac < 0.75) snow_mosaic = 1.
1576 IF(NEWSN.GT.0.) THEN
1577 !Update snow on the ground
1578 snwe=max(0.,snwe+newsn-intersn)
1579 ! Add drip to snow on the ground
1581 if (snow_mosaic==1.) then
1582 dripliq=drip*intwratio
1583 dripsn = drip - dripliq
1585 infwater=infwater+dripliq
1592 snhei=snwe*rhowater/rhosn
1593 NEWSN=NEWSN*rhowater/rhonewsn
1596 IF(SNHEI.GT.0.0) THEN
1597 !-- SNOW on the ground
1598 !--- Land-use category should be changed to snow/ice for grid points with snow>0
1600 !24nov15 - based on field exp on Pleasant View soccer fields
1601 ! if(meltfactor > 1.5) then ! all veg. types, except forests
1602 ! SNHEI_CRIT=0.01601*1.e3/rhosn
1603 ! Petzold - 1 cm of fresh snow overwrites effects from old snow.
1604 ! Need to test SNHEI_CRIT_newsn=0.01
1605 ! SNHEI_CRIT_newsn=0.01
1607 ! SNHEI_CRIT=0.02*1.e3/rhosn
1608 ! SNHEI_CRIT_newsn=0.001*1.e3/rhosn
1611 SNOWFRAC=MIN(1.,SNHEI/(2.*SNHEI_CRIT))
1612 !24nov15 - SNOWFRAC for urban category < 0.75
1613 if(ivgtyp == urban) snowfrac=min(0.75,snowfrac)
1614 ! if(meltfactor > 1.5) then
1615 ! if(isltyp > 9 .and. isltyp < 13) then
1616 !24nov15 clay soil types - SNOFRAC < 0.9
1617 ! snowfrac=min(0.9,snowfrac)
1620 !24nov15 - SNOWFRAC for forests < 0.75
1621 ! snowfrac=min(0.85,snowfrac)
1624 ! SNOWFRAC=MIN(1.,SNHEI/(2.*SNHEI_CRIT))
1625 ! elseif(snowfrac < 0.3 .and. tabs > 275.) then
1627 if(snowfrac < 0.75) snow_mosaic = 1.
1629 if(newsn > 0. ) SNOWFRACnewsn=MIN(1.,SNHEI/SNHEI_CRIT_newsn)
1631 KEEP_SNOW_ALBEDO = 0.
1632 IF (NEWSN > 0. .and. snowfracnewsn > 0.99) THEN
1634 KEEP_SNOW_ALBEDO = 1.
1635 snow_mosaic=0. ! ???
1638 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
1639 print *,'SNHEI_CRIT,SNOWFRAC,SNHEI_CRIT_newsn,SNOWFRACnewsn', &
1640 SNHEI_CRIT,SNOWFRAC,SNHEI_CRIT_newsn,SNOWFRACnewsn
1643 !-- Set znt for snow from VEGPARM table (snow/ice landuse), except for
1644 !-- land-use types with higher roughness (forests, urban).
1645 !5mar12 IF(znt.lt.0.2 .and. snowfrac.gt.0.99) znt=z0tbl(iland)
1646 ! IF(newsn==0. .and. znt.lt.0.2 .and. snowfrac.gt.0.99) znt=z0tbl(iland)
1647 IF(newsn.eq.0. .and. znt.le.0.2 .and. IVGTYP.ne.isice) then
1648 if( snhei .le. 2.*ZNT)then
1649 znt=0.55*znt+0.45*z0tbl(iland)
1650 elseif( snhei .gt. 2.*ZNT .and. snhei .le. 4.*ZNT)then
1651 znt=0.2*znt+0.8*z0tbl(iland)
1652 elseif(snhei > 4.*ZNT) then
1658 !--- GSWNEW in-coming solar for snow on land or on ice
1659 ! GSWNEW=GSWnew/(1.-ALB)
1660 !-- Time to update snow and ice albedo
1662 IF(SEAICE .LT. 0.5) THEN
1664 !-- ALB dependence on snow depth
1665 ! ALB_SNOW across Canada's forested areas is very low - 0.27-0.35, this
1666 ! causes significant warm biases. Limiting ALB in these areas to be higher than 0.4
1667 ! hwlps with these biases..
1668 if( snow_mosaic == 1.) then
1670 ! ALBsn=max(0.4,alb_snow)
1673 ALBsn = MAX(keep_snow_albedo*alb_snow, &
1674 MIN((alb_snow_free + &
1675 (alb_snow - alb_snow_free) * snowfrac), alb_snow))
1677 Emiss = MAX(keep_snow_albedo*emissn, &
1678 MIN((emiss_snowfree + &
1679 (emissn - emiss_snowfree) * snowfrac), emissn))
1681 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
1682 ! if(i.eq.279.and.j.eq.263) then
1683 print *,'Snow on soil ALBsn,emiss,snow_mosaic',i,j,ALBsn,emiss,snow_mosaic
1685 !28mar11 if canopy is covered with snow to 95% of its capacity and snow depth is
1686 ! higher than patchy snow treshold - then snow albedo is not less than 0.55
1687 ! (inspired by the flight from Fairbanks to Seatle)
1689 !test if(cst.ge.0.95*sat .and. snowfrac .gt.0.99)then
1690 ! albsn=max(alb_snow,0.55)
1693 !-- ALB dependence on snow temperature. When snow temperature is
1694 !-- below critical value of -10C - no change to albedo.
1695 !-- If temperature is higher that -10C then albedo is decreasing.
1696 !-- The minimum albedo at t=0C for snow on land is 15% less than
1697 !-- albedo of temperatures below -10C.
1698 if(albsn.lt.0.4 .or. keep_snow_albedo==1) then
1700 ! ALB=max(0.4,alb_snow)
1702 !-- change albedo when no fresh snow and snow albedo is higher than 0.5
1703 ALB = MIN(ALBSN,MAX(ALBSN - 0.1*(soilt - 263.15)/ &
1704 (273.15-263.15)*ALBSN, ALBSN - 0.05))
1708 if( snow_mosaic == 1.) then
1712 ALBsn = MAX(keep_snow_albedo*alb_snow, &
1713 MIN((albice + (alb_snow - albice) * snowfrac), alb_snow))
1714 Emiss = MAX(keep_snow_albedo*emissn, &
1715 MIN((emiss_snowfree + &
1716 (emissn - emiss_snowfree) * snowfrac), emissn))
1719 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
1720 print *,'Snow on ice snow_mosaic,ALBsn,emiss',i,j,ALBsn,emiss,snow_mosaic
1722 !-- ALB dependence on snow temperature. When snow temperature is
1723 !-- below critical value of -10C - no change to albedo.
1724 !-- If temperature is higher that -10C then albedo is decreasing.
1725 if(albsn.lt.alb_snow .or. keep_snow_albedo .eq.1.)then
1728 !-- change albedo when no fresh snow
1729 ALB = MIN(ALBSN,MAX(ALBSN - 0.15*ALBSN*(soilt - 263.15)/ &
1730 (273.15-263.15), ALBSN - 0.1))
1735 if (snow_mosaic==1.) then
1736 !may 2014 - treat separately snow-free and snow-covered areas
1738 if(SEAICE .LT. 0.5) then
1740 ! portion not covered with snow
1741 ! compute absorbed GSW for snow-free portion
1743 gswnew=GSWin*(1.-alb_snow_free)
1745 T3 = STBOLT*SOILT*SOILT*SOILT
1747 XINET = EMISS_snowfree*(GLW-UPFLUX)
1748 RNET = GSWnew + XINET
1749 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
1750 ! if(i.eq.442.and.j.eq.260) then
1751 print *,'Fractional snow - snowfrac=',snowfrac
1752 print *,'Snowfrac<1 GSWin,GSWnew -',GSWin,GSWnew,'SOILT, RNET',soilt,rnet
1755 soilm1ds(k) = soilm1d(k)
1757 smfrkeeps(k) = smfrkeep(k)
1758 keepfrs(k) = keepfr(k)
1759 soilices(k) = soilice(k)
1760 soiliqws(k) = soiliqw(k)
1774 CALL SOIL(spp_lsm,rstochcol,fieldcol_sf, &
1775 !--- input variables
1776 i,j,ilands,isoil,delt,ktau,conflx,nzs,nddzs,nroot, &
1777 PRCPMS,RAINF,PATM,QVATM,QCATM,GLW,GSWnew,gswin, &
1778 EMISS_snowfree,RNET,QKMS,TKMS,PC,csts,dripliq, &
1779 infwater,rho,vegfrac,lai,myj, &
1780 !--- soil fixed fields
1781 QWRTZ,rhocs,dqm,qmin,ref,wilt, &
1782 psis,bclh,ksat,sat,cn, &
1783 zsmain,zshalf,DTDZS,DTDZS2,tbq, &
1785 lv,CP,rovcp,G0,cw,stbolt,tabs, &
1787 !--- output variables for snow-free portion
1788 soilm1ds,ts1ds,smfrkeeps,keepfrs, &
1789 dews,soilts,qvgs,qsgs,qcgs,edir1s,ec1s, &
1790 ett1s,eetas,qfxs,hfxs,ss,evapls,prcpls,fltots,runoff1s, &
1791 runoff2s,mavails,soilices,soiliqws, &
1795 ! portion not covered with snow
1796 ! compute absorbed GSW for snow-free portion
1798 gswnew=GSWin*(1.-albice)
1800 T3 = STBOLT*SOILT*SOILT*SOILT
1802 XINET = EMISS_snowfree*(GLW-UPFLUX)
1803 RNET = GSWnew + XINET
1804 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
1805 ! if(i.eq.442.and.j.eq.260) then
1806 print *,'Fractional snow - snowfrac=',snowfrac
1807 print *,'Snowfrac<1 GSWin,GSWnew -',GSWin,GSWnew,'SOILT, RNET',soilt,rnet
1821 !--- input variables
1822 i,j,iland,isoil,delt,ktau,conflx,nzs,nddzs,nroot, &
1823 PRCPMS,RAINF,PATM,QVATM,QCATM,GLW,GSWnew, &
1824 0.98,RNET,QKMS,TKMS,rho,myj, &
1825 !--- sea ice parameters
1826 tice,rhosice,capice,thdifice, &
1827 zsmain,zshalf,DTDZS,DTDZS2,tbq, &
1829 lv,CP,rovcp,cw,stbolt,tabs, &
1830 !--- output variable
1831 ts1ds,dews,soilts,qvgs,qsgs,qcgs, &
1832 eetas,qfxs,hfxs,ss,evapls,prcpls,fltots &
1849 endif ! seaice < 0.5
1851 !return gswnew to incoming solar
1852 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
1853 ! if(i.eq.442.and.j.eq.260) then
1854 print *,'gswnew,alb_snow_free,alb',gswnew,alb_snow_free,alb
1856 ! gswnew=gswnew/(1.-alb_snow_free)
1858 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
1859 ! if(i.eq.442.and.j.eq.260) then
1860 print *,'Incoming GSWnew snowfrac<1 -',gswnew
1862 endif ! snow_mosaic=1.
1864 !--- recompute absorbed solar radiation and net radiation
1865 !--- for updated value of snow albedo - ALB
1866 gswnew=GSWin*(1.-alb)
1867 ! print *,'SNOW fraction GSWnew',gswnew,'alb=',alb
1869 T3 = STBOLT*SOILT*SOILT*SOILT
1871 XINET = EMISS*(GLW-UPFLUX)
1872 RNET = GSWnew + XINET
1873 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
1874 ! if(i.eq.442.and.j.eq.260) then
1875 ! if(i.eq.271.and.j.eq.242) then
1876 print *,'RNET=',rnet
1877 print *,'SNOW - I,J,newsn,snwe,snhei,GSW,GSWnew,GLW,UPFLUX,ALB',&
1878 i,j,newsn,snwe,snhei,GSW,GSWnew,GLW,UPFLUX,ALB
1881 if (SEAICE .LT. 0.5) then
1883 if(snow_mosaic==1.)then
1888 CALL SNOWSOIL (spp_lsm,rstochcol,fieldcol_sf, & !--- input variables
1889 i,j,isoil,delt,ktau,conflx,nzs,nddzs,nroot, &
1890 meltfactor,rhonewsn,SNHEI_CRIT, & ! new
1891 ILAND,PRCPMS,RAINF,NEWSN,snhei,SNWE,snfr, &
1892 RHOSN,PATM,QVATM,QCATM, &
1893 GLW,GSWnew,GSWin,EMISS,RNET,IVGTYP, &
1894 QKMS,TKMS,PC,CST,dripsn,infwater, &
1895 RHO,VEGFRAC,ALB,ZNT,lai, &
1897 !--- soil fixed fields
1898 QWRTZ,rhocs,dqm,qmin,ref,wilt,psis,bclh,ksat, &
1899 sat,cn,zsmain,zshalf,DTDZS,DTDZS2,tbq, &
1901 lv,CP,rovcp,G0,cw,stbolt,tabs, &
1903 !--- output variables
1904 ilnb,snweprint,snheiprint,rsm, &
1905 soilm1d,ts1d,smfrkeep,keepfr, &
1906 dew,soilt,soilt1,tsnav,qvg,qsg,qcg, &
1907 SMELT,SNOH,SNFLX,SNOM,edir1,ec1,ett1,eeta, &
1908 qfx,hfx,s,sublim,prcpl,fltot,runoff1,runoff2, &
1909 mavail,soilice,soiliqw,infiltr )
1912 if(snow_mosaic==1.)then
1919 i,j,isoil,delt,ktau,conflx,nzs,nddzs, &
1920 meltfactor,rhonewsn,SNHEI_CRIT, & ! new
1921 ILAND,PRCPMS,RAINF,NEWSN,snhei,SNWE,snfr, &
1922 RHOSN,PATM,QVATM,QCATM, &
1923 GLW,GSWnew,EMISS,RNET, &
1924 QKMS,TKMS,RHO,myj, &
1925 !--- sea ice parameters
1927 tice,rhosice,capice,thdifice, &
1928 zsmain,zshalf,DTDZS,DTDZS2,tbq, &
1930 lv,CP,rovcp,cw,stbolt,tabs, &
1931 !--- output variables
1932 ilnb,snweprint,snheiprint,rsm,ts1d, &
1933 dew,soilt,soilt1,tsnav,qvg,qsg,qcg, &
1934 SMELT,SNOH,SNFLX,SNOM,eeta, &
1935 qfx,hfx,s,sublim,prcpl,fltot &
1955 if(snhei.eq.0.) then
1956 !--- all snow is melted
1961 if (snow_mosaic==1.) then
1962 ! May 2014 - now combine snow covered and snow-free land fluxes, soil temp, moist,
1964 if(SEAICE .LT. 0.5) then
1966 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
1967 ! if(i.eq.442.and.j.eq.260) then
1968 print *,'SOILT snow on land', ktau, i,j,soilt
1969 print *,'SOILT on snow-free land', i,j,soilts
1970 print *,'ts1d,ts1ds',i,j,ts1d,ts1ds
1971 print *,' SNOW flux',i,j, snflx
1972 print *,' Ground flux on snow-covered land',i,j, s
1973 print *,' Ground flux on snow-free land', i,j,ss
1974 print *,' CSTS, CST', i,j,csts,cst
1977 soilm1d(k) = soilm1ds(k)*(1.-snowfrac) + soilm1d(k)*snowfrac
1978 ts1d(k) = ts1ds(k)*(1.-snowfrac) + ts1d(k)*snowfrac
1979 smfrkeep(k) = smfrkeeps(k)*(1.-snowfrac) + smfrkeep(k)*snowfrac
1980 if(snowfrac > 0.5) then
1981 keepfr(k) = keepfr(k)
1983 keepfr(k) = keepfrs(k)
1985 soilice(k) = soilices(k)*(1.-snowfrac) + soilice(k)*snowfrac
1986 soiliqw(k) = soiliqws(k)*(1.-snowfrac) + soiliqw(k)*snowfrac
1988 dew = dews*(1.-snowfrac) + dew*snowfrac
1989 soilt = soilts*(1.-snowfrac) + soilt*snowfrac
1990 qvg = qvgs*(1.-snowfrac) + qvg*snowfrac
1991 qsg = qsgs*(1.-snowfrac) + qsg*snowfrac
1992 qcg = qcgs*(1.-snowfrac) + qcg*snowfrac
1993 edir1 = edir1s*(1.-snowfrac) + edir1*snowfrac
1994 ec1 = ec1s*(1.-snowfrac) + ec1*snowfrac
1995 cst = csts*(1.-snowfrac) + cst*snowfrac
1996 ett1 = ett1s*(1.-snowfrac) + ett1*snowfrac
1997 eeta = eetas*(1.-snowfrac) + eeta*snowfrac
1998 qfx = qfxs*(1.-snowfrac) + qfx*snowfrac
1999 hfx = hfxs*(1.-snowfrac) + hfx*snowfrac
2000 s = ss*(1.-snowfrac) + s*snowfrac
2001 evapl = evapls*(1.-snowfrac)
2002 sublim = sublim*snowfrac
2003 prcpl = prcpls*(1.-snowfrac) + prcpl*snowfrac
2004 fltot = fltots*(1.-snowfrac) + fltot*snowfrac
2006 ALB = MAX(keep_snow_albedo*alb, &
2007 MIN((alb_snow_free + (alb - alb_snow_free) * snowfrac), alb))
2009 Emiss = MAX(keep_snow_albedo*emissn, &
2010 MIN((emiss_snowfree + &
2011 (emissn - emiss_snowfree) * snowfrac), emissn))
2013 ! alb=alb_snow_free*(1.-snowfrac) + alb*snowfrac
2014 ! emiss=emiss_snowfree*(1.-snowfrac) + emissn*snowfrac
2016 ! if(abs(fltot) > 2.) then
2017 ! print *,'i,j,fltot,snowfrac,fltots',fltot,snowfrac,fltots,i,j
2019 runoff1 = runoff1s*(1.-snowfrac) + runoff1*snowfrac
2020 runoff2 = runoff2s*(1.-snowfrac) + runoff2*snowfrac
2021 smelt = smelt * snowfrac
2022 snoh = snoh * snowfrac
2023 snflx = snflx * snowfrac
2024 snom = snom * snowfrac
2025 mavail = mavails*(1.-snowfrac) + 1.*snowfrac
2026 infiltr = infiltrs*(1.-snowfrac) + infiltr*snowfrac
2028 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
2029 print *,' Ground flux combined', i,j, s
2030 print *,'SOILT combined on land', soilt
2031 print *,'TS combined on land', ts1d
2035 ! Now combine fluxes for snow-free sea ice and snow-covered area
2036 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
2037 print *,'SOILT snow on ice', soilt
2040 ts1d(k) = ts1ds(k)*(1.-snowfrac) + ts1d(k)*snowfrac
2042 dew = dews*(1.-snowfrac) + dew*snowfrac
2043 soilt = soilts*(1.-snowfrac) + soilt*snowfrac
2044 qvg = qvgs*(1.-snowfrac) + qvg*snowfrac
2045 qsg = qsgs*(1.-snowfrac) + qsg*snowfrac
2046 qcg = qcgs*(1.-snowfrac) + qcg*snowfrac
2047 eeta = eetas*(1.-snowfrac) + eeta*snowfrac
2048 qfx = qfxs*(1.-snowfrac) + qfx*snowfrac
2049 hfx = hfxs*(1.-snowfrac) + hfx*snowfrac
2050 s = ss*(1.-snowfrac) + s*snowfrac
2052 prcpl = prcpls*(1.-snowfrac) + prcpl*snowfrac
2053 fltot = fltots*(1.-snowfrac) + fltot*snowfrac
2055 ALB = MAX(keep_snow_albedo*alb, &
2056 MIN((albice + (alb - alb_snow_free) * snowfrac), alb))
2058 Emiss = MAX(keep_snow_albedo*emissn, &
2059 MIN((emiss_snowfree + &
2060 (emissn - emiss_snowfree) * snowfrac), emissn))
2062 ! alb=alb_snow_free*(1.-snowfrac) + alb*snowfrac
2063 ! emiss=1.*(1.-snowfrac) + emissn*snowfrac
2064 runoff1 = runoff1s*(1.-snowfrac) + runoff1*snowfrac
2065 runoff2 = runoff2s*(1.-snowfrac) + runoff2*snowfrac
2066 smelt = smelt * snowfrac
2067 snoh = snoh * snowfrac
2068 snflx = snflx * snowfrac
2069 snom = snom * snowfrac
2070 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
2071 print *,'SOILT combined on ice', soilt
2074 endif ! snow_mosaic = 1.
2076 ! run-total accumulated snow based on snowfall and snowmelt in [m]
2078 snowfallac = snowfallac + max(0.,(newsn - rhowater/rhonewsn*smelt*delt*newsnowratio))
2087 T3 = STBOLT*SOILT*SOILT*SOILT
2089 XINET = EMISS*(GLW-UPFLUX)
2090 RNET = GSWnew + XINET
2091 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
2092 print *,'NO snow on the ground GSWnew -',GSWnew,'RNET=',rnet
2095 if(SEAICE .LT. 0.5) then
2097 CALL SOIL(spp_lsm,rstochcol,fieldcol_sf, &
2098 !--- input variables
2099 i,j,iland,isoil,delt,ktau,conflx,nzs,nddzs,nroot, &
2100 PRCPMS,RAINF,PATM,QVATM,QCATM,GLW,GSWnew,GSWin, &
2101 EMISS,RNET,QKMS,TKMS,PC,cst,drip,infwater, &
2102 rho,vegfrac,lai,myj, &
2103 !--- soil fixed fields
2104 QWRTZ,rhocs,dqm,qmin,ref,wilt, &
2105 psis,bclh,ksat,sat,cn, &
2106 zsmain,zshalf,DTDZS,DTDZS2,tbq, &
2108 lv,CP,rovcp,G0,cw,stbolt,tabs, &
2110 !--- output variables
2111 soilm1d,ts1d,smfrkeep,keepfr, &
2112 dew,soilt,qvg,qsg,qcg,edir1,ec1, &
2113 ett1,eeta,qfx,hfx,s,evapl,prcpl,fltot,runoff1, &
2114 runoff2,mavail,soilice,soiliqw, &
2118 ! If current ice albedo is not the same as from the previous time step, then
2119 ! update GSW, ALB and RNET for surface energy budget
2120 if(ALB.ne.ALBice) GSWnew=GSW/(1.-ALB)*(1.-ALBice)
2122 RNET = GSWnew + XINET
2125 !--- input variables
2126 i,j,iland,isoil,delt,ktau,conflx,nzs,nddzs,nroot, &
2127 PRCPMS,RAINF,PATM,QVATM,QCATM,GLW,GSWnew, &
2128 EMISS,RNET,QKMS,TKMS,rho,myj, &
2129 !--- sea ice parameters
2130 tice,rhosice,capice,thdifice, &
2131 zsmain,zshalf,DTDZS,DTDZS2,tbq, &
2133 lv,CP,rovcp,cw,stbolt,tabs, &
2134 !--- output variables
2135 ts1d,dew,soilt,qvg,qsg,qcg, &
2136 eeta,qfx,hfx,s,evapl,prcpl,fltot &
2159 !---------------------------------------------------------------
2160 END SUBROUTINE SFCTMP
2161 !---------------------------------------------------------------
2165 !****************************************************************
2166 REAL, DIMENSION(1:5001), INTENT(IN ) :: T
2167 REAL, INTENT(IN ) :: TN
2172 R=(TN-173.15)/.05+1.
2177 10 IF(I.LE.5000) GOTO 20
2182 QSN=(T(I+1)-R1)*R2 + R1
2183 ! print *,' in QSN, I,R,R1,R2,T(I+1),TN, QSN', I,R,r1,r2,t(i+1),tn,QSN
2186 !-----------------------------------------------------------------------
2188 !------------------------------------------------------------------------
2191 SUBROUTINE SOIL (spp_lsm,rstochcol, fieldcol_sf, &
2192 !--- input variables
2193 i,j,iland,isoil,delt,ktau,conflx,nzs,nddzs,nroot,&
2194 PRCPMS,RAINF,PATM,QVATM,QCATM, &
2195 GLW,GSW,GSWin,EMISS,RNET, &
2196 QKMS,TKMS,PC,cst,drip,infwater,rho,vegfrac,lai, &
2198 !--- soil fixed fields
2199 QWRTZ,rhocs,dqm,qmin,ref,wilt,psis,bclh,ksat, &
2200 sat,cn,zsmain,zshalf,DTDZS,DTDZS2,tbq, &
2202 xlv,CP,rovcp,G0_P,cw,stbolt,TABS, &
2204 !--- output variables
2205 soilmois,tso,smfrkeep,keepfr, &
2206 dew,soilt,qvg,qsg,qcg, &
2207 edir1,ec1,ett1,eeta,qfx,hfx,s,evapl, &
2208 prcpl,fltot,runoff1,runoff2,mavail,soilice, &
2209 soiliqw,infiltrp,smf)
2211 !*************************************************************
2212 ! Energy and moisture budget for vegetated surfaces
2213 ! without snow, heat diffusion and Richards eqns. in
2216 ! DELT - time step (s)
2217 ! ktau - numver of time step
2218 ! CONFLX - depth of constant flux layer (m)
2219 ! J,I - the location of grid point
2220 ! IME, JME, KME, NZS - dimensions of the domain
2221 ! NROOT - number of levels within the root zone
2222 ! PRCPMS - precipitation rate in m/s
2223 ! PATM - pressure [bar]
2224 ! QVATM,QCATM - cloud and water vapor mixing ratio (kg/kg)
2225 ! at the first atm. level
2226 ! GLW, GSW - incoming longwave and absorbed shortwave
2227 ! radiation at the surface (W/m^2)
2228 ! EMISS,RNET - emissivity of the ground surface (0-1) and net
2229 ! radiation at the surface (W/m^2)
2230 ! QKMS - exchange coefficient for water vapor in the
2231 ! surface layer (m/s)
2232 ! TKMS - exchange coefficient for heat in the surface
2234 ! PC - plant coefficient (resistance) (0-1)
2235 ! RHO - density of atmosphere near sueface (kg/m^3)
2236 ! VEGFRAC - greeness fraction
2237 ! RHOCS - volumetric heat capacity of dry soil
2238 ! DQM, QMIN - porosity minus residual soil moisture QMIN (m^3/m^3)
2239 ! REF, WILT - field capacity soil moisture and the
2240 ! wilting point (m^3/m^3)
2241 ! PSIS - matrix potential at saturation (m)
2242 ! BCLH - exponent for Clapp-Hornberger parameterization
2243 ! KSAT - saturated hydraulic conductivity (m/s)
2244 ! SAT - maximum value of water intercepted by canopy (m)
2245 ! CN - exponent for calculation of canopy water
2246 ! ZSMAIN - main levels in soil (m)
2247 ! ZSHALF - middle of the soil layers (m)
2248 ! DTDZS,DTDZS2 - dt/(2.*dzshalf*dzmain) and dt/dzshalf in soil
2249 ! TBQ - table to define saturated mixing ration
2250 ! of water vapor for given temperature and pressure
2251 ! SOILMOIS,TSO - soil moisture (m^3/m^3) and temperature (K)
2252 ! DEW - dew in kg/m^2s
2253 ! SOILT - skin temperature (K)
2254 ! QSG,QVG,QCG - saturated mixing ratio, mixing ratio of
2255 ! water vapor and cloud at the ground
2256 ! surface, respectively (kg/kg)
2257 ! EDIR1, EC1, ETT1, EETA - direct evaporation, evaporation of
2258 ! canopy water, transpiration in kg m-2 s-1 and total
2259 ! evaporation in m s-1.
2260 ! QFX, HFX - latent and sensible heat fluxes (W/m^2)
2261 ! S - soil heat flux in the top layer (W/m^2)
2262 ! RUNOFF - surface runoff (m/s)
2263 ! RUNOFF2 - underground runoff (m)
2264 ! MAVAIL - moisture availability in the top soil layer (0-1)
2265 ! INFILTRP - infiltration flux from the top of soil domain (m/s)
2267 !*****************************************************************
2269 !-----------------------------------------------------------------
2271 !--- input variables
2273 INTEGER, INTENT(IN ) :: nroot,ktau,nzs , &
2274 nddzs !nddzs=2*(nzs-2)
2275 INTEGER, INTENT(IN ) :: i,j,iland,isoil
2276 REAL, INTENT(IN ) :: DELT,CONFLX
2277 LOGICAL, INTENT(IN ) :: myj
2278 !--- 3-D Atmospheric variables
2280 INTENT(IN ) :: PATM, &
2285 INTENT(IN ) :: GLW, &
2297 !--- soil properties
2299 INTENT(IN ) :: RHOCS, &
2309 REAL, INTENT(IN ) :: CN, &
2318 REAL, DIMENSION(1:NZS), INTENT(IN) :: ZSMAIN, &
2322 REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS
2324 REAL, DIMENSION(1:5001), INTENT(IN) :: TBQ
2327 !--- input/output variables
2328 !-------- 3-d soil moisture and temperature
2329 REAL, DIMENSION( 1:nzs ) , &
2330 INTENT(INOUT) :: TSO, &
2334 REAL, DIMENSION(1:NZS), INTENT(IN) :: rstochcol
2335 REAL, DIMENSION(1:NZS), INTENT(INOUT) :: fieldcol_sf
2338 REAL, DIMENSION( 1:nzs ) , &
2339 INTENT(INOUT) :: KEEPFR
2341 !-------- 2-d variables
2343 INTENT(INOUT) :: DEW, &
2365 !-------- 1-d variables
2366 INTEGER , INTENT(IN) :: spp_lsm
2367 REAL, DIMENSION(1:NZS), INTENT(OUT) :: SOILICE, &
2370 !--- Local variables
2372 REAL :: INFILTRP, transum , &
2374 TABS, T3, UPFLUX, XINET
2375 REAL :: CP,rovcp,G0,LV,STBOLT,xlmelt,dzstop , &
2376 can,epot,fac,fltot,ft,fq,hft , &
2377 q1,ras,rhoice,sph , &
2378 trans,zn,ci,cvw,tln,tavln,pi , &
2379 DD1,CMC2MS,DRYCAN,WETCAN , &
2381 REAL, DIMENSION(1:NZS) :: transp,cap,diffu,hydro , &
2382 thdif,tranf,tav,soilmoism , &
2383 soilicem,soiliqwm,detal , &
2384 fwsat,lwsat,told,smold
2386 REAL :: soiltold,smf
2387 REAL :: soilres, alfa, fex, fex_fc, fc, psit
2389 INTEGER :: nzs1,nzs2,k
2391 !-----------------------------------------------------------------
2393 !-- define constants
2394 ! STBOLT=5.670151E-8
2409 !--- Initializing local arrays
2432 dzstop=1./(zsmain(2)-zsmain(1))
2436 !--- Computation of volumetric content of ice in soil
2440 tln=log(tso(k)/273.15)
2442 soiliqw(k)=(dqm+qmin)*(XLMELT* &
2443 (tso(k)-273.15)/tso(k)/9.81/psis) &
2445 soiliqw(k)=max(0.,soiliqw(k))
2446 soiliqw(k)=min(soiliqw(k),soilmois(k))
2447 soilice(k)=(soilmois(k)-soiliqw(k))/RIW
2449 !---- melting and freezing is balanced, soil ice cannot increase
2450 if(keepfr(k).eq.1.) then
2451 soilice(k)=min(soilice(k),smfrkeep(k))
2452 soiliqw(k)=max(0.,soilmois(k)-soilice(k)*riw)
2457 soiliqw(k)=soilmois(k)
2463 !- middle of soil layers
2464 tav(k)=0.5*(tso(k)+tso(k+1))
2465 soilmoism(k)=0.5*(soilmois(k)+soilmois(k+1))
2466 tavln=log(tav(k)/273.15)
2468 if(tavln.lt.0.) then
2469 soiliqwm(k)=(dqm+qmin)*(XLMELT* &
2470 (tav(k)-273.15)/tav(k)/9.81/psis) &
2472 fwsat(k)=dqm-soiliqwm(k)
2473 lwsat(k)=soiliqwm(k)+qmin
2474 soiliqwm(k)=max(0.,soiliqwm(k))
2475 soiliqwm(k)=min(soiliqwm(k), soilmoism(k))
2476 soilicem(k)=(soilmoism(k)-soiliqwm(k))/riw
2477 !---- melting and freezing is balanced, soil ice cannot increase
2478 if(keepfr(k).eq.1.) then
2479 soilicem(k)=min(soilicem(k), &
2480 0.5*(smfrkeep(k)+smfrkeep(k+1)))
2481 soiliqwm(k)=max(0.,soilmoism(k)-soilicem(k)*riw)
2482 fwsat(k)=dqm-soiliqwm(k)
2483 lwsat(k)=soiliqwm(k)+qmin
2488 soiliqwm(k)=soilmoism(k)
2496 if(soilice(k).gt.0.) then
2497 smfrkeep(k)=soilice(k)
2499 smfrkeep(k)=soilmois(k)/riw
2503 !******************************************************************
2504 ! SOILPROP computes thermal diffusivity, and diffusional and
2505 ! hydraulic condeuctivities
2506 !******************************************************************
2507 CALL SOILPROP(spp_lsm,rstochcol,fieldcol_sf, &
2508 !--- input variables
2509 nzs,fwsat,lwsat,tav,keepfr, &
2510 soilmois,soiliqw,soilice, &
2511 soilmoism,soiliqwm,soilicem, &
2512 !--- soil fixed fields
2513 QWRTZ,rhocs,dqm,qmin,psis,bclh,ksat, &
2515 riw,xlmelt,CP,G0_P,cvw,ci, &
2517 !--- output variables
2518 thdif,diffu,hydro,cap)
2520 !********************************************************************
2521 !--- CALCULATION OF CANOPY WATER (Smirnova et al., 1996, EQ.16) AND DEW
2528 Q1=-QKMS*RAS*(QVATM - QSG)
2531 IF(QVATM.GE.QSG)THEN
2536 ! DD1=CST+DELT*(PRCPMS +DEW*RAS)
2539 ! DELT*(PRCPMS+RAS*FQ*(QVATM-QSG) &
2543 ! DD1=CST+DELT*PRCPMS
2545 ! IF(DD1.LT.0.) DD1=0.
2546 ! if(vegfrac.eq.0.)then
2550 ! IF (vegfrac.GT.0.) THEN
2552 ! IF(CST.GT.SAT) THEN
2558 !--- WETCAN is the fraction of vegetated area covered by canopy
2559 !--- water, and DRYCAN is the fraction of vegetated area where
2560 !--- transpiration may take place.
2562 WETCAN=min(0.25,max(0.,(CST/SAT))**CN)
2563 ! if(lai > 1.) wetcan=wetcan/lai
2566 !**************************************************************
2567 ! TRANSF computes transpiration function
2568 !**************************************************************
2570 !--- input variables
2571 nzs,nroot,soiliqw,tabs,lai,gswin, &
2572 !--- soil fixed fields
2573 dqm,qmin,ref,wilt,zshalf,pc,iland, &
2574 !--- output variables
2578 !--- Save soil temp and moisture from the beginning of time step
2581 smold(k)=soilmois(k)
2584 ! Sakaguchi and Zeng (2009) - dry soil resistance to evaporation
2585 ! if (vgtype==11) then ! MODIS wetland
2588 fex=min(1.,soilmois(1)/dqm)
2590 psit=psis*fex ** (-bclh)
2591 psit = max(-1.e5, psit)
2592 alfa=min(1.,exp(g*psit/r_v/SOILT))
2596 fc=max(qmin,ref*0.5)
2598 if((soilmois(1)+qmin) > fc .or. (qvatm-qvg) > 0.) then
2601 fex_fc=min(1.,(soilmois(1)+qmin)/fc)
2602 fex_fc=max(fex_fc,0.01)
2603 soilres=0.25*(1.-cos(piconst*fex_fc))**2.
2605 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
2606 ! if (i==421.and.j==280) then
2607 print *,'fex,psit,psis,bclh,g,r_v,soilt,alfa,mavail,soilmois(1),fc,ref,soilres,fex_fc', &
2608 fex,psit,psis,bclh,g,r_v,soilt,alfa,mavail,soilmois(1),fc,ref,soilres,fex_fc
2611 !**************************************************************
2612 ! SOILTEMP soilves heat budget and diffusion eqn. in soil
2613 !**************************************************************
2616 !--- input variables
2618 delt,ktau,conflx,nzs,nddzs,nroot, &
2620 PATM,TABS,QVATM,QCATM,EMISS,RNET, &
2621 QKMS,TKMS,PC,rho,vegfrac, lai, &
2622 thdif,cap,drycan,wetcan, &
2623 transum,dew,mavail,soilres,alfa, &
2624 !--- soil fixed fields
2625 dqm,qmin,bclh,zsmain,zshalf,DTDZS,tbq, &
2627 xlv,CP,G0_P,cvw,stbolt, &
2628 !--- output variables
2629 tso,soilt,qvg,qsg,qcg,x)
2631 !************************************************************************
2633 !--- CALCULATION OF DEW USING NEW VALUE OF QSG OR TRANSP IF NO DEW
2637 IF(QVATM.GE.QSG)THEN
2638 DEW=QKMS*(QVATM-QSG)
2646 TRANSP(K)=VEGFRAC*RAS*QKMS* &
2648 TRANF(K)*DRYCAN/ZSHALF(NROOT+1)
2649 IF(TRANSP(K).GT.0.) TRANSP(K)=0.
2657 !-- Recalculate volumetric content of frozen water in soil
2660 tln=log(tso(k)/273.15)
2662 soiliqw(k)=(dqm+qmin)*(XLMELT* &
2663 (tso(k)-273.15)/tso(k)/9.81/psis) &
2665 soiliqw(k)=max(0.,soiliqw(k))
2666 soiliqw(k)=min(soiliqw(k),soilmois(k))
2667 soilice(k)=(soilmois(k)-soiliqw(k))/riw
2668 !---- melting and freezing is balanced, soil ice cannot increase
2669 if(keepfr(k).eq.1.) then
2670 soilice(k)=min(soilice(k),smfrkeep(k))
2671 soiliqw(k)=max(0.,soilmois(k)-soilice(k)*riw)
2676 soiliqw(k)=soilmois(k)
2680 !*************************************************************************
2681 ! SOILMOIST solves moisture budget (Smirnova et al., 1996, EQ.22,28)
2683 !*************************************************************************
2686 delt,nzs,nddzs,DTDZS,DTDZS2,RIW, &
2687 zsmain,zshalf,diffu,hydro, &
2688 QSG,QVG,QCG,QCATM,QVATM,-infwater, &
2689 QKMS,TRANSP,DRIP,DEW,0.,SOILICE,VEGFRAC, &
2692 DQM,QMIN,REF,KSAT,RAS,INFMAX, &
2694 SOILMOIS,SOILIQW,MAVAIL,RUNOFF1, &
2697 !--- KEEPFR is 1 when the temperature and moisture in soil
2698 !--- are both increasing. In this case soil ice should not
2699 !--- be increasing according to the freezing curve.
2700 !--- Some part of ice is melted, but additional water is
2701 !--- getting frozen. Thus, only structure of frozen soil is
2702 !--- changed, and phase changes are not affecting the heat
2703 !--- transfer. This situation may happen when it rains on the
2707 if (soilice(k).gt.0.) then
2708 if(tso(k).gt.told(k).and.soilmois(k).gt.smold(k)) then
2716 !--- THE DIAGNOSTICS OF SURFACE FLUXES
2718 T3 = STBOLT*SOILTold*SOILTold*SOILTold
2719 UPFLUX = T3 * 0.5*(SOILTold+SOILT)
2720 XINET = EMISS*(GLW-UPFLUX)
2721 ! RNET = GSW + XINET
2722 HFT=-TKMS*CP*RHO*(TABS-SOILT)
2723 HFX=-TKMS*CP*RHO*(TABS-SOILT) &
2724 *(P1000mb*0.00001/Patm)**ROVCP
2725 Q1=-QKMS*RAS*(QVATM - QSG)
2734 !-- moisture flux for coupling with MYJ PBL
2735 EETA=-QKMS*RAS*(QVATM/(1.+QVATM) - QSG/(1.+QSG))*1.E3
2736 CST= CST-EETA*DELT*vegfrac
2737 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
2738 !!! IF(i.eq.374.and.j.eq.310.or. EETA.gt.0.0004) then
2739 print *,'Cond MYJ EETA',eeta,eeta*xlv, i,j
2742 !-- actual moisture flux from RUC LSM
2744 CST=CST+DELT*DEW*RAS * vegfrac
2745 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
2746 ! IF(i.eq.374.and.j.eq.310.or. EETA.gt.0.0004) then
2747 ! IF(i.eq.440.and.j.eq.180.or. QFX.gt.1000..or.i.eq.417.and.j.eq.540) then
2748 print *,'Cond RUC LSM EETA',EETA,eeta*xlv, i,j
2755 EDIR1 =-soilres*(1.-vegfrac)*QKMS*RAS* &
2758 EC1 = Q1 * WETCAN * vegfrac
2759 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
2760 ! IF(i.eq.440.and.j.eq.180.or. QFX.gt.1000..or.i.eq.417.and.j.eq.540) then
2761 print *,'CST before update=',cst
2762 print *,'EC1=',EC1,'CMC2MS=',CMC2MS
2766 CST=max(0.,CST-EC1 * DELT)
2768 ! if (EC1 > CMC2MS) then
2769 ! EC1 = min(cmc2ms,ec1)
2774 !-- moisture flux for coupling with MYJ PBL
2775 EETA=-soilres*QKMS*RAS*(QVATM/(1.+QVATM) - QVG/(1.+QVG))*1.E3
2777 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
2778 ! IF(i.eq.440.and.j.eq.180.or. QFX.gt.1000..or.i.eq.417.and.j.eq.540) then
2779 print *,'QKMS,RAS,QVATM/(1.+QVATM),QVG/(1.+QVG),QSG ', &
2780 QKMS,RAS,QVATM/(1.+QVATM),QVG/(1.+QVG),QSG
2781 print *,'Q1*(1.-vegfrac),EDIR1',Q1*(1.-vegfrac),EDIR1
2782 print *,'CST,WETCAN,DRYCAN',CST,WETCAN,DRYCAN
2783 print *,'EC1=',EC1,'ETT1=',ETT1,'CMC2MS=',CMC2MS,'CMC2MS*ras=',CMC2MS*ras
2784 ! print *,'MYJ EETA',eeta,eeta*xlv
2786 !-- actual moisture flux from RUC LSM
2787 EETA = (EDIR1 + EC1 + ETT1)*1.E3
2788 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
2789 ! IF(i.eq.374.and.j.eq.310.or. EETA.gt.0.0004) then
2790 ! IF(i.eq.440.and.j.eq.180 .or. qfx.gt.1000..or.i.eq.417.and.j.eq.540) then
2791 print *,'RUC LSM EETA',EETA,eeta*xlv
2795 EETA = (EDIR1 + EC1 + ETT1)*1.E3
2797 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
2798 print *,'potential temp HFT ',HFT
2799 print *,'abs temp HFX ',HFX
2803 S=THDIF(1)*CAP(1)*DZSTOP*(TSO(1)-TSO(2))
2805 FLTOT=RNET-HFT-XLV*EETA-S-X
2806 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
2807 ! IF(i.eq.440.and.j.eq.180 .or. qfx.gt.1000..or.i.eq.417.and.j.eq.540) then
2808 print *,'SOIL - FLTOT,RNET,HFT,QFX,S,X=',i,j,FLTOT,RNET,HFT,XLV*EETA,s,x
2809 print *,'edir1,ec1,ett1,mavail,qkms,qvatm,qvg,qsg,vegfrac',&
2810 edir1,ec1,ett1,mavail,qkms,qvatm,qvg,qsg,vegfrac
2812 if(detal(1) .ne. 0.) then
2813 ! SMF - energy of phase change in the first soil layer
2814 ! smf=xlmelt*1.e3*(soiliqwm(1)-soiliqwmold(1))/delt
2816 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
2817 print *,'detal(1),xlmelt,soiliqwm(1),delt',detal(1),xlmelt,soiliqwm(1),delt
2818 print *,'Implicit phase change in the first layer - smf=',smf
2825 1123 FORMAT(I5,8F12.3)
2826 1133 FORMAT(I7,8E12.4)
2827 123 format(i6,f6.2,7f8.1)
2828 122 FORMAT(1X,2I3,6F8.1,F8.3,F8.2)
2829 !-------------------------------------------------------------------
2831 !-------------------------------------------------------------------
2834 !--- input variables
2835 i,j,iland,isoil,delt,ktau,conflx,nzs,nddzs,nroot, &
2836 PRCPMS,RAINF,PATM,QVATM,QCATM,GLW,GSW, &
2837 EMISS,RNET,QKMS,TKMS,rho,myj, &
2838 !--- sea ice parameters
2839 tice,rhosice,capice,thdifice, &
2840 zsmain,zshalf,DTDZS,DTDZS2,tbq, &
2842 xlv,CP,rovcp,cw,stbolt,tabs, &
2843 !--- output variables
2844 tso,dew,soilt,qvg,qsg,qcg, &
2845 eeta,qfx,hfx,s,evapl,prcpl,fltot &
2848 !*****************************************************************
2849 ! Energy budget and heat diffusion eqns. for
2851 !*************************************************************
2854 !-----------------------------------------------------------------
2856 !--- input variables
2858 INTEGER, INTENT(IN ) :: nroot,ktau,nzs , &
2859 nddzs !nddzs=2*(nzs-2)
2860 INTEGER, INTENT(IN ) :: i,j,iland,isoil
2861 REAL, INTENT(IN ) :: DELT,CONFLX
2862 LOGICAL, INTENT(IN ) :: myj
2863 !--- 3-D Atmospheric variables
2865 INTENT(IN ) :: PATM, &
2870 INTENT(IN ) :: GLW, &
2876 !--- sea ice properties
2877 REAL, DIMENSION(1:NZS) , &
2885 REAL, INTENT(IN ) :: &
2890 REAL, DIMENSION(1:NZS), INTENT(IN) :: ZSMAIN, &
2894 REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS
2896 REAL, DIMENSION(1:5001), INTENT(IN) :: TBQ
2899 !--- input/output variables
2900 !----soil temperature
2901 REAL, DIMENSION( 1:nzs ), INTENT(INOUT) :: TSO
2902 !-------- 2-d variables
2904 INTENT(INOUT) :: DEW, &
2917 !--- Local variables
2918 REAL :: x,x1,x2,x4,tn,denom
2919 REAL :: RAINF, PRCPMS , &
2920 TABS, T3, UPFLUX, XINET
2922 REAL :: CP,rovcp,G0,LV,STBOLT,xlmelt,dzstop , &
2923 epot,fltot,ft,fq,hft,ras,cvw
2925 REAL :: FKT,D1,D2,D9,D10,DID,R211,R21,R22,R6,R7,D11 , &
2926 PI,H,FKQ,R210,AA,BB,PP,Q1,QS1,TS1,TQ2,TX2 , &
2929 REAL :: AA1,RHCS, icemelt
2932 REAL, DIMENSION(1:NZS) :: cotso,rhtso
2934 INTEGER :: nzs1,nzs2,k,k1,kn,kk
2936 !-----------------------------------------------------------------
2938 !-- define constants
2939 ! STBOLT=5.670151E-8
2947 dzstop=1./(zsmain(2)-zsmain(1))
2961 X1=DTDZS(K1)*THDIFICE(KN-1)
2962 X2=DTDZS(K1+1)*THDIFICE(KN)
2963 FT=TSO(KN)+X1*(TSO(KN-1)-TSO(KN)) &
2964 -X2*(TSO(KN)-TSO(KN+1))
2965 DENOM=1.+X1+X2-X2*cotso(K)
2967 rhtso(K+1)=(FT+X2*rhtso(K))/DENOM
2970 !************************************************************************
2971 !--- THE HEAT BALANCE EQUATION (Smirnova et al., 1996, EQ. 21,26)
2978 D9=THDIFICE(1)*RHCS*dzstop
2982 R22=.5/(THDIFICE(1)*DELT*dzstop**2)
2983 R6=EMISS *STBOLT*.5*TN**4
2986 TDENOM=D9*(1.-D1+R22)+D10+R21+R7 &
2990 AA=XLS*(FKQ+R210)/TDENOM
2991 BB=(D10*TABS+R21*TN+XLS*(QVATM*FKQ &
2992 +R210*QVG)+D11+D9*(D2+R22*TN) &
2993 +RAINF*CVW*PRCPMS*max(273.15,TABS) &
2998 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
2999 PRINT *,' VILKA-SEAICE1'
3000 print *,'D10,TABS,R21,TN,QVATM,FKQ', &
3001 D10,TABS,R21,TN,QVATM,FKQ
3002 print *,'RNET, EMISS, STBOLT, SOILT',RNET, EMISS, STBOLT, SOILT
3003 print *,'R210,QVG,D11,D9,D2,R22,RAINF,CVW,PRCPMS,TDENOM', &
3004 R210,QVG,D11,D9,D2,R22,RAINF,CVW,PRCPMS,TDENOM
3005 print *,'tn,aa1,bb,pp,fkq,r210', &
3006 tn,aa1,bb,pp,fkq,r210
3009 CALL VILKA(TN,AA1,BB,PP,QS1,TS1,TBQ,KTAU,i,j,iland,isoil)
3010 !--- it is saturation over sea ice
3013 TSO(1)=min(271.4,TS1)
3015 !--- sea ice melting is not included in this simple approach
3016 !--- SOILT - skin temperature
3018 !---- Final solution for soil temperature - TSO
3021 TSO(K)=min(271.4,rhtso(KK)+cotso(KK)*TSO(K-1))
3023 !--- CALCULATION OF DEW USING NEW VALUE OF QSG OR TRANSP IF NO DEW
3026 !--- THE DIAGNOSTICS OF SURFACE FLUXES
3027 T3 = STBOLT*TN*TN*TN
3028 UPFLUX = T3 *0.5*(TN+SOILT)
3029 XINET = EMISS*(GLW-UPFLUX)
3030 ! RNET = GSW + XINET
3031 HFT=-TKMS*CP*RHO*(TABS-SOILT)
3032 HFX=-TKMS*CP*RHO*(TABS-SOILT) &
3033 *(P1000mb*0.00001/Patm)**ROVCP
3034 Q1=-QKMS*RAS*(QVATM - QSG)
3038 !-- moisture flux for coupling with MYJ PBL
3039 EETA=-QKMS*RAS*(QVATM/(1.+QVATM) - QSG/(1.+QSG))*1.E3
3040 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
3041 print *,'MYJ EETA',eeta
3044 !-- actual moisture flux from RUC LSM
3045 DEW=QKMS*(QVATM-QSG)
3047 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
3048 print *,'RUC LSM EETA',eeta
3056 !-- moisture flux for coupling with MYJ PBL
3057 EETA=-QKMS*RAS*(QVATM/(1.+QVATM) - QVG/(1.+QVG))*1.E3
3058 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
3059 print *,'MYJ EETA',eeta
3062 ! to convert from m s-1 to kg m-2 s-1: *rho water=1.e3************
3063 !-- actual moisture flux from RUC LSM
3065 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
3066 print *,'RUC LSM EETA',eeta
3074 S=THDIFICE(1)*CAPICE(1)*DZSTOP*(TSO(1)-TSO(2))
3075 ! heat storage in surface layer
3078 X= (cp*rho*r211+rhcs*zsmain(2)*0.5/delt)*(SOILT-TN) + &
3079 XLS*rho*r211*(QSG-QGOLD)
3082 -RAINF*CVW*PRCPMS*(max(273.15,TABS)-SOILT)
3084 !-- excess energy spent on sea ice melt
3085 icemelt=RNET-XLS*EETA -HFT -S -X
3086 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
3087 print *,'icemelt=',icemelt
3090 FLTOT=RNET-XLS*EETA-HFT-S-X-icemelt
3091 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
3092 print *,'SICE - FLTOT,RNET,HFT,QFX,S,SNOH,X=', &
3093 FLTOT,RNET,HFT,XLS*EETA,s,icemelt,X
3096 !-------------------------------------------------------------------
3098 !-------------------------------------------------------------------
3102 SUBROUTINE SNOWSOIL (spp_lsm,rstochcol,fieldcol_sf,&
3103 !--- input variables
3104 i,j,isoil,delt,ktau,conflx,nzs,nddzs,nroot, &
3105 meltfactor,rhonewsn,SNHEI_CRIT, & ! new
3106 ILAND,PRCPMS,RAINF,NEWSNOW,snhei,SNWE,SNOWFRAC, &
3109 GLW,GSW,GSWin,EMISS,RNET,IVGTYP, &
3110 QKMS,TKMS,PC,cst,drip,infwater, &
3111 rho,vegfrac,alb,znt,lai, &
3113 !--- soil fixed fields
3114 QWRTZ,rhocs,dqm,qmin,ref,wilt,psis,bclh,ksat, &
3115 sat,cn,zsmain,zshalf,DTDZS,DTDZS2,tbq, &
3117 xlv,CP,rovcp,G0_P,cw,stbolt,TABS, &
3119 !--- output variables
3120 ilnb,snweprint,snheiprint,rsm, &
3121 soilmois,tso,smfrkeep,keepfr, &
3122 dew,soilt,soilt1,tsnav, &
3123 qvg,qsg,qcg,SMELT,SNOH,SNFLX,SNOM, &
3124 edir1,ec1,ett1,eeta,qfx,hfx,s,sublim, &
3125 prcpl,fltot,runoff1,runoff2,mavail,soilice, &
3128 !***************************************************************
3129 ! Energy and moisture budget for snow, heat diffusion eqns.
3130 ! in snow and soil, Richards eqn. for soil covered with snow
3132 ! DELT - time step (s)
3133 ! ktau - numver of time step
3134 ! CONFLX - depth of constant flux layer (m)
3135 ! J,I - the location of grid point
3136 ! IME, JME, NZS - dimensions of the domain
3137 ! NROOT - number of levels within the root zone
3138 ! PRCPMS - precipitation rate in m/s
3139 ! NEWSNOW - pcpn in soilid form (m)
3140 ! SNHEI, SNWE - snow height and snow water equivalent (m)
3141 ! RHOSN - snow density (kg/m-3)
3142 ! PATM - pressure (bar)
3143 ! QVATM,QCATM - cloud and water vapor mixing ratio
3144 ! at the first atm. level (kg/kg)
3145 ! GLW, GSW - incoming longwave and absorbed shortwave
3146 ! radiation at the surface (W/m^2)
3147 ! EMISS,RNET - emissivity (0-1) of the ground surface and net
3148 ! radiation at the surface (W/m^2)
3149 ! QKMS - exchange coefficient for water vapor in the
3150 ! surface layer (m/s)
3151 ! TKMS - exchange coefficient for heat in the surface
3153 ! PC - plant coefficient (resistance) (0-1)
3154 ! RHO - density of atmosphere near surface (kg/m^3)
3155 ! VEGFRAC - greeness fraction (0-1)
3156 ! RHOCS - volumetric heat capacity of dry soil (J/m^3/K)
3157 ! DQM, QMIN - porosity minus residual soil moisture QMIN (m^3/m^3)
3158 ! REF, WILT - field capacity soil moisture and the
3159 ! wilting point (m^3/m^3)
3160 ! PSIS - matrix potential at saturation (m)
3161 ! BCLH - exponent for Clapp-Hornberger parameterization
3162 ! KSAT - saturated hydraulic conductivity (m/s)
3163 ! SAT - maximum value of water intercepted by canopy (m)
3164 ! CN - exponent for calculation of canopy water
3165 ! ZSMAIN - main levels in soil (m)
3166 ! ZSHALF - middle of the soil layers (m)
3167 ! DTDZS,DTDZS2 - dt/(2.*dzshalf*dzmain) and dt/dzshalf in soil
3168 ! TBQ - table to define saturated mixing ration
3169 ! of water vapor for given temperature and pressure
3170 ! ilnb - number of layers in snow
3171 ! rsm - liquid water inside snow pack (m)
3172 ! SOILMOIS,TSO - soil moisture (m^3/m^3) and temperature (K)
3173 ! DEW - dew in (kg/m^2 s)
3174 ! SOILT - skin temperature (K)
3175 ! SOILT1 - snow temperature at 7.5 cm depth (K)
3176 ! TSNAV - average temperature of snow pack (C)
3177 ! QSG,QVG,QCG - saturated mixing ratio, mixing ratio of
3178 ! water vapor and cloud at the ground
3179 ! surface, respectively (kg/kg)
3180 ! EDIR1, EC1, ETT1, EETA - direct evaporation, evaporation of
3181 ! canopy water, transpiration (kg m-2 s-1) and total
3182 ! evaporation in (m s-1).
3183 ! QFX, HFX - latent and sensible heat fluxes (W/m^2)
3184 ! S - soil heat flux in the top layer (W/m^2)
3185 ! SUBLIM - snow sublimation (kg/m^2/s)
3186 ! RUNOFF1 - surface runoff (m/s)
3187 ! RUNOFF2 - underground runoff (m)
3188 ! MAVAIL - moisture availability in the top soil layer (0-1)
3189 ! SOILICE - content of soil ice in soil layers (m^3/m^3)
3190 ! SOILIQW - lliquid water in soil layers (m^3/m^3)
3191 ! INFILTRP - infiltration flux from the top of soil domain (m/s)
3192 ! XINET - net long-wave radiation (W/m^2)
3194 !*******************************************************************
3197 !-------------------------------------------------------------------
3198 !--- input variables
3200 INTEGER, INTENT(IN ) :: nroot,ktau,nzs , &
3201 nddzs !nddzs=2*(nzs-2)
3202 INTEGER, INTENT(IN ) :: i,j,isoil
3204 REAL, INTENT(IN ) :: DELT,CONFLX,PRCPMS , &
3205 RAINF,NEWSNOW,RHONEWSN, &
3206 SNHEI_CRIT,meltfactor
3208 LOGICAL, INTENT(IN ) :: myj
3210 !--- 3-D Atmospheric variables
3212 INTENT(IN ) :: PATM, &
3217 INTENT(IN ) :: GLW, &
3228 INTEGER, INTENT(IN ) :: IVGTYP
3229 !--- soil properties
3231 INTENT(IN ) :: RHOCS, &
3242 REAL, INTENT(IN ) :: CN, &
3251 REAL, DIMENSION(1:NZS), INTENT(IN) :: ZSMAIN, &
3255 REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS
3257 REAL, DIMENSION(1:5001), INTENT(IN) :: TBQ
3259 REAL, DIMENSION(1:NZS), INTENT(IN) :: rstochcol
3260 REAL, DIMENSION(1:NZS), INTENT(INOUT) :: fieldcol_sf
3262 !--- input/output variables
3263 !-------- 3-d soil moisture and temperature
3264 REAL, DIMENSION( 1:nzs ) , &
3265 INTENT(INOUT) :: TSO, &
3269 REAL, DIMENSION( 1:nzs ) , &
3270 INTENT(INOUT) :: KEEPFR
3273 INTEGER, INTENT(INOUT) :: ILAND
3276 !-------- 2-d variables
3278 INTENT(INOUT) :: DEW, &
3311 INTEGER, INTENT(INOUT) :: ILNB
3313 !-------- 1-d variables
3314 REAL, DIMENSION(1:NZS), INTENT(OUT) :: SOILICE, &
3317 REAL, INTENT(OUT) :: RSM, &
3320 INTEGER, INTENT(IN) :: spp_lsm
3321 !--- Local variables
3324 INTEGER :: nzs1,nzs2,k
3326 REAL :: INFILTRP, TRANSUM , &
3328 TABS, T3, UPFLUX, XINET , &
3329 BETA, SNWEPR,EPDT,PP
3330 REAL :: CP,rovcp,G0,LV,xlvm,STBOLT,xlmelt,dzstop , &
3331 can,epot,fac,fltot,ft,fq,hft , &
3332 q1,ras,rhoice,sph , &
3333 trans,zn,ci,cvw,tln,tavln,pi , &
3334 DD1,CMC2MS,DRYCAN,WETCAN , &
3335 INFMAX,RIW,DELTSN,H,UMVEG
3337 REAL, DIMENSION(1:NZS) :: transp,cap,diffu,hydro , &
3338 thdif,tranf,tav,soilmoism , &
3339 soilicem,soiliqwm,detal , &
3340 fwsat,lwsat,told,smold
3341 REAL :: soiltold, qgold
3345 !-----------------------------------------------------------------
3349 !-- heat of water vapor sublimation
3351 ! STBOLT=5.670151E-8
3353 !--- SNOW flag -- ISICE
3356 !--- DELTSN - is the threshold for splitting the snow layer into 2 layers.
3357 !--- With snow density 400 kg/m^3, this threshold is equal to 7.5 cm,
3358 !--- equivalent to 0.03 m SNWE. For other snow densities the threshold is
3359 !--- computed using SNWE=0.03 m and current snow density.
3360 !--- SNTH - the threshold below which the snow layer is combined with
3361 !--- the top soil layer. SNTH is computed using snwe=0.016 m, and
3362 !--- equals 4 cm for snow density 400 kg/m^3.
3370 ! increase thinkness of top snow layer from 3 cm SWE to 5 cm SWE
3371 ! DELTSN=5.*SNHEI_CRIT
3372 ! snth=0.4*SNHEI_CRIT
3374 DELTSN=0.05*1.e3/rhosn
3375 snth=0.01*1.e3/rhosn
3376 ! snth=0.01601*1.e3/rhosn
3378 ! if(i.eq.442.and.j.eq.260) then
3379 ! print *,'deltsn,snhei,snth',i,j,deltsn,snhei,snth
3382 ! For 2-layer snow model when the snow depth is marginally higher than DELTSN,
3383 ! reset DELTSN to half of snow depth.
3384 IF(SNHEI.GE.DELTSN+SNTH) THEN
3385 if(snhei-deltsn-snth.lt.snth) deltsn=0.5*(snhei-snth)
3386 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
3387 print *,'DELTSN is changed,deltsn,snhei,snth',i,j,deltsn,snhei,snth
3421 !*** DELTSN is the depth of the top layer of snow where
3422 !*** there is a temperature gradient, the rest of the snow layer
3423 !*** is considered to have constant temperature
3428 DZSTOP=1./(zsmain(2)-zsmain(1))
3430 !----- THE CALCULATION OF THERMAL DIFFUSIVITY, DIFFUSIONAL AND ---
3431 !----- HYDRAULIC CONDUCTIVITY (SMIRNOVA ET AL. 1996, EQ.2,5,6) ---
3432 !tgs - the following loop is added to define the amount of frozen
3433 !tgs - water in soil if there is any
3436 tln=log(tso(k)/273.15)
3438 soiliqw(k)=(dqm+qmin)*(XLMELT* &
3439 (tso(k)-273.15)/tso(k)/9.81/psis) &
3441 soiliqw(k)=max(0.,soiliqw(k))
3442 soiliqw(k)=min(soiliqw(k),soilmois(k))
3443 soilice(k)=(soilmois(k)-soiliqw(k))/riw
3445 !---- melting and freezing is balanced, soil ice cannot increase
3446 if(keepfr(k).eq.1.) then
3447 soilice(k)=min(soilice(k),smfrkeep(k))
3448 soiliqw(k)=max(0.,soilmois(k)-soilice(k)*rhoice*1.e-3)
3453 soiliqw(k)=soilmois(k)
3460 tav(k)=0.5*(tso(k)+tso(k+1))
3461 soilmoism(k)=0.5*(soilmois(k)+soilmois(k+1))
3462 tavln=log(tav(k)/273.15)
3464 if(tavln.lt.0.) then
3465 soiliqwm(k)=(dqm+qmin)*(XLMELT* &
3466 (tav(k)-273.15)/tav(k)/9.81/psis) &
3468 fwsat(k)=dqm-soiliqwm(k)
3469 lwsat(k)=soiliqwm(k)+qmin
3470 soiliqwm(k)=max(0.,soiliqwm(k))
3471 soiliqwm(k)=min(soiliqwm(k), soilmoism(k))
3472 soilicem(k)=(soilmoism(k)-soiliqwm(k))/riw
3473 !---- melting and freezing is balanced, soil ice cannot increase
3474 if(keepfr(k).eq.1.) then
3475 soilicem(k)=min(soilicem(k), &
3476 0.5*(smfrkeep(k)+smfrkeep(k+1)))
3477 soiliqwm(k)=max(0.,soilmoism(k)-soilicem(k)*riw)
3478 fwsat(k)=dqm-soiliqwm(k)
3479 lwsat(k)=soiliqwm(k)+qmin
3484 soiliqwm(k)=soilmoism(k)
3492 if(soilice(k).gt.0.) then
3493 smfrkeep(k)=soilice(k)
3495 smfrkeep(k)=soilmois(k)/riw
3499 !******************************************************************
3500 ! SOILPROP computes thermal diffusivity, and diffusional and
3501 ! hydraulic condeuctivities
3502 !******************************************************************
3503 CALL SOILPROP(spp_lsm,rstochcol,fieldcol_sf, &
3504 !--- input variables
3505 nzs,fwsat,lwsat,tav,keepfr, &
3506 soilmois,soiliqw,soilice, &
3507 soilmoism,soiliqwm,soilicem, &
3508 !--- soil fixed fields
3509 QWRTZ,rhocs,dqm,qmin,psis,bclh,ksat, &
3511 riw,xlmelt,CP,G0_P,cvw,ci, &
3513 !--- output variables
3514 thdif,diffu,hydro,cap)
3516 !********************************************************************
3517 !--- CALCULATION OF CANOPY WATER (Smirnova et al., 1996, EQ.16) AND DEW
3527 !--- If vegfrac.ne.0. then part of falling snow can be
3528 !--- intercepted by the canopy.
3532 EPOT = -FQ*(QVATM-QSG)
3534 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
3535 print *,'SNWE after subtracting intercepted snow - snwe=',snwe,vegfrac,cst
3539 ! check if all snow can evaporate during DT
3541 EPDT = EPOT * RAS *DELT*UMVEG
3542 IF(EPDT.gt.0. .and. SNWEPR.LE.EPDT) THEN
3543 BETA=SNWEPR/max(1.e-8,EPDT)
3547 WETCAN=min(0.25,max(0.,(CST/SAT))**CN)
3548 ! if(lai > 1.) wetcan=wetcan/lai
3551 !**************************************************************
3552 ! TRANSF computes transpiration function
3553 !**************************************************************
3555 !--- input variables
3556 nzs,nroot,soiliqw,tabs,lai,gswin, &
3557 !--- soil fixed fields
3558 dqm,qmin,ref,wilt,zshalf,pc,iland, &
3559 !--- output variables
3562 !--- Save soil temp and moisture from the beginning of time step
3565 smold(k)=soilmois(k)
3568 !**************************************************************
3569 ! SNOWTEMP solves heat budget and diffusion eqn. in soil
3570 !**************************************************************
3572 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
3573 print *, 'TSO before calling SNOWTEMP: ', tso
3576 !--- input variables
3578 delt,ktau,conflx,nzs,nddzs,nroot, &
3579 snwe,snwepr,snhei,newsnow,snowfrac, &
3580 beta,deltsn,snth,rhosn,rhonewsn,meltfactor, & ! add meltfactor
3582 PATM,TABS,QVATM,QCATM, &
3583 GLW,GSW,EMISS,RNET, &
3584 QKMS,TKMS,PC,rho,vegfrac, &
3585 thdif,cap,drycan,wetcan,cst, &
3586 tranf,transum,dew,mavail, &
3587 !--- soil fixed fields
3588 dqm,qmin,psis,bclh, &
3589 zsmain,zshalf,DTDZS,tbq, &
3591 xlvm,CP,rovcp,G0_P,cvw,stbolt, &
3592 !--- output variables
3593 snweprint,snheiprint,rsm, &
3594 tso,soilt,soilt1,tsnav,qvg,qsg,qcg, &
3595 smelt,snoh,snflx,s,ilnb,x)
3597 !************************************************************************
3598 !--- RECALCULATION OF DEW USING NEW VALUE OF QSG OR TRANSP IF NO DEW
3602 EPOT = -FQ*(QVATM-QSG)
3606 TRANSP(K)=vegfrac*RAS*FQ*(QVATM-QSG) &
3607 *tranf(K)*DRYCAN/zshalf(NROOT+1)
3608 ! IF(TRANSP(K).GT.0.) TRANSP(K)=0.
3624 !-- recalculating of frozen water in soil
3626 tln=log(tso(k)/273.15)
3628 soiliqw(k)=(dqm+qmin)*(XLMELT* &
3629 (tso(k)-273.15)/tso(k)/9.81/psis) &
3631 soiliqw(k)=max(0.,soiliqw(k))
3632 soiliqw(k)=min(soiliqw(k),soilmois(k))
3633 soilice(k)=(soilmois(k)-soiliqw(k))/riw
3634 !---- melting and freezing is balanced, soil ice cannot increase
3635 if(keepfr(k).eq.1.) then
3636 soilice(k)=min(soilice(k),smfrkeep(k))
3637 soiliqw(k)=max(0.,soilmois(k)-soilice(k)*riw)
3642 soiliqw(k)=soilmois(k)
3646 !*************************************************************************
3647 !--- TQCAN FOR SOLUTION OF MOISTURE BALANCE (Smirnova et al. 1996, EQ.22,28)
3648 ! AND TSO,ETA PROFILES
3649 !*************************************************************************
3652 delt,nzs,nddzs,DTDZS,DTDZS2,RIW, &
3653 zsmain,zshalf,diffu,hydro, &
3654 QSG,QVG,QCG,QCATM,QVATM,-INFWATER, &
3656 0.,SMELT,soilice,vegfrac, &
3659 DQM,QMIN,REF,KSAT,RAS,INFMAX, &
3661 SOILMOIS,SOILIQW,MAVAIL,RUNOFF1, &
3666 !-- Restore land-use parameters if all snow is melted
3667 IF(SNHEI.EQ.0.) then
3672 ! SNOM [mm] goes into the passed-in ACSNOM variable in the grid derived type
3673 SNOM=SNOM+SMELT*DELT*1.e3
3675 !--- KEEPFR is 1 when the temperature and moisture in soil
3676 !--- are both increasing. In this case soil ice should not
3677 !--- be increasing according to the freezing curve.
3678 !--- Some part of ice is melted, but additional water is
3679 !--- getting frozen. Thus, only structure of frozen soil is
3680 !--- changed, and phase changes are not affecting the heat
3681 !--- transfer. This situation may happen when it rains on the
3685 if (soilice(k).gt.0.) then
3686 if(tso(k).gt.told(k).and.soilmois(k).gt.smold(k)) then
3693 !--- THE DIAGNOSTICS OF SURFACE FLUXES
3695 T3 = STBOLT*SOILTold*SOILTold*SOILTold
3696 UPFLUX = T3 *0.5*(SOILTold+SOILT)
3697 XINET = EMISS*(GLW-UPFLUX)
3698 ! RNET = GSW + XINET
3699 HFX=-TKMS*CP*RHO*(TABS-SOILT) &
3700 *(P1000mb*0.00001/Patm)**ROVCP
3701 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
3702 print *,'potential temp HFX',hfx
3704 HFT=-TKMS*CP*RHO*(TABS-SOILT)
3705 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
3706 print *,'abs temp HFX',hft
3708 Q1 = - FQ*RAS* (QVATM - QSG)
3717 !-- moisture flux for coupling with MYJ PBL
3718 EETA=-QKMS*RAS*(QVATM/(1.+QVATM) - QSG/(1.+QSG))*1.E3
3719 CST= CST-EETA*DELT*vegfrac
3720 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
3721 print *,'MYJ EETA cond', EETA
3724 !-- actual moisture flux from RUC LSM
3725 DEW=QKMS*(QVATM-QSG)
3727 CST=CST+DELT*DEW*RAS * vegfrac
3728 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
3729 print *,'RUC LSM EETA cond',EETA
3736 EDIR1 = Q1*UMVEG *BETA
3738 EC1 = Q1 * WETCAN * vegfrac
3740 CST=max(0.,CST-EC1 * DELT)
3742 ! if(EC1 > CMC2MS) then
3743 ! EC1 = min(cmc2ms,ec1)
3747 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
3748 print*,'Q1,umveg,beta',Q1,umveg,beta
3749 print *,'wetcan,vegfrac',wetcan,vegfrac
3750 print *,'EC1,CMC2MS',EC1,CMC2MS
3754 !-- moisture flux for coupling with MYJ PBL
3755 EETA=-(QKMS*RAS*(QVATM/(1.+QVATM) - QSG/(1.+QSG))*1.E3)*BETA
3756 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
3757 print *,'MYJ EETA', EETA*XLVm,EETA
3760 ! to convert from m s-1 to kg m-2 s-1: *rho water=1.e3************
3761 !-- actual moisture flux from RUC LSM
3762 EETA = (EDIR1 + EC1 + ETT1)*1.E3
3763 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
3764 print *,'RUC LSM EETA',EETA*XLVm,EETA
3768 EETA = (EDIR1 + EC1 + ETT1)*1.E3
3774 FLTOT=RNET-HFT-XLVm*EETA-S-SNOH-x
3775 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
3776 print *,'SNOWSOIL - FLTOT,RNET,HFT,QFX,S,SNOH,X=',FLTOT,RNET,HFT,XLVm*EETA,s,SNOH,X
3777 print *,'edir1,ec1,ett1,mavail,qkms,qvatm,qvg,qsg,vegfrac,beta',&
3778 edir1,ec1,ett1,mavail,qkms,qvatm,qvg,qsg,vegfrac,beta
3783 1123 FORMAT(I5,8F12.3)
3784 1133 FORMAT(I7,8E12.4)
3785 123 format(i6,f6.2,7f8.1)
3786 122 FORMAT(1X,2I3,6F8.1,F8.3,F8.2)
3788 !-------------------------------------------------------------------
3789 END SUBROUTINE SNOWSOIL
3790 !-------------------------------------------------------------------
3792 SUBROUTINE SNOWSEAICE( &
3793 i,j,isoil,delt,ktau,conflx,nzs,nddzs, &
3794 meltfactor,rhonewsn,SNHEI_CRIT, & ! new
3795 ILAND,PRCPMS,RAINF,NEWSNOW,snhei,SNWE,snowfrac, &
3796 RHOSN,PATM,QVATM,QCATM, &
3797 GLW,GSW,EMISS,RNET, &
3798 QKMS,TKMS,RHO,myj, &
3799 !--- sea ice parameters
3801 tice,rhosice,capice,thdifice, &
3802 zsmain,zshalf,DTDZS,DTDZS2,tbq, &
3804 xlv,CP,rovcp,cw,stbolt,tabs, &
3805 !--- output variables
3806 ilnb,snweprint,snheiprint,rsm,tso, &
3807 dew,soilt,soilt1,tsnav,qvg,qsg,qcg, &
3808 SMELT,SNOH,SNFLX,SNOM,eeta, &
3809 qfx,hfx,s,sublim,prcpl,fltot &
3811 !***************************************************************
3812 ! Solving energy budget for snow on sea ice and heat diffusion
3813 ! eqns. in snow and sea ice
3814 !***************************************************************
3818 !-------------------------------------------------------------------
3819 !--- input variables
3821 INTEGER, INTENT(IN ) :: ktau,nzs , &
3822 nddzs !nddzs=2*(nzs-2)
3823 INTEGER, INTENT(IN ) :: i,j,isoil
3825 REAL, INTENT(IN ) :: DELT,CONFLX,PRCPMS , &
3826 RAINF,NEWSNOW,RHONEWSN, &
3827 meltfactor, snhei_crit
3830 LOGICAL, INTENT(IN ) :: myj
3831 !--- 3-D Atmospheric variables
3833 INTENT(IN ) :: PATM, &
3838 INTENT(IN ) :: GLW, &
3844 !--- sea ice properties
3845 REAL, DIMENSION(1:NZS) , &
3852 REAL, INTENT(IN ) :: &
3856 REAL, DIMENSION(1:NZS), INTENT(IN) :: ZSMAIN, &
3860 REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS
3862 REAL, DIMENSION(1:5001), INTENT(IN) :: TBQ
3864 !--- input/output variables
3865 !-------- 3-d soil moisture and temperature
3866 REAL, DIMENSION( 1:nzs ) , &
3867 INTENT(INOUT) :: TSO
3869 INTEGER, INTENT(INOUT) :: ILAND
3872 !-------- 2-d variables
3874 INTENT(INOUT) :: DEW, &
3899 INTEGER, INTENT(INOUT) :: ILNB
3901 REAL, INTENT(OUT) :: RSM, &
3904 !--- Local variables
3907 INTEGER :: nzs1,nzs2,k,k1,kn,kk
3908 REAL :: x,x1,x2,dzstop,ft,tn,denom
3910 REAL :: SNTH, NEWSN , &
3911 TABS, T3, UPFLUX, XINET , &
3912 BETA, SNWEPR,EPDT,PP
3913 REAL :: CP,rovcp,G0,LV,xlvm,STBOLT,xlmelt , &
3914 epot,fltot,fq,hft,q1,ras,rhoice,ci,cvw , &
3917 REAL :: rhocsn,thdifsn, &
3918 xsn,ddzsn,x1sn,d1sn,d2sn,d9sn,r22sn
3920 REAL :: cotsn,rhtsn,xsn1,ddzsn1,x1sn1,ftsnow,denomsn
3922 FKT,D1,D2,D9,D10,DID,R211,R21,R22,R6,R7,D11, &
3923 FKQ,R210,AA,BB,QS1,TS1,TQ2,TX2, &
3924 TDENOM,AA1,RHCS,H1,TSOB, SNPRIM, &
3925 SNODIF,SOH,TNOLD,QGOLD,SNOHGNEW
3926 REAL, DIMENSION(1:NZS) :: cotso,rhtso
3928 REAL :: RNET,rsmfrac,soiltfrac,hsn,icemelt,rr
3932 !-----------------------------------------------------------------
3934 !-- heat of sublimation of water vapor
3936 ! STBOLT=5.670151E-8
3938 !--- SNOW flag -- ISICE
3941 !--- DELTSN - is the threshold for splitting the snow layer into 2 layers.
3942 !--- With snow density 400 kg/m^3, this threshold is equal to 7.5 cm,
3943 !--- equivalent to 0.03 m SNWE. For other snow densities the threshold is
3944 !--- computed using SNWE=0.03 m and current snow density.
3945 !--- SNTH - the threshold below which the snow layer is combined with
3946 !--- the top sea ice layer. SNTH is computed using snwe=0.016 m, and
3947 !--- equals 4 cm for snow density 400 kg/m^3.
3949 ! increase thickness of top snow layer from 3 cm SWE to 5 cm SWE
3950 ! DELTSN=5.*SNHEI_CRIT
3951 ! snth=0.4*SNHEI_CRIT
3953 DELTSN=0.05*1.e3/rhosn
3954 snth=0.01*1.e3/rhosn
3955 ! snth=0.01601*1.e3/rhosn
3957 ! For 2-layer snow model when the snow depth is marginlly higher than DELTSN,
3958 ! reset DELTSN to half of snow depth.
3959 IF(SNHEI.GE.DELTSN+SNTH) THEN
3960 if(snhei-deltsn-snth.lt.snth) deltsn=0.5*(snhei-snth)
3961 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
3962 print *,'DELTSN ICE is changed,deltsn,snhei,snth', &
3963 i,j, deltsn,snhei,snth
3975 !18apr08 - add rhonewcsn
3976 RHOnewCSN=2090.* RHOnewSN
3977 THDIFSN = 0.265/RHOCSN
3998 DZSTOP=1./(ZSMAIN(2)-ZSMAIN(1))
4004 !*** DELTSN is the depth of the top layer of snow where
4005 !*** there is a temperature gradient, the rest of the snow layer
4006 !*** is considered to have constant temperature
4013 SNHEI=SNWE*1.e3/RHOSN
4016 ! check if all snow can evaporate during DT
4018 EPOT = -FQ*(QVATM-QSG)
4019 EPDT = EPOT * RAS *DELT
4020 IF(EPDT.GT.0. .and. SNWEPR.LE.EPDT) THEN
4021 BETA=SNWEPR/max(1.e-8,EPDT)
4025 !******************************************************************************
4026 ! COEFFICIENTS FOR THOMAS ALGORITHM FOR TSO
4027 !******************************************************************************
4034 X1=DTDZS(K1)*THDIFICE(KN-1)
4035 X2=DTDZS(K1+1)*THDIFICE(KN)
4036 FT=TSO(KN)+X1*(TSO(KN-1)-TSO(KN)) &
4037 -X2*(TSO(KN)-TSO(KN+1))
4038 DENOM=1.+X1+X2-X2*cotso(K)
4040 rhtso(K+1)=(FT+X2*rhtso(K))/DENOM
4042 !--- THE NZS element in COTSO and RHTSO will be for snow
4043 !--- There will be 2 layers in snow if it is deeper than DELTSN+SNTH
4044 IF(SNHEI.GE.SNTH) then
4045 if(snhei.le.DELTSN+SNTH) then
4046 !-- 1-layer snow model
4048 snprim=max(snth,snhei)
4051 XSN = DELT/2./(zshalf(2)+0.5*SNPRIM)
4052 DDZSN = XSN / SNPRIM
4053 X1SN = DDZSN * thdifsn
4054 X2 = DTDZS(1)*THDIFICE(1)
4055 FT = TSO(1)+X1SN*(SOILT-TSO(1)) &
4057 DENOM = 1. + X1SN + X2 -X2*cotso(NZS1)
4058 cotso(NZS)=X1SN/DENOM
4059 rhtso(NZS)=(FT+X2*rhtso(NZS1))/DENOM
4062 !*** Average temperature of snow pack (C)
4063 tsnav=0.5*(soilt+tso(1)) &
4067 !-- 2 layers in snow, SOILT1 is temperasture at DELTSN depth
4071 XSN = DELT/2./(0.5*SNHEI)
4072 XSN1= DELT/2./(zshalf(2)+0.5*(SNHEI-DELTSN))
4073 DDZSN = XSN / DELTSN
4074 DDZSN1 = XSN1 / (SNHEI-DELTSN)
4075 X1SN = DDZSN * thdifsn
4076 X1SN1 = DDZSN1 * thdifsn
4077 X2 = DTDZS(1)*THDIFICE(1)
4078 FT = TSO(1)+X1SN1*(SOILT1-TSO(1)) &
4080 DENOM = 1. + X1SN1 + X2 - X2*cotso(NZS1)
4081 cotso(nzs)=x1sn1/denom
4082 rhtso(nzs)=(ft+x2*rhtso(nzs1))/denom
4083 ftsnow = soilt1+x1sn*(soilt-soilt1) &
4084 -x1sn1*(soilt1-tso(1))
4085 denomsn = 1. + X1SN + X1SN1 - X1SN1*cotso(NZS)
4087 rhtsn=(ftsnow+X1SN1*rhtso(NZS))/denomsn
4088 !*** Average temperature of snow pack (C)
4089 tsnav=0.5/snhei*((soilt+soilt1)*deltsn &
4090 +(soilt1+tso(1))*(SNHEI-DELTSN)) &
4095 IF(SNHEI.LT.SNTH.AND.SNHEI.GT.0.) then
4096 !--- snow is too thin to be treated separately, therefore it
4097 !--- is combined with the first sea ice layer.
4098 snprim=SNHEI+zsmain(2)
4103 XSN = DELT/2./((zshalf(3)-zsmain(2))+0.5*snprim)
4105 X1SN = DDZSN * (fsn*thdifsn+fso*thdifice(1))
4106 X2=DTDZS(2)*THDIFICE(2)
4107 FT=TSO(2)+X1SN*(SOILT-TSO(2))- &
4109 denom = 1. + x1sn + x2 - x2*cotso(nzs-2)
4110 cotso(nzs1) = x1sn/denom
4111 rhtso(nzs1)=(FT+X2*rhtso(NZS-2))/denom
4112 tsnav=0.5*(soilt+tso(1)) &
4114 cotso(nzs)=cotso(NZS1)
4115 rhtso(nzs)=rhtso(nzs1)
4120 !************************************************************************
4121 !--- THE HEAT BALANCE EQUATION
4122 !18apr08 nmelt is the flag for melting, and SNOH is heat of snow phase changes
4126 EPOT=-QKMS*(QVATM-QSG)
4133 D9=THDIFICE(1)*RHCS*dzstop
4137 R22=.5/(THDIFICE(1)*DELT*dzstop**2)
4138 R6=EMISS *STBOLT*.5*TN**4
4142 IF(SNHEI.GE.SNTH) THEN
4143 if(snhei.le.DELTSN+SNTH) then
4152 D9SN= THDIFSN*RHOCSN / SNPRIM
4153 R22SN = SNPRIM*SNPRIM*0.5/(THDIFSN*DELT)
4156 IF(SNHEI.LT.SNTH.AND.SNHEI.GT.0.) then
4157 !--- thin snow is combined with sea ice
4160 D9SN = (fsn*THDIFSN*RHOCSN+fso*THDIFICE(1)*RHCS)/ &
4162 R22SN = snprim*snprim*0.5 &
4163 /((fsn*THDIFSN+fso*THDIFICE(1))*delt)
4167 !--- all snow is sublimated
4175 !---- TDENOM for snow
4176 TDENOM = D9SN*(1.-D1SN +R22SN)+D10+R21+R7 &
4178 +RHOnewCSN*NEWSNOW/DELT
4182 AA=XLVM*(BETA*FKQ+R210)/TDENOM
4183 BB=(D10*TABS+R21*TN+XLVM*(QVATM* &
4185 +R210*QVG)+D11+D9SN*(D2SN+R22SN*TN) &
4186 +RAINF*CVW*PRCPMS*max(273.15,TABS) &
4187 + RHOnewCSN*NEWSNOW/DELT*min(273.15,TABS) &
4192 !18apr08 - the iteration start point
4195 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
4196 print *,'VILKA-SNOW on SEAICE'
4197 print *,'tn,aa1,bb,pp,fkq,r210', &
4198 tn,aa1,bb,pp,fkq,r210
4199 print *,'TABS,QVATM,TN,QVG=',TABS,QVATM,TN,QVG
4202 CALL VILKA(TN,AA1,BB,PP,QS1,TS1,TBQ,KTAU,i,j,iland,isoil)
4203 !--- it is saturation over snow
4208 !--- SOILT - skin temperature
4211 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
4212 print *,' AFTER VILKA-SNOW on SEAICE'
4213 print *,' TS1,QS1: ', ts1,qs1
4215 ! Solution for temperature at 7.5 cm depth and snow-seaice interface
4216 IF(SNHEI.GE.SNTH) THEN
4217 if(snhei.gt.DELTSN+SNTH) then
4218 !-- 2-layer snow model
4219 SOILT1=min(273.15,rhtsn+cotsn*SOILT)
4220 TSO(1)=min(271.4,(rhtso(NZS)+cotso(NZS)*SOILT1))
4224 TSO(1)=min(271.4,(rhtso(NZS)+cotso(NZS)*SOILT))
4228 ELSEIF (SNHEI > 0. .and. SNHEI < SNTH) THEN
4230 TSO(2)=min(271.4,(rhtso(NZS1)+cotso(NZS1)*SOILT))
4231 tso(1)=min(271.4,(tso(2)+(soilt-tso(2))*fso))
4236 TSO(1)=min(271.4,SOILT)
4237 SOILT1=min(271.4,SOILT)
4240 !---- Final solution for TSO in sea ice
4241 IF (SNHEI > 0. .and. SNHEI < SNTH) THEN
4242 ! blended or snow is melted
4245 TSO(K)=min(271.4,rhtso(KK)+cotso(KK)*TSO(K-1))
4250 TSO(K)=min(271.4,rhtso(KK)+cotso(KK)*TSO(K-1))
4253 !--- For thin snow layer combined with the top soil layer
4254 !--- TSO(i,j,1) is computed by linear interpolation between SOILT
4256 ! if(SNHEI.LT.SNTH.AND.SNHEI.GT.0.)then
4257 ! tso(1)=min(271.4,tso(2)+(soilt-tso(2))*fso)
4262 if(nmelt.eq.1) go to 220
4264 !--- IF SOILT > 273.15 F then melting of snow can happen
4265 ! IF(SOILT.GT.273.15.AND.SNWE.GT.0.) THEN
4266 ! if all snow can evaporate, then there is nothing to melt
4267 IF(SOILT.GT.273.15.AND.SNWEPR-BETA*EPOT*RAS*DELT.GT.0..AND.SNHEI.GT.0.) THEN
4271 soiltfrac=snowfrac*273.15+(1.-snowfrac)*min(271.4,SOILT)
4273 QSG= QSN(soiltfrac,TBQ)/PP
4274 T3 = STBOLT*TNold*TNold*TNold
4275 UPFLUX = T3 * 0.5*(TNold+SOILTfrac)
4276 XINET = EMISS*(GLW-UPFLUX)
4277 ! RNET = GSW + XINET
4278 EPOT = -QKMS*(QVATM-QSG)
4289 EETA = Q1 * BETA *1.E3
4290 ! to convert from kg m-2 s-1 to m s-1: 1/rho water=1.e-3************
4294 HFX=D10*(TABS-soiltfrac)
4296 IF(SNHEI.GE.SNTH)then
4297 SOH=thdifsn*RHOCSN*(soiltfrac-TSOB)/SNPRIM
4300 SOH=(fsn*thdifsn*rhocsn+fso*thdifice(1)*rhcs)* &
4301 (soiltfrac-TSOB)/snprim
4304 X= (R21+D9SN*R22SN)*(soiltfrac-TNOLD) + &
4305 XLVM*R210*(QSG-QGOLD)
4306 !-- SNOH is energy flux of snow phase change
4307 SNOH=RNET+QFX +HFX &
4308 +RHOnewCSN*NEWSNOW/DELT*(min(273.15,TABS)-soiltfrac) &
4309 -SOH-X+RAINF*CVW*PRCPMS* &
4310 (max(273.15,TABS)-soiltfrac)
4312 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
4313 print *,'SNOWSEAICE melt I,J,SNOH,RNET,QFX,HFX,SOH,X',i,j,SNOH,RNET,QFX,HFX,SOH,X
4314 print *,'RHOnewCSN*NEWSNOW/DELT*(min(273.15,TABS)-soiltfrac)', &
4315 RHOnewCSN*NEWSNOW/DELT*(min(273.15,TABS)-soiltfrac)
4316 print *,'RAINF*CVW*PRCPMS*(max(273.15,TABS)-soiltfrac)', &
4317 RAINF*CVW*PRCPMS*(max(273.15,TABS)-soiltfrac)
4320 !-- SMELT is speed of melting in M/S
4321 SMELT= SNOH /XLMELT*1.E-3
4322 SMELT=AMIN1(SMELT,SNWEPR/DELT-BETA*EPOT*RAS)
4323 SMELT=AMAX1(0.,SMELT)
4325 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
4326 print *,'1-SMELT i,j',smelt,i,j
4328 !18apr08 - Egglston limit
4329 ! SMELT= amin1 (smelt, 5.6E-7*meltfactor*max(1.,(soilt-273.15)))
4330 SMELT= amin1 (smelt, 5.6E-8*meltfactor*max(1.,(soilt-273.15)))
4331 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
4332 print *,'2-SMELT i,j',smelt,i,j
4335 ! rr - potential melting
4336 rr=SNWEPR/delt-BETA*EPOT*RAS
4338 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
4339 print *,'3- SMELT i,j,smelt,rr',i,j,smelt,rr
4341 SNOHGNEW=SMELT*XLMELT*1.E3
4342 SNODIF=AMAX1(0.,(SNOH-SNOHGNEW))
4346 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
4347 print*,'soiltfrac,soilt,SNOHGNEW,SNODIF=', &
4348 i,j,soiltfrac,soilt,snohgnew,snodif
4349 print *,'SNOH,SNODIF',SNOH,SNODIF
4352 !*** From Koren et al. (1999) 13% of snow melt stays in the snow pack
4353 rsmfrac=min(0.18,(max(0.08,snwepr/0.10*0.13)))
4354 if(snhei > 0.01) then
4355 rsm=rsmfrac*smelt*delt
4357 ! do not keep melted water if snow depth is less that 1 cm
4360 !18apr08 rsm is part of melted water that stays in snow as liquid
4361 SMELT=AMAX1(0.,SMELT-rsm/delt)
4362 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
4363 print *,'4-SMELT i,j,smelt,rsm,snwepr,rsmfrac', &
4364 i,j,smelt,rsm,snwepr,rsmfrac
4367 !-- update liquid equivalent of snow depth
4368 !-- for evaporation and snow melt
4369 SNWE = AMAX1(0.,(SNWEPR- &
4370 (SMELT+BETA*EPOT*RAS)*DELT &
4371 ! (SMELT+BETA*EPOT*RAS)*DELT*snowfrac &
4375 !--- If there is no snow melting then just evaporation
4376 !--- or condensation changes SNWE
4378 if(snhei.ne.0.) then
4379 EPOT=-QKMS*(QVATM-QSG)
4380 SNWE = AMAX1(0.,(SNWEPR- &
4381 BETA*EPOT*RAS*DELT))
4382 ! BETA*EPOT*RAS*DELT*snowfrac))
4387 ! no iteration for snow on sea ice, because it will produce
4388 ! skin temperature higher than it is possible with snow on sea ice
4389 ! if(nmelt.eq.1) goto 212 ! second iteration
4392 if(smelt > 0..and. rsm > 0.) then
4393 if(snwe.le.rsm) then
4394 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
4395 print *,'SEAICE SNWE<RSM snwe,rsm,smelt*delt,epot*ras*delt,beta', &
4396 snwe,rsm,smelt*delt,epot*ras*delt,beta
4399 !*** Update snow density on effect of snow melt, melted
4400 !*** from the top of the snow. 13% of melted water
4401 !*** remains in the pack and changes its density.
4402 !*** Eq. 9 (with my correction) in Koren et al. (1999)
4404 xsn=(rhosn*(snwe-rsm)+1.e3*rsm)/ &
4406 rhosn=MIN(MAX(58.8,XSN),500.)
4407 !13mar18 rhosn=MIN(MAX(76.9,XSN),500.)
4410 thdifsn = 0.265/RHOCSN
4416 !--- if VEGFRAC.ne.0. then some snow stays on the canopy
4417 !--- and should be added to SNWE for water conservation
4418 ! 4 Nov 07 +VEGFRAC*cst
4419 snheiprint=snweprint*1.E3 / RHOSN
4421 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
4422 print *, 'snweprint : ',snweprint
4423 print *, 'D9SN,SOILT,TSOB : ', D9SN,SOILT,TSOB
4425 IF(SNHEI.GT.0.) THEN
4427 tsnav=0.5/snhei*((soilt+soilt1)*deltsn &
4428 +(soilt1+tso(1))*(SNHEI-DELTSN)) &
4431 tsnav=0.5*(soilt+tso(1)) - 273.15
4434 !--- RECALCULATION OF DEW USING NEW VALUE OF QSG
4437 QSG= QSN(SOILT,TBQ)/PP
4438 EPOT = -FQ*(QVATM-QSG)
4444 SNOM=SNOM+SMELT*DELT*1.e3
4446 !--- THE DIAGNOSTICS OF SURFACE FLUXES
4448 T3 = STBOLT*TNold*TNold*TNold
4449 UPFLUX = T3 *0.5*(SOILT+TNold)
4450 XINET = EMISS*(GLW-UPFLUX)
4451 ! RNET = GSW + XINET
4452 HFT=-TKMS*CP*RHO*(TABS-SOILT)
4453 HFX=-TKMS*CP*RHO*(TABS-SOILT) &
4454 *(P1000mb*0.00001/Patm)**ROVCP
4455 Q1 = - FQ*RAS* (QVATM - QSG)
4459 !-- moisture flux for coupling with MYJ PBL
4460 EETA=-QKMS*RAS*(QVATM/(1.+QVATM) - QSG/(1.+QSG))*1.E3
4462 !-- actual moisture flux from RUC LSM
4463 DEW=QKMS*(QVATM-QSG)
4472 !-- moisture flux for coupling with MYJ PBL
4473 EETA=-QKMS*RAS*BETA*(QVATM/(1.+QVATM) - QVG/(1.+QVG))*1.E3
4475 ! to convert from m s-1 to kg m-2 s-1: *rho water=1.e3************
4476 !-- actual moisture flux from RUC LSM
4485 IF(SNHEI.GE.SNTH)then
4486 S=thdifsn*RHOCSN*(soilt-TSOB)/SNPRIM
4488 ELSEIF(SNHEI.lt.SNTH.and.SNHEI.GT.0.) then
4489 S=(fsn*thdifsn*rhocsn+fso*thdifice(1)*rhcs)* &
4492 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
4493 print *,'SNOW is thin, snflx',i,j,snflx
4496 SNFLX=D9SN*(SOILT-TSOB)
4497 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
4498 print *,'SNOW is GONE, snflx',i,j,snflx
4502 SNHEI=SNWE *1.E3 / RHOSN
4504 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
4505 print *,'SNHEI,SNOH',i,j,SNHEI,SNOH
4508 X= (R21+D9SN*R22SN)*(soilt-TNOLD) + &
4509 XLVM*R210*(QSG-QGOLD)
4510 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
4511 print *,'SNOWSEAICE storage ',i,j,x
4512 print *,'R21,D9sn,r22sn,soiltfrac,tnold,qsg,qgold,snprim', &
4513 R21,D9sn,r22sn,soiltfrac,tnold,qsg,qgold,snprim
4516 -RHOnewCSN*NEWSNOW/DELT*(min(273.15,TABS)-SOILT) &
4517 -RAINF*CVW*PRCPMS*(max(273.15,TABS)-SOILT)
4519 ! -- excess energy is spent on ice melt
4520 icemelt = RNET-HFT-XLVm*EETA-S-SNOH-X
4521 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
4522 print *,'SNOWSEAICE icemelt=',icemelt
4525 FLTOT=RNET-HFT-XLVm*EETA-S-SNOH-x-icemelt
4526 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
4527 print *,'i,j,snhei,qsg,soilt,soilt1,tso,TABS,QVATM', &
4528 i,j,snhei,qsg,soilt,soilt1,tso,tabs,qvatm
4529 print *,'SNOWSEAICE - FLTOT,RNET,HFT,QFX,S,SNOH,icemelt,snodif,X,SOILT=' &
4530 ,FLTOT,RNET,HFT,XLVm*EETA,s,SNOH,icemelt,snodif,X,SOILT
4532 !-- Restore sea-ice parameters if snow is less than threshold
4533 IF(SNHEI.EQ.0.) then
4540 !------------------------------------------------------------------------
4541 !------------------------------------------------------------------------
4542 END SUBROUTINE SNOWSEAICE
4543 !------------------------------------------------------------------------
4546 SUBROUTINE SOILTEMP( &
4547 !--- input variables
4549 delt,ktau,conflx,nzs,nddzs,nroot, &
4550 PRCPMS,RAINF,PATM,TABS,QVATM,QCATM, &
4552 QKMS,TKMS,PC,RHO,VEGFRAC,lai, &
4553 THDIF,CAP,DRYCAN,WETCAN, &
4554 TRANSUM,DEW,MAVAIL,soilres,alfa, &
4555 !--- soil fixed fields
4557 ZSMAIN,ZSHALF,DTDZS,TBQ, &
4559 XLV,CP,G0_P,CVW,STBOLT, &
4560 !--- output variables
4561 TSO,SOILT,QVG,QSG,QCG,X)
4563 !*************************************************************
4564 ! Energy budget equation and heat diffusion eqn are
4567 ! DELT - time step (s)
4568 ! ktau - numver of time step
4569 ! CONFLX - depth of constant flux layer (m)
4570 ! IME, JME, KME, NZS - dimensions of the domain
4571 ! NROOT - number of levels within the root zone
4572 ! PRCPMS - precipitation rate in m/s
4573 ! COTSO, RHTSO - coefficients for implicit solution of
4574 ! heat diffusion equation
4575 ! THDIF - thermal diffusivity (m^2/s)
4576 ! QSG,QVG,QCG - saturated mixing ratio, mixing ratio of
4577 ! water vapor and cloud at the ground
4578 ! surface, respectively (kg/kg)
4579 ! PATM - pressure [bar]
4580 ! QC3D,QV3D - cloud and water vapor mixing ratio
4581 ! at the first atm. level (kg/kg)
4582 ! EMISS,RNET - emissivity (0-1) of the ground surface and net
4583 ! radiation at the surface (W/m^2)
4584 ! QKMS - exchange coefficient for water vapor in the
4585 ! surface layer (m/s)
4586 ! TKMS - exchange coefficient for heat in the surface
4588 ! PC - plant coefficient (resistance)
4589 ! RHO - density of atmosphere near surface (kg/m^3)
4590 ! VEGFRAC - greeness fraction (0-1)
4591 ! CAP - volumetric heat capacity (J/m^3/K)
4592 ! DRYCAN - dry fraction of vegetated area where
4593 ! transpiration may take place (0-1)
4594 ! WETCAN - fraction of vegetated area covered by canopy
4596 ! TRANSUM - transpiration function integrated over the
4598 ! DEW - dew in kg/m^2s
4599 ! MAVAIL - fraction of maximum soil moisture in the top
4601 ! ZSMAIN - main levels in soil (m)
4602 ! ZSHALF - middle of the soil layers (m)
4603 ! DTDZS - dt/(2.*dzshalf*dzmain)
4604 ! TBQ - table to define saturated mixing ration
4605 ! of water vapor for given temperature and pressure
4606 ! TSO - soil temperature (K)
4607 ! SOILT - skin temperature (K)
4609 !****************************************************************
4612 !-----------------------------------------------------------------
4614 !--- input variables
4616 INTEGER, INTENT(IN ) :: nroot,ktau,nzs , &
4617 nddzs !nddzs=2*(nzs-2)
4618 INTEGER, INTENT(IN ) :: i,j,iland,isoil
4619 REAL, INTENT(IN ) :: DELT,CONFLX,PRCPMS, RAINF
4620 REAL, INTENT(INOUT) :: DRYCAN,WETCAN,TRANSUM
4621 !--- 3-D Atmospheric variables
4623 INTENT(IN ) :: PATM, &
4639 !--- soil properties
4650 REAL, INTENT(IN ) :: CP, &
4658 REAL, DIMENSION(1:NZS), INTENT(IN) :: ZSMAIN, &
4663 REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS
4665 REAL, DIMENSION(1:5001), INTENT(IN) :: TBQ
4668 !--- input/output variables
4669 !-------- 3-d soil moisture and temperature
4670 REAL, DIMENSION( 1:nzs ) , &
4671 INTENT(INOUT) :: TSO
4673 !-------- 2-d variables
4683 !--- Local variables
4685 REAL :: x,x1,x2,x4,dzstop,can,ft,sph , &
4686 tn,trans,umveg,denom,fex
4688 REAL :: FKT,D1,D2,D9,D10,DID,R211,R21,R22,R6,R7,D11 , &
4689 PI,H,FKQ,R210,AA,BB,PP,Q1,QS1,TS1,TQ2,TX2 , &
4692 REAL :: C,CC,AA1,RHCS,H1, QGOLD
4694 REAL, DIMENSION(1:NZS) :: cotso,rhtso
4696 INTEGER :: nzs1,nzs2,k,k1,kn,kk, iter
4699 !-----------------------------------------------------------------
4705 dzstop=1./(ZSMAIN(2)-ZSMAIN(1))
4713 !******************************************************************************
4714 ! COEFFICIENTS FOR THOMAS ALGORITHM FOR TSO
4715 !******************************************************************************
4716 ! did=2.*(ZSMAIN(nzs)-ZSHALF(nzs))
4717 ! h1=DTDZS(8)*THDIF(nzs-1)*(ZSHALF(nzs)-ZSHALF(nzs-1))/did
4718 ! cotso(1)=h1/(1.+h1)
4719 ! rhtso(1)=(tso(nzs)+h1*(tso(nzs-1)-tso(nzs)))/
4726 X1=DTDZS(K1)*THDIF(KN-1)
4727 X2=DTDZS(K1+1)*THDIF(KN)
4728 FT=TSO(KN)+X1*(TSO(KN-1)-TSO(KN)) &
4729 -X2*(TSO(KN)-TSO(KN+1))
4730 DENOM=1.+X1+X2-X2*cotso(K)
4732 rhtso(K+1)=(FT+X2*rhtso(K))/DENOM
4735 !************************************************************************
4736 !--- THE HEAT BALANCE EQUATION (Smirnova et al., 1996, EQ. 21,26)
4742 TRANS=TRANSUM*DRYCAN/ZSHALF(NROOT+1)
4744 UMVEG=(1.-VEGFRAC) * soilres
4750 D9=THDIF(1)*RHCS*dzstop
4754 R22=.5/(THDIF(1)*DELT*dzstop**2)
4755 R6=EMISS *STBOLT*.5*TN**4
4758 TDENOM=D9*(1.-D1+R22)+D10+R21+R7 &
4764 AA=XLV*(FKQ*UMVEG+R210)/TDENOM
4765 BB=(D10*TABS+R21*TN+XLV*(QVATM* &
4767 +R210*QVG)+D11+D9*(D2+R22*TN) &
4768 +RAINF*CVW*PRCPMS*max(273.15,TABS) &
4774 CALL VILKA(TN,AA1,BB,PP,QS1,TS1,TBQ,KTAU,i,j,iland,isoil)
4779 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
4780 ! if (i==421.and.j==280) then
4781 print *,'VILKA1 - TS1,QS1,TQ2,H,TX2,Q1',TS1,QS1,TQ2,H,TX2,Q1
4783 !with alfa go to 100
4784 IF(Q1.LT.QS1) GOTO 100
4785 !--- if no saturation - goto 100
4786 !--- if saturation - goto 90
4791 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
4792 ! if (i==421.and.j==280) then
4793 print *,'90 QVG,QSG,QCG,TSO(1)',QVG,QSG,QCG,TSO(1)
4798 CALL VILKA(TN,AA,BB,PP,QS1,TS1,TBQ,KTAU,i,j,iland,isoil)
4800 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
4801 ! if(i.eq.279.and.j.eq.263) then
4802 ! if (i==421.and.j==280) then
4803 print *,'VILKA2 - TS1,QS1,TQ2,H,TX2,Q1',TS1,QS1,TQ2,H,TX2,Q1
4805 IF(Q1.GE.QS1) GOTO 90
4806 !with alfa 100 continue
4809 ! if( QS1>QVATM .and. QVATM > QVG) then
4811 ! print *,'very dry soils mavail,qvg,qs1,qvatm,ts1',i,j,mavail,qvg,qs1,qvatm,ts1
4816 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
4817 ! if (i==421.and.j==280) then
4818 print *,'q1,qsg,qvg,qvatm,alfa,h',q1,qsg,qvg,qvatm,alfa,h
4821 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
4822 print *,'200 QVG,QSG,QCG,TSO(1)',QVG,QSG,QCG,TSO(1)
4825 !--- SOILT - skin temperature
4828 !---- Final solution for soil temperature - TSO
4831 TSO(K)=rhtso(KK)+cotso(KK)*TSO(K-1)
4834 X= (cp*rho*r211+rhcs*zsmain(2)*0.5/delt)*(SOILT-TN) + &
4835 XLV*rho*r211*(QVG-QGOLD)
4837 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
4838 print*,'SOILTEMP storage, i,j,x,soilt,tn,qvg,qvgold', &
4839 i,j,x,soilt,tn,qvg,qgold
4840 print *,'TEMP term (cp*rho*r211+rhcs*zsmain(2)*0.5/delt)*(SOILT-TN)',&
4841 (cp*rho*r211+rhcs*zsmain(2)*0.5/delt)*(SOILT-TN)
4842 print *,'QV term XLV*rho*r211*(QVG-QGOLD)',XLV*rho*r211*(QVG-QGOLD)
4846 -RAINF*CVW*PRCPMS*(max(273.15,TABS)-SOILT)
4848 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
4852 !--------------------------------------------------------------------
4853 END SUBROUTINE SOILTEMP
4854 !--------------------------------------------------------------------
4857 SUBROUTINE SNOWTEMP( &
4858 !--- input variables
4860 delt,ktau,conflx,nzs,nddzs,nroot, &
4861 snwe,snwepr,snhei,newsnow,snowfrac, &
4862 beta,deltsn,snth,rhosn,rhonewsn,meltfactor, & ! add meltfactor
4864 PATM,TABS,QVATM,QCATM, &
4865 GLW,GSW,EMISS,RNET, &
4866 QKMS,TKMS,PC,RHO,VEGFRAC, &
4867 THDIF,CAP,DRYCAN,WETCAN,CST, &
4868 TRANF,TRANSUM,DEW,MAVAIL, &
4869 !--- soil fixed fields
4870 DQM,QMIN,PSIS,BCLH, &
4871 ZSMAIN,ZSHALF,DTDZS,TBQ, &
4873 XLVM,CP,rovcp,G0_P,CVW,STBOLT, &
4874 !--- output variables
4875 SNWEPRINT,SNHEIPRINT,RSM, &
4876 TSO,SOILT,SOILT1,TSNAV,QVG,QSG,QCG, &
4877 SMELT,SNOH,SNFLX,S,ILNB,X)
4879 !********************************************************************
4880 ! Energy budget equation and heat diffusion eqn are
4881 ! solved here to obtain snow and soil temperatures
4883 ! DELT - time step (s)
4884 ! ktau - numver of time step
4885 ! CONFLX - depth of constant flux layer (m)
4886 ! IME, JME, KME, NZS - dimensions of the domain
4887 ! NROOT - number of levels within the root zone
4888 ! PRCPMS - precipitation rate in m/s
4889 ! COTSO, RHTSO - coefficients for implicit solution of
4890 ! heat diffusion equation
4891 ! THDIF - thermal diffusivity (W/m/K)
4892 ! QSG,QVG,QCG - saturated mixing ratio, mixing ratio of
4893 ! water vapor and cloud at the ground
4894 ! surface, respectively (kg/kg)
4895 ! PATM - pressure [bar]
4896 ! QCATM,QVATM - cloud and water vapor mixing ratio
4897 ! at the first atm. level (kg/kg)
4898 ! EMISS,RNET - emissivity (0-1) of the ground surface and net
4899 ! radiation at the surface (W/m^2)
4900 ! QKMS - exchange coefficient for water vapor in the
4901 ! surface layer (m/s)
4902 ! TKMS - exchange coefficient for heat in the surface
4904 ! PC - plant coefficient (resistance)
4905 ! RHO - density of atmosphere near surface (kg/m^3)
4906 ! VEGFRAC - greeness fraction (0-1)
4907 ! CAP - volumetric heat capacity (J/m^3/K)
4908 ! DRYCAN - dry fraction of vegetated area where
4909 ! transpiration may take place (0-1)
4910 ! WETCAN - fraction of vegetated area covered by canopy
4912 ! TRANSUM - transpiration function integrated over the
4914 ! DEW - dew in kg/m^2/s
4915 ! MAVAIL - fraction of maximum soil moisture in the top
4917 ! ZSMAIN - main levels in soil (m)
4918 ! ZSHALF - middle of the soil layers (m)
4919 ! DTDZS - dt/(2.*dzshalf*dzmain)
4920 ! TBQ - table to define saturated mixing ration
4921 ! of water vapor for given temperature and pressure
4922 ! TSO - soil temperature (K)
4923 ! SOILT - skin temperature (K)
4925 !*********************************************************************
4928 !---------------------------------------------------------------------
4929 !--- input variables
4931 INTEGER, INTENT(IN ) :: nroot,ktau,nzs , &
4932 nddzs !nddzs=2*(nzs-2)
4934 INTEGER, INTENT(IN ) :: i,j,iland,isoil
4935 REAL, INTENT(IN ) :: DELT,CONFLX,PRCPMS , &
4936 RAINF,NEWSNOW,DELTSN,SNTH , &
4937 TABS,TRANSUM,SNWEPR , &
4941 !--- 3-D Atmospheric variables
4943 INTENT(IN ) :: PATM, &
4948 INTENT(IN ) :: GLW, &
4956 !--- soil properties
4964 REAL, INTENT(IN ) :: CP, &
4972 REAL, DIMENSION(1:NZS), INTENT(IN) :: ZSMAIN, &
4978 REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS
4980 REAL, DIMENSION(1:5001), INTENT(IN) :: TBQ
4983 !--- input/output variables
4984 !-------- 3-d soil moisture and temperature
4985 REAL, DIMENSION( 1:nzs ) , &
4986 INTENT(INOUT) :: TSO
4989 !-------- 2-d variables
4991 INTENT(INOUT) :: DEW, &
5010 REAL, INTENT(INOUT) :: DRYCAN, WETCAN
5012 REAL, INTENT(OUT) :: RSM, &
5015 INTEGER, INTENT(OUT) :: ilnb
5016 !--- Local variables
5019 INTEGER :: nzs1,nzs2,k,k1,kn,kk
5021 REAL :: x,x1,x2,x4,dzstop,can,ft,sph, &
5022 tn,trans,umveg,denom
5024 REAL :: cotsn,rhtsn,xsn1,ddzsn1,x1sn1,ftsnow,denomsn
5026 REAL :: t3,upflux,xinet,ras, &
5027 xlmelt,rhocsn,thdifsn, &
5028 beta,epot,xsn,ddzsn,x1sn,d1sn,d2sn,d9sn,r22sn
5031 FKT,D1,D2,D9,D10,DID,R211,R21,R22,R6,R7,D11, &
5032 PI,H,FKQ,R210,AA,BB,PP,Q1,QS1,TS1,TQ2,TX2, &
5033 TDENOM,C,CC,AA1,RHCS,H1, &
5034 tsob, snprim, sh1, sh2, &
5035 smeltg,snohg,snodif,soh, &
5036 CMC2MS,TNOLD,QGOLD,SNOHGNEW
5038 REAL, DIMENSION(1:NZS) :: transp,cotso,rhtso
5046 REAL :: RNET,rsmfrac,soiltfrac,hsn,rr
5047 integer :: nmelt, iter
5049 !-----------------------------------------------------------------
5059 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
5060 print *, 'SNOWTEMP: SNHEI,SNTH,SOILT1: ',SNHEI,SNTH,SOILT1,soilt
5064 !18apr08 - add rhonewcsn
5065 RHOnewCSN=2090.* RHOnewSN
5066 THDIFSN = 0.265/RHOCSN
5086 DZSTOP=1./(ZSMAIN(2)-ZSMAIN(1))
5088 !******************************************************************************
5089 ! COEFFICIENTS FOR THOMAS ALGORITHM FOR TSO
5090 !******************************************************************************
5091 ! did=2.*(ZSMAIN(nzs)-ZSHALF(nzs))
5092 ! h1=DTDZS(8)*THDIF(nzs-1)*(ZSHALF(nzs)-ZSHALF(nzs-1))/did
5093 ! cotso(1)=h1/(1.+h1)
5094 ! rhtso(1)=(tso(nzs)+h1*(tso(nzs-1)-tso(nzs)))/
5102 X1=DTDZS(K1)*THDIF(KN-1)
5103 X2=DTDZS(K1+1)*THDIF(KN)
5104 FT=TSO(KN)+X1*(TSO(KN-1)-TSO(KN)) &
5105 -X2*(TSO(KN)-TSO(KN+1))
5106 DENOM=1.+X1+X2-X2*cotso(K)
5108 rhtso(K+1)=(FT+X2*rhtso(K))/DENOM
5110 !--- THE NZS element in COTSO and RHTSO will be for snow
5111 !--- There will be 2 layers in snow if it is deeper than DELTSN+SNTH
5112 IF(SNHEI.GE.SNTH) then
5113 if(snhei.le.DELTSN+SNTH) then
5114 !-- 1-layer snow model
5115 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
5116 print *,'1-layer - snth,snhei,deltsn',snth,snhei,deltsn
5119 snprim=max(snth,snhei)
5122 XSN = DELT/2./(zshalf(2)+0.5*SNPRIM)
5123 DDZSN = XSN / SNPRIM
5124 X1SN = DDZSN * thdifsn
5125 X2 = DTDZS(1)*THDIF(1)
5126 FT = TSO(1)+X1SN*(SOILT-TSO(1)) &
5128 DENOM = 1. + X1SN + X2 -X2*cotso(NZS1)
5129 cotso(NZS)=X1SN/DENOM
5130 rhtso(NZS)=(FT+X2*rhtso(NZS1))/DENOM
5133 !*** Average temperature of snow pack (C)
5134 tsnav=0.5*(soilt+tso(1)) &
5138 !-- 2 layers in snow, SOILT1 is temperasture at DELTSN depth
5139 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
5140 print *,'2-layer - snth,snhei,deltsn',snth,snhei,deltsn
5145 XSN = DELT/2./(0.5*deltsn)
5146 XSN1= DELT/2./(zshalf(2)+0.5*(SNHEI-DELTSN))
5147 DDZSN = XSN / DELTSN
5148 DDZSN1 = XSN1 / (SNHEI-DELTSN)
5149 X1SN = DDZSN * thdifsn
5150 X1SN1 = DDZSN1 * thdifsn
5151 X2 = DTDZS(1)*THDIF(1)
5152 FT = TSO(1)+X1SN1*(SOILT1-TSO(1)) &
5154 DENOM = 1. + X1SN1 + X2 - X2*cotso(NZS1)
5155 cotso(nzs)=x1sn1/denom
5156 rhtso(nzs)=(ft+x2*rhtso(nzs1))/denom
5157 ftsnow = soilt1+x1sn*(soilt-soilt1) &
5158 -x1sn1*(soilt1-tso(1))
5159 denomsn = 1. + X1SN + X1SN1 - X1SN1*cotso(NZS)
5161 rhtsn=(ftsnow+X1SN1*rhtso(NZS))/denomsn
5162 !*** Average temperature of snow pack (C)
5163 tsnav=0.5/snhei*((soilt+soilt1)*deltsn &
5164 +(soilt1+tso(1))*(SNHEI-DELTSN)) &
5168 IF(SNHEI.LT.SNTH.AND.SNHEI.GT.0.) then
5169 ! IF(SNHEI.LT.SNTH.AND.SNHEI.GE.0.) then
5170 !--- snow is too thin to be treated separately, therefore it
5171 !--- is combined with the first soil layer.
5172 snprim=SNHEI+zsmain(2)
5177 XSN = DELT/2./((zshalf(3)-zsmain(2))+0.5*snprim)
5179 X1SN = DDZSN * (fsn*thdifsn+fso*thdif(1))
5180 X2=DTDZS(2)*THDIF(2)
5181 FT=TSO(2)+X1SN*(SOILT-TSO(2))- &
5183 denom = 1. + x1sn + x2 - x2*cotso(nzs-2)
5184 cotso(nzs1) = x1sn/denom
5185 rhtso(nzs1)=(FT+X2*rhtso(NZS-2))/denom
5186 tsnav=0.5*(soilt+tso(1)) &
5188 cotso(NZS)=cotso(nzs1)
5189 rhtso(NZS)=rhtso(nzs1)
5195 !************************************************************************
5196 !--- THE HEAT BALANCE EQUATION (Smirnova et al. 1996, EQ. 21,26)
5197 !18apr08 nmelt is the flag for melting, and SNOH is heat of snow phase changes
5202 EPOT=-QKMS*(QVATM-QGOLD)
5205 TRANS=TRANSUM*DRYCAN/ZSHALF(NROOT+1)
5212 D9=THDIF(1)*RHCS*dzstop
5216 R22=.5/(THDIF(1)*DELT*dzstop**2)
5217 R6=EMISS *STBOLT*.5*TN**4
5221 IF(SNHEI.GE.SNTH) THEN
5222 if(snhei.le.DELTSN+SNTH) then
5226 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
5227 print *,'1 layer d1sn,d2sn',i,j,d1sn,d2sn
5233 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
5234 print *,'2 layers d1sn,d2sn',i,j,d1sn,d2sn
5237 D9SN= THDIFSN*RHOCSN / SNPRIM
5238 R22SN = SNPRIM*SNPRIM*0.5/(THDIFSN*DELT)
5239 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
5240 print *,'1 or 2 layers D9sn,R22sn',d9sn,r22sn
5244 IF(SNHEI.LT.SNTH.AND.SNHEI.GT.0.) then
5245 !--- thin snow is combined with soil
5248 D9SN = (fsn*THDIFSN*RHOCSN+fso*THDIF(1)*RHCS)/ &
5250 R22SN = snprim*snprim*0.5 &
5251 /((fsn*THDIFSN+fso*THDIF(1))*delt)
5252 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
5253 print *,' Combined D9SN,R22SN,D1SN,D2SN: ',D9SN,R22SN,D1SN,D2SN
5257 !--- all snow is sublimated
5262 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
5263 print *,' SNHEI = 0, D9SN,R22SN,D1SN,D2SN: ',D9SN,R22SN,D1SN,D2SN
5269 !18apr08 - the snow melt iteration start point
5272 !---- TDENOM for snow
5273 TDENOM = D9SN*(1.-D1SN +R22SN)+D10+R21+R7 &
5275 +RHOnewCSN*NEWSNOW/DELT
5281 AA=XLVM*(BETA*FKQ*UMVEG+R210)/TDENOM
5282 BB=(D10*TABS+R21*TN+XLVM*(QVATM* &
5283 (BETA*FKQ*UMVEG+C) &
5284 +R210*QGOLD)+D11+D9SN*(D2SN+R22SN*TN) &
5285 +RAINF*CVW*PRCPMS*max(273.15,TABS) &
5286 + RHOnewCSN*NEWSNOW/DELT*min(273.15,TABS) &
5293 CALL VILKA(TN,AA1,BB,PP,QS1,TS1,TBQ,KTAU,i,j,iland,isoil)
5297 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
5298 print *,'VILKA1 - TS1,QS1,TQ2,H,TX2,Q1',TS1,QS1,TQ2,H,TX2,Q1
5300 IF(Q1.LT.QS1) GOTO 100
5301 !--- if no saturation - goto 100
5302 !--- if saturation - goto 90
5306 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
5307 print *,'90 QVG,QSG,QCG,TSO(1)',QVG,QSG,QCG,TSO(1)
5312 CALL VILKA(TN,AA,BB,PP,QS1,TS1,TBQ,KTAU,i,j,iland,isoil)
5314 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
5315 print *,'VILKA2 - TS1,QS1,H,TX2,Q1',TS1,QS1,TQ2,H,TX2,Q1
5317 IF(Q1.GT.QS1) GOTO 90
5321 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
5322 print *,'No Saturation QVG,QSG,QCG,TSO(1)',QVG,QSG,QCG,TSO(1)
5326 !--- SOILT - skin temperature
5328 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
5329 IF(i.eq.266.and.j.eq.447) then
5330 print *,'snwe,snhei,soilt,soilt1,tso',i,j,snwe,snhei,soilt,soilt1,tso
5333 ! Solution for temperature at 7.5 cm depth and snow-soil interface
5334 IF(SNHEI.GE.SNTH) THEN
5335 if(snhei.gt.DELTSN+SNTH) then
5336 !-- 2-layer snow model
5337 SOILT1=min(273.15,rhtsn+cotsn*SOILT)
5338 TSO(1)=rhtso(NZS)+cotso(NZS)*SOILT1
5342 TSO(1)=rhtso(NZS)+cotso(NZS)*SOILT
5346 ELSEIF (SNHEI > 0. .and. SNHEI < SNTH) THEN
5348 TSO(2)=rhtso(NZS1)+cotso(NZS1)*SOILT
5349 tso(1)=(tso(2)+(soilt-tso(2))*fso)
5353 !-- very thin or zero snow. If snow is thin we suppose that
5354 !--- tso(i,j,1)=SOILT, and later we recompute tso(i,j,1)
5361 !---- Final solution for TSO
5362 IF (SNHEI > 0. .and. SNHEI < SNTH) THEN
5363 ! blended or snow is melted
5366 TSO(K)=rhtso(KK)+cotso(KK)*TSO(K-1)
5372 TSO(K)=rhtso(KK)+cotso(KK)*TSO(K-1)
5375 !--- For thin snow layer combined with the top soil layer
5376 !--- TSO(1) is recomputed by linear interpolation between SOILT
5378 ! if(SNHEI.LT.SNTH.AND.SNHEI.GT.0.)then
5379 ! tso(1)=tso(2)+(soilt-tso(2))*fso
5385 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
5386 ! IF(i.eq.266.and.j.eq.447) then
5387 print *,'SOILT,SOILT1,tso,TSOB,QSG',i,j,SOILT,SOILT1,tso,TSOB,QSG,'nmelt=',nmelt
5390 if(nmelt.eq.1) go to 220
5392 !--- IF SOILT > 273.15 F then melting of snow can happen
5393 ! IF(SOILT.GT.273.15.AND.SNHEI.GT.0.) THEN
5394 ! if all snow can evaporate, then there is nothing to melt
5395 IF(SOILT.GT.273.15.AND.SNWEPR-BETA*EPOT*RAS*DELT.GT.0.AND.SNHEI.GT.0.) THEN
5397 soiltfrac=snowfrac*273.15+(1.-snowfrac)*SOILT
5398 QSG=min(QSG, QSN(soiltfrac,TBQ)/PP)
5400 T3 = STBOLT*TN*TN*TN
5401 UPFLUX = T3 * 0.5*(TN + SOILTfrac)
5402 XINET = EMISS*(GLW-UPFLUX)
5403 ! RNET = GSW + XINET
5404 EPOT = -QKMS*(QVATM-QSG)
5407 IF (Q1.LE.0..or.iter==1) THEN
5419 TRANSP(K)=-VEGFRAC*q1 &
5420 *TRANF(K)*DRYCAN/zshalf(NROOT+1)
5421 ! IF(TRANSP(K).GT.0.) TRANSP(K)=0.
5428 EDIR1 = Q1*UMVEG * BETA
5429 EC1 = Q1 * WETCAN * vegfrac
5431 ! EC1=MIN(CMC2MS,EC1)
5432 EETA = (EDIR1 + EC1 + ETT1)*1.E3
5433 ! to convert from kg m-2 s-1 to m s-1: 1/rho water=1.e-3************
5437 HFX=-D10*(TABS-soiltfrac)
5439 IF(SNHEI.GE.SNTH)then
5440 SOH=thdifsn*RHOCSN*(soiltfrac-TSOB)/SNPRIM
5443 SOH=(fsn*thdifsn*rhocsn+fso*thdif(1)*rhcs)* &
5444 (soiltfrac-TSOB)/snprim
5449 X= (R21+D9SN*R22SN)*(soiltfrac-TN) + &
5450 XLVM*R210*(QVG-QGOLD)
5451 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
5452 print *,'SNOWTEMP storage ',i,j,x
5453 print *,'R21,D9sn,r22sn,soiltfrac,tn,qsg,qvg,qgold,snprim', &
5454 R21,D9sn,r22sn,soiltfrac,tn,qsg,qvg,qgold,snprim
5457 !-- SNOH is energy flux of snow phase change
5458 SNOH=RNET-QFX -HFX - SOH - X &
5459 +RHOnewCSN*NEWSNOW/DELT*(min(273.15,TABS)-soiltfrac) &
5460 +RAINF*CVW*PRCPMS*(max(273.15,TABS)-soiltfrac)
5462 !-- SMELT is speed of melting in M/S
5463 SMELT= SNOH /XLMELT*1.E-3
5464 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
5465 print *,'1- SMELT',i,j,smelt
5467 SMELT=AMIN1(SMELT,SNWEPR/DELT-BETA*EPOT*RAS)
5468 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
5469 print *,'2- SMELT',i,j,smelt
5471 SMELT=AMAX1(0.,SMELT)
5473 !18apr08 - Egglston limit
5474 ! SMELT= amin1 (smelt, 5.6E-7*meltfactor*max(1.,(soilt-273.15)))
5475 SMELT= amin1 (smelt, 5.6E-8*meltfactor*max(1.,(soilt-273.15)))
5476 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
5477 print *,'3- SMELT',i,j,smelt
5480 ! rr - potential melting
5481 rr=max(0.,SNWEPR/delt-BETA*EPOT*RAS)
5483 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
5484 print *,'4- SMELT i,j,smelt,rr',i,j,smelt,rr
5486 SNOHGNEW=SMELT*XLMELT*1.E3
5487 SNODIF=AMAX1(0.,(SNOH-SNOHGNEW))
5490 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
5491 print *,'SNOH,SNODIF',SNOH,SNODIF
5494 !*** From Koren et al. (1999) 13% of snow melt stays in the snow pack
5495 rsmfrac=min(0.18,(max(0.08,snwepr/0.10*0.13)))
5496 if(snhei > 0.01) then
5497 rsm=rsmfrac*smelt*delt
5499 ! do not keep melted water if snow depth is less that 1 cm
5502 !18apr08 rsm is part of melted water that stays in snow as liquid
5503 SMELT=max(0.,SMELT-rsm/delt)
5504 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
5505 print *,'5- SMELT i,j,smelt,rsm,snwepr,rsmfrac', &
5506 i,j,smelt,rsm,snwepr,rsmfrac
5509 !-- update of liquid equivalent of snow depth
5510 !-- due to evaporation and snow melt
5511 SNWE = AMAX1(0.,(SNWEPR- &
5512 (SMELT+BETA*EPOT*RAS)*DELT &
5513 ! (SMELT+BETA*EPOT*RAS)*DELT*snowfrac &
5514 ! (SMELT+BETA*EPOT*RAS*UMVEG)*DELT &
5516 !--- If there is no snow melting then just evaporation
5517 !--- or condensation cxhanges SNWE
5519 if(snhei.ne.0.) then
5520 EPOT=-QKMS*(QVATM-QSG)
5521 SNWE = AMAX1(0.,(SNWEPR- &
5522 BETA*EPOT*RAS*DELT))
5523 ! BETA*EPOT*RAS*DELT*snowfrac))
5527 !18apr08 - if snow melt occurred then go into iteration for energy budget
5529 if(nmelt.eq.1) goto 212 ! second interation
5532 if(smelt.gt.0..and.rsm.gt.0.) then
5533 if(snwe.le.rsm) then
5535 print *,'SNWE<RSM snwe,rsm,smelt*delt,epot*ras*delt,beta', &
5536 snwe,rsm,smelt*delt,epot*ras*delt,beta
5539 !*** Update snow density on effect of snow melt, melted
5540 !*** from the top of the snow. 13% of melted water
5541 !*** remains in the pack and changes its density.
5542 !*** Eq. 9 (with my correction) in Koren et al. (1999)
5543 xsn=(rhosn*(snwe-rsm)+1.e3*rsm)/ &
5545 rhosn=MIN(MAX(58.8,XSN),500.)
5546 !13mar18 rhosn=MIN(MAX(76.9,XSN),500.)
5549 thdifsn = 0.265/RHOCSN
5553 !--- Compute flux in the top snow layer
5554 IF(SNHEI.GE.SNTH)then
5555 S=thdifsn*RHOCSN*(soilt-TSOB)/SNPRIM
5557 S=D9*(tso(1)-tso(2))
5558 ELSEIF(SNHEI.lt.SNTH.and.SNHEI.GT.0.) then
5559 S=(fsn*thdifsn*rhocsn+fso*thdif(1)*rhcs)* &
5562 S=D9*(tso(1)-tso(2))
5566 S=D9*(tso(1)-tso(2))
5569 SNHEI=SNWE *1.E3 / RHOSN
5570 !-- If ground surface temperature
5571 !-- is above freezing snow can melt from the bottom. The following
5572 !-- piece of code will check if bottom melting is possible.
5574 IF(TSO(1).GT.273.15 .and. snhei > 0.) THEN
5575 if (snhei.GT.deltsn+snth) then
5576 hsn = snhei - deltsn
5577 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
5578 print*,'2 layer snow - snhei,hsn',snhei,hsn
5581 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
5582 print*,'1 layer snow or blended - snhei',snhei
5587 soiltfrac=snowfrac*273.15+(1.-snowfrac)*TSO(1)
5589 SNOHG=(TSO(1)-soiltfrac)*(cap(1)*zshalf(2)+ &
5590 RHOCSN*0.5*hsn) / DELT
5591 SNOHG=AMAX1(0.,SNOHG)
5593 SMELTG=SNOHG/XLMELT*1.E-3
5594 ! Egglston - empirical limit on snow melt from the bottom of snow pack
5595 SMELTG=AMIN1(SMELTG, 5.8e-9)
5597 ! rr - potential melting
5599 SMELTG=AMIN1(SMELTG, rr)
5601 SNOHGNEW=SMELTG*XLMELT*1.e3
5602 SNODIF=AMAX1(0.,(SNOHG-SNOHGNEW))
5603 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
5604 ! if(i.eq.266.and.j.eq.447) then
5605 print *,'TSO(1),soiltfrac,smeltg,SNODIF',TSO(1),soiltfrac,smeltg,SNODIF
5608 ! snwe=max(0.,snwe-smeltg*delt*snowfrac)
5609 snwe=max(0.,snwe-smeltg*delt)
5610 SNHEI=SNWE *1.E3 / RHOSN
5612 if(snhei > 0.) TSO(1) = soiltfrac
5613 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
5614 ! if(i.eq.266.and.j.eq.447) then
5615 print *,'Melt from the bottom snwe,snhei',snwe,snhei
5617 print *,'Snow is all melted on the warm ground'
5621 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
5622 print *,'SNHEI,SNOH',i,j,SNHEI,SNOH
5626 snheiprint=snweprint*1.E3 / RHOSN
5628 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
5629 print *, 'snweprint : ',snweprint
5630 print *, 'D9SN,SOILT,TSOB : ', D9SN,SOILT,TSOB
5633 X= (R21+D9SN*R22SN)*(soilt-TN) + &
5634 XLVM*R210*(QSG-QGOLD)
5635 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
5636 print *,'SNOWTEMP storage ',i,j,x
5637 print *,'R21,D9sn,r22sn,soiltfrac,soilt,tn,qsg,qgold,snprim', &
5638 R21,D9sn,r22sn,soiltfrac,soilt,tn,qsg,qgold,snprim
5642 ! "heat" from snow and rain
5643 -RHOnewCSN*NEWSNOW/DELT*(min(273.15,TABS)-SOILT) &
5644 -RAINF*CVW*PRCPMS*(max(273.15,TABS)-SOILT)
5645 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
5647 print *,'SNHEI=',snhei
5648 print *,'SNFLX=',snflx
5651 IF(SNHEI.GT.0.) THEN
5653 tsnav=0.5/snhei*((soilt+soilt1)*deltsn &
5654 +(soilt1+tso(1))*(SNHEI-DELTSN)) &
5657 tsnav=0.5*(soilt+tso(1)) - 273.15
5660 tsnav= soilt - 273.15
5663 !------------------------------------------------------------------------
5664 END SUBROUTINE SNOWTEMP
5665 !------------------------------------------------------------------------
5668 SUBROUTINE SOILMOIST ( &
5670 DELT,NZS,NDDZS,DTDZS,DTDZS2,RIW, &
5671 ZSMAIN,ZSHALF,DIFFU,HYDRO, &
5672 QSG,QVG,QCG,QCATM,QVATM,PRCP, &
5674 DEW,SMELT,SOILICE,VEGFRAC,SNOWFRAC,soilres, &
5676 DQM,QMIN,REF,KSAT,RAS,INFMAX, &
5678 SOILMOIS,SOILIQW,MAVAIL,RUNOFF,RUNOFF2,INFILTRP)
5679 !*************************************************************************
5680 ! moisture balance equation and Richards eqn.
5683 ! DELT - time step (s)
5684 ! IME,JME,NZS - dimensions of soil domain
5685 ! ZSMAIN - main levels in soil (m)
5686 ! ZSHALF - middle of the soil layers (m)
5687 ! DTDZS - dt/(2.*dzshalf*dzmain)
5688 ! DTDZS2 - dt/(2.*dzshalf)
5689 ! DIFFU - diffusional conductivity (m^2/s)
5690 ! HYDRO - hydraulic conductivity (m/s)
5691 ! QSG,QVG,QCG - saturated mixing ratio, mixing ratio of
5692 ! water vapor and cloud at the ground
5693 ! surface, respectively (kg/kg)
5694 ! QCATM,QVATM - cloud and water vapor mixing ratio
5695 ! at the first atm. level (kg/kg)
5696 ! PRCP - precipitation rate in m/s
5697 ! QKMS - exchange coefficient for water vapor in the
5698 ! surface layer (m/s)
5699 ! TRANSP - transpiration from the soil layers (m/s)
5700 ! DRIP - liquid water dripping from the canopy to soil (m)
5701 ! DEW - dew in kg/m^2s
5702 ! SMELT - melting rate in m/s
5703 ! SOILICE - volumetric content of ice in soil (m^3/m^3)
5704 ! SOILIQW - volumetric content of liquid water in soil (m^3/m^3)
5705 ! VEGFRAC - greeness fraction (0-1)
5706 ! RAS - ration of air density to soil density
5707 ! INFMAX - maximum infiltration rate (kg/m^2/s)
5709 ! SOILMOIS - volumetric soil moisture, 6 levels (m^3/m^3)
5710 ! MAVAIL - fraction of maximum soil moisture in the top
5712 ! RUNOFF - surface runoff (m/s)
5713 ! RUNOFF2 - underground runoff (m)
5714 ! INFILTRP - point infiltration flux into soil (m/s)
5715 ! /(snow bottom runoff) (mm/s)
5717 ! COSMC, RHSMC - coefficients for implicit solution of
5719 !******************************************************************
5721 !------------------------------------------------------------------
5722 !--- input variables
5723 REAL, INTENT(IN ) :: DELT
5724 INTEGER, INTENT(IN ) :: NZS,NDDZS
5728 REAL, DIMENSION(1:NZS), INTENT(IN ) :: ZSMAIN, &
5736 REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS
5738 REAL, INTENT(IN ) :: QSG,QVG,QCG,QCATM,QVATM , &
5739 QKMS,VEGFRAC,DRIP,PRCP , &
5740 DEW,SMELT,SNOWFRAC , &
5741 DQM,QMIN,REF,KSAT,RAS,RIW,SOILRES
5745 REAL, DIMENSION( 1:nzs ) , &
5747 INTENT(INOUT) :: SOILMOIS,SOILIQW
5749 REAL, INTENT(INOUT) :: MAVAIL,RUNOFF,RUNOFF2,INFILTRP, &
5754 REAL, DIMENSION( 1:nzs ) :: COSMC,RHSMC
5756 REAL :: DZS,R1,R2,R3,R4,R5,R6,R7,R8,R9,R10
5757 REAL :: REFKDT,REFDK,DELT1,F1MAX,F2MAX
5758 REAL :: F1,F2,FD,KDT,VAL,DDT,PX,FK,FKMAX
5759 REAL :: QQ,UMVEG,INFMAX1,TRANS
5760 REAL :: TOTLIQ,FLX,FLXSAT,QTOT
5761 REAL :: DID,X1,X2,X4,DENOM,Q2,Q4
5762 REAL :: dice,fcr,acrt,frzx,sum,cvfrz
5764 INTEGER :: NZS1,NZS2,K,KK,K1,KN,ialp1,jj,jk
5766 !******************************************************************************
5767 ! COEFFICIENTS FOR THOMAS ALGORITHM FOR SOILMOIS
5768 !******************************************************************************
5772 118 format(6(10Pf23.19))
5779 DID=(ZSMAIN(NZS)-ZSHALF(NZS))
5780 X1=ZSMAIN(NZS)-ZSMAIN(NZS1)
5782 !7may09 DID=(ZSMAIN(NZS)-ZSHALF(NZS))*2.
5783 ! DENOM=DID/DELT+DIFFU(NZS1)/X1
5784 ! COSMC(1)=DIFFU(NZS1)/X1/DENOM
5785 ! RHSMC(1)=(SOILMOIS(NZS)*DID/DELT
5786 ! 1 +TRANSP(NZS)-(HYDRO(NZS)*SOILMOIS(NZS)
5787 ! 1 -HYDRO(NZS1)*SOILMOIS(NZS1))*DID
5790 DENOM=(1.+DIFFU(nzs1)/X1/DID*DELT+HYDRO(NZS)/(2.*DID)*DELT)
5791 COSMC(1)=DELT*(DIFFU(nzs1)/DID/X1 &
5792 +HYDRO(NZS1)/2./DID)/DENOM
5793 RHSMC(1)=(SOILMOIS(NZS)+TRANSP(NZS)*DELT/ &
5796 ! RHSMC(1)=(SOILMOIS(NZS)*DID/DELT &
5797 ! +TRANSP(NZS)-(HYDRO(NZS)*SOILMOIS(NZS) &
5798 ! -HYDRO(NZS1)*SOILMOIS(NZS1))*DID &
5801 !12 June 2014 - low boundary condition: 1 - zero diffusion below the lowest
5802 ! level; 2 - soil moisture at the low boundary can be lost due to the root uptake.
5803 ! So far - no interaction with the water table.
5805 DENOM=1.+DIFFU(nzs1)/X1/DID*DELT
5806 !orig DENOM=(1.+DIFFU(nzs1)/X1/DID*DELT+HYDRO(NZS)/DID*DELT)
5807 !orig COSMC(1)=DELT*(DIFFU(nzs1)/DID/X1 &
5808 !orig +HYDRO(NZS1)/2./DID)/DENOM
5809 COSMC(1)=DELT*(DIFFU(nzs1)/DID/X1 &
5810 +HYDRO(NZS1)/DID)/DENOM
5812 ! RHSMC(1)=(SOILMOIS(NZS)+TRANSP(NZS)*DELT/ &
5815 RHSMC(1)=(SOILMOIS(NZS)-HYDRO(NZS)*DELT/DID*soilmois(nzs) &
5816 +TRANSP(NZS)*DELT/DID)/DENOM
5817 !test RHSMC(1)=SOILMOIS(NZS)-HYDRO(NZS)*soilmois(nzs)
5820 !this test gave smoother soil moisture, ovwerall better results
5822 RHSMC(1)=SOILMOIS(NZS)
5827 X4=2.*DTDZS(K1)*DIFFU(KN-1)
5828 X2=2.*DTDZS(K1+1)*DIFFU(KN)
5829 Q4=X4+HYDRO(KN-1)*DTDZS2(KN-1)
5830 Q2=X2-HYDRO(KN+1)*DTDZS2(KN-1)
5831 DENOM=1.+X2+X4-Q2*COSMC(K)
5833 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
5834 print *,'q2,soilmois(kn),DIFFU(KN),x2,HYDRO(KN+1),DTDZS2(KN-1),kn,k' &
5835 ,q2,soilmois(kn),DIFFU(KN),x2,HYDRO(KN+1),DTDZS2(KN-1),kn,k
5837 330 RHSMC(K+1)=(SOILMOIS(KN)+Q2*RHSMC(K) &
5839 /(ZSHALF(KN+1)-ZSHALF(KN)) &
5842 ! --- MOISTURE BALANCE BEGINS HERE
5845 UMVEG=(1.-VEGFRAC)*soilres
5856 !-- Total liquid water available on the top of soil domain
5857 !-- Without snow - 3 sources of water: precipitation,
5858 !-- water dripping from the canopy and dew
5859 !-- With snow - only one source of water - snow melt
5863 ! TOTLIQ=UMVEG*PRCP-DRIP/DELT-UMVEG*DEW*RAS-SMELT
5865 TOTLIQ=PRCP-DRIP/DELT-UMVEG*DEW*RAS-SMELT
5866 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
5867 print *,'UMVEG*PRCP,DRIP/DELT,UMVEG*DEW*RAS,SMELT', &
5868 UMVEG*PRCP,DRIP/DELT,UMVEG*DEW*RAS,SMELT
5871 !test 16 may TOTLIQ=UMVEG*PRCP-DRIP/DELT-UMVEG*DEW*RAS-SMELT
5872 !30july13 TOTLIQ=UMVEG*PRCP-DRIP/DELT-SMELT
5877 ! ----------- FROZEN GROUND VERSION -------------------------
5878 ! REFERENCE FROZEN GROUND PARAMETER, CVFRZ, IS A SHAPE PARAMETER OF
5879 ! AREAL DISTRIBUTION FUNCTION OF SOIL ICE CONTENT WHICH EQUALS 1/CV.
5880 ! CV IS A COEFFICIENT OF SPATIAL VARIATION OF SOIL ICE CONTENT.
5881 ! BASED ON FIELD DATA CV DEPENDS ON AREAL MEAN OF FROZEN DEPTH, AND IT
5882 ! CLOSE TO CONSTANT = 0.6 IF AREAL MEAN FROZEN DEPTH IS ABOVE 20 CM.
5883 ! THAT IS WHY PARAMETER CVFRZ = 3 (INT{1/0.6*0.6})
5885 ! Current logic doesn't allow CVFRZ be bigger than 3
5888 !-- SCHAAKE/KOREN EXPRESSION for calculation of max infiltration
5893 F2MAX=DQM*(ZSHALF(3)-ZSHALF(2))
5894 F1=F1MAX*(1.-SOILMOIS(1)/DQM)
5895 DICE=SOILICE(1)*ZSHALF(2)
5898 DICE=DICE+(ZSHALF(k+1)-ZSHALF(k))*SOILICE(K)
5899 FKMAX=DQM*(ZSHALF(k+1)-ZSHALF(k))
5900 FK=FKMAX*(1.-SOILMOIS(k)/DQM)
5903 KDT=REFKDT*KSAT/REFDK
5904 VAL=(1.-EXP(-KDT*DELT1))
5907 IF(PX.LT.0.0) PX = 0.0
5909 INFMAX1 = (PX*(DDT/(PX+DDT)))/DELT
5913 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
5914 print *,'INFMAX1 before frozen part',INFMAX1
5917 ! ----------- FROZEN GROUND VERSION --------------------------
5918 ! REDUCTION OF INFILTRATION BASED ON FROZEN GROUND PARAMETERS
5920 ! ------------------------------------------------------------------
5922 FRZX= 0.15*((dqm+qmin)/ref) * (0.412 / 0.468)
5924 IF ( DICE .GT. 1.E-2) THEN
5925 ACRT = CVFRZ * FRZX / DICE
5933 SUM = SUM + (ACRT ** ( CVFRZ-JK)) / FLOAT (K)
5935 FCR = 1. - EXP(-ACRT) * SUM
5937 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
5938 print *,'FCR--------',fcr
5939 print *,'DICE=',dice
5941 INFMAX1 = INFMAX1* FCR
5942 ! -------------------------------------------------------------------
5944 INFMAX = MAX(INFMAX1,HYDRO(1)*SOILMOIS(1))
5945 INFMAX = MIN(INFMAX, -TOTLIQ)
5946 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
5947 print *,'INFMAX,INFMAX1,HYDRO(1)*SOILIQW(1),-TOTLIQ', &
5948 INFMAX,INFMAX1,HYDRO(1)*SOILIQW(1),-TOTLIQ
5951 IF (-TOTLIQ.GT.INFMAX)THEN
5952 RUNOFF=-TOTLIQ-INFMAX
5954 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
5955 print *,'FLX,RUNOFF1=',flx,runoff
5958 ! INFILTRP is total infiltration flux in M/S
5960 ! Solution of moisture budget
5963 FLX=FLX-SOILMOIS(1)*R7
5964 ! R8 is for direct evaporation from soil, which occurs
5965 ! only from snow-free areas
5967 R8=UMVEG*R6*(1.-snowfrac)
5972 !-- evaporation regime
5974 QQ=(R5*R2-FLX+R9)/(R4-R5*R1-R10*R8/(REF-QMIN))
5975 FLXSAT=-DQM*(R4-R5*R1-R10*R8/(REF-QMIN)) &
5978 !-- dew formation regime
5979 QQ=(R2*R5-FLX+R8*(QTOT-QCG-QVG)+R9)/(R4-R1*R5)
5980 FLXSAT=-DQM*(R4-R1*R5)+R2*R5+R8*(QTOT-QVG-QCG)+R9
5984 ! print *,'negative QQ=',qq
5987 ELSE IF(QQ.GT.DQM) THEN
5990 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
5991 print *,'FLXSAT,FLX,DELT',FLXSAT,FLX,DELT,RUNOFF2
5993 ! RUNOFF2=(FLXSAT-FLX)
5994 RUNOFF=RUNOFF+(FLXSAT-FLX)
5996 SOILMOIS(1)=min(dqm,max(1.e-8,QQ))
5999 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
6000 print *,'SOILMOIS,SOILIQW, soilice',SOILMOIS,SOILIQW,soilice*riw
6001 print *,'COSMC,RHSMC',COSMC,RHSMC
6003 !--- FINAL SOLUTION FOR SOILMOIS
6007 QQ=COSMC(KK)*SOILMOIS(K-1)+RHSMC(KK)
6008 ! QQ=COSMC(KK)*SOILIQW(K-1)+RHSMC(KK)
6011 ! print *,'negative QQ=',qq
6014 ELSE IF(QQ.GT.DQM) THEN
6018 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
6019 print *,'hydro(k),QQ,DQM,k',hydro(k),QQ,DQM,k
6021 RUNOFF2=RUNOFF2+((QQ-DQM)*(ZSMAIN(K)-ZSHALF(K)))/DELT
6022 ! RUNOFF2=RUNOFF2+(QQ-DQM)*hydro(k)
6023 ! print *,'RUNOFF2=',RUNOFF2
6025 ! print *,'QQ,DQM,k',QQ,DQM,k
6026 RUNOFF2=RUNOFF2+((QQ-DQM)*(ZSHALF(K+1)-ZSHALF(K)))/DELT
6027 ! RUNOFF2=RUNOFF2+(QQ-DQM)*hydro(k)
6030 SOILMOIS(K)=min(dqm,max(1.e-8,QQ))
6033 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
6034 print *,'END soilmois,soiliqw,soilice',soilmois,SOILIQW,soilice*riw
6037 ! RUNOFF2=RUNOFF2+hydro(nzs)*SOILMOIS(NZS)
6038 ! MAVAIL=max(.00001,min(1.,SOILMOIS(1)/DQM))
6039 ! MAVAIL=max(.00001,min(1.,SOILMOIS(1)/(REF-QMIN)))
6040 MAVAIL=max(.00001,min(1.,(SOILMOIS(1)/(REF-QMIN)*(1.-snowfrac)+1.*snowfrac)))
6044 !-------------------------------------------------------------------
6045 END SUBROUTINE SOILMOIST
6046 !-------------------------------------------------------------------
6049 SUBROUTINE SOILPROP(spp_lsm,rstochcol,fieldcol_sf, &
6050 !--- input variables
6051 nzs,fwsat,lwsat,tav,keepfr, &
6052 soilmois,soiliqw,soilice, &
6053 soilmoism,soiliqwm,soilicem, &
6054 !--- soil fixed fields
6055 QWRTZ,rhocs,dqm,qmin,psis,bclh,ksat, &
6057 riw,xlmelt,CP,G0_P,cvw,ci, &
6059 !--- output variables
6060 thdif,diffu,hydro,cap)
6062 !******************************************************************
6063 ! SOILPROP computes thermal diffusivity, and diffusional and
6064 ! hydraulic condeuctivities
6065 !******************************************************************
6066 ! NX,NY,NZS - dimensions of soil domain
6067 ! FWSAT, LWSAT - volumetric content of frozen and liquid water
6068 ! for saturated condition at given temperatures (m^3/m^3)
6069 ! TAV - temperature averaged for soil layers (K)
6070 ! SOILMOIS - volumetric soil moisture at the main soil levels (m^3/m^3)
6071 ! SOILMOISM - volumetric soil moisture averaged for layers (m^3/m^3)
6072 ! SOILIQWM - volumetric liquid soil moisture averaged for layers (m^3/m^3)
6073 ! SOILICEM - volumetric content of soil ice averaged for layers (m^3/m^3)
6074 ! THDIF - thermal diffusivity for soil layers (W/m/K)
6075 ! DIFFU - diffusional conductivity (m^2/s)
6076 ! HYDRO - hydraulic conductivity (m/s)
6077 ! CAP - volumetric heat capacity (J/m^3/K)
6079 !******************************************************************
6082 !-----------------------------------------------------------------
6084 !--- soil properties
6085 INTEGER, INTENT(IN ) :: NZS
6087 INTENT(IN ) :: RHOCS, &
6095 REAL, DIMENSION( 1:nzs ) , &
6096 INTENT(IN ) :: SOILMOIS, &
6100 REAL, INTENT(IN ) :: CP, &
6109 REAL, DIMENSION(1:NZS), INTENT(IN) :: rstochcol
6110 REAL, DIMENSION(1:NZS), INTENT(INOUT) :: fieldcol_sf
6111 INTEGER, INTENT(IN ) :: spp_lsm
6114 !--- output variables
6115 REAL, DIMENSION(1:NZS) , &
6116 INTENT(INOUT) :: cap,diffu,hydro , &
6120 soilicem,soiliqwm , &
6123 !--- local variables
6124 REAL, DIMENSION(1:NZS) :: hk,detal,kasat,kjpl
6126 REAL :: x,x1,x2,x4,ws,wd,fact,fach,facd,psif,ci
6127 REAL :: tln,tavln,tn,pf,a,am,ame,h
6130 !-- for Johansen thermal conductivity
6131 REAL :: kzero,gamd,kdry,kas,x5,sr,ke
6136 !-- Constants for Johansen (1975) thermal conductivity
6137 kzero =2. ! if qwrtz > 0.2
6148 x1=xlmelt/(g0_p*psis)
6151 !--- Next 3 lines are for Johansen thermal conduct.
6153 kdry=(0.135*gamd+64.7)/(2700.-0.947*gamd)
6154 kas=kqwrtz**qwrtz*kzero**(1.-qwrtz)
6158 wd=ws - riw*soilicem(k)
6159 psif=psis*100.*(wd/(soiliqwm(k)+qmin))**bclh &
6161 !--- PSIF should be in [CM] to compute PF
6163 fact=1.+riw*soilicem(k)
6164 !--- HK is for McCumber thermal conductivity
6166 HK(K)=420.*EXP(-(PF+2.7))*fact
6171 IF(soilicem(k).NE.0.AND.TN.LT.0.) then
6172 !--- DETAL is taking care of energy spent on freezing or released from
6173 ! melting of soil water
6175 DETAL(K)=273.15*X2/(TAV(K)*TAV(K))* &
6176 (TAV(K)/(X1*TN))**X4
6178 if(keepfr(k).eq.1.) then
6184 !--- Next 10 lines calculate Johansen thermal conductivity KJPL
6185 kasat(k)=kas**(1.-ws)*kice**fwsat(k) &
6188 X5=(soilmoism(k)+qmin)/ws
6189 if(soilicem(k).eq.0.) then
6192 !--- next 2 lines - for coarse soils
6194 ! ke=0.7*log10(sr)+1.
6199 kjpl(k)=ke*(kasat(k)-kdry)+kdry
6201 !--- CAP -volumetric heat capacity
6202 CAP(K)=(1.-WS)*RHOCS &
6203 + (soiliqwm(K)+qmin)*CVW &
6205 + (dqm-soilmoism(k))*CP*1.2 &
6206 - DETAL(K)*1.e3*xlmelt
6210 if((ws-a).lt.0.12)then
6213 H=max(0.,(soilmoism(K)-a)/(max(1.e-8,(dqm-a))))
6215 if(a.ne.0.)facd=1.-a/max(1.e-8,soilmoism(K))
6216 ame=max(1.e-8,dqm-riw*soilicem(K))
6217 !--- DIFFU is diffusional conductivity of soil water
6218 diffu(K)=-BCLH*KSAT*PSIS/ame* &
6223 ! diffu(K)=-BCLH*KSAT*PSIS/dqm &
6227 !--- thdif - thermal diffusivity
6228 ! thdif(K)=HK(K)/CAP(K)
6229 !--- Use thermal conductivity from Johansen (1975)
6230 thdif(K)=KJPL(K)/CAP(K)
6234 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
6235 print *,'soilice*riw,soiliqw,soilmois,ws',soilice*riw,soiliqw,soilmois,ws
6239 if((ws-riw*soilice(k)).lt.0.12)then
6243 if(soilice(k).ne.0.) &
6244 fach=1.-riw*soilice(k)/max(1.e-8,soilmois(k))
6245 am=max(1.e-8,dqm-riw*soilice(k))
6246 !--- HYDRO is hydraulic conductivity of soil water
6247 hydro(K)=min(KSAT,KSAT/am* &
6251 if(hydro(k)<1.e-10)hydro(k)=0.
6256 !-----------------------------------------------------------------------
6257 END SUBROUTINE SOILPROP
6258 !-----------------------------------------------------------------------
6261 SUBROUTINE TRANSF(i,j, &
6262 !--- input variables
6263 nzs,nroot,soiliqw,tabs,lai,gswin, &
6264 !--- soil fixed fields
6265 dqm,qmin,ref,wilt,zshalf,pc,iland, &
6266 !--- output variables
6269 !-------------------------------------------------------------------
6270 !--- TRANF(K) - THE TRANSPIRATION FUNCTION (Smirnova et al. 1996, EQ. 18,19)
6271 !*******************************************************************
6272 ! NX,NY,NZS - dimensions of soil domain
6273 ! SOILIQW - volumetric liquid soil moisture at the main levels (m^3/m^3)
6274 ! TRANF - the transpiration function at levels (m)
6275 ! TRANSUM - transpiration function integrated over the rooting zone (m)
6277 !*******************************************************************
6279 !-------------------------------------------------------------------
6281 !--- input variables
6283 INTEGER, INTENT(IN ) :: i,j,nroot,nzs, iland
6286 INTENT(IN ) :: GSWin, TABS, lai
6287 !--- soil properties
6289 INTENT(IN ) :: DQM, &
6295 REAL, DIMENSION(1:NZS), INTENT(IN) :: soiliqw, &
6299 REAL, DIMENSION(1:NZS), INTENT(OUT) :: TRANF
6300 REAL, INTENT(OUT) :: TRANSUM
6306 !-- for non-linear root distribution
6307 REAL :: gx,sm1,sm2,sm3,sm4,ap0,ap1,ap2,ap3,ap4
6308 REAL :: FTEM, PCtot, fsol, f1, cmin, cmax, totcnd
6309 REAL, DIMENSION(1:NZS) :: PART
6310 !--------------------------------------------------------------------
6318 totliq=soiliqw(1)+qmin
6328 gx=ap0+ap1*sm1+ap2*sm2+ap3*sm3+ap4*sm4
6329 if(totliq.ge.ref) gx=1.
6330 if(totliq.le.0.) gx=0.
6335 IF(TOTLIQ.GT.REF) THEN
6337 ELSE IF(TOTLIQ.LE.WILT) THEN
6340 TRANF(1)=(TOTLIQ-WILT)/(REF-WILT)*DID
6342 !-- uncomment next line for non-linear root distribution
6346 totliq=soiliqw(k)+qmin
6351 gx=ap0+ap1*sm1+ap2*sm2+ap3*sm3+ap4*sm4
6352 if(totliq.ge.ref) gx=1.
6353 if(totliq.le.0.) gx=0.
6356 DID=zshalf(K+1)-zshalf(K)
6358 IF(totliq.GE.REF) THEN
6360 ELSE IF(totliq.LE.WILT) THEN
6363 TRANF(K)=(totliq-WILT) &
6366 !-- uncomment next line for non-linear root distribution
6370 ! For LAI> 3 => transpiration at potential rate (F.Tardieu, 2013)
6375 !- 26aug16- next 2 lines could lead to LH increase and higher 2-m Q during the day
6376 ! pctot=min(0.8,pc*lai)
6377 ! pctot=min(0.8,max(pc,pc*lai))
6379 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
6380 ! if (i==421.and.j==280) then
6381 print *,'i,j,pctot,lai,pc',i,j,pctot,lai,pc
6384 !--- air temperature function
6385 ! Avissar (1985) and AX 7/95
6386 IF (TABS .LE. 302.15) THEN
6387 FTEM = 1.0 / (1.0 + EXP(-0.41 * (TABS - 282.05)))
6389 FTEM = 1.0 / (1.0 + EXP(0.5 * (TABS - 314.0)))
6391 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
6392 ! if (i==421.and.j==280) then
6393 print *,'i,j,tabs,ftem',i,j,tabs,ftem
6395 !--- incoming solar function
6396 cmin = 1./rsmax_data
6397 cmax = 1./rstbl(iland)
6399 cmax = lai/rstbl(iland) ! max conductance
6401 ! Noihlan & Planton (1988)
6403 ! if(lai > 0.01) then
6404 ! f1 = 1.1/lai*gswin/rgltbl(iland)! f1=0. when GSWin=0.
6405 ! fsol = (f1+cmin/cmax)/(1.+f1)
6410 ! totcnd = max(lai/rstbl(iland), pctot * ftem * f1)
6411 ! Mahrer & Avissar (1982), Avissar et al. (1985)
6412 if (GSWin < rgltbl(iland)) then
6413 fsol = 1. / (1. + exp(-0.034 * (GSWin - 3.5)))
6417 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
6418 ! if (i==421.and.j==280) then
6419 print *,'i,j,GSWin,lai,f1,fsol',i,j,gswin,lai,f1,fsol
6421 !--- total conductance
6422 totcnd =(cmin + (cmax - cmin)*pctot*ftem*fsol)/cmax
6424 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
6425 ! if (i==421.and.j==280) then
6426 print *,'i,j,iland,RGLTBL(iland),RSTBL(iland),RSMAX_DATA,totcnd' &
6427 ,i,j,iland,RGLTBL(iland),RSTBL(iland),RSMAX_DATA,totcnd
6430 !-- TRANSUM - total for the rooting zone
6433 ! linear root distribution
6434 TRANF(k)=max(cmin,TRANF(k)*totcnd)
6435 transum=transum+tranf(k)
6437 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
6438 ! if (i==421.and.j==280) then
6439 print *,'i,j,transum,TRANF',i,j,transum,tranf
6442 !-----------------------------------------------------------------
6443 END SUBROUTINE TRANSF
6444 !-----------------------------------------------------------------
6447 SUBROUTINE VILKA(TN,D1,D2,PP,QS,TS,TT,NSTEP,ii,j,iland,isoil)
6448 !--------------------------------------------------------------
6449 !--- VILKA finds the solution of energy budget at the surface
6450 !--- using table T,QS computed from Clausius-Klapeiron
6451 !--------------------------------------------------------------
6452 REAL, DIMENSION(1:5001), INTENT(IN ) :: TT
6453 REAL, INTENT(IN ) :: TN,D1,D2,PP
6454 INTEGER, INTENT(IN ) :: NSTEP,ii,j,iland,isoil
6456 REAL, INTENT(OUT ) :: QS, TS
6461 I=(TN-1.7315E2)/.05+1
6462 T1=173.1+FLOAT(I)*.05
6464 I1=I-F1/(.05+D1*(TT(I+1)-TT(I)))
6466 IF(I.GT.5000.OR.I.LT.1) GOTO 1
6468 T1=173.1+FLOAT(I)*.05
6470 RN=F1/(.05+D1*(TT(I+1)-TT(I)))
6472 IF(I.GT.5000.OR.I.LT.1) GOTO 1
6475 QS=(TT(I)+(TT(I)-TT(I+1))*RN)/PP
6477 ! 1 PRINT *,'Crash in surface energy budget - STOP'
6478 1 PRINT *,' AVOST IN VILKA Table index= ',I
6479 ! PRINT *,TN,D1,D2,PP,NSTEP,I,TT(i),ii,j,iland,isoil
6480 print *,'I,J=',ii,j,'LU_index = ',iland, 'Psfc[hPa] = ',pp, 'Tsfc = ',tn
6481 CALL wrf_error_fatal (' Crash in surface energy budget ' )
6483 !-----------------------------------------------------------------------
6484 END SUBROUTINE VILKA
6485 !-----------------------------------------------------------------------
6487 SUBROUTINE SOILVEGIN ( mosaic_lu,mosaic_soil,soilfrac,nscat, &
6489 NLCAT,IVGTYP,ISLTYP,iswater, &
6490 IFOREST,lufrac,vegfrac,EMISS,PC,ZNT,LAI,RDLAI2D,&
6491 QWRTZ,RHOCS,BCLH,DQM,KSAT,PSIS,QMIN,REF,WILT,I,J)
6493 !************************************************************************
6494 ! Set-up soil and vegetation Parameters in the case when
6495 ! snow disappears during the forecast and snow parameters
6496 ! shold be replaced by surface parameters according to
6497 ! soil and vegetation types in this point.
6503 ! DQM: MAX soil moisture content - MIN (m^3/m^3)
6504 ! REF: Reference soil moisture (m^3/m^3)
6505 ! WILT: Wilting PT soil moisture contents (m^3/m^3)
6506 ! QMIN: Air dry soil moist content limits (m^3/m^3)
6507 ! PSIS: SAT soil potential coefs. (m)
6508 ! KSAT: SAT soil diffusivity/conductivity coefs. (m/s)
6509 ! BCLH: Soil diffusivity/conductivity exponent.
6511 ! ************************************************************************
6513 !---------------------------------------------------------------------------
6514 integer, parameter :: nsoilclas=19
6515 integer, parameter :: nvegclas=24+3
6516 integer, parameter :: ilsnow=99
6518 INTEGER, INTENT(IN ) :: nlcat, nscat, iswater, i, j
6520 !--- soiltyp classification according to STATSGO(nclasses=16)
6523 ! 2 LOAMY SAND LOAMY SAND
6524 ! 3 SANDY LOAM SANDY LOAM
6525 ! 4 SILT LOAM SILTY LOAM
6528 ! 7 SANDY CLAY LOAM SANDY CLAY LOAM
6529 ! 8 SILTY CLAY LOAM SILTY CLAY LOAM
6530 ! 9 CLAY LOAM CLAY LOAM
6531 ! 10 SANDY CLAY SANDY CLAY
6532 ! 11 SILTY CLAY SILTY CLAY
6533 ! 12 CLAY LIGHT CLAY
6534 ! 13 ORGANIC MATERIALS LOAM
6537 ! Bedrock is reclassified as class 14
6538 ! 16 OTHER (land-ice)
6543 !----------------------------------------------------------------------
6544 REAL LQMA(nsoilclas),LRHC(nsoilclas), &
6545 LPSI(nsoilclas),LQMI(nsoilclas), &
6546 LBCL(nsoilclas),LKAS(nsoilclas), &
6547 LWIL(nsoilclas),LREF(nsoilclas), &
6549 !-- LQMA Rawls et al.[1982]
6550 ! DATA LQMA /0.417, 0.437, 0.453, 0.501, 0.486, 0.463, 0.398,
6551 ! & 0.471, 0.464, 0.430, 0.479, 0.475, 0.439, 1.0, 0.20, 0.401/
6553 !-- Clapp, R. and G. Hornberger, 1978: Empirical equations for some soil
6554 ! hydraulic properties, Water Resour. Res., 14, 601-604.
6556 !-- Clapp et al. [1978]
6557 DATA LQMA /0.395, 0.410, 0.435, 0.485, 0.485, 0.451, 0.420, &
6558 0.477, 0.476, 0.426, 0.492, 0.482, 0.451, 1.0, &
6559 0.20, 0.435, 0.468, 0.200, 0.339/
6561 !-- LREF Rawls et al.[1982]
6562 ! DATA LREF /0.091, 0.125, 0.207, 0.330, 0.360, 0.270, 0.255,
6563 ! & 0.366, 0.318, 0.339, 0.387, 0.396, 0.329, 1.0, 0.108, 0.283/
6565 !-- Clapp et al. [1978]
6566 DATA LREF /0.174, 0.179, 0.249, 0.369, 0.369, 0.314, 0.299, &
6567 0.357, 0.391, 0.316, 0.409, 0.400, 0.314, 1., &
6568 0.1, 0.249, 0.454, 0.17, 0.236/
6570 !-- LWIL Rawls et al.[1982]
6571 ! DATA LWIL/0.033, 0.055, 0.095, 0.133, 0.133, 0.117, 0.148,
6572 ! & 0.208, 0.197, 0.239, 0.250, 0.272, 0.066, 0.0, 0.006, 0.029/
6574 !-- Clapp et al. [1978]
6575 DATA LWIL/0.068, 0.075, 0.114, 0.179, 0.179, 0.155, 0.175, &
6576 0.218, 0.250, 0.219, 0.283, 0.286, 0.155, 0.0, &
6577 0.006, 0.114, 0.030, 0.006, 0.01/
6579 ! DATA LQMI/0.010, 0.028, 0.047, 0.084, 0.084, 0.066, 0.067,
6580 ! & 0.120, 0.103, 0.100, 0.126, 0.138, 0.066, 0.0, 0.006, 0.028/
6582 !-- Carsel and Parrish [1988]
6583 DATA LQMI/0.045, 0.057, 0.065, 0.067, 0.034, 0.078, 0.10, &
6584 0.089, 0.095, 0.10, 0.070, 0.068, 0.078, 0.0, &
6585 0.004, 0.065, 0.020, 0.004, 0.008/
6587 !-- LPSI Cosby et al[1984]
6588 ! DATA LPSI/0.060, 0.036, 0.141, 0.759, 0.759, 0.355, 0.135,
6589 ! & 0.617, 0.263, 0.098, 0.324, 0.468, 0.355, 0.0, 0.069, 0.036/
6590 ! & 0.617, 0.263, 0.098, 0.324, 0.468, 0.355, 0.0, 0.069, 0.036/
6592 !-- Clapp et al. [1978]
6593 DATA LPSI/0.121, 0.090, 0.218, 0.786, 0.786, 0.478, 0.299, &
6594 0.356, 0.630, 0.153, 0.490, 0.405, 0.478, 0.0, &
6595 0.121, 0.218, 0.468, 0.069, 0.069/
6597 !-- LKAS Rawls et al.[1982]
6598 ! DATA LKAS/5.83E-5, 1.70E-5, 7.19E-6, 1.89E-6, 1.89E-6,
6599 ! & 3.67E-6, 1.19E-6, 4.17E-7, 6.39E-7, 3.33E-7, 2.50E-7,
6600 ! & 1.67E-7, 3.38E-6, 0.0, 1.41E-4, 1.41E-5/
6602 !-- Clapp et al. [1978]
6603 DATA LKAS/1.76E-4, 1.56E-4, 3.47E-5, 7.20E-6, 7.20E-6, &
6604 6.95E-6, 6.30E-6, 1.70E-6, 2.45E-6, 2.17E-6, &
6605 1.03E-6, 1.28E-6, 6.95E-6, 0.0, 1.41E-4, &
6606 3.47E-5, 1.28E-6, 1.41E-4, 1.76E-4/
6608 !-- LBCL Cosby et al [1984]
6609 ! DATA LBCL/2.79, 4.26, 4.74, 5.33, 5.33, 5.25, 6.66,
6610 ! & 8.72, 8.17, 10.73, 10.39, 11.55, 5.25, 0.0, 2.79, 4.26/
6612 !-- Clapp et al. [1978]
6613 DATA LBCL/4.05, 4.38, 4.90, 5.30, 5.30, 5.39, 7.12, &
6614 7.75, 8.52, 10.40, 10.40, 11.40, 5.39, 0.0, &
6615 4.05, 4.90, 11.55, 2.79, 2.79/
6617 DATA LRHC /1.47,1.41,1.34,1.27,1.27,1.21,1.18,1.32,1.23, &
6618 1.18,1.15,1.09,1.21,4.18,2.03,2.10,1.09,2.03,1.47/
6620 DATA DATQTZ/0.92,0.82,0.60,0.25,0.10,0.40,0.60,0.10,0.35, &
6621 0.52,0.10,0.25,0.00,0.,0.60,0.0,0.25,0.60,0.92/
6623 !--------------------------------------------------------------------------
6625 ! USGS Vegetation Types
6627 ! 1: Urban and Built-Up Land
6628 ! 2: Dryland Cropland and Pasture
6629 ! 3: Irrigated Cropland and Pasture
6630 ! 4: Mixed Dryland/Irrigated Cropland and Pasture
6631 ! 5: Cropland/Grassland Mosaic
6632 ! 6: Cropland/Woodland Mosaic
6635 ! 9: Mixed Shrubland/Grassland
6637 ! 11: Deciduous Broadleaf Forest
6638 ! 12: Deciduous Needleleaf Forest
6639 ! 13: Evergreen Broadleaf Forest
6640 ! 14: Evergreen Needleleaf Fores
6643 ! 17: Herbaceous Wetland
6644 ! 18: Wooded Wetland
6645 ! 19: Barren or Sparsely Vegetated
6646 ! 20: Herbaceous Tundra
6649 ! 23: Bare Ground Tundra
6656 ! MODIS vegetation categories from VEGPARM.TBL
6657 ! 1: Evergreen Needleleaf Forest
6658 ! 2: Evergreen Broadleaf Forest
6659 ! 3: Deciduous Needleleaf Forest
6660 ! 4: Deciduous Broadleaf Forest
6662 ! 6: Closed Shrublands
6663 ! 7: Open Shrublands
6667 ! 11: Permanent wetlands
6669 ! 13: Urban and Built-Up
6670 ! 14: cropland/natural vegetation mosaic
6672 ! 16: Barren or Sparsely Vegetated
6680 !---- Below are the arrays for the vegetation parameters
6681 REAL LALB(nvegclas),LMOI(nvegclas),LEMI(nvegclas), &
6682 LROU(nvegclas),LTHI(nvegclas),LSIG(nvegclas), &
6685 !************************************************************************
6686 !---- vegetation parameters
6690 DATA LALB/.18,.17,.18,.18,.18,.16,.19,.22,.20,.20,.16,.14, &
6691 .12,.12,.13,.08,.14,.14,.25,.15,.15,.15,.25,.55, &
6693 DATA LEMI/.88,4*.92,.93,.92,.88,.9,.92,.93,.94, &
6694 .95,.95,.94,.98,.95,.95,.85,.92,.93,.92,.85,.95, &
6696 !-- Roughness length is changed for forests and some others
6697 ! DATA LROU/.5,.06,.075,.065,.05,.2,.075,.1,.11,.15,.8,.85, &
6698 ! 2.0,1.0,.563,.0001,.2,.4,.05,.1,.15,.1,.065,.05/
6699 DATA LROU/.5,.06,.075,.065,.05,.2,.075,.1,.11,.15,.5,.5, &
6700 .5,.5,.5,.0001,.2,.4,.05,.1,.15,.1,.065,.05, &
6703 DATA LMOI/.1,.3,.5,.25,.25,.35,.15,.1,.15,.15,.3,.3, &
6704 .5,.3,.3,1.,.6,.35,.02,.5,.5,.5,.02,.95,.40,.50,.40/
6706 !---- still needs to be corrected
6708 ! DATA LPC/ 15*.8,0.,.8,.8,.5,.5,.5,.5,.5,.0/
6709 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, &
6710 0.3,0.3,0.4,0.4,0.3,0.,.3,0.,0./
6712 ! used in RUC DATA LPC /0.6,6*0.8,0.7,0.75,6*0.8,0.,0.8,0.8, &
6713 ! 0.5,0.7,0.6,0.7,0.5,0./
6716 !***************************************************************************
6722 INTEGER, INTENT(IN ) :: mosaic_lu, mosaic_soil
6724 REAL, INTENT(IN ) :: SHDMAX
6725 REAL, INTENT(IN ) :: SHDMIN
6726 REAL, INTENT(IN ) :: VEGFRAC
6727 REAL, DIMENSION( 1:NLCAT ), INTENT(IN):: LUFRAC
6728 REAL, DIMENSION( 1:NSCAT ), INTENT(IN):: SOILFRAC
6734 INTENT (INOUT ) :: emiss, &
6737 LOGICAL, intent(in) :: rdlai2d
6738 !--- soil properties
6740 INTENT( OUT) :: RHOCS, &
6749 INTEGER, INTENT ( OUT) :: iforest
6751 ! INTEGER, DIMENSION( 1:(lucats) ) , &
6752 ! INTENT ( OUT) :: iforest
6755 ! INTEGER, DIMENSION( 1:50 ) :: if1
6756 INTEGER :: kstart, kfin, lstart, lfin
6758 REAL :: area, factor, znt1, lb
6759 REAL, DIMENSION( 1:NLCAT ) :: ZNTtoday, LAItoday, deltalai
6761 !***********************************************************************
6762 ! DATA ZS1/0.0,0.05,0.20,0.40,1.6,3.0/ ! o - levels in soil
6763 ! DATA ZS2/0.0,0.025,0.125,0.30,1.,2.3/ ! x - levels in soil
6765 ! DATA IF1/12*0,1,1,1,12*0/
6771 iforest = IFORTBL(IVGTYP)
6773 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
6774 if(i.eq.375.and.j.eq.254)then
6775 print *,'ifortbl(ivgtyp),ivgtyp,laitbl(ivgtyp),z0tbl(ivgtyp)', &
6776 ifortbl(ivgtyp),ivgtyp,laitbl(ivgtyp),z0tbl(ivgtyp)
6782 ! 11oct2012 - seasonal correction on ZNT for crops and LAI for all veg. types
6783 ! factor = 1 with minimum greenness --> vegfrac = shdmin (cold season)
6784 ! factor = 0 with maximum greenness --> vegfrac = shdmax
6785 ! SHDMAX, SHDMIN and VEGFRAC are in % here.
6786 if((shdmax - shdmin) .lt. 1) then
6787 factor = 1. ! min greenness
6789 factor = 1. - max(0.,min(1.,(vegfrac - shdmin)/max(1.,(shdmax-shdmin))))
6792 ! 18sept18 - LAITBL and Z0TBL are the max values
6794 if(IFORTBL(k) == 1) deltalai(k)=min(0.2,0.8*LAITBL(K))
6795 if(IFORTBL(k) == 2 .or. IFORTBL(k) == 7) deltalai(k)=min(0.5,0.8*LAITBL(K))
6796 if(IFORTBL(k) == 3) deltalai(k)=min(0.45,0.8*LAITBL(K))
6797 if(IFORTBL(k) == 4) deltalai(k)=min(0.75,0.8*LAITBL(K))
6798 if(IFORTBL(k) == 5) deltalai(k)=min(0.86,0.8*LAITBL(K))
6800 if(k.ne.iswater) then
6801 !-- 20aug18 - change in LAItoday based on the greenness fraction for the current day
6802 LAItoday(k) = LAITBL(K) - deltalai(k) * factor
6804 if(IFORTBL(k) == 7) then
6805 !-- seasonal change of roughness length for crops
6806 ZNTtoday(k) = Z0TBL(K) - 0.125 * factor
6808 ZNTtoday(k) = Z0TBL(K)
6811 LAItoday(k) = LAITBL(K)
6812 ! ZNTtoday(k) = Z0TBL(K)
6813 ZNTtoday(k) = ZNT ! do not overwrite z0 over water with the table value
6817 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
6818 if(i.eq.358.and.j.eq.260)then
6819 print *,'ivgtyp,factor,vegfrac,shdmin,shdmax,deltalai,laitoday(ivgtyp),znttoday(ivgtyp)', &
6820 i,j,ivgtyp,factor,vegfrac,shdmin,shdmax,deltalai,laitoday(ivgtyp),znttoday(ivgtyp)
6828 if(.not.rdlai2d) LAI = 0.
6830 !-- mosaic approach to landuse in the grid box
6831 ! Use Mason (1988) Eq.(15) to compute effective ZNT;
6832 ! Lb - blending height = L/200., where L is the length scale
6833 ! of regions with varying Z0 (Lb = 5 if L=1000 m)
6835 if(mosaic_lu == 1) then
6837 AREA = AREA + lufrac(k)
6838 EMISS = EMISS+ LEMITBL(K)*lufrac(k)
6839 ZNT = ZNT + lufrac(k)/ALOG(LB/ZNTtoday(K))**2.
6840 ! ZNT1 - weighted average in the grid box, not used, computed for comparison
6841 ZNT1 = ZNT1 + lufrac(k)*ZNTtoday(K)
6842 if(.not.rdlai2d) LAI = LAI + LAItoday(K)*lufrac(k)
6843 PC = PC + PCTBL(K)*lufrac(k)
6846 if (area.gt.1.) area=1.
6847 if (area <= 0.) then
6848 print *,'Bad area of grid box', area
6852 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
6853 if(i.eq.358.and.j.eq.260) then
6854 print *,'area=',area,i,j,ivgtyp,nlcat,(lufrac(k),k=1,nlcat),EMISS,ZNT,ZNT1,LAI,PC
6860 ZNT = LB/EXP(SQRT(1./ZNT))
6861 if(.not.rdlai2d) LAI = LAI/AREA
6864 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
6865 if(i.eq.358.and.j.eq.260) then
6866 print *,'mosaic=',i,j,ivgtyp,nlcat,(lufrac(k),k=1,nlcat),EMISS,ZNT,ZNT1,LAI,PC
6872 EMISS = LEMITBL(IVGTYP)
6873 ZNT = ZNTtoday(IVGTYP)
6875 if(.not.rdlai2d) LAI = LAItoday(IVGTYP)
6878 ! parameters from SOILPARM.TBL
6890 if(mosaic_soil == 1 ) then
6893 !exclude watrer points from this loop
6894 AREA = AREA + soilfrac(k)
6895 RHOCS = RHOCS + HC(k)*1.E6*soilfrac(k)
6896 BCLH = BCLH + BB(K)*soilfrac(k)
6897 DQM = DQM + (MAXSMC(K)- &
6898 DRYSMC(K))*soilfrac(k)
6899 KSAT = KSAT + SATDK(K)*soilfrac(k)
6900 PSIS = PSIS - SATPSI(K)*soilfrac(k)
6901 QMIN = QMIN + DRYSMC(K)*soilfrac(k)
6902 REF = REF + REFSMC(K)*soilfrac(k)
6903 WILT = WILT + WLTSMC(K)*soilfrac(k)
6904 QWRTZ = QWRTZ + QTZ(K)*soilfrac(k)
6907 if (area.gt.1.) area=1.
6908 if (area <= 0.) then
6909 ! area = 0. for water points
6910 ! print *,'Area of a grid box', area, 'iswater = ',iswater
6911 RHOCS = HC(ISLTYP)*1.E6
6913 DQM = MAXSMC(ISLTYP)- &
6915 KSAT = SATDK(ISLTYP)
6916 PSIS = - SATPSI(ISLTYP)
6917 QMIN = DRYSMC(ISLTYP)
6918 REF = REFSMC(ISLTYP)
6919 WILT = WLTSMC(ISLTYP)
6933 ! dominant category approach
6935 if(isltyp.ne.14) then
6936 RHOCS = HC(ISLTYP)*1.E6
6938 DQM = MAXSMC(ISLTYP)- &
6940 KSAT = SATDK(ISLTYP)
6941 PSIS = - SATPSI(ISLTYP)
6942 QMIN = DRYSMC(ISLTYP)
6943 REF = REFSMC(ISLTYP)
6944 WILT = WLTSMC(ISLTYP)
6949 ! parameters from the look-up tables
6950 ! BCLH = LBCL(ISLTYP)
6951 ! DQM = LQMA(ISLTYP)- &
6953 ! KSAT = LKAS(ISLTYP)
6954 ! PSIS = - LPSI(ISLTYP)
6955 ! QMIN = LQMI(ISLTYP)
6956 ! REF = LREF(ISLTYP)
6957 ! WILT = LWIL(ISLTYP)
6958 ! QWRTZ = DATQTZ(ISLTYP)
6960 !--------------------------------------------------------------------------
6961 END SUBROUTINE SOILVEGIN
6962 !--------------------------------------------------------------------------
6964 SUBROUTINE RUCLSMINIT( SH2O,SMFR3D,TSLB,SMOIS,ISLTYP,IVGTYP, &
6965 mminlu, XICE,mavail,nzs, iswater, isice, &
6966 znt, restart, allowed_to_read , &
6967 ids,ide, jds,jde, kds,kde, &
6968 ims,ime, jms,jme, kms,kme, &
6969 its,ite, jts,jte, kts,kte )
6970 #if ( WRF_CHEM == 1 )
6971 USE module_data_gocart_dust
6976 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, &
6977 ims,ime, jms,jme, kms,kme, &
6978 its,ite, jts,jte, kts,kte, &
6980 CHARACTER(LEN=*), INTENT(IN ) :: MMINLU
6982 REAL, DIMENSION( ims:ime, 1:nzs, jms:jme ) , &
6983 INTENT(IN) :: TSLB, &
6986 INTEGER, DIMENSION( ims:ime, jms:jme ) , &
6987 INTENT(INOUT) :: ISLTYP,IVGTYP
6989 REAL, DIMENSION( ims:ime, 1:nzs, jms:jme ) , &
6990 INTENT(INOUT) :: SMFR3D, &
6993 REAL, DIMENSION( ims:ime, jms:jme ) , &
6994 INTENT(INOUT) :: XICE,MAVAIL
6996 REAL, DIMENSION( ims:ime, jms:jme ) , &
6999 REAL, DIMENSION ( 1:nzs ) :: SOILIQW
7001 LOGICAL , INTENT(IN) :: restart, allowed_to_read
7004 INTEGER :: I,J,L,itf,jtf
7005 REAL :: RIW,XLMELT,TLN,DQM,REF,PSIS,QMIN,BCLH
7007 character*8 :: MMINLURUC, MMINSL
7011 ! itf=min0(ite,ide-1)
7012 ! jtf=min0(jte,jde-1)
7018 ! initialize three LSM related tables
7019 IF ( allowed_to_read ) THEN
7020 CALL wrf_message( 'INITIALIZE THREE LSM RELATED TABLES' )
7021 if(mminlu == 'USGS') then
7022 MMINLURUC='USGS-RUC'
7023 elseif(mminlu == 'MODIS' .OR. &
7024 & mminlu == 'MODIFIED_IGBP_MODIS_NOAH') then
7025 MMINLURUC='MODI-RUC'
7028 ! CALL RUCLSM_PARM_INIT
7029 print *,'RUCLSMINIT uses ',mminluruc
7030 call RUCLSM_SOILVEGPARM( MMINLURUC, MMINSL)
7033 !#if ( WRF_CHEM == 1 )
7035 ! need this parameter for dust parameterization in wrf/chem
7038 ! porosity(i)=maxsmc(i)
7039 ! drypoint(i)=drysmc(i)
7043 IF(.not.restart)THEN
7051 IF ( ISLTYP( i,j ) .LT. 1 ) THEN
7053 WRITE(err_message,*)"module_sf_ruclsm.F: lsminit: out of range ISLTYP ",i,j,ISLTYP( i,j )
7054 CALL wrf_message(err_message)
7058 IF ( errflag .EQ. 1 ) THEN
7059 CALL wrf_error_fatal( "module_sf_ruclsm.F: lsminit: out of range value "// &
7060 "of ISLTYP. Is this field in the input?" )
7066 ZNT(I,J) = Z0TBL(IVGTYP(I,J))
7068 ! CALL SOILIN ( ISLTYP(I,J), DQM, REF, PSIS, QMIN, BCLH )
7071 !--- Computation of volumetric content of ice in soil
7072 !--- and initialize MAVAIL
7073 DQM = MAXSMC (ISLTYP(I,J)) - &
7074 DRYSMC (ISLTYP(I,J))
7075 REF = REFSMC (ISLTYP(I,J))
7076 PSIS = - SATPSI (ISLTYP(I,J))
7077 QMIN = DRYSMC (ISLTYP(I,J))
7078 BCLH = BB (ISLTYP(I,J))
7081 !!! IF (.not.restart) THEN
7083 IF(xice(i,j).gt.0.) THEN
7091 if(isltyp(i,j).ne.14 ) then
7093 mavail(i,j) = max(0.00001,min(1.,(smois(i,1,j)-qmin)/(ref-qmin)))
7095 !-- for land points initialize soil ice
7096 tln=log(TSLB(i,l,j)/273.15)
7099 soiliqw(l)=(dqm+qmin)*(XLMELT* &
7100 (tslb(i,l,j)-273.15)/tslb(i,l,j)/9.81/psis) &
7103 soiliqw(l)=max(0.,soiliqw(l))
7104 soiliqw(l)=min(soiliqw(l),smois(i,l,j))
7105 sh2o(i,l,j)=soiliqw(l)
7106 smfr3d(i,l,j)=(smois(i,l,j)-soiliqw(l))/RIW
7110 sh2o(i,l,j)=smois(i,l,j)
7115 !-- for water ISLTYP=14
7129 END SUBROUTINE ruclsminit
7131 !-----------------------------------------------------------------
7132 ! SUBROUTINE RUCLSM_PARM_INIT
7133 !-----------------------------------------------------------------
7135 ! character*9 :: MMINLU, MMINSL
7137 ! MMINLU='MODIS-RUC'
7140 ! call RUCLSM_SOILVEGPARM( MMINLU, MMINSL)
7142 !-----------------------------------------------------------------
7143 ! END SUBROUTINE RUCLSM_PARM_INIT
7144 !-----------------------------------------------------------------
7146 !-----------------------------------------------------------------
7147 SUBROUTINE RUCLSM_SOILVEGPARM( MMINLURUC, MMINSL)
7148 !-----------------------------------------------------------------
7152 integer :: LUMATCH, IINDEX, LC, NUM_SLOPE
7154 INTEGER , PARAMETER :: OPEN_OK = 0
7156 character*8 :: MMINLURUC, MMINSL
7157 character*128 :: mess , message, vege_parm_string
7158 logical, external :: wrf_dm_on_monitor
7161 !-----SPECIFY VEGETATION RELATED CHARACTERISTICS :
7162 ! ALBBCK: SFC albedo (in percentage)
7163 ! Z0: Roughness length (m)
7165 ! PC: Plant coefficient for transpiration function
7166 ! -- the rest of the parameters are read in but not used currently
7167 ! SHDFAC: Green vegetation fraction (in percentage)
7168 ! Note: The ALBEDO, Z0, and SHDFAC values read from the following table
7169 ! ALBEDO, amd Z0 are specified in LAND-USE TABLE; and SHDFAC is
7170 ! the monthly green vegetation data
7171 ! CMXTBL: MAX CNPY Capacity (m)
7172 ! RSMIN: Mimimum stomatal resistance (s m-1)
7173 ! RSMAX: Max. stomatal resistance (s m-1)
7174 ! RGL: Parameters used in radiation stress function
7175 ! HS: Parameter used in vapor pressure deficit functio
7176 ! TOPT: Optimum transpiration air temperature. (K)
7177 ! CMCMAX: Maximum canopy water capacity
7178 ! CFACTR: Parameter used in the canopy inteception calculati
7179 ! SNUP: Threshold snow depth (in water equivalent m) that
7180 ! implies 100% snow cover
7181 ! LAI: Leaf area index (dimensionless)
7182 ! MAXALB: Upper bound on maximum albedo over deep snow
7184 !-----READ IN VEGETAION PROPERTIES FROM VEGPARM.TBL
7187 IF ( wrf_dm_on_monitor() ) THEN
7189 OPEN(19, FILE='VEGPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
7190 IF(ierr .NE. OPEN_OK ) THEN
7191 WRITE(message,FMT='(A)') &
7192 'module_sf_ruclsm.F: soil_veg_gen_parm: failure opening VEGPARM.TBL'
7193 CALL wrf_error_fatal ( message )
7196 WRITE ( mess, * ) 'INPUT VEGPARM FOR ',MMINLURUC
7197 CALL wrf_message( mess )
7202 READ (19,'(A)') vege_parm_string
7204 READ (19,2000,END=2002)LUTYPE
7205 READ (19,*)LUCATS,IINDEX
7207 WRITE( mess , * ) 'VEGPARM FOR ',LUTYPE,' FOUND', LUCATS,' CATEGORIES'
7208 CALL wrf_message( mess )
7210 IF(LUTYPE.NE.MMINLURUC)THEN ! Skip over the undesired table
7211 write ( mess , * ) 'Skipping ', LUTYPE, ' table'
7212 CALL wrf_message( mess )
7216 inner : DO ! Find the next "Vegetation Parameters"
7217 READ (19,'(A)',END=2002) vege_parm_string
7218 IF (TRIM(vege_parm_string) .EQ. "Vegetation Parameters") THEN
7224 write ( mess , * ) 'Found ', LUTYPE, ' table'
7225 CALL wrf_message( mess )
7226 EXIT outer ! Found the table, read the data
7231 IF (LUMATCH == 1) then
7232 write ( mess , * ) 'Reading ',LUTYPE,' table'
7233 CALL wrf_message( mess )
7235 READ (19,*)IINDEX,ALBTBL(LC),Z0TBL(LC),LEMITBL(LC),PCTBL(LC), &
7236 SHDTBL(LC),IFORTBL(LC),RSTBL(LC),RGLTBL(LC), &
7237 HSTBL(LC),SNUPTBL(LC),LAITBL(LC),MAXALB(LC)
7241 READ (19,*)TOPT_DATA
7243 READ (19,*)CMCMAX_DATA
7245 READ (19,*)CFACTR_DATA
7247 READ (19,*)RSMAX_DATA
7255 READ (19,*,iostat=ierr)URBAN
7256 if ( ierr /= 0 ) call wrf_message ( "-------- VEGPARM.TBL READ ERROR --------")
7257 if ( ierr /= 0 ) call wrf_message ( "Problem read URBAN from VEGPARM.TBL")
7258 if ( ierr /= 0 ) call wrf_message ( " -- Use updated version of VEGPARM.TBL ")
7259 if ( ierr /= 0 ) call wrf_error_fatal ( "Problem read URBAN from VEGPARM.TBL")
7266 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
7267 print *,' LEMITBL, PCTBL, Z0TBL, LAITBL --->', LEMITBL, PCTBL, Z0TBL, LAITBL
7271 IF (LUMATCH == 0) then
7272 CALL wrf_error_fatal ("Land Use Dataset '"//MMINLURUC//"' not found in VEGPARM.TBL.")
7277 CALL wrf_dm_bcast_string ( LUTYPE , 8 )
7278 CALL wrf_dm_bcast_integer ( LUCATS , 1 )
7279 CALL wrf_dm_bcast_integer ( IINDEX , 1 )
7280 CALL wrf_dm_bcast_integer ( LUMATCH , 1 )
7281 CALL wrf_dm_bcast_real ( ALBTBL , NLUS )
7282 CALL wrf_dm_bcast_real ( Z0TBL , NLUS )
7283 CALL wrf_dm_bcast_real ( LEMITBL , NLUS )
7284 CALL wrf_dm_bcast_real ( PCTBL , NLUS )
7285 CALL wrf_dm_bcast_real ( SHDTBL , NLUS )
7286 CALL wrf_dm_bcast_real ( IFORTBL , NLUS )
7287 CALL wrf_dm_bcast_real ( RSTBL , NLUS )
7288 CALL wrf_dm_bcast_real ( RGLTBL , NLUS )
7289 CALL wrf_dm_bcast_real ( HSTBL , NLUS )
7290 CALL wrf_dm_bcast_real ( SNUPTBL , NLUS )
7291 CALL wrf_dm_bcast_real ( LAITBL , NLUS )
7292 CALL wrf_dm_bcast_real ( MAXALB , NLUS )
7293 CALL wrf_dm_bcast_real ( TOPT_DATA , 1 )
7294 CALL wrf_dm_bcast_real ( CMCMAX_DATA , 1 )
7295 CALL wrf_dm_bcast_real ( CFACTR_DATA , 1 )
7296 CALL wrf_dm_bcast_real ( RSMAX_DATA , 1 )
7297 CALL wrf_dm_bcast_integer ( BARE , 1 )
7298 CALL wrf_dm_bcast_integer ( NATURAL , 1 )
7299 CALL wrf_dm_bcast_integer ( CROP , 1 )
7300 CALL wrf_dm_bcast_integer ( URBAN , 1 )
7303 !-----READ IN SOIL PROPERTIES FROM SOILPARM.TBL
7305 IF ( wrf_dm_on_monitor() ) THEN
7306 OPEN(19, FILE='SOILPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
7307 IF(ierr .NE. OPEN_OK ) THEN
7308 WRITE(message,FMT='(A)') &
7309 'module_sf_ruclsm.F: soil_veg_gen_parm: failure opening SOILPARM.TBL'
7310 CALL wrf_error_fatal ( message )
7313 WRITE(mess,*) 'INPUT SOIL TEXTURE CLASSIFICATION = ',MMINSL
7314 CALL wrf_message( mess )
7319 READ (19,2000,END=2003)SLTYPE
7320 READ (19,*)SLCATS,IINDEX
7321 IF(SLTYPE.NE.MMINSL)THEN
7323 READ (19,*) IINDEX,BB(LC),DRYSMC(LC),HC(LC),MAXSMC(LC),&
7324 REFSMC(LC),SATPSI(LC),SATDK(LC), SATDW(LC), &
7329 READ (19,2000,END=2003)SLTYPE
7330 READ (19,*)SLCATS,IINDEX
7332 IF(SLTYPE.EQ.MMINSL)THEN
7333 WRITE( mess , * ) 'SOIL TEXTURE CLASSIFICATION = ',SLTYPE,' FOUND', &
7334 SLCATS,' CATEGORIES'
7335 CALL wrf_message ( mess )
7338 IF(SLTYPE.EQ.MMINSL)THEN
7340 READ (19,*) IINDEX,BB(LC),DRYSMC(LC),HC(LC),MAXSMC(LC),&
7341 REFSMC(LC),SATPSI(LC),SATDK(LC), SATDW(LC), &
7351 CALL wrf_dm_bcast_integer ( LUMATCH , 1 )
7352 CALL wrf_dm_bcast_string ( SLTYPE , 8 )
7353 CALL wrf_dm_bcast_string ( MMINSL , 8 ) ! since this is reset above, see oct2 ^
7354 CALL wrf_dm_bcast_integer ( SLCATS , 1 )
7355 CALL wrf_dm_bcast_integer ( IINDEX , 1 )
7356 CALL wrf_dm_bcast_real ( BB , NSLTYPE )
7357 CALL wrf_dm_bcast_real ( DRYSMC , NSLTYPE )
7358 CALL wrf_dm_bcast_real ( HC , NSLTYPE )
7359 CALL wrf_dm_bcast_real ( MAXSMC , NSLTYPE )
7360 CALL wrf_dm_bcast_real ( REFSMC , NSLTYPE )
7361 CALL wrf_dm_bcast_real ( SATPSI , NSLTYPE )
7362 CALL wrf_dm_bcast_real ( SATDK , NSLTYPE )
7363 CALL wrf_dm_bcast_real ( SATDW , NSLTYPE )
7364 CALL wrf_dm_bcast_real ( WLTSMC , NSLTYPE )
7365 CALL wrf_dm_bcast_real ( QTZ , NSLTYPE )
7367 IF(LUMATCH.EQ.0)THEN
7368 CALL wrf_message( 'SOIl TEXTURE IN INPUT FILE DOES NOT ' )
7369 CALL wrf_message( 'MATCH SOILPARM TABLE' )
7370 CALL wrf_error_fatal ( 'INCONSISTENT OR MISSING SOILPARM FILE' )
7374 !-----READ IN GENERAL PARAMETERS FROM GENPARM.TBL
7376 IF ( wrf_dm_on_monitor() ) THEN
7377 OPEN(19, FILE='GENPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
7378 IF(ierr .NE. OPEN_OK ) THEN
7379 WRITE(message,FMT='(A)') &
7380 'module_sf_ruclsm.F: soil_veg_gen_parm: failure opening GENPARM.TBL'
7381 CALL wrf_error_fatal ( message )
7386 READ (19,*) NUM_SLOPE
7391 READ (19,*)SLOPE_DATA(LC)
7395 READ (19,*)SBETA_DATA
7397 READ (19,*)FXEXP_DATA
7399 READ (19,*)CSOIL_DATA
7401 READ (19,*)SALP_DATA
7403 READ (19,*)REFDK_DATA
7405 READ (19,*)REFKDT_DATA
7407 READ (19,*)FRZK_DATA
7409 READ (19,*)ZBOT_DATA
7411 READ (19,*)CZIL_DATA
7413 READ (19,*)SMLOW_DATA
7415 READ (19,*)SMHIGH_DATA
7419 CALL wrf_dm_bcast_integer ( NUM_SLOPE , 1 )
7420 CALL wrf_dm_bcast_integer ( SLPCATS , 1 )
7421 CALL wrf_dm_bcast_real ( SLOPE_DATA , NSLOPE )
7422 CALL wrf_dm_bcast_real ( SBETA_DATA , 1 )
7423 CALL wrf_dm_bcast_real ( FXEXP_DATA , 1 )
7424 CALL wrf_dm_bcast_real ( CSOIL_DATA , 1 )
7425 CALL wrf_dm_bcast_real ( SALP_DATA , 1 )
7426 CALL wrf_dm_bcast_real ( REFDK_DATA , 1 )
7427 CALL wrf_dm_bcast_real ( REFKDT_DATA , 1 )
7428 CALL wrf_dm_bcast_real ( FRZK_DATA , 1 )
7429 CALL wrf_dm_bcast_real ( ZBOT_DATA , 1 )
7430 CALL wrf_dm_bcast_real ( CZIL_DATA , 1 )
7431 CALL wrf_dm_bcast_real ( SMLOW_DATA , 1 )
7432 CALL wrf_dm_bcast_real ( SMHIGH_DATA , 1 )
7435 !-----------------------------------------------------------------
7436 END SUBROUTINE RUCLSM_SOILVEGPARM
7437 !-----------------------------------------------------------------
7440 SUBROUTINE SOILIN (ISLTYP, DQM, REF, PSIS, QMIN, BCLH )
7442 !--- soiltyp classification according to STATSGO(nclasses=16)
7445 ! 2 LOAMY SAND LOAMY SAND
7446 ! 3 SANDY LOAM SANDY LOAM
7447 ! 4 SILT LOAM SILTY LOAM
7450 ! 7 SANDY CLAY LOAM SANDY CLAY LOAM
7451 ! 8 SILTY CLAY LOAM SILTY CLAY LOAM
7452 ! 9 CLAY LOAM CLAY LOAM
7453 ! 10 SANDY CLAY SANDY CLAY
7454 ! 11 SILTY CLAY SILTY CLAY
7455 ! 12 CLAY LIGHT CLAY
7456 ! 13 ORGANIC MATERIALS LOAM
7459 ! Bedrock is reclassified as class 14
7460 ! 16 OTHER (land-ice)
7461 ! extra classes from Fei Chen
7466 !----------------------------------------------------------------------
7467 integer, parameter :: nsoilclas=19
7469 integer, intent ( in) :: isltyp
7470 real, intent ( out) :: dqm,ref,qmin,psis
7472 REAL LQMA(nsoilclas),LREF(nsoilclas),LBCL(nsoilclas), &
7473 LPSI(nsoilclas),LQMI(nsoilclas)
7475 !-- LQMA Rawls et al.[1982]
7476 ! DATA LQMA /0.417, 0.437, 0.453, 0.501, 0.486, 0.463, 0.398,
7477 ! & 0.471, 0.464, 0.430, 0.479, 0.475, 0.439, 1.0, 0.20, 0.401/
7479 !-- Clapp, R. and G. Hornberger, Empirical equations for some soil
7480 ! hydraulic properties, Water Resour. Res., 14,601-604,1978.
7481 !-- Clapp et al. [1978]
7482 DATA LQMA /0.395, 0.410, 0.435, 0.485, 0.485, 0.451, 0.420, &
7483 0.477, 0.476, 0.426, 0.492, 0.482, 0.451, 1.0, &
7484 0.20, 0.435, 0.468, 0.200, 0.339/
7486 !-- Clapp et al. [1978]
7487 DATA LREF /0.174, 0.179, 0.249, 0.369, 0.369, 0.314, 0.299, &
7488 0.357, 0.391, 0.316, 0.409, 0.400, 0.314, 1., &
7489 0.1, 0.249, 0.454, 0.17, 0.236/
7491 !-- Carsel and Parrish [1988]
7492 DATA LQMI/0.045, 0.057, 0.065, 0.067, 0.034, 0.078, 0.10, &
7493 0.089, 0.095, 0.10, 0.070, 0.068, 0.078, 0.0, &
7494 0.004, 0.065, 0.020, 0.004, 0.008/
7496 !-- Clapp et al. [1978]
7497 DATA LPSI/0.121, 0.090, 0.218, 0.786, 0.786, 0.478, 0.299, &
7498 0.356, 0.630, 0.153, 0.490, 0.405, 0.478, 0.0, &
7499 0.121, 0.218, 0.468, 0.069, 0.069/
7501 !-- Clapp et al. [1978]
7502 DATA LBCL/4.05, 4.38, 4.90, 5.30, 5.30, 5.39, 7.12, &
7503 7.75, 8.52, 10.40, 10.40, 11.40, 5.39, 0.0, &
7504 4.05, 4.90, 11.55, 2.79, 2.79/
7507 DQM = LQMA(ISLTYP)- &
7510 PSIS = - LPSI(ISLTYP)
7514 END SUBROUTINE SOILIN
7516 END MODULE module_sf_ruclsm