Merge pull request #10 from rosiealice/megan_fates

Megan fates update

biogeochem/EDCanopyStructureMod.F90

@@ -1977,6 +1977,7 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out)
                 bc_out(s)%z0m_pa(ifp)    = EDPftvarcon_inst%z0mr(1) * bc_out(s)%htop_pa(ifp)
                 bc_out(s)%displa_pa(ifp) = EDPftvarcon_inst%displar(1) * bc_out(s)%htop_pa(ifp)
                 bc_out(s)%dleaf_pa(ifp)  = EDPftvarcon_inst%dleaf(1)
+                bc_out(s)%nocomp_MEGAN_pft_label_pa(ifp) = 1
              endif
              ! -----------------------------------------------------------------------------
 
@@ -2002,6 +2003,12 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out)
 
              bc_out(s)%nocomp_pft_label_pa(ifp) = currentPatch%nocomp_pft_label
 
+             if(currentPatch%nocomp_pft_label.gt.0)then
+                bc_out(s)%nocomp_MEGAN_pft_label_pa(ifp) = EDPftvarcon_inst%voc_pftindex(currentPatch%nocomp_pft_label)
+             else
+                bc_out(s)%nocomp_MEGAN_pft_label_pa(ifp) = 1 ! dummy for bare ground. 
+             endif
+
              ! Calculate area indices for output boundary to HLM
              ! It is assumed that cpatch%canopy_area_profile and cpat%xai_profiles
              ! have been updated (ie ed_leaf_area_profile has been called since dynamics has been called)

biogeophys/FatesPlantRespPhotosynthMod.F90

@@ -23,6 +23,7 @@ module FATESPlantRespPhotosynthMod
   use FatesGlobals,      only : endrun => fates_endrun
   use FatesGlobals,      only : fates_log
   use FatesGlobals,      only : FatesWarn,N2S,A2S,I2S
+  use FatesInterfaceTypesMod        , only : hlm_use_nocomp
   use FatesConstantsMod, only : r8 => fates_r8
   use FatesConstantsMod, only : itrue
   use FatesConstantsMod, only : nearzero
@@ -193,6 +194,9 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime)
     ! net leaf photosynthesis averaged over sun and shade leaves. [umol CO2/m**2/s]
     real(r8) :: anet_av_z(nlevleaf,maxpft,nclmax)
 
+    ! internal leaf Co2 memory variable (for passing into MEGAN) Pa
+     real(r8) :: internal_co2_z(nlevleaf,maxpft,nclmax)
+
     ! Photosynthesis [umol /m2 /s]
     real(r8) :: psn_z(nlevleaf,maxpft,nclmax)
 
@@ -355,9 +359,10 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime)
                ! ---------------------------------------------------------------------------
                bc_out(s)%rssun_pa(ifp)     = 0._r8
                bc_out(s)%rssha_pa(ifp)     = 0._r8
-
-               psn_z(:,:,:) = 0._r8
+               bc_out(s)%ci_pa(ifp)     = 0._r8
 
+               psn_z(:,:,:) = 0._r8
+               !internal_co2_z(:,:,:) = 0._r8
                g_sb_leaves = 0._r8
                patch_la    = 0._r8
 
@@ -731,6 +736,7 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime)
                                       psn_z(iv,ft,cl),                    &  ! out
                                       rs_z(iv,ft,cl),                     &  ! out
                                       anet_av_z(iv,ft,cl),                &  ! out
+                                      internal_co2_z(iv,ft,cl),           &  ! out
                                       c13disc_z(cl,ft,iv))                   ! out
 
                                  rate_mask_z(iv,ft,cl) = .true.
@@ -1080,7 +1086,13 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime)
                      ! when it comes time to calculate a flux rate per unit ground
                      bc_out(s)%rssun_pa(ifp) = r_stomata
                      bc_out(s)%rssha_pa(ifp) = r_stomata
-
+
+                     if(hlm_use_nocomp)then
+                        bc_out(s)%ci_pa(ifp) = internal_co2_z(1,currentpatch%nocomp_pft_label,1)
+                     else
+                        bc_out(s)%ci_pa(ifp) = -999
+                     endif 
+
                      ! This value is used for diagnostics, the molar form of conductance
                      ! is what is used in the field usually, so we track that form
                      currentPatch%c_stomata  = cf / r_stomata
@@ -1204,6 +1216,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl,         &  ! in
      psn_out,           &  ! out
      rstoma_out,        &  ! out
      anet_av_out,       &  ! out
+     internal_co2_out,   & ! out
      c13disc_z)            ! out
 
 
@@ -1263,6 +1276,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl,         &  ! in
   real(r8), intent(out) :: rstoma_out     ! stomatal resistance (1/gs_lsl) (s/m)
   real(r8), intent(out) :: anet_av_out    ! net leaf photosynthesis (umol CO2/m**2/s)
   ! averaged over sun and shade leaves.
+  real(r8), intent(out) :: internal_co2_out    ! internal co2 diagnostic (Pa) 
   real(r8), intent(out) :: c13disc_z      ! carbon 13 in newly assimilated carbon
 
 
@@ -1354,6 +1368,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl,         &  ! in
 
      anet_av_out = -lmr
      psn_out     = 0._r8
+     internal_co2_out =  init_co2_inter_c
 
      ! The cuticular conductance already factored in maximum resistance as a bound
      ! no need to re-bound it
@@ -1372,6 +1387,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl,         &  ! in
         psn_out     = 0._r8    ! psn is accumulated across sun and shaded leaves.
         rstoma_out  = 0._r8    ! 1/rs is accumulated across sun and shaded leaves.
         anet_av_out = 0._r8
+        internal_co2_out = 0._r8
         gstoma  = 0._r8
 
         do  sunsha = 1,2
@@ -1599,10 +1615,12 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl,         &  ! in
            if(sunsha == 1)then !sunlit
               psn_out     = psn_out + agross * f_sun_lsl
               anet_av_out = anet_av_out + anet * f_sun_lsl
+              internal_co2_out = internal_co2_out + co2_inter_c * f_sun_lsl
               gstoma  = gstoma + 1._r8/(min(1._r8/gs, rsmax0)) * f_sun_lsl
            else
               psn_out = psn_out + agross * (1.0_r8-f_sun_lsl)
               anet_av_out = anet_av_out + anet * (1.0_r8-f_sun_lsl)
+              internal_co2_out = internal_co2_out + co2_inter_c * (1.0_r8 -f_sun_lsl)              
               gstoma  = gstoma + &
                    1._r8/(min(1._r8/gs, rsmax0)) * (1.0_r8-f_sun_lsl)
            end if
@@ -1648,7 +1666,7 @@ subroutine LeafLayerPhotosynthesis(f_sun_lsl,         &  ! in
 
         psn_out     = 0._r8
         anet_av_out = 0._r8
-
+        internal_co2_out = 0._r8
         rstoma_out  = min(rsmax0,cf/(stem_cuticle_loss_frac*stomatal_intercept(ft)))
         c13disc_z = 0.0_r8
 

main/EDPftvarcon.F90

@@ -147,7 +147,8 @@ module EDPftvarcon
      real(r8), allocatable :: rhos(:, :)                 ! Stem reflectance; second dim: 1 = vis, 2 = nir
      real(r8), allocatable :: taul(:, :)                 ! Leaf transmittance; second dim: 1 = vis, 2 = nir
      real(r8), allocatable :: taus(:, :)                 ! Stem transmittance; second dim: 1 = vis, 2 = nir
-
+     real(r8),  allocatable :: voc_pftindex(:)            ! Index for MEGAN parameters 
+
      ! Fire Parameters (No PFT vector capabilities in their own routines)
      ! See fire/SFParamsMod.F90 for bulk of fire parameters
      ! -------------------------------------------------------------------------------------------
@@ -454,6 +455,10 @@ subroutine Register_PFT(this, fates_params)
     call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
          dimension_names=dim_names, lower_bounds=dim_lower_bound)
 
+    name = 'fates_voc_pftindex'
+    call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
+         dimension_names=dim_names, lower_bounds=dim_lower_bound)
+
     name = 'fates_nonhydro_smpso'
     call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
          dimension_names=dim_names, lower_bounds=dim_lower_bound)
@@ -921,6 +926,10 @@ subroutine Receive_PFT(this, fates_params)
     call fates_params%RetrieveParameterAllocate(name=name, &
          data=this%c3psn)
 
+    name = 'fates_voc_pftindex'
+    call fates_params%RetrieveParameterAllocate(name=name, &
+         data=this%voc_pftindex)    
+
     name = 'fates_nonhydro_smpso'
     call fates_params%RetrieveParameterAllocate(name=name, &
          data=this%smpso)
@@ -1728,6 +1737,7 @@ subroutine FatesReportPFTParams(is_master)
         write(fates_log(),fmt0) 'xl = ',EDPftvarcon_inst%xl
         write(fates_log(),fmt0) 'clumping_index = ',EDPftvarcon_inst%clumping_index
         write(fates_log(),fmt0) 'c3psn = ',EDPftvarcon_inst%c3psn
+        write(fates_log(),fmt0) 'voc_pftindex = ',EDPftvarcon_inst%voc_pftindex        
         write(fates_log(),fmt0) 'vcmax25top = ',EDPftvarcon_inst%vcmax25top
         write(fates_log(),fmt0) 'smpso = ',EDPftvarcon_inst%smpso
         write(fates_log(),fmt0) 'smpsc = ',EDPftvarcon_inst%smpsc
@@ -2229,7 +2239,12 @@ subroutine FatesCheckParams(is_master)
            call endrun(msg=errMsg(sourcefile, __LINE__))
 
         end if
-
+        write(*,*)  EDPftvarcon_inst%voc_pftindex(:)
+        if(EDPftvarcon_inst%voc_pftindex(ipft) .le. 0 .or. EDPftvarcon_inst%voc_pftindex(ipft) .gt. 16 ) then
+           write(fates_log(),*) 'MEGAN indices must be between 1 and 16',ipft,EDPftvarcon_inst%voc_pftindex(ipft)
+           call endrun(msg=errMsg(sourcefile, __LINE__))
+       endif
+
         if( hlm_use_fixed_biogeog .eq. itrue ) then
            ! check that the host-fates PFT map adds to one along HLM dimension so that all the HLM area
            ! goes to a FATES PFT.  Each FATES PFT can get < or > 1 of an HLM PFT.

main/FatesInterfaceMod.F90

@@ -325,6 +325,7 @@ subroutine zero_bcs(fates,s)
     ! Output boundaries
     fates%bc_out(s)%active_suction_sl(:) = .false.
     fates%bc_out(s)%fsun_pa(:)      = 0.0_r8
+    fates%bc_out(s)%ci_pa(:)        = 0.0_r8    
     fates%bc_out(s)%laisun_pa(:)    = 0.0_r8
     fates%bc_out(s)%laisha_pa(:)    = 0.0_r8
     fates%bc_out(s)%rootr_pasl(:,:) = 0.0_r8
@@ -391,6 +392,8 @@ subroutine zero_bcs(fates,s)
     fates%bc_out(s)%z0m_pa(:)    = 0.0_r8
     fates%bc_out(s)%dleaf_pa(:)   = 0.0_r8
     fates%bc_out(s)%nocomp_pft_label_pa(:) = 0
+    fates%bc_out(s)%nocomp_MEGAN_pft_label_pa(:) = 0
+
 
     fates%bc_out(s)%canopy_fraction_pa(:) = 0.0_r8
     fates%bc_out(s)%frac_veg_nosno_alb_pa(:) = 0.0_r8
@@ -605,6 +608,7 @@ subroutine allocate_bcout(bc_out, nlevsoil_in, nlevdecomp_in)
 
       ! Radiation
       allocate(bc_out%fsun_pa(maxpatch_total))
+      allocate(bc_out%ci_pa(maxpatch_total))
       allocate(bc_out%laisun_pa(maxpatch_total))
       allocate(bc_out%laisha_pa(maxpatch_total))
 
@@ -718,6 +722,7 @@ subroutine allocate_bcout(bc_out, nlevsoil_in, nlevdecomp_in)
       allocate(bc_out%frac_veg_nosno_alb_pa(maxpatch_total))
 
       allocate(bc_out%nocomp_pft_label_pa(maxpatch_total))
+      allocate(bc_out%nocomp_MEGAN_pft_label_pa(maxpatch_total))      
 
       ! Plant-Hydro BC's
       if (hlm_use_planthydro.eq.itrue) then

main/FatesInterfaceTypesMod.F90

@@ -600,6 +600,10 @@ module FatesInterfaceTypesMod
 
       ! Shaded canopy LAI
       real(r8),allocatable :: laisha_pa(:)
+
+      ! Average internal CO2 concentration
+      real(r8), allocatable :: ci_pa(:)
+
 
       ! Logical stating whether a soil layer can have water uptake by plants
       ! The only condition right now is that liquid water exists
@@ -750,6 +754,8 @@ module FatesInterfaceTypesMod
 
      integer, allocatable :: nocomp_pft_label_pa(:) ! in nocomp and SP mode, each patch has a PFT identity. 
 
+     integer, allocatable :: nocomp_MEGAN_pft_label_pa(:) ! Index to map from FATES NOCOMP PFT identity into MEGAN PFT space. 
+
       ! FATES Hydraulics
 
 

parameter_files/fates_params_default.cdl

@@ -672,6 +672,9 @@ variables:
 	double fates_woody(fates_pft) ;
 		fates_woody:units = "logical flag" ;
 		fates_woody:long_name = "Binary woody lifeform flag" ;
+	double fates_voc_pftindex(fates_pft) ;
+	       fates_voc_pftindex: units = "Integer PFT index" ;
+	       fates_voc_pftindex: long_name = "integer index for MEGAN PFT definitions" ;
 	double fates_hlm_pft_map(fates_hlm_pftno, fates_pft) ;
 		fates_hlm_pft_map:units = "area fraction" ;
 		fates_hlm_pft_map:long_name = "In fixed biogeog mode, fraction of HLM area associated with each FATES PFT" ;
@@ -1669,7 +1672,9 @@ data:
  fates_wood_density = 0.548327, 0.44235, 0.454845, 0.754336, 0.548327, 
     0.566452, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7, 0.7 ;
 
- fates_woody = 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0 ;
+ fates_woody =        1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0 ;
+
+ fates_voc_pftindex = 4, 1, 3, 5, 6, 7, 9, 10, 10, 9,  11, 12, 13, 14 ;
 
  fates_hlm_pft_map =
   0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,