Merge pull request #705 from glemieux/frag-scalar-update

Update fragmentation_scaler calculation to use HLM moisture and temperature

biogeochem/EDPatchDynamicsMod.F90

@@ -1975,6 +1975,7 @@ subroutine create_patch(currentSite, new_patch, age, areap, label)
     allocate(new_patch%fabi(hlm_numSWb))
     allocate(new_patch%sabs_dir(hlm_numSWb))
     allocate(new_patch%sabs_dif(hlm_numSWb))
+    allocate(new_patch%fragmentation_scaler(currentSite%nlevsoil))
 
 
     ! Litter
@@ -2794,6 +2795,7 @@ subroutine dealloc_patch(cpatch)
        deallocate(cpatch%fabi)
        deallocate(cpatch%sabs_dir)
        deallocate(cpatch%sabs_dif)
+       deallocate(cpatch%fragmentation_scaler)
 
     end if
 

biogeochem/EDPhysiologyMod.F90

@@ -2248,43 +2248,56 @@ subroutine fragmentation_scaler( currentPatch, bc_in)
     !
     ! !LOCAL VARIABLES:
     logical  :: use_century_tfunc = .false.
+    logical  :: use_hlm_soil_scalar = .false. ! Use hlm input decomp fraction scalars
     integer  :: j
-    integer  :: ifp                   ! Index of a FATES Patch "ifp"
-    real(r8) :: t_scalar
-    real(r8) :: w_scalar
-    real(r8) :: catanf                ! hyperbolic temperature function from CENTURY
-    real(r8) :: catanf_30             ! hyperbolic temperature function from CENTURY
-    real(r8) :: t1                    ! temperature argument
+    integer  :: ifp                          ! Index of a FATES Patch "ifp"
+    real(r8) :: t_scalar                     ! temperature scalar
+    real(r8) :: w_scalar                     ! moisture scalar
+    real(r8) :: catanf                       ! hyperbolic temperature function from CENTURY
+    real(r8) :: catanf_30                    ! hyperbolic temperature function from CENTURY
+    real(r8) :: t1                           ! temperature argument
     !----------------------------------------------------------------------
 
     catanf(t1) = 11.75_r8 +(29.7_r8 / pi) * atan( pi * 0.031_r8  * ( t1 - 15.4_r8 ))
     catanf_30 = catanf(30._r8)
 
     ifp = currentPatch%patchno 
 
-    if ( .not. use_century_tfunc ) then
-    !calculate rate constant scalar for soil temperature,assuming that the base rate constants 
-    !are assigned for non-moisture limiting conditions at 25C.
-      if (bc_in%t_veg24_pa(ifp)  >=  tfrz) then
-        t_scalar = q10_mr**((bc_in%t_veg24_pa(ifp)-(tfrz+25._r8))/10._r8)
-                 !  Q10**((t_soisno(c,j)-(tfrz+25._r8))/10._r8)
-      else
-        t_scalar = (q10_mr**(-25._r8/10._r8))*(q10_froz**((bc_in%t_veg24_pa(ifp)-tfrz)/10._r8))
-                  !Q10**(-25._r8/10._r8))*(froz_q10**((t_soisno(c,j)-tfrz)/10._r8)
-      endif
+    ! Use the hlm temp and moisture decomp fractions by default
+    if ( use_hlm_soil_scalar ) then
+
+      ! Calculate the fragmentation_scaler
+      currentPatch%fragmentation_scaler =  min(1.0_r8,max(0.0_r8,bc_in%t_scalar_sisl * bc_in%w_scalar_sisl))
+
     else
-      ! original century uses an arctangent function to calculate the 
-      ! temperature dependence of decomposition      
-      t_scalar = max(catanf(bc_in%t_veg24_pa(ifp)-tfrz)/catanf_30,0.01_r8)
-    endif    
-
-    !Moisture Limitations   
-    !BTRAN APPROACH - is quite simple, but max's out decomp at all unstressed 
-    !soil moisture values, which is not realistic.  
-    !litter decomp is proportional to water limitation on average... 
-    w_scalar = sum(currentPatch%btran_ft(1:numpft))/real(numpft,r8)
+
+      if ( .not. use_century_tfunc ) then
+      !calculate rate constant scalar for soil temperature,assuming that the base rate constants 
+      !are assigned for non-moisture limiting conditions at 25C.
+         if (bc_in%t_veg24_pa(ifp)  >=  tfrz) then
+         t_scalar = q10_mr**((bc_in%t_veg24_pa(ifp)-(tfrz+25._r8))/10._r8)
+                  !  Q10**((t_soisno(c,j)-(tfrz+25._r8))/10._r8)
+         else
+         t_scalar = (q10_mr**(-25._r8/10._r8))*(q10_froz**((bc_in%t_veg24_pa(ifp)-tfrz)/10._r8))
+                     !Q10**(-25._r8/10._r8))*(froz_q10**((t_soisno(c,j)-tfrz)/10._r8)
+         endif
+      else
+         ! original century uses an arctangent function to calculate the 
+         ! temperature dependence of decomposition      
+         t_scalar = max(catanf(bc_in%t_veg24_pa(ifp)-tfrz)/catanf_30,0.01_r8)
+      endif    
+
+      !Moisture Limitations   
+      !BTRAN APPROACH - is quite simple, but max's out decomp at all unstressed 
+      !soil moisture values, which is not realistic.  
+      !litter decomp is proportional to water limitation on average... 
+      w_scalar = sum(currentPatch%btran_ft(1:numpft))/real(numpft,r8)
 
-    currentPatch%fragmentation_scaler =  min(1.0_r8,max(0.0_r8,t_scalar * w_scalar))
+    ! Calculate the fragmentation_scaler
+      currentPatch%fragmentation_scaler(:) =  min(1.0_r8,max(0.0_r8,t_scalar * w_scalar))
+
+   endif
+
 
   end subroutine fragmentation_scaler
 
@@ -2302,8 +2315,7 @@ subroutine CWDOut( litt, fragmentation_scaler, nlev_eff_decomp )
     !
     ! !ARGUMENTS    
     type(litter_type),intent(inout),target     :: litt
-
-    real(r8),intent(in)                        :: fragmentation_scaler
+    real(r8),intent(in)                        :: fragmentation_scaler(:)
 
     ! This is not necessarily every soil layer, this is the number
     ! of effective layers that are active and can be sent
@@ -2312,21 +2324,26 @@ subroutine CWDOut( litt, fragmentation_scaler, nlev_eff_decomp )
 
     !
     ! !LOCAL VARIABLES:
-    integer :: c
-    integer :: ilyr
-    integer :: dcmpy
-    integer :: numlevsoil
+    integer :: c                       ! Fuel size class index
+    integer :: ilyr                    ! Soil layer index
+    integer :: dcmpy                   ! Decomposibility pool indexer
+    integer :: soil_layer_index = 1    ! Soil layer index associated with above ground litter
     !----------------------------------------------------------------------
 
+
+    ! Above ground litters are associated with the top soil layer temperature and
+    ! moisture scalars and fragmentation scalar associated with specified index value
+    ! is used for ag_cwd_frag and root_fines_frag calculations.
+
     do c = 1,ncwd  
 
        litt%ag_cwd_frag(c)   = litt%ag_cwd(c) * SF_val_max_decomp(c) * &
-             years_per_day * fragmentation_scaler
+             years_per_day * fragmentation_scaler(soil_layer_index)
 
        do ilyr = 1,nlev_eff_decomp
 
            litt%bg_cwd_frag(c,ilyr) = litt%bg_cwd(c,ilyr) * SF_val_max_decomp(c) * &
-                years_per_day * fragmentation_scaler
+                years_per_day * fragmentation_scaler(ilyr)
 
        enddo
     end do
@@ -2339,11 +2356,11 @@ subroutine CWDOut( litt, fragmentation_scaler, nlev_eff_decomp )
     do dcmpy = 1,ndcmpy
 
        litt%leaf_fines_frag(dcmpy) = litt%leaf_fines(dcmpy) * &
-             years_per_day * SF_val_max_decomp(dl_sf) * fragmentation_scaler
+             years_per_day * SF_val_max_decomp(dl_sf) * fragmentation_scaler(soil_layer_index)
 
        do ilyr = 1,nlev_eff_decomp
            litt%root_fines_frag(dcmpy,ilyr) = litt%root_fines(dcmpy,ilyr) * &
-                 years_per_day *  SF_val_max_decomp(dl_sf) * fragmentation_scaler
+                 years_per_day *  SF_val_max_decomp(dl_sf) * fragmentation_scaler(ilyr)
        end do
     enddo
 

main/EDTypesMod.F90

@@ -526,7 +526,7 @@ module EDTypesMod
 
      type(litter_type), pointer :: litter(:)  ! Litter (leaf,fnrt,CWD and seeds) for different elements
 
-     real(r8) :: fragmentation_scaler          ! Scale rate of litter fragmentation. 0 to 1.
+     real(r8),allocatable :: fragmentation_scaler(:)            ! Scale rate of litter fragmentation based on soil layer. 0 to 1.
 
      real(r8) ::  repro(maxpft)                                 ! allocation to reproduction per PFT : KgC/m2
 

main/FatesInterfaceMod.F90

@@ -220,6 +220,8 @@ subroutine zero_bcs(fates,s)
     fates%bc_in(s)%tot_litc            = 0.0_r8
     fates%bc_in(s)%snow_depth_si       = 0.0_r8
     fates%bc_in(s)%frac_sno_eff_si     = 0.0_r8
+    fates%bc_in(s)%w_scalar_sisl(:)    = 0.0_r8
+    fates%bc_in(s)%t_scalar_sisl(:)    = 0.0_r8
 
     if(do_fates_salinity)then
        fates%bc_in(s)%salinity_sl(:)   = 0.0_r8
@@ -397,6 +399,8 @@ subroutine allocate_bcin(bc_in, nlevsoil_in, nlevdecomp_in, num_lu_harvest_cats)
       allocate(bc_in%z_sisl(nlevsoil_in))
       allocate(bc_in%decomp_id(nlevsoil_in))
       allocate(bc_in%dz_decomp_sisl(nlevdecomp_in))
+      allocate(bc_in%w_scalar_sisl(nlevsoil_in))
+      allocate(bc_in%t_scalar_sisl(nlevsoil_in))
 
       ! Lightning (or successful ignitions) and population density
       allocate(bc_in%lightning24(maxPatchesPerSite))

main/FatesInterfaceTypesMod.F90

@@ -328,6 +328,10 @@ module FatesInterfaceTypesMod
                                                  ! soil layer maps to. This will either
                                                  ! be equivalent (ie integer ascending)
                                                  ! Or, all will be 1.
+
+      ! Decomposition fractions
+      real(r8),allocatable :: w_scalar_sisl(:)   ! fraction by which decomposition is limited by moisture availability
+      real(r8),allocatable :: t_scalar_sisl(:)   ! fraction by which decomposition is limited by temperature
 
 
       ! Vegetation Dynamics