bugfix for bb_slope, carea, and icefilter

Merge branch 'rgknox-bbfix-careafix-icefiltfix'

This PR bundles 3 fixes that address: #73, #69, #44

The fix to 73 is the only one that would be expected to have b4b
regressions. I performed baseline simulation comparsisons between
f1a14d6 and 18613d1, and tests confirmed b4b on all expected
passes. One extra step was necessary, in that I needed to update the
parameter file values of BB_slope to match what was previously hard
coded (a value of 9). The current value in the default parameter file
is 8, we can certaintly change this going forward. Not changing this
now.

The other two fixes, #69 and #44 are not supposed to generate b4b
results, and they don't. 1x1 brazil simulations were also run on eddi
to make sure that the non-b4b changes continue to generate very
similar projections of forest composition and structure, as well as
flux variables. They did.

Fixes: #73
Fixes: #69
Fixes: #44

User interface changes?: no

Code review: rgknox

Test suite: ed - lawrencium-lr3 intel, eddi (PC) gnu (visualizations);
ed - yellowstone gnu, intel, pgi

Test baseline: 18613d1
Test namelist changes: none
Test answer changes: yes, see above

Test summary: pass except for #14 known failures in f09 and f19
restart, and gnu f10 restart #43.

components/clm/src/ED/biogeophys/EDPhotosynthesisMod.F90

@@ -95,10 +95,12 @@ subroutine Photosynthesis_ED (bounds, fn, filterp, esat_tv, eair, oair, cair, &
     real(r8) :: kc( bounds%begp:bounds%endp )     ! Michaelis-Menten constant for CO2 (Pa)
     real(r8) :: ko( bounds%begp:bounds%endp )     ! Michaelis-Menten constant for O2 (Pa)
     real(r8) :: co2_cp( bounds%begp:bounds%endp ) ! CO2 compensation point (Pa)
+
+    ! ---------------------------------------------------------------
+    ! TO-DO: bbbopt is slated to be transferred to the parameter file
+    ! ----------------------------------------------------------------
     real(r8) :: bbbopt(psn_type)                  ! Ball-Berry minimum leaf conductance, unstressed (umol H2O/m**2/s)
     real(r8) :: bbb(mxpft)                        ! Ball-Berry minimum leaf conductance (umol H2O/m**2/s)
-    real(r8) :: mbbopt(psn_type)                  ! Ball-Berry slope of conductance-photosynthesis relationship, unstressed
-    real(r8) :: mbb(mxpft)                        ! Ball-Berry slope of conductance-photosynthesis relationship
 
     real(r8) :: kn(mxpft)                         ! leaf nitrogen decay coefficient
     real(r8) :: vcmax25top(mxpft)                 ! canopy top: maximum rate of carboxylation at 25C (umol CO2/m**2/s)
@@ -306,12 +308,11 @@ subroutine Photosynthesis_ED (bounds, fn, filterp, esat_tv, eair, oair, cair, &
       qe(1) = 0._r8
       theta_cj(1) = 0.98_r8
       bbbopt(1) = 10000._r8
-      mbbopt(1) = 9._r8
 
       qe(2) = 0.05_r8
       theta_cj(2) = 0.80_r8
       bbbopt(2) = 40000._r8
-      mbbopt(2) = 4._r8
+
 
       do f = 1,fn
          p = filterp(f)
@@ -355,17 +356,6 @@ subroutine Photosynthesis_ED (bounds, fn, filterp, esat_tv, eair, oair, cair, &
                enddo !ft
             enddo !CL
 
-            ! Soil water stress applied to Ball-Berry parameters
-            do FT = 1,numpft_ed
-               if (nint(c3psn(FT)) == 1)then
-                  ps = 1
-               else
-                  ps = 2
-               end if
-               bbb(FT) = max (bbbopt(ps)*currentPatch%btran_ft(FT), 1._r8)
-
-               mbb(FT) = bb_slope(ft) ! mbbopt(ps)
-            end do
 
             ! kc, ko, currentPatch, from: Bernacchi et al (2001) Plant, Cell and Environment 24:253-259
             !
@@ -410,25 +400,24 @@ subroutine Photosynthesis_ED (bounds, fn, filterp, esat_tv, eair, oair, cair, &
 
             currentPatch => map_clmpatch_to_edpatch(sites(s), p) 
 
-            do FT = 1,numpft_ed
+            NCL_p = currentPatch%NCL_p
+
+            do FT = 1,numpft_ed !calculate patch and pft specific propserties at canopy top. 
+
                if (nint(c3psn(FT)) == 1)then
                   ps = 1
                else
                   ps = 2
                end if
                bbb(FT) = max (bbbopt(ps)*currentPatch%btran_ft(FT), 1._r8)
-               mbb(FT) = mbbopt(ps)
 
+               ! THIS CALL APPEARS TO BE REDUNDANT WITH LINE 650 (RGK)
                if (nint(c3psn(FT)) == 1)then
                   ci(:,FT,:) = 0.7_r8 * cair(p)
                else
                   ci(:,FT,:) = 0.4_r8 * cair(p)
                end if
-            enddo
-
-            NCL_p = currentPatch%NCL_p
 
-            do FT = 1,numpft_ed !calculate patch and pft specific propserties at canopy top. 
 
                ! Leaf nitrogen concentration at the top of the canopy (g N leaf / m**2 leaf)
                lnc(FT) = 1._r8 / (slatop(FT) * leafcn(FT))
@@ -647,6 +636,8 @@ subroutine Photosynthesis_ED (bounds, fn, filterp, esat_tv, eair, oair, cair, &
                                  je = min(r1,r2)
 
                                  ! Iterative loop for ci beginning with initial guess
+                                 ! THIS CALL APPEARS TO BE REDUNDANT WITH LINE 423 (RGK)
+
                                  if (nint(c3psn(FT)) == 1)then
                                     ci(cl,ft,iv) = 0.7_r8 * cair(p)
                                  else
@@ -719,8 +710,8 @@ subroutine Photosynthesis_ED (bounds, fn, filterp, esat_tv, eair, oair, cair, &
                                     cs = cair(p) - 1.4_r8/gb_mol * an(cl,ft,iv) * forc_pbot(c)
                                     cs = max(cs,1.e-06_r8)
                                     aquad = cs
-                                    bquad = cs*(gb_mol - bbb(FT)) - mbb(FT)*an(cl,ft,iv)*forc_pbot(c)
-                                    cquad = -gb_mol*(cs*bbb(FT) + mbb(FT)*an(cl,ft,iv)*forc_pbot(c)*ceair/esat_tv(p))
+                                    bquad = cs*(gb_mol - bbb(FT)) - bb_slope(ft)*an(cl,ft,iv)*forc_pbot(c)
+                                    cquad = -gb_mol*(cs*bbb(FT) + bb_slope(ft)*an(cl,ft,iv)*forc_pbot(c)*ceair/esat_tv(p))
                                     call quadratic (aquad, bquad, cquad, r1, r2)
                                     gs_mol = max(r1,r2)
 
@@ -788,7 +779,7 @@ subroutine Photosynthesis_ED (bounds, fn, filterp, esat_tv, eair, oair, cair, &
 
                                  ! Compare with Ball-Berry model: gs_mol = m * an * hs/cs p + b
                                  hs = (gb_mol*ceair + gs_mol*esat_tv(p)) / ((gb_mol+gs_mol)*esat_tv(p))
-                                 gs_mol_err = mbb(FT)*max(an(cl,ft,iv), 0._r8)*hs/cs*forc_pbot(c) + bbb(FT)
+                                 gs_mol_err = bb_slope(ft)*max(an(cl,ft,iv), 0._r8)*hs/cs*forc_pbot(c) + bbb(FT)
 
                                  if (abs(gs_mol-gs_mol_err) > 1.e-01_r8) then
                                     write (iulog,*) 'CF: Ball-Berry error check - stomatal conductance error:'

components/clm/src/ED/main/EDCLMLinkMod.F90

@@ -1093,10 +1093,10 @@ subroutine ed_clm_link( this, bounds, sites, nsites, fcolumn, waterstate_inst, c
 
                   currentCohort%b = currentCohort%balive+currentCohort%bdead+currentCohort%bstore
                   currentCohort%treelai = tree_lai(currentCohort)
-                  ! Why is currentCohort%c_area used and then reset in the
-                  ! following line?
-                  canopy_leaf_area = canopy_leaf_area + currentCohort%treelai *currentCohort%c_area
+
                   currentCohort%c_area = c_area(currentCohort)
+                  canopy_leaf_area = canopy_leaf_area + currentCohort%treelai *currentCohort%c_area
+
 
                   if(currentCohort%canopy_layer==1)then
                      currentPatch%total_canopy_area = currentPatch%total_canopy_area + currentCohort%c_area

components/clm/src/main/clm_driver.F90

@@ -845,7 +845,7 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate)
                atm2lnd_inst, soilstate_inst, temperature_inst,                   &
                waterstate_inst, canopystate_inst)
 
-          call setFilters( bounds_clump, glc2lnd_inst%icemask_grc )
+          call setFilters( bounds_clump, glc2lnd_inst%icemask_grc(bounds_clump%begg:bounds_clump%endg ))
 
        end if ! use_ed branch