Rebased drydep fates

biogeochem/EDCanopyStructureMod.F90

@@ -1877,7 +1877,7 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out)
     use EDTypesMod        , only : ed_site_type, AREA
     use FatesPatchMod,      only : fates_patch_type
     use FatesInterfaceTypesMod , only : bc_out_type
-
+     use FatesInterfaceTypesMod    , only : hlm_day_of_year
     !
     ! !ARGUMENTS
     integer,            intent(in)            :: nsites
@@ -1894,7 +1894,7 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out)
     real(r8) :: total_canopy_area
     real(r8) :: total_patch_leaf_stem_area
     real(r8) :: weight  ! Weighting for cohort variables in patch
-    
+
     do s = 1,nsites
 
        ifp = 0
@@ -1904,7 +1904,8 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out)
        bc_out(s)%dleaf_pa(:) = 0._r8
        bc_out(s)%z0m_pa(:) = 0._r8
        bc_out(s)%displa_pa(:) = 0._r8
-
+       bc_out(s)%drydep_season_pa(:) = 0
+       bc_out(s)%wesley_pft_label_pa(:)=8 !for no vegetation. 
        currentPatch => sites(s)%oldest_patch
        c = fcolumn(s)
        do while(associated(currentPatch))
@@ -1979,6 +1980,46 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out)
                 bc_out(s)%dleaf_pa(ifp)  = EDPftvarcon_inst%dleaf(1)
                 bc_out(s)%nocomp_MEGAN_pft_label_pa(ifp) = 1
              endif
+
+             ! Set patch properties for the HLM dry depositioin model
+             ! This wants to know the PFT according to the 'wesley' classification scheme
+             ! and the season according to the same scheme. 
+             bc_out(s)%wesley_pft_label_pa(ifp) = EDPftvarcon_inst%wesley_pft_index_fordrydep(currentPatch%nocomp_pft_label)
+
+            !wesely seasonal "index_season"                                
+            ! 1 - midsummer with lush vegetation                           
+            ! 2 - Autumn with unharvested cropland                         
+            ! 3 - Late autumn after frost, no snow                         
+            ! 4 - Winter, snow on ground and subfreezing                   
+            ! 5 - Transitional spring with partially green short annuals
+             if(bc_out(s)%tlai_pa(ifp) .gt. 2.0_r8)then
+                bc_out(s)%drydep_season_pa(ifp) = 1 ! Summer, or something like it.
+             else ! NOT SUMMER
+                if(sites(s)%lat>0)then ! Northern HS
+                   if(hlm_day_of_year .lt. 180)then ! DOY
+                      bc_out(s)%drydep_season_pa(ifp) = 5 ! NH spring
+                   else ! autumn
+      	      	      if(bc_out(s)%tlai_pa(ifp) .gt. 1.0_r8)then
+                         bc_out(s)%drydep_season_pa(ifp) = 2 ! NH early autumn
+                      else
+                         bc_out(s)%drydep_season_pa(ifp) = 3 ! NH late autumn
+                      endif
+                   endif ! DOY
+                 else !Southern HS 
+                    if(hlm_day_of_year .gt. 180)then ! spring
+                       bc_out(s)%drydep_season_pa(ifp) = 5 ! SH spring
+                    else ! SH autumn
+                       if(bc_out(s)%tlai_pa(ifp) .gt. 1.0_r8)then
+                          bc_out(s)%drydep_season_pa(ifp) = 2 ! SH early autumn
+                       else 
+                          bc_out(s)%drydep_season_pa(ifp) = 3 ! SH late autumn
+                       endif ! autumn
+                    endif ! DOY
+                 endif ! Hemisphere
+              endif ! summer?
+
+
+
              ! -----------------------------------------------------------------------------
 
              ! We are assuming here that grass is all located underneath tree canopies.
@@ -2001,7 +2042,6 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out)
 
              total_canopy_area = total_canopy_area + bc_out(s)%canopy_fraction_pa(ifp)
 
-             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)
@@ -2061,8 +2101,9 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out)
              if(currentPatch%nocomp_pft_label.ne.nocomp_bareground)then ! for vegetated patches only
                 ifp = ifp+1
                 bc_out(s)%canopy_fraction_pa(ifp) = bc_out(s)%canopy_fraction_pa(ifp)/total_patch_area
-             endif ! veg patch
-
+                bc_out(s)%wesley_pft_label_pa(ifp)=8
+                bc_out(s)%drydep_season_pa(ifp)=4 !winter season where there is bare ground
+              endif ! veg patch
              currentPatch => currentPatch%younger
           end do
 
@@ -2082,7 +2123,6 @@ subroutine update_hlm_dynamics(nsites,sites,fcolumn,bc_out)
 
        ! Pass FATES Harvested C to bc_out.
        call UpdateHarvestC(sites(s),bc_out(s))
-
     end do
 
     ! This call to RecruitWaterStorage() makes an accounting of

biogeochem/FatesPatchMod.F90

@@ -70,6 +70,8 @@ module FatesPatchMod
     real(r8) :: age_since_anthro_disturbance ! average age for secondary forest since last anthropogenic disturbance [years]
     logical  :: changed_landuse_this_ts      ! logical flag to track patches that have just undergone land use change [only used with nocomp and land use change]
 
+
+
     !---------------------------------------------------------------------------
 
     ! RUNNING MEANS
@@ -755,7 +757,7 @@ subroutine Create(this, age, area, land_use_label, nocomp_pft, num_swb, num_pft,
         this%age_since_anthro_disturbance = fates_unset_r8
       endif
       this%nocomp_pft_label = nocomp_pft
-
+      
       this%tr_soil_dir(:) = 1.0_r8
       this%tr_soil_dif(:) = 1.0_r8
       this%NCL_p          = 1

main/EDPftvarcon.F90

@@ -287,6 +287,9 @@ module EDPftvarcon
      real(r8), allocatable :: landusechange_frac_exported(:)  ! fraction of land use change-generated wood material that is exported to wood product (the remainder is either burned or goes to litter)
      real(r8), allocatable :: landusechange_pprod10(:)        ! fraction of land use change wood product that goes to 10-year product pool (remainder goes to 100-year pool)
 
+     real(r8), allocatable :: wesley_pft_index_fordrydep(:)
+
+
    contains
      procedure, public :: Init => EDpftconInit
      procedure, public :: Register
@@ -828,6 +831,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_wesley_pft_index_fordrydep'
+    call fates_params%RegisterParameter(name=name, dimension_shape=dimension_shape_1d, &
+          dimension_names=dim_names, lower_bounds=dim_lower_bound)
+
     ! adding the hlm_pft_map variable with two dimensions - FATES PFTno and HLM PFTno
     pftmap_dim_names(1) = dimension_name_pft
     pftmap_dim_names(2) = dimension_name_hlm_pftno
@@ -1292,6 +1299,10 @@ subroutine Receive_PFT(this, fates_params)
     call fates_params%RetrieveParameterAllocate(name=name, &
          data=this%eca_lambda_ptase)
 
+    name = 'fates_wesley_pft_index_fordrydep'
+    call fates_params%RetrieveParameterAllocate(name=name, &
+         data=this%wesley_pft_index_fordrydep)
+
     name = 'fates_hlm_pft_map'
     call fates_params%RetrieveParameterAllocate(name=name, &
          data=this%hlm_pft_map)
@@ -1805,6 +1816,7 @@ subroutine FatesReportPFTParams(is_master)
         write(fates_log(),fmt0) 'hydro_pinot_node = ',EDPftvarcon_inst%hydr_pinot_node
         write(fates_log(),fmt0) 'hydro_kmax_node = ',EDPftvarcon_inst%hydr_kmax_node
         write(fates_log(),fmt0) 'hlm_pft_map = ', EDPftvarcon_inst%hlm_pft_map
+        write(fates_log(),fmt0) 'wesley_pft_index_fordrydep = ', EDPftvarcon_inst%wesley_pft_index_fordrydep
         write(fates_log(),fmt0) 'hydro_vg_alpha_node  = ',EDPftvarcon_inst%hydr_vg_alpha_node
         write(fates_log(),fmt0) 'hydro_vg_m_node  = ',EDPftvarcon_inst%hydr_vg_m_node
         write(fates_log(),fmt0) 'hydro_vg_n_node  = ',EDPftvarcon_inst%hydr_vg_n_node

main/FatesInterfaceMod.F90

@@ -393,8 +393,8 @@ subroutine zero_bcs(fates,s)
     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)%wesley_pft_label_pa(:) = 0
+    fates%bc_out(s)%drydep_season_pa(:) = 0
     fates%bc_out(s)%canopy_fraction_pa(:) = 0.0_r8
     fates%bc_out(s)%frac_veg_nosno_alb_pa(:) = 0.0_r8
 
@@ -723,6 +723,8 @@ subroutine allocate_bcout(bc_out, nlevsoil_in, nlevdecomp_in)
 
       allocate(bc_out%nocomp_pft_label_pa(maxpatch_total))
       allocate(bc_out%nocomp_MEGAN_pft_label_pa(maxpatch_total))      
+      allocate(bc_out%wesley_pft_label_pa(maxpatch_total))
+      allocate(bc_out%drydep_season_pa(maxpatch_total))
 
       ! Plant-Hydro BC's
       if (hlm_use_planthydro.eq.itrue) then

main/FatesInterfaceTypesMod.F90

@@ -753,6 +753,8 @@ module FatesInterfaceTypesMod
                                                         ! [0,1]
 
      integer, allocatable :: nocomp_pft_label_pa(:) ! in nocomp and SP mode, each patch has a PFT identity. 
+     integer, allocatable :: wesley_pft_label_pa(:) ! For dry deposition, each FATES PFT needs to correspond to a PFT from the 'wesley 1989' scheme 
+     integer, allocatable :: drydep_season_pa(:) ! For dry deposition, we need to define the season index, from 1-5, for the purposes of detrmining the deposition velocity parameters. See drydep code for details of season indices. 
 
      integer, allocatable :: nocomp_MEGAN_pft_label_pa(:) ! Index to map from FATES NOCOMP PFT identity into MEGAN PFT space. 
 

parameter_files/fates_params_default.cdl

@@ -678,6 +678,9 @@ variables:
 	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" ;
+        double fates_wesley_pft_index_fordrydep(fates_pft) ; 
+	        fates_wesley_pft_index_fordrydep:units = "integer index";
+		fates_wesley_pft_index_fordrydep:long_name = "index to map from fates pfts into the wesley pft space used to determing dry deposition veolcity." ;
 	double fates_fire_FBD(fates_litterclass) ;
 		fates_fire_FBD:units = "kg Biomass/m3" ;
 		fates_fire_FBD:long_name = "fuel bulk density" ;
@@ -1368,16 +1371,16 @@ data:
  fates_leaf_slatop = 0.012, 0.005, 0.024, 0.009, 0.03, 0.03, 0.012, 0.03, 
     0.03, 0.01, 0.032, 0.027, 0.05, 0.05 ;
 
- fates_leaf_stomatal_intercept = 10000, 10000, 10000, 10000, 10000, 10000, 
+ fates_leaf_stomatal_intercept = 50000, 10000, 10000, 10000, 10000, 10000, 
     10000, 10000, 10000, 10000, 10000, 10000, 10000, 40000 ;
 
- fates_leaf_stomatal_slope_ballberry = 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8 ;
+ fates_leaf_stomatal_slope_ballberry = 12, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8 ;
 
  fates_leaf_stomatal_slope_medlyn = 4.1, 2.3, 2.3, 4.1, 4.4, 4.4, 4.7, 4.7, 
     4.7, 4.7, 4.7, 2.2, 5.3, 1.6 ;
 
  fates_leaf_vcmax25top =
-  50, 62, 39, 61, 58, 58, 62, 54, 54, 38, 54, 86, 78, 78 ;
+  60, 62, 39, 61, 58, 58, 62, 54, 54, 38, 54, 86, 78, 78 ;
 
  fates_leaf_vcmaxha = 65330, 65330, 65330, 65330, 65330, 65330, 65330, 65330, 
     65330, 65330, 65330, 65330, 65330, 65330 ;
@@ -1653,8 +1656,8 @@ data:
  fates_turb_leaf_diameter = 0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 0.04, 
     0.04, 0.04, 0.04, 0.04, 0.04, 0.04 ;
 
- fates_turb_z0mr = 0.055, 0.055, 0.055, 0.055, 0.055, 0.055, 0.055, 0.055, 
-    0.055, 0.055, 0.055, 0.055, 0.055, 0.055 ;
+ fates_turb_z0mr = 0.075, 0.055, 0.055, 0.075, 0.055, 0.055, 0.12, 0.12,
+     0.12, 0.12, 0.12, 0.12, 0.12, 0.12 ;
 
  fates_turnover_branch = 150, 150, 150, 150, 150, 150, 150, 150, 150, 150, 
     150, 0, 0, 0 ;
@@ -1692,6 +1695,8 @@ data:
   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,
   0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1 ;
 
+ fates_wesley_pft_index_fordrydep =  4, 5, 5, 4, 4, 4, 11, 11, 11 ,3, 3, 3, 3, 3;
+
  fates_fire_FBD = 15.4, 16.8, 19.6, 999, 4, 4 ;
 
  fates_fire_low_moisture_Coeff = 1.12, 1.09, 0.98, 0.8, 1.15, 1.15 ;