diff --git a/doc/ChangeLog b/doc/ChangeLog index b8ba86ab3d..ecd9de8758 100644 --- a/doc/ChangeLog +++ b/doc/ChangeLog @@ -1,4 +1,89 @@ =============================================================== +Tag name: ctsm5.1.dev042 +Originator(s): erik (Erik Kluzek,UCAR/TSS,303-497-1326) +Date: Fri May 28 14:07:45 MDT 2021 +One-line Summary: Small answer changes for double precision constants and limit on organic soil + +Purpose and description of changes +---------------------------------- + +Change more constants to double precision. Add in limits on organic matter in +soil that @olyson added to the PPE branch (Perturbed Parameter Ensemble). + + +Significant changes to scientifically-supported configurations +-------------------------------------------------------------- + +Does this tag change answers significantly for any of the following physics configurations? +(Details of any changes will be given in the "Answer changes" section below.) + + [Put an [X] in the box for any configuration with significant answer changes.] + +[ ] clm5_1 + +[ ] clm5_0 + +[ ] ctsm5_0-nwp + +[ ] clm4_5 + + +Bugs fixed or introduced +------------------------ + +Issues fixed (include CTSM Issue #): + #1380 -- Missing NCL script + #142 --- hard coded constants aren't all double precision (some work was done on this) + +Notes of particular relevance for developers: +--------------------------------------------- +NOTE: Be sure to review the steps in README.CHECKLIST.master_tags as well as the coding style in the Developers Guide + +Caveats for developers (e.g., code that is duplicated that requires double maintenance): + + There are still more constants that should be converted to double precision. + +Testing summary: regular +---------------- + [PASS means all tests PASS; OK means tests PASS other than expected fails.] + + build-namelist tests (if CLMBuildNamelist.pm has changed): + + cheyenne - PASS + + python testing (if python code has changed; see instructions in python/README.md; document testing done): + + cheyenne - PASS + + regular tests (aux_clm: https://github.com/ESCOMP/CTSM/wiki/System-Testing-Guide#pre-merge-system-testing): + + cheyenne ---- OK + izumi ------- OK + +If the tag used for baseline comparisons was NOT the previous tag, note that here: + + +Answer changes +-------------- + +Changes answers relative to baseline: Yes! + + Summarize any changes to answers, i.e., + - what code configurations: All + - what platforms/compilers: All + - nature of change: mostly single-precision roundoff + many constants that were single precision are now double + some limits put on soil organic that weren't there previously + +Other details +------------- + +Pull Requests that document the changes (include PR ids): +(https://github.com/ESCOMP/ctsm/pull) + #1384 -- Change some constants to double precision, and add some soil limits + +=============================================================== +=============================================================== Tag name: ctsm5.1.dev041 Originator(s): Ryan Knox Date: Thu May 27 01:10:40 MDT 2021 diff --git a/doc/ChangeSum b/doc/ChangeSum index d966c4135b..f71295bc96 100644 --- a/doc/ChangeSum +++ b/doc/ChangeSum @@ -1,5 +1,6 @@ Tag Who Date Summary ============================================================================================================================ + ctsm5.1.dev042 erik 05/28/2021 Small answer changes for double precision constants and soil limits ctsm5.1.dev041 rgknox 05/27/2021 bring FATES API up to sci.1.46.0_api.16.0.0 (methane and cn hooks) ctsm5.1.dev040 slevis 05/20/2021 mksurfdata_map: replace SRC files of various masks with SRC files w no mask ctsm5.1.dev039 jedwards 05/18/2021 Add NEON sites diff --git a/src/biogeochem/CNC14DecayMod.F90 b/src/biogeochem/CNC14DecayMod.F90 index 435ec27519..2a4cbfb204 100644 --- a/src/biogeochem/CNC14DecayMod.F90 +++ b/src/biogeochem/CNC14DecayMod.F90 @@ -107,7 +107,7 @@ subroutine C14Decay( bounds, num_soilc, filter_soilc, num_soilp, filter_soilp, & spinup_term = spinup_term * get_spinup_latitude_term(grc%latdeg(col%gridcell(c))) endif else - spinup_term = 1. + spinup_term = 1._r8 endif decomp_cpools_vr(c,j,l) = decomp_cpools_vr(c,j,l) * (1._r8 - decay_const * spinup_term * dt) end do diff --git a/src/biogeochem/CNCStateUpdate1Mod.F90 b/src/biogeochem/CNCStateUpdate1Mod.F90 index c5a717f2bc..5733509539 100644 --- a/src/biogeochem/CNCStateUpdate1Mod.F90 +++ b/src/biogeochem/CNCStateUpdate1Mod.F90 @@ -166,7 +166,7 @@ subroutine CStateUpdate1( num_soilc, filter_soilc, num_soilp, filter_soilp, & real(r8) :: dt ! radiation time step (seconds) real(r8) :: check_cpool real(r8) :: cpool_delta - real(r8), parameter :: kprod05 = 1.44e-7 ! decay constant for 0.5-year product pool (1/s) (lose ~90% over a half year) + real(r8), parameter :: kprod05 = 1.44e-7_r8 ! decay constant for 0.5-year product pool (1/s) (lose ~90% over a half year) !----------------------------------------------------------------------- associate( & diff --git a/src/biogeochem/CNDVEstablishmentMod.F90 b/src/biogeochem/CNDVEstablishmentMod.F90 index 985af76a12..ffd213e3a3 100644 --- a/src/biogeochem/CNDVEstablishmentMod.F90 +++ b/src/biogeochem/CNDVEstablishmentMod.F90 @@ -82,7 +82,7 @@ subroutine Establishment(bounds, & real(r8):: bm_delta ! parameters - real(r8), parameter :: ramp_agddtw = 300.0 + real(r8), parameter :: ramp_agddtw = 300.0_r8 ! minimum individual density for persistence of PATCH (indiv/m2) real(r8), parameter :: nind_min = 1.0e-10_r8 @@ -425,8 +425,8 @@ subroutine Establishment(bounds, & greffic(p) = bm_delta / (lm_ind * slatop(ivt(p))) end if else - greffic(p) = 0. - heatstress(p) = 0. + greffic(p) = 0._r8 + heatstress(p) = 0._r8 end if end do diff --git a/src/biogeochem/CNFUNMod.F90 b/src/biogeochem/CNFUNMod.F90 index 6ab724aae2..57e8e11c86 100644 --- a/src/biogeochem/CNFUNMod.F90 +++ b/src/biogeochem/CNFUNMod.F90 @@ -1233,13 +1233,13 @@ subroutine CNFUN(bounds,num_soilc, filter_soilc,num_soilp& ! C used for uptake is reduced if the cost of N is very high frac_ideal_C_use = max(0.0_r8,1.0_r8 - (total_N_resistance-fun_cn_flex_a(ivt(p)))/fun_cn_flex_b(ivt(p)) ) ! then, if the plant is very much in need of N, the C used for uptake is increased accordingly. - if(delta_CN.lt.0.0)then + if(delta_CN.lt.0.0_r8)then frac_ideal_C_use = frac_ideal_C_use + (1.0_r8-frac_ideal_C_use)*min(1.0_r8, delta_CN/fun_cn_flex_c(ivt(p))) end if ! If we have too much N (e.g. from free N retranslocation) then make frac_ideal_c_use even lower. ! For a CN delta of fun_cn_flex_c, then we reduce C expendiure to the minimum of 0.5. ! This seems a little intense? - if(delta_CN .gt.0.and. frac_ideal_C_use.lt.1.0)then + if(delta_CN .gt.0._r8 .and. frac_ideal_C_use.lt.1.0_r8)then frac_ideal_C_use = frac_ideal_C_use + 0.5_r8*(1.0_r8*delta_CN/fun_cn_flex_c(ivt(p))) end if ! don't let this go above 1 or below an arbitrary minimum (to prevent zero N uptake). diff --git a/src/biogeochem/CNFireBaseMod.F90 b/src/biogeochem/CNFireBaseMod.F90 index b9c125716e..4c83c5a9a5 100644 --- a/src/biogeochem/CNFireBaseMod.F90 +++ b/src/biogeochem/CNFireBaseMod.F90 @@ -1042,7 +1042,7 @@ subroutine CNFireFluxes (this, bounds, num_soilc, filter_soilc, num_soilp, filte if( trotr1_col(c)+trotr2_col(c) > 0.6_r8 .and. dtrotr_col(c) > 0._r8 .and. & lfc(c) > 0._r8 .and. fbac1(c) == 0._r8) then lfc2(c) = max(0._r8, min(lfc(c), (farea_burned(c)-baf_crop(c) - & - baf_peatf(c))/2.0*dt))/(dtrotr_col(c)*dayspyr*secspday/dt)/dt + baf_peatf(c))/2.0_r8*dt))/(dtrotr_col(c)*dayspyr*secspday/dt)/dt lfc(c) = lfc(c) - max(0._r8, min(lfc(c), (farea_burned(c)-baf_crop(c) - & baf_peatf(c))*dt/2.0_r8)) end if diff --git a/src/biogeochem/CNFireLi2014Mod.F90 b/src/biogeochem/CNFireLi2014Mod.F90 index e8fd78230e..a520ba36e1 100644 --- a/src/biogeochem/CNFireLi2014Mod.F90 +++ b/src/biogeochem/CNFireLi2014Mod.F90 @@ -559,7 +559,7 @@ subroutine CNFireArea (this, bounds, num_soilc, filter_soilc, num_soilp, filter_ g = col%gridcell(c) hdmlf=this%forc_hdm(g) nfire(c) = 0._r8 - if( cropf_col(c) < 1.0 )then + if( cropf_col(c) < 1.0_r8 )then if (trotr1_col(c)+trotr2_col(c)>0.6_r8) then farea_burned(c)=min(1.0_r8,baf_crop(c)+baf_peatf(c)) else @@ -1228,7 +1228,7 @@ subroutine CNFireFluxes (this, bounds, num_soilc, filter_soilc, num_soilp, filte if( trotr1_col(c)+trotr2_col(c) > 0.6_r8 .and. dtrotr_col(c) > 0._r8 .and. & lfc(c) > 0._r8 .and. fbac1(c) == 0._r8) then lfc2(c) = max(0._r8, min(lfc(c), (farea_burned(c)-baf_crop(c) - & - baf_peatf(c))/2.0*dt))/(dtrotr_col(c)*dayspyr*secspday/dt)/dt + baf_peatf(c))/2.0_r8*dt))/(dtrotr_col(c)*dayspyr*secspday/dt)/dt lfc(c) = lfc(c) - max(0._r8, min(lfc(c), (farea_burned(c)-baf_crop(c) - & baf_peatf(c))*dt/2.0_r8)) end if diff --git a/src/biogeochem/CNGapMortalityMod.F90 b/src/biogeochem/CNGapMortalityMod.F90 index cd02221de4..abd64cecc9 100644 --- a/src/biogeochem/CNGapMortalityMod.F90 +++ b/src/biogeochem/CNGapMortalityMod.F90 @@ -112,12 +112,12 @@ subroutine CNGapMortality (bounds, num_soilc, filter_soilc, num_soilp, filter_so real(r8) , intent(in) :: stem_prof_patch(bounds%begp:,1:) ! ! !LOCAL VARIABLES: - integer :: p ! patch index - integer :: fp ! patch filter index - real(r8):: am ! rate for fractional mortality (1/yr) - real(r8):: m ! rate for fractional mortality (1/s) - real(r8):: mort_max ! asymptotic max mortality rate (/yr) - real(r8):: k_mort = 0.3 ! coeff of growth efficiency in mortality equation + integer :: p ! patch index + integer :: fp ! patch filter index + real(r8):: am ! rate for fractional mortality (1/yr) + real(r8):: m ! rate for fractional mortality (1/s) + real(r8):: mort_max ! asymptotic max mortality rate (/yr) + real(r8):: k_mort = 0.3_r8 ! coeff of growth efficiency in mortality equation !----------------------------------------------------------------------- SHR_ASSERT_ALL_FL((ubound(leaf_prof_patch) == (/bounds%endp,nlevdecomp_full/)), sourcefile, __LINE__) diff --git a/src/biogeochem/CNProductsMod.F90 b/src/biogeochem/CNProductsMod.F90 index 0586815ee1..617672d78c 100644 --- a/src/biogeochem/CNProductsMod.F90 +++ b/src/biogeochem/CNProductsMod.F90 @@ -493,9 +493,9 @@ subroutine UpdateProducts(this, bounds, & ! calculate losses from product pools ! the following (1/s) rate constants result in ~90% loss of initial state over 1, 10 and 100 years, ! respectively, using a discrete-time fractional decay algorithm. - kprod1 = 7.2e-8 - kprod10 = 7.2e-9 - kprod100 = 7.2e-10 + kprod1 = 7.2e-8_r8 + kprod10 = 7.2e-9_r8 + kprod100 = 7.2e-10_r8 do g = bounds%begg, bounds%endg ! calculate fluxes out of product pools (1/sec) diff --git a/src/biogeochem/DryDepVelocity.F90 b/src/biogeochem/DryDepVelocity.F90 index 37860e9728..01aa51d52c 100644 --- a/src/biogeochem/DryDepVelocity.F90 +++ b/src/biogeochem/DryDepVelocity.F90 @@ -400,7 +400,7 @@ subroutine depvel_compute( bounds, & endif if (index_season<0) then - if (elai(pi) < (minlai+0.05*(maxlai-minlai))) then + if (elai(pi) < (minlai+0.05_r8*(maxlai-minlai))) then index_season = 3 endif endif diff --git a/src/biogeochem/NutrientCompetitionCLM45defaultMod.F90 b/src/biogeochem/NutrientCompetitionCLM45defaultMod.F90 index f81783ad8c..d93f0a168e 100644 --- a/src/biogeochem/NutrientCompetitionCLM45defaultMod.F90 +++ b/src/biogeochem/NutrientCompetitionCLM45defaultMod.F90 @@ -272,7 +272,7 @@ subroutine calc_plant_cn_alloc (this, bounds, num_soilp, filter_soilp, & ! allocation as specified in the pft-physiology file. The value is also used ! as a trigger here: -1.0 means to use the dynamic allocation (trees). if (stem_leaf(ivt(p)) == -1._r8) then - f3 = (2.7/(1.0+exp(-0.004*(annsum_npp(p) - 300.0)))) - 0.4 + f3 = (2.7_r8/(1.0_r8+exp(-0.004_r8*(annsum_npp(p) - 300.0_r8)))) - 0.4_r8 else f3 = stem_leaf(ivt(p)) end if @@ -753,7 +753,7 @@ subroutine calc_plant_nitrogen_demand(this, bounds, num_soilp, filter_soilp, & ! as a trigger here: -1.0 means to use the dynamic allocation (trees). if (stem_leaf(ivt(p)) == -1._r8) then - f3 = (2.7/(1.0+exp(-0.004*(annsum_npp(p) - 300.0)))) - 0.4 + f3 = (2.7_r8/(1.0_r8+exp(-0.004_r8*(annsum_npp(p) - 300.0_r8)))) - 0.4_r8 else f3 = stem_leaf(ivt(p)) end if diff --git a/src/biogeochem/NutrientCompetitionFlexibleCNMod.F90 b/src/biogeochem/NutrientCompetitionFlexibleCNMod.F90 index b69c666ea4..a642a5dd2d 100644 --- a/src/biogeochem/NutrientCompetitionFlexibleCNMod.F90 +++ b/src/biogeochem/NutrientCompetitionFlexibleCNMod.F90 @@ -421,7 +421,7 @@ subroutine calc_plant_cn_alloc(this, bounds, num_soilp, filter_soilp, & ! allocation as specified in the pft-physiology file. The value is also used ! as a trigger here: -1.0 means to use the dynamic allocation (trees). if (stem_leaf(ivt(p)) == -1._r8) then - f3 = (2.7/(1.0+exp(-0.004*(annsum_npp(p) - 300.0)))) - 0.4 + f3 = (2.7_r8/(1.0_r8+exp(-0.004_r8*(annsum_npp(p) - 300.0_r8)))) - 0.4_r8 else f3 = stem_leaf(ivt(p)) end if @@ -1468,7 +1468,7 @@ subroutine calc_plant_nitrogen_demand(this, bounds, num_soilp, filter_soilp, & ! as a trigger here: -1.0 means to use the dynamic allocation (trees). if (stem_leaf(ivt(p)) == -1._r8) then - f3 = (2.7/(1.0+exp(-0.004*(annsum_npp(p) - 300.0)))) - 0.4 + f3 = (2.7_r8/(1.0_r8+exp(-0.004_r8*(annsum_npp(p) - 300.0_r8)))) - 0.4_r8 else f3 = stem_leaf(ivt(p)) end if diff --git a/src/biogeochem/VOCEmissionMod.F90 b/src/biogeochem/VOCEmissionMod.F90 index b50fe5ec77..88ae7c08cc 100644 --- a/src/biogeochem/VOCEmissionMod.F90 +++ b/src/biogeochem/VOCEmissionMod.F90 @@ -601,7 +601,7 @@ subroutine VOCEmission (bounds, num_soilp, filter_soilp, & ! Activity factor for CO2 (only for isoprene) if (trim(meg_cmp%name) == 'isoprene') then - co2_ppmv = 1.e6*forc_pco2(g)/forc_pbot(c) + co2_ppmv = 1.e6_r8*forc_pco2(g)/forc_pbot(c) gamma_c = get_gamma_C(cisun_z(p,1),cisha_z(p,1),forc_pbot(c),fsun(p), co2_ppmv) else gamma_c = 1._r8 diff --git a/src/biogeophys/ActiveLayerMod.F90 b/src/biogeophys/ActiveLayerMod.F90 index ca980fe15b..72ce965f5e 100644 --- a/src/biogeophys/ActiveLayerMod.F90 +++ b/src/biogeophys/ActiveLayerMod.F90 @@ -116,10 +116,10 @@ subroutine alt_calc(this, num_soilc, filter_soilc, & do fc = 1,num_soilc c = filter_soilc(fc) g = col%gridcell(c) - if ( grc%lat(g) > 0. ) then + if ( grc%lat(g) > 0._r8 ) then altmax_lastyear(c) = altmax(c) altmax_lastyear_indx(c) = altmax_indx(c) - altmax(c) = 0. + altmax(c) = 0._r8 altmax_indx(c) = 0 endif end do @@ -128,10 +128,10 @@ subroutine alt_calc(this, num_soilc, filter_soilc, & do fc = 1,num_soilc c = filter_soilc(fc) g = col%gridcell(c) - if ( grc%lat(g) <= 0. ) then + if ( grc%lat(g) <= 0._r8 ) then altmax_lastyear(c) = altmax(c) altmax_lastyear_indx(c) = altmax_indx(c) - altmax(c) = 0. + altmax(c) = 0._r8 altmax_indx(c) = 0 endif end do diff --git a/src/biogeophys/BalanceCheckMod.F90 b/src/biogeophys/BalanceCheckMod.F90 index 7df5b10327..8c6378d731 100644 --- a/src/biogeophys/BalanceCheckMod.F90 +++ b/src/biogeophys/BalanceCheckMod.F90 @@ -604,8 +604,8 @@ subroutine BalanceCheck( bounds, & l = col%landunit(c) if (col%itype(c) == icol_sunwall .or. col%itype(c) == icol_shadewall) then - forc_rain_col(c) = 0. - forc_snow_col(c) = 0. + forc_rain_col(c) = 0._r8 + forc_snow_col(c) = 0._r8 else forc_rain_col(c) = forc_rain(c) forc_snow_col(c) = forc_snow(c) diff --git a/src/biogeophys/BandDiagonalMod.F90 b/src/biogeophys/BandDiagonalMod.F90 index 7cd9cf204f..13529a21ef 100644 --- a/src/biogeophys/BandDiagonalMod.F90 +++ b/src/biogeophys/BandDiagonalMod.F90 @@ -176,7 +176,7 @@ subroutine BandDiagonal(bounds, lbj, ubj, jtop, jbot, numf, filter, nband, b, r, n=jbot(ci)-jtop(ci)+1 allocate(ab(m,n)) - ab=0.0 + ab=0.0_r8 ab(kl+ku-1,3:n)=b(ci,1,jtop(ci):jbot(ci)-2) ! 2nd superdiagonal ab(kl+ku+0,2:n)=b(ci,2,jtop(ci):jbot(ci)-1) ! 1st superdiagonal diff --git a/src/biogeophys/CanopyFluxesMod.F90 b/src/biogeophys/CanopyFluxesMod.F90 index f18bde2e0e..34d05f84a9 100644 --- a/src/biogeophys/CanopyFluxesMod.F90 +++ b/src/biogeophys/CanopyFluxesMod.F90 @@ -712,7 +712,7 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, ! leaf and stem surface area sa_leaf(p) = elai(p) ! double in spirit of full surface area for sensible heat - sa_leaf(p) = 2.*sa_leaf(p) + sa_leaf(p) = 2._r8*sa_leaf(p) ! Surface area for stem sa_stem(p) = nstem(patch%itype(p))*(htop(p)*shr_const_pi*dbh(p)) @@ -725,8 +725,8 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, ! (set surface area for stem, and fraction absorbed by stem to zero) if(.not.(is_tree(patch%itype(p)) .or. is_shrub(patch%itype(p))) & .or. dbh(p) < min_stem_diameter) then - frac_rad_abs_by_stem(p) = 0.0 - sa_stem(p) = 0.0 + frac_rad_abs_by_stem(p) = 0.0_r8 + sa_stem(p) = 0.0_r8 endif ! if using Satellite Phenology mode, calculate leaf and stem biomass @@ -734,12 +734,12 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, ! 2gbiomass/gC * (1/SLA) * 1e-3 = kg dry mass/m2(leaf) leaf_biomass(p) = (1.e-3_r8*c_to_b/slatop(patch%itype(p))) & * max(0.01_r8, 0.5_r8*sa_leaf(p)) & - / (1.-fbw(patch%itype(p))) + / (1._r8-fbw(patch%itype(p))) ! cross-sectional area of stems - carea_stem = shr_const_pi * (dbh(p)*0.5)**2 + carea_stem = shr_const_pi * (dbh(p)*0.5_r8)**2 stem_biomass(p) = carea_stem * htop(p) * k_cyl_vol & * nstem(patch%itype(p)) * wood_density(patch%itype(p)) & - /(1.-fbw(patch%itype(p))) + /(1._r8-fbw(patch%itype(p))) endif ! internal longwave fluxes between leaf and stem @@ -752,10 +752,10 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, ! lma_dry has units of kg dry mass/m2 here ! (Appendix B of Bonan et al., GMD, 2018) - cp_leaf(p) = leaf_biomass(p) * (c_dry_biomass*(1.-fbw(patch%itype(p))) + (fbw(patch%itype(p)))*c_water) + cp_leaf(p) = leaf_biomass(p) * (c_dry_biomass*(1._r8-fbw(patch%itype(p))) + (fbw(patch%itype(p)))*c_water) ! cp-stem will have units J/k/ground_area - cp_stem(p) = stem_biomass(p) * (c_dry_biomass*(1.-fbw(patch%itype(p))) + (fbw(patch%itype(p)))*c_water) + cp_stem(p) = stem_biomass(p) * (c_dry_biomass*(1._r8-fbw(patch%itype(p))) + (fbw(patch%itype(p)))*c_water) ! adjust for departure from cylindrical stem model cp_stem(p) = k_cyl_vol * cp_stem(p) @@ -1395,7 +1395,7 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, dt_stem(p) = (frac_rad_abs_by_stem(p)*(sabv(p) + air(p) + bir(p)*ts_ini(p)**4 & + cir(p)*lw_grnd) - eflx_sh_stem(p) & + lw_leaf(p)- lw_stem(p))/(cp_stem(p)/dtime & - - frac_rad_abs_by_stem(p)*bir(p)*4.*ts_ini(p)**3) + - frac_rad_abs_by_stem(p)*bir(p)*4._r8*ts_ini(p)**3) else dt_stem(p) = 0._r8 endif diff --git a/src/biogeophys/LakeFluxesMod.F90 b/src/biogeophys/LakeFluxesMod.F90 index 212d0ca7d1..4528d9dfb2 100644 --- a/src/biogeophys/LakeFluxesMod.F90 +++ b/src/biogeophys/LakeFluxesMod.F90 @@ -182,8 +182,8 @@ subroutine LakeFluxes(bounds, num_lakec, filter_lakec, num_lakep, filter_lakep, real(r8) :: kva0temp ! (K) temperature for kva0; will be set below real(r8), parameter :: kva0pres = 1.013e5_r8 ! (Pa) pressure for kva0 real(r8) :: kva ! kinematic viscosity of air at ground temperature and forcing pressure - real(r8), parameter :: prn = 0.713 ! Prandtl # for air at neutral stability - real(r8), parameter :: sch = 0.66 ! Schmidt # for water in air at neutral stability + real(r8), parameter :: prn = 0.713_r8 ! Prandtl # for air at neutral stability + real(r8), parameter :: sch = 0.66_r8 ! Schmidt # for water in air at neutral stability !----------------------------------------------------------------------- diff --git a/src/biogeophys/LunaMod.F90 b/src/biogeophys/LunaMod.F90 index bd4fdfb4b0..604b2b5709 100644 --- a/src/biogeophys/LunaMod.F90 +++ b/src/biogeophys/LunaMod.F90 @@ -360,11 +360,11 @@ subroutine Update_Photosynthesis_Capacity(bounds, fn, filterp, & !Implemented the nitrogen allocation model if(tlai(p) > 0.0_r8 .and. lnc(p) > 0._r8)then RadTop = par240d_z(p,1)/rabsorb - PARTop = RadTop*4.6 !conversion from w/m2 to umol/m2/s. PAR is still in umol photons, not electrons. Also the par240d_z is only for radiation at visible range. Hence 4.6 not 2.3 multiplier. + PARTop = RadTop*4.6_r8 !conversion from w/m2 to umol/m2/s. PAR is still in umol photons, not electrons. Also the par240d_z is only for radiation at visible range. Hence 4.6 not 2.3 multiplier. !------------------------------------------------------------- !the nitrogen allocation model, may need to be feed from the parameter file in CLM if (nint(c3psn(ft)) == 1)then - if(gpp_day(p)>0.0 )then !only optimize if there is growth and it is C3 plants + if(gpp_day(p)>0.0_r8 )then !only optimize if there is growth and it is C3 plants !------------------------------------------------------------- do z = 1, nrad(p) if(tlai_z(p,z)>0.0_r8)then @@ -434,7 +434,7 @@ subroutine Update_Photosynthesis_Capacity(bounds, fn, filterp, & PNlc_z(p, z)= PNlcopt - if(enzs_z(p,z)<1.0) then + if(enzs_z(p,z)<1.0_r8) then enzs_z(p,z) = enzs_z(p,z)* (1.0_r8 + max_daily_pchg) endif !nitrogen allocastion model-end @@ -581,7 +581,7 @@ subroutine Acc240_Climate_LUNA(bounds, fn, filterp, oair, cair, & if(t_veg_day(p).ne.spval) then !check whether it is the first day !--------------------------------------------------------- !calculate the 10 day running mean radiations - if(ndaysteps(p)>0.0) then + if(ndaysteps(p)>0.0_r8) then par24d_z_i=par24d_z(p,:)/(dtime * ndaysteps(p)) else par24d_z_i = 0._r8 @@ -595,7 +595,7 @@ subroutine Acc240_Climate_LUNA(bounds, fn, filterp, oair, cair, & endif !------------------------------------------------------- !calculate the 10 day running mean daytime temperature - if(ndaysteps(p)>0.0)then + if(ndaysteps(p)>0.0_r8)then t_veg_dayi = t_veg_day(p) / ndaysteps(p) else t_veg_dayi = t_veg_night(p) / nnightsteps(p) @@ -898,7 +898,7 @@ subroutine NitrogenAllocation(FNCa,forc_pbot10, relh10, CO2a10,O2a10, PARi10,PAR Nresp = PNrespold * FNCa !proportion of respirational nitrogen in functional nitrogen Ncb = PNcbold * FNCa !proportion of carboxylation nitrogen in functional nitrogen if (Nlc > FNCa * 0.5_r8) Nlc = 0.5_r8 * FNCa - chg_per_step = 0.02* FNCa + chg_per_step = 0.02_r8* FNCa PNlc = PNlcold PNlcoldi = PNlcold - 0.001_r8 PARi10c = max(PARLowLim, PARi10) diff --git a/src/biogeophys/PhotosynthesisMod.F90 b/src/biogeophys/PhotosynthesisMod.F90 index ae82794545..18aabe59b1 100644 --- a/src/biogeophys/PhotosynthesisMod.F90 +++ b/src/biogeophys/PhotosynthesisMod.F90 @@ -2851,7 +2851,7 @@ subroutine PhotosynthesisHydraulicStress ( bounds, fn, filterp, & r_soil = sqrt(1./(rpi*root_length_density)) ! length scale approach - soil_conductance = min(hksat(c,j),hk_l(c,j))/(1.e3*r_soil) + soil_conductance = min(hksat(c,j),hk_l(c,j))/(1.e3_r8*r_soil) ! use vegetation plc function to adjust root conductance fs(j)= plc(smp(c,j),p,c,root,veg) @@ -2873,7 +2873,7 @@ subroutine PhotosynthesisHydraulicStress ( bounds, fn, filterp, & if(rai(j)*rootfr(p,j) > 0._r8 .and. j > 1) then k_soil_root(p,j) = 1._r8/rs_resis else - k_soil_root(p,j) = 0. + k_soil_root(p,j) = 0._r8 endif end do @@ -4163,8 +4163,8 @@ subroutine ci_func_PHS(x,cisun, cisha, fvalsun, fvalsha, p, iv, c, bsun, bsha, b bquad = -(2.0 * (medlynintercept(patch%itype(p))*1.e-06_r8 + term) + (medlynslope(patch%itype(p)) * term)**2 / & (gb_mol*1.e-06_r8 * rh_can)) cquad = medlynintercept(patch%itype(p))*medlynintercept(patch%itype(p))*1.e-12_r8 + & - (2.0*medlynintercept(patch%itype(p))*1.e-06_r8 + term * & - (1.0 - medlynslope(patch%itype(p))* medlynslope(patch%itype(p)) / rh_can)) * term + (2.0_r8*medlynintercept(patch%itype(p))*1.e-06_r8 + term * & + (1.0_r8 - medlynslope(patch%itype(p))* medlynslope(patch%itype(p)) / rh_can)) * term call quadratic (aquad, bquad, cquad, r1, r2) gs_mol_sun = max(r1,r2) * 1.e06_r8 @@ -4278,7 +4278,7 @@ subroutine calcstress(p,c,x,bsun,bsha,gb_mol,gs_mol_sun,gs_mol_sha,qsatl,qaf, & logical :: flag ! signal that matrix was not invertible logical :: night ! signal to store vegwp within this routine, b/c it is night-time and full suite won't be called integer, parameter :: itmax=50 ! exit newton's method if iters>itmax - real(r8), parameter :: tolf=1.e-6,toldx=1.e-9 !tolerances for a satisfactory solution + real(r8), parameter :: tolf=1.e-6_r8,toldx=1.e-9_r8 !tolerances for a satisfactory solution logical :: havegs ! signals direction of calculation gs->qflx or qflx->gs real(r8) :: soilflux ! total soil column transpiration [mm/s] real(r8), parameter :: tol_lai=.001_r8 ! minimum lai where transpiration is calc'd diff --git a/src/biogeophys/RootBiophysMod.F90 b/src/biogeophys/RootBiophysMod.F90 index 6e94ddef4d..c6130c306b 100644 --- a/src/biogeophys/RootBiophysMod.F90 +++ b/src/biogeophys/RootBiophysMod.F90 @@ -269,7 +269,7 @@ function jackson1996_rootfr(bounds, ubj, varindx, water_carbon) result(rootfr) beta ** (col%zi(c,lev)*m_to_cm) ) end do else - rootfr(p,:) = 0. + rootfr(p,:) = 0._r8 endif enddo diff --git a/src/biogeophys/SnowHydrologyMod.F90 b/src/biogeophys/SnowHydrologyMod.F90 index 7127c899a8..6e5c39f8a6 100644 --- a/src/biogeophys/SnowHydrologyMod.F90 +++ b/src/biogeophys/SnowHydrologyMod.F90 @@ -3279,7 +3279,7 @@ subroutine BulkFlux_SnowCappingFluxes(bounds, num_snowc, filter_snowc, & ! Always keep at least this fraction of the bottom snow layer when doing snow capping ! This needs to be slightly greater than 0 to avoid roundoff problems - real(r8), parameter :: min_snow_to_keep = 1.e-9 ! fraction of bottom snow layer to keep with capping + real(r8), parameter :: min_snow_to_keep = 1.e-9_r8 ! fraction of bottom snow layer to keep with capping character(len=*), parameter :: subname = 'BulkFlux_SnowCappingFluxes' !----------------------------------------------------------------------- diff --git a/src/biogeophys/SnowSnicarMod.F90 b/src/biogeophys/SnowSnicarMod.F90 index 1f8e928474..196d239460 100644 --- a/src/biogeophys/SnowSnicarMod.F90 +++ b/src/biogeophys/SnowSnicarMod.F90 @@ -680,31 +680,31 @@ subroutine SNICAR_RT (flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, & ! Eddington if (APRX_TYP==1) then do i=snl_top,snl_btm,1 - gamma1(i) = (7-(omega_star(i)*(4+(3*g_star(i)))))/4 - gamma2(i) = -(1-(omega_star(i)*(4-(3*g_star(i)))))/4 - gamma3(i) = (2-(3*g_star(i)*mu_not))/4 - gamma4(i) = 1-gamma3(i) - mu_one = 0.5 + gamma1(i) = (7._r8-(omega_star(i)*(4._r8+(3._r8*g_star(i)))))/4._r8 + gamma2(i) = -(1._r8-(omega_star(i)*(4._r8-(3._r8*g_star(i)))))/4._r8 + gamma3(i) = (2._r8-(3._r8*g_star(i)*mu_not))/4._r8 + gamma4(i) = 1._r8-gamma3(i) + mu_one = 0.5_r8 enddo ! Quadrature elseif (APRX_TYP==2) then do i=snl_top,snl_btm,1 - gamma1(i) = (3**0.5)*(2-(omega_star(i)*(1+g_star(i))))/2 - gamma2(i) = omega_star(i)*(3**0.5)*(1-g_star(i))/2 - gamma3(i) = (1-((3**0.5)*g_star(i)*mu_not))/2 - gamma4(i) = 1-gamma3(i) - mu_one = 1/(3**0.5) + gamma1(i) = (3._r8**0.5)*(2._r8-(omega_star(i)*(1._r8+g_star(i))))/2._r8 + gamma2(i) = omega_star(i)*(3._r8**0.5)*(1._r8-g_star(i))/2._r8 + gamma3(i) = (1._r8-((3._r8**0.5)*g_star(i)*mu_not))/2._r8 + gamma4(i) = 1._r8-gamma3(i) + mu_one = 1._r8/(3._r8**0.5_r8) enddo ! Hemispheric Mean elseif (APRX_TYP==3) then do i=snl_top,snl_btm,1 - gamma1(i) = 2 - (omega_star(i)*(1+g_star(i))) + gamma1(i) = 2._r8 - (omega_star(i)*(1._r8+g_star(i))) gamma2(i) = omega_star(i)*(1-g_star(i)) - gamma3(i) = (1-((3**0.5)*g_star(i)*mu_not))/2 - gamma4(i) = 1-gamma3(i) - mu_one = 0.5 + gamma3(i) = (1._r8-((3._r8**0.5_r8)*g_star(i)*mu_not))/2._r8 + gamma4(i) = 1._r8-gamma3(i) + mu_one = 0.5_r8 enddo endif @@ -755,7 +755,7 @@ subroutine SNICAR_RT (flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, & !Boundary values for i=1 and i=2*snl_lcl, specifics for i=odd and i=even if (i==(2*snl_lcl+1)) then - A(i) = 0 + A(i) = 0._r8 B(i) = e1(snl_top) D(i) = -e2(snl_top) E(i) = flx_slri_lcl(bnd_idx)-C_mns_top(snl_top) @@ -763,7 +763,7 @@ subroutine SNICAR_RT (flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, & elseif(i==0) then A(i) = e1(snl_btm)-(albsfc_lcl(bnd_idx)*e3(snl_btm)) B(i) = e2(snl_btm)-(albsfc_lcl(bnd_idx)*e4(snl_btm)) - D(i) = 0 + D(i) = 0._r8 E(i) = F_direct_btm-C_pls_btm(snl_btm)+(albsfc_lcl(bnd_idx)*C_mns_btm(snl_btm)) elseif(mod(i,2)==-1) then ! If odd and i>=3 (n=1 for i=3) @@ -828,7 +828,7 @@ subroutine SNICAR_RT (flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, & ! ERROR check: negative absorption - if (flx_abs_lcl(i,bnd_idx) < -0.00001) then + if (flx_abs_lcl(i,bnd_idx) < -0.00001_r8) then trip = 1 endif enddo diff --git a/src/biogeophys/SoilHydrologyInitTimeConstMod.F90 b/src/biogeophys/SoilHydrologyInitTimeConstMod.F90 index 61b692ae37..384c843760 100644 --- a/src/biogeophys/SoilHydrologyInitTimeConstMod.F90 +++ b/src/biogeophys/SoilHydrologyInitTimeConstMod.F90 @@ -159,7 +159,7 @@ subroutine SoilHydrologyInitTimeConst(bounds, soilhydrology_inst, soilstate_inst if ( lev <= nlevsoi )then claycol(c,lev) = soilstate_inst%cellclay_col(c,lev) sandcol(c,lev) = soilstate_inst%cellsand_col(c,lev) - om_fraccol(c,lev) = soilstate_inst%cellorg_col(c,lev) / organic_max + om_fraccol(c,lev) = min(soilstate_inst%cellorg_col(c,lev) / organic_max, 1._r8) else claycol(c,lev) = soilstate_inst%cellclay_col(c,nlevsoi) sandcol(c,lev) = soilstate_inst%cellsand_col(c,nlevsoi) @@ -190,14 +190,14 @@ subroutine SoilHydrologyInitTimeConst(bounds, soilhydrology_inst, soilstate_inst soilhydrology_inst%h2osfc_thresh_col(c) = 0._r8 if (micro_sigma(c) > 1.e-6_r8 .and. (soilhydrology_inst%h2osfcflag /= 0)) then - d = 0.0 + d = 0.0_r8 do p = 1,4 - fd = 0.5*(1.0_r8+shr_spfn_erf(d/(micro_sigma(c)*sqrt(2.0)))) - pc - dfdd = exp(-d**2/(2.0*micro_sigma(c)**2))/(micro_sigma(c)*sqrt(2.0*shr_const_pi)) + fd = 0.5_r8*(1.0_r8+shr_spfn_erf(d/(micro_sigma(c)*sqrt(2.0_r8)))) - pc + dfdd = exp(-d**2/(2.0_r8*micro_sigma(c)**2))/(micro_sigma(c)*sqrt(2.0_r8*shr_const_pi)) d = d - fd/dfdd enddo - soilhydrology_inst%h2osfc_thresh_col(c) = 0.5*d*(1.0_r8+shr_spfn_erf(d/(micro_sigma(c)*sqrt(2.0)))) + & - micro_sigma(c)/sqrt(2.0*shr_const_pi)*exp(-d**2/(2.0*micro_sigma(c)**2)) + soilhydrology_inst%h2osfc_thresh_col(c) = 0.5_r8*d*(1.0_r8+shr_spfn_erf(d/(micro_sigma(c)*sqrt(2.0_r8)))) + & + micro_sigma(c)/sqrt(2.0_r8*shr_const_pi)*exp(-d**2/(2.0_r8*micro_sigma(c)**2)) soilhydrology_inst%h2osfc_thresh_col(c) = 1.e3_r8 * soilhydrology_inst%h2osfc_thresh_col(c) !convert to mm from meters else soilhydrology_inst%h2osfc_thresh_col(c) = 0._r8 @@ -298,7 +298,7 @@ subroutine initSoilParVIC(c, claycol, sandcol, om_fraccol, soilhydrology_inst) !calculate other parameters based on teh percentages soilhydrology_inst%porosity_col(c, i) = 0.489_r8 - 0.00126_r8*sandvic(i) soilhydrology_inst%expt_col(c, i) = 3._r8+ 2._r8*(2.91_r8 + 0.159_r8*clayvic(i)) - xksat = 0.0070556 *( 10.**(-0.884+0.0153*sandvic(i)) ) + xksat = 0.0070556_r8 *( 10.**(-0.884_r8+0.0153_r8*sandvic(i)) ) !consider organic matter, M.Huang soilhydrology_inst%expt_col(c, i) = & diff --git a/src/biogeophys/SoilStateInitTimeConstMod.F90 b/src/biogeophys/SoilStateInitTimeConstMod.F90 index ad2da3852f..8bacda1c99 100644 --- a/src/biogeophys/SoilStateInitTimeConstMod.F90 +++ b/src/biogeophys/SoilStateInitTimeConstMod.F90 @@ -430,20 +430,20 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) if (lev .eq. 1) then clay = clay3d(g,1) sand = sand3d(g,1) - om_frac = organic3d(g,1)/organic_max + om_frac = min(organic3d(g,1)/organic_max, 1._r8) else if (lev <= nlevsoi) then found = 0 ! reset value if (zsoi(lev) <= zisoifl(1)) then ! Search above the dataset's range of zisoifl depths clay = clay3d(g,1) sand = sand3d(g,1) - om_frac = organic3d(g,1)/organic_max + om_frac = min(organic3d(g,1)/organic_max, 1._r8) found = 1 else if (zsoi(lev) > zisoifl(nlevsoifl)) then ! Search below the dataset's range of zisoifl depths clay = clay3d(g,nlevsoifl) sand = sand3d(g,nlevsoifl) - om_frac = organic3d(g,nlevsoifl)/organic_max + om_frac = min(organic3d(g,nlevsoifl)/organic_max, 1._r8) found = 1 else ! For remaining model soil levels, search within dataset's @@ -453,7 +453,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) if (zsoi(lev) > zisoifl(j) .AND. zsoi(lev) <= zisoifl(j+1)) then clay = clay3d(g,j+1) sand = sand3d(g,j+1) - om_frac = organic3d(g,j+1)/organic_max + om_frac = min(organic3d(g,j+1)/organic_max, 1._r8) found = 1 endif if (found == 1) exit ! no need to stay in the loop @@ -503,8 +503,9 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) om_hksat = max(0.28_r8 - 0.2799_r8*(zsoi(lev)/zsapric), xksat) soilstate_inst%bd_col(c,lev) = (1._r8 - soilstate_inst%watsat_col(c,lev))*params_inst%pd - soilstate_inst%watsat_col(c,lev) = params_inst%watsat_adjustfactor * ( (1._r8 - om_frac) * & - soilstate_inst%watsat_col(c,lev) + om_watsat*om_frac ) + ! do not allow watsat_sf to push watsat above 0.93 + soilstate_inst%watsat_col(c,lev) = min(params_inst%watsat_adjustfactor * ( (1._r8 - om_frac) * & + soilstate_inst%watsat_col(c,lev) + om_watsat*om_frac ), 0.93_r8) tkm = (1._r8-om_frac) * (params_inst%tkd_sand*sand+params_inst%tkd_clay*clay)/ & (sand+clay)+params_inst%tkm_om*om_frac ! W/(m K) soilstate_inst%bsw_col(c,lev) = params_inst%bsw_adjustfactor * ( (1._r8-om_frac) * & @@ -586,9 +587,9 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) clay = soilstate_inst%cellclay_col(c,lev) sand = soilstate_inst%cellsand_col(c,lev) if ( organic_frac_squared )then - om_frac = (soilstate_inst%cellorg_col(c,lev)/organic_max)**2._r8 + om_frac = min( (soilstate_inst%cellorg_col(c,lev)/organic_max)**2._r8, 1._r8) else - om_frac = soilstate_inst%cellorg_col(c,lev)/organic_max + om_frac = min(soilstate_inst%cellorg_col(c,lev)/organic_max, 1._r8) end if else clay = soilstate_inst%cellclay_col(c,nlevsoi) @@ -604,8 +605,9 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) bd = (1._r8-soilstate_inst%watsat_col(c,lev))*params_inst%pd - soilstate_inst%watsat_col(c,lev) = params_inst%watsat_adjustfactor * ( (1._r8 - om_frac) * & - soilstate_inst%watsat_col(c,lev) + om_watsat_lake * om_frac ) + ! do not allow watsat_sf to push watsat above 0.93 + soilstate_inst%watsat_col(c,lev) = min(params_inst%watsat_adjustfactor * ( (1._r8 - om_frac) * & + soilstate_inst%watsat_col(c,lev) + om_watsat_lake * om_frac), 0.93_r8) tkm = (1._r8-om_frac)*(params_inst%tkd_sand*sand+params_inst%tkd_clay*clay)/(sand+clay) + & params_inst%tkm_om * om_frac ! W/(m K) @@ -616,7 +618,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) soilstate_inst%sucsat_col(c,lev) = params_inst%sucsat_adjustfactor * ( (1._r8-om_frac) * & soilstate_inst%sucsat_col(c,lev) + om_sucsat_lake * om_frac ) - xksat = 0.0070556 *( 10.**(-0.884+0.0153*sand) ) ! mm/s + xksat = 0.0070556_r8 *( 10._r8**(-0.884_r8+0.0153_r8*sand) ) ! mm/s ! perc_frac is zero unless perf_frac greater than percolation threshold if (om_frac > pc_lake) then @@ -631,7 +633,7 @@ subroutine SoilStateInitTimeConst(bounds, soilstate_inst, nlfilename) ! uncon_hksat is series addition of mineral/organic conductivites if (om_frac < 1._r8) then - xksat = 0.0070556 *( 10.**(-0.884+0.0153*sand) ) ! mm/s + xksat = 0.0070556_r8 *( 10._r8**(-0.884_r8+0.0153_r8*sand) ) ! mm/s uncon_hksat = uncon_frac/((1._r8-om_frac)/xksat + ((1._r8-perc_frac)*om_frac)/om_hksat_lake) else uncon_hksat = 0._r8 diff --git a/src/biogeophys/SurfaceWaterMod.F90 b/src/biogeophys/SurfaceWaterMod.F90 index 768e9a8445..031313e9e5 100644 --- a/src/biogeophys/SurfaceWaterMod.F90 +++ b/src/biogeophys/SurfaceWaterMod.F90 @@ -211,19 +211,19 @@ subroutine BulkDiag_FracH2oSfc(bounds, num_soilc, filter_soilc, & if (h2osfc(c) > min_h2osfc) then ! a cutoff is needed for numerical reasons...(nonconvergence after 5 iterations) - d=0.0 + d=0.0_r8 sigma=1.0e3 * micro_sigma(c) ! convert to mm do l=1,10 - fd = 0.5*d*(1.0_r8+erf(d/(sigma*sqrt(2.0)))) & - +sigma/sqrt(2.0*shr_const_pi)*exp(-d**2/(2.0*sigma**2)) & + fd = 0.5_r8*d*(1.0_r8+erf(d/(sigma*sqrt(2.0_r8)))) & + +sigma/sqrt(2.0_r8*shr_const_pi)*exp(-d**2/(2.0_r8*sigma**2)) & -h2osfc(c) - dfdd = 0.5*(1.0_r8+erf(d/(sigma*sqrt(2.0)))) + dfdd = 0.5_r8*(1.0_r8+erf(d/(sigma*sqrt(2.0_r8)))) d = d - fd/dfdd enddo !-- update the submerged areal fraction using the new d value - frac_h2osfc(c) = 0.5*(1.0_r8+erf(d/(sigma*sqrt(2.0)))) + frac_h2osfc(c) = 0.5_r8*(1.0_r8+erf(d/(sigma*sqrt(2.0_r8)))) qflx_too_small_h2osfc_to_soil(c) = 0._r8 @@ -453,7 +453,7 @@ subroutine QflxH2osfcSurf(bounds, num_hydrologyc, filter_hydrologyc, & ! limit runoff to value of storage above S(pc) if(h2osfc(c) > h2osfc_thresh(c) .and. h2osfcflag/=0) then ! spatially variable k_wet - k_wet=1.0e-4_r8 * sin((rpi/180.) * topo_slope(c)) + k_wet=1.0e-4_r8 * sin((rpi/180._r8) * topo_slope(c)) qflx_h2osfc_surf(c) = k_wet * frac_infclust * (h2osfc(c) - h2osfc_thresh(c)) qflx_h2osfc_surf(c)=min(qflx_h2osfc_surf(c),(h2osfc(c) - h2osfc_thresh(c))/dtime) diff --git a/src/biogeophys/TemperatureType.F90 b/src/biogeophys/TemperatureType.F90 index 56a16f1c44..f32ec83ac6 100644 --- a/src/biogeophys/TemperatureType.F90 +++ b/src/biogeophys/TemperatureType.F90 @@ -709,11 +709,11 @@ subroutine InitCold(this, bounds, & ! Set road top layer to initial air temperature and interpolate other ! layers down to 20C in bottom layer do j = 1, nlevgrnd - this%t_soisno_col(c,j) = 297.56 - (j-1) * ((297.56-293.16)/(nlevgrnd-1)) + this%t_soisno_col(c,j) = 297.56_r8 - (j-1) * ((297.56_r8-293.16_r8)/(nlevgrnd-1)) end do ! Set wall and roof layers to initial air temperature else if (col%itype(c) == icol_sunwall .or. col%itype(c) == icol_shadewall .or. col%itype(c) == icol_roof) then - this%t_soisno_col(c,1:nlevurb) = 297.56 + this%t_soisno_col(c,1:nlevurb) = 297.56_r8 else this%t_soisno_col(c,1:nlevgrnd) = 283._r8 end if @@ -722,11 +722,11 @@ subroutine InitCold(this, bounds, & ! Set road top layer to initial air temperature and interpolate other ! layers down to 22C in bottom layer do j = 1, nlevgrnd - this%t_soisno_col(c,j) = 289.46 - (j-1) * ((289.46-295.16)/(nlevgrnd-1)) + this%t_soisno_col(c,j) = 289.46_r8 - (j-1) * ((289.46_r8-295.16_r8)/(nlevgrnd-1)) end do else if (col%itype(c) == icol_sunwall .or. col%itype(c) == icol_shadewall .or. col%itype(c) == icol_roof) then ! Set wall and roof layers to initial air temperature - this%t_soisno_col(c,1:nlevurb) = 289.46 + this%t_soisno_col(c,1:nlevurb) = 289.46_r8 else this%t_soisno_col(c,1:nlevgrnd) = 283._r8 end if @@ -807,27 +807,27 @@ subroutine InitCold(this, bounds, & this%t_stem_patch(p) = this%t_veg_patch(p) if (use_vancouver) then - this%t_ref2m_patch(p) = 297.56 + this%t_ref2m_patch(p) = 297.56_r8 else if (use_mexicocity) then - this%t_ref2m_patch(p) = 289.46 + this%t_ref2m_patch(p) = 289.46_r8 else this%t_ref2m_patch(p) = 283._r8 end if if (lun%urbpoi(l)) then if (use_vancouver) then - this%t_ref2m_u_patch(p) = 297.56 + this%t_ref2m_u_patch(p) = 297.56_r8 else if (use_mexicocity) then - this%t_ref2m_u_patch(p) = 289.46 + this%t_ref2m_u_patch(p) = 289.46_r8 else this%t_ref2m_u_patch(p) = 283._r8 end if else if (.not. lun%ifspecial(l)) then if (use_vancouver) then - this%t_ref2m_r_patch(p) = 297.56 + this%t_ref2m_r_patch(p) = 297.56_r8 else if (use_mexicocity) then - this%t_ref2m_r_patch(p) = 289.46 + this%t_ref2m_r_patch(p) = 289.46_r8 else this%t_ref2m_r_patch(p) = 283._r8 end if diff --git a/src/main/FuncPedotransferMod.F90 b/src/main/FuncPedotransferMod.F90 index 41e751344e..f775b688c4 100644 --- a/src/main/FuncPedotransferMod.F90 +++ b/src/main/FuncPedotransferMod.F90 @@ -72,10 +72,10 @@ subroutine pedotransf_cosby1984_table4(sand, clay, watsat, bsw, sucsat, xksat) real(r8), intent(out):: xksat !mm/s, saturated hydraulic conductivity !Cosby et al. Table 4 - watsat = 0.505_r8-0.00142_r8*sand-0.00037*clay - bsw = 3.10+0.157*clay-0.003*sand - sucsat = 10._r8 * ( 10._r8**(1.54_r8-0.0095_r8*sand+0.0063*(100._r8-sand-clay))) - xksat = 0.0070556 *(10.**(-0.60+0.0126*sand-0.0064*clay) ) !mm/s now use table 4. + watsat = 0.505_r8-0.00142_r8*sand-0.00037_r8*clay + bsw = 3.10_r8+0.157_r8*clay-0.003_r8*sand + sucsat = 10._r8 * ( 10._r8**(1.54_r8-0.0095_r8*sand+0.0063_r8*(100._r8-sand-clay))) + xksat = 0.0070556_r8 *(10._r8**(-0.60_r8+0.0126_r8*sand-0.0064_r8*clay) ) !mm/s now use table 4. end subroutine pedotransf_cosby1984_table4 @@ -96,9 +96,9 @@ subroutine pedotransf_cosby1984_table5(sand, clay, watsat, bsw, sucsat, xksat) !Cosby et al. Table 5 watsat = 0.489_r8 - 0.00126_r8*sand - bsw = 2.91 + 0.159*clay + bsw = 2.91_r8 + 0.159_r8*clay sucsat = 10._r8 * ( 10._r8**(1.88_r8-0.0131_r8*sand) ) - xksat = 0.0070556 *( 10.**(-0.884+0.0153*sand) ) ! mm/s, from table 5 + xksat = 0.0070556_r8 *( 10._r8**(-0.884_r8+0.0153_r8*sand) ) ! mm/s, from table 5 end subroutine pedotransf_cosby1984_table5 @@ -118,10 +118,10 @@ subroutine pedotransf_noilhan_lacarrere1995(sand, clay, watsat, bsw, sucsat, xks real(r8), intent(out):: xksat !mm/s, saturated hydraulic conductivity !Noilhan and Lacarrere, 1995 - watsat = -0.00108*sand+0.494305 - bsw = 0.137*clay + 3.501 - sucsat = 10._r8**(-0.0088*sand+2.85) - xksat = 10._r8**(-0.0582*clay-0.00091*sand+0.000529*clay**2._r8-0.0001203*sand**2._r8-1.38) + watsat = -0.00108_r8*sand+0.494305_r8 + bsw = 0.137_r8*clay + 3.501_r8 + sucsat = 10._r8**(-0.0088_r8*sand+2.85_r8) + xksat = 10._r8**(-0.0582_r8*clay-0.00091_r8*sand+0.000529_r8*clay**2._r8-0.0001203_r8*sand**2._r8-1.38_r8) end subroutine pedotransf_noilhan_lacarrere1995 !------------------------------------------------------------------------------------------ function get_ipedof(soil_order)result(ipedof) diff --git a/src/main/clm_varcon.F90 b/src/main/clm_varcon.F90 index 9f66d335ad..895c7a4b80 100644 --- a/src/main/clm_varcon.F90 +++ b/src/main/clm_varcon.F90 @@ -44,8 +44,8 @@ module clm_varcon ! Initialize physical constants !------------------------------------------------------------------ - real(r8), public, parameter :: pc = 0.4 ! threshold probability - real(r8), public, parameter :: mu = 0.13889 ! connectivity exponent + real(r8), public, parameter :: pc = 0.4_r8 ! threshold probability + real(r8), public, parameter :: mu = 0.13889_r8 ! connectivity exponent real(r8), public, parameter :: secsphr = 3600._r8 ! Seconds in an hour integer, public, parameter :: isecsphr = int(secsphr) ! Integer seconds in an hour integer, public, parameter :: isecspmin= 60 ! Integer seconds in a minute @@ -126,9 +126,9 @@ module clm_varcon real(r8), public, parameter :: c_to_b = 2.0_r8 ! conversion between mass carbon and total biomass (g biomass /g C) !!! C13 - real(r8), public, parameter :: preind_atm_del13c = -6.0 ! preindustrial value for atmospheric del13C - real(r8), public, parameter :: preind_atm_ratio = SHR_CONST_PDB + (preind_atm_del13c * SHR_CONST_PDB)/1000.0 ! 13C/12C - real(r8), public :: c13ratio = preind_atm_ratio/(1.0+preind_atm_ratio) ! 13C/(12+13)C preind atmosphere + real(r8), public, parameter :: preind_atm_del13c = -6.0_r8 ! preindustrial value for atmospheric del13C + real(r8), public, parameter :: preind_atm_ratio = SHR_CONST_PDB + (preind_atm_del13c * SHR_CONST_PDB)/1000.0_r8 ! 13C/12C + real(r8), public :: c13ratio = preind_atm_ratio/(1.0_r8+preind_atm_ratio) ! 13C/(12+13)C preind atmosphere ! typical del13C for C3 photosynthesis (permil, relative to PDB) real(r8), public, parameter :: c3_del13c = -28._r8 @@ -171,7 +171,7 @@ module clm_varcon real(r8), public, parameter :: dens_floor = 2.35e3_r8 ! density of floor - concrete (Salmanca et al. 2010, TAC) (kg m-3) real(r8), public, parameter :: sh_floor = 880._r8 ! specific heat of floor - concrete (Salmanca et al. 2010, TAC) (J kg-1 K-1) real(r8), public :: cp_floor = dens_floor*sh_floor ! volumetric heat capacity of floor - concrete (Salmanca et al. 2010, TAC) (J m-3 K-1) - real(r8), public :: vent_ach = 0.3 ! ventilation rate (air exchanges per hour) + real(r8), public :: vent_ach = 0.3_r8 ! ventilation rate (air exchanges per hour) real(r8), public :: wasteheat_limit = 100._r8 ! limit on wasteheat (W/m2) diff --git a/src/soilbiogeochem/SoilBiogeochemNitrifDenitrifMod.F90 b/src/soilbiogeochem/SoilBiogeochemNitrifDenitrifMod.F90 index 784b90719b..0bdbaa2032 100644 --- a/src/soilbiogeochem/SoilBiogeochemNitrifDenitrifMod.F90 +++ b/src/soilbiogeochem/SoilBiogeochemNitrifDenitrifMod.F90 @@ -249,7 +249,7 @@ subroutine SoilBiogeochemNitrifDenitrif(bounds, num_soilc, filter_soilc, & organic_max = CNParamsShareInst%organic_max - pH(bounds%begc:bounds%endc) = 6.5 !!! set all soils with the same pH as placeholder here + pH(bounds%begc:bounds%endc) = 6.5_r8 !!! set all soils with the same pH as placeholder here co2diff_con(1) = 0.1325_r8 co2diff_con(2) = 0.0009_r8 @@ -310,7 +310,7 @@ subroutine SoilBiogeochemNitrifDenitrif(bounds, num_soilc, filter_soilc, & k_nitr_t_vr(c,j) = min(t_scalar(c,j), 1._r8) ! ph function from Parton et al., (2001, 1996) - k_nitr_ph_vr(c,j) = 0.56 + atan(rpi * 0.45 * (-5.+ pH(c)))/rpi + k_nitr_ph_vr(c,j) = 0.56_r8 + atan(rpi * 0.45_r8 * (-5._r8+ pH(c)))/rpi ! moisture function-- assume the same moisture function as limits heterotrophic respiration ! Parton et al. base their nitrification- soil moisture rate constants based on heterotrophic rates-- can we do the same? @@ -383,7 +383,7 @@ subroutine SoilBiogeochemNitrifDenitrif(bounds, num_soilc, filter_soilc, & fr_WFPS(c,j) = max(0.1_r8, 0.015_r8 * wfps_vr(c,j) - 0.32_r8) ! final ratio expression - n2_n2o_ratio_denit_vr(c,j) = max(0.16*ratio_k1(c,j), ratio_k1(c,j)*exp(-0.8 * ratio_no3_co2(c,j))) * fr_WFPS(c,j) + n2_n2o_ratio_denit_vr(c,j) = max(0.16_r8*ratio_k1(c,j), ratio_k1(c,j)*exp(-0.8_r8 * ratio_no3_co2(c,j))) * fr_WFPS(c,j) end do