Parameterizations of fine-scale topographic effects on surface radiation balance in ELM

cime_config/tests.py

@@ -101,7 +101,8 @@
             "ERS.ELM_USRDAT.I1850ELM.elm-usrdat",
             "ERS.r05_r05.IELM.elm-lnd_rof_2way",
             "ERS.r05_r05.IELM.elm-V2_ELM_MOSART_features",
-            "ERS.ELM_USRDAT.IELM.elm-surface_water_dynamics"
+            "ERS.ELM_USRDAT.IELM.elm-surface_water_dynamics",
+            "ERS.ELM_USRDAT.IELM.elm-finetop_rad"
             )
         },
 

components/elm/bld/namelist_files/namelist_defaults.xml

@@ -661,6 +661,9 @@ this mask will have smb calculated over the entire global land surface
 <!--ftopdat>same as fsurfdat </ftopdat-->
 <use_top_solar_rad>.false.</use_top_solar_rad>
 
+<!-- Use the radiation parameterization for fine-scale topographic effects -->
+<use_finetop_rad>.false.</use_finetop_rad>
+
 <!-- SNICAR (SNow, ICe, and Aerosol Radiative model) datasets -->
 <!--  ***********  Resolution independent:   ***********  -->
 <fsnowoptics >lnd/clm2/snicardata/snicar_optics_5bnd_mam_c160322.nc</fsnowoptics>

components/elm/bld/namelist_files/namelist_definition.xml

@@ -760,6 +760,13 @@ If TRUE, TOP solar radiation parameterization will be activated, which considers
 </entry>
 <!-- End TOP solar radiation parameterization                   -->
 
+<!-- fineTOP radiation parameterization -->
+<entry id="use_finetop_rad" type="logical" category="elm_physics"
+       group="elm_inparm" valid_values="" value=".false.">
+If TRUE, fineTOP radiation parameterization will be activated, which considers the effects of fine(grid)-scale topography.
+</entry>
+<!-- End fineTOP radiation parameterization -->
+
 <!-- ========================================================================================  -->
 <!-- Former CPP tokens -->
 <!-- ========================================================================================  -->

components/elm/cime_config/testdefs/testmods_dirs/elm/finetop_rad/shell_commands

@@ -0,0 +1,10 @@
+./xmlchange DATM_MODE=CLMMOSARTTEST
+
+./xmlchange LND_DOMAIN_FILE=domain.lnd.1km_SierraNevada.c240628.nc
+./xmlchange ATM_DOMAIN_FILE=domain.lnd.1km_SierraNevada.c240628.nc
+./xmlchange LND_DOMAIN_PATH="\$DIN_LOC_ROOT/share/domains/domain.clm"
+./xmlchange ATM_DOMAIN_PATH="\$DIN_LOC_ROOT/share/domains/domain.clm"
+
+./xmlchange -file env_run.xml -id DATM_CLMNCEP_YR_END -val 2000
+./xmlchange -file env_run.xml -id DATM_CLMNCEP_YR_START -val 2000
+./xmlchange -file env_run.xml -id DATM_CLMNCEP_YR_ALIGN -val 1

components/elm/cime_config/testdefs/testmods_dirs/elm/finetop_rad/user_nl_elm

@@ -0,0 +1,3 @@
+fsurdat = '$DIN_LOC_ROOT/lnd/clm2/surfdata_map/surfdata_1km_SierraNevada_simyr2010_c240628.nc'
+use_finetop_rad = .true.
+use_top_solar_rad  = .false.

components/elm/cime_config/testdefs/testmods_dirs/elm/snowveg_arctic/shell_commands

@@ -7,4 +7,4 @@
 ./xmlchange DATM_CLMNCEP_YR_END=2000
 ./xmlchange ELM_USRDAT_NAME=42_FLUXNETSITES
 ./xmlchange NTASKS=1
-./xmlchange NTHRDS=1
+./xmlchange NTHRDS=1

components/elm/src/biogeophys/CanopyFluxesMod.F90

@@ -14,7 +14,7 @@ module CanopyFluxesMod
   use shr_log_mod           , only : errMsg => shr_log_errMsg
   use abortutils            , only : endrun
   use elm_varctl            , only : iulog, use_cn, use_lch4, use_c13, use_c14, use_fates
-  use elm_varctl            , only : use_hydrstress
+  use elm_varctl            , only : use_hydrstress, use_finetop_rad
   use elm_varpar            , only : nlevgrnd, nlevsno
   use elm_varcon            , only : namep
   use elm_varcon            , only : mm_epsilon
@@ -100,6 +100,7 @@ subroutine CanopyFluxes(bounds,  num_nolakeurbanp, filter_nolakeurbanp, &
     use elm_varcon         , only : isecspday, degpsec
     use pftvarcon          , only : irrigated
     use elm_varcon         , only : c14ratio
+    use shr_const_mod      , only : SHR_CONST_PI
 
     !NEW
     use elm_varsur         , only : firrig
@@ -314,7 +315,7 @@ subroutine CanopyFluxes(bounds,  num_nolakeurbanp, filter_nolakeurbanp, &
     real(r8) :: tau_diff(bounds%begp:bounds%endp) ! Difference from previous iteration tau
     real(r8) :: prev_tau(bounds%begp:bounds%endp) ! Previous iteration tau
     real(r8) :: prev_tau_diff(bounds%begp:bounds%endp) ! Previous difference in iteration tau
-
+    real(r8) :: slope_rad, deg2rad
     character(len=64) :: event !! timing event
 
     ! Indices for raw and rah
@@ -329,6 +330,7 @@ subroutine CanopyFluxes(bounds,  num_nolakeurbanp, filter_nolakeurbanp, &
          snl                  => col_pp%snl                                   , & ! Input:  [integer  (:)   ]  number of snow layers
          dayl                 => grc_pp%dayl                                  , & ! Input:  [real(r8) (:)   ]  daylength (s)
          max_dayl             => grc_pp%max_dayl                              , & ! Input:  [real(r8) (:)   ]  maximum daylength for this grid cell (s)
+         slope_deg            => grc_pp%slope_deg                             , &
 
          forc_lwrad           => top_af%lwrad                              , & ! Input:  [real(r8) (:)   ]  downward infrared (longwave) radiation (W/m**2)                       
          forc_q               => top_as%qbot                               , & ! Input:  [real(r8) (:)   ]  atmospheric specific humidity (kg/kg)                                 
@@ -516,7 +518,7 @@ subroutine CanopyFluxes(bounds,  num_nolakeurbanp, filter_nolakeurbanp, &
       end if
 #endif
 
-
+      deg2rad = SHR_CONST_PI/180._r8
       ! Initialize
       do f = 1, fn
          p = filterp(f)
@@ -701,6 +703,11 @@ subroutine CanopyFluxes(bounds,  num_nolakeurbanp, filter_nolakeurbanp, &
          bir(p) = - (2._r8-emv(p)*(1._r8-emg(c))) * emv(p) * sb
          cir(p) =   emv(p)*emg(c)*sb
 
+         if (use_finetop_rad) then
+            slope_rad = slope_deg(g) * deg2rad
+            bir(p) = bir(p) / cos(slope_rad)
+            cir(p) = cir(p) / cos(slope_rad)
+         endif
          ! Saturated vapor pressure, specific humidity, and their derivatives
          ! at the leaf surface
 
@@ -1269,15 +1276,25 @@ subroutine CanopyFluxes(bounds,  num_nolakeurbanp, filter_nolakeurbanp, &
 
          ! Downward longwave radiation below the canopy
 
-         dlrad(p) = (1._r8-emv(p))*emg(c)*forc_lwrad(t) + &
-              emv(p)*emg(c)*sb*tlbef(p)**3*(tlbef(p) + 4._r8*dt_veg(p))
-
+         if (use_finetop_rad) then
+            slope_rad = slope_deg(g) * deg2rad
+            dlrad(p) = (1._r8-emv(p))*emg(c)*forc_lwrad(t) + &
+                  emv(p)*emg(c)*sb*tlbef(p)**3*(tlbef(p) + 4._r8*dt_veg(p))/cos(slope_rad)
+         else
+            dlrad(p) = (1._r8-emv(p))*emg(c)*forc_lwrad(t) + &
+                  emv(p)*emg(c)*sb*tlbef(p)**3*(tlbef(p) + 4._r8*dt_veg(p))
+         endif
          ! Upward longwave radiation above the canopy
-
-         ulrad(p) = ((1._r8-emg(c))*(1._r8-emv(p))*(1._r8-emv(p))*forc_lwrad(t) &
-              + emv(p)*(1._r8+(1._r8-emg(c))*(1._r8-emv(p)))*sb*tlbef(p)**3*(tlbef(p) + &
-              4._r8*dt_veg(p)) + emg(c)*(1._r8-emv(p))*sb*lw_grnd)
-
+         if (use_finetop_rad) then
+            slope_rad = slope_deg(g) * deg2rad
+            ulrad(p) = ((1._r8-emg(c))*(1._r8-emv(p))*(1._r8-emv(p))*forc_lwrad(t) &
+                + emv(p)*(1._r8+(1._r8-emg(c))*(1._r8-emv(p)))*sb*tlbef(p)**3*(tlbef(p) + &
+                4._r8*dt_veg(p))/cos(slope_rad) + emg(c)*(1._r8-emv(p))*sb*lw_grnd/cos(slope_rad))
+         else
+            ulrad(p) = ((1._r8-emg(c))*(1._r8-emv(p))*(1._r8-emv(p))*forc_lwrad(t) &
+                + emv(p)*(1._r8+(1._r8-emg(c))*(1._r8-emv(p)))*sb*tlbef(p)**3*(tlbef(p) + &
+                4._r8*dt_veg(p)) + emg(c)*(1._r8-emv(p))*sb*lw_grnd)
+         endif
          ! Derivative of soil energy flux with respect to soil temperature
 
          cgrnds(p) = cgrnds(p) + cpair*forc_rho(t)*wtg(p)*wtal(p)

components/elm/src/biogeophys/CanopyHydrologyMod.F90

@@ -29,7 +29,7 @@ module CanopyHydrologyMod
   use elm_varcon        , only : snw_rds_min
   use pftvarcon         , only : irrigated
   use GridcellType      , only : grc_pp
-  use timeinfoMod, only : dtime_mod
+  use timeinfoMod       , only : dtime_mod
   !
   ! !PUBLIC TYPES:
   implicit none
@@ -713,7 +713,7 @@ subroutine CanopyHydrology(bounds, &
        call FracH2oSfc(bounds, num_nolakec, filter_nolakec, &
              col_wf%qflx_h2osfc2topsoi, dtime)
 
-     end associate 
+     end associate
 
    end subroutine CanopyHydrology
 

components/elm/src/biogeophys/CanopyTemperatureMod.F90

@@ -13,7 +13,7 @@ module CanopyTemperatureMod
   use shr_kind_mod         , only : r8 => shr_kind_r8
   use shr_const_mod        , only : SHR_CONST_PI
   use decompMod            , only : bounds_type
-  use elm_varctl           , only : iulog, use_fates
+  use elm_varctl           , only : iulog, use_fates, use_finetop_rad
   use PhotosynthesisMod    , only : Photosynthesis, PhotosynthesisTotal, Fractionation
   use elm_instMod          , only : alm_fates
   use SurfaceResistanceMod , only : calc_soilevap_stress
@@ -29,6 +29,7 @@ module CanopyTemperatureMod
   use ColumnDataType       , only : col_es, col_ef, col_ws
   use VegetationType       , only : veg_pp
   use VegetationDataType   , only : veg_es, veg_ef, veg_wf
+  use GridcellType         , only : grc_pp
   !
   ! !PUBLIC TYPES:
   implicit none
@@ -74,6 +75,7 @@ subroutine CanopyTemperature(bounds, &
     use column_varcon      , only : icol_road_imperv, icol_road_perv
     use landunit_varcon    , only : istice, istice_mec, istwet, istsoil, istdlak, istcrop, istdlak
     use elm_varpar         , only : nlevgrnd, nlevurb, nlevsno, nlevsoi
+    use shr_const_mod   , only : SHR_CONST_PI
     !
     ! !ARGUMENTS:
     type(bounds_type)      , intent(in)    :: bounds
@@ -109,6 +111,7 @@ subroutine CanopyTemperature(bounds, &
     real(r8) :: vol_ice      ! partial volume of ice lens in layer
     real(r8) :: vol_liq      ! partial volume of liquid water in layer
     real(r8) :: fh2o_eff(bounds%begc:bounds%endc) ! effective surface water fraction (i.e. seen by atm)
+    real(r8) :: slope_rad, deg2rad
     !------------------------------------------------------------------------------
 
     associate(                                                          &
@@ -200,6 +203,7 @@ subroutine CanopyTemperature(bounds, &
          tssbef           =>    col_es%t_ssbef                          & ! Output: [real(r8) (:,:) ] soil/snow temperature before update (K)
          )
 
+      deg2rad = SHR_CONST_PI/180._r8
       do j = -nlevsno+1, nlevgrnd
          do fc = 1,num_nolakec
             c = filter_nolakec(fc)
@@ -409,6 +413,7 @@ subroutine CanopyTemperature(bounds, &
 
       do fp = 1,num_nolakep
          p = filter_nolakep(fp)
+         g = veg_pp%gridcell(p)
 
          ! Initial set (needed for history tape fields)
 
@@ -439,7 +444,12 @@ subroutine CanopyTemperature(bounds, &
          ! Vegetation Emissivity
 
          avmuir = 1._r8
-         emv(p) = 1._r8-exp(-(elai(p)+esai(p))/avmuir)
+         if (use_finetop_rad .and. (.not. lun_pp%urbpoi(l))) then
+            slope_rad = grc_pp%slope_deg(g) * deg2rad
+            emv(p) = 1._r8-exp(-(elai(p)+esai(p))*cos(slope_rad)/avmuir)
+         else
+            emv(p) = 1._r8-exp(-(elai(p)+esai(p))/avmuir)
+         endif
 
          z0mv(p)   = z0m(p)
          z0hv(p)   = z0mv(p)

components/elm/src/biogeophys/SnowSnicarMod.F90

@@ -1812,7 +1812,7 @@ subroutine SNICAR_AD_RT (flg_snw_ice, bounds, num_nourbanc, filter_nourbanc,  &
      use elm_varpar       , only : nlevsno, numrad
      use elm_time_manager , only : get_nstep
      use shr_const_mod    , only : SHR_CONST_PI
-     use elm_varctl       , only : snow_shape, snicar_atm_type, use_dust_snow_internal_mixing
+     use elm_varctl       , only : snow_shape, snicar_atm_type, use_dust_snow_internal_mixing, use_finetop_rad
      !
      ! !ARGUMENTS:
      integer           , intent(in)  :: flg_snw_ice                                        ! flag: =1 when called from CLM, =2 when called from CSIM
@@ -1868,6 +1868,7 @@ subroutine SNICAR_AD_RT (flg_snw_ice, bounds, num_nourbanc, filter_nourbanc,  &
      !integer :: trip                               ! flag: =1 to redo RT calculation if result is unrealistic
      !integer :: flg_dover                          ! defines conditions for RT redo (explained below)
 
+     real(r8):: slope_rad, deg2rad
      real(r8):: albedo                             ! temporary snow albedo [frc]
      real(r8):: flx_sum                            ! temporary summation variable for NIR weighting
      real(r8):: albout_lcl(numrad_snw)             ! snow albedo by band [frc]
@@ -2140,6 +2141,7 @@ subroutine SNICAR_AD_RT (flg_snw_ice, bounds, num_nourbanc, filter_nourbanc,  &
        ! Define constants
        pi = SHR_CONST_PI
        nint_snw_rds_min = nint(snw_rds_min)
+       deg2rad = SHR_CONST_PI/180._r8
 
        ! always use Delta approximation for snow
        DELTA = 1
@@ -2202,6 +2204,7 @@ subroutine SNICAR_AD_RT (flg_snw_ice, bounds, num_nourbanc, filter_nourbanc,  &
       ! (when called from CSIM, there is only one column)
        do fc = 1,num_nourbanc
           c_idx = filter_nourbanc(fc)
+          g_idx = col_pp%gridcell(c_idx)
 
           ! Zero absorbed radiative fluxes:
           do i=-nlevsno+1,1,1
@@ -2246,7 +2249,6 @@ subroutine SNICAR_AD_RT (flg_snw_ice, bounds, num_nourbanc, filter_nourbanc,  &
 
                 ! for debugging only
                 l_idx     = col_pp%landunit(c_idx)
-                g_idx     = col_pp%gridcell(c_idx)
                 sfctype   = lun_pp%itype(l_idx)
                 lat_coord = grc_pp%latdeg(g_idx)
                 lon_coord = grc_pp%londeg(g_idx)
@@ -2266,6 +2268,20 @@ subroutine SNICAR_AD_RT (flg_snw_ice, bounds, num_nourbanc, filter_nourbanc,  &
                 lon_coord         = 0
              endif ! end if flg_snw_ice == 1
 
+             if (use_finetop_rad) then
+                slope_rad = grc_pp%slope_deg(g_idx) * deg2rad
+
+                if ((flg_snw_ice == 1) .and. (snl(c_idx) > -1)) then
+                   h2osno_liq_lcl(0) = h2osno_liq_lcl(0) * cos(slope_rad)
+                   h2osno_ice_lcl(0) = h2osno_ice_lcl(0) * cos(slope_rad)
+                else
+                   h2osno_liq_lcl(:) = h2osno_liq_lcl(:) * cos(slope_rad)
+                   h2osno_ice_lcl(:) = h2osno_ice_lcl(:) * cos(slope_rad)
+                endif
+
+                h2osno_lcl = h2osno_lcl * cos(slope_rad)
+             endif
+
 #ifdef MODAL_AER
            !mgf++
            !

components/elm/src/biogeophys/SoilFluxesMod.F90

@@ -8,7 +8,7 @@ module SoilFluxesMod
   use shr_log_mod	, only : errMsg => shr_log_errMsg
   use decompMod		, only : bounds_type
   use abortutils	, only : endrun
-  use elm_varctl	, only : iulog, use_firn_percolation_and_compaction
+  use elm_varctl	, only : iulog, use_firn_percolation_and_compaction, use_finetop_rad
   use perfMod_GPU
   use elm_varpar	, only : nlevsno, nlevgrnd, nlevurb, max_patch_per_col
   use atm2lndType	, only : atm2lnd_type
@@ -21,6 +21,7 @@ module SoilFluxesMod
   use ColumnDataType    , only : col_es, col_ef, col_ws
   use VegetationType    , only : veg_pp
   use VegetationDataType, only : veg_ef, veg_wf
+  use GridcellType      , only : grc_pp
 
   use timeinfoMod , only : dtime_mod
   !
@@ -49,6 +50,7 @@ subroutine SoilFluxes (bounds, num_urbanl, filter_urbanl, &
     use landunit_varcon  , only : istsoil, istcrop
     use column_varcon    , only : icol_roof, icol_sunwall, icol_shadewall, icol_road_perv
     use subgridAveMod    , only : p2c
+    use shr_const_mod      , only : SHR_CONST_PI
     !
     ! !ARGUMENTS:
     type(bounds_type)      , intent(in)    :: bounds
@@ -63,7 +65,7 @@ subroutine SoilFluxes (bounds, num_urbanl, filter_urbanl, &
     type(canopystate_type) , intent(in)    :: canopystate_vars
     type(energyflux_type)  , intent(inout) :: energyflux_vars
     real(r8) :: dtime                                              ! land model time step (sec)
-
+    real(r8) :: slope_rad, deg2rad
 
     !
     ! !LOCAL VARIABLES:
@@ -159,10 +161,13 @@ subroutine SoilFluxes (bounds, num_urbanl, filter_urbanl, &
          eflx_lh_vege            => veg_ef%eflx_lh_vege      , & ! Output: [real(r8) (:)   ]  veg evaporation heat flux (W/m**2) [+ to atm]
          eflx_lh_vegt            => veg_ef%eflx_lh_vegt      , & ! Output: [real(r8) (:)   ]  veg transpiration heat flux (W/m**2) [+ to atm]
          eflx_lh_grnd            => veg_ef%eflx_lh_grnd      , & ! Output: [real(r8) (:)   ]  ground evaporation heat flux (W/m**2) [+ to atm]
-         errsoi_col              => col_ef%errsoi              , & ! Output: [real(r8) (:)   ]  column-level soil/lake energy conservation error (W/m**2)
-         errsoi_patch            => veg_ef%errsoi              & ! Output: [real(r8) (:)   ]  pft-level soil/lake energy conservation error (W/m**2)
+         errsoi_col              => col_ef%errsoi            , & ! Output: [real(r8) (:)   ]  column-level soil/lake energy conservation error (W/m**2)
+         errsoi_patch            => veg_ef%errsoi            , & ! Output: [real(r8) (:)   ]  pft-level soil/lake energy conservation error (W/m**2)
+         slope_deg               => grc_pp%slope_deg           &
          )
 
+      deg2rad = SHR_CONST_PI/180._r8
+
          dtime = dtime_mod
       event = 'bgp2_loop_1'
       call t_start_lnd(event)
@@ -282,11 +287,19 @@ subroutine SoilFluxes (bounds, num_urbanl, filter_urbanl, &
             lw_grnd=(frac_sno_eff(c)*tssbef(c,col_pp%snl(c)+1)**4 &
                  +(1._r8-frac_sno_eff(c)-frac_h2osfc(c))*tssbef(c,1)**4 &
                  +frac_h2osfc(c)*t_h2osfc_bef(c)**4)
-
-            eflx_soil_grnd(p) = ((1._r8- frac_sno_eff(c))*sabg_soil(p) + frac_sno_eff(c)*sabg_snow(p)) + dlrad(p) &
-                 + (1-frac_veg_nosno(p))*emg(c)*forc_lwrad(t) &
-                 - emg(c)*sb*lw_grnd - emg(c)*sb*t_grnd0(c)**3*(4._r8*tinc(c)) &
-                 - (eflx_sh_grnd(p)+qflx_evap_soi(p)*htvp(c))
+
+            if (use_finetop_rad) then
+               slope_rad = slope_deg(g) * deg2rad
+               eflx_soil_grnd(p) = ((1._r8- frac_sno_eff(c))*sabg_soil(p) + frac_sno_eff(c)*sabg_snow(p)) + dlrad(p) &
+                    + (1-frac_veg_nosno(p))*emg(c)*forc_lwrad(t) &
+                    - emg(c)*sb*lw_grnd/cos(slope_rad)- emg(c)*sb*t_grnd0(c)**3*(4._r8*tinc(c))/cos(slope_rad) &
+                    - (eflx_sh_grnd(p)+qflx_evap_soi(p)*htvp(c))
+            else
+               eflx_soil_grnd(p) = ((1._r8- frac_sno_eff(c))*sabg_soil(p) + frac_sno_eff(c)*sabg_snow(p)) + dlrad(p) &
+                    + (1-frac_veg_nosno(p))*emg(c)*forc_lwrad(t) &
+                    - emg(c)*sb*lw_grnd - emg(c)*sb*t_grnd0(c)**3*(4._r8*tinc(c)) &
+                    - (eflx_sh_grnd(p)+qflx_evap_soi(p)*htvp(c))
+            endif
 
             if (lun_pp%itype(l) == istsoil .or. lun_pp%itype(l) == istcrop) then
                eflx_soil_grnd_r(p) = eflx_soil_grnd(p)
@@ -301,9 +314,9 @@ subroutine SoilFluxes (bounds, num_urbanl, filter_urbanl, &
             ! Include transpiration term because needed for pervious road
             ! and wasteheat and traffic flux
             eflx_soil_grnd(p) = sabg(p) + dlrad(p) &
-                 - eflx_lwrad_net(p) - eflx_lwrad_del(p) &
-                 - (eflx_sh_grnd(p) + qflx_evap_soi(p)*htvp(c) + qflx_tran_veg(p)*hvap) &
-                 + eflx_wasteheat_patch(p) + eflx_heat_from_ac_patch(p) + eflx_traffic_patch(p)
+                - eflx_lwrad_net(p) - eflx_lwrad_del(p) &
+                - (eflx_sh_grnd(p) + qflx_evap_soi(p)*htvp(c) + qflx_tran_veg(p)*hvap) &
+                + eflx_wasteheat_patch(p) + eflx_heat_from_ac_patch(p) + eflx_traffic_patch(p)
             eflx_soil_grnd_u(p) = eflx_soil_grnd(p)
          end if
 
@@ -426,10 +439,18 @@ subroutine SoilFluxes (bounds, num_urbanl, filter_urbanl, &
                  +(1._r8-frac_sno_eff(c)-frac_h2osfc(c))*tssbef(c,1)**4 &
                  +frac_h2osfc(c)*t_h2osfc_bef(c)**4)
 
-            eflx_lwrad_out(p) = ulrad(p) &
-                 + (1-frac_veg_nosno(p))*(1.-emg(c))*forc_lwrad(t) &
-                 + (1-frac_veg_nosno(p))*emg(c)*sb*lw_grnd &
-                 + 4._r8*emg(c)*sb*t_grnd0(c)**3*tinc(c)
+            if (use_finetop_rad) then
+               slope_rad = slope_deg(g) * deg2rad
+               eflx_lwrad_out(p) = ulrad(p) &
+                    + (1-frac_veg_nosno(p))*(1.-emg(c))*forc_lwrad(t) &
+                    + (1-frac_veg_nosno(p))*emg(c)*sb*lw_grnd / cos(slope_rad) &
+                    + 4._r8*emg(c)*sb*t_grnd0(c)**3*tinc(c) / cos(slope_rad)
+            else
+               eflx_lwrad_out(p) = ulrad(p) &
+                    + (1-frac_veg_nosno(p))*(1.-emg(c))*forc_lwrad(t) &
+                    + (1-frac_veg_nosno(p))*emg(c)*sb*lw_grnd &
+                    + 4._r8*emg(c)*sb*t_grnd0(c)**3*tinc(c)
+            endif
 
             eflx_lwrad_net(p) = eflx_lwrad_out(p) - forc_lwrad(t)
             if (lun_pp%itype(l) == istsoil .or. lun_pp%itype(l) == istcrop) then