Merge branch 'master' into andre-ed-v0.1.0-buggy-science

Merge the ed-v0.1.0 buggy science branch up to master. Many runs die
at runtime (see error below). I believe that this is the error Rosie
and Charlie are expecting, so I am handing the branch off to them.

Test suite: ed on yellowstone with intel, pgi, gnu
Test baseline: n/a
Test namelist changes: none
Test status: answer changing?

Some tests run to completion, but there are many runtime failures. Example:

    SMS_D_Ld5.f10_f10.ICLM45ED.yellowstone_intel.clm-edTest

 14:forrtl: severe (408): fort: (2): Subscript #1 of the array FLOATING_CN_RATIO_DECOMP_POOLS has value 8 which is greater than the upper bound of 7
  14:
  14:Image              PC                Routine            Line        Source
  14:cesm.exe           0000000003DF3550  Unknown               Unknown  Unknow
n
  14:cesm.exe           0000000001D0CD75  soilbiogeochemdec         429  SoilBi
ogeochemDecompCascadeCNMod.F90

biogeochem/EDCohortDynamicsMod.F90

@@ -28,6 +28,8 @@ module EDCohortDynamicsMod
   public :: countCohorts
   public :: allocate_live_biomass
 
+  logical, parameter :: DEBUG  = .false. ! local debug flag
+
   ! 10/30/09: Created by Rosie Fisher
   !-------------------------------------------------------------------------------------!
 
@@ -62,7 +64,7 @@ subroutine create_cohort(patchptr, pft, nn, hite, dbh, &
     type(ed_cohort_type), pointer :: storebigcohort   
     integer :: tnull,snull                      ! are the tallest and shortest cohorts allocate
     !----------------------------------------------------------------------
- 
+
     allocate(new_cohort)
     udata%cohort_number = udata%cohort_number + 1 !give each cohort a unique number for checking cohort fusing routine.
 
@@ -88,15 +90,20 @@ subroutine create_cohort(patchptr, pft, nn, hite, dbh, &
     new_cohort%balive       = balive
     new_cohort%bstore       = bstore
 
+    if ( DEBUG ) write(iulog,*) 'EDCohortDyn I ',bstore
+
     if (new_cohort%dbh <= 0.0_r8 .or. new_cohort%n == 0._r8 .or. new_cohort%pft == 0 &
           .or. new_cohort%canopy_trim <= 0.0_r8 .or. new_cohort%balive <= 0._r8) then
-       write(iulog,*) 'ED: something is zero in create_cohort',new_cohort%indexnumber,new_cohort%dbh,new_cohort%n, &
-       new_cohort%pft,new_cohort%canopy_trim,new_cohort%balive
+             write(iulog,*) 'ED: something is zero in create_cohort', &
+                             new_cohort%indexnumber,new_cohort%dbh,new_cohort%n, &
+                             new_cohort%pft,new_cohort%canopy_trim,new_cohort%balive
     endif
-    if (new_cohort%siteptr%status==2.and.pftcon%season_decid(pft) == 1) then
+
+    if (new_cohort%siteptr%status==2 .and. pftcon%season_decid(pft) == 1) then
       new_cohort%laimemory = 0.0_r8
     endif
-    if (new_cohort%siteptr%dstatus==2.and.pftcon%stress_decid(pft) == 1) then
+
+    if (new_cohort%siteptr%dstatus==2 .and. pftcon%stress_decid(pft) == 1) then
       new_cohort%laimemory = 0.0_r8
     endif
 
@@ -187,12 +194,11 @@ subroutine allocate_live_biomass(cc_p)
     endif
 
     if (leaves_off_switch==1) then
-
-    !the purpose of this section is to figure out the root and stem biomass when the leaves are off
-    !at this point, we know the former leaf mass (laimemory) and the current alive mass
-    !because balive may decline in the off-season, we need to adjust the root and stem biomass that are predicted
-    !from the laimemory, for the fact that we now might not have enough live biomass to support the hypothesized root mass
-    !thus, we use 'ratio_balive' to adjust br and bsw. Apologies that this is so complicated! RF
+       !the purpose of this section is to figure out the root and stem biomass when the leaves are off
+       !at this point, we know the former leaf mass (laimemory) and the current alive mass
+       !because balive may decline in the off-season, we need to adjust the root and stem biomass that are predicted
+       !from the laimemory, for the fact that we now might not have enough live biomass to support the hypothesized root mass
+       !thus, we use 'ratio_balive' to adjust br and bsw. Apologies that this is so complicated! RF
        currentcohort%bl  = 0.0_r8
        ideal_balive      = currentcohort%laimemory * pftcon%froot_leaf(ft) +  &
             currentcohort%laimemory*  EDecophyscon%sapwood_ratio(ft) * currentcohort%hite
@@ -203,12 +209,13 @@ subroutine allocate_live_biomass(cc_p)
        ratio_balive           = currentcohort%balive / ideal_balive
        currentcohort%br       = currentcohort%br  * ratio_balive
        currentcohort%bsw      = currentcohort%bsw * ratio_balive               
-    endif
 
+    endif
 
     if (abs(currentcohort%balive -currentcohort%bl- currentcohort%br - currentcohort%bsw)>1e-12) then
-       write(iulog,*) 'issue with carbon allocation in create_cohort',&
-       currentcohort%balive -currentcohort%bl- currentcohort%br - currentcohort%bsw, currentcohort%status_coh,currentcohort%balive
+       write(iulog,*) 'issue with carbon allocation in create_cohort,allocate_live_biomass',&
+            currentcohort%balive -currentcohort%bl- currentcohort%br - currentcohort%bsw, &
+            currentcohort%status_coh,currentcohort%balive
        write(iulog,*) 'actual vs predicted balive',ideal_balive,currentcohort%balive ,ratio_balive,leaf_frac
        write(iulog,*) 'leaf,root,stem',currentcohort%bl,currentcohort%br,currentcohort%bsw
        write(iulog,*) 'pft',ft,pftcon%evergreen(ft),pftcon%season_decid(ft),leaves_off_switch
@@ -414,29 +421,44 @@ subroutine terminate_cohorts( patchptr )
        if (currentCohort%n/currentPatch%area <= 0.00001_r8 .or. currentCohort%dbh <  &
             0.00001_r8.and.currentCohort%bstore < 0._r8) then
           terminate = 1
-          !  write(iulog,*) 'terminating cohorts 1',currentCohort%n/currentPatch%area,currentCohort%dbh
+
+          if ( DEBUG ) then
+             write(iulog,*) 'terminating cohorts 1',currentCohort%n/currentPatch%area,currentCohort%dbh
+          endif
+
        endif
 
        ! In the third canopy layer
        if (currentCohort%canopy_layer > NCLMAX) then 
           terminate = 1
-          ! write(iulog,*) 'terminating cohorts 2', currentCohort%canopy_layer
+
+          if ( DEBUG ) then
+            write(iulog,*) 'terminating cohorts 2', currentCohort%canopy_layer
+          endif
+
        endif
 
        ! live biomass pools are terminally depleted
        if (currentCohort%balive < 1e-10_r8 .or. currentCohort%bstore < 1e-10_r8) then 
           terminate = 1  
-          ! write(iulog,*) 'terminating cohorts 3', currentCohort%balive,currentCohort%bstore
+
+          if ( DEBUG ) then
+            write(iulog,*) 'terminating cohorts 3', currentCohort%balive,currentCohort%bstore
+          endif
+
        endif
 
        ! Total cohort biomass is negative
        if (currentCohort%balive+currentCohort%bdead+currentCohort%bstore < 0._r8) then
           terminate = 1
-          ! write(iulog,*) 'terminating cohorts 4', currentCohort%balive, currentCohort%bstore, currentCohort%bdead, &
-          ! currentCohort%balive+currentCohort%bdead+&
-          ! currentCohort%bstore, currentCohort%n
-       endif
 
+          if ( DEBUG ) then
+            write(iulog,*) 'terminating cohorts 4', currentCohort%balive, currentCohort%bstore, currentCohort%bdead, &
+                           currentCohort%balive+currentCohort%bdead+&
+                           currentCohort%bstore, currentCohort%n
+          endif
+
+       endif
 
        if (terminate == 1) then 
           if (.not. associated(currentCohort%taller)) then
@@ -537,13 +559,20 @@ subroutine fuse_cohorts(patchptr)
                    if (currentCohort%pft == nextc%pft) then              
 
                       ! check cohorts in same c. layer. before fusing
-                      if (currentCohort%canopy_layer == nextc%canopy_layer) then 
+                      if (currentCohort%canopy_layer == nextc%canopy_layer) then
+
                          fusion_took_place = 1         
                          newn = currentCohort%n + nextc%n    ! sum individuals in both cohorts.     
 
                          currentCohort%balive    = (currentCohort%n*currentCohort%balive    + nextc%n*nextc%balive)/newn
                          currentCohort%bdead     = (currentCohort%n*currentCohort%bdead     + nextc%n*nextc%bdead)/newn
+
+                         if ( DEBUG ) write(iulog,*) 'EDcohortDyn I ',currentCohort%bstore
+
                          currentCohort%bstore    = (currentCohort%n*currentCohort%bstore    + nextc%n*nextc%bstore)/newn   
+
+                         if ( DEBUG ) write(iulog,*) 'EDcohortDyn II ',currentCohort%bstore
+
                          currentCohort%seed_prod = (currentCohort%n*currentCohort%seed_prod + nextc%n*nextc%seed_prod)/newn  
                          currentCohort%root_md   = (currentCohort%n*currentCohort%root_md   + nextc%n*nextc%root_md)/newn   
                          currentCohort%leaf_md   = (currentCohort%n*currentCohort%leaf_md   + nextc%n*nextc%leaf_md)/newn  
@@ -552,18 +581,34 @@ subroutine fuse_cohorts(patchptr)
 
                          currentCohort%carbon_balance = (currentCohort%n*currentCohort%carbon_balance + &
                               nextc%n*nextc%carbon_balance)/newn
+
                          currentCohort%storage_flux = (currentCohort%n*currentCohort%storage_flux + &
                               nextc%n*nextc%storage_flux)/newn               
 
                          currentCohort%b           = (currentCohort%n*currentCohort%b           + nextc%n*nextc%b)/newn
                          currentCohort%bsw         = (currentCohort%n*currentCohort%bsw         + nextc%n*nextc%bsw)/newn
                          currentCohort%bl          = (currentCohort%n*currentCohort%bl          + nextc%n*nextc%bl)/newn
+
+                         if ( DEBUG ) write(iulog,*) 'EDcohortDyn 569 ',currentCohort%br
+                         if ( DEBUG ) write(iulog,*) 'EDcohortDyn 570 ',currentCohort%n
+                         if ( DEBUG ) write(iulog,*) 'EDcohortDyn 571 ',nextc%br
+                         if ( DEBUG ) write(iulog,*) 'EDcohortDyn 572 ',nextc%n
+
                          currentCohort%br          = (currentCohort%n*currentCohort%br          + nextc%n*nextc%br)/newn
                          currentCohort%hite        = (currentCohort%n*currentCohort%hite        + nextc%n*nextc%hite)/newn         
                          currentCohort%dbh         = (currentCohort%n*currentCohort%dbh         + nextc%n*nextc%dbh)/newn
+
                          currentCohort%gpp_acc     = (currentCohort%n*currentCohort%gpp_acc     + nextc%n*nextc%gpp_acc)/newn
+
+                         if ( DEBUG ) write(iulog,*) 'EDcohortDyn III ',currentCohort%npp_acc
+                         if ( DEBUG ) write(iulog,*) 'EDcohortDyn IV ',currentCohort%resp_acc
+
                          currentCohort%npp_acc     = (currentCohort%n*currentCohort%npp_acc     + nextc%n*nextc%npp_acc)/newn
                          currentCohort%resp_acc    = (currentCohort%n*currentCohort%resp_acc    + nextc%n*nextc%resp_acc)/newn
+
+                         if ( DEBUG ) write(iulog,*) 'EDcohortDyn V ',currentCohort%npp_acc
+                         if ( DEBUG ) write(iulog,*) 'EDcohortDyn VI ',currentCohort%resp_acc
+
                          currentCohort%resp        = (currentCohort%n*currentCohort%resp        + nextc%n*nextc%resp)/newn
                          currentCohort%npp         = (currentCohort%n*currentCohort%npp         + nextc%n*nextc%npp)/newn
                          currentCohort%gpp         = (currentCohort%n*currentCohort%gpp         + nextc%n*nextc%gpp)/newn
@@ -589,9 +634,11 @@ subroutine fuse_cohorts(patchptr)
                          else
                             nextnextc%taller => nextc%taller
                          endif
+
                          if (associated(nextc)) then       
                             deallocate(nextc)            
                          endif
+
                       endif !canopy layer
                    endif !pft
                 endif  !index no. 
@@ -602,6 +649,7 @@ subroutine fuse_cohorts(patchptr)
              else
                 nextc => nextnextc !if we have removed next
              endif
+
           enddo !end checking nextc cohort loop
 
           if (associated (currentCohort%shorter)) then
@@ -621,11 +669,13 @@ subroutine fuse_cohorts(patchptr)
 
        if (nocohorts > maxcohorts) then
           iterate = 1
-          dynamic_fusion_tolerance = dynamic_fusion_tolerance * 1.1_r8
-          !write(iulog,*) 'maxcohorts exceeded',dynamic_fusion_tolerance
           !---------------------------------------------------------------------!
           ! Making profile tolerance larger means that more fusion will happen  !
           !---------------------------------------------------------------------!        
+          dynamic_fusion_tolerance = dynamic_fusion_tolerance * 1.1_r8
+
+          write(iulog,*) 'maxcohorts exceeded',dynamic_fusion_tolerance
+
        else
           iterate = 0
        endif
@@ -856,6 +906,10 @@ subroutine copy_cohort( currentCohort,copyc )
     n%gpp_clm         = o%gpp_clm
     n%npp             = o%npp
     n%npp_clm         = o%npp_clm
+
+    if ( DEBUG ) write(iulog,*) 'EDcohortDyn Ia ',o%npp_acc
+    if ( DEBUG ) write(iulog,*) 'EDcohortDyn Ib ',o%resp_acc
+
     n%npp_acc         = o%npp_acc
     n%resp_clm        = o%resp_clm
     n%resp_acc        = o%resp_acc
@@ -896,6 +950,9 @@ subroutine copy_cohort( currentCohort,copyc )
     n%dbalivedt       = o%dbalivedt
     n%dbdeaddt        = o%dbdeaddt
     n%dbstoredt       = o%dbstoredt
+
+    if ( DEBUG ) write(iulog,*) 'EDCohortDyn dpstoredt ',o%dbstoredt
+
     n%storage_flux    = o%storage_flux
 
     ! FIRE 

biogeochem/EDGrowthFunctionsMod.F90

@@ -146,7 +146,7 @@ real(r8) function tree_lai( cohort_in )
 
     if( cohort_in%status_coh  ==  2 ) then ! are the leaves on? 
        slat = 1000.0_r8 * pftcon%slatop(cohort_in%pft) ! m2/g to m2/kg
-       cohort_in%c_area = c_area(cohort_in) ! call the tree area 
+       cohort_in%c_area = c_area(cohort_in) ! call the tree area
        leafc_per_unitarea = cohort_in%bl/(cohort_in%c_area/cohort_in%n) !KgC/m2
        if(leafc_per_unitarea > 0.0_r8)then
           tree_lai = leafc_per_unitarea * slat  !kg/m2 * m2/kg = unitless LAI 

biogeochem/EDPatchDynamicsMod.F90

@@ -743,6 +743,9 @@ subroutine create_patch(currentSite, new_patch, age, areap, spread_local,cwd_ag_
     new_patch%cwd_ag_in(:)       = 0._r8
     new_patch%cwd_bg_in(:)       = 0._r8
 
+    new_patch%cwd_ag_out(:)      = 0._r8
+    new_patch%cwd_bg_out(:)      = 0._r8
+
     new_patch%f_sun              = 0._r8
     new_patch%ed_laisun_z(:,:,:) = 0._r8 
     new_patch%ed_laisha_z(:,:,:) = 0._r8 
@@ -761,6 +764,9 @@ subroutine create_patch(currentSite, new_patch, age, areap, spread_local,cwd_ag_
     new_patch%total_tree_area    = 0.0_r8  
     new_patch%NCL_p              = 1
 
+    new_patch%leaf_litter_in(:)  = 0._r8
+    new_patch%leaf_litter_out(:) = 0._r8
+
     allocate(new_patch%rootfr_ft(numpft_ed,nlevgrnd))
     allocate(new_patch%rootr_ft(numpft_ed,nlevgrnd)) 
 
@@ -802,7 +808,7 @@ subroutine zero_patch(cp_p)
     currentPatch%bare_frac_area             = nan                             
 
     currentPatch%tlai_profile(:,:,:)        = nan 
-    currentPatch%elai_profile(:,:,:)        = nan 
+    currentPatch%elai_profile(:,:,:)        = 0._r8 
     currentPatch%tsai_profile(:,:,:)        = nan 
     currentPatch%esai_profile(:,:,:)        = nan       
     currentPatch%canopy_area_profile(:,:,:) = nan       
@@ -816,7 +822,7 @@ subroutine zero_patch(cp_p)
     currentPatch%ed_laisha_z(:,:,:)         = nan 
     currentPatch%ed_parsun_z(:,:,:)         = nan 
     currentPatch%ed_parsha_z(:,:,:)         = nan 
-    currentPatch%psn_z(:,:,:)               = nan   
+    currentPatch%psn_z(:,:,:)               = 0._r8   
 
     currentPatch%f_sun(:,:,:)               = nan
     currentPatch%tr_soil_dir(:)             = nan    ! fraction of incoming direct  radiation that is transmitted to the soil as direct
@@ -847,6 +853,9 @@ subroutine zero_patch(cp_p)
     currentPatch%root_litter(:)             = 0.0_r8
     currentPatch%leaf_litter(:)             = 0.0_r8
 
+    currentPatch%leaf_litter_in(:)          = 0.0_r8
+    currentPatch%leaf_litter_out(:)         = 0.0_r8
+
     ! FIRE
     currentPatch%fuel_eff_moist             = 0.0_r8 ! average fuel moisture content of the ground fuel 
     ! (incl. live grasses. omits 1000hr fuels)

biogeochem/EDPhenologyType.F90

@@ -15,18 +15,22 @@ module EDPhenologyType
   use decompMod         , only : bounds_type
   use accumulMod        , only : update_accum_field, extract_accum_field, accumResetVal
   use clm_varctl        , only : iulog
-  use clm_time_manager  , only : get_nstep, get_step_size
+  use clm_time_manager  , only : get_nstep, get_step_size, is_restart
   !
   ! !USES:
   implicit none
   private
   !
   type, public :: ed_phenology_type
+
+     logical   ::  DEBUG = .false.
+
      !
      ! change these to allocatable
      ! add a rbuf variable that is a part of this type
      !
      real(r8), pointer :: ED_GDD_patch               (:)   ! ED Phenology growing degree days.
+
      ! This (phen_cd_status_patch?) could and should be site-level. RF
      integer , pointer :: phen_cd_status_patch       (:)   ! ED Phenology cold deciduous status
      character(10)     :: accString = 'ED_GDD0'
@@ -161,6 +165,8 @@ subroutine accumulateAndExtract( this, bounds,     &
     call update_accum_field  ( trim(this%accString), rbufslp, get_nstep() )
     call extract_accum_field ( trim(this%accString), this%ED_GDD_patch, get_nstep() )
 
+    if (this%DEBUG) write(iulog,*) 'ED_GDD accumAndExtract ', this%ED_GDD_patch
+
     deallocate(rbufslp)
 
   end subroutine accumulateAndExtract
@@ -305,6 +311,10 @@ subroutine initAccVars(this, bounds)
     call extract_accum_field (this%accString, rbufslp, get_nstep())
     this%ED_GDD_patch(bounds%begp:bounds%endp) = rbufslp(bounds%begp:bounds%endp)
 
+    if ( this%DEBUG ) then
+       write(iulog,*) 'ED_GDD initAccVars ',this%ED_GDD_patch(bounds%begp:bounds%endp)
+    endif
+
     deallocate(rbufslp)
 
   end subroutine initAccVars