MODULE module_sf_pxlsm USE module_model_constants USE module_sf_pxlsm_data INTEGER, PARAMETER :: NSOLD=20 REAL, PARAMETER :: RD = 287.04, CPD = 1004.67, & CPH2O = 4.218E+3, CPICE = 2.106E+3, & LSUBF = 3.335E+5, SIGMA = 5.67E-8, & ROVCP = RD / CPD REAL, PARAMETER :: CRANKP = 0.5 REAL, PARAMETER :: RIC = 0.25 REAL, PARAMETER :: DENW = 1000.0 REAL, PARAMETER :: TAUINV = 1.0 / 86400.0 REAL, PARAMETER :: T2TFAC = 1.0 / 10.0 REAL, PARAMETER :: PI = 3.1415926 REAL, PARAMETER :: PR0 = 0.95 REAL, PARAMETER :: CZO = 0.032 REAL, PARAMETER :: OZO = 1.E-4 CONTAINS SUBROUTINE pxlsm(U3D, V3D, DZ8W, QV3D,T3D,TH3D, RHO, & PSFC, GSW, GLW, RAINBL, EMISS, & ITIMESTEP,CURR_SECS,NSOIL,DT,ANAL_INTERVAL, & XLAND, XICE, ALBBCK, ALBEDO, & SNOALB, SMOIS, TSLB, MAVAIL, TA2, & QA2, QSFC, ZS,DZS, PSIH, & LANDUSEF,SOILCTOP,SOILCBOT,VEGFRA,VEGF_PX, & ISLTYP,RA,RS,LAI,IMPERV,CANFRA,NLCAT,NSCAT, & HFX,QFX,LH,TSK,SST,ZNT,CANWAT, & GRDFLX,SHDMIN,SHDMAX, & SNOWC,PBLH,RMOL,UST,CAPG,DTBL, & T2_NDG_OLD, T2_NDG_NEW, & Q2_NDG_OLD, Q2_NDG_NEW, & SN_NDG_OLD, SN_NDG_NEW, SNOW, SNOWH,SNOWNCV,& T2OBS, Q2OBS, PXLSM_SMOIS_INIT, & PXLSM_SOIL_NUDGE, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) IMPLICIT NONE INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte INTEGER, INTENT(IN) :: NSOIL, ITIMESTEP, NLCAT, NSCAT, & ANAL_INTERVAL, PXLSM_SMOIS_INIT, PXLSM_SOIL_NUDGE REAL, INTENT(IN ),OPTIONAL :: curr_secs REAL, INTENT(IN ) :: DT,DTBL INTEGER, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: ISLTYP REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) :: U3D, V3D, RHO, & T3D, TH3D, DZ8W, QV3D REAL, DIMENSION(1:NSOIL), INTENT(IN)::ZS,DZS REAL, DIMENSION( ims:ime , 1:NSOIL, jms:jme ), INTENT(INOUT) :: SMOIS, TSLB REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT):: RA, RS, LAI, ZNT, QSFC REAL, DIMENSION( ims:ime , jms:jme ), INTENT(OUT):: GRDFLX, TSK, TA2, QA2 REAL, DIMENSION( ims:ime , 1:NLCAT, jms:jme ), INTENT(IN):: LANDUSEF REAL, DIMENSION( ims:ime , 1:NSCAT, jms:jme ), INTENT(IN):: SOILCTOP, SOILCBOT REAL, DIMENSION( ims:ime, jms:jme ), & INTENT(IN) :: PSFC, GSW, GLW, RAINBL, & ALBBCK, SHDMIN, SHDMAX, & PBLH, RMOL, SNOWNCV, & UST, MAVAIL, SST, EMISS REAL, DIMENSION( ims:ime, jms:jme ), & INTENT(IN) :: T2_NDG_OLD, T2_NDG_NEW, & Q2_NDG_OLD, Q2_NDG_NEW, & SN_NDG_OLD, SN_NDG_NEW REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: T2OBS, Q2OBS REAL, DIMENSION( ims:ime, jms:jme ), & INTENT(INOUT) :: CAPG,CANWAT, QFX, HFX, LH, & PSIH,VEGFRA, VEGF_PX, SNOW, SNOALB, & SNOWH, SNOWC, ALBEDO, XLAND, XICE, & IMPERV, CANFRA LOGICAL :: radiation INTEGER, PARAMETER :: NSTPS = 11 REAL, PARAMETER :: DTPBLX = 40.0 INTEGER, DIMENSION( 1: NSTPS ) :: JP INTEGER:: J, I, NS, NUDGE, ISTI, WEIGHT INTEGER:: NTSPS, IT REAL, DIMENSION( ims:ime, jms:jme ) :: XLAI, XLAIMN, RSTMIN, & XVEG, XVEGMN, XSNUP, & XALB, XSNOALB, WETFRA REAL, DIMENSION( ims:ime, jms:jme ) :: RADNET, EG, ER, ETR, QST REAL:: SFCPRS,TA1,DENS1,QV1,ZLVL,SOLDN,LWDN, & EMISSI,PRECIP,THETA1,VAPPRS,QSBT, & WG,W2,WR,TG,T2,USTAR,MOLX,Z0, & RAIR,CPAIR,IFLAND,ISNOW, & ES,QSS,BETAP, & RH2_OLD, RH2_NEW, T2_OLD, T2_NEW, & CORE, CORB, TIME_BETWEEN_ANALYSIS, & G1000, ALN10,RH2OBS, HU, SNOBS, & FWSAT,FWFC,FWWLT,FB,FCGSAT,FJP,FAS, & FSEAS, T2I, HC_SNOW, SNOW_FRA,SNOWALB, & QST12,ZFUNC,ZF1,ZA2,QV2, DT_FDDA, & FC2R,FC1SAT, DTPBL, RAW CHARACTER (LEN = 6) :: LAND_USE_TYPE ALN10 = ALOG(10.0) G1000 = g*1.0E-3 WEIGHT = 0 IF (NLCAT == 50) THEN LAND_USE_TYPE = 'NLCD50' ELSE IF (NLCAT == 40) THEN LAND_USE_TYPE = 'NLCD40' ELSE IF (NLCAT == 20) THEN LAND_USE_TYPE = 'MODIS' ELSE IF (NLCAT == 21) THEN LAND_USE_TYPE = 'MODIS' ELSE IF (NLCAT == 24) THEN LAND_USE_TYPE = 'USGS' ELSE IF (NLCAT == 28) THEN LAND_USE_TYPE = 'USGS28' ELSE call wrf_error_fatal3("",349,& "Error: Unknown Land Use Category") END IF IF (ITIMESTEP .EQ. 1) THEN CALL wrf_message( 'PX LSM will use the ' // TRIM(LAND_USE_TYPE) // ' landuse tables' ) PRINT *, 'The analysis interval for surface soil and temp nudging = ',ANAL_INTERVAL,'sec.' ENDIF IF (ANAL_INTERVAL .LE. 0.0 .AND. PXLSM_SOIL_NUDGE .EQ. 1) THEN CALL wrf_message('PX LSM Error: The User specified analysis interval is zero or negative.') CALL wrf_message('If the PX LSM is used with soil nudging (pxlsm_soil_nudge=1) a wrfsfdda_d0* file is required.') CALL wrf_message('Make sure these files are present and') CALL wrf_message('Check the namelist to ensure sgfdda_interval_m is set to proper sfc analysis interval') STOP ENDIF IF(PXLSM_SOIL_NUDGE .EQ. 1) THEN DT_FDDA = ANAL_INTERVAL * 1.0 TIME_BETWEEN_ANALYSIS = MOD(CURR_SECS,DT_FDDA) IF (TIME_BETWEEN_ANALYSIS .EQ. 0.0) THEN CORB = 1.0 CORE = 0.0 ELSE CORE = TIME_BETWEEN_ANALYSIS / DT_FDDA CORB = 1.0 - CORE ENDIF ENDIF CALL VEGELAND(LANDUSEF, VEGFRA, SHDMIN, SHDMAX, & SOILCTOP, SOILCBOT, NLCAT, NSCAT, & ZNT,XLAI,XLAIMN,RSTMIN,XVEG,XVEGMN,XSNUP, & XLAND, XALB,XSNOALB,WETFRA,IMPERV,CANFRA, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte, LAND_USE_TYPE ) DO J = jts,jte DO I = its,ite IFLAND = XLAND(I,J) IF (IFLAND .LT. 1.5 ) THEN CALL SOILPROP (SOILCBOT(I,:,J), WEIGHT, & ITIMESTEP, MAVAIL(I,J), & PXLSM_SMOIS_INIT, & FWSAT,FWFC,FWWLT,FB,FCGSAT, & FJP,FAS,FC2R,FC1SAT,ISTI,SMOIS(I,1,J), & SMOIS(I,2,J)) ISLTYP(I,J) = ISTI ELSE ISLTYP(I,J) = 14 ENDIF SFCPRS = PSFC(i,j) / 1000.0 TA1 = T3D(i,1,j) DENS1 = RHO(I,1,J) QV1 = QV3D(i,1,j) QV2 = QV3D(i,2,j) ZLVL = 0.5 * DZ8W(i,1,j) ZF1 = DZ8W(i,1,j) ZA2 = ZF1 + 0.5 * DZ8W(i,2,j) LWDN = GLW(I,J) EMISSI = EMISS(I,J) PRECIP = MAX ( 1.0E-3*RAINBL(i,j)/DTBL,0.0) WR = 1.0E-3*CANWAT(I,J) THETA1 = TH3D(i,1,j) SNOBS = SNOW(I,J) IF(PXLSM_SOIL_NUDGE .EQ. 1) THEN T2_OLD = T2_NDG_OLD(I,J) T2_NEW = T2_NDG_NEW(I,J) VAPPRS = SVP1 * EXP(SVP2 * (T2_OLD - SVPT0) / ( T2_OLD - SVP3)) QSBT = EP_2 * VAPPRS / (SFCPRS - VAPPRS) RH2_OLD = Q2_NDG_OLD(I,J) / QSBT VAPPRS = SVP1 * EXP(SVP2 * (T2_NEW - SVPT0) / (T2_NEW - SVP3)) QSBT = EP_2 * VAPPRS / (SFCPRS - VAPPRS) RH2_NEW = Q2_NDG_NEW(I,J) / QSBT RH2OBS = CORB * RH2_OLD + CORE * RH2_NEW T2OBS(I,J) = CORB * T2_OLD + CORE * T2_NEW Q2OBS(I,J) = CORB * Q2_NDG_OLD(I,J) + CORE * Q2_NDG_NEW(I,J) SNOBS = CORB * SN_NDG_OLD(I,J) + CORE * SN_NDG_NEW(I,J) ENDIF USTAR = MAX(UST(I,J),0.005) IF (IFLAND .GE. 1.5) THEN ZNT(I,J) = CZO * UST(I,J) * UST(I,J) / G + OZO ENDIF Z0 = ZNT(I,J) CPAIR = CPD * (1.0 + 0.84 * QV1) SNOALB(I,J) = XSNOALB(I,J) CALL PXSNOW (ITIMESTEP, SNOBS, SNOWNCV(I,J), SNOW(I,J), & SNOWH(I,J), XSNUP(I,J), XALB(i,j), & SNOALB(I,J),VEGF_PX(I,J), SHDMIN(I,J), & HC_SNOW, SNOW_FRA, SNOWC(I,J), ALBEDO(I,J) ) IF( (XICE(I,J).GE.0.5) .OR. & (SST(I,J).LE.270.0.AND.XLAND(I,J).GE.1.50) ) THEN XLAND(I,J) = 1.0 IFLAND = 1.0 ZNT(I,J) = 0.001 SMOIS(I,1,J) = 1.0 SMOIS(I,2,J) = 1.0 XICE(I,J) = 1.0 ALBEDO(I,J) = 0.7 SNOWC(I,J) = 1.0 SNOW_FRA = 1.0 VEGF_PX(I,J) = 0.0 LAI(I,J) = 0.0 ENDIF T2I = TSLB(I,2,J) FSEAS = AMAX1(1.0 - 0.015625 * (290.0 - T2I) ** 2,0.0) IF (T2I .GE. 290.0) FSEAS = 1.0 LAI(I,J) = XLAIMN(I,J) + FSEAS*(XLAI(I,J) - XLAIMN(I,J)) VEGF_PX(I,J) = XVEGMN(I,J) + FSEAS*(XVEG(I,J) - XVEGMN(I,J)) IF (IFLAND .GE. 1.5) THEN VEGF_PX(I,J) = 0.0 LAI(I,J) = 0.0 ENDIF SOLDN = GSW(I,J) / (1.0 - ALBEDO(I,J)) ISNOW = SNOWC(I,J) NUDGE=PXLSM_SOIL_NUDGE IF ( J .LE. 2 .OR. J .GE. (jde-1) ) NUDGE=0 IF ( I .LE. 2 .OR. I .GE. (ide-1) ) NUDGE=0 IF ( RMOL(I,J) .GT. 0.0 ) THEN MOLX = AMIN1(1/RMOL(I,J),1000.0) ELSE IF ( RMOL(I,J) .LT. 0.0 ) THEN MOLX = AMAX1(1/RMOL(I,J),-1000.0) ELSE MOLX = 1000 ENDIF ZFUNC = ZF1 * (1.0 - ZF1 / MAX(100.,PBLH(I,J))) ** 2 QST12 = KARMAN * ZFUNC*(QV2-QV1) / (ZA2-ZLVL) NTSPS = INT(DT / (DTPBLX + 0.000001) + 1.0) DTPBL = DT / NTSPS DO IT=1,NTSPS IF ( TSLB(I,1,J) .LE. SVPT0 ) THEN ES = SVP1 * EXP(22.514 - 6.15E3 / TSLB(I,1,J)) ELSE ES = SVP1 * EXP(SVP2 * (TSLB(I,1,J) - SVPT0) / (TSLB(I,1,J) - SVP3)) ENDIF QSS = ES * 0.622 / (SFCPRS - ES) BETAP = 1.0 IF (IFLAND .LT. 1.5 .AND. ISNOW .LT. 0.5 .AND. SMOIS(I,1,J) .LE. FWFC) THEN BETAP = 0.25 * (1.0 - COS(SMOIS(I,1,J) / FWFC * PI)) ** 2 ENDIF CALL SURFPX (DTPBL, IFLAND, SNOWC(I,J), NUDGE, XICE(I,J), & SOLDN, GSW(I,J), LWDN, EMISSI, ZLVL, & MOLX, Z0, USTAR, & SFCPRS, DENS1, QV1, QSS, TA1, & THETA1, PRECIP, & CPAIR, PSIH(I,J), & RH2OBS,T2OBS(I,J), & VEGF_PX(I,J), ISTI, LAI(I,J), IMPERV(I,J), CANFRA(I,J), & BETAP, RSTMIN(I,J), HC_SNOW, SNOW_FRA, WETFRA(I,J), & FWWLT, FWFC, FCGSAT, FWSAT, FB, & FC1SAT,FC2R,FAS,FJP,DZS(1),DZS(2),QST12, & RADNET(I,J), GRDFLX(I,J), HFX(I,J), QFX(I,J), LH(I,J), & EG(I,J), ER(I,J), ETR(I,J), & QST(I,J), CAPG(I,J), RS(I,J), RA(I,J), & TSLB(I,1,J), TSLB(I,2,J), & SMOIS(I,1,J), SMOIS(I,2,J), WR, & TA2(I,J), QA2(I,J), LAND_USE_TYPE,I,J ) END DO TSK(I,J) = TSLB(I,1,J) CANWAT(I,J) = WR * 1000. RAW = RA(I,J) + 4.503 / USTAR QSFC(I,J) = QFX(I,J) * RAW / DENS1 + QV1 ENDDO ENDDO END SUBROUTINE pxlsm SUBROUTINE VEGELAND( landusef, vegfra, & shdmin, shdmax, & soilctop, soilcbot, nlcat, nscat,znt, xlai, & xlaimn, rstmin, xveg, xvegmn, xsnup, xland, & xalb, xsnoalb, wetfra, imperv, canfra, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte, & LAND_USE_TYPE ) IMPLICIT NONE INTEGER, INTENT(IN) :: ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte INTEGER , INTENT(IN) :: NSCAT, NLCAT REAL, DIMENSION( ims:ime , 1:NLCAT, jms:jme ), INTENT(IN) :: LANDUSEF REAL, DIMENSION( ims:ime , 1:NSCAT, jms:jme ), INTENT(IN) :: SOILCTOP, SOILCBOT REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN) :: VEGFRA, SHDMIN, SHDMAX REAL, DIMENSION( ims:ime, jms:jme ), INTENT(INOUT) :: ZNT, IMPERV, CANFRA REAL, DIMENSION( ims:ime, jms:jme ), INTENT(OUT) :: XLAI, XLAIMN, RSTMIN, XALB, & XVEG, XVEGMN, XSNUP, XLAND, & WETFRA, XSNOALB CHARACTER (LEN = 6), INTENT(IN) :: LAND_USE_TYPE INTEGER :: ITF, JTF, K, J, I REAL :: SUMLAI, SUMLMN, SUMRSI, SUMLZ0, SUMVEG, SUMVMN, & ALAI, VEGF, SUMSNUP, SUMALB, SUMSNOALB REAL :: VFMX, VFMN, VSEAS, FAREA, FWAT, ZNOTC, FCAN, FIMP, FORFRA, EXTFOR REAL, DIMENSION( NLCAT ) :: LAIMX, LAIMN, Z0, VEG, VEGMN, SNUP, ALB, SNOALB REAL, PARAMETER :: ZNOTCMN = 5.0 REAL, PARAMETER :: ZNOTCMX = 15.0 REAL, SAVE, DIMENSION(:), POINTER :: RSMIN, Z00, VEG0, VEGMN0, LAI0, & LAIMN0, SNUP0, ALBF, SNOALBF INTEGER, SAVE :: KWAT, LIMIT1, LIMIT2 IF (LAND_USE_TYPE == 'USGS') THEN KWAT = 16 RSMIN => RSMIN_USGS Z00 => Z00_USGS VEG0 => VEG0_USGS VEGMN0 => VEGMN0_USGS LAI0 => LAI0_USGS LAIMN0 => LAIMN0_USGS SNUP0 => SNUP0_USGS ALBF => ALBF_USGS SNOALBF=> SNOALB_USGS LIMIT1 = 2 LIMIT1 = 6 ELSE IF (LAND_USE_TYPE == 'USGS28') THEN KWAT = 16 RSMIN => RSMIN_USGS28 Z00 => Z00_USGS28 VEG0 => VEG0_USGS28 VEGMN0 => VEGMN0_USGS28 LAI0 => LAI0_USGS28 LAIMN0 => LAIMN0_USGS28 SNUP0 => SNUP0_USGS28 ALBF => ALBF_USGS28 SNOALBF=> SNOALB_USGS28 LIMIT1 = 2 LIMIT1 = 6 ELSE IF (LAND_USE_TYPE == 'NLCD50') THEN KWAT = 1 RSMIN => RSMIN_NLCD50 Z00 => Z00_NLCD50 VEG0 => VEG0_NLCD50 VEGMN0 => VEGMN0_NLCD50 LAI0 => LAI0_NLCD50 LAIMN0 => LAIMN0_NLCD50 SNUP0 => SNUP0_NLCD50 ALBF => ALBF_NLCD50 SNOALBF=> SNOALB_NLCD50 LIMIT1 = 20 LIMIT1 = 43 ELSE IF (LAND_USE_TYPE == 'NLCD40') THEN KWAT = 17 RSMIN => RSMIN_NLCD40 Z00 => Z00_NLCD40 VEG0 => VEG0_NLCD40 VEGMN0 => VEGMN0_NLCD40 LAI0 => LAI0_NLCD40 LAIMN0 => LAIMN0_NLCD40 SNUP0 => SNUP0_NLCD40 ALBF => ALBF_NLCD40 SNOALBF=> SNOALB_NLCD40 LIMIT1 = 20 LIMIT1 = 43 ELSE IF (LAND_USE_TYPE == 'MODIS') THEN KWAT = 17 RSMIN => RSMIN_MODIS Z00 => Z00_MODIS VEG0 => VEG0_MODIS VEGMN0 => VEGMN0_MODIS LAI0 => LAI0_MODIS LAIMN0 => LAIMN0_MODIS SNUP0 => SNUP0_MODIS ALBF => ALBF_MODIS SNOALBF=> SNOALB_MODIS LIMIT1 = 12 LIMIT1 = 14 END IF DO J = jts,jte DO I = its,ite XLAI(I,J) = 0.0 XLAIMN(I,J) = 0.0 RSTMIN(I,J) = 9999.0 XVEG(I,J) = 0.0 XVEGMN(I,J) = 0.0 XSNUP(I,J) = 0.0 XALB(I,J) = 0.0 XSNOALB(I,J)= 0.0 ENDDO ENDDO DO J = jts,jte DO I = its,ite DO K=1,NLCAT LAIMX(K) = LAI0(K) LAIMN(K) = LAIMN0(K) Z0(K) = Z00(K) VEG(K) = VEG0(K) VEGMN(K) = VEGMN0(K) SNUP(K) = SNUP0(K) ALB(K) = ALBF(K) SNOALB(K)= SNOALBF(K) ENDDO SUMLAI = 0.0 SUMLMN = 0.0 SUMRSI = 0.0 SUMLZ0 = 0.0 SUMVEG = 0.0 SUMVMN = 0.0 ALAI = 0.0 SUMSNUP = 0.0 SUMALB = 0.0 SUMSNOALB = 0.0 VFMX = SHDMAX(I,J) VFMN = SHDMIN(I,J) VEGF = VEGFRA(I,J) IF(VFMX.GT.0.0.AND.LANDUSEF(I,KWAT,J).LT.1.00) THEN VSEAS = VEGF/VFMX IF(VSEAS.GT.1.0.OR.VSEAS.LT.0.0) THEN VSEAS = MIN(VSEAS,1.0) VSEAS = MAX(VSEAS,0.0) ENDIF ZNOTC = ZNOTCMN * (1-VSEAS) + ZNOTCMX * VSEAS DO K = 1, NLCAT IF (LAND_USE_TYPE == 'MODIS') THEN IF (K.EQ.12 .OR. K.EQ.14) THEN LAIMX(K) = LAIMN0(K) * (1-VSEAS) + LAI0(K) * VSEAS LAIMN(K) = LAIMX(K) VEG(K) = VEGMN0(K) * (1-VSEAS) + VEG0(K) * VSEAS VEGMN(K) = VEG(K) IF (K .EQ. 12) THEN Z0(K) = ZNOTC ELSE IF (K .EQ.14) THEN Z0(K) = 0.5 * (ZNOTC + Z00(K)) ENDIF ENDIF ELSE IF (LAND_USE_TYPE == 'NLCD50') THEN IF (K.EQ.20 .OR. K.EQ.43 .OR. K.EQ.45) THEN LAIMX(K) = LAIMN0(K) * (1-VSEAS) + LAI0(K) * VSEAS LAIMN(K) = LAIMX(K) VEG(K) = VEGMN0(K) * (1-VSEAS) + VEG0(K) * VSEAS VEGMN(K) = VEG(K) IF (K.EQ.20 .OR. K.EQ.43) THEN Z0(K) = ZNOTC ELSE IF (K.EQ.45) THEN Z0(K) = 0.5 * (ZNOTC + Z00(K)) ENDIF ENDIF ELSE IF (LAND_USE_TYPE == 'NLCD40') THEN IF (K.EQ.12 .OR. K.EQ.14 .OR. K.EQ.38) THEN LAIMX(K) = LAIMN0(K) * (1-VSEAS) + LAI0(K) * VSEAS LAIMN(K) = LAIMX(K) VEG(K) = VEGMN0(K) * (1-VSEAS) + VEG0(K) * VSEAS VEGMN(K) = VEG(K) IF (K.EQ.12 .OR. K.EQ.38) THEN Z0(K) = ZNOTC ELSE IF (K.EQ.14) THEN Z0(K) = 0.5 * (ZNOTC + Z00(K)) ENDIF ENDIF ELSE IF (LAND_USE_TYPE == 'USGS') THEN IF (K .GE. 2 .AND. K .LE. 6) THEN LAIMX(K) = LAIMN0(K) * (1-VSEAS) + LAI0(K) * VSEAS LAIMN(K) = LAIMX(K) VEG(K) = VEGMN0(K) * (1-VSEAS) + VEG0(K) * VSEAS VEGMN(K) = VEG(K) IF (K .GE. 2 .AND. K .LE. 4) THEN Z0(K) = ZNOTC ELSE IF (K .GE.5 .AND. K .LE. 6) THEN Z0(K) = 0.5 * (ZNOTC + Z00(K)) ENDIF ENDIF ELSE IF (LAND_USE_TYPE == 'USGS28') THEN IF (K .GE. 2 .AND. K .LE. 6) THEN LAIMX(K) = LAIMN0(K) * (1-VSEAS) + LAI0(K) * VSEAS LAIMN(K) = LAIMX(K) VEG(K) = VEGMN0(K) * (1-VSEAS) + VEG0(K) * VSEAS VEGMN(K) = VEG(K) IF (K .GE. 2 .AND. K .LE. 4) THEN Z0(K) = ZNOTC ELSE IF (K .GE.5 .AND. K .LE. 6) THEN Z0(K) = 0.5 * (ZNOTC + Z00(K)) ENDIF ENDIF END IF ENDDO ENDIF DO K = 1, NLCAT FAREA = LANDUSEF(I,K,J) SUMLAI = SUMLAI + LAIMX(K) * FAREA SUMLMN = SUMLMN + LAIMN(K) * FAREA ALAI = ALAI + FAREA SUMRSI = SUMRSI + FAREA * LAIMX(K) / RSMIN(K) SUMLZ0 = SUMLZ0 + FAREA * ALOG(Z0(K)) SUMVEG = SUMVEG + FAREA * VEG(K) SUMVMN = SUMVMN + FAREA * VEGMN(K) SUMSNUP = SUMSNUP+ FAREA * SNUP(K) SUMALB = SUMALB + FAREA * ALB(K) SUMSNOALB= SUMSNOALB + FAREA * SNOALB(K) ENDDO FWAT = LANDUSEF(I,KWAT,J) IF (FWAT .GE. 0.50) THEN XLAI(I,J) = LAIMX(KWAT) XLAIMN(I,J) = LAIMN(KWAT) RSTMIN(I,J) = RSMIN(KWAT) ZNT(I,J) = Z0(KWAT) XVEG(I,J) = VEG(KWAT) XVEGMN(I,J) = VEGMN(KWAT) XSNUP(I,J) = SNUP(KWAT) XALB(I,J) = ALB(KWAT) XSNOALB(I,J)= SNOALB(KWAT) ELSE IF (FWAT .GT. 0.10) THEN ALAI = ALAI - FWAT SUMLZ0 = SUMLZ0 - FWAT * ALOG(Z0(KWAT)) ENDIF XLAI(I,J) = SUMLAI / ALAI XLAIMN(I,J) = SUMLMN / ALAI RSTMIN(I,J) = SUMLAI / SUMRSI ZNT(I,J) = EXP(SUMLZ0/ALAI) XVEG(I,J) = SUMVEG / ALAI XVEGMN(I,J) = SUMVMN / ALAI XSNUP(I,J) = SUMSNUP XALB(I,J) = SUMALB XSNOALB(I,J)= SUMSNOALB ENDIF IF (FWAT .GT. 0.50) THEN ZNT(I,J) = Z0(KWAT) XALB(I,J) = ALB(KWAT) XSNOALB(I,J)= SNOALB(KWAT) ENDIF IF (LAND_USE_TYPE == 'USGS') THEN WETFRA(I,J) = LANDUSEF(I,17,J)+LANDUSEF(I,18,J) ELSE IF (LAND_USE_TYPE == 'USGS28') THEN WETFRA(I,J) = LANDUSEF(I,17,J)+LANDUSEF(I,18,J) ELSE IF (LAND_USE_TYPE == 'NLCD50') THEN WETFRA(I,J) = LANDUSEF(I,22,J)+LANDUSEF(I,23,J)+LANDUSEF(I,27,J)+LANDUSEF(I,28,J)+LANDUSEF(I,42,J) ELSE IF (LAND_USE_TYPE == 'NLCD40') THEN WETFRA(I,J) = LANDUSEF(I,39,J)+LANDUSEF(I,40,J)+LANDUSEF(I,11,J) ELSE IF (LAND_USE_TYPE == 'MODIS') THEN WETFRA(I,J) = LANDUSEF(I,11,J) END IF ZNT(I,J) = ZNT(I,J) * 0.01 XVEG(I,J) = XVEG(I,J) * 0.01 XVEGMN(I,J) = XVEGMN(I,J) * 0.01 XLAND(I,J) = 1.0 + FWAT XALB(I,J) = XALB(I,J) * 0.01 XSNOALB(I,J)= XSNOALB(I,J) * 0.01 FIMP = IMPERV(I,J) * 0.01 FCAN = CANFRA(I,J) * 0.01 IF (LAND_USE_TYPE == 'NLCD40') THEN XVEG(I,J) = MIN(XVEG(I,J),1.0-FIMP) XVEGMN(I,J) = MIN(XVEGMN(I,J),1.0-FIMP) XVEG(I,J) = MAX(XVEG(I,J),FCAN) XVEGMN(I,J) = MAX(XVEGMN(I,J),FCAN) FORFRA = LANDUSEF(I,39,J)+LANDUSEF(I,30,J)+LANDUSEF(I,29,J)+LANDUSEF(I,28,J) EXTFOR = FCAN - FORFRA IF (EXTFOR.GE.0.01) THEN XLAI(I,J) = LAIMX(30) * EXTFOR + XLAI(I,J) * (1-EXTFOR) XLAIMN(I,J) = LAIMN(30) * EXTFOR + XLAIMN(I,J) * (1-EXTFOR) ENDIF ENDIF ENDDO ENDDO END SUBROUTINE vegeland SUBROUTINE SURFPX(DTPBL, IFLAND, ISNOW, NUDGEX, XICE1, SOLDN, GSW, & LWDN, EMISSI, Z1, MOL, ZNT, UST, PSURF, DENS1, & QV1, QSS, TA1, THETA1, PRECIP, CPAIR, PSIH, & RH2OBS, T2OBS, VEGFRC, ISTI,LAI,IMPERV,CANFRA,BETAP, & RSTMIN, HC_SNOW, SNOW_FRA, WETFRA, WWLT, WFC, & CGSAT, WSAT, B, C1SAT, C2R, AS, JP, DS1, DS2, QST12, & RADNET, GRDFLX, HFX, QFX, LH, EG, ER, ETR, & QST, CAPG, RS, RA, TG, T2, WG, W2, WR, & TA2, QA2, LAND_USE_TYPE, I, J ) IMPLICIT NONE INTEGER , INTENT(IN) :: ISTI, NUDGEX, I, J REAL , INTENT(IN) :: DTPBL, DS1, DS2 REAL , INTENT(IN) :: IFLAND, ISNOW, XICE1 REAL , INTENT(IN) :: SOLDN, GSW, LWDN, EMISSI, Z1 REAL , INTENT(IN) :: ZNT REAL , INTENT(IN) :: PSURF, DENS1, QV1, QSS, TA1, THETA1, PRECIP REAL , INTENT(IN) :: CPAIR REAL , INTENT(IN) :: VEGFRC, LAI, IMPERV, CANFRA REAL , INTENT(IN) :: RSTMIN, HC_SNOW, SNOW_FRA, WETFRA REAL , INTENT(IN) :: WWLT, WFC, CGSAT, WSAT, B, C1SAT, C2R, AS, JP REAL , INTENT(IN) :: RH2OBS,T2OBS REAL , INTENT(IN) :: QST12 REAL , INTENT(OUT) :: RADNET, EG, ER, ETR REAL , INTENT(OUT) :: QST, CAPG, RS, TA2, QA2 REAL , INTENT(INOUT) :: TG, T2, WG, W2, WR, UST, RA, BETAP REAL , INTENT(INOUT) :: GRDFLX, QFX, HFX, LH, PSIH, MOL CHARACTER (LEN = 5), INTENT(IN) :: LAND_USE_TYPE REAL :: HF, LV, CQ4, WETSAT, SM2 REAL :: RAH, RAW, ET, W2CG, CG, CT, SOILFLX, CPOT, THETAG REAL :: ZOL, ZOBOL, ZNTOL, Y, Y0, PSIH15, YNT REAL :: WGNUDG, W2NUDG, ALPH1, ALPH2, BET1, BET2, T1P5 REAL :: CQ1, CQ2, CQ3, COEFFNP1, COEFFN, TSNEW, TSHLF, T2NEW REAL :: ROFF, WRMAX, PC, DWR, PNET, TENDWR, WRNEW REAL :: COF1, CFNP1WR, CFNWR, PG, FASS REAL :: TENDW2, W2NEW, W2HLF, W2REL, C1, C2, WEQ, CFNP1, CFN, WGNEW REAL :: ALN10, TMP1, TMP2, TMP3, AA, AB, TST, RBH, CTVEG REAL :: QST1,PHIH,PSIOB REAL :: T2NUD, T2NUDF REAL :: VAPPRS, QSBT, RH2MOD, IMF, VEGF, SOILF REAL :: ZOBS, GAMAH, BETAH, SIGF, BH, CT_SNOW, CT_IMPERV REAL, PARAMETER :: CV = 1.2E-5 PARAMETER (ZOBS = 1.5) PARAMETER (BH = 15.7) PARAMETER (GAMAH = 16. ) PARAMETER (BETAH = 5.0 ) PARAMETER (SIGF = 0.5) PARAMETER (CT_SNOW = 2.0E-5) PARAMETER (CT_IMPERV = 3.268E-6) ALN10 = ALOG(10.0) RADNET = SOLDN - (EMISSI *(STBOLT *TG **4 - LWDN)) CPOT= (100.0 / PSURF) ** ROVCP THETAG = TG * CPOT ZOL = Z1/MOL ZOBOL = ZOBS/MOL ZNTOL = ZNT/MOL IF (MOL .LT. 0.0) THEN Y = ( 1.0 - GAMAH * ZOL ) ** 0.5 Y0 = ( 1.0 - GAMAH * ZOBOL ) ** 0.5 YNT = ( 1.0 - GAMAH * ZNTOL ) ** 0.5 PSIH15 = 2.0 * ALOG((Y + 1.0) / (Y0 + 1.0)) PSIH = 2.0 * ALOG((Y + 1.0) / (YNT + 1.0)) PSIOB = 2.0 * ALOG((Y0 + 1.0) / (YNT + 1.0)) PHIH = 1.0 / Y ELSE IF((ZOL-ZNTOL).LE.1.0) THEN PSIH = -BETAH*(ZOL-ZNTOL) ELSE PSIH = 1.-BETAH-(ZOL-ZNTOL) ENDIF IF ((ZOBOL - ZNTOL) .LE. 1.0) THEN PSIOB = -BETAH * (ZOBOL - ZNTOL) ELSE PSIOB = 1.0 - BETAH - (ZOBOL - ZNTOL) ENDIF PSIH15 = PSIH - PSIOB IF(ZOL.LE.1.0) THEN PHIH = 1.0 + BETAH * ZOL ELSE PHIH = BETAH + ZOL ENDIF ENDIF RA=PR0* ( ALOG(Z1/ZNT) - PSIH )/(KARMAN*UST) RAH = RA + 5.0 / UST RAW = RA + 4.503 / UST IF (IFLAND .LT. 1.5 ) THEN WETSAT = 1.00 * WSAT SM2 = (WETFRA * WETSAT) W2 = AMAX1(SM2, W2) ENDIF CALL QFLUX( DENS1, QV1, TA1, SOLDN, RAW, QSS, & VEGFRC, ISNOW, ISTI, IFLAND, LAI, BETAP, & WG, W2, WR, & RSTMIN, WWLT, WFC, & EG, ER, ETR, CQ4, RS, FASS) ET = EG + ER + ETR QST = -ET / (DENS1 * UST) LV = 2.83E6 IF (ISNOW .LT. 0.5.AND.TG.GT.273.15) & LV = (2.501 - 0.00237 * (TG - 273.15)) * 1.E6 QFX = ET LH = LV * QFX TST = (THETA1 - THETAG ) / (UST*RAH) HF = UST * TST HFX = AMAX1(-DENS1 * CPAIR * HF, -250.0) QST1 = 0.5*(QST+QST12/PHIH) TA2 = (THETAG + TST * (PR0 / KARMAN * (ALOG(ZOBS / ZNT) - PSIOB)+5.))/CPOT QA2 = QV1 - QST1 * PR0/ KARMAN * (ALOG(Z1 / ZOBS) - PSIH15) IF (QA2.LE.0.0) QA2 = QV1 VAPPRS = SVP1 * EXP(SVP2 * (TA2 - SVPT0) / (TA2 - SVP3)) QSBT = EP_2 * VAPPRS / (PSURF - VAPPRS) RH2MOD = QA2 / QSBT IF (IFLAND .LT. 1.5 ) THEN W2CG = AMAX1(W2,WWLT) CG = CGSAT * 1.0E-6 * (WSAT/ W2CG) ** & (0.5 * B / ALN10) IMF = AMAX1(0.0,IMPERV/100.0) VEGF = (1.0 - IMF) * VEGFRC SOILF= (1.0 - IMF) * (1.0 - VEGFRC) CT = 1./( IMF/CT_IMPERV + VEGF/CV + SOILF/CG) CT = 1./(SNOW_FRA/CT_SNOW + (1-SNOW_FRA)/CT) CAPG = 1.0/CT SOILFLX = 2.0 * PI * TAUINV * (TG - T2) GRDFLX = SOILFLX / CT ENDIF IF (IFLAND .LT. 1.5) THEN IF (NUDGEX .EQ. 0) THEN WGNUDG = 0.0 W2NUDG = 0.0 T2NUD = 0.0 ELSE CALL SMASS (ISTI, FASS, SOLDN, VEGFRC, RA, WWLT, WFC, & ALPH1, ALPH2, BET1, BET2, T2NUDF) WGNUDG = ALPH1 * (T2OBS - TA2) + ALPH2 * (RH2OBS - RH2MOD) * 100 W2NUDG = BET1 * (T2OBS - TA2) + BET2 * (RH2OBS - RH2MOD) * 100 IF (W2 .GE. WFC) W2NUDG = AMIN1(W2NUDG,0.0) IF (W2 .LE. WWLT) W2NUDG = AMAX1(W2NUDG,0.0) T2NUD = T2NUDF * (T2OBS - TA2) ENDIF ENDIF IF (IFLAND .LT. 1.5) THEN CQ1 = (1.0 - 0.622 * LV * CRANKP / (r_d * TG)) * QSS CQ2 = 0.622 * LV * QSS * CRANKP / (r_d * TG * TG) CQ3 = DENS1 * BETAP * (1.0 - VEGFRC) / RAW COEFFNP1 = 1.0 + DTPBL * CRANKP * (4.0 * EMISSI * STBOLT * TG ** 3 & * CT + DENS1 * CPAIR / RAH * CPOT * CT + 2.0 * PI & * TAUINV ) + DTPBL * (CT * LV * CQ2 * (CQ3 + CQ4)) COEFFN = CT * (GSW + EMISSI * (STBOLT * (4.0 * CRANKP - 1.0) & * TG*TG*TG*TG + LWDN) & + DENS1 * CPAIR / RAH * (THETA1 - (1.0 - CRANKP) * THETAG) & - LV * (CQ3 * (CQ1 - QV1) + CQ4 * (CQ1 - QV1))) & - 2.0 * PI * TAUINV * ((1.0 - CRANKP) * TG - T2) TSNEW = (TG + DTPBL * COEFFN) / COEFFNP1 IF (XICE1 .GT. 0.5) TSNEW = AMIN1(TSNEW,273.15) TSHLF = 0.5 * ( TSNEW + TG) T2NEW = (T2 + DTPBL * TAUINV * T2TFAC * (TSHLF - (1 - CRANKP) * T2) & + DTPBL*T2NUD) & / (1.0 + DTPBL * TAUINV * T2TFAC * CRANKP) TG = TSNEW T2 = T2NEW ENDIF IF (IFLAND .LT. 1.5.AND. XICE1.LT.0.5) THEN ROFF = 0.0 WRMAX = 0.2E-3 * VEGFRC * LAI IF(WRMAX.GT.0.0) THEN PC = VEGFRC * SIGF * PRECIP DWR = (WRMAX - WR) / DTPBL PNET = PC - ER/ DENW IF (PNET .GT. DWR) THEN ROFF = PNET - DWR PC = PC - ROFF ENDIF IF (QSS .LT. QV1) THEN TENDWR = PC - ER / DENW WRNEW = WR + DTPBL * TENDWR ELSE COF1 = DENS1 / DENW * VEGFRC * (QSS - QV1) / RAW CFNP1WR = 1.0 + DTPBL * COF1 * CRANKP / WRMAX CFNWR = PC - COF1 * (1.0 - CRANKP) * WR / WRMAX WRNEW = (WR + DTPBL * CFNWR) / CFNP1WR ENDIF ELSE PC=0.0 WRNEW=0.0 ENDIF PG = DENW * (PRECIP - PC) TENDW2 = 1.0 / (DENW * DS2) * (PG - EG - ETR) & + (W2NUDG + WGNUDG) / DS2 W2NEW = W2 + DTPBL * TENDW2 W2NEW = AMIN1(W2NEW,WSAT) W2NEW = AMAX1(W2NEW,0.05) W2HLF = 0.5 * (W2 + W2NEW) W2 = W2NEW WR = AMIN1(WRMAX,WRNEW) ENDIF IF (IFLAND .LT. 1.5.AND. XICE1.LT.0.5) THEN IF (ISNOW .GT.0.5) THEN WG = WSAT ELSE W2REL = W2HLF / WSAT IF (WG .GT. WWLT) THEN C1 = C1SAT * (WSAT / WG) ** (0.5 * B + 1.0) ELSE C1 = C1SAT * (WSAT / WWLT) ** (0.5 * B + 1.0) ENDIF C2 = C2R * W2HLF / (WSAT - W2HLF + 1.E-11) IF (W2HLF .GE. WSAT) THEN WEQ = WSAT ELSE WEQ = W2HLF - AS * WSAT * W2REL ** JP * & (1.0 - W2REL ** (8.0 * JP)) ENDIF CFNP1 = 1.0 + DTPBL * C2 * TAUINV * CRANKP CFN = C1 / (DENW * DS1) * (PG - EG) - C2 * TAUINV * & ((1.0 - CRANKP) * WG - WEQ) + WGNUDG/ DS1 WGNEW = AMAX1((WG + DTPBL * CFN) / CFNP1,0.001) WG = AMIN1(WGNEW,WSAT) ENDIF ENDIF END SUBROUTINE surfpx SUBROUTINE QFLUX (DENS1, QV1, TA1, RG, RAW, QSS, & VEGFRC, ISNOW, ISTI, IFLAND, LAI, BETAP, & WG, W2, WR, & RSTMIN, WWLT, WFC, & EG, ER, ETR, CQ4, RS, FASS) IMPLICIT NONE INTEGER , INTENT(IN) :: ISTI REAL , INTENT(IN) :: ISNOW, IFLAND REAL , INTENT(IN) :: DENS1, QV1, TA1, RG, RAW, QSS, & VEGFRC, LAI, & WG, W2, WR, RSTMIN REAL , INTENT(INOUT) :: BETAP REAL, INTENT(IN) :: WWLT, WFC REAL , INTENT(OUT) :: EG, ER, ETR, CQ4, RS, FASS REAL :: WRMAX, DELTA, SIGG, RADL, RADF, W2AVAIL, W2MXAV REAL :: FTOT, F1, F2, F3, F4 REAL :: FSHELT, GS, GA, FX REAL :: PAR, F1MAX REAL, PARAMETER :: RSMAX = 5000.0 REAL, PARAMETER :: FTMIN = 0.0000001 REAL, PARAMETER :: F3MIN = 0.25 ER = 0.0 ETR = 0.0 CQ4 = 0.0 IF (QSS .LT. QV1) BETAP = 1.0 EG = DENS1 * (1.0 - VEGFRC) * BETAP * (QSS - QV1) / RAW IF (IFLAND .LT. 1.5 .AND. VEGFRC .GT. 0.0) THEN WRMAX = 0.2E-3 * VEGFRC * LAI IF (WR .LE. 0.0) THEN DELTA = 0.0 ELSE DELTA = WR / WRMAX ENDIF IF (QSS .GE. QV1) THEN SIGG = DELTA ELSE SIGG = 1.0 ENDIF ER = DENS1 * VEGFRC * SIGG * (QSS - QV1) / RAW ENDIF IF (IFLAND .LT. 1.5 .AND. VEGFRC .GT. 0.0) THEN IF (RSTMIN .GT. 130.0) THEN F1MAX = 1.-0.02*LAI ELSE F1MAX = 1.-0.07*LAI ENDIF PAR = 0.45 * RG * 4.566 F1 = F1MAX*(1.0-exp(-0.0017*PAR)) F1 = AMAX1(F1,RSTMIN / RSMAX) W2AVAIL = W2 - WWLT W2MXAV = WFC - WWLT F2 = 1.0 / (1.0 + EXP(-5.0 * ( W2AVAIL / W2MXAV - & (W2MXAV / 3.0 + WWLT)))) IF (TA1 .LE. 302.15) THEN F4 = 1.0 / (1.0 + EXP(-0.41 * (TA1 - 282.05))) ELSE F4 = 1.0 / (1.0 + EXP(0.5 * (TA1 - 314.0))) ENDIF FTOT = LAI * F1 * F2 * F4 ENDIF IF (IFLAND .LT. 1.5 .AND. VEGFRC .GT. 0.0) THEN FSHELT = 1.0 GS = FTOT / (RSTMIN * FSHELT) GA = 1.0 / RAW F3 = 0.5 * (GS - GA + SQRT(GA * GA + GA * GS * & (4.0 * QV1 / QSS - 2.0) + GS * GS)) / GS F3 = AMIN1(AMAX1(F3,F3MIN),1.0) RS = 1.0 / (GS * F3) IF (RG .LT. 0.00001) THEN FX = 0.0 ELSE FX = 30.0 * F1 * F4 * LAI / (RSTMIN * FSHELT) ENDIF FASS = FX ETR = DENS1 * VEGFRC * (1.0 - SIGG) * (QSS - QV1) / (RAW + RS) CQ4 = DENS1 * VEGFRC * ((1.0 - SIGG) / (RAW + RS) + SIGG / RAW) ENDIF END SUBROUTINE qflux SUBROUTINE SMASS (ISTI, FASS, RG, VEGFRC, RA, & WWLT, WFC, & ALPH1, ALPH2, BET1, BET2, T2NUDF ) IMPLICIT NONE INTEGER, PARAMETER :: NSCAT = 16 INTEGER, INTENT(IN) :: ISTI REAL, INTENT(IN) :: FASS, RG, VEGFRC, RA REAL, INTENT(IN) :: WWLT, WFC REAL, INTENT(OUT) :: ALPH1, ALPH2, BET1, BET2, T2NUDF REAL :: FBET, FALPH, FRA, FTEXT REAL, DIMENSION( 1: NSCAT ) :: WFCX, WWLTX REAL, PARAMETER :: A1MAX = -10.E-5, A2MAX = 1.E-5 REAL, PARAMETER :: B1MAX = -10.E-3, B2MAX = 1.E-3 REAL, PARAMETER :: TASSI = 4.6296E-5 REAL, PARAMETER :: RAMIN = 10.0 DATA WFCX / 0.135, 0.150, 0.195, 0.255, 0.240, 0.255, 0.322, & 0.325, 0.310, 0.370, 0.367, 0.367, 0.367, 0.367, 0.367, 0.367 / DATA WWLTX / 0.068, 0.075, 0.114, 0.179, 0.155, 0.175, 0.218, & 0.250, 0.219, 0.283, 0.286, 0.286, 0.286, 0.286, 0.286, 0.286 / FBET = FASS FALPH = RG / 1370.0 FRA = RAMIN / RA FTEXT = TASSI * (WWLT + WFC) / (WWLTX(5) + WFCX(5)) * FRA ALPH1 = A1MAX * FALPH * (1.0 - VEGFRC) * FTEXT ALPH2 = A2MAX * FALPH * (1.0 - VEGFRC) * FTEXT BET1 = B1MAX * FBET * VEGFRC * FTEXT BET2 = B2MAX * FBET * VEGFRC * FTEXT T2NUDF = 1.0E-5 * MAX((1.0 - 5.0 * FALPH),0.0) END SUBROUTINE smass SUBROUTINE SOILPROP (SOILCBOT,WEIGHT, ITIMESTEP, MAVAIL, & PXLSM_SMOIS_INIT, & FWSAT,FWFC,FWWLT,FB,FCGSAT, & FJP,FAS,FC2R,FC1SAT,ISTI, WG, W2 ) IMPLICIT NONE INTEGER, PARAMETER :: NSCAT = 16 INTEGER, PARAMETER :: NSCATMIN = 16 INTEGER, INTENT(IN) :: WEIGHT, ITIMESTEP, PXLSM_SMOIS_INIT REAL, INTENT(IN) :: MAVAIL REAL, DIMENSION(1:NSCAT), INTENT(IN) :: SOILCBOT REAL, INTENT(OUT) :: FWSAT,FWFC,FWWLT,FB,FCGSAT, & FJP,FAS,FC2R,FC1SAT REAL, INTENT(INOUT) :: W2, WG INTEGER, INTENT(OUT) :: ISTI CHARACTER*4, AVCLASS CHARACTER*4, DIMENSION( 1: NSCAT ) :: TEXID INTEGER:: S REAL:: TFRACBOT, CFRAC, SUMSND, SUMCLY, AVS, AVC, AVSLT REAL, DIMENSION( 1: NSCAT ) :: WSAT, WFC, WWLT, B, CGSAT, AS, & JP, C2R, C1SAT REAL, DIMENSION( 1: NSCATMIN ) :: SAND, CLAY DATA SAND /92.5,80.5,61.1,19.6,4.0,40.0,57.1,11.3,26.8, & 52.0,6.5,10.2,1.0,1.0,1.0,1.0/ DATA CLAY/2.1,4.1,10.9,19.1,7.3,18.8,23.3,32.2,36.6, & 43.0,46.2,58.8,1.0,1.0,1.0,1.0/ DATA TEXID/'Sand','Lsan','Sloa','Sill','Silt','Loam','Sclo', & 'Sicl','Cllo','Sacl','Sicy','Clay','Ormt','Wate', & 'Bedr','Othe'/ DATA WSAT / 0.395, 0.410, 0.435, 0.485, 0.451, 0.420, 0.477, & 0.476, 0.426, 0.482, 0.482, 0.482, 0.482, 0.482, 0.482, 0.482 / DATA WFC / 0.135, 0.150, 0.195, 0.255, 0.240, 0.255, 0.322, & 0.325, 0.310, 0.370, 0.367, 0.367, 0.367, 0.367, 0.367, 0.367 / DATA WWLT / 0.068, 0.075, 0.114, 0.179, 0.155, 0.175, 0.218, & 0.250, 0.219, 0.283, 0.286, 0.286, 0.286, 0.286, 0.286, 0.286 / DATA B / 4.05, 4.38, 4.90, 5.30, 5.39, 7.12, 7.75, & 8.52, 10.40, 10.40, 11.40, 11.40, 11.40, 11.40, 11.40, 11.40 / DATA CGSAT / 3.222, 3.057, 3.560, 4.418, 4.111, 3.670, 3.593, & 3.995, 3.058, 3.729, 3.600, 3.600, 3.600, 3.600, 3.600, 3.600 / DATA JP / 4, 4, 4, 6, 6, 6, 8, & 10, 8, 10, 12, 12, 12, 12, 12, 12 / DATA AS / 0.387, 0.404, 0.219, 0.105, 0.148, 0.135, 0.127, & 0.084, 0.139, 0.075, 0.083, 0.083, 0.083, 0.083, 0.083, 0.083 / DATA C2R / 3.9, 3.7, 1.8, 0.8, 0.8, 0.8, 0.4, & 0.6, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3, 0.3 / DATA C1SAT / 0.082, 0.098, 0.132, 0.153, 0.191, 0.213, 0.385, & 0.227, 0.421, 0.375, 0.342, 0.342, 0.342, 0.342, 0.342, 0.342 / IF(WEIGHT.GE.1.0) THEN FWSAT =0 FWFC =0 FWWLT =0 FB =0 FCGSAT=0 FJP =0 FAS =0 FC2R =0 FC1SAT=0 TFRACBOT =0 CFRAC=0 DO S=1,NSCAT IF(SOILCBOT(S).GE.CFRAC) THEN ISTI=S CFRAC=SOILCBOT(S) ENDIF TFRACBOT=TFRACBOT+SOILCBOT(S) FWSAT =FWSAT + WSAT(S) *SOILCBOT(S) FWFC =FWFC + WFC(S) *SOILCBOT(S) FWWLT =FWWLT + WWLT(S) *SOILCBOT(S) FB =FB + B(S) *SOILCBOT(S) FCGSAT=FCGSAT + CGSAT(S) *SOILCBOT(S) FJP =FJP + JP(S) *SOILCBOT(S) FAS =FAS + AS(S) *SOILCBOT(S) FC2R =FC2R + C2R(S) *SOILCBOT(S) FC1SAT=FC1SAT + C1SAT(S) *SOILCBOT(S) ENDDO TFRACBOT = 1/TFRACBOT FWSAT =FWSAT * TFRACBOT FWFC =FWFC * TFRACBOT FWWLT =FWWLT * TFRACBOT FB =FB * TFRACBOT FCGSAT=FCGSAT * TFRACBOT FJP =FJP * TFRACBOT FAS =FAS * TFRACBOT FC2R =FC2R * TFRACBOT FC1SAT=FC1SAT * TFRACBOT ELSE CFRAC = 0.0 SUMSND = 0.0 SUMCLY = 0.0 TFRACBOT = 0.0 DO S = 1,12 TFRACBOT = TFRACBOT + SOILCBOT(S) SUMSND = SUMSND + SAND(S) * SOILCBOT(S) SUMCLY = SUMCLY + CLAY(S) * SOILCBOT(S) IF(SOILCBOT(S).GE.CFRAC) THEN ISTI=S CFRAC=SOILCBOT(S) ENDIF ENDDO IF(TFRACBOT.GT.0.001) THEN AVS = SUMSND / TFRACBOT AVC = SUMCLY / TFRACBOT AVSLT = 100 - AVS - AVC IF(AVS.GT.(85.+ 0.5*AVC)) THEN AVCLASS= 'Sand' ISTI = 1 ELSE IF(AVS.GT.(70.+ AVC)) THEN AVCLASS= 'Lsan' ISTI = 2 ELSE IF((AVC.LT.20..AND.AVS.GT.52.) & .OR.(AVC.LE.7.5.AND.AVSLT.LT.50.)) THEN AVCLASS= 'Sloa' ISTI = 3 ELSE IF(AVC.LT.35..AND.AVS.GT.45..AND.AVSLT.LT.28.) THEN AVCLASS= 'Sclo' ISTI = 6 ELSE IF(AVC.GE.35..AND.AVS.GT.45.) THEN AVCLASS = 'Sacl' ISTI = 9 ELSE IF(AVC.LT.27.0.AND.AVSLT.LT.50.) THEN AVCLASS= 'Loam' ISTI = 5 ELSE IF(AVC.LT.12..AND.AVSLT.GT.80.) THEN AVCLASS = 'Silt' ISTI = 4 ELSE IF(AVC.LT.27.) THEN AVCLASS = 'Sill' ISTI = 4 ELSE IF(AVC.LT.40..AND.AVS.GT.20.) THEN AVCLASS = 'Cllo' ISTI = 8 ELSE IF(AVC.LT.40.) THEN AVCLASS = 'Sicl' ISTI = 7 ELSE IF(AVSLT.GE.40.) THEN AVCLASS = 'Sicy' ISTI = 10 ELSE AVCLASS = 'Clay' ISTI = 11 ENDIF ELSE ISTI=5 AVCLASS = TEXID(ISTI) ENDIF FWSAT =WSAT(ISTI) FWFC =WFC(ISTI) FWWLT =WWLT(ISTI) FB =B(ISTI) FCGSAT=CGSAT(ISTI) FJP =JP(ISTI) FAS =AS(ISTI) FC2R =C2R(ISTI) FC1SAT=C1SAT(ISTI) ENDIF IF (ITIMESTEP .EQ. 1 .AND. PXLSM_SMOIS_INIT .GT. 0) THEN WG = FWWLT + (0.5*(FWSAT+FWFC) - FWWLT) * MAVAIL W2 = FWWLT + (0.5*(FWSAT+FWFC) - FWWLT) * MAVAIL ENDIF END SUBROUTINE soilprop SUBROUTINE PXSNOW (ITIMESTEP, ASNOW, CSNOW, SNOW, & SNOWH, SNUP, & ALB, SNOALB, SHDFAC, SHDMIN, & HC_SNOW, SNOW_FRA, SNOWC, SNOWALB) IMPLICIT NONE REAL, PARAMETER :: W2SN_CONV = 10.0 REAL, PARAMETER :: CS_SNOWPACK = 2092.0 REAL, PARAMETER :: SALP = 2.6 INTEGER, INTENT(IN) :: ITIMESTEP REAL, INTENT(IN) :: ASNOW, CSNOW, SNUP, ALB, SNOALB, SHDFAC, SHDMIN REAL, INTENT(INOUT) :: SNOW, SNOWH, SNOWC REAL, INTENT(OUT) :: HC_SNOW, SNOW_FRA, SNOWALB REAL:: CONV_WAT2SNOW, CSNOWW, RHO_SNOWPACK, & LIQSN_RATIO, SNEQV, RSNOW SNEQV=ASNOW*0.001 RHO_SNOWPACK = 450 LIQSN_RATIO = DENW/RHO_SNOWPACK CONV_WAT2SNOW = LIQSN_RATIO/1000 SNOW = ASNOW SNOWH= SNOW * CONV_WAT2SNOW SNOWC = 0.0 IF (SNOWH .GT. 0.005) SNOWC = 1.0 HC_SNOW = RHO_SNOWPACK * CS_SNOWPACK * SNOWH IF (SNEQV < SNUP) THEN RSNOW = SNEQV / SNUP SNOW_FRA = 1. - ( EXP ( - SALP * RSNOW) - RSNOW * EXP ( - SALP)) ELSE SNOW_FRA = 1.0 END IF SNOWC = SNOW_FRA SNOWALB = ALB + SNOWC*(SNOALB-ALB) END SUBROUTINE pxsnow END MODULE module_sf_pxlsm