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
46 !-- options for snow conductivity: 1 - constant, 2 - Sturm et al.,1997
47 ! integer, parameter :: isncond_opt = 1
49 integer, parameter :: isncond_opt=2
51 !-- Snow fraction options
52 !-- option 1: original formulation using threshold snow depth to compute snow fraction
53 !integer, parameter :: isncovr_opt = 1 (default)
54 !-- option 2: the tanh formulation from Niu,G.-Y.,and Yang,Z.-L., 2007,JGR,DOI:10.1029/2007JD008674.
55 !integer, parameter :: isncovr_opt = 2
56 !-- option 3: the tanh formulation from Niu,G.-Y.,and Yang,Z with
57 ! vegetation-dependent parameters from Noah MP (personal communication with
59 !integer, parameter :: isncovr_opt = 3
60 !-- Values of parameters are scale-dependent, have to be tuned for a given application
61 !-- Tables below are for 21-class MODI-RUC (MODIFIED_IGBP_MODIS_NOAH_15s is used in HRRR and RRFS)
62 !-- for 3-km RRFS application
63 real, dimension(30), parameter :: sncovfac = &
64 & (/ 0.030, 0.030, 0.030, 0.030, 0.030, &
65 & 0.016, 0.016, 0.020, 0.020, 0.020, &
66 & 0.020, 0.014, 0.042, 0.026, 0.030, &
67 & 0.016, 0.030, 0.030, 0.030, 0.030, &
68 & 0.000, 0.000, 0.000, 0.000, 0.000, &
69 & 0.000, 0.000, 0.000, 0.000, 0.000 /)
70 real, dimension(30), parameter :: mfsno = &
71 & (/ 1.00, 1.00, 1.00, 1.00, 2.00, 2.00, &
72 & 2.00, 2.00, 2.00, 2.00, 2.00, 2.00, &
73 & 3.00, 3.00, 2.00, 2.00, 2.00, 2.00, &
74 & 2.00, 2.00, 0.00, 0.00, 0.00, 0.00, &
75 & 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 /)
78 integer, parameter :: isncovr_opt=2
83 !-----------------------------------------------------------------
84 subroutine lsmruc(spp_lsm, &
86 pattern_spp_lsm,field_sf, &
91 graupelncv,snowncv,rainncv, &
93 zs,rainbl,snow,snowh,snowc,frzfrac,frpcpn, &
94 rhosnf,precipfr, & ! pass it out to module_diagnostics
95 z3d,p8w,t3d,qv3d,qc3d,rho3d, & !p8w in [pa]
96 glw,gsw,emiss,chklowq, chs, &
97 flqc,flhc,mavail,canwat,vegfra,alb,znt, &
98 z0,snoalb,albbck,lai, & !new
99 mminlu, landusef, nlcat, mosaic_lu, &
100 mosaic_soil, soilctop, nscat, & !new
101 qsfc,qsg,qvg,qcg,dew,soilt1,tsnav, &
102 tbot,ivgtyp,isltyp,xland, &
103 iswater,isice,xice,xice_threshold, &
104 cp,rovcp,g0,lv,stbolt, &
105 soilmois,sh2o,smavail,smmax, &
106 tso,soilt,hfx,qfx,lh, &
107 sfcrunoff,udrunoff,acrunoff,sfcexc, &
108 sfcevp,grdflx,snowfallac,acsnow,snom, &
109 smfr3d,keepfr3dflag, &
110 myj,shdmin,shdmax,rdlai2d, &
111 ids,ide, jds,jde, kds,kde, &
112 ims,ime, jms,jme, kms,kme, &
113 its,ite, jts,jte, kts,kte )
114 !-----------------------------------------------------------------
116 !-----------------------------------------------------------------
118 ! the ruc lsm model is described in:
119 ! Smirnova, T.G., J.M. Brown, and S.G. Benjamin, 1997:
120 ! performance of different soil model configurations in simulating
121 ! ground surface temperature and surface fluxes.
122 ! mon. wea. rev. 125, 1870-1884.
123 ! Smirnova, T.G., J.M. Brown, and D. Kim, 2000: parameterization of
124 ! cold-season processes in the maps land-surface scheme.
125 ! j. geophys. res. 105, 4077-4086.
126 !-----------------------------------------------------------------
127 !-- dt time step (second)
128 ! ktau - number of time step
129 ! nsl - number of soil layers
130 ! nzs - number of levels in soil
131 ! zs - depth of soil levels (m)
132 !-- rainbl - accumulated rain in [mm] between the pbl calls
133 !-- rainncv one time step grid scale precipitation (mm/step)
134 ! snow - snow water equivalent [mm]
135 ! frazfrac - fraction of frozen precipitation
136 !-- precipfr (mm) - time step frozen precipitation
137 !-- snowc flag indicating snow coverage (1 for snow cover)
139 !-- p8w 3d pressure (pa)
140 !-- t3d temperature (k)
141 !-- qv3d 3d water vapor mixing ratio (kg/kg)
142 ! qc3d - 3d cloud water mixing ratio (kg/kg)
143 ! rho3d - 3d air density (kg/m^3)
144 !-- glw downward long wave flux at ground surface (w/m^2)
145 !-- gsw absorbed short wave flux at ground surface (w/m^2)
146 !-- emiss surface emissivity (between 0 and 1)
147 ! flqc - surface exchange coefficient for moisture (kg/m^2/s)
148 ! flhc - surface exchange coefficient for heat [w/m^2/s/degreek]
149 ! sfcexc - surface exchange coefficient for heat [m/s]
150 ! canwat - canopy moisture content (mm)
151 ! vegfra - vegetation fraction (between 0 and 100)
152 ! alb - surface albedo (between 0 and 1)
153 ! snoalb - maximum snow albedo (between 0 and 1)
154 ! albbck - snow-free albedo (between 0 and 1)
155 ! znt - roughness length [m]
156 !-- tbot soil temperature at lower boundary (k)
157 ! ivgtyp - usgs vegetation type (24 classes)
158 ! isltyp - stasgo soil type (16 classes)
159 !-- xland land mask (1 for land, 2 for water)
160 !-- cp heat capacity at constant pressure for dry air (j/kg/k)
161 !-- g0 acceleration due to gravity (m/s^2)
162 !-- lv latent heat of melting (j/kg)
163 !-- stbolt stefan-boltzmann constant (w/m^2/k^4)
164 ! soilmois - soil moisture content (volumetric fraction)
165 ! tso - soil temp (k)
166 !-- soilt surface temperature (k)
167 !-- hfx upward heat flux at the surface (w/m^2)
168 !-- qfx upward moisture flux at the surface (kg/m^2/s)
169 !-- lh upward latent heat flux (w/m^2)
170 ! sfcrunoff - ground surface runoff [mm]
171 ! udrunoff - underground runoff [mm]
172 ! acrunoff - run-total surface runoff [mm]
173 ! sfcevp - total evaporation in [kg/m^2]
174 ! grdflx - soil heat flux (w/m^2: negative, if downward from surface)
175 ! snowfallac - run-total snowfall accumulation [m]
176 ! acsnow - run-toral swe of snowfall [mm]
177 !-- chklowq - is either 0 or 1 (so far set equal to 1).
178 !-- used only in myjpbl.
179 !-- tice - sea ice temperture (c)
180 !-- rhosice - sea ice density (kg m^-3)
181 !-- capice - sea ice volumetric heat capacity (j/m^3/k)
182 !-- thdifice - sea ice thermal diffusivity (m^2/s)
184 !-- ims start index for i in memory
185 !-- ime end index for i in memory
186 !-- jms start index for j in memory
187 !-- jme end index for j in memory
188 !-- kms start index for k in memory
189 !-- kme end index for k in memory
190 !-------------------------------------------------------------------------
191 ! integer, parameter :: nzss=5
192 ! integer, parameter :: nddzs=2*(nzss-2)
194 integer, parameter :: nvegclas=24+3
196 real, intent(in ) :: dt
197 logical, intent(in ) :: myj,frpcpn
198 integer, intent(in ) :: spp_lsm
199 integer, intent(in ) :: nlcat, nscat, mosaic_lu, mosaic_soil
200 integer, intent(in ) :: ktau, nsl, isice, iswater, &
201 ims,ime, jms,jme, kms,kme, &
202 ids,ide, jds,jde, kds,kde, &
203 its,ite, jts,jte, kts,kte
206 real, dimension( ims:ime, kms:kme, jms:jme ),optional:: pattern_spp_lsm
207 real, dimension( ims:ime, kms:kme, jms:jme ),optional:: field_sf
209 real, dimension( ims:ime, 1 :nsl, jms:jme ) :: field_sf_loc
211 real, dimension( ims:ime, kms:kme, jms:jme ) , &
212 intent(in ) :: qv3d, &
219 real, dimension( ims:ime , jms:jme ), &
220 intent(in ) :: rainbl, &
233 real, dimension( ims:ime , jms:jme ), &
234 intent(inout ) :: vegfra
238 real, optional, dimension( ims:ime , jms:jme ), &
239 intent(in ) :: graupelncv, &
242 real, dimension( ims:ime , jms:jme ), &
243 intent(in ) :: lakemask
244 integer, intent(in ) :: lakemodel
247 real, dimension( ims:ime , jms:jme ), intent(in ):: shdmax
248 real, dimension( ims:ime , jms:jme ), intent(in ):: shdmin
249 logical, intent(in) :: rdlai2d
251 real, dimension( 1:nsl), intent(in ) :: zs
253 real, dimension( ims:ime , jms:jme ), &
268 real, dimension( ims:ime , jms:jme ), &
272 integer, dimension( ims:ime , jms:jme ), &
273 intent(in ) :: ivgtyp, &
275 character(len=*), intent(in ) :: mminlu
276 real, dimension( ims:ime , 1:nlcat, jms:jme ), intent(in):: landusef
277 real, dimension( ims:ime , 1:nscat, jms:jme ), intent(in):: soilctop
279 real, intent(in ) :: cp,rovcp,g0,lv,stbolt,xice_threshold
281 real, dimension( ims:ime , 1:nsl, jms:jme ) , &
282 intent(inout) :: soilmois,sh2o,tso
284 real, dimension( ims:ime, jms:jme ) , &
285 intent(inout) :: soilt, &
305 real, dimension( ims:ime, jms:jme ) , &
306 intent(inout) :: smavail, &
309 real, dimension( its:ite, jts:jte ) :: &
329 ! energy and water budget variables:
330 real, dimension( its:ite, jts:jte ) :: &
340 real, dimension( ims:ime, 1:nsl, jms:jme) &
344 real, dimension( ims:ime, jms:jme ), intent(out) :: &
345 rhosnf, & !rho of snowfall
346 precipfr, & ! time-step frozen precip
348 !--- soil/snow properties
376 real, dimension(1:nsl) :: zsmain, &
380 real, dimension(1:2*(nsl-2)) :: dtdzs
382 real, dimension(1:5001) :: tbq
385 real, dimension( 1:nsl ) :: soilm1d, &
391 real, dimension( 1:nsl ) :: keepfr
393 real, dimension( 1:nlcat ) :: lufrac
394 real, dimension( 1:nscat ) :: soilfrac
422 real :: cq,r61,r273,arp,brp,x,evs,eis
423 real :: cropfr, cropsm, newsm, factor
425 real :: meltfactor, ac,as, wb
427 integer :: iland,isoil,iforest
429 integer :: i,j,k,nzs,nzs1,nddzs
430 integer :: k1,l,k2,kp,km
431 character (len=132) :: message
433 real,dimension(ims:ime,1:nsl,jms:jme) :: rstoch
435 real,dimension(ims:ime,jms:jme)::emisso,vegfrao,albo,snoalbo
436 real,dimension(its:ite,jts:jte)::emisslo
438 !-----------------------------------------------------------------
450 rstoch(i,k,j) = pattern_spp_lsm(i,k,j)
451 field_sf_loc(i,k,j)=field_sf(i,k,j)
457 !---- table tbq is for resolution of balance equation in vilka
461 arp=77455.*41.9/461.525
466 evs=exp(17.67*(cq-273.15)/(cq-29.65))
467 eis=exp(22.514-6.15e3/cq)
468 if(cq.ge.273.15) then
477 !--- initialize soil/vegetation parameters
478 #if ( NMM_CORE == 1 )
486 keepfr3dflag(i,k,j)=0.
488 !--- initializing snow fraction, thereshold = 32 mm of snow water or ~100 mm of snow height
489 if((soilt1(i,j) .lt. 170.) .or. (soilt1(i,j) .gt.400.)) then
490 if(snowc(i,j).gt.0.) then
491 soilt1(i,j)=0.5*(soilt(i,j)+tso(i,1,j))
492 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
493 write ( message , fmt='(a,f8.3,2i6)' ) &
494 'temperature inside snow is initialized in ruclsm ', soilt1(i,j),i,j
495 call wrf_debug ( 0 , message )
498 soilt1(i,j) = tso(i,1,j)
501 !-- temperature inside snow is initialized
502 tsnav(i,j) =0.5*(soilt(i,j)+tso(i,1,j))-273.15
503 patmb=p8w(i,kms,j)*1.e-2
504 qsg (i,j) = qsn(soilt(i,j),tbq)/patmb
505 if((qcg(i,j) < 0.) .or. (qcg(i,j) > 0.1)) then
506 qcg (i,j) = qc3d(i,1,j)
507 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
508 write ( message , fmt='(a,3f8.3,2i6)' ) &
509 'qvg is initialized in ruclsm ', qvg(i,j),mavail(i,j),qsg(i,j),i,j
513 if((qvg(i,j) .le. 0.) .or. (qvg(i,j) .gt.0.1)) then
514 qvg (i,j) = qsg(i,j)*mavail(i,j)
515 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
516 write ( message , fmt='(a,3f8.3,2i6)' ) &
517 'qvg is initialized in ruclsm ', qvg(i,j),mavail(i,j),qsg(i,j),i,j
518 call wrf_debug ( 0 , message )
521 qsfc(i,j) = qvg(i,j)/(1.+qvg(i,j))
526 rhosnf(i,j) = -1.e3 ! non-zero flag
539 waterbudget(i,j) = 0.
540 acwaterbudget(i,j) = 0.
544 ! for ruc lsm chklowq needed for myjpbl should
545 ! 1 because is actual specific humidity at the surface, and
546 ! not the saturation value
567 !-----------------------------------------------------------------
581 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
582 print *,' in lsmruc ','ims,ime,jms,jme,its,ite,jts,jte,nzs', &
583 ims,ime,jms,jme,its,ite,jts,jte,nzs
584 print *,' ivgtyp, isltyp ', ivgtyp(i,j),isltyp(i,j)
585 print *,' mavail ', mavail(i,j)
586 print *,' soilt,qvg,p8w',soilt(i,j),qvg(i,j),p8w(i,1,j)
587 print *, 'lsmruc, i,j,xland, qfx,hfx from sfclay',i,j,xland(i,j), &
589 print *, ' gsw, glw =',gsw(i,j),glw(i,j)
590 print *, 'soilt, tso start of time step =',soilt(i,j),(tso(i,k,j),k=1,nsl)
591 print *, 'soilmois start of time step =',(soilmois(i,k,j),k=1,nsl)
592 print *, 'smfrozen start of time step =',(smfr3d(i,k,j),k=1,nsl)
593 print *, ' i,j=, after sfclay chs,flhc ',i,j,chs(i,j),flhc(i,j)
594 print *, 'lsmruc, ivgtyp,isltyp,alb = ', ivgtyp(i,j),isltyp(i,j),alb(i,j),i,j
595 print *, 'lsmruc i,j,dt,rainbl =',i,j,dt,rainbl(i,j)
596 print *, 'xland ---->, ivgtype,isoiltyp,i,j',xland(i,j),ivgtyp(i,j),isltyp(i,j),i,j
603 qvatm = qv3d(i,kms,j)
604 qcatm = qc3d(i,kms,j)
605 patm = p8w(i,kms,j)*1.e-5
606 !-- z3d(1) is thickness between first full sigma level and the surface,
607 !-- but first mass level is at the half of the first sigma level
608 !-- (u and v are also at the half of first sigma level)
609 conflx = z3d(i,kms,j)*0.5
611 ! -- initialize snow, graupel and ice fractions in frozen precip
618 prcpncliq = rainncv(i,j)*(1.-frzfrac(i,j))
619 prcpncfr = rainncv(i,j)*frzfrac(i,j)
620 !- apply the same frozen precipitation fraction to convective precip
621 !tgs - 31 mar17 - add safety temperature check in case thompson mp produces
622 ! frozen precip at t > 273.
623 if(frzfrac(i,j) > 0..and. tabs < 273.) then
624 prcpculiq = max(0.,(rainbl(i,j)-rainncv(i,j))*(1.-frzfrac(i,j)))
625 prcpcufr = max(0.,(rainbl(i,j)-rainncv(i,j))*frzfrac(i,j))
628 prcpcufr = max(0.,(rainbl(i,j)-rainncv(i,j)))
632 prcpculiq = max(0.,(rainbl(i,j)-rainncv(i,j)))
635 !--- 1*e-3 is to convert from mm/s to m/s
636 prcpms = (prcpncliq + prcpculiq)/dt*1.e-3
637 newsnms = (prcpncfr + prcpcufr)/dt*1.e-3
639 if ( present( graupelncv ) ) then
640 graupamt = graupelncv(i,j)
645 if((prcpncfr + prcpcufr) > 0.) then
646 ! -- calculate snow, graupel and ice fractions in falling frozen precip
647 snowrat=min(1.,max(0.,snowncv(i,j)/(prcpncfr + prcpcufr)))
648 grauprat=min(1.,max(0.,graupamt/(prcpncfr + prcpcufr)))
649 icerat=min(1.,max(0.,(prcpncfr-snowncv(i,j)-graupamt) &
650 /(prcpncfr + prcpcufr)))
651 curat=min(1.,max(0.,(prcpcufr/(prcpncfr + prcpcufr))))
654 prcpms = (rainbl(i,j)/dt*1.e-3)*(1-frzfrac(i,j))
655 newsnms = (rainbl(i,j)/dt*1.e-3)*frzfrac(i,j)
656 if(newsnms == 0.) then
659 snowrat = min(1.,newsnms/(newsnms+prcpms))
664 if (tabs.le.273.15) then
666 newsnms = rainbl(i,j)/dt*1.e-3
667 !-- here no info about constituents of frozen precipitation,
668 !-- suppose it is all snow
671 prcpms = rainbl(i,j)/dt*1.e-3
676 ! -- save time-step water equivalent of frozen precipitation in precipfr array to be used in
678 precipfr(i,j) = newsnms * dt *1.e3
684 !--- convert exchange coeff qkms to [m/s]
685 qkms=flqc(i,j)/rho/mavail(i,j)
686 ! tkms=flhc(i,j)/rho/cp
687 tkms=flhc(i,j)/rho/(cp*(1.+0.84*qvatm)) ! mynnsfc uses cpm
689 !--- convert incoming snow and canwat from mm to m
692 canwatr=canwat(i,j)*1.e-3
695 rhosnfall=rhosnf(i,j)
703 zshalf(k)=0.5*(zsmain(k-1) + zsmain(k))
707 lufrac(k) = landusef(i,k,j)
710 soilfrac(k) = soilctop(i,k,j)
713 !------------------------------------------------------------
714 !----- ddzs and dsdz1 are for implicit solution of soil eqns.
715 !-------------------------------------------------------------
718 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
719 print *,' dt,nzs1, zsmain, zshalf --->', dt,nzs1,zsmain,zshalf
725 x=dt/2./(zshalf(k+1)-zshalf(k))
726 dtdzs(k1)=x/(zsmain(k)-zsmain(k-1))
728 dtdzs(k2)=x/(zsmain(k+1)-zsmain(k))
734 !--- constants used in johansen soil thermal
735 !--- conductivity method
741 !***********************************************************************
742 !--- constants for snow density calculations c1sn and c2sn
747 !***********************************************************************
749 nroot= 4 ! levels in root layer
752 if(snow(i,j).gt.0. .and. snowh(i,j).gt.0.) then
753 rhosn = snow(i,j)/snowh(i,j)
758 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
759 if(ktau.eq.1 .and.(i.eq.358.and.j.eq.260)) &
760 print *,'before soilvegin - z0,znt(195,254)',z0(i,j),znt(i,j)
762 !--- initializing soil and surface properties
763 call soilvegin ( mosaic_lu, mosaic_soil,soilfrac,nscat,shdmin(i,j),shdmax(i,j),&
764 nlcat,iland,isoil,iswater,myj,iforest,lufrac,vegfra(i,j), &
765 emissl(i,j),pc(i,j),znt(i,j),lai(i,j),rdlai2d, &
766 qwrtz,rhocs,bclh,dqm,ksat,psis,qmin,ref,wilt,i,j )
767 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
768 if(ktau.eq.1 .and.(i.eq.358.and.j.eq.260)) &
769 print *,'after soilvegin - z0,znt(375,254),lai(375,254)',z0(i,j),znt(i,j),lai(i,j)
771 if(ktau.eq.1 .and. (i.eq.358.and.j.eq.260)) then
772 print *,'nlcat,iland,lufrac,emissl(i,j),pc(i,j),znt(i,j),lai(i,j)', &
773 nlcat,iland,lufrac,emissl(i,j),pc(i,j),znt(i,j),lai(i,j),i,j
774 print *,'nscat,soilfrac,qwrtz,rhocs,bclh,dqm,ksat,psis,qmin,ref,wilt',&
775 nscat,soilfrac,qwrtz,rhocs,bclh,dqm,ksat,psis,qmin,ref,wilt,i,j
779 cn=cfactr_data ! exponent
780 sat = 5.e-4 ! units [m]
782 !-- definition of number of soil levels in the rooting zone
783 if(iforest.gt.2) then
784 !---- all vegetation types except evergreen and mixed forests
785 !18apr08 - define meltfactor for egglston melting limit:
786 ! for open areas factor is 2, and for forests - factor is 0.85
787 ! this will make limit on snow melting smaller and let snow stay
788 ! longer in the forests.
792 if(zsmain(k).ge.0.4) then
798 !---- evergreen and mixed forests
799 !18apr08 - define meltfactor
801 ! 28 march 11 - previously used value of metfactor= 1.5 needs to be further reduced
802 ! to compensate for low snow albedos in the forested areas.
803 ! melting rate in forests will reduce.
807 if(zsmain(k).ge.1.1) then
816 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
817 print *,' znt, lai, vegfra, sat, emis, pc --->', &
818 znt(i,j),lai(i,j),vegfra(i,j),sat,emissl(i,j),pc(i,j)
819 print *,' zs, zsmain, zshalf, conflx, cn, sat, --->', zs,zsmain,zshalf,conflx,cn,sat
820 print *,'nroot, meltfactor, iforest, ivgtyp, i,j ', nroot,meltfactor,iforest,ivgtyp(i,j),i,j
824 if(lakemodel==1. .and. lakemask(i,j)==1.) goto 2999
828 if((xland(i,j)-1.5).ge.0.)then
840 patmb=p8w(i,1,j)*1.e-2
841 qvg (i,j) = qsn(soilt(i,j),tbq)/patmb
842 qsfc(i,j) = qvg(i,j)/(1.+qvg(i,j))
844 q2sat=qsn(tabs,tbq)/patmb
849 tso(i,k,j)= soilt(i,j)
852 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
853 print*,' water point, i=',i, &
854 'j=',j, 'soilt=', soilt(i,j)
859 ! land point or sea ice
860 if(xice(i,j).ge.xice_threshold) then
866 if(seaice(i,j).gt.0.5)then
868 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
869 print*,' sea-ice at water point, i=',i, &
883 patmb=p8w(i,1,j)*1.e-2
884 qvg (i,j) = qsn(soilt(i,j),tbq)/patmb
886 qsfc(i,j) = qvg(i,j)/(1.+qvg(i,j))
892 keepfr3dflag(i,k,j) = 0.
893 tso(i,k,j) = min(271.4,tso(i,k,j))
897 ! attention!!!! ruc lsm uses soil moisture content minus residual (minimum
898 ! or dry soil moisture content for a given soil type) as a state variable.
901 ! soilm1d - soil moisture content minus residual [m**3/m**3]
902 soilm1d (k) = min(max(0.,soilmois(i,k,j)-qmin),dqm)
903 tso1d (k) = tso(i,k,j)
904 soiliqw (k) = min(max(0.,sh2o(i,k,j)-qmin),soilm1d(k))
905 soilice (k) =(soilm1d (k) - soiliqw (k))/0.9
909 smfrkeep(k) = smfr3d(i,k,j)
910 keepfr (k) = keepfr3dflag(i,k,j)
913 lmavail(i,j)=max(0.00001,min(1.,soilm1d(1)/(ref-qmin)))
915 #if ( NMM_CORE == 1 )
922 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
923 print *,'land, i,j,tso1d,soilm1d,patm,tabs,qvatm,qcatm,rho', &
924 i,j,tso1d,soilm1d,patm,tabs,qvatm,qcatm,rho
925 print *,'conflx =',conflx
926 print *,'smfrkeep,keepfr ',smfrkeep,keepfr
931 smtotold(i,j)=smtotold(i,j)+(qmin+soilm1d(k))* &
932 (zshalf(k+1)-zshalf(k))
935 smtotold(i,j)=smtotold(i,j)+(qmin+soilm1d(nzs))* &
936 (zsmain(nzs)-zshalf(nzs))
938 canwatold(i,j) = canwatr
939 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
940 print *,'before sfctmp, spp_lsm, rstoch, field_sf_loc', &
941 i,j,spp_lsm,(rstoch(i,k,j),k=1,nzs),(field_sf_loc(i,k,j),k=1,nzs)
943 !-----------------------------------------------------------------
944 call sfctmp (spp_lsm,rstoch(i,:,j),field_sf_loc(i,:,j), &
945 dt,ktau,conflx,i,j, &
947 nzs,nddzs,nroot,meltfactor, & !added meltfactor
948 iland,isoil,xland(i,j),ivgtyp(i,j),isltyp(i,j), &
949 prcpms, newsnms,snwe,snhei,snowfrac, &
950 rhosn,rhonewsn,rhosnfall, &
951 snowrat,grauprat,icerat,curat, &
952 patm,tabs,qvatm,qcatm,rho, &
953 glw(i,j),gsw(i,j),emissl(i,j), &
954 qkms,tkms,pc(i,j),lmavail(i,j), &
955 canwatr,vegfra(i,j),alb(i,j),znt(i,j), &
956 snoalb(i,j),albbck(i,j),lai(i,j), & !new
957 myj,seaice(i,j),isice, &
958 !--- soil fixed fields
960 rhocs,dqm,qmin,ref, &
961 wilt,psis,bclh,ksat, &
962 sat,cn,zsmain,zshalf,dtdzs,dtdzs2,tbq, &
964 cp,rovcp,g0,lv,stbolt,cw,c1sn,c2sn, &
966 !--- output variables
967 snweprint,snheiprint,rsm, &
968 soilm1d,tso1d,smfrkeep,keepfr, &
969 soilt(i,j),soilt1(i,j),tsnav(i,j),dew(i,j), &
970 qvg(i,j),qsg(i,j),qcg(i,j),smelt(i,j), &
971 snoh(i,j),snflx(i,j),snom(i,j),snowfallac(i,j), &
972 acsnow(i,j),edir(i,j),ec(i,j),ett(i,j),qfx(i,j), &
973 lh(i,j),hfx(i,j),sflx(i,j),sublim(i,j), &
974 evapl(i,j),prcpl(i,j),budget(i,j),runoff1(i,j), &
975 runoff2(i,j),soilice,soiliqw,infiltrp,smf(i,j))
976 !-----------------------------------------------------------------
978 ! irrigation: fraction of cropland category in the grid box should not have soil moisture below
979 ! wilting point during the growing season.
980 ! let's keep soil moisture 10% above wilting point for the fraction of grid box under
982 ! this change violates lsm moisture budget, but
983 ! can be considered as a compensation for irrigation not included into lsm.
985 if(mosaic_lu == 1) then
986 ! greenness factor: between 0 for min greenness and 1 for max greenness.
987 factor = max(0.,min(1.,(vegfra(i,j)-shdmin(i,j))/max(1.,(shdmax(i,j)-shdmin(i,j)))))
989 if ((lufrac(crop) > 0 .or. lufrac(natural) > 0.).and. factor > 0.75) then
990 ! cropland or grassland, apply irrigation during the growing seaspon when
993 cropsm=1.1*wilt - qmin
994 cropfr = min(1.,lufrac(crop) + 0.4*lufrac(natural)) ! assume that 40% of natural is cropland
995 newsm = cropsm*cropfr + (1.-cropfr)*soilm1d(k)
996 if(soilm1d(k) < newsm) then
997 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
998 print * ,'soil moisture is below wilting in cropland category at time step',ktau &
999 ,'i,j,lufrac(crop),k,soilm1d(k),wilt,cropsm', &
1000 i,j,lufrac(crop),k,soilm1d(k),wilt,cropsm
1003 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
1004 print * ,'added soil water to grassland category, i,j,k,soilm1d(k)',i,j,k,soilm1d(k)
1008 endif ! crop or natural
1011 ! fill in field_sf to pass perturbed field of hydraulic cond. up to model driver and output
1013 if (spp_lsm==1) then
1015 field_sf(i,k,j)=field_sf_loc(i,k,j)
1021 !--- available and maximum soil moisture content in the soil
1028 smavail(i,j)=smavail(i,j)+(qmin+soilm1d(k))* &
1029 (zshalf(k+1)-zshalf(k))
1030 smmax (i,j) =smmax (i,j)+(qmin+dqm)* &
1031 (zshalf(k+1)-zshalf(k))
1034 smavail(i,j)=smavail(i,j)+(qmin+soilm1d(nzs))* &
1035 (zsmain(nzs)-zshalf(nzs))
1036 smmax (i,j) =smmax (i,j)+(qmin+dqm)* &
1037 (zsmain(nzs)-zshalf(nzs))
1039 !--- convert the water unit into mm
1040 sfcrunoff(i,j) = sfcrunoff(i,j)+runoff1(i,j)*dt*1000.0
1041 udrunoff (i,j) = udrunoff(i,j)+runoff2(i,j)*dt*1000.0
1042 acrunoff(i,j) = acrunoff(i,j)+runoff1(i,j)*dt*1000.0
1043 smavail (i,j) = smavail(i,j) * 1000.
1044 smmax (i,j) = smmax(i,j) * 1000.
1045 smtotold (i,j) = smtotold(i,j) * 1000.
1049 soilmois(i,k,j) = soilm1d(k) + qmin
1050 sh2o (i,k,j) = min(soiliqw(k) + qmin,soilmois(i,k,j))
1051 tso(i,k,j) = tso1d(k)
1054 tso(i,nzs,j) = tbot(i,j)
1057 smfr3d(i,k,j) = smfrkeep(k)
1058 keepfr3dflag(i,k,j) = keepfr (k)
1061 z0 (i,j) = znt (i,j)
1063 patmb=p8w(i,1,j)*1.e-2
1064 q2sat=qsn(tabs,tbq)/patmb
1065 qsfc(i,j) = qvg(i,j)/(1.+qvg(i,j))
1066 ! for myj surface and pbl scheme
1068 ! myjsfc expects qsfc as actual specific humidity at the surface
1069 if((qvatm.ge.q2sat*0.95).and.qvatm.lt.qvg(i,j))then
1075 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
1076 if(chklowq(i,j).eq.0.) then
1077 print *,'i,j,chklowq', &
1082 if(snow(i,j)==0.) emissl(i,j) = lemitbl(ivgtyp(i,j))
1083 emiss (i,j) = emissl(i,j)
1084 ! snow is in [mm], snwe is in [m]; canwat is in mm, canwatr is in m
1085 snow (i,j) = snwe*1000.
1087 canwat (i,j) = canwatr*1000.
1089 infiltr(i,j) = infiltrp
1091 mavail (i,j) = lmavail(i,j)
1092 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
1093 print *,' land, i=,j=, qfx, hfx after sfctmp', i,j,lh(i,j),hfx(i,j)
1095 sfcevp (i,j) = sfcevp (i,j) + qfx (i,j) * dt
1096 grdflx (i,j) = -1. * sflx(i,j)
1098 ! if(smf(i,j) .ne.0.) then
1099 !tgs - smf.ne.0. when there is phase change in the top soil layer
1100 ! the heat of soil water freezing/thawing is not computed explicitly
1101 ! and is responsible for the residual in the energy budget.
1102 ! print *,'budget',budget(i,j),i,j,smf(i,j)
1105 !--- snowc snow cover flag
1106 if(snowfrac > 0. .and. xice(i,j).ge.xice_threshold ) then
1107 snowfrac = snowfrac*xice(i,j)
1112 !--- rhosnf - density of snowfall
1113 rhosnf(i,j)=rhosnfall
1115 ! accumulated moisture flux [kg/m^2]
1116 sfcevp (i,j) = sfcevp (i,j) + qfx (i,j) * dt
1118 ! if(smf(i,j) .ne.0.) then
1119 !tgs - smf.ne.0. when there is phase change in the top soil layer
1120 ! the heat of freezing/thawing of soil water is not computed explicitly
1121 ! and is responsible for the residual in the energy budget.
1123 ! budget(i,j)=budget(i,j)-smf(i,j)
1128 ac=max(0.,canwat(i,j)-canwatold(i,j))
1129 as=max(0.,snwe-snowold(i,j))
1130 wb =rainbl(i,j)+smelt(i,j)*dt*1.e3 & ! source
1132 -runoff1(i,j)*dt*1.e3-runoff2(i,j)*dt*1.e3 &
1133 -ac-as - (smavail(i,j)-smtotold(i,j))
1135 waterbudget(i,j)=rainbl(i,j)+smelt(i,j)*dt*1.e3 & ! source
1137 -runoff1(i,j)*dt*1.e3-runoff2(i,j)*dt*1.e3 &
1138 -ac-as - (smavail(i,j)-smtotold(i,j))
1141 acwaterbudget(i,j)=acwaterbudget(i,j)+waterbudget(i,j)
1143 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
1144 print *,'smf=',smf(i,j),i,j
1145 print *,'budget',budget(i,j),i,j
1146 print *,'runoff2= ', i,j,runoff2(i,j)
1147 print *,'water budget ', i,j,waterbudget(i,j)
1148 print *,'rainbl,qfx*dt,runoff1,smelt*dt*1.e3,smchange', &
1149 i,j,rainbl(i,j),qfx(i,j)*dt,runoff1(i,j)*dt*1.e3, &
1150 smelt(i,j)*dt*1.e3, &
1151 (smavail(i,j)-smtotold(i,j))
1153 print *,'snow,snowold',i,j,snwe,snowold(i,j)
1154 print *,'snow-snowold',i,j,max(0.,snwe-snowold(i,j))
1155 print *,'canwatold, canwat ',i,j,canwatold(i,j),canwat(i,j)
1156 print *,'canwat(i,j)-canwatold(i,j)',max(0.,canwat(i,j)-canwatold(i,j))
1160 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
1161 print *,'land, i,j,tso1d,soilm1d,soilt - end of time step', &
1162 i,j,tso1d,soilm1d,soilt(i,j)
1163 print *,'land, qfx, hfx after sfctmp', i,j,lh(i,j),hfx(i,j)
1166 !--- end of a land or sea ice point
1168 2999 continue ! lakes
1174 !-----------------------------------------------------------------
1175 end subroutine lsmruc
1176 !-----------------------------------------------------------------
1180 subroutine sfctmp (spp_lsm,rstochcol,fieldcol_sf, &
1181 delt,ktau,conflx,i,j, &
1182 !--- input variables
1183 nzs,nddzs,nroot,meltfactor, &
1184 iland,isoil,xland,ivgtyp,isltyp,prcpms, &
1185 newsnms,snwe,snhei,snowfrac, &
1186 rhosn,rhonewsn,rhosnfall, &
1187 snowrat,grauprat,icerat,curat, &
1188 patm,tabs,qvatm,qcatm,rho, &
1189 glw,gsw,emiss,qkms,tkms,pc, &
1190 mavail,cst,vegfra,alb,znt, &
1191 alb_snow,alb_snow_free,lai, &
1193 !--- soil fixed fields
1194 qwrtz,rhocs,dqm,qmin,ref,wilt,psis,bclh,ksat, &
1195 sat,cn,zsmain,zshalf,dtdzs,dtdzs2,tbq, &
1197 cp,rovcp,g0,lv,stbolt,cw,c1sn,c2sn, &
1199 !--- output variables
1200 snweprint,snheiprint,rsm, &
1201 soilm1d,ts1d,smfrkeep,keepfr,soilt,soilt1, &
1202 tsnav,dew,qvg,qsg,qcg, &
1203 smelt,snoh,snflx,snom,snowfallac,acsnow, &
1204 edir1,ec1,ett1,eeta,qfx,hfx,s,sublim, &
1205 evapl,prcpl,fltot,runoff1,runoff2,soilice, &
1206 soiliqw,infiltr,smf)
1207 !-----------------------------------------------------------------
1209 !-----------------------------------------------------------------
1211 !--- input variables
1213 integer, intent(in ) :: isice,i,j,nroot,ktau,nzs , &
1214 nddzs !nddzs=2*(nzs-2)
1216 real, intent(in ) :: delt,conflx,meltfactor
1217 real, intent(in ) :: c1sn,c2sn
1218 logical, intent(in ) :: myj
1219 !--- 3-d atmospheric variables
1221 intent(in ) :: patm, &
1226 intent(in ) :: glw, &
1238 integer, intent(in ) :: ivgtyp, isltyp
1241 intent(inout) :: emiss, &
1248 !--- soil properties
1261 real, intent(in ) :: cn, &
1272 real, dimension(1:nzs), intent(in) :: zsmain, &
1276 real, dimension(1:nzs), intent(in) :: rstochcol
1277 real, dimension(1:nzs), intent(inout) :: fieldcol_sf
1280 real, dimension(1:nddzs), intent(in) :: dtdzs
1282 real, dimension(1:5001), intent(in) :: tbq
1285 !--- input/output variables
1286 !-------- 3-d soil moisture and temperature
1287 real, dimension( 1:nzs ) , &
1288 intent(inout) :: ts1d, &
1291 real, dimension( 1:nzs ) , &
1292 intent(inout) :: keepfr
1294 real, dimension(1:nzs), intent(inout) :: soilice, &
1298 integer, intent(inout) :: iland,isoil
1301 !-------- 2-d variables
1303 intent(inout) :: dew, &
1342 real, dimension(1:nzs) :: &
1353 !-------- 1-d variables
1379 real, intent(inout) :: rsm, &
1382 integer, intent(in) :: spp_lsm
1383 !--- local variables
1387 real :: bsn, xsn , &
1388 rainf, snth, newsn, prcpms, newsnms , &
1390 real :: snhei_crit, snhei_crit_newsn, keep_snow_albedo, snowfracnewsn
1391 real :: newsnowratio, dd1, snowfrac2, m
1393 real :: rhonewgr,rhonewice
1395 real :: rnet,gswnew,gswin,emissn,zntsn,emiss_snowfree
1396 real :: vegfrac, snow_mosaic, snfr, vgfr
1397 real :: cice, albice, albsn, drip, dripsn, dripliq
1398 real :: interw, intersn, infwater, intwratio
1400 !-----------------------------------------------------------------
1401 integer, parameter :: ilsnow=99
1403 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
1404 print *,' in sfctmp',i,j,nzs,nddzs,nroot, &
1405 snwe,rhosn,snom,smelt,ts1d
1408 !-- snow fraction options
1409 !-- option 1: original formulation using critical snow depth to compute
1411 !-- option 2: the tanh formulation from niu,g.-y.,and yang,z.-l.
1412 !2007,jgr,doi:10.1029/2007jd008674.
1413 !-- option 3: the tanh formulation from niu,g.-y.,and yang,z.-l.
1414 !2007,jgr,doi:10.1029/2007jd008674.
1415 ! with vegetation dependent parameters from noah mp (personal
1416 ! communication with mike barlage)
1417 !-- snhei_crit is a threshold for fractional snow in isncovr_opt=1
1418 snhei_crit=0.01601*rhowater/rhosn
1419 snhei_crit_newsn=0.0005*rhowater/rhosn
1421 zntsn = z0tbl(isice)
1429 if(snhei == 0.) snowfrac=0.
1435 ! jul 2016 - Avissar and Pielke (1989)
1436 ! this formulation depending on lai defines relative contribution of the vegetation to
1437 ! the total heat fluxes between surface and atmosphere.
1438 ! with vegfra=100% and lai=3, vegfrac=0.86 meaning that vegetation contributes
1439 ! only 86% of the total surface fluxes.
1440 ! vgfr=0.01*vegfra ! % --> fraction
1441 ! vegfrac=2.*lai*vgfr/(1.+2.*lai*vgfr)
1451 !---initialize local arrays for sea ice
1462 albice=alb_snow_free
1465 emiss_snowfree = lemitbl(ivgtyp)
1467 !--- sea ice properties
1468 !--- n.n Zubov "arctic ice"
1469 !--- no salinity dependence because we consider the ice pack
1470 !--- to be old and to have low salinity (0.0002)
1471 if(seaice.ge.0.5) then
1473 tice(k) = ts1d(k) - 273.15
1474 rhosice(k) = 917.6/(1-0.000165*tice(k))
1475 cice = 2115.85 +7.7948*tice(k)
1476 capice(k) = cice*rhosice(k)
1477 thdifice(k) = 2.260872/capice(k)
1479 !-- sea ice alb dependence on ice temperature. when ice temperature is
1480 !-- below critical value of -10c - no change to albedo.
1481 !-- if temperature is higher that -10c then albedo is decreasing.
1482 !-- the minimum albedo at t=0c for ice is 0.1 less.
1483 albice = min(alb_snow_free,max(alb_snow_free - 0.05, &
1484 alb_snow_free - 0.1*(tice(1)+10.)/10. ))
1487 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
1488 print *,'alb_snow_free',alb_snow_free
1489 print *,'gsw,gswnew,glw,soilt,emiss,alb,albice,snwe',&
1490 gsw,gswnew,glw,soilt,emiss,alb,albice,snwe
1493 if(snhei.gt.0.0081*1.e3/rhosn) then
1494 !*** update snow density for current temperature (koren et al. 1999)
1495 bsn=delt/3600.*c1sn*exp(0.08*min(0.,tsnav)-c2sn*rhosn*1.e-3)
1496 if(bsn*snwe*100..lt.1.e-4) goto 777
1497 xsn=rhosn*(exp(bsn*snwe*100.)-1.)/(bsn*snwe*100.)
1498 rhosn=min(max(58.8,xsn),500.) ! 13mar18 - switch from 76.9 to 58.8
1503 !-- snow_mosaic from the previous time step
1504 if(snowfrac < 0.75) snow_mosaic = 1.
1508 if(newsn.gt.0.) then
1509 ! if(newsn.ge.1.e-8) then
1511 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
1512 print *, 'there is new snow, newsn', newsn
1515 newsnowratio = min(1.,newsn/(snwe+newsn))
1517 !--- 27 feb 2014 - empirical formulations from john m. brown
1518 ! rhonewsn=min(250.,rhowater/max(4.179,(13.*tanh((274.15-tabs)*0.3333))))
1519 !--- 13 mar 2018 - formulation from trevor alcott
1520 rhonewsn=min(125.,1000.0/max(8.,(17.*tanh((276.65-tabs)*0.15))))
1521 rhonewgr=min(500.,rhowater/max(2.,(3.5*tanh((274.15-tabs)*0.3333))))
1524 !--- compute density of "snowfall" from weighted contribution
1525 ! of snow, graupel and ice fractions
1527 rhosnfall = min(500.,max(58.8,(rhonewsn*snowrat + & ! 13mar18-switch from 76.9 to 58.8
1528 rhonewgr*grauprat + rhonewice*icerat + rhonewgr*curat)))
1530 ! from now on rhonewsn is the density of falling frozen precipitation
1533 !*** define average snow density of the snow pack considering
1534 !*** the amount of fresh snow (eq. 9 in koren et al.(1999)
1535 !*** without snow melt )
1536 xsn=(rhosn*snwe+rhonewsn*newsn)/ &
1538 rhosn=min(max(58.8,xsn),500.) ! 13mar18 - switch from 76.9 to 58.8
1540 endif ! end newsn > 0.
1542 if(prcpms.ne.0.) then
1544 ! prcpms is liquid precipitation rate
1545 ! rainf is a flag used for calculation of rain water
1546 ! heat content contribution into heat budget equation. rain's temperature
1547 ! is set equal to air temperature at the first atmospheric
1555 if(vegfrac > 0.01) then
1556 ! compute intercepted precipitation - eq. 1 Lawrence et al.,
1557 ! j. of hydrometeorology, 2006, clm.
1558 interw=0.25*delt*prcpms*(1.-exp(-0.5*lai))*vegfrac
1559 intersn=0.25*newsn*(1.-exp(-0.5*lai))*vegfrac
1560 infwater=prcpms - interw/delt
1561 if((interw+intersn) > 0.) then
1562 intwratio=interw/(interw+intersn)
1565 ! update water/snow intercepted by the canopy
1566 dd1=cst + interw + intersn
1578 endif ! vegfrac > 0.01
1580 if(newsn.gt.0.) then
1581 !update snow on the ground
1582 snwe=max(0.,snwe+newsn-intersn)
1583 ! add drip to snow on the ground
1585 if (snow_mosaic==1.) then
1586 dripliq=drip*intwratio
1587 dripsn = drip - dripliq
1589 infwater=infwater+dripliq
1596 snhei=snwe*rhowater/rhosn
1597 newsn=newsn*rhowater/rhonewsn
1600 if(snhei.gt.0.0) then
1601 !-- snow on the ground
1602 !--- land-use category should be changed to snow/ice for grid points with snow>0
1604 !24nov15 - based on field exp on pleasant view soccer fields
1605 ! if(meltfactor > 1.5) then ! all veg. types, except forests
1606 ! snhei_crit=0.01601*1.e3/rhosn
1607 ! petzold - 1 cm of fresh snow overwrites effects from old snow.
1608 ! need to test snhei_crit_newsn=0.01
1609 ! snhei_crit_newsn=0.01
1611 ! snhei_crit=0.02*1.e3/rhosn
1612 ! snhei_crit_newsn=0.001*1.e3/rhosn
1615 if(isncovr_opt == 1) then
1616 snowfrac=min(1.,snhei/(2.*snhei_crit))
1617 elseif(isncovr_opt == 2) then
1618 snowfrac=min(1.,snhei/(2.*snhei_crit))
1619 !if(ivgtyp == glacier .or. ivgtyp == bare) then
1620 !-- sparsely vegetated or land ice
1621 ! snowfrac2 = tanh( snhei/(2.5 * 0.2 *(rhosn/rhonewsn)**1.))
1623 !-- Niu&Yang: znt=0.01 m for 1 degree (100km) resolution tests
1624 ! on 3-km scale use actual roughness, but not higher than 0.2 m.
1625 ! the factor is 20 for forests (~100/dx = 33., dx=3 km)
1626 snowfrac2 = tanh( snhei/(2.5 * min(0.2,znt) *(rhosn/rhonewsn)**1.))
1628 !-- snow fraction is average between method 1 and 2
1629 snowfrac = 0.5*(snowfrac+snowfrac2)
1632 !m = mfsno(ivgtyp) ! vegetation dependent facsnf/msnf from noah mp
1634 !-- for rrfs a factor 10. was added to 'facsnf' to get reasonal values of
1635 ! snow cover fractions on the 3-km scale.
1636 ! this factor is scale dependent.
1637 snowfrac = tanh( snhei/(10. * sncovfac(ivgtyp)*(rhosn/rhonewsn)**m))
1640 if(newsn > 0. ) then
1641 snowfracnewsn=min(1.,snowfallac*1.e-3/snhei_crit_newsn)
1644 !24nov15 - snowfrac for urban category < 0.75
1645 if(ivgtyp == urban) snowfrac=min(0.75,snowfrac)
1646 ! if(meltfactor > 1.5) then
1647 ! if(isltyp > 9 .and. isltyp < 13) then
1648 !24nov15 clay soil types - snofrac < 0.9
1649 ! snowfrac=min(0.9,snowfrac)
1652 !24nov15 - snowfrac for forests < 0.75
1653 ! snowfrac=min(0.85,snowfrac)
1656 if(snowfrac < 0.75) snow_mosaic = 1.
1658 keep_snow_albedo = 0.
1659 if (snowfracnewsn > 0.99 .and. rhosnfall < 450.) then
1661 keep_snow_albedo = 1.
1665 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
1666 print *,'snhei_crit,snowfrac,snhei_crit_newsn,snowfracnewsn', &
1667 snhei_crit,snowfrac,snhei_crit_newsn,snowfracnewsn
1670 !-- set znt for snow from VEGPARM table (snow/ice landuse), except for
1671 !-- land-use types with higher roughness (forests, urban).
1672 if(newsn.eq.0. .and. znt.le.0.2 .and. ivgtyp.ne.isice) then
1673 if( snhei .le. 2.*znt)then
1674 znt=0.55*znt+0.45*z0tbl(iland)
1675 elseif( snhei .gt. 2.*znt .and. snhei .le. 4.*znt)then
1676 znt=0.2*znt+0.8*z0tbl(iland)
1677 elseif(snhei > 4.*znt) then
1682 if(seaice .lt. 0.5) then
1684 !-- alb dependence on snow depth
1685 ! alb_snow across canada's forested areas is very low - 0.27-0.35, this
1686 ! causes significant warm biases. limiting alb in these areas to be higher than 0.4
1687 ! hwlps with these biases..
1688 if( snow_mosaic == 1.) then
1690 if(keep_snow_albedo > 0.9 .and. albsn < 0.4) then
1691 !-- Albedo correction with fresh snow and deep snow pack
1692 !-- will reduce warm bias in western Canada
1693 !-- and US West coast, where max snow albedo is low (0.3-0.5).
1694 !print *,'ALB increase to 0.7',alb_snow,snhei,snhei_crit,albsn,i,j
1700 albsn = max(keep_snow_albedo*alb_snow, &
1701 min((alb_snow_free + &
1702 (alb_snow - alb_snow_free) * snowfrac), alb_snow))
1703 if(newsn > 0. .and. keep_snow_albedo > 0.9 .and. albsn < 0.4) then
1704 !-- Albedo correction with fresh snow and deep snow pack
1705 !-- will reduce warm bias in western Canada
1706 !-- and US West coast, where max snow albedo is low (0.3-0.5).
1707 !print *,'ALB increase to 0.7',alb_snow,snhei,snhei_crit,albsn,i,j
1709 !print *,'NO mosaic ALB increase to 0.7',alb_snow,snhei,snhei_crit,alb,i,j
1712 emiss = max(keep_snow_albedo*emissn, &
1713 min((emiss_snowfree + &
1714 (emissn - emiss_snowfree) * snowfrac), emissn))
1716 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
1717 print *,'snow on soil albsn,emiss,snow_mosaic',i,j,albsn,emiss,snow_mosaic
1719 !28mar11 if canopy is covered with snow to 95% of its capacity and snow depth is
1720 ! higher than patchy snow treshold - then snow albedo is not less than 0.55
1721 ! (inspired by the flight from fairbanks to seatle)
1722 !test if(cst.ge.0.95*sat .and. snowfrac .gt.0.99)then
1723 ! albsn=max(alb_snow,0.55)
1726 !-- alb dependence on snow temperature. when snow temperature is
1727 !-- below critical value of -10c - no change to albedo.
1728 !-- if temperature is higher that -10c then albedo is decreasing.
1729 !-- the minimum albedo at t=0c for snow on land is 15% less than
1730 !-- albedo of temperatures below -10c.
1731 if(albsn.lt.0.4 .or. keep_snow_albedo==1) then
1734 !-- change albedo when no fresh snow and snow albedo is higher than 0.5
1735 alb = min(albsn,max(albsn - 0.1*(soilt - 263.15)/ &
1736 (273.15-263.15)*albsn, albsn - 0.05))
1740 if( snow_mosaic == 1.) then
1744 albsn = max(keep_snow_albedo*alb_snow, &
1745 min((albice + (alb_snow - albice) * snowfrac), alb_snow))
1746 emiss = max(keep_snow_albedo*emissn, &
1747 min((emiss_snowfree + &
1748 (emissn - emiss_snowfree) * snowfrac), emissn))
1751 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
1752 print *,'snow on ice snow_mosaic,albsn,emiss',i,j,albsn,emiss,snow_mosaic
1754 !-- alb dependence on snow temperature. when snow temperature is
1755 !-- below critical value of -10c - no change to albedo.
1756 !-- if temperature is higher that -10c then albedo is decreasing.
1757 if(albsn.lt.alb_snow .or. keep_snow_albedo .eq.1.)then
1760 !-- change albedo when no fresh snow
1761 alb = min(albsn,max(albsn - 0.15*albsn*(soilt - 263.15)/ &
1762 (273.15-263.15), albsn - 0.1))
1767 if (snow_mosaic==1.) then
1768 !may 2014 - treat separately snow-free and snow-covered areas
1770 if(seaice .lt. 0.5) then
1772 ! portion not covered with snow
1773 ! compute absorbed gsw for snow-free portion
1775 gswnew=gswin*(1.-alb_snow_free)
1777 t3 = stbolt*soilt*soilt*soilt
1779 xinet = emiss_snowfree*(glw-upflux)
1780 rnet = gswnew + xinet
1781 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
1782 print *,'fractional snow - snowfrac=',snowfrac
1783 print *,'snowfrac<1 gswin,gswnew -',gswin,gswnew,'soilt, rnet',soilt,rnet
1786 soilm1ds(k) = soilm1d(k)
1788 smfrkeeps(k) = smfrkeep(k)
1789 keepfrs(k) = keepfr(k)
1790 soilices(k) = soilice(k)
1791 soiliqws(k) = soiliqw(k)
1804 call soil(spp_lsm,rstochcol,fieldcol_sf, &
1805 !--- input variables
1806 i,j,ilands,isoil,delt,ktau,conflx,nzs,nddzs,nroot, &
1807 prcpms,rainf,patm,qvatm,qcatm,glw,gswnew,gswin, &
1808 emiss_snowfree,rnet,qkms,tkms,pc,csts,dripliq, &
1809 infwater,rho,vegfrac,lai,myj, &
1810 !--- soil fixed fields
1811 qwrtz,rhocs,dqm,qmin,ref,wilt, &
1812 psis,bclh,ksat,sat,cn, &
1813 zsmain,zshalf,dtdzs,dtdzs2,tbq, &
1815 lv,cp,rovcp,g0,cw,stbolt,tabs, &
1817 !--- output variables for snow-free portion
1818 soilm1ds,ts1ds,smfrkeeps,keepfrs, &
1819 dews,soilts,qvgs,qsgs,qcgs,edir1s,ec1s, &
1820 ett1s,eetas,qfxs,hfxs,ss,evapls,prcpls,fltots,runoff1s, &
1821 runoff2s,mavails,soilices,soiliqws, &
1825 ! portion not covered with snow
1826 ! compute absorbed gsw for snow-free portion
1828 gswnew=gswin*(1.-albice)
1830 t3 = stbolt*soilt*soilt*soilt
1832 xinet = emiss_snowfree*(glw-upflux)
1833 rnet = gswnew + xinet
1834 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
1835 print *,'fractional snow - snowfrac=',snowfrac
1836 print *,'snowfrac<1 gswin,gswnew -',gswin,gswnew,'soilt, rnet',soilt,rnet
1850 !--- input variables
1851 i,j,iland,isoil,delt,ktau,conflx,nzs,nddzs,nroot, &
1852 prcpms,rainf,patm,qvatm,qcatm,glw,gswnew, &
1853 0.98,rnet,qkms,tkms,rho,myj, &
1854 !--- sea ice parameters
1855 tice,rhosice,capice,thdifice, &
1856 zsmain,zshalf,dtdzs,dtdzs2,tbq, &
1858 lv,cp,rovcp,cw,stbolt,tabs, &
1859 !--- output variable
1860 ts1ds,dews,soilts,qvgs,qsgs,qcgs, &
1861 eetas,qfxs,hfxs,ss,evapls,prcpls,fltots &
1878 endif ! seaice < 0.5
1880 !return gswnew to incoming solar
1881 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
1882 print *,'gswnew,alb_snow_free,alb',gswnew,alb_snow_free,alb
1885 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
1886 print *,'incoming gswnew snowfrac<1 -',gswnew
1888 endif ! snow_mosaic=1.
1890 !--- recompute absorbed solar radiation and net radiation
1891 !--- for updated value of snow albedo - alb
1892 gswnew=gswin*(1.-alb)
1893 ! print *,'snow fraction gswnew',gswnew,'alb=',alb
1895 t3 = stbolt*soilt*soilt*soilt
1897 xinet = emiss*(glw-upflux)
1898 rnet = gswnew + xinet
1899 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
1900 print *,'rnet=',rnet
1901 print *,'snow - i,j,newsn,snwe,snhei,gsw,gswnew,glw,upflux,alb',&
1902 i,j,newsn,snwe,snhei,gsw,gswnew,glw,upflux,alb
1905 if (seaice .lt. 0.5) then
1907 if(snow_mosaic==1.)then
1912 call snowsoil (spp_lsm,rstochcol,fieldcol_sf, & !--- input variables
1913 i,j,isoil,delt,ktau,conflx,nzs,nddzs,nroot, &
1914 meltfactor,rhonewsn,snhei_crit, & ! new
1915 iland,prcpms,rainf,newsn,snhei,snwe,snfr, &
1916 rhosn,patm,qvatm,qcatm, &
1917 glw,gswnew,gswin,emiss,rnet,ivgtyp, &
1918 qkms,tkms,pc,cst,dripsn,infwater, &
1919 rho,vegfrac,alb,znt,lai, &
1921 !--- soil fixed fields
1922 qwrtz,rhocs,dqm,qmin,ref,wilt,psis,bclh,ksat, &
1923 sat,cn,zsmain,zshalf,dtdzs,dtdzs2,tbq, &
1925 lv,cp,rovcp,g0,cw,stbolt,tabs, &
1927 !--- output variables
1928 ilnb,snweprint,snheiprint,rsm, &
1929 soilm1d,ts1d,smfrkeep,keepfr, &
1930 dew,soilt,soilt1,tsnav,qvg,qsg,qcg, &
1931 smelt,snoh,snflx,snom,edir1,ec1,ett1,eeta, &
1932 qfx,hfx,s,sublim,prcpl,fltot,runoff1,runoff2, &
1933 mavail,soilice,soiliqw,infiltr )
1936 if(snow_mosaic==1.)then
1943 i,j,isoil,delt,ktau,conflx,nzs,nddzs, &
1944 meltfactor,rhonewsn,snhei_crit, & ! new
1945 iland,prcpms,rainf,newsn,snhei,snwe,snfr, &
1946 rhosn,patm,qvatm,qcatm, &
1947 glw,gswnew,emiss,rnet, &
1948 qkms,tkms,rho,myj, &
1949 !--- sea ice parameters
1951 tice,rhosice,capice,thdifice, &
1952 zsmain,zshalf,dtdzs,dtdzs2,tbq, &
1954 lv,cp,rovcp,cw,stbolt,tabs, &
1955 !--- output variables
1956 ilnb,snweprint,snheiprint,rsm,ts1d, &
1957 dew,soilt,soilt1,tsnav,qvg,qsg,qcg, &
1958 smelt,snoh,snflx,snom,eeta, &
1959 qfx,hfx,s,sublim,prcpl,fltot &
1979 if (snow_mosaic==1.) then
1980 ! may 2014 - now combine snow covered and snow-free land fluxes, soil temp, moist,
1982 if(seaice .lt. 0.5) then
1984 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
1985 print *,'soilt snow on land', ktau, i,j,soilt
1986 print *,'soilt on snow-free land', i,j,soilts
1987 print *,'ts1d,ts1ds',i,j,ts1d,ts1ds
1988 print *,' snow flux',i,j, snflx
1989 print *,' ground flux on snow-covered land',i,j, s
1990 print *,' ground flux on snow-free land', i,j,ss
1991 print *,' csts, cst', i,j,csts,cst
1994 soilm1d(k) = soilm1ds(k)*(1.-snowfrac) + soilm1d(k)*snowfrac
1995 ts1d(k) = ts1ds(k)*(1.-snowfrac) + ts1d(k)*snowfrac
1996 smfrkeep(k) = smfrkeeps(k)*(1.-snowfrac) + smfrkeep(k)*snowfrac
1997 if(snowfrac > 0.5) then
1998 keepfr(k) = keepfr(k)
2000 keepfr(k) = keepfrs(k)
2002 soilice(k) = soilices(k)*(1.-snowfrac) + soilice(k)*snowfrac
2003 soiliqw(k) = soiliqws(k)*(1.-snowfrac) + soiliqw(k)*snowfrac
2005 dew = dews*(1.-snowfrac) + dew*snowfrac
2006 soilt = soilts*(1.-snowfrac) + soilt*snowfrac
2007 qvg = qvgs*(1.-snowfrac) + qvg*snowfrac
2008 qsg = qsgs*(1.-snowfrac) + qsg*snowfrac
2009 qcg = qcgs*(1.-snowfrac) + qcg*snowfrac
2010 edir1 = edir1s*(1.-snowfrac) + edir1*snowfrac
2011 ec1 = ec1s*(1.-snowfrac) + ec1*snowfrac
2012 cst = csts*(1.-snowfrac) + cst*snowfrac
2013 ett1 = ett1s*(1.-snowfrac) + ett1*snowfrac
2014 eeta = eetas*(1.-snowfrac) + eeta*snowfrac
2015 qfx = qfxs*(1.-snowfrac) + qfx*snowfrac
2016 hfx = hfxs*(1.-snowfrac) + hfx*snowfrac
2017 s = ss*(1.-snowfrac) + s*snowfrac
2018 evapl = evapls*(1.-snowfrac)
2019 sublim = sublim*snowfrac
2020 prcpl = prcpls*(1.-snowfrac) + prcpl*snowfrac
2021 fltot = fltots*(1.-snowfrac) + fltot*snowfrac
2023 alb = max(keep_snow_albedo*alb, &
2024 min((alb_snow_free + (alb - alb_snow_free) * snowfrac), alb))
2026 emiss = max(keep_snow_albedo*emissn, &
2027 min((emiss_snowfree + &
2028 (emissn - emiss_snowfree) * snowfrac), emissn))
2030 runoff1 = runoff1s*(1.-snowfrac) + runoff1*snowfrac
2031 runoff2 = runoff2s*(1.-snowfrac) + runoff2*snowfrac
2032 mavail = mavails*(1.-snowfrac) + 1.*snowfrac
2033 infiltr = infiltrs*(1.-snowfrac) + infiltr*snowfrac
2035 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
2036 print *,' ground flux combined', i,j, s
2037 print *,'soilt combined on land', soilt
2038 print *,'ts combined on land', ts1d
2042 ! now combine fluxes for snow-free sea ice and snow-covered area
2043 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
2044 print *,'soilt snow on ice', soilt
2047 ts1d(k) = ts1ds(k)*(1.-snowfrac) + ts1d(k)*snowfrac
2049 dew = dews*(1.-snowfrac) + dew*snowfrac
2050 soilt = soilts*(1.-snowfrac) + soilt*snowfrac
2051 qvg = qvgs*(1.-snowfrac) + qvg*snowfrac
2052 qsg = qsgs*(1.-snowfrac) + qsg*snowfrac
2053 qcg = qcgs*(1.-snowfrac) + qcg*snowfrac
2054 eeta = eetas*(1.-snowfrac) + eeta*snowfrac
2055 qfx = qfxs*(1.-snowfrac) + qfx*snowfrac
2056 hfx = hfxs*(1.-snowfrac) + hfx*snowfrac
2057 s = ss*(1.-snowfrac) + s*snowfrac
2059 prcpl = prcpls*(1.-snowfrac) + prcpl*snowfrac
2060 fltot = fltots*(1.-snowfrac) + fltot*snowfrac
2062 alb = max(keep_snow_albedo*alb, &
2063 min((albice + (alb - alb_snow_free) * snowfrac), alb))
2065 emiss = max(keep_snow_albedo*emissn, &
2066 min((emiss_snowfree + &
2067 (emissn - emiss_snowfree) * snowfrac), emissn))
2069 runoff1 = runoff1s*(1.-snowfrac) + runoff1*snowfrac
2070 runoff2 = runoff2s*(1.-snowfrac) + runoff2*snowfrac
2071 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
2072 print *,'soilt combined on ice', soilt
2075 endif ! snow_mosaic = 1.
2077 if(snhei.eq.0.) then
2078 !-- all snow is melted
2082 !-- snow on the ground
2083 if(isncovr_opt == 1) then
2084 snowfrac=min(1.,snhei/(2.*snhei_crit))
2085 elseif(isncovr_opt == 2) then
2086 snowfrac=min(1.,snhei/(2.*snhei_crit))
2087 !if(ivgtyp == glacier .or. ivgtyp == bare) then
2088 !-- sparsely vegetated or land ice
2089 ! snowfrac2 = tanh( snhei/(2.5 * 0.2 *(rhosn/rhonewsn)**1.))
2091 !-- niu&yang: znt=0.01 m for 1 degree (100km) resolution tests
2092 ! on 3-km scale use actual roughness, but not higher than 0.2 m.
2093 ! the factor is 20 for forests (~100/dx = 33.)
2094 snowfrac2 = tanh( snhei/(2.5 * min(0.2,znt) *(rhosn/rhonewsn)**1.))
2096 !-- snow fraction is average between method 1 and 2
2097 snowfrac = 0.5*(snowfrac+snowfrac2)
2100 !m = mfsno(ivgtyp) ! vegetation dependent facsnf/msnf from noah mp
2102 !-- for rrfs a factor 10. was added to 'facsnf' to get reasonal values
2104 ! snow cover fractions on the 3-km scale.
2105 ! this factor is scale dependent.
2106 snowfrac = tanh( snhei/(10. * sncovfac(ivgtyp)*(rhosn/rhonewsn)**m))
2111 if(ivgtyp == urban) snowfrac=min(0.75,snowfrac)
2113 ! run-total accumulated snow based on snowfall and snowmelt in [m]
2115 snowfallac = snowfallac + newsn * 1.e3 ! accumulated snow depth [mm], using variable snow den
2116 !snowfallac = snowfallac + max(0.,(newsn - rhowater/rhonewsn*smelt*delt*newsnowratio))
2125 t3 = stbolt*soilt*soilt*soilt
2127 xinet = emiss*(glw-upflux)
2128 rnet = gswnew + xinet
2129 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
2130 print *,'no snow on the ground gswnew -',gswnew,'rnet=',rnet
2133 if(seaice .lt. 0.5) then
2135 call soil(spp_lsm,rstochcol,fieldcol_sf, &
2136 !--- input variables
2137 i,j,iland,isoil,delt,ktau,conflx,nzs,nddzs,nroot, &
2138 prcpms,rainf,patm,qvatm,qcatm,glw,gswnew,gswin, &
2139 emiss,rnet,qkms,tkms,pc,cst,drip,infwater, &
2140 rho,vegfrac,lai,myj, &
2141 !--- soil fixed fields
2142 qwrtz,rhocs,dqm,qmin,ref,wilt, &
2143 psis,bclh,ksat,sat,cn, &
2144 zsmain,zshalf,dtdzs,dtdzs2,tbq, &
2146 lv,cp,rovcp,g0,cw,stbolt,tabs, &
2148 !--- output variables
2149 soilm1d,ts1d,smfrkeep,keepfr, &
2150 dew,soilt,qvg,qsg,qcg,edir1,ec1, &
2151 ett1,eeta,qfx,hfx,s,evapl,prcpl,fltot,runoff1, &
2152 runoff2,mavail,soilice,soiliqw, &
2156 ! if current ice albedo is not the same as from the previous time step, then
2157 ! update gsw, alb and rnet for surface energy budget
2158 if(alb.ne.albice) gswnew=gsw/(1.-alb)*(1.-albice)
2160 rnet = gswnew + xinet
2163 !--- input variables
2164 i,j,iland,isoil,delt,ktau,conflx,nzs,nddzs,nroot, &
2165 prcpms,rainf,patm,qvatm,qcatm,glw,gswnew, &
2166 emiss,rnet,qkms,tkms,rho,myj, &
2167 !--- sea ice parameters
2168 tice,rhosice,capice,thdifice, &
2169 zsmain,zshalf,dtdzs,dtdzs2,tbq, &
2171 lv,cp,rovcp,cw,stbolt,tabs, &
2172 !--- output variables
2173 ts1d,dew,soilt,qvg,qsg,qcg, &
2174 eeta,qfx,hfx,s,evapl,prcpl,fltot &
2197 !---------------------------------------------------------------
2198 end subroutine sfctmp
2199 !---------------------------------------------------------------
2203 !****************************************************************
2204 real, dimension(1:5001), intent(in ) :: t
2205 real, intent(in ) :: tn
2210 r=(tn-173.15)/.05+1.
2215 10 if(i.le.5000) goto 20
2220 qsn=(t(i+1)-r1)*r2 + r1
2221 ! print *,' in qsn, i,r,r1,r2,t(i+1),tn, qsn', i,r,r1,r2,t(i+1),tn,qsn
2224 !-----------------------------------------------------------------------
2226 !------------------------------------------------------------------------
2229 subroutine soil (spp_lsm,rstochcol, fieldcol_sf, &
2230 !--- input variables
2231 i,j,iland,isoil,delt,ktau,conflx,nzs,nddzs,nroot,&
2232 prcpms,rainf,patm,qvatm,qcatm, &
2233 glw,gsw,gswin,emiss,rnet, &
2234 qkms,tkms,pc,cst,drip,infwater,rho,vegfrac,lai, &
2236 !--- soil fixed fields
2237 qwrtz,rhocs,dqm,qmin,ref,wilt,psis,bclh,ksat, &
2238 sat,cn,zsmain,zshalf,dtdzs,dtdzs2,tbq, &
2240 xlv,cp,rovcp,g0_p,cw,stbolt,tabs, &
2242 !--- output variables
2243 soilmois,tso,smfrkeep,keepfr, &
2244 dew,soilt,qvg,qsg,qcg, &
2245 edir1,ec1,ett1,eeta,qfx,hfx,s,evapl, &
2246 prcpl,fltot,runoff1,runoff2,mavail,soilice, &
2247 soiliqw,infiltrp,smf)
2249 !*************************************************************
2250 ! energy and moisture budget for vegetated surfaces
2251 ! without snow, heat diffusion and richards eqns. in
2254 ! delt - time step (s)
2255 ! ktau - numver of time step
2256 ! conflx - depth of constant flux layer (m)
2257 ! j,i - the location of grid point
2258 ! ime, jme, kme, nzs - dimensions of the domain
2259 ! nroot - number of levels within the root zone
2260 ! prcpms - precipitation rate in m/s
2261 ! patm - pressure [bar]
2262 ! qvatm,qcatm - cloud and water vapor mixing ratio (kg/kg)
2263 ! at the first atm. level
2264 ! glw, gsw - incoming longwave and absorbed shortwave
2265 ! radiation at the surface (w/m^2)
2266 ! emiss,rnet - emissivity of the ground surface (0-1) and net
2267 ! radiation at the surface (w/m^2)
2268 ! qkms - exchange coefficient for water vapor in the
2269 ! surface layer (m/s)
2270 ! tkms - exchange coefficient for heat in the surface
2272 ! pc - plant coefficient (resistance) (0-1)
2273 ! rho - density of atmosphere near sueface (kg/m^3)
2274 ! vegfrac - greeness fraction
2275 ! rhocs - volumetric heat capacity of dry soil
2276 ! dqm, qmin - porosity minus residual soil moisture qmin (m^3/m^3)
2277 ! ref, wilt - field capacity soil moisture and the
2278 ! wilting point (m^3/m^3)
2279 ! psis - matrix potential at saturation (m)
2280 ! bclh - exponent for clapp-hornberger parameterization
2281 ! ksat - saturated hydraulic conductivity (m/s)
2282 ! sat - maximum value of water intercepted by canopy (m)
2283 ! cn - exponent for calculation of canopy water
2284 ! zsmain - main levels in soil (m)
2285 ! zshalf - middle of the soil layers (m)
2286 ! dtdzs,dtdzs2 - dt/(2.*dzshalf*dzmain) and dt/dzshalf in soil
2287 ! tbq - table to define saturated mixing ration
2288 ! of water vapor for given temperature and pressure
2289 ! soilmois,tso - soil moisture (m^3/m^3) and temperature (k)
2290 ! dew - dew in kg/m^2s
2291 ! soilt - skin temperature (k)
2292 ! qsg,qvg,qcg - saturated mixing ratio, mixing ratio of
2293 ! water vapor and cloud at the ground
2294 ! surface, respectively (kg/kg)
2295 ! edir1, ec1, ett1, eeta - direct evaporation, evaporation of
2296 ! canopy water, transpiration in kg m-2 s-1 and total
2297 ! evaporation in m s-1.
2298 ! qfx, hfx - latent and sensible heat fluxes (w/m^2)
2299 ! s - soil heat flux in the top layer (w/m^2)
2300 ! runoff - surface runoff (m/s)
2301 ! runoff2 - underground runoff (m)
2302 ! mavail - moisture availability in the top soil layer (0-1)
2303 ! infiltrp - infiltration flux from the top of soil domain (m/s)
2305 !*****************************************************************
2307 !-----------------------------------------------------------------
2309 !--- input variables
2311 integer, intent(in ) :: nroot,ktau,nzs , &
2312 nddzs !nddzs=2*(nzs-2)
2313 integer, intent(in ) :: i,j,iland,isoil
2314 real, intent(in ) :: delt,conflx
2315 logical, intent(in ) :: myj
2316 !--- 3-d atmospheric variables
2318 intent(in ) :: patm, &
2323 intent(in ) :: glw, &
2335 !--- soil properties
2337 intent(in ) :: rhocs, &
2347 real, intent(in ) :: cn, &
2356 real, dimension(1:nzs), intent(in) :: zsmain, &
2360 real, dimension(1:nddzs), intent(in) :: dtdzs
2362 real, dimension(1:5001), intent(in) :: tbq
2365 !--- input/output variables
2366 !-------- 3-d soil moisture and temperature
2367 real, dimension( 1:nzs ) , &
2368 intent(inout) :: tso, &
2372 real, dimension(1:nzs), intent(in) :: rstochcol
2373 real, dimension(1:nzs), intent(inout) :: fieldcol_sf
2376 real, dimension( 1:nzs ) , &
2377 intent(inout) :: keepfr
2379 !-------- 2-d variables
2381 intent(inout) :: dew, &
2403 !-------- 1-d variables
2404 integer , intent(in) :: spp_lsm
2405 real, dimension(1:nzs), intent(out) :: soilice, &
2408 !--- local variables
2410 real :: infiltrp, transum , &
2412 tabs, t3, upflux, xinet
2413 real :: cp,rovcp,g0,lv,stbolt,xlmelt,dzstop , &
2414 can,epot,fac,fltot,ft,fq,hft , &
2415 q1,ras,rhoice,sph , &
2416 trans,zn,ci,cvw,tln,tavln,pi , &
2417 dd1,cmc2ms,drycan,wetcan , &
2419 real, dimension(1:nzs) :: transp,cap,diffu,hydro , &
2420 thdif,tranf,tav,soilmoism , &
2421 soilicem,soiliqwm,detal , &
2422 fwsat,lwsat,told,smold
2424 real :: soiltold,smf
2425 real :: soilres, alfa, fex, fex_fc, fc, psit
2427 integer :: nzs1,nzs2,k
2429 !-----------------------------------------------------------------
2431 !-- define constants
2445 !--- initializing local arrays
2468 dzstop=1./(zsmain(2)-zsmain(1))
2472 !--- computation of volumetric content of ice in soil
2476 tln=log(tso(k)/273.15)
2478 soiliqw(k)=(dqm+qmin)*(xlmelt* &
2479 (tso(k)-273.15)/tso(k)/9.81/psis) &
2481 soiliqw(k)=max(0.,soiliqw(k))
2482 soiliqw(k)=min(soiliqw(k),soilmois(k))
2483 soilice(k)=(soilmois(k)-soiliqw(k))/riw
2485 !---- melting and freezing is balanced, soil ice cannot increase
2486 if(keepfr(k).eq.1.) then
2487 soilice(k)=min(soilice(k),smfrkeep(k))
2488 soiliqw(k)=max(0.,soilmois(k)-soilice(k)*riw)
2493 soiliqw(k)=soilmois(k)
2499 !- middle of soil layers
2500 tav(k)=0.5*(tso(k)+tso(k+1))
2501 soilmoism(k)=0.5*(soilmois(k)+soilmois(k+1))
2502 tavln=log(tav(k)/273.15)
2504 if(tavln.lt.0.) then
2505 soiliqwm(k)=(dqm+qmin)*(xlmelt* &
2506 (tav(k)-273.15)/tav(k)/9.81/psis) &
2508 fwsat(k)=dqm-soiliqwm(k)
2509 lwsat(k)=soiliqwm(k)+qmin
2510 soiliqwm(k)=max(0.,soiliqwm(k))
2511 soiliqwm(k)=min(soiliqwm(k), soilmoism(k))
2512 soilicem(k)=(soilmoism(k)-soiliqwm(k))/riw
2513 !---- melting and freezing is balanced, soil ice cannot increase
2514 if(keepfr(k).eq.1.) then
2515 soilicem(k)=min(soilicem(k), &
2516 0.5*(smfrkeep(k)+smfrkeep(k+1)))
2517 soiliqwm(k)=max(0.,soilmoism(k)-soilicem(k)*riw)
2518 fwsat(k)=dqm-soiliqwm(k)
2519 lwsat(k)=soiliqwm(k)+qmin
2524 soiliqwm(k)=soilmoism(k)
2532 if(soilice(k).gt.0.) then
2533 smfrkeep(k)=soilice(k)
2535 smfrkeep(k)=soilmois(k)/riw
2538 !******************************************************************
2539 ! soilprop computes thermal diffusivity, and diffusional and
2540 ! hydraulic condeuctivities
2541 !******************************************************************
2542 call soilprop(spp_lsm,rstochcol,fieldcol_sf, &
2543 !--- input variables
2544 nzs,fwsat,lwsat,tav,keepfr, &
2545 soilmois,soiliqw,soilice, &
2546 soilmoism,soiliqwm,soilicem, &
2547 !--- soil fixed fields
2548 qwrtz,rhocs,dqm,qmin,psis,bclh,ksat, &
2550 riw,xlmelt,cp,g0_p,cvw,ci, &
2552 !--- output variables
2553 thdif,diffu,hydro,cap)
2555 !********************************************************************
2556 !--- calculation of canopy water (Smirnova et al., 1996, eq.16) and dew
2563 q1=-qkms*ras*(qvatm - qsg)
2566 if(qvatm.ge.qsg)then
2571 ! dd1=cst+delt*(prcpms +dew*ras)
2574 ! delt*(prcpms+ras*fq*(qvatm-qsg) &
2578 ! dd1=cst+delt*prcpms
2580 ! if(dd1.lt.0.) dd1=0.
2581 ! if(vegfrac.eq.0.)then
2585 ! if (vegfrac.gt.0.) then
2587 ! if(cst.gt.sat) then
2593 !--- wetcan is the fraction of vegetated area covered by canopy
2594 !--- water, and drycan is the fraction of vegetated area where
2595 !--- transpiration may take place.
2597 wetcan=min(0.25,max(0.,(cst/sat))**cn)
2598 ! if(lai > 1.) wetcan=wetcan/lai
2601 !**************************************************************
2602 ! transf computes transpiration function
2603 !**************************************************************
2605 !--- input variables
2606 nzs,nroot,soiliqw,tabs,lai,gswin, &
2607 !--- soil fixed fields
2608 dqm,qmin,ref,wilt,zshalf,pc,iland, &
2609 !--- output variables
2613 !--- save soil temp and moisture from the beginning of time step
2616 smold(k)=soilmois(k)
2619 ! sakaguchi and zeng (2009) - dry soil resistance to evaporation
2620 ! if (vgtype==11) then ! modis wetland
2623 fex=min(1.,soilmois(1)/dqm)
2625 psit=psis*fex ** (-bclh)
2626 psit = max(-1.e5, psit)
2627 alfa=min(1.,exp(g*psit/r_v/soilt))
2631 fc=max(qmin,ref*0.5)
2633 if((soilmois(1)+qmin) > fc .or. (qvatm-qvg) > 0.) then
2636 fex_fc=min(1.,(soilmois(1)+qmin)/fc)
2637 fex_fc=max(fex_fc,0.01)
2638 soilres=0.25*(1.-cos(piconst*fex_fc))**2.
2640 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
2641 ! if (i==421.and.j==280) then
2642 print *,'fex,psit,psis,bclh,g,r_v,soilt,alfa,mavail,soilmois(1),fc,ref,soilres,fex_fc', &
2643 fex,psit,psis,bclh,g,r_v,soilt,alfa,mavail,soilmois(1),fc,ref,soilres,fex_fc
2646 !**************************************************************
2647 ! soiltemp soilves heat budget and diffusion eqn. in soil
2648 !**************************************************************
2651 !--- input variables
2653 delt,ktau,conflx,nzs,nddzs,nroot, &
2655 patm,tabs,qvatm,qcatm,emiss,rnet, &
2656 qkms,tkms,pc,rho,vegfrac, lai, &
2657 thdif,cap,drycan,wetcan, &
2658 transum,dew,mavail,soilres,alfa, &
2659 !--- soil fixed fields
2660 dqm,qmin,bclh,zsmain,zshalf,dtdzs,tbq, &
2662 xlv,cp,g0_p,cvw,stbolt, &
2663 !--- output variables
2664 tso,soilt,qvg,qsg,qcg,x)
2666 !************************************************************************
2668 !--- calculation of dew using new value of qsg or transp if no dew
2672 if(qvatm.ge.qsg)then
2673 dew=qkms*(qvatm-qsg)
2681 transp(k)=vegfrac*ras*qkms* &
2683 tranf(k)*drycan/zshalf(nroot+1)
2684 if(transp(k).gt.0.) transp(k)=0.
2692 !-- recalculate volumetric content of frozen water in soil
2695 tln=log(tso(k)/273.15)
2697 soiliqw(k)=(dqm+qmin)*(xlmelt* &
2698 (tso(k)-273.15)/tso(k)/9.81/psis) &
2700 soiliqw(k)=max(0.,soiliqw(k))
2701 soiliqw(k)=min(soiliqw(k),soilmois(k))
2702 soilice(k)=(soilmois(k)-soiliqw(k))/riw
2703 !---- melting and freezing is balanced, soil ice cannot increase
2704 if(keepfr(k).eq.1.) then
2705 soilice(k)=min(soilice(k),smfrkeep(k))
2706 soiliqw(k)=max(0.,soilmois(k)-soilice(k)*riw)
2711 soiliqw(k)=soilmois(k)
2715 !*************************************************************************
2716 ! soilmoist solves moisture budget (Smirnova et al., 1996, eq.22,28)
2718 !*************************************************************************
2721 delt,nzs,nddzs,dtdzs,dtdzs2,riw, &
2722 zsmain,zshalf,diffu,hydro, &
2723 qsg,qvg,qcg,qcatm,qvatm,-infwater, &
2724 qkms,transp,drip,dew,0.,soilice,vegfrac, &
2727 dqm,qmin,ref,ksat,ras,infmax, &
2729 soilmois,soiliqw,mavail,runoff1, &
2732 !--- keepfr is 1 when the temperature and moisture in soil
2733 !--- are both increasing. in this case soil ice should not
2734 !--- be increasing according to the freezing curve.
2735 !--- some part of ice is melted, but additional water is
2736 !--- getting frozen. thus, only structure of frozen soil is
2737 !--- changed, and phase changes are not affecting the heat
2738 !--- transfer. this situation may happen when it rains on the
2742 if (soilice(k).gt.0.) then
2743 if(tso(k).gt.told(k).and.soilmois(k).gt.smold(k)) then
2751 !--- the diagnostics of surface fluxes
2753 t3 = stbolt*soiltold*soiltold*soiltold
2754 upflux = t3 * 0.5*(soiltold+soilt)
2755 xinet = emiss*(glw-upflux)
2756 hft=-tkms*cp*rho*(tabs-soilt)
2757 hfx=-tkms*cp*rho*(tabs-soilt) &
2758 *(p1000mb*0.00001/patm)**rovcp
2759 q1=-qkms*ras*(qvatm - qsg)
2768 !-- moisture flux for coupling with myj pbl
2769 eeta=-qkms*ras*(qvatm/(1.+qvatm) - qsg/(1.+qsg))*1.e3
2770 cst= cst-eeta*delt*vegfrac
2771 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
2772 print *,'cond myj eeta',eeta,eeta*xlv, i,j
2775 !-- actual moisture flux from ruc lsm
2777 cst=cst+delt*dew*ras * vegfrac
2778 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
2779 print *,'cond ruc lsm eeta',eeta,eeta*xlv, i,j
2786 edir1 =-soilres*(1.-vegfrac)*qkms*ras* &
2789 ec1 = q1 * wetcan * vegfrac
2790 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
2791 print *,'cst before update=',cst
2792 print *,'ec1=',ec1,'cmc2ms=',cmc2ms
2796 cst=max(0.,cst-ec1 * delt)
2799 !-- moisture flux for coupling with myj pbl
2800 eeta=-soilres*qkms*ras*(qvatm/(1.+qvatm) - qvg/(1.+qvg))*1.e3
2802 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
2803 print *,'qkms,ras,qvatm/(1.+qvatm),qvg/(1.+qvg),qsg ', &
2804 qkms,ras,qvatm/(1.+qvatm),qvg/(1.+qvg),qsg
2805 print *,'q1*(1.-vegfrac),edir1',q1*(1.-vegfrac),edir1
2806 print *,'cst,wetcan,drycan',cst,wetcan,drycan
2807 print *,'ec1=',ec1,'ett1=',ett1,'cmc2ms=',cmc2ms,'cmc2ms*ras=',cmc2ms*ras
2809 !-- actual moisture flux from ruc lsm
2810 eeta = (edir1 + ec1 + ett1)*1.e3
2811 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
2812 print *,'ruc lsm eeta',eeta,eeta*xlv
2816 eeta = (edir1 + ec1 + ett1)*1.e3
2818 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
2819 print *,'potential temp hft ',hft
2820 print *,'abs temp hfx ',hfx
2824 s=thdif(1)*cap(1)*dzstop*(tso(1)-tso(2))
2826 fltot=rnet-hft-xlv*eeta-s-x
2827 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
2828 print *,'soil - fltot,rnet,hft,qfx,s,x=',i,j,fltot,rnet,hft,xlv*eeta,s,x
2829 print *,'edir1,ec1,ett1,mavail,qkms,qvatm,qvg,qsg,vegfrac',&
2830 edir1,ec1,ett1,mavail,qkms,qvatm,qvg,qsg,vegfrac
2832 if(detal(1) .ne. 0.) then
2833 ! smf - energy of phase change in the first soil layer
2834 ! smf=xlmelt*1.e3*(soiliqwm(1)-soiliqwmold(1))/delt
2836 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
2837 print *,'detal(1),xlmelt,soiliqwm(1),delt',detal(1),xlmelt,soiliqwm(1),delt
2838 print *,'implicit phase change in the first layer - smf=',smf
2845 1123 format(i5,8f12.3)
2846 1133 format(i7,8e12.4)
2847 123 format(i6,f6.2,7f8.1)
2848 122 format(1x,2i3,6f8.1,f8.3,f8.2)
2849 !-------------------------------------------------------------------
2851 !-------------------------------------------------------------------
2854 !--- input variables
2855 i,j,iland,isoil,delt,ktau,conflx,nzs,nddzs,nroot, &
2856 prcpms,rainf,patm,qvatm,qcatm,glw,gsw, &
2857 emiss,rnet,qkms,tkms,rho,myj, &
2858 !--- sea ice parameters
2859 tice,rhosice,capice,thdifice, &
2860 zsmain,zshalf,dtdzs,dtdzs2,tbq, &
2862 xlv,cp,rovcp,cw,stbolt,tabs, &
2863 !--- output variables
2864 tso,dew,soilt,qvg,qsg,qcg, &
2865 eeta,qfx,hfx,s,evapl,prcpl,fltot &
2868 !*****************************************************************
2869 ! energy budget and heat diffusion eqns. for
2871 !*************************************************************
2874 !-----------------------------------------------------------------
2876 !--- input variables
2878 integer, intent(in ) :: nroot,ktau,nzs , &
2879 nddzs !nddzs=2*(nzs-2)
2880 integer, intent(in ) :: i,j,iland,isoil
2881 real, intent(in ) :: delt,conflx
2882 logical, intent(in ) :: myj
2883 !--- 3-d atmospheric variables
2885 intent(in ) :: patm, &
2890 intent(in ) :: glw, &
2896 !--- sea ice properties
2897 real, dimension(1:nzs) , &
2905 real, intent(in ) :: &
2910 real, dimension(1:nzs), intent(in) :: zsmain, &
2914 real, dimension(1:nddzs), intent(in) :: dtdzs
2916 real, dimension(1:5001), intent(in) :: tbq
2919 !--- input/output variables
2920 !----soil temperature
2921 real, dimension( 1:nzs ), intent(inout) :: tso
2922 !-------- 2-d variables
2924 intent(inout) :: dew, &
2937 !--- local variables
2938 real :: x,x1,x2,x4,tn,denom
2939 real :: rainf, prcpms , &
2940 tabs, t3, upflux, xinet
2942 real :: cp,rovcp,g0,lv,stbolt,xlmelt,dzstop , &
2943 epot,fltot,ft,fq,hft,ras,cvw
2945 real :: fkt,d1,d2,d9,d10,did,r211,r21,r22,r6,r7,d11 , &
2946 pi,h,fkq,r210,aa,bb,pp,q1,qs1,ts1,tq2,tx2 , &
2949 real :: aa1,rhcs, icemelt
2952 real, dimension(1:nzs) :: cotso,rhtso
2954 integer :: nzs1,nzs2,k,k1,kn,kk
2956 !-----------------------------------------------------------------
2958 !-- define constants
2966 dzstop=1./(zsmain(2)-zsmain(1))
2980 x1=dtdzs(k1)*thdifice(kn-1)
2981 x2=dtdzs(k1+1)*thdifice(kn)
2982 ft=tso(kn)+x1*(tso(kn-1)-tso(kn)) &
2983 -x2*(tso(kn)-tso(kn+1))
2984 denom=1.+x1+x2-x2*cotso(k)
2986 rhtso(k+1)=(ft+x2*rhtso(k))/denom
2989 !************************************************************************
2990 !--- the heat balance equation (Smirnova et al., 1996, eq. 21,26)
2997 d9=thdifice(1)*rhcs*dzstop
3001 r22=.5/(thdifice(1)*delt*dzstop**2)
3002 r6=emiss *stbolt*.5*tn**4
3005 tdenom=d9*(1.-d1+r22)+d10+r21+r7 &
3009 aa=xls*(fkq+r210)/tdenom
3010 bb=(d10*tabs+r21*tn+xls*(qvatm*fkq &
3011 +r210*qvg)+d11+d9*(d2+r22*tn) &
3012 +rainf*cvw*prcpms*max(273.15,tabs) &
3017 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
3018 print *,' vilka-seaice1'
3019 print *,'d10,tabs,r21,tn,qvatm,fkq', &
3020 d10,tabs,r21,tn,qvatm,fkq
3021 print *,'rnet, emiss, stbolt, soilt',rnet, emiss, stbolt, soilt
3022 print *,'r210,qvg,d11,d9,d2,r22,rainf,cvw,prcpms,tdenom', &
3023 r210,qvg,d11,d9,d2,r22,rainf,cvw,prcpms,tdenom
3024 print *,'tn,aa1,bb,pp,fkq,r210', &
3025 tn,aa1,bb,pp,fkq,r210
3028 call vilka(tn,aa1,bb,pp,qs1,ts1,tbq,ktau,i,j,iland,isoil)
3029 !--- it is saturation over sea ice
3032 tso(1)=min(271.4,ts1)
3034 !--- sea ice melting is not included in this simple approach
3035 !--- soilt - skin temperature
3037 !---- final solution for soil temperature - tso
3040 tso(k)=min(271.4,rhtso(kk)+cotso(kk)*tso(k-1))
3042 !--- calculation of dew using new value of qsg or transp if no dew
3045 !--- the diagnostics of surface fluxes
3046 t3 = stbolt*tn*tn*tn
3047 upflux = t3 *0.5*(tn+soilt)
3048 xinet = emiss*(glw-upflux)
3049 hft=-tkms*cp*rho*(tabs-soilt)
3050 hfx=-tkms*cp*rho*(tabs-soilt) &
3051 *(p1000mb*0.00001/patm)**rovcp
3052 q1=-qkms*ras*(qvatm - qsg)
3056 !-- moisture flux for coupling with myj pbl
3057 eeta=-qkms*ras*(qvatm/(1.+qvatm) - qsg/(1.+qsg))*1.e3
3058 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
3059 print *,'myj eeta',eeta
3062 !-- actual moisture flux from ruc lsm
3063 dew=qkms*(qvatm-qsg)
3065 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
3066 print *,'ruc lsm eeta',eeta
3074 !-- moisture flux for coupling with myj pbl
3075 eeta=-qkms*ras*(qvatm/(1.+qvatm) - qvg/(1.+qvg))*1.e3
3076 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
3077 print *,'myj eeta',eeta
3080 ! to convert from m s-1 to kg m-2 s-1: *rho water=1.e3************
3081 !-- actual moisture flux from ruc lsm
3083 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
3084 print *,'ruc lsm eeta',eeta
3092 s=thdifice(1)*capice(1)*dzstop*(tso(1)-tso(2))
3093 ! heat storage in surface layer
3096 x= (cp*rho*r211+rhcs*zsmain(2)*0.5/delt)*(soilt-tn) + &
3097 xls*rho*r211*(qsg-qgold)
3100 -rainf*cvw*prcpms*(max(273.15,tabs)-soilt)
3102 !-- excess energy spent on sea ice melt
3103 icemelt=rnet-xls*eeta -hft -s -x
3104 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
3105 print *,'icemelt=',icemelt
3108 fltot=rnet-xls*eeta-hft-s-x-icemelt
3109 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
3110 print *,'sice - fltot,rnet,hft,qfx,s,snoh,x=', &
3111 fltot,rnet,hft,xls*eeta,s,icemelt,x
3114 !-------------------------------------------------------------------
3116 !-------------------------------------------------------------------
3120 subroutine snowsoil (spp_lsm,rstochcol,fieldcol_sf,&
3121 !--- input variables
3122 i,j,isoil,delt,ktau,conflx,nzs,nddzs,nroot, &
3123 meltfactor,rhonewsn,snhei_crit, & ! new
3124 iland,prcpms,rainf,newsnow,snhei,snwe,snowfrac, &
3127 glw,gsw,gswin,emiss,rnet,ivgtyp, &
3128 qkms,tkms,pc,cst,drip,infwater, &
3129 rho,vegfrac,alb,znt,lai, &
3131 !--- soil fixed fields
3132 qwrtz,rhocs,dqm,qmin,ref,wilt,psis,bclh,ksat, &
3133 sat,cn,zsmain,zshalf,dtdzs,dtdzs2,tbq, &
3135 xlv,cp,rovcp,g0_p,cw,stbolt,tabs, &
3137 !--- output variables
3138 ilnb,snweprint,snheiprint,rsm, &
3139 soilmois,tso,smfrkeep,keepfr, &
3140 dew,soilt,soilt1,tsnav, &
3141 qvg,qsg,qcg,smelt,snoh,snflx,snom, &
3142 edir1,ec1,ett1,eeta,qfx,hfx,s,sublim, &
3143 prcpl,fltot,runoff1,runoff2,mavail,soilice, &
3146 !***************************************************************
3147 ! energy and moisture budget for snow, heat diffusion eqns.
3148 ! in snow and soil, richards eqn. for soil covered with snow
3150 ! delt - time step (s)
3151 ! ktau - numver of time step
3152 ! conflx - depth of constant flux layer (m)
3153 ! j,i - the location of grid point
3154 ! ime, jme, nzs - dimensions of the domain
3155 ! nroot - number of levels within the root zone
3156 ! prcpms - precipitation rate in m/s
3157 ! newsnow - pcpn in soilid form (m)
3158 ! snhei, snwe - snow height and snow water equivalent (m)
3159 ! rhosn - snow density (kg/m-3)
3160 ! patm - pressure (bar)
3161 ! qvatm,qcatm - cloud and water vapor mixing ratio
3162 ! at the first atm. level (kg/kg)
3163 ! glw, gsw - incoming longwave and absorbed shortwave
3164 ! radiation at the surface (w/m^2)
3165 ! emiss,rnet - emissivity (0-1) of the ground surface and net
3166 ! radiation at the surface (w/m^2)
3167 ! qkms - exchange coefficient for water vapor in the
3168 ! surface layer (m/s)
3169 ! tkms - exchange coefficient for heat in the surface
3171 ! pc - plant coefficient (resistance) (0-1)
3172 ! rho - density of atmosphere near surface (kg/m^3)
3173 ! vegfrac - greeness fraction (0-1)
3174 ! rhocs - volumetric heat capacity of dry soil (j/m^3/k)
3175 ! dqm, qmin - porosity minus residual soil moisture qmin (m^3/m^3)
3176 ! ref, wilt - field capacity soil moisture and the
3177 ! wilting point (m^3/m^3)
3178 ! psis - matrix potential at saturation (m)
3179 ! bclh - exponent for clapp-hornberger parameterization
3180 ! ksat - saturated hydraulic conductivity (m/s)
3181 ! sat - maximum value of water intercepted by canopy (m)
3182 ! cn - exponent for calculation of canopy water
3183 ! zsmain - main levels in soil (m)
3184 ! zshalf - middle of the soil layers (m)
3185 ! dtdzs,dtdzs2 - dt/(2.*dzshalf*dzmain) and dt/dzshalf in soil
3186 ! tbq - table to define saturated mixing ration
3187 ! of water vapor for given temperature and pressure
3188 ! ilnb - number of layers in snow
3189 ! rsm - liquid water inside snow pack (m)
3190 ! soilmois,tso - soil moisture (m^3/m^3) and temperature (k)
3191 ! dew - dew in (kg/m^2 s)
3192 ! soilt - skin temperature (k)
3193 ! soilt1 - snow temperature at 7.5 cm depth (k)
3194 ! tsnav - average temperature of snow pack (c)
3195 ! qsg,qvg,qcg - saturated mixing ratio, mixing ratio of
3196 ! water vapor and cloud at the ground
3197 ! surface, respectively (kg/kg)
3198 ! edir1, ec1, ett1, eeta - direct evaporation, evaporation of
3199 ! canopy water, transpiration (kg m-2 s-1) and total
3200 ! evaporation in (m s-1).
3201 ! qfx, hfx - latent and sensible heat fluxes (w/m^2)
3202 ! s - soil heat flux in the top layer (w/m^2)
3203 ! sublim - snow sublimation (kg/m^2/s)
3204 ! runoff1 - surface runoff (m/s)
3205 ! runoff2 - underground runoff (m)
3206 ! mavail - moisture availability in the top soil layer (0-1)
3207 ! soilice - content of soil ice in soil layers (m^3/m^3)
3208 ! soiliqw - lliquid water in soil layers (m^3/m^3)
3209 ! infiltrp - infiltration flux from the top of soil domain (m/s)
3210 ! xinet - net long-wave radiation (w/m^2)
3212 !*******************************************************************
3215 !-------------------------------------------------------------------
3216 !--- input variables
3218 integer, intent(in ) :: nroot,ktau,nzs , &
3219 nddzs !nddzs=2*(nzs-2)
3220 integer, intent(in ) :: i,j,isoil
3222 real, intent(in ) :: delt,conflx,prcpms , &
3223 rainf,newsnow,rhonewsn, &
3224 snhei_crit,meltfactor
3226 logical, intent(in ) :: myj
3228 !--- 3-d atmospheric variables
3230 intent(in ) :: patm, &
3235 intent(in ) :: glw, &
3246 integer, intent(in ) :: ivgtyp
3247 !--- soil properties
3249 intent(in ) :: rhocs, &
3260 real, intent(in ) :: cn, &
3269 real, dimension(1:nzs), intent(in) :: zsmain, &
3273 real, dimension(1:nddzs), intent(in) :: dtdzs
3275 real, dimension(1:5001), intent(in) :: tbq
3277 real, dimension(1:nzs), intent(in) :: rstochcol
3278 real, dimension(1:nzs), intent(inout) :: fieldcol_sf
3280 !--- input/output variables
3281 !-------- 3-d soil moisture and temperature
3282 real, dimension( 1:nzs ) , &
3283 intent(inout) :: tso, &
3287 real, dimension( 1:nzs ) , &
3288 intent(inout) :: keepfr
3291 integer, intent(inout) :: iland
3294 !-------- 2-d variables
3296 intent(inout) :: dew, &
3329 integer, intent(inout) :: ilnb
3331 !-------- 1-d variables
3332 real, dimension(1:nzs), intent(out) :: soilice, &
3335 real, intent(out) :: rsm, &
3338 integer, intent(in) :: spp_lsm
3339 !--- local variables
3342 integer :: nzs1,nzs2,k
3344 real :: infiltrp, transum , &
3346 tabs, t3, upflux, xinet , &
3347 beta, snwepr,epdt,pp
3348 real :: cp,rovcp,g0,lv,xlvm,stbolt,xlmelt,dzstop , &
3349 can,epot,fac,fltot,ft,fq,hft , &
3350 q1,ras,rhoice,sph , &
3351 trans,zn,ci,cvw,tln,tavln,pi , &
3352 dd1,cmc2ms,drycan,wetcan , &
3353 infmax,riw,deltsn,h,umveg
3355 real, dimension(1:nzs) :: transp,cap,diffu,hydro , &
3356 thdif,tranf,tav,soilmoism , &
3357 soilicem,soiliqwm,detal , &
3358 fwsat,lwsat,told,smold
3359 real :: soiltold, qgold
3363 !-----------------------------------------------------------------
3367 !-- heat of water vapor sublimation
3370 !--- snow flag -- isice
3373 !--- deltsn - is the threshold for splitting the snow layer into 2 layers.
3374 !--- with snow density 400 kg/m^3, this threshold is equal to 7.5 cm,
3375 !--- equivalent to 0.03 m snwe. for other snow densities the threshold is
3376 !--- computed using snwe=0.03 m and current snow density.
3377 !--- snth - the threshold below which the snow layer is combined with
3378 !--- the top soil layer. snth is computed using snwe=0.016 m, and
3379 !--- equals 4 cm for snow density 400 kg/m^3.
3387 deltsn=0.05*1.e3/rhosn
3388 snth=0.01*1.e3/rhosn
3389 ! print *,'deltsn,snhei,snth',i,j,deltsn,snhei,snth
3391 ! for 2-layer snow model when the snow depth is marginally higher than deltsn,
3392 ! reset deltsn to half of snow depth.
3393 if(snhei.ge.deltsn+snth) then
3394 if(snhei-deltsn-snth.lt.snth) deltsn=0.5*(snhei-snth)
3395 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
3396 print *,'deltsn is changed,deltsn,snhei,snth',i,j,deltsn,snhei,snth
3429 !*** deltsn is the depth of the top layer of snow where
3430 !*** there is a temperature gradient, the rest of the snow layer
3431 !*** is considered to have constant temperature
3436 dzstop=1./(zsmain(2)-zsmain(1))
3438 !----- the calculation of thermal diffusivity, diffusional and ---
3439 !----- hydraulic conductivity (Smirnova et al. 1996, eq.2,5,6) ---
3440 !tgs - the following loop is added to define the amount of frozen
3441 !tgs - water in soil if there is any
3444 tln=log(tso(k)/273.15)
3446 soiliqw(k)=(dqm+qmin)*(xlmelt* &
3447 (tso(k)-273.15)/tso(k)/9.81/psis) &
3449 soiliqw(k)=max(0.,soiliqw(k))
3450 soiliqw(k)=min(soiliqw(k),soilmois(k))
3451 soilice(k)=(soilmois(k)-soiliqw(k))/riw
3453 !---- melting and freezing is balanced, soil ice cannot increase
3454 if(keepfr(k).eq.1.) then
3455 soilice(k)=min(soilice(k),smfrkeep(k))
3456 soiliqw(k)=max(0.,soilmois(k)-soilice(k)*rhoice*1.e-3)
3461 soiliqw(k)=soilmois(k)
3468 tav(k)=0.5*(tso(k)+tso(k+1))
3469 soilmoism(k)=0.5*(soilmois(k)+soilmois(k+1))
3470 tavln=log(tav(k)/273.15)
3472 if(tavln.lt.0.) then
3473 soiliqwm(k)=(dqm+qmin)*(xlmelt* &
3474 (tav(k)-273.15)/tav(k)/9.81/psis) &
3476 fwsat(k)=dqm-soiliqwm(k)
3477 lwsat(k)=soiliqwm(k)+qmin
3478 soiliqwm(k)=max(0.,soiliqwm(k))
3479 soiliqwm(k)=min(soiliqwm(k), soilmoism(k))
3480 soilicem(k)=(soilmoism(k)-soiliqwm(k))/riw
3481 !---- melting and freezing is balanced, soil ice cannot increase
3482 if(keepfr(k).eq.1.) then
3483 soilicem(k)=min(soilicem(k), &
3484 0.5*(smfrkeep(k)+smfrkeep(k+1)))
3485 soiliqwm(k)=max(0.,soilmoism(k)-soilicem(k)*riw)
3486 fwsat(k)=dqm-soiliqwm(k)
3487 lwsat(k)=soiliqwm(k)+qmin
3492 soiliqwm(k)=soilmoism(k)
3500 if(soilice(k).gt.0.) then
3501 smfrkeep(k)=soilice(k)
3503 smfrkeep(k)=soilmois(k)/riw
3506 !******************************************************************
3507 ! soilprop computes thermal diffusivity, and diffusional and
3508 ! hydraulic condeuctivities
3509 !******************************************************************
3510 call soilprop(spp_lsm,rstochcol,fieldcol_sf, &
3511 !--- input variables
3512 nzs,fwsat,lwsat,tav,keepfr, &
3513 soilmois,soiliqw,soilice, &
3514 soilmoism,soiliqwm,soilicem, &
3515 !--- soil fixed fields
3516 qwrtz,rhocs,dqm,qmin,psis,bclh,ksat, &
3518 riw,xlmelt,cp,g0_p,cvw,ci, &
3520 !--- output variables
3521 thdif,diffu,hydro,cap)
3523 !********************************************************************
3524 !--- calculation of canopy water (Smirnova et al., 1996, eq.16) and dew
3532 !--- if vegfrac.ne.0. then part of falling snow can be
3533 !--- intercepted by the canopy.
3537 epot = -fq*(qvatm-qsg)
3539 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
3540 print *,'snwe after subtracting intercepted snow - snwe=',snwe,vegfrac,cst
3544 ! check if all snow can evaporate during dt
3546 epdt = epot * ras *delt*umveg
3547 if(epdt.gt.0. .and. snwepr.le.epdt) then
3548 beta=snwepr/max(1.e-8,epdt)
3552 wetcan=min(0.25,max(0.,(cst/sat))**cn)
3553 ! if(lai > 1.) wetcan=wetcan/lai
3556 !**************************************************************
3557 ! transf computes transpiration function
3558 !**************************************************************
3560 !--- input variables
3561 nzs,nroot,soiliqw,tabs,lai,gswin, &
3562 !--- soil fixed fields
3563 dqm,qmin,ref,wilt,zshalf,pc,iland, &
3564 !--- output variables
3567 !--- save soil temp and moisture from the beginning of time step
3570 smold(k)=soilmois(k)
3573 !**************************************************************
3574 ! snowtemp solves heat budget and diffusion eqn. in soil
3575 !**************************************************************
3577 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
3578 print *, 'tso before calling snowtemp: ', tso
3581 !--- input variables
3583 delt,ktau,conflx,nzs,nddzs,nroot, &
3584 snwe,snwepr,snhei,newsnow,snowfrac, &
3585 beta,deltsn,snth,rhosn,rhonewsn,meltfactor, & ! add meltfactor
3587 patm,tabs,qvatm,qcatm, &
3588 glw,gsw,emiss,rnet, &
3589 qkms,tkms,pc,rho,vegfrac, &
3590 thdif,cap,drycan,wetcan,cst, &
3591 tranf,transum,dew,mavail, &
3592 !--- soil fixed fields
3593 dqm,qmin,psis,bclh, &
3594 zsmain,zshalf,dtdzs,tbq, &
3596 xlvm,cp,rovcp,g0_p,cvw,stbolt, &
3597 !--- output variables
3598 snweprint,snheiprint,rsm, &
3599 tso,soilt,soilt1,tsnav,qvg,qsg,qcg, &
3600 smelt,snoh,snflx,s,ilnb,x)
3602 !************************************************************************
3603 !--- recalculation of dew using new value of qsg or transp if no dew
3607 epot = -fq*(qvatm-qsg)
3611 transp(k)=vegfrac*ras*fq*(qvatm-qsg) &
3612 *tranf(k)*drycan/zshalf(nroot+1)
3628 !-- recalculating of frozen water in soil
3630 tln=log(tso(k)/273.15)
3632 soiliqw(k)=(dqm+qmin)*(xlmelt* &
3633 (tso(k)-273.15)/tso(k)/9.81/psis) &
3635 soiliqw(k)=max(0.,soiliqw(k))
3636 soiliqw(k)=min(soiliqw(k),soilmois(k))
3637 soilice(k)=(soilmois(k)-soiliqw(k))/riw
3638 !---- melting and freezing is balanced, soil ice cannot increase
3639 if(keepfr(k).eq.1.) then
3640 soilice(k)=min(soilice(k),smfrkeep(k))
3641 soiliqw(k)=max(0.,soilmois(k)-soilice(k)*riw)
3646 soiliqw(k)=soilmois(k)
3650 !*************************************************************************
3651 !--- tqcan for solution of moisture balance (Smirnova et al. 1996, eq.22,28)
3652 ! and tso,eta profiles
3653 !*************************************************************************
3656 delt,nzs,nddzs,dtdzs,dtdzs2,riw, &
3657 zsmain,zshalf,diffu,hydro, &
3658 qsg,qvg,qcg,qcatm,qvatm,-infwater, &
3660 0.,smelt,soilice,vegfrac, &
3663 dqm,qmin,ref,ksat,ras,infmax, &
3665 soilmois,soiliqw,mavail,runoff1, &
3670 !-- restore land-use parameters if all snow is melted
3671 if(snhei.eq.0.) then
3676 ! snom [mm] goes into the passed-in acsnom variable in the grid derived type
3677 snom=snom+smelt*delt*1.e3
3679 !--- keepfr is 1 when the temperature and moisture in soil
3680 !--- are both increasing. in this case soil ice should not
3681 !--- be increasing according to the freezing curve.
3682 !--- some part of ice is melted, but additional water is
3683 !--- getting frozen. thus, only structure of frozen soil is
3684 !--- changed, and phase changes are not affecting the heat
3685 !--- transfer. this situation may happen when it rains on the
3689 if (soilice(k).gt.0.) then
3690 if(tso(k).gt.told(k).and.soilmois(k).gt.smold(k)) then
3697 !--- the diagnostics of surface fluxes
3699 t3 = stbolt*soiltold*soiltold*soiltold
3700 upflux = t3 *0.5*(soiltold+soilt)
3701 xinet = emiss*(glw-upflux)
3702 hfx=-tkms*cp*rho*(tabs-soilt) &
3703 *(p1000mb*0.00001/patm)**rovcp
3704 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
3705 print *,'potential temp hfx',hfx
3707 hft=-tkms*cp*rho*(tabs-soilt)
3708 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
3709 print *,'abs temp hfx',hft
3711 q1 = - fq*ras* (qvatm - qsg)
3720 !-- moisture flux for coupling with myj pbl
3721 eeta=-qkms*ras*(qvatm/(1.+qvatm) - qsg/(1.+qsg))*1.e3
3722 cst= cst-eeta*delt*vegfrac
3723 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
3724 print *,'myj eeta cond', eeta
3727 !-- actual moisture flux from ruc lsm
3728 dew=qkms*(qvatm-qsg)
3730 cst=cst+delt*dew*ras * vegfrac
3731 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
3732 print *,'ruc lsm eeta cond',eeta
3739 edir1 = q1*umveg *beta
3741 ec1 = q1 * wetcan * vegfrac
3743 cst=max(0.,cst-ec1 * delt)
3745 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
3746 print*,'q1,umveg,beta',q1,umveg,beta
3747 print *,'wetcan,vegfrac',wetcan,vegfrac
3748 print *,'ec1,cmc2ms',ec1,cmc2ms
3752 !-- moisture flux for coupling with myj pbl
3753 eeta=-(qkms*ras*(qvatm/(1.+qvatm) - qsg/(1.+qsg))*1.e3)*beta
3754 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
3755 print *,'myj eeta', eeta*xlvm,eeta
3758 ! to convert from m s-1 to kg m-2 s-1: *rho water=1.e3************
3759 !-- actual moisture flux from ruc lsm
3760 eeta = (edir1 + ec1 + ett1)*1.e3
3761 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
3762 print *,'ruc lsm eeta',eeta*xlvm,eeta
3766 eeta = (edir1 + ec1 + ett1)*1.e3
3771 fltot=rnet-hft-xlvm*eeta-s-snoh-x
3772 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
3773 print *,'snowsoil - fltot,rnet,hft,qfx,s,snoh,x=',fltot,rnet,hft,xlvm*eeta,s,snoh,x
3774 print *,'edir1,ec1,ett1,mavail,qkms,qvatm,qvg,qsg,vegfrac,beta',&
3775 edir1,ec1,ett1,mavail,qkms,qvatm,qvg,qsg,vegfrac,beta
3780 1123 format(i5,8f12.3)
3781 1133 format(i7,8e12.4)
3782 123 format(i6,f6.2,7f8.1)
3783 122 format(1x,2i3,6f8.1,f8.3,f8.2)
3785 !-------------------------------------------------------------------
3786 end subroutine snowsoil
3787 !-------------------------------------------------------------------
3789 subroutine snowseaice( &
3790 i,j,isoil,delt,ktau,conflx,nzs,nddzs, &
3791 meltfactor,rhonewsn,snhei_crit, & ! new
3792 iland,prcpms,rainf,newsnow,snhei,snwe,snowfrac, &
3793 rhosn,patm,qvatm,qcatm, &
3794 glw,gsw,emiss,rnet, &
3795 qkms,tkms,rho,myj, &
3796 !--- sea ice parameters
3798 tice,rhosice,capice,thdifice, &
3799 zsmain,zshalf,dtdzs,dtdzs2,tbq, &
3801 xlv,cp,rovcp,cw,stbolt,tabs, &
3802 !--- output variables
3803 ilnb,snweprint,snheiprint,rsm,tso, &
3804 dew,soilt,soilt1,tsnav,qvg,qsg,qcg, &
3805 smelt,snoh,snflx,snom,eeta, &
3806 qfx,hfx,s,sublim,prcpl,fltot &
3808 !***************************************************************
3809 ! solving energy budget for snow on sea ice and heat diffusion
3810 ! eqns. in snow and sea ice
3811 !***************************************************************
3815 !-------------------------------------------------------------------
3816 !--- input variables
3818 integer, intent(in ) :: ktau,nzs , &
3819 nddzs !nddzs=2*(nzs-2)
3820 integer, intent(in ) :: i,j,isoil
3822 real, intent(in ) :: delt,conflx,prcpms , &
3823 rainf,newsnow,rhonewsn, &
3824 meltfactor, snhei_crit
3827 logical, intent(in ) :: myj
3828 !--- 3-d atmospheric variables
3830 intent(in ) :: patm, &
3835 intent(in ) :: glw, &
3841 !--- sea ice properties
3842 real, dimension(1:nzs) , &
3849 real, intent(in ) :: &
3853 real, dimension(1:nzs), intent(in) :: zsmain, &
3857 real, dimension(1:nddzs), intent(in) :: dtdzs
3859 real, dimension(1:5001), intent(in) :: tbq
3861 !--- input/output variables
3862 !-------- 3-d soil moisture and temperature
3863 real, dimension( 1:nzs ) , &
3864 intent(inout) :: tso
3866 integer, intent(inout) :: iland
3869 !-------- 2-d variables
3871 intent(inout) :: dew, &
3896 integer, intent(inout) :: ilnb
3898 real, intent(out) :: rsm, &
3901 !--- local variables
3904 integer :: nzs1,nzs2,k,k1,kn,kk
3905 real :: x,x1,x2,dzstop,ft,tn,denom
3907 real :: snth, newsn , &
3908 tabs, t3, upflux, xinet , &
3909 beta, snwepr,epdt,pp
3910 real :: cp,rovcp,g0,lv,xlvm,stbolt,xlmelt , &
3911 epot,fltot,fq,hft,q1,ras,rhoice,ci,cvw , &
3914 real :: rhocsn,thdifsn, &
3915 xsn,ddzsn,x1sn,d1sn,d2sn,d9sn,r22sn
3917 real :: cotsn,rhtsn,xsn1,ddzsn1,x1sn1,ftsnow,denomsn
3919 fkt,d1,d2,d9,d10,did,r211,r21,r22,r6,r7,d11, &
3920 fkq,r210,aa,bb,qs1,ts1,tq2,tx2, &
3921 tdenom,aa1,rhcs,h1,tsob, snprim, &
3922 snodif,soh,tnold,qgold,snohgnew
3923 real, dimension(1:nzs) :: cotso,rhtso
3925 real :: rnet,rsmfrac,soiltfrac,hsn,icemelt,rr
3929 !-----------------------------------------------------------------
3931 !-- heat of sublimation of water vapor
3934 !--- snow flag -- isice
3937 !--- deltsn - is the threshold for splitting the snow layer into 2 layers.
3938 !--- with snow density 400 kg/m^3, this threshold is equal to 7.5 cm,
3939 !--- equivalent to 0.03 m snwe. for other snow densities the threshold is
3940 !--- computed using snwe=0.03 m and current snow density.
3941 !--- snth - the threshold below which the snow layer is combined with
3942 !--- the top sea ice layer. snth is computed using snwe=0.016 m, and
3943 !--- equals 4 cm for snow density 400 kg/m^3.
3945 deltsn=0.05*1.e3/rhosn
3946 snth=0.01*1.e3/rhosn
3948 ! for 2-layer snow model when the snow depth is marginlly higher than deltsn,
3949 ! reset deltsn to half of snow depth.
3950 if(snhei.ge.deltsn+snth) then
3951 if(snhei-deltsn-snth.lt.snth) deltsn=0.5*(snhei-snth)
3952 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
3953 print *,'deltsn ice is changed,deltsn,snhei,snth', &
3954 i,j, deltsn,snhei,snth
3966 !18apr08 - add rhonewcsn
3967 rhonewcsn=2090.* rhonewsn
3968 thdifsn = 0.265/rhocsn
3989 dzstop=1./(zsmain(2)-zsmain(1))
3995 !*** deltsn is the depth of the top layer of snow where
3996 !*** there is a temperature gradient, the rest of the snow layer
3997 !*** is considered to have constant temperature
4004 snhei=snwe*1.e3/rhosn
4007 ! check if all snow can evaporate during dt
4009 epot = -fq*(qvatm-qsg)
4010 epdt = epot * ras *delt
4011 if(epdt.gt.0. .and. snwepr.le.epdt) then
4012 beta=snwepr/max(1.e-8,epdt)
4016 !******************************************************************************
4017 ! coefficients for thomas algorithm for tso
4018 !******************************************************************************
4025 x1=dtdzs(k1)*thdifice(kn-1)
4026 x2=dtdzs(k1+1)*thdifice(kn)
4027 ft=tso(kn)+x1*(tso(kn-1)-tso(kn)) &
4028 -x2*(tso(kn)-tso(kn+1))
4029 denom=1.+x1+x2-x2*cotso(k)
4031 rhtso(k+1)=(ft+x2*rhtso(k))/denom
4033 !--- the nzs element in cotso and rhtso will be for snow
4034 !--- there will be 2 layers in snow if it is deeper than deltsn+snth
4035 if(snhei.ge.snth) then
4036 if(snhei.le.deltsn+snth) then
4037 !-- 1-layer snow model
4039 snprim=max(snth,snhei)
4042 xsn = delt/2./(zshalf(2)+0.5*snprim)
4043 ddzsn = xsn / snprim
4044 x1sn = ddzsn * thdifsn
4045 x2 = dtdzs(1)*thdifice(1)
4046 ft = tso(1)+x1sn*(soilt-tso(1)) &
4048 denom = 1. + x1sn + x2 -x2*cotso(nzs1)
4049 cotso(nzs)=x1sn/denom
4050 rhtso(nzs)=(ft+x2*rhtso(nzs1))/denom
4053 !*** average temperature of snow pack (c)
4054 tsnav=0.5*(soilt+tso(1)) &
4058 !-- 2 layers in snow, soilt1 is temperasture at deltsn depth
4062 xsn = delt/2./(0.5*snhei)
4063 xsn1= delt/2./(zshalf(2)+0.5*(snhei-deltsn))
4064 ddzsn = xsn / deltsn
4065 ddzsn1 = xsn1 / (snhei-deltsn)
4066 x1sn = ddzsn * thdifsn
4067 x1sn1 = ddzsn1 * thdifsn
4068 x2 = dtdzs(1)*thdifice(1)
4069 ft = tso(1)+x1sn1*(soilt1-tso(1)) &
4071 denom = 1. + x1sn1 + x2 - x2*cotso(nzs1)
4072 cotso(nzs)=x1sn1/denom
4073 rhtso(nzs)=(ft+x2*rhtso(nzs1))/denom
4074 ftsnow = soilt1+x1sn*(soilt-soilt1) &
4075 -x1sn1*(soilt1-tso(1))
4076 denomsn = 1. + x1sn + x1sn1 - x1sn1*cotso(nzs)
4078 rhtsn=(ftsnow+x1sn1*rhtso(nzs))/denomsn
4079 !*** average temperature of snow pack (c)
4080 tsnav=0.5/snhei*((soilt+soilt1)*deltsn &
4081 +(soilt1+tso(1))*(snhei-deltsn)) &
4086 if(snhei.lt.snth.and.snhei.gt.0.) then
4087 !--- snow is too thin to be treated separately, therefore it
4088 !--- is combined with the first sea ice layer.
4089 snprim=snhei+zsmain(2)
4094 xsn = delt/2./((zshalf(3)-zsmain(2))+0.5*snprim)
4096 x1sn = ddzsn * (fsn*thdifsn+fso*thdifice(1))
4097 x2=dtdzs(2)*thdifice(2)
4098 ft=tso(2)+x1sn*(soilt-tso(2))- &
4100 denom = 1. + x1sn + x2 - x2*cotso(nzs-2)
4101 cotso(nzs1) = x1sn/denom
4102 rhtso(nzs1)=(ft+x2*rhtso(nzs-2))/denom
4103 tsnav=0.5*(soilt+tso(1)) &
4105 cotso(nzs)=cotso(nzs1)
4106 rhtso(nzs)=rhtso(nzs1)
4111 !************************************************************************
4112 !--- the heat balance equation
4113 !18apr08 nmelt is the flag for melting, and snoh is heat of snow phase changes
4117 epot=-qkms*(qvatm-qsg)
4124 d9=thdifice(1)*rhcs*dzstop
4128 r22=.5/(thdifice(1)*delt*dzstop**2)
4129 r6=emiss *stbolt*.5*tn**4
4133 if(snhei.ge.snth) then
4134 if(snhei.le.deltsn+snth) then
4143 d9sn= thdifsn*rhocsn / snprim
4144 r22sn = snprim*snprim*0.5/(thdifsn*delt)
4147 if(snhei.lt.snth.and.snhei.gt.0.) then
4148 !--- thin snow is combined with sea ice
4151 d9sn = (fsn*thdifsn*rhocsn+fso*thdifice(1)*rhcs)/ &
4153 r22sn = snprim*snprim*0.5 &
4154 /((fsn*thdifsn+fso*thdifice(1))*delt)
4158 !--- all snow is sublimated
4166 !---- tdenom for snow
4167 tdenom = d9sn*(1.-d1sn +r22sn)+d10+r21+r7 &
4169 +rhonewcsn*newsnow/delt
4173 aa=xlvm*(beta*fkq+r210)/tdenom
4174 bb=(d10*tabs+r21*tn+xlvm*(qvatm* &
4176 +r210*qvg)+d11+d9sn*(d2sn+r22sn*tn) &
4177 +rainf*cvw*prcpms*max(273.15,tabs) &
4178 + rhonewcsn*newsnow/delt*min(273.15,tabs) &
4183 !18apr08 - the iteration start point
4186 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
4187 print *,'vilka-snow on seaice'
4188 print *,'tn,aa1,bb,pp,fkq,r210', &
4189 tn,aa1,bb,pp,fkq,r210
4190 print *,'tabs,qvatm,tn,qvg=',tabs,qvatm,tn,qvg
4193 call vilka(tn,aa1,bb,pp,qs1,ts1,tbq,ktau,i,j,iland,isoil)
4194 !--- it is saturation over snow
4199 !--- soilt - skin temperature
4202 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
4203 print *,' after vilka-snow on seaice'
4204 print *,' ts1,qs1: ', ts1,qs1
4206 ! solution for temperature at 7.5 cm depth and snow-seaice interface
4207 if(snhei.ge.snth) then
4208 if(snhei.gt.deltsn+snth) then
4209 !-- 2-layer snow model
4210 soilt1=min(273.15,rhtsn+cotsn*soilt)
4211 tso(1)=min(271.4,(rhtso(nzs)+cotso(nzs)*soilt1))
4215 tso(1)=min(271.4,(rhtso(nzs)+cotso(nzs)*soilt))
4219 elseif (snhei > 0. .and. snhei < snth) then
4221 tso(2)=min(271.4,(rhtso(nzs1)+cotso(nzs1)*soilt))
4222 tso(1)=min(271.4,(tso(2)+(soilt-tso(2))*fso))
4227 tso(1)=min(271.4,soilt)
4228 soilt1=min(271.4,soilt)
4231 !---- final solution for tso in sea ice
4232 if (snhei > 0. .and. snhei < snth) then
4233 ! blended or snow is melted
4236 tso(k)=min(271.4,rhtso(kk)+cotso(kk)*tso(k-1))
4241 tso(k)=min(271.4,rhtso(kk)+cotso(kk)*tso(k-1))
4244 !--- for thin snow layer combined with the top soil layer
4245 !--- tso(i,j,1) is computed by linear interpolation between soilt
4247 ! if(snhei.lt.snth.and.snhei.gt.0.)then
4248 ! tso(1)=min(271.4,tso(2)+(soilt-tso(2))*fso)
4253 if(nmelt.eq.1) go to 220
4255 !--- if soilt > 273.15 f then melting of snow can happen
4256 ! if all snow can evaporate, then there is nothing to melt
4257 if(soilt.gt.273.15.and.snwepr-beta*epot*ras*delt.gt.0..and.snhei.gt.0.) then
4260 soiltfrac=snowfrac*273.15+(1.-snowfrac)*min(271.4,soilt)
4262 qsg= qsn(soiltfrac,tbq)/pp
4263 t3 = stbolt*tnold*tnold*tnold
4264 upflux = t3 * 0.5*(tnold+soiltfrac)
4265 xinet = emiss*(glw-upflux)
4266 epot = -qkms*(qvatm-qsg)
4277 eeta = q1 * beta *1.e3
4278 ! to convert from kg m-2 s-1 to m s-1: 1/rho water=1.e-3************
4282 hfx=d10*(tabs-soiltfrac)
4284 if(snhei.ge.snth)then
4285 soh=thdifsn*rhocsn*(soiltfrac-tsob)/snprim
4288 soh=(fsn*thdifsn*rhocsn+fso*thdifice(1)*rhcs)* &
4289 (soiltfrac-tsob)/snprim
4292 x= (r21+d9sn*r22sn)*(soiltfrac-tnold) + &
4293 xlvm*r210*(qsg-qgold)
4294 !-- snoh is energy flux of snow phase change
4295 snoh=rnet+qfx +hfx &
4296 +rhonewcsn*newsnow/delt*(min(273.15,tabs)-soiltfrac) &
4297 -soh-x+rainf*cvw*prcpms* &
4298 (max(273.15,tabs)-soiltfrac)
4300 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
4301 print *,'snowseaice melt i,j,snoh,rnet,qfx,hfx,soh,x',i,j,snoh,rnet,qfx,hfx,soh,x
4302 print *,'rhonewcsn*newsnow/delt*(min(273.15,tabs)-soiltfrac)', &
4303 rhonewcsn*newsnow/delt*(min(273.15,tabs)-soiltfrac)
4304 print *,'rainf*cvw*prcpms*(max(273.15,tabs)-soiltfrac)', &
4305 rainf*cvw*prcpms*(max(273.15,tabs)-soiltfrac)
4308 !-- smelt is speed of melting in m/s
4309 smelt= snoh /xlmelt*1.e-3
4310 smelt=amin1(smelt,snwepr/delt-beta*epot*ras)
4311 smelt=amax1(0.,smelt)
4313 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
4314 print *,'1-smelt i,j',smelt,i,j
4316 !18apr08 - egglston limit
4317 ! smelt= amin1 (smelt, 5.6e-7*meltfactor*max(1.,(soilt-273.15)))
4318 smelt= amin1 (smelt, 5.6e-8*meltfactor*max(1.,(soilt-273.15)))
4319 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
4320 print *,'2-smelt i,j',smelt,i,j
4323 ! rr - potential melting
4324 rr=snwepr/delt-beta*epot*ras
4326 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
4327 print *,'3- smelt i,j,smelt,rr',i,j,smelt,rr
4329 snohgnew=smelt*xlmelt*1.e3
4330 snodif=amax1(0.,(snoh-snohgnew))
4334 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
4335 print*,'soiltfrac,soilt,snohgnew,snodif=', &
4336 i,j,soiltfrac,soilt,snohgnew,snodif
4337 print *,'snoh,snodif',snoh,snodif
4340 !*** from koren et al. (1999) 13% of snow melt stays in the snow pack
4341 rsmfrac=min(0.18,(max(0.08,snwepr/0.10*0.13)))
4342 if(snhei > 0.01) then
4343 rsm=rsmfrac*smelt*delt
4345 ! do not keep melted water if snow depth is less that 1 cm
4348 !18apr08 rsm is part of melted water that stays in snow as liquid
4349 smelt=amax1(0.,smelt-rsm/delt)
4350 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
4351 print *,'4-smelt i,j,smelt,rsm,snwepr,rsmfrac', &
4352 i,j,smelt,rsm,snwepr,rsmfrac
4355 !-- update liquid equivalent of snow depth
4356 !-- for evaporation and snow melt
4357 snwe = amax1(0.,(snwepr- &
4358 (smelt+beta*epot*ras)*delt &
4361 !--- if there is no snow melting then just evaporation
4362 !--- or condensation changes snwe
4364 if(snhei.ne.0.) then
4365 epot=-qkms*(qvatm-qsg)
4366 snwe = amax1(0.,(snwepr- &
4367 beta*epot*ras*delt))
4372 ! no iteration for snow on sea ice, because it will produce
4373 ! skin temperature higher than it is possible with snow on sea ice
4374 ! if(nmelt.eq.1) goto 212 ! second iteration
4377 if(smelt > 0..and. rsm > 0.) then
4378 if(snwe.le.rsm) then
4379 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
4380 print *,'seaice snwe<rsm snwe,rsm,smelt*delt,epot*ras*delt,beta', &
4381 snwe,rsm,smelt*delt,epot*ras*delt,beta
4384 !*** update snow density on effect of snow melt, melted
4385 !*** from the top of the snow. 13% of melted water
4386 !*** remains in the pack and changes its density.
4387 !*** eq. 9 (with my correction) in koren et al. (1999)
4389 xsn=(rhosn*(snwe-rsm)+1.e3*rsm)/ &
4391 rhosn=min(max(58.8,xsn),500.)
4394 thdifsn = 0.265/rhocsn
4400 !--- if vegfrac.ne.0. then some snow stays on the canopy
4401 !--- and should be added to snwe for water conservation
4402 ! 4 nov 07 +vegfrac*cst
4403 snheiprint=snweprint*1.e3 / rhosn
4405 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
4406 print *, 'snweprint : ',snweprint
4407 print *, 'd9sn,soilt,tsob : ', d9sn,soilt,tsob
4409 if(snhei.gt.0.) then
4411 tsnav=0.5/snhei*((soilt+soilt1)*deltsn &
4412 +(soilt1+tso(1))*(snhei-deltsn)) &
4415 tsnav=0.5*(soilt+tso(1)) - 273.15
4418 !--- recalculation of dew using new value of qsg
4421 qsg= qsn(soilt,tbq)/pp
4422 epot = -fq*(qvatm-qsg)
4428 snom=snom+smelt*delt*1.e3
4430 !--- the diagnostics of surface fluxes
4432 t3 = stbolt*tnold*tnold*tnold
4433 upflux = t3 *0.5*(soilt+tnold)
4434 xinet = emiss*(glw-upflux)
4435 ! rnet = gsw + xinet
4436 hft=-tkms*cp*rho*(tabs-soilt)
4437 hfx=-tkms*cp*rho*(tabs-soilt) &
4438 *(p1000mb*0.00001/patm)**rovcp
4439 q1 = - fq*ras* (qvatm - qsg)
4443 !-- moisture flux for coupling with myj pbl
4444 eeta=-qkms*ras*(qvatm/(1.+qvatm) - qsg/(1.+qsg))*1.e3
4446 !-- actual moisture flux from ruc lsm
4447 dew=qkms*(qvatm-qsg)
4456 !-- moisture flux for coupling with myj pbl
4457 eeta=-qkms*ras*beta*(qvatm/(1.+qvatm) - qvg/(1.+qvg))*1.e3
4459 ! to convert from m s-1 to kg m-2 s-1: *rho water=1.e3************
4460 !-- actual moisture flux from ruc lsm
4469 if(snhei.ge.snth)then
4470 s=thdifsn*rhocsn*(soilt-tsob)/snprim
4472 elseif(snhei.lt.snth.and.snhei.gt.0.) then
4473 s=(fsn*thdifsn*rhocsn+fso*thdifice(1)*rhcs)* &
4476 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
4477 print *,'snow is thin, snflx',i,j,snflx
4480 snflx=d9sn*(soilt-tsob)
4481 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
4482 print *,'snow is gone, snflx',i,j,snflx
4486 snhei=snwe *1.e3 / rhosn
4488 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
4489 print *,'snhei,snoh',i,j,snhei,snoh
4492 x= (r21+d9sn*r22sn)*(soilt-tnold) + &
4493 xlvm*r210*(qsg-qgold)
4494 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
4495 print *,'snowseaice storage ',i,j,x
4496 print *,'r21,d9sn,r22sn,soiltfrac,tnold,qsg,qgold,snprim', &
4497 r21,d9sn,r22sn,soiltfrac,tnold,qsg,qgold,snprim
4500 -rhonewcsn*newsnow/delt*(min(273.15,tabs)-soilt) &
4501 -rainf*cvw*prcpms*(max(273.15,tabs)-soilt)
4503 ! -- excess energy is spent on ice melt
4504 icemelt = rnet-hft-xlvm*eeta-s-snoh-x
4505 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
4506 print *,'snowseaice icemelt=',icemelt
4509 fltot=rnet-hft-xlvm*eeta-s-snoh-x-icemelt
4510 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
4511 print *,'i,j,snhei,qsg,soilt,soilt1,tso,tabs,qvatm', &
4512 i,j,snhei,qsg,soilt,soilt1,tso,tabs,qvatm
4513 print *,'snowseaice - fltot,rnet,hft,qfx,s,snoh,icemelt,snodif,x,soilt=' &
4514 ,fltot,rnet,hft,xlvm*eeta,s,snoh,icemelt,snodif,x,soilt
4516 !-- restore sea-ice parameters if snow is less than threshold
4517 if(snhei.eq.0.) then
4524 !------------------------------------------------------------------------
4525 !------------------------------------------------------------------------
4526 end subroutine snowseaice
4527 !------------------------------------------------------------------------
4530 subroutine soiltemp( &
4531 !--- input variables
4533 delt,ktau,conflx,nzs,nddzs,nroot, &
4534 prcpms,rainf,patm,tabs,qvatm,qcatm, &
4536 qkms,tkms,pc,rho,vegfrac,lai, &
4537 thdif,cap,drycan,wetcan, &
4538 transum,dew,mavail,soilres,alfa, &
4539 !--- soil fixed fields
4541 zsmain,zshalf,dtdzs,tbq, &
4543 xlv,cp,g0_p,cvw,stbolt, &
4544 !--- output variables
4545 tso,soilt,qvg,qsg,qcg,x)
4547 !*************************************************************
4548 ! energy budget equation and heat diffusion eqn are
4551 ! delt - time step (s)
4552 ! ktau - numver of time step
4553 ! conflx - depth of constant flux layer (m)
4554 ! ime, jme, kme, nzs - dimensions of the domain
4555 ! nroot - number of levels within the root zone
4556 ! prcpms - precipitation rate in m/s
4557 ! cotso, rhtso - coefficients for implicit solution of
4558 ! heat diffusion equation
4559 ! thdif - thermal diffusivity (m^2/s)
4560 ! qsg,qvg,qcg - saturated mixing ratio, mixing ratio of
4561 ! water vapor and cloud at the ground
4562 ! surface, respectively (kg/kg)
4563 ! patm - pressure [bar]
4564 ! qc3d,qv3d - cloud and water vapor mixing ratio
4565 ! at the first atm. level (kg/kg)
4566 ! emiss,rnet - emissivity (0-1) of the ground surface and net
4567 ! radiation at the surface (w/m^2)
4568 ! qkms - exchange coefficient for water vapor in the
4569 ! surface layer (m/s)
4570 ! tkms - exchange coefficient for heat in the surface
4572 ! pc - plant coefficient (resistance)
4573 ! rho - density of atmosphere near surface (kg/m^3)
4574 ! vegfrac - greeness fraction (0-1)
4575 ! cap - volumetric heat capacity (j/m^3/k)
4576 ! drycan - dry fraction of vegetated area where
4577 ! transpiration may take place (0-1)
4578 ! wetcan - fraction of vegetated area covered by canopy
4580 ! transum - transpiration function integrated over the
4582 ! dew - dew in kg/m^2s
4583 ! mavail - fraction of maximum soil moisture in the top
4585 ! zsmain - main levels in soil (m)
4586 ! zshalf - middle of the soil layers (m)
4587 ! dtdzs - dt/(2.*dzshalf*dzmain)
4588 ! tbq - table to define saturated mixing ration
4589 ! of water vapor for given temperature and pressure
4590 ! tso - soil temperature (k)
4591 ! soilt - skin temperature (k)
4593 !****************************************************************
4596 !-----------------------------------------------------------------
4598 !--- input variables
4600 integer, intent(in ) :: nroot,ktau,nzs , &
4601 nddzs !nddzs=2*(nzs-2)
4602 integer, intent(in ) :: i,j,iland,isoil
4603 real, intent(in ) :: delt,conflx,prcpms, rainf
4604 real, intent(inout) :: drycan,wetcan,transum
4605 !--- 3-d atmospheric variables
4607 intent(in ) :: patm, &
4623 !--- soil properties
4634 real, intent(in ) :: cp, &
4642 real, dimension(1:nzs), intent(in) :: zsmain, &
4647 real, dimension(1:nddzs), intent(in) :: dtdzs
4649 real, dimension(1:5001), intent(in) :: tbq
4652 !--- input/output variables
4653 !-------- 3-d soil moisture and temperature
4654 real, dimension( 1:nzs ) , &
4655 intent(inout) :: tso
4657 !-------- 2-d variables
4667 !--- local variables
4669 real :: x,x1,x2,x4,dzstop,can,ft,sph , &
4670 tn,trans,umveg,denom,fex
4672 real :: fkt,d1,d2,d9,d10,did,r211,r21,r22,r6,r7,d11 , &
4673 pi,h,fkq,r210,aa,bb,pp,q1,qs1,ts1,tq2,tx2 , &
4676 real :: c,cc,aa1,rhcs,h1, qgold
4678 real, dimension(1:nzs) :: cotso,rhtso
4680 integer :: nzs1,nzs2,k,k1,kn,kk, iter
4683 !-----------------------------------------------------------------
4689 dzstop=1./(zsmain(2)-zsmain(1))
4697 !******************************************************************************
4698 ! coefficients for thomas algorithm for tso
4699 !******************************************************************************
4700 ! did=2.*(zsmain(nzs)-zshalf(nzs))
4701 ! h1=dtdzs(8)*thdif(nzs-1)*(zshalf(nzs)-zshalf(nzs-1))/did
4702 ! cotso(1)=h1/(1.+h1)
4703 ! rhtso(1)=(tso(nzs)+h1*(tso(nzs-1)-tso(nzs)))/
4710 x1=dtdzs(k1)*thdif(kn-1)
4711 x2=dtdzs(k1+1)*thdif(kn)
4712 ft=tso(kn)+x1*(tso(kn-1)-tso(kn)) &
4713 -x2*(tso(kn)-tso(kn+1))
4714 denom=1.+x1+x2-x2*cotso(k)
4716 rhtso(k+1)=(ft+x2*rhtso(k))/denom
4719 !************************************************************************
4720 !--- the heat balance equation (Smirnova et al., 1996, eq. 21,26)
4726 trans=transum*drycan/zshalf(nroot+1)
4728 umveg=(1.-vegfrac) * soilres
4734 d9=thdif(1)*rhcs*dzstop
4738 r22=.5/(thdif(1)*delt*dzstop**2)
4739 r6=emiss *stbolt*.5*tn**4
4742 tdenom=d9*(1.-d1+r22)+d10+r21+r7 &
4748 aa=xlv*(fkq*umveg+r210)/tdenom
4749 bb=(d10*tabs+r21*tn+xlv*(qvatm* &
4751 +r210*qvg)+d11+d9*(d2+r22*tn) &
4752 +rainf*cvw*prcpms*max(273.15,tabs) &
4758 call vilka(tn,aa1,bb,pp,qs1,ts1,tbq,ktau,i,j,iland,isoil)
4763 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
4764 print *,'vilka1 - ts1,qs1,tq2,h,tx2,q1',ts1,qs1,tq2,h,tx2,q1
4766 !with alfa go to 100
4767 if(q1.lt.qs1) goto 100
4768 !--- if no saturation - goto 100
4769 !--- if saturation - goto 90
4774 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
4775 print *,'90 qvg,qsg,qcg,tso(1)',qvg,qsg,qcg,tso(1)
4780 call vilka(tn,aa,bb,pp,qs1,ts1,tbq,ktau,i,j,iland,isoil)
4782 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
4783 print *,'vilka2 - ts1,qs1,tq2,h,tx2,q1',ts1,qs1,tq2,h,tx2,q1
4785 if(q1.ge.qs1) goto 90
4786 !with alfa 100 continue
4789 ! if( qs1>qvatm .and. qvatm > qvg) then
4791 ! print *,'very dry soils mavail,qvg,qs1,qvatm,ts1',i,j,mavail,qvg,qs1,qvatm,ts1
4796 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
4797 print *,'q1,qsg,qvg,qvatm,alfa,h',q1,qsg,qvg,qvatm,alfa,h
4800 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
4801 print *,'200 qvg,qsg,qcg,tso(1)',qvg,qsg,qcg,tso(1)
4804 !--- soilt - skin temperature
4807 !---- final solution for soil temperature - tso
4810 tso(k)=rhtso(kk)+cotso(kk)*tso(k-1)
4813 x= (cp*rho*r211+rhcs*zsmain(2)*0.5/delt)*(soilt-tn) + &
4814 xlv*rho*r211*(qvg-qgold)
4816 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
4817 print*,'soiltemp storage, i,j,x,soilt,tn,qvg,qvgold', &
4818 i,j,x,soilt,tn,qvg,qgold
4819 print *,'temp term (cp*rho*r211+rhcs*zsmain(2)*0.5/delt)*(soilt-tn)',&
4820 (cp*rho*r211+rhcs*zsmain(2)*0.5/delt)*(soilt-tn)
4821 print *,'qv term xlv*rho*r211*(qvg-qgold)',xlv*rho*r211*(qvg-qgold)
4825 -rainf*cvw*prcpms*(max(273.15,tabs)-soilt)
4827 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
4831 !--------------------------------------------------------------------
4832 end subroutine soiltemp
4833 !--------------------------------------------------------------------
4836 subroutine snowtemp( &
4837 !--- input variables
4839 delt,ktau,conflx,nzs,nddzs,nroot, &
4840 snwe,snwepr,snhei,newsnow,snowfrac, &
4841 beta,deltsn,snth,rhosn,rhonewsn,meltfactor, & ! add meltfactor
4843 patm,tabs,qvatm,qcatm, &
4844 glw,gsw,emiss,rnet, &
4845 qkms,tkms,pc,rho,vegfrac, &
4846 thdif,cap,drycan,wetcan,cst, &
4847 tranf,transum,dew,mavail, &
4848 !--- soil fixed fields
4849 dqm,qmin,psis,bclh, &
4850 zsmain,zshalf,dtdzs,tbq, &
4852 xlvm,cp,rovcp,g0_p,cvw,stbolt, &
4853 !--- output variables
4854 snweprint,snheiprint,rsm, &
4855 tso,soilt,soilt1,tsnav,qvg,qsg,qcg, &
4856 smelt,snoh,snflx,s,ilnb,x)
4858 !********************************************************************
4859 ! energy budget equation and heat diffusion eqn are
4860 ! solved here to obtain snow and soil temperatures
4862 ! delt - time step (s)
4863 ! ktau - numver of time step
4864 ! conflx - depth of constant flux layer (m)
4865 ! ime, jme, kme, nzs - dimensions of the domain
4866 ! nroot - number of levels within the root zone
4867 ! prcpms - precipitation rate in m/s
4868 ! cotso, rhtso - coefficients for implicit solution of
4869 ! heat diffusion equation
4870 ! thdif - thermal diffusivity (w/m/k)
4871 ! qsg,qvg,qcg - saturated mixing ratio, mixing ratio of
4872 ! water vapor and cloud at the ground
4873 ! surface, respectively (kg/kg)
4874 ! patm - pressure [bar]
4875 ! qcatm,qvatm - cloud and water vapor mixing ratio
4876 ! at the first atm. level (kg/kg)
4877 ! emiss,rnet - emissivity (0-1) of the ground surface and net
4878 ! radiation at the surface (w/m^2)
4879 ! qkms - exchange coefficient for water vapor in the
4880 ! surface layer (m/s)
4881 ! tkms - exchange coefficient for heat in the surface
4883 ! pc - plant coefficient (resistance)
4884 ! rho - density of atmosphere near surface (kg/m^3)
4885 ! vegfrac - greeness fraction (0-1)
4886 ! cap - volumetric heat capacity (j/m^3/k)
4887 ! drycan - dry fraction of vegetated area where
4888 ! transpiration may take place (0-1)
4889 ! wetcan - fraction of vegetated area covered by canopy
4891 ! transum - transpiration function integrated over the
4893 ! dew - dew in kg/m^2/s
4894 ! mavail - fraction of maximum soil moisture in the top
4896 ! zsmain - main levels in soil (m)
4897 ! zshalf - middle of the soil layers (m)
4898 ! dtdzs - dt/(2.*dzshalf*dzmain)
4899 ! tbq - table to define saturated mixing ration
4900 ! of water vapor for given temperature and pressure
4901 ! tso - soil temperature (k)
4902 ! soilt - skin temperature (k)
4904 !*********************************************************************
4907 !---------------------------------------------------------------------
4908 !--- input variables
4910 integer, intent(in ) :: nroot,ktau,nzs , &
4911 nddzs !nddzs=2*(nzs-2)
4913 integer, intent(in ) :: i,j,iland,isoil
4914 real, intent(in ) :: delt,conflx,prcpms , &
4915 rainf,newsnow,deltsn,snth , &
4916 tabs,transum,snwepr , &
4920 !--- 3-d atmospheric variables
4922 intent(in ) :: patm, &
4927 intent(in ) :: glw, &
4935 !--- soil properties
4943 real, intent(in ) :: cp, &
4951 real, dimension(1:nzs), intent(in) :: zsmain, &
4957 real, dimension(1:nddzs), intent(in) :: dtdzs
4959 real, dimension(1:5001), intent(in) :: tbq
4962 !--- input/output variables
4963 !-------- 3-d soil moisture and temperature
4964 real, dimension( 1:nzs ) , &
4965 intent(inout) :: tso
4968 !-------- 2-d variables
4970 intent(inout) :: dew, &
4989 real, intent(inout) :: drycan, wetcan
4991 real, intent(out) :: rsm, &
4994 integer, intent(out) :: ilnb
4995 !--- local variables
4998 integer :: nzs1,nzs2,k,k1,kn,kk
5000 real :: x,x1,x2,x4,dzstop,can,ft,sph, &
5001 tn,trans,umveg,denom
5003 real :: cotsn,rhtsn,xsn1,ddzsn1,x1sn1,ftsnow,denomsn
5005 real :: t3,upflux,xinet,ras, &
5006 xlmelt,rhocsn,thdifsn, &
5007 beta,epot,xsn,ddzsn,x1sn,d1sn,d2sn,d9sn,r22sn
5010 fkt,d1,d2,d9,d10,did,r211,r21,r22,r6,r7,d11, &
5011 pi,h,fkq,r210,aa,bb,pp,q1,qs1,ts1,tq2,tx2, &
5012 tdenom,c,cc,aa1,rhcs,h1, &
5013 tsob, snprim, sh1, sh2, &
5014 smeltg,snohg,snodif,soh, &
5015 cmc2ms,tnold,qgold,snohgnew
5017 real, dimension(1:nzs) :: transp,cotso,rhtso
5025 real :: rnet,rsmfrac,soiltfrac,hsn,rr,keff,fact
5026 integer :: nmelt, iter
5028 !-----------------------------------------------------------------
5037 !-- options for snow conductivity:
5039 !-- opt 2 - Sturm et al., 1997
5042 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
5043 print *, 'snowtemp: snhei,snth,soilt1: ',snhei,snth,soilt1,soilt
5047 !18apr08 - add rhonewcsn
5048 rhonewcsn=2090.* rhonewsn
5050 if(isncond_opt == 1) then
5051 !-- old version thdifsn = 0.265/rhocsn
5052 thdifsn = 0.265/rhocsn
5054 !-- 07jun19 - thermal conductivity (k_eff) from Sturm et al.(1997)
5055 !-- keff = 10. ** (2.650 * rhosn*1.e-3 - 1.652)
5057 if(rhosn < 156. .or. (newsnow > 0. .and. rhonewsn < 156.)) then
5058 keff = 0.023 + 0.234 * rhosn * 1.e-3
5060 keff = 0.138 - 1.01 * rhosn*1.e-3 + 3.233 * rhosn**2 * 1.e-6
5062 if(newsnow <= 0. .and. snhei > 1. .and. rhosn > 250.) then
5063 !-- some areas with large snow depth have unrealistically
5064 !-- low snow density (in the rockie's with snow depth > 1 m).
5065 !-- based on Sturm et al. keff=0.452 typical for hard snow slabs
5066 !-- with rhosn=488 kg/m^3. thdifsn = 0.452/(2090*488)=4.431718e-7
5067 !-- in future a better compaction scheme is needed for these areas.
5068 thdifsn = 4.431718e-7
5070 thdifsn = keff/rhocsn * fact
5092 dzstop=1./(zsmain(2)-zsmain(1))
5094 !******************************************************************************
5095 ! coefficients for thomas algorithm for tso
5096 !******************************************************************************
5097 ! did=2.*(zsmain(nzs)-zshalf(nzs))
5098 ! h1=dtdzs(8)*thdif(nzs-1)*(zshalf(nzs)-zshalf(nzs-1))/did
5099 ! cotso(1)=h1/(1.+h1)
5100 ! rhtso(1)=(tso(nzs)+h1*(tso(nzs-1)-tso(nzs)))/
5108 x1=dtdzs(k1)*thdif(kn-1)
5109 x2=dtdzs(k1+1)*thdif(kn)
5110 ft=tso(kn)+x1*(tso(kn-1)-tso(kn)) &
5111 -x2*(tso(kn)-tso(kn+1))
5112 denom=1.+x1+x2-x2*cotso(k)
5114 rhtso(k+1)=(ft+x2*rhtso(k))/denom
5116 !--- the nzs element in cotso and rhtso will be for snow
5117 !--- there will be 2 layers in snow if it is deeper than deltsn+snth
5118 if(snhei.ge.snth) then
5119 if(snhei.le.deltsn+snth) then
5120 !-- 1-layer snow model
5121 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
5122 print *,'1-layer - snth,snhei,deltsn',snth,snhei,deltsn
5125 snprim=max(snth,snhei)
5128 xsn = delt/2./(zshalf(2)+0.5*snprim)
5129 ddzsn = xsn / snprim
5130 x1sn = ddzsn * thdifsn
5131 x2 = dtdzs(1)*thdif(1)
5132 ft = tso(1)+x1sn*(soilt-tso(1)) &
5134 denom = 1. + x1sn + x2 -x2*cotso(nzs1)
5135 cotso(nzs)=x1sn/denom
5136 rhtso(nzs)=(ft+x2*rhtso(nzs1))/denom
5139 !*** average temperature of snow pack (c)
5140 tsnav=0.5*(soilt+tso(1)) &
5144 !-- 2 layers in snow, soilt1 is temperasture at deltsn depth
5145 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
5146 print *,'2-layer - snth,snhei,deltsn',snth,snhei,deltsn
5151 xsn = delt/2./(0.5*deltsn)
5152 xsn1= delt/2./(zshalf(2)+0.5*(snhei-deltsn))
5153 ddzsn = xsn / deltsn
5154 ddzsn1 = xsn1 / (snhei-deltsn)
5155 x1sn = ddzsn * thdifsn
5156 x1sn1 = ddzsn1 * thdifsn
5157 x2 = dtdzs(1)*thdif(1)
5158 ft = tso(1)+x1sn1*(soilt1-tso(1)) &
5160 denom = 1. + x1sn1 + x2 - x2*cotso(nzs1)
5161 cotso(nzs)=x1sn1/denom
5162 rhtso(nzs)=(ft+x2*rhtso(nzs1))/denom
5163 ftsnow = soilt1+x1sn*(soilt-soilt1) &
5164 -x1sn1*(soilt1-tso(1))
5165 denomsn = 1. + x1sn + x1sn1 - x1sn1*cotso(nzs)
5167 rhtsn=(ftsnow+x1sn1*rhtso(nzs))/denomsn
5168 !*** average temperature of snow pack (c)
5169 tsnav=0.5/snhei*((soilt+soilt1)*deltsn &
5170 +(soilt1+tso(1))*(snhei-deltsn)) &
5174 if(snhei.lt.snth.and.snhei.gt.0.) then
5175 !--- snow is too thin to be treated separately, therefore it
5176 !--- is combined with the first soil layer.
5177 snprim=snhei+zsmain(2)
5182 xsn = delt/2./((zshalf(3)-zsmain(2))+0.5*snprim)
5184 x1sn = ddzsn * (fsn*thdifsn+fso*thdif(1))
5185 x2=dtdzs(2)*thdif(2)
5186 ft=tso(2)+x1sn*(soilt-tso(2))- &
5188 denom = 1. + x1sn + x2 - x2*cotso(nzs-2)
5189 cotso(nzs1) = x1sn/denom
5190 rhtso(nzs1)=(ft+x2*rhtso(nzs-2))/denom
5191 tsnav=0.5*(soilt+tso(1)) &
5193 cotso(nzs)=cotso(nzs1)
5194 rhtso(nzs)=rhtso(nzs1)
5200 !************************************************************************
5201 !--- the heat balance equation (Smirnova et al. 1996, eq. 21,26)
5202 !18apr08 nmelt is the flag for melting, and snoh is heat of snow phase changes
5207 epot=-qkms*(qvatm-qgold)
5210 trans=transum*drycan/zshalf(nroot+1)
5217 d9=thdif(1)*rhcs*dzstop
5221 r22=.5/(thdif(1)*delt*dzstop**2)
5222 r6=emiss *stbolt*.5*tn**4
5226 if(snhei.ge.snth) then
5227 if(snhei.le.deltsn+snth) then
5231 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
5232 print *,'1 layer d1sn,d2sn',i,j,d1sn,d2sn
5238 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
5239 print *,'2 layers d1sn,d2sn',i,j,d1sn,d2sn
5242 d9sn= thdifsn*rhocsn / snprim
5243 r22sn = snprim*snprim*0.5/(thdifsn*delt)
5244 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
5245 print *,'1 or 2 layers d9sn,r22sn',d9sn,r22sn
5249 if(snhei.lt.snth.and.snhei.gt.0.) then
5250 !--- thin snow is combined with soil
5253 d9sn = (fsn*thdifsn*rhocsn+fso*thdif(1)*rhcs)/ &
5255 r22sn = snprim*snprim*0.5 &
5256 /((fsn*thdifsn+fso*thdif(1))*delt)
5257 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
5258 print *,' combined d9sn,r22sn,d1sn,d2sn: ',d9sn,r22sn,d1sn,d2sn
5262 !--- all snow is sublimated
5267 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
5268 print *,' snhei = 0, d9sn,r22sn,d1sn,d2sn: ',d9sn,r22sn,d1sn,d2sn
5274 !18apr08 - the snow melt iteration start point
5277 !---- tdenom for snow
5278 tdenom = d9sn*(1.-d1sn +r22sn)+d10+r21+r7 &
5280 +rhonewcsn*newsnow/delt
5286 aa=xlvm*(beta*fkq*umveg+r210)/tdenom
5287 bb=(d10*tabs+r21*tn+xlvm*(qvatm* &
5288 (beta*fkq*umveg+c) &
5289 +r210*qgold)+d11+d9sn*(d2sn+r22sn*tn) &
5290 +rainf*cvw*prcpms*max(273.15,tabs) &
5291 + rhonewcsn*newsnow/delt*min(273.15,tabs) &
5298 call vilka(tn,aa1,bb,pp,qs1,ts1,tbq,ktau,i,j,iland,isoil)
5302 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
5303 print *,'vilka1 - ts1,qs1,tq2,h,tx2,q1',ts1,qs1,tq2,h,tx2,q1
5305 if(q1.lt.qs1) goto 100
5306 !--- if no saturation - goto 100
5307 !--- if saturation - goto 90
5311 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
5312 print *,'90 qvg,qsg,qcg,tso(1)',qvg,qsg,qcg,tso(1)
5317 call vilka(tn,aa,bb,pp,qs1,ts1,tbq,ktau,i,j,iland,isoil)
5319 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
5320 print *,'vilka2 - ts1,qs1,h,tx2,q1',ts1,qs1,tq2,h,tx2,q1
5322 if(q1.gt.qs1) goto 90
5326 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
5327 print *,'no saturation qvg,qsg,qcg,tso(1)',qvg,qsg,qcg,tso(1)
5331 !--- soilt - skin temperature
5334 if(nmelt==1 .and. snowfrac==1. .and. snwe > 0. .and. soilt > 273.15) then
5335 !--7feb22 on the second iteration when snoh is known and snwe > 0. after melting,
5336 !-- check if the snow skin temperature is =<tfrzk
5337 !-- when a grid cell is fully covered with snow (snowfrac=1)
5338 !-- or with partial snow cover and snow_mosaic=1 (snowfrac=1).
5339 soilt = min(273.15,soilt)
5342 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
5343 if(i.eq.266.and.j.eq.447) then
5344 print *,'snwe,snhei,soilt,soilt1,tso',i,j,snwe,snhei,soilt,soilt1,tso
5347 ! solution for temperature at 7.5 cm depth and snow-soil interface
5348 if(snhei.ge.snth) then
5349 if(snhei.gt.deltsn+snth) then
5350 !-- 2-layer snow model
5351 soilt1=min(273.15,rhtsn+cotsn*soilt)
5352 tso(1)=rhtso(nzs)+cotso(nzs)*soilt1
5356 tso(1)=rhtso(nzs)+cotso(nzs)*soilt
5360 elseif (snhei > 0. .and. snhei < snth) then
5362 tso(2)=rhtso(nzs1)+cotso(nzs1)*soilt
5363 tso(1)=(tso(2)+(soilt-tso(2))*fso)
5367 !-- very thin or zero snow. if snow is thin we suppose that
5368 !--- tso(i,j,1)=soilt, and later we recompute tso(i,j,1)
5373 if(nmelt==1.and.snowfrac==1.) then
5374 !-- second iteration with full snow cover
5375 soilt1= min(273.15,soilt1)
5376 tso(1)= min(273.15,tso(1))
5377 tsob = min(273.15,tsob)
5381 !---- final solution for tso
5382 if (snhei > 0. .and. snhei < snth) then
5383 ! blended or snow is melted
5386 tso(k)=rhtso(kk)+cotso(kk)*tso(k-1)
5392 tso(k)=rhtso(kk)+cotso(kk)*tso(k-1)
5395 !--- for thin snow layer combined with the top soil layer
5396 !--- tso(1) is recomputed by linear interpolation between soilt
5398 ! if(snhei.lt.snth.and.snhei.gt.0.)then
5399 ! tso(1)=tso(2)+(soilt-tso(2))*fso
5405 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
5406 print *,'soilt,soilt1,tso,tsob,qsg',i,j,soilt,soilt1,tso,tsob,qsg,'nmelt=',nmelt
5409 if(nmelt.eq.1) go to 220
5411 !--- if soilt > 273.15 f then melting of snow can happen
5412 ! if(soilt.gt.273.15.and.snhei.gt.0.) then
5413 ! if all snow can evaporate, then there is nothing to melt
5414 if(soilt.gt.273.15.and.beta==1..and.snhei.gt.0.) then
5416 soiltfrac=snowfrac*273.15+(1.-snowfrac)*soilt
5417 qsg=min(qsg, qsn(soiltfrac,tbq)/pp)
5418 qvg=snowfrac*qsg+(1.-snowfrac)*qvg
5419 t3 = stbolt*tn*tn*tn
5420 upflux = t3 * 0.5*(tn + soiltfrac)
5421 xinet = emiss*(glw-upflux)
5422 epot = -qkms*(qvatm-qsg)
5425 if (q1.le.0..or.iter==1) then
5437 transp(k)=-vegfrac*q1 &
5438 *tranf(k)*drycan/zshalf(nroot+1)
5445 edir1 = q1*umveg * beta
5446 ec1 = q1 * wetcan * vegfrac
5448 eeta = (edir1 + ec1 + ett1)*1.e3
5449 ! to convert from kg m-2 s-1 to m s-1: 1/rho water=1.e-3************
5453 hfx=-d10*(tabs-soiltfrac)
5455 if(snhei.ge.snth)then
5456 soh=thdifsn*rhocsn*(soiltfrac-tsob)/snprim
5459 soh=(fsn*thdifsn*rhocsn+fso*thdif(1)*rhcs)* &
5460 (soiltfrac-tsob)/snprim
5465 x= (r21+d9sn*r22sn)*(soiltfrac-tn) + &
5466 xlvm*r210*(qvg-qgold)
5467 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
5468 print *,'snowtemp storage ',i,j,x
5469 print *,'r21,d9sn,r22sn,soiltfrac,tn,qsg,qvg,qgold,snprim', &
5470 r21,d9sn,r22sn,soiltfrac,tn,qsg,qvg,qgold,snprim
5473 !-- snoh is energy flux of snow phase change
5474 snoh=rnet-qfx -hfx - soh - x &
5475 +rhonewcsn*newsnow/delt*(min(273.15,tabs)-soiltfrac) &
5476 +rainf*cvw*prcpms*(max(273.15,tabs)-soiltfrac)
5478 !-- smelt is speed of melting in m/s
5479 smelt= snoh /xlmelt*1.e-3
5480 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
5481 print *,'1- smelt',i,j,smelt
5483 if(epot.gt.0. .and. snwepr.le.epot*ras*delt) then
5484 !-- all snow can evaporate
5485 beta=snwepr/(epot*ras*delt)
5486 smelt=amin1(smelt,snwepr/delt-beta*epot*ras)
5488 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
5489 print *,'2- smelt',i,j,smelt
5494 smelt=amax1(0.,smelt)
5496 !18apr08 - egglston limit
5497 !-- 22apr22 do not limit snow melting for hail (rhonewsn > 450), or dense snow
5498 !-- (rhosn > 350.) with very warm surface temperatures (>10c)
5499 if( (rhosn < 350. .or. (newsnow > 0. .and. rhonewsn < 450.)) .and. soilt < 283. ) then
5500 ! smelt= amin1 (smelt, 5.6e-7*meltfactor*max(1.,(soilt-273.15)))
5501 smelt= amin1 (smelt, delt/60.*5.6e-8*meltfactor*max(1.,(soilt-273.15)))
5503 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
5504 print *,'3- smelt',i,j,smelt
5508 ! rr - potential melting
5509 rr=max(0.,snwepr/delt-beta*epot*ras)
5513 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
5514 print *,'4- smelt i,j,smelt,rr',i,j,smelt,rr
5519 snohgnew=smelt*xlmelt*1.e3
5520 snodif=amax1(0.,(snoh-snohgnew))
5523 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
5524 print *,'snoh,snodif',snoh,snodif
5527 if( smelt > 0.) then
5528 !*** from koren et al. (1999) 13% of snow melt stays in the snow pack
5529 rsmfrac=min(0.18,(max(0.08,snwepr/0.10*0.13)))
5530 if(snhei > 0.01 .and. rhosn < 350.) then
5531 rsm=rsmfrac*smelt*delt
5533 ! do not keep melted water if snow depth is less that 1 cm
5536 !18apr08 rsm is part of melted water that stays in snow as liquid
5538 smelt=max(0.,smelt-rsm/delt)
5539 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
5540 print *,'5- smelt i,j,smelt,rsm,snwepr,rsmfrac', &
5541 i,j,smelt,rsm,snwepr,rsmfrac
5547 !-- update of liquid equivalent of snow depth
5548 !-- due to evaporation and snow melt
5550 snwe = amax1(0.,(snwepr- &
5551 (smelt+beta*epot*ras)*delt &
5555 !--- if there is no snow melting then just evaporation
5556 !--- or condensation cxhanges snwe
5558 if(snhei.ne.0..and. beta == 1.) then
5559 epot=-qkms*(qvatm-qsg)
5560 snwe = amax1(0.,(snwepr- &
5561 beta*epot*ras*delt))
5563 !-- all snow is sublibated
5568 !18apr08 - if snow melt occurred then go into iteration for energy budget solution
5569 if(nmelt.eq.1) goto 212 ! second interation
5572 if(smelt.gt.0..and.rsm.gt.0.) then
5573 if(snwe.le.rsm) then
5575 print *,'snwe<rsm snwe,rsm,smelt*delt,epot*ras*delt,beta', &
5576 snwe,rsm,smelt*delt,epot*ras*delt,beta
5579 !*** update snow density on effect of snow melt, melted
5580 !*** from the top of the snow. 13% of melted water
5581 !*** remains in the pack and changes its density.
5582 !*** eq. 9 (with my correction) in koren et al. (1999)
5583 xsn=(rhosn*(snwe-rsm)+1.e3*rsm)/ &
5585 rhosn=min(max(58.8,xsn),500.)
5588 if(isncond_opt == 1) then
5589 !-- old version thdifsn = 0.265/rhocsn
5590 thdifsn = 0.265/rhocsn
5592 !-- 07jun19 - thermal conductivity (k_eff) from Sturm et al.(1997)
5593 !-- keff = 10. ** (2.650 * rhosn*1.e-3 - 1.652)
5595 if(rhosn < 156. .or. (newsnow > 0. .and. rhonewsn < 156.)) then
5596 keff = 0.023 + 0.234 * rhosn * 1.e-3
5598 keff = 0.138 - 1.01 * rhosn*1.e-3 + 3.233 * rhosn**2 * 1.e-6
5600 if(newsnow <= 0. .and. snhei > 1. .and. rhosn > 250.) then
5601 !-- some areas with large snow depth have unrealistically
5602 !-- low snow density (in the rockie's with snow depth > 1 m).
5603 !-- based on Sturm et al. keff=0.452 typical for hard snow slabs
5604 !-- with rhosn=488 kg/m^3. thdifsn = 0.452/(2090*488)=4.431718e-7
5605 !-- in future a better compaction scheme is needed for these areas.
5606 thdifsn = 4.431718e-7
5608 thdifsn = keff/rhocsn * fact
5615 !--- compute flux in the top snow layer
5616 if(snhei.ge.snth)then
5617 s=thdifsn*rhocsn*(soilt-tsob)/snprim
5619 s=d9*(tso(1)-tso(2))
5620 elseif(snhei.lt.snth.and.snhei.gt.0.) then
5621 s=(fsn*thdifsn*rhocsn+fso*thdif(1)*rhcs)* &
5624 s=d9*(tso(1)-tso(2))
5628 s=d9*(tso(1)-tso(2))
5631 snhei=snwe *1.e3 / rhosn
5632 !-- if ground surface temperature
5633 !-- is above freezing snow can melt from the bottom. the following
5634 !-- piece of code will check if bottom melting is possible.
5636 if(tso(1).gt.273.15 .and. snhei > 0.) then
5637 if (snhei.gt.deltsn+snth) then
5638 hsn = snhei - deltsn
5639 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
5640 print*,'2 layer snow - snhei,hsn',snhei,hsn
5643 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
5644 print*,'1 layer snow or blended - snhei',snhei
5649 soiltfrac=snowfrac*273.15+(1.-snowfrac)*tso(1)
5651 snohg=(tso(1)-soiltfrac)*(cap(1)*zshalf(2)+ &
5652 rhocsn*0.5*hsn) / delt
5653 snohg=amax1(0.,snohg)
5655 smeltg=snohg/xlmelt*1.e-3
5656 ! egglston - empirical limit on snow melt from the bottom of snow pack
5657 !9jun22-- the next line excludeis cases of summer hail from snowmelt limiting
5658 if( (rhosn < 350. .or. (newsnow > 0. .and. rhonewsn < 450.)) .and. soilt < 283. ) then
5659 smeltg=amin1(smeltg, 5.8e-9)
5662 ! rr - potential melting
5664 smeltg=amin1(smeltg, rr)
5666 snohgnew=smeltg*xlmelt*1.e3
5667 snodif=amax1(0.,(snohg-snohgnew))
5668 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
5669 print *,'tso(1),soiltfrac,smeltg,snodif',tso(1),soiltfrac,smeltg,snodif
5672 snwe=max(0.,snwe-smeltg*delt)
5673 snhei=snwe *1.e3 / rhosn
5674 !-- add up all snow melt
5675 smelt = smelt + smeltg
5677 if(snhei > 0.) tso(1) = soiltfrac
5678 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
5679 print *,'melt from the bottom snwe,snhei',snwe,snhei
5681 print *,'snow is all melted on the warm ground'
5685 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
5686 print *,'snhei,snoh',i,j,snhei,snoh
5690 snheiprint=snweprint*1.e3 / rhosn
5692 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
5693 print *, 'snweprint : ',snweprint
5694 print *, 'd9sn,soilt,tsob : ', d9sn,soilt,tsob
5697 x= (r21+d9sn*r22sn)*(soilt-tn) + &
5698 xlvm*r210*(qsg-qgold)
5699 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
5700 print *,'snowtemp storage ',i,j,x
5701 print *,'r21,d9sn,r22sn,soiltfrac,soilt,tn,qsg,qgold,snprim', &
5702 r21,d9sn,r22sn,soiltfrac,soilt,tn,qsg,qgold,snprim
5706 ! "heat" from snow and rain
5707 -rhonewcsn*newsnow/delt*(min(273.15,tabs)-soilt) &
5708 -rainf*cvw*prcpms*(max(273.15,tabs)-soilt)
5709 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
5711 print *,'snhei=',snhei
5712 print *,'snflx=',snflx
5715 if(snhei.gt.0.) then
5717 tsnav=0.5/snhei*((soilt+soilt1)*deltsn &
5718 +(soilt1+tso(1))*(snhei-deltsn)) &
5721 tsnav=0.5*(soilt+tso(1)) - 273.15
5724 tsnav= soilt - 273.15
5727 !------------------------------------------------------------------------
5728 end subroutine snowtemp
5729 !------------------------------------------------------------------------
5732 subroutine soilmoist ( &
5734 delt,nzs,nddzs,dtdzs,dtdzs2,riw, &
5735 zsmain,zshalf,diffu,hydro, &
5736 qsg,qvg,qcg,qcatm,qvatm,prcp, &
5738 dew,smelt,soilice,vegfrac,snowfrac,soilres, &
5740 dqm,qmin,ref,ksat,ras,infmax, &
5742 soilmois,soiliqw,mavail,runoff,runoff2,infiltrp)
5743 !*************************************************************************
5744 ! moisture balance equation and richards eqn.
5747 ! delt - time step (s)
5748 ! ime,jme,nzs - dimensions of soil domain
5749 ! zsmain - main levels in soil (m)
5750 ! zshalf - middle of the soil layers (m)
5751 ! dtdzs - dt/(2.*dzshalf*dzmain)
5752 ! dtdzs2 - dt/(2.*dzshalf)
5753 ! diffu - diffusional conductivity (m^2/s)
5754 ! hydro - hydraulic conductivity (m/s)
5755 ! qsg,qvg,qcg - saturated mixing ratio, mixing ratio of
5756 ! water vapor and cloud at the ground
5757 ! surface, respectively (kg/kg)
5758 ! qcatm,qvatm - cloud and water vapor mixing ratio
5759 ! at the first atm. level (kg/kg)
5760 ! prcp - precipitation rate in m/s
5761 ! qkms - exchange coefficient for water vapor in the
5762 ! surface layer (m/s)
5763 ! transp - transpiration from the soil layers (m/s)
5764 ! drip - liquid water dripping from the canopy to soil (m)
5765 ! dew - dew in kg/m^2s
5766 ! smelt - melting rate in m/s
5767 ! soilice - volumetric content of ice in soil (m^3/m^3)
5768 ! soiliqw - volumetric content of liquid water in soil (m^3/m^3)
5769 ! vegfrac - greeness fraction (0-1)
5770 ! ras - ration of air density to soil density
5771 ! infmax - maximum infiltration rate (kg/m^2/s)
5773 ! soilmois - volumetric soil moisture, 6 levels (m^3/m^3)
5774 ! mavail - fraction of maximum soil moisture in the top
5776 ! runoff - surface runoff (m/s)
5777 ! runoff2 - underground runoff (m)
5778 ! infiltrp - point infiltration flux into soil (m/s)
5779 ! /(snow bottom runoff) (mm/s)
5781 ! cosmc, rhsmc - coefficients for implicit solution of
5783 !******************************************************************
5785 !------------------------------------------------------------------
5786 !--- input variables
5787 real, intent(in ) :: delt
5788 integer, intent(in ) :: nzs,nddzs
5792 real, dimension(1:nzs), intent(in ) :: zsmain, &
5800 real, dimension(1:nddzs), intent(in) :: dtdzs
5802 real, intent(in ) :: qsg,qvg,qcg,qcatm,qvatm , &
5803 qkms,vegfrac,drip,prcp , &
5804 dew,smelt,snowfrac , &
5805 dqm,qmin,ref,ksat,ras,riw,soilres
5809 real, dimension( 1:nzs ) , &
5811 intent(inout) :: soilmois,soiliqw
5813 real, intent(inout) :: mavail,runoff,runoff2,infiltrp, &
5818 real, dimension( 1:nzs ) :: cosmc,rhsmc
5820 real :: dzs,r1,r2,r3,r4,r5,r6,r7,r8,r9,r10
5821 real :: refkdt,refdk,delt1,f1max,f2max
5822 real :: f1,f2,fd,kdt,val,ddt,px,fk,fkmax
5823 real :: qq,umveg,infmax1,trans
5824 real :: totliq,flx,flxsat,qtot
5825 real :: did,x1,x2,x4,denom,q2,q4
5826 real :: dice,fcr,acrt,frzx,sum,cvfrz
5828 integer :: nzs1,nzs2,k,kk,k1,kn,ialp1,jj,jk
5830 !******************************************************************************
5831 ! coefficients for thomas algorithm for soilmois
5832 !******************************************************************************
5836 118 format(6(10pf23.19))
5843 did=(zsmain(nzs)-zshalf(nzs))
5844 x1=zsmain(nzs)-zsmain(nzs1)
5846 !7may09 did=(zsmain(nzs)-zshalf(nzs))*2.
5847 ! denom=did/delt+diffu(nzs1)/x1
5848 ! cosmc(1)=diffu(nzs1)/x1/denom
5849 ! rhsmc(1)=(soilmois(nzs)*did/delt
5850 ! 1 +transp(nzs)-(hydro(nzs)*soilmois(nzs)
5851 ! 1 -hydro(nzs1)*soilmois(nzs1))*did
5854 denom=(1.+diffu(nzs1)/x1/did*delt+hydro(nzs)/(2.*did)*delt)
5855 cosmc(1)=delt*(diffu(nzs1)/did/x1 &
5856 +hydro(nzs1)/2./did)/denom
5857 rhsmc(1)=(soilmois(nzs)+transp(nzs)*delt/ &
5860 ! rhsmc(1)=(soilmois(nzs)*did/delt &
5861 ! +transp(nzs)-(hydro(nzs)*soilmois(nzs) &
5862 ! -hydro(nzs1)*soilmois(nzs1))*did &
5865 !12 june 2014 - low boundary condition: 1 - zero diffusion below the lowest
5866 ! level; 2 - soil moisture at the low boundary can be lost due to the root uptake.
5867 ! so far - no interaction with the water table.
5869 denom=1.+diffu(nzs1)/x1/did*delt
5870 cosmc(1)=delt*(diffu(nzs1)/did/x1 &
5871 +hydro(nzs1)/did)/denom
5872 rhsmc(1)=(soilmois(nzs)-hydro(nzs)*delt/did*soilmois(nzs) &
5873 +transp(nzs)*delt/did)/denom
5875 rhsmc(1)=soilmois(nzs)
5880 x4=2.*dtdzs(k1)*diffu(kn-1)
5881 x2=2.*dtdzs(k1+1)*diffu(kn)
5882 q4=x4+hydro(kn-1)*dtdzs2(kn-1)
5883 q2=x2-hydro(kn+1)*dtdzs2(kn-1)
5884 denom=1.+x2+x4-q2*cosmc(k)
5886 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
5887 print *,'q2,soilmois(kn),diffu(kn),x2,hydro(kn+1),dtdzs2(kn-1),kn,k' &
5888 ,q2,soilmois(kn),diffu(kn),x2,hydro(kn+1),dtdzs2(kn-1),kn,k
5890 330 rhsmc(k+1)=(soilmois(kn)+q2*rhsmc(k) &
5892 /(zshalf(kn+1)-zshalf(kn)) &
5895 ! --- moisture balance begins here
5898 umveg=(1.-vegfrac)*soilres
5909 !-- total liquid water available on the top of soil domain
5910 !-- without snow - 3 sources of water: precipitation,
5911 !-- water dripping from the canopy and dew
5912 !-- with snow - only one source of water - snow melt
5916 ! totliq=umveg*prcp-drip/delt-umveg*dew*ras-smelt
5918 totliq=prcp-drip/delt-umveg*dew*ras-smelt
5919 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
5920 print *,'umveg*prcp,drip/delt,umveg*dew*ras,smelt', &
5921 umveg*prcp,drip/delt,umveg*dew*ras,smelt
5927 ! ----------- frozen ground version -------------------------
5928 ! reference frozen ground parameter, cvfrz, is a shape parameter of
5929 ! areal distribution function of soil ice content which equals 1/cv.
5930 ! cv is a coefficient of spatial variation of soil ice content.
5931 ! based on field data cv depends on areal mean of frozen depth, and it
5932 ! close to constant = 0.6 if areal mean frozen depth is above 20 cm.
5933 ! that is why parameter cvfrz = 3 (int{1/0.6*0.6})
5935 ! current logic doesn't allow cvfrz be bigger than 3
5938 !-- schaake/koren expression for calculation of max infiltration
5943 f2max=dqm*(zshalf(3)-zshalf(2))
5944 f1=f1max*(1.-soilmois(1)/dqm)
5945 dice=soilice(1)*zshalf(2)
5948 dice=dice+(zshalf(k+1)-zshalf(k))*soilice(k)
5949 fkmax=dqm*(zshalf(k+1)-zshalf(k))
5950 fk=fkmax*(1.-soilmois(k)/dqm)
5953 kdt=refkdt*ksat/refdk
5954 val=(1.-exp(-kdt*delt1))
5957 if(px.lt.0.0) px = 0.0
5959 infmax1 = (px*(ddt/(px+ddt)))/delt
5963 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
5964 print *,'infmax1 before frozen part',infmax1
5967 ! ----------- frozen ground version --------------------------
5968 ! reduction of infiltration based on frozen ground parameters
5970 ! ------------------------------------------------------------------
5972 frzx= 0.15*((dqm+qmin)/ref) * (0.412 / 0.468)
5974 if ( dice .gt. 1.e-2) then
5975 acrt = cvfrz * frzx / dice
5983 sum = sum + (acrt ** ( cvfrz-jk)) / float (k)
5985 fcr = 1. - exp(-acrt) * sum
5987 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
5988 print *,'fcr--------',fcr
5989 print *,'dice=',dice
5991 infmax1 = infmax1* fcr
5992 ! -------------------------------------------------------------------
5994 infmax = max(infmax1,hydro(1)*soilmois(1))
5995 infmax = min(infmax, -totliq)
5996 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
5997 print *,'infmax,infmax1,hydro(1)*soiliqw(1),-totliq', &
5998 infmax,infmax1,hydro(1)*soiliqw(1),-totliq
6001 if (-totliq.gt.infmax)then
6002 runoff=-totliq-infmax
6004 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
6005 print *,'flx,runoff1=',flx,runoff
6008 ! infiltrp is total infiltration flux in m/s
6010 ! solution of moisture budget
6013 flx=flx-soilmois(1)*r7
6014 ! r8 is for direct evaporation from soil, which occurs
6015 ! only from snow-free areas
6017 r8=umveg*r6*(1.-snowfrac)
6022 !-- evaporation regime
6024 qq=(r5*r2-flx+r9)/(r4-r5*r1-r10*r8/(ref-qmin))
6025 flxsat=-dqm*(r4-r5*r1-r10*r8/(ref-qmin)) &
6028 !-- dew formation regime
6029 qq=(r2*r5-flx+r8*(qtot-qcg-qvg)+r9)/(r4-r1*r5)
6030 flxsat=-dqm*(r4-r1*r5)+r2*r5+r8*(qtot-qvg-qcg)+r9
6034 ! print *,'negative qq=',qq
6037 else if(qq.gt.dqm) then
6040 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
6041 print *,'flxsat,flx,delt',flxsat,flx,delt,runoff2
6043 ! runoff2=(flxsat-flx)
6044 runoff=runoff+(flxsat-flx)
6046 soilmois(1)=min(dqm,max(1.e-8,qq))
6049 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
6050 print *,'soilmois,soiliqw, soilice',soilmois,soiliqw,soilice*riw
6051 print *,'cosmc,rhsmc',cosmc,rhsmc
6053 !--- final solution for soilmois
6057 qq=cosmc(kk)*soilmois(k-1)+rhsmc(kk)
6058 ! qq=cosmc(kk)*soiliqw(k-1)+rhsmc(kk)
6061 ! print *,'negative qq=',qq
6064 else if(qq.gt.dqm) then
6068 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
6069 print *,'hydro(k),qq,dqm,k',hydro(k),qq,dqm,k
6071 runoff2=runoff2+((qq-dqm)*(zsmain(k)-zshalf(k)))/delt
6073 runoff2=runoff2+((qq-dqm)*(zshalf(k+1)-zshalf(k)))/delt
6076 soilmois(k)=min(dqm,max(1.e-8,qq))
6079 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
6080 print *,'end soilmois,soiliqw,soilice',soilmois,soiliqw,soilice*riw
6083 mavail=max(.00001,min(1.,(soilmois(1)/(ref-qmin)*(1.-snowfrac)+1.*snowfrac)))
6087 !-------------------------------------------------------------------
6088 end subroutine soilmoist
6089 !-------------------------------------------------------------------
6092 subroutine soilprop(spp_lsm,rstochcol,fieldcol_sf, &
6093 !--- input variables
6094 nzs,fwsat,lwsat,tav,keepfr, &
6095 soilmois,soiliqw,soilice, &
6096 soilmoism,soiliqwm,soilicem, &
6097 !--- soil fixed fields
6098 qwrtz,rhocs,dqm,qmin,psis,bclh,ksat, &
6100 riw,xlmelt,cp,g0_p,cvw,ci, &
6102 !--- output variables
6103 thdif,diffu,hydro,cap)
6105 !******************************************************************
6106 ! soilprop computes thermal diffusivity, and diffusional and
6107 ! hydraulic condeuctivities
6108 !******************************************************************
6109 ! nx,ny,nzs - dimensions of soil domain
6110 ! fwsat, lwsat - volumetric content of frozen and liquid water
6111 ! for saturated condition at given temperatures (m^3/m^3)
6112 ! tav - temperature averaged for soil layers (k)
6113 ! soilmois - volumetric soil moisture at the main soil levels (m^3/m^3)
6114 ! soilmoism - volumetric soil moisture averaged for layers (m^3/m^3)
6115 ! soiliqwm - volumetric liquid soil moisture averaged for layers (m^3/m^3)
6116 ! soilicem - volumetric content of soil ice averaged for layers (m^3/m^3)
6117 ! thdif - thermal diffusivity for soil layers (w/m/k)
6118 ! diffu - diffusional conductivity (m^2/s)
6119 ! hydro - hydraulic conductivity (m/s)
6120 ! cap - volumetric heat capacity (j/m^3/k)
6122 !******************************************************************
6125 !-----------------------------------------------------------------
6127 !--- soil properties
6128 integer, intent(in ) :: nzs
6130 intent(in ) :: rhocs, &
6138 real, dimension( 1:nzs ) , &
6139 intent(in ) :: soilmois, &
6143 real, intent(in ) :: cp, &
6152 real, dimension(1:nzs), intent(in) :: rstochcol
6153 real, dimension(1:nzs), intent(inout) :: fieldcol_sf
6154 integer, intent(in ) :: spp_lsm
6157 !--- output variables
6158 real, dimension(1:nzs) , &
6159 intent(inout) :: cap,diffu,hydro , &
6163 soilicem,soiliqwm , &
6166 !--- local variables
6167 real, dimension(1:nzs) :: hk,detal,kasat,kjpl
6169 real :: x,x1,x2,x4,ws,wd,fact,fach,facd,psif,ci
6170 real :: tln,tavln,tn,pf,a,am,ame,h
6173 !-- for johansen thermal conductivity
6174 real :: kzero,gamd,kdry,kas,x5,sr,ke
6179 !-- constants for johansen (1975) thermal conductivity
6180 kzero =2. ! if qwrtz > 0.2
6191 x1=xlmelt/(g0_p*psis)
6194 !--- next 3 lines are for johansen thermal conduct.
6196 kdry=(0.135*gamd+64.7)/(2700.-0.947*gamd)
6197 !-- one more option from christa's paper
6198 if(qwrtz > 0.2) then
6199 kas=kqwrtz**qwrtz*kzero**(1.-qwrtz)
6201 kas=kqwrtz**qwrtz*3.**(1.-qwrtz)
6206 wd=ws - riw*soilicem(k)
6207 psif=psis*100.*(wd/(soiliqwm(k)+qmin))**bclh &
6209 !--- psif should be in [cm] to compute pf
6211 fact=1.+riw*soilicem(k)
6212 !--- hk is for mccumber thermal conductivity
6214 hk(k)=420.*exp(-(pf+2.7))*fact
6219 if(soilicem(k).ne.0.and.tn.lt.0.) then
6220 !--- detal is taking care of energy spent on freezing or released from
6221 ! melting of soil water
6223 detal(k)=273.15*x2/(tav(k)*tav(k))* &
6224 (tav(k)/(x1*tn))**x4
6226 if(keepfr(k).eq.1.) then
6232 !--- next 10 lines calculate johansen thermal conductivity kjpl
6233 kasat(k)=kas**(1.-ws)*kice**fwsat(k) &
6236 x5=(soilmoism(k)+qmin)/ws
6237 if(soilicem(k).eq.0.) then
6240 !--- next 2 lines - for coarse soils
6242 ! ke=0.7*log10(sr)+1.
6247 kjpl(k)=ke*(kasat(k)-kdry)+kdry
6249 !--- cap -volumetric heat capacity
6250 cap(k)=(1.-ws)*rhocs &
6251 + (soiliqwm(k)+qmin)*cvw &
6253 + (dqm-soilmoism(k))*cp*1.2 &
6254 - detal(k)*1.e3*xlmelt
6258 if((ws-a).lt.0.12)then
6261 h=max(0.,(soilmoism(k)+qmin-a)/(max(1.e-8,(ws-a))))
6263 if(a.ne.0.)facd=1.-a/max(1.e-8,soilmoism(k))
6264 ame=max(1.e-8,ws-riw*soilicem(k))
6265 !--- diffu is diffusional conductivity of soil water
6266 diffu(k)=-bclh*ksat*psis/ame* &
6271 !--- thdif - thermal diffusivity
6272 ! thdif(k)=hk(k)/cap(k)
6273 !--- use thermal conductivity from johansen (1975)
6274 thdif(k)=kjpl(k)/cap(k)
6278 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
6279 print *,'soilice*riw,soiliqw,soilmois,ws',soilice*riw,soiliqw,soilmois,ws
6283 if((ws-riw*soilice(k)).lt.0.12)then
6287 if(soilice(k).ne.0.) &
6288 fach=1.-riw*soilice(k)/max(1.e-8,soilmois(k))
6289 am=max(1.e-8,ws-riw*soilice(k))
6290 !--- hydro is hydraulic conductivity of soil water
6291 hydro(k)=min(ksat,ksat/am* &
6295 if(hydro(k)<1.e-10)hydro(k)=0.
6300 !-----------------------------------------------------------------------
6301 end subroutine soilprop
6302 !-----------------------------------------------------------------------
6305 subroutine transf(i,j, &
6306 !--- input variables
6307 nzs,nroot,soiliqw,tabs,lai,gswin, &
6308 !--- soil fixed fields
6309 dqm,qmin,ref,wilt,zshalf,pc,iland, &
6310 !--- output variables
6313 !-------------------------------------------------------------------
6314 !--- tranf(k) - the transpiration function (Smirnova et al. 1996, eq. 18,19)
6315 !*******************************************************************
6316 ! nx,ny,nzs - dimensions of soil domain
6317 ! soiliqw - volumetric liquid soil moisture at the main levels (m^3/m^3)
6318 ! tranf - the transpiration function at levels (m)
6319 ! transum - transpiration function integrated over the rooting zone (m)
6321 !*******************************************************************
6323 !-------------------------------------------------------------------
6325 !--- input variables
6327 integer, intent(in ) :: i,j,nroot,nzs, iland
6330 intent(in ) :: gswin, tabs, lai
6331 !--- soil properties
6333 intent(in ) :: dqm, &
6339 real, dimension(1:nzs), intent(in) :: soiliqw, &
6343 real, dimension(1:nzs), intent(out) :: tranf
6344 real, intent(out) :: transum
6350 !-- for non-linear root distribution
6351 real :: gx,sm1,sm2,sm3,sm4,ap0,ap1,ap2,ap3,ap4
6352 real :: ftem, pctot, fsol, f1, cmin, cmax, totcnd
6353 real, dimension(1:nzs) :: part
6354 !--------------------------------------------------------------------
6362 totliq=soiliqw(1)+qmin
6372 gx=ap0+ap1*sm1+ap2*sm2+ap3*sm3+ap4*sm4
6373 if(totliq.ge.ref) gx=1.
6374 if(totliq.le.0.) gx=0.
6379 if(totliq.gt.ref) then
6381 else if(totliq.le.wilt) then
6384 tranf(1)=(totliq-wilt)/(ref-wilt)*did
6386 !-- uncomment next line for non-linear root distribution
6390 totliq=soiliqw(k)+qmin
6395 gx=ap0+ap1*sm1+ap2*sm2+ap3*sm3+ap4*sm4
6396 if(totliq.ge.ref) gx=1.
6397 if(totliq.le.0.) gx=0.
6400 did=zshalf(k+1)-zshalf(k)
6402 if(totliq.ge.ref) then
6404 else if(totliq.le.wilt) then
6407 tranf(k)=(totliq-wilt) &
6410 !-- uncomment next line for non-linear root distribution
6414 ! for lai> 3 => transpiration at potential rate (f.tardieu, 2013)
6419 !- 26aug16- next 2 lines could lead to lh increase and higher 2-m q during the day
6420 ! pctot=min(0.8,pc*lai)
6421 ! pctot=min(0.8,max(pc,pc*lai))
6423 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
6424 print *,'i,j,pctot,lai,pc',i,j,pctot,lai,pc
6427 !--- air temperature function
6428 ! Avissar (1985) and Ax 7/95
6429 if (tabs .le. 302.15) then
6430 ftem = 1.0 / (1.0 + exp(-0.41 * (tabs - 282.05)))
6432 ftem = 1.0 / (1.0 + exp(0.5 * (tabs - 314.0)))
6434 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
6435 print *,'i,j,tabs,ftem',i,j,tabs,ftem
6437 !--- incoming solar function
6438 cmin = 1./rsmax_data
6439 cmax = 1./rstbl(iland)
6441 cmax = lai/rstbl(iland) ! max conductance
6443 ! noihlan & planton (1988)
6445 ! if(lai > 0.01) then
6446 ! f1 = 1.1/lai*gswin/rgltbl(iland)! f1=0. when gswin=0.
6447 ! fsol = (f1+cmin/cmax)/(1.+f1)
6452 ! totcnd = max(lai/rstbl(iland), pctot * ftem * f1)
6453 ! Mahrer & Avissar (1982), Avissar et al. (1985)
6454 if (gswin < rgltbl(iland)) then
6455 fsol = 1. / (1. + exp(-0.034 * (gswin - 3.5)))
6459 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
6460 print *,'i,j,gswin,lai,f1,fsol',i,j,gswin,lai,f1,fsol
6462 !--- total conductance
6463 totcnd =(cmin + (cmax - cmin)*pctot*ftem*fsol)/cmax
6465 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
6466 print *,'i,j,iland,rgltbl(iland),rstbl(iland),rsmax_data,totcnd' &
6467 ,i,j,iland,rgltbl(iland),rstbl(iland),rsmax_data,totcnd
6470 !-- transum - total for the rooting zone
6473 ! linear root distribution
6474 tranf(k)=max(cmin,tranf(k)*totcnd)
6475 transum=transum+tranf(k)
6477 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
6478 print *,'i,j,transum,tranf',i,j,transum,tranf
6481 !-----------------------------------------------------------------
6482 end subroutine transf
6483 !-----------------------------------------------------------------
6486 subroutine vilka(tn,d1,d2,pp,qs,ts,tt,nstep,ii,j,iland,isoil)
6487 !--------------------------------------------------------------
6488 !--- vilka finds the solution of energy budget at the surface
6489 !--- using table t,qs computed from clausius-klapeiron
6490 !--------------------------------------------------------------
6491 real, dimension(1:5001), intent(in ) :: tt
6492 real, intent(in ) :: tn,d1,d2,pp
6493 integer, intent(in ) :: nstep,ii,j,iland,isoil
6495 real, intent(out ) :: qs, ts
6500 i=(tn-1.7315e2)/.05+1
6501 t1=173.1+float(i)*.05
6503 i1=i-f1/(.05+d1*(tt(i+1)-tt(i)))
6505 if(i.gt.5000.or.i.lt.1) goto 1
6507 t1=173.1+float(i)*.05
6509 rn=f1/(.05+d1*(tt(i+1)-tt(i)))
6511 if(i.gt.5000.or.i.lt.1) goto 1
6514 qs=(tt(i)+(tt(i)-tt(i+1))*rn)/pp
6516 ! 1 print *,'crash in surface energy budget - stop'
6517 1 print *,' avost in vilka table index= ',i
6518 ! print *,tn,d1,d2,pp,nstep,i,tt(i),ii,j,iland,isoil
6519 print *,'i,j=',ii,j,'lu_index = ',iland, 'psfc[hpa] = ',pp, 'tsfc = ',tn
6520 call wrf_error_fatal (' crash in surface energy budget ' )
6522 !-----------------------------------------------------------------------
6523 end subroutine vilka
6524 !-----------------------------------------------------------------------
6526 subroutine soilvegin ( mosaic_lu,mosaic_soil,soilfrac,nscat, &
6528 nlcat,ivgtyp,isltyp,iswater,myj, &
6529 iforest,lufrac,vegfrac,emiss,pc,znt,lai,rdlai2d,&
6530 qwrtz,rhocs,bclh,dqm,ksat,psis,qmin,ref,wilt,i,j)
6532 !************************************************************************
6533 ! set-up soil and vegetation parameters in the case when
6534 ! snow disappears during the forecast and snow parameters
6535 ! shold be replaced by surface parameters according to
6536 ! soil and vegetation types in this point.
6542 ! dqm: max soil moisture content - min (m^3/m^3)
6543 ! ref: reference soil moisture (m^3/m^3)
6544 ! wilt: wilting pt soil moisture contents (m^3/m^3)
6545 ! qmin: air dry soil moist content limits (m^3/m^3)
6546 ! psis: sat soil potential coefs. (m)
6547 ! ksat: sat soil diffusivity/conductivity coefs. (m/s)
6548 ! bclh: soil diffusivity/conductivity exponent.
6550 ! ************************************************************************
6552 !---------------------------------------------------------------------------
6553 integer, parameter :: nsoilclas=19
6554 integer, parameter :: nvegclas=24+3
6555 integer, parameter :: ilsnow=99
6557 integer, intent(in ) :: nlcat, nscat, iswater, i, j
6559 !--- soiltyp classification according to statsgo(nclasses=16)
6562 ! 2 loamy sand loamy sand
6563 ! 3 sandy loam sandy loam
6564 ! 4 silt loam silty loam
6567 ! 7 sandy clay loam sandy clay loam
6568 ! 8 silty clay loam silty clay loam
6569 ! 9 clay loam clay loam
6570 ! 10 sandy clay sandy clay
6571 ! 11 silty clay silty clay
6572 ! 12 clay light clay
6573 ! 13 organic materials loam
6576 ! bedrock is reclassified as class 14
6577 ! 16 other (land-ice)
6582 !----------------------------------------------------------------------
6583 real lqma(nsoilclas),lrhc(nsoilclas), &
6584 lpsi(nsoilclas),lqmi(nsoilclas), &
6585 lbcl(nsoilclas),lkas(nsoilclas), &
6586 lwil(nsoilclas),lref(nsoilclas), &
6588 !-- lqma rawls et al.[1982]
6589 ! data lqma /0.417, 0.437, 0.453, 0.501, 0.486, 0.463, 0.398,
6590 ! & 0.471, 0.464, 0.430, 0.479, 0.475, 0.439, 1.0, 0.20, 0.401/
6592 !-- clapp, r. and g. hornberger, 1978: empirical equations for some soil
6593 ! hydraulic properties, water resour. res., 14, 601-604.
6595 !-- clapp et al. [1978]
6596 data lqma /0.395, 0.410, 0.435, 0.485, 0.485, 0.451, 0.420, &
6597 0.477, 0.476, 0.426, 0.492, 0.482, 0.451, 1.0, &
6598 0.20, 0.435, 0.468, 0.200, 0.339/
6600 !-- lref rawls et al.[1982]
6601 ! data lref /0.091, 0.125, 0.207, 0.330, 0.360, 0.270, 0.255,
6602 ! & 0.366, 0.318, 0.339, 0.387, 0.396, 0.329, 1.0, 0.108, 0.283/
6604 !-- clapp et al. [1978]
6605 data lref /0.174, 0.179, 0.249, 0.369, 0.369, 0.314, 0.299, &
6606 0.357, 0.391, 0.316, 0.409, 0.400, 0.314, 1., &
6607 0.1, 0.249, 0.454, 0.17, 0.236/
6609 !-- lwil rawls et al.[1982]
6610 ! data lwil/0.033, 0.055, 0.095, 0.133, 0.133, 0.117, 0.148,
6611 ! & 0.208, 0.197, 0.239, 0.250, 0.272, 0.066, 0.0, 0.006, 0.029/
6613 !-- clapp et al. [1978]
6614 data lwil/0.068, 0.075, 0.114, 0.179, 0.179, 0.155, 0.175, &
6615 0.218, 0.250, 0.219, 0.283, 0.286, 0.155, 0.0, &
6616 0.006, 0.114, 0.030, 0.006, 0.01/
6618 ! data lqmi/0.010, 0.028, 0.047, 0.084, 0.084, 0.066, 0.067,
6619 ! & 0.120, 0.103, 0.100, 0.126, 0.138, 0.066, 0.0, 0.006, 0.028/
6621 !-- carsel and parrish [1988]
6622 data lqmi/0.045, 0.057, 0.065, 0.067, 0.034, 0.078, 0.10, &
6623 0.089, 0.095, 0.10, 0.070, 0.068, 0.078, 0.0, &
6624 0.004, 0.065, 0.020, 0.004, 0.008/
6626 !-- lpsi cosby et al[1984]
6627 ! data lpsi/0.060, 0.036, 0.141, 0.759, 0.759, 0.355, 0.135,
6628 ! & 0.617, 0.263, 0.098, 0.324, 0.468, 0.355, 0.0, 0.069, 0.036/
6629 ! & 0.617, 0.263, 0.098, 0.324, 0.468, 0.355, 0.0, 0.069, 0.036/
6631 !-- clapp et al. [1978]
6632 data lpsi/0.121, 0.090, 0.218, 0.786, 0.786, 0.478, 0.299, &
6633 0.356, 0.630, 0.153, 0.490, 0.405, 0.478, 0.0, &
6634 0.121, 0.218, 0.468, 0.069, 0.069/
6636 !-- lkas rawls et al.[1982]
6637 ! data lkas/5.83e-5, 1.70e-5, 7.19e-6, 1.89e-6, 1.89e-6,
6638 ! & 3.67e-6, 1.19e-6, 4.17e-7, 6.39e-7, 3.33e-7, 2.50e-7,
6639 ! & 1.67e-7, 3.38e-6, 0.0, 1.41e-4, 1.41e-5/
6641 !-- clapp et al. [1978]
6642 data lkas/1.76e-4, 1.56e-4, 3.47e-5, 7.20e-6, 7.20e-6, &
6643 6.95e-6, 6.30e-6, 1.70e-6, 2.45e-6, 2.17e-6, &
6644 1.03e-6, 1.28e-6, 6.95e-6, 0.0, 1.41e-4, &
6645 3.47e-5, 1.28e-6, 1.41e-4, 1.76e-4/
6647 !-- lbcl cosby et al [1984]
6648 ! data lbcl/2.79, 4.26, 4.74, 5.33, 5.33, 5.25, 6.66,
6649 ! & 8.72, 8.17, 10.73, 10.39, 11.55, 5.25, 0.0, 2.79, 4.26/
6651 !-- clapp et al. [1978]
6652 data lbcl/4.05, 4.38, 4.90, 5.30, 5.30, 5.39, 7.12, &
6653 7.75, 8.52, 10.40, 10.40, 11.40, 5.39, 0.0, &
6654 4.05, 4.90, 11.55, 2.79, 2.79/
6656 data lrhc /1.47,1.41,1.34,1.27,1.27,1.21,1.18,1.32,1.23, &
6657 1.18,1.15,1.09,1.21,4.18,2.03,2.10,1.09,2.03,1.47/
6659 data datqtz/0.92,0.82,0.60,0.25,0.10,0.40,0.60,0.10,0.35, &
6660 0.52,0.10,0.25,0.00,0.,0.60,0.0,0.25,0.60,0.92/
6662 !--------------------------------------------------------------------------
6664 ! usgs vegetation types
6666 ! 1: urban and built-up land
6667 ! 2: dryland cropland and pasture
6668 ! 3: irrigated cropland and pasture
6669 ! 4: mixed dryland/irrigated cropland and pasture
6670 ! 5: cropland/grassland mosaic
6671 ! 6: cropland/woodland mosaic
6674 ! 9: mixed shrubland/grassland
6676 ! 11: deciduous broadleaf forest
6677 ! 12: deciduous needleleaf forest
6678 ! 13: evergreen broadleaf forest
6679 ! 14: evergreen needleleaf fores
6682 ! 17: herbaceous wetland
6683 ! 18: wooded wetland
6684 ! 19: barren or sparsely vegetated
6685 ! 20: herbaceous tundra
6688 ! 23: bare ground tundra
6695 ! modis vegetation categories from VEGPARM.TBL
6696 ! 1: evergreen needleleaf forest
6697 ! 2: evergreen broadleaf forest
6698 ! 3: deciduous needleleaf forest
6699 ! 4: deciduous broadleaf forest
6701 ! 6: closed shrublands
6702 ! 7: open shrublands
6706 ! 11: permanent wetlands
6708 ! 13: urban and built-up
6709 ! 14: cropland/natural vegetation mosaic
6711 ! 16: barren or sparsely vegetated
6719 !---- below are the arrays for the vegetation parameters
6720 real lalb(nvegclas),lmoi(nvegclas),lemi(nvegclas), &
6721 lrou(nvegclas),lthi(nvegclas),lsig(nvegclas), &
6724 !************************************************************************
6725 !---- vegetation parameters
6729 data lalb/.18,.17,.18,.18,.18,.16,.19,.22,.20,.20,.16,.14, &
6730 .12,.12,.13,.08,.14,.14,.25,.15,.15,.15,.25,.55, &
6732 data lemi/.88,4*.92,.93,.92,.88,.9,.92,.93,.94, &
6733 .95,.95,.94,.98,.95,.95,.85,.92,.93,.92,.85,.95, &
6735 !-- roughness length is changed for forests and some others
6736 ! data lrou/.5,.06,.075,.065,.05,.2,.075,.1,.11,.15,.8,.85, &
6737 ! 2.0,1.0,.563,.0001,.2,.4,.05,.1,.15,.1,.065,.05/
6738 data lrou/.5,.06,.075,.065,.05,.2,.075,.1,.11,.15,.5,.5, &
6739 .5,.5,.5,.0001,.2,.4,.05,.1,.15,.1,.065,.05, &
6742 data lmoi/.1,.3,.5,.25,.25,.35,.15,.1,.15,.15,.3,.3, &
6743 .5,.3,.3,1.,.6,.35,.02,.5,.5,.5,.02,.95,.40,.50,.40/
6745 !---- still needs to be corrected
6747 ! data lpc/ 15*.8,0.,.8,.8,.5,.5,.5,.5,.5,.0/
6748 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, &
6749 0.3,0.3,0.4,0.4,0.3,0.,.3,0.,0./
6751 ! used in ruc data lpc /0.6,6*0.8,0.7,0.75,6*0.8,0.,0.8,0.8, &
6752 ! 0.5,0.7,0.6,0.7,0.5,0./
6755 !***************************************************************************
6761 integer, intent(in ) :: mosaic_lu, mosaic_soil
6763 logical, intent(in ) :: myj
6764 real, intent(in ) :: shdmax
6765 real, intent(in ) :: shdmin
6766 real, intent(in ) :: vegfrac
6767 real, dimension( 1:nlcat ), intent(in):: lufrac
6768 real, dimension( 1:nscat ), intent(in):: soilfrac
6774 intent (inout ) :: emiss, &
6777 logical, intent(in) :: rdlai2d
6778 !--- soil properties
6780 intent( out) :: rhocs, &
6789 integer, intent ( out) :: iforest
6791 ! integer, dimension( 1:(lucats) ) , &
6792 ! intent ( out) :: iforest
6795 ! integer, dimension( 1:50 ) :: if1
6796 integer :: kstart, kfin, lstart, lfin
6798 real :: area, factor, znt1, lb
6799 real, dimension( 1:nlcat ) :: znttoday, laitoday, deltalai
6801 !***********************************************************************
6802 ! data zs1/0.0,0.05,0.20,0.40,1.6,3.0/ ! o - levels in soil
6803 ! data zs2/0.0,0.025,0.125,0.30,1.,2.3/ ! x - levels in soil
6805 ! data if1/12*0,1,1,1,12*0/
6811 iforest = ifortbl(ivgtyp)
6813 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
6814 if(i.eq.375.and.j.eq.254)then
6815 print *,'ifortbl(ivgtyp),ivgtyp,laitbl(ivgtyp),z0tbl(ivgtyp)', &
6816 ifortbl(ivgtyp),ivgtyp,laitbl(ivgtyp),z0tbl(ivgtyp)
6822 ! 11oct2012 - seasonal correction on znt for crops and lai for all veg. types
6823 ! factor = 1 with minimum greenness --> vegfrac = shdmin (cold season)
6824 ! factor = 0 with maximum greenness --> vegfrac = shdmax
6825 ! shdmax, shdmin and vegfrac are in % here.
6826 if((shdmax - shdmin) .lt. 1) then
6827 factor = 1. ! min greenness
6829 factor = 1. - max(0.,min(1.,(vegfrac - shdmin)/max(1.,(shdmax-shdmin))))
6832 ! 18sept18 - laitbl and z0tbl are the max values
6834 if(ifortbl(k) == 1) deltalai(k)=min(0.2,0.8*laitbl(k))
6835 if(ifortbl(k) == 2 .or. ifortbl(k) == 7) deltalai(k)=min(0.5,0.8*laitbl(k))
6836 if(ifortbl(k) == 3) deltalai(k)=min(0.45,0.8*laitbl(k))
6837 if(ifortbl(k) == 4) deltalai(k)=min(0.75,0.8*laitbl(k))
6838 if(ifortbl(k) == 5) deltalai(k)=min(0.86,0.8*laitbl(k))
6840 if(k.ne.iswater) then
6841 !-- 20aug18 - change in laitoday based on the greenness fraction for the current day
6842 laitoday(k) = laitbl(k) - deltalai(k) * factor
6844 if(ifortbl(k) == 7) then
6845 !-- seasonal change of roughness length for crops
6846 znttoday(k) = z0tbl(k) - 0.125 * factor
6848 znttoday(k) = z0tbl(k)
6851 laitoday(k) = laitbl(k)
6852 znttoday(k) = znt ! do not overwrite z0 over water with the table value
6856 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
6857 if(i.eq.358.and.j.eq.260)then
6858 print *,'ivgtyp,factor,vegfrac,shdmin,shdmax,deltalai,laitoday(ivgtyp),znttoday(ivgtyp)', &
6859 i,j,ivgtyp,factor,vegfrac,shdmin,shdmax,deltalai,laitoday(ivgtyp),znttoday(ivgtyp)
6867 if(.not.rdlai2d) lai = 0.
6869 !-- mosaic approach to landuse in the grid box
6870 ! use mason (1988) eq.(15) to compute effective znt;
6871 ! lb - blending height = l/200., where l is the length scale of regions with varying z0 (lb = 5 if l=1000 m)
6873 if(mosaic_lu == 1) then
6875 area = area + lufrac(k)
6876 emiss = emiss+ lemitbl(k)*lufrac(k)
6877 znt = znt + lufrac(k)/alog(lb/znttoday(k))**2.
6878 ! znt1 - weighted average in the grid box, not used, computed for comparison
6879 znt1 = znt1 + lufrac(k)*znttoday(k)
6880 if(.not.rdlai2d) lai = lai + laitoday(k)*lufrac(k)
6881 pc = pc + pctbl(k)*lufrac(k)
6884 if (area.gt.1.) area=1.
6885 if (area <= 0.) then
6886 print *,'bad area of grid box', area
6890 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
6891 if(i.eq.358.and.j.eq.260) then
6892 print *,'area=',area,i,j,ivgtyp,nlcat,(lufrac(k),k=1,nlcat),emiss,znt,znt1,lai,pc
6898 znt = lb/exp(sqrt(1./znt))
6899 if(.not.rdlai2d) lai = lai/area
6902 if ( wrf_at_debug_level(lsmruc_dbg_lvl) ) then
6903 if(i.eq.358.and.j.eq.260) then
6904 print *,'mosaic=',i,j,ivgtyp,nlcat,(lufrac(k),k=1,nlcat),emiss,znt,znt1,lai,pc
6910 emiss = lemitbl(ivgtyp)
6911 znt = znttoday(ivgtyp)
6913 if(.not.rdlai2d) lai = laitoday(ivgtyp)
6916 ! parameters from soilparm.tbl
6928 if(mosaic_soil == 1 ) then
6930 if(k.ne.14) then ! statsgo value for water
6931 !exclude watrer points from this loop
6932 area = area + soilfrac(k)
6933 rhocs = rhocs + hc(k)*1.e6*soilfrac(k)
6934 bclh = bclh + bb(k)*soilfrac(k)
6935 dqm = dqm + (maxsmc(k)- &
6936 drysmc(k))*soilfrac(k)
6937 ksat = ksat + satdk(k)*soilfrac(k)
6938 psis = psis - satpsi(k)*soilfrac(k)
6939 qmin = qmin + drysmc(k)*soilfrac(k)
6940 ref = ref + refsmc(k)*soilfrac(k)
6941 wilt = wilt + wltsmc(k)*soilfrac(k)
6942 qwrtz = qwrtz + qtz(k)*soilfrac(k)
6945 if (area.gt.1.) area=1.
6946 if (area <= 0.) then
6947 ! area = 0. for water points
6948 ! print *,'area of a grid box', area, 'iswater = ',iswater
6949 rhocs = hc(isltyp)*1.e6
6951 dqm = maxsmc(isltyp)- &
6953 ksat = satdk(isltyp)
6954 psis = - satpsi(isltyp)
6955 qmin = drysmc(isltyp)
6956 ref = refsmc(isltyp)
6957 wilt = wltsmc(isltyp)
6971 ! dominant category approach
6973 if(isltyp.ne.14) then
6974 rhocs = hc(isltyp)*1.e6
6976 dqm = maxsmc(isltyp)- &
6978 ksat = satdk(isltyp)
6979 psis = - satpsi(isltyp)
6980 qmin = drysmc(isltyp)
6981 ref = refsmc(isltyp)
6982 wilt = wltsmc(isltyp)
6987 !--------------------------------------------------------------------------
6988 end subroutine soilvegin
6989 !--------------------------------------------------------------------------
6991 subroutine ruclsminit( sh2o,smfr3d,tslb,smois,isltyp,ivgtyp, &
6992 mminlu, xice,mavail,nzs, iswater, isice, &
6993 znt, restart, allowed_to_read , &
6994 ids,ide, jds,jde, kds,kde, &
6995 ims,ime, jms,jme, kms,kme, &
6996 its,ite, jts,jte, kts,kte )
6997 #if ( wrf_chem == 1 )
6998 use module_data_gocart_dust
7003 integer, intent(in ) :: ids,ide, jds,jde, kds,kde, &
7004 ims,ime, jms,jme, kms,kme, &
7005 its,ite, jts,jte, kts,kte, &
7007 character(len=*), intent(in ) :: mminlu
7009 real, dimension( ims:ime, 1:nzs, jms:jme ) , &
7010 intent(in) :: tslb, &
7013 integer, dimension( ims:ime, jms:jme ) , &
7014 intent(inout) :: isltyp,ivgtyp
7016 real, dimension( ims:ime, 1:nzs, jms:jme ) , &
7017 intent(inout) :: smfr3d, &
7020 real, dimension( ims:ime, jms:jme ) , &
7021 intent(inout) :: xice,mavail
7023 real, dimension( ims:ime, jms:jme ) , &
7026 real, dimension ( 1:nzs ) :: soiliqw
7028 logical , intent(in) :: restart, allowed_to_read
7031 integer :: i,j,l,itf,jtf
7032 real :: riw,xlmelt,tln,dqm,ref,psis,qmin,bclh
7034 character*8 :: mminluruc, mminsl
7041 ! initialize three lsm related tables
7042 if ( allowed_to_read ) then
7043 call wrf_message( 'initialize three lsm related tables' )
7044 if(mminlu == 'USGS') then
7045 mminluruc='USGS-RUC'
7046 elseif(mminlu == 'MODIS' .or. &
7047 & mminlu == 'MODIFIED_IGBP_MODIS_NOAH') then
7048 mminluruc='MODI-RUC'
7051 print *,'ruclsminit uses ',mminluruc
7052 call ruclsm_soilvegparm( mminluruc, mminsl)
7055 !#if ( wrf_chem == 1 )
7057 ! need this parameter for dust parameterization in wrf/chem
7060 ! porosity(i)=maxsmc(i)
7061 ! drypoint(i)=drysmc(i)
7065 if(.not.restart)then
7073 if ( isltyp( i,j ) .lt. 1 ) then
7075 write(err_message,*)"module_sf_ruclsm.f: lsminit: out of range isltyp ",i,j,isltyp( i,j )
7076 call wrf_message(err_message)
7080 if ( errflag .eq. 1 ) then
7081 call wrf_error_fatal( "module_sf_ruclsm.f: lsminit: out of range value "// &
7082 "of isltyp. is this field in the input?" )
7088 znt(i,j) = z0tbl(ivgtyp(i,j))
7090 !--- computation of volumetric content of ice in soil
7091 !--- and initialize mavail
7092 dqm = maxsmc (isltyp(i,j)) - &
7093 drysmc (isltyp(i,j))
7094 ref = refsmc (isltyp(i,j))
7095 psis = - satpsi (isltyp(i,j))
7096 qmin = drysmc (isltyp(i,j))
7097 bclh = bb (isltyp(i,j))
7100 if(xice(i,j).gt.0.) then
7108 if(isltyp(i,j).ne.14 ) then
7110 mavail(i,j) = max(0.00001,min(1.,(smois(i,1,j)-qmin)/(ref-qmin)))
7112 !-- for land points initialize soil ice
7113 tln=log(tslb(i,l,j)/273.15)
7116 soiliqw(l)=(dqm+qmin)*(xlmelt* &
7117 (tslb(i,l,j)-273.15)/tslb(i,l,j)/9.81/psis) &
7119 soiliqw(l)=max(0.,soiliqw(l))
7120 soiliqw(l)=min(soiliqw(l),smois(i,l,j))
7121 sh2o(i,l,j)=soiliqw(l)
7122 smfr3d(i,l,j)=(smois(i,l,j)-soiliqw(l))/riw
7126 sh2o(i,l,j)=smois(i,l,j)
7131 !-- for water isltyp=14
7145 !-----------------------------------------------------------------
7146 end subroutine ruclsminit
7147 !-----------------------------------------------------------------
7149 subroutine ruclsm_soilvegparm( mminluruc, mminsl)
7150 !-----------------------------------------------------------------
7154 integer :: LUMATCH, IINDEX, LC, NUM_SLOPE
7156 INTEGER , PARAMETER :: OPEN_OK = 0
7158 character*8 :: MMINLURUC, MMINSL
7159 character*128 :: mess , message, vege_parm_string
7160 logical, external :: wrf_dm_on_monitor
7163 !-----specify vegetation related characteristics :
7164 ! albbck: sfc albedo (in percentage)
7165 ! z0: roughness length (m)
7167 ! pc: plant coefficient for transpiration function
7168 ! -- the rest of the parameters are read in but not used currently
7169 ! shdfac: green vegetation fraction (in percentage)
7170 ! note: the albedo, z0, and shdfac values read from the following table
7171 ! albedo, amd z0 are specified in land-use table; and shdfac is
7172 ! the monthly green vegetation data
7173 ! cmxtbl: max cnpy capacity (m)
7174 ! rsmin: mimimum stomatal resistance (s m-1)
7175 ! rsmax: max. stomatal resistance (s m-1)
7176 ! rgl: parameters used in radiation stress function
7177 ! hs: parameter used in vapor pressure deficit functio
7178 ! topt: optimum transpiration air temperature. (k)
7179 ! cmcmax: maximum canopy water capacity
7180 ! cfactr: parameter used in the canopy inteception calculati
7181 ! snup: threshold snow depth (in water equivalent m) that
7182 ! implies 100% snow cover
7183 ! lai: leaf area index (dimensionless)
7184 ! maxalb: upper bound on maximum albedo over deep snow
7186 !-----read in vegetaion properties from VEGPARM.TBL
7189 IF ( wrf_dm_on_monitor() ) THEN
7191 OPEN(19, FILE='VEGPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
7192 IF(ierr .NE. OPEN_OK ) THEN
7193 WRITE(message,FMT='(A)') &
7194 'module_sf_ruclsm.F: soil_veg_gen_parm: failure opening VEGPARM.TBL'
7195 CALL wrf_error_fatal ( message )
7198 WRITE ( mess, * ) 'INPUT VEGPARM FOR ',MMINLURUC
7199 CALL wrf_message( mess )
7204 READ (19,'(A)') vege_parm_string
7206 READ (19,2000,END=2002)LUTYPE
7207 READ (19,*)LUCATS,IINDEX
7209 WRITE( mess , * ) 'VEGPARM FOR ',LUTYPE,' FOUND', LUCATS,' CATEGORIES'
7210 CALL wrf_message( mess )
7212 IF(LUTYPE.NE.MMINLURUC)THEN ! Skip over the undesired table
7213 write ( mess , * ) 'Skipping ', LUTYPE, ' table'
7214 CALL wrf_message( mess )
7218 inner : DO ! Find the next "Vegetation Parameters"
7219 READ (19,'(A)',END=2002) vege_parm_string
7220 IF (TRIM(vege_parm_string) .EQ. "Vegetation Parameters") THEN
7226 write ( mess , * ) 'Found ', LUTYPE, ' table'
7227 CALL wrf_message( mess )
7228 EXIT outer ! Found the table, read the data
7233 IF (LUMATCH == 1) then
7234 write ( mess , * ) 'Reading ',LUTYPE,' table'
7235 CALL wrf_message( mess )
7237 READ (19,*)IINDEX,ALBTBL(LC),Z0TBL(LC),LEMITBL(LC),PCTBL(LC), &
7238 SHDTBL(LC),IFORTBL(LC),RSTBL(LC),RGLTBL(LC), &
7239 HSTBL(LC),SNUPTBL(LC),LAITBL(LC),MAXALB(LC)
7243 READ (19,*)TOPT_DATA
7245 READ (19,*)CMCMAX_DATA
7247 READ (19,*)CFACTR_DATA
7249 READ (19,*)RSMAX_DATA
7257 READ (19,*,iostat=ierr)URBAN
7258 if ( ierr /= 0 ) call wrf_message ( "-------- VEGPARM.TBL READ ERROR --------")
7259 if ( ierr /= 0 ) call wrf_message ( "Problem read URBAN from VEGPARM.TBL")
7260 if ( ierr /= 0 ) call wrf_message ( " -- Use updated version of VEGPARM.TBL ")
7261 if ( ierr /= 0 ) call wrf_error_fatal ( "Problem read URBAN from VEGPARM.TBL")
7268 IF ( wrf_at_debug_level(lsmruc_dbg_lvl) ) THEN
7269 print *,' LEMITBL, PCTBL, Z0TBL, LAITBL --->', LEMITBL, PCTBL, Z0TBL, LAITBL
7272 IF (LUMATCH == 0) then
7273 CALL wrf_error_fatal ("Land Use Dataset '"//MMINLURUC//"' not found in VEGPARM.TBL.")
7278 CALL wrf_dm_bcast_string ( LUTYPE , 8 )
7279 CALL wrf_dm_bcast_integer ( LUCATS , 1 )
7280 CALL wrf_dm_bcast_integer ( IINDEX , 1 )
7281 CALL wrf_dm_bcast_integer ( LUMATCH , 1 )
7282 CALL wrf_dm_bcast_real ( ALBTBL , NLUS )
7283 CALL wrf_dm_bcast_real ( Z0TBL , NLUS )
7284 CALL wrf_dm_bcast_real ( LEMITBL , NLUS )
7285 CALL wrf_dm_bcast_real ( PCTBL , NLUS )
7286 CALL wrf_dm_bcast_real ( SHDTBL , NLUS )
7287 CALL wrf_dm_bcast_real ( IFORTBL , NLUS )
7288 CALL wrf_dm_bcast_real ( RSTBL , NLUS )
7289 CALL wrf_dm_bcast_real ( RGLTBL , NLUS )
7290 CALL wrf_dm_bcast_real ( HSTBL , NLUS )
7291 CALL wrf_dm_bcast_real ( SNUPTBL , NLUS )
7292 CALL wrf_dm_bcast_real ( LAITBL , NLUS )
7293 CALL wrf_dm_bcast_real ( MAXALB , NLUS )
7294 CALL wrf_dm_bcast_real ( TOPT_DATA , 1 )
7295 CALL wrf_dm_bcast_real ( CMCMAX_DATA , 1 )
7296 CALL wrf_dm_bcast_real ( CFACTR_DATA , 1 )
7297 CALL wrf_dm_bcast_real ( RSMAX_DATA , 1 )
7298 CALL wrf_dm_bcast_integer ( BARE , 1 )
7299 CALL wrf_dm_bcast_integer ( NATURAL , 1 )
7300 CALL wrf_dm_bcast_integer ( CROP , 1 )
7301 CALL wrf_dm_bcast_integer ( URBAN , 1 )
7304 !-----READ IN SOIL PROPERTIES FROM SOILPARM.TBL
7306 IF ( wrf_dm_on_monitor() ) THEN
7307 OPEN(19, FILE='SOILPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr)
7308 IF(ierr .NE. OPEN_OK ) THEN
7309 WRITE(message,FMT='(A)') &
7310 'module_sf_ruclsm.F: soil_veg_gen_parm: failure opening SOILPARM.TBL'
7311 CALL wrf_error_fatal ( message )
7314 WRITE(mess,*) 'INPUT SOIL TEXTURE CLASSIFICATION = ',MMINSL
7315 CALL wrf_message( mess )
7320 READ (19,2000,END=2003)SLTYPE
7321 READ (19,*)SLCATS,IINDEX
7322 IF(SLTYPE.NE.MMINSL)THEN
7324 READ (19,*) IINDEX,BB(LC),DRYSMC(LC),HC(LC),MAXSMC(LC),&
7325 REFSMC(LC),SATPSI(LC),SATDK(LC), SATDW(LC), &
7330 READ (19,2000,END=2003)SLTYPE
7331 READ (19,*)SLCATS,IINDEX
7333 IF(SLTYPE.EQ.MMINSL)THEN
7334 WRITE( mess , * ) 'SOIL TEXTURE CLASSIFICATION = ',SLTYPE,' FOUND', &
7335 SLCATS,' CATEGORIES'
7336 CALL wrf_message ( mess )
7339 IF(SLTYPE.EQ.MMINSL)THEN
7341 READ (19,*) IINDEX,BB(LC),DRYSMC(LC),HC(LC),MAXSMC(LC),&
7342 REFSMC(LC),SATPSI(LC),SATDK(LC), SATDW(LC), &
7352 CALL wrf_dm_bcast_integer ( LUMATCH , 1 )
7353 CALL wrf_dm_bcast_string ( SLTYPE , 8 )
7354 CALL wrf_dm_bcast_string ( MMINSL , 8 ) ! since this is reset above, see oct2 ^
7355 CALL wrf_dm_bcast_integer ( SLCATS , 1 )
7356 CALL wrf_dm_bcast_integer ( IINDEX , 1 )
7357 CALL wrf_dm_bcast_real ( BB , NSLTYPE )
7358 CALL wrf_dm_bcast_real ( DRYSMC , NSLTYPE )
7359 CALL wrf_dm_bcast_real ( HC , NSLTYPE )
7360 CALL wrf_dm_bcast_real ( MAXSMC , NSLTYPE )
7361 CALL wrf_dm_bcast_real ( REFSMC , NSLTYPE )
7362 CALL wrf_dm_bcast_real ( SATPSI , NSLTYPE )
7363 CALL wrf_dm_bcast_real ( SATDK , NSLTYPE )
7364 CALL wrf_dm_bcast_real ( SATDW , NSLTYPE )
7365 CALL wrf_dm_bcast_real ( WLTSMC , NSLTYPE )
7366 CALL wrf_dm_bcast_real ( QTZ , NSLTYPE )
7368 IF(LUMATCH.EQ.0)THEN
7369 CALL wrf_message( 'SOIl TEXTURE IN INPUT FILE DOES NOT ' )
7370 CALL wrf_message( 'MATCH SOILPARM TABLE' )
7371 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