module noahmp_globals use module_sf_noahlsm, only: & & SLCATS, & & LUCATS, & & CSOIL_DATA, & & BB, & & SATDK, & & SATDW, & & F11, & & SATPSI, & & QTZ, & & DRYSMC, & & MAXSMC, & & REFSMC, & & WLTSMC, & & RSTBL, & & RGLTBL, & & HSTBL, & & NROTBL, & & TOPT_DATA, & & RSMAX_DATA, & & ZBOT_DATA, & & CZIL_DATA, & & FRZK_DATA, & & SLOPE_DATA, & & REFDK_DATA, & & REFKDT_DATA implicit none REAL, PARAMETER :: GRAV = 9.80616 REAL, PARAMETER :: SB = 5.67E-08 REAL, PARAMETER :: VKC = 0.40 REAL, PARAMETER :: TFRZ = 273.16 REAL, PARAMETER :: HSUB = 2.8440E06 REAL, PARAMETER :: HVAP = 2.5104E06 REAL, PARAMETER :: HFUS = 0.3336E06 REAL, PARAMETER :: CWAT = 4.188E06 REAL, PARAMETER :: CICE = 2.094E06 REAL, PARAMETER :: CPAIR = 1004.64 REAL, PARAMETER :: TKWAT = 0.6 REAL, PARAMETER :: TKICE = 2.2 REAL, PARAMETER :: TKAIR = 0.023 REAL, PARAMETER :: RAIR = 287.04 REAL, PARAMETER :: RW = 461.269 REAL, PARAMETER :: DENH2O = 1000. REAL, PARAMETER :: DENICE = 917. INTEGER :: NROOT REAL :: RGL REAL :: RSMIN REAL :: HS REAL :: RSMAX REAL :: TOPT REAL :: BEXP REAL :: SMCDRY REAL :: F1 REAL :: SMCMAX REAL :: SMCREF REAL :: PSISAT REAL :: DKSAT REAL :: DWSAT REAL :: SMCWLT REAL :: QUARTZ REAL :: SLOPE REAL :: CSOIL REAL :: ZBOT REAL :: CZIL REAL :: KDT REAL :: FRZX INTEGER :: DVEG INTEGER :: OPT_CRS INTEGER :: OPT_BTR INTEGER :: OPT_RUN INTEGER :: OPT_SFC INTEGER :: OPT_FRZ INTEGER :: OPT_INF INTEGER :: OPT_RAD INTEGER :: OPT_ALB INTEGER :: OPT_SNF INTEGER :: OPT_TBOT INTEGER :: OPT_STC REAL, PARAMETER :: TIMEAN = 10.5 REAL, PARAMETER :: FSATMX = 0.38 REAL, PARAMETER :: M = 1.0 REAL, PARAMETER :: Z0SNO = 0.002 REAL, PARAMETER :: SSI = 0.03 REAL, PARAMETER :: SWEMX = 1.00 END MODULE NOAHMP_GLOBALS MODULE NOAHMP_VEG_PARAMETERS IMPLICIT NONE INTEGER, PARAMETER :: MAX_VEG_PARAMS = 33 INTEGER, PARAMETER :: MVT = 27 INTEGER, PARAMETER :: MBAND = 2 INTEGER, PRIVATE :: ISURBAN INTEGER :: ISWATER INTEGER :: ISBARREN INTEGER :: ISSNOW INTEGER :: EBLFOREST REAL :: CH2OP(MVT) REAL :: DLEAF(MVT) REAL :: Z0MVT(MVT) REAL :: HVT(MVT) REAL :: HVB(MVT) REAL :: DEN(MVT) REAL :: RC(MVT) REAL :: SAIM(MVT,12) REAL :: LAIM(MVT,12) REAL :: SLA(MVT) REAL :: DILEFC(MVT) REAL :: DILEFW(MVT) REAL :: FRAGR(MVT) REAL :: LTOVRC(MVT) REAL :: C3PSN(MVT) REAL :: KC25(MVT) REAL :: AKC(MVT) REAL :: KO25(MVT) REAL :: AKO(MVT) REAL :: VCMX25(MVT) REAL :: AVCMX(MVT) REAL :: BP(MVT) REAL :: MP(MVT) REAL :: QE25(MVT) REAL :: AQE(MVT) REAL :: RMF25(MVT) REAL :: RMS25(MVT) REAL :: RMR25(MVT) REAL :: ARM(MVT) REAL :: FOLNMX(MVT) REAL :: TMIN(MVT) REAL :: XL(MVT) REAL :: RHOL(MVT,MBAND) REAL :: RHOS(MVT,MBAND) REAL :: TAUL(MVT,MBAND) REAL :: TAUS(MVT,MBAND) REAL :: MRP(MVT) REAL :: CWPVT(MVT) REAL :: WRRAT(MVT) REAL :: WDPOOL(MVT) REAL :: TDLEF(MVT) INTEGER :: IK,IM REAL :: TMP10(MVT*MBAND) REAL :: TMP11(MVT*MBAND) REAL :: TMP12(MVT*MBAND) REAL :: TMP13(MVT*MBAND) REAL :: TMP14(MVT*12) REAL :: TMP15(MVT*12) REAL :: TMP16(MVT*5) real slarea(MVT) real eps(MVT,5) CONTAINS subroutine read_mp_veg_parameters(DATASET_IDENTIFIER) implicit none character(len=*), intent(in) :: DATASET_IDENTIFIER integer :: ierr REAL :: TMP10(MVT*MBAND) REAL :: TMP11(MVT*MBAND) REAL :: TMP12(MVT*MBAND) REAL :: TMP13(MVT*MBAND) REAL :: TMP14(MVT*12) REAL :: TMP15(MVT*12) REAL :: TMP16(MVT*5) integer :: NVEG character(len=256) :: VEG_DATASET_DESCRIPTION NAMELIST / noah_mp_usgs_veg_categories / VEG_DATASET_DESCRIPTION, NVEG NAMELIST / noah_mp_usgs_parameters / ISURBAN, ISWATER, ISBARREN, ISSNOW, EBLFOREST, & CH2OP, DLEAF, Z0MVT, HVT, HVB, DEN, RC, RHOL, RHOS, TAUL, TAUS, XL, CWPVT, C3PSN, KC25, AKC, KO25, AKO, AVCMX, AQE, & LTOVRC, DILEFC, DILEFW, RMF25 , SLA , FRAGR , TMIN , VCMX25, TDLEF , BP, MP, QE25, RMS25, RMR25, ARM, FOLNMX, WDPOOL, WRRAT, MRP, & SAIM, LAIM, SLAREA, EPS NAMELIST / noah_mp_modis_veg_categories / VEG_DATASET_DESCRIPTION, NVEG NAMELIST / noah_mp_modis_parameters / ISURBAN, ISWATER, ISBARREN, ISSNOW, EBLFOREST, & CH2OP, DLEAF, Z0MVT, HVT, HVB, DEN, RC, RHOL, RHOS, TAUL, TAUS, XL, CWPVT, C3PSN, KC25, AKC, KO25, AKO, AVCMX, AQE, & LTOVRC, DILEFC, DILEFW, RMF25 , SLA , FRAGR , TMIN , VCMX25, TDLEF , BP, MP, QE25, RMS25, RMR25, ARM, FOLNMX, WDPOOL, WRRAT, MRP, & SAIM, LAIM, SLAREA, EPS CH2OP = -1.E36 DLEAF = -1.E36 Z0MVT = -1.E36 HVT = -1.E36 HVB = -1.E36 DEN = -1.E36 RC = -1.E36 RHOL = -1.E36 RHOS = -1.E36 TAUL = -1.E36 TAUS = -1.E36 XL = -1.E36 CWPVT = -1.E36 C3PSN = -1.E36 KC25 = -1.E36 AKC = -1.E36 KO25 = -1.E36 AKO = -1.E36 AVCMX = -1.E36 AQE = -1.E36 LTOVRC = -1.E36 DILEFC = -1.E36 DILEFW = -1.E36 RMF25 = -1.E36 SLA = -1.E36 FRAGR = -1.E36 TMIN = -1.E36 VCMX25 = -1.E36 TDLEF = -1.E36 BP = -1.E36 MP = -1.E36 QE25 = -1.E36 RMS25 = -1.E36 RMR25 = -1.E36 ARM = -1.E36 FOLNMX = -1.E36 WDPOOL = -1.E36 WRRAT = -1.E36 MRP = -1.E36 SAIM = -1.E36 LAIM = -1.E36 SLAREA = -1.E36 EPS = -1.E36 open(15, file="MPTABLE.TBL", status='old', form='formatted', action='read', iostat=ierr) if (ierr /= 0) then write(*,'("****** Error ******************************************************")') write(*,'("Cannot find file MPTABLE.TBL")') write(*,'("STOP")') write(*,'("*******************************************************************")') call wrf_error_fatal3("",388,& "STOP in Noah-MP read_mp_veg_parameters") endif if ( trim(DATASET_IDENTIFIER) == "USGS" ) then read(15,noah_mp_usgs_veg_categories) read(15,noah_mp_usgs_parameters) else if ( trim(DATASET_IDENTIFIER) == "MODIFIED_IGBP_MODIS_NOAH" ) then read(15,noah_mp_modis_veg_categories) read(15,noah_mp_modis_parameters) else write(*,'("Unrecognized DATASET_IDENTIFIER in subroutine READ_MP_VEG_PARAMETERS")') write(*,'("DATASET_IDENTIFIER = ''", A, "''")') trim(DATASET_IDENTIFIER) call wrf_error_fatal3("",401,& "STOP in Noah-MP read_mp_veg_parameters") endif close(15) if ( MVT > NVEG ) then TMP10 = reshape( RHOL, (/ MVT*size(RHOL,2) /)) TMP11 = reshape( RHOS, (/ MVT*size(RHOS,2) /)) TMP12 = reshape( TAUL, (/ MVT*size(TAUL,2) /)) TMP13 = reshape( TAUS, (/ MVT*size(TAUS,2) /)) TMP14 = reshape( SAIM, (/ MVT*size(SAIM,2) /)) TMP15 = reshape( LAIM, (/ MVT*size(LAIM,2) /)) TMP16 = reshape( EPS, (/ MVT*size(EPS ,2) /)) RHOL(1:NVEG,:) = reshape( TMP10, (/ NVEG, size(RHOL,2) /)) RHOS(1:NVEG,:) = reshape( TMP11, (/ NVEG, size(RHOS,2) /)) TAUL(1:NVEG,:) = reshape( TMP12, (/ NVEG, size(TAUL,2) /)) TAUS(1:NVEG,:) = reshape( TMP13, (/ NVEG, size(TAUS,2) /)) SAIM(1:NVEG,:) = reshape( TMP14, (/ NVEG, size(SAIM,2) /)) LAIM(1:NVEG,:) = reshape( TMP15, (/ NVEG, size(LAIM,2) /)) EPS(1:NVEG,:) = reshape( TMP16, (/ NVEG, size(EPS,2) /)) RHOL(NVEG+1:MVT,:) = -1.E36 RHOS(NVEG+1:MVT,:) = -1.E36 TAUL(NVEG+1:MVT,:) = -1.E36 TAUS(NVEG+1:MVT,:) = -1.E36 SAIM(NVEG+1:MVT,:) = -1.E36 LAIM(NVEG+1:MVT,:) = -1.E36 EPS( NVEG+1:MVT,:) = -1.E36 endif end subroutine read_mp_veg_parameters END MODULE NOAHMP_VEG_PARAMETERS MODULE NOAHMP_RAD_PARAMETERS IMPLICIT NONE INTEGER I INTEGER, PARAMETER :: MSC = 9 INTEGER, PARAMETER :: MBAND = 2 REAL :: ALBSAT(MSC,MBAND) REAL :: ALBDRY(MSC,MBAND) REAL :: ALBICE(MBAND) REAL :: ALBLAK(MBAND) REAL :: OMEGAS(MBAND) REAL :: BETADS REAL :: BETAIS REAL :: EG(2) DATA(ALBSAT(I,1),I=1,8)/0.15,0.11,0.10,0.09,0.08,0.07,0.06,0.05/ DATA(ALBSAT(I,2),I=1,8)/0.30,0.22,0.20,0.18,0.16,0.14,0.12,0.10/ DATA(ALBDRY(I,1),I=1,8)/0.27,0.22,0.20,0.18,0.16,0.14,0.12,0.10/ DATA(ALBDRY(I,2),I=1,8)/0.54,0.44,0.40,0.36,0.32,0.28,0.24,0.20/ DATA (ALBICE(I),I=1,MBAND) /0.80, 0.55/ DATA (ALBLAK(I),I=1,MBAND) /0.60, 0.40/ DATA (OMEGAS(I),I=1,MBAND) /0.8, 0.4/ DATA BETADS, BETAIS /0.5, 0.5/ DATA EG /0.97, 0.98/ END MODULE NOAHMP_RAD_PARAMETERS MODULE NOAHMP_ROUTINES USE NOAHMP_GLOBALS IMPLICIT NONE public :: noahmp_options public :: NOAHMP_SFLX public :: REDPRM public :: FRH2O private :: ATM private :: PHENOLOGY private :: ENERGY private :: THERMOPROP private :: CSNOW private :: TDFCND private :: RADIATION private :: ALBEDO private :: SNOW_AGE private :: SNOWALB_BATS private :: SNOWALB_CLASS private :: GROUNDALB private :: TWOSTREAM private :: SURRAD private :: VEGE_FLUX private :: SFCDIF1 private :: SFCDIF2 private :: STOMATA private :: CANRES private :: ESAT private :: RAGRB private :: BARE_FLUX private :: TSNOSOI private :: HRT private :: HSTEP private :: ROSR12 private :: PHASECHANGE private :: WATER private :: CANWATER private :: SNOWWATER private :: SNOWFALL private :: COMBINE private :: DIVIDE private :: COMBO private :: COMPACT private :: SNOWH2O private :: SOILWATER private :: ZWTEQ private :: INFIL private :: SRT private :: WDFCND1 private :: WDFCND2 private :: SSTEP private :: GROUNDWATER private :: CARBON private :: CO2FLUX private :: ERROR contains SUBROUTINE NOAHMP_SFLX (& & ILOC , JLOC , LAT , YEARLEN , JULIAN , COSZ , & & DT , DX , DZ8W , NSOIL , ZSOIL , NSNOW , & & SHDFAC , SHDMAX , VEGTYP , ISURBAN , ICE , IST , & & ISC , & & IZ0TLND , & & SFCTMP , SFCPRS , PSFC , UU , VV , Q2 , & & QC , SOLDN , LWDN , PRCP , TBOT , CO2AIR , & & O2AIR , FOLN , FICEOLD , PBLH , & & ZLVL , ALBOLD , SNEQVO , & & STC , SH2O , SMC , TAH , EAH , FWET , & & CANLIQ , CANICE , TV , TG , QSFC , QSNOW , & & ISNOW , ZSNSO , SNOWH , SNEQV , SNICE , SNLIQ , & & ZWT , WA , WT , WSLAKE , LFMASS , RTMASS , & & STMASS , WOOD , STBLCP , FASTCP , LAI , SAI , & & CM , CH , CHSTAR , & & FSA , FSR , FIRA , FSH , SSOIL , FCEV , & & FGEV , FCTR , ECAN , ETRAN , EDIR , TRAD , & & TS , TGB , TGV , T2MV , T2MB , TSTAR , & & Q1 , Q2V , Q2B , RUNSRF , RUNSUB , APAR , & & PSN , SAV , SAG , FSNO , NEE , GPP , & & NPP , FVEG , ALBEDO , QMELT , PONDING , PONDING1, & & PONDING2, RSSUN , RSSHA , BGAP , WGAP , GAP , & & ERRWAT , CHV , CHB , EMISSI) USE NOAHMP_VEG_PARAMETERS USE NOAHMP_RAD_PARAMETERS implicit none INTEGER , INTENT(IN) :: ICE INTEGER , INTENT(IN) :: IST INTEGER , INTENT(IN) :: VEGTYP INTEGER , INTENT(IN) :: ISC INTEGER , INTENT(IN) :: NSNOW INTEGER , INTENT(IN) :: NSOIL INTEGER , INTENT(IN) :: ILOC INTEGER , INTENT(IN) :: JLOC REAL , INTENT(IN) :: DT REAL, DIMENSION( 1:NSOIL), INTENT(IN) :: ZSOIL REAL , INTENT(IN) :: Q2 REAL , INTENT(IN) :: SFCTMP REAL , INTENT(IN) :: UU REAL , INTENT(IN) :: VV REAL , INTENT(IN) :: SOLDN REAL , INTENT(IN) :: PRCP REAL , INTENT(IN) :: LWDN REAL , INTENT(IN) :: SFCPRS REAL , INTENT(INOUT) :: ZLVL REAL , INTENT(IN) :: COSZ REAL , INTENT(IN) :: TBOT REAL , INTENT(IN) :: FOLN REAL , INTENT(IN) :: SHDFAC INTEGER , INTENT(IN) :: YEARLEN REAL , INTENT(IN) :: JULIAN REAL , INTENT(IN) :: LAT REAL, DIMENSION(-NSNOW+1: 0), INTENT(IN) :: FICEOLD INTEGER , INTENT(IN) :: ISURBAN INTEGER , INTENT(IN) :: IZ0TLND REAL , INTENT(IN) :: QC REAL , INTENT(IN) :: PBLH REAL , INTENT(INOUT) :: QSFC REAL , INTENT(IN) :: PSFC REAL , INTENT(IN) :: DZ8W REAL , INTENT(IN) :: DX REAL , INTENT(IN) :: SHDMAX REAL , INTENT(INOUT) :: QSNOW REAL , INTENT(INOUT) :: FWET REAL , INTENT(INOUT) :: SNEQVO REAL , INTENT(INOUT) :: EAH REAL , INTENT(INOUT) :: TAH REAL , INTENT(INOUT) :: ALBOLD REAL , INTENT(INOUT) :: CM REAL , INTENT(INOUT) :: CH INTEGER , INTENT(INOUT) :: ISNOW REAL , INTENT(INOUT) :: CANLIQ REAL , INTENT(INOUT) :: CANICE REAL , INTENT(INOUT) :: SNEQV REAL, DIMENSION( 1:NSOIL), INTENT(INOUT) :: SMC REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(INOUT) :: ZSNSO REAL , INTENT(INOUT) :: SNOWH REAL, DIMENSION(-NSNOW+1: 0), INTENT(INOUT) :: SNICE REAL, DIMENSION(-NSNOW+1: 0), INTENT(INOUT) :: SNLIQ REAL , INTENT(INOUT) :: TV REAL , INTENT(INOUT) :: TG REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(INOUT) :: STC REAL, DIMENSION( 1:NSOIL), INTENT(INOUT) :: SH2O REAL , INTENT(INOUT) :: ZWT REAL , INTENT(INOUT) :: WA REAL , INTENT(INOUT) :: WT REAL , INTENT(INOUT) :: WSLAKE REAL , INTENT(OUT) :: FSA REAL , INTENT(OUT) :: FSR REAL , INTENT(OUT) :: FIRA REAL , INTENT(OUT) :: FSH REAL , INTENT(OUT) :: FCEV REAL , INTENT(OUT) :: FGEV REAL , INTENT(OUT) :: FCTR REAL , INTENT(OUT) :: SSOIL REAL , INTENT(OUT) :: TRAD REAL , INTENT(OUT) :: TS REAL , INTENT(OUT) :: ECAN REAL , INTENT(OUT) :: ETRAN REAL , INTENT(OUT) :: EDIR REAL , INTENT(OUT) :: RUNSRF REAL , INTENT(OUT) :: RUNSUB REAL , INTENT(OUT) :: PSN REAL , INTENT(OUT) :: APAR REAL , INTENT(OUT) :: SAV REAL , INTENT(OUT) :: SAG REAL , INTENT(OUT) :: FSNO REAL , INTENT(OUT) :: FVEG REAL , INTENT(OUT) :: ALBEDO REAL , INTENT(OUT) :: ERRWAT REAL , INTENT(OUT) :: QMELT REAL , INTENT(OUT) :: PONDING REAL , INTENT(OUT) :: PONDING1 REAL , INTENT(OUT) :: PONDING2 REAL , INTENT(OUT) :: CHSTAR REAL , INTENT(OUT) :: TSTAR REAL , INTENT(OUT) :: T2MV REAL , INTENT(OUT) :: T2MB REAL, INTENT(OUT) :: RSSUN REAL, INTENT(OUT) :: RSSHA REAL, INTENT(OUT) :: BGAP REAL, INTENT(OUT) :: WGAP REAL, INTENT(OUT) :: GAP REAL, INTENT(OUT) :: TGV REAL, INTENT(OUT) :: TGB REAL, INTENT(OUT) :: Q1 REAL, INTENT(OUT) :: EMISSI INTEGER :: IZ INTEGER, DIMENSION(-NSNOW+1:NSOIL) :: IMELT REAL :: CMC REAL :: TAUX REAL :: TAUY REAL :: RHOAIR REAL, DIMENSION(-NSNOW+1:NSOIL) :: DZSNSO REAL :: THAIR REAL :: QAIR REAL :: EAIR REAL, DIMENSION( 1: 2) :: SOLAD REAL, DIMENSION( 1: 2) :: SOLAI REAL :: QPRECC REAL :: QPRECL REAL :: IGS REAL :: ELAI REAL :: ESAI REAL :: BEVAP REAL, DIMENSION( 1:NSOIL) :: BTRANI REAL :: BTRAN REAL :: HTOP REAL :: QIN REAL :: QDIS REAL, DIMENSION( 1:NSOIL) :: SICE REAL, DIMENSION(-NSNOW+1: 0) :: SNICEV REAL, DIMENSION(-NSNOW+1: 0) :: SNLIQV REAL, DIMENSION(-NSNOW+1: 0) :: EPORE REAL :: TOTSC REAL :: TOTLB REAL :: T2M REAL :: QDEW REAL :: QVAP REAL :: LATHEA REAL :: SWDOWN REAL :: BEG_WB REAL :: FSRV REAL :: FSRG REAL,INTENT(OUT) :: Q2V REAL,INTENT(OUT) :: Q2B REAL :: Q2E REAL :: QFX REAL,INTENT(OUT) :: CHV REAL,INTENT(OUT) :: CHB REAL , INTENT(IN) :: CO2AIR REAL , INTENT(IN) :: O2AIR REAL , INTENT(INOUT) :: LFMASS REAL , INTENT(INOUT) :: RTMASS REAL , INTENT(INOUT) :: STMASS REAL , INTENT(INOUT) :: WOOD REAL , INTENT(INOUT) :: STBLCP REAL , INTENT(INOUT) :: FASTCP REAL , INTENT(INOUT) :: LAI REAL , INTENT(INOUT) :: SAI REAL , INTENT(OUT) :: NEE REAL , INTENT(OUT) :: GPP REAL , INTENT(OUT) :: NPP REAL :: AUTORS REAL :: HETERS REAL :: TROOT nee = 0.0 npp = 0.0 gpp = 0.0 CALL ATM (SFCPRS ,SFCTMP ,Q2 ,PRCP ,SOLDN ,COSZ ,THAIR , & QAIR ,EAIR ,RHOAIR ,QPRECC ,QPRECL ,SOLAD ,SOLAI , & SWDOWN ) DO IZ = ISNOW+1, NSOIL IF(IZ == ISNOW+1) THEN DZSNSO(IZ) = - ZSNSO(IZ) ELSE DZSNSO(IZ) = ZSNSO(IZ-1) - ZSNSO(IZ) END IF END DO TROOT = 0. DO IZ=1,NROOT TROOT = TROOT + STC(IZ)*DZSNSO(IZ)/(-ZSOIL(NROOT)) ENDDO IF(IST == 1) THEN BEG_WB = CANLIQ + CANICE + SNEQV + WA DO IZ = 1,NSOIL BEG_WB = BEG_WB + SMC(IZ) * DZSNSO(IZ) * 1000. END DO END IF CALL PHENOLOGY (VEGTYP , SNOWH , TV , LAT , YEARLEN , JULIAN , & LAI , SAI , TROOT , HTOP , ELAI , ESAI ,IGS) IF(DVEG == 1) THEN FVEG = SHDFAC IF(FVEG <= 0.05) FVEG = 0.05 ELSE IF (DVEG == 2 .or. DVEG == 3) THEN FVEG = 1.-EXP(-0.52*(LAI+SAI)) IF(FVEG <= 0.05) FVEG = 0.05 ELSE IF (DVEG == 4) THEN FVEG = SHDMAX IF(FVEG <= 0.05) FVEG = 0.05 ELSE WRITE(*,*) "-------- FATAL CALLED IN SFLX -----------" CALL wrf_error_fatal3("",829,& "Namelist parameter DVEG unknown") ENDIF CALL ENERGY (ICE ,VEGTYP ,IST ,ISC ,NSNOW ,NSOIL , & ISNOW ,NROOT ,DT ,RHOAIR ,SFCPRS ,QAIR , & SFCTMP ,THAIR ,LWDN ,UU ,VV ,ZLVL , & CO2AIR ,O2AIR ,SOLAD ,SOLAI ,COSZ ,IGS , & EAIR ,HTOP ,TBOT ,ZBOT ,ZSNSO ,ZSOIL , & ELAI ,ESAI ,CSOIL ,FWET ,FOLN , & FVEG , & QSNOW ,DZSNSO ,LAT ,CANLIQ ,CANICE ,iloc, jloc , & IMELT ,SNICEV ,SNLIQV ,EPORE ,T2M ,FSNO , & SAV ,SAG ,QMELT ,FSA ,FSR ,TAUX , & TAUY ,FIRA ,FSH ,FCEV ,FGEV ,FCTR , & TRAD ,PSN ,APAR ,SSOIL ,BTRANI ,BTRAN , & PONDING,TS ,LATHEA , & TV ,TG ,STC ,SNOWH ,EAH ,TAH , & SNEQVO ,SNEQV ,SH2O ,SMC ,SNICE ,SNLIQ , & ALBOLD ,CM ,CH ,DX ,DZ8W ,Q2 , & QC ,PBLH ,QSFC ,PSFC ,ISURBAN,IZ0TLND, & CHSTAR ,TSTAR ,T2MV ,T2MB ,FSRV , & FSRG ,RSSUN ,RSSHA ,BGAP ,WGAP, GAP,TGV,TGB,& Q1 ,Q2V ,Q2B ,Q2E ,CHV ,CHB , & EMISSI) SICE(:) = MAX(0.0, SMC(:) - SH2O(:)) SNEQVO = SNEQV QVAP = MAX( FGEV/LATHEA, 0.) QDEW = ABS( MIN(FGEV/LATHEA, 0.)) EDIR = QVAP - QDEW CALL WATER (VEGTYP ,NSNOW ,NSOIL ,IMELT ,DT ,UU , & VV ,FCEV ,FCTR ,QPRECC ,QPRECL ,ELAI , & ESAI ,SFCTMP ,QVAP ,QDEW ,ZSOIL ,BTRANI , & FICEOLD,PONDING,TG ,IST ,FVEG ,iloc,jloc , & ISNOW ,CANLIQ ,CANICE ,TV ,SNOWH ,SNEQV , & SNICE ,SNLIQ ,STC ,ZSNSO ,SH2O ,SMC , & SICE ,ZWT ,WA ,WT ,DZSNSO ,WSLAKE , & CMC ,ECAN ,ETRAN ,FWET ,RUNSRF ,RUNSUB , & QIN ,QDIS ,QSNOW ,PONDING1 ,PONDING2,& ISURBAN) IF (DVEG == 2) THEN CALL CARBON (NSNOW ,NSOIL ,VEGTYP ,NROOT ,DT ,ZSOIL , & DZSNSO ,STC ,SMC ,TV ,TG ,PSN , & FOLN ,SMCMAX ,BTRAN ,APAR ,FVEG ,IGS , & TROOT ,IST ,LAT ,iloc ,jloc , & LFMASS ,RTMASS ,STMASS ,WOOD ,STBLCP ,FASTCP , & GPP ,NPP ,NEE ,AUTORS ,HETERS ,TOTSC , & TOTLB ,LAI ,SAI ) END IF CALL ERROR (SWDOWN ,FSA ,FSR ,FIRA ,FSH ,FCEV , & FGEV ,FCTR ,SSOIL ,BEG_WB ,CANLIQ ,CANICE , & SNEQV ,WA ,SMC ,DZSNSO ,PRCP ,ECAN , & ETRAN ,EDIR ,RUNSRF ,RUNSUB ,DT ,NSOIL , & NSNOW ,IST ,ERRWAT ,ILOC , JLOC ,FVEG , & SAV ,SAG ,FSRV ,FSRG) QFX = ETRAN + ECAN + EDIR IF ( VEGTYP == ISURBAN ) THEN QSFC = (QFX/RHOAIR*CHSTAR) + QAIR Q2B = QSFC END IF IF(SWDOWN.NE.0.) THEN ALBEDO = FSR / SWDOWN ELSE ALBEDO = -999.9 END IF END SUBROUTINE NOAHMP_SFLX SUBROUTINE ATM (SFCPRS ,SFCTMP ,Q2 ,PRCP ,SOLDN ,COSZ ,THAIR , & QAIR ,EAIR ,RHOAIR ,QPRECC ,QPRECL ,SOLAD ,SOLAI , & SWDOWN ) IMPLICIT NONE REAL , INTENT(IN) :: SFCPRS REAL , INTENT(IN) :: SFCTMP REAL , INTENT(IN) :: Q2 REAL , INTENT(IN) :: SOLDN REAL , INTENT(IN) :: PRCP REAL , INTENT(IN) :: COSZ REAL , INTENT(OUT) :: THAIR REAL , INTENT(OUT) :: QAIR REAL , INTENT(OUT) :: EAIR REAL, DIMENSION( 1: 2), INTENT(OUT) :: SOLAD REAL, DIMENSION( 1: 2), INTENT(OUT) :: SOLAI REAL , INTENT(OUT) :: QPRECC REAL , INTENT(OUT) :: QPRECL REAL , INTENT(OUT) :: RHOAIR REAL , INTENT(OUT) :: SWDOWN REAL :: PAIR PAIR = SFCPRS THAIR = SFCTMP * (SFCPRS/PAIR)**(RAIR/CPAIR) QAIR = Q2 / (1.0+Q2) EAIR = QAIR*SFCPRS / (0.622+0.378*QAIR) RHOAIR = (SFCPRS-0.378*EAIR) / (RAIR*SFCTMP) QPRECC = 0.10 * PRCP QPRECL = 0.90 * PRCP IF(COSZ <= 0.) THEN SWDOWN = 0. ELSE SWDOWN = SOLDN END IF SOLAD(1) = SWDOWN*0.7*0.5 SOLAD(2) = SWDOWN*0.7*0.5 SOLAI(1) = SWDOWN*0.3*0.5 SOLAI(2) = SWDOWN*0.3*0.5 END SUBROUTINE ATM SUBROUTINE PHENOLOGY (VEGTYP , SNOWH , TV , LAT , YEARLEN , JULIAN , & LAI , SAI , TROOT , HTOP , ELAI , ESAI , IGS) USE NOAHMP_VEG_PARAMETERS IMPLICIT NONE INTEGER , INTENT(IN ) :: VEGTYP REAL , INTENT(IN ) :: SNOWH REAL , INTENT(IN ) :: TV REAL , INTENT(IN ) :: LAT INTEGER , INTENT(IN ) :: YEARLEN REAL , INTENT(IN ) :: JULIAN real , INTENT(IN ) :: TROOT REAL , INTENT(INOUT) :: LAI REAL , INTENT(INOUT) :: SAI REAL , INTENT(OUT ) :: HTOP REAL , INTENT(OUT ) :: ELAI REAL , INTENT(OUT ) :: ESAI REAL , INTENT(OUT ) :: IGS REAL :: DB REAL :: FB REAL :: SNOWHC INTEGER :: K INTEGER :: IT1,IT2 REAL :: DAY REAL :: WT1,WT2 REAL :: T IF ( DVEG == 1 .or. DVEG == 3 .or. DVEG == 4 ) THEN IF (LAT >= 0.) THEN DAY = JULIAN ELSE DAY = MOD ( JULIAN + ( 0.5 * YEARLEN ) , REAL(YEARLEN) ) ENDIF T = 12. * DAY / REAL(YEARLEN) IT1 = T + 0.5 IT2 = IT1 + 1 WT1 = (IT1+0.5) - T WT2 = 1.-WT1 IF (IT1 .LT. 1) IT1 = 12 IF (IT2 .GT. 12) IT2 = 1 LAI = WT1*LAIM(VEGTYP,IT1) + WT2*LAIM(VEGTYP,IT2) SAI = WT1*SAIM(VEGTYP,IT1) + WT2*SAIM(VEGTYP,IT2) ENDIF IF ( ( VEGTYP == ISWATER ) .OR. ( VEGTYP == ISBARREN ) .OR. ( VEGTYP == ISSNOW ) ) THEN LAI = 0. SAI = 0. ENDIF DB = MIN( MAX(SNOWH - HVB(VEGTYP),0.), HVT(VEGTYP)-HVB(VEGTYP) ) FB = DB / MAX(1.E-06,HVT(VEGTYP)-HVB(VEGTYP)) IF(HVT(VEGTYP)> 0. .AND. HVT(VEGTYP) <= 0.5) THEN SNOWHC = HVT(VEGTYP)*EXP(-SNOWH/0.1) FB = MIN(SNOWH,SNOWHC)/SNOWHC ENDIF ELAI = LAI*(1.-FB) ESAI = SAI*(1.-FB) IF (TV .GT. TMIN(VEGTYP)) THEN IGS = 1. ELSE IGS = 0. ENDIF HTOP = HVT(VEGTYP) END SUBROUTINE PHENOLOGY SUBROUTINE ERROR (SWDOWN ,FSA ,FSR ,FIRA ,FSH ,FCEV , & FGEV ,FCTR ,SSOIL ,BEG_WB ,CANLIQ ,CANICE , & SNEQV ,WA ,SMC ,DZSNSO ,PRCP ,ECAN , & ETRAN ,EDIR ,RUNSRF ,RUNSUB ,DT ,NSOIL , & NSNOW ,IST ,ERRWAT, ILOC ,JLOC ,FVEG , & SAV ,SAG ,FSRV ,FSRG) IMPLICIT NONE INTEGER , INTENT(IN) :: NSNOW INTEGER , INTENT(IN) :: NSOIL INTEGER , INTENT(IN) :: IST INTEGER , INTENT(IN) :: ILOC INTEGER , INTENT(IN) :: JLOC REAL , INTENT(IN) :: SWDOWN REAL , INTENT(IN) :: FSA REAL , INTENT(IN) :: FSR REAL , INTENT(IN) :: FIRA REAL , INTENT(IN) :: FSH REAL , INTENT(IN) :: FCEV REAL , INTENT(IN) :: FGEV REAL , INTENT(IN) :: FCTR REAL , INTENT(IN) :: SSOIL REAL , INTENT(IN) :: FVEG REAL , INTENT(IN) :: SAV REAL , INTENT(IN) :: SAG REAL , INTENT(IN) :: FSRV REAL , INTENT(IN) :: FSRG REAL , INTENT(IN) :: PRCP REAL , INTENT(IN) :: ECAN REAL , INTENT(IN) :: ETRAN REAL , INTENT(IN) :: EDIR REAL , INTENT(IN) :: RUNSRF REAL , INTENT(IN) :: RUNSUB REAL , INTENT(IN) :: CANLIQ REAL , INTENT(IN) :: CANICE REAL , INTENT(IN) :: SNEQV REAL, DIMENSION( 1:NSOIL), INTENT(IN) :: SMC REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(IN) :: DZSNSO REAL , INTENT(IN) :: WA REAL , INTENT(IN) :: DT REAL , INTENT(IN) :: BEG_WB REAL , INTENT(OUT) :: ERRWAT INTEGER :: IZ REAL :: END_WB REAL :: ERRENG REAL :: ERRSW REAL :: FSRVG CHARACTER(len=256) :: message ERRSW = SWDOWN - (SAV+SAG + FSRV+FSRG) IF (ERRSW > 0.01) THEN WRITE(*,*) "VEGETATION!" WRITE(*,*) "SWDOWN*FVEG =",SWDOWN*FVEG WRITE(*,*) "FVEG*(SAV+SAG) =",FVEG*SAV + SAG WRITE(*,*) "FVEG*(FSRV +FSRG)=",FVEG*FSRV + FSRG WRITE(*,*) "GROUND!" WRITE(*,*) "(1-.FVEG)*SWDOWN =",(1.-FVEG)*SWDOWN WRITE(*,*) "(1.-FVEG)*SAG =",(1.-FVEG)*SAG WRITE(*,*) "(1.-FVEG)*FSRG=",(1.-FVEG)*FSRG WRITE(*,*) "FSRV =",FSRV WRITE(*,*) "FSRG =",FSRG WRITE(*,*) "FSR =",FSR WRITE(*,*) "SAV =",SAV WRITE(*,*) "SAG =",SAG WRITE(*,*) "FSA =",FSA WRITE(message,*) 'ERRSW =',ERRSW call wrf_message(trim(message)) call wrf_error_fatal3("",1152,& "Stop in Noah-MP") END IF ERRENG = FVEG*SAV+SAG-(FIRA+FSH+FCEV+FGEV+FCTR+SSOIL) IF(ERRENG > 0.01) THEN write(message,*) 'ERRENG =',ERRENG call wrf_message(trim(message)) WRITE(message,'(i6,1x,i6,1x,7F10.4)')ILOC,JLOC,FSA,FIRA,FSH,FCEV,FGEV,FCTR,SSOIL call wrf_message(trim(message)) call wrf_error_fatal3("",1164,& "Energy budget problem in NOAHMP LSM") END IF IF (IST == 1) THEN END_WB = CANLIQ + CANICE + SNEQV + WA DO IZ = 1,NSOIL END_WB = END_WB + SMC(IZ) * DZSNSO(IZ) * 1000. END DO ERRWAT = END_WB-BEG_WB-(PRCP-ECAN-ETRAN-EDIR-RUNSRF-RUNSUB)*DT IF(ABS(ERRWAT) > 0.1) THEN if (ERRWAT > 0) then call wrf_message ('The model is gaining water (ERRWAT is positive)') else call wrf_message('The model is losing water (ERRWAT is negative)') endif write(message, *) 'ERRWAT =',ERRWAT, "kg m{-2} timestep{-1}" call wrf_message(trim(message)) WRITE(message,'(" I J END_WB BEG_WB PRCP ECAN EDIR ETRAN RUNSRF RUNSUB")') call wrf_message(trim(message)) WRITE(message,'(i6,1x,i6,1x,2f15.3,8f11.5)')ILOC,JLOC,END_WB,BEG_WB,PRCP*DT,ECAN*DT,& EDIR*DT,ETRAN*DT,RUNSRF*DT,RUNSUB*DT call wrf_message(trim(message)) call wrf_error_fatal3("",1188,& "Water budget problem in NOAHMP LSM") END IF ELSE ERRWAT = 0.0 ENDIF END SUBROUTINE ERROR SUBROUTINE ENERGY (ICE ,VEGTYP ,IST ,ISC ,NSNOW ,NSOIL , & ISNOW ,NROOT ,DT ,RHOAIR ,SFCPRS ,QAIR , & SFCTMP ,THAIR ,LWDN ,UU ,VV ,ZREF , & CO2AIR ,O2AIR ,SOLAD ,SOLAI ,COSZ ,IGS , & EAIR ,HTOP ,TBOT ,ZBOT ,ZSNSO ,ZSOIL , & ELAI ,ESAI ,CSOIL ,FWET ,FOLN , & FVEG , & QSNOW ,DZSNSO ,LAT ,CANLIQ ,CANICE ,ILOC , JLOC, & IMELT ,SNICEV ,SNLIQV ,EPORE ,T2M ,FSNO , & SAV ,SAG ,QMELT ,FSA ,FSR ,TAUX , & TAUY ,FIRA ,FSH ,FCEV ,FGEV ,FCTR , & TRAD ,PSN ,APAR ,SSOIL ,BTRANI ,BTRAN , & PONDING,TS ,LATHEA , & TV ,TG ,STC ,SNOWH ,EAH ,TAH , & SNEQVO ,SNEQV ,SH2O ,SMC ,SNICE ,SNLIQ , & ALBOLD ,CM ,CH ,DX ,DZ8W ,Q2 , & QC ,PBLH ,QSFC ,PSFC ,ISURBAN,IZ0TLND, & CHSTAR ,TSTAR ,T2MV ,T2MB ,FSRV , & FSRG ,RSSUN ,RSSHA ,BGAP ,WGAP,GAP,TGV,TGB,& Q1 ,Q2V ,Q2B ,Q2E ,CHV ,CHB, EMISSI ) USE NOAHMP_VEG_PARAMETERS USE NOAHMP_RAD_PARAMETERS IMPLICIT NONE integer , INTENT(IN) :: ILOC integer , INTENT(IN) :: JLOC INTEGER , INTENT(IN) :: ICE INTEGER , INTENT(IN) :: VEGTYP INTEGER , INTENT(IN) :: IST INTEGER , INTENT(IN) :: ISC INTEGER , INTENT(IN) :: NSNOW INTEGER , INTENT(IN) :: NSOIL INTEGER , INTENT(IN) :: NROOT INTEGER , INTENT(IN) :: ISNOW REAL , INTENT(IN) :: DT REAL , INTENT(IN) :: QSNOW REAL , INTENT(IN) :: RHOAIR REAL , INTENT(IN) :: EAIR REAL , INTENT(IN) :: SFCPRS REAL , INTENT(IN) :: QAIR REAL , INTENT(IN) :: SFCTMP REAL , INTENT(IN) :: THAIR REAL , INTENT(IN) :: LWDN REAL , INTENT(IN) :: UU REAL , INTENT(IN) :: VV REAL , DIMENSION( 1: 2), INTENT(IN) :: SOLAD REAL , DIMENSION( 1: 2), INTENT(IN) :: SOLAI REAL , INTENT(IN) :: COSZ REAL , INTENT(IN) :: ELAI REAL , INTENT(IN) :: ESAI REAL , INTENT(IN) :: CSOIL REAL , INTENT(IN) :: FWET REAL , INTENT(IN) :: HTOP REAL , INTENT(IN) :: FVEG REAL , INTENT(IN) :: LAT REAL , INTENT(IN) :: CANLIQ REAL , INTENT(IN) :: CANICE REAL , INTENT(IN) :: FOLN REAL , INTENT(IN) :: CO2AIR REAL , INTENT(IN) :: O2AIR REAL , INTENT(IN) :: IGS REAL , INTENT(IN) :: ZREF REAL , INTENT(IN) :: TBOT REAL , INTENT(IN) :: ZBOT REAL , DIMENSION(-NSNOW+1:NSOIL), INTENT(IN) :: ZSNSO REAL , DIMENSION( 1:NSOIL), INTENT(IN) :: ZSOIL REAL , DIMENSION(-NSNOW+1:NSOIL), INTENT(IN) :: DZSNSO INTEGER , INTENT(IN) :: ISURBAN INTEGER , INTENT(IN) :: IZ0TLND REAL , INTENT(IN) :: QC REAL , INTENT(IN) :: PBLH REAL , INTENT(INOUT) :: QSFC REAL , INTENT(IN) :: PSFC REAL , INTENT(IN) :: DX REAL , INTENT(IN) :: DZ8W REAL , INTENT(IN) :: Q2 INTEGER, DIMENSION(-NSNOW+1:NSOIL), INTENT(OUT) :: IMELT REAL , DIMENSION(-NSNOW+1: 0), INTENT(OUT) :: SNICEV REAL , DIMENSION(-NSNOW+1: 0), INTENT(OUT) :: SNLIQV REAL , DIMENSION(-NSNOW+1: 0), INTENT(OUT) :: EPORE REAL , INTENT(OUT) :: FSNO REAL , INTENT(OUT) :: QMELT REAL , INTENT(OUT) :: PONDING REAL , INTENT(OUT) :: SAV REAL , INTENT(OUT) :: SAG REAL , INTENT(OUT) :: FSA REAL , INTENT(OUT) :: FSR REAL , INTENT(OUT) :: TAUX REAL , INTENT(OUT) :: TAUY REAL , INTENT(OUT) :: FIRA REAL , INTENT(OUT) :: FSH REAL , INTENT(OUT) :: FCEV REAL , INTENT(OUT) :: FGEV REAL , INTENT(OUT) :: FCTR REAL , INTENT(OUT) :: TRAD REAL , INTENT(OUT) :: T2M REAL , INTENT(OUT) :: PSN REAL , INTENT(OUT) :: APAR REAL , INTENT(OUT) :: SSOIL REAL , DIMENSION( 1:NSOIL), INTENT(OUT) :: BTRANI REAL , INTENT(OUT) :: BTRAN REAL , INTENT(OUT) :: LATHEA REAL , INTENT(OUT) :: FSRV REAL , INTENT(OUT) :: FSRG REAL, INTENT(OUT) :: RSSUN REAL, INTENT(OUT) :: RSSHA REAL , INTENT(OUT) :: CHSTAR REAL , INTENT(OUT) :: TSTAR REAL , INTENT(OUT) :: T2MV REAL , INTENT(OUT) :: T2MB REAL , INTENT(OUT) :: BGAP REAL , INTENT(OUT) :: WGAP REAL , INTENT(OUT) :: GAP REAL , INTENT(INOUT) :: TS REAL , INTENT(INOUT) :: TV REAL , INTENT(INOUT) :: TG REAL , DIMENSION(-NSNOW+1:NSOIL), INTENT(INOUT) :: STC REAL , INTENT(INOUT) :: SNOWH REAL , INTENT(INOUT) :: SNEQV REAL , INTENT(INOUT) :: SNEQVO REAL , DIMENSION( 1:NSOIL), INTENT(INOUT) :: SH2O REAL , DIMENSION( 1:NSOIL), INTENT(INOUT) :: SMC REAL , DIMENSION(-NSNOW+1: 0), INTENT(INOUT) :: SNICE REAL , DIMENSION(-NSNOW+1: 0), INTENT(INOUT) :: SNLIQ REAL , INTENT(INOUT) :: EAH REAL , INTENT(INOUT) :: TAH REAL , INTENT(INOUT) :: ALBOLD REAL , INTENT(INOUT) :: CM REAL , INTENT(INOUT) :: CH REAL , INTENT(INOUT) :: Q1 REAL, INTENT(OUT) :: EMISSI INTEGER :: IZ LOGICAL :: VEG REAL :: UR REAL :: ZLVL REAL :: FSUN REAL :: RB REAL :: RSURF REAL :: L_RSURF REAL :: D_RSURF REAL :: BEVAP REAL :: MOL REAL :: VAI REAL :: CWP REAL :: ZPD REAL :: Z0M REAL :: ZPDG REAL :: Z0MG REAL :: EMV REAL :: EMG REAL :: FIRE REAL :: LAISUN REAL :: LAISHA REAL :: PSNSUN REAL :: PSNSHA REAL :: PARSUN REAL :: PARSHA REAL, DIMENSION(-NSNOW+1:NSOIL) :: FACT REAL, DIMENSION(-NSNOW+1:NSOIL) :: DF REAL, DIMENSION(-NSNOW+1:NSOIL) :: HCPCT REAL :: BDSNO REAL :: FMELT REAL :: GX REAL, DIMENSION(-NSNOW+1:NSOIL) :: PHI REAL :: GAMMA REAL :: PSI REAL :: RHSUR REAL :: TAUXV REAL :: TAUYV REAL :: IRC REAL :: IRG REAL :: SHC REAL :: SHG REAL,INTENT(OUT) :: Q2V REAL,INTENT(OUT) :: Q2B REAL,INTENT(OUT) :: Q2E REAL :: EVC REAL :: EVG REAL :: TR REAL :: GHV REAL,INTENT(OUT) :: TGV REAL :: CMV REAL,INTENT(OUT) :: CHV REAL :: TAUXB REAL :: TAUYB REAL :: IRB REAL :: SHB REAL :: EVB REAL :: GHB REAL,INTENT(OUT) :: TGB REAL :: CMB REAL,INTENT(OUT) :: CHB REAL :: CAH2 REAL :: EHB2 REAL :: noahmpres REAL, PARAMETER :: MPE = 1.E-6 REAL, PARAMETER :: PSIWLT = -150. REAL, PARAMETER :: Z0 = 0.01 TAUXV = 0. TAUYV = 0. IRC = 0. SHC = 0. IRG = 0. SHG = 0. EVG = 0. EVC = 0. TR = 0. GHV = 0. PSNSUN = 0. PSNSHA = 0. UR = MAX( SQRT(UU**2.+VV**2.), 1. ) VAI = ELAI + ESAI VEG = .FALSE. IF(VAI > 0.) VEG = .TRUE. FSNO = 0. IF(SNOWH.GT.0.) THEN BDSNO = SNEQV / SNOWH FMELT = (BDSNO/100.)**M FSNO = TANH( SNOWH /(2.5* Z0 * FMELT)) ENDIF IF(IST == 2) THEN IF(TG .LE. TFRZ) THEN Z0MG = 0.01 * (1.0-FSNO) + FSNO * Z0SNO ELSE Z0MG = 0.01 END IF ELSE Z0MG = Z0 * (1.0-FSNO) + FSNO * Z0SNO END IF ZPDG = SNOWH IF(VEG) THEN Z0M = Z0MVT(VEGTYP) ZPD = 0.65 * HTOP IF(SNOWH.GT.ZPD) ZPD = SNOWH ELSE Z0M = Z0MG ZPD = ZPDG END IF ZLVL = MAX(ZPD,HTOP) + ZREF IF(ZPDG >= ZLVL) ZLVL = ZPDG + ZREF CWP = CWPVT(VEGTYP) CALL THERMOPROP (NSOIL ,NSNOW ,ISNOW ,IST ,DZSNSO , & DT ,SNOWH ,SNICE ,SNLIQ ,CSOIL , & SMC ,SH2O ,TG ,STC ,UR , & LAT ,Z0M ,ZLVL ,VEGTYP ,ISURBAN , & DF ,HCPCT ,SNICEV ,SNLIQV ,EPORE , & FACT ) CALL RADIATION (VEGTYP ,IST ,ISC ,ICE ,NSOIL , & SNEQVO ,SNEQV ,DT ,COSZ ,SNOWH , & TG ,TV ,FSNO ,QSNOW ,FWET , & ELAI ,ESAI ,SMC ,SOLAD ,SOLAI , & FVEG ,ILOC ,JLOC , & ALBOLD , & FSUN ,LAISUN ,LAISHA ,PARSUN ,PARSHA , & SAV ,SAG ,FSR ,FSA ,FSRV , & FSRG ,BGAP ,WGAP ,GAP) EMV = 1. - EXP(-(ELAI+ESAI)/1.0) IF (ICE == 1) THEN EMG = 0.98*(1.-FSNO) + 1.0*FSNO ELSE EMG = EG(IST)*(1.-FSNO) + 1.0*FSNO END IF BTRAN = 0. IF(IST ==1 ) THEN DO IZ = 1, NROOT IF(OPT_BTR == 1) then GX = (SH2O(IZ)-SMCWLT) / (SMCREF-SMCWLT) END IF IF(OPT_BTR == 2) then PSI = MAX(PSIWLT,-PSISAT*(MAX(0.01,SH2O(IZ))/SMCMAX)**(-BEXP) ) GX = (1.-PSI/PSIWLT)/(1.+PSISAT/PSIWLT) END IF IF(OPT_BTR == 3) then PSI = MAX(PSIWLT,-PSISAT*(MAX(0.01,SH2O(IZ))/SMCMAX)**(-BEXP) ) GX = 1.-EXP(-5.8*(LOG(PSIWLT/PSI))) END IF GX = MIN(1.,MAX(0.,GX)) BTRANI(IZ) = MAX(MPE,DZSNSO(IZ) / (-ZSOIL(NROOT)) * GX) BTRAN = BTRAN + BTRANI(IZ) END DO BTRAN = MAX(MPE,BTRAN) BTRANI(1:NROOT) = BTRANI(1:NROOT)/BTRAN END IF BEVAP = MAX(0.0,SH2O(1)/SMCMAX) IF(IST == 2) THEN RSURF = 1. RHSUR = 1.0 ELSE L_RSURF = (-ZSOIL(1)) * ( exp ( (1.0 - MIN(1.0,SH2O(1)/SMCMAX)) ** 5 ) - 1.0 ) / ( 2.71828 - 1.0 ) D_RSURF = 2.2E-5 * SMCMAX * SMCMAX * ( 1.0 - SMCWLT / SMCMAX ) ** (2.0+3.0/BEXP) RSURF = L_RSURF / D_RSURF IF(SH2O(1) < 0.01 .and. SNOWH == 0.) RSURF = 1.E6 PSI = -PSISAT*(MAX(0.01,SH2O(1))/SMCMAX)**(-BEXP) RHSUR = FSNO + (1.-FSNO) * EXP(PSI*GRAV/(RW*TG)) END IF IF (VEGTYP == ISURBAN .and. SNOWH == 0. ) THEN RSURF = 1.E6 ENDIF IF (SFCTMP .GT. TFRZ) THEN LATHEA = HVAP ELSE LATHEA = HSUB END IF GAMMA = CPAIR*SFCPRS/(0.622*LATHEA) IF (VEG) THEN TGV = TG CMV = CM CHV = CH CALL VEGE_FLUX (NSNOW ,NSOIL ,ISNOW ,VEGTYP ,VEG , & DT ,SAV ,SAG ,LWDN ,UR , & UU ,VV ,SFCTMP ,THAIR ,QAIR , & EAIR ,RHOAIR ,SNOWH ,VAI ,GAMMA , & FWET ,LAISUN ,LAISHA ,CWP ,DZSNSO , & HTOP ,ZLVL ,ZPD ,Z0M ,FVEG , & Z0MG ,EMV ,EMG ,CANLIQ , & CANICE ,STC ,DF ,RSSUN ,RSSHA , & RSURF ,LATHEA ,PARSUN ,PARSHA ,IGS , & FOLN ,CO2AIR ,O2AIR ,BTRAN ,SFCPRS , & RHSUR ,ILOC ,JLOC ,Q2 , & EAH ,TAH ,TV ,TGV ,CMV , & CHV ,DX ,DZ8W , & TAUXV ,TAUYV ,IRG ,IRC ,SHG , & SHC ,EVG ,EVC ,TR ,GHV , & T2MV ,PSNSUN ,PSNSHA , & QC ,PBLH ,QSFC ,PSFC ,ISURBAN , & IZ0TLND ,Q2V ,CAH2) END IF TGB = TG CMB = CM CHB = CH CALL BARE_FLUX (NSNOW ,NSOIL ,ISNOW ,DT ,SAG , & LWDN ,UR ,UU ,VV ,SFCTMP , & THAIR ,QAIR ,EAIR ,RHOAIR ,SNOWH , & DZSNSO ,ZLVL ,ZPDG ,Z0MG , & EMG ,STC ,DF ,RSURF ,LATHEA , & GAMMA ,RHSUR ,ILOC ,JLOC ,Q2 , & TGB ,CMB ,CHB , & TAUXB ,TAUYB ,IRB ,SHB ,EVB , & GHB ,T2MB ,DX ,DZ8W ,VEGTYP , & QC ,PBLH ,QSFC ,PSFC ,ISURBAN , & IZ0TLND ,SFCPRS ,Q2B, EHB2) IF (VEG) THEN TAUX = FVEG * TAUXV + (1.0 - FVEG) * TAUXB TAUY = FVEG * TAUYV + (1.0 - FVEG) * TAUYB FIRA = FVEG * IRG + (1.0 - FVEG) * IRB + FVEG * IRC FSH = FVEG * SHG + (1.0 - FVEG) * SHB + FVEG * SHC FGEV = FVEG * EVG + (1.0 - FVEG) * EVB SSOIL = FVEG * GHV + (1.0 - FVEG) * GHB FCEV = FVEG * EVC FCTR = FVEG * TR TG = FVEG * TGV + (1.0 - FVEG) * TGB T2M = FVEG * T2MV + (1.0 - FVEG) * T2MB TS = FVEG * TV + (1.0 - FVEG) * TGB CM = FVEG * CMV + (1.0 - FVEG) * CMB CH = FVEG * CHV + (1.0 - FVEG) * CHB Q1 = FVEG * (EAH*0.622/(SFCPRS - 0.378*EAH)) + (1.0 - FVEG)*QSFC Q2E = FVEG * Q2V + (1.0 - FVEG) * Q2B ELSE TAUX = TAUXB TAUY = TAUYB FIRA = IRB FSH = SHB FGEV = EVB SSOIL = GHB TG = TGB T2M = T2MB FCEV = 0. FCTR = 0. TS = TG CM = CMB CH = CHB Q1 = QSFC Q2E = Q2B END IF FIRE = LWDN + FIRA IF(FIRE <=0.) THEN WRITE(6,*) 'emitted longwave <0; skin T may be wrong due to inconsistent' WRITE(6,*) 'input of SHDFAC with LAI' WRITE(6,*) ILOC, JLOC, 'SHDFAC=',FVEG,'VAI=',VAI,'TV=',TV,'TG=',TG WRITE(6,*) 'LWDN=',LWDN,'FIRA=',FIRA,'SNOWH=',SNOWH call wrf_error_fatal3("",1723,& "STOP in Noah-MP") END IF EMISSI = FVEG * ( EMG*(1-EMV) + EMV + EMV*(1-EMV)*(1-EMG) ) + & (1-FVEG) * EMG TRAD = ( ( FIRE - (1-EMISSI)*LWDN ) / (EMISSI*SB) ) ** 0.25 APAR = PARSUN*LAISUN + PARSHA*LAISHA PSN = PSNSUN*LAISUN + PSNSHA*LAISHA CALL EPARM(ILOC ,JLOC ,TAH ,TGB ,FVEG ,& CHV ,CHB ,VEG ,CHSTAR ,TSTAR) CALL TSNOSOI (ICE ,NSOIL ,NSNOW ,ISNOW ,IST , & TBOT ,ZSNSO ,SSOIL ,DF ,HCPCT , & ZBOT ,SAG ,DT ,SNOWH ,DZSNSO , & TG ,ILOC ,JLOC , & STC ) IF(OPT_STC == 2) THEN IF (SNOWH > 0.05 .AND. TG > TFRZ) THEN TGV = TFRZ TGB = TFRZ IF (VEG) THEN TG = FVEG * TGV + (1.0 - FVEG) * TGB TS = FVEG * TV + (1.0 - FVEG) * TGB ELSE TG = TGB TS = TGB END IF END IF END IF CALL PHASECHANGE (NSNOW ,NSOIL ,ISNOW ,DT ,FACT , & DZSNSO ,HCPCT ,IST ,ILOC ,JLOC , & STC ,SNICE ,SNLIQ ,SNEQV ,SNOWH , & SMC ,SH2O , & QMELT ,IMELT ,PONDING ) END SUBROUTINE ENERGY SUBROUTINE THERMOPROP (NSOIL ,NSNOW ,ISNOW ,IST ,DZSNSO , & DT ,SNOWH ,SNICE ,SNLIQ ,CSOIL , & SMC ,SH2O ,TG ,STC ,UR , & LAT ,Z0M ,ZLVL ,VEGTYP ,ISURBAN , & DF ,HCPCT ,SNICEV ,SNLIQV ,EPORE , & FACT ) IMPLICIT NONE INTEGER , INTENT(IN) :: NSOIL INTEGER , INTENT(IN) :: NSNOW INTEGER , INTENT(IN) :: ISNOW INTEGER , INTENT(IN) :: IST REAL , INTENT(IN) :: DT REAL, DIMENSION(-NSNOW+1: 0), INTENT(IN) :: SNICE REAL, DIMENSION(-NSNOW+1: 0), INTENT(IN) :: SNLIQ REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(IN) :: DZSNSO REAL, DIMENSION( 1:NSOIL), INTENT(IN) :: SMC REAL, DIMENSION( 1:NSOIL), INTENT(IN) :: SH2O REAL , INTENT(IN) :: SNOWH REAL , INTENT(IN) :: CSOIL REAL, INTENT(IN) :: TG REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(IN) :: STC REAL, INTENT(IN) :: UR REAL, INTENT(IN) :: LAT REAL, INTENT(IN) :: Z0M REAL, INTENT(IN) :: ZLVL INTEGER , INTENT(IN) :: VEGTYP INTEGER , INTENT(IN) :: ISURBAN REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(OUT) :: DF REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(OUT) :: HCPCT REAL, DIMENSION(-NSNOW+1: 0), INTENT(OUT) :: SNICEV REAL, DIMENSION(-NSNOW+1: 0), INTENT(OUT) :: SNLIQV REAL, DIMENSION(-NSNOW+1: 0), INTENT(OUT) :: EPORE REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(OUT) :: FACT INTEGER :: IZ REAL, DIMENSION(-NSNOW+1: 0) :: CVSNO REAL, DIMENSION(-NSNOW+1: 0) :: TKSNO REAL, DIMENSION( 1:NSOIL) :: SICE CALL CSNOW (ISNOW ,NSNOW ,NSOIL ,SNICE ,SNLIQ ,DZSNSO , & TKSNO ,CVSNO ,SNICEV ,SNLIQV ,EPORE ) DO IZ = ISNOW+1, 0 DF (IZ) = TKSNO(IZ) HCPCT(IZ) = CVSNO(IZ) END DO DO IZ = 1, NSOIL SICE(IZ) = SMC(IZ) - SH2O(IZ) HCPCT(IZ) = SH2O(IZ)*CWAT + (1.0-SMCMAX)*CSOIL & + (SMCMAX-SMC(IZ))*CPAIR + SICE(IZ)*CICE CALL TDFCND (DF(IZ), SMC(IZ), SH2O(IZ)) END DO IF ( VEGTYP == ISURBAN ) THEN DO IZ = 1,NSOIL DF(IZ) = 3.24 END DO ENDIF IF(IST == 2) THEN DO IZ = 1, NSOIL IF(STC(IZ) > TFRZ) THEN HCPCT(IZ) = CWAT DF(IZ) = TKWAT ELSE HCPCT(IZ) = CICE DF(IZ) = TKICE END IF END DO END IF DO IZ = ISNOW+1,NSOIL FACT(IZ) = DT/(HCPCT(IZ)*DZSNSO(IZ)) END DO IF(ISNOW == 0) THEN DF(1) = (DF(1)*DZSNSO(1)+0.35*SNOWH) / (SNOWH +DZSNSO(1)) ELSE DF(1) = (DF(1)*DZSNSO(1)+DF(0)*DZSNSO(0)) / (DZSNSO(0)+DZSNSO(1)) END IF END SUBROUTINE THERMOPROP SUBROUTINE CSNOW (ISNOW ,NSNOW ,NSOIL ,SNICE ,SNLIQ ,DZSNSO , & TKSNO ,CVSNO ,SNICEV ,SNLIQV ,EPORE ) IMPLICIT NONE INTEGER, INTENT(IN) :: ISNOW INTEGER , INTENT(IN) :: NSNOW INTEGER , INTENT(IN) :: NSOIL REAL, DIMENSION(-NSNOW+1: 0), INTENT(IN) :: SNICE REAL, DIMENSION(-NSNOW+1: 0), INTENT(IN) :: SNLIQ REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(IN) :: DZSNSO REAL, DIMENSION(-NSNOW+1: 0), INTENT(OUT) :: CVSNO REAL, DIMENSION(-NSNOW+1: 0), INTENT(OUT) :: TKSNO REAL, DIMENSION(-NSNOW+1: 0), INTENT(OUT) :: SNICEV REAL, DIMENSION(-NSNOW+1: 0), INTENT(OUT) :: SNLIQV REAL, DIMENSION(-NSNOW+1: 0), INTENT(OUT) :: EPORE INTEGER :: IZ REAL, DIMENSION(-NSNOW+1: 0) :: BDSNOI DO IZ = ISNOW+1, 0 SNICEV(IZ) = MIN(1., SNICE(IZ)/(DZSNSO(IZ)*DENICE) ) EPORE(IZ) = 1. - SNICEV(IZ) SNLIQV(IZ) = MIN(EPORE(IZ),SNLIQ(IZ)/(DZSNSO(IZ)*DENH2O)) ENDDO DO IZ = ISNOW+1, 0 BDSNOI(IZ) = (SNICE(IZ)+SNLIQ(IZ))/DZSNSO(IZ) CVSNO(IZ) = CICE*SNICEV(IZ)+CWAT*SNLIQV(IZ) enddo DO IZ = ISNOW+1, 0 TKSNO(IZ) = 3.2217E-6*BDSNOI(IZ)**2. ENDDO END SUBROUTINE CSNOW SUBROUTINE TDFCND ( DF, SMC, SH2O) IMPLICIT NONE REAL, INTENT(IN) :: SMC REAL, INTENT(IN) :: SH2O REAL, INTENT(OUT) :: DF REAL :: AKE REAL :: GAMMD REAL :: THKDRY REAL :: THKO REAL :: THKQTZ REAL :: THKSAT REAL :: THKS REAL :: THKW REAL :: SATRATIO REAL :: XU REAL :: XUNFROZ SATRATIO = SMC / SMCMAX THKW = 0.57 THKO = 2.0 THKQTZ = 7.7 THKS = (THKQTZ ** QUARTZ)* (THKO ** (1. - QUARTZ)) XUNFROZ = SH2O / SMC XU = XUNFROZ * SMCMAX THKSAT = THKS ** (1. - SMCMAX)* TKICE ** (SMCMAX - XU)* THKW ** & (XU) GAMMD = (1. - SMCMAX)*2700. THKDRY = (0.135* GAMMD+ 64.7)/ (2700. - 0.947* GAMMD) IF ( (SH2O + 0.0005) < SMC ) THEN AKE = SATRATIO ELSE IF ( SATRATIO > 0.1 ) THEN AKE = LOG10 (SATRATIO) + 1.0 ELSE AKE = 0.0 END IF END IF DF = AKE * (THKSAT - THKDRY) + THKDRY end subroutine TDFCND SUBROUTINE RADIATION (VEGTYP ,IST ,ISC ,ICE ,NSOIL , & SNEQVO ,SNEQV ,DT ,COSZ ,SNOWH , & TG ,TV ,FSNO ,QSNOW ,FWET , & ELAI ,ESAI ,SMC ,SOLAD ,SOLAI , & FVEG ,ILOC ,JLOC , & ALBOLD , & FSUN ,LAISUN ,LAISHA ,PARSUN ,PARSHA , & SAV ,SAG ,FSR ,FSA ,FSRV , & FSRG ,BGAP ,WGAP,GAP) IMPLICIT NONE INTEGER, INTENT(IN) :: ILOC INTEGER, INTENT(IN) :: JLOC INTEGER, INTENT(IN) :: VEGTYP INTEGER, INTENT(IN) :: IST INTEGER, INTENT(IN) :: ISC INTEGER, INTENT(IN) :: ICE INTEGER, INTENT(IN) :: NSOIL REAL, INTENT(IN) :: DT REAL, INTENT(IN) :: QSNOW REAL, INTENT(IN) :: SNEQVO REAL, INTENT(IN) :: SNEQV REAL, INTENT(IN) :: SNOWH REAL, INTENT(IN) :: COSZ REAL, INTENT(IN) :: TG REAL, INTENT(IN) :: TV REAL, INTENT(IN) :: ELAI REAL, INTENT(IN) :: ESAI REAL, INTENT(IN) :: FWET REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SMC REAL, DIMENSION(1:2) , INTENT(IN) :: SOLAD REAL, DIMENSION(1:2) , INTENT(IN) :: SOLAI REAL, INTENT(IN) :: FSNO REAL, INTENT(IN) :: FVEG REAL, INTENT(INOUT) :: ALBOLD REAL, INTENT(OUT) :: FSUN REAL, INTENT(OUT) :: LAISUN REAL, INTENT(OUT) :: LAISHA REAL, INTENT(OUT) :: PARSUN REAL, INTENT(OUT) :: PARSHA REAL, INTENT(OUT) :: SAV REAL, INTENT(OUT) :: SAG REAL, INTENT(OUT) :: FSA REAL, INTENT(OUT) :: FSR REAL, INTENT(OUT) :: FSRV REAL, INTENT(OUT) :: FSRG REAL, INTENT(OUT) :: BGAP REAL, INTENT(OUT) :: WGAP REAL, INTENT(OUT) :: GAP REAL :: FAGE REAL, DIMENSION(1:2) :: ALBGRD REAL, DIMENSION(1:2) :: ALBGRI REAL, DIMENSION(1:2) :: ALBD REAL, DIMENSION(1:2) :: ALBI REAL, DIMENSION(1:2) :: FABD REAL, DIMENSION(1:2) :: FABI REAL, DIMENSION(1:2) :: FTDD REAL, DIMENSION(1:2) :: FTID REAL, DIMENSION(1:2) :: FTII REAL, DIMENSION(1:2) :: FREVI REAL, DIMENSION(1:2) :: FREVD REAL, DIMENSION(1:2) :: FREGI REAL, DIMENSION(1:2) :: FREGD REAL :: FSHA REAL :: VAI REAL,PARAMETER :: MPE = 1.E-6 LOGICAL VEG CALL ALBEDO (VEGTYP ,IST ,ISC ,ICE ,NSOIL , & DT ,COSZ ,FAGE ,ELAI ,ESAI , & TG ,TV ,SNOWH ,FSNO ,FWET , & SMC ,SNEQVO ,SNEQV ,QSNOW ,FVEG , & ILOC ,JLOC , & ALBOLD , & ALBGRD ,ALBGRI ,ALBD ,ALBI ,FABD , & FABI ,FTDD ,FTID ,FTII ,FSUN , & FREVI ,FREVD ,FREGD ,FREGI ,BGAP , & WGAP ,GAP) FSHA = 1.-FSUN LAISUN = ELAI*FSUN LAISHA = ELAI*FSHA VAI = ELAI+ ESAI IF (VAI .GT. 0.) THEN VEG = .TRUE. ELSE VEG = .FALSE. END IF CALL SURRAD (MPE ,FSUN ,FSHA ,ELAI ,VAI , & LAISUN ,LAISHA ,SOLAD ,SOLAI ,FABD , & FABI ,FTDD ,FTID ,FTII ,ALBGRD , & ALBGRI ,ALBD ,ALBI ,ILOC ,JLOC , & PARSUN ,PARSHA ,SAV ,SAG ,FSA , & FSR , & FREVI ,FREVD ,FREGD ,FREGI ,FSRV , & FSRG) END SUBROUTINE RADIATION SUBROUTINE ALBEDO (VEGTYP ,IST ,ISC ,ICE ,NSOIL , & DT ,COSZ ,FAGE ,ELAI ,ESAI , & TG ,TV ,SNOWH ,FSNO ,FWET , & SMC ,SNEQVO ,SNEQV ,QSNOW ,FVEG , & ILOC ,JLOC , & ALBOLD , & ALBGRD ,ALBGRI ,ALBD ,ALBI ,FABD , & FABI ,FTDD ,FTID ,FTII ,FSUN , & FREVI ,FREVD ,FREGD ,FREGI ,BGAP , & WGAP ,GAP) USE NOAHMP_VEG_PARAMETERS IMPLICIT NONE INTEGER, INTENT(IN) :: ILOC INTEGER, INTENT(IN) :: JLOC INTEGER, INTENT(IN) :: NSOIL INTEGER, INTENT(IN) :: VEGTYP INTEGER, INTENT(IN) :: IST INTEGER, INTENT(IN) :: ISC INTEGER, INTENT(IN) :: ICE REAL, INTENT(IN) :: DT REAL, INTENT(IN) :: QSNOW REAL, INTENT(IN) :: COSZ REAL, INTENT(IN) :: SNOWH REAL, INTENT(IN) :: TG REAL, INTENT(IN) :: TV REAL, INTENT(IN) :: ELAI REAL, INTENT(IN) :: ESAI REAL, INTENT(IN) :: FSNO REAL, INTENT(IN) :: FWET REAL, INTENT(IN) :: SNEQVO REAL, INTENT(IN) :: SNEQV REAL, INTENT(IN) :: FVEG REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SMC REAL, INTENT(INOUT) :: ALBOLD REAL, DIMENSION(1: 2), INTENT(OUT) :: ALBGRD REAL, DIMENSION(1: 2), INTENT(OUT) :: ALBGRI REAL, DIMENSION(1: 2), INTENT(OUT) :: ALBD REAL, DIMENSION(1: 2), INTENT(OUT) :: ALBI REAL, DIMENSION(1: 2), INTENT(OUT) :: FABD REAL, DIMENSION(1: 2), INTENT(OUT) :: FABI REAL, DIMENSION(1: 2), INTENT(OUT) :: FTDD REAL, DIMENSION(1: 2), INTENT(OUT) :: FTID REAL, DIMENSION(1: 2), INTENT(OUT) :: FTII REAL, INTENT(OUT) :: FSUN REAL, DIMENSION(1: 2), INTENT(OUT) :: FREVD REAL, DIMENSION(1: 2), INTENT(OUT) :: FREVI REAL, DIMENSION(1: 2), INTENT(OUT) :: FREGD REAL, DIMENSION(1: 2), INTENT(OUT) :: FREGI REAL, INTENT(OUT) :: BGAP REAL, INTENT(OUT) :: WGAP REAL, INTENT(OUT) :: GAP REAL :: FAGE REAL :: ALB INTEGER :: IB INTEGER :: NBAND INTEGER :: IC REAL :: WL REAL :: WS REAL :: MPE REAL, DIMENSION(1:2) :: RHO REAL, DIMENSION(1:2) :: TAU REAL, DIMENSION(1:2) :: FTDI REAL, DIMENSION(1:2) :: ALBSND REAL, DIMENSION(1:2) :: ALBSNI REAL :: VAI REAL :: GDIR REAL :: EXT NBAND = 2 MPE = 1.E-06 BGAP = 0. WGAP = 0. GAP = 0. DO IB = 1, NBAND ALBD(IB) = 0. ALBI(IB) = 0. ALBGRD(IB) = 0. ALBGRI(IB) = 0. FABD(IB) = 0. FABI(IB) = 0. FTDD(IB) = 0. FTID(IB) = 0. FTII(IB) = 0. IF (IB.EQ.1) FSUN = 0. END DO IF(COSZ <= 0) GOTO 100 DO IB = 1, NBAND VAI = ELAI + ESAI WL = ELAI / MAX(VAI,MPE) WS = ESAI / MAX(VAI,MPE) RHO(IB) = MAX(RHOL(VEGTYP,IB)*WL+RHOS(VEGTYP,IB)*WS, MPE) TAU(IB) = MAX(TAUL(VEGTYP,IB)*WL+TAUS(VEGTYP,IB)*WS, MPE) END DO CALL SNOW_AGE (DT,TG,SNEQVO,SNEQV,FAGE) IF(OPT_ALB == 1) & CALL SNOWALB_BATS (NBAND, FSNO,COSZ,FAGE,ALBSND,ALBSNI) IF(OPT_ALB == 2) THEN CALL SNOWALB_CLASS (NBAND,QSNOW,DT,ALB,ALBOLD,ALBSND,ALBSNI,ILOC,JLOC) ALBOLD = ALB END IF CALL GROUNDALB (NSOIL ,NBAND ,ICE ,IST ,ISC , & FSNO ,SMC ,ALBSND ,ALBSNI ,COSZ , & TG ,ILOC ,JLOC , & ALBGRD ,ALBGRI ) DO IB = 1, NBAND IC = 0 CALL TWOSTREAM (IB ,IC ,VEGTYP ,COSZ ,VAI , & FWET ,TV ,ALBGRD ,ALBGRI ,RHO , & TAU ,FVEG ,IST ,ILOC ,JLOC , & FABD ,ALBD ,FTDD ,FTID ,GDIR , & FREVD ,FREGD ,BGAP ,WGAP ,GAP) IC = 1 CALL TWOSTREAM (IB ,IC ,VEGTYP ,COSZ ,VAI , & FWET ,TV ,ALBGRD ,ALBGRI ,RHO , & TAU ,FVEG ,IST ,ILOC ,JLOC , & FABI ,ALBI ,FTDI ,FTII ,GDIR , & FREVI ,FREGI ,BGAP ,WGAP ,GAP) END DO EXT = GDIR/COSZ * SQRT(1.-RHO(1)-TAU(1)) FSUN = (1.-EXP(-EXT*VAI)) / MAX(EXT*VAI,MPE) EXT = FSUN IF (EXT .LT. 0.01) THEN WL = 0. ELSE WL = EXT END IF FSUN = WL 100 CONTINUE END SUBROUTINE ALBEDO SUBROUTINE SURRAD (MPE ,FSUN ,FSHA ,ELAI ,VAI , & LAISUN ,LAISHA ,SOLAD ,SOLAI ,FABD , & FABI ,FTDD ,FTID ,FTII ,ALBGRD , & ALBGRI ,ALBD ,ALBI ,ILOC ,JLOC , & PARSUN ,PARSHA ,SAV ,SAG ,FSA , & FSR , & FREVI ,FREVD ,FREGD ,FREGI ,FSRV , & FSRG) IMPLICIT NONE INTEGER, INTENT(IN) :: ILOC INTEGER, INTENT(IN) :: JLOC REAL, INTENT(IN) :: MPE REAL, INTENT(IN) :: FSUN REAL, INTENT(IN) :: FSHA REAL, INTENT(IN) :: ELAI REAL, INTENT(IN) :: VAI REAL, INTENT(IN) :: LAISUN REAL, INTENT(IN) :: LAISHA REAL, DIMENSION(1:2), INTENT(IN) :: SOLAD REAL, DIMENSION(1:2), INTENT(IN) :: SOLAI REAL, DIMENSION(1:2), INTENT(IN) :: FABD REAL, DIMENSION(1:2), INTENT(IN) :: FABI REAL, DIMENSION(1:2), INTENT(IN) :: FTDD REAL, DIMENSION(1:2), INTENT(IN) :: FTID REAL, DIMENSION(1:2), INTENT(IN) :: FTII REAL, DIMENSION(1:2), INTENT(IN) :: ALBGRD REAL, DIMENSION(1:2), INTENT(IN) :: ALBGRI REAL, DIMENSION(1:2), INTENT(IN) :: ALBD REAL, DIMENSION(1:2), INTENT(IN) :: ALBI REAL, DIMENSION(1:2), INTENT(IN) :: FREVD REAL, DIMENSION(1:2), INTENT(IN) :: FREVI REAL, DIMENSION(1:2), INTENT(IN) :: FREGD REAL, DIMENSION(1:2), INTENT(IN) :: FREGI REAL, INTENT(OUT) :: PARSUN REAL, INTENT(OUT) :: PARSHA REAL, INTENT(OUT) :: SAV REAL, INTENT(OUT) :: SAG REAL, INTENT(OUT) :: FSA REAL, INTENT(OUT) :: FSR REAL, INTENT(OUT) :: FSRV REAL, INTENT(OUT) :: FSRG INTEGER :: IB INTEGER :: NBAND REAL :: ABS REAL :: RNIR REAL :: RVIS REAL :: LAIFRA REAL :: TRD REAL :: TRI REAL, DIMENSION(1:2) :: CAD REAL, DIMENSION(1:2) :: CAI NBAND = 2 SAG = 0. SAV = 0. FSA = 0. DO IB = 1, NBAND CAD(IB) = SOLAD(IB)*FABD(IB) CAI(IB) = SOLAI(IB)*FABI(IB) SAV = SAV + CAD(IB) + CAI(IB) FSA = FSA + CAD(IB) + CAI(IB) TRD = SOLAD(IB)*FTDD(IB) TRI = SOLAD(IB)*FTID(IB) + SOLAI(IB)*FTII(IB) ABS = TRD*(1.-ALBGRD(IB)) + TRI*(1.-ALBGRI(IB)) SAG = SAG + ABS FSA = FSA + ABS END DO LAIFRA = ELAI / MAX(VAI,MPE) IF (FSUN .GT. 0.) THEN PARSUN = (CAD(1)+FSUN*CAI(1)) * LAIFRA / MAX(LAISUN,MPE) PARSHA = (FSHA*CAI(1))*LAIFRA / MAX(LAISHA,MPE) ELSE PARSUN = 0. PARSHA = (CAD(1)+CAI(1))*LAIFRA /MAX(LAISHA,MPE) ENDIF RVIS = ALBD(1)*SOLAD(1) + ALBI(1)*SOLAI(1) RNIR = ALBD(2)*SOLAD(2) + ALBI(2)*SOLAI(2) FSR = RVIS + RNIR FSRV = FREVD(1)*SOLAD(1)+FREVI(1)*SOLAI(1)+FREVD(2)*SOLAD(2)+FREVI(2)*SOLAI(2) FSRG = FREGD(1)*SOLAD(1)+FREGI(1)*SOLAI(1)+FREGD(2)*SOLAD(2)+FREGI(2)*SOLAI(2) END SUBROUTINE SURRAD SUBROUTINE SNOW_AGE (DT,TG,SNEQVO,SNEQV,FAGE) IMPLICIT NONE REAL, INTENT(IN) :: DT REAL, INTENT(IN) :: TG REAL, INTENT(IN) :: SNEQVO REAL, INTENT(IN) :: SNEQV REAL, INTENT(OUT) :: FAGE REAL :: TAUSS REAL :: TAGE REAL :: AGE1 REAL :: AGE2 REAL :: AGE3 REAL :: DELA REAL :: SGE REAL :: DELS REAL :: DELA0 REAL :: ARG IF(SNEQV.LE.0.0) THEN TAUSS = 0. ELSE IF (SNEQV.GT.800.) THEN TAUSS = 0. ELSE TAUSS = 0. DELA0 = 1.E-6*DT ARG = 5.E3*(1./TFRZ-1./TG) AGE1 = EXP(ARG) AGE2 = EXP(AMIN1(0.,10.*ARG)) AGE3 = 0.3 TAGE = AGE1+AGE2+AGE3 DELA = DELA0*TAGE DELS = AMAX1(0.0,SNEQV-SNEQVO) / SWEMX SGE = (TAUSS+DELA)*(1.0-DELS) TAUSS = AMAX1(0.,SGE) ENDIF FAGE= TAUSS/(TAUSS+1.) END SUBROUTINE SNOW_AGE SUBROUTINE SNOWALB_BATS (NBAND,FSNO,COSZ,FAGE,ALBSND,ALBSNI) IMPLICIT NONE INTEGER,INTENT(IN) :: NBAND REAL,INTENT(IN) :: COSZ REAL,INTENT(IN) :: FSNO REAL,INTENT(IN) :: FAGE REAL, DIMENSION(1:2),INTENT(OUT) :: ALBSND REAL, DIMENSION(1:2),INTENT(OUT) :: ALBSNI INTEGER :: IB REAL :: FZEN REAL :: CF1 REAL :: SL2 REAL :: SL1 REAL :: SL REAL, PARAMETER :: C1 = 0.2 REAL, PARAMETER :: C2 = 0.5 ALBSND(1: NBAND) = 0. ALBSNI(1: NBAND) = 0. SL=2.0 SL1=1./SL SL2=2.*SL CF1=((1.+SL1)/(1.+SL2*COSZ)-SL1) FZEN=AMAX1(CF1,0.) ALBSNI(1)=0.95*(1.-C1*FAGE) ALBSNI(2)=0.65*(1.-C2*FAGE) ALBSND(1)=ALBSNI(1)+0.4*FZEN*(1.-ALBSNI(1)) ALBSND(2)=ALBSNI(2)+0.4*FZEN*(1.-ALBSNI(2)) END SUBROUTINE SNOWALB_BATS SUBROUTINE SNOWALB_CLASS (NBAND,QSNOW,DT,ALB,ALBOLD,ALBSND,ALBSNI,ILOC,JLOC) IMPLICIT NONE INTEGER,INTENT(IN) :: ILOC INTEGER,INTENT(IN) :: JLOC INTEGER,INTENT(IN) :: NBAND REAL,INTENT(IN) :: QSNOW REAL,INTENT(IN) :: DT REAL,INTENT(IN) :: ALBOLD REAL, INTENT(INOUT) :: ALB REAL, DIMENSION(1:2),INTENT(OUT) :: ALBSND REAL, DIMENSION(1:2),INTENT(OUT) :: ALBSNI INTEGER :: IB ALBSND(1: NBAND) = 0. ALBSNI(1: NBAND) = 0. ALB = 0.55 + (ALBOLD-0.55) * EXP(-0.01*DT/3600.) IF (QSNOW > 0.) then ALB = ALB + MIN(QSNOW*DT,SWEMX) * (0.84-ALB)/(SWEMX) ENDIF ALBSNI(1)= ALB ALBSNI(2)= ALB ALBSND(1)= ALB ALBSND(2)= ALB END SUBROUTINE SNOWALB_CLASS SUBROUTINE GROUNDALB (NSOIL ,NBAND ,ICE ,IST ,ISC , & FSNO ,SMC ,ALBSND ,ALBSNI ,COSZ , & TG ,ILOC ,JLOC , & ALBGRD ,ALBGRI ) USE NOAHMP_RAD_PARAMETERS IMPLICIT NONE INTEGER, INTENT(IN) :: ILOC INTEGER, INTENT(IN) :: JLOC INTEGER, INTENT(IN) :: NSOIL INTEGER, INTENT(IN) :: NBAND INTEGER, INTENT(IN) :: ICE INTEGER, INTENT(IN) :: IST INTEGER, INTENT(IN) :: ISC REAL, INTENT(IN) :: FSNO REAL, INTENT(IN) :: TG REAL, INTENT(IN) :: COSZ REAL, DIMENSION(1:NSOIL), INTENT(IN) :: SMC REAL, DIMENSION(1: 2), INTENT(IN) :: ALBSND REAL, DIMENSION(1: 2), INTENT(IN) :: ALBSNI REAL, DIMENSION(1: 2), INTENT(OUT) :: ALBGRD REAL, DIMENSION(1: 2), INTENT(OUT) :: ALBGRI INTEGER :: IB REAL :: INC REAL :: ALBSOD REAL :: ALBSOI DO IB = 1, NBAND INC = MAX(0.11-0.40*SMC(1), 0.) IF (IST .EQ. 1) THEN ALBSOD = MIN(ALBSAT(ISC,IB)+INC,ALBDRY(ISC,IB)) ALBSOI = ALBSOD ELSE IF (TG .GT. TFRZ) THEN ALBSOD = 0.06/(MAX(0.01,COSZ)**1.7 + 0.15) ALBSOI = 0.06 ELSE ALBSOD = ALBLAK(IB) ALBSOI = ALBSOD END IF IF (IST .EQ. 1 .AND. ISC .EQ. 9) THEN ALBSOD = ALBSOD + 0.10 ALBSOI = ALBSOI + 0.10 end if ALBGRD(IB) = ALBSOD*(1.-FSNO) + ALBSND(IB)*FSNO ALBGRI(IB) = ALBSOI*(1.-FSNO) + ALBSNI(IB)*FSNO END DO END SUBROUTINE GROUNDALB SUBROUTINE TWOSTREAM (IB ,IC ,VEGTYP ,COSZ ,VAI , & FWET ,T ,ALBGRD ,ALBGRI ,RHO , & TAU ,FVEG ,IST ,ILOC ,JLOC , & FAB ,FRE ,FTD ,FTI ,GDIR , & FREV ,FREG ,BGAP ,WGAP ,GAP) USE NOAHMP_VEG_PARAMETERS USE NOAHMP_RAD_PARAMETERS IMPLICIT NONE INTEGER, INTENT(IN) :: ILOC INTEGER, INTENT(IN) :: JLOC INTEGER, INTENT(IN) :: IST INTEGER, INTENT(IN) :: IB INTEGER, INTENT(IN) :: IC INTEGER, INTENT(IN) :: VEGTYP REAL, INTENT(IN) :: COSZ REAL, INTENT(IN) :: VAI REAL, INTENT(IN) :: FWET REAL, INTENT(IN) :: T REAL, DIMENSION(1:2), INTENT(IN) :: ALBGRD REAL, DIMENSION(1:2), INTENT(IN) :: ALBGRI REAL, DIMENSION(1:2), INTENT(IN) :: RHO REAL, DIMENSION(1:2), INTENT(IN) :: TAU REAL, INTENT(IN) :: FVEG REAL, DIMENSION(1:2), INTENT(OUT) :: FAB REAL, DIMENSION(1:2), INTENT(OUT) :: FRE REAL, DIMENSION(1:2), INTENT(OUT) :: FTD REAL, DIMENSION(1:2), INTENT(OUT) :: FTI REAL, INTENT(OUT) :: GDIR REAL, DIMENSION(1:2), INTENT(OUT) :: FREV REAL, DIMENSION(1:2), INTENT(OUT) :: FREG REAL :: OMEGA REAL :: OMEGAL REAL :: BETAI REAL :: BETAIL REAL :: BETAD REAL :: BETADL REAL :: EXT REAL :: AVMU REAL :: COSZI REAL :: ASU REAL :: CHIL REAL :: TMP0,TMP1,TMP2,TMP3,TMP4,TMP5,TMP6,TMP7,TMP8,TMP9 REAL :: P1,P2,P3,P4,S1,S2,U1,U2,U3 REAL :: B,C,D,D1,D2,F,H,H1,H2,H3,H4,H5,H6,H7,H8,H9,H10 REAL :: PHI1,PHI2,SIGMA REAL :: FTDS,FTIS,FRES REAL :: FREVEG,FREBAR,FTDVEG,FTIVEG,FTDBAR,FTIBAR REAL :: THETAZ REAL, PARAMETER :: PAI = 3.14159265 REAL :: HD REAL :: BB REAL :: THETAP REAL :: FA REAL :: NEWVAI REAL,INTENT(INOUT) :: BGAP REAL,INTENT(INOUT) :: WGAP REAL :: KOPEN REAL, INTENT(OUT) :: GAP if(VAI == 0.0) THEN GAP = 1.0 KOPEN = 1.0 ELSE IF(OPT_RAD == 1) THEN HD = HVT(VEGTYP) - HVB(VEGTYP) BB = 0.5 * HD THETAP = ATAN(BB/RC(VEGTYP) * TAN(ACOS(MAX(0.01,COSZ))) ) BGAP = EXP(-DEN(VEGTYP) * PAI * RC(VEGTYP)**2/COS(THETAP) ) FA = VAI/(1.33 * PAI * RC(VEGTYP)**3.0 *(BB/RC(VEGTYP))*DEN(VEGTYP)) NEWVAI = HD*FA WGAP = (1.0-BGAP) * EXP(-0.5*NEWVAI/COSZ) BGAP = (1.0-FVEG)*BGAP WGAP = FVEG*WGAP GAP = MIN(1.0-FVEG, BGAP+WGAP) KOPEN = 0.05 END IF IF(OPT_RAD == 2) THEN GAP = 0.0 KOPEN = 0.0 END IF IF(OPT_RAD == 3) THEN GAP = 1.0-FVEG KOPEN = 0.0 END IF end if COSZI = MAX(0.001, COSZ) CHIL = MIN( MAX(XL(VEGTYP), -0.4), 0.6) IF (ABS(CHIL) .LE. 0.01) CHIL = 0.01 PHI1 = 0.5 - 0.633*CHIL - 0.330*CHIL*CHIL PHI2 = 0.877 * (1.-2.*PHI1) GDIR = PHI1 + PHI2*COSZI EXT = GDIR/COSZI AVMU = ( 1. - PHI1/PHI2 * LOG((PHI1+PHI2)/PHI1) ) / PHI2 OMEGAL = RHO(IB) + TAU(IB) TMP0 = GDIR + PHI2*COSZI TMP1 = PHI1*COSZI ASU = 0.5*OMEGAL*GDIR/TMP0 * ( 1.-TMP1/TMP0*LOG((TMP1+TMP0)/TMP1) ) BETADL = (1.+AVMU*EXT)/(OMEGAL*AVMU*EXT)*ASU BETAIL = 0.5 * ( RHO(IB)+TAU(IB) + (RHO(IB)-TAU(IB)) & * ((1.+CHIL)/2.)**2 ) / OMEGAL IF (T .GT. TFRZ) THEN TMP0 = OMEGAL TMP1 = BETADL TMP2 = BETAIL ELSE TMP0 = (1.-FWET)*OMEGAL + FWET*OMEGAS(IB) TMP1 = ( (1.-FWET)*OMEGAL*BETADL + FWET*OMEGAS(IB)*BETADS ) / TMP0 TMP2 = ( (1.-FWET)*OMEGAL*BETAIL + FWET*OMEGAS(IB)*BETAIS ) / TMP0 END IF OMEGA = TMP0 BETAD = TMP1 BETAI = TMP2 B = 1. - OMEGA + OMEGA*BETAI C = OMEGA*BETAI TMP0 = AVMU*EXT D = TMP0 * OMEGA*BETAD F = TMP0 * OMEGA*(1.-BETAD) TMP1 = B*B - C*C H = SQRT(TMP1) / AVMU SIGMA = TMP0*TMP0 - TMP1 if(SIGMA == 0.) SIGMA = 1.e-6 P1 = B + AVMU*H P2 = B - AVMU*H P3 = B + TMP0 P4 = B - TMP0 S1 = EXP(-H*VAI) S2 = EXP(-EXT*VAI) IF (IC .EQ. 0) THEN U1 = B - C/ALBGRD(IB) U2 = B - C*ALBGRD(IB) U3 = F + C*ALBGRD(IB) ELSE U1 = B - C/ALBGRI(IB) U2 = B - C*ALBGRI(IB) U3 = F + C*ALBGRI(IB) END IF TMP2 = U1 - AVMU*H TMP3 = U1 + AVMU*H D1 = P1*TMP2/S1 - P2*TMP3*S1 TMP4 = U2 + AVMU*H TMP5 = U2 - AVMU*H D2 = TMP4/S1 - TMP5*S1 H1 = -D*P4 - C*F TMP6 = D - H1*P3/SIGMA TMP7 = ( D - C - H1/SIGMA*(U1+TMP0) ) * S2 H2 = ( TMP6*TMP2/S1 - P2*TMP7 ) / D1 H3 = - ( TMP6*TMP3*S1 - P1*TMP7 ) / D1 H4 = -F*P3 - C*D TMP8 = H4/SIGMA TMP9 = ( U3 - TMP8*(U2-TMP0) ) * S2 H5 = - ( TMP8*TMP4/S1 + TMP9 ) / D2 H6 = ( TMP8*TMP5*S1 + TMP9 ) / D2 H7 = (C*TMP2) / (D1*S1) H8 = (-C*TMP3*S1) / D1 H9 = TMP4 / (D2*S1) H10 = (-TMP5*S1) / D2 IF (IC .EQ. 0) THEN FTDS = S2 *(1.0-GAP) + GAP FTIS = (H4*S2/SIGMA + H5*S1 + H6/S1)*(1.0-GAP) ELSE FTDS = 0. FTIS = (H9*S1 + H10/S1)*(1.0-KOPEN) + KOPEN END IF FTD(IB) = FTDS FTI(IB) = FTIS IF (IC .EQ. 0) THEN FRES = (H1/SIGMA + H2 + H3)*(1.0-GAP ) + ALBGRD(IB)*GAP FREVEG = (H1/SIGMA + H2 + H3)*(1.0-GAP ) FREBAR = ALBGRD(IB)*GAP ELSE FRES = (H7 + H8) *(1.0-KOPEN) + ALBGRI(IB)*KOPEN FREVEG = (H7 + H8) *(1.0-KOPEN)+ALBGRI(IB)*KOPEN FREBAR = 0 END IF FRE(IB) = FRES FREV(IB) = FREVEG FREG(IB) = FREBAR FAB(IB) = 1. - FRE(IB) - (1.-ALBGRD(IB))*FTD(IB) & - (1.-ALBGRI(IB))*FTI(IB) END SUBROUTINE TWOSTREAM SUBROUTINE VEGE_FLUX(NSNOW ,NSOIL ,ISNOW ,VEGTYP ,VEG , & DT ,SAV ,SAG ,LWDN ,UR , & UU ,VV ,SFCTMP ,THAIR ,QAIR , & EAIR ,RHOAIR ,SNOWH ,VAI ,GAMMA , & FWET ,LAISUN ,LAISHA ,CWP ,DZSNSO , & HTOP ,ZLVL ,ZPD ,Z0M ,FVEG , & Z0MG ,EMV ,EMG ,CANLIQ , & CANICE ,STC ,DF ,RSSUN ,RSSHA , & RSURF ,LATHEA ,PARSUN ,PARSHA ,IGS , & FOLN ,CO2AIR ,O2AIR ,BTRAN ,SFCPRS , & RHSUR ,ILOC ,JLOC ,Q2 , & EAH ,TAH ,TV ,TG ,CM , & CH ,DX ,DZ8W , & TAUXV ,TAUYV ,IRG ,IRC ,SHG , & SHC ,EVG ,EVC ,TR ,GH , & T2MV ,PSNSUN ,PSNSHA , & QC ,PBLH ,QSFC ,PSFC ,ISURBAN , & IZ0TLND ,Q2V ,CAH2) USE NOAHMP_VEG_PARAMETERS USE MODULE_MODEL_CONSTANTS IMPLICIT NONE INTEGER, INTENT(IN) :: ILOC INTEGER, INTENT(IN) :: JLOC LOGICAL, INTENT(IN) :: VEG INTEGER, INTENT(IN) :: NSNOW INTEGER, INTENT(IN) :: NSOIL INTEGER, INTENT(IN) :: ISNOW INTEGER, INTENT(IN) :: VEGTYP REAL, INTENT(IN) :: FVEG REAL, INTENT(IN) :: SAV REAL, INTENT(IN) :: SAG REAL, INTENT(IN) :: LWDN REAL, INTENT(IN) :: UR REAL, INTENT(IN) :: UU REAL, INTENT(IN) :: VV REAL, INTENT(IN) :: SFCTMP REAL, INTENT(IN) :: THAIR REAL, INTENT(IN) :: EAIR REAL, INTENT(IN) :: QAIR REAL, INTENT(IN) :: RHOAIR REAL, INTENT(IN) :: DT REAL, INTENT(IN) :: SNOWH REAL, INTENT(IN) :: FWET REAL, INTENT(IN) :: HTOP REAL, INTENT(IN) :: CWP REAL, INTENT(IN) :: VAI REAL, INTENT(IN) :: LAISUN REAL, INTENT(IN) :: LAISHA REAL, INTENT(IN) :: ZLVL REAL, INTENT(IN) :: ZPD REAL, INTENT(IN) :: Z0M REAL, INTENT(IN) :: Z0MG REAL, INTENT(IN) :: EMV REAL, INTENT(IN) :: EMG REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(IN) :: STC REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(IN) :: DF REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(IN) :: DZSNSO REAL, INTENT(IN) :: CANLIQ REAL, INTENT(IN) :: CANICE REAL, INTENT(IN) :: RSURF REAL, INTENT(IN) :: GAMMA REAL, INTENT(IN) :: LATHEA REAL, INTENT(IN) :: PARSUN REAL, INTENT(IN) :: PARSHA REAL, INTENT(IN) :: FOLN REAL, INTENT(IN) :: CO2AIR REAL, INTENT(IN) :: O2AIR REAL, INTENT(IN) :: IGS REAL, INTENT(IN) :: SFCPRS REAL, INTENT(IN) :: BTRAN REAL, INTENT(IN) :: RHSUR INTEGER , INTENT(IN) :: ISURBAN INTEGER , INTENT(IN) :: IZ0TLND REAL , INTENT(IN) :: QC REAL , INTENT(IN) :: PBLH REAL , INTENT(IN) :: PSFC REAL , INTENT(IN) :: DX REAL , INTENT(IN) :: Q2 REAL , INTENT(IN) :: DZ8W REAL , INTENT(INOUT) :: QSFC REAL, INTENT(INOUT) :: EAH REAL, INTENT(INOUT) :: TAH REAL, INTENT(INOUT) :: TV REAL, INTENT(INOUT) :: TG REAL, INTENT(INOUT) :: CM REAL, INTENT(INOUT) :: CH REAL, INTENT(OUT) :: TAUXV REAL, INTENT(OUT) :: TAUYV REAL, INTENT(OUT) :: IRC REAL, INTENT(OUT) :: SHC REAL, INTENT(OUT) :: EVC REAL, INTENT(OUT) :: IRG REAL, INTENT(OUT) :: SHG REAL, INTENT(OUT) :: EVG REAL, INTENT(OUT) :: TR REAL, INTENT(OUT) :: GH REAL, INTENT(OUT) :: T2MV REAL, INTENT(OUT) :: PSNSUN REAL, INTENT(OUT) :: PSNSHA REAL, INTENT(OUT) :: Q2V REAL :: CAH REAL :: U10V REAL :: V10V REAL :: WSPD REAL :: CW REAL :: FV REAL :: WSTAR REAL :: Z0H REAL :: Z0HG REAL :: RB REAL :: RAMC REAL :: RAHC REAL :: RAWC REAL :: RAMG REAL :: RAHG REAL :: RAWG REAL, INTENT(OUT) :: RSSUN REAL, INTENT(OUT) :: RSSHA REAL :: MOL REAL :: DTV REAL :: DTG REAL :: AIR,CIR REAL :: CSH REAL :: CEV REAL :: CGH REAL :: ATR,CTR REAL :: ATA,BTA REAL :: AEA,BEA REAL :: ESTV REAL :: ESTG REAL :: DESTV REAL :: DESTG REAL :: ESATW REAL :: ESATI REAL :: DSATW REAL :: DSATI REAL :: FM REAL :: FH REAL :: FHG REAL :: HCAN REAL :: A REAL :: B REAL :: CVH REAL :: CAW REAL :: CTW REAL :: CEW REAL :: CGW REAL :: COND REAL :: UC REAL :: KH REAL :: H REAL :: HG REAL :: MOZ REAL :: MOZG REAL :: MOZOLD REAL :: THVAIR REAL :: THAH REAL :: RAHC2 REAL :: RAWC2 REAL, INTENT(OUT):: CAH2 REAL :: CH2V REAL :: CQ2V REAL :: EAH2 REAL :: QFX REAL :: E1 REAL :: VAIE REAL :: LAISUNE REAL :: LAISHAE INTEGER :: K INTEGER :: ITER INTEGER, PARAMETER :: NITERC = 20 INTEGER, PARAMETER :: NITERG = 5 INTEGER :: MOZSGN REAL :: MPE INTEGER :: LITER REAL :: T, TDC character(len=80) :: message TDC(T) = MIN( 50., MAX(-50.,(T-TFRZ)) ) MPE = 1E-6 LITER = 0 FV = 0.1 DTV = 0. DTG = 0. MOZSGN = 0 MOZOLD = 0. HG = 0. H = 0. QFX = 0. VAIE = MIN(6.,VAI / FVEG) LAISUNE = MIN(6.,LAISUN / FVEG) LAISHAE = MIN(6.,LAISHA / FVEG) T = TDC(TG) CALL ESAT(T, ESATW, ESATI, DSATW, DSATI) IF (T .GT. 0.) THEN ESTG = ESATW ELSE ESTG = ESATI END IF QSFC = 0.622*EAIR/(PSFC-0.378*EAIR) HCAN = HTOP UC = UR*LOG(HCAN/Z0M)/LOG(ZLVL/Z0M) IF((HCAN-ZPD) <= 0.) THEN WRITE(message,*) "CRITICAL PROBLEM: HCAN <= ZPD" call wrf_message ( message ) WRITE(message,*) 'i,j point=',ILOC, JLOC call wrf_message ( message ) WRITE(message,*) 'HCAN =',HCAN call wrf_message ( message ) WRITE(message,*) 'ZPD =',ZPD call wrf_message ( message ) write (message, *) 'SNOWH =',SNOWH call wrf_message ( message ) call wrf_error_fatal3("",3227,& "CRITICAL PROBLEM IN MODULE_SF_NOAHMPLSM:VEGEFLUX" ) END IF AIR = -EMV*(1.+(1.-EMV)*(1.-EMG))*LWDN - EMV*EMG*SB*TG**4 CIR = (2.-EMV*(1.-EMG))*EMV*SB loop1: DO ITER = 1, NITERC IF(ITER == 1) THEN Z0H = Z0M Z0HG = Z0MG ELSE Z0H = Z0M Z0HG = Z0MG END IF IF(OPT_SFC == 1) THEN CALL SFCDIF1(ITER ,SFCTMP ,RHOAIR ,H ,QAIR , & ZLVL ,ZPD ,Z0M ,Z0H ,UR , & MPE ,ILOC ,JLOC , & MOZ ,MOZSGN ,FM ,FH , & CM ,CH ,FV ) ENDIF IF(OPT_SFC == 2) THEN CALL SFCDIF2(ITER ,Z0M ,TAH ,THAIR ,UR , & CZIL ,ZLVL ,ILOC ,JLOC , & CM ,CH ,MOZ ,WSTAR , & FV ) CH = CH / UR CM = CM / UR ENDIF IF(OPT_SFC == 3) THEN CALL SFCDIF3(ILOC ,JLOC ,TAH ,QSFC ,PSFC ,& PBLH ,Z0M ,Z0MG ,VEGTYP ,ISURBAN,& IZ0TLND,UC ,ITER ,NITERC ,SFCTMP ,& THAIR ,QAIR ,QC ,ZLVL , & SFCPRS ,FV ,CM ,CH ,CH2V ,& CQ2V ,MOZ) CH = CH / UR CM = CM / UR CH2V = CH2V / UR ENDIF IF(OPT_SFC == 4) THEN CALL SFCDIF4(ILOC ,JLOC ,UU ,VV ,SFCTMP ,& SFCPRS ,PSFC ,PBLH ,DX ,Z0M ,& TAH ,QAIR ,ZLVL ,IZ0TLND,QSFC ,& H ,QFX ,CM ,CH ,CH2V ,& CQ2V ,MOZ ,FV ,U10V ,V10V) CH = CH / UR CM = CM / UR CH2V = CH2V / UR ENDIF RAMC = MAX(1.,1./(CM*UR)) RAHC = MAX(1.,1./(CH*UR)) RAWC = RAHC IF (OPT_SFC == 3 .OR. OPT_SFC == 4 ) THEN RAHC2 = MAX(1.,1./(CH2V*UR)) RAWC2 = RAHC2 CAH2 = 1./RAHC2 CQ2V = CAH2 ENDIF CALL RAGRB(ITER ,VAIE ,RHOAIR ,HG ,TAH , & ZPD ,Z0MG ,Z0HG ,HCAN ,UC , & Z0H ,FV ,CWP ,VEGTYP ,MPE , & TV ,MOZG ,FHG ,ILOC ,JLOC , & RAMG ,RAHG ,RAWG ,RB ) T = TDC(TV) CALL ESAT(T, ESATW, ESATI, DSATW, DSATI) IF (T .GT. 0.) THEN ESTV = ESATW DESTV = DSATW ELSE ESTV = ESATI DESTV = DSATI END IF IF(ITER == 1) THEN IF (OPT_CRS == 1) then CALL STOMATA (VEGTYP,MPE ,PARSUN ,FOLN ,ILOC , JLOC , & TV ,ESTV ,EAH ,SFCTMP,SFCPRS, & O2AIR ,CO2AIR,IGS ,BTRAN ,RB , & RSSUN ,PSNSUN) CALL STOMATA (VEGTYP,MPE ,PARSHA ,FOLN ,ILOC , JLOC , & TV ,ESTV ,EAH ,SFCTMP,SFCPRS, & O2AIR ,CO2AIR,IGS ,BTRAN ,RB , & RSSHA ,PSNSHA) END IF IF (OPT_CRS == 2) then CALL CANRES (PARSUN,TV ,BTRAN ,EAH ,SFCPRS, & RSSUN ,PSNSUN,ILOC ,JLOC ) CALL CANRES (PARSHA,TV ,BTRAN ,EAH ,SFCPRS, & RSSHA ,PSNSHA,ILOC ,JLOC ) END IF END IF CAH = 1./RAHC CVH = 2.*VAIE/RB CGH = 1./RAHG COND = CAH + CVH + CGH ATA = (SFCTMP*CAH + TG*CGH) / COND BTA = CVH/COND CSH = (1.-BTA)*RHOAIR*CPAIR*CVH CAW = 1./RAWC CEW = FWET*VAIE/RB CTW = (1.-FWET)*(LAISUNE/(RB+RSSUN) + LAISHAE/(RB+RSSHA)) CGW = 1./(RAWG+RSURF) COND = CAW + CEW + CTW + CGW AEA = (EAIR*CAW + ESTG*CGW) / COND BEA = (CEW+CTW)/COND CEV = (1.-BEA)*CEW*RHOAIR*CPAIR/GAMMA CTR = (1.-BEA)*CTW*RHOAIR*CPAIR/GAMMA TAH = ATA + BTA*TV EAH = AEA + BEA*ESTV IRC = AIR + CIR*TV**4 SHC = RHOAIR*CPAIR*CVH * ( TV-TAH) EVC = RHOAIR*CPAIR*CEW * (ESTV-EAH) / GAMMA TR = RHOAIR*CPAIR*CTW * (ESTV-EAH) / GAMMA EVC = MIN(CANLIQ*LATHEA/DT,EVC) B = SAV-IRC-SHC-EVC-TR A = 4.*CIR*TV**3 + CSH + (CEV+CTR)*DESTV DTV = B/A IRC = IRC + 4.*CIR*TV**3*DTV SHC = SHC + CSH*DTV EVC = EVC + CEV*DESTV*DTV TR = TR + CTR*DESTV*DTV TV = TV + DTV H = RHOAIR*CPAIR*(TAH - SFCTMP) /RAHC HG = RHOAIR*CPAIR*(TG - TAH) /RAHG QSFC = (0.622*EAH)/(SFCPRS-0.378*EAH) IF ( OPT_SFC == 4 ) THEN QFX = (QSFC-QAIR)*RHOAIR*CAW ENDIF IF (LITER == 1) THEN exit loop1 ENDIF IF (ITER >= 5 .AND. ABS(DTV) <= 0.01 .AND. LITER == 0) THEN LITER = 1 ENDIF END DO loop1 AIR = - EMG*(1.-EMV)*LWDN - EMG*EMV*SB*TV**4 CIR = EMG*SB CSH = RHOAIR*CPAIR/RAHG CEV = RHOAIR*CPAIR / (GAMMA*(RAWG+RSURF)) CGH = 2.*DF(ISNOW+1)/DZSNSO(ISNOW+1) loop2: DO ITER = 1, NITERG T = TDC(TG) CALL ESAT(T, ESATW, ESATI, DSATW, DSATI) IF (T .GT. 0.) THEN ESTG = ESATW DESTG = DSATW ELSE ESTG = ESATI DESTG = DSATI END IF IRG = CIR*TG**4 + AIR SHG = CSH * (TG - TAH ) EVG = CEV * (ESTG*RHSUR - EAH ) GH = CGH * (TG - STC(ISNOW+1)) B = SAG-IRG-SHG-EVG-GH A = 4.*CIR*TG**3+CSH+CEV*DESTG+CGH DTG = B/A IRG = IRG + 4.*CIR*TG**3*DTG SHG = SHG + CSH*DTG EVG = EVG + CEV*DESTG*DTG GH = GH + CGH*DTG TG = TG + DTG END DO loop2 IF(OPT_STC == 1) THEN IF (SNOWH > 0.05 .AND. TG > TFRZ) THEN TG = TFRZ IRG = CIR*TG**4 - EMG*(1.-EMV)*LWDN - EMG*EMV*SB*TV**4 SHG = CSH * (TG - TAH) EVG = CEV * (ESTG*RHSUR - EAH) GH = SAG - (IRG+SHG+EVG) END IF END IF TAUXV = -RHOAIR*CM*UR*UU TAUYV = -RHOAIR*CM*UR*VV TAH = SFCTMP + (SHG+SHC)/(RHOAIR*CPAIR*CAH) EAH = EAIR + (EVC+TR+EVG)/(RHOAIR*CAW*CPAIR/GAMMA ) QFX = (QSFC-QAIR)*RHOAIR*CAW IF (OPT_SFC == 1 .OR. OPT_SFC == 2) THEN CAH2 = FV*1./VKC*LOG((2.+Z0H)/Z0H) CQ2V = CAH2 IF (CAH2 .LT. 1.E-5 ) THEN T2MV = TAH Q2V = (EAH*0.622/(SFCPRS - 0.378*EAH)) ELSE T2MV = TAH - (SHG+SHC)/(RHOAIR*CPAIR*FV) * 1./VKC * LOG((2.+Z0H)/Z0H) Q2V = (EAH*0.622/(SFCPRS - 0.378*EAH))- QFX/(RHOAIR*FV)* 1./VKC * LOG((2.+Z0H)/Z0H) ENDIF ENDIF IF (OPT_SFC == 3 .OR. OPT_SFC == 4 ) THEN IF (CAH2 .LT. 1.E-5 ) THEN T2MV = TAH Q2V = (EAH*0.622/(SFCPRS - 0.378*EAH)) ELSE T2MV = TAH - (SHG+SHC)/(RHOAIR*CPAIR*CAH2) Q2V = (EAH*0.622/(SFCPRS - 0.378*EAH)) - QFX/(RHOAIR*CQ2V) ENDIF ENDIF CH = CAH END SUBROUTINE VEGE_FLUX SUBROUTINE BARE_FLUX (NSNOW ,NSOIL ,ISNOW ,DT ,SAG , & LWDN ,UR ,UU ,VV ,SFCTMP , & THAIR ,QAIR ,EAIR ,RHOAIR ,SNOWH , & DZSNSO ,ZLVL ,ZPD ,Z0M , & EMG ,STC ,DF ,RSURF ,LATHEA , & GAMMA ,RHSUR ,ILOC ,JLOC ,Q2 , & TGB ,CM ,CH , & TAUXB ,TAUYB ,IRB ,SHB ,EVB , & GHB ,T2MB ,DX ,DZ8W ,IVGTYP , & QC ,PBLH ,QSFC ,PSFC ,ISURBAN , & IZ0TLND ,SFCPRS ,Q2B ,EHB2) USE NOAHMP_VEG_PARAMETERS USE MODULE_MODEL_CONSTANTS IMPLICIT NONE integer , INTENT(IN) :: ILOC integer , INTENT(IN) :: JLOC INTEGER, INTENT(IN) :: NSNOW INTEGER, INTENT(IN) :: NSOIL INTEGER, INTENT(IN) :: ISNOW REAL, INTENT(IN) :: DT REAL, INTENT(IN) :: SAG REAL, INTENT(IN) :: LWDN REAL, INTENT(IN) :: UR REAL, INTENT(IN) :: UU REAL, INTENT(IN) :: VV REAL, INTENT(IN) :: SFCTMP REAL, INTENT(IN) :: THAIR REAL, INTENT(IN) :: QAIR REAL, INTENT(IN) :: EAIR REAL, INTENT(IN) :: RHOAIR REAL, INTENT(IN) :: SNOWH REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(IN) :: DZSNSO REAL, INTENT(IN) :: ZLVL REAL, INTENT(IN) :: ZPD REAL, INTENT(IN) :: Z0M REAL, INTENT(IN) :: EMG REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(IN) :: STC REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(IN) :: DF REAL, INTENT(IN) :: RSURF REAL, INTENT(IN) :: LATHEA REAL, INTENT(IN) :: GAMMA REAL, INTENT(IN) :: RHSUR INTEGER , INTENT(IN) :: ISURBAN INTEGER , INTENT(IN) :: IVGTYP INTEGER , INTENT(IN) :: IZ0TLND REAL , INTENT(IN) :: QC REAL , INTENT(IN) :: PBLH REAL , INTENT(INOUT) :: QSFC REAL , INTENT(IN) :: PSFC REAL , INTENT(IN) :: SFCPRS REAL , INTENT(IN) :: DX REAL , INTENT(IN) :: Q2 REAL , INTENT(IN) :: DZ8W REAL, INTENT(INOUT) :: TGB REAL, INTENT(INOUT) :: CM REAL, INTENT(INOUT) :: CH REAL, INTENT(OUT) :: TAUXB REAL, INTENT(OUT) :: TAUYB REAL, INTENT(OUT) :: IRB REAL, INTENT(OUT) :: SHB REAL, INTENT(OUT) :: EVB REAL, INTENT(OUT) :: GHB REAL, INTENT(OUT) :: T2MB REAL, INTENT(OUT) :: Q2B REAL :: EHB REAL :: U10B REAL :: V10B REAL :: WSPD REAL :: TAUX REAL :: TAUY REAL :: FIRA REAL :: FSH REAL :: FGEV REAL :: SSOIL REAL :: FIRE REAL :: TRAD REAL :: TAH REAL :: CW REAL :: FV REAL :: WSTAR REAL :: Z0H REAL :: RB REAL :: RAMB REAL :: RAHB REAL :: RAWB REAL :: MOL REAL :: DTG REAL :: CIR REAL :: CSH REAL :: CEV REAL :: CGH REAL :: RAHB2 REAL :: RAWB2 REAL,INTENT(OUT) :: EHB2 REAL :: CH2B REAL :: CQ2B REAL :: THVAIR REAL :: THGH REAL :: EMB REAL :: QFX REAL :: ESTG2 INTEGER :: VEGTYP REAL :: E1 REAL :: ESTG REAL :: DESTG REAL :: ESATW REAL :: ESATI REAL :: DSATW REAL :: DSATI REAL :: A REAL :: B REAL :: H REAL :: MOZ REAL :: MOZOLD REAL :: FM REAL :: FH INTEGER :: MOZSGN INTEGER :: ITER INTEGER :: NITERB REAL :: MPE DATA NITERB /5/ SAVE NITERB REAL :: T, TDC TDC(T) = MIN( 50., MAX(-50.,(T-TFRZ)) ) MPE = 1E-6 DTG = 0. MOZSGN = 0 MOZOLD = 0. H = 0. QFX = 0. FV = 0.1 CIR = EMG*SB CGH = 2.*DF(ISNOW+1)/DZSNSO(ISNOW+1) loop3: DO ITER = 1, NITERB IF(ITER == 1) THEN Z0H = Z0M ELSE Z0H = Z0M END IF IF(OPT_SFC == 1) THEN CALL SFCDIF1(ITER ,SFCTMP ,RHOAIR ,H ,QAIR , & ZLVL ,ZPD ,Z0M ,Z0H ,UR , & MPE ,ILOC ,JLOC , & MOZ ,MOZSGN ,FM ,FH , & CM ,CH ,FV ) ENDIF IF(OPT_SFC == 2) THEN CALL SFCDIF2(ITER ,Z0M ,TGB ,THAIR ,UR , & CZIL ,ZLVL ,ILOC ,JLOC , & CM ,CH ,MOZ ,WSTAR , & FV ) CH = CH / UR CM = CM / UR IF(SNOWH > 0.) THEN CM = MIN(0.01,CM) CH = MIN(0.01,CH) END IF ENDIF IF(OPT_SFC == 3) THEN VEGTYP = ISBARREN CALL SFCDIF3(ILOC ,JLOC ,TGB ,QSFC ,PSFC ,& PBLH ,Z0M ,Z0M ,VEGTYP ,ISURBAN,& IZ0TLND,UR ,ITER ,NITERB ,SFCTMP ,& THAIR ,QAIR ,QC ,ZLVL , & SFCPRS ,FV ,CM ,CH ,CH2B ,& CQ2B ,MOZ) CH = CH / UR CM = CM / UR CH2B = CH2B / UR IF(SNOWH > 0.) THEN CM = MIN(0.01,CM) CH = MIN(0.01,CH) CH2B = MIN(0.01,CH2B) CQ2B = MIN(0.01,CQ2B) END IF ENDIF IF(OPT_SFC == 4) THEN CALL SFCDIF4(ILOC ,JLOC ,UU ,VV ,SFCTMP ,& SFCPRS ,PSFC ,PBLH ,DX ,Z0M ,& TGB ,QAIR ,ZLVL ,IZ0TLND,QSFC ,& H ,QFX ,CM ,CH ,CH2B ,& CQ2B ,MOZ ,FV ,U10B ,V10B) CH = CH / UR CM = CM / UR CH2B = CH2B / UR IF(SNOWH > 0.) THEN CM = MIN(0.01,CM) CH = MIN(0.01,CH) CH2B = MIN(0.01,CH2B) CQ2B = MIN(0.01,CQ2B) END IF ENDIF RAMB = MAX(1.,1./(CM*UR)) RAHB = MAX(1.,1./(CH*UR)) RAWB = RAHB EMB = 1./RAMB EHB = 1./RAHB IF (OPT_SFC == 3 .OR. OPT_SFC == 4) THEN RAHB2 = MAX(1.,1./(CH2B*UR)) EHB2 = 1./RAHB2 CQ2B = EHB2 END IF T = TDC(TGB) CALL ESAT(T, ESATW, ESATI, DSATW, DSATI) IF (T .GT. 0.) THEN ESTG = ESATW DESTG = DSATW ELSE ESTG = ESATI DESTG = DSATI END IF CSH = RHOAIR*CPAIR/RAHB CEV = RHOAIR*CPAIR/GAMMA/(RSURF+RAWB) IRB = CIR * TGB**4 - EMG*LWDN SHB = CSH * (TGB - SFCTMP ) EVB = CEV * (ESTG*RHSUR - EAIR ) GHB = CGH * (TGB - STC(ISNOW+1)) B = SAG-IRB-SHB-EVB-GHB A = 4.*CIR*TGB**3 + CSH + CEV*DESTG + CGH DTG = B/A IRB = IRB + 4.*CIR*TGB**3*DTG SHB = SHB + CSH*DTG EVB = EVB + CEV*DESTG*DTG GHB = GHB + CGH*DTG TGB = TGB + DTG H = CSH * (TGB - SFCTMP) T = TDC(TGB) CALL ESAT(T, ESATW, ESATI, DSATW, DSATI) IF (T .GT. 0.) THEN ESTG = ESATW ELSE ESTG = ESATI END IF QSFC = 0.622*(ESTG*RHSUR)/(PSFC-0.378*(ESTG*RHSUR)) QFX = (QSFC-QAIR)*CEV*GAMMA/CPAIR END DO loop3 IF(OPT_STC == 1) THEN IF (SNOWH > 0.05 .AND. TGB > TFRZ) THEN TGB = TFRZ IRB = CIR * TGB**4 - EMG*LWDN SHB = CSH * (TGB - SFCTMP) EVB = CEV * (ESTG*RHSUR - EAIR ) GHB = SAG - (IRB+SHB+EVB) END IF END IF TAUXB = -RHOAIR*CM*UR*UU TAUYB = -RHOAIR*CM*UR*VV IF(OPT_SFC == 1 .OR. OPT_SFC ==2) THEN EHB2 = FV*VKC/LOG((2.+Z0H)/Z0H) CQ2B = EHB2 IF (EHB2.lt.1.E-5 ) THEN T2MB = TGB Q2B = QSFC ELSE T2MB = TGB - SHB/(RHOAIR*CPAIR*FV) * 1./VKC * LOG((2.+Z0H)/Z0H) Q2B = QSFC - EVB/(LATHEA*RHOAIR)*(1./CQ2B + RSURF) ENDIF IF (IVGTYP == ISURBAN) Q2B = QSFC END IF IF(OPT_SFC ==3 .OR. OPT_SFC == 4) THEN IF (EHB2.lt.1.E-5 ) THEN T2MB = TGB Q2B = QSFC ELSE T2MB = TGB - SHB/(RHOAIR*CPAIR*EHB2) Q2B = QSFC - QFX/(RHOAIR*CQ2B) END IF END IF CH = EHB END SUBROUTINE BARE_FLUX SUBROUTINE RAGRB(ITER ,VAI ,RHOAIR ,HG ,TAH , & ZPD ,Z0MG ,Z0HG ,HCAN ,UC , & Z0H ,FV ,CWP ,VEGTYP ,MPE , & TV ,MOZG ,FHG ,ILOC ,JLOC , & RAMG ,RAHG ,RAWG ,RB ) USE NOAHMP_VEG_PARAMETERS IMPLICIT NONE INTEGER, INTENT(IN) :: ILOC INTEGER, INTENT(IN) :: JLOC INTEGER, INTENT(IN) :: ITER INTEGER, INTENT(IN) :: VEGTYP REAL, INTENT(IN) :: VAI REAL, INTENT(IN) :: RHOAIR REAL, INTENT(IN) :: HG REAL, INTENT(IN) :: TV REAL, INTENT(IN) :: TAH REAL, INTENT(IN) :: ZPD REAL, INTENT(IN) :: Z0MG REAL, INTENT(IN) :: HCAN REAL, INTENT(IN) :: UC REAL, INTENT(IN) :: Z0H REAL, INTENT(IN) :: Z0HG REAL, INTENT(IN) :: FV REAL, INTENT(IN) :: CWP REAL, INTENT(IN) :: MPE REAL, INTENT(INOUT) :: MOZG REAL, INTENT(INOUT) :: FHG REAL :: RAMG REAL :: RAHG REAL :: RAWG REAL :: RB REAL :: KH REAL :: TMP1 REAL :: TMP2 REAL :: TMPRAH2 REAL :: TMPRB real :: MOLG,FHGNEW,CWPC MOZG = 0. MOLG = 0. IF(ITER > 1) THEN TMP1 = VKC * (GRAV/TAH) * HG/(RHOAIR*CPAIR) IF (ABS(TMP1) .LE. MPE) TMP1 = MPE MOLG = -1. * FV**3 / TMP1 MOZG = MIN( (ZPD-Z0MG)/MOLG, 1.) END IF IF (MOZG < 0.) THEN FHGNEW = (1. - 15.*MOZG)**(-0.25) ELSE FHGNEW = 1.+ 4.7*MOZG ENDIF IF (ITER == 1) THEN FHG = FHGNEW ELSE FHG = 0.5 * (FHG+FHGNEW) ENDIF CWPC = (CWP * VAI * HCAN * FHG)**0.5 TMP1 = EXP( -CWPC*Z0HG/HCAN ) TMP2 = EXP( -CWPC*(Z0H+ZPD)/HCAN ) TMPRAH2 = HCAN*EXP(CWPC) / CWPC * (TMP1-TMP2) KH = MAX ( VKC*FV*(HCAN-ZPD), MPE ) RAMG = 0. RAHG = TMPRAH2 / KH RAWG = RAHG TMPRB = CWPC*50. / (1. - EXP(-CWPC/2.)) RB = TMPRB * SQRT(DLEAF(VEGTYP)/UC) END SUBROUTINE RAGRB SUBROUTINE SFCDIF1(ITER ,SFCTMP ,RHOAIR ,H ,QAIR , & & ZLVL ,ZPD ,Z0M ,Z0H ,UR , & & MPE ,ILOC ,JLOC , & & MOZ ,MOZSGN ,FM ,FH , & & CM ,CH ,FV ) IMPLICIT NONE INTEGER, INTENT(IN) :: ILOC INTEGER, INTENT(IN) :: JLOC INTEGER, INTENT(IN) :: ITER REAL, INTENT(IN) :: SFCTMP REAL, INTENT(IN) :: RHOAIR REAL, INTENT(IN) :: H REAL, INTENT(IN) :: QAIR REAL, INTENT(IN) :: ZLVL REAL, INTENT(IN) :: ZPD REAL, INTENT(IN) :: Z0H REAL, INTENT(IN) :: Z0M REAL, INTENT(IN) :: UR REAL, INTENT(IN) :: MPE INTEGER, INTENT(INOUT) :: MOZSGN REAL, INTENT(INOUT) :: MOZ REAL, INTENT(INOUT) :: FM REAL, INTENT(INOUT) :: FH REAL, INTENT(OUT) :: CM REAL, INTENT(OUT) :: CH REAL, INTENT(OUT) :: FV REAL :: MOL REAL :: TMPCM REAL :: TMPCH REAL :: FMNEW REAL :: FHNEW REAL :: MOZOLD REAL :: TMP1,TMP2,TMP3,TMP4,TMP5 REAL :: TVIR REAL :: CMFM, CHFH MOZOLD = MOZ IF(ZLVL <= ZPD) THEN write(*,*) 'critical problem: ZLVL <= ZPD; model stops' call wrf_error_fatal3("",4026,& "STOP in Noah-MP") ENDIF TMPCM = LOG((ZLVL-ZPD) / Z0M) TMPCH = LOG((ZLVL-ZPD) / Z0H) IF(ITER == 1) THEN FV = 0.0 MOZ = 0.0 MOL = 0.0 ELSE TVIR = (1. + 0.61*QAIR) * SFCTMP TMP1 = VKC * (GRAV/TVIR) * H/(RHOAIR*CPAIR) IF (ABS(TMP1) .LE. MPE) TMP1 = MPE MOL = -1. * FV**3 / TMP1 MOZ = MIN( (ZLVL-ZPD)/MOL, 1.) ENDIF IF (MOZOLD*MOZ .LT. 0.) MOZSGN = MOZSGN+1 IF (MOZSGN .GE. 2) THEN MOZ = 0. FM = 0. FH = 0. ENDIF IF (MOZ .LT. 0.) THEN TMP1 = (1. - 16.*MOZ)**0.25 TMP2 = LOG((1.+TMP1*TMP1)/2.) TMP3 = LOG((1.+TMP1)/2.) FMNEW = 2.*TMP3 + TMP2 - 2.*ATAN(TMP1) + 1.5707963 FHNEW = 2*TMP2 ELSE FMNEW = -5.*MOZ FHNEW = FMNEW ENDIF IF (ITER == 1) THEN FM = FMNEW FH = FHNEW ELSE FM = 0.5 * (FM+FMNEW) FH = 0.5 * (FH+FHNEW) ENDIF CMFM = TMPCM-FM CHFH = TMPCH-FH IF(ABS(CMFM) <= MPE) CMFM = MPE IF(ABS(CHFH) <= MPE) CHFH = MPE CM = VKC*VKC/(CMFM*CMFM) CH = VKC*VKC/(CMFM*CHFH) FV = UR * SQRT(CM) END SUBROUTINE SFCDIF1 SUBROUTINE SFCDIF2(ITER ,Z0 ,THZ0 ,THLM ,SFCSPD , & CZIL ,ZLM ,ILOC ,JLOC , & AKMS ,AKHS ,RLMO ,WSTAR2 , & USTAR ) IMPLICIT NONE INTEGER, INTENT(IN) :: ILOC INTEGER, INTENT(IN) :: JLOC INTEGER, INTENT(IN) :: ITER REAL, INTENT(IN) :: ZLM, Z0, THZ0, THLM, SFCSPD, CZIL REAL, intent(INOUT) :: AKMS REAL, intent(INOUT) :: AKHS REAL, intent(INOUT) :: RLMO REAL, intent(INOUT) :: WSTAR2 REAL, intent(OUT) :: USTAR REAL ZZ, PSLMU, PSLMS, PSLHU, PSLHS REAL XX, PSPMU, YY, PSPMS, PSPHU, PSPHS REAL ZILFC, ZU, ZT, RDZ, CXCH REAL DTHV, DU2, BTGH, ZSLU, ZSLT, RLOGU, RLOGT REAL ZETALT, ZETALU, ZETAU, ZETAT, XLU4, XLT4, XU4, XT4 REAL XLU, XLT, XU, XT, PSMZ, SIMM, PSHZ, SIMH, USTARK, RLMN, & & RLMA INTEGER ILECH, ITR INTEGER, PARAMETER :: ITRMX = 5 REAL, PARAMETER :: WWST = 1.2 REAL, PARAMETER :: WWST2 = WWST * WWST REAL, PARAMETER :: VKRM = 0.40 REAL, PARAMETER :: EXCM = 0.001 REAL, PARAMETER :: BETA = 1.0 / 270.0 REAL, PARAMETER :: BTG = BETA * GRAV REAL, PARAMETER :: ELFC = VKRM * BTG REAL, PARAMETER :: WOLD = 0.15 REAL, PARAMETER :: WNEW = 1.0 - WOLD REAL, PARAMETER :: PIHF = 3.14159265 / 2. REAL, PARAMETER :: EPSU2 = 1.E-4 REAL, PARAMETER :: EPSUST = 0.07 REAL, PARAMETER :: EPSIT = 1.E-4 REAL, PARAMETER :: EPSA = 1.E-8 REAL, PARAMETER :: ZTMIN = -5.0 REAL, PARAMETER :: ZTMAX = 1.0 REAL, PARAMETER :: HPBL = 1000.0 REAL, PARAMETER :: SQVISC = 258.2 REAL, PARAMETER :: RIC = 0.183 REAL, PARAMETER :: RRIC = 1.0 / RIC REAL, PARAMETER :: FHNEU = 0.8 REAL, PARAMETER :: RFC = 0.191 REAL, PARAMETER :: RFAC = RIC / ( FHNEU * RFC * RFC ) PSLMU (ZZ)= -0.96* log (1.0-4.5* ZZ) PSLMS (ZZ)= ZZ * RRIC -2.076* (1. -1./ (ZZ +1.)) PSLHU (ZZ)= -0.96* log (1.0-4.5* ZZ) PSLHS (ZZ)= ZZ * RFAC -2.076* (1. -1./ (ZZ +1.)) PSPMU (XX)= -2.* log ( (XX +1.)*0.5) - log ( (XX * XX +1.)*0.5) & & +2.* ATAN (XX) & &- PIHF PSPMS (YY)= 5.* YY PSPHU (XX)= -2.* log ( (XX * XX +1.)*0.5) PSPHS (YY)= 5.* YY ILECH = 0 ZILFC = - CZIL * VKRM * SQVISC ZU = Z0 RDZ = 1./ ZLM CXCH = EXCM * RDZ DTHV = THLM - THZ0 DU2 = MAX (SFCSPD * SFCSPD,EPSU2) BTGH = BTG * HPBL IF(ITER == 1) THEN IF (BTGH * AKHS * DTHV .ne. 0.0) THEN WSTAR2 = WWST2* ABS (BTGH * AKHS * DTHV)** (2./3.) ELSE WSTAR2 = 0.0 END IF USTAR = MAX (SQRT (AKMS * SQRT (DU2+ WSTAR2)),EPSUST) RLMO = ELFC * AKHS * DTHV / USTAR **3 END IF ZT = MAX(1.E-6,EXP (ZILFC * SQRT (USTAR * Z0))* Z0) ZSLU = ZLM + ZU ZSLT = ZLM + ZT RLOGU = log (ZSLU / ZU) RLOGT = log (ZSLT / ZT) ZETALT = MAX (ZSLT * RLMO,ZTMIN) RLMO = ZETALT / ZSLT ZETALU = ZSLU * RLMO ZETAU = ZU * RLMO ZETAT = ZT * RLMO IF (ILECH .eq. 0) THEN IF (RLMO .lt. 0.)THEN XLU4 = 1. -16.* ZETALU XLT4 = 1. -16.* ZETALT XU4 = 1. -16.* ZETAU XT4 = 1. -16.* ZETAT XLU = SQRT (SQRT (XLU4)) XLT = SQRT (SQRT (XLT4)) XU = SQRT (SQRT (XU4)) XT = SQRT (SQRT (XT4)) PSMZ = PSPMU (XU) SIMM = PSPMU (XLU) - PSMZ + RLOGU PSHZ = PSPHU (XT) SIMH = PSPHU (XLT) - PSHZ + RLOGT ELSE ZETALU = MIN (ZETALU,ZTMAX) ZETALT = MIN (ZETALT,ZTMAX) PSMZ = PSPMS (ZETAU) SIMM = PSPMS (ZETALU) - PSMZ + RLOGU PSHZ = PSPHS (ZETAT) SIMH = PSPHS (ZETALT) - PSHZ + RLOGT END IF ELSE IF (RLMO .lt. 0.)THEN PSMZ = PSLMU (ZETAU) SIMM = PSLMU (ZETALU) - PSMZ + RLOGU PSHZ = PSLHU (ZETAT) SIMH = PSLHU (ZETALT) - PSHZ + RLOGT ELSE ZETALU = MIN (ZETALU,ZTMAX) ZETALT = MIN (ZETALT,ZTMAX) PSMZ = PSLMS (ZETAU) SIMM = PSLMS (ZETALU) - PSMZ + RLOGU PSHZ = PSLHS (ZETAT) SIMH = PSLHS (ZETALT) - PSHZ + RLOGT END IF END IF USTAR = MAX (SQRT (AKMS * SQRT (DU2+ WSTAR2)),EPSUST) ZT = MAX(1.E-6,EXP (ZILFC * SQRT (USTAR * Z0))* Z0) ZSLT = ZLM + ZT RLOGT = log (ZSLT / ZT) USTARK = USTAR * VKRM AKMS = MAX (USTARK / SIMM,CXCH) AKHS = MAX (USTARK / SIMH,CXCH) IF (BTGH * AKHS * DTHV .ne. 0.0) THEN WSTAR2 = WWST2* ABS (BTGH * AKHS * DTHV)** (2./3.) ELSE WSTAR2 = 0.0 END IF RLMN = ELFC * AKHS * DTHV / USTAR **3 RLMA = RLMO * WOLD+ RLMN * WNEW RLMO = RLMA END SUBROUTINE SFCDIF2 SUBROUTINE SFCDIF3(ILOC ,JLOC ,TSK ,QS ,PSFC ,& PBLH ,Z0 ,Z0BASE ,VEGTYP ,ISURBAN,& IZ0TLND,SFCSPD ,ITER ,ITRMX ,TLOW ,& THLOW ,QLOW ,CWMLOW ,ZSL , & PLOW ,USTAR ,AKMS ,AKHS ,CHS2 ,& CQS2 ,RLMO ) USE MODULE_SF_MYJSFC, ONLY : & & EPSU2 , & & EPSUST , & & EPSZT , & & BETA , & & EXCML , & & RIC , & & SQVISC , & & ZTFC , & & BTG , & & CZIV , & & PI , & & PIHF , & & KZTM , & & KZTM2 , & & DZETA1 , & & DZETA2 , & & FH01 , & & FH02 , & & WWST2 , & & WWST , & & ZTMAX1 , & & ZTMAX2 , & & ZTMIN1 , & & ZTMIN2 , & & PSIH1 , & & PSIH2 , & & PSIM1 , & & PSIM2 USE MODULE_MODEL_CONSTANTS IMPLICIT NONE INTEGER,INTENT(IN) :: ILOC INTEGER,INTENT(IN) :: JLOC REAL ,INTENT(IN) :: TSK REAL ,INTENT(IN) :: PSFC REAL ,INTENT(IN) :: PBLH INTEGER,INTENT(IN) :: VEGTYP INTEGER,INTENT(IN) :: ISURBAN INTEGER,INTENT(IN) :: IZ0TLND REAL ,INTENT(IN) :: QLOW REAL ,INTENT(IN) :: THLOW REAL ,INTENT(IN) :: TLOW REAL ,INTENT(IN) :: CWMLOW REAL ,INTENT(IN) :: SFCSPD REAL ,INTENT(IN) :: PLOW REAL ,INTENT(IN) :: ZSL REAL ,INTENT(IN) :: Z0BASE INTEGER,INTENT(IN) :: ITER INTEGER,INTENT(IN) :: ITRMX REAL ,INTENT(OUT) :: CHS2 REAL ,INTENT(OUT) :: CQS2 REAL ,INTENT(OUT) :: RLMO REAL ,INTENT(INOUT) :: AKHS REAL ,INTENT(INOUT) :: AKMS REAL :: QZ0 REAL ,INTENT(INOUT) :: USTAR REAL ,INTENT(IN) :: Z0 REAL ,INTENT(INOUT):: QS REAL :: RIB INTEGER :: ITR,K REAL :: THZ0 REAL :: THVLOW REAL :: CT REAL :: BTGH REAL :: BTGX REAL :: CXCHL REAL :: DTHV REAL :: DU2 REAL :: ELFC REAL :: PSH02 REAL :: PSH10 REAL :: PSHZ REAL :: PSHZL REAL :: PSM10 REAL :: PSMZ REAL :: PSMZL REAL :: RDZ REAL :: RDZT REAL :: RLMA REAL :: RLMN REAL :: RLOGT REAL :: RLOGU REAL :: RZ REAL :: SIMH REAL :: SIMM REAL :: USTARK REAL :: WSTAR2 REAL :: WSTAR REAL :: CHS REAL :: RZSU REAL :: RZST REAL :: X,XLT,XLT4,XLU,XLU4,XT,XT4,XU,XU4,ZETALT,ZETALU , & ZETAT,ZETAU,ZQ,ZSLT,ZSLU,ZT,ZU,TOPOTERM,ZZIL REAL :: AKHS02,AKHS10,AKMS02,AKMS10 REAL :: ZU10 REAL :: ZT02 REAL :: ZT10 REAL :: RLNU10 REAL :: RLNT02 REAL :: RLNT10 REAL :: ZTAU10 REAL :: ZTAT02 REAL :: ZTAT10 REAL :: SIMM10 REAL :: SIMH02 REAL :: SIMH10 REAL :: ZUUZ REAL :: EKMS10 REAL :: test REAL :: E1 REAL, PARAMETER :: VKRM = 0.40 REAL, PARAMETER :: CZETMAX = 10. REAL :: CZIL REAL :: ZILFC THVLOW = THLOW*(1.+EP_1*QLOW) THZ0 = TSK*(P1000mb/PSFC)**RCP ZU = Z0 ZT = ZU*ZTFC ZQ = ZT QZ0 = QS RDZ = 1./ZSL CXCHL = EXCML*RDZ DTHV = THVLOW-THZ0*(0.608*QZ0+1.) BTGX=GRAV/THLOW ELFC=VKRM*BTGX IF(PBLH > 1000.)THEN BTGH = BTGX*PBLH ELSE BTGH = BTGX*1000. ENDIF DU2 = MAX(SFCSPD*SFCSPD,EPSU2) RIB = BTGX*DTHV*ZSL/DU2 ZSLU = ZSL+ZU RZSU = ZSLU/ZU RLOGU = LOG(RZSU) ZSLT = ZSL + ZU IF ( (IZ0TLND==0) .or. (VEGTYP == ISURBAN) ) THEN CZIL = 0.1 ELSE CZIL = 10.0 ** ( -0.40 * ( Z0 / 0.07 ) ) ENDIF ZILFC=-CZIL*VKRM*SQVISC IF(DTHV>0.)THEN IF (RIB",8723,& 'REDPRM: Error: too many input soil types') END IF IF (VEGTYP .gt. LUCATS) THEN call wrf_message('VEGTYP must be less than LUCATS:') write(message, '("VEGTYP = ", I6, "; LUCATS = ", I6)') VEGTYP, LUCATS call wrf_message(trim(message)) call wrf_error_fatal3("",8730,& 'Error: too many input landuse types') END IF CSOIL = CSOIL_DATA BEXP = BB (SOILTYP) DKSAT = SATDK (SOILTYP) DWSAT = SATDW (SOILTYP) F1 = F11 (SOILTYP) PSISAT = SATPSI (SOILTYP) QUARTZ = QTZ (SOILTYP) SMCDRY = DRYSMC (SOILTYP) SMCMAX = MAXSMC (SOILTYP) SMCREF = REFSMC (SOILTYP) SMCWLT = WLTSMC (SOILTYP) IF(VEGTYP==ISURBAN)THEN SMCMAX = 0.45 SMCREF = 0.42 SMCWLT = 0.40 SMCDRY = 0.40 CSOIL = 3.E6 ENDIF ZBOT = ZBOT_DATA CZIL = CZIL_DATA FRZK = FRZK_DATA REFDK = REFDK_DATA REFKDT = REFKDT_DATA KDT = REFKDT * DKSAT / REFDK SLOPE = SLOPE_DATA (SLOPETYP) if(SOILTYP /= 14) then FRZFACT = (SMCMAX / SMCREF) * (0.412 / 0.468) FRZX = FRZK * FRZFACT end if TOPT = TOPT_DATA RGL = RGLTBL (VEGTYP) RSMAX = RSMAX_DATA RSMIN = RSTBL (VEGTYP) HS = HSTBL (VEGTYP) NROOT = NROTBL (VEGTYP) IF(VEGTYP==ISURBAN)THEN RSMIN=400.0 ENDIF IF (NROOT .gt. NSOIL) THEN WRITE (*,*) 'Warning: too many root layers' write (*,*) 'NROOT = ', nroot write (*,*) 'NSOIL = ', nsoil call wrf_error_fatal3("",8799,& "STOP in Noah-MP") END IF END SUBROUTINE REDPRM SUBROUTINE EPARM(ILOC ,JLOC ,TAH ,TGB ,FVEG , & CHV ,CHB ,VEG ,CHSTAR ,TSTAR) IMPLICIT NONE INTEGER, INTENT(IN) :: ILOC INTEGER, INTENT(IN) :: JLOC REAL , INTENT(IN) :: TAH REAL , INTENT(IN) :: TGB REAL , INTENT(IN) :: FVEG REAL , INTENT(IN) :: CHV REAL , INTENT(IN) :: CHB LOGICAL, INTENT(IN) :: VEG REAL, INTENT(OUT) :: CHSTAR REAL, INTENT(OUT) :: TSTAR REAL :: W IF (VEG) THEN CHSTAR = FVEG*CHV +(1.-FVEG)*CHB W = FVEG*CHV*TAH+(1.-FVEG)*CHB*TGB TSTAR = W/CHSTAR ELSE CHSTAR = CHB TSTAR = TGB ENDIF END SUBROUTINE EPARM subroutine noahmp_options(idveg ,iopt_crs ,iopt_btr ,iopt_run ,iopt_sfc ,iopt_frz , & iopt_inf ,iopt_rad ,iopt_alb ,iopt_snf ,iopt_tbot, iopt_stc ) implicit none INTEGER, INTENT(IN) :: idveg INTEGER, INTENT(IN) :: iopt_crs INTEGER, INTENT(IN) :: iopt_btr INTEGER, INTENT(IN) :: iopt_run INTEGER, INTENT(IN) :: iopt_sfc INTEGER, INTENT(IN) :: iopt_frz INTEGER, INTENT(IN) :: iopt_inf INTEGER, INTENT(IN) :: iopt_rad INTEGER, INTENT(IN) :: iopt_alb INTEGER, INTENT(IN) :: iopt_snf INTEGER, INTENT(IN) :: iopt_tbot INTEGER, INTENT(IN) :: iopt_stc dveg = idveg opt_crs = iopt_crs opt_btr = iopt_btr opt_run = iopt_run opt_sfc = iopt_sfc opt_frz = iopt_frz opt_inf = iopt_inf opt_rad = iopt_rad opt_alb = iopt_alb opt_snf = iopt_snf opt_tbot = iopt_tbot opt_stc = iopt_stc end subroutine noahmp_options END MODULE NOAHMP_ROUTINES MODULE MODULE_SF_NOAHMPLSM USE NOAHMP_ROUTINES USE NOAHMP_GLOBALS USE NOAHMP_VEG_PARAMETERS END MODULE MODULE_SF_NOAHMPLSM