Merge pull request #632 from rosiealice/fates_sp

Fates Fixed Biogeography

Author: glemieux

biogeochem/EDPhysiologyMod.F90

@@ -1205,6 +1205,7 @@ subroutine SeedIn( currentSite, bc_in )
     ! !USES:
     use EDTypesMod, only : area
     use EDTypesMod, only : homogenize_seed_pfts
+    !use FatesInterfaceTypesMod,  only : hlm_use_fixed_biogeog    ! For future reduced complexity?
     !
     ! !ARGUMENTS    
     type(ed_site_type), intent(inout), target  :: currentSite
@@ -1308,18 +1309,17 @@ subroutine SeedIn( currentSite, bc_in )
 
           litt => currentPatch%litter(el)
           do pft = 1,numpft
-             
+             if(currentSite%use_this_pft(pft).eq.itrue)then
              ! Seed input from local sources (within site)
              litt%seed_in_local(pft) = litt%seed_in_local(pft) + site_seed_rain(pft)/area
 
              ! Seed input from external sources (user param seed rain, or dispersal model)
              seed_in_external =  seed_stoich*EDPftvarcon_inst%seed_suppl(pft)*years_per_day
-             
              litt%seed_in_extern(pft) = litt%seed_in_extern(pft) + seed_in_external
 
              ! Seeds entering externally [kg/site/day]
              site_mass%seed_in = site_mass%seed_in + seed_in_external*currentPatch%area
-
+            end if !use this pft  
           enddo
 
 
@@ -1411,6 +1411,7 @@ subroutine SeedGermination( litt, cold_stat, drought_stat )
            litt%seed_germ_in(pft) = 0.0_r8
        end if
 
+
     enddo
 
   end subroutine SeedGermination
@@ -1467,14 +1468,15 @@ subroutine recruitment( currentSite, currentPatch, bc_in )
     real(r8) :: mass_demand ! Total mass demanded by the plant to achieve the stoichiometric targets
                             ! of all the organs in the recruits. Used for both [kg per plant] and [kg per cohort]
     real(r8) :: stem_drop_fraction 
-                              
+                             
     !----------------------------------------------------------------------
 
     allocate(temp_cohort) ! create temporary cohort
     call zero_cohort(temp_cohort)
 
-    do ft = 1,numpft
 
+    do ft = 1,numpft
+      if(currentSite%use_this_pft(ft).eq.itrue)then
        temp_cohort%canopy_trim = 0.8_r8  !starting with the canopy not fully expanded 
        temp_cohort%pft         = ft
        temp_cohort%hite        = EDPftvarcon_inst%hgt_min(ft)
@@ -1711,6 +1713,7 @@ subroutine recruitment( currentSite, currentPatch, bc_in )
 
 
        endif
+      endif !use_this_pft
      enddo  !pft loop
 
      deallocate(temp_cohort) ! delete temporary cohort

main/EDInitMod.F90

@@ -35,10 +35,10 @@ module EDInitMod
   use EDTypesMod                , only : phen_dstat_moistoff
   use EDTypesMod                , only : phen_cstat_notcold
   use EDTypesMod                , only : phen_dstat_moiston
-  use EDTypesMod                , only : element_pos
   use FatesInterfaceTypesMod         , only : bc_in_type
   use FatesInterfaceTypesMod         , only : hlm_use_planthydro
   use FatesInterfaceTypesMod         , only : hlm_use_inventory_init
+  use FatesInterfaceTypesMod         , only : hlm_use_fixed_biogeog
   use FatesInterfaceTypesMod         , only : numpft
   use FatesInterfaceTypesMod         , only : nleafage
   use FatesInterfaceTypesMod         , only : nlevsclass
@@ -124,6 +124,9 @@ subroutine init_site_vars( site_in, bc_in )
     allocate(site_in%dz_soil(site_in%nlevsoil))
     allocate(site_in%z_soil(site_in%nlevsoil))
 
+    allocate(site_in%area_pft(1:numpft))
+    allocate(site_in%use_this_pft(1:numpft))
+
     do el=1,num_elements
         allocate(site_in%flux_diags(el)%leaf_litter_input(1:numpft))
         allocate(site_in%flux_diags(el)%root_litter_input(1:numpft))
@@ -136,7 +139,7 @@ subroutine init_site_vars( site_in, bc_in )
     site_in%zi_soil(:) = bc_in%zi_sisl(:)
     site_in%dz_soil(:) = bc_in%dz_sisl(:)
     site_in%z_soil(:)  = bc_in%z_sisl(:)
-
+    
 
     !
     end subroutine init_site_vars
@@ -218,10 +221,12 @@ subroutine zero_site( site_in )
     ! canopy spread
     site_in%spread = 0._r8
 
+    site_in%area_pft(:) = 0._r8
+    site_in%use_this_pft(:) = fates_unset_int
   end subroutine zero_site
 
   ! ============================================================================
-  subroutine set_site_properties( nsites, sites )
+  subroutine set_site_properties( nsites, sites,bc_in )
     !
     ! !DESCRIPTION:
     !
@@ -231,7 +236,7 @@ subroutine set_site_properties( nsites, sites )
 
     integer, intent(in)                        :: nsites
     type(ed_site_type) , intent(inout), target :: sites(nsites)
-
+    type(bc_in_type), intent(in)               :: bc_in(nsites)
     !
     ! !LOCAL VARIABLES:
     integer  :: s
@@ -244,6 +249,7 @@ subroutine set_site_properties( nsites, sites )
     integer  :: cleafoff   ! DOY for cold-decid leaf-off, initial guess
     integer  :: dleafoff   ! DOY for drought-decid leaf-off, initial guess
     integer  :: dleafon    ! DOY for drought-decid leaf-on, initial guess
+    integer  :: ft         ! PFT loop
     !----------------------------------------------------------------------
 
 
@@ -286,7 +292,24 @@ subroutine set_site_properties( nsites, sites )
           sites(s)%acc_NI     = acc_NI
           sites(s)%NF         = 0.0_r8         
           sites(s)%frac_burnt = 0.0_r8
-
+         
+         ! PLACEHOLDER FOR PFT AREA DATA MOVED ACROSS INTERFACE                                                                                   
+          if(hlm_use_fixed_biogeog.eq.itrue)then
+            do ft =  1,numpft
+              sites(s)%area_pft(ft) = bc_in(s)%pft_areafrac(ft)
+            end do
+          end if
+
+          do ft = 1,numpft
+           sites(s)%use_this_pft(ft) = itrue
+           if(hlm_use_fixed_biogeog.eq.itrue)then
+             if(sites(s)%area_pft(ft).gt.0.0_r8)then
+                sites(s)%use_this_pft(ft) = itrue
+             else
+                sites(s)%use_this_pft(ft) = ifalse
+             end if !area
+           end if !SBG
+          end do !ft
 
        end do
 
@@ -466,7 +489,7 @@ subroutine init_cohorts( site_in, patch_in, bc_in)
     patch_in%shortest => null()
 
     do pft =  1,numpft
-
+     if(site_in%use_this_pft(pft).eq.itrue)then
        if(EDPftvarcon_inst%initd(pft)>1.0E-7) then
 
        allocate(temp_cohort) ! temporary cohort
@@ -610,7 +633,7 @@ subroutine init_cohorts( site_in, patch_in, bc_in)
        deallocate(temp_cohort) ! get rid of temporary cohort
 
        endif
-
+     endif !use_this_pft
     enddo !numpft
 
     ! Zero the mass flux pools of the new cohorts

main/EDMainMod.F90

@@ -673,6 +673,7 @@ subroutine TotalBalanceCheck (currentSite, call_index )
           write(fates_log(),*) 'biomass', biomass_stock
           write(fates_log(),*) 'litter',litter_stock
           write(fates_log(),*) 'seeds',seed_stock
+          write(fates_log(),*) 'total stock', total_stock
           write(fates_log(),*) 'previous total',site_mass%old_stock  
           write(fates_log(),*) 'lat lon',currentSite%lat,currentSite%lon
 
@@ -690,6 +691,7 @@ subroutine TotalBalanceCheck (currentSite, call_index )
                 write(fates_log(),*) 'root litter (by layer): ',sum(litt%root_fines,dim=1)
                 write(fates_log(),*) 'dist mode: ',currentPatch%disturbance_mode
                 write(fates_log(),*) 'anthro_disturbance_label: ',currentPatch%anthro_disturbance_label
+                write(fates_log(),*) 'use_this_pft: ', currentSite%use_this_pft(:)
                 if(print_cohorts)then
                     write(fates_log(),*) '---- Biomass by cohort and organ -----'
                     currentCohort => currentPatch%tallest

main/EDTypesMod.F90

@@ -32,7 +32,7 @@ module EDTypesMod
                                                           ! are not the top canopy layer)
 
   integer, parameter, public :: nlevleaf = 30             ! number of leaf layers in canopy layer
-  integer, parameter, public :: maxpft = 15               ! maximum number of PFTs allowed
+  integer, parameter, public :: maxpft = 16               ! maximum number of PFTs allowed
                                                           ! the parameter file may determine that fewer
                                                           ! are used, but this helps allocate scratch
                                                           ! space and output arrays.
@@ -657,7 +657,11 @@ module EDTypesMod
      ! INDICES 
      real(r8) ::  lat                                          ! latitude:  degrees 
      real(r8) ::  lon                                          ! longitude: degrees 
-     
+
+     ! Fixed Biogeography mode inputs
+     real(r8), allocatable :: area_PFT(:)                      ! Area allocated to individual PFTs    
+     integer, allocatable  :: use_this_pft(:)                  ! Is area_PFT > 0 ? (1=yes, 0=no)
+ 
      ! Mass Balance (allocation for each element)
 
      type(site_massbal_type), pointer :: mass_balance(:)

main/FatesInterfaceMod.F90

@@ -76,6 +76,7 @@ module FatesInterfaceMod
    public :: allocate_bcin
    public :: allocate_bcout
    public :: zero_bcs
+   public :: set_bcs
 
 contains
 
@@ -354,6 +355,9 @@ subroutine allocate_bcin(bc_in, nlevsoil_in, nlevdecomp_in)
          allocate(bc_in%h2o_liq_sisl(nlevsoil_in)); bc_in%h2o_liq_sisl = nan
       end if
 
+         allocate(bc_in%pft_areafrac(maxpft))
+
+
       return
    end subroutine allocate_bcin
 
@@ -986,6 +990,8 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval)
          hlm_use_logging   = unset_int
          hlm_use_ed_st3    = unset_int
          hlm_use_ed_prescribed_phys = unset_int
+         hlm_use_fixed_biogeog = unset_int
+         !hlm_use_nocomp = unset_int    ! future reduced complexity mode
          hlm_use_inventory_init = unset_int
          hlm_inventory_ctrl_file = 'unset'
 
@@ -1186,6 +1192,21 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval)
             call endrun(msg=errMsg(sourcefile, __LINE__))
          end if
 
+        if(hlm_use_fixed_biogeog.eq.unset_int) then
+           if(fates_global_verbose()) then
+             write(fates_log(), *) 'switch for fixed biogeog unset: him_use_fixed_biogeog, exiting'
+           end if
+           call endrun(msg=errMsg(sourcefile, __LINE__))
+         end if
+
+        ! Future reduced complexity mode   
+        !if(hlm_use_nocomp.eq.unset_int) then
+        !      if(fates_global_verbose()) then
+        !     write(fates_log(), *) 'switch for no competition mode. '
+        !    end if
+        !   call endrun(msg=errMsg(sourcefile, __LINE__))
+        ! end if
+
          if(hlm_use_cohort_age_tracking .eq. unset_int) then
             if (fates_global_verbose()) then
                write(fates_log(), *) 'switch for cohort_age_tracking  unset: hlm_use_cohort_age_tracking, exiting'
@@ -1269,7 +1290,21 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval)
                if (fates_global_verbose()) then
                   write(fates_log(),*) 'Transfering hlm_use_spitfire= ',ival,' to FATES'
                end if
-               
+
+            case('use_fixed_biogeog')
+                hlm_use_fixed_biogeog = ival
+               if (fates_global_verbose()) then
+                   write(fates_log(),*) 'Transfering hlm_use_fixed_biogeog= ',ival,' to FATES'
+               end if
+            
+            ! Future reduced complexity mode   
+            !case('use_nocomp')
+            !    hlm_use_nocomp = ival
+            !   if (fates_global_verbose()) then
+            !       write(fates_log(),*) 'Transfering hlm_use_nocomp= ',ival,' to FATES'
+            !   end if
+
+
             case('use_planthydro')
                hlm_use_planthydro = ival
                if (fates_global_verbose()) then

main/FatesInterfaceTypesMod.F90

@@ -133,6 +133,9 @@ module FatesInterfaceTypesMod
                                                                     ! This need only be defined when
                                                                     ! hlm_use_inventory_init = 1
 
+  integer, public ::  hlm_use_fixed_biogeog                         !  Flag to use FATES fixed biogeography mode
+                                                                    !  1 = TRUE, 0 = FALSE 
+
    ! -------------------------------------------------------------------------------------
    ! Parameters that are dictated by FATES and known to be required knowledge
    !  needed by the HLMs
@@ -423,7 +426,10 @@ module FatesInterfaceTypesMod
       real(r8),allocatable :: hksat_sisl(:)        ! hydraulic conductivity at saturation (mm H2O /s)
       real(r8),allocatable :: h2o_liq_sisl(:)      ! Liquid water mass in each layer (kg/m2)
       real(r8) :: smpmin_si                        ! restriction for min of soil potential (mm)
-      
+   
+      ! Fixed biogeography mode 
+      real(r8), allocatable :: pft_areafrac(:)     ! Fractional area of the FATES column occupied by each PFT  
+    
    end type bc_in_type
 
 

main/FatesRestartInterfaceMod.F90

@@ -152,6 +152,8 @@ module FatesRestartInterfaceMod
   integer :: ir_seed_bank_sift
   integer :: ir_spread_si
   integer :: ir_recrate_sift
+  integer :: ir_use_this_pft_sift
+  integer :: ir_area_pft_sift
   integer :: ir_fmortrate_cano_siscpf
   integer :: ir_fmortrate_usto_siscpf
   integer :: ir_imortrate_siscpf
@@ -177,7 +179,6 @@ module FatesRestartInterfaceMod
   integer :: ir_errfates_mbal
   integer :: ir_prt_base     ! Base index for all PRT variables
 
-
   ! Hydraulic indices
   integer :: ir_hydro_th_ag_covec
   integer :: ir_hydro_th_troot
@@ -1005,6 +1006,17 @@ subroutine define_restart_vars(this, initialize_variables)
          units='indiv/ha/day', flushval = flushzero, &
          hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_recrate_sift)
 
+    call this%set_restart_var(vname='fates_use_this_pft', vtype=cohort_int, & !should this be cohort_int as above?
+         long_name='in fixed biogeog mode, is pft in gridcell?', &
+         units='0/1', flushval = flushzero, &
+         hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_use_this_pft_sift)
+
+    call this%set_restart_var(vname='fates_area_pft', vtype=cohort_r8, &
+         long_name='in fixed biogeog mode, what is pft area in gridcell?', &
+         units='0/1', flushval = flushzero, &
+         hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_area_pft_sift)
+
+
     call this%set_restart_var(vname='fates_fmortrate_canopy', vtype=cohort_r8, &
          long_name='fates diagnostics on fire mortality canopy', &
          units='indiv/ha/year', flushval = flushzero, &
@@ -1533,6 +1545,8 @@ subroutine set_restart_vectors(this,nc,nsites,sites)
            rio_watermem_siwm           => this%rvars(ir_watermem_siwm)%r81d, &
            rio_vegtempmem_sitm         => this%rvars(ir_vegtempmem_sitm)%r81d, &
            rio_recrate_sift            => this%rvars(ir_recrate_sift)%r81d, &
+           rio_use_this_pft_sift       => this%rvars(ir_use_this_pft_sift)%int1d, &
+           rio_area_pft_sift           => this%rvars(ir_area_pft_sift)%r81d, &
            rio_fmortrate_cano_siscpf   => this%rvars(ir_fmortrate_cano_siscpf)%r81d, &
            rio_fmortrate_usto_siscpf   => this%rvars(ir_fmortrate_usto_siscpf)%r81d, &
            rio_imortrate_siscpf        => this%rvars(ir_imortrate_siscpf)%r81d, &
@@ -1588,6 +1602,14 @@ subroutine set_restart_vectors(this,nc,nsites,sites)
              rio_recrate_sift(io_idx_co_1st+i_pft-1)   = sites(s)%recruitment_rate(i_pft)
           end do
 
+         do i_pft = 1,numpft
+             rio_use_this_pft_sift(io_idx_co_1st+i_pft-1)   = sites(s)%use_this_pft(i_pft)        
+         end do
+
+         do i_pft = 1,numpft
+              rio_area_pft_sift(io_idx_co_1st+i_pft-1)      = sites(s)%area_pft(i_pft)
+         end do
+
           do el = 1, num_elements
 
              io_idx_si_cwd = io_idx_co_1st
@@ -2270,6 +2292,8 @@ subroutine get_restart_vectors(this, nc, nsites, sites)
           rio_watermem_siwm           => this%rvars(ir_watermem_siwm)%r81d, &
           rio_vegtempmem_sitm         => this%rvars(ir_vegtempmem_sitm)%r81d, &
           rio_recrate_sift            => this%rvars(ir_recrate_sift)%r81d, &
+          rio_use_this_pft_sift       => this%rvars(ir_use_this_pft_sift)%int1d, &
+          rio_area_pft_sift           => this%rvars(ir_area_pft_sift)%r81d,&
           rio_fmortrate_cano_siscpf   => this%rvars(ir_fmortrate_cano_siscpf)%r81d, &
           rio_fmortrate_usto_siscpf   => this%rvars(ir_fmortrate_usto_siscpf)%r81d, &
           rio_imortrate_siscpf        => this%rvars(ir_imortrate_siscpf)%r81d, &
@@ -2313,7 +2337,12 @@ subroutine get_restart_vectors(this, nc, nsites, sites)
           do i_pft = 1,numpft 
              sites(s)%recruitment_rate(i_pft) = rio_recrate_sift(io_idx_co_1st+i_pft-1)
           enddo
-
+          
+         !variables for fixed biogeography mode. These are currently used in restart even when this is off. 
+          do i_pft = 1,numpft
+             sites(s)%use_this_pft(i_pft) = rio_use_this_pft_sift(io_idx_co_1st+i_pft-1)
+             sites(s)%area_pft(i_pft)     = rio_area_pft_sift(io_idx_co_1st+i_pft-1)
+          enddo
 
           ! Mass balance and diagnostics across elements at the site level
           do el = 1, num_elements