Merge pull request #11 from rgknox/mpaiao-pr-allom-modes-moreparams

added parameters and merged up

Author: mpaiao

biogeophys/FatesPlantHydraulicsMod.F90

@@ -937,6 +937,7 @@ subroutine UpdatePlantHydrLenVol(ccohort,csite_hydr)
     real(r8), parameter :: min_trim      = 0.1_r8   ! The lower cap on trimming function used
     ! to estimate maximum leaf carbon
 
+
     ccohort_hydr => ccohort%co_hydr
     ft           = ccohort%pft
     nlevrhiz     = csite_hydr%nlevrhiz
@@ -980,13 +981,13 @@ subroutine UpdatePlantHydrLenVol(ccohort,csite_hydr)
     ! Get the target, or rather, maximum leaf carrying capacity of plant
     ! Lets also avoid super-low targets that have very low trimming functions
 
+    ! efleaf_coh hard-coded to 1 in the call below to avoid zero leaf volume
     call bleaf(ccohort%dbh,ccohort%pft,ccohort%crowndamage, &
-         max(ccohort%canopy_trim,min_trim),ccohort%efleaf_coh, leaf_c_target)
+         max(ccohort%canopy_trim,min_trim),1.0_r8, leaf_c_target)
+
+    ccohort_hydr%v_ag(1:n_hypool_leaf) = max(leaf_c,min_leaf_frac*leaf_c_target) * &
+         prt_params%c2b(ft) / denleaf/ real(n_hypool_leaf,r8)
 
-    if( (ccohort%status_coh == leaves_on) .or. ccohort_hydr%is_newly_recruited ) then
-       ccohort_hydr%v_ag(1:n_hypool_leaf) = max(leaf_c,min_leaf_frac*leaf_c_target) * &
-            prt_params%c2b(ft) / denleaf/ real(n_hypool_leaf,r8)
-    end if
 
     ! Step sapwood volume
     ! -----------------------------------------------------------------------------------

biogeophys/FatesPlantRespPhotosynthMod.F90

@@ -38,6 +38,7 @@ module FATESPlantRespPhotosynthMod
   use FatesInterfaceTypesMod, only : hlm_parteh_mode
   use FatesInterfaceTypesMod, only : numpft
   use FatesInterfaceTypesMod, only : nleafage
+  use FatesUtilsMod,          only : QuadraticRoots => QuadraticRootsSridharachary
   use EDParamsMod,           only : maxpft
   use EDParamsMod,       only : nlevleaf
   use EDParamsMod,       only : nclmax
@@ -1379,7 +1380,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl,         &  ! in
            aquad = theta_psii
            bquad = -(qabs + jmax)
            cquad = qabs * jmax
-           call quadratic_f (aquad, bquad, cquad, r1, r2)
+           call QuadraticRoots(aquad, bquad, cquad, r1, r2)
            je = min(r1,r2)
 
            ! Initialize intercellular co2
@@ -1409,7 +1410,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl,         &  ! in
                     aquad = theta_cj_c3
                     bquad = -(ac + aj)
                     cquad = ac * aj
-                    call quadratic_f (aquad, bquad, cquad, r1, r2)
+                    call QuadraticRoots(aquad, bquad, cquad, r1, r2)
                     agross = min(r1,r2)
 
               else
@@ -1440,13 +1441,13 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl,         &  ! in
                  aquad = theta_cj_c4
                  bquad = -(ac + aj)
                  cquad = ac * aj
-                 call quadratic_f (aquad, bquad, cquad, r1, r2)
+                 call QuadraticRoots(aquad, bquad, cquad, r1, r2)
                  ai = min(r1,r2)
 
                  aquad = theta_ip
                  bquad = -(ai + ap)
                  cquad = ai * ap
-                 call quadratic_f (aquad, bquad, cquad, r1, r2)
+                 call QuadraticRoots(aquad, bquad, cquad, r1, r2)
                  agross = min(r1,r2)
 
               end if
@@ -1485,7 +1486,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl,         &  ! in
                       (2.0*stomatal_intercept_btran + term * &
                       (1.0 - medlyn_slope(ft)* medlyn_slope(ft) / vpd)) * term
 
-                 call quadratic_f (aquad, bquad, cquad, r1, r2)
+                 call QuadraticRoots(aquad, bquad, cquad, r1, r2)
                  gs_mol = max(r1,r2)
 
               else if ( stomatal_model == ballberry_model ) then         !stomatal conductance calculated from Ball et al. (1987)
@@ -1494,7 +1495,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl,         &  ! in
                  cquad = -gb_mol*(leaf_co2_ppress*stomatal_intercept_btran + &
                       bb_slope(ft)*anet*can_press * ceair/ veg_esat )
 
-                 call quadratic_f (aquad, bquad, cquad, r1, r2)
+                 call QuadraticRoots(aquad, bquad, cquad, r1, r2)
                  gs_mol = max(r1,r2)
               end if
 
@@ -1922,145 +1923,6 @@ end function fth25_f
 
   ! =====================================================================================
 
-  subroutine quadratic_f (a, b, c, r1, r2)
-    !
-    ! !DESCRIPTION:
-    !==============================================================================!
-    !----------------- Solve quadratic equation for its two roots -----------------!
-    !==============================================================================!
-    ! This solution is mostly derived from:
-    ! Press WH, Teukolsky SA, Vetterling WT, Flannery BP. 1992. Numerical Recipes
-    !    in Fortran77: The Art of Scientific Computing. 2nd edn. Cambridge 
-    !    University Press, Cambridge UK, ISBN 0-521-43064-X.
-    !    Available at: http://numerical.recipes/oldverswitcher.html, section 5.6.
-    !
-    ! !REVISION HISTORY:
-    ! 4/5/10: Adapted from /home/bonan/ecm/psn/An_gs_iterative.f90 by Keith Oleson
-    ! 7/23/16: Copied over from CLM by Ryan Knox
-    ! 12/30/23: Instead of issuing errors when a=0, solve the trivial cases too.
-    !           Check determinant sign, and stop the run when it is negative.
-    !
-    ! !USES:
-    !
-    ! !ARGUMENTS:
-    real(r8), intent(in)  :: a,b,c       ! Terms for quadratic equation
-    real(r8), intent(out) :: r1,r2       ! Roots of quadratic equation
-    !
-    ! !LOCAL VARIABLES:
-    real(r8) :: discriminant             ! Discriminant
-    real(r8) :: q                        ! Temporary term for quadratic solution
-    logical  :: a_offzero                ! Is a close to zero?
-    logical  :: b_offzero                ! Is b close to zero?
-    logical  :: c_offzero                ! Is c close to zero?
-    ! ! Local constants:
-    real(r8), parameter :: discard = 1.e36_r8 ! Large number for discarding answer
-    !------------------------------------------------------------------------------
-
-    ! Save logical tests.
-    a_offzero = abs(a) > nearzero
-    b_offzero = abs(b) > nearzero
-    c_offzero = abs(c) > nearzero
-
-    if (a_offzero .and. ( b_offzero .or. c_offzero ) ) then
-       ! Quadratic equation with two non-zero solutions (but may be complex solutions)
-       discriminant = b*b - 4._r8 * a * c
-
-       ! Proceed only when the discriminant is non-negative or only tiny negative
-       if (discriminant >= - nearzero) then
-          ! Coerce discriminant to non-negative
-          discriminant = max(0._r8,discriminant)
-
-          ! Find q as in the numerical recipes. If b or c are non-zero, q cannot
-          ! be zero, no need for additional checks.
-          q  = - 0.5_r8 * (b + sign(sqrt(discriminant),b))
-          r1 = q / a
-          r2 = c / q
-       else
-          ! Negative discriminant, stop the run.
-          write (fates_log(),'(a)') '---~---'
-          write (fates_log(),'(a)') ' Fatal error!'
-          write (fates_log(),'(a)') ' Quadratic equation discriminant is negative.'
-          write (fates_log(),'(a)') '---~---'
-          write (fates_log(),'(a,1x,es12.5)') ' a            = ',a
-          write (fates_log(),'(a,1x,es12.5)') ' b            = ',b
-          write (fates_log(),'(a,1x,es12.5)') ' c            = ',c
-          write (fates_log(),'(a,1x,es12.5)') ' discriminant = ',discriminant
-          call endrun(msg=errMsg(sourcefile, __LINE__))
-       end if
-    else if (a_offzero) then
-       ! b and c are nearly zero. Both roots must be zero.
-       r1 = 0._r8
-       r2 = 0._r8
-    else if (b_offzero) then
-       ! "a" is not zero, not a true quadratic equation. Single root.
-       r1 = - c / b
-       r2 = discard
-    else
-       ! Both a and b are zero, this really doesn't make any sense and should never
-       ! happen. If it does, issue an error and stop the run.
-       write (fates_log(),'(a)') '---~---'
-       write (fates_log(),'(a)') ' Fatal error!'
-       write (fates_log(),'(a)') ' This solver expects ''a'' and/or ''b'' to be non-zero.'
-       write (fates_log(),'(a)') '---~---'
-       write (fates_log(),'(a,1x,es12.5)') ' a = ',a
-       write (fates_log(),'(a,1x,es12.5)') ' b = ',b
-       write (fates_log(),'(a,1x,es12.5)') ' c = ',c
-       write (fates_log(),'(a)') '---~---'
-       call endrun(msg=errMsg(sourcefile, __LINE__))
-    end if
-
-    return
-  end subroutine quadratic_f
-
-  ! ====================================================================================
-
-  subroutine quadratic_fast (a, b, c, r1, r2)
-    !
-    ! !DESCRIPTION:
-    !==============================================================================!
-    !----------------- Solve quadratic equation for its two roots -----------------!
-    ! THIS METHOD SIMPLY REMOVES THE DIV0 CHECK AND ERROR REPORTING                !
-    !==============================================================================!
-    ! Solution from Press et al (1986) Numerical Recipes: The Art of Scientific
-    ! Computing (Cambridge University Press, Cambridge), pp. 145.
-    !
-    ! !REVISION HISTORY:
-    ! 4/5/10: Adapted from /home/bonan/ecm/psn/An_gs_iterative.f90 by Keith Oleson
-    ! 7/23/16: Copied over from CLM by Ryan Knox
-    !
-    ! !USES:
-    !
-    ! !ARGUMENTS:
-    real(r8), intent(in)  :: a,b,c       ! Terms for quadratic equation
-    real(r8), intent(out) :: r1,r2       ! Roots of quadratic equation
-    !
-    ! !LOCAL VARIABLES:
-    real(r8) :: q                        ! Temporary term for quadratic solution
-    !------------------------------------------------------------------------------
-
-    !  if (a == 0._r8) then
-    !     write (fates_log(),*) 'Quadratic solution error: a = ',a
-    !     call endrun(msg=errMsg(sourcefile, __LINE__))
-    !  end if
-
-    if (b >= 0._r8) then
-       q = -0.5_r8 * (b + sqrt(b*b - 4._r8*a*c))
-    else
-       q = -0.5_r8 * (b - sqrt(b*b - 4._r8*a*c))
-    end if
-
-    r1 = q / a
-    !  if (q /= 0._r8) then
-    r2 = c / q
-    !  else
-    !     r2 = 1.e36_r8
-    !  end if
-
-  end subroutine quadratic_fast
-
-
-  ! ====================================================================================
-
   subroutine UpdateCanopyNCanNRadPresent(currentPatch)
 
     ! ---------------------------------------------------------------------------------