Merge pull request #1384 from ekluzek/dblconstantsandlimit

Change some constants to double precision, and add some soil limits

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

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

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

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(                                                               & 

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

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).

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

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

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__)

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)

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

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

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

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

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

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)

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

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,21 +725,21 @@ 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
             if(.not. use_cn) then
                ! 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

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
 
     !-----------------------------------------------------------------------