Refactoring: use function do generate derived variables rather using …

…variables in type.

Note: the output changed slightly as the value of the previous year for BA (BA_ys) was getting a value which was not updated after the annual diagnostic (while DBH was). In other words, both _ys values are now coherent

src/datatypes_biomee.mod.f90

@@ -93,9 +93,6 @@ module datatypes_biomee
     real :: height    = 0.0              ! vegetation height, m
     real :: crownarea = 1.0              ! crown area, m2 tree-1
     real :: leafarea  = 0.0              ! total area of leaves, m2 tree-1
-    real :: lai       = 0.0              ! crown leaf area index, m2/m2
-    real :: BA        = 0.0              ! tree basal area, m2 tree-1
-    real :: Volume    = 0.0              ! tree basal volume, m3 tree-1
 
     !===== Biological prognostic variables
     integer :: species    = 1            ! vegetation species
@@ -156,8 +153,8 @@ module datatypes_biomee
     real    :: resg               = 0.0           ! growth respiration, kg C day-1 tree-1
 
     !===== Memory variables used for computing deltas
-    real    :: DBH_ys            = dummy          ! DBH at the begining of a year (growing season)
-    real    :: BA_ys             = dummy          ! Basal area at the beginning og a year
+    real    :: DBH_ys            = 0.0            ! DBH at the begining of a year (growing season)
+    real    :: BA_ys             = 0.0            ! Basal area at the beginning og a year
 
   end type cohort_type
 
@@ -277,6 +274,29 @@ module datatypes_biomee
   end type vegn_tile_type
 
 contains
+  function lai(cohort) result(res)
+    ! Leaf area index: surface of leaves per m2 of crown
+    real :: res
+    type(cohort_type) :: cohort
+
+    res = cohort%leafarea / cohort%crownarea !(cohort%crownarea *(1.0-sp%internal_gap_frac))
+  end function lai
+
+  function basal_area(cohort) result(ba)
+    ! Tree basal area, m2 tree-1
+    real :: ba
+    type(cohort_type) :: cohort
+
+    ba = pi/4 * cohort%dbh * cohort%dbh
+  end function basal_area
+
+  function volume(cohort) result(vol)
+    ! Tree basal volume, m3 tree-1
+    real :: vol
+    type(cohort_type) :: cohort
+
+    vol = (cohort%psapw%c%c12 + cohort%pwood%c%c12) / myinterface%params_species(cohort%species)%rho_wood
+  end function volume
 
   subroutine update_fluxes(fluxes, delta)
     ! Add delta quantities to partial fluxes (accounting)
@@ -299,7 +319,6 @@ subroutine Zero_diagnostics(vegn)
     type(vegn_tile_type), intent(inout) :: vegn
 
     ! local variables
-    type(common_fluxes) :: cf
     type(cohort_type), pointer :: cc
     integer :: i
 
@@ -326,6 +345,10 @@ subroutine Zero_diagnostics(vegn)
     do i = 1, vegn%n_cohorts
       cc => vegn%cohorts(i)
 
+      ! Save last year's values
+      cc%DBH_ys       = cc%dbh
+      cc%BA_ys        = basal_area(cc)
+
       cc%C_growth     = 0.0
       cc%N_growth     = 0.0
       cc%gpp          = 0.0
@@ -336,16 +359,14 @@ subroutine Zero_diagnostics(vegn)
       cc%resg         = 0.0
       cc%transp       = 0.0
 
-      ! daily
-      cc%daily_fluxes = cf
+      ! Reset daily
+      cc%daily_fluxes = common_fluxes()
 
       ! annual
-      cc%annual_fluxes = cf
+      cc%annual_fluxes = common_fluxes()
       cc%NPPleaf      = 0.0
       cc%NPProot      = 0.0
       cc%NPPwood      = 0.0
-      cc%DBH_ys       = cc%dbh
-      cc%BA_ys        = cc%BA
       cc%n_deadtrees  = 0.0
       cc%c_deadtrees  = 0.0
       cc%m_turnover   = 0.0
@@ -426,15 +447,9 @@ subroutine summarize_tile( vegn )
         vegn%DBH12pow2  = vegn%DBH12pow2 + cc%dbh      * cc%dbh     * cc%nindivs
       endif
 
-      if (cc%age    > vegn%MaxAge)       vegn%MaxAge    = cc%age
-      if (cc%Volume > vegn%MaxVolume)    vegn%MaxVolume = cc%Volume ! maxloc(cc%age)
-      if (cc%dbh    > vegn%MaxDBH)       vegn%MaxDBH    = cc%dbh    ! maxloc(cc%age)
-
-      ! vegn%NPPL      = vegn%NPPL   + fleaf * cc%nindivs
-      ! vegn%NPPW      = vegn%NPPW   + fwood * cc%nindivs
-
-      ! vegn%n_deadtrees = vegn%n_deadtrees + cc%n_deadtrees
-      ! vegn%c_deadtrees = vegn%c_deadtrees + cc%c_deadtrees
+      vegn%MaxAge    = MAX(cc%age, vegn%MaxAge)
+      vegn%MaxVolume = MAX(volume(cc), vegn%MaxVolume)
+      vegn%MaxDBH    = MAX(cc%dbh, vegn%MaxDBH)
 
     enddo
 
@@ -513,7 +528,6 @@ subroutine daily_diagnostics( vegn , iyears, idoy, state, out_daily_tile )
 
     ! local variables
     type(cohort_type), pointer :: cc    ! current cohort
-    type(common_fluxes) :: cf
     integer :: i
 
     ! cohorts output
@@ -523,8 +537,8 @@ subroutine daily_diagnostics( vegn , iyears, idoy, state, out_daily_tile )
       ! running annual sum
       call update_fluxes(cc%annual_fluxes, cc%daily_fluxes)
 
-      ! Zero Daily variables
-      cc%daily_fluxes = cf
+      ! Reset daily variables
+      cc%daily_fluxes = common_fluxes()
     enddo
 
     ! Tile level, daily
@@ -599,7 +613,7 @@ subroutine annual_diagnostics(vegn, iyears, out_annual_cohorts, out_annual_tile)
 
     ! local variables
     type(cohort_type), pointer :: cc
-    real :: treeG, fseed, fleaf=0, froot, fwood=0, dDBH, dBA
+    real :: treeG, fseed, fleaf=0, froot, fwood=0, dDBH, BA, dBA
     real :: plantC, plantN, soilC, soilN
     integer :: i
 
@@ -648,9 +662,8 @@ subroutine annual_diagnostics(vegn, iyears, out_annual_cohorts, out_annual_tile)
       froot     = cc%NPProot / treeG
       fwood     = cc%NPPwood / treeG
       dDBH      = cc%dbh - cc%DBH_ys !in m
-      cc%BA     = pi/4 * cc%dbh * cc%dbh
-      dBA       = cc%BA - cc%BA_ys
-      cc%Volume = (cc%psapw%c%c12 + cc%pwood%c%c12) / myinterface%params_species(cc%species)%rho_wood
+      BA        = basal_area(cc)
+      dBA       = BA - cc%BA_ys
 
       out_annual_cohorts(i)%year        = iyears
       out_annual_cohorts(i)%cID         = cc%ccID
@@ -662,7 +675,7 @@ subroutine annual_diagnostics(vegn, iyears, out_annual_cohorts, out_annual_tile)
       out_annual_cohorts(i)%dDBH        = dDBH * 100     ! *100 to convert m in cm
       out_annual_cohorts(i)%height      = cc%height
       out_annual_cohorts(i)%age         = cc%age
-      out_annual_cohorts(i)%BA          = cc%BA
+      out_annual_cohorts(i)%BA          = BA
       out_annual_cohorts(i)%dBA         = dBA
       out_annual_cohorts(i)%Acrown      = cc%crownarea
       out_annual_cohorts(i)%Aleaf       = cc%leafarea

src/gpp_biomee.mod.f90

@@ -163,7 +163,7 @@ subroutine gpp( forcing, vegn, first_simu_step )
           fw = 0.0
           fs = 0.0
 
-          call gs_leuning(rad_top, rad_net, TairK, cana_q, cc%lai, &
+          call gs_leuning(rad_top, rad_net, TairK, cana_q, lai(cc), &
             p_surf, water_supply, cc%species, sp%pt, &
             cana_co2, extinct, fs+fw, &
             psyn, resp, w_scale2, transp )

src/vegetation_biomee.mod.f90

@@ -365,8 +365,6 @@ subroutine vegn_growth_EW( vegn )
         cc%height    = cc%height    + dHeight
         cc%crownarea = cc%crownarea + dCA
         cc%leafarea  = leaf_area_from_biomass(cc%pleaf%c%c12, cc%species)
-        ! cc%lai is surface of leaves per m2 of crown
-        cc%lai       = cc%leafarea/cc%crownarea !(cc%crownarea *(1.0-sp%internal_gap_frac))
         ! vegn%LAI is the surface of leaves per m2 of ground/tile
         vegn%LAI     = vegn%LAI + cc%leafarea * cc%nindivs
 
@@ -462,8 +460,6 @@ subroutine vegn_phenology( vegn )
     ! local variables
     type(cohort_type), pointer :: cc
     integer :: i
-    ! real    :: grassdensity   ! for grasses only
-    ! real    :: BL_u,BL_c
     real    :: ccNSC, ccNSN
     logical :: cc_firstday = .false.
     logical :: TURN_ON_life = .false., TURN_OFF_life
@@ -482,14 +478,14 @@ subroutine vegn_phenology( vegn )
       associate (sp => myinterface%params_species(cc%species) )
 
       ! for evergreen
-      if (sp%phenotype==1 .and. cc%status==LEAF_OFF) cc%status=LEAF_ON
+      if (sp%phenotype == 1 .and. cc%status == LEAF_OFF) cc%status=LEAF_ON
 
       ! for deciduous and grasses
       TURN_ON_life = (sp%phenotype == 0 .and. &
         cc%status    == LEAF_OFF       .and. &
         cc%gdd        > sp%gdd_crit    .and. &
         vegn%tc_pheno > sp%tc_crit_on) .and. &
-        (sp%lifeform .ne. 0 .OR.(sp%lifeform .eq. 0 .and. cc%layer==1))
+        (sp%lifeform /= 0 .OR.(sp%lifeform == 0 .and. cc%layer == 1))
 
       cc_firstday = .false.
       if (TURN_ON_life) then
@@ -1277,8 +1273,7 @@ subroutine update_plant_pools( cc, vegn, dBL, dBR, dBStem, dNL, dNR, dNStem)
       cc%proot%n%n14 = cc%proot%n%n14 - dNR
 
       ! update leaf area and LAI
-      cc%leafarea= leaf_area_from_biomass(cc%pleaf%c%c12, cc%species)
-      cc%lai     = cc%leafarea / (cc%crownarea *(1.0-sp%internal_gap_frac))
+      cc%leafarea = leaf_area_from_biomass(cc%pleaf%c%c12, cc%species)
 
       ! update NPP for leaves, fine roots, and wood
       cc%NPPleaf = cc%NPPleaf - myinterface%params_tile%l_fract * dBL
@@ -1842,10 +1837,6 @@ subroutine merge_cohorts(c1, c2)
 
       x1 = c1%nindivs/(c1%nindivs+c2%nindivs)
       x2 = c2%nindivs/(c1%nindivs+c2%nindivs)
-      !else
-      !   x1 = 0.5
-      !   x2 = 0.5
-      !endif
 
       ! update number of individuals in merged cohort
       c2%nindivs = c1%nindivs + c2%nindivs
@@ -1858,25 +1849,20 @@ subroutine merge_cohorts(c1, c2)
       c2%pseed%c%c12 = x1*c1%pseed%c%c12 + x2*c2%pseed%c%c12
       c2%plabl%c%c12 = x1*c1%plabl%c%c12 + x2*c2%plabl%c%c12
 
-      ! Allometry
-      c2%dbh = x1*c1%dbh + x2*c2%dbh
-      c2%height = x1*c1%height + x2*c2%height
-      c2%crownarea = x1*c1%crownarea + x2*c2%crownarea
-      c2%age = x1*c1%age + x2*c2%age
-      c2%C_growth = x1*c1%C_growth + x2*c2%C_growth
-      c2%topyear = x1*c1%topyear + x2*c2%topyear
-
       ! Nitrogen
       c2%pleaf%n%n14 = x1*c1%pleaf%n%n14 + x2*c2%pleaf%n%n14
       c2%proot%n%n14 = x1*c1%proot%n%n14 + x2*c2%proot%n%n14
       c2%psapw%n%n14 = x1*c1%psapw%n%n14 + x2*c2%psapw%n%n14
       c2%pwood%n%n14 = x1*c1%pwood%n%n14 + x2*c2%pwood%n%n14
       c2%pseed%n%n14 = x1*c1%pseed%n%n14 + x2*c2%pseed%n%n14
       c2%plabl%n%n14 = x1*c1%plabl%n%n14 + x2*c2%plabl%n%n14
-
-      !  calculate the resulting dry heat capacity
-      c2%leafarea = leaf_area_from_biomass(c2%pleaf%c%c12, c2%species)
-      call init_cohort_allometry(c2) !Enseng comments
+
+      ! Cohort age
+      c2%age = x1*c1%age + x2*c2%age
+      c2%topyear = x1*c1%topyear + x2*c2%topyear
+
+      ! Allometry
+      call init_cohort_allometry(c2)
     endif
 
   end subroutine merge_cohorts
@@ -1956,6 +1942,8 @@ subroutine init_cohort_allometry( cc )
     integer :: layer
     real    :: btot ! total biomass per individual, kg C
 
+    cc%leafarea = leaf_area_from_biomass(cc%pleaf%c%c12, cc%species)
+
     associate(sp=>myinterface%params_species(cc%species))
       btot = max(0.0001, cc%pwood%c%c12 + cc%psapw%c%c12)
       layer = max(1, cc%layer)
@@ -1983,14 +1971,8 @@ function leaf_area_from_biomass(bl,species) result (area)
     real :: area ! returned value
     real,    intent(in) :: bl      ! biomass of leaves, kg C/individual
     integer, intent(in) :: species ! species
-    ! integer, intent(in) :: layer, firstlayer
-    ! modified by Weng (2014-01-09), 07-18-2017
+
     area = bl/myinterface%params_species(species)%LMA
-    !if (layer > 1.AND. firstlayer == 0) then
-    !   area = bl/(0.5*myinterface%params_species(species)%LMA) ! half thickness for leaves in understory
-    !else
-    !   area = bl/myinterface%params_species(species)%LMA
-    !endif
 
   end function