Megan fates pr follow up

bld/namelist_files/namelist_defaults_drydep.xml

@@ -25,6 +25,7 @@ attributes from the config_cache.xml file (with keys converted to upper-case).
 
 <megan_specifier>'ISOP = isoprene', 'C10H16 = pinene_a + carene_3 + thujene_a', 'CH3OH = methanol', 'C2H5OH = ethanol', 'CH2O = formaldehyde', 'CH3CHO = acetaldehyde', 'CH3COOH = acetic_acid', 'CH3COCH3 = acetone'</megan_specifier>
 
+<megan_factors_file phys="clm6_0" use_fates=".true." >atm/cam/chem/trop_mozart/emis/megan21_emis_factors_78pft_c20161108.nc</megan_factors_file>
 <megan_factors_file phys="clm6_0" >atm/cam/chem/trop_mozart/emis/megan21_emis_factors_78pft_c20161108.nc</megan_factors_file>
 <megan_factors_file phys="clm5_1" >atm/cam/chem/trop_mozart/emis/megan21_emis_factors_78pft_c20161108.nc</megan_factors_file>
 <megan_factors_file phys="clm5_0" >atm/cam/chem/trop_mozart/emis/megan21_emis_factors_78pft_c20161108.nc</megan_factors_file>

cime_config/testdefs/testmods_dirs/clm/FatesColdMeganSatPhen/user_nl_clm

@@ -2,3 +2,4 @@
 hist_mfilt        = 365
 hist_nhtfrq       = -24
 
+hist_fincl1 += 'GAMMA' 'GAMMAC' 'GAMMAA' 'GAMMAL' 'GAMMAT' 'GAMMAP' 'GAMMAS' 'CISHA' 'CISUN' 'FSUN24' 'FSUN' 'FSUN240' 'PCO2'

src/biogeochem/VOCEmissionMod.F90

@@ -92,7 +92,7 @@ subroutine Init(this, bounds)
 
     if ( shr_megan_mechcomps_n > 0) then
        if (use_fates) then
-          if (.not. use_fates_nocomp) then
+          if (.not. use_fates_nocomp) then ! SP implies NOCOMP is on.
              call endrun( msg='ERROR: MEGAN currently only works with when FATES is in SP and/or NOCOMP mode '//&
                   errMsg(sourcefile, __LINE__))
           end if
@@ -417,6 +417,8 @@ subroutine VOCEmission (bounds, num_soilp, filter_soilp, &
     ! !USES:
     use subgridAveMod        , only : p2g
     use clm_varctl           , only : use_fates
+    use clm_varcon           , only : smallValue
+    use decompMod            , only : subgrid_level_patch
     !
     ! !ARGUMENTS:
     type(bounds_type)      , intent(in)    :: bounds                  
@@ -452,7 +454,9 @@ 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.
+    integer                            :: l_pft_itype(bounds%begp:bounds%endp) ! local index of pft type 
+                                                      ! that corresponds to pfts on megan factors 
+                                                      ! for BGC it will be 1 to 1 with pftcon%itype(p)
     character(len=16)                  :: mech_name
     type(shr_megan_megcomp_t), pointer :: meg_cmp
     real(r8)                           :: cp, alpha,  Eopt, topt  ! for history output
@@ -499,7 +503,6 @@ 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)
 
@@ -531,11 +534,29 @@ subroutine VOCEmission (bounds, num_soilp, filter_soilp, &
     ! initialize variables which get passed to the atmosphere
     vocflx(bounds%begp:bounds%endp,:)   = 0._r8
     vocflx_tot(bounds%begp:bounds%endp) = 0._r8
-
+    
     do imeg=1,shr_megan_megcomps_n
       meg_out(imeg)%flux_out(bounds%begp:bounds%endp) = 0._r8
     enddo
-
+
+    ! Get local pft types:
+    ! this has to be done earlier, so if use_fates, we locally know what is not bare ground
+    ! voc_pft_index comes from fates-internal mapping between pft's in megan_factors_file and fates pfts
+    l_pft_itype(bounds%begp:bounds%endp) = 0
+    if (use_fates) then
+       do fp = 1,num_soilp
+          p = filter_soilp(fp)
+          if (patch%is_fates(p)) then
+            l_pft_itype(p) = canopystate_inst%voc_pftindex_patch(p)
+          endif
+       end do
+    else
+       do fp = 1,num_soilp
+          p = filter_soilp(fp)
+          l_pft_itype(p) = patch%itype(p)
+       end do
+    end if
+
     ! Begin loop over points
     !_______________________________________________________________________________
     do fp = 1,num_soilp
@@ -551,7 +572,7 @@ subroutine VOCEmission (bounds, num_soilp, filter_soilp, &
        vocflx_meg(:) = 0._r8
 
        ! calculate VOC emissions for non-bare ground Patches
-       if (patch%itype(p) > 0) then 
+       if (l_pft_itype(p) > 0) then 
           gamma=0._r8
 
           ! Calculate PAR: multiply w/m2 by 4.6 to get umol/m2/s for par (added 8/14/02)
@@ -583,16 +604,10 @@ 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(patchpft,g, vocemis_inst)
+                epsilon = get_map_EF(l_pft_itype(p),g, vocemis_inst)
              else
-                epsilon = meg_cmp%emis_factors(patchpft)
+                epsilon = meg_cmp%emis_factors(l_pft_itype(p))
              end if
 
 
@@ -604,19 +619,20 @@ subroutine VOCEmission (bounds, num_soilp, filter_soilp, &
 
              ! 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, patchpft)
+                                   betaT(class_num),LDF(class_num), Ceo(class_num), Eopt, topt, l_pft_itype(p))
 
              ! Activity factor for Leaf Age
-             gamma_a = get_gamma_A(patchpft, elai240(p),elai(p),class_num)
+             gamma_a = get_gamma_A(l_pft_itype(p), 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)
-                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)
+                ! Check of valid intercellular co2 pressure values.
+                if (cisha_z(p,1) < 0.0_r8 .or. cisun_z(p,1) < 0.0_r8) then
+                   write(iulog,*) 'Invalid intercellular co2 pressure (patch itype, sunlit, shaded): ', l_pft_itype(p),cisun_z(p,1),cisha_z(p,1)
+                   call endrun(subgrid_index=p, subgrid_level=subgrid_level_patch, msg=errMsg(sourcefile, __LINE__))
                 endif
+                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)
 
              else
                 gamma_c = 1._r8

src/biogeophys/CanopyFluxesMod.F90

@@ -1111,9 +1111,9 @@ subroutine CanopyFluxes(bounds,  num_exposedvegp, filter_exposedvegp,
             vpd(p)  = max((svpts(p) - eah(p)), 50._r8) * 0.001_r8
 
          end do
-
          if ( use_fates ) then      
-
+            ! zeroed out here for megan isoprene calculations, in case patches are not transfered properly.
+
             call clm_fates%wrap_photosynthesis(nc, bounds, fn, filterp(1:fn), &
                  svpts(begp:endp), eah(begp:endp), o2(begp:endp), &
                  co2(begp:endp), rb(begp:endp), dayl_factor(begp:endp), &

src/biogeophys/CanopyStateType.F90

@@ -46,7 +46,6 @@ 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
@@ -141,7 +140,6 @@ 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

src/biogeophys/PhotosynthesisMod.F90

@@ -709,6 +709,30 @@ subroutine InitHistory(this, bounds)
 
     endif
 
+    ! megan variables
+    this%cisun_z_patch(begp:endp,:) = spval
+    this%cisha_z_patch(begp:endp,:) = spval
+    if (nlevcan>1) then 
+       call hist_addfld2d (fname='CISUN', units='(Pa)', type2d='nlevcan', &
+          avgflag='A', long_name='sunlit leaf intracellular CO2 pressure', &
+          ptr_patch=this%cisun_z_patch, set_lake=spval, set_urb=spval, default='inactive')
+
+       call hist_addfld2d (fname='CISHA', units='(Pa)', type2d='nlevcan', &
+          avgflag='A', long_name='shaded leaf intracellular CO2 pressure', &
+          ptr_patch=this%cisha_z_patch, set_lake=spval, set_urb=spval, default='inactive')
+    else
+       ptr_1d => this%cisun_z_patch(begp:endp,1)
+       call hist_addfld1d (fname='CISUN', units='(Pa)', &
+          avgflag='A', long_name='sunlit leaf intracellular CO2 pressure', &
+          ptr_patch=ptr_1d, default='inactive')
+
+       ptr_1d => this%cisha_z_patch(begp:endp,1)
+       call hist_addfld1d (fname='CISHA', units='(Pa)', &
+          avgflag='A', long_name='shaded leaf intracellular CO2 pressure', &
+          ptr_patch=ptr_1d, default='inactive')
+
+    endif
+
   end subroutine InitHistory
 
   !-----------------------------------------------------------------------

src/utils/clmfates_interfaceMod.F90

@@ -1616,6 +1616,10 @@ subroutine wrap_update_hlmfates_dyn(this, nc, bounds_clump,      &
           patch%is_bareground(col%patchi(c)) = .true.
           npatch = this%fates(nc)%sites(s)%youngest_patch%patchno
 
+          ! set voc_pft_index of bareground to 0 explicitly, so the bare ground is properly ignored in VOCEmissionMod
+          if (patch%is_bareground(col%patchi(c))) then
+             voc_pftindex(col%patchi(c)) = 0
+          endif
           ! Precision errors on the canopy_fraction_pa sum, even small (e-12)
           ! do exist, and can create potentially negetive bare-soil fractions
           ! (ie -1e-12 or smaller). Even though this is effectively zero,
@@ -2578,7 +2582,9 @@ subroutine wrap_photosynthesis(this, nc, bounds, fn, filterp, &
           rssha     => photosyns_inst%rssha_patch,   &
           psnsun    => photosyns_inst%psnsun_patch,  &
           psnsha    => photosyns_inst%psnsha_patch,  &
-          ci        => canopystate_inst%ci_patch)
+          cisun_z   => photosyns_inst%cisun_z_patch, & 
+          cisha_z   => photosyns_inst%cisha_z_patch & 
+         )
       do s = 1, this%fates(nc)%nsites
 
          c = this%f2hmap(nc)%fcolumn(s)
@@ -2645,7 +2651,14 @@ 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)
+            ! this is needed for MEGAN to work with FATES
+            cisun_z(p,:) = this%fates(nc)%bc_out(s)%ci_pa(ifp)
+            cisha_z(p,:) = this%fates(nc)%bc_out(s)%ci_pa(ifp)
+            if (this%fates(nc)%bc_out(s)%ci_pa(ifp) <0.0_r8) then
+               write(iulog,*) 'ci_pa is less than 0: ', this%fates(nc)%bc_out(s)%ci_pa(ifp)
+               write(iulog,*) 'filter ran photosynthesis s p icp ifp ilter',s,p,icp,ifp
+               call endrun(msg=errMsg(sourcefile, __LINE__))
+            endif
             ! These fields are marked with a bad-value flag
             photosyns_inst%psnsun_patch(p)   = spval
             photosyns_inst%psnsha_patch(p)   = spval