Merge pull request #347 from rgknox/rgknox-parameter-updates

updates to the parameter file

Author: rgknox

biogeochem/EDCanopyStructureMod.F90

@@ -148,20 +148,6 @@ subroutine canopy_structure( currentSite , bc_in )
          enddo
 
 
-         ! ------------------------------------------------------------------------------
-         ! Check patch area to prevent numerical weirdness
-         ! ------------------------------------------------------------------------------
-         
-         if (currentPatch%area .lt. min_patch_area) then  
-            
-            write(fates_log(),*) 'An incredibly small patch exists that should'
-            write(fates_log(),*) 'had been fused or culled already'
-            write(fates_log(),*) 'currentPatch%area: ',currentPatch%area
-            write(fates_log(),*) 'min_patch_area: ',min_patch_area
-            call endrun(msg=errMsg(sourcefile, __LINE__))
-            
-         end if
-
          ! Does any layer have excess area in it? Keep going until it does not...
          patch_area_counter = 0
          area_not_balanced = .true.
@@ -1098,6 +1084,13 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si)
        currentPatch%canopy_area_profile(:,:,:)  = 0._r8       
        currentPatch%canopy_mask(:,:)            = 0
 
+       ! ------------------------------------------------------------------------------
+       ! It is remotely possible that in deserts we will not have any canopy
+       ! area, ie not plants at all...
+       ! ------------------------------------------------------------------------------
+       
+       if (currentPatch%total_canopy_area > tiny(currentPatch%total_canopy_area)) then
+
        currentCohort => currentPatch%shortest
        do while(associated(currentCohort)) 
 
@@ -1226,12 +1219,6 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si)
           ! and canopy area to the accumulators. 
           ! -----------------------------------------------------------------------------
 
-          ! ------------------------------------------------------------------------------
-          ! It is remotely possible that in deserts we will not have any canopy
-          ! area, ie not plants at all...
-          ! ------------------------------------------------------------------------------
-          
-          if (currentPatch%total_canopy_area > tiny(currentPatch%total_canopy_area)) then
 
              currentCohort => currentPatch%shortest
              do while(associated(currentCohort))   
@@ -1389,7 +1376,7 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si)
              do cl = 1,currentPatch%NCL_p
                 do iv = 1,currentPatch%ncan(cl,ft)
 
-                   if( sum(currentPatch%canopy_area_profile(cl,:,iv)) > 1.0001_r8 ) then
+                   if( DEBUG .and. sum(currentPatch%canopy_area_profile(cl,:,iv)) > 1.0001_r8 ) then
 
                       write(fates_log(), *) 'FATES: A canopy_area_profile exceeded 1.0'
                       write(fates_log(), *) 'cl: ',cl

biogeochem/EDLoggingMortalityMod.F90

@@ -31,6 +31,7 @@ module EDLoggingMortalityMod
    use EDParamsMod       , only : logging_direct_frac
    use EDParamsMod       , only : logging_mechanical_frac 
    use EDParamsMod       , only : logging_coll_under_frac 
+   use EDParamsMod       , only : logging_dbhmax_infra
    use FatesInterfaceMod , only : hlm_current_year
    use FatesInterfaceMod , only : hlm_current_month
    use FatesInterfaceMod , only : hlm_current_day
@@ -154,7 +155,6 @@ subroutine LoggingMortality_frac( pft_i, dbh, lmort_direct,lmort_collateral,lmor
 
       ! Parameters
       real(r8), parameter   :: adjustment = 1.0 ! adjustment for mortality rates
-      real(r8), parameter   :: logging_dbhmax_infra = 35 !(cm), based on Feldpaush et al. (2005) and Ferry et al. (2010)
 
       if (logging_time) then 
          if(EDPftvarcon_inst%woody(pft_i) == 1)then ! only set logging rates for trees

biogeochem/EDMortalityFunctionsMod.F90

@@ -12,7 +12,6 @@ module EDMortalityFunctionsMod
    use EDTypesMod            , only : ed_patch_type
    use FatesConstantsMod     , only : itrue,ifalse
    use FatesAllometryMod     , only : bleaf
-   use EDParamsMod           , only : ED_val_stress_mort
    use FatesInterfaceMod     , only : bc_in_type
    use FatesInterfaceMod     , only : hlm_use_ed_prescribed_phys
    use FatesInterfaceMod     , only : hlm_freq_day
@@ -62,7 +61,6 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort )
     real(r8) :: hf_sm_threshold    ! hydraulic failure soil moisture threshold 
     real(r8) :: temp_dep           ! Temp. function (freezing mortality)
     real(r8) :: temp_in_C          ! Daily averaged temperature in Celcius
-    real(r8),parameter :: frost_mort_scaler = 3.0_r8  ! Scaling factor for freezing mortality
     real(r8),parameter :: frost_mort_buffer = 5.0_r8  ! 5deg buffer for freezing mortality
 
     logical, parameter :: test_zero_mortality = .false. ! Developer test which
@@ -79,7 +77,7 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort )
     hf_sm_threshold = EDPftvarcon_inst%hf_sm_threshold(cohort_in%pft)
 
     if(cohort_in%patchptr%btran_ft(cohort_in%pft) <= hf_sm_threshold)then 
-       hmort = ED_val_stress_mort
+       hmort = EDPftvarcon_inst%mort_scalar_hydrfailure(cohort_in%pft)
      else
        hmort = 0.0_r8
      endif 
@@ -89,7 +87,8 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort )
        call bleaf(cohort_in%dbh,cohort_in%pft,cohort_in%canopy_trim,b_leaf)
        if( b_leaf > 0._r8 .and. cohort_in%bstore <= b_leaf )then
           frac = cohort_in%bstore/ b_leaf
-          cmort = max(0.0_r8,ED_val_stress_mort*(1.0_r8 - frac))
+          cmort = max(0.0_r8,EDPftvarcon_inst%mort_scalar_cstarvation(cohort_in%pft) * &
+               (1.0_r8 - frac))
         else
           cmort = 0.0_r8
        endif
@@ -109,7 +108,7 @@ subroutine mortality_rates( cohort_in,bc_in,cmort,hmort,bmort,frmort )
     temp_in_C = bc_in%t_veg24_si - tfrz
     temp_dep  = max(0.0,min(1.0,1.0 - (temp_in_C - &
                 EDPftvarcon_inst%freezetol(cohort_in%pft))/frost_mort_buffer) )
-    frmort    = frost_mort_scaler * temp_dep
+    frmort    = EDPftvarcon_inst%mort_scalar_coldstress(cohort_in%pft) * temp_dep
 
 
     !mortality_rates = bmort + hmort + cmort

biogeochem/EDPhysiologyMod.F90

@@ -872,15 +872,9 @@ subroutine PlantGrowth( currentSite, currentCohort, bc_in )
     integer , parameter :: max_substeps = 300
     real(r8), parameter :: max_trunc_error = 1.0_r8
     integer,  parameter :: ODESolve = 2    ! 1=RKF45,  2=Euler
-    real(r8), parameter :: global_branch_turnover = 0.0_r8 ! Temporary branch turnover setting
-                                                           ! Branch-turnover control will be 
-                                                           ! introduced in a later PR
-
 
     ipft = currentCohort%pft
 
-    EDPftvarcon_inst%branch_turnover(ipft) = global_branch_turnover
-
     ! Initialize seed production
     currentCohort%seed_prod  = 0.0_r8
 
@@ -973,35 +967,6 @@ subroutine PlantGrowth( currentSite, currentCohort, bc_in )
              currentCohort%canopy_trim, currentCohort%dbh, currentCohort%hite )
     end if
 
-    ! -----------------------------------------------------------------------------------
-    ! III(a). Calculate the maintenance turnover demands 
-    !       Pre-check, make sure phenology is mutually exclusive and at least one chosen
-    !       (MOVE THIS TO THE PARAMETER READ-IN SECTION)
-    ! -----------------------------------------------------------------------------------
-
-    if (EDPftvarcon_inst%evergreen(ipft) == 1) then
-       if (EDPftvarcon_inst%season_decid(ipft) == 1)then 
-          write(fates_log(),*) 'PFT # ',ipft,' was specified as being both evergreen'
-          write(fates_log(),*) '       and seasonally deciduous, impossible, aborting'
-          call endrun(msg=errMsg(sourcefile, __LINE__))
-       end if
-       if (EDPftvarcon_inst%stress_decid(ipft) == 1)then 
-          write(fates_log(),*) 'PFT # ',ipft,' was specified as being both evergreen'
-          write(fates_log(),*) '       and stress deciduous, impossible, aborting'
-          call endrun(msg=errMsg(sourcefile, __LINE__))
-       end if
-    end if
-    if (EDPftvarcon_inst%stress_decid(ipft) /= 1 .and. &
-        EDPftvarcon_inst%season_decid(ipft) /= 1 .and. &
-        EDPftvarcon_inst%evergreen(ipft)    /= 1) then
-       write(fates_log(),*) 'PFT # ',ipft,' must be defined as having one of three'
-       write(fates_log(),*) 'phenology habits, ie == 1'
-       write(fates_log(),*) 'stress_decid: ',EDPftvarcon_inst%stress_decid(ipft)
-       write(fates_log(),*) 'season_decid: ',EDPftvarcon_inst%season_decid(ipft)
-       write(fates_log(),*) 'evergreen: ',EDPftvarcon_inst%evergreen(ipft)
-       call endrun(msg=errMsg(sourcefile, __LINE__))
-    endif
-    
 
     ! -----------------------------------------------------------------------------------
     ! III(b). Calculate the maintenance turnover demands 
@@ -1460,7 +1425,7 @@ function AllomCGrowthDeriv(c_pools,c_mask,cbalance,currentCohort) result(dCdx)
         if (dbh <= EDPftvarcon_inst%dbh_repro_threshold(ipft)) then ! cap on leaf biomass
            repro_fraction = EDPftvarcon_inst%seed_alloc(ipft)
         else
-           repro_fraction = EDPftvarcon_inst%seed_alloc(ipft) + EDPftvarcon_inst%clone_alloc(ipft)
+           repro_fraction = EDPftvarcon_inst%seed_alloc(ipft) + EDPftvarcon_inst%seed_alloc_mature(ipft)
         end if
 
         dCdx = 0.0_r8

biogeophys/FatesPlantRespPhotosynthMod.F90

@@ -191,7 +191,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime)
     ! Ball-Berry minimum leaf conductance, unstressed (umol H2O/m**2/s)
     ! For C3 and C4 plants
     ! -----------------------------------------------------------------------------------
-    real(r8), dimension(2) :: bbbopt 
+    real(r8), dimension(0:1) :: bbbopt 
 
     associate(  &
          c3psn     => EDPftvarcon_inst%c3psn  , &
@@ -202,9 +202,9 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime)
          frootcn   => EDPftvarcon_inst%frootcn, & ! froot C:N (gc/gN)   ! slope of BB relationship
          q10       => FatesSynchronizedParamsInst%Q10 )
 
+      bbbopt(0) = ED_val_bbopt_c4
       bbbopt(1) = ED_val_bbopt_c3
-      bbbopt(2) = ED_val_bbopt_c4
-
+      
       do s = 1,nsites
 
          ! Multi-layer parameters scaled by leaf nitrogen profile.
@@ -721,7 +721,8 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl,         &  ! in
 
    ! Locals
    ! ------------------------------------------------------------------------
-   integer :: pp_type            ! Index for the different photosynthetic pathways C3,C4
+   integer :: c3c4_path_index    ! Index for which photosynthetic pathway 
+                                 ! is active.  C4 = 0,  C3 = 1
    integer :: sunsha             ! Index for differentiating sun and shade
    real(r8) :: gstoma            ! Stomatal Conductance of this leaf layer (m/s)
    real(r8) :: agross            ! co-limited gross leaf photosynthesis (umol CO2/m**2/s)
@@ -758,22 +759,23 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl,         &  ! in
    real(r8),parameter :: init_a2l_co2_c3 = 0.7_r8
    real(r8),parameter :: init_a2l_co2_c4 = 0.4_r8
 
-   ! quantum efficiency, used only for C4 (mol CO2 / mol photons)
-   real(r8),parameter,dimension(2) :: quant_eff = [0.0_r8,0.05_r8]
+   ! quantum efficiency, used only for C4 (mol CO2 / mol photons) (index 0)
+   real(r8),parameter,dimension(0:1) :: quant_eff = [0.05_r8,0.0_r8]
 
    ! empirical curvature parameter for ac, aj photosynthesis co-limitation
-   real(r8),parameter,dimension(2) :: theta_cj  = [0.98_r8,0.80_r8]
+   real(r8),parameter,dimension(0:1) :: theta_cj  = [0.80_r8,0.98_r8]
 
    ! empirical curvature parameter for ap photosynthesis co-limitation
    real(r8),parameter :: theta_ip = 0.95_r8
 
-   associate( bb_slope  => EDPftvarcon_inst%BB_slope ) ! slope of BB relationship
+   associate( bb_slope  => EDPftvarcon_inst%BB_slope)    ! slope of BB relationship
+
+     ! photosynthetic pathway: 0. = c4, 1. = c3
+     c3c4_path_index = nint(EDPftvarcon_inst%c3psn(ft))
 
-     if (nint(EDPftvarcon_inst%c3psn(ft)) == 1) then! photosynthetic pathway: 0. = c4, 1. = c3
-        pp_type = 1
+     if (c3c4_path_index == 1) then
         init_co2_intra_c = init_a2l_co2_c3 * can_co2_ppress
      else
-        pp_type = 2
         init_co2_intra_c = init_a2l_co2_c4 * can_co2_ppress
      end if
 
@@ -845,7 +847,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl,         &  ! in
                  co2_intra_c_old = co2_intra_c
 
                  ! Photosynthesis limitation rate calculations 
-                 if (pp_type == 1)then    
+                 if (c3c4_path_index == 1)then    
 
                     ! C3: Rubisco-limited photosynthesis
                     ac = vcmax * max(co2_intra_c-co2_cpoint, 0._r8) / &
@@ -867,14 +869,14 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl,         &  ! in
                     if(sunsha == 1)then !sunlit
                        !guard against /0's in the night.
                        if((laisun_lsl * canopy_area_lsl) > 0.0000000001_r8) then   
-                          aj = quant_eff(pp_type) * parsun_lsl * 4.6_r8
+                          aj = quant_eff(c3c4_path_index) * parsun_lsl * 4.6_r8
                           !convert from per cohort to per m2 of leaf)
                           aj = aj / (laisun_lsl * canopy_area_lsl)
                        else
                           aj = 0._r8
                        end if
                     else
-                       aj = quant_eff(pp_type) * parsha_lsl * 4.6_r8
+                       aj = quant_eff(c3c4_path_index) * parsha_lsl * 4.6_r8
                        aj = aj / (laisha_lsl * canopy_area_lsl)
                     end if
 
@@ -884,7 +886,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl,         &  ! in
                  end if
 
                  ! Gross photosynthesis smoothing calculations. First co-limit ac and aj. Then co-limit ap
-                 aquad = theta_cj(pp_type)
+                 aquad = theta_cj(c3c4_path_index)
                  bquad = -(ac + aj)
                  cquad = ac * aj
                  call quadratic_f (aquad, bquad, cquad, r1, r2)