add finetop

debug & define cosinc_col

add model outputs; & modify surface area correction for longwave radiation & debug

correct surface area for thermal radiation for soilT module

update snowsnicar_AD model

change COSINC output name to avoid duplicate

fix a small bug in variable name

fix small bugs

debug

debug

fix bugs

set SZA<85 limit

fix small bugs

fix one bug and add initialization for few variables

fix small bugs

delete log output and add control to adjust albedo

delete log output

delete write output

add more outputs like MODIS daytime and nighttime ST

debug

fix a small bug

add MODIS FSNO output

fix a bug

clean up the code in biogeophys

clean up the code in main

clean up the code

fix urban bug

add code to calculate apparent surface albedo for EAM

fix a small bug

debug

debug & modify code to not change for urban landunits

debug

fix typo

debug

consider the top effects on upward radaiton to sky

debug

Conflicts:
	components/elm/src/biogeophys/SoilFluxesMod.F90
	components/elm/src/biogeophys/SurfaceAlbedoMod.F90
	components/elm/src/main/atm2lndType.F90
	components/elm/src/main/elm_initializeMod.F90

fix a LAI bug

delete ;

debug & add test

add restart variables

fix a bug

debug

debug restart issues

remove log output

Author: Dalei Hao

cime_config/tests.py

@@ -100,7 +100,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
@@ -2266,6 +2269,13 @@ 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
+                h2osno_liq_lcl = h2osno_liq_lcl * cos(slope_rad)
+                h2osno_ice_lcl = h2osno_ice_lcl * cos(slope_rad)
+                h2osno_lcl = h2osno_lcl * cos(slope_rad)
+             endif
+
 #ifdef MODAL_AER
            !mgf++
            !