Merge pull request #1424 from rgknox/nclmax-3-sortfixes

promotion/demotion refactor and nclmax = 3

Author: glemieux

biogeochem/EDCanopyStructureMod.F90

@@ -340,12 +340,13 @@ subroutine terminate_cohorts( currentSite, currentPatch, level , call_index, bc_
           terminate = itrue
           termination_type = i_term_mort_type_numdens
           if ( debug ) then
-             write(fates_log(),*) 'terminating cohorts 0',currentCohort%n/currentPatch%area,currentCohort%dbh,currentCohort%pft,call_index
+             write(fates_log(),*) 'terminating cohorts 0',currentCohort%n/currentPatch%area, &
+                  currentCohort%dbh,currentCohort%pft,call_index
           endif
        endif
 
        ! The rest of these are only allowed if we are not dealing with a recruit (level 2)
-       if (.not.currentCohort%isnew .and. level == 2) then
+       if_level_2: if (.not.currentCohort%isnew .and. level == 2) then
 
          ! Not enough n or dbh
          if  (currentCohort%n/currentPatch%area <= min_npm2 .or.	&  !
@@ -354,18 +355,13 @@ subroutine terminate_cohorts( currentSite, currentPatch, level , call_index, bc_
             terminate = itrue
             termination_type = i_term_mort_type_numdens
             if ( debug ) then
-               write(fates_log(),*) 'terminating cohorts 1',currentCohort%n/currentPatch%area,currentCohort%dbh,currentCohort%pft,call_index
+               write(fates_log(),*) 'terminating cohorts 1', &
+                    currentCohort%n/currentPatch%area,currentCohort%dbh, &
+                    currentCohort%pft,call_index
             endif
          endif
 
-         ! Outside the maximum canopy layer
-         if (currentCohort%canopy_layer > nclmax ) then
-           terminate = itrue
-           termination_type = i_term_mort_type_canlev
-           if ( debug ) then
-             write(fates_log(),*) 'terminating cohorts 2', currentCohort%canopy_layer,currentCohort%pft,call_index
-           endif
-         endif
+        
 
          ! live biomass pools are terminally depleted
          if ( ( sapw_c+leaf_c+fnrt_c ) < 1e-10_r8  .or.  &
@@ -387,8 +383,18 @@ subroutine terminate_cohorts( currentSite, currentPatch, level , call_index, bc_
                     struct_c,sapw_c,leaf_c,fnrt_c,store_c,currentCohort%pft,call_index
             endif
 
-        endif
-      endif    !  if (.not.currentCohort%isnew .and. level == 2) then
+         endif
+
+      end if if_level_2
+      
+      ! Outside the maximum canopy layer
+      if (currentCohort%canopy_layer > nclmax .and. level == 3) then
+         terminate = itrue
+         termination_type = i_term_mort_type_canlev
+         if ( debug ) then
+            write(fates_log(),*) 'terminating cohorts 2', currentCohort%canopy_layer,currentCohort%pft,call_index
+         endif
+      endif
 
       if (terminate == itrue) then
          call terminate_cohort(currentSite, currentPatch, currentCohort, bc_in, termination_type)

biogeochem/EDCohortDynamicsMod.F90

@@ -21,7 +21,7 @@ module FatesPatchMod
   use PRTGenericMod,          only : struct_organ, leaf_organ, sapw_organ
   use PRTParametersMod,       only : prt_params
   use FatesConstantsMod,      only : nocomp_bareground
-  use EDParamsMod,            only : nlevleaf, nclmax, maxpft
+  use EDParamsMod,            only : nlevleaf, nclmax, maxpft,max_cohort_per_patch
   use FatesConstantsMod,      only : n_dbh_bins, n_dist_types
   use FatesConstantsMod,      only : t_water_freeze_k_1atm
   use FatesRunningMeanMod,    only : ema_24hr, fixed_24hr, ema_lpa, ema_longterm
@@ -41,14 +41,35 @@ module FatesPatchMod
   ! for error message writing
   character(len=*), parameter :: sourcefile = __FILE__
 
+  type :: fates_cohort_vec_type
+
+     ! This is a scratch array for cohort pointers
+     ! this is useful if you want to loop over a sparse subset
+     ! of fates cohorts over and over again, allowing
+     ! you to iterate them in a do loop
+     
+     type(fates_cohort_type), pointer :: p => null()
+
+     ! This is the area of the cohort (less than or equal to cohort%carea)
+     ! that will be promoted or demoted, ie promoted/demoted crown area
+     ! units [m2/site] or [m2/ha] (same as the patch area and crown area)
+     ! We track it here because we construct the cohort list for specific
+     ! canopy layers
+     
+     real(r8) :: pd_area
+     
+  end type fates_cohort_vec_type
+
+  
   type, public :: fates_patch_type
 
     ! POINTERS
     type (fates_cohort_type), pointer :: tallest => null()  ! pointer to patch's tallest cohort    
     type (fates_cohort_type), pointer :: shortest => null() ! pointer to patch's shortest cohort
     type (fates_patch_type),  pointer :: older => null()    ! pointer to next older patch   
     type (fates_patch_type),  pointer :: younger => null()  ! pointer to next younger patch
-  
+    type (fates_cohort_vec_type), pointer :: co_scr(:)      ! Scratch vector of cohort properties
+    
     !---------------------------------------------------------------------------
 
     ! INDICES
@@ -279,7 +300,8 @@ subroutine Init(this, num_swb, num_levsoil)
       allocate(this%sabs_dir(num_swb))
       allocate(this%sabs_dif(num_swb))
       allocate(this%fragmentation_scaler(num_levsoil))
-
+      allocate(this%co_scr(max_cohort_per_patch))
+      
       ! initialize all values to nan
       call this%NanValues()
 
@@ -898,6 +920,7 @@ subroutine FreeMemory(this, regeneration_model, numpft)
                  this%sabs_dir,                 &
                  this%sabs_dif,                 &
                  this%fragmentation_scaler,     &
+                 this%co_scr,                   &
                  stat=istat, errmsg=smsg)
 
       ! These arrays are allocated via a call from EDCanopyStructureMod
@@ -1161,19 +1184,29 @@ end subroutine CountCohorts
 
     !===========================================================================
 
-    subroutine SortCohorts(this)
+    subroutine SortCohorts(this,check_order)
       !
       ! DESCRIPTION: sort cohorts in patch's linked list
       ! uses insertion sort to build a new list
       !
 
       ! ARGUMENTS:
       class(fates_patch_type), intent(inout), target :: this ! patch
+
+      logical, optional, intent(in) :: check_order
 
       ! LOCALS:
       type(fates_cohort_type), pointer :: currentCohort
       type(fates_cohort_type), pointer :: nextCohort
+
+      logical :: check_order_present
 
+      if (present(check_order)) then
+        check_order_present = check_order
+      else
+        check_order_present = .false.
+      end if
+
       ! check for inconsistent list state
       if (.not. associated(this%shortest) .and. .not. associated(this%tallest)) then
           ! empty list
@@ -1186,6 +1219,23 @@ subroutine SortCohorts(this)
 
       ! hold on to current linked list so we don't lose it
       currentCohort => this%shortest
+
+      if(check_order_present)then
+         do while (associated(currentCohort))
+            if( associated(currentCohort%taller)) then
+               if(currentCohort%height > currentCohort%taller%height)then
+                  write(fates_log(),*) 'Cohort sort checking has failed,'
+                  write(fates_log(),*) 'they are not in height order:'
+                  write(fates_log(),*) 'current: ',currentCohort%height
+                  write(fates_log(),*) 'taller: ',currentCohort%taller%height
+                  call endrun(msg=errMsg(sourcefile, __LINE__))
+               end if
+            end if
+            currentCohort => currentCohort%taller
+         end do
+         return
+      end if
+
 
       ! reset the current list: we'll build it incrementally
       this%shortest => null()

biogeochem/FatesPatchMod.F90

@@ -44,7 +44,7 @@ module EDParamsMod
                                                                       !moving average of par at the seedling layer used to 
                                                                       !calculate seedling to sapling transition rates
    real(r8),protected, public :: fates_mortality_disturbance_fraction ! the fraction of canopy mortality that results in disturbance
-   real(r8),protected, public :: ED_val_comp_excln                    ! weighting factor for canopy layer exclusion and promotion
+   real(r8),protected, public :: comp_excln_exp                       ! weighting factor (exponent) for canopy layer exclusion and promotion
    real(r8),protected, public :: ED_val_vai_top_bin_width             ! width in VAI units of uppermost leaf+stem layer scattering element
    real(r8),protected, public :: ED_val_vai_width_increase_factor     ! factor by which each leaf+stem scattering element increases in VAI width
    real(r8),protected, public :: ED_val_nignitions                    ! number of annual ignitions per square km
@@ -81,10 +81,10 @@ module EDParamsMod
 
    real(r8), parameter, public :: soil_tfrz_thresh = -2.0_r8          ! Soil temperature threshold below which hydraulic failure mortality is off (non-hydro only) in degrees C
 
-   integer, parameter, public :: nclmax = 2                           ! Maximum number of canopy layers (used only for scratch arrays)
-                                                                      ! We would make this even higher, but making this
-                                                                      ! a little lower keeps the size down on some output arrays
-                                                                      ! For large arrays at patch level we use dynamic allocation
+   integer, parameter, public :: nclmax = 3   ! Maximum number of canopy layers allowed
+                                              ! We would make this even higher, but making this
+                                              ! a little lower keeps the size down on some output arrays
+                                              ! For large arrays at patch level we use dynamic allocation
 
                                                                       ! parameters that govern the VAI (LAI+SAI) bins used in radiative transfer code
    integer, parameter, public :: nlevleaf = 30                        ! number of leaf+stem layers in each canopy layer
@@ -303,7 +303,7 @@ subroutine FatesParamsInit()
     sdlng2sap_par_timescale               = nan
     photo_temp_acclim_thome_time          = nan
     fates_mortality_disturbance_fraction  = nan
-    ED_val_comp_excln                     = nan
+    comp_excln_exp                        = nan
     ED_val_vai_top_bin_width              = nan
     ED_val_vai_width_increase_factor      = nan
     ED_val_nignitions                     = nan
@@ -610,7 +610,7 @@ subroutine FatesReceiveParams(fates_params)
           data=fates_mortality_disturbance_fraction)
 
     call fates_params%RetrieveParameter(name=ED_name_comp_excln, &
-         data=ED_val_comp_excln)
+         data=comp_excln_exp)
 
     call fates_params%RetrieveParameter(name=ED_name_vai_top_bin_width, &
          data=ED_val_vai_top_bin_width)
@@ -825,7 +825,7 @@ subroutine FatesReportParams(is_master)
         write(fates_log(),fmt0) 'photo_temp_acclim_thome_time (years) = ',photo_temp_acclim_thome_time
         write(fates_log(),fmti) 'hydr_htftype_node = ',hydr_htftype_node
         write(fates_log(),fmt0) 'fates_mortality_disturbance_fraction = ',fates_mortality_disturbance_fraction
-        write(fates_log(),fmt0) 'ED_val_comp_excln = ',ED_val_comp_excln
+        write(fates_log(),fmt0) 'comp_excln_exp = ',comp_excln_exp
         write(fates_log(),fmt0) 'ED_val_vai_top_bin_width = ',ED_val_vai_top_bin_width
         write(fates_log(),fmt0) 'ED_val_vai_width_increase_factor = ',ED_val_vai_width_increase_factor
         write(fates_log(),fmt0) 'ED_val_nignitions = ',ED_val_nignitions

main/EDParamsMod.F90

@@ -55,7 +55,7 @@ module FatesHistoryInterfaceMod
   use FatesInterfaceTypesMod        , only : hlm_freq_day
   use FatesInterfaceTypesMod        , only : hlm_parteh_mode
   use FatesInterfaceTypesMod        , only : hlm_use_sp
-  use EDParamsMod              , only : ED_val_comp_excln
+  use EDParamsMod              , only : comp_excln_exp
   use EDParamsMod              , only : ED_val_phen_coldtemp
   use EDParamsMod                   , only : nlevleaf
   use EDParamsMod               , only : ED_val_history_height_bin_edges
@@ -2693,7 +2693,7 @@ subroutine update_history_dyn_sitelevel(this,nc,nsites,sites)
             hio_ncl_si(io_si) = hio_ncl_si(io_si) + cpatch%ncl_p * cpatch%area * AREA_INV
 
             ! only valid when "strict ppa" enabled
-            if ( ED_val_comp_excln .lt. 0._r8 ) then
+            if ( comp_excln_exp .lt. 0._r8 ) then
                hio_zstar_si(io_si) = hio_zstar_si(io_si) &
                     + cpatch%zstar * cpatch%area * AREA_INV
             end if
@@ -4861,7 +4861,7 @@ subroutine update_history_dyn_subsite_ageclass(this,nc,nsites,sites)
           end do
 
           ! only valid when "strict ppa" enabled
-          if ( ED_val_comp_excln .lt. 0._r8 ) then
+          if ( comp_excln_exp .lt. 0._r8 ) then
              hio_zstar_si_age(io_si,cpatch%age_class) = hio_zstar_si_age(io_si,cpatch%age_class) &
                   + cpatch%zstar * patch_area_div_site_area
           end if
@@ -7347,7 +7347,7 @@ subroutine define_history_vars(this, initialize_variables)
                upfreq=group_dyna_complx, ivar=ivar, initialize=initialize_variables,                 &
                index=ih_npatches_si_age)
 
-          if ( ED_val_comp_excln .lt. 0._r8 ) then ! only valid when "strict ppa" enabled
+          if ( comp_excln_exp .lt. 0._r8 ) then ! only valid when "strict ppa" enabled
              tempstring = 'active'
           else
              tempstring = 'inactive'

main/FatesHistoryInterfaceMod.F90

@@ -24,7 +24,7 @@ module test_SortCohorts
       class(TestSortCohorts),  intent(inout) :: this  ! test object
       type(fates_patch_type)                 :: patch ! patch object
 
-      ! sort cohorts - should pass
+      ! sort cohorts - should pass - the argument
       call patch%SortCohorts()
 
     end subroutine EmptyList_SortCohorts_Passes