Stem area index AND canopy snow radiation fixes.

biogeochem/EDPatchDynamicsMod.F90

@@ -430,7 +430,7 @@ subroutine disturbance_rates( site_in, bc_in)
     enddo !patch loop 
 
   end subroutine disturbance_rates
-  
+
     ! ============================================================================
 
   subroutine spawn_patches( currentSite, bc_in)
@@ -451,7 +451,7 @@ subroutine spawn_patches( currentSite, bc_in)
     ! !USES:
 
     use EDParamsMod         , only : ED_val_understorey_death, logging_coll_under_frac
-    use EDCohortDynamicsMod , only : zero_cohort, copy_cohort, terminate_cohorts 
+    use EDCohortDynamicsMod , only : zero_cohort, copy_cohort, terminate_cohorts
     use FatesConstantsMod   , only : rsnbl_math_prec
 
     !
@@ -2201,6 +2201,7 @@ subroutine zero_patch(cp_p)
     currentPatch%c_stomata                  = 0.0_r8 ! This is calculated immediately before use
     currentPatch%c_lblayer                  = 0.0_r8
     currentPatch%fragmentation_scaler(:)    = 0.0_r8
+    currentPatch%radiation_error            = 0.0_r8
 
     ! diagnostic radiation profiles
     currentPatch%nrmlzd_parprof_pft_dir_z(:,:,:,:) = 0._r8
@@ -2210,6 +2211,7 @@ subroutine zero_patch(cp_p)
 
     currentPatch%solar_zenith_flag          = .false.
     currentPatch%solar_zenith_angle         = nan
+    currentPatch%fcansno                    = nan
 
     currentPatch%gnd_alb_dir(:)             = nan
     currentPatch%gnd_alb_dif(:)             = nan
@@ -2397,7 +2399,7 @@ subroutine fuse_patches( csite, bc_in )
                          !-------------------------------------------------------------------------!
 
                          if(fuse_flag  ==  1)then
-
+                            
                             !-----------------------!
                             ! fuse the two patches  !
                             !-----------------------!
@@ -2552,7 +2554,8 @@ subroutine fuse_2_patches(csite, dp, rp)
     rp%zstar                = (dp%zstar*dp%area + rp%zstar*rp%area) * inv_sum_area
     rp%c_stomata            = (dp%c_stomata*dp%area + rp%c_stomata*rp%area) * inv_sum_area
     rp%c_lblayer            = (dp%c_lblayer*dp%area + rp%c_lblayer*rp%area) * inv_sum_area
-
+    rp%radiation_error      = (dp%radiation_error*dp%area + rp%radiation_error*rp%area) * inv_sum_area
+
     rp%area = rp%area + dp%area !THIS MUST COME AT THE END!
 
     !insert donor cohorts into recipient patch
@@ -2674,9 +2677,9 @@ subroutine terminate_patches(currentSite)
 
     currentPatch => currentSite%youngest_patch
     do while(associated(currentPatch)) 
-
+       
        if(currentPatch%area <= min_patch_area)then
-
+          
 
           ! Even if the patch area is small, avoid fusing it into its neighbor
           ! if it is the youngest of all patches. We do this in attempts to maintain

main/EDTypesMod.F90

@@ -430,7 +430,7 @@ module EDTypesMod
      real(r8) ::  elai_profile(nclmax,maxpft,nlevleaf)          ! exposed leaf area in each canopy layer, pft, and leaf layer
      real(r8) ::  tsai_profile(nclmax,maxpft,nlevleaf)          ! total   stem area in each canopy layer, pft, and leaf layer
      real(r8) ::  esai_profile(nclmax,maxpft,nlevleaf)          ! exposed stem area in each canopy layer, pft, and leaf layer
-
+     real(r8) ::  radiation_error                               ! radiation error (w/m2)
      real(r8) ::  layer_height_profile(nclmax,maxpft,nlevleaf)
      real(r8) ::  canopy_area_profile(nclmax,maxpft,nlevleaf)   ! fraction of crown area per canopy area in each layer
                                                                 ! they will sum to 1.0 in the fully closed canopy layers
@@ -444,6 +444,7 @@ module EDTypesMod
      integer  ::  ncan(nclmax,maxpft)                           ! number of total   leaf layers for each canopy layer and pft
 
      !RADIATION FLUXES      
+     real(r8) :: fcansno                                        ! Fraction of canopy covered in snow
 
      logical  ::  solar_zenith_flag                             ! integer flag specifying daylight (based on zenith angle)
      real(r8) ::  solar_zenith_angle                            ! solar zenith angle (radians)
@@ -661,7 +662,7 @@ module EDTypesMod
      procedure :: ZeroMassBalFlux
 
   end type site_massbal_type
-
+  
 
   !************************************
   !** Site type structure           **
@@ -685,7 +686,7 @@ module EDTypesMod
      ! Global index of this site in the history output file
      integer :: h_gid
 
-
+     
      ! INDICES 
      real(r8) ::  lat                                          ! latitude:  degrees 
      real(r8) ::  lon                                          ! longitude: degrees 

main/FatesHistoryInterfaceMod.F90

@@ -304,6 +304,7 @@ module FatesHistoryInterfaceMod
 
   integer :: ih_c_stomata_si
   integer :: ih_c_lblayer_si
+  integer :: ih_rad_error_si
 
   integer :: ih_fire_c_to_atm_si
 
@@ -3431,6 +3432,7 @@ subroutine update_history_hifrq(this,nc,nsites,sites,bc_in,dt_tstep)
                hio_growth_resp_si => this%hvars(ih_growth_resp_si)%r81d, &
                hio_c_stomata_si   => this%hvars(ih_c_stomata_si)%r81d, &
                hio_c_lblayer_si   => this%hvars(ih_c_lblayer_si)%r81d, &
+               hio_rad_error_si   => this%hvars(ih_rad_error_si)%r81d, &
                hio_nep_si         => this%hvars(ih_nep_si)%r81d, &
                hio_hr_si          => this%hvars(ih_hr_si)%r81d, &
                hio_ar_si_scpf     => this%hvars(ih_ar_si_scpf)%r82d, &
@@ -3537,6 +3539,9 @@ subroutine update_history_hifrq(this,nc,nsites,sites,bc_in,dt_tstep)
             hio_c_lblayer_si(io_si) = hio_c_lblayer_si(io_si) + &
                  cpatch%c_lblayer * cpatch%total_canopy_area
 
+          hio_rad_error_si(io_si) = hio_rad_error_si(io_si) + &
+                 cpatch%radiation_error * cpatch%area * AREA_INV
+
             ccohort => cpatch%shortest
             do while(associated(ccohort))
 
@@ -4910,6 +4915,14 @@ subroutine define_history_vars(this, initialize_variables)
          avgflag='A', vtype=site_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=2,   &
          ivar=ivar, initialize=initialize_variables, index = ih_c_lblayer_si )
 
+! radiation error 
+
+    call this%set_history_var(vname='RAD_ERROR', units='W m-2 ',                   &
+         long='radiation error in FATES RTM', use_default='active',                   &
+         avgflag='A', vtype=site_r8, hlms='CLM:ALM', flushval=0.0_r8, upfreq=2,   &
+         ivar=ivar, initialize=initialize_variables, index = ih_rad_error_si )
+
+
 
     ! Ecosystem Carbon Fluxes (updated rapidly, upfreq=2)
 

main/FatesInterfaceMod.F90

@@ -250,6 +250,7 @@ subroutine zero_bcs(fates,s)
     fates%bc_in(s)%h2o_liqvol_sl(:)    = 0.0_r8
     fates%bc_in(s)%filter_vegzen_pa(:) = .false.
     fates%bc_in(s)%coszen_pa(:)        = 0.0_r8
+    fates%bc_in(s)%fcansno_pa(:)       = 0.0_r8
     fates%bc_in(s)%albgr_dir_rb(:)     = 0.0_r8
     fates%bc_in(s)%albgr_dif_rb(:)     = 0.0_r8
     fates%bc_in(s)%max_rooting_depth_index_col = 0
@@ -337,7 +338,7 @@ subroutine zero_bcs(fates,s)
     fates%bc_out(s)%z0m_pa(:)    = 0.0_r8
     fates%bc_out(s)%dleaf_pa(:)   = 0.0_r8
     fates%bc_out(s)%nocomp_pft_label_pa(:) = 0
-
+    
     fates%bc_out(s)%canopy_fraction_pa(:) = 0.0_r8
     fates%bc_out(s)%frac_veg_nosno_alb_pa(:) = 0.0_r8
 
@@ -346,7 +347,7 @@ subroutine zero_bcs(fates,s)
        fates%bc_out(s)%qflx_ro_sisl(:)        = 0.0_r8
     end if
     fates%bc_out(s)%plant_stored_h2o_si = 0.0_r8
-    
+
     return
   end subroutine zero_bcs
 
@@ -488,6 +489,7 @@ subroutine allocate_bcin(bc_in, nlevsoil_in, nlevdecomp_in, num_lu_harvest_cats)
       ! Canopy Radiation
       allocate(bc_in%filter_vegzen_pa(maxPatchesPerSite))
       allocate(bc_in%coszen_pa(maxPatchesPerSite))
+      allocate(bc_in%fcansno_pa(maxPatchesPerSite))
       allocate(bc_in%albgr_dir_rb(hlm_numSWb))
       allocate(bc_in%albgr_dif_rb(hlm_numSWb))
 
@@ -606,7 +608,7 @@ subroutine allocate_bcout(bc_out, nlevsoil_in, nlevdecomp_in)
          bc_out%rootfr_pa(0,1:nlevsoil_in)=1._r8/real(nlevsoil_in,r8)
       end if
 
-
+         
       ! Fates -> BGC fragmentation mass fluxes
       select case(hlm_parteh_mode) 
       case(prt_carbon_allom_hyp)
@@ -650,7 +652,7 @@ subroutine allocate_bcout(bc_out, nlevsoil_in, nlevdecomp_in)
 
       allocate(bc_out%canopy_fraction_pa(maxPatchesPerSite))
       allocate(bc_out%frac_veg_nosno_alb_pa(maxPatchesPerSite))
-     
+
       allocate(bc_out%nocomp_pft_label_pa(maxPatchesPerSite))
 
       ! Plant-Hydro BC's
@@ -1521,9 +1523,9 @@ subroutine set_fates_ctrlparms(tag,ival,rval,cval)
                call endrun(msg=errMsg(sourcefile, __LINE__))
             end if
          end if
-
 
-         if(hlm_use_fixed_biogeog.eq.unset_int) then
+
+        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
@@ -1824,7 +1826,7 @@ subroutine FatesReportParameters(masterproc)
       call PRTDerivedParams()              ! Update PARTEH derived constants
       call PRTCheckParams(masterproc)      ! Check PARTEH parameters
       call SpitFireCheckParams(masterproc)
-
+      
 
 
       return

main/FatesInterfaceTypesMod.F90

@@ -404,7 +404,7 @@ module FatesInterfaceTypesMod
       ! 2 = patch is currently marked for photosynthesis
       ! 3 = patch has been called for photosynthesis at least once
       integer, allocatable  :: filter_photo_pa(:)
-     
+
       ! atmospheric pressure (Pa)
       real(r8)              :: forc_pbot             
 
@@ -447,7 +447,10 @@ module FatesInterfaceTypesMod
       !           I am leaving it at this scale for simplicity.  Patches should
       !           have no spacially variable information
       real(r8), allocatable :: coszen_pa(:)
-
+
+      ! fraction of canopy that is covered in snow
+      real(r8), allocatable :: fcansno_pa(:)
+
       ! Abledo of the ground for direct radiation, by site broadband (0-1)
       real(r8), allocatable :: albgr_dir_rb(:)
 
@@ -681,7 +684,7 @@ module FatesInterfaceTypesMod
 
      integer, allocatable :: nocomp_pft_label_pa(:) ! in nocomp and SP mode, each patch has a PFT identity. 
 
-       ! FATES Hydraulics
+      ! FATES Hydraulics
 
 
 
@@ -735,7 +738,7 @@ module FatesInterfaceTypesMod
                                                 ! increasing, or all 1s)
 
    end type bc_pconst_type
-   
+
 
 
  contains

main/FatesRestartInterfaceMod.F90

@@ -138,6 +138,7 @@ module FatesRestartInterfaceMod
   integer :: ir_lmort_infra_co
 
   ! Radiation
+  integer :: ir_fcansno_pa
   integer :: ir_solar_zenith_flag_pa
   integer :: ir_solar_zenith_angle_pa
   integer :: ir_gnd_alb_dif_pasb
@@ -654,6 +655,10 @@ subroutine define_restart_vars(this, initialize_variables)
          long_name='the number of cohorts per patch', units='unitless', flushval = flushinvalid, &
          hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_ncohort_pa )
 
+    call this%set_restart_var(vname='fates_fcansno_pa', vtype=cohort_r8, &
+         long_name='Fraction of canopy covered in snow', units='unitless', flushval = flushinvalid, &
+         hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_fcansno_pa )
+
     call this%set_restart_var(vname='fates_solar_zenith_flag_pa', vtype=cohort_int, &
          long_name='switch specifying if zenith is positive', units='unitless', flushval = flushinvalid, &
          hlms='CLM:ALM', initialize=initialize_variables, ivar=ivar, index = ir_solar_zenith_flag_pa )
@@ -1651,6 +1656,7 @@ subroutine set_restart_vectors(this,nc,nsites,sites)
            rio_snow_depth_si           => this%rvars(ir_snow_depth_si)%r81d, &
            rio_trunk_product_si        => this%rvars(ir_trunk_product_si)%r81d, &
            rio_ncohort_pa              => this%rvars(ir_ncohort_pa)%int1d, &
+           rio_fcansno_pa              => this%rvars(ir_fcansno_pa)%r81d, &
            rio_solar_zenith_flag_pa    => this%rvars(ir_solar_zenith_flag_pa)%int1d, &
            rio_solar_zenith_angle_pa   => this%rvars(ir_solar_zenith_angle_pa)%r81d, &
            rio_canopy_layer_co         => this%rvars(ir_canopy_layer_co)%int1d, &
@@ -1971,6 +1977,8 @@ subroutine set_restart_vectors(this,nc,nsites,sites)
              ! set cohorts per patch for IO
              rio_ncohort_pa( io_idx_co_1st )   = cohortsperpatch
 
+             rio_fcansno_pa( io_idx_co_1st )   = cpatch%fcansno
+
              ! Set zenith angle info
              if ( cpatch%solar_zenith_flag ) then
                 rio_solar_zenith_flag_pa(io_idx_co_1st)     = itrue
@@ -2466,6 +2474,7 @@ subroutine get_restart_vectors(this, nc, nsites, sites)
           rio_snow_depth_si           => this%rvars(ir_snow_depth_si)%r81d, &
           rio_trunk_product_si        => this%rvars(ir_trunk_product_si)%r81d, &
           rio_ncohort_pa              => this%rvars(ir_ncohort_pa)%int1d, &
+          rio_fcansno_pa              => this%rvars(ir_fcansno_pa)%r81d, &
           rio_solar_zenith_flag_pa    => this%rvars(ir_solar_zenith_flag_pa)%int1d, &
           rio_solar_zenith_angle_pa   => this%rvars(ir_solar_zenith_angle_pa)%r81d, &
           rio_canopy_layer_co         => this%rvars(ir_canopy_layer_co)%int1d, &
@@ -2770,6 +2779,7 @@ subroutine get_restart_vectors(this, nc, nsites, sites)
              cpatch%nocomp_pft_label               = rio_nocomp_pft_label_pa(io_idx_co_1st)
              cpatch%area               = rio_area_pa(io_idx_co_1st)
              cpatch%age_class          = get_age_class_index(cpatch%age)
+             cpatch%fcansno            = rio_fcansno_pa(io_idx_co_1st)
 
              ! Set zenith angle info
              cpatch%solar_zenith_flag  = ( rio_solar_zenith_flag_pa(io_idx_co_1st) .eq. itrue )