#define LSMRUC_DBG_LVL 3000 !WRF:MODEL_LAYER:PHYSICS ! MODULE module_sf_ruclsm ! Notes for perturbations of soil properties (Judith Berner) ! Perturbations are applied in subroutine soilprob to array hydro; ! soilprop is called from subroutine SFCTMP which is called from subroutine LSMRUC; ! subroutine LSMRUC had two new 3D fields: pattern_spp_lsm (in) and field_sf(inout); ! their vertical dimension is number of atmospheric levels (kms:kme) - (suboptimal, but easiest hack) ! field_sf is used to pass perturbed fields of hydrop up to model (and output) driver; ! in argument list to SFCTMP the arrays are passed as pattern_spp_lsm(i,1:nzs,j), and exist henceforth as ! column arrays; ! in the subroutines below SFCTMP (SNOW and SNOWSOIL) the fields are called rstochcol,fieldcol_sf ! to reflect their dimension rstochcol (1:nzs) USE module_model_constants USE module_wrf_error ! VEGETATION PARAMETERS INTEGER :: LUCATS , BARE, NATURAL, CROP, URBAN integer, PARAMETER :: NLUS=50 CHARACTER*8 LUTYPE INTEGER, DIMENSION(1:NLUS) :: IFORTBL real, dimension(1:NLUS) :: SNUPTBL, RSTBL, RGLTBL, HSTBL, LAITBL, & ALBTBL, Z0TBL, LEMITBL, PCTBL, SHDTBL, MAXALB REAL :: TOPT_DATA,CMCMAX_DATA,CFACTR_DATA,RSMAX_DATA ! SOIL PARAMETERS INTEGER :: SLCATS INTEGER, PARAMETER :: NSLTYPE=30 CHARACTER*8 SLTYPE REAL, DIMENSION (1:NSLTYPE) :: BB,DRYSMC,HC, & MAXSMC, REFSMC,SATPSI,SATDK,SATDW, WLTSMC,QTZ ! LSM GENERAL PARAMETERS INTEGER :: SLPCATS INTEGER, PARAMETER :: NSLOPE=30 REAL, DIMENSION (1:NSLOPE) :: SLOPE_DATA REAL :: SBETA_DATA,FXEXP_DATA,CSOIL_DATA,SALP_DATA,REFDK_DATA, & REFKDT_DATA,FRZK_DATA,ZBOT_DATA, SMLOW_DATA,SMHIGH_DATA, & CZIL_DATA CHARACTER*256 :: err_message CONTAINS !----------------------------------------------------------------- SUBROUTINE LSMRUC(spp_lsm, & #if (EM_CORE==1) pattern_spp_lsm,field_sf, & #endif DT,KTAU,NSL, & #if (EM_CORE==1) lakemodel,lakemask, & graupelncv,snowncv,rainncv, & #endif ZS,RAINBL,SNOW,SNOWH,SNOWC,FRZFRAC,frpcpn, & rhosnf,precipfr, & ! pass it out to module_diagnostics Z3D,P8W,T3D,QV3D,QC3D,RHO3D, & !p8W in [PA] GLW,GSW,EMISS,CHKLOWQ, CHS, & FLQC,FLHC,MAVAIL,CANWAT,VEGFRA,ALB,ZNT, & Z0,SNOALB,ALBBCK,LAI, & !new mminlu, landusef, nlcat, mosaic_lu, & mosaic_soil, soilctop, nscat, & !new QSFC,QSG,QVG,QCG,DEW,SOILT1,TSNAV, & TBOT,IVGTYP,ISLTYP,XLAND, & ISWATER,ISICE,XICE,XICE_THRESHOLD, & CP,ROVCP,G0,LV,STBOLT, & SOILMOIS,SH2O,SMAVAIL,SMMAX, & TSO,SOILT,HFX,QFX,LH, & SFCRUNOFF,UDRUNOFF,ACRUNOFF,SFCEXC, & SFCEVP,GRDFLX,SNOWFALLAC,ACSNOW,SNOM, & SMFR3D,KEEPFR3DFLAG, & myjpbl,shdmin,shdmax,rdlai2d, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) !----------------------------------------------------------------- IMPLICIT NONE !----------------------------------------------------------------- ! ! The RUC LSM model is described in: ! Smirnova, T.G., J.M. Brown, and S.G. Benjamin, 1997: ! Performance of different soil model configurations in simulating ! ground surface temperature and surface fluxes. ! Mon. Wea. Rev. 125, 1870-1884. ! Smirnova, T.G., J.M. Brown, and D. Kim, 2000: Parameterization of ! cold-season processes in the MAPS land-surface scheme. ! J. Geophys. Res. 105, 4077-4086. !----------------------------------------------------------------- !-- DT time step (second) ! ktau - number of time step ! NSL - number of soil layers ! NZS - number of levels in soil ! ZS - depth of soil levels (m) !-- RAINBL - accumulated rain in [mm] between the PBL calls !-- RAINNCV one time step grid scale precipitation (mm/step) ! SNOW - snow water equivalent [mm] ! FRAZFRAC - fraction of frozen precipitation !-- PRECIPFR (mm) - time step frozen precipitation !-- SNOWC flag indicating snow coverage (1 for snow cover) !-- Z3D heights (m) !-- P8W 3D pressure (Pa) !-- T3D temperature (K) !-- QV3D 3D water vapor mixing ratio (Kg/Kg) ! QC3D - 3D cloud water mixing ratio (Kg/Kg) ! RHO3D - 3D air density (kg/m^3) !-- GLW downward long wave flux at ground surface (W/m^2) !-- GSW absorbed short wave flux at ground surface (W/m^2) !-- EMISS surface emissivity (between 0 and 1) ! FLQC - surface exchange coefficient for moisture (kg/m^2/s) ! FLHC - surface exchange coefficient for heat [W/m^2/s/degreeK] ! SFCEXC - surface exchange coefficient for heat [m/s] ! CANWAT - CANOPY MOISTURE CONTENT (mm) ! VEGFRA - vegetation fraction (between 0 and 100) ! ALB - surface albedo (between 0 and 1) ! SNOALB - maximum snow albedo (between 0 and 1) ! ALBBCK - snow-free albedo (between 0 and 1) ! ZNT - roughness length [m] !-- TBOT soil temperature at lower boundary (K) ! IVGTYP - USGS vegetation type (24 classes) ! ISLTYP - STASGO soil type (16 classes) !-- XLAND land mask (1 for land, 2 for water) !-- CP heat capacity at constant pressure for dry air (J/kg/K) !-- G0 acceleration due to gravity (m/s^2) !-- LV latent heat of melting (J/kg) !-- STBOLT Stefan-Boltzmann constant (W/m^2/K^4) ! SOILMOIS - soil moisture content (volumetric fraction) ! TSO - soil temp (K) !-- SOILT surface temperature (K) !-- HFX upward heat flux at the surface (W/m^2) !-- QFX upward moisture flux at the surface (kg/m^2/s) !-- LH upward latent heat flux (W/m^2) ! SFCRUNOFF - ground surface runoff [mm] ! UDRUNOFF - underground runoff [mm] ! ACRUNOFF - run-total surface runoff [mm] ! SFCEVP - total evaporation in [kg/m^2] ! GRDFLX - soil heat flux (W/m^2: negative, if downward from surface) ! SNOWFALLAC - run-total snowfall accumulation [m] ! ACSNOW - run-toral SWE of snowfall [mm] !-- CHKLOWQ - is either 0 or 1 (so far set equal to 1). !-- used only in MYJPBL. !-- tice - sea ice temperture (C) !-- rhosice - sea ice density (kg m^-3) !-- capice - sea ice volumetric heat capacity (J/m^3/K) !-- thdifice - sea ice thermal diffusivity (m^2/s) !-- !-- ims start index for i in memory !-- ime end index for i in memory !-- jms start index for j in memory !-- jme end index for j in memory !-- kms start index for k in memory !-- kme end index for k in memory !------------------------------------------------------------------------- ! INTEGER, PARAMETER :: nzss=5 ! INTEGER, PARAMETER :: nddzs=2*(nzss-2) INTEGER, PARAMETER :: nvegclas=24+3 REAL, INTENT(IN ) :: DT LOGICAL, INTENT(IN ) :: myjpbl,frpcpn INTEGER, INTENT(IN ) :: spp_lsm INTEGER, INTENT(IN ) :: NLCAT, NSCAT, mosaic_lu, mosaic_soil INTEGER, INTENT(IN ) :: ktau, nsl, isice, iswater, & ims,ime, jms,jme, kms,kme, & ids,ide, jds,jde, kds,kde, & its,ite, jts,jte, kts,kte #if (EM_CORE==1) REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),OPTIONAL:: pattern_spp_lsm REAL, DIMENSION( ims:ime, kms:kme, jms:jme ),OPTIONAL:: field_sf #endif REAL, DIMENSION( ims:ime, 1 :nsl, jms:jme ) :: field_sf_loc REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , & INTENT(IN ) :: QV3D, & QC3D, & p8w, & rho3D, & T3D, & z3D REAL, DIMENSION( ims:ime , jms:jme ), & INTENT(IN ) :: RAINBL, & GLW, & GSW, & ALBBCK, & CHS , & XICE, & XLAND, & VEGFRA, & TBOT !beka REAL, DIMENSION( ims:ime , jms:jme ), & INTENT(INOUT ) :: FLQC,FLHC #if (EM_CORE==1) REAL, OPTIONAL, DIMENSION( ims:ime , jms:jme ), & INTENT(IN ) :: GRAUPELNCV, & SNOWNCV, & RAINNCV REAL, DIMENSION( ims:ime , jms:jme ), & INTENT(IN ) :: lakemask INTEGER, INTENT(IN ) :: LakeModel #endif REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN ):: SHDMAX REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN ):: SHDMIN LOGICAL, intent(in) :: rdlai2d REAL, DIMENSION( 1:nsl), INTENT(IN ) :: ZS REAL, DIMENSION( ims:ime , jms:jme ), & INTENT(INOUT) :: & SNOW, & SNOWH, & SNOWC, & CANWAT, & ! new SNOALB, & ALB, & EMISS, & LAI, & MAVAIL, & SFCEXC, & Z0 , & ZNT REAL, DIMENSION( ims:ime , jms:jme ), & INTENT(IN ) :: & FRZFRAC INTEGER, DIMENSION( ims:ime , jms:jme ), & INTENT(IN ) :: IVGTYP, & ISLTYP CHARACTER(LEN=*), INTENT(IN ) :: MMINLU REAL, DIMENSION( ims:ime , 1:nlcat, jms:jme ), INTENT(IN):: LANDUSEF REAL, DIMENSION( ims:ime , 1:nscat, jms:jme ), INTENT(IN):: SOILCTOP REAL, INTENT(IN ) :: CP,ROVCP,G0,LV,STBOLT,XICE_threshold REAL, DIMENSION( ims:ime , 1:nsl, jms:jme ) , & INTENT(INOUT) :: SOILMOIS,SH2O,TSO REAL, DIMENSION( ims:ime, jms:jme ) , & INTENT(INOUT) :: SOILT, & HFX, & QFX, & LH, & SFCEVP, & SFCRUNOFF, & UDRUNOFF, & ACRUNOFF, & GRDFLX, & ACSNOW, & SNOM, & QVG, & QCG, & DEW, & QSFC, & QSG, & CHKLOWQ, & SOILT1, & TSNAV REAL, DIMENSION( ims:ime, jms:jme ) , & INTENT(INOUT) :: SMAVAIL, & SMMAX REAL, DIMENSION( its:ite, jts:jte ) :: & PC, & RUNOFF1, & RUNOFF2, & EMISSL, & ZNTL, & LMAVAIL, & SMELT, & SNOH, & SNFLX, & EDIR, & EC, & ETT, & SUBLIM, & sflx, & smf, & EVAPL, & PRCPL, & SEAICE, & INFILTR ! Energy and water budget variables: REAL, DIMENSION( its:ite, jts:jte ) :: & budget, & acbudget, & waterbudget, & acwaterbudget, & smtotold, & snowold, & canwatold REAL, DIMENSION( ims:ime, 1:nsl, jms:jme) & :: KEEPFR3DFLAG, & SMFR3D REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: & RHOSNF, & !RHO of snowfall PRECIPFR, & ! time-step frozen precip SNOWFALLAC !--- soil/snow properties REAL & :: RHOCS, & RHONEWSN, & RHOSN, & RHOSNFALL, & BCLH, & DQM, & KSAT, & PSIS, & QMIN, & QWRTZ, & REF, & WILT, & CANWATR, & SNOWFRAC, & SNHEI, & SNWE REAL :: CN, & SAT,CW, & C1SN, & C2SN, & KQWRTZ, & KICE, & KWT REAL, DIMENSION(1:NSL) :: ZSMAIN, & ZSHALF, & DTDZS2 REAL, DIMENSION(1:2*(nsl-2)) :: DTDZS REAL, DIMENSION(1:5001) :: TBQ REAL, DIMENSION( 1:nsl ) :: SOILM1D, & TSO1D, & SOILICE, & SOILIQW, & SMFRKEEP REAL, DIMENSION( 1:nsl ) :: KEEPFR REAL, DIMENSION( 1:nlcat ) :: lufrac REAL, DIMENSION( 1:nscat ) :: soilfrac REAL :: RSM, & SNWEPRINT, & SNHEIPRINT REAL :: PRCPMS, & NEWSNMS, & prcpncliq, & prcpncfr, & prcpculiq, & prcpcufr, & PATM, & PATMB, & TABS, & QVATM, & QCATM, & Q2SAT, & CONFLX, & RHO, & QKMS, & TKMS, & snowrat, & grauprat, & graupamt, & icerat, & curat, & INFILTRP REAL :: cq,r61,r273,arp,brp,x,evs,eis REAL :: cropsm REAL :: meltfactor, ac,as, wb INTEGER :: NROOT INTEGER :: ILAND,ISOIL,IFOREST INTEGER :: I,J,K,NZS,NZS1,NDDZS INTEGER :: k1,l,k2,kp,km CHARACTER (LEN=132) :: message REAL,DIMENSION(ims:ime,1:nsl,jms:jme) :: rstoch LOGICAL :: print_flag print_flag=wrf_at_debug_level(LSMRUC_DBG_LVL) !----------------------------------------------------------------- NZS=NSL NDDZS=2*(nzs-2) rstoch(ims:ime,1:nsl,jms:jme)=0.0 field_sf_loc(ims:ime,1:nsl,jms:jme)=0.0 !beka added #if (EM_CORE==1) if (spp_lsm==1) then do J=jts,jte do i=its,ite do k=1,nsl rstoch(i,k,j) = pattern_spp_lsm(i,k,j) field_sf_loc(i,k,j)=field_sf(i,k,j) enddo enddo enddo if(ktau.eq.10 .and. ktau.eq.20 .and. ktau.eq.30) then print *,'ktau, min/max rstoch',ktau,minval(rstoch(:,1,:)),maxval(rstoch(:,1,:)) endif endif #endif !---- table TBQ is for resolution of balance equation in VILKA CQ=173.15-.05 R273=1./273.15 R61=6.1153*0.62198 ARP=77455.*41.9/461.525 BRP=64.*41.9/461.525 DO K=1,5001 CQ=CQ+.05 ! TBQ(K)=R61*EXP(ARP*(R273-1./CQ)-BRP*LOG(CQ*R273)) EVS=EXP(17.67*(CQ-273.15)/(CQ-29.65)) EIS=EXP(22.514-6.15E3/CQ) if(CQ.ge.273.15) then ! tbq is in mb tbq(k) = R61*evs else tbq(k) = R61*eis endif END DO !--- Initialize soil/vegetation parameters !--- This is temporary until SI is added to mass coordinate ---!!!!! #if ( NMM_CORE == 1 ) if(ktau+1.eq.1) then #else if(ktau.eq.1) then #endif DO J=jts,jte DO i=its,ite do k=1,nsl keepfr3dflag(i,k,j)=0. enddo !--- initializing snow fraction (thereshold = 32 mm of snow water ! or ~100 mm of snow height) if it is not defined yet if(snow(i,j) > 0. .and. snowc(i,j) <= 0.) then snowc(i,j) = min(1.,snow(i,j)/32.) endif !--- initializing inside snow temp if it is not defined IF((soilt1(i,j) .LT. 170.) .or. (soilt1(i,j) .GT.400.)) THEN IF(snow(i,j).gt.32.) THEN soilt1(i,j)=0.5*(soilt(i,j)+tso(i,1,j)) IF ( print_flag ) THEN WRITE ( message , FMT='(A,F8.3,2I6)' ) & 'Temperature inside snow is initialized in RUCLSM ', soilt1(i,j),i,j CALL wrf_debug ( 0 , message ) ENDIF ELSE soilt1(i,j) = tso(i,1,j) ENDIF ENDIF tsnav(i,j) =0.5*(soilt(i,j)+tso(i,1,j))-273.15 patmb=P8w(i,kms,j)*1.e-2 QSG (i,j) = QSN(SOILT(i,j),TBQ)/PATMB IF((qvg(i,j) .LE. 0.) .or. (qvg(i,j) .GT.0.1)) THEN !17sept19 - bad approximation with very low mavail. !qvg(i,j) = QSG(i,j)*mavail(i,j) qvg (i,j) = qv3d(i,1,j) qcg (i,j) =0. IF ( print_flag ) THEN WRITE ( message , FMT='(A,3F8.3,2I6)' ) & 'QVG is initialized in RUCLSM ', qvg(i,j),mavail(i,j),qsg(i,j),i,j CALL wrf_debug ( 0 , message ) ENDIF ENDIF qsfc(i,j) = qvg(i,j)/(1.+qvg(i,j)) SMELT(i,j) = 0. SNOM (i,j) = 0. ! ACSNOW(i,j) = 0. SNOWFALLAC(i,j) = 0. PRECIPFR(i,j) = 0. RHOSNF(i,j) = -1.e3 ! non-zero flag SNFLX(i,j) = 0. DEW (i,j) = 0. PC (i,j) = 0. zntl (i,j) = 0. RUNOFF1(i,j) = 0. RUNOFF2(i,j) = 0. SFCRUNOFF(i,j) = 0. UDRUNOFF(i,j) = 0. ACRUNOFF(i,j) = 0. emissl (i,j) = 0. budget(i,j) = 0. acbudget(i,j) = 0. waterbudget(i,j) = 0. acwaterbudget(i,j) = 0. smtotold(i,j)=0. canwatold(i,j)=0. ! Temporarily!!! ! canwat(i,j)=0. ! For RUC LSM CHKLOWQ needed for MYJPBL should ! 1 because is actual specific humidity at the surface, and ! not the saturation value chklowq(i,j) = 1. infiltr(i,j) = 0. snoh (i,j) = 0. edir (i,j) = 0. ec (i,j) = 0. ett (i,j) = 0. sublim(i,j) = 0. sflx (i,j) = 0. smf (i,j) = 0. evapl (i,j) = 0. prcpl (i,j) = 0. ENDDO ENDDO do k=1,nsl soilice(k)=0. soiliqw(k)=0. enddo endif !----------------------------------------------------------------- PRCPMS = 0. newsnms = 0. prcpncliq = 0. prcpculiq = 0. prcpncfr = 0. prcpcufr = 0. DO J=jts,jte DO i=its,ite IF ( print_flag ) THEN print *,' IN LSMRUC ','ims,ime,jms,jme,its,ite,jts,jte,nzs', & ims,ime,jms,jme,its,ite,jts,jte,nzs print *,' IVGTYP, ISLTYP ', ivgtyp(i,j),isltyp(i,j) print *,' MAVAIL ', mavail(i,j) print *,' SOILT,QVG,P8w',soilt(i,j),qvg(i,j),p8w(i,1,j) print *, 'LSMRUC, I,J,xland, QFX,HFX from SFCLAY',i,j,xland(i,j), & qfx(i,j),hfx(i,j) print *, ' GSW, GLW =',gsw(i,j),glw(i,j) print *, 'SOILT, TSO start of time step =',soilt(i,j),(tso(i,k,j),k=1,nsl) print *, 'SOILMOIS start of time step =',(soilmois(i,k,j),k=1,nsl) print *, 'SMFROZEN start of time step =',(smfr3d(i,k,j),k=1,nsl) print *, ' I,J=, after SFCLAY CHS,FLHC ',i,j,chs(i,j),flhc(i,j) print *, 'LSMRUC, IVGTYP,ISLTYP,ALB = ', ivgtyp(i,j),isltyp(i,j),alb(i,j),i,j print *, 'LSMRUC I,J,DT,RAINBL =',I,J,dt,RAINBL(i,j) print *, 'XLAND ---->, ivgtype,isoiltyp,i,j',xland(i,j),ivgtyp(i,j),isltyp(i,j),i,j ENDIF ILAND = IVGTYP(i,j) ISOIL = ISLTYP(I,J) TABS = T3D(i,kms,j) QVATM = QV3D(i,kms,j) QCATM = QC3D(i,kms,j) PATM = P8w(i,kms,j)*1.e-5 !-- Z3D(1) is thickness between first full sigma level and the surface, !-- but first mass level is at the half of the first sigma level !-- (u and v are also at the half of first sigma level) CONFLX = Z3D(i,kms,j)*0.5 RHO = RHO3D(I,kms,J) ! -- initialize snow, graupel and ice fractions in frozen precip snowrat = 0. grauprat = 0. icerat = 0. curat = 0. IF(FRPCPN) THEN #if (EM_CORE==1) prcpncliq = rainncv(i,j)*(1.-frzfrac(i,j)) prcpncfr = rainncv(i,j)*frzfrac(i,j) !- apply the same frozen precipitation fraction to convective precip !tgs - 31 mar17 - add safety temperature check in case Thompson MP produces ! frozen precip at T > 273. if(frzfrac(i,j) > 0..and. tabs < 273.) then prcpculiq = max(0.,(rainbl(i,j)-rainncv(i,j))*(1.-frzfrac(i,j))) prcpcufr = max(0.,(rainbl(i,j)-rainncv(i,j))*frzfrac(i,j)) else if(tabs < 273.) then prcpcufr = max(0.,(rainbl(i,j)-rainncv(i,j))) prcpculiq = 0. else prcpcufr = 0. prcpculiq = max(0.,(rainbl(i,j)-rainncv(i,j))) endif ! tabs < 273. endif ! frzfrac > 0. !--- 1*e-3 is to convert from mm/s to m/s PRCPMS = (prcpncliq + prcpculiq)/DT*1.e-3 NEWSNMS = (prcpncfr + prcpcufr)/DT*1.e-3 IF ( PRESENT( graupelncv ) ) THEN graupamt = graupelncv(i,j) ELSE graupamt = 0. ENDIF if((prcpncfr + prcpcufr) > 0.) then ! -- calculate snow, graupel and ice fractions in falling frozen precip snowrat=min(1.,max(0.,snowncv(i,j)/(prcpncfr + prcpcufr))) grauprat=min(1.,max(0.,graupamt/(prcpncfr + prcpcufr))) icerat=min(1.,max(0.,(prcpncfr-snowncv(i,j)-graupamt) & /(prcpncfr + prcpcufr))) curat=min(1.,max(0.,(prcpcufr/(prcpncfr + prcpcufr)))) endif #else PRCPMS = (RAINBL(i,j)/DT*1.e-3)*(1-FRZFRAC(I,J)) NEWSNMS = (RAINBL(i,j)/DT*1.e-3)*FRZFRAC(I,J) if(newsnms == 0.) then snowrat = 0. else snowrat = min(1.,newsnms/(newsnms+prcpms)) endif #endif ELSE ! .not. FRPCPN if (tabs.le.273.15) then PRCPMS = 0. NEWSNMS = RAINBL(i,j)/DT*1.e-3 !-- here no info about constituents of frozen precipitation, !-- suppose it is all snow snowrat = 1. else PRCPMS = RAINBL(i,j)/DT*1.e-3 NEWSNMS = 0. endif ENDIF ! -- save time-step water equivalent of frozen precipitation in PRECIPFR array to be used in ! module_diagnostics precipfr(i,j) = NEWSNMS * DT *1.e3 !--- convert exchange coeff QKMS to [m/s] QKMS=FLQC(I,J)/RHO/MAVAIL(I,J) ! TKMS=FLHC(I,J)/RHO/CP TKMS=FLHC(I,J)/RHO/(CP*(1.+0.84*QVATM)) ! mynnsfc uses CPM !--- convert incoming snow and canwat from mm to m SNWE=SNOW(I,J)*1.E-3 SNHEI=SNOWH(I,J) CANWATR=CANWAT(I,J)*1.E-3 SNOWFRAC=SNOWC(I,J) RHOSNFALL=RHOSNF(I,J) snowold(i,j)=snwe !----- zsmain(1)=0. zshalf(1)=0. do k=2,nzs zsmain(k)= zs(k) zshalf(k)=0.5*(zsmain(k-1) + zsmain(k)) enddo do k=1,nlcat lufrac(k) = landusef(i,k,j) enddo do k=1,nscat soilfrac(k) = soilctop(i,k,j) enddo !------------------------------------------------------------ !----- DDZS and DSDZ1 are for implicit solution of soil eqns. !------------------------------------------------------------- NZS1=NZS-1 !----- IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,' DT,NZS1, ZSMAIN, ZSHALF --->', dt,nzs1,zsmain,zshalf ENDIF DO K=2,NZS1 K1=2*K-3 K2=K1+1 X=DT/2./(ZSHALF(K+1)-ZSHALF(K)) DTDZS(K1)=X/(ZSMAIN(K)-ZSMAIN(K-1)) DTDZS2(K-1)=X DTDZS(K2)=X/(ZSMAIN(K+1)-ZSMAIN(K)) END DO !27jul2011 - CN and SAT are defined in VEGPARM.TBL ! CN=0.5 ! exponent ! SAT=0.0004 ! canopy water saturated CW =4.183E6 !--- Constants used in Johansen soil thermal !--- conductivity method KQWRTZ=7.7 KICE=2.2 KWT=0.57 !*********************************************************************** !--- Constants for snow density calculations C1SN and C2SN c1sn=0.026 ! c1sn=0.01 c2sn=21. !*********************************************************************** NROOT= 4 ! ! rooting depth RHONEWSN = 200. if(SNOW(i,j).gt.0. .and. SNOWH(i,j).gt.0.) then RHOSN = SNOW(i,j)/SNOWH(i,j) else RHOSN = 300. endif IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN if(ktau.eq.1 .and.(i.eq.358.and.j.eq.260)) & print *,'before SOILVEGIN - z0,znt(195,254)',z0(i,j),znt(i,j) ENDIF !--- initializing soil and surface properties CALL SOILVEGIN ( mosaic_lu, mosaic_soil,soilfrac,nscat,shdmin(i,j),shdmax(i,j),& NLCAT,ILAND,ISOIL,iswater,IFOREST,lufrac,VEGFRA(I,J), & EMISSL(I,J),PC(I,J),ZNT(I,J),LAI(I,J),RDLAI2D, & QWRTZ,RHOCS,BCLH,DQM,KSAT,PSIS,QMIN,REF,WILT,i,j ) IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN if(ktau.eq.1 .and.(i.eq.358.and.j.eq.260)) & print *,'after SOILVEGIN - z0,znt(375,254),lai(375,254)',z0(i,j),znt(i,j),lai(i,j) if(ktau.eq.1 .and. (i.eq.358.and.j.eq.260)) then print *,'NLCAT,iland,lufrac,EMISSL(I,J),PC(I,J),ZNT(I,J),LAI(I,J)', & NLCAT,iland,lufrac,EMISSL(I,J),PC(I,J),ZNT(I,J),LAI(I,J),i,j print *,'NSCAT,soilfrac,QWRTZ,RHOCS,BCLH,DQM,KSAT,PSIS,QMIN,REF,WILT',& NSCAT,soilfrac,QWRTZ,RHOCS,BCLH,DQM,KSAT,PSIS,QMIN,REF,WILT,i,j endif ENDIF CN=CFACTR_DATA ! exponent ! SAT=max(1.e-5,(min(5.e-4,(CMCMAX_DATA * (1.-exp(-0.5*lai(i,j))) * 0.01*VEGFRA(I,J))))) ! canopy water saturated SAT = 5.e-4 ! units [m] ! if(i==666.and.j==282) print *,'second 666,282 - sat',sat !-- definition of number of soil levels in the rooting zone ! IF(iforest(ivgtyp(i,j)).ne.1) THEN IF(iforest.gt.2) THEN !---- all vegetation types except evergreen and mixed forests !18apr08 - define meltfactor for Egglston melting limit: ! for open areas factor is 2, and for forests - factor is 0.85 ! This will make limit on snow melting smaller and let snow stay ! longer in the forests. meltfactor = 2.0 do k=2,nzs if(zsmain(k).ge.0.4) then NROOT=K goto 111 endif enddo ELSE !---- evergreen and mixed forests !18apr08 - define meltfactor ! meltfactor = 1.5 ! 28 March 11 - Previously used value of metfactor= 1.5 needs to be further reduced ! to compensate for low snow albedos in the forested areas. ! Melting rate in forests will reduce. meltfactor = 0.85 do k=2,nzs if(zsmain(k).ge.1.1) then NROOT=K goto 111 endif enddo ENDIF 111 continue !----- IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,' ZNT, LAI, VEGFRA, SAT, EMIS, PC --->', & ZNT(I,J),LAI(I,J),VEGFRA(I,J),SAT,EMISSL(I,J),PC(I,J) print *,' ZS, ZSMAIN, ZSHALF, CONFLX, CN, SAT, --->', zs,zsmain,zshalf,conflx,cn,sat print *,'NROOT, meltfactor, iforest, ivgtyp, i,j ', nroot,meltfactor,iforest,ivgtyp(I,J),I,J ! print *,'NROOT, iforest, ivgtyp, i,j ', nroot,iforest(ivgtyp(i,j)),ivgtyp(I,J),I,J ENDIF !!*** SET ZERO-VALUE FOR SOME OUTPUT DIAGNOSTIC ARRAYS ! if(i.eq.397.and.j.eq.562) then ! print *,'RUC LSM - xland(i,j),xice(i,j),snow(i,j)',i,j,xland(i,j),xice(i,j),snow(i,j) ! endif #if (EM_CORE==1) if(lakemodel==1. .and. lakemask(i,j)==1.) goto 2999 !Lakes #endif IF((XLAND(I,J)-1.5).GE.0.)THEN !-- Water SMAVAIL(I,J)=1.0 SMMAX(I,J)=1.0 SNOW(I,J)=0.0 SNOWH(I,J)=0.0 SNOWC(I,J)=0.0 LMAVAIL(I,J)=1.0 ! accumulated water equivalent of frozen precipitation over water [mm] ! acsnow(i,j)=acsnow(i,j)+precipfr(i,j) ILAND=iswater ! ILAND=16 ISOIL=14 patmb=P8w(i,1,j)*1.e-2 qvg (i,j) = QSN(SOILT(i,j),TBQ)/PATMB qsfc(i,j) = qvg(i,j)/(1.+qvg(i,j)) CHKLOWQ(I,J)=1. Q2SAT=QSN(TABS,TBQ)/PATMB DO K=1,NZS SOILMOIS(I,K,J)=1.0 SH2O (I,K,J)=1.0 TSO(I,K,J)= SOILT(I,J) ENDDO IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN PRINT*,' water point, I=',I, & 'J=',J, 'SOILT=', SOILT(i,j) ENDIF ELSE ! LAND POINT OR SEA ICE if(xice(i,j).ge.xice_threshold) then ! if(IVGTYP(i,j).eq.isice) then SEAICE(i,j)=1. else SEAICE(i,j)=0. endif IF(SEAICE(I,J).GT.0.5)THEN !-- Sea-ice case IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN PRINT*,' sea-ice at water point, I=',I, & 'J=',J ENDIF ! ILAND = 24 ILAND = isice ISOIL = 16 ZNT(I,J) = 0.011 snoalb(i,j) = 0.75 dqm = 1. ref = 1. qmin = 0. wilt = 0. emissl(i,j) = 0.98 patmb=P8w(i,1,j)*1.e-2 qvg (i,j) = QSN(SOILT(i,j),TBQ)/PATMB qsg (i,j) = qvg(i,j) qsfc(i,j) = qvg(i,j)/(1.+qvg(i,j)) DO K=1,NZS soilmois(i,k,j) = 1. smfr3d(i,k,j) = 1. sh2o(i,k,j) = 0. keepfr3dflag(i,k,j) = 0. tso(i,k,j) = min(271.4,tso(i,k,j)) ENDDO ENDIF !beka !stochastic perturbation at time step 1 if (spp_lsm == 1) then IF (KTAU == 2) THEN do k=1,nzs SOILMOIS(i,k,j)=SOILMOIS(i,k,j)*(1+1.5*rstoch(i,k,j)) enddo ENDIF endif !beka ! Attention!!!! RUC LSM uses soil moisture content minus residual (minimum ! or dry soil moisture content for a given soil type) as a state variable. DO k=1,nzs ! soilm1d - soil moisture content minus residual [m**3/m**3] soilm1d (k) = min(max(0.,soilmois(i,k,j)-qmin),dqm) ! soilm1d (k) = min(max(0.,soilmois(i,k,j)),dqm) tso1d (k) = tso(i,k,j) soiliqw (k) = min(max(0.,sh2o(i,k,j)-qmin),soilm1d(k)) soilice (k) =(soilm1d (k) - soiliqw (k))/0.9 ENDDO do k=1,nzs smfrkeep(k) = smfr3d(i,k,j) keepfr (k) = keepfr3dflag(i,k,j) enddo LMAVAIL(I,J)=max(0.00001,min(1.,soilm1d(1)/(REF-QMIN))) ! LMAVAIL(I,J)=max(0.00001,min(1.,soilm1d(1)/dqm)) #if ( NMM_CORE == 1 ) if(ktau+1.gt.1) then #else if(ktau.gt.1) then #endif ! extract dew from the cloud water at the surface !30july13 QCG(I,J)=QCG(I,J)-DEW(I,J)/QKMS endif IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'LAND, i,j,tso1d,soilm1d,PATM,TABS,QVATM,QCATM,RHO', & i,j,tso1d,soilm1d,PATM,TABS,QVATM,QCATM,RHO print *,'CONFLX =',CONFLX print *,'SMFRKEEP,KEEPFR ',SMFRKEEP,KEEPFR ENDIF smtotold(i,j)=0. do k=1,nzs-1 smtotold(i,j)=smtotold(i,j)+(qmin+soilm1d(k))* & (zshalf(k+1)-zshalf(k)) enddo smtotold(i,j)=smtotold(i,j)+(qmin+soilm1d(nzs))* & (zsmain(nzs)-zshalf(nzs)) canwatold(i,j) = canwatr IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'before SFCTMP, spp_lsm, rstoch, field_sf_loc', & i,j,spp_lsm,(rstoch(i,k,j),k=1,nzs),(field_sf_loc(i,k,j),k=1,nzs) ENDIF !----------------------------------------------------------------- CALL SFCTMP (spp_lsm,rstoch(i,:,j),field_sf_loc(i,:,j), & dt,ktau,conflx,i,j, & !--- input variables nzs,nddzs,nroot,meltfactor, & !added meltfactor iland,isoil,xland(i,j),ivgtyp(i,j),isltyp(i,j), & PRCPMS, NEWSNMS,SNWE,SNHEI,SNOWFRAC, & RHOSN,RHONEWSN,RHOSNFALL, & snowrat,grauprat,icerat,curat, & PATM,TABS,QVATM,QCATM,RHO, & GLW(I,J),GSW(I,J),EMISSL(I,J), & QKMS,TKMS,PC(I,J),LMAVAIL(I,J), & canwatr,vegfra(I,J),alb(I,J),znt(I,J), & snoalb(i,j),albbck(i,j),lai(i,j), & !new myjpbl,seaice(i,j),isice, & !--- soil fixed fields QWRTZ, & rhocs,dqm,qmin,ref, & wilt,psis,bclh,ksat, & sat,cn,zsmain,zshalf,DTDZS,DTDZS2,tbq, & !--- constants cp,rovcp,g0,lv,stbolt,cw,c1sn,c2sn, & KQWRTZ,KICE,KWT, & !--- output variables snweprint,snheiprint,rsm, & soilm1d,tso1d,smfrkeep,keepfr, & soilt(I,J),soilt1(i,j),tsnav(i,j),dew(I,J), & qvg(I,J),qsg(I,J),qcg(I,J),SMELT(I,J), & SNOH(I,J),SNFLX(I,J),SNOM(I,J),SNOWFALLAC(I,J), & ACSNOW(I,J),edir(I,J),ec(I,J),ett(I,J),qfx(I,J), & lh(I,J),hfx(I,J),sflx(I,J),sublim(I,J), & evapl(I,J),prcpl(I,J),budget(i,j),runoff1(i,j), & runoff2(I,J),soilice,soiliqw,infiltrp,smf(i,j)) !----------------------------------------------------------------- ! Fraction of cropland category in the grid box should not have soil moisture below ! wilting point during the growing season. ! Let's keep soil moisture 20% above wilting point for the fraction of grid box under ! croplands. ! This change violates LSM moisture budget, but ! can be considered as a compensation for irrigation not included into LSM. IF (lufrac(crop) > 0 .and. lai(i,j) > 1.1) THEN ! IF (ivgtyp(i,j) == crop .and. lai(i,j) > 1.1) THEN ! cropland do k=1,nroot cropsm=1.1*wilt - qmin if(soilm1d(k) < cropsm*lufrac(crop)) then IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print * ,'Soil moisture is below wilting in cropland category at time step',ktau & ,'i,j,lufrac(crop),k,soilm1d(k),wilt,cropsm', & i,j,lufrac(crop),k,soilm1d(k),wilt,cropsm ENDIF soilm1d(k) = cropsm*lufrac(crop) IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print * ,'Added soil water to cropland category, i,j,k,soilm1d(k)',i,j,k,soilm1d(k) ENDIF endif enddo ELSEIF (ivgtyp(i,j) == natural .and. lai(i,j) > 0.7) THEN ! grassland: assume that 40% of grassland is irrigated cropland do k=1,nroot cropsm=1.2*wilt - qmin if(soilm1d(k) < cropsm*lufrac(natural)*0.4) then IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print * ,'Soil moisture is below wilting in mixed grassland/cropland category at time step',ktau & ,'i,j,lufrac(natural),k,soilm1d(k),wilt', & i,j,lufrac(natural),k,soilm1d(k),wilt ENDIF soilm1d(k) = cropsm * lufrac(natural)*0.4 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print * ,'Added soil water to grassland category, i,j,k,soilm1d(k)',i,j,k,soilm1d(k) ENDIF endif enddo ENDIF ! Fill in field_sf to pass perturbed field of hydraulic cond. up to model driver and output #if (EM_CORE==1) if (spp_lsm==1) then do k=1,nsl field_sf(i,k,j)=field_sf_loc(i,k,j) enddo endif #endif !*** DIAGNOSTICS !--- available and maximum soil moisture content in the soil !--- domain smavail(i,j) = 0. smmax (i,j) = 0. do k=1,nzs-1 smavail(i,j)=smavail(i,j)+(qmin+soilm1d(k))* & (zshalf(k+1)-zshalf(k)) smmax (i,j) =smmax (i,j)+(qmin+dqm)* & (zshalf(k+1)-zshalf(k)) enddo smavail(i,j)=smavail(i,j)+(qmin+soilm1d(nzs))* & (zsmain(nzs)-zshalf(nzs)) smmax (i,j) =smmax (i,j)+(qmin+dqm)* & (zsmain(nzs)-zshalf(nzs)) !--- Convert the water unit into mm SFCRUNOFF(I,J) = SFCRUNOFF(I,J)+RUNOFF1(I,J)*DT*1000.0 UDRUNOFF (I,J) = UDRUNOFF(I,J)+RUNOFF2(I,J)*DT*1000.0 ACRUNOFF(I,J) = ACRUNOFF(I,J)+RUNOFF1(I,J)*DT*1000.0 SMAVAIL (I,J) = SMAVAIL(I,J) * 1000. SMMAX (I,J) = SMMAX(I,J) * 1000. smtotold (I,J) = smtotold(I,J) * 1000. do k=1,nzs ! soilmois(i,k,j) = soilm1d(k) soilmois(i,k,j) = soilm1d(k) + qmin sh2o (i,k,j) = min(soiliqw(k) + qmin,soilmois(i,k,j)) tso(i,k,j) = tso1d(k) enddo tso(i,nzs,j) = tbot(i,j) do k=1,nzs smfr3d(i,k,j) = smfrkeep(k) keepfr3dflag(i,k,j) = keepfr (k) enddo !tgs add together dew and cloud at the ground surface !30july13 qcg(i,j)=qcg(i,j)+dew(i,j)/qkms Z0 (I,J) = ZNT (I,J) SFCEXC (I,J) = TKMS patmb=P8w(i,1,j)*1.e-2 Q2SAT=QSN(TABS,TBQ)/PATMB QSFC(I,J) = QVG(I,J)/(1.+QVG(I,J)) ! for MYJ PBL scheme IF((myjpbl).AND.(QVATM.GE.Q2SAT*0.95).AND.QVATM.LT.qvg(I,J))THEN CHKLOWQ(I,J)=0. ELSE CHKLOWQ(I,J)=1. ENDIF IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN if(CHKLOWQ(I,J).eq.0.) then print *,'i,j,CHKLOWQ', & i,j,CHKLOWQ(I,J) endif ENDIF if(snow(i,j)==0.) EMISSL(i,j) = LEMITBL(IVGTYP(i,j)) ! evan if (spp_lsm == 1) then EMISSL(i,j) = min(EMISSL(i,j) * (1. + 0.1*rstoch(i,1,j)), 1.) endif EMISS (I,J) = EMISSL(I,J) ! SNOW is in [mm], SNWE is in [m]; CANWAT is in mm, CANWATR is in m SNOW (i,j) = SNWE*1000. SNOWH (I,J) = SNHEI CANWAT (I,J) = CANWATR*1000. INFILTR(I,J) = INFILTRP MAVAIL (i,j) = LMAVAIL(I,J) IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,' LAND, I=,J=, QFX, HFX after SFCTMP', i,j,lh(i,j),hfx(i,j) ENDIF !!! QFX (I,J) = LH(I,J)/LV SFCEVP (I,J) = SFCEVP (I,J) + QFX (I,J) * DT GRDFLX (I,J) = -1. * sflx(I,J) ! if(smf(i,j) .ne.0.) then !tgs - SMF.NE.0. when there is phase change in the top soil layer ! The heat of soil water freezing/thawing is not computed explicitly ! and is responsible for the residual in the energy budget. ! print *,'Budget',budget(i,j),i,j,smf(i,j) ! endif !--- SNOWC snow cover flag if(snowfrac > 0. .and. xice(i,j).ge.xice_threshold ) then SNOWFRAC = SNOWFRAC*XICE(I,J) endif SNOWC(I,J)=SNOWFRAC !--- RHOSNF - density of snowfall RHOSNF(I,J)=RHOSNFALL ! Accumulated moisture flux [kg/m^2] SFCEVP (I,J) = SFCEVP (I,J) + QFX (I,J) * DT !TEST!!!! for test put heat budget term in GRDFLX ! acbudget(i,j)=acbudget(i,j)+budget(i,j)-smf(i,j) ! GRDFLX (I,J) = acbudget(i,j) ! if(smf(i,j) .ne.0.) then !tgs - SMF.NE.0. when there is phase change in the top soil layer ! The heat of freezing/thawing of soil water is not computed explicitly ! and is responsible for the residual in the energy budget. ! endif ! budget(i,j)=budget(i,j)-smf(i,j) ac=0. as=0. ac=max(0.,canwat(i,j)-canwatold(i,j)) as=max(0.,snwe-snowold(i,j)) wb =rainbl(i,j)+smelt(i,j)*dt*1.e3 & ! source -qfx(i,j)*dt & -runoff1(i,j)*dt*1.e3-runoff2(i,j)*dt*1.e3 & -ac-as - (smavail(i,j)-smtotold(i,j)) waterbudget(i,j)=rainbl(i,j)+smelt(i,j)*dt*1.e3 & ! source -qfx(i,j)*dt & -runoff1(i,j)*dt*1.e3-runoff2(i,j)*dt*1.e3 & -ac-as - (smavail(i,j)-smtotold(i,j)) ! waterbudget(i,j)=rainbl(i,j)-qfx(i,j)*dt-(smavail(i,j)-smtotold(i,j)) & acwaterbudget(i,j)=acwaterbudget(i,j)+waterbudget(i,j) !!!!TEST use LH to check water budget ! GRDFLX (I,J) = waterbudget(i,j) IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'Smf=',smf(i,j),i,j print *,'Budget',budget(i,j),i,j print *,'RUNOFF2= ', i,j,runoff2(i,j) print *,'Water budget ', i,j,waterbudget(i,j) print *,'rainbl,qfx*dt,runoff1,smelt*dt*1.e3,smchange', & i,j,rainbl(i,j),qfx(i,j)*dt,runoff1(i,j)*dt*1.e3, & smelt(i,j)*dt*1.e3, & (smavail(i,j)-smtotold(i,j)) print *,'SNOW,SNOWold',i,j,snwe,snowold(i,j) print *,'SNOW-SNOWold',i,j,max(0.,snwe-snowold(i,j)) print *,'CANWATold, canwat ',i,j,canwatold(i,j),canwat(i,j) print *,'canwat(i,j)-canwatold(i,j)',max(0.,canwat(i,j)-canwatold(i,j)) ENDIF IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'LAND, i,j,tso1d,soilm1d,soilt - end of time step', & i,j,tso1d,soilm1d,soilt(i,j) print *,'LAND, QFX, HFX after SFCTMP', i,j,lh(i,j),hfx(i,j) ENDIF !--- end of a land or sea ice point ENDIF 2999 continue ! lakes ENDDO ENDDO !----------------------------------------------------------------- END SUBROUTINE LSMRUC !----------------------------------------------------------------- SUBROUTINE SFCTMP (spp_lsm,rstochcol,fieldcol_sf, & delt,ktau,conflx,i,j, & !--- input variables nzs,nddzs,nroot,meltfactor, & ILAND,ISOIL,XLAND,IVGTYP,ISLTYP,PRCPMS, & NEWSNMS,SNWE,SNHEI,SNOWFRAC, & RHOSN,RHONEWSN,RHOSNFALL, & snowrat,grauprat,icerat,curat, & PATM,TABS,QVATM,QCATM,rho, & GLW,GSW,EMISS,QKMS,TKMS,PC, & MAVAIL,CST,VEGFRA,ALB,ZNT, & ALB_SNOW,ALB_SNOW_FREE,lai, & MYJ,SEAICE,ISICE, & !--- soil fixed fields QWRTZ,rhocs,dqm,qmin,ref,wilt,psis,bclh,ksat, & sat,cn,zsmain,zshalf,DTDZS,DTDZS2,tbq, & !--- constants cp,rovcp,g0,lv,stbolt,cw,c1sn,c2sn, & KQWRTZ,KICE,KWT, & !--- output variables snweprint,snheiprint,rsm, & soilm1d,ts1d,smfrkeep,keepfr,soilt,soilt1, & tsnav,dew,qvg,qsg,qcg, & SMELT,SNOH,SNFLX,SNOM,SNOWFALLAC,ACSNOW, & edir1,ec1,ett1,eeta,qfx,hfx,s,sublim, & evapl,prcpl,fltot,runoff1,runoff2,soilice, & soiliqw,infiltr,smf) !----------------------------------------------------------------- IMPLICIT NONE !----------------------------------------------------------------- !--- input variables INTEGER, INTENT(IN ) :: isice,i,j,nroot,ktau,nzs , & nddzs !nddzs=2*(nzs-2) REAL, INTENT(IN ) :: DELT,CONFLX,meltfactor REAL, INTENT(IN ) :: C1SN,C2SN LOGICAL, INTENT(IN ) :: myj !--- 3-D Atmospheric variables REAL , & INTENT(IN ) :: PATM, & TABS, & QVATM, & QCATM REAL , & INTENT(IN ) :: GLW, & GSW, & PC, & VEGFRA, & ALB_SNOW_FREE, & lai, & SEAICE, & XLAND, & RHO, & QKMS, & TKMS INTEGER, INTENT(IN ) :: IVGTYP, ISLTYP !--- 2-D variables REAL , & INTENT(INOUT) :: EMISS, & MAVAIL, & SNOWFRAC, & ALB_SNOW, & ALB, & CST !--- soil properties REAL :: & RHOCS, & BCLH, & DQM, & KSAT, & PSIS, & QMIN, & QWRTZ, & REF, & SAT, & WILT REAL, INTENT(IN ) :: CN, & CW, & CP, & ROVCP, & G0, & LV, & STBOLT, & KQWRTZ, & KICE, & KWT REAL, DIMENSION(1:NZS), INTENT(IN) :: ZSMAIN, & ZSHALF, & DTDZS2 REAL, DIMENSION(1:NZS), INTENT(IN) :: rstochcol REAL, DIMENSION(1:NZS), INTENT(INOUT) :: fieldcol_sf REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS REAL, DIMENSION(1:5001), INTENT(IN) :: TBQ !--- input/output variables !-------- 3-d soil moisture and temperature REAL, DIMENSION( 1:nzs ) , & INTENT(INOUT) :: TS1D, & SOILM1D, & SMFRKEEP REAL, DIMENSION( 1:nzs ) , & INTENT(INOUT) :: KEEPFR REAL, DIMENSION(1:NZS), INTENT(INOUT) :: SOILICE, & SOILIQW INTEGER, INTENT(INOUT) :: ILAND,ISOIL INTEGER :: ILANDs !-------- 2-d variables REAL , & INTENT(INOUT) :: DEW, & EDIR1, & EC1, & ETT1, & EETA, & EVAPL, & INFILTR, & RHOSN, & RHONEWSN, & rhosnfall, & snowrat, & grauprat, & icerat, & curat, & SUBLIM, & PRCPL, & QVG, & QSG, & QCG, & QFX, & HFX, & fltot, & smf, & S, & RUNOFF1, & RUNOFF2, & ACSNOW, & SNOWFALLAC, & SNWE, & SNHEI, & SMELT, & SNOM, & SNOH, & SNFLX, & SOILT, & SOILT1, & TSNAV, & ZNT REAL, DIMENSION(1:NZS) :: & tice, & rhosice, & capice, & thdifice, & TS1DS, & SOILM1DS, & SMFRKEEPS, & SOILIQWS, & SOILICES, & KEEPFRS !-------- 1-d variables REAL :: & DEWS, & MAVAILS, & EDIR1s, & EC1s, & csts, & ETT1s, & EETAs, & EVAPLs, & INFILTRs, & PRCPLS, & QVGS, & QSGS, & QCGS, & QFXS, & HFXS, & fltots, & RUNOFF1S, & RUNOFF2s, & SS, & SOILTs REAL, INTENT(INOUT) :: RSM, & SNWEPRINT, & SNHEIPRINT INTEGER, INTENT(IN) :: spp_lsm !--- Local variables INTEGER :: K,ILNB REAL :: BSN, XSN , & RAINF, SNTH, NEWSN, PRCPMS, NEWSNMS , & T3, UPFLUX, XINET REAL :: snhei_crit, snhei_crit_newsn, keep_snow_albedo, SNOWFRACnewsn REAL :: newsnowratio, dd1 REAL :: rhonewgr,rhonewice REAL :: RNET,GSWNEW,GSWIN,EMISSN,ZNTSN,EMISS_snowfree REAL :: VEGFRAC, snow_mosaic, snfr, vgfr real :: cice, albice, albsn, drip, dripsn, dripliq real :: interw, intersn, infwater, intwratio !----------------------------------------------------------------- integer, parameter :: ilsnow=99 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,' in SFCTMP',i,j,nzs,nddzs,nroot, & SNWE,RHOSN,SNOM,SMELT,TS1D ENDIF snow_mosaic=0. snfr = 1. NEWSN=0. newsnowratio = 0. snowfracnewsn=0. if(snhei == 0.) snowfrac=0. smelt = 0. RAINF = 0. RSM=0. DD1=0. INFILTR=0. ! Jul 2016 - Avissar and Pielke (1989) ! This formulation depending on LAI defines relative contribution of the vegetation to ! the total heat fluxes between surface and atmosphere. ! With VEGFRA=100% and LAI=3, VEGFRAC=0.86 meaning that vegetation contributes ! only 86% of the total surface fluxes. ! VGFR=0.01*VEGFRA ! % --> fraction ! VEGFRAC=2.*lai*vgfr/(1.+2.*lai*vgfr) VEGFRAC=0.01*VEGFRA !beka if (spp_lsm == 1) then VEGFRAC = min(VEGFRAC * (1. + 0.33*rstochcol(1)), 1.) endif drip = 0. dripsn = 0. dripliq = 0. smf = 0. interw=0. intersn=0. infwater=0. !---initialize local arrays for sea ice do k=1,nzs tice(k) = 0. rhosice(k) = 0. cice = 0. capice(k) = 0. thdifice(k) = 0. enddo GSWnew=GSW ! evan if (spp_lsm == 1) then alb = min(ALBTBL(IVGTYP) * (1. + 0.4*rstochcol(1)),1.) endif GSWin=GSW/(1.-alb) ALBice=ALB_SNOW_FREE ALBsn=alb_snow EMISSN = 0.98 EMISS_snowfree = LEMITBL(IVGTYP) ! evan if (spp_lsm == 1) then EMISS_snowfree = min(EMISS_snowfree * (1. + 0.1*rstochcol(1)), 1.) endif !beka stochastic pert of soil mositure ! if(spp_lsm == 1) then !! stochastic perturbations for soil moisture ! do k=1,nzs ! soilm1d(k)=max(0.,min(soilm1d(k)*(1+rstochcol(k)),dqm)) ! enddo ! endif !--- sea ice properties !--- N.N Zubov "Arctic Ice" !--- no salinity dependence because we consider the ice pack !--- to be old and to have low salinity (0.0002) if(SEAICE.ge.0.5) then do k=1,nzs tice(k) = ts1d(k) - 273.15 rhosice(k) = 917.6/(1-0.000165*tice(k)) cice = 2115.85 +7.7948*tice(k) capice(k) = cice*rhosice(k) thdifice(k) = 2.260872/capice(k) enddo !-- SEA ICE ALB dependence on ice temperature. When ice temperature is !-- below critical value of -10C - no change to albedo. !-- If temperature is higher that -10C then albedo is decreasing. !-- The minimum albedo at t=0C for ice is 0.1 less. ALBice = MIN(ALB_SNOW_FREE,MAX(ALB_SNOW_FREE - 0.05, & ALB_SNOW_FREE - 0.1*(tice(1)+10.)/10. )) endif IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN ! print *,'I,J,KTAU,QKMS,TKMS', i,j,ktau,qkms,tkms print *,'alb_snow_free',ALB_SNOW_FREE print *,'GSW,GSWnew,GLW,SOILT,EMISS,ALB,ALBice,SNWE',& GSW,GSWnew,GLW,SOILT,EMISS,ALB,ALBice,SNWE ENDIF if(snhei.gt.0.0081*1.e3/rhosn) then !*** Update snow density for current temperature (Koren et al. 1999) BSN=delt/3600.*c1sn*exp(0.08*min(0.,tsnav)-c2sn*rhosn*1.e-3) if(bsn*snwe*100..lt.1.e-4) goto 777 XSN=rhosn*(exp(bsn*snwe*100.)-1.)/(bsn*snwe*100.) rhosn=MIN(MAX(58.8,XSN),500.) !13mar18 rhosn=MIN(MAX(76.9,XSN),500.) ! rhosn=MIN(MAX(62.5,XSN),890.) ! rhosn=MIN(MAX(100.,XSN),400.) ! rhosn=MIN(MAX(50.,XSN),400.) 777 continue ! else ! rhosn =200. ! rhonewsn =200. endif newsn=newsnms*delt !---- ACSNOW - run-total snowfall water [mm] ! acsnow=acsnow+newsn*1.e3 IF(NEWSN.GT.0.) THEN ! IF(NEWSN.GE.1.E-8) THEN IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *, 'THERE IS NEW SNOW, newsn', newsn ENDIF newsnowratio = min(1.,newsn/(snwe+newsn)) !*** Calculate fresh snow density (t > -15C, else MIN value) !*** Eq. 10 from Koren et al. (1999) !--- old formulation from Koren (1999) ! if(tabs.lt.258.15) then ! rhonewsn=50. ! rhonewsn=100. ! rhonewsn=62.5 ! else ! rhonewsn=MIN(rhonewsn,400.) ! endif !--- end of old formulation !--- 27 Feb 2014 - empirical formulations from John M. Brown ! rhonewsn=min(250.,rhowater/max(4.179,(13.*tanh((274.15-Tabs)*0.3333)))) !--- 13 Mar 2018 - formulation from Trevor Alcott rhonewsn=min(125.,1000.0/max(8.,(17.*tanh((276.65-Tabs)*0.15)))) rhonewgr=min(500.,rhowater/max(2.,(3.5*tanh((274.15-Tabs)*0.3333)))) rhonewice=rhonewsn !--- compute density of "snowfall" from weighted contribution ! of snow, graupel and ice fractions rhosnfall = min(500.,max(58.8,(rhonewsn*snowrat + & !13mar18 rhosnfall = min(500.,max(76.9,(rhonewsn*snowrat + & rhonewgr*grauprat + rhonewice*icerat + rhonewgr*curat))) ! from now on rhonewsn is the density of falling frozen precipitation rhonewsn=rhosnfall !*** Define average snow density of the snow pack considering !*** the amount of fresh snow (eq. 9 in Koren et al.(1999) !*** without snow melt ) xsn=(rhosn*snwe+rhonewsn*newsn)/ & (snwe+newsn) rhosn=MIN(MAX(58.8,XSN),500.) !13mar18 rhosn=MIN(MAX(76.9,XSN),500.) ! rhosn=MIN(MAX(100.,XSN),500.) ! rhosn=MIN(MAX(50.,XSN),400.) !Update snow on the ground ! snwe=snwe+newsn ! newsnowratio = min(1.,newsn/snwe) ! snhei=snwe*rhowater/rhosn ! NEWSN=NEWSN*rhowater/rhonewsn ENDIF ! end NEWSN > 0. IF(PRCPMS.NE.0.) THEN ! PRCPMS is liquid precipitation rate ! RAINF is a flag used for calculation of rain water ! heat content contribution into heat budget equation. Rain's temperature ! is set equal to air temperature at the first atmospheric ! level. RAINF=1. ENDIF drip = 0. intwratio=0. if(vegfrac > 0.01) then ! compute intercepted precipitation - Eq. 1 Lawrence et al., ! J. of Hydrometeorology, 2006, CLM. interw=0.25*DELT*PRCPMS*(1.-exp(-0.5*lai))*vegfrac intersn=0.25*NEWSN*(1.-exp(-0.5*lai))*vegfrac !original - next 2 lines ! interw=DELT*PRCPMS*vegfrac ! intersn=NEWSN*vegfrac infwater=PRCPMS - interw/delt if((interw+intersn) > 0.) then intwratio=interw/(interw+intersn) endif ! Update water/snow intercepted by the canopy dd1=CST + interw + intersn CST=DD1 ! if(i==666.and.j==282) print *,'666,282 - cst,sat,interw,intersn',cst,sat,interw,intersn IF(CST.GT.SAT) THEN CST=SAT DRIP=DD1-SAT ENDIF else CST=0. DRIP=0. interw=0. intersn=0. infwater=PRCPMS endif ! vegfrac > 0.01 ! SNHEI_CRIT is a threshold for fractional snow SNHEI_CRIT=0.01601*1.e3/rhosn SNHEI_CRIT_newsn=0.0005*1.e3/rhosn ! snowfrac from the previous time step SNOWFRAC=MIN(1.,SNHEI/(2.*SNHEI_CRIT)) if(snowfrac < 0.75) snow_mosaic = 1. IF(NEWSN.GT.0.) THEN !Update snow on the ground snwe=max(0.,snwe+newsn-intersn) ! Add drip to snow on the ground if(drip > 0.) then if (snow_mosaic==1.) then dripliq=drip*intwratio dripsn = drip - dripliq snwe=snwe+dripsn infwater=infwater+dripliq dripliq=0. dripsn = 0. else snwe=snwe+drip endif endif snhei=snwe*rhowater/rhosn NEWSN=NEWSN*rhowater/rhonewsn ENDIF IF(SNHEI.GT.0.0) THEN !-- SNOW on the ground !--- Land-use category should be changed to snow/ice for grid points with snow>0 ILAND=ISICE !24nov15 - based on field exp on Pleasant View soccer fields ! if(meltfactor > 1.5) then ! all veg. types, except forests ! SNHEI_CRIT=0.01601*1.e3/rhosn ! Petzold - 1 cm of fresh snow overwrites effects from old snow. ! Need to test SNHEI_CRIT_newsn=0.01 ! SNHEI_CRIT_newsn=0.01 ! else ! forests ! SNHEI_CRIT=0.02*1.e3/rhosn ! SNHEI_CRIT_newsn=0.001*1.e3/rhosn ! endif SNOWFRAC=MIN(1.,SNHEI/(2.*SNHEI_CRIT)) !24nov15 - SNOWFRAC for urban category < 0.75 if(ivgtyp == urban) snowfrac=min(0.75,snowfrac) ! if(meltfactor > 1.5) then ! if(isltyp > 9 .and. isltyp < 13) then !24nov15 clay soil types - SNOFRAC < 0.9 ! snowfrac=min(0.9,snowfrac) ! endif ! else !24nov15 - SNOWFRAC for forests < 0.75 ! snowfrac=min(0.85,snowfrac) ! endif ! SNOWFRAC=MIN(1.,SNHEI/(2.*SNHEI_CRIT)) ! elseif(snowfrac < 0.3 .and. tabs > 275.) then ! if(snowfrac < 0.3.and. tabs > 275.) snow_mosaic = 1. if(snowfrac < 0.75) snow_mosaic = 1. if(newsn > 0. ) SNOWFRACnewsn=MIN(1.,SNHEI/SNHEI_CRIT_newsn) KEEP_SNOW_ALBEDO = 0. IF (NEWSN > 0. .and. snowfracnewsn > 0.99) THEN ! new snow KEEP_SNOW_ALBEDO = 1. snow_mosaic=0. ! ??? ENDIF !7Mar18 - turn off snow mosaic for T<271K to prevent from too warm ! temperature and loss of low-level clouds in HRRR (case 2 Feb. 2018, 15z) ! IF (TABS < 271.) then ! snow_mosaic=0. ! ENDIF IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'SNHEI_CRIT,SNOWFRAC,SNHEI_CRIT_newsn,SNOWFRACnewsn', & SNHEI_CRIT,SNOWFRAC,SNHEI_CRIT_newsn,SNOWFRACnewsn ENDIF !-- Set znt for snow from VEGPARM table (snow/ice landuse), except for !-- land-use types with higher roughness (forests, urban). !5mar12 IF(znt.lt.0.2 .and. snowfrac.gt.0.99) znt=z0tbl(iland) ! IF(newsn==0. .and. znt.lt.0.2 .and. snowfrac.gt.0.99) znt=z0tbl(iland) IF(newsn.eq.0. .and. znt.le.0.2 .and. IVGTYP.ne.isice) then if( snhei .le. 2.*ZNT)then znt=0.55*znt+0.45*z0tbl(iland) elseif( snhei .gt. 2.*ZNT .and. snhei .le. 4.*ZNT)then znt=0.2*znt+0.8*z0tbl(iland) elseif(snhei > 4.*ZNT) then znt=z0tbl(iland) endif ENDIF !--- GSWNEW in-coming solar for snow on land or on ice ! GSWNEW=GSWnew/(1.-ALB) !-- Time to update snow and ice albedo IF(SEAICE .LT. 0.5) THEN !----- SNOW on soil !-- ALB dependence on snow depth ! ALB_SNOW across Canada's forested areas is very low - 0.27-0.35, this ! causes significant warm biases. Limiting ALB in these areas to be higher than 0.4 ! hwlps with these biases.. if( snow_mosaic == 1.) then ALBsn=alb_snow ! ALBsn=max(0.4,alb_snow) Emiss= emissn else ALBsn = MAX(keep_snow_albedo*alb_snow, & MIN((alb_snow_free + & (alb_snow - alb_snow_free) * snowfrac), alb_snow)) Emiss = MAX(keep_snow_albedo*emissn, & MIN((emiss_snowfree + & (emissn - emiss_snowfree) * snowfrac), emissn)) endif IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN ! if(i.eq.279.and.j.eq.263) then print *,'Snow on soil ALBsn,emiss,snow_mosaic',i,j,ALBsn,emiss,snow_mosaic ENDIF !28mar11 if canopy is covered with snow to 95% of its capacity and snow depth is ! higher than patchy snow treshold - then snow albedo is not less than 0.55 ! (inspired by the flight from Fairbanks to Seatle) !test if(cst.ge.0.95*sat .and. snowfrac .gt.0.99)then ! albsn=max(alb_snow,0.55) ! endif !-- ALB dependence on snow temperature. When snow temperature is !-- below critical value of -10C - no change to albedo. !-- If temperature is higher that -10C then albedo is decreasing. !-- The minimum albedo at t=0C for snow on land is 15% less than !-- albedo of temperatures below -10C. if(albsn.lt.0.4 .or. keep_snow_albedo==1) then ALB=ALBsn ! ALB=max(0.4,alb_snow) else !-- change albedo when no fresh snow and snow albedo is higher than 0.5 ALB = MIN(ALBSN,MAX(ALBSN - 0.1*(soilt - 263.15)/ & (273.15-263.15)*ALBSN, ALBSN - 0.05)) endif ELSE !----- SNOW on ice if( snow_mosaic == 1.) then ALBsn=alb_snow Emiss= emissn else ALBsn = MAX(keep_snow_albedo*alb_snow, & MIN((albice + (alb_snow - albice) * snowfrac), alb_snow)) Emiss = MAX(keep_snow_albedo*emissn, & MIN((emiss_snowfree + & (emissn - emiss_snowfree) * snowfrac), emissn)) endif IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'Snow on ice snow_mosaic,ALBsn,emiss',i,j,ALBsn,emiss,snow_mosaic ENDIF !-- ALB dependence on snow temperature. When snow temperature is !-- below critical value of -10C - no change to albedo. !-- If temperature is higher that -10C then albedo is decreasing. if(albsn.lt.alb_snow .or. keep_snow_albedo .eq.1.)then ALB=ALBsn else !-- change albedo when no fresh snow ALB = MIN(ALBSN,MAX(ALBSN - 0.15*ALBSN*(soilt - 263.15)/ & (273.15-263.15), ALBSN - 0.1)) endif ENDIF if (snow_mosaic==1.) then !may 2014 - treat separately snow-free and snow-covered areas if(SEAICE .LT. 0.5) then ! LAND ! portion not covered with snow ! compute absorbed GSW for snow-free portion gswnew=GSWin*(1.-alb_snow_free) !-------------- T3 = STBOLT*SOILT*SOILT*SOILT UPFLUX = T3 *SOILT XINET = EMISS_snowfree*(GLW-UPFLUX) RNET = GSWnew + XINET IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN ! if(i.eq.442.and.j.eq.260) then print *,'Fractional snow - snowfrac=',snowfrac print *,'Snowfrac<1 GSWin,GSWnew -',GSWin,GSWnew,'SOILT, RNET',soilt,rnet ENDIF do k=1,nzs soilm1ds(k) = soilm1d(k) ts1ds(k) = ts1d(k) smfrkeeps(k) = smfrkeep(k) keepfrs(k) = keepfr(k) soilices(k) = soilice(k) soiliqws(k) = soiliqw(k) enddo soilts = soilt qvgs = qvg qsgs = qsg qcgs = qcg csts = cst mavails = mavail smelt=0. runoff1s=0. runoff2s=0. ilands = ivgtyp CALL SOIL(spp_lsm,rstochcol,fieldcol_sf, & !--- input variables i,j,ilands,isoil,delt,ktau,conflx,nzs,nddzs,nroot, & PRCPMS,RAINF,PATM,QVATM,QCATM,GLW,GSWnew,gswin, & EMISS_snowfree,RNET,QKMS,TKMS,PC,csts,dripliq, & infwater,rho,vegfrac,lai,myj, & !--- soil fixed fields QWRTZ,rhocs,dqm,qmin,ref,wilt, & psis,bclh,ksat,sat,cn, & zsmain,zshalf,DTDZS,DTDZS2,tbq, & !--- constants lv,CP,rovcp,G0,cw,stbolt,tabs, & KQWRTZ,KICE,KWT, & !--- output variables for snow-free portion soilm1ds,ts1ds,smfrkeeps,keepfrs, & dews,soilts,qvgs,qsgs,qcgs,edir1s,ec1s, & ett1s,eetas,qfxs,hfxs,ss,evapls,prcpls,fltots,runoff1s, & runoff2s,mavails,soilices,soiliqws, & infiltrs,smf) else ! SEA ICE ! portion not covered with snow ! compute absorbed GSW for snow-free portion gswnew=GSWin*(1.-albice) !-------------- T3 = STBOLT*SOILT*SOILT*SOILT UPFLUX = T3 *SOILT XINET = EMISS_snowfree*(GLW-UPFLUX) RNET = GSWnew + XINET IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN ! if(i.eq.442.and.j.eq.260) then print *,'Fractional snow - snowfrac=',snowfrac print *,'Snowfrac<1 GSWin,GSWnew -',GSWin,GSWnew,'SOILT, RNET',soilt,rnet ENDIF do k=1,nzs ts1ds(k) = ts1d(k) enddo soilts = soilt qvgs = qvg qsgs = qsg qcgs = qcg smelt=0. runoff1s=0. runoff2s=0. CALL SICE( & !--- input variables i,j,iland,isoil,delt,ktau,conflx,nzs,nddzs,nroot, & PRCPMS,RAINF,PATM,QVATM,QCATM,GLW,GSWnew, & 0.98,RNET,QKMS,TKMS,rho,myj, & !--- sea ice parameters tice,rhosice,capice,thdifice, & zsmain,zshalf,DTDZS,DTDZS2,tbq, & !--- constants lv,CP,rovcp,cw,stbolt,tabs, & !--- output variable ts1ds,dews,soilts,qvgs,qsgs,qcgs, & eetas,qfxs,hfxs,ss,evapls,prcpls,fltots & ) edir1 = eeta*1.e-3 ec1 = 0. ett1 = 0. runoff1 = prcpms runoff2 = 0. mavail = 1. infiltr=0. cst=0. do k=1,nzs soilm1d(k)=1. soiliqw(k)=0. soilice(k)=1. smfrkeep(k)=1. keepfr(k)=0. enddo endif ! seaice < 0.5 !return gswnew to incoming solar IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN ! if(i.eq.442.and.j.eq.260) then print *,'gswnew,alb_snow_free,alb',gswnew,alb_snow_free,alb ENDIF ! gswnew=gswnew/(1.-alb_snow_free) IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN ! if(i.eq.442.and.j.eq.260) then print *,'Incoming GSWnew snowfrac<1 -',gswnew ENDIF endif ! snow_mosaic=1. !--- recompute absorbed solar radiation and net radiation !--- for updated value of snow albedo - ALB gswnew=GSWin*(1.-alb) ! print *,'SNOW fraction GSWnew',gswnew,'alb=',alb !-------------- T3 = STBOLT*SOILT*SOILT*SOILT UPFLUX = T3 *SOILT XINET = EMISS*(GLW-UPFLUX) RNET = GSWnew + XINET IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN ! if(i.eq.442.and.j.eq.260) then ! if(i.eq.271.and.j.eq.242) then print *,'RNET=',rnet print *,'SNOW - I,J,newsn,snwe,snhei,GSW,GSWnew,GLW,UPFLUX,ALB',& i,j,newsn,snwe,snhei,GSW,GSWnew,GLW,UPFLUX,ALB ENDIF if (SEAICE .LT. 0.5) then ! LAND if(snow_mosaic==1.)then snfr=1. else snfr=snowfrac endif CALL SNOWSOIL (spp_lsm,rstochcol,fieldcol_sf, & !--- input variables i,j,isoil,delt,ktau,conflx,nzs,nddzs,nroot, & meltfactor,rhonewsn,SNHEI_CRIT, & ! new ILAND,PRCPMS,RAINF,NEWSN,snhei,SNWE,snfr, & RHOSN,PATM,QVATM,QCATM, & GLW,GSWnew,GSWin,EMISS,RNET,IVGTYP, & QKMS,TKMS,PC,CST,dripsn,infwater, & RHO,VEGFRAC,ALB,ZNT,lai, & MYJ, & !--- soil fixed fields QWRTZ,rhocs,dqm,qmin,ref,wilt,psis,bclh,ksat, & sat,cn,zsmain,zshalf,DTDZS,DTDZS2,tbq, & !--- constants lv,CP,rovcp,G0,cw,stbolt,tabs, & KQWRTZ,KICE,KWT, & !--- output variables ilnb,snweprint,snheiprint,rsm, & soilm1d,ts1d,smfrkeep,keepfr, & dew,soilt,soilt1,tsnav,qvg,qsg,qcg, & SMELT,SNOH,SNFLX,SNOM,edir1,ec1,ett1,eeta, & qfx,hfx,s,sublim,prcpl,fltot,runoff1,runoff2, & mavail,soilice,soiliqw,infiltr ) else ! SEA ICE if(snow_mosaic==1.)then snfr=1. else snfr=snowfrac endif CALL SNOWSEAICE ( & i,j,isoil,delt,ktau,conflx,nzs,nddzs, & meltfactor,rhonewsn,SNHEI_CRIT, & ! new ILAND,PRCPMS,RAINF,NEWSN,snhei,SNWE,snfr, & RHOSN,PATM,QVATM,QCATM, & GLW,GSWnew,EMISS,RNET, & QKMS,TKMS,RHO,myj, & !--- sea ice parameters ALB,ZNT, & tice,rhosice,capice,thdifice, & zsmain,zshalf,DTDZS,DTDZS2,tbq, & !--- constants lv,CP,rovcp,cw,stbolt,tabs, & !--- output variables ilnb,snweprint,snheiprint,rsm,ts1d, & dew,soilt,soilt1,tsnav,qvg,qsg,qcg, & SMELT,SNOH,SNFLX,SNOM,eeta, & qfx,hfx,s,sublim,prcpl,fltot & ) edir1 = eeta*1.e-3 ec1 = 0. ett1 = 0. runoff1 = smelt runoff2 = 0. mavail = 1. infiltr=0. cst=0. do k=1,nzs soilm1d(k)=1. soiliqw(k)=0. soilice(k)=1. smfrkeep(k)=1. keepfr(k)=0. enddo endif if(snhei.eq.0.) then !--- all snow is melted alb=alb_snow_free iland=ivgtyp endif if (snow_mosaic==1.) then ! May 2014 - now combine snow covered and snow-free land fluxes, soil temp, moist, ! etc. if(SEAICE .LT. 0.5) then ! LAND IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN ! if(i.eq.442.and.j.eq.260) then print *,'SOILT snow on land', ktau, i,j,soilt print *,'SOILT on snow-free land', i,j,soilts print *,'ts1d,ts1ds',i,j,ts1d,ts1ds print *,' SNOW flux',i,j, snflx print *,' Ground flux on snow-covered land',i,j, s print *,' Ground flux on snow-free land', i,j,ss print *,' CSTS, CST', i,j,csts,cst ENDIF do k=1,nzs soilm1d(k) = soilm1ds(k)*(1.-snowfrac) + soilm1d(k)*snowfrac ts1d(k) = ts1ds(k)*(1.-snowfrac) + ts1d(k)*snowfrac smfrkeep(k) = smfrkeeps(k)*(1.-snowfrac) + smfrkeep(k)*snowfrac if(snowfrac > 0.5) then keepfr(k) = keepfr(k) else keepfr(k) = keepfrs(k) endif soilice(k) = soilices(k)*(1.-snowfrac) + soilice(k)*snowfrac soiliqw(k) = soiliqws(k)*(1.-snowfrac) + soiliqw(k)*snowfrac enddo dew = dews*(1.-snowfrac) + dew*snowfrac soilt = soilts*(1.-snowfrac) + soilt*snowfrac qvg = qvgs*(1.-snowfrac) + qvg*snowfrac qsg = qsgs*(1.-snowfrac) + qsg*snowfrac qcg = qcgs*(1.-snowfrac) + qcg*snowfrac edir1 = edir1s*(1.-snowfrac) + edir1*snowfrac ec1 = ec1s*(1.-snowfrac) + ec1*snowfrac cst = csts*(1.-snowfrac) + cst*snowfrac ett1 = ett1s*(1.-snowfrac) + ett1*snowfrac eeta = eetas*(1.-snowfrac) + eeta*snowfrac qfx = qfxs*(1.-snowfrac) + qfx*snowfrac hfx = hfxs*(1.-snowfrac) + hfx*snowfrac s = ss*(1.-snowfrac) + s*snowfrac evapl = evapls*(1.-snowfrac) sublim = sublim*snowfrac prcpl = prcpls*(1.-snowfrac) + prcpl*snowfrac fltot = fltots*(1.-snowfrac) + fltot*snowfrac !alb ALB = MAX(keep_snow_albedo*alb, & MIN((alb_snow_free + (alb - alb_snow_free) * snowfrac), alb)) Emiss = MAX(keep_snow_albedo*emissn, & MIN((emiss_snowfree + & (emissn - emiss_snowfree) * snowfrac), emissn)) ! alb=alb_snow_free*(1.-snowfrac) + alb*snowfrac ! emiss=emiss_snowfree*(1.-snowfrac) + emissn*snowfrac ! if(abs(fltot) > 2.) then ! print *,'i,j,fltot,snowfrac,fltots',fltot,snowfrac,fltots,i,j ! endif runoff1 = runoff1s*(1.-snowfrac) + runoff1*snowfrac runoff2 = runoff2s*(1.-snowfrac) + runoff2*snowfrac smelt = smelt * snowfrac snoh = snoh * snowfrac snflx = snflx * snowfrac snom = snom * snowfrac mavail = mavails*(1.-snowfrac) + 1.*snowfrac infiltr = infiltrs*(1.-snowfrac) + infiltr*snowfrac IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,' Ground flux combined', i,j, s print *,'SOILT combined on land', soilt print *,'TS combined on land', ts1d ENDIF else ! SEA ICE ! Now combine fluxes for snow-free sea ice and snow-covered area IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'SOILT snow on ice', soilt ENDIF do k=1,nzs ts1d(k) = ts1ds(k)*(1.-snowfrac) + ts1d(k)*snowfrac enddo dew = dews*(1.-snowfrac) + dew*snowfrac soilt = soilts*(1.-snowfrac) + soilt*snowfrac qvg = qvgs*(1.-snowfrac) + qvg*snowfrac qsg = qsgs*(1.-snowfrac) + qsg*snowfrac qcg = qcgs*(1.-snowfrac) + qcg*snowfrac eeta = eetas*(1.-snowfrac) + eeta*snowfrac qfx = qfxs*(1.-snowfrac) + qfx*snowfrac hfx = hfxs*(1.-snowfrac) + hfx*snowfrac s = ss*(1.-snowfrac) + s*snowfrac sublim = eeta prcpl = prcpls*(1.-snowfrac) + prcpl*snowfrac fltot = fltots*(1.-snowfrac) + fltot*snowfrac !alb ALB = MAX(keep_snow_albedo*alb, & MIN((albice + (alb - alb_snow_free) * snowfrac), alb)) Emiss = MAX(keep_snow_albedo*emissn, & MIN((emiss_snowfree + & (emissn - emiss_snowfree) * snowfrac), emissn)) ! alb=alb_snow_free*(1.-snowfrac) + alb*snowfrac ! emiss=1.*(1.-snowfrac) + emissn*snowfrac runoff1 = runoff1s*(1.-snowfrac) + runoff1*snowfrac runoff2 = runoff2s*(1.-snowfrac) + runoff2*snowfrac smelt = smelt * snowfrac snoh = snoh * snowfrac snflx = snflx * snowfrac snom = snom * snowfrac IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'SOILT combined on ice', soilt ENDIF endif endif ! snow_mosaic = 1. ! run-total accumulated snow based on snowfall and snowmelt in [m] snowfallac = snowfallac + max(0.,(newsn - rhowater/rhonewsn*smelt*delt*newsnowratio)) ELSE !--- no snow snheiprint=0. snweprint=0. smelt=0. !-------------- T3 = STBOLT*SOILT*SOILT*SOILT UPFLUX = T3 *SOILT XINET = EMISS*(GLW-UPFLUX) RNET = GSWnew + XINET IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'NO snow on the ground GSWnew -',GSWnew,'RNET=',rnet ENDIF if(SEAICE .LT. 0.5) then ! LAND CALL SOIL(spp_lsm,rstochcol,fieldcol_sf, & !--- input variables i,j,iland,isoil,delt,ktau,conflx,nzs,nddzs,nroot, & PRCPMS,RAINF,PATM,QVATM,QCATM,GLW,GSWnew,GSWin, & EMISS,RNET,QKMS,TKMS,PC,cst,drip,infwater, & rho,vegfrac,lai,myj, & !--- soil fixed fields QWRTZ,rhocs,dqm,qmin,ref,wilt, & psis,bclh,ksat,sat,cn, & zsmain,zshalf,DTDZS,DTDZS2,tbq, & !--- constants lv,CP,rovcp,G0,cw,stbolt,tabs, & KQWRTZ,KICE,KWT, & !--- output variables soilm1d,ts1d,smfrkeep,keepfr, & dew,soilt,qvg,qsg,qcg,edir1,ec1, & ett1,eeta,qfx,hfx,s,evapl,prcpl,fltot,runoff1, & runoff2,mavail,soilice,soiliqw, & infiltr,smf) else ! SEA ICE ! If current ice albedo is not the same as from the previous time step, then ! update GSW, ALB and RNET for surface energy budget if(ALB.ne.ALBice) GSWnew=GSW/(1.-ALB)*(1.-ALBice) alb=albice RNET = GSWnew + XINET CALL SICE( & !--- input variables i,j,iland,isoil,delt,ktau,conflx,nzs,nddzs,nroot, & PRCPMS,RAINF,PATM,QVATM,QCATM,GLW,GSWnew, & EMISS,RNET,QKMS,TKMS,rho,myj, & !--- sea ice parameters tice,rhosice,capice,thdifice, & zsmain,zshalf,DTDZS,DTDZS2,tbq, & !--- constants lv,CP,rovcp,cw,stbolt,tabs, & !--- output variables ts1d,dew,soilt,qvg,qsg,qcg, & eeta,qfx,hfx,s,evapl,prcpl,fltot & ) edir1 = eeta*1.e-3 ec1 = 0. ett1 = 0. runoff1 = prcpms runoff2 = 0. mavail = 1. infiltr=0. cst=0. do k=1,nzs soilm1d(k)=1. soiliqw(k)=0. soilice(k)=1. smfrkeep(k)=1. keepfr(k)=0. enddo endif ENDIF ! RETURN ! END !--------------------------------------------------------------- END SUBROUTINE SFCTMP !--------------------------------------------------------------- FUNCTION QSN(TN,T) !**************************************************************** REAL, DIMENSION(1:5001), INTENT(IN ) :: T REAL, INTENT(IN ) :: TN REAL QSN, R,R1,R2 INTEGER I R=(TN-173.15)/.05+1. I=INT(R) IF(I.GE.1) goto 10 I=1 R=1. 10 IF(I.LE.5000) GOTO 20 I=5000 R=5001. 20 R1=T(I) R2=R-I QSN=(T(I+1)-R1)*R2 + R1 ! print *,' in QSN, I,R,R1,R2,T(I+1),TN, QSN', I,R,r1,r2,t(i+1),tn,QSN ! RETURN ! END !----------------------------------------------------------------------- END FUNCTION QSN !------------------------------------------------------------------------ SUBROUTINE SOIL (spp_lsm,rstochcol, fieldcol_sf, & !--- input variables i,j,iland,isoil,delt,ktau,conflx,nzs,nddzs,nroot,& PRCPMS,RAINF,PATM,QVATM,QCATM, & GLW,GSW,GSWin,EMISS,RNET, & QKMS,TKMS,PC,cst,drip,infwater,rho,vegfrac,lai, & myj, & !--- soil fixed fields QWRTZ,rhocs,dqm,qmin,ref,wilt,psis,bclh,ksat, & sat,cn,zsmain,zshalf,DTDZS,DTDZS2,tbq, & !--- constants xlv,CP,rovcp,G0_P,cw,stbolt,TABS, & KQWRTZ,KICE,KWT, & !--- output variables soilmois,tso,smfrkeep,keepfr, & dew,soilt,qvg,qsg,qcg, & edir1,ec1,ett1,eeta,qfx,hfx,s,evapl, & prcpl,fltot,runoff1,runoff2,mavail,soilice, & soiliqw,infiltrp,smf) !************************************************************* ! Energy and moisture budget for vegetated surfaces ! without snow, heat diffusion and Richards eqns. in ! soil ! ! DELT - time step (s) ! ktau - numver of time step ! CONFLX - depth of constant flux layer (m) ! J,I - the location of grid point ! IME, JME, KME, NZS - dimensions of the domain ! NROOT - number of levels within the root zone ! PRCPMS - precipitation rate in m/s ! PATM - pressure [bar] ! QVATM,QCATM - cloud and water vapor mixing ratio (kg/kg) ! at the first atm. level ! GLW, GSW - incoming longwave and absorbed shortwave ! radiation at the surface (W/m^2) ! EMISS,RNET - emissivity of the ground surface (0-1) and net ! radiation at the surface (W/m^2) ! QKMS - exchange coefficient for water vapor in the ! surface layer (m/s) ! TKMS - exchange coefficient for heat in the surface ! layer (m/s) ! PC - plant coefficient (resistance) (0-1) ! RHO - density of atmosphere near sueface (kg/m^3) ! VEGFRAC - greeness fraction ! RHOCS - volumetric heat capacity of dry soil ! DQM, QMIN - porosity minus residual soil moisture QMIN (m^3/m^3) ! REF, WILT - field capacity soil moisture and the ! wilting point (m^3/m^3) ! PSIS - matrix potential at saturation (m) ! BCLH - exponent for Clapp-Hornberger parameterization ! KSAT - saturated hydraulic conductivity (m/s) ! SAT - maximum value of water intercepted by canopy (m) ! CN - exponent for calculation of canopy water ! ZSMAIN - main levels in soil (m) ! ZSHALF - middle of the soil layers (m) ! DTDZS,DTDZS2 - dt/(2.*dzshalf*dzmain) and dt/dzshalf in soil ! TBQ - table to define saturated mixing ration ! of water vapor for given temperature and pressure ! SOILMOIS,TSO - soil moisture (m^3/m^3) and temperature (K) ! DEW - dew in kg/m^2s ! SOILT - skin temperature (K) ! QSG,QVG,QCG - saturated mixing ratio, mixing ratio of ! water vapor and cloud at the ground ! surface, respectively (kg/kg) ! EDIR1, EC1, ETT1, EETA - direct evaporation, evaporation of ! canopy water, transpiration in kg m-2 s-1 and total ! evaporation in m s-1. ! QFX, HFX - latent and sensible heat fluxes (W/m^2) ! S - soil heat flux in the top layer (W/m^2) ! RUNOFF - surface runoff (m/s) ! RUNOFF2 - underground runoff (m) ! MAVAIL - moisture availability in the top soil layer (0-1) ! INFILTRP - infiltration flux from the top of soil domain (m/s) ! !***************************************************************** IMPLICIT NONE !----------------------------------------------------------------- !--- input variables INTEGER, INTENT(IN ) :: nroot,ktau,nzs , & nddzs !nddzs=2*(nzs-2) INTEGER, INTENT(IN ) :: i,j,iland,isoil REAL, INTENT(IN ) :: DELT,CONFLX LOGICAL, INTENT(IN ) :: myj !--- 3-D Atmospheric variables REAL, & INTENT(IN ) :: PATM, & QVATM, & QCATM !--- 2-D variables REAL, & INTENT(IN ) :: GLW, & GSW, & GSWin, & EMISS, & RHO, & PC, & VEGFRAC, & lai, & infwater, & QKMS, & TKMS !--- soil properties REAL, & INTENT(IN ) :: RHOCS, & BCLH, & DQM, & KSAT, & PSIS, & QMIN, & QWRTZ, & REF, & WILT REAL, INTENT(IN ) :: CN, & CW, & KQWRTZ, & KICE, & KWT, & XLV, & g0_p REAL, DIMENSION(1:NZS), INTENT(IN) :: ZSMAIN, & ZSHALF, & DTDZS2 REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS REAL, DIMENSION(1:5001), INTENT(IN) :: TBQ !--- input/output variables !-------- 3-d soil moisture and temperature REAL, DIMENSION( 1:nzs ) , & INTENT(INOUT) :: TSO, & SOILMOIS, & SMFRKEEP REAL, DIMENSION(1:NZS), INTENT(IN) :: rstochcol REAL, DIMENSION(1:NZS), INTENT(INOUT) :: fieldcol_sf REAL, DIMENSION( 1:nzs ) , & INTENT(INOUT) :: KEEPFR !-------- 2-d variables REAL, & INTENT(INOUT) :: DEW, & CST, & DRIP, & EDIR1, & EC1, & ETT1, & EETA, & EVAPL, & PRCPL, & MAVAIL, & QVG, & QSG, & QCG, & RNET, & QFX, & HFX, & S, & SAT, & RUNOFF1, & RUNOFF2, & SOILT !-------- 1-d variables INTEGER , INTENT(IN) :: spp_lsm REAL, DIMENSION(1:NZS), INTENT(OUT) :: SOILICE, & SOILIQW !--- Local variables REAL :: INFILTRP, transum , & RAINF, PRCPMS , & TABS, T3, UPFLUX, XINET REAL :: CP,rovcp,G0,LV,STBOLT,xlmelt,dzstop , & can,epot,fac,fltot,ft,fq,hft , & q1,ras,rhoice,sph , & trans,zn,ci,cvw,tln,tavln,pi , & DD1,CMC2MS,DRYCAN,WETCAN , & INFMAX,RIW, X REAL, DIMENSION(1:NZS) :: transp,cap,diffu,hydro , & thdif,tranf,tav,soilmoism , & soilicem,soiliqwm,detal , & fwsat,lwsat,told,smold REAL :: soiltold,smf REAL :: soilres, alfa, fex, fex_fc, fc, psit INTEGER :: nzs1,nzs2,k !----------------------------------------------------------------- !-- define constants ! STBOLT=5.670151E-8 RHOICE=900. CI=RHOICE*2100. XLMELT=3.35E+5 cvw=cw ! SAT=0.0004 prcpl=prcpms smf=0. soiltold = soilt wetcan=0. drycan=1. !--- Initializing local arrays DO K=1,NZS TRANSP (K)=0. soilmoism(k)=0. soilice (k)=0. soiliqw (k)=0. soilicem (k)=0. soiliqwm (k)=0. lwsat (k)=0. fwsat (k)=0. tav (k)=0. cap (k)=0. thdif (k)=0. diffu (k)=0. hydro (k)=0. tranf (k)=0. detal (k)=0. told (k)=0. smold (k)=0. ENDDO NZS1=NZS-1 NZS2=NZS-2 dzstop=1./(zsmain(2)-zsmain(1)) RAS=RHO*1.E-3 RIW=rhoice*1.e-3 !--- Computation of volumetric content of ice in soil DO K=1,NZS !- main levels tln=log(tso(k)/273.15) if(tln.lt.0.) then soiliqw(k)=(dqm+qmin)*(XLMELT* & (tso(k)-273.15)/tso(k)/9.81/psis) & **(-1./bclh)-qmin soiliqw(k)=max(0.,soiliqw(k)) soiliqw(k)=min(soiliqw(k),soilmois(k)) soilice(k)=(soilmois(k)-soiliqw(k))/RIW !---- melting and freezing is balanced, soil ice cannot increase if(keepfr(k).eq.1.) then soilice(k)=min(soilice(k),smfrkeep(k)) soiliqw(k)=max(0.,soilmois(k)-soilice(k)*riw) endif else soilice(k)=0. soiliqw(k)=soilmois(k) endif ENDDO DO K=1,NZS1 !- middle of soil layers tav(k)=0.5*(tso(k)+tso(k+1)) soilmoism(k)=0.5*(soilmois(k)+soilmois(k+1)) tavln=log(tav(k)/273.15) if(tavln.lt.0.) then soiliqwm(k)=(dqm+qmin)*(XLMELT* & (tav(k)-273.15)/tav(k)/9.81/psis) & **(-1./bclh)-qmin fwsat(k)=dqm-soiliqwm(k) lwsat(k)=soiliqwm(k)+qmin soiliqwm(k)=max(0.,soiliqwm(k)) soiliqwm(k)=min(soiliqwm(k), soilmoism(k)) soilicem(k)=(soilmoism(k)-soiliqwm(k))/riw !---- melting and freezing is balanced, soil ice cannot increase if(keepfr(k).eq.1.) then soilicem(k)=min(soilicem(k), & 0.5*(smfrkeep(k)+smfrkeep(k+1))) soiliqwm(k)=max(0.,soilmoism(k)-soilicem(k)*riw) fwsat(k)=dqm-soiliqwm(k) lwsat(k)=soiliqwm(k)+qmin endif else soilicem(k)=0. soiliqwm(k)=soilmoism(k) lwsat(k)=dqm+qmin fwsat(k)=0. endif ENDDO do k=1,nzs if(soilice(k).gt.0.) then smfrkeep(k)=soilice(k) else smfrkeep(k)=soilmois(k)/riw endif enddo !****************************************************************** ! SOILPROP computes thermal diffusivity, and diffusional and ! hydraulic condeuctivities !****************************************************************** CALL SOILPROP(spp_lsm,rstochcol,fieldcol_sf, & !--- input variables nzs,fwsat,lwsat,tav,keepfr, & soilmois,soiliqw,soilice, & soilmoism,soiliqwm,soilicem, & !--- soil fixed fields QWRTZ,rhocs,dqm,qmin,psis,bclh,ksat, & !--- constants riw,xlmelt,CP,G0_P,cvw,ci, & kqwrtz,kice,kwt, & !--- output variables thdif,diffu,hydro,cap) !******************************************************************** !--- CALCULATION OF CANOPY WATER (Smirnova et al., 1996, EQ.16) AND DEW ! DRIP=0. ! DD1=0. FQ=QKMS Q1=-QKMS*RAS*(QVATM - QSG) DEW=0. IF(QVATM.GE.QSG)THEN DEW=FQ*(QVATM-QSG) ENDIF ! IF(DEW.NE.0.)THEN ! DD1=CST+DELT*(PRCPMS +DEW*RAS) ! ELSE ! DD1=CST+ & ! DELT*(PRCPMS+RAS*FQ*(QVATM-QSG) & ! *(CST/SAT)**CN) ! ENDIF ! DD1=CST+DELT*PRCPMS ! IF(DD1.LT.0.) DD1=0. ! if(vegfrac.eq.0.)then ! cst=0. ! drip=0. ! endif ! IF (vegfrac.GT.0.) THEN ! CST=DD1 ! IF(CST.GT.SAT) THEN ! CST=SAT ! DRIP=DD1-SAT ! ENDIF ! ENDIF ! !--- WETCAN is the fraction of vegetated area covered by canopy !--- water, and DRYCAN is the fraction of vegetated area where !--- transpiration may take place. WETCAN=min(0.25,(CST/SAT)**CN) ! if(lai > 1.) wetcan=wetcan/lai DRYCAN=1.-WETCAN !************************************************************** ! TRANSF computes transpiration function !************************************************************** CALL TRANSF(i,j, & !--- input variables nzs,nroot,soiliqw,tabs,lai,gswin, & !--- soil fixed fields dqm,qmin,ref,wilt,zshalf,pc,iland, & !--- output variables tranf,transum) !--- Save soil temp and moisture from the beginning of time step do k=1,nzs told(k)=tso(k) smold(k)=soilmois(k) enddo ! Sakaguchi and Zeng (2009) - dry soil resistance to evaporation ! 20jun18 - alfa is the ratio of surface Q to the saturation RH ! qvg = alfa * qsg ! if (vgtype==11) then ! MODIS wetland ! alfa=1. ! else fex=min(1.,soilmois(1)/dqm) fex=max(fex,0.01) psit=psis*fex ** (-bclh) psit = max(-1.e5, psit) alfa=min(1.,exp(g*psit/r_v/SOILT)) ! endif ! 20jun18 - at this time alfa is not used, and qvg is computed from the energy budget. alfa=1. ! field capacity ! 20jun18 - beta in Eq. (4) is called soilres here - it limits soil evaporation ! when soil moisture is below field capacity. [Lee and Pielke, 1992] ! This formulation agrees with obsevations when top layer is < 2 cm thick. ! Soilres = 1 for snow, glaciers and wetland. ! fc=ref - suggested in the paper ! fc=max(qmin,ref*0.5) ! used prior to 20jun18 change ! Switch from ref*0.5 to ref*0.25 will reduce soil resistance, increase direct ! evaporation, effects sparsely vegetated areas--> cooler during the day ! fc=max(qmin,ref*0.25) ! ! For now we'll go back to ref*0.5 fc=max(qmin,ref*0.5) fex_fc=1. if((soilmois(1)+qmin) > fc .or. (qvatm-qvg) > 0.) then ! 20jun18 - soil moisture is higher than field capacity or condensation process. soilres = 1. else fex_fc=min(1.,(soilmois(1)+qmin)/fc) fex_fc=max(fex_fc,0.01) soilres=0.25*(1.-cos(piconst*fex_fc))**2. endif IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN ! if (i==421.and.j==280) then print *,'fex,psit,psis,bclh,g,r_v,soilt,alfa,mavail,soilmois(1),fc,ref,soilres,fex_fc', & fex,psit,psis,bclh,g,r_v,soilt,alfa,mavail,soilmois(1),fc,ref,soilres,fex_fc endif !************************************************************** ! SOILTEMP soilves heat budget and diffusion eqn. in soil !************************************************************** CALL SOILTEMP( & !--- input variables i,j,iland,isoil, & delt,ktau,conflx,nzs,nddzs,nroot, & PRCPMS,RAINF, & PATM,TABS,QVATM,QCATM,EMISS,RNET, & QKMS,TKMS,PC,rho,vegfrac, lai, & thdif,cap,drycan,wetcan, & transum,dew,mavail,soilres,alfa, & !--- soil fixed fields dqm,qmin,bclh,zsmain,zshalf,DTDZS,tbq, & !--- constants xlv,CP,G0_P,cvw,stbolt, & !--- output variables tso,soilt,qvg,qsg,qcg,x) !************************************************************************ !--- CALCULATION OF DEW USING NEW VALUE OF QSG OR TRANSP IF NO DEW ETT1=0. DEW=0. IF(QVATM.GE.QSG)THEN DEW=QKMS*(QVATM-QSG) ETT1=0. DO K=1,NZS TRANSP(K)=0. ENDDO ELSE DO K=1,NROOT TRANSP(K)=VEGFRAC*RAS*QKMS* & (QVATM-QSG)* & TRANF(K)*DRYCAN/ZSHALF(NROOT+1) IF(TRANSP(K).GT.0.) TRANSP(K)=0. ETT1=ETT1-TRANSP(K) ENDDO DO k=nroot+1,nzs transp(k)=0. enddo ENDIF !-- Recalculate volumetric content of frozen water in soil DO K=1,NZS !- main levels tln=log(tso(k)/273.15) if(tln.lt.0.) then soiliqw(k)=(dqm+qmin)*(XLMELT* & (tso(k)-273.15)/tso(k)/9.81/psis) & **(-1./bclh)-qmin soiliqw(k)=max(0.,soiliqw(k)) soiliqw(k)=min(soiliqw(k),soilmois(k)) soilice(k)=(soilmois(k)-soiliqw(k))/riw !---- melting and freezing is balanced, soil ice cannot increase if(keepfr(k).eq.1.) then soilice(k)=min(soilice(k),smfrkeep(k)) soiliqw(k)=max(0.,soilmois(k)-soilice(k)*riw) endif else soilice(k)=0. soiliqw(k)=soilmois(k) endif ENDDO !************************************************************************* ! SOILMOIST solves moisture budget (Smirnova et al., 1996, EQ.22,28) ! and Richards eqn. !************************************************************************* CALL SOILMOIST ( & !-- input delt,nzs,nddzs,DTDZS,DTDZS2,RIW, & zsmain,zshalf,diffu,hydro, & QSG,QVG,QCG,QCATM,QVATM,-infwater, & QKMS,TRANSP,DRIP,DEW,0.,SOILICE,VEGFRAC, & 0.,soilres, & !-- soil properties DQM,QMIN,REF,KSAT,RAS,INFMAX, & !-- output SOILMOIS,SOILIQW,MAVAIL,RUNOFF1, & RUNOFF2,INFILTRP) !--- KEEPFR is 1 when the temperature and moisture in soil !--- are both increasing. In this case soil ice should not !--- be increasing according to the freezing curve. !--- Some part of ice is melted, but additional water is !--- getting frozen. Thus, only structure of frozen soil is !--- changed, and phase changes are not affecting the heat !--- transfer. This situation may happen when it rains on the !--- frozen soil. do k=1,nzs if (soilice(k).gt.0.) then if(tso(k).gt.told(k).and.soilmois(k).gt.smold(k)) then keepfr(k)=1. else keepfr(k)=0. endif endif enddo !--- THE DIAGNOSTICS OF SURFACE FLUXES T3 = STBOLT*SOILTold*SOILTold*SOILTold UPFLUX = T3 * 0.5*(SOILTold+SOILT) XINET = EMISS*(GLW-UPFLUX) ! RNET = GSW + XINET HFT=-TKMS*CP*RHO*(TABS-SOILT) HFX=-TKMS*CP*RHO*(TABS-SOILT) & *(P1000mb*0.00001/Patm)**ROVCP Q1=-QKMS*RAS*(QVATM - QSG) CMC2MS = 0. IF (Q1.LE.0.) THEN ! --- condensation EC1=0. EDIR1=0. ETT1=0. if(myj) then !-- moisture flux for coupling with MYJ PBL EETA=-QKMS*RAS*(QVATM/(1.+QVATM) - QSG/(1.+QSG))*1.E3 CST= CST-EETA*DELT*vegfrac IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN !!! IF(i.eq.374.and.j.eq.310.or. EETA.gt.0.0004) then print *,'Cond MYJ EETA',eeta,eeta*xlv, i,j ENDIF else ! myj !-- actual moisture flux from RUC LSM EETA= - RHO*DEW CST=CST+DELT*DEW*RAS * vegfrac IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN ! IF(i.eq.374.and.j.eq.310.or. EETA.gt.0.0004) then ! IF(i.eq.440.and.j.eq.180.or. QFX.gt.1000..or.i.eq.417.and.j.eq.540) then print *,'Cond RUC LSM EETA',EETA,eeta*xlv, i,j ENDIF endif ! myj QFX= XLV*EETA EETA= - RHO*DEW ELSE ! --- evaporation EDIR1 =-soilres*(1.-vegfrac)*QKMS*RAS* & (QVATM-QVG) CMC2MS=CST/DELT*RAS EC1 = Q1 * WETCAN * vegfrac IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN ! IF(i.eq.440.and.j.eq.180.or. QFX.gt.1000..or.i.eq.417.and.j.eq.540) then print *,'CST before update=',cst print *,'EC1=',EC1,'CMC2MS=',CMC2MS ENDIF ! ENDIF CST=max(0.,CST-EC1 * DELT) ! if (EC1 > CMC2MS) then ! EC1 = min(cmc2ms,ec1) ! CST = 0. ! endif if (myj) then !-- moisture flux for coupling with MYJ PBL EETA=-soilres*QKMS*RAS*(QVATM/(1.+QVATM) - QVG/(1.+QVG))*1.E3 else ! myj IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN ! IF(i.eq.440.and.j.eq.180.or. QFX.gt.1000..or.i.eq.417.and.j.eq.540) then print *,'QKMS,RAS,QVATM/(1.+QVATM),QVG/(1.+QVG),QSG ', & QKMS,RAS,QVATM/(1.+QVATM),QVG/(1.+QVG),QSG print *,'Q1*(1.-vegfrac),EDIR1',Q1*(1.-vegfrac),EDIR1 print *,'CST,WETCAN,DRYCAN',CST,WETCAN,DRYCAN print *,'EC1=',EC1,'ETT1=',ETT1,'CMC2MS=',CMC2MS,'CMC2MS*ras=',CMC2MS*ras ! print *,'MYJ EETA',eeta,eeta*xlv ENDIF !-- actual moisture flux from RUC LSM EETA = (EDIR1 + EC1 + ETT1)*1.E3 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN ! IF(i.eq.374.and.j.eq.310.or. EETA.gt.0.0004) then ! IF(i.eq.440.and.j.eq.180 .or. qfx.gt.1000..or.i.eq.417.and.j.eq.540) then print *,'RUC LSM EETA',EETA,eeta*xlv ENDIF endif ! myj QFX= XLV * EETA EETA = (EDIR1 + EC1 + ETT1)*1.E3 ENDIF IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'potential temp HFT ',HFT print *,'abs temp HFX ',HFX ENDIF EVAPL=EETA S=THDIF(1)*CAP(1)*DZSTOP*(TSO(1)-TSO(2)) ! Energy budget FLTOT=RNET-HFT-XLV*EETA-S-X IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN ! IF(i.eq.440.and.j.eq.180 .or. qfx.gt.1000..or.i.eq.417.and.j.eq.540) then print *,'SOIL - FLTOT,RNET,HFT,QFX,S,X=',i,j,FLTOT,RNET,HFT,XLV*EETA,s,x print *,'edir1,ec1,ett1,mavail,qkms,qvatm,qvg,qsg,vegfrac',& edir1,ec1,ett1,mavail,qkms,qvatm,qvg,qsg,vegfrac ENDIF if(detal(1) .ne. 0.) then ! SMF - energy of phase change in the first soil layer ! smf=xlmelt*1.e3*(soiliqwm(1)-soiliqwmold(1))/delt smf=fltot IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'detal(1),xlmelt,soiliqwm(1),delt',detal(1),xlmelt,soiliqwm(1),delt print *,'Implicit phase change in the first layer - smf=',smf ENDIF endif 222 CONTINUE 1123 FORMAT(I5,8F12.3) 1133 FORMAT(I7,8E12.4) 123 format(i6,f6.2,7f8.1) 122 FORMAT(1X,2I3,6F8.1,F8.3,F8.2) !------------------------------------------------------------------- END SUBROUTINE SOIL !------------------------------------------------------------------- SUBROUTINE SICE ( & !--- input variables i,j,iland,isoil,delt,ktau,conflx,nzs,nddzs,nroot, & PRCPMS,RAINF,PATM,QVATM,QCATM,GLW,GSW, & EMISS,RNET,QKMS,TKMS,rho,myj, & !--- sea ice parameters tice,rhosice,capice,thdifice, & zsmain,zshalf,DTDZS,DTDZS2,tbq, & !--- constants xlv,CP,rovcp,cw,stbolt,tabs, & !--- output variables tso,dew,soilt,qvg,qsg,qcg, & eeta,qfx,hfx,s,evapl,prcpl,fltot & ) !***************************************************************** ! Energy budget and heat diffusion eqns. for ! sea ice !************************************************************* IMPLICIT NONE !----------------------------------------------------------------- !--- input variables INTEGER, INTENT(IN ) :: nroot,ktau,nzs , & nddzs !nddzs=2*(nzs-2) INTEGER, INTENT(IN ) :: i,j,iland,isoil REAL, INTENT(IN ) :: DELT,CONFLX LOGICAL, INTENT(IN ) :: myj !--- 3-D Atmospheric variables REAL, & INTENT(IN ) :: PATM, & QVATM, & QCATM !--- 2-D variables REAL, & INTENT(IN ) :: GLW, & GSW, & EMISS, & RHO, & QKMS, & TKMS !--- sea ice properties REAL, DIMENSION(1:NZS) , & INTENT(IN ) :: & tice, & rhosice, & capice, & thdifice REAL, INTENT(IN ) :: & CW, & XLV REAL, DIMENSION(1:NZS), INTENT(IN) :: ZSMAIN, & ZSHALF, & DTDZS2 REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS REAL, DIMENSION(1:5001), INTENT(IN) :: TBQ !--- input/output variables !----soil temperature REAL, DIMENSION( 1:nzs ), INTENT(INOUT) :: TSO !-------- 2-d variables REAL, & INTENT(INOUT) :: DEW, & EETA, & EVAPL, & PRCPL, & QVG, & QSG, & QCG, & RNET, & QFX, & HFX, & S, & SOILT !--- Local variables REAL :: x,x1,x2,x4,tn,denom REAL :: RAINF, PRCPMS , & TABS, T3, UPFLUX, XINET REAL :: CP,rovcp,G0,LV,STBOLT,xlmelt,dzstop , & epot,fltot,ft,fq,hft,ras,cvw REAL :: FKT,D1,D2,D9,D10,DID,R211,R21,R22,R6,R7,D11 , & PI,H,FKQ,R210,AA,BB,PP,Q1,QS1,TS1,TQ2,TX2 , & TDENOM,QGOLD,SNOH REAL :: AA1,RHCS, icemelt REAL, DIMENSION(1:NZS) :: cotso,rhtso INTEGER :: nzs1,nzs2,k,k1,kn,kk !----------------------------------------------------------------- !-- define constants ! STBOLT=5.670151E-8 XLMELT=3.35E+5 cvw=cw prcpl=prcpms NZS1=NZS-1 NZS2=NZS-2 dzstop=1./(zsmain(2)-zsmain(1)) RAS=RHO*1.E-3 do k=1,nzs cotso(k)=0. rhtso(k)=0. enddo cotso(1)=0. rhtso(1)=TSO(NZS) DO 33 K=1,NZS2 KN=NZS-K K1=2*KN-3 X1=DTDZS(K1)*THDIFICE(KN-1) X2=DTDZS(K1+1)*THDIFICE(KN) FT=TSO(KN)+X1*(TSO(KN-1)-TSO(KN)) & -X2*(TSO(KN)-TSO(KN+1)) DENOM=1.+X1+X2-X2*cotso(K) cotso(K+1)=X1/DENOM rhtso(K+1)=(FT+X2*rhtso(K))/DENOM 33 CONTINUE !************************************************************************ !--- THE HEAT BALANCE EQUATION (Smirnova et al., 1996, EQ. 21,26) RHCS=CAPICE(1) H=1. FKT=TKMS D1=cotso(NZS1) D2=rhtso(NZS1) TN=SOILT D9=THDIFICE(1)*RHCS*dzstop D10=TKMS*CP*RHO R211=.5*CONFLX/DELT R21=R211*CP*RHO R22=.5/(THDIFICE(1)*DELT*dzstop**2) R6=EMISS *STBOLT*.5*TN**4 R7=R6/TN D11=RNET+R6 TDENOM=D9*(1.-D1+R22)+D10+R21+R7 & +RAINF*CVW*PRCPMS FKQ=QKMS*RHO R210=R211*RHO AA=XLS*(FKQ+R210)/TDENOM BB=(D10*TABS+R21*TN+XLS*(QVATM*FKQ & +R210*QVG)+D11+D9*(D2+R22*TN) & +RAINF*CVW*PRCPMS*max(273.15,TABS) & )/TDENOM AA1=AA PP=PATM*1.E3 AA1=AA1/PP IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN PRINT *,' VILKA-SEAICE1' print *,'D10,TABS,R21,TN,QVATM,FKQ', & D10,TABS,R21,TN,QVATM,FKQ print *,'RNET, EMISS, STBOLT, SOILT',RNET, EMISS, STBOLT, SOILT print *,'R210,QVG,D11,D9,D2,R22,RAINF,CVW,PRCPMS,TDENOM', & R210,QVG,D11,D9,D2,R22,RAINF,CVW,PRCPMS,TDENOM print *,'tn,aa1,bb,pp,fkq,r210', & tn,aa1,bb,pp,fkq,r210 ENDIF QGOLD=QSG CALL VILKA(TN,AA1,BB,PP,QS1,TS1,TBQ,KTAU,i,j,iland,isoil) !--- it is saturation over sea ice QVG=QS1 QSG=QS1 TSO(1)=min(271.4,TS1) QCG=0. !--- sea ice melting is not included in this simple approach !--- SOILT - skin temperature SOILT=TSO(1) !---- Final solution for soil temperature - TSO DO K=2,NZS KK=NZS-K+1 TSO(K)=min(271.4,rhtso(KK)+cotso(KK)*TSO(K-1)) END DO !--- CALCULATION OF DEW USING NEW VALUE OF QSG OR TRANSP IF NO DEW DEW=0. !--- THE DIAGNOSTICS OF SURFACE FLUXES T3 = STBOLT*TN*TN*TN UPFLUX = T3 *0.5*(TN+SOILT) XINET = EMISS*(GLW-UPFLUX) ! RNET = GSW + XINET HFT=-TKMS*CP*RHO*(TABS-SOILT) HFX=-TKMS*CP*RHO*(TABS-SOILT) & *(P1000mb*0.00001/Patm)**ROVCP Q1=-QKMS*RAS*(QVATM - QSG) IF (Q1.LE.0.) THEN ! --- condensation if(myj) then !-- moisture flux for coupling with MYJ PBL EETA=-QKMS*RAS*(QVATM/(1.+QVATM) - QSG/(1.+QSG))*1.E3 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'MYJ EETA',eeta ENDIF else ! myj !-- actual moisture flux from RUC LSM DEW=QKMS*(QVATM-QSG) EETA= - RHO*DEW IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'RUC LSM EETA',eeta ENDIF endif ! myj QFX= XLS*EETA EETA= - RHO*DEW ELSE ! --- evaporation if(myj) then !-- moisture flux for coupling with MYJ PBL EETA=-QKMS*RAS*(QVATM/(1.+QVATM) - QVG/(1.+QVG))*1.E3 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'MYJ EETA',eeta ENDIF else ! myj ! to convert from m s-1 to kg m-2 s-1: *rho water=1.e3************ !-- actual moisture flux from RUC LSM EETA = Q1*1.E3 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'RUC LSM EETA',eeta ENDIF endif ! myj QFX= XLS * EETA EETA = Q1*1.E3 ENDIF EVAPL=EETA S=THDIFICE(1)*CAPICE(1)*DZSTOP*(TSO(1)-TSO(2)) ! heat storage in surface layer SNOH=0. ! There is ice melt X= (cp*rho*r211+rhcs*zsmain(2)*0.5/delt)*(SOILT-TN) + & XLS*rho*r211*(QSG-QGOLD) X=X & ! "heat" from rain -RAINF*CVW*PRCPMS*(max(273.15,TABS)-SOILT) !-- excess energy spent on sea ice melt icemelt=RNET-XLS*EETA -HFT -S -X IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'icemelt=',icemelt ENDIF FLTOT=RNET-XLS*EETA-HFT-S-X-icemelt IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'SICE - FLTOT,RNET,HFT,QFX,S,SNOH,X=', & FLTOT,RNET,HFT,XLS*EETA,s,icemelt,X ENDIF !------------------------------------------------------------------- END SUBROUTINE SICE !------------------------------------------------------------------- SUBROUTINE SNOWSOIL (spp_lsm,rstochcol,fieldcol_sf,& !--- input variables i,j,isoil,delt,ktau,conflx,nzs,nddzs,nroot, & meltfactor,rhonewsn,SNHEI_CRIT, & ! new ILAND,PRCPMS,RAINF,NEWSNOW,snhei,SNWE,SNOWFRAC, & RHOSN, & PATM,QVATM,QCATM, & GLW,GSW,GSWin,EMISS,RNET,IVGTYP, & QKMS,TKMS,PC,cst,drip,infwater, & rho,vegfrac,alb,znt,lai, & MYJ, & !--- soil fixed fields QWRTZ,rhocs,dqm,qmin,ref,wilt,psis,bclh,ksat, & sat,cn,zsmain,zshalf,DTDZS,DTDZS2,tbq, & !--- constants xlv,CP,rovcp,G0_P,cw,stbolt,TABS, & KQWRTZ,KICE,KWT, & !--- output variables ilnb,snweprint,snheiprint,rsm, & soilmois,tso,smfrkeep,keepfr, & dew,soilt,soilt1,tsnav, & qvg,qsg,qcg,SMELT,SNOH,SNFLX,SNOM, & edir1,ec1,ett1,eeta,qfx,hfx,s,sublim, & prcpl,fltot,runoff1,runoff2,mavail,soilice, & soiliqw,infiltrp ) !*************************************************************** ! Energy and moisture budget for snow, heat diffusion eqns. ! in snow and soil, Richards eqn. for soil covered with snow ! ! DELT - time step (s) ! ktau - numver of time step ! CONFLX - depth of constant flux layer (m) ! J,I - the location of grid point ! IME, JME, NZS - dimensions of the domain ! NROOT - number of levels within the root zone ! PRCPMS - precipitation rate in m/s ! NEWSNOW - pcpn in soilid form (m) ! SNHEI, SNWE - snow height and snow water equivalent (m) ! RHOSN - snow density (kg/m-3) ! PATM - pressure (bar) ! QVATM,QCATM - cloud and water vapor mixing ratio ! at the first atm. level (kg/kg) ! GLW, GSW - incoming longwave and absorbed shortwave ! radiation at the surface (W/m^2) ! EMISS,RNET - emissivity (0-1) of the ground surface and net ! radiation at the surface (W/m^2) ! QKMS - exchange coefficient for water vapor in the ! surface layer (m/s) ! TKMS - exchange coefficient for heat in the surface ! layer (m/s) ! PC - plant coefficient (resistance) (0-1) ! RHO - density of atmosphere near surface (kg/m^3) ! VEGFRAC - greeness fraction (0-1) ! RHOCS - volumetric heat capacity of dry soil (J/m^3/K) ! DQM, QMIN - porosity minus residual soil moisture QMIN (m^3/m^3) ! REF, WILT - field capacity soil moisture and the ! wilting point (m^3/m^3) ! PSIS - matrix potential at saturation (m) ! BCLH - exponent for Clapp-Hornberger parameterization ! KSAT - saturated hydraulic conductivity (m/s) ! SAT - maximum value of water intercepted by canopy (m) ! CN - exponent for calculation of canopy water ! ZSMAIN - main levels in soil (m) ! ZSHALF - middle of the soil layers (m) ! DTDZS,DTDZS2 - dt/(2.*dzshalf*dzmain) and dt/dzshalf in soil ! TBQ - table to define saturated mixing ration ! of water vapor for given temperature and pressure ! ilnb - number of layers in snow ! rsm - liquid water inside snow pack (m) ! SOILMOIS,TSO - soil moisture (m^3/m^3) and temperature (K) ! DEW - dew in (kg/m^2 s) ! SOILT - skin temperature (K) ! SOILT1 - snow temperature at 7.5 cm depth (K) ! TSNAV - average temperature of snow pack (C) ! QSG,QVG,QCG - saturated mixing ratio, mixing ratio of ! water vapor and cloud at the ground ! surface, respectively (kg/kg) ! EDIR1, EC1, ETT1, EETA - direct evaporation, evaporation of ! canopy water, transpiration (kg m-2 s-1) and total ! evaporation in (m s-1). ! QFX, HFX - latent and sensible heat fluxes (W/m^2) ! S - soil heat flux in the top layer (W/m^2) ! SUBLIM - snow sublimation (kg/m^2/s) ! RUNOFF1 - surface runoff (m/s) ! RUNOFF2 - underground runoff (m) ! MAVAIL - moisture availability in the top soil layer (0-1) ! SOILICE - content of soil ice in soil layers (m^3/m^3) ! SOILIQW - lliquid water in soil layers (m^3/m^3) ! INFILTRP - infiltration flux from the top of soil domain (m/s) ! XINET - net long-wave radiation (W/m^2) ! !******************************************************************* IMPLICIT NONE !------------------------------------------------------------------- !--- input variables INTEGER, INTENT(IN ) :: nroot,ktau,nzs , & nddzs !nddzs=2*(nzs-2) INTEGER, INTENT(IN ) :: i,j,isoil REAL, INTENT(IN ) :: DELT,CONFLX,PRCPMS , & RAINF,NEWSNOW,RHONEWSN, & SNHEI_CRIT,meltfactor LOGICAL, INTENT(IN ) :: myj !--- 3-D Atmospheric variables REAL, & INTENT(IN ) :: PATM, & QVATM, & QCATM !--- 2-D variables REAL , & INTENT(IN ) :: GLW, & GSW, & GSWin, & RHO, & PC, & VEGFRAC, & lai, & infwater, & QKMS, & TKMS INTEGER, INTENT(IN ) :: IVGTYP !--- soil properties REAL , & INTENT(IN ) :: RHOCS, & BCLH, & DQM, & KSAT, & PSIS, & QMIN, & QWRTZ, & REF, & SAT, & WILT REAL, INTENT(IN ) :: CN, & CW, & XLV, & G0_P, & KQWRTZ, & KICE, & KWT REAL, DIMENSION(1:NZS), INTENT(IN) :: ZSMAIN, & ZSHALF, & DTDZS2 REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS REAL, DIMENSION(1:5001), INTENT(IN) :: TBQ REAL, DIMENSION(1:NZS), INTENT(IN) :: rstochcol REAL, DIMENSION(1:NZS), INTENT(INOUT) :: fieldcol_sf !--- input/output variables !-------- 3-d soil moisture and temperature REAL, DIMENSION( 1:nzs ) , & INTENT(INOUT) :: TSO, & SOILMOIS, & SMFRKEEP REAL, DIMENSION( 1:nzs ) , & INTENT(INOUT) :: KEEPFR INTEGER, INTENT(INOUT) :: ILAND !-------- 2-d variables REAL , & INTENT(INOUT) :: DEW, & CST, & DRIP, & EDIR1, & EC1, & ETT1, & EETA, & RHOSN, & SUBLIM, & PRCPL, & ALB, & EMISS, & ZNT, & MAVAIL, & QVG, & QSG, & QCG, & QFX, & HFX, & S, & RUNOFF1, & RUNOFF2, & SNWE, & SNHEI, & SMELT, & SNOM, & SNOH, & SNFLX, & SOILT, & SOILT1, & SNOWFRAC, & TSNAV INTEGER, INTENT(INOUT) :: ILNB !-------- 1-d variables REAL, DIMENSION(1:NZS), INTENT(OUT) :: SOILICE, & SOILIQW REAL, INTENT(OUT) :: RSM, & SNWEPRINT, & SNHEIPRINT INTEGER, INTENT(IN) :: spp_lsm !--- Local variables INTEGER :: nzs1,nzs2,k REAL :: INFILTRP, TRANSUM , & SNTH, NEWSN , & TABS, T3, UPFLUX, XINET , & BETA, SNWEPR,EPDT,PP REAL :: CP,rovcp,G0,LV,xlvm,STBOLT,xlmelt,dzstop , & can,epot,fac,fltot,ft,fq,hft , & q1,ras,rhoice,sph , & trans,zn,ci,cvw,tln,tavln,pi , & DD1,CMC2MS,DRYCAN,WETCAN , & INFMAX,RIW,DELTSN,H,UMVEG REAL, DIMENSION(1:NZS) :: transp,cap,diffu,hydro , & thdif,tranf,tav,soilmoism , & soilicem,soiliqwm,detal , & fwsat,lwsat,told,smold REAL :: soiltold, qgold REAL :: RNET, X !----------------------------------------------------------------- cvw=cw XLMELT=3.35E+5 !-- heat of water vapor sublimation XLVm=XLV+XLMELT ! STBOLT=5.670151E-8 !--- SNOW flag -- ISICE ! ILAND=isice !--- DELTSN - is the threshold for splitting the snow layer into 2 layers. !--- With snow density 400 kg/m^3, this threshold is equal to 7.5 cm, !--- equivalent to 0.03 m SNWE. For other snow densities the threshold is !--- computed using SNWE=0.03 m and current snow density. !--- SNTH - the threshold below which the snow layer is combined with !--- the top soil layer. SNTH is computed using snwe=0.016 m, and !--- equals 4 cm for snow density 400 kg/m^3. !save SOILT and QVG soiltold=soilt qgold=qvg x=0. ! increase thinkness of top snow layer from 3 cm SWE to 5 cm SWE ! DELTSN=5.*SNHEI_CRIT ! snth=0.4*SNHEI_CRIT DELTSN=0.05*1.e3/rhosn snth=0.01*1.e3/rhosn ! snth=0.01601*1.e3/rhosn ! if(i.eq.442.and.j.eq.260) then ! print *,'deltsn,snhei,snth',i,j,deltsn,snhei,snth ! ENDIF ! For 2-layer snow model when the snow depth is marginally higher than DELTSN, ! reset DELTSN to half of snow depth. IF(SNHEI.GE.DELTSN+SNTH) THEN if(snhei-deltsn-snth.lt.snth) deltsn=0.5*(snhei-snth) IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'DELTSN is changed,deltsn,snhei,snth',i,j,deltsn,snhei,snth ENDIF ENDIF RHOICE=900. CI=RHOICE*2100. RAS=RHO*1.E-3 RIW=rhoice*1.e-3 ! MAVAIL=1. RSM=0. DO K=1,NZS TRANSP (K)=0. soilmoism (k)=0. soiliqwm (k)=0. soilice (k)=0. soilicem (k)=0. lwsat (k)=0. fwsat (k)=0. tav (k)=0. cap (k)=0. diffu (k)=0. hydro (k)=0. thdif (k)=0. tranf (k)=0. detal (k)=0. told (k)=0. smold (k)=0. ENDDO snweprint=0. snheiprint=0. prcpl=prcpms !*** DELTSN is the depth of the top layer of snow where !*** there is a temperature gradient, the rest of the snow layer !*** is considered to have constant temperature NZS1=NZS-1 NZS2=NZS-2 DZSTOP=1./(zsmain(2)-zsmain(1)) !----- THE CALCULATION OF THERMAL DIFFUSIVITY, DIFFUSIONAL AND --- !----- HYDRAULIC CONDUCTIVITY (SMIRNOVA ET AL. 1996, EQ.2,5,6) --- !tgs - the following loop is added to define the amount of frozen !tgs - water in soil if there is any DO K=1,NZS tln=log(tso(k)/273.15) if(tln.lt.0.) then soiliqw(k)=(dqm+qmin)*(XLMELT* & (tso(k)-273.15)/tso(k)/9.81/psis) & **(-1./bclh)-qmin soiliqw(k)=max(0.,soiliqw(k)) soiliqw(k)=min(soiliqw(k),soilmois(k)) soilice(k)=(soilmois(k)-soiliqw(k))/riw !---- melting and freezing is balanced, soil ice cannot increase if(keepfr(k).eq.1.) then soilice(k)=min(soilice(k),smfrkeep(k)) soiliqw(k)=max(0.,soilmois(k)-soilice(k)*rhoice*1.e-3) endif else soilice(k)=0. soiliqw(k)=soilmois(k) endif ENDDO DO K=1,NZS1 tav(k)=0.5*(tso(k)+tso(k+1)) soilmoism(k)=0.5*(soilmois(k)+soilmois(k+1)) tavln=log(tav(k)/273.15) if(tavln.lt.0.) then soiliqwm(k)=(dqm+qmin)*(XLMELT* & (tav(k)-273.15)/tav(k)/9.81/psis) & **(-1./bclh)-qmin fwsat(k)=dqm-soiliqwm(k) lwsat(k)=soiliqwm(k)+qmin soiliqwm(k)=max(0.,soiliqwm(k)) soiliqwm(k)=min(soiliqwm(k), soilmoism(k)) soilicem(k)=(soilmoism(k)-soiliqwm(k))/riw !---- melting and freezing is balanced, soil ice cannot increase if(keepfr(k).eq.1.) then soilicem(k)=min(soilicem(k), & 0.5*(smfrkeep(k)+smfrkeep(k+1))) soiliqwm(k)=max(0.,soilmoism(k)-soilicem(k)*riw) fwsat(k)=dqm-soiliqwm(k) lwsat(k)=soiliqwm(k)+qmin endif else soilicem(k)=0. soiliqwm(k)=soilmoism(k) lwsat(k)=dqm+qmin fwsat(k)=0. endif ENDDO do k=1,nzs if(soilice(k).gt.0.) then smfrkeep(k)=soilice(k) else smfrkeep(k)=soilmois(k)/riw endif enddo !****************************************************************** ! SOILPROP computes thermal diffusivity, and diffusional and ! hydraulic condeuctivities !****************************************************************** CALL SOILPROP(spp_lsm,rstochcol,fieldcol_sf, & !--- input variables nzs,fwsat,lwsat,tav,keepfr, & soilmois,soiliqw,soilice, & soilmoism,soiliqwm,soilicem, & !--- soil fixed fields QWRTZ,rhocs,dqm,qmin,psis,bclh,ksat, & !--- constants riw,xlmelt,CP,G0_P,cvw,ci, & kqwrtz,kice,kwt, & !--- output variables thdif,diffu,hydro,cap) !******************************************************************** !--- CALCULATION OF CANOPY WATER (Smirnova et al., 1996, EQ.16) AND DEW ! DRIP=0. SMELT=0. ! DD1=0. H=1. FQ=QKMS !--- If vegfrac.ne.0. then part of falling snow can be !--- intercepted by the canopy. DEW=0. UMVEG=1.-vegfrac EPOT = -FQ*(QVATM-QSG) ! IF(vegfrac.EQ.0.) then ! cst=0. ! drip=0. ! ELSE ! IF(EPOT.GE.0.) THEN !! DD1=CST+ NEWSNOW*RHOnewSN*1.E-3 & !! -DELT*RAS*EPOT*(CST/SAT)**CN ! IF ( 1==2 ) THEN ! print *,'DD1=',dd1,cst,sat ! ENDIF ! ELSE !! Sublimation ! DEW = - EPOT !! DD1=CST+(NEWSNOW*RHOnewSN*1.E-3+delt*DEW*RAS) ! IF ( 1==2 ) THEN ! print *,'Sublimation DD1=',dd1,cst,sat ! ENDIF ! ENDIF ! DD1=CST+NEWSNOW*RHOnewSN*1.E-3 !! IF(DD1.LT.0.) DD1=0. ! IF (vegfrac.GT.0.) THEN ! CST=DD1 ! IF(CST.GT.SAT) THEN ! CST=SAT ! DRIP=DD1-SAT ! ENDIF ! ENDIF !--- With vegetation part of NEWSNOW can be intercepted by canopy until !--- the saturation is reached. After the canopy saturation is reached !--- DRIP in the solid form will be added to SNOW cover. ! SNWE=SNWE-vegfrac*NEWSNOW*RHOnewSN*1.E-3 & ! SNWE=SNHEI*RHOSN*1.e-3-vegfrac*NEWSNOW*RHOnewSN*1.E-3 & ! + DRIP IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'SNWE after subtracting intercepted snow - snwe=',snwe,vegfrac,cst ENDIF ! ENDIF ! vegfrac=0. ! DRIP=0. ! SNHEI=SNWE*1.e3/RHOSN SNWEPR=SNWE ! check if all snow can evaporate during DT BETA=1. EPDT = EPOT * RAS *DELT*UMVEG IF(EPDT.gt.0. .and. SNWEPR.LE.EPDT) THEN BETA=SNWEPR/max(1.e-8,EPDT) SNWE=0. ENDIF WETCAN=min(0.25,(CST/SAT)**CN) ! if(lai > 1.) wetcan=wetcan/lai DRYCAN=1.-WETCAN !************************************************************** ! TRANSF computes transpiration function !************************************************************** CALL TRANSF(i,j, & !--- input variables nzs,nroot,soiliqw,tabs,lai,gswin, & !--- soil fixed fields dqm,qmin,ref,wilt,zshalf,pc,iland, & !--- output variables tranf,transum) !--- Save soil temp and moisture from the beginning of time step do k=1,nzs told(k)=tso(k) smold(k)=soilmois(k) enddo !************************************************************** ! SNOWTEMP solves heat budget and diffusion eqn. in soil !************************************************************** IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *, 'TSO before calling SNOWTEMP: ', tso ENDIF CALL SNOWTEMP( & !--- input variables i,j,iland,isoil, & delt,ktau,conflx,nzs,nddzs,nroot, & snwe,snwepr,snhei,newsnow,snowfrac, & beta,deltsn,snth,rhosn,rhonewsn,meltfactor, & ! add meltfactor PRCPMS,RAINF, & PATM,TABS,QVATM,QCATM, & GLW,GSW,EMISS,RNET, & QKMS,TKMS,PC,rho,vegfrac, & thdif,cap,drycan,wetcan,cst, & tranf,transum,dew,mavail, & !--- soil fixed fields dqm,qmin,psis,bclh, & zsmain,zshalf,DTDZS,tbq, & !--- constants xlvm,CP,rovcp,G0_P,cvw,stbolt, & !--- output variables snweprint,snheiprint,rsm, & tso,soilt,soilt1,tsnav,qvg,qsg,qcg, & smelt,snoh,snflx,s,ilnb,x) !************************************************************************ !--- RECALCULATION OF DEW USING NEW VALUE OF QSG OR TRANSP IF NO DEW DEW=0. ETT1=0. PP=PATM*1.E3 EPOT = -FQ*(QVATM-QSG) IF(EPOT.GT.0.) THEN ! Evaporation DO K=1,NROOT TRANSP(K)=vegfrac*RAS*FQ*(QVATM-QSG) & *tranf(K)*DRYCAN/zshalf(NROOT+1) ! IF(TRANSP(K).GT.0.) TRANSP(K)=0. ETT1=ETT1-TRANSP(K) ENDDO DO k=nroot+1,nzs transp(k)=0. enddo ELSE ! Sublimation DEW=-EPOT DO K=1,NZS TRANSP(K)=0. ENDDO ETT1=0. ENDIF !-- recalculating of frozen water in soil DO K=1,NZS tln=log(tso(k)/273.15) if(tln.lt.0.) then soiliqw(k)=(dqm+qmin)*(XLMELT* & (tso(k)-273.15)/tso(k)/9.81/psis) & **(-1./bclh)-qmin soiliqw(k)=max(0.,soiliqw(k)) soiliqw(k)=min(soiliqw(k),soilmois(k)) soilice(k)=(soilmois(k)-soiliqw(k))/riw !---- melting and freezing is balanced, soil ice cannot increase if(keepfr(k).eq.1.) then soilice(k)=min(soilice(k),smfrkeep(k)) soiliqw(k)=max(0.,soilmois(k)-soilice(k)*riw) endif else soilice(k)=0. soiliqw(k)=soilmois(k) endif ENDDO !************************************************************************* !--- TQCAN FOR SOLUTION OF MOISTURE BALANCE (Smirnova et al. 1996, EQ.22,28) ! AND TSO,ETA PROFILES !************************************************************************* CALL SOILMOIST ( & !-- input delt,nzs,nddzs,DTDZS,DTDZS2,RIW, & zsmain,zshalf,diffu,hydro, & QSG,QVG,QCG,QCATM,QVATM,-INFWATER, & QKMS,TRANSP,0., & 0.,SMELT,soilice,vegfrac, & snowfrac,1., & !-- soil properties DQM,QMIN,REF,KSAT,RAS,INFMAX, & !-- output SOILMOIS,SOILIQW,MAVAIL,RUNOFF1, & RUNOFF2,infiltrp) ! endif !-- Restore land-use parameters if all snow is melted IF(SNHEI.EQ.0.) then tsnav=soilt-273.15 ENDIF ! 21apr2009 ! SNOM [mm] goes into the passed-in ACSNOM variable in the grid derived type SNOM=SNOM+SMELT*DELT*1.e3 ! !--- KEEPFR is 1 when the temperature and moisture in soil !--- are both increasing. In this case soil ice should not !--- be increasing according to the freezing curve. !--- Some part of ice is melted, but additional water is !--- getting frozen. Thus, only structure of frozen soil is !--- changed, and phase changes are not affecting the heat !--- transfer. This situation may happen when it rains on the !--- frozen soil. do k=1,nzs if (soilice(k).gt.0.) then if(tso(k).gt.told(k).and.soilmois(k).gt.smold(k)) then keepfr(k)=1. else keepfr(k)=0. endif endif enddo !--- THE DIAGNOSTICS OF SURFACE FLUXES T3 = STBOLT*SOILTold*SOILTold*SOILTold UPFLUX = T3 *0.5*(SOILTold+SOILT) XINET = EMISS*(GLW-UPFLUX) ! RNET = GSW + XINET HFX=-TKMS*CP*RHO*(TABS-SOILT) & *(P1000mb*0.00001/Patm)**ROVCP IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'potential temp HFX',hfx ENDIF HFT=-TKMS*CP*RHO*(TABS-SOILT) IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'abs temp HFX',hft ENDIF Q1 = - FQ*RAS* (QVATM - QSG) CMC2MS=0. IF (Q1.LT.0.) THEN ! --- condensation EDIR1=0. EC1=0. ETT1=0. ! --- condensation if(myj) then !-- moisture flux for coupling with MYJ PBL EETA=-QKMS*RAS*(QVATM/(1.+QVATM) - QSG/(1.+QSG))*1.E3 CST= CST-EETA*DELT*vegfrac IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'MYJ EETA cond', EETA ENDIF else ! myj !-- actual moisture flux from RUC LSM DEW=QKMS*(QVATM-QSG) EETA= - RHO*DEW CST=CST+DELT*DEW*RAS * vegfrac IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'RUC LSM EETA cond',EETA ENDIF endif ! myj QFX= XLVm*EETA EETA= - RHO*DEW ELSE ! --- evaporation EDIR1 = Q1*UMVEG *BETA CMC2MS=CST/DELT*RAS EC1 = Q1 * WETCAN * vegfrac CST=max(0.,CST-EC1 * DELT) ! if(EC1 > CMC2MS) then ! EC1 = min(cmc2ms,ec1) ! CST = 0. ! endif IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print*,'Q1,umveg,beta',Q1,umveg,beta print *,'wetcan,vegfrac',wetcan,vegfrac print *,'EC1,CMC2MS',EC1,CMC2MS ENDIF if(myj) then !-- moisture flux for coupling with MYJ PBL EETA=-(QKMS*RAS*(QVATM/(1.+QVATM) - QSG/(1.+QSG))*1.E3)*BETA IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'MYJ EETA', EETA*XLVm,EETA ENDIF else ! myj ! to convert from m s-1 to kg m-2 s-1: *rho water=1.e3************ !-- actual moisture flux from RUC LSM EETA = (EDIR1 + EC1 + ETT1)*1.E3 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'RUC LSM EETA',EETA*XLVm,EETA ENDIF endif ! myj QFX= XLVm * EETA EETA = (EDIR1 + EC1 + ETT1)*1.E3 ENDIF S=SNFLX ! sublim=eeta sublim=EDIR1*1.E3 ! Energy budget FLTOT=RNET-HFT-XLVm*EETA-S-SNOH-x IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'SNOWSOIL - FLTOT,RNET,HFT,QFX,S,SNOH,X=',FLTOT,RNET,HFT,XLVm*EETA,s,SNOH,X print *,'edir1,ec1,ett1,mavail,qkms,qvatm,qvg,qsg,vegfrac,beta',& edir1,ec1,ett1,mavail,qkms,qvatm,qvg,qsg,vegfrac,beta ENDIF 222 CONTINUE 1123 FORMAT(I5,8F12.3) 1133 FORMAT(I7,8E12.4) 123 format(i6,f6.2,7f8.1) 122 FORMAT(1X,2I3,6F8.1,F8.3,F8.2) !------------------------------------------------------------------- END SUBROUTINE SNOWSOIL !------------------------------------------------------------------- SUBROUTINE SNOWSEAICE( & i,j,isoil,delt,ktau,conflx,nzs,nddzs, & meltfactor,rhonewsn,SNHEI_CRIT, & ! new ILAND,PRCPMS,RAINF,NEWSNOW,snhei,SNWE,snowfrac, & RHOSN,PATM,QVATM,QCATM, & GLW,GSW,EMISS,RNET, & QKMS,TKMS,RHO,myj, & !--- sea ice parameters ALB,ZNT, & tice,rhosice,capice,thdifice, & zsmain,zshalf,DTDZS,DTDZS2,tbq, & !--- constants xlv,CP,rovcp,cw,stbolt,tabs, & !--- output variables ilnb,snweprint,snheiprint,rsm,tso, & dew,soilt,soilt1,tsnav,qvg,qsg,qcg, & SMELT,SNOH,SNFLX,SNOM,eeta, & qfx,hfx,s,sublim,prcpl,fltot & ) !*************************************************************** ! Solving energy budget for snow on sea ice and heat diffusion ! eqns. in snow and sea ice !*************************************************************** IMPLICIT NONE !------------------------------------------------------------------- !--- input variables INTEGER, INTENT(IN ) :: ktau,nzs , & nddzs !nddzs=2*(nzs-2) INTEGER, INTENT(IN ) :: i,j,isoil REAL, INTENT(IN ) :: DELT,CONFLX,PRCPMS , & RAINF,NEWSNOW,RHONEWSN, & meltfactor, snhei_crit real :: rhonewcsn LOGICAL, INTENT(IN ) :: myj !--- 3-D Atmospheric variables REAL, & INTENT(IN ) :: PATM, & QVATM, & QCATM !--- 2-D variables REAL , & INTENT(IN ) :: GLW, & GSW, & RHO, & QKMS, & TKMS !--- sea ice properties REAL, DIMENSION(1:NZS) , & INTENT(IN ) :: & tice, & rhosice, & capice, & thdifice REAL, INTENT(IN ) :: & CW, & XLV REAL, DIMENSION(1:NZS), INTENT(IN) :: ZSMAIN, & ZSHALF, & DTDZS2 REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS REAL, DIMENSION(1:5001), INTENT(IN) :: TBQ !--- input/output variables !-------- 3-d soil moisture and temperature REAL, DIMENSION( 1:nzs ) , & INTENT(INOUT) :: TSO INTEGER, INTENT(INOUT) :: ILAND !-------- 2-d variables REAL , & INTENT(INOUT) :: DEW, & EETA, & RHOSN, & SUBLIM, & PRCPL, & ALB, & EMISS, & ZNT, & QVG, & QSG, & QCG, & QFX, & HFX, & S, & SNWE, & SNHEI, & SMELT, & SNOM, & SNOH, & SNFLX, & SOILT, & SOILT1, & SNOWFRAC, & TSNAV INTEGER, INTENT(INOUT) :: ILNB REAL, INTENT(OUT) :: RSM, & SNWEPRINT, & SNHEIPRINT !--- Local variables INTEGER :: nzs1,nzs2,k,k1,kn,kk REAL :: x,x1,x2,dzstop,ft,tn,denom REAL :: SNTH, NEWSN , & TABS, T3, UPFLUX, XINET , & BETA, SNWEPR,EPDT,PP REAL :: CP,rovcp,G0,LV,xlvm,STBOLT,xlmelt , & epot,fltot,fq,hft,q1,ras,rhoice,ci,cvw , & RIW,DELTSN,H REAL :: rhocsn,thdifsn, & xsn,ddzsn,x1sn,d1sn,d2sn,d9sn,r22sn REAL :: cotsn,rhtsn,xsn1,ddzsn1,x1sn1,ftsnow,denomsn REAL :: fso,fsn, & FKT,D1,D2,D9,D10,DID,R211,R21,R22,R6,R7,D11, & FKQ,R210,AA,BB,QS1,TS1,TQ2,TX2, & TDENOM,AA1,RHCS,H1,TSOB, SNPRIM, & SNODIF,SOH,TNOLD,QGOLD,SNOHGNEW REAL, DIMENSION(1:NZS) :: cotso,rhtso REAL :: RNET,rsmfrac,soiltfrac,hsn,icemelt,rr integer :: nmelt !----------------------------------------------------------------- XLMELT=3.35E+5 !-- heat of sublimation of water vapor XLVm=XLV+XLMELT ! STBOLT=5.670151E-8 !--- SNOW flag -- ISICE ! ILAND=isice !--- DELTSN - is the threshold for splitting the snow layer into 2 layers. !--- With snow density 400 kg/m^3, this threshold is equal to 7.5 cm, !--- equivalent to 0.03 m SNWE. For other snow densities the threshold is !--- computed using SNWE=0.03 m and current snow density. !--- SNTH - the threshold below which the snow layer is combined with !--- the top sea ice layer. SNTH is computed using snwe=0.016 m, and !--- equals 4 cm for snow density 400 kg/m^3. ! increase thickness of top snow layer from 3 cm SWE to 5 cm SWE ! DELTSN=5.*SNHEI_CRIT ! snth=0.4*SNHEI_CRIT DELTSN=0.05*1.e3/rhosn snth=0.01*1.e3/rhosn ! snth=0.01601*1.e3/rhosn ! For 2-layer snow model when the snow depth is marginlly higher than DELTSN, ! reset DELTSN to half of snow depth. IF(SNHEI.GE.DELTSN+SNTH) THEN if(snhei-deltsn-snth.lt.snth) deltsn=0.5*(snhei-snth) IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'DELTSN ICE is changed,deltsn,snhei,snth', & i,j, deltsn,snhei,snth ENDIF ENDIF RHOICE=900. CI=RHOICE*2100. RAS=RHO*1.E-3 RIW=rhoice*1.e-3 RSM=0. XLMELT=3.35E+5 RHOCSN=2090.* RHOSN !18apr08 - add rhonewcsn RHOnewCSN=2090.* RHOnewSN THDIFSN = 0.265/RHOCSN RAS=RHO*1.E-3 SOILTFRAC=SOILT SMELT=0. SOH=0. SNODIF=0. SNOH=0. SNOHGNEW=0. RSM = 0. RSMFRAC = 0. fsn=1. fso=0. cvw=cw NZS1=NZS-1 NZS2=NZS-2 QGOLD=QSG TNOLD=SOILT DZSTOP=1./(ZSMAIN(2)-ZSMAIN(1)) snweprint=0. snheiprint=0. prcpl=prcpms !*** DELTSN is the depth of the top layer of snow where !*** there is a temperature gradient, the rest of the snow layer !*** is considered to have constant temperature H=1. SMELT=0. FQ=QKMS SNHEI=SNWE*1.e3/RHOSN SNWEPR=SNWE ! check if all snow can evaporate during DT BETA=1. EPOT = -FQ*(QVATM-QSG) EPDT = EPOT * RAS *DELT IF(EPDT.GT.0. .and. SNWEPR.LE.EPDT) THEN BETA=SNWEPR/max(1.e-8,EPDT) SNWE=0. ENDIF !****************************************************************************** ! COEFFICIENTS FOR THOMAS ALGORITHM FOR TSO !****************************************************************************** cotso(1)=0. rhtso(1)=TSO(NZS) DO 33 K=1,NZS2 KN=NZS-K K1=2*KN-3 X1=DTDZS(K1)*THDIFICE(KN-1) X2=DTDZS(K1+1)*THDIFICE(KN) FT=TSO(KN)+X1*(TSO(KN-1)-TSO(KN)) & -X2*(TSO(KN)-TSO(KN+1)) DENOM=1.+X1+X2-X2*cotso(K) cotso(K+1)=X1/DENOM rhtso(K+1)=(FT+X2*rhtso(K))/DENOM 33 CONTINUE !--- THE NZS element in COTSO and RHTSO will be for snow !--- There will be 2 layers in snow if it is deeper than DELTSN+SNTH IF(SNHEI.GE.SNTH) then if(snhei.le.DELTSN+SNTH) then !-- 1-layer snow model ilnb=1 snprim=max(snth,snhei) soilt1=tso(1) tsob=tso(1) XSN = DELT/2./(zshalf(2)+0.5*SNPRIM) DDZSN = XSN / SNPRIM X1SN = DDZSN * thdifsn X2 = DTDZS(1)*THDIFICE(1) FT = TSO(1)+X1SN*(SOILT-TSO(1)) & -X2*(TSO(1)-TSO(2)) DENOM = 1. + X1SN + X2 -X2*cotso(NZS1) cotso(NZS)=X1SN/DENOM rhtso(NZS)=(FT+X2*rhtso(NZS1))/DENOM cotsn=cotso(NZS) rhtsn=rhtso(NZS) !*** Average temperature of snow pack (C) tsnav=0.5*(soilt+tso(1)) & -273.15 else !-- 2 layers in snow, SOILT1 is temperasture at DELTSN depth ilnb=2 snprim=deltsn tsob=soilt1 XSN = DELT/2./(0.5*SNHEI) XSN1= DELT/2./(zshalf(2)+0.5*(SNHEI-DELTSN)) DDZSN = XSN / DELTSN DDZSN1 = XSN1 / (SNHEI-DELTSN) X1SN = DDZSN * thdifsn X1SN1 = DDZSN1 * thdifsn X2 = DTDZS(1)*THDIFICE(1) FT = TSO(1)+X1SN1*(SOILT1-TSO(1)) & -X2*(TSO(1)-TSO(2)) DENOM = 1. + X1SN1 + X2 - X2*cotso(NZS1) cotso(nzs)=x1sn1/denom rhtso(nzs)=(ft+x2*rhtso(nzs1))/denom ftsnow = soilt1+x1sn*(soilt-soilt1) & -x1sn1*(soilt1-tso(1)) denomsn = 1. + X1SN + X1SN1 - X1SN1*cotso(NZS) cotsn=x1sn/denomsn rhtsn=(ftsnow+X1SN1*rhtso(NZS))/denomsn !*** Average temperature of snow pack (C) tsnav=0.5/snhei*((soilt+soilt1)*deltsn & +(soilt1+tso(1))*(SNHEI-DELTSN)) & -273.15 endif ENDIF IF(SNHEI.LT.SNTH.AND.SNHEI.GT.0.) then !--- snow is too thin to be treated separately, therefore it !--- is combined with the first sea ice layer. snprim=SNHEI+zsmain(2) fsn=SNHEI/snprim fso=1.-fsn soilt1=tso(1) tsob=tso(2) XSN = DELT/2./((zshalf(3)-zsmain(2))+0.5*snprim) DDZSN = XSN /snprim X1SN = DDZSN * (fsn*thdifsn+fso*thdifice(1)) X2=DTDZS(2)*THDIFICE(2) FT=TSO(2)+X1SN*(SOILT-TSO(2))- & X2*(TSO(2)-TSO(3)) denom = 1. + x1sn + x2 - x2*cotso(nzs-2) cotso(nzs1) = x1sn/denom rhtso(nzs1)=(FT+X2*rhtso(NZS-2))/denom tsnav=0.5*(soilt+tso(1)) & -273.15 cotso(nzs)=cotso(NZS1) rhtso(nzs)=rhtso(nzs1) cotsn=cotso(NZS) rhtsn=rhtso(NZS) ENDIF !************************************************************************ !--- THE HEAT BALANCE EQUATION !18apr08 nmelt is the flag for melting, and SNOH is heat of snow phase changes nmelt=0 SNOH=0. EPOT=-QKMS*(QVATM-QSG) RHCS=CAPICE(1) H=1. FKT=TKMS D1=cotso(NZS1) D2=rhtso(NZS1) TN=SOILT D9=THDIFICE(1)*RHCS*dzstop D10=TKMS*CP*RHO R211=.5*CONFLX/DELT R21=R211*CP*RHO R22=.5/(THDIFICE(1)*DELT*dzstop**2) R6=EMISS *STBOLT*.5*TN**4 R7=R6/TN D11=RNET+R6 IF(SNHEI.GE.SNTH) THEN if(snhei.le.DELTSN+SNTH) then !--- 1-layer snow D1SN = cotso(NZS) D2SN = rhtso(NZS) else !--- 2-layer snow D1SN = cotsn D2SN = rhtsn endif D9SN= THDIFSN*RHOCSN / SNPRIM R22SN = SNPRIM*SNPRIM*0.5/(THDIFSN*DELT) ENDIF IF(SNHEI.LT.SNTH.AND.SNHEI.GT.0.) then !--- thin snow is combined with sea ice D1SN = D1 D2SN = D2 D9SN = (fsn*THDIFSN*RHOCSN+fso*THDIFICE(1)*RHCS)/ & snprim R22SN = snprim*snprim*0.5 & /((fsn*THDIFSN+fso*THDIFICE(1))*delt) ENDIF IF(SNHEI.eq.0.)then !--- all snow is sublimated D9SN = D9 R22SN = R22 D1SN = D1 D2SN = D2 ENDIF !---- TDENOM for snow TDENOM = D9SN*(1.-D1SN +R22SN)+D10+R21+R7 & +RAINF*CVW*PRCPMS & +RHOnewCSN*NEWSNOW/DELT FKQ=QKMS*RHO R210=R211*RHO AA=XLVM*(BETA*FKQ+R210)/TDENOM BB=(D10*TABS+R21*TN+XLVM*(QVATM* & (BETA*FKQ) & +R210*QVG)+D11+D9SN*(D2SN+R22SN*TN) & +RAINF*CVW*PRCPMS*max(273.15,TABS) & + RHOnewCSN*NEWSNOW/DELT*min(273.15,TABS) & )/TDENOM AA1=AA PP=PATM*1.E3 AA1=AA1/PP !18apr08 - the iteration start point 212 continue BB=BB-SNOH/TDENOM IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'VILKA-SNOW on SEAICE' print *,'tn,aa1,bb,pp,fkq,r210', & tn,aa1,bb,pp,fkq,r210 print *,'TABS,QVATM,TN,QVG=',TABS,QVATM,TN,QVG ENDIF CALL VILKA(TN,AA1,BB,PP,QS1,TS1,TBQ,KTAU,i,j,iland,isoil) !--- it is saturation over snow QVG=QS1 QSG=QS1 QCG=0. !--- SOILT - skin temperature SOILT=TS1 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,' AFTER VILKA-SNOW on SEAICE' print *,' TS1,QS1: ', ts1,qs1 ENDIF ! Solution for temperature at 7.5 cm depth and snow-seaice interface IF(SNHEI.GE.SNTH) THEN if(snhei.gt.DELTSN+SNTH) then !-- 2-layer snow model SOILT1=min(273.15,rhtsn+cotsn*SOILT) TSO(1)=min(271.4,(rhtso(NZS)+cotso(NZS)*SOILT1)) tsob=soilt1 else !-- 1 layer in snow TSO(1)=min(271.4,(rhtso(NZS)+cotso(NZS)*SOILT)) SOILT1=TSO(1) tsob=tso(1) endif ELSEIF (SNHEI > 0. .and. SNHEI < SNTH) THEN ! blended TSO(2)=min(271.4,(rhtso(NZS1)+cotso(NZS1)*SOILT)) tso(1)=min(271.4,(tso(2)+(soilt-tso(2))*fso)) SOILT1=TSO(1) tsob=TSO(2) ELSE ! snow is melted TSO(1)=min(271.4,SOILT) SOILT1=min(271.4,SOILT) tsob=tso(1) ENDIF !---- Final solution for TSO in sea ice IF (SNHEI > 0. .and. SNHEI < SNTH) THEN ! blended or snow is melted DO K=3,NZS KK=NZS-K+1 TSO(K)=min(271.4,rhtso(KK)+cotso(KK)*TSO(K-1)) END DO ELSE DO K=2,NZS KK=NZS-K+1 TSO(K)=min(271.4,rhtso(KK)+cotso(KK)*TSO(K-1)) END DO ENDIF !--- For thin snow layer combined with the top soil layer !--- TSO(i,j,1) is computed by linear interpolation between SOILT !--- and TSO(i,j,2) ! if(SNHEI.LT.SNTH.AND.SNHEI.GT.0.)then ! tso(1)=min(271.4,tso(2)+(soilt-tso(2))*fso) ! soilt1=tso(1) ! tsob = tso(2) ! endif if(nmelt.eq.1) go to 220 !--- IF SOILT > 273.15 F then melting of snow can happen ! IF(SOILT.GT.273.15.AND.SNWE.GT.0.) THEN ! if all snow can evaporate, then there is nothing to melt IF(SOILT.GT.273.15.AND.SNWEPR-BETA*EPOT*RAS*DELT.GT.0..AND.SNHEI.GT.0.) THEN ! nmelt = 1 ! soiltfrac=273.15 soiltfrac=snowfrac*273.15+(1.-snowfrac)*min(271.4,SOILT) QSG= QSN(soiltfrac,TBQ)/PP T3 = STBOLT*TNold*TNold*TNold UPFLUX = T3 * 0.5*(TNold+SOILTfrac) XINET = EMISS*(GLW-UPFLUX) ! RNET = GSW + XINET EPOT = -QKMS*(QVATM-QSG) Q1=EPOT*RAS IF (Q1.LE.0.) THEN ! --- condensation DEW=-EPOT QFX= XLVM*RHO*DEW EETA=QFX/XLVM ELSE ! --- evaporation EETA = Q1 * BETA *1.E3 ! to convert from kg m-2 s-1 to m s-1: 1/rho water=1.e-3************ QFX= - XLVM * EETA ENDIF HFX=D10*(TABS-soiltfrac) IF(SNHEI.GE.SNTH)then SOH=thdifsn*RHOCSN*(soiltfrac-TSOB)/SNPRIM SNFLX=SOH ELSE SOH=(fsn*thdifsn*rhocsn+fso*thdifice(1)*rhcs)* & (soiltfrac-TSOB)/snprim SNFLX=SOH ENDIF X= (R21+D9SN*R22SN)*(soiltfrac-TNOLD) + & XLVM*R210*(QSG-QGOLD) !-- SNOH is energy flux of snow phase change SNOH=RNET+QFX +HFX & +RHOnewCSN*NEWSNOW/DELT*(min(273.15,TABS)-soiltfrac) & -SOH-X+RAINF*CVW*PRCPMS* & (max(273.15,TABS)-soiltfrac) IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'SNOWSEAICE melt I,J,SNOH,RNET,QFX,HFX,SOH,X',i,j,SNOH,RNET,QFX,HFX,SOH,X print *,'RHOnewCSN*NEWSNOW/DELT*(min(273.15,TABS)-soiltfrac)', & RHOnewCSN*NEWSNOW/DELT*(min(273.15,TABS)-soiltfrac) print *,'RAINF*CVW*PRCPMS*(max(273.15,TABS)-soiltfrac)', & RAINF*CVW*PRCPMS*(max(273.15,TABS)-soiltfrac) ENDIF SNOH=AMAX1(0.,SNOH) !-- SMELT is speed of melting in M/S SMELT= SNOH /XLMELT*1.E-3 SMELT=AMIN1(SMELT,SNWEPR/DELT-BETA*EPOT*RAS) SMELT=AMAX1(0.,SMELT) IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'1-SMELT i,j',smelt,i,j ENDIF !18apr08 - Egglston limit ! SMELT= amin1 (smelt, 5.6E-7*meltfactor*max(1.,(soilt-273.15))) SMELT= amin1 (smelt, 5.6E-8*meltfactor*max(1.,(soilt-273.15))) IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'2-SMELT i,j',smelt,i,j ENDIF ! rr - potential melting rr=SNWEPR/delt-BETA*EPOT*RAS SMELT=min(SMELT,rr) IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'3- SMELT i,j,smelt,rr',i,j,smelt,rr ENDIF SNOHGNEW=SMELT*XLMELT*1.E3 SNODIF=AMAX1(0.,(SNOH-SNOHGNEW)) SNOH=SNOHGNEW IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print*,'soiltfrac,soilt,SNOHGNEW,SNODIF=', & i,j,soiltfrac,soilt,snohgnew,snodif print *,'SNOH,SNODIF',SNOH,SNODIF ENDIF !*** From Koren et al. (1999) 13% of snow melt stays in the snow pack rsmfrac=min(0.18,(max(0.08,snwepr/0.10*0.13))) if(snhei > 0.01) then rsm=rsmfrac*smelt*delt else ! do not keep melted water if snow depth is less that 1 cm rsm=0. endif !18apr08 rsm is part of melted water that stays in snow as liquid SMELT=AMAX1(0.,SMELT-rsm/delt) IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'4-SMELT i,j,smelt,rsm,snwepr,rsmfrac', & i,j,smelt,rsm,snwepr,rsmfrac ENDIF !-- update liquid equivalent of snow depth !-- for evaporation and snow melt SNWE = AMAX1(0.,(SNWEPR- & (SMELT+BETA*EPOT*RAS)*DELT & ! (SMELT+BETA*EPOT*RAS)*DELT*snowfrac & ) ) !!!! soilt=soiltfrac !--- If there is no snow melting then just evaporation !--- or condensation changes SNWE ELSE if(snhei.ne.0.) then EPOT=-QKMS*(QVATM-QSG) SNWE = AMAX1(0.,(SNWEPR- & BETA*EPOT*RAS*DELT)) ! BETA*EPOT*RAS*DELT*snowfrac)) endif ENDIF ! no iteration for snow on sea ice, because it will produce ! skin temperature higher than it is possible with snow on sea ice ! if(nmelt.eq.1) goto 212 ! second iteration 220 continue if(smelt > 0..and. rsm > 0.) then if(snwe.le.rsm) then IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'SEAICE SNWEQVATM .and. QVATM > QVG) then ! very dry soil ! print *,'very dry soils mavail,qvg,qs1,qvatm,ts1',i,j,mavail,qvg,qs1,qvatm,ts1 ! QVG = QVATM ! endif TSO(1)=TS1 QCG=0. IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN ! if (i==421.and.j==280) then print *,'q1,qsg,qvg,qvatm,alfa,h',q1,qsg,qvg,qvatm,alfa,h endif 200 CONTINUE IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'200 QVG,QSG,QCG,TSO(1)',QVG,QSG,QCG,TSO(1) ENDIF if(1==2) then !--11Aug18 - do not use this iteration, it warms up skin temperature too much in summer time at night. if(qvatm > QSG .and. iter==0) then !condensation regime IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'turn off canopy evaporation and transpiration' print *,' QVATM,QVG,QSG,TS1',QVATM,QVG,QSG,TS1 ENDIF can=0. umveg=1. iter=1 goto 2111 endif IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN if(iter == 1) then print *,'QVATM,QVG,QSG,QCG,TS1',QVATM,QVG,QSG,QCG,TS1 endif ENDIF endif ! 1==2 !--- SOILT - skin temperature SOILT=TS1 !---- Final solution for soil temperature - TSO DO K=2,NZS KK=NZS-K+1 TSO(K)=rhtso(KK)+cotso(KK)*TSO(K-1) END DO X= (cp*rho*r211+rhcs*zsmain(2)*0.5/delt)*(SOILT-TN) + & XLV*rho*r211*(QVG-QGOLD) ! IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print*,'SOILTEMP storage, i,j,x,soilt,tn,qvg,qvgold', & i,j,x,soilt,tn,qvg,qgold print *,'TEMP term (cp*rho*r211+rhcs*zsmain(2)*0.5/delt)*(SOILT-TN)',& (cp*rho*r211+rhcs*zsmain(2)*0.5/delt)*(SOILT-TN) print *,'QV term XLV*rho*r211*(QVG-QGOLD)',XLV*rho*r211*(QVG-QGOLD) ENDIF X=X & ! "heat" from rain -RAINF*CVW*PRCPMS*(max(273.15,TABS)-SOILT) IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'x=',x ENDIF !-------------------------------------------------------------------- END SUBROUTINE SOILTEMP !-------------------------------------------------------------------- SUBROUTINE SNOWTEMP( & !--- input variables i,j,iland,isoil, & delt,ktau,conflx,nzs,nddzs,nroot, & snwe,snwepr,snhei,newsnow,snowfrac, & beta,deltsn,snth,rhosn,rhonewsn,meltfactor, & ! add meltfactor PRCPMS,RAINF, & PATM,TABS,QVATM,QCATM, & GLW,GSW,EMISS,RNET, & QKMS,TKMS,PC,RHO,VEGFRAC, & THDIF,CAP,DRYCAN,WETCAN,CST, & TRANF,TRANSUM,DEW,MAVAIL, & !--- soil fixed fields DQM,QMIN,PSIS,BCLH, & ZSMAIN,ZSHALF,DTDZS,TBQ, & !--- constants XLVM,CP,rovcp,G0_P,CVW,STBOLT, & !--- output variables SNWEPRINT,SNHEIPRINT,RSM, & TSO,SOILT,SOILT1,TSNAV,QVG,QSG,QCG, & SMELT,SNOH,SNFLX,S,ILNB,X) !******************************************************************** ! Energy budget equation and heat diffusion eqn are ! solved here to obtain snow and soil temperatures ! ! DELT - time step (s) ! ktau - numver of time step ! CONFLX - depth of constant flux layer (m) ! IME, JME, KME, NZS - dimensions of the domain ! NROOT - number of levels within the root zone ! PRCPMS - precipitation rate in m/s ! COTSO, RHTSO - coefficients for implicit solution of ! heat diffusion equation ! THDIF - thermal diffusivity (W/m/K) ! QSG,QVG,QCG - saturated mixing ratio, mixing ratio of ! water vapor and cloud at the ground ! surface, respectively (kg/kg) ! PATM - pressure [bar] ! QCATM,QVATM - cloud and water vapor mixing ratio ! at the first atm. level (kg/kg) ! EMISS,RNET - emissivity (0-1) of the ground surface and net ! radiation at the surface (W/m^2) ! QKMS - exchange coefficient for water vapor in the ! surface layer (m/s) ! TKMS - exchange coefficient for heat in the surface ! layer (m/s) ! PC - plant coefficient (resistance) ! RHO - density of atmosphere near surface (kg/m^3) ! VEGFRAC - greeness fraction (0-1) ! CAP - volumetric heat capacity (J/m^3/K) ! DRYCAN - dry fraction of vegetated area where ! transpiration may take place (0-1) ! WETCAN - fraction of vegetated area covered by canopy ! water (0-1) ! TRANSUM - transpiration function integrated over the ! rooting zone (m) ! DEW - dew in kg/m^2/s ! MAVAIL - fraction of maximum soil moisture in the top ! layer (0-1) ! ZSMAIN - main levels in soil (m) ! ZSHALF - middle of the soil layers (m) ! DTDZS - dt/(2.*dzshalf*dzmain) ! TBQ - table to define saturated mixing ration ! of water vapor for given temperature and pressure ! TSO - soil temperature (K) ! SOILT - skin temperature (K) ! !********************************************************************* IMPLICIT NONE !--------------------------------------------------------------------- !--- input variables INTEGER, INTENT(IN ) :: nroot,ktau,nzs , & nddzs !nddzs=2*(nzs-2) INTEGER, INTENT(IN ) :: i,j,iland,isoil REAL, INTENT(IN ) :: DELT,CONFLX,PRCPMS , & RAINF,NEWSNOW,DELTSN,SNTH , & TABS,TRANSUM,SNWEPR , & rhonewsn,meltfactor real :: rhonewcsn !--- 3-D Atmospheric variables REAL, & INTENT(IN ) :: PATM, & QVATM, & QCATM !--- 2-D variables REAL , & INTENT(IN ) :: GLW, & GSW, & RHO, & PC, & VEGFRAC, & QKMS, & TKMS !--- soil properties REAL , & INTENT(IN ) :: & BCLH, & DQM, & PSIS, & QMIN REAL, INTENT(IN ) :: CP, & ROVCP, & CVW, & STBOLT, & XLVM, & G0_P REAL, DIMENSION(1:NZS), INTENT(IN) :: ZSMAIN, & ZSHALF, & THDIF, & CAP, & TRANF REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS REAL, DIMENSION(1:5001), INTENT(IN) :: TBQ !--- input/output variables !-------- 3-d soil moisture and temperature REAL, DIMENSION( 1:nzs ) , & INTENT(INOUT) :: TSO !-------- 2-d variables REAL , & INTENT(INOUT) :: DEW, & CST, & RHOSN, & EMISS, & MAVAIL, & QVG, & QSG, & QCG, & SNWE, & SNHEI, & SNOWFRAC, & SMELT, & SNOH, & SNFLX, & S, & SOILT, & SOILT1, & TSNAV REAL, INTENT(INOUT) :: DRYCAN, WETCAN REAL, INTENT(OUT) :: RSM, & SNWEPRINT, & SNHEIPRINT INTEGER, INTENT(OUT) :: ilnb !--- Local variables INTEGER :: nzs1,nzs2,k,k1,kn,kk REAL :: x,x1,x2,x4,dzstop,can,ft,sph, & tn,trans,umveg,denom REAL :: cotsn,rhtsn,xsn1,ddzsn1,x1sn1,ftsnow,denomsn REAL :: t3,upflux,xinet,ras, & xlmelt,rhocsn,thdifsn, & beta,epot,xsn,ddzsn,x1sn,d1sn,d2sn,d9sn,r22sn REAL :: fso,fsn, & FKT,D1,D2,D9,D10,DID,R211,R21,R22,R6,R7,D11, & PI,H,FKQ,R210,AA,BB,PP,Q1,QS1,TS1,TQ2,TX2, & TDENOM,C,CC,AA1,RHCS,H1, & tsob, snprim, sh1, sh2, & smeltg,snohg,snodif,soh, & CMC2MS,TNOLD,QGOLD,SNOHGNEW REAL, DIMENSION(1:NZS) :: transp,cotso,rhtso REAL :: edir1, & ec1, & ett1, & eeta, & qfx, & hfx REAL :: RNET,rsmfrac,soiltfrac,hsn,rr integer :: nmelt, iter !----------------------------------------------------------------- iter = 0 do k=1,nzs transp (k)=0. cotso (k)=0. rhtso (k)=0. enddo IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *, 'SNOWTEMP: SNHEI,SNTH,SOILT1: ',SNHEI,SNTH,SOILT1,soilt ENDIF XLMELT=3.35E+5 RHOCSN=2090.* RHOSN !18apr08 - add rhonewcsn RHOnewCSN=2090.* RHOnewSN THDIFSN = 0.265/RHOCSN RAS=RHO*1.E-3 SOILTFRAC=SOILT SMELT=0. SOH=0. SMELTG=0. SNOHG=0. SNODIF=0. RSM = 0. RSMFRAC = 0. fsn=1. fso=0. ! hsn=snhei NZS1=NZS-1 NZS2=NZS-2 QGOLD=QVG DZSTOP=1./(ZSMAIN(2)-ZSMAIN(1)) !****************************************************************************** ! COEFFICIENTS FOR THOMAS ALGORITHM FOR TSO !****************************************************************************** ! did=2.*(ZSMAIN(nzs)-ZSHALF(nzs)) ! h1=DTDZS(8)*THDIF(nzs-1)*(ZSHALF(nzs)-ZSHALF(nzs-1))/did ! cotso(1)=h1/(1.+h1) ! rhtso(1)=(tso(nzs)+h1*(tso(nzs-1)-tso(nzs)))/ ! 1 (1.+h1) cotso(1)=0. rhtso(1)=TSO(NZS) DO 33 K=1,NZS2 KN=NZS-K K1=2*KN-3 X1=DTDZS(K1)*THDIF(KN-1) X2=DTDZS(K1+1)*THDIF(KN) FT=TSO(KN)+X1*(TSO(KN-1)-TSO(KN)) & -X2*(TSO(KN)-TSO(KN+1)) DENOM=1.+X1+X2-X2*cotso(K) cotso(K+1)=X1/DENOM rhtso(K+1)=(FT+X2*rhtso(K))/DENOM 33 CONTINUE !--- THE NZS element in COTSO and RHTSO will be for snow !--- There will be 2 layers in snow if it is deeper than DELTSN+SNTH IF(SNHEI.GE.SNTH) then if(snhei.le.DELTSN+SNTH) then !-- 1-layer snow model IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'1-layer - snth,snhei,deltsn',snth,snhei,deltsn ENDIF ilnb=1 snprim=max(snth,snhei) tsob=tso(1) soilt1=tso(1) XSN = DELT/2./(zshalf(2)+0.5*SNPRIM) DDZSN = XSN / SNPRIM X1SN = DDZSN * thdifsn X2 = DTDZS(1)*THDIF(1) FT = TSO(1)+X1SN*(SOILT-TSO(1)) & -X2*(TSO(1)-TSO(2)) DENOM = 1. + X1SN + X2 -X2*cotso(NZS1) cotso(NZS)=X1SN/DENOM rhtso(NZS)=(FT+X2*rhtso(NZS1))/DENOM cotsn=cotso(NZS) rhtsn=rhtso(NZS) !*** Average temperature of snow pack (C) tsnav=0.5*(soilt+tso(1)) & -273.15 else !-- 2 layers in snow, SOILT1 is temperasture at DELTSN depth IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'2-layer - snth,snhei,deltsn',snth,snhei,deltsn ENDIF ilnb=2 snprim=deltsn tsob=soilt1 XSN = DELT/2./(0.5*deltsn) XSN1= DELT/2./(zshalf(2)+0.5*(SNHEI-DELTSN)) DDZSN = XSN / DELTSN DDZSN1 = XSN1 / (SNHEI-DELTSN) X1SN = DDZSN * thdifsn X1SN1 = DDZSN1 * thdifsn X2 = DTDZS(1)*THDIF(1) FT = TSO(1)+X1SN1*(SOILT1-TSO(1)) & -X2*(TSO(1)-TSO(2)) DENOM = 1. + X1SN1 + X2 - X2*cotso(NZS1) cotso(nzs)=x1sn1/denom rhtso(nzs)=(ft+x2*rhtso(nzs1))/denom ftsnow = soilt1+x1sn*(soilt-soilt1) & -x1sn1*(soilt1-tso(1)) denomsn = 1. + X1SN + X1SN1 - X1SN1*cotso(NZS) cotsn=x1sn/denomsn rhtsn=(ftsnow+X1SN1*rhtso(NZS))/denomsn !*** Average temperature of snow pack (C) tsnav=0.5/snhei*((soilt+soilt1)*deltsn & +(soilt1+tso(1))*(SNHEI-DELTSN)) & -273.15 endif ENDIF IF(SNHEI.LT.SNTH.AND.SNHEI.GT.0.) then ! IF(SNHEI.LT.SNTH.AND.SNHEI.GE.0.) then !--- snow is too thin to be treated separately, therefore it !--- is combined with the first soil layer. snprim=SNHEI+zsmain(2) fsn=SNHEI/snprim fso=1.-fsn soilt1=tso(1) tsob=tso(2) XSN = DELT/2./((zshalf(3)-zsmain(2))+0.5*snprim) DDZSN = XSN /snprim X1SN = DDZSN * (fsn*thdifsn+fso*thdif(1)) X2=DTDZS(2)*THDIF(2) FT=TSO(2)+X1SN*(SOILT-TSO(2))- & X2*(TSO(2)-TSO(3)) denom = 1. + x1sn + x2 - x2*cotso(nzs-2) cotso(nzs1) = x1sn/denom rhtso(nzs1)=(FT+X2*rhtso(NZS-2))/denom tsnav=0.5*(soilt+tso(1)) & -273.15 cotso(NZS)=cotso(nzs1) rhtso(NZS)=rhtso(nzs1) cotsn=cotso(NZS) rhtsn=rhtso(NZS) ENDIF !************************************************************************ !--- THE HEAT BALANCE EQUATION (Smirnova et al. 1996, EQ. 21,26) !18apr08 nmelt is the flag for melting, and SNOH is heat of snow phase changes nmelt=0 SNOH=0. ETT1=0. EPOT=-QKMS*(QVATM-QGOLD) RHCS=CAP(1) H=1. TRANS=TRANSUM*DRYCAN/ZSHALF(NROOT+1) CAN=WETCAN+TRANS UMVEG=1.-VEGFRAC FKT=TKMS D1=cotso(NZS1) D2=rhtso(NZS1) TN=SOILT D9=THDIF(1)*RHCS*dzstop D10=TKMS*CP*RHO R211=.5*CONFLX/DELT R21=R211*CP*RHO R22=.5/(THDIF(1)*DELT*dzstop**2) R6=EMISS *STBOLT*.5*TN**4 R7=R6/TN D11=RNET+R6 IF(SNHEI.GE.SNTH) THEN if(snhei.le.DELTSN+SNTH) then !--- 1-layer snow D1SN = cotso(NZS) D2SN = rhtso(NZS) IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'1 layer d1sn,d2sn',i,j,d1sn,d2sn ENDIF else !--- 2-layer snow D1SN = cotsn D2SN = rhtsn IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'2 layers d1sn,d2sn',i,j,d1sn,d2sn ENDIF endif D9SN= THDIFSN*RHOCSN / SNPRIM R22SN = SNPRIM*SNPRIM*0.5/(THDIFSN*DELT) IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'1 or 2 layers D9sn,R22sn',d9sn,r22sn ENDIF ENDIF IF(SNHEI.LT.SNTH.AND.SNHEI.GT.0.) then !--- thin snow is combined with soil D1SN = D1 D2SN = D2 D9SN = (fsn*THDIFSN*RHOCSN+fso*THDIF(1)*RHCS)/ & snprim R22SN = snprim*snprim*0.5 & /((fsn*THDIFSN+fso*THDIF(1))*delt) IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,' Combined D9SN,R22SN,D1SN,D2SN: ',D9SN,R22SN,D1SN,D2SN ENDIF ENDIF IF(SNHEI.eq.0.)then !--- all snow is sublimated D9SN = D9 R22SN = R22 D1SN = D1 D2SN = D2 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,' SNHEI = 0, D9SN,R22SN,D1SN,D2SN: ',D9SN,R22SN,D1SN,D2SN ENDIF ENDIF 2211 continue !18apr08 - the snow melt iteration start point 212 continue !---- TDENOM for snow TDENOM = D9SN*(1.-D1SN +R22SN)+D10+R21+R7 & +RAINF*CVW*PRCPMS & +RHOnewCSN*NEWSNOW/DELT FKQ=QKMS*RHO R210=R211*RHO C=VEGFRAC*FKQ*CAN CC=C*XLVM/TDENOM AA=XLVM*(BETA*FKQ*UMVEG+R210)/TDENOM BB=(D10*TABS+R21*TN+XLVM*(QVATM* & (BETA*FKQ*UMVEG+C) & +R210*QGOLD)+D11+D9SN*(D2SN+R22SN*TN) & +RAINF*CVW*PRCPMS*max(273.15,TABS) & + RHOnewCSN*NEWSNOW/DELT*min(273.15,TABS) & )/TDENOM AA1=AA+CC PP=PATM*1.E3 AA1=AA1/PP BB=BB-SNOH/TDENOM CALL VILKA(TN,AA1,BB,PP,QS1,TS1,TBQ,KTAU,i,j,iland,isoil) TQ2=QVATM TX2=TQ2*(1.-H) Q1=TX2+H*QS1 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'VILKA1 - TS1,QS1,TQ2,H,TX2,Q1',TS1,QS1,TQ2,H,TX2,Q1 ENDIF IF(Q1.LT.QS1) GOTO 100 !--- if no saturation - goto 100 !--- if saturation - goto 90 90 QVG=QS1 QSG=QS1 QCG=max(0.,Q1-QS1) IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'90 QVG,QSG,QCG,TSO(1)',QVG,QSG,QCG,TSO(1) ENDIF GOTO 200 100 BB=BB-AA*TX2 AA=(AA*H+CC)/PP CALL VILKA(TN,AA,BB,PP,QS1,TS1,TBQ,KTAU,i,j,iland,isoil) Q1=TX2+H*QS1 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'VILKA2 - TS1,QS1,H,TX2,Q1',TS1,QS1,TQ2,H,TX2,Q1 ENDIF IF(Q1.GT.QS1) GOTO 90 QSG=QS1 QVG=Q1 QCG=0. IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'No Saturation QVG,QSG,QCG,TSO(1)',QVG,QSG,QCG,TSO(1) ENDIF 200 CONTINUE if(1==2) then !--11Aug18 - do not use this iteration to be consistent with SOILTEMP. if(qvatm > QSG .and. iter==0) then !condensation regime IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'SNOW turn off canopy evaporation and transpiration' print *,' QVATM,QVG,QSG,TS1',QVATM,QVG,QSG,TS1 ENDIF can=0. umveg=1. iter=1 goto 2211 endif IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN if(iter==1) then print *,'SNOW - QVATM,QVG,QSG,QCG,TS1',QVATM,QVG,QSG,QCG,TS1 endif ENDIF endif ! 1==2 !--- SOILT - skin temperature SOILT=TS1 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN IF(i.eq.266.and.j.eq.447) then print *,'snwe,snhei,soilt,soilt1,tso',i,j,snwe,snhei,soilt,soilt1,tso endif ENDIF ! Solution for temperature at 7.5 cm depth and snow-soil interface IF(SNHEI.GE.SNTH) THEN if(snhei.gt.DELTSN+SNTH) then !-- 2-layer snow model SOILT1=min(273.15,rhtsn+cotsn*SOILT) TSO(1)=rhtso(NZS)+cotso(NZS)*SOILT1 tsob=soilt1 else !-- 1 layer in snow TSO(1)=rhtso(NZS)+cotso(NZS)*SOILT SOILT1=TSO(1) tsob=tso(1) endif ELSEIF (SNHEI > 0. .and. SNHEI < SNTH) THEN ! blended TSO(2)=rhtso(NZS1)+cotso(NZS1)*SOILT tso(1)=(tso(2)+(soilt-tso(2))*fso) SOILT1=TSO(1) tsob=TSO(2) ELSE !-- very thin or zero snow. If snow is thin we suppose that !--- tso(i,j,1)=SOILT, and later we recompute tso(i,j,1) TSO(1)=SOILT SOILT1=SOILT tsob=TSO(1) !new tsob=tso(2) ENDIF !---- Final solution for TSO IF (SNHEI > 0. .and. SNHEI < SNTH) THEN ! blended or snow is melted DO K=3,NZS KK=NZS-K+1 TSO(K)=rhtso(KK)+cotso(KK)*TSO(K-1) END DO ELSE DO K=2,NZS KK=NZS-K+1 TSO(K)=rhtso(KK)+cotso(KK)*TSO(K-1) END DO ENDIF !--- For thin snow layer combined with the top soil layer !--- TSO(1) is recomputed by linear interpolation between SOILT !--- and TSO(i,j,2) ! if(SNHEI.LT.SNTH.AND.SNHEI.GT.0.)then ! tso(1)=tso(2)+(soilt-tso(2))*fso ! soilt1=tso(1) ! tsob = tso(2) ! endif IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN ! IF(i.eq.266.and.j.eq.447) then print *,'SOILT,SOILT1,tso,TSOB,QSG',i,j,SOILT,SOILT1,tso,TSOB,QSG,'nmelt=',nmelt ENDIF if(nmelt.eq.1) go to 220 !--- IF SOILT > 273.15 F then melting of snow can happen ! IF(SOILT.GT.273.15.AND.SNHEI.GT.0.) THEN ! if all snow can evaporate, then there is nothing to melt IF(SOILT.GT.273.15.AND.SNWEPR-BETA*EPOT*RAS*DELT.GT.0.AND.SNHEI.GT.0.) THEN nmelt = 1 soiltfrac=snowfrac*273.15+(1.-snowfrac)*SOILT QSG=min(QSG, QSN(soiltfrac,TBQ)/PP) qvg=qsg T3 = STBOLT*TN*TN*TN UPFLUX = T3 * 0.5*(TN + SOILTfrac) XINET = EMISS*(GLW-UPFLUX) ! RNET = GSW + XINET EPOT = -QKMS*(QVATM-QSG) Q1=EPOT*RAS IF (Q1.LE.0..or.iter==1) THEN ! --- condensation DEW=-EPOT DO K=1,NZS TRANSP(K)=0. ENDDO QFX = -XLVM*RHO*DEW EETA = QFX/XLVM ELSE ! --- evaporation DO K=1,NROOT TRANSP(K)=-VEGFRAC*q1 & *TRANF(K)*DRYCAN/zshalf(NROOT+1) ! IF(TRANSP(K).GT.0.) TRANSP(K)=0. ETT1=ETT1-TRANSP(K) ENDDO DO k=nroot+1,nzs transp(k)=0. enddo EDIR1 = Q1*UMVEG * BETA EC1 = Q1 * WETCAN * vegfrac CMC2MS=CST/DELT*RAS ! EC1=MIN(CMC2MS,EC1) EETA = (EDIR1 + EC1 + ETT1)*1.E3 ! to convert from kg m-2 s-1 to m s-1: 1/rho water=1.e-3************ QFX= XLVM * EETA ENDIF HFX=-D10*(TABS-soiltfrac) IF(SNHEI.GE.SNTH)then SOH=thdifsn*RHOCSN*(soiltfrac-TSOB)/SNPRIM SNFLX=SOH ELSE SOH=(fsn*thdifsn*rhocsn+fso*thdif(1)*rhcs)* & (soiltfrac-TSOB)/snprim SNFLX=SOH ENDIF ! X= (R21+D9SN*R22SN)*(soiltfrac-TN) + & XLVM*R210*(QVG-QGOLD) IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'SNOWTEMP storage ',i,j,x print *,'R21,D9sn,r22sn,soiltfrac,tn,qsg,qvg,qgold,snprim', & R21,D9sn,r22sn,soiltfrac,tn,qsg,qvg,qgold,snprim ENDIF !-- SNOH is energy flux of snow phase change SNOH=RNET-QFX -HFX - SOH - X & +RHOnewCSN*NEWSNOW/DELT*(min(273.15,TABS)-soiltfrac) & +RAINF*CVW*PRCPMS*(max(273.15,TABS)-soiltfrac) SNOH=AMAX1(0.,SNOH) !-- SMELT is speed of melting in M/S SMELT= SNOH /XLMELT*1.E-3 IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'1- SMELT',i,j,smelt ENDIF SMELT=AMIN1(SMELT,SNWEPR/DELT-BETA*EPOT*RAS) IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'2- SMELT',i,j,smelt ENDIF SMELT=AMAX1(0.,SMELT) !18apr08 - Egglston limit ! SMELT= amin1 (smelt, 5.6E-7*meltfactor*max(1.,(soilt-273.15))) SMELT= amin1 (smelt, 5.6E-8*meltfactor*max(1.,(soilt-273.15))) IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'3- SMELT',i,j,smelt ENDIF ! rr - potential melting rr=max(0.,SNWEPR/delt-BETA*EPOT*RAS) SMELT=min(SMELT,rr) IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'4- SMELT i,j,smelt,rr',i,j,smelt,rr ENDIF SNOHGNEW=SMELT*XLMELT*1.E3 SNODIF=AMAX1(0.,(SNOH-SNOHGNEW)) SNOH=SNOHGNEW IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'SNOH,SNODIF',SNOH,SNODIF ENDIF !*** From Koren et al. (1999) 13% of snow melt stays in the snow pack rsmfrac=min(0.18,(max(0.08,snwepr/0.10*0.13))) if(snhei > 0.01) then rsm=rsmfrac*smelt*delt else ! do not keep melted water if snow depth is less that 1 cm rsm=0. endif !18apr08 rsm is part of melted water that stays in snow as liquid SMELT=max(0.,SMELT-rsm/delt) IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'5- SMELT i,j,smelt,rsm,snwepr,rsmfrac', & i,j,smelt,rsm,snwepr,rsmfrac ENDIF !-- update of liquid equivalent of snow depth !-- due to evaporation and snow melt SNWE = AMAX1(0.,(SNWEPR- & (SMELT+BETA*EPOT*RAS)*DELT & ! (SMELT+BETA*EPOT*RAS)*DELT*snowfrac & ! (SMELT+BETA*EPOT*RAS*UMVEG)*DELT & ) ) !--- If there is no snow melting then just evaporation !--- or condensation cxhanges SNWE ELSE if(snhei.ne.0.) then EPOT=-QKMS*(QVATM-QSG) SNWE = AMAX1(0.,(SNWEPR- & BETA*EPOT*RAS*DELT)) ! BETA*EPOT*RAS*DELT*snowfrac)) endif ENDIF !18apr08 - if snow melt occurred then go into iteration for energy budget ! solution if(nmelt.eq.1) goto 212 ! second interation 220 continue if(smelt.gt.0..and.rsm.gt.0.) then if(snwe.le.rsm) then IF ( 1==1 ) THEN print *,'SNWE 0.) THEN if (snhei.GT.deltsn+snth) then hsn = snhei - deltsn IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print*,'2 layer snow - snhei,hsn',snhei,hsn ENDIF else IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print*,'1 layer snow or blended - snhei',snhei ENDIF hsn = snhei endif soiltfrac=snowfrac*273.15+(1.-snowfrac)*TSO(1) SNOHG=(TSO(1)-soiltfrac)*(cap(1)*zshalf(2)+ & RHOCSN*0.5*hsn) / DELT SNOHG=AMAX1(0.,SNOHG) SNODIF=0. SMELTG=SNOHG/XLMELT*1.E-3 ! Egglston - empirical limit on snow melt from the bottom of snow pack SMELTG=AMIN1(SMELTG, 5.8e-9) ! rr - potential melting rr=SNWE/delt SMELTG=AMIN1(SMELTG, rr) SNOHGNEW=SMELTG*XLMELT*1.e3 SNODIF=AMAX1(0.,(SNOHG-SNOHGNEW)) IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN ! if(i.eq.266.and.j.eq.447) then print *,'TSO(1),soiltfrac,smeltg,SNODIF',TSO(1),soiltfrac,smeltg,SNODIF ENDIF ! snwe=max(0.,snwe-smeltg*delt*snowfrac) snwe=max(0.,snwe-smeltg*delt) SNHEI=SNWE *1.E3 / RHOSN if(snhei > 0.) TSO(1) = soiltfrac IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN ! if(i.eq.266.and.j.eq.447) then print *,'Melt from the bottom snwe,snhei',snwe,snhei if (snhei==0.) & print *,'Snow is all melted on the warm ground' ENDIF ENDIF IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'SNHEI,SNOH',i,j,SNHEI,SNOH ENDIF ! & snweprint=snwe snheiprint=snweprint*1.E3 / RHOSN IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *, 'snweprint : ',snweprint print *, 'D9SN,SOILT,TSOB : ', D9SN,SOILT,TSOB ENDIF X= (R21+D9SN*R22SN)*(soilt-TN) + & XLVM*R210*(QSG-QGOLD) IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'SNOWTEMP storage ',i,j,x print *,'R21,D9sn,r22sn,soiltfrac,soilt,tn,qsg,qgold,snprim', & R21,D9sn,r22sn,soiltfrac,soilt,tn,qsg,qgold,snprim ENDIF X=X & ! "heat" from snow and rain -RHOnewCSN*NEWSNOW/DELT*(min(273.15,TABS)-SOILT) & -RAINF*CVW*PRCPMS*(max(273.15,TABS)-SOILT) IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'x=',x print *,'SNHEI=',snhei print *,'SNFLX=',snflx ENDIF IF(SNHEI.GT.0.) THEN if(ilnb.gt.1) then tsnav=0.5/snhei*((soilt+soilt1)*deltsn & +(soilt1+tso(1))*(SNHEI-DELTSN)) & -273.15 else tsnav=0.5*(soilt+tso(1)) - 273.15 endif ELSE tsnav= soilt - 273.15 ENDIF !------------------------------------------------------------------------ END SUBROUTINE SNOWTEMP !------------------------------------------------------------------------ SUBROUTINE SOILMOIST ( & !--input parameters DELT,NZS,NDDZS,DTDZS,DTDZS2,RIW, & ZSMAIN,ZSHALF,DIFFU,HYDRO, & QSG,QVG,QCG,QCATM,QVATM,PRCP, & QKMS,TRANSP,DRIP, & DEW,SMELT,SOILICE,VEGFRAC,SNOWFRAC,soilres, & !--soil properties DQM,QMIN,REF,KSAT,RAS,INFMAX, & !--output SOILMOIS,SOILIQW,MAVAIL,RUNOFF,RUNOFF2,INFILTRP) !************************************************************************* ! moisture balance equation and Richards eqn. ! are solved here ! ! DELT - time step (s) ! IME,JME,NZS - dimensions of soil domain ! ZSMAIN - main levels in soil (m) ! ZSHALF - middle of the soil layers (m) ! DTDZS - dt/(2.*dzshalf*dzmain) ! DTDZS2 - dt/(2.*dzshalf) ! DIFFU - diffusional conductivity (m^2/s) ! HYDRO - hydraulic conductivity (m/s) ! QSG,QVG,QCG - saturated mixing ratio, mixing ratio of ! water vapor and cloud at the ground ! surface, respectively (kg/kg) ! QCATM,QVATM - cloud and water vapor mixing ratio ! at the first atm. level (kg/kg) ! PRCP - precipitation rate in m/s ! QKMS - exchange coefficient for water vapor in the ! surface layer (m/s) ! TRANSP - transpiration from the soil layers (m/s) ! DRIP - liquid water dripping from the canopy to soil (m) ! DEW - dew in kg/m^2s ! SMELT - melting rate in m/s ! SOILICE - volumetric content of ice in soil (m^3/m^3) ! SOILIQW - volumetric content of liquid water in soil (m^3/m^3) ! VEGFRAC - greeness fraction (0-1) ! RAS - ration of air density to soil density ! INFMAX - maximum infiltration rate (kg/m^2/s) ! ! SOILMOIS - volumetric soil moisture, 6 levels (m^3/m^3) ! MAVAIL - fraction of maximum soil moisture in the top ! layer (0-1) ! RUNOFF - surface runoff (m/s) ! RUNOFF2 - underground runoff (m) ! INFILTRP - point infiltration flux into soil (m/s) ! /(snow bottom runoff) (mm/s) ! ! COSMC, RHSMC - coefficients for implicit solution of ! Richards equation !****************************************************************** IMPLICIT NONE !------------------------------------------------------------------ !--- input variables REAL, INTENT(IN ) :: DELT INTEGER, INTENT(IN ) :: NZS,NDDZS ! input variables REAL, DIMENSION(1:NZS), INTENT(IN ) :: ZSMAIN, & ZSHALF, & DIFFU, & HYDRO, & TRANSP, & SOILICE, & DTDZS2 REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS REAL, INTENT(IN ) :: QSG,QVG,QCG,QCATM,QVATM , & QKMS,VEGFRAC,DRIP,PRCP , & DEW,SMELT,SNOWFRAC , & DQM,QMIN,REF,KSAT,RAS,RIW,SOILRES ! output REAL, DIMENSION( 1:nzs ) , & INTENT(INOUT) :: SOILMOIS,SOILIQW REAL, INTENT(INOUT) :: MAVAIL,RUNOFF,RUNOFF2,INFILTRP, & INFMAX ! local variables REAL, DIMENSION( 1:nzs ) :: COSMC,RHSMC REAL :: DZS,R1,R2,R3,R4,R5,R6,R7,R8,R9,R10 REAL :: REFKDT,REFDK,DELT1,F1MAX,F2MAX REAL :: F1,F2,FD,KDT,VAL,DDT,PX,FK,FKMAX REAL :: QQ,UMVEG,INFMAX1,TRANS REAL :: TOTLIQ,FLX,FLXSAT,QTOT REAL :: DID,X1,X2,X4,DENOM,Q2,Q4 REAL :: dice,fcr,acrt,frzx,sum,cvfrz INTEGER :: NZS1,NZS2,K,KK,K1,KN,ialp1,jj,jk !****************************************************************************** ! COEFFICIENTS FOR THOMAS ALGORITHM FOR SOILMOIS !****************************************************************************** NZS1=NZS-1 NZS2=NZS-2 118 format(6(10Pf23.19)) do k=1,nzs cosmc(k)=0. rhsmc(k)=0. enddo DID=(ZSMAIN(NZS)-ZSHALF(NZS)) X1=ZSMAIN(NZS)-ZSMAIN(NZS1) !7may09 DID=(ZSMAIN(NZS)-ZSHALF(NZS))*2. ! DENOM=DID/DELT+DIFFU(NZS1)/X1 ! COSMC(1)=DIFFU(NZS1)/X1/DENOM ! RHSMC(1)=(SOILMOIS(NZS)*DID/DELT ! 1 +TRANSP(NZS)-(HYDRO(NZS)*SOILMOIS(NZS) ! 1 -HYDRO(NZS1)*SOILMOIS(NZS1))*DID ! 1 /X1) /DENOM DENOM=(1.+DIFFU(nzs1)/X1/DID*DELT+HYDRO(NZS)/(2.*DID)*DELT) COSMC(1)=DELT*(DIFFU(nzs1)/DID/X1 & +HYDRO(NZS1)/2./DID)/DENOM RHSMC(1)=(SOILMOIS(NZS)+TRANSP(NZS)*DELT/ & DID)/DENOM ! RHSMC(1)=(SOILMOIS(NZS)*DID/DELT & ! +TRANSP(NZS)-(HYDRO(NZS)*SOILMOIS(NZS) & ! -HYDRO(NZS1)*SOILMOIS(NZS1))*DID & ! /X1) /DENOM !12 June 2014 - low boundary condition: 1 - zero diffusion below the lowest ! level; 2 - soil moisture at the low boundary can be lost due to the root uptake. ! So far - no interaction with the water table. DENOM=1.+DIFFU(nzs1)/X1/DID*DELT !orig DENOM=(1.+DIFFU(nzs1)/X1/DID*DELT+HYDRO(NZS)/DID*DELT) !orig COSMC(1)=DELT*(DIFFU(nzs1)/DID/X1 & !orig +HYDRO(NZS1)/2./DID)/DENOM COSMC(1)=DELT*(DIFFU(nzs1)/DID/X1 & +HYDRO(NZS1)/DID)/DENOM ! RHSMC(1)=(SOILMOIS(NZS)+TRANSP(NZS)*DELT/ & ! DID)/DENOM RHSMC(1)=(SOILMOIS(NZS)-HYDRO(NZS)*DELT/DID*soilmois(nzs) & +TRANSP(NZS)*DELT/DID)/DENOM !test RHSMC(1)=SOILMOIS(NZS)-HYDRO(NZS)*soilmois(nzs) !test!!! !this test gave smoother soil moisture, ovwerall better results COSMC(1)=0. RHSMC(1)=SOILMOIS(NZS) ! DO 330 K=1,NZS2 KN=NZS-K K1=2*KN-3 X4=2.*DTDZS(K1)*DIFFU(KN-1) X2=2.*DTDZS(K1+1)*DIFFU(KN) Q4=X4+HYDRO(KN-1)*DTDZS2(KN-1) Q2=X2-HYDRO(KN+1)*DTDZS2(KN-1) DENOM=1.+X2+X4-Q2*COSMC(K) COSMC(K+1)=Q4/DENOM IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'q2,soilmois(kn),DIFFU(KN),x2,HYDRO(KN+1),DTDZS2(KN-1),kn,k' & ,q2,soilmois(kn),DIFFU(KN),x2,HYDRO(KN+1),DTDZS2(KN-1),kn,k ENDIF 330 RHSMC(K+1)=(SOILMOIS(KN)+Q2*RHSMC(K) & +TRANSP(KN) & /(ZSHALF(KN+1)-ZSHALF(KN)) & *DELT)/DENOM ! --- MOISTURE BALANCE BEGINS HERE TRANS=TRANSP(1) UMVEG=(1.-VEGFRAC)*soilres RUNOFF=0. RUNOFF2=0. DZS=ZSMAIN(2) R1=COSMC(NZS1) R2= RHSMC(NZS1) R3=DIFFU(1)/DZS R4=R3+HYDRO(1)*.5 R5=R3-HYDRO(2)*.5 R6=QKMS*RAS !-- Total liquid water available on the top of soil domain !-- Without snow - 3 sources of water: precipitation, !-- water dripping from the canopy and dew !-- With snow - only one source of water - snow melt 191 format (f23.19) ! TOTLIQ=UMVEG*PRCP-DRIP/DELT-UMVEG*DEW*RAS-SMELT TOTLIQ=PRCP-DRIP/DELT-UMVEG*DEW*RAS-SMELT IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'UMVEG*PRCP,DRIP/DELT,UMVEG*DEW*RAS,SMELT', & UMVEG*PRCP,DRIP/DELT,UMVEG*DEW*RAS,SMELT ENDIF !test 16 may TOTLIQ=UMVEG*PRCP-DRIP/DELT-UMVEG*DEW*RAS-SMELT !30july13 TOTLIQ=UMVEG*PRCP-DRIP/DELT-SMELT FLX=TOTLIQ INFILTRP=TOTLIQ ! ----------- FROZEN GROUND VERSION ------------------------- ! REFERENCE FROZEN GROUND PARAMETER, CVFRZ, IS A SHAPE PARAMETER OF ! AREAL DISTRIBUTION FUNCTION OF SOIL ICE CONTENT WHICH EQUALS 1/CV. ! CV IS A COEFFICIENT OF SPATIAL VARIATION OF SOIL ICE CONTENT. ! BASED ON FIELD DATA CV DEPENDS ON AREAL MEAN OF FROZEN DEPTH, AND IT ! CLOSE TO CONSTANT = 0.6 IF AREAL MEAN FROZEN DEPTH IS ABOVE 20 CM. ! THAT IS WHY PARAMETER CVFRZ = 3 (INT{1/0.6*0.6}) ! ! Current logic doesn't allow CVFRZ be bigger than 3 CVFRZ = 3. !-- SCHAAKE/KOREN EXPRESSION for calculation of max infiltration REFKDT=3. REFDK=3.4341E-6 DELT1=DELT/86400. F1MAX=DQM*ZSHALF(2) F2MAX=DQM*(ZSHALF(3)-ZSHALF(2)) F1=F1MAX*(1.-SOILMOIS(1)/DQM) DICE=SOILICE(1)*ZSHALF(2) FD=F1 do k=2,nzs1 DICE=DICE+(ZSHALF(k+1)-ZSHALF(k))*SOILICE(K) FKMAX=DQM*(ZSHALF(k+1)-ZSHALF(k)) FK=FKMAX*(1.-SOILMOIS(k)/DQM) FD=FD+FK enddo KDT=REFKDT*KSAT/REFDK VAL=(1.-EXP(-KDT*DELT1)) DDT = FD*VAL PX= - TOTLIQ * DELT IF(PX.LT.0.0) PX = 0.0 IF(PX.gt.0.0) THEN INFMAX1 = (PX*(DDT/(PX+DDT)))/DELT ELSE INFMAX1 = 0. ENDIF IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'INFMAX1 before frozen part',INFMAX1 ENDIF ! ----------- FROZEN GROUND VERSION -------------------------- ! REDUCTION OF INFILTRATION BASED ON FROZEN GROUND PARAMETERS ! ! ------------------------------------------------------------------ FRZX= 0.15*((dqm+qmin)/ref) * (0.412 / 0.468) FCR = 1. IF ( DICE .GT. 1.E-2) THEN ACRT = CVFRZ * FRZX / DICE SUM = 1. IALP1 = CVFRZ - 1 DO JK = 1,IALP1 K = 1 DO JJ = JK+1, IALP1 K = K * JJ END DO SUM = SUM + (ACRT ** ( CVFRZ-JK)) / FLOAT (K) END DO FCR = 1. - EXP(-ACRT) * SUM END IF IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'FCR--------',fcr print *,'DICE=',dice ENDIF INFMAX1 = INFMAX1* FCR ! ------------------------------------------------------------------- INFMAX = MAX(INFMAX1,HYDRO(1)*SOILMOIS(1)) INFMAX = MIN(INFMAX, -TOTLIQ) IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'INFMAX,INFMAX1,HYDRO(1)*SOILIQW(1),-TOTLIQ', & INFMAX,INFMAX1,HYDRO(1)*SOILIQW(1),-TOTLIQ ENDIF !---- IF (-TOTLIQ.GT.INFMAX)THEN RUNOFF=-TOTLIQ-INFMAX FLX=-INFMAX IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'FLX,RUNOFF1=',flx,runoff ENDIF ENDIF ! INFILTRP is total infiltration flux in M/S INFILTRP=FLX ! Solution of moisture budget R7=.5*DZS/DELT R4=R4+R7 FLX=FLX-SOILMOIS(1)*R7 ! R8 is for direct evaporation from soil, which occurs ! only from snow-free areas ! R8=UMVEG*R6 R8=UMVEG*R6*(1.-snowfrac) QTOT=QVATM+QCATM R9=TRANS R10=QTOT-QSG !-- evaporation regime IF(R10.LE.0.) THEN QQ=(R5*R2-FLX+R9)/(R4-R5*R1-R10*R8/(REF-QMIN)) FLXSAT=-DQM*(R4-R5*R1-R10*R8/(REF-QMIN)) & +R5*R2+R9 ELSE !-- dew formation regime QQ=(R2*R5-FLX+R8*(QTOT-QCG-QVG)+R9)/(R4-R1*R5) FLXSAT=-DQM*(R4-R1*R5)+R2*R5+R8*(QTOT-QVG-QCG)+R9 END IF IF(QQ.LT.0.) THEN ! print *,'negative QQ=',qq SOILMOIS(1)=1.e-8 ELSE IF(QQ.GT.DQM) THEN !-- saturation SOILMOIS(1)=DQM IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'FLXSAT,FLX,DELT',FLXSAT,FLX,DELT,RUNOFF2 ENDIF ! RUNOFF2=(FLXSAT-FLX) RUNOFF=RUNOFF+(FLXSAT-FLX) ELSE SOILMOIS(1)=min(dqm,max(1.e-8,QQ)) END IF IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'SOILMOIS,SOILIQW, soilice',SOILMOIS,SOILIQW,soilice*riw print *,'COSMC,RHSMC',COSMC,RHSMC ENDIF !--- FINAL SOLUTION FOR SOILMOIS ! DO K=2,NZS1 DO K=2,NZS KK=NZS-K+1 QQ=COSMC(KK)*SOILMOIS(K-1)+RHSMC(KK) ! QQ=COSMC(KK)*SOILIQW(K-1)+RHSMC(KK) IF (QQ.LT.0.) THEN ! print *,'negative QQ=',qq SOILMOIS(K)=1.e-8 ELSE IF(QQ.GT.DQM) THEN !-- saturation SOILMOIS(K)=DQM IF(K.EQ.NZS)THEN IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'hydro(k),QQ,DQM,k',hydro(k),QQ,DQM,k ENDIF RUNOFF2=RUNOFF2+((QQ-DQM)*(ZSMAIN(K)-ZSHALF(K)))/DELT ! RUNOFF2=RUNOFF2+(QQ-DQM)*hydro(k) ! print *,'RUNOFF2=',RUNOFF2 ELSE ! print *,'QQ,DQM,k',QQ,DQM,k RUNOFF2=RUNOFF2+((QQ-DQM)*(ZSHALF(K+1)-ZSHALF(K)))/DELT ! RUNOFF2=RUNOFF2+(QQ-DQM)*hydro(k) ENDIF ELSE SOILMOIS(K)=min(dqm,max(1.e-8,QQ)) END IF END DO IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'END soilmois,soiliqw,soilice',soilmois,SOILIQW,soilice*riw ENDIF ! RUNOFF2=RUNOFF2+hydro(nzs)*SOILMOIS(NZS) ! MAVAIL=max(.00001,min(1.,SOILMOIS(1)/DQM)) ! MAVAIL=max(.00001,min(1.,SOILMOIS(1)/(REF-QMIN))) MAVAIL=max(.00001,min(1.,(SOILMOIS(1)/(REF-QMIN)*(1.-snowfrac)+1.*snowfrac))) ! RETURN ! END !------------------------------------------------------------------- END SUBROUTINE SOILMOIST !------------------------------------------------------------------- SUBROUTINE SOILPROP(spp_lsm,rstochcol,fieldcol_sf, & !--- input variables nzs,fwsat,lwsat,tav,keepfr, & soilmois,soiliqw,soilice, & soilmoism,soiliqwm,soilicem, & !--- soil fixed fields QWRTZ,rhocs,dqm,qmin,psis,bclh,ksat, & !--- constants riw,xlmelt,CP,G0_P,cvw,ci, & kqwrtz,kice,kwt, & !--- output variables thdif,diffu,hydro,cap) !****************************************************************** ! SOILPROP computes thermal diffusivity, and diffusional and ! hydraulic condeuctivities !****************************************************************** ! NX,NY,NZS - dimensions of soil domain ! FWSAT, LWSAT - volumetric content of frozen and liquid water ! for saturated condition at given temperatures (m^3/m^3) ! TAV - temperature averaged for soil layers (K) ! SOILMOIS - volumetric soil moisture at the main soil levels (m^3/m^3) ! SOILMOISM - volumetric soil moisture averaged for layers (m^3/m^3) ! SOILIQWM - volumetric liquid soil moisture averaged for layers (m^3/m^3) ! SOILICEM - volumetric content of soil ice averaged for layers (m^3/m^3) ! THDIF - thermal diffusivity for soil layers (W/m/K) ! DIFFU - diffusional conductivity (m^2/s) ! HYDRO - hydraulic conductivity (m/s) ! CAP - volumetric heat capacity (J/m^3/K) ! !****************************************************************** IMPLICIT NONE !----------------------------------------------------------------- !--- soil properties INTEGER, INTENT(IN ) :: NZS REAL , & INTENT(IN ) :: RHOCS, & BCLH, & DQM, & KSAT, & PSIS, & QWRTZ, & QMIN REAL, DIMENSION( 1:nzs ) , & INTENT(IN ) :: SOILMOIS, & keepfr REAL, INTENT(IN ) :: CP, & CVW, & RIW, & kqwrtz, & kice, & kwt, & XLMELT, & G0_P REAL, DIMENSION(1:NZS), INTENT(IN) :: rstochcol REAL, DIMENSION(1:NZS), INTENT(INOUT) :: fieldcol_sf INTEGER, INTENT(IN ) :: spp_lsm !--- output variables REAL, DIMENSION(1:NZS) , & INTENT(INOUT) :: cap,diffu,hydro , & thdif,tav , & soilmoism , & soiliqw,soilice , & soilicem,soiliqwm , & fwsat,lwsat !--- local variables REAL, DIMENSION(1:NZS) :: hk,detal,kasat,kjpl REAL :: x,x1,x2,x4,ws,wd,fact,fach,facd,psif,ci REAL :: tln,tavln,tn,pf,a,am,ame,h INTEGER :: nzs1,k !-- for Johansen thermal conductivity REAL :: kzero,gamd,kdry,kas,x5,sr,ke nzs1=nzs-1 !-- Constants for Johansen (1975) thermal conductivity kzero =2. ! if qwrtz > 0.2 do k=1,nzs detal (k)=0. kasat (k)=0. kjpl (k)=0. hk (k)=0. enddo ws=dqm+qmin x1=xlmelt/(g0_p*psis) x2=x1/bclh*ws x4=(bclh+1.)/bclh !--- Next 3 lines are for Johansen thermal conduct. gamd=(1.-ws)*2700. kdry=(0.135*gamd+64.7)/(2700.-0.947*gamd) kas=kqwrtz**qwrtz*kzero**(1.-qwrtz) DO K=1,NZS1 tn=tav(k) - 273.15 wd=ws - riw*soilicem(k) psif=psis*100.*(wd/(soiliqwm(k)+qmin))**bclh & * (ws/wd)**3. !--- PSIF should be in [CM] to compute PF pf=log10(abs(psif)) fact=1.+riw*soilicem(k) !--- HK is for McCumber thermal conductivity IF(PF.LE.5.2) THEN HK(K)=420.*EXP(-(PF+2.7))*fact ELSE HK(K)=.1744*fact END IF IF(soilicem(k).NE.0.AND.TN.LT.0.) then !--- DETAL is taking care of energy spent on freezing or released from ! melting of soil water DETAL(K)=273.15*X2/(TAV(K)*TAV(K))* & (TAV(K)/(X1*TN))**X4 if(keepfr(k).eq.1.) then detal(k)=0. endif ENDIF !--- Next 10 lines calculate Johansen thermal conductivity KJPL kasat(k)=kas**(1.-ws)*kice**fwsat(k) & *kwt**lwsat(k) X5=(soilmoism(k)+qmin)/ws if(soilicem(k).eq.0.) then sr=max(0.101,x5) ke=log10(sr)+1. !--- next 2 lines - for coarse soils ! sr=max(0.0501,x5) ! ke=0.7*log10(sr)+1. else ke=x5 endif kjpl(k)=ke*(kasat(k)-kdry)+kdry !--- CAP -volumetric heat capacity CAP(K)=(1.-WS)*RHOCS & + (soiliqwm(K)+qmin)*CVW & + soilicem(K)*CI & + (dqm-soilmoism(k))*CP*1.2 & - DETAL(K)*1.e3*xlmelt a=RIW*soilicem(K) if((ws-a).lt.0.12)then diffu(K)=0. else H=max(0.,(soilmoism(K)-a)/(max(1.e-8,(dqm-a)))) facd=1. if(a.ne.0.)facd=1.-a/max(1.e-8,soilmoism(K)) ame=max(1.e-8,dqm-riw*soilicem(K)) !--- DIFFU is diffusional conductivity of soil water diffu(K)=-BCLH*KSAT*PSIS/ame* & (dqm/ame)**3. & *H**(BCLH+2.)*facd endif ! diffu(K)=-BCLH*KSAT*PSIS/dqm & ! *H**(BCLH+2.) !--- thdif - thermal diffusivity ! thdif(K)=HK(K)/CAP(K) !--- Use thermal conductivity from Johansen (1975) thdif(K)=KJPL(K)/CAP(K) END DO IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,'soilice*riw,soiliqw,soilmois,ws',soilice*riw,soiliqw,soilmois,ws ENDIF DO K=1,NZS if((ws-riw*soilice(k)).lt.0.12)then hydro(k)=0. else fach=1. if(soilice(k).ne.0.) & fach=1.-riw*soilice(k)/max(1.e-8,soilmois(k)) am=max(1.e-8,dqm-riw*soilice(k)) !--- HYDRO is hydraulic conductivity of soil water hydro(K)=min(KSAT,KSAT/am* & (soiliqw(K)/am) & **(2.*BCLH+2.) & * fach) if(hydro(k)<1.e-10)hydro(k)=0. endif ENDDO !beka perturb hydrolic conductivity by 10-30%, not more than 50% ! if (spp_lsm==1) then ! DO K=1,NZS !lala ! fieldcol_sf(k)=hydro(k)*rstochcol(k) ! hydro(k)=hydro(k)*(1+rstochcol(k)) ! ENDDO ! ENDIF ! RETURN ! END !----------------------------------------------------------------------- END SUBROUTINE SOILPROP !----------------------------------------------------------------------- SUBROUTINE TRANSF(i,j, & !--- input variables nzs,nroot,soiliqw,tabs,lai,gswin, & !--- soil fixed fields dqm,qmin,ref,wilt,zshalf,pc,iland, & !--- output variables tranf,transum) !------------------------------------------------------------------- !--- TRANF(K) - THE TRANSPIRATION FUNCTION (Smirnova et al. 1996, EQ. 18,19) !******************************************************************* ! NX,NY,NZS - dimensions of soil domain ! SOILIQW - volumetric liquid soil moisture at the main levels (m^3/m^3) ! TRANF - the transpiration function at levels (m) ! TRANSUM - transpiration function integrated over the rooting zone (m) ! !******************************************************************* IMPLICIT NONE !------------------------------------------------------------------- !--- input variables INTEGER, INTENT(IN ) :: i,j,nroot,nzs, iland REAL , & INTENT(IN ) :: GSWin, TABS, lai !--- soil properties REAL , & INTENT(IN ) :: DQM, & QMIN, & REF, & PC, & WILT REAL, DIMENSION(1:NZS), INTENT(IN) :: soiliqw, & ZSHALF !-- output REAL, DIMENSION(1:NZS), INTENT(OUT) :: TRANF REAL, INTENT(OUT) :: TRANSUM !-- local variables REAL :: totliq, did INTEGER :: k !-- for non-linear root distribution REAL :: gx,sm1,sm2,sm3,sm4,ap0,ap1,ap2,ap3,ap4 REAL :: FTEM, PCtot, fsol, f1, cmin, cmax, totcnd REAL, DIMENSION(1:NZS) :: PART !-------------------------------------------------------------------- do k=1,nzs part(k)=0. tranf(k)=0. enddo transum=0. totliq=soiliqw(1)+qmin sm1=totliq sm2=sm1*sm1 sm3=sm2*sm1 sm4=sm3*sm1 ap0=0.299 ap1=-8.152 ap2=61.653 ap3=-115.876 ap4=59.656 gx=ap0+ap1*sm1+ap2*sm2+ap3*sm3+ap4*sm4 if(totliq.ge.ref) gx=1. if(totliq.le.0.) gx=0. if(gx.gt.1.) gx=1. if(gx.lt.0.) gx=0. DID=zshalf(2) part(1)=DID*gx IF(TOTLIQ.GT.REF) THEN TRANF(1)=DID ELSE IF(TOTLIQ.LE.WILT) THEN TRANF(1)=0. ELSE TRANF(1)=(TOTLIQ-WILT)/(REF-WILT)*DID ENDIF !-- uncomment next line for non-linear root distribution ! TRANF(1)=part(1) DO K=2,NROOT totliq=soiliqw(k)+qmin sm1=totliq sm2=sm1*sm1 sm3=sm2*sm1 sm4=sm3*sm1 gx=ap0+ap1*sm1+ap2*sm2+ap3*sm3+ap4*sm4 if(totliq.ge.ref) gx=1. if(totliq.le.0.) gx=0. if(gx.gt.1.) gx=1. if(gx.lt.0.) gx=0. DID=zshalf(K+1)-zshalf(K) part(k)=did*gx IF(totliq.GE.REF) THEN TRANF(K)=DID ELSE IF(totliq.LE.WILT) THEN TRANF(K)=0. ELSE TRANF(K)=(totliq-WILT) & /(REF-WILT)*DID ENDIF !-- uncomment next line for non-linear root distribution ! TRANF(k)=part(k) END DO ! For LAI> 3 => transpiration at potential rate (F.Tardieu, 2013) if(lai > 4.) then pctot=0.8 else pctot=pc !- 26aug16- next 2 lines could lead to LH increase and higher 2-m Q during the day ! pctot=min(0.8,pc*lai) ! pctot=min(0.8,max(pc,pc*lai)) endif IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN ! if (i==421.and.j==280) then print *,'i,j,pctot,lai,pc',i,j,pctot,lai,pc ENDIF !--- !--- air temperature function ! Avissar (1985) and AX 7/95 IF (TABS .LE. 302.15) THEN FTEM = 1.0 / (1.0 + EXP(-0.41 * (TABS - 282.05))) ELSE FTEM = 1.0 / (1.0 + EXP(0.5 * (TABS - 314.0))) ENDIF IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN ! if (i==421.and.j==280) then print *,'i,j,tabs,ftem',i,j,tabs,ftem ENDIF !--- incoming solar function cmin = 1./rsmax_data cmax = 1./rstbl(iland) if(lai > 1.) then cmax = lai/rstbl(iland) ! max conductance endif ! Noihlan & Planton (1988) f1=0. ! if(lai > 0.01) then ! f1 = 1.1/lai*gswin/rgltbl(iland)! f1=0. when GSWin=0. ! fsol = (f1+cmin/cmax)/(1.+f1) ! fsol=min(1.,fsol) ! else ! fsol=cmin/cmax ! endif ! totcnd = max(lai/rstbl(iland), pctot * ftem * f1) ! Mahrer & Avissar (1982), Avissar et al. (1985) if (GSWin < rgltbl(iland)) then fsol = 1. / (1. + exp(-0.034 * (GSWin - 3.5))) else fsol = 1. endif IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN ! if (i==421.and.j==280) then print *,'i,j,GSWin,lai,f1,fsol',i,j,gswin,lai,f1,fsol ENDIF !--- total conductance totcnd =(cmin + (cmax - cmin)*pctot*ftem*fsol)/cmax IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN ! if (i==421.and.j==280) then print *,'i,j,iland,RGLTBL(iland),RSTBL(iland),RSMAX_DATA,totcnd' & ,i,j,iland,RGLTBL(iland),RSTBL(iland),RSMAX_DATA,totcnd ENDIF !-- TRANSUM - total for the rooting zone transum=0. DO K=1,NROOT ! linear root distribution TRANF(k)=max(cmin,TRANF(k)*totcnd) transum=transum+tranf(k) END DO IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN ! if (i==421.and.j==280) then print *,'i,j,transum,TRANF',i,j,transum,tranf endif !----------------------------------------------------------------- END SUBROUTINE TRANSF !----------------------------------------------------------------- SUBROUTINE VILKA(TN,D1,D2,PP,QS,TS,TT,NSTEP,ii,j,iland,isoil) !-------------------------------------------------------------- !--- VILKA finds the solution of energy budget at the surface !--- using table T,QS computed from Clausius-Klapeiron !-------------------------------------------------------------- REAL, DIMENSION(1:5001), INTENT(IN ) :: TT REAL, INTENT(IN ) :: TN,D1,D2,PP INTEGER, INTENT(IN ) :: NSTEP,ii,j,iland,isoil REAL, INTENT(OUT ) :: QS, TS REAL :: F1,T1,T2,RN INTEGER :: I,I1 I=(TN-1.7315E2)/.05+1 T1=173.1+FLOAT(I)*.05 F1=T1+D1*TT(I)-D2 I1=I-F1/(.05+D1*(TT(I+1)-TT(I))) I=I1 IF(I.GT.5000.OR.I.LT.1) GOTO 1 10 I1=I T1=173.1+FLOAT(I)*.05 F1=T1+D1*TT(I)-D2 RN=F1/(.05+D1*(TT(I+1)-TT(I))) I=I-INT(RN) IF(I.GT.5000.OR.I.LT.1) GOTO 1 IF(I1.NE.I) GOTO 10 TS=T1-.05*RN QS=(TT(I)+(TT(I)-TT(I+1))*RN)/PP GOTO 20 ! 1 PRINT *,'Crash in surface energy budget - STOP' 1 PRINT *,' AVOST IN VILKA Table index= ',I ! PRINT *,TN,D1,D2,PP,NSTEP,I,TT(i),ii,j,iland,isoil print *,'I,J=',ii,j,'LU_index = ',iland, 'Psfc[hPa] = ',pp, 'Tsfc = ',tn CALL wrf_error_fatal (' Crash in surface energy budget ' ) 20 CONTINUE !----------------------------------------------------------------------- END SUBROUTINE VILKA !----------------------------------------------------------------------- SUBROUTINE SOILVEGIN ( mosaic_lu,mosaic_soil,soilfrac,nscat, & shdmin, shdmax, & NLCAT,IVGTYP,ISLTYP,iswater, & IFOREST,lufrac,vegfrac,EMISS,PC,ZNT,LAI,RDLAI2D,& QWRTZ,RHOCS,BCLH,DQM,KSAT,PSIS,QMIN,REF,WILT,I,J) !************************************************************************ ! Set-up soil and vegetation Parameters in the case when ! snow disappears during the forecast and snow parameters ! shold be replaced by surface parameters according to ! soil and vegetation types in this point. ! ! Output: ! ! ! Soil parameters: ! DQM: MAX soil moisture content - MIN (m^3/m^3) ! REF: Reference soil moisture (m^3/m^3) ! WILT: Wilting PT soil moisture contents (m^3/m^3) ! QMIN: Air dry soil moist content limits (m^3/m^3) ! PSIS: SAT soil potential coefs. (m) ! KSAT: SAT soil diffusivity/conductivity coefs. (m/s) ! BCLH: Soil diffusivity/conductivity exponent. ! ! ************************************************************************ IMPLICIT NONE !--------------------------------------------------------------------------- integer, parameter :: nsoilclas=19 integer, parameter :: nvegclas=24+3 integer, parameter :: ilsnow=99 INTEGER, INTENT(IN ) :: nlcat, nscat, iswater, i, j !--- soiltyp classification according to STATSGO(nclasses=16) ! ! 1 SAND SAND ! 2 LOAMY SAND LOAMY SAND ! 3 SANDY LOAM SANDY LOAM ! 4 SILT LOAM SILTY LOAM ! 5 SILT SILTY LOAM ! 6 LOAM LOAM ! 7 SANDY CLAY LOAM SANDY CLAY LOAM ! 8 SILTY CLAY LOAM SILTY CLAY LOAM ! 9 CLAY LOAM CLAY LOAM ! 10 SANDY CLAY SANDY CLAY ! 11 SILTY CLAY SILTY CLAY ! 12 CLAY LIGHT CLAY ! 13 ORGANIC MATERIALS LOAM ! 14 WATER ! 15 BEDROCK ! Bedrock is reclassified as class 14 ! 16 OTHER (land-ice) ! 17 Playa ! 18 Lava ! 19 White Sand ! !---------------------------------------------------------------------- REAL LQMA(nsoilclas),LRHC(nsoilclas), & LPSI(nsoilclas),LQMI(nsoilclas), & LBCL(nsoilclas),LKAS(nsoilclas), & LWIL(nsoilclas),LREF(nsoilclas), & DATQTZ(nsoilclas) !-- LQMA Rawls et al.[1982] ! DATA LQMA /0.417, 0.437, 0.453, 0.501, 0.486, 0.463, 0.398, ! & 0.471, 0.464, 0.430, 0.479, 0.475, 0.439, 1.0, 0.20, 0.401/ !--- !-- Clapp, R. and G. Hornberger, 1978: Empirical equations for some soil ! hydraulic properties, Water Resour. Res., 14, 601-604. !-- Clapp et al. [1978] DATA LQMA /0.395, 0.410, 0.435, 0.485, 0.485, 0.451, 0.420, & 0.477, 0.476, 0.426, 0.492, 0.482, 0.451, 1.0, & 0.20, 0.435, 0.468, 0.200, 0.339/ !-- LREF Rawls et al.[1982] ! DATA LREF /0.091, 0.125, 0.207, 0.330, 0.360, 0.270, 0.255, ! & 0.366, 0.318, 0.339, 0.387, 0.396, 0.329, 1.0, 0.108, 0.283/ !-- Clapp et al. [1978] DATA LREF /0.174, 0.179, 0.249, 0.369, 0.369, 0.314, 0.299, & 0.357, 0.391, 0.316, 0.409, 0.400, 0.314, 1., & 0.1, 0.249, 0.454, 0.17, 0.236/ !-- LWIL Rawls et al.[1982] ! DATA LWIL/0.033, 0.055, 0.095, 0.133, 0.133, 0.117, 0.148, ! & 0.208, 0.197, 0.239, 0.250, 0.272, 0.066, 0.0, 0.006, 0.029/ !-- Clapp et al. [1978] DATA LWIL/0.068, 0.075, 0.114, 0.179, 0.179, 0.155, 0.175, & 0.218, 0.250, 0.219, 0.283, 0.286, 0.155, 0.0, & 0.006, 0.114, 0.030, 0.006, 0.01/ ! DATA LQMI/0.010, 0.028, 0.047, 0.084, 0.084, 0.066, 0.067, ! & 0.120, 0.103, 0.100, 0.126, 0.138, 0.066, 0.0, 0.006, 0.028/ !-- Carsel and Parrish [1988] DATA LQMI/0.045, 0.057, 0.065, 0.067, 0.034, 0.078, 0.10, & 0.089, 0.095, 0.10, 0.070, 0.068, 0.078, 0.0, & 0.004, 0.065, 0.020, 0.004, 0.008/ !-- LPSI Cosby et al[1984] ! DATA LPSI/0.060, 0.036, 0.141, 0.759, 0.759, 0.355, 0.135, ! & 0.617, 0.263, 0.098, 0.324, 0.468, 0.355, 0.0, 0.069, 0.036/ ! & 0.617, 0.263, 0.098, 0.324, 0.468, 0.355, 0.0, 0.069, 0.036/ !-- Clapp et al. [1978] DATA LPSI/0.121, 0.090, 0.218, 0.786, 0.786, 0.478, 0.299, & 0.356, 0.630, 0.153, 0.490, 0.405, 0.478, 0.0, & 0.121, 0.218, 0.468, 0.069, 0.069/ !-- LKAS Rawls et al.[1982] ! DATA LKAS/5.83E-5, 1.70E-5, 7.19E-6, 1.89E-6, 1.89E-6, ! & 3.67E-6, 1.19E-6, 4.17E-7, 6.39E-7, 3.33E-7, 2.50E-7, ! & 1.67E-7, 3.38E-6, 0.0, 1.41E-4, 1.41E-5/ !-- Clapp et al. [1978] DATA LKAS/1.76E-4, 1.56E-4, 3.47E-5, 7.20E-6, 7.20E-6, & 6.95E-6, 6.30E-6, 1.70E-6, 2.45E-6, 2.17E-6, & 1.03E-6, 1.28E-6, 6.95E-6, 0.0, 1.41E-4, & 3.47E-5, 1.28E-6, 1.41E-4, 1.76E-4/ !-- LBCL Cosby et al [1984] ! DATA LBCL/2.79, 4.26, 4.74, 5.33, 5.33, 5.25, 6.66, ! & 8.72, 8.17, 10.73, 10.39, 11.55, 5.25, 0.0, 2.79, 4.26/ !-- Clapp et al. [1978] DATA LBCL/4.05, 4.38, 4.90, 5.30, 5.30, 5.39, 7.12, & 7.75, 8.52, 10.40, 10.40, 11.40, 5.39, 0.0, & 4.05, 4.90, 11.55, 2.79, 2.79/ DATA LRHC /1.47,1.41,1.34,1.27,1.27,1.21,1.18,1.32,1.23, & 1.18,1.15,1.09,1.21,4.18,2.03,2.10,1.09,2.03,1.47/ DATA DATQTZ/0.92,0.82,0.60,0.25,0.10,0.40,0.60,0.10,0.35, & 0.52,0.10,0.25,0.00,0.,0.60,0.0,0.25,0.60,0.92/ !-------------------------------------------------------------------------- ! ! USGS Vegetation Types ! ! 1: Urban and Built-Up Land ! 2: Dryland Cropland and Pasture ! 3: Irrigated Cropland and Pasture ! 4: Mixed Dryland/Irrigated Cropland and Pasture ! 5: Cropland/Grassland Mosaic ! 6: Cropland/Woodland Mosaic ! 7: Grassland ! 8: Shrubland ! 9: Mixed Shrubland/Grassland ! 10: Savanna ! 11: Deciduous Broadleaf Forest ! 12: Deciduous Needleleaf Forest ! 13: Evergreen Broadleaf Forest ! 14: Evergreen Needleleaf Fores ! 15: Mixed Forest ! 16: Water Bodies ! 17: Herbaceous Wetland ! 18: Wooded Wetland ! 19: Barren or Sparsely Vegetated ! 20: Herbaceous Tundra ! 21: Wooded Tundra ! 22: Mixed Tundra ! 23: Bare Ground Tundra ! 24: Snow or Ice ! ! 25: Playa ! 26: Lava ! 27: White Sand ! MODIS vegetation categories from VEGPARM.TBL ! 1: Evergreen Needleleaf Forest ! 2: Evergreen Broadleaf Forest ! 3: Deciduous Needleleaf Forest ! 4: Deciduous Broadleaf Forest ! 5: Mixed Forests ! 6: Closed Shrublands ! 7: Open Shrublands ! 8: Woody Savannas ! 9: Savannas ! 10: Grasslands ! 11: Permanent wetlands ! 12: Croplands ! 13: Urban and Built-Up ! 14: cropland/natural vegetation mosaic ! 15: Snow and Ice ! 16: Barren or Sparsely Vegetated ! 17: Water ! 18: Wooded Tundra ! 19: Mixed Tundra ! 20: Barren Tundra ! 21: Lakes !---- Below are the arrays for the vegetation parameters REAL LALB(nvegclas),LMOI(nvegclas),LEMI(nvegclas), & LROU(nvegclas),LTHI(nvegclas),LSIG(nvegclas), & LPC(nvegclas) !************************************************************************ !---- vegetation parameters ! !-- USGS model ! DATA LALB/.18,.17,.18,.18,.18,.16,.19,.22,.20,.20,.16,.14, & .12,.12,.13,.08,.14,.14,.25,.15,.15,.15,.25,.55, & .30,.16,.60 / DATA LEMI/.88,4*.92,.93,.92,.88,.9,.92,.93,.94, & .95,.95,.94,.98,.95,.95,.85,.92,.93,.92,.85,.95, & .85,.85,.90 / !-- Roughness length is changed for forests and some others ! DATA LROU/.5,.06,.075,.065,.05,.2,.075,.1,.11,.15,.8,.85, & ! 2.0,1.0,.563,.0001,.2,.4,.05,.1,.15,.1,.065,.05/ DATA LROU/.5,.06,.075,.065,.05,.2,.075,.1,.11,.15,.5,.5, & .5,.5,.5,.0001,.2,.4,.05,.1,.15,.1,.065,.05, & .01,.15,.01 / DATA LMOI/.1,.3,.5,.25,.25,.35,.15,.1,.15,.15,.3,.3, & .5,.3,.3,1.,.6,.35,.02,.5,.5,.5,.02,.95,.40,.50,.40/ ! !---- still needs to be corrected ! ! DATA LPC/ 15*.8,0.,.8,.8,.5,.5,.5,.5,.5,.0/ 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, & 0.3,0.3,0.4,0.4,0.3,0.,.3,0.,0./ ! used in RUC DATA LPC /0.6,6*0.8,0.7,0.75,6*0.8,0.,0.8,0.8, & ! 0.5,0.7,0.6,0.7,0.5,0./ !*************************************************************************** INTEGER :: & IVGTYP, & ISLTYP INTEGER, INTENT(IN ) :: mosaic_lu, mosaic_soil REAL, INTENT(IN ) :: SHDMAX REAL, INTENT(IN ) :: SHDMIN REAL, INTENT(IN ) :: VEGFRAC REAL, DIMENSION( 1:NLCAT ), INTENT(IN):: LUFRAC REAL, DIMENSION( 1:NSCAT ), INTENT(IN):: SOILFRAC REAL , & INTENT ( OUT) :: pc REAL , & INTENT (INOUT ) :: emiss, & lai, & znt LOGICAL, intent(in) :: rdlai2d !--- soil properties REAL , & INTENT( OUT) :: RHOCS, & BCLH, & DQM, & KSAT, & PSIS, & QMIN, & QWRTZ, & REF, & WILT INTEGER, INTENT ( OUT) :: iforest ! INTEGER, DIMENSION( 1:(lucats) ) , & ! INTENT ( OUT) :: iforest ! INTEGER, DIMENSION( 1:50 ) :: if1 INTEGER :: kstart, kfin, lstart, lfin INTEGER :: k REAL :: area, factor, znt1, lb REAL, DIMENSION( 1:NLCAT ) :: ZNTtoday, LAItoday, deltalai !*********************************************************************** ! DATA ZS1/0.0,0.05,0.20,0.40,1.6,3.0/ ! o - levels in soil ! DATA ZS2/0.0,0.025,0.125,0.30,1.,2.3/ ! x - levels in soil ! DATA IF1/12*0,1,1,1,12*0/ ! do k=1,LUCATS ! iforest(k)=if1(k) ! enddo iforest = IFORTBL(IVGTYP) IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN if(i.eq.375.and.j.eq.254)then print *,'ifortbl(ivgtyp),ivgtyp,laitbl(ivgtyp),z0tbl(ivgtyp)', & ifortbl(ivgtyp),ivgtyp,laitbl(ivgtyp),z0tbl(ivgtyp) endif ENDIF deltalai(:) = 0. ! 11oct2012 - seasonal correction on ZNT for crops and LAI for all veg. types ! factor = 1 with minimum greenness --> vegfrac = shdmin (cold season) ! factor = 0 with maximum greenness --> vegfrac = shdmax ! SHDMAX, SHDMIN and VEGFRAC are in % here. if((shdmax - shdmin) .lt. 1) then factor = 1. ! min greenness else factor = 1. - max(0.,min(1.,(vegfrac - shdmin)/max(1.,(shdmax-shdmin)))) endif ! 18sept18 - LAITBL and Z0TBL are the max values do k = 1,nlcat if(IFORTBL(k) == 1) deltalai(k)=min(0.2,0.8*LAITBL(K)) if(IFORTBL(k) == 2 .or. IFORTBL(k) == 7) deltalai(k)=min(0.5,0.8*LAITBL(K)) if(IFORTBL(k) == 3) deltalai(k)=min(0.45,0.8*LAITBL(K)) if(IFORTBL(k) == 4) deltalai(k)=min(0.75,0.8*LAITBL(K)) if(IFORTBL(k) == 5) deltalai(k)=min(0.86,0.8*LAITBL(K)) if(k.ne.iswater) then !-- 20aug18 - change in LAItoday based on the greenness fraction for the current day LAItoday(k) = LAITBL(K) - deltalai(k) * factor if(IFORTBL(k) == 7) then !-- seasonal change of roughness length for crops ZNTtoday(k) = Z0TBL(K) - 0.125 * factor else ZNTtoday(k) = Z0TBL(K) endif else LAItoday(k) = LAITBL(K) ! ZNTtoday(k) = Z0TBL(K) ZNTtoday(k) = ZNT ! do not overwrite z0 over water with the table value endif enddo IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN if(i.eq.539.and.j.eq.334)then print *,'ivgtyp,factor,vegfrac,shdmin,shdmax,deltalai(ivgtyp),laitoday(ivgtyp),znttoday(ivgtyp)', & i,j,ivgtyp,factor,vegfrac,shdmin,shdmax,deltalai(ivgtyp),laitoday(ivgtyp),znttoday(ivgtyp) endif ENDIF EMISS = 0. ZNT = 0. ZNT1 = 0. PC = 0. if(.not.rdlai2d) LAI = 0. AREA = 0. !-- mosaic approach to landuse in the grid box ! Use Mason (1988) Eq.(15) to compute effective ZNT; ! Lb - blending height = L/200., where L is the length scale ! of regions with varying Z0 (Lb = 5 if L=1000 m) LB = 5. if(mosaic_lu == 1) then do k = 1,nlcat AREA = AREA + lufrac(k) EMISS = EMISS+ LEMITBL(K)*lufrac(k) ZNT = ZNT + lufrac(k)/ALOG(LB/ZNTtoday(K))**2. ! ZNT1 - weighted average in the grid box, not used, computed for comparison ZNT1 = ZNT1 + lufrac(k)*ZNTtoday(K) if(.not.rdlai2d) LAI = LAI + LAItoday(K)*lufrac(k) PC = PC + PCTBL(K)*lufrac(k) enddo if (area.gt.1.) area=1. if (area <= 0.) then print *,'Bad area of grid box', area stop endif IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN if(i.eq.358.and.j.eq.260) then print *,'area=',area,i,j,ivgtyp,nlcat,(lufrac(k),k=1,nlcat),EMISS,ZNT,ZNT1,LAI,PC endif ENDIF EMISS = EMISS/AREA ZNT1 = ZNT1/AREA ZNT = LB/EXP(SQRT(1./ZNT)) if(.not.rdlai2d) LAI = LAI/AREA PC = PC /AREA IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN if(i.eq.358.and.j.eq.260) then print *,'mosaic=',i,j,ivgtyp,nlcat,(lufrac(k),k=1,nlcat),EMISS,ZNT,ZNT1,LAI,PC endif ENDIF else EMISS = LEMITBL(IVGTYP) ZNT = ZNTtoday(IVGTYP) PC = PCTBL(IVGTYP) if(.not.rdlai2d) LAI = LAItoday(IVGTYP) endif ! parameters from SOILPARM.TBL RHOCS = 0. BCLH = 0. DQM = 0. KSAT = 0. PSIS = 0. QMIN = 0. REF = 0. WILT = 0. QWRTZ = 0. AREA = 0. ! mosaic approach if(mosaic_soil == 1 ) then do k = 1, nscat if(k.ne.14) then ! STATSGO value for water !exclude watrer points from this loop AREA = AREA + soilfrac(k) RHOCS = RHOCS + HC(k)*1.E6*soilfrac(k) BCLH = BCLH + BB(K)*soilfrac(k) DQM = DQM + (MAXSMC(K)- & DRYSMC(K))*soilfrac(k) KSAT = KSAT + SATDK(K)*soilfrac(k) PSIS = PSIS - SATPSI(K)*soilfrac(k) QMIN = QMIN + DRYSMC(K)*soilfrac(k) REF = REF + REFSMC(K)*soilfrac(k) WILT = WILT + WLTSMC(K)*soilfrac(k) QWRTZ = QWRTZ + QTZ(K)*soilfrac(k) endif enddo if (area.gt.1.) area=1. if (area <= 0.) then ! area = 0. for water points ! print *,'Area of a grid box', area, 'iswater = ',iswater RHOCS = HC(ISLTYP)*1.E6 BCLH = BB(ISLTYP) DQM = MAXSMC(ISLTYP)- & DRYSMC(ISLTYP) KSAT = SATDK(ISLTYP) PSIS = - SATPSI(ISLTYP) QMIN = DRYSMC(ISLTYP) REF = REFSMC(ISLTYP) WILT = WLTSMC(ISLTYP) QWRTZ = QTZ(ISLTYP) else RHOCS = RHOCS/AREA BCLH = BCLH/AREA DQM = DQM/AREA KSAT = KSAT/AREA PSIS = PSIS/AREA QMIN = QMIN/AREA REF = REF/AREA WILT = WILT/AREA QWRTZ = QWRTZ/AREA endif ! dominant category approach else if(isltyp.ne.14) then RHOCS = HC(ISLTYP)*1.E6 BCLH = BB(ISLTYP) DQM = MAXSMC(ISLTYP)- & DRYSMC(ISLTYP) KSAT = SATDK(ISLTYP) PSIS = - SATPSI(ISLTYP) QMIN = DRYSMC(ISLTYP) REF = REFSMC(ISLTYP) WILT = WLTSMC(ISLTYP) QWRTZ = QTZ(ISLTYP) endif endif ! parameters from the look-up tables ! BCLH = LBCL(ISLTYP) ! DQM = LQMA(ISLTYP)- & ! LQMI(ISLTYP) ! KSAT = LKAS(ISLTYP) ! PSIS = - LPSI(ISLTYP) ! QMIN = LQMI(ISLTYP) ! REF = LREF(ISLTYP) ! WILT = LWIL(ISLTYP) ! QWRTZ = DATQTZ(ISLTYP) !-------------------------------------------------------------------------- END SUBROUTINE SOILVEGIN !-------------------------------------------------------------------------- SUBROUTINE RUCLSMINIT( SH2O,SMFR3D,TSLB,SMOIS,ISLTYP,IVGTYP, & mminlu, XICE,mavail,nzs, iswater, isice, & znt, restart, allowed_to_read , & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) #if ( WRF_CHEM == 1 ) USE module_data_gocart_dust #endif IMPLICIT NONE INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte, & nzs, iswater, isice CHARACTER(LEN=*), INTENT(IN ) :: MMINLU REAL, DIMENSION( ims:ime, 1:nzs, jms:jme ) , & INTENT(IN) :: TSLB, & SMOIS INTEGER, DIMENSION( ims:ime, jms:jme ) , & INTENT(INOUT) :: ISLTYP,IVGTYP REAL, DIMENSION( ims:ime, 1:nzs, jms:jme ) , & INTENT(INOUT) :: SMFR3D, & SH2O REAL, DIMENSION( ims:ime, jms:jme ) , & INTENT(INOUT) :: XICE,MAVAIL REAL, DIMENSION( ims:ime, jms:jme ) , & INTENT( OUT) :: znt REAL, DIMENSION ( 1:nzs ) :: SOILIQW LOGICAL , INTENT(IN) :: restart, allowed_to_read ! INTEGER :: I,J,L,itf,jtf REAL :: RIW,XLMELT,TLN,DQM,REF,PSIS,QMIN,BCLH character*8 :: MMINLURUC, MMINSL INTEGER :: errflag ! itf=min0(ite,ide-1) ! jtf=min0(jte,jde-1) RIW=900.*1.e-3 XLMELT=3.35E+5 ! initialize three LSM related tables IF ( allowed_to_read ) THEN CALL wrf_message( 'INITIALIZE THREE LSM RELATED TABLES' ) if(mminlu == 'USGS') then MMINLURUC='USGS-RUC' elseif(mminlu == 'MODIS' .OR. & & mminlu == 'MODIFIED_IGBP_MODIS_NOAH') then MMINLURUC='MODI-RUC' endif MMINSL='STAS-RUC' ! CALL RUCLSM_PARM_INIT print *,'RUCLSMINIT uses ',mminluruc call RUCLSM_SOILVEGPARM( MMINLURUC, MMINSL) ENDIF !#if ( WRF_CHEM == 1 ) ! ! need this parameter for dust parameterization in wrf/chem ! ! do I=1,NSLTYPE ! porosity(i)=maxsmc(i) ! drypoint(i)=drysmc(i) ! enddo !#endif ! IF(.not.restart)THEN itf=min0(ite,ide-1) jtf=min0(jte,jde-1) errflag = 0 DO j = jts,jtf DO i = its,itf IF ( ISLTYP( i,j ) .LT. 1 ) THEN errflag = 1 WRITE(err_message,*)"module_sf_ruclsm.F: lsminit: out of range ISLTYP ",i,j,ISLTYP( i,j ) CALL wrf_message(err_message) ENDIF ENDDO ENDDO IF ( errflag .EQ. 1 ) THEN CALL wrf_error_fatal( "module_sf_ruclsm.F: lsminit: out of range value "// & "of ISLTYP. Is this field in the input?" ) ENDIF DO J=jts,jtf DO I=its,itf ZNT(I,J) = Z0TBL(IVGTYP(I,J)) ! CALL SOILIN ( ISLTYP(I,J), DQM, REF, PSIS, QMIN, BCLH ) !--- Computation of volumetric content of ice in soil !--- and initialize MAVAIL DQM = MAXSMC (ISLTYP(I,J)) - & DRYSMC (ISLTYP(I,J)) REF = REFSMC (ISLTYP(I,J)) PSIS = - SATPSI (ISLTYP(I,J)) QMIN = DRYSMC (ISLTYP(I,J)) BCLH = BB (ISLTYP(I,J)) !!! IF (.not.restart) THEN IF(xice(i,j).gt.0.) THEN !-- for ice DO L=1,NZS smfr3d(i,l,j)=1. sh2o(i,l,j)=0. mavail(i,j) = 1. ENDDO ELSE if(isltyp(i,j).ne.14 ) then !-- land mavail(i,j) = max(0.00001,min(1.,(smois(i,1,j)-qmin)/(ref-qmin))) DO L=1,NZS !-- for land points initialize soil ice tln=log(TSLB(i,l,j)/273.15) if(tln.lt.0.) then soiliqw(l)=(dqm+qmin)*(XLMELT* & (tslb(i,l,j)-273.15)/tslb(i,l,j)/9.81/psis) & **(-1./bclh) ! **(-1./bclh)-qmin soiliqw(l)=max(0.,soiliqw(l)) soiliqw(l)=min(soiliqw(l),smois(i,l,j)) sh2o(i,l,j)=soiliqw(l) smfr3d(i,l,j)=(smois(i,l,j)-soiliqw(l))/RIW else smfr3d(i,l,j)=0. sh2o(i,l,j)=smois(i,l,j) endif ENDDO else !-- for water ISLTYP=14 DO L=1,NZS smfr3d(i,l,j)=0. sh2o(i,l,j)=1. mavail(i,j) = 1. ENDDO endif ENDIF ENDDO ENDDO ENDIF END SUBROUTINE ruclsminit ! !----------------------------------------------------------------- ! SUBROUTINE RUCLSM_PARM_INIT !----------------------------------------------------------------- ! character*9 :: MMINLU, MMINSL ! MMINLU='MODIS-RUC' ! MMINLU='USGS-RUC' ! MMINSL='STAS-RUC' ! call RUCLSM_SOILVEGPARM( MMINLU, MMINSL) !----------------------------------------------------------------- ! END SUBROUTINE RUCLSM_PARM_INIT !----------------------------------------------------------------- !----------------------------------------------------------------- SUBROUTINE RUCLSM_SOILVEGPARM( MMINLURUC, MMINSL) !----------------------------------------------------------------- IMPLICIT NONE integer :: LUMATCH, IINDEX, LC, NUM_SLOPE integer :: ierr INTEGER , PARAMETER :: OPEN_OK = 0 character*8 :: MMINLURUC, MMINSL character*128 :: mess , message, vege_parm_string logical, external :: wrf_dm_on_monitor !-----SPECIFY VEGETATION RELATED CHARACTERISTICS : ! ALBBCK: SFC albedo (in percentage) ! Z0: Roughness length (m) ! LEMI: Emissivity ! PC: Plant coefficient for transpiration function ! -- the rest of the parameters are read in but not used currently ! SHDFAC: Green vegetation fraction (in percentage) ! Note: The ALBEDO, Z0, and SHDFAC values read from the following table ! ALBEDO, amd Z0 are specified in LAND-USE TABLE; and SHDFAC is ! the monthly green vegetation data ! CMXTBL: MAX CNPY Capacity (m) ! RSMIN: Mimimum stomatal resistance (s m-1) ! RSMAX: Max. stomatal resistance (s m-1) ! RGL: Parameters used in radiation stress function ! HS: Parameter used in vapor pressure deficit functio ! TOPT: Optimum transpiration air temperature. (K) ! CMCMAX: Maximum canopy water capacity ! CFACTR: Parameter used in the canopy inteception calculati ! SNUP: Threshold snow depth (in water equivalent m) that ! implies 100% snow cover ! LAI: Leaf area index (dimensionless) ! MAXALB: Upper bound on maximum albedo over deep snow ! !-----READ IN VEGETAION PROPERTIES FROM VEGPARM.TBL ! IF ( wrf_dm_on_monitor() ) THEN OPEN(19, FILE='VEGPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr) IF(ierr .NE. OPEN_OK ) THEN WRITE(message,FMT='(A)') & 'module_sf_ruclsm.F: soil_veg_gen_parm: failure opening VEGPARM.TBL' CALL wrf_error_fatal ( message ) END IF WRITE ( mess, * ) 'INPUT VEGPARM FOR ',MMINLURUC CALL wrf_message( mess ) LUMATCH=0 2000 FORMAT (A8) READ (19,'(A)') vege_parm_string outer : DO READ (19,2000,END=2002)LUTYPE READ (19,*)LUCATS,IINDEX WRITE( mess , * ) 'VEGPARM FOR ',LUTYPE,' FOUND', LUCATS,' CATEGORIES' CALL wrf_message( mess ) IF(LUTYPE.NE.MMINLURUC)THEN ! Skip over the undesired table write ( mess , * ) 'Skipping ', LUTYPE, ' table' CALL wrf_message( mess ) DO LC=1,LUCATS READ (19,*) ENDDO inner : DO ! Find the next "Vegetation Parameters" READ (19,'(A)',END=2002) vege_parm_string IF (TRIM(vege_parm_string) .EQ. "Vegetation Parameters") THEN EXIT inner END IF ENDDO inner ELSE LUMATCH=1 write ( mess , * ) 'Found ', LUTYPE, ' table' CALL wrf_message( mess ) EXIT outer ! Found the table, read the data END IF ENDDO outer IF (LUMATCH == 1) then write ( mess , * ) 'Reading ',LUTYPE,' table' CALL wrf_message( mess ) DO LC=1,LUCATS READ (19,*)IINDEX,ALBTBL(LC),Z0TBL(LC),LEMITBL(LC),PCTBL(LC), & SHDTBL(LC),IFORTBL(LC),RSTBL(LC),RGLTBL(LC), & HSTBL(LC),SNUPTBL(LC),LAITBL(LC),MAXALB(LC) ENDDO ! READ (19,*) READ (19,*)TOPT_DATA READ (19,*) READ (19,*)CMCMAX_DATA READ (19,*) READ (19,*)CFACTR_DATA READ (19,*) READ (19,*)RSMAX_DATA READ (19,*) READ (19,*)BARE READ (19,*) READ (19,*)NATURAL READ (19,*) READ (19,*)CROP READ (19,*) READ (19,*,iostat=ierr)URBAN if ( ierr /= 0 ) call wrf_message ( "-------- VEGPARM.TBL READ ERROR --------") if ( ierr /= 0 ) call wrf_message ( "Problem read URBAN from VEGPARM.TBL") if ( ierr /= 0 ) call wrf_message ( " -- Use updated version of VEGPARM.TBL ") if ( ierr /= 0 ) call wrf_error_fatal ( "Problem read URBAN from VEGPARM.TBL") ENDIF 2002 CONTINUE CLOSE (19) !----- IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN print *,' LEMITBL, PCTBL, Z0TBL, LAITBL --->', LEMITBL, PCTBL, Z0TBL, LAITBL ENDIF IF (LUMATCH == 0) then CALL wrf_error_fatal ("Land Use Dataset '"//MMINLURUC//"' not found in VEGPARM.TBL.") ENDIF END IF CALL wrf_dm_bcast_string ( LUTYPE , 8 ) CALL wrf_dm_bcast_integer ( LUCATS , 1 ) CALL wrf_dm_bcast_integer ( IINDEX , 1 ) CALL wrf_dm_bcast_integer ( LUMATCH , 1 ) CALL wrf_dm_bcast_real ( ALBTBL , NLUS ) CALL wrf_dm_bcast_real ( Z0TBL , NLUS ) CALL wrf_dm_bcast_real ( LEMITBL , NLUS ) CALL wrf_dm_bcast_real ( PCTBL , NLUS ) CALL wrf_dm_bcast_real ( SHDTBL , NLUS ) CALL wrf_dm_bcast_real ( IFORTBL , NLUS ) CALL wrf_dm_bcast_real ( RSTBL , NLUS ) CALL wrf_dm_bcast_real ( RGLTBL , NLUS ) CALL wrf_dm_bcast_real ( HSTBL , NLUS ) CALL wrf_dm_bcast_real ( SNUPTBL , NLUS ) CALL wrf_dm_bcast_real ( LAITBL , NLUS ) CALL wrf_dm_bcast_real ( MAXALB , NLUS ) CALL wrf_dm_bcast_real ( TOPT_DATA , 1 ) CALL wrf_dm_bcast_real ( CMCMAX_DATA , 1 ) CALL wrf_dm_bcast_real ( CFACTR_DATA , 1 ) CALL wrf_dm_bcast_real ( RSMAX_DATA , 1 ) CALL wrf_dm_bcast_integer ( BARE , 1 ) CALL wrf_dm_bcast_integer ( NATURAL , 1 ) CALL wrf_dm_bcast_integer ( CROP , 1 ) CALL wrf_dm_bcast_integer ( URBAN , 1 ) ! !-----READ IN SOIL PROPERTIES FROM SOILPARM.TBL ! IF ( wrf_dm_on_monitor() ) THEN OPEN(19, FILE='SOILPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr) IF(ierr .NE. OPEN_OK ) THEN WRITE(message,FMT='(A)') & 'module_sf_ruclsm.F: soil_veg_gen_parm: failure opening SOILPARM.TBL' CALL wrf_error_fatal ( message ) END IF WRITE(mess,*) 'INPUT SOIL TEXTURE CLASSIFICATION = ',MMINSL CALL wrf_message( mess ) LUMATCH=0 READ (19,*) READ (19,2000,END=2003)SLTYPE READ (19,*)SLCATS,IINDEX IF(SLTYPE.NE.MMINSL)THEN DO LC=1,SLCATS READ (19,*) IINDEX,BB(LC),DRYSMC(LC),HC(LC),MAXSMC(LC),& REFSMC(LC),SATPSI(LC),SATDK(LC), SATDW(LC), & WLTSMC(LC), QTZ(LC) ENDDO ENDIF READ (19,*) READ (19,2000,END=2003)SLTYPE READ (19,*)SLCATS,IINDEX IF(SLTYPE.EQ.MMINSL)THEN WRITE( mess , * ) 'SOIL TEXTURE CLASSIFICATION = ',SLTYPE,' FOUND', & SLCATS,' CATEGORIES' CALL wrf_message ( mess ) LUMATCH=1 ENDIF IF(SLTYPE.EQ.MMINSL)THEN DO LC=1,SLCATS READ (19,*) IINDEX,BB(LC),DRYSMC(LC),HC(LC),MAXSMC(LC),& REFSMC(LC),SATPSI(LC),SATDK(LC), SATDW(LC), & WLTSMC(LC), QTZ(LC) ENDDO ENDIF 2003 CONTINUE CLOSE (19) ENDIF CALL wrf_dm_bcast_integer ( LUMATCH , 1 ) CALL wrf_dm_bcast_string ( SLTYPE , 8 ) CALL wrf_dm_bcast_string ( MMINSL , 8 ) ! since this is reset above, see oct2 ^ CALL wrf_dm_bcast_integer ( SLCATS , 1 ) CALL wrf_dm_bcast_integer ( IINDEX , 1 ) CALL wrf_dm_bcast_real ( BB , NSLTYPE ) CALL wrf_dm_bcast_real ( DRYSMC , NSLTYPE ) CALL wrf_dm_bcast_real ( HC , NSLTYPE ) CALL wrf_dm_bcast_real ( MAXSMC , NSLTYPE ) CALL wrf_dm_bcast_real ( REFSMC , NSLTYPE ) CALL wrf_dm_bcast_real ( SATPSI , NSLTYPE ) CALL wrf_dm_bcast_real ( SATDK , NSLTYPE ) CALL wrf_dm_bcast_real ( SATDW , NSLTYPE ) CALL wrf_dm_bcast_real ( WLTSMC , NSLTYPE ) CALL wrf_dm_bcast_real ( QTZ , NSLTYPE ) IF(LUMATCH.EQ.0)THEN CALL wrf_message( 'SOIl TEXTURE IN INPUT FILE DOES NOT ' ) CALL wrf_message( 'MATCH SOILPARM TABLE' ) CALL wrf_error_fatal ( 'INCONSISTENT OR MISSING SOILPARM FILE' ) ENDIF ! !-----READ IN GENERAL PARAMETERS FROM GENPARM.TBL ! IF ( wrf_dm_on_monitor() ) THEN OPEN(19, FILE='GENPARM.TBL',FORM='FORMATTED',STATUS='OLD',IOSTAT=ierr) IF(ierr .NE. OPEN_OK ) THEN WRITE(message,FMT='(A)') & 'module_sf_ruclsm.F: soil_veg_gen_parm: failure opening GENPARM.TBL' CALL wrf_error_fatal ( message ) END IF READ (19,*) READ (19,*) READ (19,*) NUM_SLOPE SLPCATS=NUM_SLOPE DO LC=1,SLPCATS READ (19,*)SLOPE_DATA(LC) ENDDO READ (19,*) READ (19,*)SBETA_DATA READ (19,*) READ (19,*)FXEXP_DATA READ (19,*) READ (19,*)CSOIL_DATA READ (19,*) READ (19,*)SALP_DATA READ (19,*) READ (19,*)REFDK_DATA READ (19,*) READ (19,*)REFKDT_DATA READ (19,*) READ (19,*)FRZK_DATA READ (19,*) READ (19,*)ZBOT_DATA READ (19,*) READ (19,*)CZIL_DATA READ (19,*) READ (19,*)SMLOW_DATA READ (19,*) READ (19,*)SMHIGH_DATA CLOSE (19) ENDIF CALL wrf_dm_bcast_integer ( NUM_SLOPE , 1 ) CALL wrf_dm_bcast_integer ( SLPCATS , 1 ) CALL wrf_dm_bcast_real ( SLOPE_DATA , NSLOPE ) CALL wrf_dm_bcast_real ( SBETA_DATA , 1 ) CALL wrf_dm_bcast_real ( FXEXP_DATA , 1 ) CALL wrf_dm_bcast_real ( CSOIL_DATA , 1 ) CALL wrf_dm_bcast_real ( SALP_DATA , 1 ) CALL wrf_dm_bcast_real ( REFDK_DATA , 1 ) CALL wrf_dm_bcast_real ( REFKDT_DATA , 1 ) CALL wrf_dm_bcast_real ( FRZK_DATA , 1 ) CALL wrf_dm_bcast_real ( ZBOT_DATA , 1 ) CALL wrf_dm_bcast_real ( CZIL_DATA , 1 ) CALL wrf_dm_bcast_real ( SMLOW_DATA , 1 ) CALL wrf_dm_bcast_real ( SMHIGH_DATA , 1 ) !----------------------------------------------------------------- END SUBROUTINE RUCLSM_SOILVEGPARM !----------------------------------------------------------------- SUBROUTINE SOILIN (ISLTYP, DQM, REF, PSIS, QMIN, BCLH ) !--- soiltyp classification according to STATSGO(nclasses=16) ! ! 1 SAND SAND ! 2 LOAMY SAND LOAMY SAND ! 3 SANDY LOAM SANDY LOAM ! 4 SILT LOAM SILTY LOAM ! 5 SILT SILTY LOAM ! 6 LOAM LOAM ! 7 SANDY CLAY LOAM SANDY CLAY LOAM ! 8 SILTY CLAY LOAM SILTY CLAY LOAM ! 9 CLAY LOAM CLAY LOAM ! 10 SANDY CLAY SANDY CLAY ! 11 SILTY CLAY SILTY CLAY ! 12 CLAY LIGHT CLAY ! 13 ORGANIC MATERIALS LOAM ! 14 WATER ! 15 BEDROCK ! Bedrock is reclassified as class 14 ! 16 OTHER (land-ice) ! extra classes from Fei Chen ! 17 Playa ! 18 Lava ! 19 White Sand ! !---------------------------------------------------------------------- integer, parameter :: nsoilclas=19 integer, intent ( in) :: isltyp real, intent ( out) :: dqm,ref,qmin,psis REAL LQMA(nsoilclas),LREF(nsoilclas),LBCL(nsoilclas), & LPSI(nsoilclas),LQMI(nsoilclas) !-- LQMA Rawls et al.[1982] ! DATA LQMA /0.417, 0.437, 0.453, 0.501, 0.486, 0.463, 0.398, ! & 0.471, 0.464, 0.430, 0.479, 0.475, 0.439, 1.0, 0.20, 0.401/ !--- !-- Clapp, R. and G. Hornberger, Empirical equations for some soil ! hydraulic properties, Water Resour. Res., 14,601-604,1978. !-- Clapp et al. [1978] DATA LQMA /0.395, 0.410, 0.435, 0.485, 0.485, 0.451, 0.420, & 0.477, 0.476, 0.426, 0.492, 0.482, 0.451, 1.0, & 0.20, 0.435, 0.468, 0.200, 0.339/ !-- Clapp et al. [1978] DATA LREF /0.174, 0.179, 0.249, 0.369, 0.369, 0.314, 0.299, & 0.357, 0.391, 0.316, 0.409, 0.400, 0.314, 1., & 0.1, 0.249, 0.454, 0.17, 0.236/ !-- Carsel and Parrish [1988] DATA LQMI/0.045, 0.057, 0.065, 0.067, 0.034, 0.078, 0.10, & 0.089, 0.095, 0.10, 0.070, 0.068, 0.078, 0.0, & 0.004, 0.065, 0.020, 0.004, 0.008/ !-- Clapp et al. [1978] DATA LPSI/0.121, 0.090, 0.218, 0.786, 0.786, 0.478, 0.299, & 0.356, 0.630, 0.153, 0.490, 0.405, 0.478, 0.0, & 0.121, 0.218, 0.468, 0.069, 0.069/ !-- Clapp et al. [1978] DATA LBCL/4.05, 4.38, 4.90, 5.30, 5.30, 5.39, 7.12, & 7.75, 8.52, 10.40, 10.40, 11.40, 5.39, 0.0, & 4.05, 4.90, 11.55, 2.79, 2.79/ DQM = LQMA(ISLTYP)- & LQMI(ISLTYP) REF = LREF(ISLTYP) PSIS = - LPSI(ISLTYP) QMIN = LQMI(ISLTYP) BCLH = LBCL(ISLTYP) END SUBROUTINE SOILIN END MODULE module_sf_ruclsm