diff --git a/README.md b/README.md index 18a0e752..912568f7 100644 --- a/README.md +++ b/README.md @@ -1,39 +1,5 @@ -# noahmp -Noah-MP Community Repository +# Noah-MP version 4.7 release -This is the official Noah-MP unified Github repository for code downloading and contribution. Note that Noah-MP is a community model contributed by the whole Noah-MP scientific community. For maintenance and release of this GitHub, please contact: Cenlin He (cenlinhe@ucar.edu) and Fei Chen (feichen@ucar.edu). +This is the official Noah-MP code version 4.7 consistent with that released in WRF v4.7. Note that WRF v4.7 GitHub code is directly connected to this Noah-MP GitHub through the submodule mechanism. -Some changes have been made to the structure of archiving the stand-alone version of Noah-MP/HRLDAS code in the Github repository. Now, it separately archives the core Noah-MP source code, Noah-MP driver, and parameter tables in this unified Noah-MP Github repository and the rest of the HRLDAS related files (e.g., module_sf_urban.F, etc.) in another unified HRLDAS Github repository (https://github.com/NCAR/hrldas). The HRLDAS Github repo is already linked to this unified core Noah-MP code repository. This new archiving structure will allow different host model platforms/systems (e.g., HRLDAS, WRF, UFS, NWM, LIS, etc.) to connect to the core Noah-MP source code and develop their own host model drivers. - -Model developers can make code development based on the develop branch and create pull request to the develop branch. The pull request will be reviewed by Noah-MP model release team and if everything looks good, the new code development will be merged to the develop branch. Eventually, the updates in the develop branch will be merged to the master branch for official Noah-MP model release. - -Branch structures of this Noah-MP repository: - -1. "master" branch: most stable & recent version, updated whenever there are bug fixes or major model update/release; - -2. "develop" branch: used for continuous NoahMP development, keep being updated by including bug fixes and code updates (e.g., new physics options, processes, etc.); - -3. "develop_no_gecros" branch: same as the "develop" branch, except for excluding the gecros crop section (to be consistent with recent Noah-MP changes for National Water Model (NWM)); - -4. "release-v4.3-WRF" branch: version archive of the stable version consistent with the WRFv4.3 release; - -5. "release-v4.3-NWM" branch: version archive, same as the "release-v4.3-WRF" branch, except for excluding the gecros crop section (to be consistent with recent Noah-MP changes for National Water Model (NWM)); - -Some suggestions for model developers to contribute to Noah-MP code through the Github repository (typical procedures): - -1. Step (1) Create a fork of this official NoahMP repository to your own Github account; - -2. Step (2) Make code updates based on the "develop" branch of the forked repository under your own account; - -3. Step (3) Finalize and test the code updates you make; - -4. Step (4) Submit a pull request for your model updates from your own forked Github repository to the "develop" branch of the official repository; - -5. Step (5) The Noah-MP release team reviews and tests the model updates in the submitted pull request and discusses with the developer if there is any problem; - -6. Step (6) The Noah-MP release team confirms the pull request and merges the updated code to the "develop" and "develop_no_gecros" branches in the official repository; - -7. Step (7) The Noah-MP release team will merge the updated "develop" branch to the master branch and version-release branch during the annual model release. - -Note: If updates are made to both the core NoahMP source codes and other HRLDAS files, two separate pull requests need to be submitted to this NoahMP repository and the HRLDAS repository (https://github.com/NCAR/hrldas), respectively, regarding the changes in the Noah-MP code files and other HRLDAS-related files. This could be done by using the same titles for the two pull requests to ensure that the submitted code changes are handled together by the release team in the two repositories. diff --git a/drivers/wrf/module_sf_noahmpdrv.F b/drivers/wrf/module_sf_noahmpdrv.F index 46e61a17..ab47e023 100644 --- a/drivers/wrf/module_sf_noahmpdrv.F +++ b/drivers/wrf/module_sf_noahmpdrv.F @@ -645,7 +645,6 @@ SUBROUTINE noahmplsm(ITIMESTEP, YR, JULIAN, COSZIN,XLAT,XLONG, & ! IN IPRINT = .false. ! debug printout ! for using soil update timestep difference from noahmp main timestep - calculate_soil = .false. soil_update_steps = nint(soiltstep/DT) ! 3600 = 1 hour soil_update_steps = max(soil_update_steps,1) if ( soil_update_steps == 1 ) then @@ -673,7 +672,7 @@ SUBROUTINE noahmplsm(ITIMESTEP, YR, JULIAN, COSZIN,XLAT,XLONG, & ! IN end if end if - if (mod(itimestep,soil_update_steps) == 0) calculate_soil = .true. + calculate_soil = mod(itimestep,soil_update_steps) == 0 ! end soil timestep YEARLEN = 365 ! find length of year for phenology (also S Hemisphere) @@ -1056,7 +1055,7 @@ SUBROUTINE noahmplsm(ITIMESTEP, YR, JULIAN, COSZIN,XLAT,XLONG, & ! IN FSA, FSR, FIRA, FSH, FGEV, SSOIL, & ! OUT : TRAD, ESOIL, RUNSF, RUNSB, SAG, SALB, & ! OUT : QSNBOT,PONDING,PONDING1,PONDING2, T2MB, Q2MB, & ! OUT : - EMISSI, FPICE, CHB2, QMELT & ! OUT : + EMISSI, FPICE, CHB2, QMELT, HCPCT & ! OUT : #ifdef WRF_HYDRO , sfcheadrt(i,j) & #endif @@ -2004,7 +2003,7 @@ SUBROUTINE NOAHMP_INIT ( MMINLU, SNOW , SNOWH , CANWAT , ISLTYP , IVGTYP, XLAT call read_mp_crop_parameters() call read_tiledrain_parameters() call read_mp_optional_parameters() - if(iopt_irr >= 1) call read_mp_irrigation_parameters() + call read_mp_irrigation_parameters(iopt_irr) IF( .NOT. restart ) THEN @@ -2181,7 +2180,11 @@ SUBROUTINE NOAHMP_INIT ( MMINLU, SNOW , SNOWH , CANWAT , ISLTYP , IVGTYP, XLAT lai (I,J) = max(lai(i,j),0.05) ! at least start with 0.05 for arbitrary initialization (v3.7) xsaixy (I,J) = max(0.1*lai(I,J),0.05) ! MB: arbitrarily initialize SAI using input LAI (v3.7) - masslai = 1000. / max(SLA_TABLE(IVGTYP(I,J)),1.0) ! conversion from lai to mass (v3.7) + if (urbanpt_flag) then + masslai = 1000. / max(SLA_TABLE(NATURAL_TABLE),1.0) ! conversion from lai to mass (v3.7) + else + masslai = 1000. / max(SLA_TABLE(IVGTYP(I,J)),1.0) ! conversion from lai to mass (v3.7) + endif lfmassxy (I,J) = lai(i,j)*masslai ! use LAI to initialize (v3.7) masssai = 1000. / 3.0 ! conversion from lai to mass (v3.7) stmassxy (I,J) = xsaixy(i,j)*masssai ! use SAI to initialize (v3.7) @@ -3035,6 +3038,9 @@ SUBROUTINE noahmp_urban(sf_urban_physics, NSOIL, IVGTYP, ITIMESTEP, dgr_urb3d, dg_urb3d, lfr_urb3d, lfg_urb3d, & !RMS lp_urb2d, hi_urb2d, lb_urb2d, hgt_urb2d, & !H multi-layer urban mh_urb2d, stdh_urb2d, lf_urb2d, & !SLUCM + lf_urb2d_s, z0_urb2d, vegfra, & !SLUCM + albr_urb2d, albb_urb2d, albg_urb2d, & !I urban explicit radiation PK2025 + epsr_urb2d, epsb_urb2d, epsg_urb2d, & !I urban explicit radiation PK2025 th_phy, rho, p_phy, ust, & !I multi-layer urban gmt, julday, xlong, xlat, & !I multi-layer urban a_u_bep, a_v_bep, a_t_bep, a_q_bep, & !O multi-layer urban @@ -3172,6 +3178,15 @@ SUBROUTINE noahmp_urban(sf_urban_physics, NSOIL, IVGTYP, ITIMESTEP, REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_roof_layers, jms:jme ), INTENT(INOUT) :: TGRL_URB3D REAL, OPTIONAL, DIMENSION( ims:ime, 1:num_roof_layers, jms:jme ), INTENT(INOUT) :: SMR_URB3D +! Explicit radiation parameters PK2025 + real :: albr_urb,albb_urb,albg_urb,epsr_urb,epsb_urb,epsg_urb + real, optional, dimension( ims:ime, jms:jme ), intent(in) :: albr_urb2d + real, optional, dimension( ims:ime, jms:jme ), intent(in) :: albb_urb2d + real, optional, dimension( ims:ime, jms:jme ), intent(in) :: albg_urb2d + real, optional, dimension( ims:ime, jms:jme ), intent(in) :: epsr_urb2d + real, optional, dimension( ims:ime, jms:jme ), intent(in) :: epsb_urb2d + real, optional, dimension( ims:ime, jms:jme ), intent(in) :: epsg_urb2d + ! state variable surface_driver <--> lsm <--> urban @@ -3239,6 +3254,11 @@ SUBROUTINE noahmp_urban(sf_urban_physics, NSOIL, IVGTYP, ITIMESTEP, REAL :: hgt_urb REAL, DIMENSION(4) :: lf_urb +! Distributed aerodynamics parameters + REAL :: lf_urb_s + REAL :: z0_urb + REAL :: vegfrac + ! Local variables INTEGER :: I,J,K @@ -3277,6 +3297,9 @@ SUBROUTINE noahmp_urban(sf_urban_physics, NSOIL, IVGTYP, ITIMESTEP, REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: mh_urb2d REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: stdh_urb2d REAL, OPTIONAL, DIMENSION( ims:ime, 4, jms:jme ), INTENT(IN ) :: lf_urb2d + REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: lf_urb2d_s + REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: z0_urb2d + REAL, OPTIONAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: vegfra REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zrd, jms:jme ), INTENT(INOUT) :: trb_urb4d REAL, OPTIONAL, DIMENSION( ims:ime, 1:urban_map_zwd, jms:jme ), INTENT(INOUT) :: tw1_urb4d @@ -3472,6 +3495,18 @@ SUBROUTINE noahmp_urban(sf_urban_physics, NSOIL, IVGTYP, ITIMESTEP, IF (I.EQ.73.AND.J.EQ.125)THEN CHECK = 1 END IF +! Distributed aerodynamics + lf_urb_s = lf_urb2d_s(I, J) + z0_urb = z0_urb2d(I, J) + vegfrac = vegfra(I, J) / 100 + +! Explicit radiation parameters PK2025 + albr_urb = albr_urb2d(i,j) + albb_urb = albb_urb2d(i,j) + albg_urb = albg_urb2d(i,j) + epsr_urb = epsr_urb2d(i,j) + epsb_urb = epsb_urb2d(i,j) + epsg_urb = epsg_urb2d(i,j) ! Call urban @@ -3493,10 +3528,14 @@ SUBROUTINE noahmp_urban(sf_urban_physics, NSOIL, IVGTYP, ITIMESTEP, CMR_URB, CHR_URB, CMC_URB, CHC_URB, & U10_URB, V10_URB, TH2_URB, Q2_URB, & ! O UST_URB, mh_urb, stdh_urb, lf_urb, lp_urb, & ! 0 - hgt_urb, frc_urb, lb_urb, check, CMCR_URB,TGR_URB, & ! H + hgt_urb, frc_urb, lb_urb, & ! H + albr_urb, albb_urb, albg_urb, & ! I + epsr_urb, epsb_urb, epsg_urb, & ! I + check, CMCR_URB, TGR_URB, & ! H TGRL_URB, SMR_URB, CMGR_URB, CHGR_URB, jmonth, & ! H DRELR_URB, DRELB_URB, & ! H - DRELG_URB,FLXHUMR_URB,FLXHUMB_URB,FLXHUMG_URB ) + DRELG_URB,FLXHUMR_URB,FLXHUMB_URB,FLXHUMG_URB, & + lf_urb_s, z0_urb, vegfrac) TS_URB2D(I,J) = TS_URB @@ -3505,7 +3544,7 @@ SUBROUTINE noahmp_urban(sf_urban_physics, NSOIL, IVGTYP, ITIMESTEP, QFX(I,J) = FRC_URB2D(I,J) * LH_KINEMATIC_URB & + (1-FRC_URB2D(I,J))* QFX(I,J) ![kg/m/m/s] LH(I,J) = FRC_URB2D(I,J) * LH_URB + (1-FRC_URB2D(I,J)) * LH(I,J) ![W/m/m] - GRDFLX(I,J) = FRC_URB2D(I,J) * (G_URB) + (1-FRC_URB2D(I,J)) * GRDFLX(I,J) ![W/m/m] + GRDFLX(I,J) = FRC_URB2D(I,J) * (G_URB*(-1.0)) + (1-FRC_URB2D(I,J)) * GRDFLX(I,J) ![W/m/m] positive: downward direction TSK(I,J) = FRC_URB2D(I,J) * TS_URB + (1-FRC_URB2D(I,J)) * TSK(I,J) ![K] ! Q1 = QSFC(I,J)/(1.0+QSFC(I,J)) ! Q1 = FRC_URB2D(I,J) * QS_URB + (1-FRC_URB2D(I,J)) * Q1 ![-] @@ -3734,7 +3773,7 @@ SUBROUTINE noahmp_urban(sf_urban_physics, NSOIL, IVGTYP, ITIMESTEP, ! rename *_urb to sh_urb2d,lh_urb2d,g_urb2d,rn_urb2d - grdflx(i,j) = (1.-frc_urb2d(i,j))*grdflx_rural(i,j)+ frc_urb2d(i,j)*grdflx_urb(i,j) + grdflx(i,j) = (1.-frc_urb2d(i,j))*grdflx_rural(i,j)+ frc_urb2d(i,j)*grdflx_urb(i,j)*(-1.0) ! positive: downward direction qfx(i,j) = (1.-frc_urb2d(i,j))*qfx_rural(i,j) + qfx_urb(i,j) lh(i,j) = qfx(i,j)*xlv hfx(i,j) = hfx_urb(i,j) + (1-frc_urb2d(i,j))*hfx_rural(i,j) ![W/m/m] diff --git a/src/module_sf_noahmp_glacier.F b/src/module_sf_noahmp_glacier.F index 5ec65e58..df1b09d2 100644 --- a/src/module_sf_noahmp_glacier.F +++ b/src/module_sf_noahmp_glacier.F @@ -112,7 +112,7 @@ SUBROUTINE NOAHMP_GLACIER (& FSA ,FSR ,FIRA ,FSH ,FGEV ,SSOIL , & ! OUT : TRAD ,EDIR ,RUNSRF ,RUNSUB ,SAG ,ALBEDO , & ! OUT : QSNBOT ,PONDING ,PONDING1,PONDING2,T2M ,Q2E , & ! OUT : - EMISSI, FPICE, CH2B ,QMELT & ! OUT : + EMISSI, FPICE, CH2B ,QMELT ,HCPCT & ! OUT : #ifdef WRF_HYDRO , sfcheadrt & #endif @@ -191,6 +191,8 @@ SUBROUTINE NOAHMP_GLACIER (& REAL , INTENT(OUT) :: EMISSI REAL , INTENT(OUT) :: FPICE REAL , INTENT(OUT) :: CH2B + REAL, DIMENSION(-NSNOW+1:NSOIL), INTENT(OUT) :: HCPCT !heat capacity [j/m3/k] + REAL , INTENT(OUT) :: QMELT !internal pack melt due to phase change [mm/s] ! local INTEGER :: IZ !do-loop index @@ -209,7 +211,6 @@ SUBROUTINE NOAHMP_GLACIER (& REAL :: QDEW !ground surface dew rate [mm/s] REAL :: QVAP !ground surface evap. rate [mm/s] REAL :: LATHEA !latent heat [j/kg] - REAL, INTENT(OUT) :: QMELT !internal pack melt due to phase change [mm/s] REAL :: SWDOWN !downward solar [w/m2] REAL :: BEG_WB !beginning water for error check REAL :: ZBOT = -8.0 @@ -245,7 +246,8 @@ SUBROUTINE NOAHMP_GLACIER (& TAUSS ,QSFC , & !inout IMELT ,SNICEV ,SNLIQV ,EPORE ,QMELT ,PONDING, & !out SAG ,FSA ,FSR ,FIRA ,FSH ,FGEV , & !out - TRAD ,T2M ,SSOIL ,LATHEA ,Q2E ,EMISSI, CH2B ) !out + TRAD ,T2M ,SSOIL ,LATHEA ,Q2E ,EMISSI, & !out + CH2B ,HCPCT ) !out SICE = MAX(0.0, SMC - SH2O) SNEQVO = SNEQV @@ -358,7 +360,7 @@ SUBROUTINE ENERGY_GLACIER (NSNOW ,NSOIL ,ISNOW ,DT ,QSNOW ,RHOAIR , & !i TAUSS ,QSFC , & !inout IMELT ,SNICEV ,SNLIQV ,EPORE ,QMELT ,PONDING, & !out SAG ,FSA ,FSR ,FIRA ,FSH ,FGEV , & !out - TRAD ,T2M ,SSOIL ,LATHEA ,Q2E ,EMISSI, CH2B ) !out + TRAD ,T2M ,SSOIL ,LATHEA ,Q2E,EMISSI, CH2B, HCPCT ) !out ! -------------------------------------------------------------------------------------------------- ! -------------------------------------------------------------------------------------------------- @@ -426,6 +428,7 @@ SUBROUTINE ENERGY_GLACIER (NSNOW ,NSOIL ,ISNOW ,DT ,QSNOW ,RHOAIR , & !i REAL , INTENT(OUT) :: Q2E REAL , INTENT(OUT) :: EMISSI REAL , INTENT(OUT) :: CH2B !sensible heat conductance, canopy air to ZLVL air (m/s) + REAL, DIMENSION(-NSNOW+1:NSOIL) , INTENT(OUT) :: HCPCT !heat capacity [j/m3/k] ! local @@ -438,7 +441,6 @@ SUBROUTINE ENERGY_GLACIER (NSNOW ,NSOIL ,ISNOW ,DT ,QSNOW ,RHOAIR , & !i REAL :: FIRE !emitted IR (w/m2) REAL, DIMENSION(-NSNOW+1:NSOIL) :: FACT !temporary used in phase change REAL, DIMENSION(-NSNOW+1:NSOIL) :: DF !thermal conductivity [w/m/k] - REAL, DIMENSION(-NSNOW+1:NSOIL) :: HCPCT !heat capacity [j/m3/k] REAL :: GAMMA !psychrometric constant (pa/k) REAL :: RHSUR !raltive humidity in surface soil/snow air space (-) @@ -711,26 +713,25 @@ SUBROUTINE RADIATION_GLACIER (DT ,TG ,SNEQVO ,SNEQV ,COSZ , & !i ! compute snow albedo only if cosz >0 ! to be consistent with the main NoahMP code, currently do not include snow aging when sun is not present ! this needs more future work - if (COSZ > 0.0) then - ! snow age + ! snow age (allow nighttime aging) - CALL SNOW_AGE_GLACIER (DT,TG,SNEQVO,SNEQV,TAUSS,FAGE) + CALL SNOW_AGE_GLACIER (DT,TG,SNEQVO,SNEQV,TAUSS,FAGE) + if (cosz > 0.0) then ! only compute albedo when there is sunlight IF(OPT_ALB == 1) & CALL SNOWALB_BATS_GLACIER (NBAND,COSZ,FAGE,ALBSND,ALBSNI) IF(OPT_ALB == 2) THEN CALL SNOWALB_CLASS_GLACIER(NBAND,QSNOW,DT,ALB,ALBOLD,ALBSND,ALBSNI) ALBOLD = ALB END IF - endif ! zero summed solar fluxes - SAG = 0. - FSA = 0. - FSR = 0. + SAG = 0.0 + FSA = 0.0 + FSR = 0.0 FSNO = 0.0 IF(SNEQV > 0.0) FSNO = 1.0 diff --git a/src/module_sf_noahmplsm.F b/src/module_sf_noahmplsm.F index b2588c81..4e9659a0 100644 --- a/src/module_sf_noahmplsm.F +++ b/src/module_sf_noahmplsm.F @@ -1447,7 +1447,7 @@ SUBROUTINE PRECIP_HEAT (parameters,ILOC ,JLOC ,VEGTYP ,DT ,UU ,VV ! --------------------------- liquid water ------------------------------ ! maximum canopy water - MAXLIQ = parameters%CH2OP * (ELAI+ ESAI) + MAXLIQ = FVEG * parameters%CH2OP * (ELAI+ ESAI) ! average interception and throughfall @@ -1477,7 +1477,7 @@ SUBROUTINE PRECIP_HEAT (parameters,ILOC ,JLOC ,VEGTYP ,DT ,UU ,VV ! --------------------------- canopy ice ------------------------------ ! for canopy ice - MAXSNO = 6.6*(0.27+46.0/BDFALL) * (ELAI+ ESAI) + MAXSNO = FVEG * 6.6*(0.27+46.0/BDFALL) * (ELAI+ ESAI) IF((ELAI+ ESAI) .GT. 0.0) THEN QINTS = FVEG * SNOW * FP @@ -2920,7 +2920,11 @@ SUBROUTINE ALBEDO (parameters,VEGTYP ,IST ,ICE ,NSOIL , & !in IF (IB.EQ.1) FSUN = 0.0 END DO - IF(COSZ <= 0) GOTO 100 +! snow age (allow nighttime aging) + + CALL SNOW_AGE (parameters,DT,TG,SNEQVO,SNEQV,TAUSS,FAGE) + + IF (cosz > 0) THEN ! weight reflectance/transmittance by LAI and SAI @@ -2932,10 +2936,6 @@ SUBROUTINE ALBEDO (parameters,VEGTYP ,IST ,ICE ,NSOIL , & !in TAU(IB) = MAX(parameters%TAUL(IB)*WL+parameters%TAUS(IB)*WS, MPE) END DO -! snow age - - CALL SNOW_AGE (parameters,DT,TG,SNEQVO,SNEQV,TAUSS,FAGE) - ! snow albedos: only if COSZ > 0 and FSNO > 0 IF(OPT_ALB == 1) & @@ -2985,7 +2985,7 @@ SUBROUTINE ALBEDO (parameters,VEGTYP ,IST ,ICE ,NSOIL , & !in END IF FSUN = WL -100 CONTINUE + END IF ! cosz>0 END SUBROUTINE ALBEDO @@ -3935,21 +3935,21 @@ SUBROUTINE VEGE_FLUX(parameters,NSNOW ,NSOIL ,ISNOW ,VEGTYP ,VEG , & IF(ITER == 1) THEN IF (OPT_CRS == 1) then ! Ball-Berry CALL STOMATA (parameters,VEGTYP,MPE ,PARSUN ,FOLN ,ILOC , JLOC , & !in - TV ,ESTV ,EAH ,SFCTMP,SFCPRS, & !in + TV ,ESTV ,EAH ,SFCTMP,SFCPRS, FVEG, & !in O2AIR ,CO2AIR,IGS ,BTRAN ,RB , & !in RSSUN ,PSNSUN) !out CALL STOMATA (parameters,VEGTYP,MPE ,PARSHA ,FOLN ,ILOC , JLOC , & !in - TV ,ESTV ,EAH ,SFCTMP,SFCPRS, & !in + TV ,ESTV ,EAH ,SFCTMP,SFCPRS, FVEG, & !in O2AIR ,CO2AIR,IGS ,BTRAN ,RB , & !in RSSHA ,PSNSHA) !out END IF IF (OPT_CRS == 2) then ! Jarvis - CALL CANRES (parameters,PARSUN,TV ,BTRAN ,EAH ,SFCPRS, & !in + CALL CANRES (parameters,PARSUN,TV ,BTRAN ,EAH ,SFCPRS, FVEG, & !in RSSUN ,PSNSUN,ILOC ,JLOC ) !out - CALL CANRES (parameters,PARSHA,TV ,BTRAN ,EAH ,SFCPRS, & !in + CALL CANRES (parameters,PARSHA,TV ,BTRAN ,EAH ,SFCPRS, FVEG, & !in RSSHA ,PSNSHA,ILOC ,JLOC ) !out END IF @@ -4052,14 +4052,14 @@ SUBROUTINE VEGE_FLUX(parameters,NSNOW ,NSOIL ,ISNOW ,VEGTYP ,VEG , & B = SAV-IRC-SHC-EVC-TR+PAHV !additional w/m2 ! A = FVEG*(4.0*CIR*TV**3 + CSH + (CEV+CTR)*DESTV) !volumetric heat capacity - A = FVEG*(4.0*CIR*TV**3 + CSH + (CEV+CTR)*DESTV) + HCV/DT ! total volumetric heat capacity, add canopy heat capacity, more stable + A = FVEG*(4.0*CIR*TV**3 + CSH + (CEV+CTR)*DESTV + HCV/DT) ! total volumetric heat capacity, add canopy heat capacity, more stable DTV = B/A IRC = IRC + FVEG*4.0*CIR*TV**3*DTV SHC = SHC + FVEG*CSH*DTV EVC = EVC + FVEG*CEV*DESTV*DTV TR = TR + FVEG*CTR*DESTV*DTV - CANHS = DTV * HCV/DT ! w/m2 canopy heat storage change + CANHS = DTV * FVEG*HCV/DT ! w/m2 canopy heat storage change ! update vegetation surface temperature TV = TV + DTV @@ -4808,9 +4808,9 @@ SUBROUTINE SFCDIF2(parameters,ITER ,Z0 ,THZ0 ,THLM ,SFCSPD , & !in ! ---------------------------------------------------------------------- ! LECH'S SURFACE FUNCTIONS PSLMU (ZZ)= -0.96* log (1.0-4.5* ZZ) - PSLMS (ZZ)= ZZ * RRIC -2.076* (1.0 -1.0/ (ZZ +1.0)) + PSLMS (ZZ)= (ZZ / RFC) -2.076* (1.0 -1.0/ (ZZ +1.0)) PSLHU (ZZ)= -0.96* log (1.0-4.5* ZZ) - PSLHS (ZZ)= ZZ * RFAC -2.076* (1.0 -1.0/ (ZZ +1.0)) + PSLHS (ZZ)= ZZ * RFAC -2.076* (1.0 - exp(-1.2 * ZZ)) ! PAULSON'S SURFACE FUNCTIONS PSPMU (XX)= -2.0* log ( (XX +1.0)*0.5) - log ( (XX * XX +1.0)*0.5) & & +2.0* ATAN (XX) & @@ -5003,7 +5003,7 @@ END SUBROUTINE ESAT !== begin stomata ================================================================================== SUBROUTINE STOMATA (parameters,VEGTYP ,MPE ,APAR ,FOLN ,ILOC , JLOC, & !in - TV ,EI ,EA ,SFCTMP ,SFCPRS , & !in + TV ,EI ,EA ,SFCTMP ,SFCPRS , FVEG, & !in O2 ,CO2 ,IGS ,BTRAN ,RB , & !in RS ,PSN ) !out ! -------------------------------------------------------------------------------------------------- @@ -5017,7 +5017,7 @@ SUBROUTINE STOMATA (parameters,VEGTYP ,MPE ,APAR ,FOLN ,ILOC , JLO REAL, INTENT(IN) :: IGS !growing season index (0=off, 1=on) REAL, INTENT(IN) :: MPE !prevents division by zero errors - + REAL, INTENT(IN) :: FVEG !vegetation greeness fraction REAL, INTENT(IN) :: TV !foliage temperature (k) REAL, INTENT(IN) :: EI !vapor pressure inside leaf (sat vapor press at tv) (pa) REAL, INTENT(IN) :: EA !vapor pressure of canopy air (pa) @@ -5067,6 +5067,8 @@ SUBROUTINE STOMATA (parameters,VEGTYP ,MPE ,APAR ,FOLN ,ILOC , JLO REAL :: J !electron transport (umol co2/m2/s) REAL :: CEA !constrain ea or else model blows up REAL :: CF !s m2/umol -> s/m + REAL :: APAR_scale !APAR scaled by FVEG from grid level to stand scale + F1(AB,BC) = AB**((BC-25.0)/10.0) F2(AB) = 1.0 + EXP((-2.2E05+710.0*(AB+273.16))/(8.314*(AB+273.16))) @@ -5076,15 +5078,16 @@ SUBROUTINE STOMATA (parameters,VEGTYP ,MPE ,APAR ,FOLN ,ILOC , JLO ! initialize RS=RSMAX and PSN=0 because will only do calculations ! for APAR > 0, in which case RS <= RSMAX and PSN >= 0 + APAR_scale = APAR / max(FVEG,1.0e-6) ! scaling APAR back to stand scale CF = SFCPRS/(8.314*SFCTMP)*1.0e06 RS = 1.0/parameters%BP * CF PSN = 0.0 - IF (APAR .LE. 0.0) RETURN + IF (APAR_scale .LE. 0.0) RETURN FNF = MIN( FOLN/MAX(MPE,parameters%FOLNMX), 1.0 ) TC = TV-TFRZ - PPF = 4.6*APAR + PPF = 4.6*APAR_scale J = PPF*parameters%QE25 KC = parameters%KC25 * F1(parameters%AKC,TC) KO = parameters%KO25 * F1(parameters%AKO,TC) @@ -5135,7 +5138,7 @@ END SUBROUTINE STOMATA !== begin canres =================================================================================== - SUBROUTINE CANRES (parameters,PAR ,SFCTMP,RCSOIL ,EAH ,SFCPRS , & !in + SUBROUTINE CANRES (parameters,PAR ,SFCTMP,RCSOIL ,EAH ,SFCPRS , FVEG, & !in RC ,PSN ,ILOC ,JLOC ) !out ! -------------------------------------------------------------------------------------------------- @@ -5162,7 +5165,7 @@ SUBROUTINE CANRES (parameters,PAR ,SFCTMP,RCSOIL ,EAH ,SFCPRS , & !in REAL, INTENT(IN) :: SFCPRS !surface pressure (pa) REAL, INTENT(IN) :: EAH !water vapor pressure (pa) REAL, INTENT(IN) :: RCSOIL !soil moisture stress factor - + REAL, INTENT(IN) :: FVEG !vegetation greeness fraction !outputs REAL, INTENT(OUT) :: RC !canopy resistance per unit LAI @@ -5174,9 +5177,10 @@ SUBROUTINE CANRES (parameters,PAR ,SFCTMP,RCSOIL ,EAH ,SFCPRS , & !in REAL :: RCS REAL :: RCT REAL :: FF - REAL :: Q2 !water vapor mixing ratio (kg/kg) - REAL :: Q2SAT !saturation Q2 - REAL :: DQSDT2 !d(Q2SAT)/d(T) + REAL :: Q2 !water vapor mixing ratio (kg/kg) + REAL :: Q2SAT !saturation Q2 + REAL :: DQSDT2 !d(Q2SAT)/d(T) + REAL :: PAR_scale !PAR scaled by FVEG from grid level to stand scale ! RSMIN, RSMAX, TOPT, RGL, HS are canopy stress parameters set in REDPRM ! ---------------------------------------------------------------------- @@ -5186,6 +5190,7 @@ SUBROUTINE CANRES (parameters,PAR ,SFCTMP,RCSOIL ,EAH ,SFCPRS , & !in RCS = 0.0 RCT = 0.0 RCQ = 0.0 + PAR_scale = PAR / max(FVEG,1.0e-6) ! scaling PAR back to stand scale ! compute Q2 and Q2SAT @@ -5196,7 +5201,7 @@ SUBROUTINE CANRES (parameters,PAR ,SFCTMP,RCSOIL ,EAH ,SFCPRS , & !in ! contribution due to incoming solar radiation - FF = 2.0 * PAR / parameters%RGL + FF = 2.0 * PAR_scale / parameters%RGL RCS = (FF + parameters%RSMIN / parameters%RSMAX) / (1.0+ FF) RCS = MAX (RCS,0.0001) @@ -6315,7 +6320,7 @@ SUBROUTINE CANWATER (parameters,VEGTYP ,DT , & !in ! --------------------------- liquid water ------------------------------ ! maximum canopy water - MAXLIQ = parameters%CH2OP * (ELAI+ ESAI) + MAXLIQ = FVEG * parameters%CH2OP * (ELAI+ ESAI) ! evaporation, transpiration, and dew @@ -6343,7 +6348,7 @@ SUBROUTINE CANWATER (parameters,VEGTYP ,DT , & !in ! --------------------------- canopy ice ------------------------------ ! for canopy ice - MAXSNO = 6.6*(0.27+46.0/BDFALL) * (ELAI+ ESAI) + MAXSNO = FVEG * 6.6*(0.27+46.0/BDFALL) * (ELAI+ ESAI) QSUBC = MIN(CANICE/DT,QSUBC) CANICE= MAX(0.0,CANICE + (QFROC-QSUBC)*DT) @@ -6654,7 +6659,7 @@ SUBROUTINE COMBINE (parameters,NSNOW ,NSOIL ,ILOC ,JLOC , & !in SNICE(J+1) = SNICE(J+1) + SNICE(J) DZSNSO(J+1) = DZSNSO(J+1) + DZSNSO(J) ELSE - IF (ISNOW_OLD < -1) THEN ! MB/KM: change to ISNOW + IF (ISNOW < -1) THEN ! MB/KM: change to ISNOW SNLIQ(J-1) = SNLIQ(J-1) + SNLIQ(J) SNICE(J-1) = SNICE(J-1) + SNICE(J) DZSNSO(J-1) = DZSNSO(J-1) + DZSNSO(J) @@ -12183,8 +12188,10 @@ subroutine read_mp_crop_parameters() end subroutine read_mp_crop_parameters - subroutine read_mp_irrigation_parameters() + subroutine read_mp_irrigation_parameters(iopt_irr) implicit none + + INTEGER, INTENT(in) :: iopt_irr integer :: ierr logical :: file_named @@ -12211,12 +12218,18 @@ subroutine read_mp_irrigation_parameters() FIRTFAC_TABLE = -1.0E36 ! flood application rate factor IR_RAIN_TABLE = -1.0E36 ! maximum precipitation to stop irrigation trigger - inquire( file='MPTABLE.TBL', exist=file_named ) - if ( file_named ) then - open(15, file="MPTABLE.TBL", status='old', form='formatted', action='read', iostat=ierr) - else - open(15, status='old', form='formatted', action='read', iostat=ierr) - end if + IF (iopt_irr >= 1) THEN + inquire( file='MPTABLE.TBL', exist=file_named ) + if ( file_named ) then + open(15, file="MPTABLE.TBL", status='old', form='formatted', action='read', iostat=ierr) + else + open(15, status='old', form='formatted', action='read', iostat=ierr) + end if + ELSE + ! Nothing to read, nothing to do + RETURN + + END IF if (ierr /= 0) then write(*,'("WARNING: Cannot find file MPTABLE.TBL")')