From 898148ed4c09a9a076de0594877ec0da882292b1 Mon Sep 17 00:00:00 2001 From: Dalei Hao Date: Wed, 19 Jun 2024 13:09:23 -0700 Subject: [PATCH 1/2] add finetop debug & define cosinc_col add model outputs; & modify surface area correction for longwave radiation & debug correct surface area for thermal radiation for soilT module update snowsnicar_AD model change COSINC output name to avoid duplicate fix a small bug in variable name fix small bugs debug debug fix bugs set SZA<85 limit fix small bugs fix one bug and add initialization for few variables fix small bugs delete log output and add control to adjust albedo delete log output delete write output add more outputs like MODIS daytime and nighttime ST debug fix a small bug add MODIS FSNO output fix a bug clean up the code in biogeophys clean up the code in main clean up the code fix urban bug add code to calculate apparent surface albedo for EAM fix a small bug debug debug & modify code to not change for urban landunits debug fix typo debug consider the top effects on upward radaiton to sky debug Conflicts: components/elm/src/biogeophys/SoilFluxesMod.F90 components/elm/src/biogeophys/SurfaceAlbedoMod.F90 components/elm/src/main/atm2lndType.F90 components/elm/src/main/elm_initializeMod.F90 fix a LAI bug delete ; debug & add test add restart variables fix a bug debug debug restart issues remove log output --- cime_config/tests.py | 3 +- .../bld/namelist_files/namelist_defaults.xml | 3 + .../namelist_files/namelist_definition.xml | 7 + .../elm/finetop_rad/shell_commands | 10 + .../testmods_dirs/elm/finetop_rad/user_nl_elm | 3 + .../elm/snowveg_arctic/shell_commands | 2 +- .../elm/src/biogeophys/CanopyFluxesMod.F90 | 39 ++-- .../elm/src/biogeophys/CanopyHydrologyMod.F90 | 4 +- .../src/biogeophys/CanopyTemperatureMod.F90 | 14 +- .../elm/src/biogeophys/SnowSnicarMod.F90 | 12 +- .../elm/src/biogeophys/SoilFluxesMod.F90 | 53 ++++-- .../elm/src/biogeophys/SoilTemperatureMod.F90 | 20 +- .../elm/src/biogeophys/SolarAbsorbedType.F90 | 36 +++- .../elm/src/biogeophys/SurfaceAlbedoMod.F90 | 112 ++++++++---- .../elm/src/biogeophys/SurfaceAlbedoType.F90 | 16 +- .../src/biogeophys/SurfaceRadiationMod.F90 | 38 ++-- .../elm/src/biogeophys/UrbanRadiationMod.F90 | 17 +- components/elm/src/cpl/lnd_comp_mct.F90 | 19 +- components/elm/src/cpl/lnd_import_export.F90 | 19 +- .../elm/src/data_types/GridcellType.F90 | 20 +- .../elm/src/data_types/TopounitDataType.F90 | 12 ++ components/elm/src/main/atm2lndMod.F90 | 171 +++++++++++++++++- components/elm/src/main/atm2lndType.F90 | 99 +++++++++- components/elm/src/main/controlMod.F90 | 16 +- components/elm/src/main/elm_driver.F90 | 9 +- components/elm/src/main/elm_initializeMod.F90 | 25 ++- components/elm/src/main/elm_varctl.F90 | 3 +- components/elm/src/main/elm_varpar.F90 | 2 + components/elm/src/main/lnd2atmMod.F90 | 95 +++++++++- components/elm/src/main/lnd2atmType.F90 | 28 ++- components/elm/src/main/surfrdMod.F90 | 106 +++++++++++ share/util/shr_orb_mod.F90 | 37 ++++ 32 files changed, 925 insertions(+), 125 deletions(-) create mode 100644 components/elm/cime_config/testdefs/testmods_dirs/elm/finetop_rad/shell_commands create mode 100644 components/elm/cime_config/testdefs/testmods_dirs/elm/finetop_rad/user_nl_elm 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..61397261e496 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 @@ -2266,6 +2269,13 @@ subroutine SNICAR_AD_RT (flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, & lon_coord = 0 endif ! end if flg_snw_ice == 1 + if (use_finetop_rad) then + slope_rad = grc_pp%slope_deg(g_idx) * deg2rad + h2osno_liq_lcl = h2osno_liq_lcl * cos(slope_rad) + h2osno_ice_lcl = h2osno_ice_lcl * cos(slope_rad) + h2osno_lcl = h2osno_lcl * cos(slope_rad) + endif + #ifdef MODAL_AER !mgf++ ! 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..2d3834656178 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 + deg2rad = SHR_CONST_PI/180._r8 do g = bounds%begg,bounds%endg coszen_gcell(g) = shr_orb_cosz (nextsw_cday, grc_pp%lat(g), grc_pp%lon(g), declinp1) + + if (use_finetop_rad) then + ! solar zenith angle + sza = acos(coszen_gcell(g)) + ! solar azimuth angle + saa = shr_orb_azimuth(nextsw_cday, grc_pp%lat(g), grc_pp%lon(g), declinp1, sza) + + 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)) + !write(iulog,*) 'cosinc_gcell',cosinc_gcell(g) + !write(iulog,*) 'coszen_gcell',coszen_gcell(g) + 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 + else + cosinc_gcell(g) = coszen_gcell(g) + endif end do 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..239972003af1 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) = nan + saa(g) = nan + cosinc(g) = nan + ! 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..2a650e1f321f 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 (:) = nan + allocate(this%f_short_dif (begg:endg)) ; this%f_short_dif (:) = nan + allocate(this%f_short_refl (begg:endg,numrad)) ; this%f_short_refl (:,:) = nan + allocate(this%f_long_dif (begg:endg)) ; this%f_long_dif (:) = nan + allocate(this%f_long_refl (begg:endg)) ; this%f_long_refl (:) = nan + allocate(this%sza (begg:endg)) ; this%sza (:) = nan + allocate(this%saa (begg:endg)) ; this%saa (:) = nan + allocate(this%cosinc (begg:endg)) ; this%cosinc (:) = nan + 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,58 @@ 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') + endif + + if ( use_finetop_rad ) then + this%forc_solad_pp_grc(begg:endg, :) = 0._r8 + this%forc_solai_pp_grc(begg:endg, :) = 0._r8 + endif + end subroutine InitHistory !----------------------------------------------------------------------- @@ -760,6 +841,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..9de22e72659a 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) @@ -1123,7 +1124,13 @@ subroutine control_print () else write(iulog,*) ' use_top_solar_rad is False, so do not run TOP solar radiation parameterization' end if - + + if (use_finetop_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 +1249,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) From 093e6e46ea9bcbe63058ae298bf5b8bed0f7ba48 Mon Sep 17 00:00:00 2001 From: daleihao Date: Fri, 31 Oct 2025 12:24:01 -0700 Subject: [PATCH 2/2] reformat and clean up the code fix a typo add error message fix a bug to avoid the initilization issue in the debug mode --- .../elm/src/biogeophys/SnowSnicarMod.F90 | 12 ++++-- .../elm/src/biogeophys/SurfaceAlbedoMod.F90 | 40 +++++++++---------- components/elm/src/main/atm2lndMod.F90 | 6 +-- components/elm/src/main/atm2lndType.F90 | 22 +++++----- components/elm/src/main/controlMod.F90 | 10 +++-- 5 files changed, 49 insertions(+), 41 deletions(-) diff --git a/components/elm/src/biogeophys/SnowSnicarMod.F90 b/components/elm/src/biogeophys/SnowSnicarMod.F90 index 61397261e496..157f540c8256 100644 --- a/components/elm/src/biogeophys/SnowSnicarMod.F90 +++ b/components/elm/src/biogeophys/SnowSnicarMod.F90 @@ -2249,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) @@ -2271,8 +2270,15 @@ subroutine SNICAR_AD_RT (flg_snw_ice, bounds, num_nourbanc, filter_nourbanc, & if (use_finetop_rad) then slope_rad = grc_pp%slope_deg(g_idx) * deg2rad - h2osno_liq_lcl = h2osno_liq_lcl * cos(slope_rad) - h2osno_ice_lcl = h2osno_ice_lcl * cos(slope_rad) + + 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 diff --git a/components/elm/src/biogeophys/SurfaceAlbedoMod.F90 b/components/elm/src/biogeophys/SurfaceAlbedoMod.F90 index 2d3834656178..95faf476a117 100644 --- a/components/elm/src/biogeophys/SurfaceAlbedoMod.F90 +++ b/components/elm/src/biogeophys/SurfaceAlbedoMod.F90 @@ -242,27 +242,27 @@ subroutine SurfaceAlbedo(bounds, & ! Cosine solar zenith angle for next time step deg2rad = SHR_CONST_PI/180._r8 - do g = bounds%begg,bounds%endg - coszen_gcell(g) = shr_orb_cosz (nextsw_cday, grc_pp%lat(g), grc_pp%lon(g), declinp1) + 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 - if (use_finetop_rad) then - ! solar zenith angle - sza = acos(coszen_gcell(g)) - ! solar azimuth angle - saa = shr_orb_azimuth(nextsw_cday, grc_pp%lat(g), grc_pp%lon(g), declinp1, sza) - - 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)) - !write(iulog,*) 'cosinc_gcell',cosinc_gcell(g) - !write(iulog,*) 'coszen_gcell',coszen_gcell(g) - 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 - else - cosinc_gcell(g) = coszen_gcell(g) - endif - end do do c = bounds%begc,bounds%endc g = col_pp%gridcell(c) coszen_col(c) = coszen_gcell(g) diff --git a/components/elm/src/main/atm2lndMod.F90 b/components/elm/src/main/atm2lndMod.F90 index 239972003af1..92b721ee9c70 100644 --- a/components/elm/src/main/atm2lndMod.F90 +++ b/components/elm/src/main/atm2lndMod.F90 @@ -540,9 +540,9 @@ subroutine topographic_effects_on_radiation(bounds, atm2lnd_vars, nextsw_cday, d f_short_dir(g) = 1._r8 f_short_dif(g) = 1._r8 f_short_refl(g,:) = 0._r8 - sza(g) = nan - saa(g) = nan - cosinc(g) = nan + 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 diff --git a/components/elm/src/main/atm2lndType.F90 b/components/elm/src/main/atm2lndType.F90 index 2a650e1f321f..d663fe085d75 100644 --- a/components/elm/src/main/atm2lndType.F90 +++ b/components/elm/src/main/atm2lndType.F90 @@ -331,14 +331,14 @@ subroutine InitAllocate(this, bounds) 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 (:) = nan - allocate(this%f_short_dif (begg:endg)) ; this%f_short_dif (:) = nan - allocate(this%f_short_refl (begg:endg,numrad)) ; this%f_short_refl (:,:) = nan - allocate(this%f_long_dif (begg:endg)) ; this%f_long_dif (:) = nan - allocate(this%f_long_refl (begg:endg)) ; this%f_long_refl (:) = nan - allocate(this%sza (begg:endg)) ; this%sza (:) = nan - allocate(this%saa (begg:endg)) ; this%saa (:) = nan - allocate(this%cosinc (begg:endg)) ; this%cosinc (:) = nan + 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 @@ -578,11 +578,9 @@ subroutine InitHistory(this, bounds) 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') - endif - if ( use_finetop_rad ) then - this%forc_solad_pp_grc(begg:endg, :) = 0._r8 - this%forc_solai_pp_grc(begg:endg, :) = 0._r8 + this%forc_solad_pp_grc(begg:endg, :) = 0._r8 + this%forc_solai_pp_grc(begg:endg, :) = 0._r8 endif end subroutine InitHistory diff --git a/components/elm/src/main/controlMod.F90 b/components/elm/src/main/controlMod.F90 index 9de22e72659a..ba23a9ab57b8 100755 --- a/components/elm/src/main/controlMod.F90 +++ b/components/elm/src/main/controlMod.F90 @@ -1122,13 +1122,17 @@ 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) then + 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' + write(iulog,*) ' use_finetop_rad is False, so do not run fineTOP radiation parameterization' end if if (use_cn) then