/wrfv2_fire/phys/module_sf_ruclsm.F
FORTRAN Legacy | 6031 lines | 3585 code | 793 blank | 1653 comment | 96 complexity | 926e2faf7df8889ef7fefb2ef52fc9ed MD5 | raw file
Possible License(s): AGPL-1.0
- #define LSMRUC_DBG_LVL 3000
- !WRF:MODEL_LAYER:PHYSICS
- !
- MODULE module_sf_ruclsm
- USE module_model_constants
- USE module_wrf_error
- ! VEGETATION PARAMETERS
- INTEGER :: LUCATS , BARE, NATURAL
- 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( &
- DT,KTAU,NSL,ZS, &
- RAINBL,SNOW,SNOWH,SNOWC,FRZFRAC,frpcpn, &
- 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,SFCEXC, &
- SFCEVP,GRDFLX,ACSNOW,SNOM, &
- SMFR3D,KEEPFR3DFLAG, &
- myj, &
- 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
- !-- 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 1)
- ! 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]
- ! SFCEVP - total evaporation in [kg/m^2]
- ! GRDFLX - soil heat flux (W/m^2: negative, if downward from surface)
- ! ACSNOW - accumulation of snow water [m]
- !-- 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 ) :: myj,frpcpn
- 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
- 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, &
- FLHC, &
- FLQC, &
- CHS , &
- EMISS, &
- XICE, &
- XLAND, &
- ALBBCK, &
- VEGFRA, &
- TBOT
- REAL, DIMENSION( 1:nsl), INTENT(IN ) :: ZS
- REAL, DIMENSION( ims:ime , jms:jme ), &
- INTENT(INOUT) :: &
- SNOW, & !new
- SNOWH, &
- SNOWC, &
- CANWAT, & ! new
- SNOALB, &
- ALB, &
- 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, &
- 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, &
- EVAPL, &
- PRCPL, &
- SEAICE, &
- INFILTR
- !--- soil/snow properties
- REAL, DIMENSION( ims:ime, 1:nsl, jms:jme) &
- :: KEEPFR3DFLAG, &
- SMFR3D
- REAL &
- :: RHOCS, &
- RHOSN, &
- RHONEWSN, &
- 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:4001) :: 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, &
- PATM, &
- PATMB, &
- TABS, &
- QVATM, &
- QCATM, &
- Q2SAT, &
- CONFLX, &
- RHO, &
- QKMS, &
- TKMS, &
- INFILTRP
- REAL :: cq,r61,r273,arp,brp,x,evs,eis
- REAL :: meltfactor
- INTEGER :: NROOT
- INTEGER :: ILAND,ISOIL,IFOREST
-
- INTEGER :: I,J,K,NZS,NZS1,NDDZS
- INTEGER :: k1,l,k2,kp,km
- CHARACTER (LEN=132) :: message
- !-----------------------------------------------------------------
- NZS=NSL
- NDDZS=2*(nzs-2)
- !---- 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,4001
- 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
- ! smfr3d (i,k,j)=soilmois(i,k,j)/900.*1.e3
- ! sh2o (i,k,j)=soilmois(i,k,j)-smfr3d(i,k,j)/1.e3*900.
- keepfr3dflag(i,k,j)=0.
- enddo
- !--- initializing to zero snow fraction
- snowc(i,j) = min(1.,snowh(i,j)/0.05)
- !--- initializing inside snow temp if it is not defined
- IF((soilt1(i,j) .LT. 170.) .or. (soilt1(i,j) .GT.400.)) THEN
- IF(snowc(i,j).gt.0.1) THEN
- soilt1(i,j)=0.5*(soilt(i,j)+tso(i,1,j))
- IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) 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) = soilt(i,j)
- ENDIF
- ENDIF
- tsnav(i,j) =0.5*(soilt(i,j)+tso(i,1,j))-273.
- qcg (i,j) =0.
- 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
- qvg (i,j) = QSG(i,j)*mavail(i,j)
- IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
- WRITE ( message , FMT='(A,F8.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
- ! qvg (i,j) =qv3d(i,1,j)
- ! qsfc(i,j) = qsg(i,j)/(1.+qsg(i,j))
- qsfc(i,j) = qvg(i,j)/(1.+qvg(i,j))
- SMELT(i,j) = 0.
- SNOM (i,j) = 0.
- 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.
- emissl (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.
- evapl (i,j) = 0.
- prcpl (i,j) = 0.
- ENDDO
- ENDDO
- do k=1,nsl
- soilice(k)=0.
- soiliqw(k)=0.
- enddo
- endif
- !-----------------------------------------------------------------
- PRCPMS = 0.
- DO J=jts,jte
- DO i=its,ite
- IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) 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)
- !--- 1*e-3 is to convert from mm/s to m/s
- IF(FRPCPN) THEN
- PRCPMS = (RAINBL(i,j)/DT*1.e-3)*(1-FRZFRAC(I,J))
- NEWSNMS = (RAINBL(i,j)/DT*1.e-3)*FRZFRAC(I,J)
- ELSE
- if (tabs.le.273.15) then
- PRCPMS = 0.
- NEWSNMS = RAINBL(i,j)/DT*1.e-3
- else
- PRCPMS = RAINBL(i,j)/DT*1.e-3
- NEWSNMS = 0.
- endif
- ENDIF
- ! if (myj) then
- QKMS=CHS(i,j)
- TKMS=CHS(i,j)
- ! else
- !--- convert exchange coeff QKMS to [m/s]
- ! QKMS=FLQC(I,J)/RHO/MAVAIL(I,J)
- ! TKMS=FLHC(I,J)/RHO/CP
- ! endif
- !--- 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)
- !-----
- 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.195.and.j.eq.254)) &
- 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, &
- NLCAT,ILAND,ISOIL,iswater,MYJ,IFOREST,lufrac, &
- EMISSL(I,J),PC(I,J),ZNT(I,J),LAI(I,J), &
- 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.195.and.j.eq.254)) &
- print *,'after SOILVEGIN - z0,znt(195,254)',z0(i,j),znt(i,j)
- if(ktau.eq.1 .and. (i.eq.195.and.j.eq.254)) 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-4,CMCMAX_DATA * LAI(I,J) * 0.01*VEGFRA(I,J)) ! canopy water saturated
- !-- definition of number of soil levels in the rooting zone
- ! IF(iforest(ivgtyp(i,j)).ne.1) THEN
- IF(iforest.ne.1) 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 1.5
- ! 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((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
- 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) = 1.0
- DO K=1,NZS
- soilmois(i,k,j) = 1.
- smfr3d(i,k,j) = 1.
- sh2o(i,k,j) = 0.
- keepfr3dflag(i,k,j) = 0.
- ENDDO
- ENDIF
- ! 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))
- 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.,soilmois(i,1,j)/(REF-QMIN)))
- ! LMAVAIL(I,J)=max(0.00001,min(1.,soilmois(i,1,j)/dqm))
- 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
- QCG(I,J)=QCG(I,J)-DEW(I,J)
- 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
- !-----------------------------------------------------------------
- CALL SFCTMP (dt,ktau,conflx,i,j, &
- !--- input variables
- nzs,nddzs,nroot,meltfactor, & !added meltfactor
- iland,isoil,xland(i,j),ivgtyp(i,j),PRCPMS, &
- NEWSNMS,SNWE,SNHEI,SNOWFRAC,RHOSN,RHONEWSN, &
- 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), & !new
- myj,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),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),runoff1(I,J), &
- runoff2(I,J),soilice,soiliqw,infiltrp)
- !-----------------------------------------------------------------
- !*** 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)*1000.0
- SMAVAIL (I,J) = SMAVAIL(I,J) * 1000.
- SMMAX (I,J) = SMMAX(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
- qcg(i,j)=qcg(i,j)+dew(i,j)
- 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 surface and PBL scheme
- ! if (myj) then
- ! MYJSFC expects QSFC as actual specific humidity at the surface
- IF((QVATM.GE.Q2SAT*0.95).AND.QVATM.LT.qvg(I,J))THEN
- CHKLOWQ(I,J)=0.
- ELSE
- CHKLOWQ(I,J)=1.
- ENDIF
- ! 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
- ! 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 ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
- print *,' QFX after change, LH ', i,j, QFX(i,j),LH(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
- !--- get 3d soil fields
- IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
- print *,'LAND, i,j,tso1d,soilm1d - end of time step', &
- i,j,tso1d,soilm1d
- ENDIF
- !--- end of a land or sea ice point
- ENDIF
- ENDDO
- ENDDO
- !-----------------------------------------------------------------
- END SUBROUTINE LSMRUC
- !-----------------------------------------------------------------
- SUBROUTINE SFCTMP (delt,ktau,conflx,i,j, &
- !--- input variables
- nzs,nddzs,nroot,meltfactor, &
- ILAND,ISOIL,XLAND,IVGTYP,PRCPMS, &
- NEWSNMS,SNWE,SNHEI,SNOWFRAC,RHOSN,RHONEWSN, &
- PATM,TABS,QVATM,QCATM,rho, &
- GLW,GSW,EMISS,QKMS,TKMS,PC, &
- MAVAIL,CST,VEGFRA,ALB,ZNT, &
- ALB_SNOW,ALB_SNOW_FREE, &
- 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,ACSNOW, &
- edir1,ec1,ett1,eeta,qfx,hfx,s,sublim, &
- evapl,prcpl,runoff1,runoff2,soilice, &
- soiliqw,infiltr)
- !-----------------------------------------------------------------
- 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, &
- SEAICE, &
- XLAND, &
- RHO, &
- QKMS, &
- TKMS
-
- INTEGER, INTENT(IN ) :: IVGTYP
- !--- 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:NDDZS), INTENT(IN) :: DTDZS
- REAL, DIMENSION(1:4001), 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
-
- INTEGER, INTENT(INOUT) :: ILAND,ISOIL
- !-------- 2-d variables
- REAL , &
- INTENT(INOUT) :: DEW, &
- EDIR1, &
- EC1, &
- ETT1, &
- EETA, &
- EVAPL, &
- INFILTR, &
- RHOSN, &
- RHONEWSN, &
- SUBLIM, &
- PRCPL, &
- QVG, &
- QSG, &
- QCG, &
- QFX, &
- HFX, &
- S, &
- RUNOFF1, &
- RUNOFF2, &
- ACSNOW, &
- SNWE, &
- SNHEI, &
- SMELT, &
- SNOM, &
- SNOH, &
- SNFLX, &
- SOILT, &
- SOILT1, &
- TSNAV, &
- ZNT
- REAL, DIMENSION(1:NZS) :: &
- tice, &
- rhosice, &
- capice, &
- thdifice
- !-------- 1-d variables
- REAL, DIMENSION(1:NZS), INTENT(OUT) :: SOILICE, &
- SOILIQW
-
-
- REAL, INTENT(OUT) :: RSM, &
- SNWEPRINT, &
- SNHEIPRINT
- !--- Local variables
- INTEGER :: K,ILNB
- REAL :: BSN, XSN , &
- RAINF, SNTH, NEWSN, PRCPMS, NEWSNMS , &
- T3, UPFLUX, XINET
- REAL :: snhei_crit, keep_snow_albedo
- REAL :: RNET,GSWNEW,EMISSN,ZNTSN
- REAL :: VEGFRAC
- real :: cice, albice, albsn
- !-----------------------------------------------------------------
- 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
- NEWSN=0.
- RAINF = 0.
- RSM=0.
- INFILTR=0.
- VEGFRAC=0.01*VEGFRA
- ! if(VEGFRAC.le.0.01) VEGFRAC=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
- ALBice=ALB_SNOW_FREE
- !--- 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.
- GSWNEW=GSW/(1.-ALB)
- ALBice = MIN(ALB_SNOW_FREE,MAX(ALB_SNOW_FREE - 0.05, &
- ALB_SNOW_FREE - 0.1*(tice(1)+10.)/10. ))
- GSWNEW=GSW*(1.-ALBice)
- endif
- IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
- print *,'I,J,KTAU,QKMS,TKMS', i,j,ktau,qkms,tkms
- print *,'GSW,GSWnew,GLW,SOILT,EMISS,ALB,ALBice,SNWE',&
- GSW,GSWnew,GLW,SOILT,EMISS,ALB,ALBice,SNWE
- ENDIF
- SNHEI = SNWE * 1000. / RHOSN
- !--------------
- T3 = STBOLT*SOILT*SOILT*SOILT
- UPFLUX = T3 *SOILT
- XINET = EMISS*(GLW-UPFLUX)
- RNET = GSWnew + XINET
- !Calculate the amount (m) of fresh snow
- if(snhei.gt.0.0081*1.e3/rhosn) then
- !*** Correct snow density for current temperature (Koren et al. 1999)
- BSN=delt/3600.*c1sn*exp(0.08*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(100.,XSN),400.)
- ! rhosn=MIN(MAX(50.,XSN),400.)
- 777 continue
- ! else
- ! rhosn =200.
- ! rhonewsn =200.
- endif
- newsn=newsnms*delt
- !---- ACSNOW - accumulated snow water [kg m-2]
- acsnow=acsnow+newsn*1.e3
- IF(NEWSN.GT.0.) THEN
- ! IF(NEWSN.GE.1.E-8) THEN
- !*** Calculate fresh snow density (t > -15C, else MIN value)
- !*** Eq. 10 from Koren et al. (1999)
- IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
- print *, 'THERE IS NEW SNOW, newsn', newsn
- ENDIF
- if(tabs.lt.258.15) then
- ! rhonewsn=50.
- rhonewsn=100.
- else
- rhonewsn=1.e3*max((0.10+0.0017*(Tabs-273.15+15.)**1.5) &
- , 0.10)
- ! rhonewsn=1.e3*max((0.05+0.0017*(Tabs-273.15+15.)**1.5) &
- ! , 0.05)
- rhonewsn=MIN(rhonewsn,400.)
- ! rhonewsn=100.
- endif
- !*** 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(100.,XSN),400.)
- ! rhosn=MIN(MAX(50.,XSN),400.)
- snwe=snwe+newsn
- snhei=snwe*1.E3/rhosn
- NEWSN=NEWSN*1.E3/rhonewsn
- ENDIF
- 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
- IF(SNHEI.GT.0.0) THEN
- !--- Set of surface parameters should be changed to snow values for grid
- !--- points where the snow cover exceeds snow threshold of 2 cm
- ILAND=ISICE
- SNHEI_CRIT=0.01601*1.e3/rhosn
- ! SNOWFRAC=MIN(1.,SNHEI/SNHEI_CRIT)
- SNOWFRAC=MIN(1.,SNHEI/(2.*SNHEI_CRIT))
- !tgs - low limit on snow fraction
- ! if(SNOWFRAC.lt.0.01) snowfrac=0.
- !--- EMISS = 0.98 for snow
- EMISS = EMISS*(1.-snowfrac)+0.98*snowfrac
- !-- Set znt for snow from VEGPARM table (snow/ice landuse), except for
- !-- land-use types with higher roughness (forests, urban).
- IF(znt.lt.0.2 .and. snowfrac.gt.0.99) znt=z0tbl(iland)
- KEEP_SNOW_ALBEDO = 0.
- IF (NEWSN.GT.0.) KEEP_SNOW_ALBEDO = 1.
- !--- GSW in-coming solar
- GSWNEW=GSW/(1.-ALB)
- IF(SEAICE .LT. 0.5) THEN
- !----- SNOW on soil
- !-- ALB dependence on snow depth
- ALBsn = MAX(keep_snow_albedo*alb_snow, &
- MIN((alb_snow_free + &
- (alb_snow - alb_snow_free) * snowfrac), alb_snow))
- !28mar11 if canopy is covered with snow to its capacity and snow depth is
- ! higher than patchy snow treshold - then snow albedo is not less than 0.7
- if(cst.ge.sat .and. snowfrac .eq.1.) albsn=max(alb_snow,0.7)
- !-- 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.5) then
- ALB=ALBsn
- 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
- ALBsn = MAX(keep_snow_albedo*alb_snow, &
- MIN((albice + (alb_snow - albice) * snowfrac), alb_snow))
- !-- 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
- !--- recompute absorbed solar radiation and net radiation
- !--- for new value of albedo
- gswnew=gswnew*(1.-alb)
- XINET = EMISS*(GLW-UPFLUX)
- RNET = GSWnew + XINET
- IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
- print *,'I,J,GSW,GSWnew,GLW,UPFLUX,ALB',&
- i,j,GSW,GSWnew,GLW,UPFLUX,ALB
- ENDIF
- if (SEAICE .LT. 0.5) then
- ! LAND
- CALL SNOWSOIL ( & !--- input variables
- i,j,isoil,delt,ktau,conflx,nzs,nddzs,nroot, &
- meltfactor,rhonewsn, & ! new
- ILAND,PRCPMS,RAINF,NEWSN,snhei,SNWE,snowfrac, &
- RHOSN,PATM,QVATM,QCATM, &
- GLW,GSWnew,EMISS,RNET,IVGTYP, &
- QKMS,TKMS,PC,CST, &
- RHO,VEGFRAC,ALB,ZNT, &
- 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,runoff1,runoff2, &
- mavail,soilice,soiliqw,infiltr )
- else
- ! SEA ICE
- CALL SNOWSEAICE ( &
- i,j,isoil,delt,ktau,conflx,nzs,nddzs, &
- meltfactor,rhonewsn, & ! new
- ILAND,PRCPMS,RAINF,NEWSN,snhei,SNWE,snowfrac, &
- RHOSN,PATM,QVATM,QCATM, &
- GLW,GSWnew,EMISS,RNET, &
- QKMS,TKMS,RHO, &
- !--- 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 &
- )
- edir1 = eeta
- 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
- ELSE
- !--- no snow
- snheiprint=0.
- snweprint=0.
- if(SEAICE .LT. 0.5) then
- ! LAND
- CALL SOIL( &
- !--- input variables
- i,j,iland,isoil,delt,ktau,conflx,nzs,nddzs,nroot, &
- PRCPMS,RAINF,PATM,QVATM,QCATM,GLW,GSWnew, &
- EMISS,RNET,QKMS,TKMS,PC,cst,rho,vegfrac, &
- !--- 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,runoff1, &
- runoff2,mavail,soilice,soiliqw, &
- infiltr)
- else
- ! SEA ICE
- 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, &
- !--- 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 &
- )
- edir1 = eeta
- 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
- ! ENDIF
- !
- ! RETURN
- ! END
- !---------------------------------------------------------------
- END SUBROUTINE SFCTMP
- !---------------------------------------------------------------
- FUNCTION QSN(TN,T)
- !****************************************************************
- REAL, DIMENSION(1:4001), 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.4000) GOTO 20
- I=4000
- R=4001.
- 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 ( &
- !--- input variables
- i,j,iland,isoil,delt,ktau,conflx,nzs,nddzs,nroot,&
- PRCPMS,RAINF,PATM,QVATM,QCATM, &
- GLW,GSW,EMISS,RNET, &
- QKMS,TKMS,PC,cst,rho,vegfrac, &
- !--- 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,runoff1,runoff2,mavail,soilice, &
- soiliqw,infiltrp)
- !*************************************************************
- ! 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
- !--- 3-D Atmospheric variables
- REAL, &
- INTENT(IN ) :: PATM, &
- QVATM, &
- QCATM
- !--- 2-D variables
- REAL, &
- INTENT(IN ) :: GLW, &
- GSW, &
- EMISS, &
- RHO, &
- PC, &
- VEGFRAC, &
- 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:4001), 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(INOUT) :: KEEPFR
- !-------- 2-d variables
- REAL, &
- INTENT(INOUT) :: DEW, &
- CST, &
- EDIR1, &
- EC1, &
- ETT1, &
- EETA, &
- EVAPL, &
- PRCPL, &
- MAVAIL, &
- QVG, &
- QSG, &
- QCG, &
- RNET, &
- QFX, &
- HFX, &
- S, &
- SAT, &
- RUNOFF1, &
- RUNOFF2, &
- SOILT
- !-------- 1-d variables
- 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
- REAL, DIMENSION(1:NZS) :: transp,cap,diffu,hydro , &
- thdif,tranf,tav,soilmoism , &
- soilicem,soiliqwm,detal , &
- fwsat,lwsat,told,smold
- REAL :: drip
- 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
- !--- 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( &
- !--- 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)*vegfrac
- ELSE
- DD1=CST+ &
- DELT*(PRCPMS+RAS*FQ*(QVATM-QSG) &
- *(CST/SAT)**CN)*vegfrac
- ENDIF
- 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=(CST/SAT)**CN
- DRYCAN=1.-WETCAN
- !**************************************************************
- ! TRANSF computes transpiration function
- !**************************************************************
- CALL TRANSF( &
- !--- input variables
- nzs,nroot,soiliqw,tabs, &
- !--- soil fixed fields
- dqm,qmin,ref,wilt,zshalf, &
- !--- 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
- !**************************************************************
- ! 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, &
- thdif,cap,drycan,wetcan, &
- transum,dew,mavail, &
- !--- soil fixed fields
- dqm,qmin,bclh,zsmain,zshalf,DTDZS,tbq, &
- !--- constants
- xlv,CP,G0_P,cvw,stbolt, &
- !--- output variables
- tso,soilt,qvg,qsg,qcg)
- !************************************************************************
- !--- 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)
- DO K=1,NZS
- TRANSP(K)=0.
- ENDDO
- ELSE
- DO K=1,NROOT
- TRANSP(K)=VEGFRAC*RAS*QKMS* &
- (QVATM-QSG)* &
- PC*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
- !-- Recalculating of 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,-PRCPMS, &
- QKMS,TRANSP,DRIP,DEW,0.,SOILICE,VEGFRAC, &
- ! QKMS,TRANSP,DRIP,DEW,0.,SOILICE,VEGFRAC, &
- !-- 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*SOILT*SOILT*SOILT
- UPFLUX = T3 *SOILT
- XINET = EMISS*(GLW-UPFLUX)
- RNET = GSW + XINET
- HFT=-TKMS*CP*RHO*(TABS-SOILT) &
- *(P1000mb*0.00001/Patm)**ROVCP
- Q1=-QKMS*RAS*(QVATM - QSG)
- IF (Q1.LE.0.) THEN
- ! --- condensation
- EC1=0.
- EDIR1=0.
- ETT1=0.
- !-- moisture flux for coupling with PBL
- EETA=-QKMS*RAS*(QVATM/(1.+QVATM) - QSG/(1.+QSG))*1.E3
- QFX= XLV*EETA
- !-- actual moisture flux from RUC LSM
- EETA= RHO*DEW
- ELSE
- ! --- evaporation
- EDIR1 =-(1.-vegfrac)*QKMS*RAS* &
- (QVATM-QVG)
- EC1 = Q1 * WETCAN
- CMC2MS=CST/DELT
- if(EC1.gt.CMC2MS) cst=0.
- EC1=MIN(CMC2MS,EC1)*vegfrac
- !-- moisture flux for coupling with PBL
- EETA=-QKMS*RAS*(QVATM/(1.+QVATM) - QVG/(1.+QVG))*1.E3
- ! to convert from kg m-2 s-1 to m s-1: 1/rho water=1.e-3************
- QFX= XLV * EETA
- !-- actual moisture flux from RUC LSM
- EETA = (EDIR1 + EC1 + ETT1)*1.E3
- ENDIF
- EVAPL=EETA
- S=THDIF(1)*CAP(1)*DZSTOP*(TSO(1)-TSO(2))
- HFX=HFT
- FLTOT=RNET-HFT-QFX-S
- 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)
- ! RETURN
- ! END
- !-------------------------------------------------------------------
- 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, &
- !--- 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 &
- )
- !*****************************************************************
- ! 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
- !--- 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:4001), 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
- REAL :: AA1,RHCS
- 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
- 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*SOILT*SOILT*SOILT
- UPFLUX = T3 *SOILT
- XINET = EMISS*(GLW-UPFLUX)
- RNET = GSW + XINET
- HFT=-TKMS*CP*RHO*(TABS-SOILT) &
- *(P1000mb*0.00001/Patm)**ROVCP
- Q1=-QKMS*RAS*(QVATM - QSG)
- IF (Q1.LE.0.) THEN
- ! --- condensation
- !-- moisture flux for coupling with PBL
- EETA=-QKMS*RAS*(QVATM/(1.+QVATM) - QSG/(1.+QSG))*1.E3
- QFX= XLS*EETA
- !-- actual moisture flux from RUC LSM
- DEW=QKMS*(QVATM-QSG)
- EETA= RHO*DEW
- ELSE
- ! --- evaporation
- !-- moisture flux for coupling with PBL
- EETA=-QKMS*RAS*(QVATM/(1.+QVATM) - QVG/(1.+QVG))*1.E3
- ! EETA=-QKMS*RAS*(QVATM/(1.+QVATM) - QSG/(1.+QSG))*1.E3
- ! to convert from kg m-2 s-1 to m s-1: 1/rho water=1.e-3************
- QFX= XLS * EETA
- !-- actual moisture flux from RUC LSM
- EETA = Q1*1.E3
- ENDIF
- EVAPL=EETA
- S=THDIFICE(1)*CAPICE(1)*DZSTOP*(TSO(1)-TSO(2))
- HFX=HFT
- FLTOT=RNET-HFT-QFX-S
- !-------------------------------------------------------------------
- END SUBROUTINE SICE
- !-------------------------------------------------------------------
- SUBROUTINE SNOWSOIL ( &
- !--- input variables
- i,j,isoil,delt,ktau,conflx,nzs,nddzs,nroot, &
- meltfactor,rhonewsn, & ! new
- ILAND,PRCPMS,RAINF,NEWSNOW,snhei,SNWE,SNOWFRAC, &
- RHOSN, &
- PATM,QVATM,QCATM, &
- GLW,GSW,EMISS,RNET,IVGTYP, &
- QKMS,TKMS,PC,cst,rho,vegfrac,alb,znt, &
- 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,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,meltfactor
- LOGICAL, INTENT(IN ) :: myj
- !--- 3-D Atmospheric variables
- REAL, &
- INTENT(IN ) :: PATM, &
- QVATM, &
- QCATM
- !--- 2-D variables
- REAL , &
- INTENT(IN ) :: GLW, &
- GSW, &
- RHO, &
- PC, &
- VEGFRAC, &
- 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:4001), 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(INOUT) :: KEEPFR
- INTEGER, INTENT(INOUT) :: ILAND
- !-------- 2-d variables
- REAL , &
- INTENT(INOUT) :: DEW, &
- CST, &
- 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
- !--- 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 :: drip
- REAL :: RNET
- !-----------------------------------------------------------------
- cvw=cw
- XLMELT=3.35E+5
- !-- the next line calculates heat of sublimation of water vapor
- XLVm=XLV+XLMELT
- ! STBOLT=5.670151E-8
- !--- SNOW flag -- 99
- ILAND=99
- !--- 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.
- DELTSN=0.0301*1.e3/rhosn
- snth=0.01601*1.e3/rhosn
- !tgs when the snow depth is marginlly higher than DELTSN,
- ! reset DELTSN to half of snow depth.
- IF(SNHEI.GE.DELTSN+SNTH) THEN
- ! 2-layer model
- if(snhei-deltsn-snth.lt.snth) deltsn=0.5*(snhei-snth)
- 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( &
- !--- 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
- ! Evaporation
- ! DD1=CST+(NEWSNOW*RHOSN*1.E-3 &
- DD1=CST+(NEWSNOW*RHOnewSN*1.E-3 &
- !-- this change will not let liquid waer be intercepted by the canopy....
- -DELT*(RAS*EPOT &
- ! -DELT*(-PRCPMS+RAS*EPOT &
- *(CST/SAT)**CN)) *vegfrac
- ELSE
- ! Sublimation
- DEW = - EPOT
- ! DD1=CST+(NEWSNOW*RHOSN*1.E-3+delt*( &
- DD1=CST+(NEWSNOW*RHOnewSN*1.E-3+delt*( &
- ! DD1=CST+(NEWSNOW*RHOSN*1.E-3+delt*(PRCPMS &
- +DEW*RAS)) *vegfrac
- ENDIF
- IF(DD1.LT.0.) DD1=0.
- IF (vegfrac.GT.0.) THEN
- CST=DD1
- IF(CST.GT.SAT*vegfrac) THEN
- CST=SAT*vegfrac
- DRIP=DD1-SAT*vegfrac
- 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=SNHEI*RHOSN*1.e-3-vegfrac*NEWSNOW*RHOnewSN*1.E-3 &
- + DRIP
- ENDIF
-
- 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(SNWEPR.LE.EPDT) THEN
- BETA=SNWEPR/max(1.e-8,EPDT)
- SNWE=0.
- SNHEI=0.
- ENDIF
- WETCAN=(CST/SAT)**CN
- DRYCAN=1.-WETCAN
- !**************************************************************
- ! TRANSF computes transpiration function
- !**************************************************************
- CALL TRANSF( &
- !--- input variables
- nzs,nroot,soiliqw,tabs, &
- !--- soil fixed fields
- dqm,qmin,ref,wilt,zshalf, &
- !--- 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,ilnb)
- !************************************************************************
- !--- RECALCULATION OF DEW USING NEW VALUE OF QSG OR TRANSP IF NO DEW
- DEW=0.
- ETT1=0.
- PP=PATM*1.E3
- QSG= QSN(SOILT,TBQ)/PP
- EPOT = -FQ*(QVATM-QSG)
- IF(EPOT.GE.0.) THEN
- ! Evaporation
- DO K=1,NROOT
- TRANSP(K)=vegfrac*RAS*FQ*(QVATM-QSG) &
- *PC*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,-PRCPMS, &
- 0.,TRANSP,0., &
- 0.,SMELT,soilice,vegfrac, &
- !-- 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
- smelt=smelt+snwe/delt
- rsm=0.
- ! snwe=0.
- 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*SOILT*SOILT*SOILT
- UPFLUX = T3 *SOILT
- XINET = EMISS*(GLW-UPFLUX)
- RNET = GSW + XINET
- HFT=-TKMS*CP*RHO*(TABS-SOILT) &
- *(P1000mb*0.00001/Patm)**ROVCP
- Q1 = - FQ*RAS* (QVATM - QSG)
- IF (Q1.LT.0.) THEN
- ! --- condensation
- EDIR1=0.
- EC1=0.
- ETT1=0.
- ! --- condensation
- !-- moisture flux for coupling with PBL
- EETA=-QKMS*RAS*(QVATM/(1.+QVATM) - QSG/(1.+QSG))*1.E3
- QFX= XLVm*EETA
- !-- actual moisture flux from RUC LSM
- DEW=QKMS*(QVATM-QSG)
- EETA= RHO*DEW
- ELSE
- ! --- evaporation
- EDIR1 = Q1*UMVEG *BETA
- EC1 = Q1 * WETCAN
- CMC2MS=CST/DELT
- if(EC1.gt.CMC2MS) cst=0.
- EC1=MIN(CMC2MS,EC1)*vegfrac
- !-- moisture flux for coupling with PBL
- EETA=-QKMS*RAS*(QVATM/(1.+QVATM) - QVG/(1.+QVG))*1.E3
- ! EETA=-QKMS*RAS*(QVATM/(1.+QVATM) - QSG/(1.+QSG))*1.E3
- ! to convert from kg m-2 s-1 to m s-1: 1/rho water=1.e-3************
- QFX= XLVm * EETA
- !-- actual moisture flux from RUC LSM
- EETA = (EDIR1 + EC1 + ETT1)*1.E3
- ENDIF
- s=THDIF(1)*CAP(1)*dzstop*(tso(1)-tso(2))
- HFX=HFT
- FLTOT=RNET-HFT-QFX-S
- 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)
- ! RETURN
- ! END
- !-------------------------------------------------------------------
- END SUBROUTINE SNOWSOIL
- !-------------------------------------------------------------------
- SUBROUTINE SNOWSEAICE( &
- i,j,isoil,delt,ktau,conflx,nzs,nddzs, &
- meltfactor,rhonewsn, & ! new
- ILAND,PRCPMS,RAINF,NEWSNOW,snhei,SNWE,snowfrac, &
- RHOSN,PATM,QVATM,QCATM, &
- GLW,GSW,EMISS,RNET, &
- QKMS,TKMS,RHO, &
- !--- 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 &
- )
- !***************************************************************
- ! 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
- real :: rhonewcsn
- !--- 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:4001), 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
- integer :: nmelt
- !-----------------------------------------------------------------
- XLMELT=3.35E+5
- !-- the next line calculates heat of sublimation of water vapor
- XLVm=XLV+XLMELT
- ! STBOLT=5.670151E-8
- !--- SNOW flag -- 99
- ILAND=99
- !--- 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.
- DELTSN=0.0301*1.e3/rhosn
- snth=0.01601*1.e3/rhosn
- !tgs when the snow depth is marginlly higher than DELTSN,
- ! reset DELTSN to half of snow depth.
- IF(SNHEI.GE.DELTSN+SNTH) THEN
- ! 2-layer model
- if(snhei-deltsn-snth.lt.snth) deltsn=0.5*(snhei-snth)
- 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.
- hsn=snhei
- ! cvw=cw
- NZS1=NZS-1
- NZS2=NZS-2
- QGOLD=QVG
- 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(SNWEPR.LE.EPDT) THEN
- BETA=SNWEPR/max(1.e-8,EPDT)
- SNWE=0.
- SNHEI=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=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.
- fsn=SNHEI/(SNHEI+zsmain(2))
- fso=1.-fsn
- soilt1=tso(1)
- tsob=tso(2)
- snprim=SNHEI+zsmain(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
- 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
- IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
- print *,' SNHEI = 0, D9SN,R22SN,D1SN,D2SN: ',D9SN,R22SN,D1SN,D2SN
- ENDIF
- ENDIF
- !---- TDENOM for snow
- !18apr08 - the iteration start point
- 212 continue
- 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) &
- !18apr08 - add heat of snow phase change
- -SNOH &
- )/TDENOM
- AA1=AA
- PP=PATM*1.E3
- AA1=AA1/PP
- 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
- 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=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
- ELSE
- TSO(1)=SOILT
- SOILT1=SOILT
- tsob=SOILT
- ENDIF
- !---- Final solution for TSO in sea ice
- DO K=2,NZS
- KK=NZS-K+1
- TSO(K)=min(271.4,(rhtso(KK)+cotso(KK)*TSO(K-1)))
- END DO
- !--- For thin snow layer combined with the top sea ice layer
- !--- TSO is computed by linear inmterpolation between SOILT
- !--- and TSO(2)
- 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
- nmelt = 1
- soiltfrac=snowfrac*273.15+(1.-snowfrac)*SOILT
- QSG= QSN(soiltfrac,TBQ)/PP
- QVG=QSG
- T3 = STBOLT*SOILTfrac*SOILTfrac*SOILTfrac
- UPFLUX = T3 * 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)-TN) &
- -SOH-X+RAINF*CVW*PRCPMS* &
- (max(273.15,TABS)-TN)
- 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)
- !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)))
- !*** From Koren et al. (1999) 13% of snow melt stays in the snow pack
- !!! rsm=0.13*smelt*delt
- if(snwepr.gt.0.) then
- rsmfrac=min(0.18,(max(0.08,0.10/snwepr*0.13)))
- ! else
- ! rsmfrac=0.13
- endif
- rsm=rsmfrac*smelt*delt
- !18apr08 rsm is part of melted water that stays in snow as liquid
- SMELT=AMAX1(0.,SMELT-rsm/delt)
- SNOHGNEW=SMELT*XLMELT*1.E3
- SNODIF=AMAX1(0.,(SNOH-SNOHGNEW))
- SNOH=SNOHGNEW
- !18apr08 - if snow melt occurred then go into iteration for energy budget
- ! solution
- !-- correction of liquid equivalent of snow depth
- !-- due to evaporation and snow melt
- SNWE = AMAX1(0.,(SNWEPR- &
- (SMELT+BETA*EPOT*RAS)*DELT &
- ) )
- !--- If all snow melts, then 13% of snow melt we kept in the
- !--- snow pack should be added back to snow melt and infiltrate
- !--- into soil.
- if(snwe.le.rsm) then
- smelt=smelt+rsm/delt
- snwe=0.
- rsm=0.
- else
- !*** Correct snow density on effect of snow melt, melted
- !*** from the top of the snow. 13% of melted water
- !*** remains in the pack and changes its density.
- !*** Eq. 9 (with my correction) in Koren et al. (1999)
- if(snwe.gt.0.) then
- xsn=(rhosn*(snwe-rsm)+1.e3*rsm)/ &
- snwe
- rhosn=MIN(XSN,400.)
- RHOCSN=2090.* RHOSN
- thdifsn = 0.265/RHOCSN
- endif
- endif
- !--- If there is no snow melting then just evaporation
- !--- or condensation cxhanges SNWE
- ELSE
- EPOT=-QKMS*(QVATM-QSG)
- SNWE = AMAX1(0.,(SNWEPR- &
- BETA*EPOT*RAS*DELT))
- ENDIF
- !*** Correct snow density on effect of snow melt, melted
- !*** from the top of the snow. 13% of melted water
- !*** remains in the pack and changes its density.
- !*** Eq. 9 (with my correction) in Koren et al. (1999)
- SNHEI=SNWE *1.E3 / RHOSN
- snweprint=snwe
- ! &
- !--- if VEGFRAC.ne.0. then some snow stays on the canopy
- !--- and should be added to SNWE for water conservation
- ! 4 Nov 07 +VEGFRAC*cst
- snheiprint=snweprint*1.E3 / RHOSN
- if(nmelt.eq.1) goto 212
- 220 continue
- IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
- print *, 'snweprint : ',snweprint
- print *, 'D9SN,SOILT,TSOB : ', D9SN,SOILT,TSOB
- ENDIF
- !--- Compute flux in the top snow layer
- SNFLX=D9SN*(SOILT-TSOB)
- 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
- ENDIF
- !--- RECALCULATION OF DEW USING NEW VALUE OF QSG
- DEW=0.
- PP=PATM*1.E3
- QSG= QSN(SOILT,TBQ)/PP
- EPOT = -FQ*(QVATM-QSG)
- IF(EPOT.LT.0.) THEN
- ! Sublimation
- DEW=-EPOT
- ENDIF
- !-- Restore sea-ice parameters if snow is less than threshold
- IF(SNHEI.EQ.0.) then
- tsnav=soilt-273.15
- smelt=smelt+snwe/delt
- rsm=0.
- emiss=1.
- znt=0.011
- alb=0.55
- ENDIF
- SNOM=SNOM+SMELT*DELT*1.e3
- !--- THE DIAGNOSTICS OF SURFACE FLUXES
- T3 = STBOLT*SOILT*SOILT*SOILT
- UPFLUX = T3 *SOILT
- XINET = EMISS*(GLW-UPFLUX)
- RNET = GSW + XINET
- HFT=-TKMS*CP*RHO*(TABS-SOILT) &
- *(P1000mb*0.00001/Patm)**ROVCP
- Q1 = - FQ*RAS* (QVATM - QSG)
- IF (Q1.LT.0.) THEN
- ! --- condensation
- !-- moisture flux for coupling with PBL
- EETA=-QKMS*RAS*(QVATM/(1.+QVATM) - QSG/(1.+QSG))*1.E3
- QFX= XLVm*EETA
- !-- actual moisture flux from RUC LSM
- DEW=QKMS*(QVATM-QSG)
- EETA= RHO*DEW
- sublim = EETA
- ELSE
- ! --- evaporation
- !-- moisture flux for coupling with PBL
- EETA=-QKMS*RAS*(QVATM/(1.+QVATM) - QVG/(1.+QVG))*1.E3
- ! EETA=-QKMS*RAS*(QVATM/(1.+QVATM) - QSG/(1.+QSG))*1.E3
- ! to convert from kg m-2 s-1 to m s-1: 1/rho water=1.e-3************
- QFX= XLVm * EETA
- !-- actual moisture flux from RUC LSM
- EETA = Q1*1.E3
- sublim = EETA
- ENDIF
- s=THDIFICE(1)*CAPICE(1)*dzstop*(tso(1)-tso(2))
- ! s=D9SN*(SOILT-TSOB)
- HFX=HFT
- FLTOT=RNET-HFT-QFX-S
- !------------------------------------------------------------------------
- !------------------------------------------------------------------------
- END SUBROUTINE SNOWSEAICE
- !------------------------------------------------------------------------
- SUBROUTINE 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, &
- THDIF,CAP,DRYCAN,WETCAN, &
- TRANSUM,DEW,MAVAIL, &
- !--- soil fixed fields
- DQM,QMIN,BCLH, &
- ZSMAIN,ZSHALF,DTDZS,TBQ, &
- !--- constants
- XLV,CP,G0_P,CVW,STBOLT, &
- !--- output variables
- TSO,SOILT,QVG,QSG,QCG)
- !*************************************************************
- ! Energy budget equation and heat diffusion eqn are
- ! solved here and
- !
- ! 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 (m^2/s)
- ! QSG,QVG,QCG - saturated mixing ratio, mixing ratio of
- ! water vapor and cloud at the ground
- ! surface, respectively (kg/kg)
- ! PATM - pressure [bar]
- ! QC3D,QV3D - 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^2s
- ! 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
- REAL, INTENT(INOUT) :: DRYCAN,WETCAN,TRANSUM
- !--- 3-D Atmospheric variables
- REAL, &
- INTENT(IN ) :: PATM, &
- QVATM, &
- QCATM
- !--- 2-D variables
- REAL , &
- INTENT(IN ) :: &
- EMISS, &
- RHO, &
- RNET, &
- PC, &
- VEGFRAC, &
- DEW, &
- QKMS, &
- TKMS
- !--- soil properties
- REAL , &
- INTENT(IN ) :: &
- BCLH, &
- DQM, &
- QMIN
- REAL, INTENT(IN ) :: CP, &
- CVW, &
- XLV, &
- STBOLT, &
- TABS, &
- G0_P
- REAL, DIMENSION(1:NZS), INTENT(IN) :: ZSMAIN, &
- ZSHALF, &
- THDIF, &
- CAP
- REAL, DIMENSION(1:NDDZS), INTENT(IN) :: DTDZS
- REAL, DIMENSION(1:4001), 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) :: &
- MAVAIL, &
- QVG, &
- QSG, &
- QCG, &
- SOILT
- !--- Local variables
- REAL :: x,x1,x2,x4,dzstop,can,ft,sph , &
- tn,trans,umveg,denom
- 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
- REAL :: C,CC,AA1,RHCS,H1
- REAL, DIMENSION(1:NZS) :: cotso,rhtso
- INTEGER :: nzs1,nzs2,k,k1,kn,kk
- !-----------------------------------------------------------------
- NZS1=NZS-1
- NZS2=NZS-2
- dzstop=1./(ZSMAIN(2)-ZSMAIN(1))
- do k=1,nzs
- cotso(k)=0.
- rhtso(k)=0.
- enddo
- !******************************************************************************
- ! 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 HEAT BALANCE EQUATION (Smirnova et al., 1996, EQ. 21,26)
- RHCS=CAP(1)
- H=MAVAIL
- IF(DEW.NE.0.)THEN
- DRYCAN=0.
- WETCAN=1.
- ENDIf
- TRANS=PC*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
- TDENOM=D9*(1.-D1+R22)+D10+R21+R7 &
- +RAINF*CVW*PRCPMS
- FKQ=QKMS*RHO
- R210=R211*RHO
- C=VEGFRAC*FKQ*CAN
- CC=C*XLV/TDENOM
- AA=XLV*(FKQ*UMVEG+R210)/TDENOM
- BB=(D10*TABS+R21*TN+XLV*(QVATM* &
- (FKQ*UMVEG+C) &
- +R210*QVG)+D11+D9*(D2+R22*TN) &
- +RAINF*CVW*PRCPMS*max(273.15,TABS) &
- )/TDENOM
- AA1=AA+CC
- PP=PATM*1.E3
- AA1=AA1/PP
- IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
- PRINT *,' VILKA-1'
- print *,'D10,TABS,R21,TN,QVATM,FKQ,UMVEG,VEGFRAC,CAN', &
- D10,TABS,R21,TN,QVATM,FKQ,UMVEG,VEGFRAC,CAN
- 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,umveg,fkq,r210,vegfrac', &
- tn,aa1,bb,pp,umveg,fkq,r210,vegfrac
- ENDIF
- 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(Q1.LT.QS1) GOTO 100
- !--- if no saturation - goto 100
- !--- if saturation - goto 90
- 90 QVG=QS1
- QSG=QS1
- TSO(1)=TS1
- QCG=max(0.,Q1-QS1)
- GOTO 200
- 100 BB=BB-AA*TX2
- AA=(AA*H+CC)/PP
- IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
- PRINT *,' VILKA-2'
- print *,'D10,TABS,R21,TN,QVATM,FKQ,UMVEG,VEGFRAC,CAN', &
- D10,TABS,R21,TN,QVATM,FKQ,UMVEG,VEGFRAC,CAN
- 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,umveg,fkq,r210,vegfrac', &
- tn,aa1,bb,pp,umveg,fkq,r210,vegfrac
- ENDIF
- CALL VILKA(TN,AA,BB,PP,QS1,TS1,TBQ,KTAU,i,j,iland,isoil)
- Q1=TX2+H*QS1
- IF(Q1.GT.QS1) GOTO 90
- QSG=QS1
- QVG=Q1
- TSO(1)=TS1
- QCG=0.
- 200 CONTINUE
- !--- 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
- !--------------------------------------------------------------------
- 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,ILNB)
- !********************************************************************
- ! 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:4001), 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, &
- 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, &
- s, &
- qfx, &
- hfx
- REAL :: RNET,rsmfrac,soiltfrac,hsn
- integer :: nmelt
- !-----------------------------------------------------------------
- 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
- TNOLD=SOILT
- 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
- ilnb=1
- snprim=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
- ilnb=2
- snprim=deltsn
- tsob=soilt1
- XSN = DELT/2./(0.5*SNPRIM)
- 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
- !--- snow is too thin to be treated separately, therefore it
- !--- is combined with the first soil layer.
- fsn=SNHEI/(SNHEI+zsmain(2))
- fso=1.-fsn
- soilt1=tso(1)
- tsob=tso(2)
- snprim=SNHEI+zsmain(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
- 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-QSG)
- RHCS=CAP(1)
- H=MAVAIL
- IF(DEW.NE.0.)THEN
- DRYCAN=0.
- WETCAN=1.
- ENDIF
- TRANS=PC*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)
- 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 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)
- 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
- !---- TDENOM for snow
- !18apr08 - the iteration start point
- 212 continue
- 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*QVG)+D11+D9SN*(D2SN+R22SN*TN) &
- +RAINF*CVW*PRCPMS*max(273.15,TABS) &
- + RHOnewCSN*NEWSNOW/DELT*min(273.15,TABS) &
- !18apr08 - added heat of snow phase change computed in the first iteration
- -SNOH &
- )/TDENOM
- AA1=AA+CC
- PP=PATM*1.E3
- AA1=AA1/PP
- IF ( wrf_at_debug_level(LSMRUC_DBG_LVL) ) THEN
- print *,'VILKA-SNOW'
- print *,'tn,aa1,bb,pp,umveg,fkq,r210,vegfrac', &
- tn,aa1,bb,pp,umveg,fkq,r210,vegfrac
- 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'
- print *,' TS1,QS1: ', ts1,qs1
- 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.,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
- ELSE
- TSO(1)=SOILT
- SOILT1=SOILT
- tsob=SOILT
- ENDIF
- !---- Final solution for TSO
- DO K=2,NZS
- KK=NZS-K+1
- TSO(K)=rhtso(KK)+cotso(KK)*TSO(K-1)
- END DO
- !--- For thin snow layer combined with the top soil layer
- !--- TSO is computed by linear inmterpolation between SOILT
- !--- and TSO(2)
- if(SNHEI.LT.SNTH.AND.SNHEI.GT.0.)then
- tso(1)=tso(2)+(soilt-tso(2))*fso
- SOILT1=TSO(1)
- tsob=tso(2)
- !!! tsob=tso(1)
- 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
- nmelt = 1
- soiltfrac=snowfrac*273.15+(1.-snowfrac)*SOILT
- QSG= QSN(soiltfrac,TBQ)/PP
- QVG=QSG
- T3 = STBOLT*SOILTfrac*SOILTfrac*SOILTfrac
- UPFLUX = T3 * SOILTfrac
- XINET = EMISS*(GLW-UPFLUX)
- RNET = GSW + XINET
- EPOT = -QKMS*(QVATM-QSG)
- Q1=EPOT*RAS
- IF (Q1.LE.0.) 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 &
- *PC*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
- 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-TNOLD) + &
- XLVM*R210*(QSG-QGOLD)
- !-- SNOH is energy flux of snow phase change
- SNOH=RNET+QFX +HFX &
- +RHOnewCSN*NEWSNOW/DELT*(min(273.15,TABS)-TN) &
- -SOH-X+RAINF*CVW*PRCPMS* &
- (max(273.15,TABS)-TN)
- 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*UMVEG)
- SMELT=AMIN1(SMELT,SNWEPR/DELT-BETA*EPOT*RAS)
- 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)))
- !*** From Koren et al. (1999) 13% of snow melt stays in the snow pack
- !!! rsm=0.13*smelt*delt
- if(snwepr.gt.0.) then
- rsmfrac=min(0.18,(max(0.08,0.10/snwepr*0.13)))
- ! else
- ! rsmfrac=0.13
- endif
- rsm=rsmfrac*smelt*delt
- !18apr08 rsm is part of melted water that stays in snow as liquid
- SMELT=AMAX1(0.,SMELT-rsm/delt)
- SNOHGNEW=SMELT*XLMELT*1.E3
- SNODIF=AMAX1(0.,(SNOH-SNOHGNEW))
- SNOH=SNOHGNEW
- !-- correction 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*UMVEG)*DELT &
- ) )
- if(snwe.le.rsm) then
- smelt=smelt+rsm/delt
- snwe=0.
- rsm=0.
- else
- !*** Correct snow density on effect of snow melt, melted
- !*** from the top of the snow. 13% of melted water
- !*** remains in the pack and changes its density.
- !*** Eq. 9 (with my correction) in Koren et al. (1999)
- if(snwe.gt.0.) then
- xsn=(rhosn*(snwe-rsm)+917.*rsm)/ &
- snwe
- rhosn=MIN(XSN,400.)
- RHOCSN=2090.* RHOSN
- thdifsn = 0.265/RHOCSN
- endif
- endif
- !--- If there is no snow melting then just evaporation
- !--- or condensation cxhanges SNWE
- ELSE
- EPOT=-QKMS*(QVATM-QSG)
- SNWE = AMAX1(0.,(SNWEPR- &
- BETA*EPOT*RAS*DELT))
- ! BETA*EPOT*RAS*UMVEG*DELT))
- ENDIF
- !*** Correct snow density on effect of snow melt, melted
- !*** from the top of the snow. 13% of melted water
- !*** remains in the pack and changes its density.
- !*** Eq. 9 (with my correction) in Koren et al. (1999)
- SNHEI=SNWE *1.E3 / RHOSN
- !18apr08 - if snow melt occurred then go into iteration for energy budget
- ! solution
- if(nmelt.eq.1) goto 212
- 220 continue
- !-- Snow melt from the top is done. But if ground surface temperature
- !-- is above freezing snow can melt from the bottom. The following
- !-- piece of code will check if bottom melting is possible.
- IF(TSO(1).GT.273.15.AND.SNHEI.GT.0.) THEN
- if (snhei.GE.deltsn+snth) then
- hsn = snhei - deltsn
- else
- hsn = snhei
- endif
- soiltfrac=snowfrac*273.15+(1.-snowfrac)*TSO(1)
- SNOHG=(TSO(1)-soiltfrac)*(RHCS*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)
- if(SNWE-SMELTG*DELT.ge.rsm) then
- SNWE = AMAX1(0.,SNWE-SMELTG*DELT)
- else
- smeltg=snwe/delt
- snwe=0.
- rsm=0.
- hsn=0.
- endif
- SNOHGNEW=SMELTG*XLMELT*1.e3
- SNODIF=AMAX1(0.,(SNOHG-SNOHGNEW))
- TSO(1)=soiltfrac &
- + SNODIF/(RHCS*zshalf(2)+ RHOCSN*0.5*hsn)* DELT
- SMELT=SMELT+SMELTG
- SNOH=SNOH+SNOHGNEW
- ENDIF
- SNHEI=SNWE *1.E3 / RHOSN
- snweprint=snwe
- ! &
- !--- if VEGFRAC.ne.0. then some snow stays on the canopy
- !--- and should be added to SNWE for water conservation
- ! 4 Nov 07 +cst
- 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
- !--- Compute flux in the top snow layer
- SNFLX=D9SN*(SOILT-TSOB)
- 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
- ENDIF
- ! return
- ! end
- !------------------------------------------------------------------------
- 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, &
- !--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 , &
- DQM,QMIN,REF,KSAT,RAS,RIW
-
- ! 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)=(SOILIQW(NZS)+TRANSP(NZS)*DELT/ &
- DID)/DENOM
- 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
- 330 RHSMC(K+1)=(SOILIQW(KN)+Q2*RHSMC(K) &
- +TRANSP(KN) &
- /(ZSHALF(KN+1)-ZSHALF(KN)) &
- *DELT)/DENOM
- ! --- MOISTURE BALANCE BEGINS HERE
- TRANS=TRANSP(1)
- UMVEG=1.-VEGFRAC
- 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
- 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
- INFMAX1 = (PX*(DDT/(PX+DDT)))/DELT
- ! print *,'INFMAX1=,ksat',infmax1,ksat,f1,f2,kdt,val,ddt,px
- ! ----------- 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
- ! print *,'FCR--------',fcr
- INFMAX1 = INFMAX1* FCR
- ! -------------------------------------------------------------------
- INFMAX = MAX(INFMAX1,HYDRO(1)*SOILIQW(1))
- INFMAX = MIN(INFMAX, -TOTLIQ)
-
- !----
- IF (-TOTLIQ.GT.INFMAX)THEN
- RUNOFF=-TOTLIQ-INFMAX
- FLX=-INFMAX
- ENDIF
- ! INFILTRP is total infiltration flux in M/S
- INFILTRP=FLX
- ! Solution of moisture budget
- R7=.5*DZS/DELT
- R4=R4+R7
- FLX=FLX-SOILIQW(1)*R7
- R8=UMVEG*R6
- 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
- SOILIQW (1)=1.e-8
- SOILMOIS(1)=1.e-8+soilice(1)*riw
- ELSE IF(QQ.GT.DQM) THEN
- !-- saturation
- SOILIQW (1)=DQM
- SOILMOIS(1)=DQM
- RUNOFF2=(FLXSAT-FLX)*DELT
- RUNOFF=RUNOFF+(FLXSAT-FLX)
- ELSE
- SOILIQW (1)=min(dqm,max(1.e-8,QQ))
- SOILMOIS(1)=max(1.e-8,QQ)+soilice(1)*riw
- END IF
- !--- FINAL SOLUTION FOR SOILMOIS
- DO K=2,NZS
- KK=NZS-K+1
- QQ=COSMC(KK)*SOILIQW(K-1)+RHSMC(KK)
- IF (QQ.LT.0.) THEN
- SOILIQW(K) =1.e-8
- SOILMOIS(K)=1.e-8 + soilice(k)*riw
- ELSE IF(QQ.GT.DQM) THEN
- !-- saturation
- SOILIQW (K)=DQM
- SOILMOIS(K)=DQM
- IF(K.EQ.NZS)THEN
- RUNOFF2=RUNOFF2+(QQ-DQM)*(ZSMAIN(K)-ZSHALF(K))
- ELSE
- RUNOFF2=RUNOFF2+(QQ-DQM)*(ZSHALF(K+1)-ZSHALF(K))
- ENDIF
- ELSE
- SOILIQW (K)=min(dqm,max(1.e-8,QQ))
- SOILMOIS(K)=max(1.e-8,QQ)+soilice(k)*riw
- END IF
- END DO
- ! MAVAIL=min(1.,SOILMOIS(1)/(REF-QMIN))
- MAVAIL=min(1.,SOILMOIS(1)/DQM)
- if (MAVAIL.EQ.0.) MAVAIL=.00001
- ! RETURN
- ! END
- !-------------------------------------------------------------------
- END SUBROUTINE SOILMOIST
- !-------------------------------------------------------------------
- SUBROUTINE SOILPROP( &
- !--- 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
- !--- 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
- 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)=KSAT/am* &
- (soiliqw(K)/am) &
- **(2.*BCLH+2.) &
- * fach
- endif
- ENDDO
- ! RETURN
- ! END
- !-----------------------------------------------------------------------
- END SUBROUTINE SOILPROP
- !-----------------------------------------------------------------------
- SUBROUTINE TRANSF( &
- !--- input variables
- nzs,nroot,soiliqw,tabs, &
- !--- soil fixed fields
- dqm,qmin,ref,wilt,zshalf, &
- !--- 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 ) :: nroot,nzs
- REAL , &
- INTENT(IN ) :: TABS
- !--- soil properties
- REAL , &
- INTENT(IN ) :: DQM, &
- QMIN, &
- REF, &
- 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
- REAL, DIMENSION(1:NZS) :: PART
- !--------------------------------------------------------------------
- do k=1,nzs
- part(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
- !--- 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(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
- !cc TRANF(1)=part(1)
- TRANF(1)=TRANF(1)*FTEM
- 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
- !cc TRANF(k)=part(k)
- TRANF(k)=TRANF(k)*FTEM
- END DO
- !-- TRANSUM - total for the rooting zone
- transum=0.
- DO K=1,NROOT
- transum=transum+tranf(k)
- END DO
- ! RETURN
- ! END
- !-----------------------------------------------------------------
- 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:4001), 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.4000.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.4000.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 *,' AVOST IN VILKA '
- ! WRITE(12,*)'AVOST',TN,D1,D2,PP,NSTEP
- PRINT *,TN,D1,D2,PP,NSTEP,I,TT(i),ii,j,iland,isoil
- CALL wrf_error_fatal (' AVOST IN VILKA ' )
- 20 CONTINUE
- ! RETURN
- ! END
- !-----------------------------------------------------------------------
- END SUBROUTINE VILKA
- !-----------------------------------------------------------------------
- SUBROUTINE SOILVEGIN ( mosaic_lu,mosaic_soil,soilfrac,nscat, &
- NLCAT,IVGTYP,ISLTYP,iswater,MYJ, &
- IFOREST,lufrac,EMISS,PC,ZNT,LAI,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
- !---- 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
- LOGICAL, INTENT(IN ) :: myj
- REAL, DIMENSION( 1:NLCAT ), INTENT(IN):: LUFRAC
- REAL, DIMENSION( 1:NSCAT ), INTENT(IN):: SOILFRAC
- REAL , &
- INTENT ( OUT) :: pc
- REAL , &
- INTENT (INOUT ) :: emiss, &
- lai, &
- znt
- !--- 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
- !***********************************************************************
- ! 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)
- EMISS = 0.
- ZNT = 0.
- PC = 0.
- LAI = 0.
- AREA = 0.
- !-- mosaic approach to landuse in the grid box
- if(mosaic_lu == 1) then
- do k = 1,nlcat
- AREA = AREA + lufrac(k)
- EMISS = EMISS+ LEMITBL(K)*lufrac(k)
- ZNT = ZNT + Z0TBL(K)*lufrac(k)
- LAI = LAI + LAITBL(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(i.eq.222.and.j.eq.335) then
- ! print *,'area=',area,i,j,ivgtyp,nlcat,(lufrac(k),k=1,nlcat),EMISS,ZNT,LAI,PC
- ! endif
- EMISS = EMISS/AREA
- ZNT = ZNT/AREA
- LAI = LAI/AREA
- PC = PC /AREA
- ! EMISS = LEMI(IVGTYP)
- else
- EMISS = LEMITBL(IVGTYP)
- ZNT = Z0TBL(IVGTYP)
- PC = PCTBL(IVGTYP)
- LAI = LAITBL(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
- !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
- 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
- ! 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 )
- #ifdef WRF_CHEM
- 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') then
- MMINLURUC='MODI-RUC'
- endif
- MMINSL='STAS-RUC'
- ! CALL RUCLSM_PARM_INIT
- print *,'RUCLSMINIT uses ',mminluruc
- call RUCLSM_SOILVEGPARM( MMINLURUC, MMINSL)
- ENDIF
- #ifdef WRF_CHEM
- !
- ! 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)/(dqm+qmin)))
- ! mavail(i,j) = max(0.00001,min(1.,smois(i,1,j)/(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
- 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 )
- !
- !-----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