Merge pull request #6 from rgknox/xuchongang/fates_hydro_improve_read…

…y_to_merge

merge resolution (trivial)

biogeochem/EDPatchDynamicsMod.F90

@@ -1145,14 +1145,20 @@ subroutine mortality_litter_fluxes(currentSite, cp_target, new_patch_target, pat
              !not right to recalcualte dmort here.
              canopy_dead = currentCohort%n * min(1.0_r8,currentCohort%dmort * hlm_freq_day * fates_mortality_disturbance_fraction)
 
-
-
              canopy_mortality_woody_litter(p)= canopy_mortality_woody_litter(p) + &
                   canopy_dead*(struct_c + sapw_c)
              canopy_mortality_leaf_litter(p) = canopy_mortality_leaf_litter(p) + &
                   canopy_dead*leaf_c
+
+             ! Some plants upon death will transfer storage carbon to seed production
+             ! Storage carbon that is not transferred to seeds goes to root litter flux
+
              canopy_mortality_root_litter(p) = canopy_mortality_root_litter(p) + &
-                  canopy_dead*(fnrt_c + store_c)
+                  canopy_dead*(fnrt_c + store_c*(1.0_r8-EDPftvarcon_inst%allom_frbstor_repro(p)) )
+
+             currentSite%seed_bank(p) = currentSite%seed_bank(p) + &
+                   canopy_dead * store_c * EDPftvarcon_inst%allom_frbstor_repro(p)/AREA
+
 
              if( hlm_use_planthydro == itrue ) then
                 call AccumulateMortalityWaterStorage(currentSite,currentCohort, canopy_dead)

biogeochem/EDPhysiologyMod.F90

@@ -180,7 +180,7 @@ subroutine non_canopy_derivs( currentSite, currentPatch, bc_in )
     ! update fragmenting pool fluxes
     call cwd_input( currentSite, currentPatch)
     call cwd_out( currentSite, currentPatch, bc_in)
-
+  
     do p = 1,numpft
        currentSite%dseed_dt(p) = currentSite%dseed_dt(p) + &
             (currentPatch%seeds_in(p) - currentPatch%seed_decay(p) - &
@@ -800,6 +800,7 @@ subroutine seeds_in( currentSite, cp_pnt )
     type(ed_cohort_type), pointer :: currentCohort
     integer :: p
     logical :: pft_present(maxpft)
+    real(r8) :: store_c_to_repro   ! carbon sent from storage to reproduction upon death [kg/plant]
     real(r8) :: npfts_present
     !----------------------------------------------------------------------
 
@@ -830,10 +831,18 @@ subroutine seeds_in( currentSite, cp_pnt )
        currentPatch => cp_pnt
        currentCohort => currentPatch%tallest
        do while (associated(currentCohort))
+
+          ! a certain fraction of bstore goes to clonal reproduction when plants die
+          store_c_to_repro = currentCohort%prt%GetState(store_organ,all_carbon_elements) * &
+                EDPftvarcon_inst%allom_frbstor_repro(currentCohort%pft)
+
           do p = 1, numpft
              if (pft_present(p)) then
-                currentPatch%seeds_in(p) = currentPatch%seeds_in(p) +  currentCohort%seed_prod * currentCohort%n / &
-                     (currentPatch%area * npfts_present)
+
+                  currentPatch%seeds_in(p) = currentPatch%seeds_in(p) + &
+                        (currentCohort%seed_prod * currentCohort%n - &
+                        currentCohort%dndt*store_c_to_repro) &
+                        /(currentPatch%area * npfts_present)		  
              endif
           end do
           currentCohort => currentCohort%shorter
@@ -844,8 +853,15 @@ subroutine seeds_in( currentSite, cp_pnt )
     currentCohort => currentPatch%tallest
     do while (associated(currentCohort))
        p = currentCohort%pft
+
+       ! a certain fraction of bstore goes to clonal reproduction when plants die
+       store_c_to_repro = currentCohort%prt%GetState(store_organ,all_carbon_elements) * &
+             EDPftvarcon_inst%allom_frbstor_repro(p)
+
        currentPatch%seeds_in(p) = currentPatch%seeds_in(p) +  &
-             currentCohort%seed_prod * currentCohort%n/currentPatch%area
+           (currentCohort%seed_prod * currentCohort%n - &
+	   currentCohort%dndt*store_c_to_repro)/currentPatch%area
+
        currentCohort => currentCohort%shorter
     enddo !cohort loop
 
@@ -921,13 +937,18 @@ subroutine seed_germination( currentSite, currentPatch )
 
     do p = 1,numpft
        currentPatch%seed_germination(p) =  min(currentSite%seed_bank(p) * &
-             EDPftvarcon_inst%germination_timescale(p),max_germination)
+             EDPftvarcon_inst%germination_timescale(p),max_germination)     
+       !set the germination only under the growing season...c.xu
+       if (EDPftvarcon_inst%season_decid(p) == 1.and.currentSite%status == 1)then 
+             currentPatch%seed_germination(p) = 0.0_r8
+       endif
+       if (EDPftvarcon_inst%stress_decid(p) == 1.and.currentSite%dstatus == 1)then
+             currentPatch%seed_germination(p) = 0.0_r8
+       endif
     enddo
 
   end subroutine seed_germination
 
-  ! ============================================================================
-
   subroutine recruitment( currentSite, currentPatch, bc_in )
     !
     ! !DESCRIPTION:
@@ -1182,9 +1203,8 @@ subroutine CWD_Input( currentSite, currentPatch)
                 ! the litter flux has already been counted since it captured
                 ! the losses of live trees and those flagged for death
 
-
           currentPatch%root_litter_in(pft) = currentPatch%root_litter_in(pft) + &
-               (fnrt_c + store_c ) * dead_n
+               (fnrt_c + store_c*(1._r8-EDPftvarcon_inst%allom_frbstor_repro(pft)) ) * dead_n
 
           ! Update diagnostics that track resource management
           currentSite%resources_management%delta_litter_stock  = &

biogeochem/FatesAllometryMod.F90

@@ -2184,7 +2184,7 @@ subroutine StructureResetOfDH( bdead, ipft, canopy_trim, d, h )
      ! of the diameter increment
      counter = 0
      step_frac = step_frac0
-     do while( (bdead-bt_dead) > calloc_abs_error )
+     do while( (bdead-bt_dead) > calloc_abs_error .and. dbt_dead_dd>0.0_r8)
 
         ! vulnerable to div0
         dd    = step_frac*(bdead-bt_dead)/dbt_dead_dd

biogeophys/EDBtranMod.F90

@@ -25,6 +25,7 @@ module EDBtranMod
 
    public :: btran_ed
    public :: get_active_suction_layers
+   public :: check_layer_water
 
 contains 
 

biogeophys/FatesBstressMod.F90

@@ -0,0 +1,99 @@
+module FatesBstressMod
+
+   !-------------------------------------------------------------------------------------
+   ! Description: calculate the stress impact on transpiration from salinity and sulphide in soils
+   !              note that water stress is calculated in EDBtranMod or HYDRO
+   ! ------------------------------------------------------------------------------------
+   !
+   use FatesConstantsMod , only : tfrz => t_water_freeze_k_1atm 
+   use FatesConstantsMod , only : itrue,ifalse
+   use EDTypesMod        , only : ed_site_type,       &
+                                  ed_patch_type,      &
+                                  ed_cohort_type,     &
+                                  maxpft
+   use shr_kind_mod      , only : r8 => shr_kind_r8
+   use FatesInterfaceMod , only : bc_in_type, &
+                                  bc_out_type, &
+                                  numpft
+   use FatesInterfaceMod , only : hlm_use_planthydro
+   use FatesGlobals      , only : fates_log
+   use EDBtranMod        , only : check_layer_water
+
+   implicit none
+   private
+
+   public :: btran_sal_stress_fates
+
+contains
+ ! =====================================================================================
+
+  subroutine btran_sal_stress_fates( nsites, sites, bc_in)
+
+
+      ! ---------------------------------------------------------------------------------
+      ! Calculate the transpiration wetness function (BTRAN) and the root uptake
+      ! distribution (ROOTR).
+      ! Boundary conditions in: bc_in(s)%salinity_sl(j)        salinity concontration[ppm]
+      ! Boundary conditions in: bc_in(s)%sulphide_sl(j)        sulphide concontration[ppm]
+      ! Output cpatch%bstress_sal_ft(ft)
+      ! Output cpatch%bstress_sul_ft(ft)
+      ! ---------------------------------------------------------------------------------
+
+      ! Arguments
+
+      integer,intent(in)                      :: nsites
+      type(ed_site_type),intent(inout),target :: sites(nsites)
+      type(bc_in_type),intent(in)             :: bc_in(nsites)
+
+      !
+      ! !LOCAL VARIABLES:
+      type(ed_patch_type),pointer             :: cpatch ! Current Patch Pointer
+      type(ed_cohort_type),pointer            :: ccohort ! Current cohort pointer
+      integer  :: s                 ! site
+      integer  :: j                 ! soil layer
+      integer  :: ft                ! plant functional type index
+      real(r8) :: salinity_node     ! salinity in the soil water [ppt]
+      real(r8) :: rresis            ! salinity limitation to transpiration independent
+      !------------------------------------------------------------------------------
+
+        do s = 1,nsites
+
+           cpatch => sites(s)%oldest_patch
+           do while (associated(cpatch))                 
+
+              ! THIS SHOULD REALLY BE A COHORT LOOP ONCE WE HAVE rootfr_ft FOR COHORTS (RGK)
+
+              do ft = 1,numpft
+                 cpatch%bstress_sal_ft(ft) = 0.0_r8
+                 do j = 1,bc_in(s)%nlevsoil
+
+                    ! Calculations are only relevant where liquid water exists
+                    ! see clm_fates%wrap_btran for calculation with CLM/ELM
+
+                    if ( check_layer_water(bc_in(s)%h2o_liqvol_sl(j),bc_in(s)%tempk_sl(j)) )  then
+
+                       salinity_node =  bc_in(s)%salinity_sl(j)
+
+                       rresis  = min( 1.244_r8/(1+exp((0.186_r8-salinity_node)/(-0.132_r8))), 1._r8)
+
+                       cpatch%bstress_sal_ft(ft) = cpatch%bstress_sal_ft(ft)+ &
+		                                     cpatch%rootfr_ft(ft,j)*rresis
+
+                    end if
+
+                 end do !j
+
+              end do !PFT              
+
+              cpatch => cpatch%younger
+
+           end do
+
+        end do
+
+
+    end subroutine btran_sal_stress_fates
+
+  ! ====================================================================================
+
+end module FatesBstressMod

biogeophys/FatesPlantRespPhotosynthMod.F90

@@ -20,27 +20,27 @@ module FATESPlantRespPhotosynthMod
 
    ! !USES:
 
-   use FatesGlobals, only      : endrun => fates_endrun
-   use FatesGlobals, only      : fates_log
+   use FatesGlobals,      only : endrun => fates_endrun
+   use FatesGlobals,      only : fates_log
    use FatesConstantsMod, only : r8 => fates_r8
    use FatesConstantsMod, only : itrue
    use FatesInterfaceMod, only : hlm_use_planthydro
    use FatesInterfaceMod, only : hlm_parteh_mode
    use FatesInterfaceMod, only : numpft
-   use EDTypesMod, only        : maxpft
-   use EDTypesMod, only        : nlevleaf
-   use EDTypesMod, only        : nclmax
-
-   use PRTGenericMod,          only : prt_carbon_allom_hyp
-   use PRTGenericMod,          only : prt_cnp_flex_allom_hyp 
-   use PRTGenericMod,          only : all_carbon_elements
-   use PRTGenericMod,          only : nitrogen_element
-   use PRTGenericMod,          only : leaf_organ
-   use PRTGenericMod,          only : fnrt_organ
-   use PRTGenericMod,          only : sapw_organ
-   use PRTGenericMod,          only : store_organ
-   use PRTGenericMod,          only : repro_organ
-   use PRTGenericMod,          only : struct_organ
+   use EDTypesMod,        only : maxpft
+   use EDTypesMod,        only : nlevleaf
+   use EDTypesMod,        only : nclmax
+   use EDTypesMod,        only : do_fates_salinity 
+   use PRTGenericMod,     only : prt_carbon_allom_hyp
+   use PRTGenericMod,     only : prt_cnp_flex_allom_hyp 
+   use PRTGenericMod,     only : all_carbon_elements
+   use PRTGenericMod,     only : nitrogen_element
+   use PRTGenericMod,     only : leaf_organ
+   use PRTGenericMod,     only : fnrt_organ
+   use PRTGenericMod,     only : sapw_organ
+   use PRTGenericMod,     only : store_organ
+   use PRTGenericMod,     only : repro_organ
+   use PRTGenericMod,     only : struct_organ
 
    ! CIME Globals
    use shr_log_mod , only      : errMsg => shr_log_errMsg
@@ -423,7 +423,11 @@ subroutine FatesPlantRespPhotosynthDrive (nsites, sites,bc_in,bc_out,dtime)
 
 
                               end if
-
+
+                              if(do_fates_salinity)then
+                                btran_eff = btran_eff*currentPatch%bstress_sal_ft(ft)
+                              endif 
+
 
                               ! Scale for leaf nitrogen profile
                               nscaler = exp(-kn(ft) * cumulative_lai)

main/EDParamsMod.F90

@@ -43,6 +43,7 @@ module EDParamsMod
    real(r8),protected :: ED_val_patch_fusion_tol
    real(r8),protected :: ED_val_canopy_closure_thresh ! site-level canopy closure point where trees take on forest (narrow) versus savannah (wide) crown allometry
 
+
    ! two special parameters whose size is defined in the parameter file
    real(r8),protected,allocatable :: ED_val_history_sizeclass_bin_edges(:)
    real(r8),protected,allocatable :: ED_val_history_ageclass_bin_edges(:)
@@ -87,7 +88,11 @@ module EDParamsMod
    real(r8),protected :: hydr_psicap        !  sapwood water potential at which capillary reserves exhausted (MPa)
    character(len=param_string_length),parameter :: hydr_name_psicap = "fates_hydr_psicap"
 
-
+   !Soil BGC parameters, mostly used for testing FATES when not coupled to the dynamics bgc hlm
+   ! ----------------------------------------------------------------------------------------------
+   real(r8),protected :: bgc_soil_salinity ! site-level soil salinity for FATES when not coupled to dynamic soil BGC of salinity
+   character(len=param_string_length),parameter :: bgc_name_soil_salinity= "fates_soil_salinity"      
+
    ! Logging Control Parameters (ONLY RELEVANT WHEN USE_FATES_LOGGING = TRUE)
    ! ----------------------------------------------------------------------------------------------
 
@@ -155,6 +160,8 @@ subroutine FatesParamsInit()
     hydr_kmax_rsurf                       = nan
     hydr_psi0                             = nan
     hydr_psicap                           = nan
+
+    bgc_soil_salinity                     = nan
 
     logging_dbhmin                        = nan
     logging_collateral_frac               = nan
@@ -261,6 +268,9 @@ subroutine FatesRegisterParams(fates_params)
     call fates_params%RegisterParameter(name=hydr_name_psicap, dimension_shape=dimension_shape_1d, &
          dimension_names=dim_names)
 
+    call fates_params%RegisterParameter(name=bgc_name_soil_salinity, dimension_shape=dimension_shape_1d, &
+         dimension_names=dim_names) 
+
     call fates_params%RegisterParameter(name=logging_name_dbhmin, dimension_shape=dimension_shape_1d, &
          dimension_names=dim_names)
 
@@ -378,6 +388,9 @@ subroutine FatesReceiveParams(fates_params)
 
     call fates_params%RetreiveParameter(name=hydr_name_psicap, &
           data=hydr_psicap)
+
+    call fates_params%RetreiveParameter(name=bgc_name_soil_salinity, &
+          data=bgc_soil_salinity)	  
 
     call fates_params%RetreiveParameter(name=logging_name_dbhmin, &
           data=logging_dbhmin)
@@ -450,6 +463,7 @@ subroutine FatesReportParams(is_master)
 	write(fates_log(),fmt0) 'hydr_kmax_rsurf = ',hydr_kmax_rsurf  
         write(fates_log(),fmt0) 'hydr_psi0 = ',hydr_psi0
         write(fates_log(),fmt0) 'hydr_psicap = ',hydr_psicap
+        write(fates_log(),fmt0) 'bgc_soil_salinity = ', bgc_soil_salinity
         write(fates_log(),fmt0) 'logging_dbhmin = ',logging_dbhmin
         write(fates_log(),fmt0) 'logging_collateral_frac = ',logging_collateral_frac
         write(fates_log(),fmt0) 'logging_coll_under_frac = ',logging_coll_under_frac