debugged bisection, leaf unit tests generating stable results

NGEET · Oct 22, 2024 · c7af3ef · c7af3ef
1 parent bd89109
commit c7af3ef
Show file tree

Hide file tree

Showing 2 changed files with 24 additions and 33 deletions.
diff --git a/biogeophys/LeafBiophysicsMod.F90 b/biogeophys/LeafBiophysicsMod.F90
@@ -202,6 +202,16 @@ module LeafBiophysicsMod
 
   subroutine StomatalCondMedlyn(anet,ft,veg_esat,can_vpress,stomatal_intercept_btran,leaf_co2_ppress,can_press,gb,gs)
 
+    ! -----------------------------------------------------------------------------------
+    ! Calculate stomatal conductance (gs [umol/m2/s]) via the Medlyn approach, described in:
+    ! Medlyn et al. Reconciling the optimal and empirical approaches to modelling stomatal
+    ! conductance.  Global Change Biology (2011) 17, 2134–2144, doi: 10.1111/j.1365-2486.2010.02375.x
+    !
+    ! THe implementation, and adaptation to include a leaf boundary layer in serial
+    ! resitance was taken from CLM5.0 [kPa]
+    ! -----------------------------------------------------------------------------------
+
+
     ! Input
     real(r8), intent(in) :: anet                     ! net leaf photosynthesis (umol CO2/m**2/s)
     integer, intent(in)  :: ft                       ! plant functional type index
@@ -223,20 +233,16 @@ subroutine StomatalCondMedlyn(anet,ft,veg_esat,can_vpress,stomatal_intercept_btr
     real(r8) :: aquad,bquad,cquad                    ! quadradic solve terms
     real(r8) :: r1,r2                                ! quadradic solve roots
 
+    real(r8),parameter :: min_vpd_pa = 50._r8        ! Minimum allowable vpd [Pa]
+
     ! Evaluate trival solution, if there is no positive net assimiolation
     ! the stomatal conductance is the intercept conductance
     if (anet <= nearzero) then
        gs = stomatal_intercept_btran
        return
     end if
 
-
-
-
-    ! stomatal conductance calculated from Medlyn et al. (2011), the numerical &
-    ! implementation was adapted from the equations in CLM5.0 [kPa]
-
-    vpd =  max((veg_esat - can_vpress), 50._r8) * kpa_per_pa       !addapted from CLM5. Put some constraint on VPD
+    vpd =  max((veg_esat - can_vpress), min_vpd_pa) * kpa_per_pa
 
     term = h2o_co2_stoma_diffuse_ratio * anet / (leaf_co2_ppress / can_press)
     aquad = 1.0_r8
@@ -283,6 +289,7 @@ subroutine StomatalCondBallBerry(a_gs,ft,veg_esat,can_vpress,stomatal_intercept_
 
     ! Apply a constraint to the vapor pressure
     ceair = GetConstrainedVPress(can_vpress,veg_esat)
+    ceair = can_vpress
     aquad = leaf_co2_ppress
     bquad = leaf_co2_ppress*(gb - stomatal_intercept_btran) - lb_params%bb_slope(ft) * a_gs * can_press
     cquad = -gb*(leaf_co2_ppress * stomatal_intercept_btran + &
@@ -291,20 +298,6 @@ subroutine StomatalCondBallBerry(a_gs,ft,veg_esat,can_vpress,stomatal_intercept_
     call QuadraticRoots(aquad, bquad, cquad, r1, r2)
     gs = max(r1,r2)
 
-    ! Testing poor convergence
-    !if( ceair/veg_esat < 0.04_r8 ) then
-    !   print*,"weird constrained:",ceair,veg_esat
-    !   stop
-    !end if
-
-    !if(leaf_co2_ppress < 5._r8) then
-    !   print*,"low leaf_co2_ppress:",leaf_co2_ppress
-    !   stop
-    !end if
-
-    !gs = lb_params%bb_slope(ft)*a_gs*(ceair/veg_esat)/leaf_co2_ppress*can_press + stomatal_intercept_btran
-
-
     return
   end subroutine StomatalCondBallBerry
 
@@ -647,7 +640,7 @@ subroutine CiBisection(ft,vcmax,jmax,kp,co2_cpoint,mm_kco2,mm_ko2, &
             stomatal_intercept_btran, &
             anet,agross,gs,fval_b)
 
-       if(fval_b<ci_tol) then
+       if(abs(fval_b)<ci_tol) then
           loop_continue = .false.
        end if
 
@@ -1531,26 +1524,24 @@ subroutine LowstorageMainRespReduction(frac, ft, maintresp_reduction_factor)
 
   end subroutine LowstorageMainRespReduction
 
-
-
   ! =====================================================================================
 
   real(r8) function GetConstrainedVPress(air_vpress,veg_esat) result(ceair)
 
     ! Return a constrained vapor pressure [Pa]
+    ! This function is used to make sure that Ball-Berry conductance is well
+    ! behaved.
 
-    real(r8) :: air_vpress   ! vapor pressure of the air (unconstrained) [Pa]
-    real(r8) :: veg_esat     ! saturated vapor pressure [Pa]
-
-    ! Constrain eair >= 0.05*esat_tv so that solution does not blow up. This ensures
-    ! that hs (ie air_vpress / veg_esat)  in the Ball-Berry equation does not go to zero.
-    ! Also eair <= veg_esat so that hs <= 1
+    real(r8) :: air_vpress              ! vapor pressure of the air (unconstrained) [Pa]
+    real(r8) :: veg_esat                ! saturated vapor pressure [Pa]
+    real(r8) :: min_frac_esat = 0.01_r8 ! We don't allow vapor pressures
+                                        ! below this fraction amount of saturation vapor pressure
 
-    ceair = min( max(air_vpress, 0.05_r8*veg_esat ),veg_esat )
+    ceair = min( max(air_vpress, min_frac_esat*veg_esat ),veg_esat )
 
   end function GetConstrainedVPress
 
-   ! ====================================================================================
+  ! ====================================================================================
 
   real(r8) function GetStomatalInterceptBtran(ft,btran) result(stomatal_intercept_btran)
 

diff --git a/functional_unit_testing/leaf_biophys/leaf_biophys_controls.xml b/functional_unit_testing/leaf_biophys/leaf_biophys_controls.xml
@@ -9,7 +9,7 @@
       <param name='fates_daylength_factor_switch'>1</param>
       <param name='fates_leaf_theta_cj_c3'>0.999</param>
       <param name='fates_leaf_theta_cj_c4'>0.999</param>
-      <param name='fates_leaf_stomatal_model'>1</param>
+      <param name='fates_leaf_stomatal_model'>2</param>
       <param name='fates_leaf_stomatal_assim_model'>1</param>
       <param name='fates_leaf_photo_tempsens_model'>1</param>
     </scalar_dim>