Introduce fates-host parameter interface

Merge remote-tracking branch 'pr/andre-ed-params'

Introduce an interface to pass parameter information from the host to
fates. Fates registers a set of parameters that it needs read in, and
indicates if they are fates only parameters, or need to be synced with
values from the host. The host reads the parameters and returns them
to fates.

This refactor attempted to be as minimally invasive as possible to the
fates science code. All existing storage and conventions for fates
parameters were left in place. The only exception was the
consolidation of all pft dimensioned parameters into the EDpftvarcon
type. Fates no longer uses variables from the host pftcon type.

This introduces dynamic allocation of pft level parameter in
preparation for setting the number of pfts at run time, but still
requires a hard coded number of pfts until the host code can be
modified.

Note that the default clm and old clm-ed parameter files have diverged
before this work began. To do this development in a bit for bit way,
the clm-ed parameter file was updated to agree with the clm default
parameter file. This is answer changing compared to the fates master
branch, but code was refactored in a bit for bit way. No netcdf
variables were added or removed in this PR.

Fixes: #155, #156

User interface changes?: Yes.

1. Users will need to update custom parameter files. This introduces a
new namelist variable, fates_paramfile. The fates parameters are
**always** read from the netcdf file pointed to by
fates_paramfile. All host parameters are **always** read from the
netcdf file pointed to by paramfile. The host paramfile and fates
paramfile **may** point to the same netcdf file.

2. All fates parameters and dimensions are now name spaced by
'fates_'. The variable names have remained the same, but the dimension
'param' is now 'fates_scalar'. A new default parameter file with the
required changes is available. See fates_params.c170308.nc in the
input data repo.

Code review: self. Discussion with clm-cmt, code walk through with
ckoven, rgknox, rosie

Testing:

  andre: 2017-03-14

    Test suite: ed - yellowstone gnu, intel, pgi
                     hobart nag
    Test baseline: 0ea3fed
    Test namelist changes: addition of fates_paramfile
    Test answer changes: bit for bit
    Test status: all tests pass

    Test suite: clm_short - yellowstone gnu, intel, pgi
    Test baseline: clm4_5_12_r195
    Test namelist changes: none
    Test answer changes: bit for bit
    Test status: all tests pass

components/clm/bld/CLMBuildNamelist.pm

@@ -1976,7 +1976,7 @@ sub setup_logic_params_file {
 
   if ( $physv->as_long() >= $physv->as_long("clm4_5") ) {
     add_default($test_files, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'paramfile', 
-                'use_ed'=>$nl_flags->{'use_ed'}, 'phys'=>$nl_flags->{'phys'},
+                'phys'=>$nl_flags->{'phys'},
                 'use_flexibleCN'=>$nl_flags->{'use_flexibleCN'} );
   } else {
     add_default($test_files, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'fpftcon');
@@ -3367,6 +3367,7 @@ sub setup_logic_ed {
 
     if ($physv->as_long() >= $physv->as_long("clm4_5") && value_is_true( $nl_flags->{'use_ed'})  ) {
  	add_default($test_files, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'use_ed_spit_fire', 'use_ed'=>$nl_flags->{'use_ed'} );
+        add_default($test_files, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl, 'fates_paramfile', 'phys'=>$nl_flags->{'phys'});
     }
 }
 

components/clm/bld/namelist_files/namelist_defaults_clm4_5.xml

@@ -240,9 +240,10 @@ attributes from the config_cache.xml file (with keys converted to upper-case).
 <!-- The default filenames are given relative to the root directory
      for the CLM2 data in the CESM distribution -->
 <!-- Plant function types (relative to {csmdata}) -->
-<paramfile phys="clm5_0" use_ed=".false.">lnd/clm2/paramdata/clm5_params.c160713.nc</paramfile>
-<paramfile phys="clm4_5" use_ed=".false.">lnd/clm2/paramdata/clm_params.c160713.nc</paramfile>
-<paramfile phys="clm4_5" use_ed=".true." >lnd/clm2/paramdata/clm_params_ed.c160808.nc</paramfile>
+<paramfile phys="clm5_0">lnd/clm2/paramdata/clm5_params.c160713.nc</paramfile>
+<paramfile phys="clm4_5">lnd/clm2/paramdata/clm_params.c160713.nc</paramfile>
+
+<fates_paramfile phys="clm4_5">lnd/clm2/paramdata/fates_params.c170308.nc</fates_paramfile>
 
 <!-- ========================================================================================  -->
 <!-- clm 5.0 BGC nitrogen model                                                                -->

components/clm/bld/namelist_files/namelist_definition_clm4_5.xml

@@ -510,6 +510,11 @@ Full pathname datafile with plant function type (PFT) constants combined with
 constants for biogeochem modules
 </entry>
 
+<entry id="fates_paramfile" type="char*256" category="datasets"
+       input_pathname="abs" group="clm_inparm" valid_values="" >
+Full pathname datafile with fates parameters
+</entry>
+
 <entry id="fsurdat" type="char*256"     category="datasets"
        input_pathname="abs" group="clm_inparm" valid_values="" >
 Full pathname of surface data file.

components/clm/src/ED/biogeochem/EDCanopyStructureMod.F90

@@ -7,7 +7,7 @@ module EDCanopyStructureMod
 
   use FatesConstantsMod     , only : r8 => fates_r8
   use FatesGlobals          , only : fates_log
-  use pftconMod             , only : pftcon
+  use EDPftvarcon             , only : EDPftvarcon_inst
   use EDGrowthFunctionsMod  , only : c_area
   use EDCohortDynamicsMod   , only : copy_cohort, terminate_cohorts, fuse_cohorts
   use EDtypesMod            , only : ed_site_type, ed_patch_type, ed_cohort_type, ncwd
@@ -651,7 +651,7 @@ subroutine canopy_spread( currentSite )
        currentCohort => currentPatch%tallest
        do while (associated(currentCohort))
           currentCohort%c_area = c_area(currentCohort) 
-          if(pftcon%woody(currentCohort%pft) == 1)then
+          if(EDPftvarcon_inst%woody(currentCohort%pft) == 1)then
              arealayer(currentCohort%canopy_layer) = arealayer(currentCohort%canopy_layer) + currentCohort%c_area
           endif
           currentCohort => currentCohort%shorter
@@ -698,7 +698,7 @@ subroutine canopy_summarization( nsites, sites, bc_in )
     use EDGrowthFunctionsMod , only : tree_lai, c_area
     use EDEcophysConType     , only : EDecophyscon
     use EDtypesMod           , only : area
-    use pftconMod            , only : pftcon
+    use EDPftvarcon            , only : EDPftvarcon_inst
 
     ! !ARGUMENTS    
     integer                 , intent(in)            :: nsites
@@ -762,7 +762,7 @@ subroutine canopy_summarization( nsites, sites, bc_in )
 
              if(currentCohort%canopy_layer==1)then
                 currentPatch%total_canopy_area = currentPatch%total_canopy_area + currentCohort%c_area
-                if(pftcon%woody(ft)==1)then
+                if(EDPftvarcon_inst%woody(ft)==1)then
                    currentPatch%total_tree_area = currentPatch%total_tree_area + currentCohort%c_area
                 endif
              endif
@@ -1024,11 +1024,11 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si)
              do iv = 1,currentCohort%NV-1 
 
                 ! what is the height of this layer? (for snow burial purposes...)  
-                ! pftcon%vertical_canopy_frac(ft))! fudge - this should be pft specific but i cant get it to compile. 
+                ! EDPftvarcon_inst%vertical_canopy_frac(ft))! fudge - this should be pft specific but i cant get it to compile. 
                 layer_top_hite = currentCohort%hite-((iv/currentCohort%NV) * currentCohort%hite * &
                       EDecophyscon%crown(currentCohort%pft) )
                 layer_bottom_hite = currentCohort%hite-(((iv+1)/currentCohort%NV) * currentCohort%hite * &
-                      EDecophyscon%crown(currentCohort%pft)) ! pftcon%vertical_canopy_frac(ft))
+                      EDecophyscon%crown(currentCohort%pft)) ! EDPftvarcon_inst%vertical_canopy_frac(ft))
 
                 fraction_exposed =1.0_r8
 
@@ -1057,10 +1057,10 @@ subroutine leaf_area_profile( currentSite , snow_depth_si, frac_sno_eff_si)
 
              !Bottom layer
              iv = currentCohort%NV
-             ! pftcon%vertical_canopy_frac(ft))! fudge - this should be pft specific but i cant get it to compile.
+             ! EDPftvarcon_inst%vertical_canopy_frac(ft))! fudge - this should be pft specific but i cant get it to compile.
              layer_top_hite = currentCohort%hite-((iv/currentCohort%NV) * currentCohort%hite * &
                    EDecophyscon%crown(currentCohort%pft) )
-             ! pftcon%vertical_canopy_frac(ft))
+             ! EDPftvarcon_inst%vertical_canopy_frac(ft))
              layer_bottom_hite = currentCohort%hite-(((iv+1)/currentCohort%NV) * currentCohort%hite * &
                    EDecophyscon%crown(currentCohort%pft))
 

components/clm/src/ED/biogeochem/EDCohortDynamicsMod.F90

@@ -9,7 +9,7 @@ module EDCohortDynamicsMod
   use FatesInterfaceMod     , only : hlm_freq_day
   use FatesConstantsMod     , only : r8 => fates_r8
   use FatesConstantsMod     , only : fates_unset_int
-  use pftconMod             , only : pftcon
+  use EDPftvarcon           , only : EDPftvarcon_inst
   use EDEcophysContype      , only : EDecophyscon
   use EDGrowthFunctionsMod  , only : c_area, tree_lai
   use EDTypesMod            , only : ed_site_type, ed_patch_type, ed_cohort_type
@@ -121,11 +121,11 @@ subroutine create_cohort(patchptr, pft, nn, hite, dbh, &
              call endrun(msg=errMsg(sourcefile, __LINE__))
     endif
 
-    if (new_cohort%siteptr%status==2 .and. pftcon%season_decid(pft) == 1) then
+    if (new_cohort%siteptr%status==2 .and. EDPftvarcon_inst%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. EDPftvarcon_inst%stress_decid(pft) == 1) then
       new_cohort%laimemory = 0.0_r8
     endif
 
@@ -201,27 +201,27 @@ subroutine allocate_live_biomass(cc_p,mode)
 
     currentCohort => cc_p
     ft = currentcohort%pft
-    leaf_frac = 1.0_r8/(1.0_r8 + EDecophyscon%sapwood_ratio(ft) * currentcohort%hite + pftcon%froot_leaf(ft))     
+    leaf_frac = 1.0_r8/(1.0_r8 + EDecophyscon%sapwood_ratio(ft) * currentcohort%hite + EDPftvarcon_inst%froot_leaf(ft))     
 
     !currentcohort%bl = currentcohort%balive*leaf_frac    
     !for deciduous trees, there are no leaves  
 
-    if (pftcon%evergreen(ft) == 1) then
+    if (EDPftvarcon_inst%evergreen(ft) == 1) then
        currentcohort%laimemory = 0._r8
        currentcohort%status_coh    = 2      
     endif    
 
     ! iagnore the root and stem biomass from the functional balance hypothesis. This is used when the leaves are 
     !fully on. 
-    !currentcohort%br  = pftcon%froot_leaf(ft) * (currentcohort%balive + currentcohort%laimemory) * leaf_frac
+    !currentcohort%br  = EDPftvarcon_inst%froot_leaf(ft) * (currentcohort%balive + currentcohort%laimemory) * leaf_frac
     !currentcohort%bsw = EDecophyscon%sapwood_ratio(ft) * currentcohort%hite *(currentcohort%balive + &
     !     currentcohort%laimemory)*leaf_frac 
 
     leaves_off_switch = 0
-    if (currentcohort%status_coh == 1.and.pftcon%stress_decid(ft) == 1.and.currentcohort%siteptr%dstatus==1) then !no leaves 
+    if (currentcohort%status_coh == 1.and.EDPftvarcon_inst%stress_decid(ft) == 1.and.currentcohort%siteptr%dstatus==1) then !no leaves 
       leaves_off_switch = 1 !drought decid
     endif
-    if (currentcohort%status_coh == 1.and.pftcon%season_decid(ft) == 1.and.currentcohort%siteptr%status==1) then !no leaves
+    if (currentcohort%status_coh == 1.and.EDPftvarcon_inst%season_decid(ft) == 1.and.currentcohort%siteptr%status==1) then !no leaves
       leaves_off_switch = 1 !cold decid
     endif
 
@@ -230,7 +230,7 @@ subroutine allocate_live_biomass(cc_p,mode)
 
        new_bl = currentcohort%balive*leaf_frac
 
-       new_br = pftcon%froot_leaf(ft) * (currentcohort%balive + currentcohort%laimemory) * leaf_frac
+       new_br = EDpftvarcon_inst%froot_leaf(ft) * (currentcohort%balive + currentcohort%laimemory) * leaf_frac
 
        new_bsw = EDecophyscon%sapwood_ratio(ft) * currentcohort%hite *(currentcohort%balive + &
             currentcohort%laimemory)*leaf_frac
@@ -255,7 +255,6 @@ subroutine allocate_live_biomass(cc_p,mode)
        currentcohort%br = new_br
        currentcohort%bsw = new_bsw
 
-
     else ! Leaves are off (leaves_off_switch==1)
 
        !the purpose of this section is to figure out the root and stem biomass when the leaves are off
@@ -265,11 +264,11 @@ subroutine allocate_live_biomass(cc_p,mode)
        !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
 
-       ideal_balive      = currentcohort%laimemory * pftcon%froot_leaf(ft) +  &
+       ideal_balive      = currentcohort%laimemory * EDPftvarcon_inst%froot_leaf(ft) +  &
             currentcohort%laimemory*  EDecophyscon%sapwood_ratio(ft) * currentcohort%hite
        ratio_balive      = currentcohort%balive / ideal_balive
 
-       new_br  = pftcon%froot_leaf(ft) * (ideal_balive + currentcohort%laimemory) * &
+       new_br  = EDpftvarcon_inst%froot_leaf(ft) * (ideal_balive + currentcohort%laimemory) * &
              leaf_frac *  ratio_balive
        new_bsw = EDecophyscon%sapwood_ratio(ft) * currentcohort%hite * &
              (ideal_balive + currentcohort%laimemory) * leaf_frac * ratio_balive
@@ -298,7 +297,7 @@ subroutine allocate_live_biomass(cc_p,mode)
             currentcohort%status_coh,currentcohort%balive
        write(fates_log(),*) 'actual vs predicted balive',ideal_balive,currentcohort%balive ,ratio_balive,leaf_frac
        write(fates_log(),*) 'leaf,root,stem',currentcohort%bl,currentcohort%br,currentcohort%bsw
-       write(fates_log(),*) 'pft',ft,pftcon%evergreen(ft),pftcon%season_decid(ft),leaves_off_switch
+       write(fates_log(),*) 'pft',ft,EDPftvarcon_inst%evergreen(ft),EDPftvarcon_inst%season_decid(ft),leaves_off_switch
     endif
     currentCohort%b   = currentCohort%bdead + currentCohort%balive
 

components/clm/src/ED/biogeochem/EDGrowthFunctionsMod.F90

@@ -8,7 +8,7 @@ module EDGrowthFunctionsMod
 
   use FatesConstantsMod, only : r8 => fates_r8
   use FatesGlobals     , only : fates_log
-  use pftconMod        , only : pftcon
+  use EDPftvarcon        , only : EDPftvarcon_inst
   use EDEcophysContype , only : EDecophyscon
   use EDTypesMod       , only : ed_cohort_type, nlevcan, dinc_ed
 
@@ -114,7 +114,7 @@ real(r8) function Bleaf( cohort_in )
     else  
        bleaf = 0.0419_r8 * (EDecophyscon%max_dbh(cohort_in%pft)**1.56) * EDecophyscon%wood_density(cohort_in%pft)**0.55_r8      
     endif  
-    slascaler = 0.03_r8/pftcon%slatop(cohort_in%pft)
+    slascaler = 0.03_r8/EDPftvarcon_inst%slatop(cohort_in%pft)
     bleaf = bleaf * slascaler
 
     !write(fates_log(),*) 'bleaf',bleaf, slascaler,cohort_in%pft
@@ -145,7 +145,7 @@ real(r8) function tree_lai( cohort_in )
     endif
 
     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
+       slat = 1000.0_r8 * EDPftvarcon_inst%slatop(cohort_in%pft) ! m2/g to m2/kg
        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
@@ -225,15 +225,15 @@ real(r8) function c_area( cohort_in )
     if (DEBUG_growth) then
        write(fates_log(),*) 'z_area 1',cohort_in%dbh,cohort_in%pft
        write(fates_log(),*) 'z_area 2',EDecophyscon%max_dbh
-       write(fates_log(),*) 'z_area 3',pftcon%woody
+       write(fates_log(),*) 'z_area 3',EDPftvarcon_inst%woody
        write(fates_log(),*) 'z_area 4',cohort_in%n
        write(fates_log(),*) 'z_area 5',cohort_in%patchptr%spread
        write(fates_log(),*) 'z_area 6',cohort_in%canopy_layer
        write(fates_log(),*) 'z_area 7',ED_val_grass_spread
     end if
 
     dbh = min(cohort_in%dbh,EDecophyscon%max_dbh(cohort_in%pft))
-    if(pftcon%woody(cohort_in%pft) == 1)then 
+    if(EDPftvarcon_inst%woody(cohort_in%pft) == 1)then 
        c_area = 3.142_r8 * cohort_in%n * &
             (cohort_in%patchptr%spread(cohort_in%canopy_layer)*dbh)**1.56_r8
     else

components/clm/src/ED/biogeochem/EDPatchDynamicsMod.F90

@@ -3,10 +3,9 @@ module EDPatchDynamicsMod
   ! ============================================================================
   ! Controls formation, creation, fusing and termination of patch level processes. 
   ! ============================================================================
-
   use FatesGlobals         , only : fates_log 
   use FatesInterfaceMod    , only : hlm_freq_day
-  use pftconMod            , only : pftcon
+  use EDPftvarcon          , only : EDPftvarcon_inst
   use EDCohortDynamicsMod  , only : fuse_cohorts, sort_cohorts, insert_cohort
   use EDtypesMod           , only : ncwd, n_dbh_bins, ntol, area, dbhmax
   use EDTypesMod           , only : numpft_ed
@@ -293,7 +292,7 @@ subroutine spawn_patches( currentSite )
                    nc%imort = nan
                 else
                    ! small trees 
-                   if(pftcon%woody(currentCohort%pft) == 1)then
+                   if(EDPftvarcon_inst%woody(currentCohort%pft) == 1)then
 
                       ! Number of trees in the understory of new patch, before we impose impact mortality and survivorship
                       nc%n = currentCohort%n * patch_site_areadis/currentPatch%area
@@ -578,7 +577,7 @@ subroutine fire_litter_fluxes(cp_target, new_patch_target, patch_site_areadis)
        currentCohort => currentPatch%shortest
        do while(associated(currentCohort))
           p = currentCohort%pft
-          if(pftcon%woody(p) == 1)then !DEAD (FROM FIRE) TREES
+          if(EDPftvarcon_inst%woody(p) == 1)then !DEAD (FROM FIRE) TREES
              !************************************/ 
              ! Number of trees that died because of the fire, per m2 of ground. 
              ! Divide their litter into the four litter streams, and spread evenly across ground surface. 
@@ -661,7 +660,7 @@ subroutine fire_litter_fluxes(cp_target, new_patch_target, patch_site_areadis)
        do while(associated(currentCohort))
 
           currentCohort%c_area = c_area(currentCohort) 
-          if(pftcon%woody(currentCohort%pft) == 1)then
+          if(EDPftvarcon_inst%woody(currentCohort%pft) == 1)then
              burned_leaves = (currentCohort%bl+currentCohort%bsw) * currentCohort%cfa
           else
              burned_leaves = (currentCohort%bl+currentCohort%bsw) * currentPatch%burnt_frac_litter(6)
@@ -738,7 +737,7 @@ subroutine mortality_litter_fluxes(cp_target, new_patch_target, patch_site_aread
                   canopy_dead*(currentCohort%br+currentCohort%bstore)
 
          else 
-             if(pftcon%woody(currentCohort%pft) == 1)then
+             if(EDPftvarcon_inst%woody(currentCohort%pft) == 1)then
 
                 understorey_dead = ED_val_understorey_death * currentCohort%n * (patch_site_areadis/currentPatch%area)  !kgC/site/day
                 currentPatch%canopy_mortality_woody_litter  = currentPatch%canopy_mortality_woody_litter  + &
@@ -1523,7 +1522,7 @@ subroutine set_root_fraction( cpatch , depth_gl )
     !  Calculates the fractions of the root biomass in each layer for each pft. 
     !
     ! !USES:
-    use pftconMod   , only : pftcon
+
     !
     ! !ARGUMENTS
     type(ed_patch_type),intent(inout), target :: cpatch
@@ -1540,10 +1539,10 @@ subroutine set_root_fraction( cpatch , depth_gl )
 
        do lev = 1, hlm_numlevsoil-1
           cpatch%rootfr_ft(ft,lev) = .5_r8*( &
-                 exp(-pftcon%roota_par(ft) * depth_gl(lev-1))  &
-               + exp(-pftcon%rootb_par(ft) * depth_gl(lev-1))  &
-               - exp(-pftcon%roota_par(ft) * depth_gl(lev))    &
-               - exp(-pftcon%rootb_par(ft) * depth_gl(lev)))
+                 exp(-EDPftvarcon_inst%roota_par(ft) * depth_gl(lev-1))  &
+               + exp(-EDPftvarcon_inst%rootb_par(ft) * depth_gl(lev-1))  &
+               - exp(-EDPftvarcon_inst%roota_par(ft) * depth_gl(lev))    &
+               - exp(-EDPftvarcon_inst%rootb_par(ft) * depth_gl(lev)))
        end do
     end do