Merge pull request #10 from glemieux/fates_harvest_offmaster-deconflict

Fates harvest offmaster deconflict

biogeochem/EDPatchDynamicsMod.F90

@@ -1185,7 +1185,7 @@ subroutine check_patch_area( currentSite )
     ! !USES:
     !
     ! !ARGUMENTS:
-    type(ed_site_type), intent(in), target  :: currentSite 
+    type(ed_site_type), intent(inout), target  :: currentSite
     !
     ! !LOCAL VARIABLES:
     real(r8)                     :: areatot
@@ -1579,6 +1579,7 @@ subroutine fire_litter_fluxes(currentSite, currentPatch, newPatch, patch_site_ar
           num_dead_trees  = (currentCohort%fire_mort * &
                 currentCohort%n*patch_site_areadis/currentPatch%area)
           call AccumulateMortalityWaterStorage(currentSite,currentCohort,num_dead_trees)
+          currentCohort => currentCohort%taller
        end do
     end if
 

biogeochem/EDPhysiologyMod.F90

@@ -370,36 +370,90 @@ subroutine trim_canopy( currentSite )
     type (ed_cohort_type) , pointer :: currentCohort
     type (ed_patch_type)  , pointer :: currentPatch
 
-    integer  :: z                ! leaf layer
-    integer  :: ipft             ! pft index
-    logical  :: trimmed          ! was this layer trimmed in this year? If not expand the canopy. 
-    real(r8) :: tar_bl           ! target leaf biomass       (leaves flushed, trimmed)
-    real(r8) :: tar_bfr          ! target fine-root biomass  (leaves flushed, trimmed)
-    real(r8) :: bfr_per_bleaf    ! ratio of fine root per leaf biomass
-    real(r8) :: sla_levleaf      ! sla at leaf level z
-    real(r8) :: nscaler_levleaf  ! nscaler value at leaf level z
-    integer  :: cl               ! canopy layer index
-    real(r8) :: kn               ! nitrogen decay coefficient
-    real(r8) :: sla_max          ! Observational constraint on how large sla (m2/gC) can become
-    real(r8) :: leaf_c           ! leaf carbon [kg]
-    real(r8) :: sapw_c           ! sapwood carbon [kg]
-    real(r8) :: store_c          ! storage carbon [kg]
-    real(r8) :: struct_c         ! structure carbon [kg]
-    real(r8) :: leaf_inc         ! LAI-only portion of the vegetation increment of dinc_ed
-    real(r8) :: lai_canopy_above ! the LAI in the canopy layers above the layer of interest
-    real(r8) :: lai_layers_above ! the LAI in the leaf layers, within the current canopy, 
-                                 ! above the leaf layer of interest
-    real(r8) :: lai_current      ! the LAI in the current leaf layer
-    real(r8) :: cumulative_lai   ! the cumulative LAI, top down, to the leaf layer of interest
+    integer  :: z                     ! leaf layer
+    integer  :: ipft                  ! pft index
+    logical  :: trimmed               ! was this layer trimmed in this year? If not expand the canopy. 
+    real(r8) :: tar_bl                ! target leaf biomass       (leaves flushed, trimmed)
+    real(r8) :: tar_bfr               ! target fine-root biomass  (leaves flushed, trimmed)
+    real(r8) :: bfr_per_bleaf         ! ratio of fine root per leaf biomass
+    real(r8) :: sla_levleaf           ! sla at leaf level z
+    real(r8) :: nscaler_levleaf       ! nscaler value at leaf level z
+    integer  :: cl                    ! canopy layer index
+    real(r8) :: kn                    ! nitrogen decay coefficient
+    real(r8) :: sla_max               ! Observational constraint on how large sla (m2/gC) can become
+    real(r8) :: leaf_c                ! leaf carbon [kg]
+    real(r8) :: sapw_c                ! sapwood carbon [kg]
+    real(r8) :: store_c               ! storage carbon [kg]
+    real(r8) :: struct_c              ! structure carbon [kg]
+    real(r8) :: leaf_inc              ! LAI-only portion of the vegetation increment of dinc_ed
+    real(r8) :: lai_canopy_above      ! the LAI in the canopy layers above the layer of interest
+    real(r8) :: lai_layers_above      ! the LAI in the leaf layers, within the current canopy, 
+                                      ! above the leaf layer of interest
+    real(r8) :: lai_current           ! the LAI in the current leaf layer
+    real(r8) :: cumulative_lai        ! whole canopy cumulative LAI, top down, to the leaf layer of interest
+    real(r8) :: cumulative_lai_cohort ! cumulative LAI within the current cohort only
+
+    ! Temporary diagnostic ouptut
+    integer :: ipatch
+    integer :: icohort
+
+    ! LAPACK linear least squares fit variables
+    ! The standard equation for a linear fit, y = mx + b, is converted to a linear system, AX=B and has
+    ! the form: [n  sum(x); sum(x)  sum(x^2)] * [b; m]  = [sum(y); sum(x*y)] where 
+    ! n is the number of leaf layers 
+    ! x is yearly_net_uptake minus the leaf cost aka the net-net uptake
+    ! y is the cumulative lai for the current cohort
+    ! b is the y-intercept i.e. the cumulative lai that has zero net-net uptake
+    ! m is the slope of the linear fit
+    integer :: nll = 3                    ! Number of leaf layers to fit a regression to for calculating the optimum lai
+    character(1) :: trans = 'N'           ! Input matrix is not transposed
+
+    integer, parameter :: m = 2, n = 2    ! Number of rows and columns, respectively, in matrix A
+    integer, parameter :: nrhs = 1        ! Number of columns in matrix B and X
+    integer, parameter :: workmax = 100   ! Maximum iterations to minimize work
+
+    integer :: lda = m, ldb = n           ! Leading dimension of A and B, respectively
+    integer :: lwork                      ! Dimension of work array
+    integer :: info                       ! Procedure diagnostic ouput
+
+    real(r8) :: nnu_clai_a(m,n)           ! LHS of linear least squares fit, A matrix
+    real(r8) :: nnu_clai_b(m,nrhs)        ! RHS of linear least squares fit, B matrix
+    real(r8) :: work(workmax)             ! work array
+
+    real(r8) :: initial_trim              ! Initial trim
+    real(r8) :: optimum_trim              ! Optimum trim value 
+    real(r8) :: initial_laimem            ! Initial laimemory
+    real(r8) :: optimum_laimem            ! Optimum laimemory
 
     !----------------------------------------------------------------------
 
+    ipatch = 1 ! Start counting patches
+
     currentPatch => currentSite%youngest_patch
     do while(associated(currentPatch))
-
+
+       ! Add debug diagnstic output to determine which patch
+       if (debug) then
+          write(fates_log(),*) 'Current patch:', ipatch
+          write(fates_log(),*) 'Current patch cohorts:', currentPatch%countcohorts
+       endif
+
+       icohort = 1
+
        currentCohort => currentPatch%tallest
        do while (associated(currentCohort)) 
 
+         ! Save off the incoming trim and laimemory
+         initial_trim = currentCohort%canopy_trim
+         initial_laimem = currentCohort%laimemory
+
+          ! Add debug diagnstic output to determine which cohort
+          if (debug) then
+            write(fates_log(),*) 'Current cohort:', icohort
+            write(fates_log(),*) 'Starting canopy trim:', initial_trim
+            write(fates_log(),*) 'Starting laimemory:', currentCohort%laimemory
+          endif   
+
           trimmed = .false.
           ipft = currentCohort%pft
           call carea_allom(currentCohort%dbh,currentCohort%n,currentSite%spread,currentCohort%pft,currentCohort%c_area)
@@ -438,6 +492,10 @@ subroutine trim_canopy( currentSite )
           ! PFT-level maximum SLA value, even if under a thick canopy (same units as slatop)
           sla_max = EDPftvarcon_inst%slamax(ipft)
 
+          ! Initialize nnu_clai_a
+          nnu_clai_a(:,:) = 0._r8
+          nnu_clai_b(:,:) = 0._r8
+
           !Leaf cost vs netuptake for each leaf layer. 
           do z = 1, currentCohort%nv
 
@@ -446,14 +504,17 @@ subroutine trim_canopy( currentSite )
              leaf_inc    = dinc_ed * &
                    currentCohort%treelai/(currentCohort%treelai+currentCohort%treesai)
 
-             ! Now calculate the cumulative top-down lai of the current layer's midpoint
+             ! Now calculate the cumulative top-down lai of the current layer's midpoint within the current cohort
+             lai_layers_above      = leaf_inc * (z-1)
+             lai_current           = min(leaf_inc, currentCohort%treelai - lai_layers_above)
+             cumulative_lai_cohort = lai_layers_above + 0.5*lai_current
+
+             ! Now add in the lai above the current cohort for calculating the sla leaf level
              lai_canopy_above  = sum(currentPatch%canopy_layer_tlai(1:cl-1)) 
-             lai_layers_above  = leaf_inc * (z-1)
-             lai_current       = min(leaf_inc, currentCohort%treelai - lai_layers_above)
-             cumulative_lai    = lai_canopy_above + lai_layers_above + 0.5*lai_current
-
-             if (currentCohort%year_net_uptake(z) /= 999._r8)then !there was activity this year in this leaf layer.
-
+             cumulative_lai    = lai_canopy_above + cumulative_lai_cohort
+
+             ! There was activity this year in this leaf layer.  This should only occur for bottom most leaf layer
+             if (currentCohort%year_net_uptake(z) /= 999._r8)then 
 
                 ! Calculate sla_levleaf following the sla profile with overlying leaf area
                 ! Scale for leaf nitrogen profile
@@ -502,42 +563,125 @@ subroutine trim_canopy( currentSite )
                    currentCohort%leaf_cost = currentCohort%leaf_cost * &
                          (EDPftvarcon_inst%grperc(ipft) + 1._r8)
                 endif
-                if (currentCohort%year_net_uptake(z) < currentCohort%leaf_cost)then
-                   if (currentCohort%canopy_trim > EDPftvarcon_inst%trim_limit(ipft))then
 
-                      if ( debug ) then
-                         write(fates_log(),*) 'trimming leaves', &
-                               currentCohort%canopy_trim,currentCohort%leaf_cost
-                      endif
+                ! Construct the arrays for a least square fit of the net_net_uptake versus the cumulative lai
+                ! if at least nll leaf layers are present in the current cohort and only for the bottom nll 
+                ! leaf layers.
+                if (currentCohort%nv > nll .and. currentCohort%nv - z < nll) then 
+
+                  ! Build the A matrix for the LHS of the linear system. A = [n  sum(x); sum(x)  sum(x^2)]
+                  ! where n = nll and x = yearly_net_uptake-leafcost
+                  nnu_clai_a(1,1) = nnu_clai_a(1,1) + 1 ! Increment for each layer used
+                  nnu_clai_a(1,2) = nnu_clai_a(1,2) + currentCohort%year_net_uptake(z) - currentCohort%leaf_cost
+                  nnu_clai_a(2,1) = nnu_clai_a(1,2)
+                  nnu_clai_a(2,2) = nnu_clai_a(2,2) + (currentCohort%year_net_uptake(z) - currentCohort%leaf_cost)**2
+
+                  ! Build the B matrix for the RHS of the linear system. B = [sum(y); sum(x*y)]
+                  ! where x = yearly_net_uptake-leafcost and y = cumulative_lai_cohort
+                  nnu_clai_b(1,1) = nnu_clai_b(1,1) + cumulative_lai_cohort
+                  nnu_clai_b(2,1) = nnu_clai_b(2,1) + (cumulative_lai_cohort * & 
+                                    (currentCohort%year_net_uptake(z) - currentCohort%leaf_cost))
+                end if
+
+                ! Check leaf cost against the yearly net uptake for that cohort leaf layer
+                if (currentCohort%year_net_uptake(z) < currentCohort%leaf_cost) then
+                   ! Make sure the cohort trim fraction is great than the pft trim limit
+                   if (currentCohort%canopy_trim > EDPftvarcon_inst%trim_limit(ipft)) then
+
+                     !  if ( debug ) then
+                     !     write(fates_log(),*) 'trimming leaves', &
+                     !           currentCohort%canopy_trim,currentCohort%leaf_cost
+                     !  endif
 
                       ! keep trimming until none of the canopy is in negative carbon balance.              
-                      if (currentCohort%hite > EDPftvarcon_inst%hgt_min(ipft))then
+                      if (currentCohort%hite > EDPftvarcon_inst%hgt_min(ipft)) then
                          currentCohort%canopy_trim = currentCohort%canopy_trim - &
                                EDPftvarcon_inst%trim_inc(ipft)
-                         if (EDPftvarcon_inst%evergreen(ipft) /= 1)then
+
+                         if (EDPftvarcon_inst%evergreen(ipft) /= 1) then
                             currentCohort%laimemory = currentCohort%laimemory * &
                                   (1.0_r8 - EDPftvarcon_inst%trim_inc(ipft)) 
                          endif
+
                          trimmed = .true.
-                      endif
-                   endif
-                endif
-             endif !leaf activity? 
-          enddo !z
 
+                      endif ! hite check
+                   endif ! trim limit check
+                endif ! net uptake check
+             endif ! leaf activity check 
+          enddo ! z, leaf layer loop
+
+          ! Compute the optimal cumulative lai based on the cohort net-net uptake profile if at least 2 leaf layers
+          if (nnu_clai_a(1,1) > 1) then
+
+            ! Compute the optimum size of the work array
+            lwork = -1 ! Ask sgels to compute optimal number of entries for work
+            call dgels(trans, m, n, nrhs, nnu_clai_a, lda, nnu_clai_b, ldb, work, lwork, info)
+            lwork = int(work(1)) ! Pick the optimum.  TBD, can work(1) come back with greater than work size?
+
+            ! if (debug) then
+            !    write(fates_log(),*) 'LLSF lwork output (info, lwork):', info, lwork
+            ! endif
+
+            ! Compute the minimum of 2-norm of of the least squares fit to solve for X
+            ! Note that dgels returns the solution by overwriting the nnu_clai_b array.  
+            ! The result has the form: X = [b; m]
+            ! where b = y-intercept (i.e. the cohort lai that has zero yearly net-net uptake)
+            ! and m is the slope of the linear fit
+            call dgels(trans, m, n, nrhs, nnu_clai_a, lda, nnu_clai_b, ldb, work, lwork, info)
+
+            if (info < 0) then
+               write(fates_log(),*) 'LLSF optimium LAI calculation returned illegal value'
+               call endrun(msg=errMsg(sourcefile, __LINE__))
+            endif
+
+            if (debug) then
+               write(fates_log(),*) 'LLSF optimium LAI (intercept,slope):', nnu_clai_b
+               write(fates_log(),*) 'LLSF optimium LAI:', nnu_clai_b(1,1)
+               write(fates_log(),*) 'LLSF optimium LAI info:', info
+               write(fates_log(),*) 'LAI fraction (optimum_lai/cumulative_lai):', nnu_clai_b(1,1) / cumulative_lai_cohort
+            endif
+
+            ! Calculate the optimum trim based on the initial canopy trim value
+            if (cumulative_lai_cohort > 0._r8) then  ! Sometime cumulative_lai comes in at 0.0?
+
+               ! 
+               optimum_trim = (nnu_clai_b(1,1) / cumulative_lai_cohort) * initial_trim
+               optimum_laimem = (nnu_clai_b(1,1) / cumulative_lai_cohort) * initial_laimem
+
+               ! Determine if the optimum trim value makes sense.  The smallest cohorts tend to have unrealistic fits.
+               if (optimum_trim > 0. .and. optimum_trim < 1.) then
+                  currentCohort%canopy_trim = optimum_trim
+
+                  ! If the cohort pft is not evergreen we reduce the laimemory as well
+                  if (EDPftvarcon_inst%evergreen(ipft) /= 1) then
+                     currentCohort%laimemory = optimum_laimem
+                  endif
+
+                  trimmed = .true.
+
+               endif
+            endif
+         endif
+
+          ! Reset activity for the cohort for the start of the next year
           currentCohort%year_net_uptake(:) = 999.0_r8
+
+          ! Add to trim fraction if cohort not trimmed at all 
           if ( (.not.trimmed) .and.currentCohort%canopy_trim < 1.0_r8)then
              currentCohort%canopy_trim = currentCohort%canopy_trim + EDPftvarcon_inst%trim_inc(ipft)
           endif 
 
           if ( debug ) then
-             write(fates_log(),*) 'trimming',currentCohort%canopy_trim
+             write(fates_log(),*) 'trimming:',currentCohort%canopy_trim
           endif
 
           ! currentCohort%canopy_trim = 1.0_r8 !FIX(RF,032414) this turns off ctrim for now. 
           currentCohort => currentCohort%shorter
+          icohort = icohort + 1
        enddo
        currentPatch => currentPatch%older
+       ipatch = ipatch + 1
     enddo
 
   end subroutine trim_canopy

biogeochem/FatesAllometryMod.F90

@@ -1333,6 +1333,9 @@ subroutine d2h_chave2014(d,p1,p2,p3,dbh_maxh,h,dhdd)
     ! Chave et al. Improved allometric models to estimate the abovegroud
     ! biomass of tropical trees.  Global Change Biology. V20, p3177-3190. 2015.
     !
+    ! p1 =  0.893 - E
+    ! p2 =  0.76
+    ! p3 = -0.034
     ! =========================================================================
 
     real(r8),intent(in)  :: d        ! plant diameter [cm]
@@ -1565,6 +1568,9 @@ subroutine dh2bagw_chave2014(d,h,dhdd,p1,p2,wood_density,c2b,bagw,dbagwdd)
     ! Output:
     ! bagw:   Total above ground biomass [kgC]
     !
+    ! Chave's Paper has p1 = 0.0673, p2 = 0.976
+    !
+
     ! =========================================================================
 
 

biogeophys/FatesHydroWTFMod.F90

@@ -210,6 +210,7 @@ function th_from_psi_base(this,psi) result(th)
     class(wrf_type)     :: this
     real(r8),intent(in) :: psi
     real(r8)            :: th
+    th = 0._r8
     write(fates_log(),*) 'The base water retention function'
     write(fates_log(),*) 'should never be actualized'
     write(fates_log(),*) 'check how the class pointer was setup'
@@ -219,6 +220,7 @@ function psi_from_th_base(this,th) result(psi)
     class(wrf_type)     :: this
     real(r8),intent(in) :: th
     real(r8)            :: psi
+    psi = 0._r8
     write(fates_log(),*) 'The base water retention function'
     write(fates_log(),*) 'should never be actualized'
     write(fates_log(),*) 'check how the class pointer was setup'
@@ -228,6 +230,7 @@ function dpsidth_from_th_base(this,th) result(dpsidth)
     class(wrf_type)     :: this
     real(r8),intent(in) :: th
     real(r8)            :: dpsidth
+    dpsidth = 0._r8
     write(fates_log(),*) 'The base water retention function'
     write(fates_log(),*) 'should never be actualized'
     write(fates_log(),*) 'check how the class pointer was setup'
@@ -237,6 +240,7 @@ function ftc_from_psi_base(this,psi) result(ftc)
     class(wkf_type)     :: this
     real(r8),intent(in) :: psi
     real(r8)            :: ftc
+    ftc = 0._r8
     write(fates_log(),*) 'The base water retention function'
     write(fates_log(),*) 'should never be actualized'
     write(fates_log(),*) 'check how the class pointer was setup'
@@ -246,6 +250,7 @@ function dftcdpsi_from_psi_base(this,psi) result(dftcdpsi)
     class(wkf_type)     :: this
     real(r8),intent(in) :: psi
     real(r8)            :: dftcdpsi
+    dftcdpsi = 0._r8
     write(fates_log(),*) 'The base water retention function'
     write(fates_log(),*) 'should never be actualized'
     write(fates_log(),*) 'check how the class pointer was setup'

biogeophys/FatesPlantHydraulicsMod.F90

@@ -3604,9 +3604,9 @@ subroutine ImTaylorSolve1D(site_hydr,cohort,cohort_hydr,dtime,q_top, &
         end if
 
         ! Save the number of times we refined our sub-step counts (iterh1)
-        cohort_hydr%iterh1 = max(cohort_hydr%iterh1,real(iter))
+        cohort_hydr%iterh1 = max(cohort_hydr%iterh1,real(iter,r8))
         ! Save the number of sub-steps we ultimately used
-        cohort_hydr%iterh2 = max(cohort_hydr%iterh2,real(nsteps))
+        cohort_hydr%iterh2 = max(cohort_hydr%iterh2,real(nsteps,r8))
 
         ! Update water contents in the relevant plant compartments [m3/m3]
         ! -------------------------------------------------------------------------------
@@ -4967,7 +4967,7 @@ subroutine MatSolve2D(site_hydr,cohort,cohort_hydr, &
          cohort_hydr%iterh1 = cohort_hydr%iterh1 + 1
 
          ! Save the  max number of Newton iterations needed
-         cohort_hydr%iterh2 = max(cohort_hydr%iterh2,real(nwtn_iter))
+         cohort_hydr%iterh2 = max(cohort_hydr%iterh2,real(nwtn_iter,r8))
 
          print*,'Completed a newton solve'
          print*,psi_node(:)