Megan ctsm

bld/CLMBuildNamelist.pm

@@ -3257,6 +3257,7 @@ sub setup_logic_supplemental_nitrogen {
       add_default($opts, $nl_flags->{'inputdata_rootdir'}, $definition, $defaults, $nl,
 		  'suplnitro', 'use_fates'=>$nl_flags->{'use_fates'});
   }
+
   #
   # Error checking for suplnitro
   #
@@ -4112,6 +4113,7 @@ sub setup_logic_dust_emis {
          }
       }
   # Otherwise make sure dust settings are NOT being set in CLM
+
   } else {
       my @vars = ( "dust_emis_method", "zender_soil_erod_source" );
       foreach my $option ( @vars ) {
@@ -4154,10 +4156,6 @@ sub setup_logic_megan {
   if ( defined($nl->get_value('megan_specifier')) ||
          defined($nl->get_value('megan_factors_file')) ) {
     check_megan_spec( $opts, $nl, $definition );
-    if ( &value_is_true( $nl_flags->{'use_fates'} ) ) {
-      $log->warning("MEGAN can NOT be on when FATES is also on.\n" .
-                    "   Use the '-no-megan' option when '-bgc fates' is activated");
-    }
   }
 }
 

src/biogeochem/MEGANFactorsMod.F90

@@ -58,9 +58,9 @@ subroutine megan_factors_get( comp_name, factors, class_n, molecwght )
     !
     ! !ARGUMENTS:
     character(len=*),intent(in)  :: comp_name      ! MEGAN compound name
-    real(r8),        intent(out) :: factors(npfts) ! vegitation type factors for the compound of intrest
-    integer,         intent(out) :: class_n        ! MEGAN class number for the compound of intrest
-    real(r8),        intent(out) :: molecwght      ! molecular weight of the compound of intrest
+    real(r8),        intent(out) :: factors(npfts) ! vegetation type factors for the compound of interest
+    integer,         intent(out) :: class_n        ! MEGAN class number for the compound of interest
+    real(r8),        intent(out) :: molecwght      ! molecular weight of the compound of interest
     !
     ! LOCAL VARS:
     integer :: hashkey, ndx

src/biogeochem/VOCEmissionMod.F90

@@ -416,6 +416,7 @@ subroutine VOCEmission (bounds, num_soilp, filter_soilp, &
     !
     ! !USES:
     use subgridAveMod        , only : p2g
+    use clm_varctl           , only : use_fates
     !
     ! !ARGUMENTS:
     type(bounds_type)      , intent(in)    :: bounds                  
@@ -451,6 +452,7 @@ subroutine VOCEmission (bounds, num_soilp, filter_soilp, &
     real(r8) :: par240_sha              ! temporary
 
     integer                            :: class_num, n_meg_comps, imech, imeg, ii
+    integer                            :: patchpft ! to transfer FATES PFT space into CLM PFT space.
     character(len=16)                  :: mech_name
     type(shr_megan_megcomp_t), pointer :: meg_cmp
     real(r8)                           :: cp, alpha,  Eopt, topt  ! for history output
@@ -497,7 +499,7 @@ subroutine VOCEmission (bounds, num_soilp, filter_soilp, &
          fsun240       => canopystate_inst%fsun240_patch        , & ! Input:  [real(r8) (:)   ]  sunlit fraction of canopy last 240 hrs            
          elai          => canopystate_inst%elai_patch           , & ! Input:  [real(r8) (:)   ]  one-sided leaf area index with burying by snow
          elai240       => canopystate_inst%elai240_patch        , & ! Input:  [real(r8) (:)   ]  one-sided leaf area index with burying by snow last 240 hrs
-
+         ci_fates      => canopystate_inst%ci_patch            , & !Input:  [real(r8) (:)   ]  FATES-calculated internalleaf ci    
          cisun_z       => photosyns_inst%cisun_z_patch          , & ! Input:  [real(r8) (:,:) ]  sunlit intracellular CO2 (Pa)
          cisha_z       => photosyns_inst%cisha_z_patch          , & ! Input:  [real(r8) (:,:) ]  shaded intracellular CO2 (Pa)
 
@@ -581,33 +583,45 @@ subroutine VOCEmission (bounds, num_soilp, filter_soilp, &
 
              ! set emis factor
              ! if specified, set EF for isoprene with mapped values
+             if(use_fates)then
+                patchpft = canopystate_inst%voc_pftindex_patch(p)
+             else
+                patchpft = patch%itype(p)
+             endif
+
              if ( trim(meg_cmp%name) == 'isoprene' .and. shr_megan_mapped_emisfctrs) then
-                epsilon = get_map_EF(patch%itype(p),g, vocemis_inst)
+                epsilon = get_map_EF(patchpft,g, vocemis_inst)
              else
-                epsilon = meg_cmp%emis_factors(patch%itype(p))
+                epsilon = meg_cmp%emis_factors(patchpft)
              end if
 
+
              class_num = meg_cmp%class_number
 
              ! Activity factor for PPFD
              gamma_p = get_gamma_P(par_sun, par24_sun, par240_sun, par_sha, par24_sha, par240_sha, &
                   fsun(p), fsun240(p), forc_solad240(p),forc_solai240(p), LDF(class_num), cp, alpha)
-
+             
              ! Activity factor for T
              gamma_t = get_gamma_T(t_veg240(p), t_veg24(p),t_veg(p), ct1(class_num), ct2(class_num),&
-                                   betaT(class_num),LDF(class_num), Ceo(class_num), Eopt, topt, patch%itype(p))
+                                   betaT(class_num),LDF(class_num), Ceo(class_num), Eopt, topt, patchpft)
 
              ! Activity factor for Leaf Age
-             gamma_a = get_gamma_A(patch%itype(p), elai240(p),elai(p),class_num)
+             gamma_a = get_gamma_A(patchpft, elai240(p),elai(p),class_num)
 
              ! Activity factor for CO2 (only for isoprene)
              if (trim(meg_cmp%name) == 'isoprene') then 
                 co2_ppmv = 1.e6_r8*forc_pco2(g)/forc_pbot(c)
-                gamma_c = get_gamma_C(cisun_z(p,1),cisha_z(p,1),forc_pbot(c),fsun(p), co2_ppmv)
+                if(use_fates)then
+                   gamma_c = get_gamma_C(ci_fates(p),ci_fates(p),forc_pbot(c),fsun(p), co2_ppmv)
+                else
+                   gamma_c = get_gamma_C(cisun_z(p,1),cisha_z(p,1),forc_pbot(c),fsun(p), co2_ppmv)
+                endif
+
              else
                 gamma_c = 1._r8
              end if
-
+             
              ! Calculate total scaling factor
              gamma = gamma_l * gamma_sm * gamma_a * gamma_p * gamma_T * gamma_c
 
@@ -618,7 +632,6 @@ subroutine VOCEmission (bounds, num_soilp, filter_soilp, &
                 ! assign to arrays for history file output (not weighted by landfrac)
                 meg_out(imeg)%flux_out(p) = meg_out(imeg)%flux_out(p) &
                                           + epsilon * gamma * megemis_units_factor*1.e-3_r8 ! Kg/m2/sec
-
                 if (imeg==1) then 
                    ! 
                    gamma_out(p)=gamma
@@ -692,7 +705,6 @@ function get_map_EF(ivt_in, g_in, vocemis_inst)
     ! vocemis_inst%efisop_patch ! Output: [real(r8) (:,:)]  emission factors for isoprene for each patch [ug m-2 h-1]
 
     get_map_EF = 0._r8
-
     if (     ivt_in == ndllf_evr_tmp_tree  &
          .or.     ivt_in == ndllf_evr_brl_tree) then   !fineleaf evergreen
        get_map_EF = vocemis_inst%efisop_grc(2,g_in)
@@ -710,7 +722,7 @@ function get_map_EF(ivt_in, g_in, vocemis_inst)
     else if (ivt_in >= nc3crop) then                   !crops
        get_map_EF = vocemis_inst%efisop_grc(6,g_in)
     end if
-
+    
   end function get_map_EF
 
   !-----------------------------------------------------------------------

src/biogeophys/CanopyStateType.F90

@@ -46,12 +46,14 @@ module CanopyStateType
      real(r8) , pointer :: hbot_patch               (:)   ! patch canopy bottom (m)
      real(r8) , pointer :: z0m_patch                (:)   ! patch momentum roughness length (m)
      real(r8) , pointer :: displa_patch             (:)   ! patch displacement height (m)
+     real(r8) , pointer :: ci_patch                 (:)   ! Internal leaf CO2 concentration for MEGAN
      real(r8) , pointer :: fsun_patch               (:)   ! patch sunlit fraction of canopy
      real(r8) , pointer :: fsun24_patch             (:)   ! patch 24hr average of sunlit fraction of canopy
      real(r8) , pointer :: fsun240_patch            (:)   ! patch 240hr average of sunlit fraction of canopy
 
      real(r8) , pointer :: dleaf_patch              (:)   ! patch characteristic leaf width (diameter) [m]
-                                                          ! for non-ED/FATES this is the same as pftcon%dleaf()
+     ! for non-ED/FATES this is the same as pftcon%dleaf()
+     integer ,  pointer :: voc_pftindex_patch      (:)   ! FATES specific MEGAN pft index. 
      real(r8) , pointer :: rscanopy_patch           (:)   ! patch canopy stomatal resistance (s/m) (ED specific)
 
      real(r8) , pointer :: vegwp_patch              (:,:) ! patch vegetation water matric potential (mm)
@@ -139,11 +141,13 @@ subroutine InitAllocate(this, bounds)
     allocate(this%hbot_patch               (begp:endp))           ; this%hbot_patch               (:)   = nan
     allocate(this%z0m_patch                (begp:endp))           ; this%z0m_patch                (:)   = nan
     allocate(this%displa_patch             (begp:endp))           ; this%displa_patch             (:)   = nan
+    allocate(this%ci_patch                 (begp:endp))           ; this%ci_patch                 (:)   = nan
     allocate(this%fsun_patch               (begp:endp))           ; this%fsun_patch               (:)   = nan
     allocate(this%fsun24_patch             (begp:endp))           ; this%fsun24_patch             (:)   = nan
     allocate(this%fsun240_patch            (begp:endp))           ; this%fsun240_patch            (:)   = nan
 
     allocate(this%dleaf_patch              (begp:endp))           ; this%dleaf_patch              (:)   = nan
+    allocate(this%voc_pftindex_patch       (begp:endp))           ; this%voc_pftindex_patch       (:)   = 0    
     allocate(this%rscanopy_patch           (begp:endp))           ; this%rscanopy_patch           (:)   = nan
 !    allocate(this%gccanopy_patch           (begp:endp))           ; this%gccanopy_patch           (:)   = 0.0_r8
     allocate(this%vegwp_patch              (begp:endp,1:nvegwcs)) ; this%vegwp_patch              (:,:) = nan

src/utils/clmfates_interfaceMod.F90

@@ -1485,6 +1485,7 @@ subroutine wrap_update_hlmfates_dyn(this, nc, bounds_clump,      &
          z0m  => canopystate_inst%z0m_patch  , & ! Output: [real(r8) (:)   ] momentum roughness length (m)
          displa => canopystate_inst%displa_patch, &
          dleaf_patch => canopystate_inst%dleaf_patch, &
+         voc_pftindex => canopystate_inst%voc_pftindex_patch, &
          snow_depth => waterdiagnosticbulk_inst%snow_depth_col, &
          frac_sno_eff => waterdiagnosticbulk_inst%frac_sno_eff_col, &
          frac_veg_nosno_alb => canopystate_inst%frac_veg_nosno_alb_patch)
@@ -1608,7 +1609,7 @@ subroutine wrap_update_hlmfates_dyn(this, nc, bounds_clump,      &
           z0m(col%patchi(c)+1:col%patchf(c)) = 0.0_r8
           displa(col%patchi(c)+1:col%patchf(c)) = 0.0_r8
           dleaf_patch(col%patchi(c)+1:col%patchf(c)) = 0.0_r8
-
+          voc_pftindex(col%patchi(c)+1:col%patchf(c)) = 0
           frac_veg_nosno_alb(col%patchi(c):col%patchf(c)) = 0.0_r8
 
           ! Set the bareground patch indicator
@@ -1670,6 +1671,7 @@ subroutine wrap_update_hlmfates_dyn(this, nc, bounds_clump,      &
              z0m(p)    = this%fates(nc)%bc_out(s)%z0m_pa(ifp)
              displa(p) = this%fates(nc)%bc_out(s)%displa_pa(ifp)
              dleaf_patch(p) = this%fates(nc)%bc_out(s)%dleaf_pa(ifp)
+             voc_pftindex(p) = this%fates(nc)%bc_out(s)%nocomp_MEGAN_pft_label_pa(ifp)
           end do ! veg pach
 
           if(abs(areacheck - 1.0_r8).gt.1.e-9_r8)then
@@ -2575,8 +2577,8 @@ subroutine wrap_photosynthesis(this, nc, bounds, fn, filterp, &
           rssun     => photosyns_inst%rssun_patch  , &
           rssha     => photosyns_inst%rssha_patch,   &
           psnsun    => photosyns_inst%psnsun_patch,  &
-          psnsha    => photosyns_inst%psnsha_patch)
-
+          psnsha    => photosyns_inst%psnsha_patch,  &
+          ci        => canopystate_inst%ci_patch)
       do s = 1, this%fates(nc)%nsites
 
          c = this%f2hmap(nc)%fcolumn(s)
@@ -2643,7 +2645,7 @@ subroutine wrap_photosynthesis(this, nc, bounds, fn, filterp, &
             this%fates(nc)%bc_in(s)%filter_photo_pa(ifp) = 3
             rssun(p) = this%fates(nc)%bc_out(s)%rssun_pa(ifp)
             rssha(p) = this%fates(nc)%bc_out(s)%rssha_pa(ifp)
-
+            ci(p) = this%fates(nc)%bc_out(s)%ci_pa(ifp)
             ! These fields are marked with a bad-value flag
             photosyns_inst%psnsun_patch(p)   = spval
             photosyns_inst%psnsha_patch(p)   = spval