diff --git a/cime_config/tests.py b/cime_config/tests.py
index e06aa029a6aa..5152404b68db 100644
--- a/cime_config/tests.py
+++ b/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"
)
},
diff --git a/components/elm/bld/namelist_files/namelist_defaults.xml b/components/elm/bld/namelist_files/namelist_defaults.xml
index bbd68b2c0843..aebc38a9c978 100644
--- a/components/elm/bld/namelist_files/namelist_defaults.xml
+++ b/components/elm/bld/namelist_files/namelist_defaults.xml
@@ -661,6 +661,9 @@ this mask will have smb calculated over the entire global land surface
.false.
+
+.false.
+
lnd/clm2/snicardata/snicar_optics_5bnd_mam_c160322.nc
diff --git a/components/elm/bld/namelist_files/namelist_definition.xml b/components/elm/bld/namelist_files/namelist_definition.xml
index 81948f0a2821..fc4967b46876 100644
--- a/components/elm/bld/namelist_files/namelist_definition.xml
+++ b/components/elm/bld/namelist_files/namelist_definition.xml
@@ -760,6 +760,13 @@ If TRUE, TOP solar radiation parameterization will be activated, which considers
+
+
+If TRUE, fineTOP radiation parameterization will be activated, which considers the effects of fine(grid)-scale topography.
+
+
+
diff --git a/components/elm/cime_config/testdefs/testmods_dirs/elm/finetop_rad/shell_commands b/components/elm/cime_config/testdefs/testmods_dirs/elm/finetop_rad/shell_commands
new file mode 100644
index 000000000000..f27bbc407386
--- /dev/null
+++ b/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
\ No newline at end of file
diff --git a/components/elm/cime_config/testdefs/testmods_dirs/elm/finetop_rad/user_nl_elm b/components/elm/cime_config/testdefs/testmods_dirs/elm/finetop_rad/user_nl_elm
new file mode 100644
index 000000000000..f91db8e1da33
--- /dev/null
+++ b/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.
diff --git a/components/elm/cime_config/testdefs/testmods_dirs/elm/snowveg_arctic/shell_commands b/components/elm/cime_config/testdefs/testmods_dirs/elm/snowveg_arctic/shell_commands
index 396d8cea8d66..4b467e35f9ac 100644
--- a/components/elm/cime_config/testdefs/testmods_dirs/elm/snowveg_arctic/shell_commands
+++ b/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
\ No newline at end of file
diff --git a/components/elm/src/biogeophys/CanopyFluxesMod.F90 b/components/elm/src/biogeophys/CanopyFluxesMod.F90
index 47d4558990f9..1c278c75ce97 100755
--- a/components/elm/src/biogeophys/CanopyFluxesMod.F90
+++ b/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)
diff --git a/components/elm/src/biogeophys/CanopyHydrologyMod.F90 b/components/elm/src/biogeophys/CanopyHydrologyMod.F90
index 343380fa18fb..a6622b6cbd3e 100755
--- a/components/elm/src/biogeophys/CanopyHydrologyMod.F90
+++ b/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
diff --git a/components/elm/src/biogeophys/CanopyTemperatureMod.F90 b/components/elm/src/biogeophys/CanopyTemperatureMod.F90
index 017c3770a91d..22e4e1b2d6da 100644
--- a/components/elm/src/biogeophys/CanopyTemperatureMod.F90
+++ b/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)
diff --git a/components/elm/src/biogeophys/SnowSnicarMod.F90 b/components/elm/src/biogeophys/SnowSnicarMod.F90
index de161c1d4c9d..157f540c8256 100644
--- a/components/elm/src/biogeophys/SnowSnicarMod.F90
+++ b/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++
!
diff --git a/components/elm/src/biogeophys/SoilFluxesMod.F90 b/components/elm/src/biogeophys/SoilFluxesMod.F90
index 1b221a617a83..14ed1b4cac06 100644
--- a/components/elm/src/biogeophys/SoilFluxesMod.F90
+++ b/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
diff --git a/components/elm/src/biogeophys/SoilTemperatureMod.F90 b/components/elm/src/biogeophys/SoilTemperatureMod.F90
index 6826b3e22e63..79b452af41b2 100644
--- a/components/elm/src/biogeophys/SoilTemperatureMod.F90
+++ b/components/elm/src/biogeophys/SoilTemperatureMod.F90
@@ -12,7 +12,7 @@ module SoilTemperatureMod
use shr_infnan_mod , only : nan => shr_infnan_nan
use decompMod , only : bounds_type
use abortutils , only : endrun
- use elm_varctl , only : iulog
+ use elm_varctl , only : iulog, use_finetop_rad
use elm_varcon , only : spval
use UrbanParamsType , only : urbanparams_type
use atm2lndType , only : atm2lnd_type
@@ -32,6 +32,8 @@ module SoilTemperatureMod
use ExternalModelConstants , only : EM_ID_PTM
use ExternalModelConstants , only : EM_PTM_TBASED_SOLVE_STAGE
use ExternalModelInterfaceMod, only : EMI_Driver
+ use shr_const_mod , only : SHR_CONST_PI
+ use GridcellType , only : grc_pp
!! Needed beacuse EMI is still using them as arguments
use WaterstateType , only : waterstate_type
@@ -1750,6 +1752,7 @@ subroutine ComputeGroundHeatFluxAndDeriv(bounds, num_nolakec, filter_nolakec, &
real(r8) :: eflx_gnet_snow !
real(r8) :: eflx_gnet_soil !
real(r8) :: eflx_gnet_h2osfc !
+ real(r8) :: slope_rad, deg2rad
!-----------------------------------------------------------------------
! Enforce expected array sizes
@@ -1800,7 +1803,8 @@ subroutine ComputeGroundHeatFluxAndDeriv(bounds, num_nolakec, filter_nolakec, &
sabg_lyr => solarabs_vars%sabg_lyr_patch , & ! Output: [real(r8) (:,:) ] absorbed solar radiation (pft,lyr) [W/m2]
begc => bounds%begc , & ! Input: [integer ] beginning column index
- endc => bounds%endc & ! Input: [integer ] ending column index
+ endc => bounds%endc , & ! Input: [integer ] ending column index
+ slope_deg => grc_pp%slope_deg &
)
! Net ground heat flux into the surface and its temperature derivative
@@ -1809,6 +1813,8 @@ subroutine ComputeGroundHeatFluxAndDeriv(bounds, num_nolakec, filter_nolakec, &
do fc = 1,num_nolakec
c = filter_nolakec(fc)
+ g = col_pp%gridcell(c)
+ l = col_pp%landunit(c)
lwrad_emit(c) = emg(c) * sb * t_grnd(c)**4
dlwrad_emit(c) = 4._r8*emg(c) * sb * t_grnd(c)**3
@@ -1816,6 +1822,16 @@ subroutine ComputeGroundHeatFluxAndDeriv(bounds, num_nolakec, filter_nolakec, &
lwrad_emit_snow(c) = emg(c) * sb * t_soisno(c,snl(c)+1)**4
lwrad_emit_soil(c) = emg(c) * sb * t_soisno(c,1)**4
lwrad_emit_h2osfc(c) = emg(c) * sb * t_h2osfc(c)**4
+
+ if (use_finetop_rad .and. (.not. lun_pp%urbpoi(l))) then
+ deg2rad = SHR_CONST_PI/180._r8
+ slope_rad = slope_deg(g) * deg2rad
+ lwrad_emit(c) = lwrad_emit(c) / cos(slope_rad)
+ dlwrad_emit(c) = dlwrad_emit(c) / cos(slope_rad)
+ lwrad_emit_snow(c) = lwrad_emit_snow(c) / cos(slope_rad)
+ lwrad_emit_soil(c) = lwrad_emit_soil(c) / cos(slope_rad)
+ lwrad_emit_h2osfc(c) = lwrad_emit_h2osfc(c) / cos(slope_rad)
+ endif
end do
hs_soil(begc:endc) = 0._r8
diff --git a/components/elm/src/biogeophys/SolarAbsorbedType.F90 b/components/elm/src/biogeophys/SolarAbsorbedType.F90
index 75e985b858cd..03f26403401d 100644
--- a/components/elm/src/biogeophys/SolarAbsorbedType.F90
+++ b/components/elm/src/biogeophys/SolarAbsorbedType.F90
@@ -54,6 +54,8 @@ module SolarAbsorbedType
real(r8), pointer :: fsr_nir_d_patch (:) => null() ! patch reflected direct beam nir solar radiation (W/m**2)
real(r8), pointer :: fsr_nir_i_patch (:) => null() ! patch reflected diffuse nir solar radiation (W/m**2)
real(r8), pointer :: fsr_nir_d_ln_patch (:) => null() ! patch reflected direct beam nir solar radiation at local noon (W/m**2)
+ real(r8), pointer :: fsr_vis_d_patch (:) => null() ! patch reflected direct beam vis solar radiation (W/m**2)
+ real(r8), pointer :: fsr_vis_i_patch (:) => null() ! patch reflected diffuse vis solar radiation (W/m**2)
contains
@@ -131,6 +133,8 @@ subroutine InitAllocate(this, bounds)
allocate(this%fsds_nir_d_patch (begp:endp)) ; this%fsds_nir_d_patch (:) = spval
allocate(this%fsds_nir_i_patch (begp:endp)) ; this%fsds_nir_i_patch (:) = spval
allocate(this%fsds_nir_d_ln_patch (begp:endp)) ; this%fsds_nir_d_ln_patch (:) = spval
+ allocate(this%fsr_vis_d_patch (begp:endp)) ; this%fsr_vis_d_patch (:) = spval
+ allocate(this%fsr_vis_i_patch (begp:endp)) ; this%fsr_vis_i_patch (:) = spval
end subroutine InitAllocate
@@ -256,9 +260,11 @@ subroutine InitCold(this, bounds)
!
! !LOCAL VARIABLES:
integer :: begl, endl
+ integer :: begp, endp
!-----------------------------------------------------------------------
-
+
begl = bounds%begl; endl = bounds%endl
+ begp = bounds%begp; endp= bounds%endp
this%sabs_roof_dir_lun (begl:endl, :) = 0._r8
this%sabs_roof_dif_lun (begl:endl, :) = 0._r8
@@ -270,6 +276,10 @@ subroutine InitCold(this, bounds)
this%sabs_improad_dif_lun (begl:endl, :) = 0._r8
this%sabs_perroad_dir_lun (begl:endl, :) = 0._r8
this%sabs_perroad_dif_lun (begl:endl, :) = 0._r8
+ this%fsr_nir_d_patch (begp:endp) = 0._r8
+ this%fsr_nir_i_patch (begp:endp) = 0._r8
+ this%fsr_vis_d_patch (begp:endp) = 0._r8
+ this%fsr_vis_i_patch (begp:endp) = 0._r8
end subroutine InitCold
@@ -281,7 +291,7 @@ subroutine Restart(this, bounds, ncid, flag)
!
! !USES:
use shr_infnan_mod , only : shr_infnan_isnan
- use elm_varctl , only : use_snicar_frc, iulog
+ use elm_varctl , only : use_snicar_frc, iulog, use_finetop_rad
use spmdMod , only : masterproc
use abortutils , only : endrun
use ncdio_pio , only : file_desc_t, ncd_defvar, ncd_io, ncd_double, ncd_int, ncd_inqvdlen
@@ -347,6 +357,28 @@ subroutine Restart(this, bounds, ncid, flag)
dim2name='numrad', switchdim=.true., &
long_name='diffuse solar absorbed by pervious road per unit ground area per unit incident flux', units='', &
interpinic_flag='interp', readvar=readvar, data=this%sabs_perroad_dif_lun)
+
+ if (use_finetop_rad) then
+ call restartvar(ncid=ncid, flag=flag, varname='fsr_nir_d', xtype=ncd_double, &
+ dim1name='pft', &
+ long_name='direct nir reflected solar radiation', units='W/m^2', &
+ interpinic_flag='interp', readvar=readvar, data=this%fsr_nir_d_patch)
+
+ call restartvar(ncid=ncid, flag=flag, varname='fsr_nir_i', xtype=ncd_double, &
+ dim1name='pft', &
+ long_name='diffuse nir reflected solar radiation', units='W/m^2', &
+ interpinic_flag='interp', readvar=readvar, data=this%fsr_nir_i_patch)
+
+ call restartvar(ncid=ncid, flag=flag, varname='fsr_vis_d', xtype=ncd_double, &
+ dim1name='pft', &
+ long_name='direct vis reflected solar radiation', units='W/m^2', &
+ interpinic_flag='interp', readvar=readvar, data=this%fsr_vis_d_patch)
+
+ call restartvar(ncid=ncid, flag=flag, varname='fsr_vis_i', xtype=ncd_double, &
+ dim1name='pft', &
+ long_name='diffuse vis reflected solar radiation', units='W/m^2', &
+ interpinic_flag='interp', readvar=readvar, data=this%fsr_vis_i_patch)
+ end if
end subroutine Restart
diff --git a/components/elm/src/biogeophys/SurfaceAlbedoMod.F90 b/components/elm/src/biogeophys/SurfaceAlbedoMod.F90
index 4aed983b0ae0..95faf476a117 100644
--- a/components/elm/src/biogeophys/SurfaceAlbedoMod.F90
+++ b/components/elm/src/biogeophys/SurfaceAlbedoMod.F90
@@ -14,7 +14,7 @@ module SurfaceAlbedoMod
use landunit_varcon , only : istsoil, istcrop, istdlak
use elm_varcon , only : grlnd, namep, namet
use elm_varpar , only : numrad, nlevcan, nlevsno, nlevcan
- use elm_varctl , only : fsurdat, iulog, subgridflag, use_snicar_frc, use_fates, use_snicar_ad, use_top_solar_rad
+ use elm_varctl , only : fsurdat, iulog, subgridflag, use_snicar_frc, use_fates, use_snicar_ad, use_top_solar_rad, use_finetop_rad
use VegetationPropertiesType , only : veg_vp
use SnowSnicarMod , only : sno_nbr_aer, SNICAR_RT, SNICAR_AD_RT, DO_SNO_AER, DO_SNO_OC
use AerosolType , only : aerosol_type
@@ -93,6 +93,7 @@ subroutine SurfaceAlbedo(bounds, &
!$acc routine seq
use elm_varctl , only : iulog, subgridflag, use_snicar_frc, use_fates, use_snicar_ad, use_top_solar_rad
use shr_orb_mod
+ use shr_const_mod , only : SHR_CONST_PI
!
! !ARGUMENTS:
@@ -134,8 +135,12 @@ subroutine SurfaceAlbedo(bounds, &
real(r8) :: ws (bounds%begp:bounds%endp) ! fraction of LAI+SAI that is SAI
real(r8) :: blai(bounds%begp:bounds%endp) ! lai buried by snow: tlai - elai
real(r8) :: bsai(bounds%begp:bounds%endp) ! sai buried by snow: tsai - esai
+ real(r8) :: sza,saa,deg2rad,slope_rad,aspect_rad
+
real(r8) :: coszen_gcell (bounds%begg:bounds%endg) ! cosine solar zenith angle for next time step (grc)
real(r8) :: coszen_patch (bounds%begp:bounds%endp) ! cosine solar zenith angle for next time step (pft)
+ real(r8) :: cosinc_gcell (bounds%begg:bounds%endg) ! cosine solar incidence angle to local surface for next time step (grc)
+ real(r8) :: cosinc_patch (bounds%begp:bounds%endp) ! cosine solar incidence angle to local surface for next time step (pft)
real(r8) :: rho(bounds%begp:bounds%endp,numrad) ! leaf/stem refl weighted by fraction LAI and SAI
real(r8) :: tau(bounds%begp:bounds%endp,numrad) ! leaf/stem tran weighted by fraction LAI and SAI
real(r8) :: albsfc (bounds%begc:bounds%endc,numrad) ! albedo of surface underneath snow (col,bnd)
@@ -198,6 +203,7 @@ subroutine SurfaceAlbedo(bounds, &
ncan => surfalb_vars%ncan_patch , & ! Output: [integer (:) ] number of canopy layers
nrad => surfalb_vars%nrad_patch , & ! Output: [integer (:) ] number of canopy layers, above snow for radiative transfer
coszen_col => surfalb_vars%coszen_col , & ! Output: [real(r8) (:) ] cosine of solar zenith angle
+ cosinc_col => surfalb_vars%cosinc_col , & ! Output: [real(r8) (:) ] cosine of solar zenith angle
albgrd => surfalb_vars%albgrd_col , & ! Output: [real(r8) (:,:) ] ground albedo (direct)
albgri => surfalb_vars%albgri_col , & ! Output: [real(r8) (:,:) ] ground albedo (diffuse)
albsod => surfalb_vars%albsod_col , & ! Output: [real(r8) (:,:) ] direct-beam soil albedo (col,bnd) [frc]
@@ -235,17 +241,38 @@ subroutine SurfaceAlbedo(bounds, &
! Cosine solar zenith angle for next time step
- do g = bounds%begg,bounds%endg
- coszen_gcell(g) = shr_orb_cosz (nextsw_cday, grc_pp%lat(g), grc_pp%lon(g), declinp1)
- end do
+ deg2rad = SHR_CONST_PI/180._r8
+ if (.not. use_finetop_rad) then
+ do g = bounds%begg,bounds%endg
+ coszen_gcell(g) = shr_orb_cosz (nextsw_cday, grc_pp%lat(g), grc_pp%lon(g), declinp1)
+ cosinc_gcell(g) = coszen_gcell(g)
+ end do
+ else
+ do g = bounds%begg,bounds%endg
+ coszen_gcell(g) = shr_orb_cosz (nextsw_cday, grc_pp%lat(g), grc_pp%lon(g), declinp1)
+ sza = acos(coszen_gcell(g)) ! solar zenith angle
+ saa = shr_orb_azimuth(nextsw_cday, grc_pp%lat(g), grc_pp%lon(g), declinp1, sza) ! solar azimuth angle
+
+ slope_rad = grc_pp%slope_deg(g) * deg2rad
+ aspect_rad = grc_pp%aspect_deg(g) * deg2rad
+
+ cosinc_gcell(g) = cos(slope_rad) * coszen_gcell(g) + sin(slope_rad) * sin(sza) * cos(aspect_rad - saa)
+ cosinc_gcell(g) = max(-1._r8, min(cosinc_gcell(g), 1._r8))
+
+ if (cosinc_gcell(g) <= 0._r8) cosinc_gcell(g) = 0.1_r8 ! although direct solar radiation is zero, we need to calculate diffuse albedo in this case
+ end do
+ endif
+
do c = bounds%begc,bounds%endc
g = col_pp%gridcell(c)
coszen_col(c) = coszen_gcell(g)
+ cosinc_col(c) = cosinc_gcell(g)
end do
do fp = 1,num_nourbanp
p = filter_nourbanp(fp)
g = veg_pp%gridcell(p)
- coszen_patch(p) = coszen_gcell(g)
+ coszen_patch(p) = coszen_gcell(g)
+ cosinc_patch(p) = cosinc_gcell(g)
end do
! Initialize output because solar radiation only done if coszen > 0
@@ -305,7 +332,7 @@ subroutine SurfaceAlbedo(bounds, &
call SoilAlbedo(bounds, &
num_nourbanc, filter_nourbanc, &
- coszen_col(bounds%begc:bounds%endc), &
+ cosinc_col(bounds%begc:bounds%endc), &
albsnd(bounds%begc:bounds%endc, :), &
albsni(bounds%begc:bounds%endc, :), &
lakestate_vars, surfalb_vars)
@@ -373,7 +400,7 @@ subroutine SurfaceAlbedo(bounds, &
flg_slr = 1; ! direct-beam
if (use_snicar_ad) then
call SNICAR_AD_RT(flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, &
- coszen_col(bounds%begc:bounds%endc), &
+ cosinc_col(bounds%begc:bounds%endc), &
flg_slr, &
h2osno_liq(bounds%begc:bounds%endc, :), &
h2osno_ice(bounds%begc:bounds%endc, :), &
@@ -384,7 +411,7 @@ subroutine SurfaceAlbedo(bounds, &
foo_snw(bounds%begc:bounds%endc, :, :) )
else
call SNICAR_RT(flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, &
- coszen_col(bounds%begc:bounds%endc), &
+ cosinc_col(bounds%begc:bounds%endc), &
flg_slr, &
h2osno_liq(bounds%begc:bounds%endc, :), &
h2osno_ice(bounds%begc:bounds%endc, :), &
@@ -399,7 +426,7 @@ subroutine SurfaceAlbedo(bounds, &
flg_slr = 2; ! diffuse
if (use_snicar_ad) then
call SNICAR_AD_RT(flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, &
- coszen_col(bounds%begc:bounds%endc), &
+ cosinc_col(bounds%begc:bounds%endc), &
flg_slr, &
h2osno_liq(bounds%begc:bounds%endc, :), &
h2osno_ice(bounds%begc:bounds%endc, :), &
@@ -410,7 +437,7 @@ subroutine SurfaceAlbedo(bounds, &
foo_snw(bounds%begc:bounds%endc, :, :) )
else
call SNICAR_RT(flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, &
- coszen_col(bounds%begc:bounds%endc), &
+ cosinc_col(bounds%begc:bounds%endc), &
flg_slr, &
h2osno_liq(bounds%begc:bounds%endc, :), &
h2osno_ice(bounds%begc:bounds%endc, :), &
@@ -435,7 +462,7 @@ subroutine SurfaceAlbedo(bounds, &
flg_slr = 1; ! direct-beam
if (use_snicar_ad) then
call SNICAR_AD_RT(flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, &
- coszen_col(bounds%begc:bounds%endc), &
+ cosinc_col(bounds%begc:bounds%endc), &
flg_slr, &
h2osno_liq(bounds%begc:bounds%endc, :), &
h2osno_ice(bounds%begc:bounds%endc, :), &
@@ -446,7 +473,7 @@ subroutine SurfaceAlbedo(bounds, &
foo_snw(bounds%begc:bounds%endc, :, :) )
else
call SNICAR_RT(flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, &
- coszen_col(bounds%begc:bounds%endc), &
+ cosinc_col(bounds%begc:bounds%endc), &
flg_slr, &
h2osno_liq(bounds%begc:bounds%endc, :), &
h2osno_ice(bounds%begc:bounds%endc, :), &
@@ -460,7 +487,7 @@ subroutine SurfaceAlbedo(bounds, &
flg_slr = 2; ! diffuse
if (use_snicar_ad) then
call SNICAR_AD_RT(flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, &
- coszen_col(bounds%begc:bounds%endc), &
+ cosinc_col(bounds%begc:bounds%endc), &
flg_slr, &
h2osno_liq(bounds%begc:bounds%endc, :), &
h2osno_ice(bounds%begc:bounds%endc, :), &
@@ -471,7 +498,7 @@ subroutine SurfaceAlbedo(bounds, &
foo_snw(bounds%begc:bounds%endc, :, :) )
else
call SNICAR_RT(flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, &
- coszen_col(bounds%begc:bounds%endc), &
+ cosinc_col(bounds%begc:bounds%endc), &
flg_slr, &
h2osno_liq(bounds%begc:bounds%endc, :), &
h2osno_ice(bounds%begc:bounds%endc, :), &
@@ -496,7 +523,7 @@ subroutine SurfaceAlbedo(bounds, &
flg_slr = 1; ! direct-beam
if (use_snicar_ad) then
call SNICAR_AD_RT(flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, &
- coszen_col(bounds%begc:bounds%endc), &
+ cosinc_col(bounds%begc:bounds%endc), &
flg_slr, &
h2osno_liq(bounds%begc:bounds%endc, :), &
h2osno_ice(bounds%begc:bounds%endc, :), &
@@ -507,7 +534,7 @@ subroutine SurfaceAlbedo(bounds, &
foo_snw(bounds%begc:bounds%endc, :, :) )
else
call SNICAR_RT(flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, &
- coszen_col(bounds%begc:bounds%endc), &
+ cosinc_col(bounds%begc:bounds%endc), &
flg_slr, &
h2osno_liq(bounds%begc:bounds%endc, :), &
h2osno_ice(bounds%begc:bounds%endc, :), &
@@ -521,7 +548,7 @@ subroutine SurfaceAlbedo(bounds, &
flg_slr = 2; ! diffuse
if (use_snicar_ad) then
call SNICAR_AD_RT(flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, &
- coszen_col(bounds%begc:bounds%endc), &
+ cosinc_col(bounds%begc:bounds%endc), &
flg_slr, &
h2osno_liq(bounds%begc:bounds%endc, :), &
h2osno_ice(bounds%begc:bounds%endc, :), &
@@ -532,7 +559,7 @@ subroutine SurfaceAlbedo(bounds, &
foo_snw(bounds%begc:bounds%endc, :, :) )
else
call SNICAR_RT(flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, &
- coszen_col(bounds%begc:bounds%endc), &
+ cosinc_col(bounds%begc:bounds%endc), &
flg_slr, &
h2osno_liq(bounds%begc:bounds%endc, :), &
h2osno_ice(bounds%begc:bounds%endc, :), &
@@ -548,7 +575,7 @@ subroutine SurfaceAlbedo(bounds, &
flg_slr = 1; ! direct-beam
if (use_snicar_ad) then
call SNICAR_AD_RT(flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, &
- coszen_col(bounds%begc:bounds%endc), &
+ cosinc_col(bounds%begc:bounds%endc), &
flg_slr, &
h2osno_liq(bounds%begc:bounds%endc, :), &
h2osno_ice(bounds%begc:bounds%endc, :), &
@@ -559,7 +586,7 @@ subroutine SurfaceAlbedo(bounds, &
foo_snw(bounds%begc:bounds%endc, :, :) )
else
call SNICAR_RT(flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, &
- coszen_col(bounds%begc:bounds%endc), &
+ cosinc_col(bounds%begc:bounds%endc), &
flg_slr, &
h2osno_liq(bounds%begc:bounds%endc, :), &
h2osno_ice(bounds%begc:bounds%endc, :), &
@@ -573,7 +600,7 @@ subroutine SurfaceAlbedo(bounds, &
flg_slr = 2; ! diffuse
if (use_snicar_ad) then
call SNICAR_AD_RT(flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, &
- coszen_col(bounds%begc:bounds%endc), &
+ cosinc_col(bounds%begc:bounds%endc), &
flg_slr, &
h2osno_liq(bounds%begc:bounds%endc, :), &
h2osno_ice(bounds%begc:bounds%endc, :), &
@@ -584,7 +611,7 @@ subroutine SurfaceAlbedo(bounds, &
foo_snw(bounds%begc:bounds%endc, :, :) )
else
call SNICAR_RT(flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, &
- coszen_col(bounds%begc:bounds%endc), &
+ cosinc_col(bounds%begc:bounds%endc), &
flg_slr, &
h2osno_liq(bounds%begc:bounds%endc, :), &
h2osno_ice(bounds%begc:bounds%endc, :), &
@@ -601,7 +628,7 @@ subroutine SurfaceAlbedo(bounds, &
flg_slr = 1; ! direct-beam
if (use_snicar_ad) then
call SNICAR_AD_RT(flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, &
- coszen_col(bounds%begc:bounds%endc), &
+ cosinc_col(bounds%begc:bounds%endc), &
flg_slr, &
h2osno_liq(bounds%begc:bounds%endc, :), &
h2osno_ice(bounds%begc:bounds%endc, :), &
@@ -612,7 +639,7 @@ subroutine SurfaceAlbedo(bounds, &
flx_absd_snw(bounds%begc:bounds%endc, :, :) )
else
call SNICAR_RT(flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, &
- coszen_col(bounds%begc:bounds%endc), &
+ cosinc_col(bounds%begc:bounds%endc), &
flg_slr, &
h2osno_liq(bounds%begc:bounds%endc, :), &
h2osno_ice(bounds%begc:bounds%endc, :), &
@@ -626,7 +653,7 @@ subroutine SurfaceAlbedo(bounds, &
flg_slr = 2; ! diffuse
if (use_snicar_ad) then
call SNICAR_AD_RT(flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, &
- coszen_col(bounds%begc:bounds%endc), &
+ cosinc_col(bounds%begc:bounds%endc), &
flg_slr, &
h2osno_liq(bounds%begc:bounds%endc, :), &
h2osno_ice(bounds%begc:bounds%endc, :), &
@@ -637,7 +664,7 @@ subroutine SurfaceAlbedo(bounds, &
flx_absi_snw(bounds%begc:bounds%endc, :, :) )
else
call SNICAR_RT(flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, &
- coszen_col(bounds%begc:bounds%endc), &
+ cosinc_col(bounds%begc:bounds%endc), &
flg_slr, &
h2osno_liq(bounds%begc:bounds%endc, :), &
h2osno_ice(bounds%begc:bounds%endc, :), &
@@ -947,12 +974,11 @@ subroutine SurfaceAlbedo(bounds, &
endif
call TwoStream (bounds, filter_vegsol, num_vegsol, &
- coszen_patch(bounds%begp:bounds%endp), &
+ cosinc_patch(bounds%begp:bounds%endp), &
rho(bounds%begp:bounds%endp, :), &
tau(bounds%begp:bounds%endp, :), &
canopystate_vars, surfalb_vars, &
nextsw_cday, declinp1)
-
endif
! Determine values for non-vegetated patches where coszen > 0
@@ -980,7 +1006,7 @@ subroutine SurfaceAlbedo(bounds, &
coszen_patch(bounds%begp:bounds%endp), declinp1, surfalb_vars, .false.)
endif
- end associate
+ end associate
end subroutine SurfaceAlbedo
@@ -1138,6 +1164,7 @@ subroutine TwoStream(bounds, &
use elm_varpar, only : numrad, nlevcan
use elm_varcon, only : omegas, tfrz, betads, betais
use elm_varctl, only : iulog, use_top_solar_rad
+ use shr_const_mod , only : SHR_CONST_PI
!
! !ARGUMENTS:
type(bounds_type) , intent(in) :: bounds
@@ -1153,7 +1180,7 @@ subroutine TwoStream(bounds, &
!
! !LOCAL VARIABLES:
- integer :: fp,p,c,iv ! array indices
+ integer :: fp,p,c,iv,g ! array indices
integer :: ib ! waveband number
real(r8) :: cosz ! 0.001 <= coszen <= 1.000
real(r8) :: asu ! single scattering albedo
@@ -1187,6 +1214,9 @@ subroutine TwoStream(bounds, &
real(r8) :: laisum ! cumulative lai+sai for canopy layer (at middle of layer)
real(r8) :: extkb ! direct beam extinction coefficient
real(r8) :: extkn ! nitrogen allocation coefficient
+ real(r8) :: deg2rad,slope_rad
+ real(r8) :: elaislope(bounds%begp:bounds%endp) ! elai projected to slope
+ real(r8) :: esaislope(bounds%begp:bounds%endp) ! esai projected to slope
!-----------------------------------------------------------------------
! Enforce expected array sizes
@@ -1234,11 +1264,23 @@ subroutine TwoStream(bounds, &
ftii => surfalb_vars%ftii_patch & ! Output: [real(r8) (:,:) ] down diffuse flux below canopy per unit diffuse flx
)
+ deg2rad = SHR_CONST_PI/180._r8
! Calculate two-stream parameters that are independent of waveband:
! chil, gdir, twostext, avmu, and temp0 and temp2 (used for asu)
do fp = 1,num_vegsol
p = filter_vegsol(fp)
+ g = veg_pp%gridcell(p)
+
+ slope_rad = grc_pp%slope_deg(g) * deg2rad
+
+ if (use_finetop_rad) then
+ elaislope(p) = elai(p) * cos(slope_rad)
+ esaislope(p) = esai(p) * cos(slope_rad)
+ else
+ elaislope(p) = elai(p)
+ esaislope(p) = esai(p)
+ endif
! note that the following limit only acts on cosz values > 0 and less than
! 0.001, not on values cosz = 0, since these zero have already been filtered
@@ -1343,9 +1385,9 @@ subroutine TwoStream(bounds, &
! Absorbed, reflected, transmitted fluxes per unit incoming radiation
! for full canopy
- t1 = min(h*(elai(p)+esai(p)), 40._r8)
+ t1 = min(h*(elaislope(p)+esaislope(p)), 40._r8)
s1 = exp(-t1)
- t1 = min(twostext(p)*(elai(p)+esai(p)), 40._r8)
+ t1 = min(twostext(p)*(elaislope(p)+esaislope(p)), 40._r8)
s2 = exp(-t1)
! Direct beam
@@ -1468,7 +1510,7 @@ subroutine TwoStream(bounds, &
fsun_z(p,1) = (1._r8 - s2) / t1
! absorbed PAR (per unit sun/shade lai+sai)
- laisum = elai(p)+esai(p)
+ laisum = elaislope(p)+esaislope(p)
fabd_sun_z(p,1) = fabd_sun(p,ib) / (fsun_z(p,1)*laisum)
fabi_sun_z(p,1) = fabi_sun(p,ib) / (fsun_z(p,1)*laisum)
fabd_sha_z(p,1) = fabd_sha(p,ib) / ((1._r8 - fsun_z(p,1))*laisum)
@@ -1477,11 +1519,11 @@ subroutine TwoStream(bounds, &
! leaf to canopy scaling coefficients
extkn = 0.30_r8
extkb = twostext(p)
- vcmaxcintsun(p) = (1._r8 - exp(-(extkn+extkb)*elai(p))) / (extkn + extkb)
- vcmaxcintsha(p) = (1._r8 - exp(-extkn*elai(p))) / extkn - vcmaxcintsun(p)
- if (elai(p) > 0._r8) then
- vcmaxcintsun(p) = vcmaxcintsun(p) / (fsun_z(p,1)*elai(p))
- vcmaxcintsha(p) = vcmaxcintsha(p) / ((1._r8 - fsun_z(p,1))*elai(p))
+ vcmaxcintsun(p) = (1._r8 - exp(-(extkn+extkb)*elaislope(p))) / (extkn + extkb)
+ vcmaxcintsha(p) = (1._r8 - exp(-extkn*elaislope(p))) / extkn - vcmaxcintsun(p)
+ if (elaislope(p) > 0._r8) then
+ vcmaxcintsun(p) = vcmaxcintsun(p) / (fsun_z(p,1)*elaislope(p))
+ vcmaxcintsha(p) = vcmaxcintsha(p) / ((1._r8 - fsun_z(p,1))*elaislope(p))
else
vcmaxcintsun(p) = 0._r8
vcmaxcintsha(p) = 0._r8
diff --git a/components/elm/src/biogeophys/SurfaceAlbedoType.F90 b/components/elm/src/biogeophys/SurfaceAlbedoType.F90
index 10c148df7bff..c57a0efdaf9b 100644
--- a/components/elm/src/biogeophys/SurfaceAlbedoType.F90
+++ b/components/elm/src/biogeophys/SurfaceAlbedoType.F90
@@ -8,7 +8,7 @@ module SurfaceAlbedoType
use decompMod , only : bounds_type
use elm_varpar , only : numrad, nlevcan, nlevsno
use abortutils , only : endrun
- use elm_varctl , only : fsurdat, iulog
+ use elm_varctl , only : fsurdat, iulog, use_finetop_rad
use elm_varcon , only : grlnd
use ColumnType , only : col_pp
@@ -59,6 +59,7 @@ module SurfaceAlbedoType
type, public :: surfalb_type
real(r8), pointer :: coszen_col (:) => null() ! col cosine of solar zenith angle
+ real(r8), pointer :: cosinc_col (:) => null() ! col cosine of solar incident angle (local)
real(r8), pointer :: albd_patch (:,:) => null() ! patch surface albedo (direct) (numrad)
real(r8), pointer :: albi_patch (:,:) => null() ! patch surface albedo (diffuse) (numrad)
real(r8), pointer :: albgrd_pur_col (:,:) => null() ! col pure snow ground direct albedo (numrad)
@@ -255,6 +256,7 @@ subroutine InitAllocate(this, bounds)
begc = bounds%begc; endc = bounds%endc
allocate(this%coszen_col (begc:endc)) ; this%coszen_col (:) = spval
+ allocate(this%cosinc_col (begc:endc)) ; this%cosinc_col (:) = spval
allocate(this%albgrd_col (begc:endc,numrad)) ; this%albgrd_col (:,:) =spval
allocate(this%albgri_col (begc:endc,numrad)) ; this%albgri_col (:,:) =spval
allocate(this%albsnd_hst_col (begc:endc,numrad)) ; this%albsnd_hst_col (:,:) = spval
@@ -334,6 +336,11 @@ subroutine InitHistory(this, bounds)
avgflag='A', long_name='cosine of solar zenith angle', &
ptr_col=this%coszen_col, default='inactive')
+ this%cosinc_col(begc:endc) = spval
+ call hist_addfld1d (fname='COSINC', units='1', &
+ avgflag='A', long_name='cosine of solar incident angle (local surface)', &
+ ptr_col=this%cosinc_col, default='inactive')
+
this%albgri_col(begc:endc,:) = spval
call hist_addfld2d (fname='ALBGRD', units='proportion', type2d='numrad', &
avgflag='A', long_name='ground albedo (direct)', &
@@ -450,6 +457,13 @@ subroutine Restart(this, bounds, ncid, flag, &
long_name='cosine of solar zenith angle', units='1', &
interpinic_flag='interp', readvar=readvar, data=this%coszen_col)
+ if (use_finetop_rad) then
+ call restartvar(ncid=ncid, flag=flag, varname='cosinc', xtype=ncd_double, &
+ dim1name='column', &
+ long_name='cosine of solar incident angle (local surface)', units='1', &
+ interpinic_flag='interp', readvar=readvar, data=this%cosinc_col)
+ end if
+
call restartvar(ncid=ncid, flag=flag, varname='albd', xtype=ncd_double, &
dim1name='pft', dim2name='numrad', switchdim=.true., &
long_name='surface albedo (direct) (0 to 1)', units='', &
diff --git a/components/elm/src/biogeophys/SurfaceRadiationMod.F90 b/components/elm/src/biogeophys/SurfaceRadiationMod.F90
index 47fe244c530e..68e3c2e917af 100644
--- a/components/elm/src/biogeophys/SurfaceRadiationMod.F90
+++ b/components/elm/src/biogeophys/SurfaceRadiationMod.F90
@@ -328,7 +328,7 @@ subroutine SurfaceRadiation(bounds, num_nourbanp, filter_nourbanp, &
use elm_varpar , only : numrad, nlevsno
use elm_varcon , only : spval, degpsec, isecspday
use landunit_varcon , only : istsoil, istcrop
- use elm_varctl , only : subgridflag, use_snicar_frc
+ use elm_varctl , only : subgridflag, use_snicar_frc, use_finetop_rad
use SnowSnicarMod , only : DO_SNO_OC
!
! !ARGUMENTS:
@@ -388,6 +388,8 @@ subroutine SurfaceRadiation(bounds, num_nourbanp, filter_nourbanp, &
forc_solad => top_af%solad , & ! Input: [real(r8) (:,:) ] direct beam radiation (W/m**2)
forc_solai => top_af%solai , & ! Input: [real(r8) (:,:) ] diffuse radiation (W/m**2)
+ forc_solad_pp => top_af%solad_pp , & ! Input: [real(r8) (:,:) ] direct beam radiation under PP (W/m**2)
+ forc_solai_pp => top_af%solai_pp , & ! Input: [real(r8) (:,:) ] diffuse radiation under PP (W/m**2)
snow_depth => col_ws%snow_depth , & ! Input: [real(r8) (:) ] snow height (m)
frac_sno => col_ws%frac_sno , & ! Input: [real(r8) (:) ] fraction of ground covered by snow (0 to 1)
@@ -475,6 +477,12 @@ subroutine SurfaceRadiation(bounds, num_nourbanp, filter_nourbanp, &
dtime = dtime_mod
secs = secs_curr
+
+ if (.not. use_finetop_rad) then
+ forc_solad_pp(:,:) = forc_solad(:,:)
+ forc_solai_pp(:,:) = forc_solai(:,:)
+ end if
+
! Initialize fluxes
do fp = 1,num_nourbanp
p = filter_nourbanp(fp)
@@ -736,6 +744,9 @@ subroutine SurfaceRadiation(bounds, num_nourbanp, filter_nourbanp, &
fsr_vis_i(p) = albi(p,1)*forc_solai(t,1)
fsr_nir_i(p) = albi(p,2)*forc_solai(t,2)
+ solarabs_vars%fsr_vis_d_patch(p) = fsr_vis_d(p)
+ solarabs_vars%fsr_vis_i_patch(p) = fsr_vis_i(p)
+
local_secp1 = secs + nint((grc_pp%londeg(g)/degpsec)/dtime)*dtime
local_secp1 = mod(local_secp1,isecspday)
if (local_secp1 == isecspday/2) then
@@ -794,16 +805,16 @@ subroutine SurfaceRadiation(bounds, num_nourbanp, filter_nourbanp, &
endif
! Solar incident
- fsds_vis_d(p) = forc_solad(t,1)
- fsds_nir_d(p) = forc_solad(t,2)
- fsds_vis_i(p) = forc_solai(t,1)
- fsds_nir_i(p) = forc_solai(t,2)
+ fsds_vis_d(p) = forc_solad_pp(t,1)
+ fsds_nir_d(p) = forc_solad_pp(t,2)
+ fsds_vis_i(p) = forc_solai_pp(t,1)
+ fsds_nir_i(p) = forc_solai_pp(t,2)
! Determine local noon incident solar
if (local_secp1 == noonsec) then
- fsds_vis_d_ln(p) = forc_solad(t,1)
- fsds_nir_d_ln(p) = forc_solad(t,2)
- fsds_vis_i_ln(p) = forc_solai(t,1)
+ fsds_vis_d_ln(p) = forc_solad_pp(t,1)
+ fsds_nir_d_ln(p) = forc_solad_pp(t,2)
+ fsds_vis_i_ln(p) = forc_solai_pp(t,1)
parveg_ln(p) = 0._r8
else
fsds_vis_d_ln(p) = spval
@@ -815,10 +826,13 @@ subroutine SurfaceRadiation(bounds, num_nourbanp, filter_nourbanp, &
! Solar reflected
! per unit ground area (roof, road) and per unit wall area (sunwall, shadewall)
- fsr_vis_d(p) = albd(p,1) * forc_solad(t,1)
- fsr_nir_d(p) = albd(p,2) * forc_solad(t,2)
- fsr_vis_i(p) = albi(p,1) * forc_solai(t,1)
- fsr_nir_i(p) = albi(p,2) * forc_solai(t,2)
+ fsr_vis_d(p) = albd(p,1) * forc_solad_pp(t,1)
+ fsr_nir_d(p) = albd(p,2) * forc_solad_pp(t,2)
+ fsr_vis_i(p) = albi(p,1) * forc_solai_pp(t,1)
+ fsr_nir_i(p) = albi(p,2) * forc_solai_pp(t,2)
+
+ solarabs_vars%fsr_vis_d_patch(p) = fsr_vis_d(p)
+ solarabs_vars%fsr_vis_i_patch(p) = fsr_vis_i(p)
! Determine local noon reflected solar
if (local_secp1 == noonsec) then
diff --git a/components/elm/src/biogeophys/UrbanRadiationMod.F90 b/components/elm/src/biogeophys/UrbanRadiationMod.F90
index 8cccc4e9ac32..22313fc85f29 100644
--- a/components/elm/src/biogeophys/UrbanRadiationMod.F90
+++ b/components/elm/src/biogeophys/UrbanRadiationMod.F90
@@ -13,7 +13,7 @@ module UrbanRadiationMod
use decompMod , only : bounds_type
use elm_varpar , only : numrad
use elm_varcon , only : isecspday, degpsec, namel
- use elm_varctl , only : iulog
+ use elm_varctl , only : iulog, use_finetop_rad
use abortutils , only : endrun
use UrbanParamsType , only : urbanparams_type
use atm2lndType , only : atm2lnd_type
@@ -116,10 +116,10 @@ subroutine UrbanRadiation (bounds , &
canyon_hwr => lun_pp%canyon_hwr , & ! Input: [real(r8) (:) ] ratio of building height to street width
wtroad_perv => lun_pp%wtroad_perv , & ! Input: [real(r8) (:) ] weight of pervious road wrt total road
- forc_solad => top_af%solad , & ! Input: [real(r8) (:,:) ] direct beam radiation (vis=forc_sols , nir=forc_soll ) (W/m**2)
- forc_solai => top_af%solai , & ! Input: [real(r8) (:,:) ] diffuse beam radiation (vis=forc_sols , nir=forc_soll ) (W/m**2)
- forc_solar => top_af%solar , & ! Input: [real(r8) (:) ] incident solar radiation (W/m**2)
- forc_lwrad => top_af%lwrad , & ! Input: [real(r8) (:) ] downward infrared (longwave) radiation (W/m**2)
+ forc_solad => top_af%solad_pp , & ! Input: [real(r8) (:,:) ] direct beam radiation under PP (vis=forc_sols , nir=forc_soll ) (W/m**2)
+ forc_solai => top_af%solai_pp , & ! Input: [real(r8) (:,:) ] diffuse beam radiation under PP (vis=forc_sols , nir=forc_soll ) (W/m**2)
+ forc_solar => top_af%solar_pp , & ! Input: [real(r8) (:) ] incident solar radiation under PP (W/m**2)
+ forc_lwrad => top_af%lwrad_pp , & ! Input: [real(r8) (:) ] downward infrared (longwave) radiation under PP (W/m**2)
frac_sno => col_ws%frac_sno , & ! Input: [real(r8) (:) ] fraction of ground covered by snow (0 to 1)
@@ -157,6 +157,13 @@ subroutine UrbanRadiation (bounds , &
endl => bounds%endl &
)
+ if (.not. use_finetop_rad) then
+ forc_solad(:,:) = top_af%solad(:,:)
+ forc_solai(:,:) = top_af%solai(:,:)
+ forc_solar(:) = top_af%solar(:)
+ forc_lwrad(:) = top_af%lwrad(:)
+ end if
+
! Define fields that appear on the restart file for non-urban landunits
do fl = 1,num_nourbanl
diff --git a/components/elm/src/cpl/lnd_comp_mct.F90 b/components/elm/src/cpl/lnd_comp_mct.F90
index 6b5a6dbe6b08..0e6d41c4f559 100644
--- a/components/elm/src/cpl/lnd_comp_mct.F90
+++ b/components/elm/src/cpl/lnd_comp_mct.F90
@@ -1046,7 +1046,7 @@ subroutine lnd_export_moab(EClock, bounds, lnd2atm_vars, lnd2glc_vars)
! !USES:
use shr_kind_mod , only : r8 => shr_kind_r8
use shr_kind_mod , only : CXX => SHR_KIND_CXX
- use elm_varctl , only : iulog, create_glacier_mec_landunit
+ use elm_varctl , only : iulog, create_glacier_mec_landunit, use_finetop_rad
use elm_time_manager , only : get_nstep, get_step_size
use domainMod , only : ldomain
use seq_drydep_mod , only : n_drydep
@@ -1083,10 +1083,19 @@ subroutine lnd_export_moab(EClock, bounds, lnd2atm_vars, lnd2glc_vars)
i = 1 + (g-bounds%begg)
l2x_lm(i,index_l2x_Sl_t) = lnd2atm_vars%t_rad_grc(g)
l2x_lm(i,index_l2x_Sl_snowh) = lnd2atm_vars%h2osno_grc(g)
- l2x_lm(i,index_l2x_Sl_avsdr) = lnd2atm_vars%albd_grc(g,1)
- l2x_lm(i,index_l2x_Sl_anidr) = lnd2atm_vars%albd_grc(g,2)
- l2x_lm(i,index_l2x_Sl_avsdf) = lnd2atm_vars%albi_grc(g,1)
- l2x_lm(i,index_l2x_Sl_anidf) = lnd2atm_vars%albi_grc(g,2)
+
+ if (use_finetop_rad) then
+ l2x_lm(i,index_l2x_Sl_avsdr) = lnd2atm_vars%apparent_albd_grc(g,1)
+ l2x_lm(i,index_l2x_Sl_anidr) = lnd2atm_vars%apparent_albd_grc(g,2)
+ l2x_lm(i,index_l2x_Sl_avsdf) = lnd2atm_vars%apparent_albi_grc(g,1)
+ l2x_lm(i,index_l2x_Sl_anidf) = lnd2atm_vars%apparent_albi_grc(g,2)
+ else
+ l2x_lm(i,index_l2x_Sl_avsdr) = lnd2atm_vars%albd_grc(g,1)
+ l2x_lm(i,index_l2x_Sl_anidr) = lnd2atm_vars%albd_grc(g,2)
+ l2x_lm(i,index_l2x_Sl_avsdf) = lnd2atm_vars%albi_grc(g,1)
+ l2x_lm(i,index_l2x_Sl_anidf) = lnd2atm_vars%albi_grc(g,2)
+ end if
+
l2x_lm(i,index_l2x_Sl_tref) = lnd2atm_vars%t_ref2m_grc(g)
l2x_lm(i,index_l2x_Sl_qref) = lnd2atm_vars%q_ref2m_grc(g)
l2x_lm(i,index_l2x_Sl_u10) = lnd2atm_vars%u_ref10m_grc(g)
diff --git a/components/elm/src/cpl/lnd_import_export.F90 b/components/elm/src/cpl/lnd_import_export.F90
index 75c6111da77f..20df1510a137 100644
--- a/components/elm/src/cpl/lnd_import_export.F90
+++ b/components/elm/src/cpl/lnd_import_export.F90
@@ -1403,7 +1403,7 @@ subroutine lnd_export( bounds, lnd2atm_vars, lnd2glc_vars, lnd2iac_vars, l2x)
!
! !USES:
use shr_kind_mod , only : r8 => shr_kind_r8
- use elm_varctl , only : iulog, create_glacier_mec_landunit, iac_present
+ use elm_varctl , only : iulog, create_glacier_mec_landunit, iac_present, use_finetop_rad
use elm_time_manager , only : get_nstep, get_step_size
use domainMod , only : ldomain
use seq_drydep_mod , only : n_drydep
@@ -1437,10 +1437,19 @@ subroutine lnd_export( bounds, lnd2atm_vars, lnd2glc_vars, lnd2iac_vars, l2x)
i = 1 + (g-bounds%begg)
l2x(index_l2x_Sl_t,i) = lnd2atm_vars%t_rad_grc(g)
l2x(index_l2x_Sl_snowh,i) = lnd2atm_vars%h2osno_grc(g)
- l2x(index_l2x_Sl_avsdr,i) = lnd2atm_vars%albd_grc(g,1)
- l2x(index_l2x_Sl_anidr,i) = lnd2atm_vars%albd_grc(g,2)
- l2x(index_l2x_Sl_avsdf,i) = lnd2atm_vars%albi_grc(g,1)
- l2x(index_l2x_Sl_anidf,i) = lnd2atm_vars%albi_grc(g,2)
+
+ if (use_finetop_rad) then
+ l2x(index_l2x_Sl_avsdr,i) = lnd2atm_vars%apparent_albd_grc(g,1)
+ l2x(index_l2x_Sl_anidr,i) = lnd2atm_vars%apparent_albd_grc(g,2)
+ l2x(index_l2x_Sl_avsdf,i) = lnd2atm_vars%apparent_albi_grc(g,1)
+ l2x(index_l2x_Sl_anidf,i) = lnd2atm_vars%apparent_albi_grc(g,2)
+ else
+ l2x(index_l2x_Sl_avsdr,i) = lnd2atm_vars%albd_grc(g,1)
+ l2x(index_l2x_Sl_anidr,i) = lnd2atm_vars%albd_grc(g,2)
+ l2x(index_l2x_Sl_avsdf,i) = lnd2atm_vars%albi_grc(g,1)
+ l2x(index_l2x_Sl_anidf,i) = lnd2atm_vars%albi_grc(g,2)
+ end if
+
l2x(index_l2x_Sl_tref,i) = lnd2atm_vars%t_ref2m_grc(g)
l2x(index_l2x_Sl_qref,i) = lnd2atm_vars%q_ref2m_grc(g)
l2x(index_l2x_Sl_u10,i) = lnd2atm_vars%u_ref10m_grc(g)
diff --git a/components/elm/src/data_types/GridcellType.F90 b/components/elm/src/data_types/GridcellType.F90
index e33800256e60..3ecabfb0cb89 100644
--- a/components/elm/src/data_types/GridcellType.F90
+++ b/components/elm/src/data_types/GridcellType.F90
@@ -15,6 +15,7 @@ module GridcellType
use landunit_varcon, only : max_lunit
use elm_varcon , only : ispval, spval
use topounit_varcon, only : max_topounits
+ use elm_varpar , only : ndir_horizon_angle
!
! !PUBLIC TYPES:
implicit none
@@ -50,6 +51,12 @@ module GridcellType
real(r8), pointer :: terrain_config (:) => null() ! mean of (terrain configuration factor / cos(slope))
real(r8), pointer :: sinsl_cosas (:) => null() ! sin(slope)*cos(aspect) / cos(slope)
real(r8), pointer :: sinsl_sinas (:) => null() ! sin(slope)*sin(aspect) / cos(slope)
+
+ real(r8), pointer :: slope_deg (:) => null() ! gridcell slope in degree (0-90)
+ real(r8), pointer :: aspect_deg (:) => null() ! gridcell aspect in degree (0-360)
+ real(r8), pointer :: horizon_angle_deg (:,:) => null() ! horizon angles in degree (0-90)
+ real(r8), pointer :: sky_view_factor (:) => null() ! sky view factor (unitless, 0-1)
+ real(r8), pointer :: terrain_config_factor (:) => null() ! terrain configuration factor (unitless, 0-1)
! Daylength
real(r8) , pointer :: max_dayl (:) => null() ! maximum daylength for this grid cell (seconds)
@@ -115,6 +122,12 @@ subroutine grc_pp_init(this, begg, endg)
allocate(this%sinsl_cosas(begg:endg)) ; this%sinsl_cosas(:) = ispval ! sin(slope)*cos(aspect) / cos(slope)
allocate(this%sinsl_sinas(begg:endg)) ; this%sinsl_sinas(:) = ispval ! sin(slope)*sin(aspect) / cos(slope)
+ allocate(this%slope_deg (begg:endg)) ; this%slope_deg (:) = ispval
+ allocate(this%aspect_deg (begg:endg)) ; this%aspect_deg (:) = ispval
+ allocate(this%horizon_angle_deg (begg:endg,1:ndir_horizon_angle)) ; this%horizon_angle_deg(:,:) = ispval
+ allocate(this%sky_view_factor (begg:endg)) ; this%sky_view_factor (:) = ispval
+ allocate(this%terrain_config_factor(begg:endg)) ; this%terrain_config_factor(:) = ispval
+
! This is initiailized in module DayLength
allocate(this%max_dayl (begg:endg)) ; this%max_dayl (:) = spval
allocate(this%dayl (begg:endg)) ; this%dayl (:) = spval
@@ -167,7 +180,12 @@ subroutine grc_pp_clean(this)
deallocate(this%terrain_config )
deallocate(this%sinsl_cosas )
deallocate(this%sinsl_sinas )
-
+ deallocate(this%slope_deg )
+ deallocate(this%aspect_deg )
+ deallocate(this%horizon_angle_deg )
+ deallocate(this%sky_view_factor )
+ deallocate(this%terrain_config_factor)
+
end subroutine grc_pp_clean
end module GridcellType
diff --git a/components/elm/src/data_types/TopounitDataType.F90 b/components/elm/src/data_types/TopounitDataType.F90
index 33de4bb73542..7683b62f376d 100644
--- a/components/elm/src/data_types/TopounitDataType.F90
+++ b/components/elm/src/data_types/TopounitDataType.F90
@@ -63,6 +63,10 @@ module TopounitDataType
real(r8), pointer :: solai (:,:) => null() ! diffuse radiation (numrad) (vis=forc_solsd, nir=forc_solld) (W/m**2)
real(r8), pointer :: solar (:) => null() ! incident solar radiation (W/m**2)
real(r8), pointer :: lwrad (:) => null() ! atm downwrd IR longwave radiation (W/m**2)
+ real(r8), pointer :: solad_pp (:,:) => null() ! direct beam radiation under PP (numrad)(vis=forc_sols , nir=forc_soll) (W/m**2)
+ real(r8), pointer :: solai_pp (:,:) => null() ! diffuse radiation under PP (numrad) (vis=forc_solsd, nir=forc_solld) (W/m**2)
+ real(r8), pointer :: solar_pp (:) => null() ! incident solar radiation under PP (W/m**2)
+ real(r8), pointer :: lwrad_pp (:) => null() ! atm downwrd IR longwave radiation under PP (W/m**2)
! Accumulated fields
real(r8), pointer :: prec24h (:) => null() ! 24-hour mean precip rate (kg H2O/m**2/s, equivalent to mm liquid H2O/s)
real(r8), pointer :: prec10d (:) => null() ! 10-day mean precip rate (kg H2O/m**2/s, equivalent to mm liquid H2O/s)
@@ -367,6 +371,10 @@ subroutine init_top_af(this, begt, endt)
allocate(this%solai (begt:endt, numrad)) ; this%solai (:,:) = spval
allocate(this%solar (begt:endt)) ; this%solar (:) = spval
allocate(this%lwrad (begt:endt)) ; this%lwrad (:) = spval
+ allocate(this%solad_pp (begt:endt, numrad)) ; this%solad_pp (:,:) = spval
+ allocate(this%solai_pp (begt:endt, numrad)) ; this%solai_pp (:,:) = spval
+ allocate(this%solar_pp (begt:endt)) ; this%solar_pp (:) = spval
+ allocate(this%lwrad_pp (begt:endt)) ; this%lwrad_pp (:) = spval
if (use_fates) then
allocate(this%prec24h (begt:endt)) ; this%prec24h (:) = spval
end if
@@ -416,6 +424,10 @@ subroutine clean_top_af(this, begt, endt)
deallocate(this%solai)
deallocate(this%solar)
deallocate(this%lwrad)
+ deallocate(this%solad_pp)
+ deallocate(this%solai_pp)
+ deallocate(this%solar_pp)
+ deallocate(this%lwrad_pp)
if (use_fates) then
deallocate(this%prec24h)
end if
diff --git a/components/elm/src/main/atm2lndMod.F90 b/components/elm/src/main/atm2lndMod.F90
index 27d18c89e867..92b721ee9c70 100644
--- a/components/elm/src/main/atm2lndMod.F90
+++ b/components/elm/src/main/atm2lndMod.F90
@@ -19,7 +19,10 @@ module atm2lndMod
use decompMod , only : bounds_type
use atm2lndType , only : atm2lnd_type
use LandunitType , only : lun_pp
- use ColumnType , only : col_pp
+ use ColumnType , only : col_pp
+ use GridCellType , only : grc_pp
+ use TopounitDataType, only: top_af
+ use lnd2atmType , only : lnd2atm_type
!
! !PUBLIC TYPES:
implicit none
@@ -28,6 +31,7 @@ module atm2lndMod
!
! !PUBLIC MEMBER FUNCTIONS:
public :: downscale_forcings ! Downscale atm forcing fields from gridcell to column
+ public :: topographic_effects_on_radiation ! Topographic effects on shortwave/longwave radiation
!
! !PRIVATE MEMBER FUNCTIONS:
private :: downscale_longwave ! Downscale longwave radiation from gridcell to column
@@ -447,4 +451,169 @@ subroutine check_downscale_consistency(bounds, atm2lnd_vars)
end subroutine check_downscale_consistency
+ !-----------------------------------------------------------------------
+ subroutine topographic_effects_on_radiation(bounds, atm2lnd_vars, nextsw_cday, declin, lnd2atm_vars)
+ !
+ ! !DESCRIPTION:
+ ! Calculate fine-scale topographic effects on shortwave and longwave radiation
+ !
+ ! !USES:
+ use domainMod , only : ldomain
+ use shr_orb_mod
+ use shr_const_mod , only : SHR_CONST_PI
+ use elm_varpar , only : numrad, ndir_horizon_angle
+ !
+ ! !ARGUMENTS:
+ type(bounds_type) , intent(in) :: bounds
+ type(atm2lnd_type) , intent(inout) :: atm2lnd_vars
+ type(lnd2atm_type) , intent(in) :: lnd2atm_vars
+ real(r8) , intent(in) :: nextsw_cday ! calendar day at Greenwich (1.00, ..., days/year)
+ real(r8) , intent(in) :: declin ! declination angle (radians) for next time step
+
+ !
+ ! !LOCAL VARIABLES:
+ integer :: g,topo ! indices
+ integer :: dd,ib
+ real(r8) :: cossza
+ real(r8) :: slope_rad
+ real(r8) :: aspect_rad
+ real(r8) :: deg2rad
+ real(r8) :: horizon_angle_twd_sun_rad
+
+ character(len=*), parameter :: subname = 'topographic_effects_on_radiation'
+ !-----------------------------------------------------------------------
+
+ associate(&
+ ! Gridcell-level fields:
+ lat => grc_pp%lat , & ! Input: latitude
+ lon => grc_pp%lon , & ! Input: longitude
+ slope_deg => grc_pp%slope_deg , &
+ aspect_deg => grc_pp%aspect_deg , &
+ horizon_angle_deg => grc_pp%horizon_angle_deg , &
+ sky_view_factor => grc_pp%sky_view_factor , &
+ terrain_config_factor => grc_pp%terrain_config_factor , &
+ f_short_dir => atm2lnd_vars%f_short_dir , &
+ f_short_dif => atm2lnd_vars%f_short_dif , &
+ f_short_refl => atm2lnd_vars%f_short_refl , &
+ f_long_dif => atm2lnd_vars%f_long_dif , &
+ f_long_refl => atm2lnd_vars%f_long_refl , &
+ sza => atm2lnd_vars%sza , &
+ saa => atm2lnd_vars%saa , &
+ cosinc => atm2lnd_vars%cosinc , &
+ albd => lnd2atm_vars%albd_grc , & ! Input: [real(r8) (:,:) ] surface albedo (direct)
+ albi => lnd2atm_vars%albi_grc , & ! Input: [real(r8) (:,:) ] surface albedo (diffuse)
+ eflx_lwrad_out_grc => lnd2atm_vars%eflx_lwrad_out_grc , &
+ forc_solad_grc => atm2lnd_vars%forc_solad_grc , &
+ forc_solai_grc => atm2lnd_vars%forc_solai_grc , &
+ forc_solar_grc => atm2lnd_vars%forc_solar_grc , &
+ forc_lwrad_grc => atm2lnd_vars%forc_lwrad_not_downscaled_grc , &
+ forc_solad_pp_grc => atm2lnd_vars%forc_solad_pp_grc , &
+ forc_solai_pp_grc => atm2lnd_vars%forc_solai_pp_grc , &
+ forc_solar_pp_grc => atm2lnd_vars%forc_solar_pp_grc , &
+ forc_lwrad_pp_grc => atm2lnd_vars%forc_lwrad_not_downscaled_pp_grc &
+ )
+
+ deg2rad = SHR_CONST_PI/180._r8
+
+ ! Initialize column forcing (needs to be done for ALL active columns)
+ do g = bounds%begg, bounds%endg
+
+ forc_solad_pp_grc(g,:) = forc_solad_grc(g,:)
+ forc_solai_pp_grc(g,:) = forc_solai_grc(g,:)
+ forc_solar_pp_grc(g) = forc_solar_grc(g)
+ forc_lwrad_pp_grc(g) = forc_lwrad_grc(g)
+
+ ! copy radiation values from gridcell to topounit
+ do topo = grc_pp%topi(g), grc_pp%topf(g)
+ top_af%solad_pp(topo,:) = forc_solad_pp_grc(g,:)
+ top_af%solai_pp(topo,:) = forc_solai_pp_grc(g,:)
+ top_af%solar_pp(topo) = forc_solar_pp_grc(g)
+ top_af%lwrad_pp(topo) = forc_lwrad_pp_grc(g)
+ end do
+
+ ! calculate cosine of solar zenith angle
+ cossza = shr_orb_cosz(nextsw_cday, lat(g), lon(g), declin)
+
+ slope_rad = slope_deg(g) * deg2rad
+ aspect_rad = aspect_deg(g) * deg2rad
+
+ f_short_dir(g) = 1._r8
+ f_short_dif(g) = 1._r8
+ f_short_refl(g,:) = 0._r8
+ sza(g) = spval
+ saa(g) = spval
+ cosinc(g) = spval
+ ! scale shortwave radiation
+ if (cossza > 0.0872_r8) then ! just modify when SZA > 85 degree 0.0872_r8
+
+ ! solar zenith angle
+ sza(g) = acos(cossza)
+
+ ! solar azimuth angle
+ saa(g) = shr_orb_azimuth(nextsw_cday, lat(g), lon(g), declin,sza(g))
+
+ ! calculat local solar zenith angle
+ cosinc(g) = cos(slope_rad) * cossza + sin(slope_rad) * sin(sza(g)) * cos(aspect_rad - saa(g))
+ cosinc(g) = max(-1._r8, min(cosinc(g), 1._r8))
+
+ if (cosinc(g) < 0._r8) then
+ f_short_dir(g) = 0._r8
+ else
+ f_short_dir(g) = cosinc(g) / cossza / cos(slope_rad)
+ endif
+
+ ! Find horizion angle towards the sun
+ dd = nint(saa(g) / (2._r8*SHR_CONST_PI) * ndir_horizon_angle) + 1
+ if (dd > ndir_horizon_angle) dd = 1
+ horizon_angle_twd_sun_rad = horizon_angle_deg(g,dd) * deg2rad
+
+ ! Check if sun is above horizon angle
+ if (cossza < sin(horizon_angle_twd_sun_rad)) then
+ f_short_dir(g) = 0._r8
+ endif
+
+ f_short_dif(g) = sky_view_factor(g) / cos(slope_rad)
+ if (f_short_dif(g) < 0._r8) f_short_dif(g) = 0._r8
+
+ forc_solar_grc(g) = 0._r8
+ do ib = 1, numrad
+ ! Calculate reflected radiation from adjacent terrain
+ f_short_refl(g,ib) = terrain_config_factor(g) / cos(slope_rad) * (albd(g,ib) * forc_solad_grc(g,ib) + albi(g,ib) * forc_solai_grc(g,ib))
+
+ if (f_short_refl(g,ib) < 0._r8) f_short_refl(g,ib) = 0._r8
+
+ ! scale direct solar radiation: vis & nir
+ forc_solad_grc(g,ib) = forc_solad_grc(g,ib) * f_short_dir(g)
+ ! scale diffuse solar radiation: vis & nir
+ forc_solai_grc(g,ib) = forc_solai_grc(g,ib) * f_short_dif(g) + f_short_refl(g,ib)
+
+ forc_solar_grc(g) = forc_solar_grc(g) + forc_solad_grc(g,ib) + forc_solai_grc(g,ib)
+ end do
+
+ end if
+
+ f_long_dif(g) = 1._r8
+ f_long_refl(g) = 0._r8
+ ! scale longwave radiation
+ f_long_dif(g) = sky_view_factor(g) / cos(slope_rad)
+ if (f_long_dif(g) < 0._r8) f_long_dif(g) = 0._r8
+ f_long_refl(g) = terrain_config_factor(g) * eflx_lwrad_out_grc(g)
+ if (f_long_refl(g) < 0._r8) f_long_refl(g) = 0._r8
+
+ forc_lwrad_grc(g) = forc_lwrad_grc(g) * f_long_dif(g) + f_long_refl(g)
+
+ ! copy radiation values from gridcell to topounit
+ do topo = grc_pp%topi(g), grc_pp%topf(g)
+ top_af%solad(topo,:) = forc_solad_grc(g,:)
+ top_af%solai(topo,:) = forc_solai_grc(g,:)
+ top_af%solar(topo) = forc_solar_grc(g)
+ top_af%lwrad(topo) = forc_lwrad_grc(g)
+ end do
+
+ end do
+
+ end associate
+
+ end subroutine topographic_effects_on_radiation
+
end module atm2lndMod
diff --git a/components/elm/src/main/atm2lndType.F90 b/components/elm/src/main/atm2lndType.F90
index f2facc776224..d663fe085d75 100644
--- a/components/elm/src/main/atm2lndType.F90
+++ b/components/elm/src/main/atm2lndType.F90
@@ -11,7 +11,7 @@ module atm2lndType
use shr_megan_mod , only : shr_megan_mechcomps_n
use elm_varpar , only : numrad, ndst, nlevgrnd !ndst = number of dust bins.
use elm_varcon , only : rair, grav, cpair, hfus, tfrz, spval
- use elm_varctl , only : iulog, use_c13, use_cn, use_lch4, use_fates, use_fan
+ use elm_varctl , only : iulog, use_c13, use_cn, use_lch4, use_fates, use_fan, use_finetop_rad
use seq_drydep_mod, only : n_drydep, drydep_method, DD_XLND
use decompMod , only : bounds_type
use abortutils , only : endrun
@@ -148,7 +148,20 @@ module atm2lndType
! Needed for FNM precip downscaling, when used within CPL_BYPASS
real(r8), pointer :: forc_uovern (:) => null() ! Froude number (dimensionless)
-
+
+ real(r8), pointer :: sza (:) => null() ! solar zenith angle
+ real(r8), pointer :: saa (:) => null() ! solar azmith angle
+ real(r8), pointer :: cosinc (:) => null() ! cosine of the local incident angle
+ real(r8), pointer :: f_short_dir (:) => null() ! adjust factor for direct shortwave radiation
+ real(r8), pointer :: f_short_dif (:) => null() ! adjust factor for diffuse shortwave radiation
+ real(r8), pointer :: f_short_refl (:,:) => null() ! adjust factor for reflected shortwave radiation
+ real(r8), pointer :: f_long_dif (:) => null() ! adjust factor for diffuse longwave radiation
+ real(r8), pointer :: f_long_refl (:) => null() ! adjust factor for reflected longwave radiation
+
+ real(r8), pointer :: forc_solad_pp_grc (:,:) => null() ! direct beam radiation (numrad) (vis=forc_sols , nir=forc_soll ) under PP
+ real(r8), pointer :: forc_solai_pp_grc (:,:) => null() ! diffuse radiation (numrad) (vis=forc_solsd, nir=forc_solld) under PP
+ real(r8), pointer :: forc_solar_pp_grc (:) => null() ! incident solar radiation under PP
+ real(r8), pointer :: forc_lwrad_not_downscaled_pp_grc(:) => null() ! not downscaled atm downwrd IR longwave radiation (W/m**2) under PP
contains
@@ -314,7 +327,23 @@ subroutine InitAllocate(this, bounds)
allocate(this%forc_ndep_nitr_grc (begg:endg)) ; this%forc_ndep_nitr_grc (:) = ival
allocate(this%forc_soilph_grc (begg:endg)) ; this%forc_soilph_grc (:) = ival
end if
+
allocate(this%forc_uovern (begg:endg)) ; this%forc_uovern (:) = ival
+
+ if ( use_finetop_rad ) then
+ allocate(this%f_short_dir (begg:endg)) ; this%f_short_dir (:) = spval
+ allocate(this%f_short_dif (begg:endg)) ; this%f_short_dif (:) = spval
+ allocate(this%f_short_refl (begg:endg,numrad)) ; this%f_short_refl (:,:) = spval
+ allocate(this%f_long_dif (begg:endg)) ; this%f_long_dif (:) = spval
+ allocate(this%f_long_refl (begg:endg)) ; this%f_long_refl (:) = spval
+ allocate(this%sza (begg:endg)) ; this%sza (:) = spval
+ allocate(this%saa (begg:endg)) ; this%saa (:) = spval
+ allocate(this%cosinc (begg:endg)) ; this%cosinc (:) = spval
+ allocate(this%forc_solad_pp_grc (begg:endg,numrad)) ; this%forc_solad_pp_grc (:,:) = ival
+ allocate(this%forc_solai_pp_grc (begg:endg,numrad)) ; this%forc_solai_pp_grc (:,:) = ival
+ allocate(this%forc_solar_pp_grc (begg:endg)) ; this%forc_solar_pp_grc (:) = ival
+ allocate(this%forc_lwrad_not_downscaled_pp_grc (begg:endg)) ; this%forc_lwrad_not_downscaled_pp_grc (:)= ival
+ endif
end subroutine InitAllocate
@@ -322,7 +351,7 @@ end subroutine InitAllocate
subroutine InitHistory(this, bounds)
!
! !USES:
- use histFileMod, only : hist_addfld1d
+ use histFileMod, only : hist_addfld1d, hist_addfld2d
!
! !ARGUMENTS:
class(atm2lnd_type) :: this
@@ -504,6 +533,56 @@ subroutine InitHistory(this, bounds)
avgflag='A', long_name='direct radiation (last 240hrs)', &
ptr_patch=this%fsd240_patch, default='inactive')
+ if ( use_finetop_rad ) then
+ this%f_short_dir(begg:endg) = spval
+ call hist_addfld1d (fname='F_SHORT_DIR', units='unitless', &
+ avgflag='A', long_name='f_short_dir', &
+ ptr_gcell=this%f_short_dir, default='inactive')
+
+ this%f_short_dif(begg:endg) = spval
+ call hist_addfld1d (fname='F_SHORT_DIF', units='unitless', &
+ avgflag='A', long_name='f_short_dif', &
+ ptr_gcell=this%f_short_dif, default='inactive')
+
+ this%f_short_refl(begg:endg,:) = spval
+ call hist_addfld2d (fname='F_SHORT_REFL', units='unitless', type2d='numrad', &
+ avgflag='A', long_name='f_short_refl', &
+ ptr_gcell=this%f_short_refl, default='inactive')
+
+ this%f_long_dif(begg:endg) = spval
+ call hist_addfld1d (fname='F_LONG_DIF', units='unitless', &
+ avgflag='A', long_name='f_long_dif', &
+ ptr_gcell=this%f_long_dif, default='inactive')
+
+ this%f_long_refl(begg:endg) = spval
+ call hist_addfld1d (fname='F_LONG_REFL', units='unitless', &
+ avgflag='A', long_name='f_long_refl', &
+ ptr_gcell=this%f_long_refl, default='inactive')
+
+ this%forc_solad_pp_grc(begg:endg,:) = spval
+ call hist_addfld2d (fname='forc_solad_pp_grc', units='unitless', type2d='numrad', &
+ avgflag='A', long_name='forc_solad_pp_grc', &
+ ptr_gcell=this%forc_solad_pp_grc, default='inactive')
+
+ this%forc_solai_pp_grc(begg:endg,:) = spval
+ call hist_addfld2d (fname='forc_solai_pp_grc', units='unitless', type2d='numrad', &
+ avgflag='A', long_name='forc_solai_pp_grc', &
+ ptr_gcell=this%forc_solai_pp_grc, default='inactive')
+
+ this%forc_solar_pp_grc(begg:endg) = spval
+ call hist_addfld1d (fname='SWdown_PP', units='W/m^2', &
+ avgflag='A', long_name='atmospheric incident solar radiation (PP)', &
+ ptr_gcell=this%forc_solar_pp_grc, default='inactive')
+
+ this%forc_lwrad_not_downscaled_pp_grc(begg:endg) = spval
+ call hist_addfld1d (fname='LWdown_PP', units='W/m^2', &
+ avgflag='A', long_name='atmospheric longwave radiation (PP)', &
+ ptr_gcell=this%forc_lwrad_not_downscaled_pp_grc, default='inactive')
+
+ this%forc_solad_pp_grc(begg:endg, :) = 0._r8
+ this%forc_solai_pp_grc(begg:endg, :) = 0._r8
+ endif
+
end subroutine InitHistory
!-----------------------------------------------------------------------
@@ -760,6 +839,18 @@ subroutine Restart(this, bounds, ncid, flag)
this%forc_flood_grc = 0._r8
endif
+ if ( use_finetop_rad ) then
+ call restartvar(ncid=ncid, flag=flag, varname='forc_solad_pp', xtype=ncd_double, &
+ dim1name='gridcell', dim2name='numrad', switchdim=.true., &
+ long_name='direct beam radiation (numrad)', units='W/m^2', &
+ interpinic_flag='interp', readvar=readvar, data=this%forc_solad_pp_grc)
+
+ call restartvar(ncid=ncid, flag=flag, varname='forc_solai_pp', xtype=ncd_double, &
+ dim1name='gridcell', dim2name='numrad', switchdim=.true., &
+ long_name='diffuse radiation (numrad)', units='W/m^2', &
+ interpinic_flag='interp', readvar=readvar, data=this%forc_solai_pp_grc)
+ endif
+
end subroutine Restart
end module atm2lndType
diff --git a/components/elm/src/main/controlMod.F90 b/components/elm/src/main/controlMod.F90
index edb50f316959..ba23a9ab57b8 100755
--- a/components/elm/src/main/controlMod.F90
+++ b/components/elm/src/main/controlMod.F90
@@ -53,7 +53,7 @@ module controlMod
use elm_varctl , only: startdate_add_temperature, startdate_add_co2
use elm_varctl , only: add_temperature, add_co2
use elm_varctl , only: const_climate_hist
- use elm_varctl , only: use_top_solar_rad
+ use elm_varctl , only: use_top_solar_rad, use_finetop_rad
use elm_varctl , only: snow_shape, snicar_atm_type, use_dust_snow_internal_mixing
use EcosystemBalanceCheckMod, only: bgc_balance_check_tolerance => balance_check_tolerance
@@ -342,7 +342,7 @@ subroutine control_init( )
use_erosion, ero_ccycle
namelist /elm_inparm/ &
- use_top_solar_rad
+ use_top_solar_rad, use_finetop_rad
namelist /elm_mosart/ &
lnd_rof_coupling_nstep
@@ -925,6 +925,7 @@ subroutine control_spmd()
call mpi_bcast (more_vertlayers,1, MPI_LOGICAL, 0, mpicom, ier)
call mpi_bcast (const_climate_hist, 1, MPI_LOGICAL, 0, mpicom, ier)
call mpi_bcast (use_top_solar_rad, 1, MPI_LOGICAL, 0, mpicom, ier) ! TOP solar radiation parameterization
+ call mpi_bcast (use_finetop_rad, 1, MPI_LOGICAL, 0, mpicom, ier) ! fineTOP radiation parameterization
! glacier_mec variables
call mpi_bcast (create_glacier_mec_landunit, 1, MPI_LOGICAL, 0, mpicom, ier)
@@ -1121,9 +1122,19 @@ subroutine control_print ()
if (use_top_solar_rad) then
write(iulog,*) ' use TOP solar radiation parameterization instead of PP'
else
- write(iulog,*) ' use_top_solar_rad is False, so do not run TOP solar radiation parameterization'
+ write(iulog,*) ' use_top_solar_rad is False, so do not run TOP solar radiation parameterization'
end if
-
+
+ if (use_finetop_rad .and. use_top_solar_rad) then
+ write(iulog,*) ' cannot use both TOP and fineTOP radiation parameterizations simultaneously'
+ call endrun(msg=' ERROR: use_finetop_rad and use_top_solar_rad cannot both be set to true.'//&
+ errMsg(__FILE__, __LINE__))
+ else if (use_finetop_rad .and. (.not. use_top_solar_rad)) then
+ write(iulog,*) ' use fineTOP radiation parameterization instead of PP'
+ else
+ write(iulog,*) ' use_finetop_rad is False, so do not run fineTOP radiation parameterization'
+ end if
+
if (use_cn) then
if (suplnitro /= suplnNon)then
write(iulog,*) ' Supplemental Nitrogen mode is set to run over Patches: ', &
@@ -1242,7 +1253,8 @@ subroutine control_print ()
write(iulog,*) ' more vertical layers = ', more_vertlayers
write(iulog,*) ' Sub-grid topographic effects on solar radiation = ', use_top_solar_rad ! TOP solar radiation parameterization
-
+ write(iulog,*) ' Grid-scale topographic effects on radiation (fineTOP) = ', use_finetop_rad ! fineTOP radiation parameterization
+
if (nsrest == nsrContinue) then
write(iulog,*) 'restart warning:'
write(iulog,*) ' Namelist not checked for agreement with initial run.'
diff --git a/components/elm/src/main/elm_driver.F90 b/components/elm/src/main/elm_driver.F90
index 896095a12c92..f1257b7b5e27 100644
--- a/components/elm/src/main/elm_driver.F90
+++ b/components/elm/src/main/elm_driver.F90
@@ -83,7 +83,7 @@ module elm_driver
!
use filterMod , only : setFilters
!
- use atm2lndMod , only : downscale_forcings
+ use atm2lndMod , only : downscale_forcings, topographic_effects_on_radiation
use lnd2atmMod , only : lnd2atm
use lnd2glcMod , only : lnd2glc_type
use lnd2iacMod , only : lnd2iac_type
@@ -181,6 +181,7 @@ module elm_driver
use CNPBudgetMod , only : CNPBudget_SetBeginningMonthlyStates, CNPBudget_SetEndingMonthlyStates
use elm_varctl , only : do_budgets, budget_inst, budget_daily, budget_month
use elm_varctl , only : budget_ann, budget_ltann, budget_ltend
+ use elm_varctl , only : use_finetop_rad
use timeinfoMod
!
@@ -687,6 +688,12 @@ subroutine elm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate)
filter(nc)%num_soilp , filter(nc)%soilp, &
canopystate_vars, energyflux_vars)
+ if (use_finetop_rad) then
+ call topographic_effects_on_radiation(bounds_clump, &
+ atm2lnd_vars, nextsw_cday, declinp1, &
+ lnd2atm_vars)
+ endif
+
call downscale_forcings(bounds_clump, &
filter(nc)%num_do_smb_c, filter(nc)%do_smb_c, &
atm2lnd_vars)
diff --git a/components/elm/src/main/elm_initializeMod.F90 b/components/elm/src/main/elm_initializeMod.F90
index 43e124372c83..b061e481884d 100755
--- a/components/elm/src/main/elm_initializeMod.F90
+++ b/components/elm/src/main/elm_initializeMod.F90
@@ -12,7 +12,7 @@ module elm_initializeMod
use elm_varctl , only : nsrest, nsrStartup, nsrContinue, nsrBranch
use elm_varctl , only : create_glacier_mec_landunit, iulog, iac_present
use elm_varctl , only : use_lch4, use_cn, use_voc, use_c13, use_c14
- use elm_varctl , only : use_fates, use_betr, use_fates_sp, use_fan, use_fates_luh
+ use elm_varctl , only : use_fates, use_betr, use_fates_sp, use_fan, use_fates_luh, use_finetop_rad
use elm_varsur , only : wt_lunit, urban_valid, wt_nat_patch, wt_cft, wt_glc_mec, topo_glc_mec,firrig,f_surf,f_grd
use elm_varsur , only : fert_cft, fert_p_cft, wt_polygon
use elm_varsur , only : wt_tunit, elv_tunit, slp_tunit,asp_tunit,num_tunit_per_grd
@@ -82,7 +82,7 @@ subroutine initialize1( )
use decompInitMod , only: decompInit_moab
#endif
use domainMod , only: domain_check, ldomain, domain_init
- use surfrdMod , only: surfrd_get_globmask, surfrd_get_grid, surfrd_get_topo, surfrd_get_data,surfrd_get_topo_for_solar_rad
+ use surfrdMod , only: surfrd_get_globmask, surfrd_get_grid, surfrd_get_topo, surfrd_get_data, surfrd_get_topo_for_solar_rad, surfrd_finetop_data
use controlMod , only: control_init, control_print, NLFilename
use ncdio_pio , only: ncd_pio_init
use initGridCellsMod , only: initGridCells, initGhostGridCells
@@ -276,9 +276,10 @@ subroutine initialize1( )
write(iulog,*) 'Attempting to read topo parameters for TOP solar radiation parameterization from ',trim(fsurdat)
call shr_sys_flush(iulog)
endif
- call surfrd_get_topo_for_solar_rad(ldomain, fsurdat)
- endif
+ call surfrd_get_topo_for_solar_rad(ldomain, fsurdat)
+ endif
+
!-------------------------------------------------------------------------
! Topounit
!-------------------------------------------------------------------------
@@ -422,6 +423,14 @@ subroutine initialize1( )
call initGridCells()
+ if (fsurdat /= " " .and. use_finetop_rad) then
+ if (masterproc) then
+ write(iulog,*) 'Attempting to read topo parameters for fineTOP parameterization from ',trim(fsurdat)
+ call shr_sys_flush(iulog)
+ endif
+ call surfrd_finetop_data(ldomain, fsurdat)
+ endif
+
! Set global seg maps for gridcells, topounits, landlunits, columns and patches
!if(max_topounits > 1) then
! if (create_glacier_mec_landunit) then
@@ -1035,10 +1044,16 @@ subroutine initialize2( )
if (nsrest == nsrStartup) then
call t_startf('init_map2gc')
- call lnd2atm_minimal(bounds_proc, surfalb_vars, energyflux_vars, lnd2atm_vars)
+ call lnd2atm_minimal(bounds_proc, surfalb_vars, solarabs_vars, energyflux_vars, atm2lnd_vars, lnd2atm_vars)
+ call t_stopf('init_map2gc')
+ else if ( use_finetop_rad .and. ((nsrest == nsrContinue) .or. (nsrest == nsrBranch))) then
+ call t_startf('init_map2gc')
+ call lnd2atm_minimal(bounds_proc, surfalb_vars, solarabs_vars, energyflux_vars, atm2lnd_vars, lnd2atm_vars)
call t_stopf('init_map2gc')
end if
+
+
!------------------------------------------------------------
! Initialize sno export state to send to glc
!------------------------------------------------------------
diff --git a/components/elm/src/main/elm_varctl.F90 b/components/elm/src/main/elm_varctl.F90
index 66fb18bc23b2..c9c44e06154b 100644
--- a/components/elm/src/main/elm_varctl.F90
+++ b/components/elm/src/main/elm_varctl.F90
@@ -401,7 +401,8 @@ module elm_varctl
character(len = SHR_KIND_CS), public :: precip_downscaling_method = 'ERMM' ! Precip downscaling method values can be ERMM or FNM
logical, public :: use_lake_wat_storage = .false.
logical, public :: use_top_solar_rad = .false. ! TOP : sub-grid topographic effect on surface solar radiation
-
+ logical, public :: use_finetop_rad = .false. ! fineTOP : fine(grid)-scale topographic effect on surface radiation balance (longwave + shortwave)
+
!----------------------------------------------------------
! Fan controls (use_fan)
!----------------------------------------------------------
diff --git a/components/elm/src/main/elm_varpar.F90 b/components/elm/src/main/elm_varpar.F90
index 36be64825abc..fa9c38005471 100644
--- a/components/elm/src/main/elm_varpar.F90
+++ b/components/elm/src/main/elm_varpar.F90
@@ -83,6 +83,8 @@ module elm_varpar
real(r8), parameter :: scalez = 0.025_r8
real(r8), parameter :: zecoeff = 0.50_r8
+ integer, parameter :: ndir_horizon_angle = 8 ! number of directions for horizon angle
+
! constants for decomposition cascade
integer :: i_met_lit
diff --git a/components/elm/src/main/lnd2atmMod.F90 b/components/elm/src/main/lnd2atmMod.F90
index 978be5a1b6ff..49d052d1e4d3 100644
--- a/components/elm/src/main/lnd2atmMod.F90
+++ b/components/elm/src/main/lnd2atmMod.F90
@@ -13,7 +13,7 @@ module lnd2atmMod
use elm_varpar , only : numrad, ndst, nlevgrnd, nlevsno, nlevsoi !ndst = number of dust bins.
use elm_varcon , only : rair, grav, cpair, hfus, tfrz, spval
use elm_varctl , only : iulog, use_c13, use_cn, use_lch4, use_voc, use_fates, use_atm_downscaling_to_topunit, use_fan
- use elm_varctl , only : use_lnd_rof_two_way
+ use elm_varctl , only : use_lnd_rof_two_way, use_finetop_rad
use tracer_varcon , only : is_active_betr_bgc
use seq_drydep_mod , only : n_drydep, drydep_method, DD_XLND
use decompMod , only : bounds_type
@@ -55,7 +55,7 @@ module lnd2atmMod
!------------------------------------------------------------------------
subroutine lnd2atm_minimal(bounds, &
- surfalb_vars, energyflux_vars, lnd2atm_vars)
+ surfalb_vars, solarabs_vars, energyflux_vars, atm2lnd_vars, lnd2atm_vars)
!
! !DESCRIPTION:
! Compute elm_l2a_vars component of gridcell derived type. This routine computes
@@ -70,11 +70,13 @@ subroutine lnd2atm_minimal(bounds, &
! !ARGUMENTS:
type(bounds_type) , intent(in) :: bounds
type(surfalb_type) , intent(in) :: surfalb_vars
+ type(solarabs_type) , intent(in) :: solarabs_vars
type(energyflux_type) , intent(in) :: energyflux_vars
+ type(atm2lnd_type) , intent(in) :: atm2lnd_vars
type(lnd2atm_type) , intent(inout) :: lnd2atm_vars
!
! !LOCAL VARIABLES:
- integer :: g, t ! index
+ integer :: g, t, ib ! index
!------------------------------------------------------------------------
associate( &
h2osno => col_ws%h2osno , &
@@ -83,10 +85,23 @@ subroutine lnd2atm_minimal(bounds, &
h2osoi_vol_grc => lnd2atm_vars%h2osoi_vol_grc , &
albd_patch => surfalb_vars%albd_patch , &
albd_grc => lnd2atm_vars%albd_grc , &
+ apparent_albd_grc => lnd2atm_vars%apparent_albd_grc , &
albi_patch => surfalb_vars%albi_patch , &
albi_grc => lnd2atm_vars%albi_grc , &
+ apparent_albi_grc => lnd2atm_vars%apparent_albi_grc , &
eflx_lwrad_out => veg_ef%eflx_lwrad_out , &
- eflx_lwrad_out_grc => lnd2atm_vars%eflx_lwrad_out_grc &
+ eflx_lwrad_out_grc => lnd2atm_vars%eflx_lwrad_out_grc , &
+ forc_solad_pp_grc => atm2lnd_vars%forc_solad_pp_grc , &
+ forc_solai_pp_grc => atm2lnd_vars%forc_solai_pp_grc , &
+ fsr_vis_d_patch => solarabs_vars%fsr_vis_d_patch , &
+ fsr_vis_i_patch => solarabs_vars%fsr_vis_i_patch , &
+ fsr_nir_d_patch => solarabs_vars%fsr_nir_d_patch , &
+ fsr_nir_i_patch => solarabs_vars%fsr_nir_i_patch , &
+ fsr_vis_d_grc => lnd2atm_vars%fsr_vis_d_grc , &
+ fsr_vis_i_grc => lnd2atm_vars%fsr_vis_i_grc , &
+ fsr_nir_d_grc => lnd2atm_vars%fsr_nir_d_grc , &
+ fsr_nir_i_grc => lnd2atm_vars%fsr_nir_i_grc , &
+ sky_view_factor => grc_pp%sky_view_factor &
)
call c2g(bounds, &
@@ -118,10 +133,16 @@ subroutine lnd2atm_minimal(bounds, &
eflx_lwrad_out_grc (bounds%begg:bounds%endg), &
p2c_scale_type=unity, c2l_scale_type= urbanf, l2g_scale_type=unity)
+ if (use_finetop_rad) then
+ do g = bounds%begg,bounds%endg
+ eflx_lwrad_out_grc(g) = eflx_lwrad_out_grc(g)*sky_view_factor(g)
+ end do
+ end if
+
do g = bounds%begg,bounds%endg
lnd2atm_vars%t_rad_grc(g) = sqrt(sqrt(eflx_lwrad_out_grc(g)/sb))
end do
-
+
! Calculate topounit level eflx_lwrad_out_topo for downscaling purpose
if (use_atm_downscaling_to_topunit) then
call p2t(bounds, &
@@ -134,6 +155,68 @@ subroutine lnd2atm_minimal(bounds, &
end do
end if
+ if (use_finetop_rad) then
+ ! calculate gridcell level reflected radiation
+ call p2g(bounds, &
+ fsr_vis_d_patch(bounds%begp:bounds%endp) , &
+ fsr_vis_d_grc (bounds%begg:bounds%endg) , &
+ p2c_scale_type=unity, c2l_scale_type= urbanf, l2g_scale_type=unity)
+
+ call p2g(bounds, &
+ fsr_vis_i_patch(bounds%begp:bounds%endp) , &
+ fsr_vis_i_grc (bounds%begg:bounds%endg) , &
+ p2c_scale_type=unity, c2l_scale_type= urbanf, l2g_scale_type=unity)
+
+ call p2g(bounds, &
+ fsr_nir_d_patch(bounds%begp:bounds%endp) , &
+ fsr_nir_d_grc (bounds%begg:bounds%endg) , &
+ p2c_scale_type=unity, c2l_scale_type= urbanf, l2g_scale_type=unity)
+
+ call p2g(bounds, &
+ fsr_nir_i_patch(bounds%begp:bounds%endp) , &
+ fsr_nir_i_grc (bounds%begg:bounds%endg) , &
+ p2c_scale_type=unity, c2l_scale_type= urbanf, l2g_scale_type=unity)
+
+ ! calculate gridcell level apparent albedo
+ do g = bounds%begg,bounds%endg
+ do ib = 1, numrad
+ apparent_albd_grc(g,ib) = 1._r8
+ apparent_albi_grc(g,ib) = 1._r8
+ end do
+
+ ! calculate direct albedo
+ ib = 1
+ if (forc_solad_pp_grc(g,ib) > 0._r8) then
+ apparent_albd_grc(g,ib) = fsr_vis_d_grc(g)/forc_solad_pp_grc(g,ib)*sky_view_factor(g)
+ end if
+ ib = 2
+ if (forc_solad_pp_grc(g,ib) > 0._r8) then
+ apparent_albd_grc(g,ib) = fsr_nir_d_grc(g)/forc_solad_pp_grc(g,ib)*sky_view_factor(g)
+ end if
+
+ ! calculate diffuse albedo
+ ib = 1
+ if (forc_solai_pp_grc(g,ib) > 0._r8) then
+ apparent_albi_grc(g,ib) = fsr_vis_i_grc(g)/forc_solai_pp_grc(g,ib)*sky_view_factor(g)
+ end if
+ ib = 2
+ if (forc_solai_pp_grc(g,ib) > 0._r8) then
+ apparent_albi_grc(g,ib) = fsr_nir_i_grc(g)/forc_solai_pp_grc(g,ib)*sky_view_factor(g)
+ end if
+
+ ! limit albedo to be <= 1
+ do ib = 1, numrad
+ if (apparent_albd_grc(g,ib) > 1._r8) then
+ apparent_albd_grc(g,ib) = 1._r8
+ end if
+
+ if (apparent_albi_grc(g,ib) > 1._r8) then
+ apparent_albi_grc(g,ib) = 1._r8
+ end if
+ end do
+ end do
+ end if
+
end associate
end subroutine lnd2atm_minimal
@@ -246,7 +329,7 @@ subroutine lnd2atm(bounds, &
! First, compute the "minimal" set of fields.
call lnd2atm_minimal(bounds, &
- surfalb_vars, energyflux_vars, lnd2atm_vars)
+ surfalb_vars, solarabs_vars, energyflux_vars, atm2lnd_vars, lnd2atm_vars)
call p2g(bounds, &
t_ref2m (bounds%begp:bounds%endp), &
diff --git a/components/elm/src/main/lnd2atmType.F90 b/components/elm/src/main/lnd2atmType.F90
index 64bffe47925a..4396256a0a88 100644
--- a/components/elm/src/main/lnd2atmType.F90
+++ b/components/elm/src/main/lnd2atmType.F90
@@ -13,7 +13,7 @@ module lnd2atmType
use abortutils , only : endrun
use elm_varpar , only : numrad, ndst, nlevsno, nlevgrnd !ndst = number of dust bins.
use elm_varcon , only : rair, grav, cpair, hfus, tfrz, spval
- use elm_varctl , only : iulog, use_c13, use_cn, use_lch4, use_fan
+ use elm_varctl , only : iulog, use_c13, use_cn, use_lch4, use_fan, use_finetop_rad
use seq_drydep_mod, only : n_drydep, drydep_method, DD_XLND
use decompMod , only : bounds_type
!
@@ -37,6 +37,8 @@ module lnd2atmType
real(r8), pointer :: h2osoi_vol_grc (:,:) => null() ! volumetric soil water (0~watsat, m3/m3, nlevgrnd) (for dust model)
real(r8), pointer :: albd_grc (:,:) => null() ! (numrad) surface albedo (direct)
real(r8), pointer :: albi_grc (:,:) => null() ! (numrad) surface albedo (diffuse)
+ real(r8), pointer :: apparent_albd_grc (:,:) => null() ! (numrad) apparent surface albedo (direct)
+ real(r8), pointer :: apparent_albi_grc (:,:) => null() ! (numrad) apparent surface albedo (diffuse)
real(r8), pointer :: taux_grc (:) => null() ! wind stress: e-w (kg/m/s**2)
real(r8), pointer :: tauy_grc (:) => null() ! wind stress: n-s (kg/m/s**2)
real(r8), pointer :: eflx_lh_tot_grc (:) => null() ! total latent HF (W/m**2) [+ to atm]
@@ -44,6 +46,10 @@ module lnd2atmType
real(r8), pointer :: eflx_lwrad_out_grc (:) => null() ! IR (longwave) radiation (W/m**2)
real(r8), pointer :: qflx_evap_tot_grc (:) => null() ! qflx_evap_soi + qflx_evap_can + qflx_tran_veg
real(r8), pointer :: fsa_grc (:) => null() ! solar rad absorbed (total) (W/m**2)
+ real(r8), pointer :: fsr_vis_d_grc (:) => null() ! reflected direct beam vis solar radiation (W/m**2)
+ real(r8), pointer :: fsr_vis_i_grc (:) => null() ! reflected diffuse vis solar radiation (W/m**2)
+ real(r8), pointer :: fsr_nir_d_grc (:) => null() ! reflected direct beam nir solar radiation (W/m**2)
+ real(r8), pointer :: fsr_nir_i_grc (:) => null() ! reflected diffuse nir solar radiation (W/m**2)
real(r8), pointer :: nee_grc (:) => null() ! net CO2 flux (kg CO2/m**2/s) [+ to atm]
real(r8), pointer :: nem_grc (:) => null() ! gridcell average net methane correction to CO2 flux (g C/m^2/s)
real(r8), pointer :: ram1_grc (:) => null() ! aerodynamical resistance (s/m)
@@ -127,6 +133,8 @@ subroutine InitAllocate(this, bounds)
allocate(this%h2osoi_vol_grc (begg:endg,1:nlevgrnd)) ; this%h2osoi_vol_grc (:,:) =ival
allocate(this%albd_grc (begg:endg,1:numrad)) ; this%albd_grc (:,:) =ival
allocate(this%albi_grc (begg:endg,1:numrad)) ; this%albi_grc (:,:) =ival
+ allocate(this%apparent_albd_grc (begg:endg,1:numrad)) ; this%apparent_albd_grc (:,:) =ival
+ allocate(this%apparent_albi_grc (begg:endg,1:numrad)) ; this%apparent_albi_grc (:,:) =ival
allocate(this%taux_grc (begg:endg)) ; this%taux_grc (:) =ival
allocate(this%tauy_grc (begg:endg)) ; this%tauy_grc (:) =ival
allocate(this%eflx_lwrad_out_grc (begg:endg)) ; this%eflx_lwrad_out_grc (:) =ival
@@ -134,6 +142,10 @@ subroutine InitAllocate(this, bounds)
allocate(this%eflx_lh_tot_grc (begg:endg)) ; this%eflx_lh_tot_grc (:) =ival
allocate(this%qflx_evap_tot_grc (begg:endg)) ; this%qflx_evap_tot_grc (:) =ival
allocate(this%fsa_grc (begg:endg)) ; this%fsa_grc (:) =ival
+ allocate(this%fsr_vis_d_grc (begg:endg)) ; this%fsr_vis_d_grc (:) =ival
+ allocate(this%fsr_vis_i_grc (begg:endg)) ; this%fsr_vis_i_grc (:) =ival
+ allocate(this%fsr_nir_d_grc (begg:endg)) ; this%fsr_nir_d_grc (:) =ival
+ allocate(this%fsr_nir_i_grc (begg:endg)) ; this%fsr_nir_i_grc (:) =ival
allocate(this%nee_grc (begg:endg)) ; this%nee_grc (:) =ival
allocate(this%nem_grc (begg:endg)) ; this%nem_grc (:) =ival
allocate(this%ram1_grc (begg:endg)) ; this%ram1_grc (:) =ival
@@ -185,7 +197,7 @@ end subroutine InitAllocate
subroutine InitHistory(this, bounds)
!
! !USES:
- use histFileMod, only : hist_addfld1d
+ use histFileMod, only : hist_addfld1d, hist_addfld2d
!
! !ARGUMENTS:
class(lnd2atm_type) :: this
@@ -231,6 +243,18 @@ subroutine InitHistory(this, bounds)
ptr_lnd=this%flux_nh3_grc)
end if
+ if (use_finetop_rad) then
+ this%apparent_albd_grc(begg:endg,:) = spval
+ call hist_addfld2d (fname='APPARENT_ALBD', units='proportion', type2d='numrad', &
+ avgflag='A', long_name='apparent surface albedo (direct)', &
+ ptr_gcell=this%apparent_albd_grc, default='inactive', c2l_scale_type='urbanf')
+
+ this%apparent_albi_grc(begg:endg,:) = spval
+ call hist_addfld2d (fname='APPARENT_ALBI', units='proportion', type2d='numrad', &
+ avgflag='A', long_name='apparent surface albedo (indirect)', &
+ ptr_gcell=this%apparent_albi_grc, default='inactive', c2l_scale_type='urbanf')
+ end if
+
end subroutine InitHistory
end module lnd2atmType
diff --git a/components/elm/src/main/surfrdMod.F90 b/components/elm/src/main/surfrdMod.F90
index 48cd79c91f36..a48b2003fbda 100755
--- a/components/elm/src/main/surfrdMod.F90
+++ b/components/elm/src/main/surfrdMod.F90
@@ -36,6 +36,7 @@ module surfrdMod
public :: surfrd_get_grid_conn ! Reads grid connectivity information from domain file
public :: surfrd_topounit_data ! Read topounit physical properties
public :: surfrd_get_topo_for_solar_rad ! Read topography dataset for TOP solar radiation parameterization
+ public :: surfrd_finetop_data ! Read topography dataset for fineTOP parameterization
!
! !PRIVATE MEMBER FUNCTIONS:
private :: surfrd_special ! Read the special landunits
@@ -1705,6 +1706,111 @@ subroutine surfrd_get_topo_for_solar_rad(domain,filename)
end subroutine surfrd_get_topo_for_solar_rad
+!-----------------------------------------------------------------------
+ subroutine surfrd_finetop_data(domain,filename)
+! !DESCRIPTION:
+! Read the topography parameters for fineTOP parameterization:
+! Assume domain has already been initialized and read
+
+! !USES:
+ use domainMod , only : domain_type
+ use fileutils , only : getfil
+ use GridcellType, only : grc_pp
+ use elm_varpar , only : ndir_horizon_angle
+
+! !ARGUMENTS:
+ implicit none
+ type(domain_type),intent(in) :: domain ! domain to init
+ character(len=*) ,intent(in) :: filename ! grid filename
+!
+! !CALLED FROM:
+! subroutine initialize
+!
+! !REVISION HISTORY:
+! Created by Dalei Hao
+!
+! !LOCAL VARIABLES:
+!EOP
+ type(file_desc_t) :: ncid ! netcdf file id
+ integer :: n ! indices
+ integer :: ni,nj,ns ! size of grid on file
+ integer :: dimid,varid ! netCDF id's
+ integer :: ier ! error status
+ real(r8) :: eps = 1.0e-12_r8 ! lat/lon error tolerance
+ integer :: beg,end ! local beg,end indices
+ logical :: isgrid2d ! true => file is 2d lat/lon
+ real(r8),pointer :: lonc(:),latc(:) ! local lat/lon
+ character(len=256) :: locfn ! local file name
+ logical :: readvar ! is variable on file
+ character(len=32) :: subname = 'surfrd_finetop_data' ! subroutine name
+!-----------------------------------------------------------------------
+
+ if (masterproc) then
+ if (filename == ' ') then
+ write(iulog,*) trim(subname),' ERROR: filename must be specified '
+ call endrun()
+ else
+ write(iulog,*) 'Attempting to read topography parameters from fsurdat ',trim(filename)
+ endif
+ end if
+
+ call getfil( filename, locfn, 0 )
+ call ncd_pio_openfile (ncid, trim(locfn), 0)
+ call ncd_inqfdims(ncid, isgrid2d, ni, nj, ns)
+
+ if (domain%ns /= ns) then
+ write(iulog,*) trim(subname),' ERROR: fsurdat file mismatch ns',&
+ domain%ns,ns
+ call endrun()
+ endif
+
+ beg = domain%nbeg
+ end = domain%nend
+
+ allocate(latc(beg:end),lonc(beg:end))
+
+ call ncd_io(ncid=ncid, varname='LONGXY', flag='read', data=lonc, &
+ dim1name=grlnd, readvar=readvar)
+ if (.not. readvar) call endrun( trim(subname)//' ERROR: LONGXY NOT on fsurdat file' )
+
+ call ncd_io(ncid=ncid, varname='LATIXY', flag='read', data=latc, &
+ dim1name=grlnd, readvar=readvar)
+ if (.not. readvar) call endrun( trim(subname)//' ERROR: LATIXY NOT on fsurdat file' )
+
+ do n = beg,end
+ if (abs(latc(n)-domain%latc(n)) > eps .or. &
+ abs(lonc(n)-domain%lonc(n)) > eps) then
+ write(iulog,*) trim(subname),' ERROR: fsurdat file mismatch lat,lon',latc(n),&
+ domain%latc(n),lonc(n),domain%lonc(n),eps
+ call endrun()
+ endif
+ enddo
+
+ call check_dim(ncid, 'ndir_horizon_angle', ndir_horizon_angle)
+
+ call ncd_io(ncid=ncid, varname='SLOPE_DEG', flag='read', data=grc_pp%slope_deg, &
+ dim1name=grlnd, readvar=readvar)
+ if (.not. readvar) call endrun( trim(subname)//' ERROR: slope_deg NOT on fsurdat file' )
+ call ncd_io(ncid=ncid, varname='ASPECT_DEG', flag='read', data=grc_pp%aspect_deg, &
+ dim1name=grlnd, readvar=readvar)
+ if (.not. readvar) call endrun( trim(subname)//' ERROR: aspect_deg NOT on fsurdat file' )
+ call ncd_io(ncid=ncid, varname='SKY_VIEW_FACTOR', flag='read', data=grc_pp%sky_view_factor, &
+ dim1name=grlnd, readvar=readvar)
+ if (.not. readvar) call endrun( trim(subname)//' ERROR: sky_view_factor NOT on fsurdat file' )
+ call ncd_io(ncid=ncid, varname='TERRAIN_CONFIG_FACTOR', flag='read', data=grc_pp%terrain_config_factor, &
+ dim1name=grlnd, readvar=readvar)
+ if (.not. readvar) call endrun( trim(subname)//' ERROR: terrain_config_factor NOT on fsurdat file' )
+ call ncd_io(ncid=ncid, varname='HORIZON_ANGLE_DEG', flag='read', data=grc_pp%horizon_angle_deg, &
+ dim1name=grlnd, readvar=readvar)
+ If (.not. readvar) call endrun( trim(subname)//' ERROR: horizon_angle_deg NOT on fsurdat file' )
+
+ deallocate(latc,lonc)
+
+ call ncd_pio_closefile(ncid)
+
+ end subroutine surfrd_finetop_data
+
+
subroutine surfrd_fates_nocropmod( ncid, begg, endg )
!--------------------------------------------------------------------------
diff --git a/share/util/shr_orb_mod.F90 b/share/util/shr_orb_mod.F90
index 18a5a91d028f..272662700ba5 100644
--- a/share/util/shr_orb_mod.F90
+++ b/share/util/shr_orb_mod.F90
@@ -11,6 +11,7 @@ MODULE shr_orb_mod
!----------------------------------------------------------------------------
! PUBLIC: Interfaces and global data
!----------------------------------------------------------------------------
+ public :: shr_orb_azimuth
public :: shr_orb_cosz
public :: shr_orb_params
public :: shr_orb_decl
@@ -47,6 +48,42 @@ SUBROUTINE set_constant_zenith_angle_deg(angle_deg)
END SUBROUTINE set_constant_zenith_angle_deg
!=======================================================================
+ !=======================================================================
+ real(SHR_KIND_R8) pure function shr_orb_azimuth(jday,lat,lon,declin,z)
+
+ !----------------------------------------------------------------------------
+ !
+ ! function returns the solar azimuth angle.
+ ! azimuth angle is defined with respect to north, positive to east
+ ! based on Sproul, Renewable Energy, 2007
+ !
+ !----------------------------------------------------------------------------
+ real (SHR_KIND_R8),intent(in) :: jday ! Julian cal day (1.xx to 365.xx)
+ real (SHR_KIND_R8),intent(in) :: lat ! Centered latitude (radians)
+ real (SHR_KIND_R8),intent(in) :: lon ! Centered longitude (radians)
+ real (SHR_KIND_R8),intent(in) :: declin ! Solar declination (radians)
+ real (SHR_KIND_R8),intent(in) :: z ! Solar zenith angle (radians)
+
+ real(SHR_KIND_R8) :: hour_angle
+ !----------------------------------------------------------------------------
+
+ hour_angle = 2.0_SHR_KIND_R8*pi*((jday-floor(jday)) - 0.5_SHR_KIND_R8) + lon
+
+ ! constrain hour_angle to [-pi,pi] to determine east/west below
+ if(hour_angle > pi) hour_angle = hour_angle - 2.0_SHR_KIND_R8*pi
+
+ shr_orb_azimuth = (sin(declin)*cos(lat) - cos(declin)*sin(lat)*cos(hour_angle))/ sin(z)
+
+ shr_orb_azimuth = max(-1._SHR_KIND_R8, min(shr_orb_azimuth, 1._SHR_KIND_R8))
+
+ shr_orb_azimuth = acos(shr_orb_azimuth)
+
+ ! azimuth is east for times between midnight and noon (h < 0)
+ ! azimuth is west for times between noon and midnight (h > 0)
+ if(hour_angle > 0.) shr_orb_azimuth = 2.0_SHR_KIND_R8*pi - shr_orb_azimuth
+
+ end function shr_orb_azimuth
+
!=======================================================================
real(SHR_KIND_R8) pure FUNCTION shr_orb_cosz(jday,lat,lon,declin,dt_avg,uniform_angle)