module UrbanMod !----------------------------------------------------------------------- !BOP ! ! !MODULE: UrbanMod ! ! !DESCRIPTION: ! Calculate solar and longwave radiation, and turbulent fluxes for urban landunit ! ! !USES: use shr_kind_mod, only : r8 => shr_kind_r8 use clm_varpar , only : numrad use clm_varcon , only : isecspday, degpsec use clm_varctl , only : iulog use abortutils , only : endrun use shr_sys_mod , only : shr_sys_flush ! ! !PUBLIC TYPES: implicit none save ! ! !PUBLIC MEMBER FUNCTIONS: public :: UrbanClumpInit ! Initialization of urban clump data structure public :: UrbanRadiation ! Urban radiative fluxes public :: UrbanAlbedo ! Urban albedos public :: UrbanSnowAlbedo ! Urban snow albedos public :: UrbanFluxes ! Urban turbulent fluxes ! !Urban control variables character(len= *), parameter, public :: urban_hac_off = 'OFF' ! character(len= *), parameter, public :: urban_hac_on = 'ON' ! character(len= *), parameter, public :: urban_wasteheat_on = 'ON_WASTEHEAT' ! character(len= 16), public :: urban_hac = urban_hac_off logical, public :: urban_traffic = .false. ! urban traffic fluxes ! ! !REVISION HISTORY: ! Created by Gordon Bonan and Mariana Vertenstein and Keith Oleson 04/2003 ! !EOP ! ! PRIVATE MEMBER FUNCTIONS private :: view_factor ! View factors for road and one wall private :: incident_direct ! Direct beam solar rad incident on walls and road in urban canyon private :: incident_diffuse ! Diffuse solar rad incident on walls and road in urban canyon private :: net_solar ! Solar radiation absorbed by road and both walls in urban canyon private :: net_longwave ! Net longwave radiation for road and both walls in urban canyon ! PRIVATE TYPES private type urban_clump_t real(r8), pointer :: canyon_hwr(:) ! ratio of building height to street width real(r8), pointer :: wtroad_perv(:) ! weight of pervious road wrt total road real(r8), pointer :: ht_roof(:) ! height of urban roof (m) real(r8), pointer :: wtlunit_roof(:) ! weight of roof with respect to landunit real(r8), pointer :: wind_hgt_canyon(:) ! height above road at which wind in canyon is to be computed (m) real(r8), pointer :: em_roof(:) ! roof emissivity real(r8), pointer :: em_improad(:) ! impervious road emissivity real(r8), pointer :: em_perroad(:) ! pervious road emissivity real(r8), pointer :: em_wall(:) ! wall emissivity real(r8), pointer :: alb_roof_dir(:,:) ! direct roof albedo real(r8), pointer :: alb_roof_dif(:,:) ! diffuse roof albedo real(r8), pointer :: alb_improad_dir(:,:) ! direct impervious road albedo real(r8), pointer :: alb_improad_dif(:,:) ! diffuse impervious road albedo real(r8), pointer :: alb_perroad_dir(:,:) ! direct pervious road albedo real(r8), pointer :: alb_perroad_dif(:,:) ! diffuse pervious road albedo real(r8), pointer :: alb_wall_dir(:,:) ! direct wall albedo real(r8), pointer :: alb_wall_dif(:,:) ! diffuse wall albedo end type urban_clump_t type (urban_clump_t), private, pointer :: urban_clump(:) ! array of urban clumps for this processor integer, private, parameter :: noonsec = isecspday / 2 ! seconds at local noon !----------------------------------------------------------------------- contains !----------------------------------------------------------------------- !BOP ! ! !IROUTINE: UrbanAlbedo ! ! !INTERFACE: subroutine UrbanAlbedo (nc, lbl, ubl, lbc, ubc, lbp, ubp, & num_urbanl, filter_urbanl, & num_urbanc, filter_urbanc, & num_urbanp, filter_urbanp) ! ! !DESCRIPTION: ! Determine urban landunit component albedos ! ! !USES: use clmtype use shr_orb_mod , only : shr_orb_decl, shr_orb_cosz use clm_varcon , only : icol_roof, icol_sunwall, icol_shadewall, icol_road_perv, icol_road_imperv, & sb ! ! !ARGUMENTS: implicit none integer , intent(in) :: nc ! clump index integer, intent(in) :: lbl, ubl ! landunit-index bounds integer, intent(in) :: lbc, ubc ! column-index bounds integer, intent(in) :: lbp, ubp ! pft-index bounds integer , intent(in) :: num_urbanl ! number of urban landunits in clump integer , intent(in) :: filter_urbanl(ubl-lbl+1) ! urban landunit filter integer , intent(in) :: num_urbanc ! number of urban columns in clump integer , intent(in) :: filter_urbanc(ubc-lbc+1) ! urban column filter integer , intent(in) :: num_urbanp ! number of urban pfts in clump integer , intent(in) :: filter_urbanp(ubp-lbp+1) ! urban pft filter ! ! !CALLED FROM: ! subroutine clm_driver1 ! ! !REVISION HISTORY: ! Author: Gordon Bonan ! 03/2003, Mariana Vertenstein: Migrated to clm2.2 ! 01/2008, Erik Kluzek: Migrated to clm3.5.15 ! ! !LOCAL VARIABLES: ! ! local pointers to original implicit in arguments ! integer , pointer :: pgridcell(:) ! gridcell of corresponding pft integer , pointer :: lgridcell(:) ! gridcell of corresponding landunit integer , pointer :: clandunit(:) ! column's landunit integer , pointer :: cgridcell(:) ! gridcell of corresponding column integer , pointer :: coli(:) ! beginning column index for landunit integer , pointer :: colf(:) ! ending column index for landunit integer , pointer :: ctype(:) ! column type integer , pointer :: pcolumn(:) ! column of corresponding pft real(r8), pointer :: czen(:) ! cosine of solar zenith angle for each column real(r8), pointer :: lat(:) ! latitude (radians) real(r8), pointer :: lon(:) ! longitude (radians) real(r8), pointer :: frac_sno(:) ! fraction of ground covered by snow (0 to 1) ! ! local pointers to original implicit out arguments ! real(r8), pointer :: albgrd(:,:) ! ground albedo (direct) real(r8), pointer :: albgri(:,:) ! ground albedo (diffuse) real(r8), pointer :: albd(:,:) ! surface albedo (direct) real(r8), pointer :: albi(:,:) ! surface albedo (diffuse) real(r8), pointer :: fabd(:,:) ! flux absorbed by veg per unit direct flux real(r8), pointer :: fabi(:,:) ! flux absorbed by veg per unit diffuse flux real(r8), pointer :: ftdd(:,:) ! down direct flux below veg per unit dir flx real(r8), pointer :: ftid(:,:) ! down diffuse flux below veg per unit dir flx real(r8), pointer :: ftii(:,:) ! down diffuse flux below veg per unit dif flx real(r8), pointer :: fsun(:) ! sunlit fraction of canopy real(r8), pointer :: gdir(:) ! leaf projection in solar direction (0 to 1) real(r8), pointer :: omega(:,:) ! fraction of intercepted radiation that is scattered (0 to 1) real(r8), pointer :: vf_sr(:) ! view factor of sky for road real(r8), pointer :: vf_wr(:) ! view factor of one wall for road real(r8), pointer :: vf_sw(:) ! view factor of sky for one wall real(r8), pointer :: vf_rw(:) ! view factor of road for one wall real(r8), pointer :: vf_ww(:) ! view factor of opposing wall for one wall real(r8), pointer :: sabs_roof_dir(:,:) ! direct solar absorbed by roof per unit ground area per unit incident flux real(r8), pointer :: sabs_roof_dif(:,:) ! diffuse solar absorbed by roof per unit ground area per unit incident flux real(r8), pointer :: sabs_sunwall_dir(:,:) ! direct solar absorbed by sunwall per unit wall area per unit incident flux real(r8), pointer :: sabs_sunwall_dif(:,:) ! diffuse solar absorbed by sunwall per unit wall area per unit incident flux real(r8), pointer :: sabs_shadewall_dir(:,:) ! direct solar absorbed by shadewall per unit wall area per unit incident flux real(r8), pointer :: sabs_shadewall_dif(:,:) ! diffuse solar absorbed by shadewall per unit wall area per unit incident flux real(r8), pointer :: sabs_improad_dir(:,:) ! direct solar absorbed by impervious road per unit ground area per unit incident flux real(r8), pointer :: sabs_improad_dif(:,:) ! diffuse solar absorbed by impervious road per unit ground area per unit incident flux real(r8), pointer :: sabs_perroad_dir(:,:) ! direct solar absorbed by pervious road per unit ground area per unit incident flux real(r8), pointer :: sabs_perroad_dif(:,:) ! diffuse solar absorbed by pervious road per unit ground area per unit incident flux ! ! ! !OTHER LOCAL VARIABLES !EOP ! real(r8) :: coszen(num_urbanl) ! cosine solar zenith angle real(r8) :: coszen_pft(num_urbanp) ! cosine solar zenith angle for next time step (pft level) real(r8) :: zen(num_urbanl) ! solar zenith angle (radians) real(r8) :: sdir(num_urbanl, numrad) ! direct beam solar radiation on horizontal surface real(r8) :: sdif(num_urbanl, numrad) ! diffuse solar radiation on horizontal surface real(r8) :: sdir_road(num_urbanl, numrad) ! direct beam solar radiation incident on road real(r8) :: sdif_road(num_urbanl, numrad) ! diffuse solar radiation incident on road real(r8) :: sdir_sunwall(num_urbanl, numrad) ! direct beam solar radiation (per unit wall area) incident on sunlit wall per unit incident flux real(r8) :: sdif_sunwall(num_urbanl, numrad) ! diffuse solar radiation (per unit wall area) incident on sunlit wall per unit incident flux real(r8) :: sdir_shadewall(num_urbanl, numrad) ! direct beam solar radiation (per unit wall area) incident on shaded wall per unit incident flux real(r8) :: sdif_shadewall(num_urbanl, numrad) ! diffuse solar radiation (per unit wall area) incident on shaded wall per unit incident flux real(r8) :: albsnd_roof(num_urbanl,numrad) ! snow albedo for roof (direct) real(r8) :: albsni_roof(num_urbanl,numrad) ! snow albedo for roof (diffuse) real(r8) :: albsnd_improad(num_urbanl,numrad) ! snow albedo for impervious road (direct) real(r8) :: albsni_improad(num_urbanl,numrad) ! snow albedo for impervious road (diffuse) real(r8) :: albsnd_perroad(num_urbanl,numrad) ! snow albedo for pervious road (direct) real(r8) :: albsni_perroad(num_urbanl,numrad) ! snow albedo for pervious road (diffuse) integer :: fl,fp,fc,g,l,p,c,ib ! indices integer :: ic ! 0=unit incoming direct; 1=unit incoming diffuse integer :: num_solar ! counter real(r8) :: alb_roof_dir_s(num_urbanl,numrad) ! direct roof albedo with snow effects real(r8) :: alb_roof_dif_s(num_urbanl,numrad) ! diffuse roof albedo with snow effects real(r8) :: alb_improad_dir_s(num_urbanl,numrad) ! direct impervious road albedo with snow effects real(r8) :: alb_perroad_dir_s(num_urbanl,numrad) ! direct pervious road albedo with snow effects real(r8) :: alb_improad_dif_s(num_urbanl,numrad) ! diffuse impervious road albedo with snow effects real(r8) :: alb_perroad_dif_s(num_urbanl,numrad) ! diffuse pervious road albedo with snow effects real(r8) :: sref_roof_dir(num_urbanl,numrad) ! direct solar reflected by roof per unit ground area per unit incident flux real(r8) :: sref_roof_dif(num_urbanl,numrad) ! diffuse solar reflected by roof per unit ground area per unit incident flux real(r8) :: sref_sunwall_dir(num_urbanl,numrad) ! direct solar reflected by sunwall per unit wall area per unit incident flux real(r8) :: sref_sunwall_dif(num_urbanl,numrad) ! diffuse solar reflected by sunwall per unit wall area per unit incident flux real(r8) :: sref_shadewall_dir(num_urbanl,numrad) ! direct solar reflected by shadewall per unit wall area per unit incident flux real(r8) :: sref_shadewall_dif(num_urbanl,numrad) ! diffuse solar reflected by shadewall per unit wall area per unit incident flux real(r8) :: sref_improad_dir(num_urbanl,numrad) ! direct solar reflected by impervious road per unit ground area per unit incident flux real(r8) :: sref_improad_dif(num_urbanl,numrad) ! diffuse solar reflected by impervious road per unit ground area per unit incident flux real(r8) :: sref_perroad_dir(num_urbanl,numrad) ! direct solar reflected by pervious road per unit ground area per unit incident flux real(r8) :: sref_perroad_dif(num_urbanl,numrad) ! diffuse solar reflected by pervious road per unit ground area per unit incident flux real(r8), pointer :: canyon_hwr(:) ! ratio of building height to street width real(r8), pointer :: wtroad_perv(:) ! weight of pervious road wrt total road real(r8), pointer :: alb_roof_dir(:,:) ! direct roof albedo real(r8), pointer :: alb_roof_dif(:,:) ! diffuse roof albedo real(r8), pointer :: alb_improad_dir(:,:) ! direct impervious road albedo real(r8), pointer :: alb_perroad_dir(:,:) ! direct pervious road albedo real(r8), pointer :: alb_improad_dif(:,:) ! diffuse imprevious road albedo real(r8), pointer :: alb_perroad_dif(:,:) ! diffuse pervious road albedo real(r8), pointer :: alb_wall_dir(:,:) ! direct wall albedo real(r8), pointer :: alb_wall_dif(:,:) ! diffuse wall albedo !----------------------------------------------------------------------- ! Assign pointers into module urban clumps canyon_hwr => urban_clump(nc)%canyon_hwr wtroad_perv => urban_clump(nc)%wtroad_perv alb_roof_dir => urban_clump(nc)%alb_roof_dir alb_roof_dif => urban_clump(nc)%alb_roof_dif alb_improad_dir => urban_clump(nc)%alb_improad_dir alb_improad_dif => urban_clump(nc)%alb_improad_dif alb_perroad_dir => urban_clump(nc)%alb_perroad_dir alb_perroad_dif => urban_clump(nc)%alb_perroad_dif alb_wall_dir => urban_clump(nc)%alb_wall_dir alb_wall_dif => urban_clump(nc)%alb_wall_dif ! Assign gridcell level pointers lat => clm3%g%lat lon => clm3%g%lon ! Assign landunit level pointer lgridcell => clm3%g%l%gridcell coli => clm3%g%l%coli colf => clm3%g%l%colf vf_sr => clm3%g%l%lps%vf_sr vf_wr => clm3%g%l%lps%vf_wr vf_sw => clm3%g%l%lps%vf_sw vf_rw => clm3%g%l%lps%vf_rw vf_ww => clm3%g%l%lps%vf_ww sabs_roof_dir => clm3%g%l%lps%sabs_roof_dir sabs_roof_dif => clm3%g%l%lps%sabs_roof_dif sabs_sunwall_dir => clm3%g%l%lps%sabs_sunwall_dir sabs_sunwall_dif => clm3%g%l%lps%sabs_sunwall_dif sabs_shadewall_dir => clm3%g%l%lps%sabs_shadewall_dir sabs_shadewall_dif => clm3%g%l%lps%sabs_shadewall_dif sabs_improad_dir => clm3%g%l%lps%sabs_improad_dir sabs_improad_dif => clm3%g%l%lps%sabs_improad_dif sabs_perroad_dir => clm3%g%l%lps%sabs_perroad_dir sabs_perroad_dif => clm3%g%l%lps%sabs_perroad_dif ! Assign column level pointers ctype => clm3%g%l%c%itype albgrd => clm3%g%l%c%cps%albgrd albgri => clm3%g%l%c%cps%albgri frac_sno => clm3%g%l%c%cps%frac_sno clandunit => clm3%g%l%c%landunit cgridcell => clm3%g%l%c%gridcell czen => clm3%g%l%c%cps%coszen ! Assign pft level pointers pgridcell => clm3%g%l%c%p%gridcell pcolumn => clm3%g%l%c%p%column albd => clm3%g%l%c%p%pps%albd albi => clm3%g%l%c%p%pps%albi fabd => clm3%g%l%c%p%pps%fabd fabi => clm3%g%l%c%p%pps%fabi ftdd => clm3%g%l%c%p%pps%ftdd ftid => clm3%g%l%c%p%pps%ftid ftii => clm3%g%l%c%p%pps%ftii fsun => clm3%g%l%c%p%pps%fsun gdir => clm3%g%l%c%p%pps%gdir omega => clm3%g%l%c%p%pps%omega ! ---------------------------------------------------------------------------- ! Solar declination and cosine solar zenith angle and zenith angle for ! next time step ! ---------------------------------------------------------------------------- do fl = 1,num_urbanl l = filter_urbanl(fl) g = lgridcell(l) coszen(fl) = czen(coli(l)) ! Assumes coszen for each column are the same zen(fl) = acos(coszen(fl)) end do do fp = 1,num_urbanp p = filter_urbanp(fp) g = pgridcell(p) c = pcolumn(p) coszen_pft(fp) = czen(c) end do ! ---------------------------------------------------------------------------- ! Initialize clmtype output since solar radiation is only done if coszen > 0 ! ---------------------------------------------------------------------------- do ib = 1,numrad do fc = 1,num_urbanc c = filter_urbanc(fc) albgrd(c,ib) = 0._r8 albgri(c,ib) = 0._r8 end do do fp = 1,num_urbanp p = filter_urbanp(fp) g = pgridcell(p) albd(p,ib) = 1._r8 albi(p,ib) = 1._r8 fabd(p,ib) = 0._r8 fabi(p,ib) = 0._r8 if (coszen_pft(fp) > 0._r8) then ftdd(p,ib) = 1._r8 else ftdd(p,ib) = 0._r8 end if ftid(p,ib) = 0._r8 if (coszen_pft(fp) > 0._r8) then ftii(p,ib) = 1._r8 else ftii(p,ib) = 0._r8 end if omega(p,ib) = 0._r8 if (ib == 1) then gdir(p) = 0._r8 fsun(p) = 0._r8 end if end do end do ! ---------------------------------------------------------------------------- ! Urban Code ! ---------------------------------------------------------------------------- num_solar = 0 do fl = 1,num_urbanl if (coszen(fl) > 0._r8) num_solar = num_solar + 1 end do ! Initialize urban clump components do ib = 1,numrad do fl = 1,num_urbanl l = filter_urbanl(fl) sabs_roof_dir(l,ib) = 0._r8 sabs_roof_dif(l,ib) = 0._r8 sabs_sunwall_dir(l,ib) = 0._r8 sabs_sunwall_dif(l,ib) = 0._r8 sabs_shadewall_dir(l,ib) = 0._r8 sabs_shadewall_dif(l,ib) = 0._r8 sabs_improad_dir(l,ib) = 0._r8 sabs_improad_dif(l,ib) = 0._r8 sabs_perroad_dir(l,ib) = 0._r8 sabs_perroad_dif(l,ib) = 0._r8 sref_roof_dir(fl,ib) = 1._r8 sref_roof_dif(fl,ib) = 1._r8 sref_sunwall_dir(fl,ib) = 1._r8 sref_sunwall_dif(fl,ib) = 1._r8 sref_shadewall_dir(fl,ib) = 1._r8 sref_shadewall_dif(fl,ib) = 1._r8 sref_improad_dir(fl,ib) = 1._r8 sref_improad_dif(fl,ib) = 1._r8 sref_perroad_dir(fl,ib) = 1._r8 sref_perroad_dif(fl,ib) = 1._r8 end do end do ! View factors for road and one wall in urban canyon (depends only on canyon_hwr) if (num_urbanl .gt. 0) then call view_factor (lbl, ubl, num_urbanl, filter_urbanl, canyon_hwr) end if ! ---------------------------------------------------------------------------- ! Only do the rest if all coszen are positive ! ---------------------------------------------------------------------------- if (num_solar > 0)then ! Set constants - solar fluxes are per unit incoming flux do ib = 1,numrad do fl = 1,num_urbanl sdir(fl,ib) = 1._r8 sdif(fl,ib) = 1._r8 end do end do ! Incident direct beam radiation for ! (a) roof and (b) road and both walls in urban canyon if (num_urbanl .gt. 0) then call incident_direct (lbl, ubl, num_urbanl, canyon_hwr, coszen, zen, sdir, sdir_road, sdir_sunwall, sdir_shadewall) end if ! Incident diffuse radiation for ! (a) roof and (b) road and both walls in urban canyon. if (num_urbanl .gt. 0) then call incident_diffuse (lbl, ubl, num_urbanl, filter_urbanl, canyon_hwr, sdif, sdif_road, & sdif_sunwall, sdif_shadewall) end if ! Get snow albedos for roof and impervious and pervious road if (num_urbanl .gt. 0) then ic = 0; call UrbanSnowAlbedo(lbl, ubl, num_urbanl, filter_urbanl, coszen, ic, albsnd_roof, albsnd_improad, albsnd_perroad) ic = 1; call UrbanSnowAlbedo(lbl, ubl, num_urbanl, filter_urbanl, coszen, ic, albsni_roof, albsni_improad, albsni_perroad) end if ! Combine snow-free and snow albedos do ib = 1,numrad do fl = 1,num_urbanl l = filter_urbanl(fl) do c = coli(l),colf(l) if (ctype(c) == icol_roof) then alb_roof_dir_s(fl,ib) = alb_roof_dir(fl,ib)*(1._r8-frac_sno(c)) & + albsnd_roof(fl,ib)*frac_sno(c) alb_roof_dif_s(fl,ib) = alb_roof_dif(fl,ib)*(1._r8-frac_sno(c)) & + albsni_roof(fl,ib)*frac_sno(c) else if (ctype(c) == icol_road_imperv) then alb_improad_dir_s(fl,ib) = alb_improad_dir(fl,ib)*(1._r8-frac_sno(c)) & + albsnd_improad(fl,ib)*frac_sno(c) alb_improad_dif_s(fl,ib) = alb_improad_dif(fl,ib)*(1._r8-frac_sno(c)) & + albsni_improad(fl,ib)*frac_sno(c) else if (ctype(c) == icol_road_perv) then alb_perroad_dir_s(fl,ib) = alb_perroad_dir(fl,ib)*(1._r8-frac_sno(c)) & + albsnd_perroad(fl,ib)*frac_sno(c) alb_perroad_dif_s(fl,ib) = alb_perroad_dif(fl,ib)*(1._r8-frac_sno(c)) & + albsni_perroad(fl,ib)*frac_sno(c) end if end do end do end do ! Reflected and absorbed solar radiation per unit incident radiation ! for road and both walls in urban canyon allowing for multiple reflection ! Reflected and absorbed solar radiation per unit incident radiation for roof if (num_urbanl .gt. 0) then call net_solar (lbl, ubl, num_urbanl, filter_urbanl, coszen, canyon_hwr, wtroad_perv, sdir, sdif, & alb_improad_dir_s, alb_perroad_dir_s, alb_wall_dir, alb_roof_dir_s, & alb_improad_dif_s, alb_perroad_dif_s, alb_wall_dif, alb_roof_dif_s, & sdir_road, sdir_sunwall, sdir_shadewall, & sdif_road, sdif_sunwall, sdif_shadewall, & sref_improad_dir, sref_perroad_dir, sref_sunwall_dir, sref_shadewall_dir, sref_roof_dir, & sref_improad_dif, sref_perroad_dif, sref_sunwall_dif, sref_shadewall_dif, sref_roof_dif) end if ! ---------------------------------------------------------------------------- ! Map urban output to clmtype components ! ---------------------------------------------------------------------------- ! Set albgrd and albgri (ground albedos) and albd and albi (surface albedos) do ib = 1,numrad do fl = 1,num_urbanl l = filter_urbanl(fl) do c = coli(l),colf(l) if (ctype(c) == icol_roof) then albgrd(c,ib) = sref_roof_dir(fl,ib) albgri(c,ib) = sref_roof_dif(fl,ib) else if (ctype(c) == icol_sunwall) then albgrd(c,ib) = sref_sunwall_dir(fl,ib) albgri(c,ib) = sref_sunwall_dif(fl,ib) else if (ctype(c) == icol_shadewall) then albgrd(c,ib) = sref_shadewall_dir(fl,ib) albgri(c,ib) = sref_shadewall_dif(fl,ib) else if (ctype(c) == icol_road_perv) then albgrd(c,ib) = sref_perroad_dir(fl,ib) albgri(c,ib) = sref_perroad_dif(fl,ib) else if (ctype(c) == icol_road_imperv) then albgrd(c,ib) = sref_improad_dir(fl,ib) albgri(c,ib) = sref_improad_dif(fl,ib) endif end do end do do fp = 1,num_urbanp p = filter_urbanp(fp) c = pcolumn(p) albd(p,ib) = albgrd(c,ib) albi(p,ib) = albgri(c,ib) end do end do end if end subroutine UrbanAlbedo !----------------------------------------------------------------------- !BOP ! ! !IROUTINE: UrbanSnowAlbedo ! ! !INTERFACE: subroutine UrbanSnowAlbedo (lbl, ubl, num_urbanl, filter_urbanl, coszen, ind, & albsn_roof, albsn_improad, albsn_perroad) ! ! !DESCRIPTION: ! Determine urban snow albedos ! ! !USES: use shr_kind_mod, only: r8 => shr_kind_r8 use clmtype use clm_varcon , only : icol_roof, icol_road_perv, icol_road_imperv ! ! !ARGUMENTS: implicit none integer, intent(in) :: lbl, ubl ! landunit-index bounds integer , intent(in) :: num_urbanl ! number of urban landunits in clump integer , intent(in) :: filter_urbanl(ubl-lbl+1) ! urban landunit filter integer , intent(in) :: ind ! 0=direct beam, 1=diffuse radiation real(r8), intent(in) :: coszen(num_urbanl) ! cosine solar zenith angle real(r8), intent(out):: albsn_roof(num_urbanl,2) ! roof snow albedo by waveband (assume 2 wavebands) real(r8), intent(out):: albsn_improad(num_urbanl,2) ! impervious road snow albedo by waveband (assume 2 wavebands) real(r8), intent(out):: albsn_perroad(num_urbanl,2) ! pervious road snow albedo by waveband (assume 2 wavebands) ! ! !CALLED FROM: ! subroutine UrbanAlbedo in this module ! ! !REVISION HISTORY: ! Author: Keith Oleson 9/2005 ! ! !LOCAL VARIABLES: ! ! local pointers to implicit in arguments integer , pointer :: coli(:) ! beginning column index for landunit integer , pointer :: colf(:) ! ending column index for landunit real(r8), pointer :: h2osno(:) ! snow water (mm H2O) integer , pointer :: ctype(:) ! column type ! ! ! !OTHER LOCAL VARIABLES: !EOP integer :: fl,c,l ! indices ! ! variables and constants for snow albedo calculation ! ! These values are derived from Marshall (1989) assuming soot content of 1.5e-5 ! (three times what LSM uses globally). Note that snow age effects are ignored here. real(r8), parameter :: snal0 = 0.66_r8 ! vis albedo of urban snow real(r8), parameter :: snal1 = 0.56_r8 ! nir albedo of urban snow !----------------------------------------------------------------------- ! Assign local pointers to derived type members (landunit level) coli => clm3%g%l%coli colf => clm3%g%l%colf ! Assign local pointers to derived subtypes components (column-level) ctype => clm3%g%l%c%itype h2osno => clm3%g%l%c%cws%h2osno ! this code assumes that numrad = 2 , with the following ! index values: 1 = visible, 2 = NIR do fl = 1,num_urbanl l = filter_urbanl(fl) do c = coli(l),colf(l) if (coszen(fl) > 0._r8 .and. h2osno(c) > 0._r8) then if (ctype(c) == icol_roof) then albsn_roof(fl,1) = snal0 albsn_roof(fl,2) = snal1 else if (ctype(c) == icol_road_imperv) then albsn_improad(fl,1) = snal0 albsn_improad(fl,2) = snal1 else if (ctype(c) == icol_road_perv) then albsn_perroad(fl,1) = snal0 albsn_perroad(fl,2) = snal1 end if else if (ctype(c) == icol_roof) then albsn_roof(fl,1) = 0._r8 albsn_roof(fl,2) = 0._r8 else if (ctype(c) == icol_road_imperv) then albsn_improad(fl,1) = 0._r8 albsn_improad(fl,2) = 0._r8 else if (ctype(c) == icol_road_perv) then albsn_perroad(fl,1) = 0._r8 albsn_perroad(fl,2) = 0._r8 end if end if end do end do end subroutine UrbanSnowAlbedo !----------------------------------------------------------------------- !BOP ! ! !IROUTINE: UrbanRadiation ! ! !INTERFACE: subroutine UrbanRadiation (nc, lbl, ubl, lbc, ubc, lbp, ubp, & num_nourbanl, filter_nourbanl, & num_urbanl, filter_urbanl, & num_urbanc, filter_urbanc, & num_urbanp, filter_urbanp) ! ! !DESCRIPTION: ! Solar fluxes absorbed and reflected by roof and canyon (walls, road). ! Also net and upward longwave fluxes. ! !USES: use clmtype use clm_varcon , only : spval, icol_roof, icol_sunwall, icol_shadewall, & icol_road_perv, icol_road_imperv, sb use clm_varcon , only : tfrz ! To use new constant.. use clm_time_manager , only : get_curr_date, get_step_size use clm_atmlnd , only : clm_a2l ! ! !ARGUMENTS: implicit none integer , intent(in) :: nc ! clump index integer, intent(in) :: lbl, ubl ! landunit-index bounds integer, intent(in) :: lbc, ubc ! column-index bounds integer, intent(in) :: lbp, ubp ! pft-index bounds integer , intent(in) :: num_nourbanl ! number of non-urban landunits in clump integer , intent(in) :: filter_nourbanl(ubl-lbl+1) ! non-urban landunit filter integer , intent(in) :: num_urbanl ! number of urban landunits in clump integer , intent(in) :: filter_urbanl(ubl-lbl+1) ! urban landunit filter integer , intent(in) :: num_urbanc ! number of urban columns in clump integer , intent(in) :: filter_urbanc(ubc-lbc+1) ! urban column filter integer , intent(in) :: num_urbanp ! number of urban pfts in clump integer , intent(in) :: filter_urbanp(ubp-lbp+1) ! urban pft filter ! ! !CALLED FROM: ! subroutine clm_driver1 ! ! !REVISION HISTORY: ! Author: Gordon Bonan ! 03/2003, Mariana Vertenstein: Migrated to clm2.2 ! 07/2004, Mariana Vertenstein: Migrated to clm3.0 ! 01/2008, Erik Kluzek: Migrated to clm3.5.15 ! ! !LOCAL VARIABLES: ! ! local pointers to original implicit in arguments (urban clump) ! real(r8), pointer :: canyon_hwr(:) ! ratio of building height to street width real(r8), pointer :: wtroad_perv(:) ! weight of pervious road wrt total road real(r8), pointer :: em_roof(:) ! roof emissivity real(r8), pointer :: em_improad(:) ! impervious road emissivity real(r8), pointer :: em_perroad(:) ! pervious road emissivity real(r8), pointer :: em_wall(:) ! wall emissivity ! ! local pointers to original implicit in arguments (clmtype) ! integer , pointer :: pgridcell(:) ! gridcell of corresponding pft integer , pointer :: pcolumn(:) ! column of corresponding pft integer , pointer :: lgridcell(:) ! gridcell of corresponding landunit integer , pointer :: ctype(:) ! column type integer , pointer :: coli(:) ! beginning column index for landunit integer , pointer :: colf(:) ! ending column index for landunit integer , pointer :: pfti(:) ! beginning pfti index for landunit integer , pointer :: pftf(:) ! ending pftf index for landunit real(r8), pointer :: londeg(:) ! longitude (degrees) real(r8), pointer :: forc_lwrad(:) ! downward infrared (longwave) radiation (W/m**2) real(r8), pointer :: forc_solad(:,:) ! direct beam radiation (vis=forc_sols , nir=forc_soll ) (W/m**2) real(r8), pointer :: forc_solai(:,:) ! diffuse beam radiation (vis=forc_sols , nir=forc_soll ) (W/m**2) real(r8), pointer :: forc_solar(:) ! incident solar radiation (W/m**2) real(r8), pointer :: albd(:,:) ! surface albedo (direct) real(r8), pointer :: albi(:,:) ! surface albedo (diffuse) real(r8), pointer :: t_grnd(:) ! ground temperature (K) real(r8), pointer :: frac_sno(:) ! fraction of ground covered by snow (0 to 1) real(r8), pointer :: t_ref2m(:) ! 2 m height surface air temperature (K) real(r8), pointer :: vf_sr(:) ! view factor of sky for road real(r8), pointer :: vf_wr(:) ! view factor of one wall for road real(r8), pointer :: vf_sw(:) ! view factor of sky for one wall real(r8), pointer :: vf_rw(:) ! view factor of road for one wall real(r8), pointer :: vf_ww(:) ! view factor of opposing wall for one wall real(r8), pointer :: sabs_roof_dir(:,:) ! direct solar absorbed by roof per unit ground area per unit incident flux real(r8), pointer :: sabs_roof_dif(:,:) ! diffuse solar absorbed by roof per unit ground area per unit incident flux real(r8), pointer :: sabs_sunwall_dir(:,:) ! direct solar absorbed by sunwall per unit wall area per unit incident flux real(r8), pointer :: sabs_sunwall_dif(:,:) ! diffuse solar absorbed by sunwall per unit wall area per unit incident flux real(r8), pointer :: sabs_shadewall_dir(:,:) ! direct solar absorbed by shadewall per unit wall area per unit incident flux real(r8), pointer :: sabs_shadewall_dif(:,:) ! diffuse solar absorbed by shadewall per unit wall area per unit incident flux real(r8), pointer :: sabs_improad_dir(:,:) ! direct solar absorbed by impervious road per unit ground area per unit incident flux real(r8), pointer :: sabs_improad_dif(:,:) ! diffuse solar absorbed by impervious road per unit ground area per unit incident flux real(r8), pointer :: sabs_perroad_dir(:,:) ! direct solar absorbed by pervious road per unit ground area per unit incident flux real(r8), pointer :: sabs_perroad_dif(:,:) ! diffuse solar absorbed by pervious road per unit ground area per unit incident flux ! ! local pointers to original implicit out arguments (clmtype) ! real(r8), pointer :: parsun(:) ! average absorbed PAR for sunlit leaves (W/m**2) real(r8), pointer :: parsha(:) ! average absorbed PAR for shaded leaves (W/m**2) real(r8), pointer :: sabg(:) ! solar radiation absorbed by ground (W/m**2) real(r8), pointer :: sabv(:) ! solar radiation absorbed by vegetation (W/m**2) real(r8), pointer :: fsa(:) ! solar radiation absorbed (total) (W/m**2) real(r8), pointer :: fsa_u(:) ! urban solar radiation absorbed (total) (W/m**2) real(r8), pointer :: fsr(:) ! solar radiation reflected (total) (W/m**2) real(r8), pointer :: fsds_vis_d(:) ! incident direct beam vis solar radiation (W/m**2) real(r8), pointer :: fsds_nir_d(:) ! incident direct beam nir solar radiation (W/m**2) real(r8), pointer :: fsds_vis_i(:) ! incident diffuse vis solar radiation (W/m**2) real(r8), pointer :: fsds_nir_i(:) ! incident diffuse nir solar radiation (W/m**2) real(r8), pointer :: fsr_vis_d(:) ! reflected direct beam vis solar radiation (W/m**2) real(r8), pointer :: fsr_nir_d(:) ! reflected direct beam nir solar radiation (W/m**2) real(r8), pointer :: fsr_vis_i(:) ! reflected diffuse vis solar radiation (W/m**2) real(r8), pointer :: fsr_nir_i(:) ! reflected diffuse nir solar radiation (W/m**2) real(r8), pointer :: fsds_vis_d_ln(:) ! incident direct beam vis solar rad at local noon (W/m**2) real(r8), pointer :: fsds_nir_d_ln(:) ! incident direct beam nir solar rad at local noon (W/m**2) real(r8), pointer :: fsr_vis_d_ln(:) ! reflected direct beam vis solar rad at local noon (W/m**2) real(r8), pointer :: fsr_nir_d_ln(:) ! reflected direct beam nir solar rad at local noon (W/m**2) real(r8), pointer :: eflx_lwrad_out(:) ! emitted infrared (longwave) radiation (W/m**2) real(r8), pointer :: eflx_lwrad_net(:) ! net infrared (longwave) rad (W/m**2) [+ = to atm] real(r8), pointer :: eflx_lwrad_net_u(:) ! urban net infrared (longwave) rad (W/m**2) [+ = to atm] ! ! ! !OTHER LOCAL VARIABLES !EOP ! integer :: fp,fl,p,c,l,g ! indices integer :: local_secp1 ! seconds into current date in local time real(r8) :: dtime ! land model time step (sec) integer :: year,month,day ! temporaries (not used) integer :: secs ! seconds into current date real(r8), parameter :: mpe = 1.e-06_r8 ! prevents overflow for division by zero real(r8), parameter :: snoem = 0.97_r8 ! snow emissivity (should use value from Biogeophysics1) real(r8) :: lwnet_roof(num_urbanl) ! net (outgoing-incoming) longwave radiation (per unit ground area), roof (W/m**2) real(r8) :: lwnet_improad(num_urbanl) ! net (outgoing-incoming) longwave radiation (per unit ground area), impervious road (W/m**2) real(r8) :: lwnet_perroad(num_urbanl) ! net (outgoing-incoming) longwave radiation (per unit ground area), pervious road (W/m**2) real(r8) :: lwnet_sunwall(num_urbanl) ! net (outgoing-incoming) longwave radiation (per unit wall area), sunlit wall (W/m**2) real(r8) :: lwnet_shadewall(num_urbanl)! net (outgoing-incoming) longwave radiation (per unit wall area), shaded wall (W/m**2) real(r8) :: lwnet_canyon(num_urbanl) ! net (outgoing-incoming) longwave radiation for canyon, per unit ground area (W/m**2) real(r8) :: lwup_roof(num_urbanl) ! upward longwave radiation (per unit ground area), roof (W/m**2) real(r8) :: lwup_improad(num_urbanl) ! upward longwave radiation (per unit ground area), impervious road (W/m**2) real(r8) :: lwup_perroad(num_urbanl) ! upward longwave radiation (per unit ground area), pervious road (W/m**2) real(r8) :: lwup_sunwall(num_urbanl) ! upward longwave radiation, (per unit wall area), sunlit wall (W/m**2) real(r8) :: lwup_shadewall(num_urbanl) ! upward longwave radiation, (per unit wall area), shaded wall (W/m**2) real(r8) :: lwup_canyon(num_urbanl) ! upward longwave radiation for canyon, per unit ground area (W/m**2) real(r8) :: t_roof(num_urbanl) ! roof temperature (K) real(r8) :: t_improad(num_urbanl) ! imppervious road temperature (K) real(r8) :: t_perroad(num_urbanl) ! pervious road temperature (K) real(r8) :: t_sunwall(num_urbanl) ! sunlit wall temperature (K) real(r8) :: t_shadewall(num_urbanl) ! shaded wall temperature (K) real(r8) :: lwdown(num_urbanl) ! atmospheric downward longwave radiation (W/m**2) real(r8) :: em_roof_s(num_urbanl) ! roof emissivity with snow effects real(r8) :: em_improad_s(num_urbanl) ! impervious road emissivity with snow effects real(r8) :: em_perroad_s(num_urbanl) ! pervious road emissivity with snow effects !----------------------------------------------------------------------- ! Assign pointers into module urban clumps if( num_urbanl > 0 )then canyon_hwr => urban_clump(nc)%canyon_hwr wtroad_perv => urban_clump(nc)%wtroad_perv em_roof => urban_clump(nc)%em_roof em_improad => urban_clump(nc)%em_improad em_perroad => urban_clump(nc)%em_perroad em_wall => urban_clump(nc)%em_wall end if ! Assign local pointers to multi-level derived type members (gridcell level) londeg => clm3%g%londeg forc_solad => clm_a2l%forc_solad forc_solai => clm_a2l%forc_solai forc_solar => clm_a2l%forc_solar forc_lwrad => clm_a2l%forc_lwrad ! Assign local pointers to derived type members (landunit level) pfti => clm3%g%l%pfti pftf => clm3%g%l%pftf coli => clm3%g%l%coli colf => clm3%g%l%colf lgridcell => clm3%g%l%gridcell vf_sr => clm3%g%l%lps%vf_sr vf_wr => clm3%g%l%lps%vf_wr vf_sw => clm3%g%l%lps%vf_sw vf_rw => clm3%g%l%lps%vf_rw vf_ww => clm3%g%l%lps%vf_ww sabs_roof_dir => clm3%g%l%lps%sabs_roof_dir sabs_roof_dif => clm3%g%l%lps%sabs_roof_dif sabs_sunwall_dir => clm3%g%l%lps%sabs_sunwall_dir sabs_sunwall_dif => clm3%g%l%lps%sabs_sunwall_dif sabs_shadewall_dir => clm3%g%l%lps%sabs_shadewall_dir sabs_shadewall_dif => clm3%g%l%lps%sabs_shadewall_dif sabs_improad_dir => clm3%g%l%lps%sabs_improad_dir sabs_improad_dif => clm3%g%l%lps%sabs_improad_dif sabs_perroad_dir => clm3%g%l%lps%sabs_perroad_dir sabs_perroad_dif => clm3%g%l%lps%sabs_perroad_dif ! Assign local pointers to derived type members (column level) ctype => clm3%g%l%c%itype t_grnd => clm3%g%l%c%ces%t_grnd frac_sno => clm3%g%l%c%cps%frac_sno ! Assign local pointers to derived type members (pft level) pgridcell => clm3%g%l%c%p%gridcell pcolumn => clm3%g%l%c%p%column albd => clm3%g%l%c%p%pps%albd albi => clm3%g%l%c%p%pps%albi sabg => clm3%g%l%c%p%pef%sabg sabv => clm3%g%l%c%p%pef%sabv fsa => clm3%g%l%c%p%pef%fsa fsa_u => clm3%g%l%c%p%pef%fsa_u fsr => clm3%g%l%c%p%pef%fsr fsds_vis_d => clm3%g%l%c%p%pef%fsds_vis_d fsds_nir_d => clm3%g%l%c%p%pef%fsds_nir_d fsds_vis_i => clm3%g%l%c%p%pef%fsds_vis_i fsds_nir_i => clm3%g%l%c%p%pef%fsds_nir_i fsr_vis_d => clm3%g%l%c%p%pef%fsr_vis_d fsr_nir_d => clm3%g%l%c%p%pef%fsr_nir_d fsr_vis_i => clm3%g%l%c%p%pef%fsr_vis_i fsr_nir_i => clm3%g%l%c%p%pef%fsr_nir_i fsds_vis_d_ln => clm3%g%l%c%p%pef%fsds_vis_d_ln fsds_nir_d_ln => clm3%g%l%c%p%pef%fsds_nir_d_ln fsr_vis_d_ln => clm3%g%l%c%p%pef%fsr_vis_d_ln fsr_nir_d_ln => clm3%g%l%c%p%pef%fsr_nir_d_ln eflx_lwrad_out => clm3%g%l%c%p%pef%eflx_lwrad_out eflx_lwrad_net => clm3%g%l%c%p%pef%eflx_lwrad_net eflx_lwrad_net_u => clm3%g%l%c%p%pef%eflx_lwrad_net_u parsun => clm3%g%l%c%p%pef%parsun parsha => clm3%g%l%c%p%pef%parsha t_ref2m => clm3%g%l%c%p%pes%t_ref2m ! Define fields that appear on the restart file for non-urban landunits do fl = 1,num_nourbanl l = filter_nourbanl(fl) sabs_roof_dir(l,:) = spval sabs_roof_dif(l,:) = spval sabs_sunwall_dir(l,:) = spval sabs_sunwall_dif(l,:) = spval sabs_shadewall_dir(l,:) = spval sabs_shadewall_dif(l,:) = spval sabs_improad_dir(l,:) = spval sabs_improad_dif(l,:) = spval sabs_perroad_dir(l,:) = spval sabs_perroad_dif(l,:) = spval vf_sr(l) = spval vf_wr(l) = spval vf_sw(l) = spval vf_rw(l) = spval vf_ww(l) = spval end do ! Set input forcing fields do fl = 1,num_urbanl l = filter_urbanl(fl) g = lgridcell(l) ! Need to set the following temperatures to some defined value even if it ! does not appear in the urban landunit for the net_longwave computation t_roof(fl) = 19._r8 + tfrz t_sunwall(fl) = 19._r8 + tfrz t_shadewall(fl) = 19._r8 + tfrz t_improad(fl) = 19._r8 + tfrz t_perroad(fl) = 19._r8 + tfrz ! Initial assignment of emissivity em_roof_s(fl) = em_roof(fl) em_improad_s(fl) = em_improad(fl) em_perroad_s(fl) = em_perroad(fl) ! Set urban temperatures and emissivity including snow effects. do c = coli(l),colf(l) if (ctype(c) == icol_roof ) then t_roof(fl) = t_grnd(c) em_roof_s(fl) = em_roof(fl)*(1._r8-frac_sno(c)) + snoem*frac_sno(c) else if (ctype(c) == icol_road_imperv) then t_improad(fl) = t_grnd(c) em_improad_s(fl) = em_improad(fl)*(1._r8-frac_sno(c)) + snoem*frac_sno(c) else if (ctype(c) == icol_road_perv ) then t_perroad(fl) = t_grnd(c) em_perroad_s(fl) = em_perroad(fl)*(1._r8-frac_sno(c)) + snoem*frac_sno(c) else if (ctype(c) == icol_sunwall ) then t_sunwall(fl) = t_grnd(c) else if (ctype(c) == icol_shadewall ) then t_shadewall(fl) = t_grnd(c) end if end do lwdown(fl) = forc_lwrad(g) end do ! Net longwave radiation for road and both walls in urban canyon allowing for multiple re-emission if (num_urbanl .gt. 0) then call net_longwave (lbl, ubl, num_urbanl, filter_urbanl, canyon_hwr, wtroad_perv, & lwdown, em_roof_s, em_improad_s, em_perroad_s, em_wall, & t_roof, t_improad, t_perroad, t_sunwall, t_shadewall, & lwnet_roof, lwnet_improad, lwnet_perroad, lwnet_sunwall, lwnet_shadewall, lwnet_canyon, & lwup_roof, lwup_improad, lwup_perroad, lwup_sunwall, lwup_shadewall, lwup_canyon) end if dtime = get_step_size() call get_curr_date (year, month, day, secs) ! Determine clmtype variables needed for history output and communication with atm ! Loop over urban pfts in clump do fp = 1,num_urbanp p = filter_urbanp(fp) g = pgridcell(p) local_secp1 = secs + nint((londeg(g)/degpsec)/dtime)*dtime local_secp1 = mod(local_secp1,isecspday) ! Solar incident fsds_vis_d(p) = forc_solad(g,1) fsds_nir_d(p) = forc_solad(g,2) fsds_vis_i(p) = forc_solai(g,1) fsds_nir_i(p) = forc_solai(g,2) ! Determine local noon incident solar if (local_secp1 == noonsec) then fsds_vis_d_ln(p) = forc_solad(g,1) fsds_nir_d_ln(p) = forc_solad(g,2) else fsds_vis_d_ln(p) = spval fsds_nir_d_ln(p) = spval endif ! Solar reflected ! per unit ground area (roof, road) and per unit wall area (sunwall, shadewall) fsr_vis_d(p) = albd(p,1) * forc_solad(g,1) fsr_nir_d(p) = albd(p,2) * forc_solad(g,2) fsr_vis_i(p) = albi(p,1) * forc_solai(g,1) fsr_nir_i(p) = albi(p,2) * forc_solai(g,2) ! Determine local noon reflected solar if (local_secp1 == noonsec) then fsr_vis_d_ln(p) = fsr_vis_d(p) fsr_nir_d_ln(p) = fsr_nir_d(p) else fsr_vis_d_ln(p) = spval fsr_nir_d_ln(p) = spval endif fsr(p) = fsr_vis_d(p) + fsr_nir_d(p) + fsr_vis_i(p) + fsr_nir_i(p) end do ! Loop over urban landunits in clump do fl = 1,num_urbanl l = filter_urbanl(fl) g = lgridcell(l) ! Solar absorbed and longwave out and net ! per unit ground area (roof, road) and per unit wall area (sunwall, shadewall) ! Each urban pft has its own column - this is used in the logic below do p = pfti(l), pftf(l) c = pcolumn(p) if (ctype(c) == icol_roof) then eflx_lwrad_out(p) = lwup_roof(fl) eflx_lwrad_net(p) = lwnet_roof(fl) eflx_lwrad_net_u(p) = lwnet_roof(fl) sabg(p) = sabs_roof_dir(l,1)*forc_solad(g,1) + & sabs_roof_dif(l,1)*forc_solai(g,1) + & sabs_roof_dir(l,2)*forc_solad(g,2) + & sabs_roof_dif(l,2)*forc_solai(g,2) else if (ctype(c) == icol_sunwall) then eflx_lwrad_out(p) = lwup_sunwall(fl) eflx_lwrad_net(p) = lwnet_sunwall(fl) eflx_lwrad_net_u(p) = lwnet_sunwall(fl) sabg(p) = sabs_sunwall_dir(l,1)*forc_solad(g,1) + & sabs_sunwall_dif(l,1)*forc_solai(g,1) + & sabs_sunwall_dir(l,2)*forc_solad(g,2) + & sabs_sunwall_dif(l,2)*forc_solai(g,2) else if (ctype(c) == icol_shadewall) then eflx_lwrad_out(p) = lwup_shadewall(fl) eflx_lwrad_net(p) = lwnet_shadewall(fl) eflx_lwrad_net_u(p) = lwnet_shadewall(fl) sabg(p) = sabs_shadewall_dir(l,1)*forc_solad(g,1) + & sabs_shadewall_dif(l,1)*forc_solai(g,1) + & sabs_shadewall_dir(l,2)*forc_solad(g,2) + & sabs_shadewall_dif(l,2)*forc_solai(g,2) else if (ctype(c) == icol_road_perv) then eflx_lwrad_out(p) = lwup_perroad(fl) eflx_lwrad_net(p) = lwnet_perroad(fl) eflx_lwrad_net_u(p) = lwnet_perroad(fl) sabg(p) = sabs_perroad_dir(l,1)*forc_solad(g,1) + & sabs_perroad_dif(l,1)*forc_solai(g,1) + & sabs_perroad_dir(l,2)*forc_solad(g,2) + & sabs_perroad_dif(l,2)*forc_solai(g,2) else if (ctype(c) == icol_road_imperv) then eflx_lwrad_out(p) = lwup_improad(fl) eflx_lwrad_net(p) = lwnet_improad(fl) eflx_lwrad_net_u(p) = lwnet_improad(fl) sabg(p) = sabs_improad_dir(l,1)*forc_solad(g,1) + & sabs_improad_dif(l,1)*forc_solai(g,1) + & sabs_improad_dir(l,2)*forc_solad(g,2) + & sabs_improad_dif(l,2)*forc_solai(g,2) end if sabv(p) = 0._r8 fsa(p) = sabv(p) + sabg(p) fsa_u(p) = fsa(p) parsun(p) = 0._r8 parsha(p) = 0._r8 end do ! end loop over urban pfts end do ! end loop over urban landunits end subroutine UrbanRadiation !----------------------------------------------------------------------- !BOP ! ! !IROUTINE: view_factor ! ! !INTERFACE: subroutine view_factor (lbl, ubl, num_urbanl, filter_urbanl, canyon_hwr) ! ! !DESCRIPTION: ! View factors for road and one wall ! WALL | ! ROAD | ! wall | ! -----\ /----- - - |\----------/ ! | \ vsr / | | r | | \ vww / s ! | \ / | h o w | \ / k ! wall | \ / | wall | a | | \ / y ! |vwr \ / vwr| | d | |vrw \ / vsw ! ------\/------ - - |-----\/----- ! road wall | ! <----- w ----> | ! <---- h --->| ! ! vsr = view factor of sky for road vrw = view factor of road for wall ! vwr = view factor of one wall for road vww = view factor of opposing wall for wall ! vsw = view factor of sky for wall ! vsr + vwr + vwr = 1 vrw + vww + vsw = 1 ! ! Source: Masson, V. (2000) A physically-based scheme for the urban energy budget in ! atmospheric models. Boundary-Layer Meteorology 94:357-397 ! ! !USES: use shr_kind_mod, only: r8 => shr_kind_r8 use clmtype ! ! !ARGUMENTS: implicit none integer, intent(in) :: lbl, ubl ! landunit-index bounds integer , intent(in) :: num_urbanl ! number of urban landunits integer , intent(in) :: filter_urbanl(ubl-lbl+1) ! urban landunit filter real(r8), intent(in) :: canyon_hwr(num_urbanl) ! ratio of building height to street width ! ! local pointers to original implicit out arguments (clmtype) ! real(r8), pointer :: vf_sr(:) ! view factor of sky for road real(r8), pointer :: vf_wr(:) ! view factor of one wall for road real(r8), pointer :: vf_sw(:) ! view factor of sky for one wall real(r8), pointer :: vf_rw(:) ! view factor of road for one wall real(r8), pointer :: vf_ww(:) ! view factor of opposing wall for one wall ! ! !CALLED FROM: ! subroutine UrbanAlbedo in this module ! ! !REVISION HISTORY: ! Author: Gordon Bonan ! 03/2003, Mariana Vertenstein: Migrated to clm2.2 ! 01/2008, Erik Kluzek: Migrated to clm3.5.15 ! ! ! !LOCAL VARIABLES: !EOP integer :: l, fl ! indices real(r8) :: sum ! sum of view factors for wall or road !----------------------------------------------------------------------- ! Assign landunit level pointer vf_sr => clm3%g%l%lps%vf_sr vf_wr => clm3%g%l%lps%vf_wr vf_sw => clm3%g%l%lps%vf_sw vf_rw => clm3%g%l%lps%vf_rw vf_ww => clm3%g%l%lps%vf_ww do fl = 1,num_urbanl l = filter_urbanl(fl) ! road -- sky view factor -> 1 as building height -> 0 ! and -> 0 as building height -> infinity vf_sr(l) = sqrt(canyon_hwr(fl)**2 + 1._r8) - canyon_hwr(fl) vf_wr(l) = 0.5_r8 * (1._r8 - vf_sr(l)) ! one wall -- sky view factor -> 0.5 as building height -> 0 ! and -> 0 as building height -> infinity vf_sw(l) = 0.5_r8 * (canyon_hwr(fl) + 1._r8 - sqrt(canyon_hwr(fl)**2+1._r8)) / canyon_hwr(fl) vf_rw(l) = vf_sw(l) vf_ww(l) = 1._r8 - vf_sw(l) - vf_rw(l) end do ! error check -- make sure view factor sums to one for road and wall do fl = 1,num_urbanl l = filter_urbanl(fl) sum = vf_sr(l) + 2._r8*vf_wr(l) if (abs(sum-1._r8) > 1.e-06_r8 ) then write (iulog,*) 'urban road view factor error',sum write (iulog,*) 'clm model is stopping' call endrun() endif sum = vf_sw(l) + vf_rw(l) + vf_ww(l) if (abs(sum-1._r8) > 1.e-06_r8 ) then write (iulog,*) 'urban wall view factor error',sum write (iulog,*) 'clm model is stopping' call endrun() endif end do end subroutine view_factor !----------------------------------------------------------------------- !BOP ! ! !IROUTINE: incident_direct ! ! !INTERFACE: subroutine incident_direct (lbl, ubl, num_urbanl, canyon_hwr, coszen, zen, sdir, sdir_road, sdir_sunwall, sdir_shadewall) ! ! !DESCRIPTION: ! Direct beam solar radiation incident on walls and road in urban canyon ! ! Sun ! / ! roof / ! ------ /--- - ! | / | | ! sunlit wall | / | shaded wall h ! | / | | ! -----/----- - ! road ! <--- w ---> ! ! Method: ! Road = Horizontal surface. Account for shading by wall. Integrate over all canyon orientations ! Wall (sunlit) = Adjust horizontal radiation for 90 degree surface. Account for shading by opposing wall. ! Integrate over all canyon orientations ! Wall (shaded) = 0 ! ! Conservation check: Total incoming direct beam (sdir) = sdir_road + (sdir_shadewall + sdir_sunwall)*canyon_hwr ! Multiplication by canyon_hwr scales wall fluxes (per unit wall area) to per unit ground area ! ! Source: Masson, V. (2000) A physically-based scheme for the urban energy budget in ! atmospheric models. Boundary-Layer Meteorology 94:357-397 ! ! This analytical solution from Masson (2000) agrees with the numerical solution to ! within 0.6 W/m**2 for sdir = 1000 W/m**2 and for all H/W from 0.1 to 10 by 0.1 ! and all solar zenith angles from 1 to 90 deg by 1 ! ! !USES: use shr_kind_mod , only : r8 => shr_kind_r8 use clm_varcon , only : rpi implicit none ! ! !ARGUMENTS: integer, intent(in) :: lbl, ubl ! landunit-index bounds integer , intent(in) :: num_urbanl ! number of urban landunits real(r8), intent(in) :: canyon_hwr(num_urbanl) ! ratio of building height to street width real(r8), intent(in) :: coszen(num_urbanl) ! cosine solar zenith angle real(r8), intent(in) :: zen(num_urbanl) ! solar zenith angle (radians) real(r8), intent(in) :: sdir(num_urbanl, numrad) ! direct beam solar radiation incident on horizontal surface real(r8), intent(out) :: sdir_road(num_urbanl, numrad) ! direct beam solar radiation incident on road per unit incident flux real(r8), intent(out) :: sdir_sunwall(num_urbanl, numrad) ! direct beam solar radiation (per unit wall area) incident on sunlit wall per unit incident flux real(r8), intent(out) :: sdir_shadewall(num_urbanl, numrad) ! direct beam solar radiation (per unit wall area) incident on shaded wall per unit incident flux ! ! !CALLED FROM: ! subroutine UrbanAlbedo in this module ! ! !REVISION HISTORY: ! Author: Gordon Bonan ! ! ! !LOCAL VARIABLES: !EOP integer :: l,i,ib ! indices !KO logical :: numchk = .true. ! true => perform numerical check of analytical solution logical :: numchk = .false. ! true => perform numerical check of analytical solution real(r8) :: theta0(num_urbanl) ! critical canyon orientation for which road is no longer illuminated real(r8) :: tanzen(num_urbanl) ! tan(zenith angle) real(r8) :: swall_projected ! direct beam solar radiation (per unit ground area) incident on wall real(r8) :: err1(num_urbanl) ! energy conservation error real(r8) :: err2(num_urbanl) ! energy conservation error real(r8) :: err3(num_urbanl) ! energy conservation error real(r8) :: sumr ! sum of sroad for each orientation (0 <= theta <= pi/2) real(r8) :: sumw ! sum of swall for each orientation (0 <= theta <= pi/2) real(r8) :: num ! number of orientations real(r8) :: theta ! canyon orientation relative to sun (0 <= theta <= pi/2) real(r8) :: zen0 ! critical solar zenith angle for which sun begins to illuminate road !----------------------------------------------------------------------- do l = 1,num_urbanl if (coszen(l) > 0._r8) then theta0(l) = asin(min( (1._r8/(canyon_hwr(l)*tan(max(zen(l),0.000001_r8)))), 1._r8 )) tanzen(l) = tan(zen(l)) end if end do do ib = 1,numrad do l = 1,num_urbanl if (coszen(l) > 0._r8) then sdir_shadewall(l,ib) = 0._r8 ! incident solar radiation on wall and road integrated over all canyon orientations (0 <= theta <= pi/2) sdir_road(l,ib) = sdir(l,ib) * & (2._r8*theta0(l)/rpi - 2./rpi*canyon_hwr(l)*tanzen(l)*(1._r8-cos(theta0(l)))) sdir_sunwall(l,ib) = 2._r8 * sdir(l,ib) * ((1._r8/canyon_hwr(l))* & (0.5_r8-theta0(l)/rpi) + (1._r8/rpi)*tanzen(l)*(1._r8-cos(theta0(l)))) ! conservation check for road and wall. need to use wall fluxes converted to ground area swall_projected = (sdir_shadewall(l,ib) + sdir_sunwall(l,ib)) * canyon_hwr(l) err1(l) = sdir(l,ib) - (sdir_road(l,ib) + swall_projected) else sdir_road(l,ib) = 0._r8 sdir_sunwall(l,ib) = 0._r8 sdir_shadewall(l,ib) = 0._r8 endif end do do l = 1,num_urbanl if (coszen(l) > 0._r8) then if (abs(err1(l)) > 0.001_r8) then write (iulog,*) 'urban direct beam solar radiation balance error',err1(l) write (iulog,*) 'clm model is stopping' call endrun() endif endif end do ! numerical check of analytical solution ! sum sroad and swall over all canyon orientations (0 <= theta <= pi/2) if (numchk) then do l = 1,num_urbanl if (coszen(l) > 0._r8) then sumr = 0._r8 sumw = 0._r8 num = 0._r8 do i = 1, 9000 theta = i/100._r8 * rpi/180._r8 zen0 = atan(1._r8/(canyon_hwr(l)*sin(theta))) if (zen(l) >= zen0) then sumr = sumr + 0._r8 sumw = sumw + sdir(l,ib) / canyon_hwr(l) else sumr = sumr + sdir(l,ib) * (1._r8-canyon_hwr(l)*sin(theta)*tanzen(l)) sumw = sumw + sdir(l,ib) * sin(theta)*tanzen(l) end if num = num + 1._r8 end do err2(l) = sumr/num - sdir_road(l,ib) err3(l) = sumw/num - sdir_sunwall(l,ib) endif end do do l = 1,num_urbanl if (coszen(l) > 0._r8) then if (abs(err2(l)) > 0.0006_r8 ) then write (iulog,*) 'urban road incident direct beam solar radiation error',err2(l) write (iulog,*) 'clm model is stopping' call endrun endif if (abs(err3(l)) > 0.0006_r8 ) then write (iulog,*) 'urban wall incident direct beam solar radiation error',err3(l) write (iulog,*) 'clm model is stopping' call endrun end if end if end do end if end do end subroutine incident_direct !----------------------------------------------------------------------- !BOP ! ! !IROUTINE: incident_diffuse ! ! !INTERFACE: subroutine incident_diffuse (lbl, ubl, num_urbanl, filter_urbanl, canyon_hwr, sdif, sdif_road, sdif_sunwall, sdif_shadewall) ! ! !DESCRIPTION: ! Diffuse solar radiation incident on walls and road in urban canyon ! Conservation check: Total incoming diffuse ! (sdif) = sdif_road + (sdif_shadewall + sdif_sunwall)*canyon_hwr ! Multiplication by canyon_hwr scales wall fluxes (per unit wall area) to per unit ground area ! ! !USES: use shr_kind_mod, only: r8 => shr_kind_r8 use clmtype ! ! ! !ARGUMENTS: implicit none integer, intent(in) :: lbl, ubl ! landunit-index bounds integer , intent(in) :: num_urbanl ! number of urban landunits integer , intent(in) :: filter_urbanl(ubl-lbl+1) ! urban landunit filter real(r8), intent(in) :: canyon_hwr(num_urbanl) ! ratio of building height to street width real(r8), intent(in) :: sdif(num_urbanl, numrad) ! diffuse solar radiation incident on horizontal surface real(r8), intent(out) :: sdif_road(num_urbanl, numrad) ! diffuse solar radiation incident on road real(r8), intent(out) :: sdif_sunwall(num_urbanl, numrad) ! diffuse solar radiation (per unit wall area) incident on sunlit wall real(r8), intent(out) :: sdif_shadewall(num_urbanl, numrad) ! diffuse solar radiation (per unit wall area) incident on shaded wall ! ! local pointers to original implicit in arguments (clmtype) ! real(r8), pointer :: vf_sr(:) ! view factor of sky for road real(r8), pointer :: vf_sw(:) ! view factor of sky for one wall ! ! !CALLED FROM: ! subroutine UrbanAlbedo in this module ! ! !REVISION HISTORY: ! Author: Gordon Bonan ! ! ! !LOCAL VARIABLES: !EOP integer :: l, fl, ib ! indices real(r8) :: err(num_urbanl) ! energy conservation error (W/m**2) real(r8) :: swall_projected ! diffuse solar radiation (per unit ground area) incident on wall (W/m**2) !----------------------------------------------------------------------- ! Assign landunit level pointer vf_sr => clm3%g%l%lps%vf_sr vf_sw => clm3%g%l%lps%vf_sw do ib = 1, numrad ! diffuse solar and conservation check. need to convert wall fluxes to ground area do fl = 1,num_urbanl l = filter_urbanl(fl) sdif_road(fl,ib) = sdif(fl,ib) * vf_sr(l) sdif_sunwall(fl,ib) = sdif(fl,ib) * vf_sw(l) sdif_shadewall(fl,ib) = sdif(fl,ib) * vf_sw(l) swall_projected = (sdif_shadewall(fl,ib) + sdif_sunwall(fl,ib)) * canyon_hwr(fl) err(fl) = sdif(fl,ib) - (sdif_road(fl,ib) + swall_projected) end do ! error check do l = 1, num_urbanl if (abs(err(l)) > 0.001_r8) then write (iulog,*) 'urban diffuse solar radiation balance error',err(l) write (iulog,*) 'clm model is stopping' call endrun endif end do end do end subroutine incident_diffuse !----------------------------------------------------------------------- !BOP ! ! !IROUTINE: net_solar ! ! !INTERFACE: subroutine net_solar (lbl, ubl, num_urbanl, filter_urbanl, coszen, canyon_hwr, wtroad_perv, sdir, sdif, & alb_improad_dir, alb_perroad_dir, alb_wall_dir, alb_roof_dir, & alb_improad_dif, alb_perroad_dif, alb_wall_dif, alb_roof_dif, & sdir_road, sdir_sunwall, sdir_shadewall, & sdif_road, sdif_sunwall, sdif_shadewall, & sref_improad_dir, sref_perroad_dir, sref_sunwall_dir, sref_shadewall_dir, sref_roof_dir, & sref_improad_dif, sref_perroad_dif, sref_sunwall_dif, sref_shadewall_dif, sref_roof_dif) ! ! !DESCRIPTION: ! Solar radiation absorbed by road and both walls in urban canyon allowing ! for multiple reflection. ! ! !USES: use shr_kind_mod, only : r8 => shr_kind_r8 use clmtype ! ! ! !ARGUMENTS: implicit none integer, intent(in) :: lbl, ubl ! landunit-index bounds integer , intent(in) :: num_urbanl ! number of urban landunits integer , intent(in) :: filter_urbanl(ubl-lbl+1) ! urban landunit filter real(r8), intent(in) :: coszen(num_urbanl) ! cosine solar zenith angle real(r8), intent(in) :: canyon_hwr(num_urbanl) ! ratio of building height to street width real(r8), intent(in) :: wtroad_perv(num_urbanl) ! weight of pervious road wrt total road real(r8), intent(in) :: sdir(num_urbanl, numrad) ! direct beam solar radiation incident on horizontal surface real(r8), intent(in) :: sdif(num_urbanl, numrad) ! diffuse solar radiation on horizontal surface real(r8), intent(in) :: alb_improad_dir(num_urbanl, numrad) ! direct impervious road albedo real(r8), intent(in) :: alb_perroad_dir(num_urbanl, numrad) ! direct pervious road albedo real(r8), intent(in) :: alb_wall_dir(num_urbanl, numrad) ! direct wall albedo real(r8), intent(in) :: alb_roof_dir(num_urbanl, numrad) ! direct roof albedo real(r8), intent(in) :: alb_improad_dif(num_urbanl, numrad) ! diffuse impervious road albedo real(r8), intent(in) :: alb_perroad_dif(num_urbanl, numrad) ! diffuse pervious road albedo real(r8), intent(in) :: alb_wall_dif(num_urbanl, numrad) ! diffuse wall albedo real(r8), intent(in) :: alb_roof_dif(num_urbanl, numrad) ! diffuse roof albedo real(r8), intent(in) :: sdir_road(num_urbanl, numrad) ! direct beam solar radiation incident on road per unit incident flux real(r8), intent(in) :: sdir_sunwall(num_urbanl, numrad) ! direct beam solar radiation (per unit wall area) incident on sunlit wall per unit incident flux real(r8), intent(in) :: sdir_shadewall(num_urbanl, numrad) ! direct beam solar radiation (per unit wall area) incident on shaded wall per unit incident flux real(r8), intent(in) :: sdif_road(num_urbanl, numrad) ! diffuse solar radiation incident on road per unit incident flux real(r8), intent(in) :: sdif_sunwall(num_urbanl, numrad) ! diffuse solar radiation (per unit wall area) incident on sunlit wall per unit incident flux real(r8), intent(in) :: sdif_shadewall(num_urbanl, numrad) ! diffuse solar radiation (per unit wall area) incident on shaded wall per unit incident flux real(r8), intent(inout) :: sref_improad_dir(num_urbanl, numrad) ! direct solar rad reflected by impervious road (per unit ground area) per unit incident flux real(r8), intent(inout) :: sref_perroad_dir(num_urbanl, numrad) ! direct solar rad reflected by pervious road (per unit ground area) per unit incident flux real(r8), intent(inout) :: sref_improad_dif(num_urbanl, numrad) ! diffuse solar rad reflected by impervious road (per unit ground area) per unit incident flux real(r8), intent(inout) :: sref_perroad_dif(num_urbanl, numrad) ! diffuse solar rad reflected by pervious road (per unit ground area) per unit incident flux real(r8), intent(inout) :: sref_sunwall_dir(num_urbanl, numrad) ! direct solar rad reflected by sunwall (per unit wall area) per unit incident flux real(r8), intent(inout) :: sref_sunwall_dif(num_urbanl, numrad) ! diffuse solar rad reflected by sunwall (per unit wall area) per unit incident flux real(r8), intent(inout) :: sref_shadewall_dir(num_urbanl, numrad) ! direct solar rad reflected by shadewall (per unit wall area) per unit incident flux real(r8), intent(inout) :: sref_shadewall_dif(num_urbanl, numrad) ! diffuse solar rad reflected by shadewall (per unit wall area) per unit incident flux real(r8), intent(inout) :: sref_roof_dir(num_urbanl, numrad) ! direct solar rad reflected by roof (per unit ground area) per unit incident flux real(r8), intent(inout) :: sref_roof_dif(num_urbanl, numrad) ! diffuse solar rad reflected by roof (per unit ground area) per unit incident flux ! ! local pointers to original implicit in arguments (clmtype) ! real(r8), pointer :: vf_sr(:) ! view factor of sky for road real(r8), pointer :: vf_wr(:) ! view factor of one wall for road real(r8), pointer :: vf_sw(:) ! view factor of sky for one wall real(r8), pointer :: vf_rw(:) ! view factor of road for one wall real(r8), pointer :: vf_ww(:) ! view factor of opposing wall for one wall real(r8), pointer :: sabs_roof_dir(:,:) ! direct solar absorbed by roof per unit ground area per unit incident flux real(r8), pointer :: sabs_roof_dif(:,:) ! diffuse solar absorbed by roof per unit ground area per unit incident flux real(r8), pointer :: sabs_sunwall_dir(:,:) ! direct solar absorbed by sunwall per unit wall area per unit incident flux real(r8), pointer :: sabs_sunwall_dif(:,:) ! diffuse solar absorbed by sunwall per unit wall area per unit incident flux real(r8), pointer :: sabs_shadewall_dir(:,:) ! direct solar absorbed by shadewall per unit wall area per unit incident flux real(r8), pointer :: sabs_shadewall_dif(:,:) ! diffuse solar absorbed by shadewall per unit wall area per unit incident flux real(r8), pointer :: sabs_improad_dir(:,:) ! direct solar absorbed by impervious road per unit ground area per unit incident flux real(r8), pointer :: sabs_improad_dif(:,:) ! diffuse solar absorbed by impervious road per unit ground area per unit incident flux real(r8), pointer :: sabs_perroad_dir(:,:) ! direct solar absorbed by pervious road per unit ground area per unit incident flux real(r8), pointer :: sabs_perroad_dif(:,:) ! diffuse solar absorbed by pervious road per unit ground area per unit incident flux ! ! !CALLED FROM: ! subroutine UrbanAlbedo in this module ! ! !REVISION HISTORY: ! Author: Gordon Bonan ! ! ! !LOCAL VARIABLES !EOP ! real(r8) :: wtroad_imperv(num_urbanl) ! weight of impervious road wrt total road real(r8) :: sabs_canyon_dir(num_urbanl) ! direct solar rad absorbed by canyon per unit incident flux real(r8) :: sabs_canyon_dif(num_urbanl) ! diffuse solar rad absorbed by canyon per unit incident flux real(r8) :: sref_canyon_dir(num_urbanl) ! direct solar reflected by canyon per unit incident flux real(r8) :: sref_canyon_dif(num_urbanl) ! diffuse solar reflected by canyon per unit incident flux real(r8) :: improad_a_dir(num_urbanl) ! absorbed direct solar for impervious road after "n" reflections per unit incident flux real(r8) :: improad_a_dif(num_urbanl) ! absorbed diffuse solar for impervious road after "n" reflections per unit incident flux real(r8) :: improad_r_dir(num_urbanl) ! reflected direct solar for impervious road after "n" reflections per unit incident flux real(r8) :: improad_r_dif(num_urbanl) ! reflected diffuse solar for impervious road after "n" reflections per unit incident flux real(r8) :: improad_r_sky_dir(num_urbanl) ! improad_r_dir to sky per unit incident flux real(r8) :: improad_r_sunwall_dir(num_urbanl) ! improad_r_dir to sunlit wall per unit incident flux real(r8) :: improad_r_shadewall_dir(num_urbanl) ! improad_r_dir to shaded wall per unit incident flux real(r8) :: improad_r_sky_dif(num_urbanl) ! improad_r_dif to sky per unit incident flux real(r8) :: improad_r_sunwall_dif(num_urbanl) ! improad_r_dif to sunlit wall per unit incident flux real(r8) :: improad_r_shadewall_dif(num_urbanl) ! improad_r_dif to shaded wall per unit incident flux real(r8) :: perroad_a_dir(num_urbanl) ! absorbed direct solar for pervious road after "n" reflections per unit incident flux real(r8) :: perroad_a_dif(num_urbanl) ! absorbed diffuse solar for pervious road after "n" reflections per unit incident flux real(r8) :: perroad_r_dir(num_urbanl) ! reflected direct solar for pervious road after "n" reflections per unit incident flux real(r8) :: perroad_r_dif(num_urbanl) ! reflected diffuse solar for pervious road after "n" reflections per unit incident flux real(r8) :: perroad_r_sky_dir(num_urbanl) ! perroad_r_dir to sky per unit incident flux real(r8) :: perroad_r_sunwall_dir(num_urbanl) ! perroad_r_dir to sunlit wall per unit incident flux real(r8) :: perroad_r_shadewall_dir(num_urbanl) ! perroad_r_dir to shaded wall per unit incident flux real(r8) :: perroad_r_sky_dif(num_urbanl) ! perroad_r_dif to sky per unit incident flux real(r8) :: perroad_r_sunwall_dif(num_urbanl) ! perroad_r_dif to sunlit wall per unit incident flux real(r8) :: perroad_r_shadewall_dif(num_urbanl) ! perroad_r_dif to shaded wall per unit incident flux real(r8) :: road_a_dir(num_urbanl) ! absorbed direct solar for total road after "n" reflections per unit incident flux real(r8) :: road_a_dif(num_urbanl) ! absorbed diffuse solar for total road after "n" reflections per unit incident flux real(r8) :: road_r_dir(num_urbanl) ! reflected direct solar for total road after "n" reflections per unit incident flux real(r8) :: road_r_dif(num_urbanl) ! reflected diffuse solar for total road after "n" reflections per unit incident flux real(r8) :: road_r_sky_dir(num_urbanl) ! road_r_dir to sky per unit incident flux real(r8) :: road_r_sunwall_dir(num_urbanl) ! road_r_dir to sunlit wall per unit incident flux real(r8) :: road_r_shadewall_dir(num_urbanl) ! road_r_dir to shaded wall per unit incident flux real(r8) :: road_r_sky_dif(num_urbanl) ! road_r_dif to sky per unit incident flux real(r8) :: road_r_sunwall_dif(num_urbanl) ! road_r_dif to sunlit wall per unit incident flux real(r8) :: road_r_shadewall_dif(num_urbanl) ! road_r_dif to shaded wall per unit incident flux real(r8) :: sunwall_a_dir(num_urbanl) ! absorbed direct solar for sunlit wall (per unit wall area) after "n" reflections per unit incident flux real(r8) :: sunwall_a_dif(num_urbanl) ! absorbed diffuse solar for sunlit wall (per unit wall area) after "n" reflections per unit incident flux real(r8) :: sunwall_r_dir(num_urbanl) ! reflected direct solar for sunlit wall (per unit wall area) after "n" reflections per unit incident flux real(r8) :: sunwall_r_dif(num_urbanl) ! reflected diffuse solar for sunlit wall (per unit wall area) after "n" reflections per unit incident flux real(r8) :: sunwall_r_sky_dir(num_urbanl) ! sunwall_r_dir to sky per unit incident flux real(r8) :: sunwall_r_road_dir(num_urbanl) ! sunwall_r_dir to road per unit incident flux real(r8) :: sunwall_r_shadewall_dir(num_urbanl) ! sunwall_r_dir to opposing (shaded) wall per unit incident flux real(r8) :: sunwall_r_sky_dif(num_urbanl) ! sunwall_r_dif to sky per unit incident flux real(r8) :: sunwall_r_road_dif(num_urbanl) ! sunwall_r_dif to road per unit incident flux real(r8) :: sunwall_r_shadewall_dif(num_urbanl) ! sunwall_r_dif to opposing (shaded) wall per unit incident flux real(r8) :: shadewall_a_dir(num_urbanl) ! absorbed direct solar for shaded wall (per unit wall area) after "n" reflections per unit incident flux real(r8) :: shadewall_a_dif(num_urbanl) ! absorbed diffuse solar for shaded wall (per unit wall area) after "n" reflections per unit incident flux real(r8) :: shadewall_r_dir(num_urbanl) ! reflected direct solar for shaded wall (per unit wall area) after "n" reflections per unit incident flux real(r8) :: shadewall_r_dif(num_urbanl) ! reflected diffuse solar for shaded wall (per unit wall area) after "n" reflections per unit incident flux real(r8) :: shadewall_r_sky_dir(num_urbanl) ! shadewall_r_dir to sky per unit incident flux real(r8) :: shadewall_r_road_dir(num_urbanl) ! shadewall_r_dir to road per unit incident flux real(r8) :: shadewall_r_sunwall_dir(num_urbanl) ! shadewall_r_dir to opposing (sunlit) wall per unit incident flux real(r8) :: shadewall_r_sky_dif(num_urbanl) ! shadewall_r_dif to sky per unit incident flux real(r8) :: shadewall_r_road_dif(num_urbanl) ! shadewall_r_dif to road per unit incident flux real(r8) :: shadewall_r_sunwall_dif(num_urbanl) ! shadewall_r_dif to opposing (sunlit) wall per unit incident flux real(r8) :: canyon_alb_dir(num_urbanl) ! direct canyon albedo real(r8) :: canyon_alb_dif(num_urbanl) ! diffuse canyon albedo real(r8) :: stot(num_urbanl) ! sum of radiative terms real(r8) :: stot_dir(num_urbanl) ! sum of direct radiative terms real(r8) :: stot_dif(num_urbanl) ! sum of diffuse radiative terms integer :: l,fl,ib ! indices integer :: iter_dir,iter_dif ! iteration counter real(r8) :: crit ! convergence criterion real(r8) :: err ! energy conservation error integer :: pass integer, parameter :: n = 50 ! number of interations real(r8) :: sabs_road ! temporary for absorption over road real(r8) :: sref_road ! temporary for reflected over road real(r8), parameter :: errcrit = .00001_r8 ! error criteria !----------------------------------------------------------------------- ! Assign landunit level pointer vf_sr => clm3%g%l%lps%vf_sr vf_wr => clm3%g%l%lps%vf_wr vf_sw => clm3%g%l%lps%vf_sw vf_rw => clm3%g%l%lps%vf_rw vf_ww => clm3%g%l%lps%vf_ww sabs_roof_dir => clm3%g%l%lps%sabs_roof_dir sabs_roof_dif => clm3%g%l%lps%sabs_roof_dif sabs_sunwall_dir => clm3%g%l%lps%sabs_sunwall_dir sabs_sunwall_dif => clm3%g%l%lps%sabs_sunwall_dif sabs_shadewall_dir => clm3%g%l%lps%sabs_shadewall_dir sabs_shadewall_dif => clm3%g%l%lps%sabs_shadewall_dif sabs_improad_dir => clm3%g%l%lps%sabs_improad_dir sabs_improad_dif => clm3%g%l%lps%sabs_improad_dif sabs_perroad_dir => clm3%g%l%lps%sabs_perroad_dir sabs_perroad_dif => clm3%g%l%lps%sabs_perroad_dif ! Calculate impervious road do l = 1,num_urbanl wtroad_imperv(l) = 1._r8 - wtroad_perv(l) end do do ib = 1,numrad do fl = 1,num_urbanl if (coszen(fl) .gt. 0._r8) then l = filter_urbanl(fl) ! initial absorption and reflection for road and both walls. ! distribute reflected radiation to sky, road, and walls ! according to appropriate view factor. radiation reflected to ! road and walls will undergo multiple reflections within the canyon. ! do separately for direct beam and diffuse radiation. ! direct beam road_a_dir(fl) = 0.0_r8 road_r_dir(fl) = 0.0_r8 if ( wtroad_imperv(fl) > 0.0_r8 ) then improad_a_dir(fl) = (1._r8-alb_improad_dir(fl,ib)) * sdir_road(fl,ib) improad_r_dir(fl) = alb_improad_dir(fl,ib) * sdir_road(fl,ib) improad_r_sky_dir(fl) = improad_r_dir(fl) * vf_sr(l) improad_r_sunwall_dir(fl) = improad_r_dir(fl) * vf_wr(l) improad_r_shadewall_dir(fl) = improad_r_dir(fl) * vf_wr(l) road_a_dir(fl) = road_a_dir(fl) + improad_a_dir(fl)*wtroad_imperv(fl) road_r_dir(fl) = road_r_dir(fl) + improad_r_dir(fl)*wtroad_imperv(fl) end if if ( wtroad_perv(fl) > 0.0_r8 ) then perroad_a_dir(fl) = (1._r8-alb_perroad_dir(fl,ib)) * sdir_road(fl,ib) perroad_r_dir(fl) = alb_perroad_dir(fl,ib) * sdir_road(fl,ib) perroad_r_sky_dir(fl) = perroad_r_dir(fl) * vf_sr(l) perroad_r_sunwall_dir(fl) = perroad_r_dir(fl) * vf_wr(l) perroad_r_shadewall_dir(fl) = perroad_r_dir(fl) * vf_wr(l) road_a_dir(fl) = road_a_dir(fl) + perroad_a_dir(fl)*wtroad_perv(fl) road_r_dir(fl) = road_r_dir(fl) + perroad_r_dir(fl)*wtroad_perv(fl) end if road_r_sky_dir(fl) = road_r_dir(fl) * vf_sr(l) road_r_sunwall_dir(fl) = road_r_dir(fl) * vf_wr(l) road_r_shadewall_dir(fl) = road_r_dir(fl) * vf_wr(l) sunwall_a_dir(fl) = (1._r8-alb_wall_dir(fl,ib)) * sdir_sunwall(fl,ib) sunwall_r_dir(fl) = alb_wall_dir(fl,ib) * sdir_sunwall(fl,ib) sunwall_r_sky_dir(fl) = sunwall_r_dir(fl) * vf_sw(l) sunwall_r_road_dir(fl) = sunwall_r_dir(fl) * vf_rw(l) sunwall_r_shadewall_dir(fl) = sunwall_r_dir(fl) * vf_ww(l) shadewall_a_dir(fl) = (1._r8-alb_wall_dir(fl,ib)) * sdir_shadewall(fl,ib) shadewall_r_dir(fl) = alb_wall_dir(fl,ib) * sdir_shadewall(fl,ib) shadewall_r_sky_dir(fl) = shadewall_r_dir(fl) * vf_sw(l) shadewall_r_road_dir(fl) = shadewall_r_dir(fl) * vf_rw(l) shadewall_r_sunwall_dir(fl) = shadewall_r_dir(fl) * vf_ww(l) ! diffuse road_a_dif(fl) = 0.0_r8 road_r_dif(fl) = 0.0_r8 if ( wtroad_imperv(fl) > 0.0_r8 ) then improad_a_dif(fl) = (1._r8-alb_improad_dif(fl,ib)) * sdif_road(fl,ib) improad_r_dif(fl) = alb_improad_dif(fl,ib) * sdif_road(fl,ib) improad_r_sky_dif(fl) = improad_r_dif(fl) * vf_sr(l) improad_r_sunwall_dif(fl) = improad_r_dif(fl) * vf_wr(l) improad_r_shadewall_dif(fl) = improad_r_dif(fl) * vf_wr(l) road_a_dif(fl) = road_a_dif(fl) + improad_a_dif(fl)*wtroad_imperv(fl) road_r_dif(fl) = road_r_dif(fl) + improad_r_dif(fl)*wtroad_imperv(fl) end if if ( wtroad_perv(fl) > 0.0_r8 ) then perroad_a_dif(fl) = (1._r8-alb_perroad_dif(fl,ib)) * sdif_road(fl,ib) perroad_r_dif(fl) = alb_perroad_dif(fl,ib) * sdif_road(fl,ib) perroad_r_sky_dif(fl) = perroad_r_dif(fl) * vf_sr(l) perroad_r_sunwall_dif(fl) = perroad_r_dif(fl) * vf_wr(l) perroad_r_shadewall_dif(fl) = perroad_r_dif(fl) * vf_wr(l) road_a_dif(fl) = road_a_dif(fl) + perroad_a_dif(fl)*wtroad_perv(fl) road_r_dif(fl) = road_r_dif(fl) + perroad_r_dif(fl)*wtroad_perv(fl) end if road_r_sky_dif(fl) = road_r_dif(fl) * vf_sr(l) road_r_sunwall_dif(fl) = road_r_dif(fl) * vf_wr(l) road_r_shadewall_dif(fl) = road_r_dif(fl) * vf_wr(l) sunwall_a_dif(fl) = (1._r8-alb_wall_dif(fl,ib)) * sdif_sunwall(fl,ib) sunwall_r_dif(fl) = alb_wall_dif(fl,ib) * sdif_sunwall(fl,ib) sunwall_r_sky_dif(fl) = sunwall_r_dif(fl) * vf_sw(l) sunwall_r_road_dif(fl) = sunwall_r_dif(fl) * vf_rw(l) sunwall_r_shadewall_dif(fl) = sunwall_r_dif(fl) * vf_ww(l) shadewall_a_dif(fl) = (1._r8-alb_wall_dif(fl,ib)) * sdif_shadewall(fl,ib) shadewall_r_dif(fl) = alb_wall_dif(fl,ib) * sdif_shadewall(fl,ib) shadewall_r_sky_dif(fl) = shadewall_r_dif(fl) * vf_sw(l) shadewall_r_road_dif(fl) = shadewall_r_dif(fl) * vf_rw(l) shadewall_r_sunwall_dif(fl) = shadewall_r_dif(fl) * vf_ww(l) ! initialize sum of direct and diffuse solar absorption and reflection for road and both walls if ( wtroad_imperv(fl) > 0.0_r8 ) sabs_improad_dir(l,ib) = improad_a_dir(fl) if ( wtroad_perv(fl) > 0.0_r8 ) sabs_perroad_dir(l,ib) = perroad_a_dir(fl) sabs_sunwall_dir(l,ib) = sunwall_a_dir(fl) sabs_shadewall_dir(l,ib) = shadewall_a_dir(fl) if ( wtroad_imperv(fl) > 0.0_r8 ) sabs_improad_dif(l,ib) = improad_a_dif(fl) if ( wtroad_perv(fl) > 0.0_r8 ) sabs_perroad_dif(l,ib) = perroad_a_dif(fl) sabs_sunwall_dif(l,ib) = sunwall_a_dif(fl) sabs_shadewall_dif(l,ib) = shadewall_a_dif(fl) if ( wtroad_imperv(fl) > 0.0_r8 ) sref_improad_dir(fl,ib) = improad_r_sky_dir(fl) if ( wtroad_perv(fl) > 0.0_r8 ) sref_perroad_dir(fl,ib) = perroad_r_sky_dir(fl) sref_sunwall_dir(fl,ib) = sunwall_r_sky_dir(fl) sref_shadewall_dir(fl,ib) = shadewall_r_sky_dir(fl) if ( wtroad_imperv(fl) > 0.0_r8 ) sref_improad_dif(fl,ib) = improad_r_sky_dif(fl) if ( wtroad_perv(fl) > 0.0_r8 ) sref_perroad_dif(fl,ib) = perroad_r_sky_dif(fl) sref_sunwall_dif(fl,ib) = sunwall_r_sky_dif(fl) sref_shadewall_dif(fl,ib) = shadewall_r_sky_dif(fl) endif end do ! absorption and reflection for walls and road with multiple reflections ! (i.e., absorb and reflect initial reflection in canyon and allow for ! subsequent scattering) ! ! (1) absorption and reflection of scattered solar radiation ! road: reflected fluxes from walls need to be projected to ground area ! wall: reflected flux from road needs to be projected to wall area ! ! (2) add absorbed radiation for ith reflection to total absorbed ! ! (3) distribute reflected radiation to sky, road, and walls according to view factors ! ! (4) add solar reflection to sky for ith reflection to total reflection ! ! (5) stop iteration when absorption for ith reflection is less than some nominal amount. ! small convergence criteria is required to ensure solar radiation is conserved ! ! do separately for direct beam and diffuse do fl = 1,num_urbanl if (coszen(fl) .gt. 0._r8) then l = filter_urbanl(fl) ! reflected direct beam do iter_dir = 1, n ! step (1) stot(fl) = (sunwall_r_road_dir(fl) + shadewall_r_road_dir(fl))*canyon_hwr(fl) road_a_dir(fl) = 0.0_r8 road_r_dir(fl) = 0.0_r8 if ( wtroad_imperv(fl) > 0.0_r8 ) then improad_a_dir(fl) = (1._r8-alb_improad_dir(fl,ib)) * stot(fl) improad_r_dir(fl) = alb_improad_dir(fl,ib) * stot(fl) road_a_dir(fl) = road_a_dir(fl) + improad_a_dir(fl)*wtroad_imperv(fl) road_r_dir(fl) = road_r_dir(fl) + improad_r_dir(fl)*wtroad_imperv(fl) end if if ( wtroad_perv(fl) > 0.0_r8 ) then perroad_a_dir(fl) = (1._r8-alb_perroad_dir(fl,ib)) * stot(fl) perroad_r_dir(fl) = alb_perroad_dir(fl,ib) * stot(fl) road_a_dir(fl) = road_a_dir(fl) + perroad_a_dir(fl)*wtroad_perv(fl) road_r_dir(fl) = road_r_dir(fl) + perroad_r_dir(fl)*wtroad_perv(fl) end if stot(fl) = road_r_sunwall_dir(fl)/canyon_hwr(fl) + shadewall_r_sunwall_dir(fl) sunwall_a_dir(fl) = (1._r8-alb_wall_dir(fl,ib)) * stot(fl) sunwall_r_dir(fl) = alb_wall_dir(fl,ib) * stot(fl) stot(fl) = road_r_shadewall_dir(fl)/canyon_hwr(fl) + sunwall_r_shadewall_dir(fl) shadewall_a_dir(fl) = (1._r8-alb_wall_dir(fl,ib)) * stot(fl) shadewall_r_dir(fl) = alb_wall_dir(fl,ib) * stot(fl) ! step (2) if ( wtroad_imperv(fl) > 0.0_r8 ) sabs_improad_dir(l,ib) = sabs_improad_dir(l,ib) + improad_a_dir(fl) if ( wtroad_perv(fl) > 0.0_r8 ) sabs_perroad_dir(l,ib) = sabs_perroad_dir(l,ib) + perroad_a_dir(fl) sabs_sunwall_dir(l,ib) = sabs_sunwall_dir(l,ib) + sunwall_a_dir(fl) sabs_shadewall_dir(l,ib) = sabs_shadewall_dir(l,ib) + shadewall_a_dir(fl) ! step (3) if ( wtroad_imperv(fl) > 0.0_r8 ) then improad_r_sky_dir(fl) = improad_r_dir(fl) * vf_sr(l) improad_r_sunwall_dir(fl) = improad_r_dir(fl) * vf_wr(l) improad_r_shadewall_dir(fl) = improad_r_dir(fl) * vf_wr(l) end if if ( wtroad_perv(fl) > 0.0_r8 ) then perroad_r_sky_dir(fl) = perroad_r_dir(fl) * vf_sr(l) perroad_r_sunwall_dir(fl) = perroad_r_dir(fl) * vf_wr(l) perroad_r_shadewall_dir(fl) = perroad_r_dir(fl) * vf_wr(l) end if road_r_sky_dir(fl) = road_r_dir(fl) * vf_sr(l) road_r_sunwall_dir(fl) = road_r_dir(fl) * vf_wr(l) road_r_shadewall_dir(fl) = road_r_dir(fl) * vf_wr(l) sunwall_r_sky_dir(fl) = sunwall_r_dir(fl) * vf_sw(l) sunwall_r_road_dir(fl) = sunwall_r_dir(fl) * vf_rw(l) sunwall_r_shadewall_dir(fl) = sunwall_r_dir(fl) * vf_ww(l) shadewall_r_sky_dir(fl) = shadewall_r_dir(fl) * vf_sw(l) shadewall_r_road_dir(fl) = shadewall_r_dir(fl) * vf_rw(l) shadewall_r_sunwall_dir(fl) = shadewall_r_dir(fl) * vf_ww(l) ! step (4) if ( wtroad_imperv(fl) > 0.0_r8 ) sref_improad_dir(fl,ib) = sref_improad_dir(fl,ib) + improad_r_sky_dir(fl) if ( wtroad_perv(fl) > 0.0_r8 ) sref_perroad_dir(fl,ib) = sref_perroad_dir(fl,ib) + perroad_r_sky_dir(fl) sref_sunwall_dir(fl,ib) = sref_sunwall_dir(fl,ib) + sunwall_r_sky_dir(fl) sref_shadewall_dir(fl,ib) = sref_shadewall_dir(fl,ib) + shadewall_r_sky_dir(fl) ! step (5) crit = max(road_a_dir(fl), sunwall_a_dir(fl), shadewall_a_dir(fl)) if (crit < errcrit) exit end do if (iter_dir >= n) then write (iulog,*) 'urban net solar radiation error: no convergence, direct beam' write (iulog,*) 'clm model is stopping' call endrun endif ! reflected diffuse do iter_dif = 1, n ! step (1) stot(fl) = (sunwall_r_road_dif(fl) + shadewall_r_road_dif(fl))*canyon_hwr(fl) road_a_dif(fl) = 0.0_r8 road_r_dif(fl) = 0.0_r8 if ( wtroad_imperv(fl) > 0.0_r8 ) then improad_a_dif(fl) = (1._r8-alb_improad_dif(fl,ib)) * stot(fl) improad_r_dif(fl) = alb_improad_dif(fl,ib) * stot(fl) road_a_dif(fl) = road_a_dif(fl) + improad_a_dif(fl)*wtroad_imperv(fl) road_r_dif(fl) = road_r_dif(fl) + improad_r_dif(fl)*wtroad_imperv(fl) end if if ( wtroad_perv(fl) > 0.0_r8 ) then perroad_a_dif(fl) = (1._r8-alb_perroad_dif(fl,ib)) * stot(fl) perroad_r_dif(fl) = alb_perroad_dif(fl,ib) * stot(fl) road_a_dif(fl) = road_a_dif(fl) + perroad_a_dif(fl)*wtroad_perv(fl) road_r_dif(fl) = road_r_dif(fl) + perroad_r_dif(fl)*wtroad_perv(fl) end if stot(fl) = road_r_sunwall_dif(fl)/canyon_hwr(fl) + shadewall_r_sunwall_dif(fl) sunwall_a_dif(fl) = (1._r8-alb_wall_dif(fl,ib)) * stot(fl) sunwall_r_dif(fl) = alb_wall_dif(fl,ib) * stot(fl) stot(fl) = road_r_shadewall_dif(fl)/canyon_hwr(fl) + sunwall_r_shadewall_dif(fl) shadewall_a_dif(fl) = (1._r8-alb_wall_dif(fl,ib)) * stot(fl) shadewall_r_dif(fl) = alb_wall_dif(fl,ib) * stot(fl) ! step (2) if ( wtroad_imperv(fl) > 0.0_r8 ) sabs_improad_dif(l,ib) = sabs_improad_dif(l,ib) + improad_a_dif(fl) if ( wtroad_perv(fl) > 0.0_r8 ) sabs_perroad_dif(l,ib) = sabs_perroad_dif(l,ib) + perroad_a_dif(fl) sabs_sunwall_dif(l,ib) = sabs_sunwall_dif(l,ib) + sunwall_a_dif(fl) sabs_shadewall_dif(l,ib) = sabs_shadewall_dif(l,ib) + shadewall_a_dif(fl) ! step (3) if ( wtroad_imperv(fl) > 0.0_r8 ) then improad_r_sky_dif(fl) = improad_r_dif(fl) * vf_sr(l) improad_r_sunwall_dif(fl) = improad_r_dif(fl) * vf_wr(l) improad_r_shadewall_dif(fl) = improad_r_dif(fl) * vf_wr(l) end if if ( wtroad_perv(fl) > 0.0_r8 ) then perroad_r_sky_dif(fl) = perroad_r_dif(fl) * vf_sr(l) perroad_r_sunwall_dif(fl) = perroad_r_dif(fl) * vf_wr(l) perroad_r_shadewall_dif(fl) = perroad_r_dif(fl) * vf_wr(l) end if road_r_sky_dif(fl) = road_r_dif(fl) * vf_sr(l) road_r_sunwall_dif(fl) = road_r_dif(fl) * vf_wr(l) road_r_shadewall_dif(fl) = road_r_dif(fl) * vf_wr(l) sunwall_r_sky_dif(fl) = sunwall_r_dif(fl) * vf_sw(l) sunwall_r_road_dif(fl) = sunwall_r_dif(fl) * vf_rw(l) sunwall_r_shadewall_dif(fl) = sunwall_r_dif(fl) * vf_ww(l) shadewall_r_sky_dif(fl) = shadewall_r_dif(fl) * vf_sw(l) shadewall_r_road_dif(fl) = shadewall_r_dif(fl) * vf_rw(l) shadewall_r_sunwall_dif(fl) = shadewall_r_dif(fl) * vf_ww(l) ! step (4) if ( wtroad_imperv(fl) > 0.0_r8 ) sref_improad_dif(fl,ib) = sref_improad_dif(fl,ib) + improad_r_sky_dif(fl) if ( wtroad_perv(fl) > 0.0_r8 ) sref_perroad_dif(fl,ib) = sref_perroad_dif(fl,ib) + perroad_r_sky_dif(fl) sref_sunwall_dif(fl,ib) = sref_sunwall_dif(fl,ib) + sunwall_r_sky_dif(fl) sref_shadewall_dif(fl,ib) = sref_shadewall_dif(fl,ib) + shadewall_r_sky_dif(fl) ! step (5) crit = max(road_a_dif(fl), sunwall_a_dif(fl), shadewall_a_dif(fl)) if (crit < errcrit) exit end do if (iter_dif >= n) then write (iulog,*) 'urban net solar radiation error: no convergence, diffuse' write (iulog,*) 'clm model is stopping' call endrun() endif ! total reflected by canyon - sum of solar reflection to sky from canyon. ! project wall fluxes to horizontal surface sref_canyon_dir(fl) = 0.0_r8 sref_canyon_dif(fl) = 0.0_r8 if ( wtroad_imperv(fl) > 0.0_r8 ) then sref_canyon_dir(fl) = sref_canyon_dir(fl) + sref_improad_dir(fl,ib)*wtroad_imperv(fl) sref_canyon_dif(fl) = sref_canyon_dif(fl) + sref_improad_dif(fl,ib)*wtroad_imperv(fl) end if if ( wtroad_perv(fl) > 0.0_r8 ) then sref_canyon_dir(fl) = sref_canyon_dir(fl) + sref_perroad_dir(fl,ib)*wtroad_perv(fl) sref_canyon_dif(fl) = sref_canyon_dif(fl) + sref_perroad_dif(fl,ib)*wtroad_perv(fl) end if sref_canyon_dir(fl) = sref_canyon_dir(fl) + (sref_sunwall_dir(fl,ib) + sref_shadewall_dir(fl,ib))*canyon_hwr(fl) sref_canyon_dif(fl) = sref_canyon_dif(fl) + (sref_sunwall_dif(fl,ib) + sref_shadewall_dif(fl,ib))*canyon_hwr(fl) ! total absorbed by canyon. project wall fluxes to horizontal surface sabs_canyon_dir(fl) = 0.0_r8 sabs_canyon_dif(fl) = 0.0_r8 if ( wtroad_imperv(fl) > 0.0_r8 ) then sabs_canyon_dir(fl) = sabs_canyon_dir(fl) + sabs_improad_dir(l,ib)*wtroad_imperv(fl) sabs_canyon_dif(fl) = sabs_canyon_dif(fl) + sabs_improad_dif(l,ib)*wtroad_imperv(fl) end if if ( wtroad_perv(fl) > 0.0_r8 ) then sabs_canyon_dir(fl) = sabs_canyon_dir(fl) + sabs_perroad_dir(l,ib)*wtroad_perv(fl) sabs_canyon_dif(fl) = sabs_canyon_dif(fl) + sabs_perroad_dif(l,ib)*wtroad_perv(fl) end if sabs_canyon_dir(fl) = sabs_canyon_dir(fl) + (sabs_sunwall_dir(l,ib) + sabs_shadewall_dir(l,ib))*canyon_hwr(fl) sabs_canyon_dif(fl) = sabs_canyon_dif(fl) + (sabs_sunwall_dif(l,ib) + sabs_shadewall_dif(l,ib))*canyon_hwr(fl) ! conservation check. note: previous conservation checks confirm partioning of total direct ! beam and diffuse radiation from atmosphere to road and walls is conserved as ! sdir (from atmosphere) = sdir_road + (sdir_sunwall + sdir_shadewall)*canyon_hwr ! sdif (from atmosphere) = sdif_road + (sdif_sunwall + sdif_shadewall)*canyon_hwr stot_dir(fl) = sdir_road(fl,ib) + (sdir_sunwall(fl,ib) + sdir_shadewall(fl,ib))*canyon_hwr(fl) stot_dif(fl) = sdif_road(fl,ib) + (sdif_sunwall(fl,ib) + sdif_shadewall(fl,ib))*canyon_hwr(fl) err = stot_dir(fl) + stot_dif(fl) & - (sabs_canyon_dir(fl) + sabs_canyon_dif(fl) + sref_canyon_dir(fl) + sref_canyon_dif(fl)) if (abs(err) > 0.001_r8 ) then write(iulog,*)'urban net solar radiation balance error for ib=',ib,' err= ',err write(iulog,*)' l= ',l,' ib= ',ib write(iulog,*)' stot_dir = ',stot_dir(fl) write(iulog,*)' stot_dif = ',stot_dif(fl) write(iulog,*)' sabs_canyon_dir = ',sabs_canyon_dir(fl) write(iulog,*)' sabs_canyon_dif = ',sabs_canyon_dif(fl) write(iulog,*)' sref_canyon_dir = ',sref_canyon_dir(fl) write(iulog,*)' sref_canyon_dif = ',sref_canyon_dir(fl) write(iulog,*) 'clm model is stopping' call endrun() endif ! canyon albedo canyon_alb_dif(fl) = sref_canyon_dif(fl) / max(stot_dif(fl), 1.e-06_r8) canyon_alb_dir(fl) = sref_canyon_dir(fl) / max(stot_dir(fl), 1.e-06_r8) end if end do ! end of landunit loop ! Refected and absorbed solar radiation per unit incident radiation for roof do fl = 1,num_urbanl if (coszen(fl) .gt. 0._r8) then l = filter_urbanl(fl) sref_roof_dir(fl,ib) = alb_roof_dir(fl,ib) * sdir(fl,ib) sref_roof_dif(fl,ib) = alb_roof_dif(fl,ib) * sdif(fl,ib) sabs_roof_dir(l,ib) = sdir(fl,ib) - sref_roof_dir(fl,ib) sabs_roof_dif(l,ib) = sdif(fl,ib) - sref_roof_dif(fl,ib) end if end do end do ! end of radiation band loop end subroutine net_solar !----------------------------------------------------------------------- !BOP ! ! !IROUTINE: net_longwave ! ! !INTERFACE: subroutine net_longwave (lbl, ubl, num_urbanl, filter_urbanl, canyon_hwr, wtroad_perv, & lwdown, em_roof, em_improad, em_perroad, em_wall, & t_roof, t_improad, t_perroad, t_sunwall, t_shadewall, & lwnet_roof, lwnet_improad, lwnet_perroad, lwnet_sunwall, lwnet_shadewall, lwnet_canyon, & lwup_roof, lwup_improad, lwup_perroad, lwup_sunwall, lwup_shadewall, lwup_canyon) ! ! !DESCRIPTION: ! Net longwave radiation for road and both walls in urban canyon allowing for ! multiple reflection. Also net longwave radiation for urban roof. ! ! !USES: use shr_kind_mod , only : r8 => shr_kind_r8 use clm_varcon , only : sb use clmtype ! ! !ARGUMENTS: implicit none integer , intent(in) :: num_urbanl ! number of urban landunits integer, intent(in) :: lbl, ubl ! landunit-index bounds integer , intent(in) :: filter_urbanl(ubl-lbl+1) ! urban landunit filter real(r8), intent(in) :: canyon_hwr(num_urbanl) ! ratio of building height to street width real(r8), intent(in) :: wtroad_perv(num_urbanl) ! weight of pervious road wrt total road real(r8), intent(in) :: lwdown(num_urbanl) ! atmospheric longwave radiation (W/m**2) real(r8), intent(in) :: em_roof(num_urbanl) ! roof emissivity real(r8), intent(in) :: em_improad(num_urbanl) ! impervious road emissivity real(r8), intent(in) :: em_perroad(num_urbanl) ! pervious road emissivity real(r8), intent(in) :: em_wall(num_urbanl) ! wall emissivity real(r8), intent(in) :: t_roof(num_urbanl) ! roof temperature (K) real(r8), intent(in) :: t_improad(num_urbanl) ! impervious road temperature (K) real(r8), intent(in) :: t_perroad(num_urbanl) ! ervious road temperature (K) real(r8), intent(in) :: t_sunwall(num_urbanl) ! sunlit wall temperature (K) real(r8), intent(in) :: t_shadewall(num_urbanl) ! shaded wall temperature (K) real(r8), intent(out) :: lwnet_roof(num_urbanl) ! net (outgoing-incoming) longwave radiation, roof (W/m**2) real(r8), intent(out) :: lwnet_improad(num_urbanl) ! net (outgoing-incoming) longwave radiation, impervious road (W/m**2) real(r8), intent(out) :: lwnet_perroad(num_urbanl) ! net (outgoing-incoming) longwave radiation, pervious road (W/m**2) real(r8), intent(out) :: lwnet_sunwall(num_urbanl) ! net (outgoing-incoming) longwave radiation (per unit wall area), sunlit wall (W/m**2) real(r8), intent(out) :: lwnet_shadewall(num_urbanl) ! net (outgoing-incoming) longwave radiation (per unit wall area), shaded wall (W/m**2) real(r8), intent(out) :: lwnet_canyon(num_urbanl) ! net (outgoing-incoming) longwave radiation for canyon, per unit ground area (W/m**2) real(r8), intent(out) :: lwup_roof(num_urbanl) ! upward longwave radiation, roof (W/m**2) real(r8), intent(out) :: lwup_improad(num_urbanl) ! upward longwave radiation, impervious road (W/m**2) real(r8), intent(out) :: lwup_perroad(num_urbanl) ! upward longwave radiation, pervious road (W/m**2) real(r8), intent(out) :: lwup_sunwall(num_urbanl) ! upward longwave radiation (per unit wall area), sunlit wall (W/m**2) real(r8), intent(out) :: lwup_shadewall(num_urbanl) ! upward longwave radiation (per unit wall area), shaded wall (W/m**2) real(r8), intent(out) :: lwup_canyon(num_urbanl) ! upward longwave radiation for canyon, per unit ground area (W/m**2) ! ! local pointers to original implicit in arguments (clmtype) ! real(r8), pointer :: vf_sr(:) ! view factor of sky for road real(r8), pointer :: vf_wr(:) ! view factor of one wall for road real(r8), pointer :: vf_sw(:) ! view factor of sky for one wall real(r8), pointer :: vf_rw(:) ! view factor of road for one wall real(r8), pointer :: vf_ww(:) ! view factor of opposing wall for one wall ! ! !CALLED FROM: ! subroutine UrbanRadiation in this module ! ! !REVISION HISTORY: ! Author: Gordon Bonan ! ! ! !LOCAL VARIABLES: !EOP real(r8) :: lwdown_road(num_urbanl) ! atmospheric longwave radiation for total road (W/m**2) real(r8) :: lwdown_sunwall(num_urbanl) ! atmospheric longwave radiation (per unit wall area) for sunlit wall (W/m**2) real(r8) :: lwdown_shadewall(num_urbanl) ! atmospheric longwave radiation (per unit wall area) for shaded wall (W/m**2) real(r8) :: lwtot(num_urbanl) ! incoming longwave radiation (W/m**2) real(r8) :: improad_a(num_urbanl) ! absorbed longwave for improad (W/m**2) real(r8) :: improad_r(num_urbanl) ! reflected longwave for improad (W/m**2) real(r8) :: improad_r_sky(num_urbanl) ! improad_r to sky (W/m**2) real(r8) :: improad_r_sunwall(num_urbanl) ! improad_r to sunlit wall (W/m**2) real(r8) :: improad_r_shadewall(num_urbanl) ! improad_r to shaded wall (W/m**2) real(r8) :: improad_e(num_urbanl) ! emitted longwave for improad (W/m**2) real(r8) :: improad_e_sky(num_urbanl) ! improad_e to sky (W/m**2) real(r8) :: improad_e_sunwall(num_urbanl) ! improad_e to sunlit wall (W/m**2) real(r8) :: improad_e_shadewall(num_urbanl) ! improad_e to shaded wall (W/m**2) real(r8) :: perroad_a(num_urbanl) ! absorbed longwave for perroad (W/m**2) real(r8) :: perroad_r(num_urbanl) ! reflected longwave for perroad (W/m**2) real(r8) :: perroad_r_sky(num_urbanl) ! perroad_r to sky (W/m**2) real(r8) :: perroad_r_sunwall(num_urbanl) ! perroad_r to sunlit wall (W/m**2) real(r8) :: perroad_r_shadewall(num_urbanl) ! perroad_r to shaded wall (W/m**2) real(r8) :: perroad_e(num_urbanl) ! emitted longwave for perroad (W/m**2) real(r8) :: perroad_e_sky(num_urbanl) ! perroad_e to sky (W/m**2) real(r8) :: perroad_e_sunwall(num_urbanl) ! perroad_e to sunlit wall (W/m**2) real(r8) :: perroad_e_shadewall(num_urbanl) ! perroad_e to shaded wall (W/m**2) real(r8) :: road_a(num_urbanl) ! absorbed longwave for total road (W/m**2) real(r8) :: road_r(num_urbanl) ! reflected longwave for total road (W/m**2) real(r8) :: road_r_sky(num_urbanl) ! total road_r to sky (W/m**2) real(r8) :: road_r_sunwall(num_urbanl) ! total road_r to sunlit wall (W/m**2) real(r8) :: road_r_shadewall(num_urbanl) ! total road_r to shaded wall (W/m**2) real(r8) :: road_e(num_urbanl) ! emitted longwave for total road (W/m**2) real(r8) :: road_e_sky(num_urbanl) ! total road_e to sky (W/m**2) real(r8) :: road_e_sunwall(num_urbanl) ! total road_e to sunlit wall (W/m**2) real(r8) :: road_e_shadewall(num_urbanl) ! total road_e to shaded wall (W/m**2) real(r8) :: sunwall_a(num_urbanl) ! absorbed longwave (per unit wall area) for sunlit wall (W/m**2) real(r8) :: sunwall_r(num_urbanl) ! reflected longwave (per unit wall area) for sunlit wall (W/m**2) real(r8) :: sunwall_r_sky(num_urbanl) ! sunwall_r to sky (W/m**2) real(r8) :: sunwall_r_road(num_urbanl) ! sunwall_r to road (W/m**2) real(r8) :: sunwall_r_shadewall(num_urbanl) ! sunwall_r to opposing (shaded) wall (W/m**2) real(r8) :: sunwall_e(num_urbanl) ! emitted longwave (per unit wall area) for sunlit wall (W/m**2) real(r8) :: sunwall_e_sky(num_urbanl) ! sunwall_e to sky (W/m**2) real(r8) :: sunwall_e_road(num_urbanl) ! sunwall_e to road (W/m**2) real(r8) :: sunwall_e_shadewall(num_urbanl) ! sunwall_e to opposing (shaded) wall (W/m**2) real(r8) :: shadewall_a(num_urbanl) ! absorbed longwave (per unit wall area) for shaded wall (W/m**2) real(r8) :: shadewall_r(num_urbanl) ! reflected longwave (per unit wall area) for shaded wall (W/m**2) real(r8) :: shadewall_r_sky(num_urbanl) ! shadewall_r to sky (W/m**2) real(r8) :: shadewall_r_road(num_urbanl) ! shadewall_r to road (W/m**2) real(r8) :: shadewall_r_sunwall(num_urbanl) ! shadewall_r to opposing (sunlit) wall (W/m**2) real(r8) :: shadewall_e(num_urbanl) ! emitted longwave (per unit wall area) for shaded wall (W/m**2) real(r8) :: shadewall_e_sky(num_urbanl) ! shadewall_e to sky (W/m**2) real(r8) :: shadewall_e_road(num_urbanl) ! shadewall_e to road (W/m**2) real(r8) :: shadewall_e_sunwall(num_urbanl) ! shadewall_e to opposing (sunlit) wall (W/m**2) integer :: l,fl,iter ! indices integer, parameter :: n = 50 ! number of interations real(r8) :: crit ! convergence criterion (W/m**2) real(r8) :: err ! energy conservation error (W/m**2) real(r8) :: wtroad_imperv(num_urbanl) ! weight of impervious road wrt total road !----------------------------------------------------------------------- ! Assign landunit level pointer vf_sr => clm3%g%l%lps%vf_sr vf_wr => clm3%g%l%lps%vf_wr vf_sw => clm3%g%l%lps%vf_sw vf_rw => clm3%g%l%lps%vf_rw vf_ww => clm3%g%l%lps%vf_ww ! Calculate impervious road do l = 1,num_urbanl wtroad_imperv(l) = 1._r8 - wtroad_perv(l) end do do fl = 1,num_urbanl l = filter_urbanl(fl) ! atmospheric longwave radiation incident on walls and road in urban canyon. ! check for conservation (need to convert wall fluxes to ground area). ! lwdown (from atmosphere) = lwdown_road + (lwdown_sunwall + lwdown_shadewall)*canyon_hwr lwdown_road(fl) = lwdown(fl) * vf_sr(l) lwdown_sunwall(fl) = lwdown(fl) * vf_sw(l) lwdown_shadewall(fl) = lwdown(fl) * vf_sw(l) err = lwdown(fl) - (lwdown_road(fl) + (lwdown_shadewall(fl) + lwdown_sunwall(fl))*canyon_hwr(fl)) if (abs(err) > 0.10_r8 ) then write (iulog,*) 'urban incident atmospheric longwave radiation balance error',err write (iulog,*) 'clm model is stopping' call endrun endif end do do fl = 1,num_urbanl l = filter_urbanl(fl) ! initial absorption, reflection, and emission for road and both walls. ! distribute reflected and emitted radiation to sky, road, and walls according ! to appropriate view factor. radiation reflected to road and walls will ! undergo multiple reflections within the canyon. road_a(fl) = 0.0_r8 road_r(fl) = 0.0_r8 road_e(fl) = 0.0_r8 if ( wtroad_imperv(fl) > 0.0_r8 ) then improad_a(fl) = em_improad(fl) * lwdown_road(fl) improad_r(fl) = (1._r8-em_improad(fl)) * lwdown_road(fl) improad_r_sky(fl) = improad_r(fl) * vf_sr(l) improad_r_sunwall(fl) = improad_r(fl) * vf_wr(l) improad_r_shadewall(fl) = improad_r(fl) * vf_wr(l) improad_e(fl) = em_improad(fl) * sb * (t_improad(fl)**4) improad_e_sky(fl) = improad_e(fl) * vf_sr(l) improad_e_sunwall(fl) = improad_e(fl) * vf_wr(l) improad_e_shadewall(fl) = improad_e(fl) * vf_wr(l) road_a(fl) = road_a(fl) + improad_a(fl)*wtroad_imperv(fl) road_r(fl) = road_r(fl) + improad_r(fl)*wtroad_imperv(fl) road_e(fl) = road_e(fl) + improad_e(fl)*wtroad_imperv(fl) end if if ( wtroad_perv(fl) > 0.0_r8 ) then perroad_a(fl) = em_perroad(fl) * lwdown_road(fl) perroad_r(fl) = (1._r8-em_perroad(fl)) * lwdown_road(fl) perroad_r_sky(fl) = perroad_r(fl) * vf_sr(l) perroad_r_sunwall(fl) = perroad_r(fl) * vf_wr(l) perroad_r_shadewall(fl) = perroad_r(fl) * vf_wr(l) perroad_e(fl) = em_perroad(fl) * sb * (t_perroad(fl)**4) perroad_e_sky(fl) = perroad_e(fl) * vf_sr(l) perroad_e_sunwall(fl) = perroad_e(fl) * vf_wr(l) perroad_e_shadewall(fl) = perroad_e(fl) * vf_wr(l) road_a(fl) = road_a(fl) + perroad_a(fl)*wtroad_perv(fl) road_r(fl) = road_r(fl) + perroad_r(fl)*wtroad_perv(fl) road_e(fl) = road_e(fl) + perroad_e(fl)*wtroad_perv(fl) end if road_r_sky(fl) = road_r(fl) * vf_sr(l) road_r_sunwall(fl) = road_r(fl) * vf_wr(l) road_r_shadewall(fl) = road_r(fl) * vf_wr(l) road_e_sky(fl) = road_e(fl) * vf_sr(l) road_e_sunwall(fl) = road_e(fl) * vf_wr(l) road_e_shadewall(fl) = road_e(fl) * vf_wr(l) sunwall_a(fl) = em_wall(fl) * lwdown_sunwall(fl) sunwall_r(fl) = (1._r8-em_wall(fl)) * lwdown_sunwall(fl) sunwall_r_sky(fl) = sunwall_r(fl) * vf_sw(l) sunwall_r_road(fl) = sunwall_r(fl) * vf_rw(l) sunwall_r_shadewall(fl) = sunwall_r(fl) * vf_ww(l) sunwall_e(fl) = em_wall(fl) * sb * (t_sunwall(fl)**4) sunwall_e_sky(fl) = sunwall_e(fl) * vf_sw(l) sunwall_e_road(fl) = sunwall_e(fl) * vf_rw(l) sunwall_e_shadewall(fl) = sunwall_e(fl) * vf_ww(l) shadewall_a(fl) = em_wall(fl) * lwdown_shadewall(fl) shadewall_r(fl) = (1._r8-em_wall(fl)) * lwdown_shadewall(fl) shadewall_r_sky(fl) = shadewall_r(fl) * vf_sw(l) shadewall_r_road(fl) = shadewall_r(fl) * vf_rw(l) shadewall_r_sunwall(fl) = shadewall_r(fl) * vf_ww(l) shadewall_e(fl) = em_wall(fl) * sb * (t_shadewall(fl)**4) shadewall_e_sky(fl) = shadewall_e(fl) * vf_sw(l) shadewall_e_road(fl) = shadewall_e(fl) * vf_rw(l) shadewall_e_sunwall(fl) = shadewall_e(fl) * vf_ww(l) ! initialize sum of net and upward longwave radiation for road and both walls if ( wtroad_imperv(fl) > 0.0_r8 ) lwnet_improad(fl) = improad_e(fl) - improad_a(fl) if ( wtroad_perv(fl) > 0.0_r8 ) lwnet_perroad(fl) = perroad_e(fl) - perroad_a(fl) lwnet_sunwall(fl) = sunwall_e(fl) - sunwall_a(fl) lwnet_shadewall(fl) = shadewall_e(fl) - shadewall_a(fl) if ( wtroad_imperv(fl) > 0.0_r8 ) lwup_improad(fl) = improad_r_sky(fl) + improad_e_sky(fl) if ( wtroad_perv(fl) > 0.0_r8 ) lwup_perroad(fl) = perroad_r_sky(fl) + perroad_e_sky(fl) lwup_sunwall(fl) = sunwall_r_sky(fl) + sunwall_e_sky(fl) lwup_shadewall(fl) = shadewall_r_sky(fl) + shadewall_e_sky(fl) end do ! now account for absorption and reflection within canyon of fluxes from road and walls ! allowing for multiple reflections ! ! (1) absorption and reflection. note: emission from road and walls absorbed by walls and roads ! only occurs in first iteration. zero out for later iterations. ! ! road: fluxes from walls need to be projected to ground area ! wall: fluxes from road need to be projected to wall area ! ! (2) add net longwave for ith reflection to total net longwave ! ! (3) distribute reflected radiation to sky, road, and walls according to view factors ! ! (4) add upward longwave radiation to sky from road and walls for ith reflection to total ! ! (5) stop iteration when absorption for ith reflection is less than some nominal amount. ! small convergence criteria is required to ensure radiation is conserved do fl = 1,num_urbanl l = filter_urbanl(fl) do iter = 1, n ! step (1) lwtot(fl) = (sunwall_r_road(fl) + sunwall_e_road(fl) & + shadewall_r_road(fl) + shadewall_e_road(fl))*canyon_hwr(fl) road_a(fl) = 0.0_r8 road_r(fl) = 0.0_r8 if ( wtroad_imperv(fl) > 0.0_r8 ) then improad_r(fl) = (1._r8-em_improad(fl)) * lwtot(fl) improad_a(fl) = em_improad(fl) * lwtot(fl) road_a(fl) = road_a(fl) + improad_a(fl)*wtroad_imperv(fl) road_r(fl) = road_r(fl) + improad_r(fl)*wtroad_imperv(fl) end if if ( wtroad_perv(fl) > 0.0_r8 ) then perroad_r(fl) = (1._r8-em_perroad(fl)) * lwtot(fl) perroad_a(fl) = em_perroad(fl) * lwtot(fl) road_a(fl) = road_a(fl) + perroad_a(fl)*wtroad_perv(fl) road_r(fl) = road_r(fl) + perroad_r(fl)*wtroad_perv(fl) end if lwtot(fl) = (road_r_sunwall(fl) + road_e_sunwall(fl))/canyon_hwr(fl) & + (shadewall_r_sunwall(fl) + shadewall_e_sunwall(fl)) sunwall_a(fl) = em_wall(fl) * lwtot(fl) sunwall_r(fl) = (1._r8-em_wall(fl)) * lwtot(fl) lwtot(fl) = (road_r_shadewall(fl) + road_e_shadewall(fl))/canyon_hwr(fl) & + (sunwall_r_shadewall(fl) + sunwall_e_shadewall(fl)) shadewall_a(fl) = em_wall(fl) * lwtot(fl) shadewall_r(fl) = (1._r8-em_wall(fl)) * lwtot(fl) sunwall_e_road(fl) = 0._r8 shadewall_e_road(fl) = 0._r8 road_e_sunwall(fl) = 0._r8 shadewall_e_sunwall(fl) = 0._r8 road_e_shadewall(fl) = 0._r8 sunwall_e_shadewall(fl) = 0._r8 ! step (2) if ( wtroad_imperv(fl) > 0.0_r8 ) lwnet_improad(fl) = lwnet_improad(fl) - improad_a(fl) if ( wtroad_perv(fl) > 0.0_r8 ) lwnet_perroad(fl) = lwnet_perroad(fl) - perroad_a(fl) lwnet_sunwall(fl) = lwnet_sunwall(fl) - sunwall_a(fl) lwnet_shadewall(fl) = lwnet_shadewall(fl) - shadewall_a(fl) ! step (3) if ( wtroad_imperv(fl) > 0.0_r8 ) then improad_r_sky(fl) = improad_r(fl) * vf_sr(l) improad_r_sunwall(fl) = improad_r(fl) * vf_wr(l) improad_r_shadewall(fl) = improad_r(fl) * vf_wr(l) end if if ( wtroad_perv(fl) > 0.0_r8 ) then perroad_r_sky(fl) = perroad_r(fl) * vf_sr(l) perroad_r_sunwall(fl) = perroad_r(fl) * vf_wr(l) perroad_r_shadewall(fl) = perroad_r(fl) * vf_wr(l) end if road_r_sky(fl) = road_r(fl) * vf_sr(l) road_r_sunwall(fl) = road_r(fl) * vf_wr(l) road_r_shadewall(fl) = road_r(fl) * vf_wr(l) sunwall_r_sky(fl) = sunwall_r(fl) * vf_sw(l) sunwall_r_road(fl) = sunwall_r(fl) * vf_rw(l) sunwall_r_shadewall(fl) = sunwall_r(fl) * vf_ww(l) shadewall_r_sky(fl) = shadewall_r(fl) * vf_sw(l) shadewall_r_road(fl) = shadewall_r(fl) * vf_rw(l) shadewall_r_sunwall(fl) = shadewall_r(fl) * vf_ww(l) ! step (4) if ( wtroad_imperv(fl) > 0.0_r8 ) lwup_improad(fl) = lwup_improad(fl) + improad_r_sky(fl) if ( wtroad_perv(fl) > 0.0_r8 ) lwup_perroad(fl) = lwup_perroad(fl) + perroad_r_sky(fl) lwup_sunwall(fl) = lwup_sunwall(fl) + sunwall_r_sky(fl) lwup_shadewall(fl) = lwup_shadewall(fl) + shadewall_r_sky(fl) ! step (5) crit = max(road_a(fl), sunwall_a(fl), shadewall_a(fl)) if (crit < .001_r8) exit end do if (iter >= n) then write (iulog,*) 'urban net longwave radiation error: no convergence' write (iulog,*) 'clm model is stopping' call endrun endif ! total net longwave radiation for canyon. project wall fluxes to horizontal surface lwnet_canyon(fl) = 0.0_r8 if ( wtroad_imperv(fl) > 0.0_r8 ) lwnet_canyon(fl) = lwnet_canyon(fl) + lwnet_improad(fl)*wtroad_imperv(fl) if ( wtroad_perv(fl) > 0.0_r8 ) lwnet_canyon(fl) = lwnet_canyon(fl) + lwnet_perroad(fl)*wtroad_perv(fl) lwnet_canyon(fl) = lwnet_canyon(fl) + (lwnet_sunwall(fl) + lwnet_shadewall(fl))*canyon_hwr(fl) ! total emitted longwave for canyon. project wall fluxes to horizontal lwup_canyon(fl) = 0.0_r8 if( wtroad_imperv(fl) > 0.0_r8 ) lwup_canyon(fl) = lwup_canyon(fl) + lwup_improad(fl)*wtroad_imperv(fl) if( wtroad_perv(fl) > 0.0_r8 ) lwup_canyon(fl) = lwup_canyon(fl) + lwup_perroad(fl)*wtroad_perv(fl) lwup_canyon(fl) = lwup_canyon(fl) + (lwup_sunwall(fl) + lwup_shadewall(fl))*canyon_hwr(fl) ! conservation check. note: previous conservation check confirms partioning of incident ! atmospheric longwave radiation to road and walls is conserved as ! lwdown (from atmosphere) = lwdown_improad + lwdown_perroad + (lwdown_sunwall + lwdown_shadewall)*canyon_hwr err = lwnet_canyon(fl) - (lwup_canyon(fl) - lwdown(fl)) if (abs(err) > .10_r8 ) then write (iulog,*) 'urban net longwave radiation balance error',err write (iulog,*) 'clm model is stopping' call endrun() end if end do ! Net longwave radiation for roof do l = 1,num_urbanl lwup_roof(l) = em_roof(l)*sb*(t_roof(l)**4) + (1._r8-em_roof(l))*lwdown(l) lwnet_roof(l) = lwup_roof(l) - lwdown(l) end do end subroutine net_longwave !----------------------------------------------------------------------- !BOP ! ! !IROUTINE: UrbanClumpInit ! ! !INTERFACE: subroutine UrbanClumpInit() ! ! !DESCRIPTION: ! Initialize urban radiation module ! ! !USES: use clmtype use clm_varcon , only : spval, icol_roof, icol_sunwall, icol_shadewall, & icol_road_perv, icol_road_imperv use decompMod , only : get_proc_clumps, ldecomp use filterMod , only : filter use UrbanInputMod, only : urbinp ! ! !ARGUMENTS: implicit none ! ! !CALLED FROM: ! subroutine initialize ! ! !REVISION HISTORY: ! Author: Mariana Vertenstein 04/2003 ! ! !LOCAL VARIABLES: ! ! local pointers to original implicit in arguments ! integer , pointer :: coli(:) ! beginning column index for landunit integer , pointer :: colf(:) ! ending column index for landunit integer , pointer :: lgridcell(:) ! gridcell of corresponding landunit integer , pointer :: ctype(:) ! column type ! ! ! !OTHER LOCAL VARIABLES !EOP ! integer :: nc,fl,ib,l,c,p,g ! indices integer :: nclumps ! number of clumps on processor integer :: num_urbanl ! number of per-clump urban landunits integer :: ier ! error status !----------------------------------------------------------------------- ! Assign local pointers to derived type members (landunit-level) coli => clm3%g%l%coli colf => clm3%g%l%colf lgridcell => clm3%g%l%gridcell ! Assign local pointers to derived type members (column-level) ctype => clm3%g%l%c%itype ! Allocate memory nclumps = get_proc_clumps() allocate(urban_clump(nclumps), stat=ier) if (ier /= 0) then write (iulog,*) 'UrbanInit: allocation error for urban clumps'; call endrun() end if ! Loop over all clumps on this processor do nc = 1, nclumps ! Determine number of unrban landunits in clump num_urbanl = filter(nc)%num_urbanl ! Consistency check for urban columns do fl = 1,num_urbanl l = filter(nc)%urbanl(fl) do c = coli(l),colf(l) if ( ctype(c) /= icol_roof .and. & ctype(c) /= icol_sunwall .and. ctype(c) /= icol_shadewall .and. & ctype(c) /= icol_road_perv .and. ctype(c) /= icol_road_imperv) then write(iulog,*)'error in urban column types for landunit = ',l write(iulog,*)'ctype= ',ctype(c) call endrun() endif end do end do ! Allocate memory for urban clump clumponents if (num_urbanl > 0) then allocate( urban_clump(nc)%canyon_hwr (num_urbanl), & urban_clump(nc)%wtroad_perv (num_urbanl), & urban_clump(nc)%ht_roof (num_urbanl), & urban_clump(nc)%wtlunit_roof (num_urbanl), & urban_clump(nc)%wind_hgt_canyon (num_urbanl), & urban_clump(nc)%em_roof (num_urbanl), & urban_clump(nc)%em_improad (num_urbanl), & urban_clump(nc)%em_perroad (num_urbanl), & urban_clump(nc)%em_wall (num_urbanl), & urban_clump(nc)%alb_roof_dir (num_urbanl,numrad), & urban_clump(nc)%alb_roof_dif (num_urbanl,numrad), & urban_clump(nc)%alb_improad_dir (num_urbanl,numrad), & urban_clump(nc)%alb_perroad_dir (num_urbanl,numrad), & urban_clump(nc)%alb_improad_dif (num_urbanl,numrad), & urban_clump(nc)%alb_perroad_dif (num_urbanl,numrad), & urban_clump(nc)%alb_wall_dir (num_urbanl,numrad), & urban_clump(nc)%alb_wall_dif (num_urbanl,numrad), stat=ier ) if (ier /= 0) then write(iulog,*)'UrbanRadInit: allocation error for urban derived type'; call endrun() endif end if ! Set constants in derived type values for urban clump do fl = 1,num_urbanl l = filter(nc)%urbanl(fl) g = clm3%g%l%gridcell(l) urban_clump(nc)%canyon_hwr (fl) = urbinp%canyon_hwr (g) urban_clump(nc)%wtroad_perv (fl) = urbinp%wtroad_perv (g) urban_clump(nc)%ht_roof (fl) = urbinp%ht_roof (g) urban_clump(nc)%wtlunit_roof (fl) = urbinp%wtlunit_roof (g) urban_clump(nc)%wind_hgt_canyon(fl) = urbinp%wind_hgt_canyon(g) do ib = 1,numrad urban_clump(nc)%alb_roof_dir (fl,ib) = urbinp%alb_roof_dir (g,ib) urban_clump(nc)%alb_roof_dif (fl,ib) = urbinp%alb_roof_dif (g,ib) urban_clump(nc)%alb_improad_dir(fl,ib) = urbinp%alb_improad_dir(g,ib) urban_clump(nc)%alb_perroad_dir(fl,ib) = urbinp%alb_perroad_dir(g,ib) urban_clump(nc)%alb_improad_dif(fl,ib) = urbinp%alb_improad_dif(g,ib) urban_clump(nc)%alb_perroad_dif(fl,ib) = urbinp%alb_perroad_dif(g,ib) urban_clump(nc)%alb_wall_dir (fl,ib) = urbinp%alb_wall_dir (g,ib) urban_clump(nc)%alb_wall_dif (fl,ib) = urbinp%alb_wall_dif (g,ib) end do urban_clump(nc)%em_roof (fl) = urbinp%em_roof (g) urban_clump(nc)%em_improad(fl) = urbinp%em_improad(g) urban_clump(nc)%em_perroad(fl) = urbinp%em_perroad(g) urban_clump(nc)%em_wall (fl) = urbinp%em_wall (g) ! write(iulog,*)'g: ',g ! write(iulog,*)'l: ',l ! write(iulog,*)'fl: ',fl ! write(iulog,*)'urban_clump(nc)%canyon_hwr: ',urban_clump(nc)%canyon_hwr(fl) ! write(iulog,*)'urban_clump(nc)%wtroad_perv: ',urban_clump(nc)%wtroad_perv(fl) ! write(iulog,*)'urban_clump(nc)%ht_roof: ',urban_clump(nc)%ht_roof(fl) ! write(iulog,*)'urban_clump(nc)%wtlunit_roof: ',urban_clump(nc)%wtlunit_roof(fl) ! write(iulog,*)'urban_clump(nc)%wind_hgt_canyon: ',urban_clump(nc)%wind_hgt_canyon(fl) ! write(iulog,*)'urban_clump(nc)%alb_roof_dir: ',urban_clump(nc)%alb_roof_dir(fl,:) ! write(iulog,*)'urban_clump(nc)%alb_roof_dif: ',urban_clump(nc)%alb_roof_dif(fl,:) ! write(iulog,*)'urban_clump(nc)%alb_improad_dir: ',urban_clump(nc)%alb_improad_dir(fl,:) ! write(iulog,*)'urban_clump(nc)%alb_improad_dif: ',urban_clump(nc)%alb_improad_dif(fl,:) ! write(iulog,*)'urban_clump(nc)%alb_perroad_dir: ',urban_clump(nc)%alb_perroad_dir(fl,:) ! write(iulog,*)'urban_clump(nc)%alb_perroad_dif: ',urban_clump(nc)%alb_perroad_dif(fl,:) ! write(iulog,*)'urban_clump(nc)%alb_wall_dir: ',urban_clump(nc)%alb_wall_dir(fl,:) ! write(iulog,*)'urban_clump(nc)%alb_wall_dif: ',urban_clump(nc)%alb_wall_dif(fl,:) ! write(iulog,*)'urban_clump(nc)%em_roof: ',urban_clump(nc)%em_roof(fl) ! write(iulog,*)'urban_clump(nc)%em_improad: ',urban_clump(nc)%em_improad(fl) ! write(iulog,*)'urban_clump(nc)%em_perroad: ',urban_clump(nc)%em_perroad(fl) ! write(iulog,*)'urban_clump(nc)%em_wall: ',urban_clump(nc)%em_wall(fl) end do end do ! end of loop over clumps end subroutine UrbanClumpInit !----------------------------------------------------------------------- !BOP ! ! !IROUTINE: UrbanFluxes ! ! !INTERFACE: subroutine UrbanFluxes (nc, lbp, ubp, lbl, ubl, lbc, ubc, & num_nourbanl, filter_nourbanl, & num_urbanl, filter_urbanl, & num_urbanc, filter_urbanc, & num_urbanp, filter_urbanp) ! ! !DESCRIPTION: ! Turbulent and momentum fluxes from urban canyon (consisting of roof, sunwall, ! shadewall, pervious and impervious road). ! !USES: use clmtype use clm_varcon , only : cpair, vkc, spval, icol_roof, icol_sunwall, & icol_shadewall, icol_road_perv, icol_road_imperv, & grav, pondmx_urban, rpi, rgas, & ht_wasteheat_factor, ac_wasteheat_factor, & wasteheat_limit use filterMod , only : filter use FrictionVelocityMod, only : FrictionVelocity, MoninObukIni use QSatMod , only : QSat use clm_varpar , only : maxpatch_urb, nlevurb use clm_time_manager , only : get_curr_date, get_step_size, get_nstep use clm_atmlnd , only : clm_a2l ! ! !ARGUMENTS: implicit none integer , intent(in) :: nc ! clump index integer, intent(in) :: lbp, ubp ! pft-index bounds integer, intent(in) :: lbl, ubl ! landunit-index bounds integer, intent(in) :: lbc, ubc ! column-index bounds integer , intent(in) :: num_nourbanl ! number of non-urban landunits in clump integer , intent(in) :: filter_nourbanl(ubl-lbl+1) ! non-urban landunit filter integer , intent(in) :: num_urbanl ! number of urban landunits in clump integer , intent(in) :: filter_urbanl(ubl-lbl+1) ! urban landunit filter integer , intent(in) :: num_urbanc ! number of urban columns in clump integer , intent(in) :: filter_urbanc(ubc-lbc+1) ! urban column filter integer , intent(in) :: num_urbanp ! number of urban pfts in clump integer , intent(in) :: filter_urbanp(ubp-lbp+1) ! urban pft filter ! ! !CALLED FROM: ! subroutine clm_driver1 ! ! !REVISION HISTORY: ! Author: Keith Oleson 10/2005 ! ! !LOCAL VARIABLES: ! ! local pointers to original implicit in arguments (urban clump) ! real(r8), pointer :: ht_roof(:) ! height of urban roof (m) real(r8), pointer :: wtlunit_roof(:) ! weight of roof with respect to landunit real(r8), pointer :: canyon_hwr(:) ! ratio of building height to street width real(r8), pointer :: wtroad_perv(:) ! weight of pervious road wrt total road real(r8), pointer :: wind_hgt_canyon(:) ! height above road at which wind in canyon is to be computed (m) ! ! local pointers to original implicit in arguments (clmtype) ! real(r8), pointer :: forc_u(:) ! atmospheric wind speed in east direction (m/s) real(r8), pointer :: forc_v(:) ! atmospheric wind speed in north direction (m/s) real(r8), pointer :: forc_rho(:) ! density (kg/m**3) real(r8), pointer :: forc_hgt_u_pft(:) ! observational height of wind at pft-level (m) real(r8), pointer :: forc_hgt_t_pft(:) ! observational height of temperature at pft-level (m) real(r8), pointer :: forc_q(:) ! atmospheric specific humidity (kg/kg) real(r8), pointer :: forc_t(:) ! atmospheric temperature (K) real(r8), pointer :: forc_th(:) ! atmospheric potential temperature (K) real(r8), pointer :: forc_pbot(:) ! atmospheric pressure (Pa) real(r8), pointer :: z_0_town(:) ! momentum roughness length of urban landunit (m) real(r8), pointer :: z_d_town(:) ! displacement height of urban landunit (m) integer , pointer :: pgridcell(:) ! gridcell of corresponding pft integer , pointer :: pcolumn(:) ! column of corresponding pft integer , pointer :: lgridcell(:) ! gridcell of corresponding landunit integer , pointer :: plandunit(:) ! pft's landunit index integer , pointer :: ctype(:) ! column type integer , pointer :: coli(:) ! beginning column index for landunit integer , pointer :: colf(:) ! ending column index for landunit integer , pointer :: pfti(:) ! beginning pft index for landunit integer , pointer :: pftf(:) ! ending pft index for landunit real(r8), pointer :: taf(:) ! urban canopy air temperature (K) real(r8), pointer :: qaf(:) ! urban canopy air specific humidity (kg/kg) integer , pointer :: npfts(:) ! landunit's number of pfts (columns) real(r8), pointer :: t_grnd(:) ! ground surface temperature (K) real(r8), pointer :: qg(:) ! specific humidity at ground surface (kg/kg) real(r8), pointer :: htvp(:) ! latent heat of evaporation (/sublimation) (J/kg) real(r8), pointer :: dqgdT(:) ! temperature derivative of "qg" real(r8), pointer :: eflx_traffic(:) ! traffic sensible heat flux (W/m**2) real(r8), pointer :: eflx_traffic_factor(:) ! multiplicative urban traffic factor for sensible heat flux real(r8), pointer :: eflx_wasteheat(:) ! sensible heat flux from urban heating/cooling sources of waste heat (W/m**2) real(r8), pointer :: eflx_heat_from_ac(:) ! sensible heat flux put back into canyon due to removal by AC (W/m**2) real(r8), pointer :: t_soisno(:,:) ! soil temperature (K) real(r8), pointer :: eflx_urban_ac(:) ! urban air conditioning flux (W/m**2) real(r8), pointer :: eflx_urban_heat(:) ! urban heating flux (W/m**2) real(r8), pointer :: londeg(:) ! longitude (degrees) real(r8), pointer :: h2osoi_ice(:,:) ! ice lens (kg/m2) real(r8), pointer :: h2osoi_liq(:,:) ! liquid water (kg/m2) real(r8), pointer :: frac_sno(:) ! fraction of ground covered by snow (0 to 1) real(r8), pointer :: snowdp(:) ! snow height (m) real(r8), pointer :: h2osno(:) ! snow water (mm H2O) integer , pointer :: snl(:) ! number of snow layers real(r8), pointer :: rootr_road_perv(:,:) ! effective fraction of roots in each soil layer for urban pervious road real(r8), pointer :: soilalpha_u(:) ! Urban factor that reduces ground saturated specific humidity (-) ! ! local pointers to original implicit out arguments ! real(r8), pointer :: dlrad(:) ! downward longwave radiation below the canopy (W/m**2) real(r8), pointer :: ulrad(:) ! upward longwave radiation above the canopy (W/m**2) real(r8), pointer :: cgrnds(:) ! deriv, of soil sensible heat flux wrt soil temp (W/m**2/K) real(r8), pointer :: cgrndl(:) ! deriv of soil latent heat flux wrt soil temp (W/m**2/K) real(r8), pointer :: cgrnd(:) ! deriv. of soil energy flux wrt to soil temp (W/m**2/K) real(r8), pointer :: taux(:) ! wind (shear) stress: e-w (kg/m/s**2) real(r8), pointer :: tauy(:) ! wind (shear) stress: n-s (kg/m/s**2) real(r8), pointer :: eflx_sh_grnd(:) ! sensible heat flux from ground (W/m**2) [+ to atm] real(r8), pointer :: eflx_sh_tot(:) ! total sensible heat flux (W/m**2) [+ to atm] real(r8), pointer :: eflx_sh_tot_u(:) ! urban total sensible heat flux (W/m**2) [+ to atm] real(r8), pointer :: qflx_evap_soi(:) ! soil evaporation (mm H2O/s) (+ = to atm) real(r8), pointer :: qflx_tran_veg(:) ! vegetation transpiration (mm H2O/s) (+ = to atm) real(r8), pointer :: qflx_evap_veg(:) ! vegetation evaporation (mm H2O/s) (+ = to atm) real(r8), pointer :: qflx_evap_tot(:) ! qflx_evap_soi + qflx_evap_can + qflx_tran_veg real(r8), pointer :: t_ref2m(:) ! 2 m height surface air temperature (K) real(r8), pointer :: q_ref2m(:) ! 2 m height surface specific humidity (kg/kg) real(r8), pointer :: t_ref2m_u(:) ! Urban 2 m height surface air temperature (K) real(r8), pointer :: t_veg(:) ! vegetation temperature (K) real(r8), pointer :: ram1(:) ! aerodynamical resistance (s/m) real(r8), pointer :: rootr(:,:) ! effective fraction of roots in each soil layer real(r8), pointer :: psnsun(:) ! sunlit leaf photosynthesis (umol CO2 /m**2/ s) real(r8), pointer :: psnsha(:) ! shaded leaf photosynthesis (umol CO2 /m**2/ s) real(r8), pointer :: t_building(:) ! internal building temperature (K) real(r8), pointer :: rh_ref2m(:) ! 2 m height surface relative humidity (%) real(r8), pointer :: rh_ref2m_u(:) ! Urban 2 m height surface relative humidity (%) ! ! ! !OTHER LOCAL VARIABLES !EOP ! character(len=*), parameter :: sub="UrbanFluxes" integer :: fp,fc,fl,f,p,c,l,g,j,pi,i ! indices real(r8) :: canyontop_wind(num_urbanl) ! wind at canyon top (m/s) real(r8) :: canyon_u_wind(num_urbanl) ! u-component of wind speed inside canyon (m/s) real(r8) :: canyon_wind(num_urbanl) ! net wind speed inside canyon (m/s) real(r8) :: canyon_resistance(num_urbanl) ! resistance to heat and moisture transfer from canyon road/walls to canyon air (s/m) real(r8) :: ur(lbl:ubl) ! wind speed at reference height (m/s) real(r8) :: ustar(lbl:ubl) ! friction velocity (m/s) real(r8) :: ramu(lbl:ubl) ! aerodynamic resistance (s/m) real(r8) :: rahu(lbl:ubl) ! thermal resistance (s/m) real(r8) :: rawu(lbl:ubl) ! moisture resistance (s/m) real(r8) :: temp1(lbl:ubl) ! relation for potential temperature profile real(r8) :: temp12m(lbl:ubl) ! relation for potential temperature profile applied at 2-m real(r8) :: temp2(lbl:ubl) ! relation for specific humidity profile real(r8) :: temp22m(lbl:ubl) ! relation for specific humidity profile applied at 2-m real(r8) :: thm_g(lbl:ubl) ! intermediate variable (forc_t+0.0098*forc_hgt_t) real(r8) :: thv_g(lbl:ubl) ! virtual potential temperature (K) real(r8) :: dth(lbl:ubl) ! diff of virtual temp. between ref. height and surface real(r8) :: dqh(lbl:ubl) ! diff of humidity between ref. height and surface real(r8) :: zldis(lbl:ubl) ! reference height "minus" zero displacement height (m) real(r8) :: um(lbl:ubl) ! wind speed including the stablity effect (m/s) real(r8) :: obu(lbl:ubl) ! Monin-Obukhov length (m) real(r8) :: taf_numer(lbl:ubl) ! numerator of taf equation (K m/s) real(r8) :: taf_denom(lbl:ubl) ! denominator of taf equation (m/s) real(r8) :: qaf_numer(lbl:ubl) ! numerator of qaf equation (kg m/kg s) real(r8) :: qaf_denom(lbl:ubl) ! denominator of qaf equation (m/s) real(r8) :: wtas(lbl:ubl) ! sensible heat conductance for urban air to atmospheric air (m/s) real(r8) :: wtaq(lbl:ubl) ! latent heat conductance for urban air to atmospheric air (m/s) real(r8) :: wts_sum(lbl:ubl) ! sum of wtas, wtus_roof, wtus_road_perv, wtus_road_imperv, wtus_sunwall, wtus_shadewall real(r8) :: wtq_sum(lbl:ubl) ! sum of wtaq, wtuq_roof, wtuq_road_perv, wtuq_road_imperv, wtuq_sunwall, wtuq_shadewall real(r8) :: beta(lbl:ubl) ! coefficient of convective velocity real(r8) :: zii(lbl:ubl) ! convective boundary layer height (m) real(r8) :: fm(lbl:ubl) ! needed for BGC only to diagnose 10m wind speed real(r8) :: wtus(lbc:ubc) ! sensible heat conductance for urban columns (m/s) real(r8) :: wtuq(lbc:ubc) ! latent heat conductance for urban columns (m/s) integer :: iter ! iteration index real(r8) :: dthv ! diff of vir. poten. temp. between ref. height and surface real(r8) :: tstar ! temperature scaling parameter real(r8) :: qstar ! moisture scaling parameter real(r8) :: thvstar ! virtual potential temperature scaling parameter real(r8) :: wtus_roof(lbl:ubl) ! sensible heat conductance for roof (not scaled) (m/s) real(r8) :: wtuq_roof(lbl:ubl) ! latent heat conductance for roof (not scaled) (m/s) real(r8) :: wtus_road_perv(lbl:ubl) ! sensible heat conductance for pervious road (not scaled) (m/s) real(r8) :: wtuq_road_perv(lbl:ubl) ! latent heat conductance for pervious road (not scaled) (m/s) real(r8) :: wtus_road_imperv(lbl:ubl) ! sensible heat conductance for impervious road (not scaled) (m/s) real(r8) :: wtuq_road_imperv(lbl:ubl) ! latent heat conductance for impervious road (not scaled) (m/s) real(r8) :: wtus_sunwall(lbl:ubl) ! sensible heat conductance for sunwall (not scaled) (m/s) real(r8) :: wtuq_sunwall(lbl:ubl) ! latent heat conductance for sunwall (not scaled) (m/s) real(r8) :: wtus_shadewall(lbl:ubl) ! sensible heat conductance for shadewall (not scaled) (m/s) real(r8) :: wtuq_shadewall(lbl:ubl) ! latent heat conductance for shadewall (not scaled) (m/s) real(r8) :: t_sunwall_innerl(lbl:ubl) ! temperature of inner layer of sunwall (K) real(r8) :: t_shadewall_innerl(lbl:ubl) ! temperature of inner layer of shadewall (K) real(r8) :: t_roof_innerl(lbl:ubl) ! temperature of inner layer of roof (K) real(r8) :: lngth_roof ! length of roof (m) real(r8) :: wc ! convective velocity (m/s) real(r8) :: zeta ! dimensionless height used in Monin-Obukhov theory real(r8) :: eflx_sh_grnd_scale(lbp:ubp) ! scaled sensible heat flux from ground (W/m**2) [+ to atm] real(r8) :: qflx_evap_soi_scale(lbp:ubp) ! scaled soil evaporation (mm H2O/s) (+ = to atm) real(r8) :: eflx_wasteheat_roof(lbl:ubl) ! sensible heat flux from urban heating/cooling sources of waste heat for roof (W/m**2) real(r8) :: eflx_wasteheat_sunwall(lbl:ubl) ! sensible heat flux from urban heating/cooling sources of waste heat for sunwall (W/m**2) real(r8) :: eflx_wasteheat_shadewall(lbl:ubl) ! sensible heat flux from urban heating/cooling sources of waste heat for shadewall (W/m**2) real(r8) :: eflx_heat_from_ac_roof(lbl:ubl) ! sensible heat flux put back into canyon due to heat removal by AC for roof (W/m**2) real(r8) :: eflx_heat_from_ac_sunwall(lbl:ubl) ! sensible heat flux put back into canyon due to heat removal by AC for sunwall (W/m**2) real(r8) :: eflx_heat_from_ac_shadewall(lbl:ubl) ! sensible heat flux put back into canyon due to heat removal by AC for shadewall (W/m**2) real(r8) :: eflx(lbl:ubl) ! total sensible heat flux for error check (W/m**2) real(r8) :: qflx(lbl:ubl) ! total water vapor flux for error check (kg/m**2/s) real(r8) :: eflx_scale(lbl:ubl) ! sum of scaled sensible heat fluxes for urban columns for error check (W/m**2) real(r8) :: qflx_scale(lbl:ubl) ! sum of scaled water vapor fluxes for urban columns for error check (kg/m**2/s) real(r8) :: eflx_err(lbl:ubl) ! sensible heat flux error (W/m**2) real(r8) :: qflx_err(lbl:ubl) ! water vapor flux error (kg/m**2/s) real(r8) :: fwet_roof ! fraction of roof surface that is wet (-) real(r8) :: fwet_road_imperv ! fraction of impervious road surface that is wet (-) integer, parameter :: niters = 3 ! maximum number of iterations for surface temperature integer :: local_secp1(lbl:ubl) ! seconds into current date in local time (sec) real(r8) :: dtime ! land model time step (sec) integer :: year,month,day,secs ! calendar info for current time step logical :: found ! flag in search loop integer :: indexl ! index of first found in search loop integer :: nstep ! time step number real(r8) :: z_d_town_loc(lbl:ubl) ! temporary copy real(r8) :: z_0_town_loc(lbl:ubl) ! temporary copy real(r8), parameter :: lapse_rate = 0.0098_r8 ! Dry adiabatic lapse rate (K/m) real(r8) :: e_ref2m ! 2 m height surface saturated vapor pressure [Pa] real(r8) :: de2mdT ! derivative of 2 m height surface saturated vapor pressure on t_ref2m real(r8) :: qsat_ref2m ! 2 m height surface saturated specific humidity [kg/kg] real(r8) :: dqsat2mdT ! derivative of 2 m height surface saturated specific humidity on t_ref2m !----------------------------------------------------------------------- ! Assign pointers into module urban clumps if ( num_urbanl > 0 )then ht_roof => urban_clump(nc)%ht_roof wtlunit_roof => urban_clump(nc)%wtlunit_roof canyon_hwr => urban_clump(nc)%canyon_hwr wtroad_perv => urban_clump(nc)%wtroad_perv wind_hgt_canyon => urban_clump(nc)%wind_hgt_canyon end if ! Assign local pointers to multi-level derived type members (gridcell level) forc_t => clm_a2l%forc_t forc_th => clm_a2l%forc_th forc_u => clm_a2l%forc_u forc_v => clm_a2l%forc_v forc_rho => clm_a2l%forc_rho forc_q => clm_a2l%forc_q forc_pbot => clm_a2l%forc_pbot londeg => clm3%g%londeg ! Assign local pointers to derived type members (landunit level) pfti => clm3%g%l%pfti pftf => clm3%g%l%pftf coli => clm3%g%l%coli colf => clm3%g%l%colf lgridcell => clm3%g%l%gridcell z_0_town => clm3%g%l%z_0_town z_d_town => clm3%g%l%z_d_town taf => clm3%g%l%lps%taf qaf => clm3%g%l%lps%qaf npfts => clm3%g%l%npfts eflx_traffic => clm3%g%l%lef%eflx_traffic eflx_traffic_factor => clm3%g%l%lef%eflx_traffic_factor eflx_wasteheat => clm3%g%l%lef%eflx_wasteheat eflx_heat_from_ac => clm3%g%l%lef%eflx_heat_from_ac t_building => clm3%g%l%lps%t_building ! Assign local pointers to derived type members (column level) ctype => clm3%g%l%c%itype t_grnd => clm3%g%l%c%ces%t_grnd qg => clm3%g%l%c%cws%qg htvp => clm3%g%l%c%cps%htvp dqgdT => clm3%g%l%c%cws%dqgdT t_soisno => clm3%g%l%c%ces%t_soisno eflx_urban_ac => clm3%g%l%c%cef%eflx_urban_ac eflx_urban_heat => clm3%g%l%c%cef%eflx_urban_heat h2osoi_ice => clm3%g%l%c%cws%h2osoi_ice h2osoi_liq => clm3%g%l%c%cws%h2osoi_liq frac_sno => clm3%g%l%c%cps%frac_sno snowdp => clm3%g%l%c%cps%snowdp h2osno => clm3%g%l%c%cws%h2osno snl => clm3%g%l%c%cps%snl rootr_road_perv => clm3%g%l%c%cps%rootr_road_perv soilalpha_u => clm3%g%l%c%cws%soilalpha_u ! Assign local pointers to derived type members (pft level) pgridcell => clm3%g%l%c%p%gridcell pcolumn => clm3%g%l%c%p%column plandunit => clm3%g%l%c%p%landunit ram1 => clm3%g%l%c%p%pps%ram1 dlrad => clm3%g%l%c%p%pef%dlrad ulrad => clm3%g%l%c%p%pef%ulrad cgrnds => clm3%g%l%c%p%pef%cgrnds cgrndl => clm3%g%l%c%p%pef%cgrndl cgrnd => clm3%g%l%c%p%pef%cgrnd taux => clm3%g%l%c%p%pmf%taux tauy => clm3%g%l%c%p%pmf%tauy eflx_sh_grnd => clm3%g%l%c%p%pef%eflx_sh_grnd eflx_sh_tot => clm3%g%l%c%p%pef%eflx_sh_tot eflx_sh_tot_u => clm3%g%l%c%p%pef%eflx_sh_tot_u qflx_evap_soi => clm3%g%l%c%p%pwf%qflx_evap_soi qflx_tran_veg => clm3%g%l%c%p%pwf%qflx_tran_veg qflx_evap_veg => clm3%g%l%c%p%pwf%qflx_evap_veg qflx_evap_tot => clm3%g%l%c%p%pwf%qflx_evap_tot t_ref2m => clm3%g%l%c%p%pes%t_ref2m q_ref2m => clm3%g%l%c%p%pes%q_ref2m t_ref2m_u => clm3%g%l%c%p%pes%t_ref2m_u t_veg => clm3%g%l%c%p%pes%t_veg rootr => clm3%g%l%c%p%pps%rootr psnsun => clm3%g%l%c%p%pcf%psnsun psnsha => clm3%g%l%c%p%pcf%psnsha forc_hgt_u_pft => clm3%g%l%c%p%pps%forc_hgt_u_pft forc_hgt_t_pft => clm3%g%l%c%p%pps%forc_hgt_t_pft forc_hgt_u_pft => clm3%g%l%c%p%pps%forc_hgt_u_pft forc_hgt_t_pft => clm3%g%l%c%p%pps%forc_hgt_t_pft rh_ref2m => clm3%g%l%c%p%pes%rh_ref2m rh_ref2m_u => clm3%g%l%c%p%pes%rh_ref2m_u ! Define fields that appear on the restart file for non-urban landunits do fl = 1,num_nourbanl l = filter_nourbanl(fl) taf(l) = spval qaf(l) = spval end do ! Get time step nstep = get_nstep() ! Set constants (same as in Biogeophysics1Mod) beta(:) = 1._r8 ! Should be set to the same values as in Biogeophysics1Mod zii(:) = 1000._r8 ! Should be set to the same values as in Biogeophysics1Mod ! Get current date dtime = get_step_size() call get_curr_date (year, month, day, secs) ! Compute canyontop wind using Masson (2000) do fl = 1, num_urbanl l = filter_urbanl(fl) g = lgridcell(l) local_secp1(l) = secs + nint((londeg(g)/degpsec)/dtime)*dtime local_secp1(l) = mod(local_secp1(l),isecspday) ! Error checks if (ht_roof(fl) - z_d_town(l) <= z_0_town(l)) then write (iulog,*) 'aerodynamic parameter error in UrbanFluxes' write (iulog,*) 'h_r - z_d <= z_0' write (iulog,*) 'ht_roof, z_d_town, z_0_town: ', ht_roof(fl), z_d_town(l), & z_0_town(l) write (iulog,*) 'clm model is stopping' call endrun() end if if (forc_hgt_u_pft(pfti(l)) - z_d_town(l) <= z_0_town(l)) then write (iulog,*) 'aerodynamic parameter error in UrbanFluxes' write (iulog,*) 'h_u - z_d <= z_0' write (iulog,*) 'forc_hgt_u_pft, z_d_town, z_0_town: ', forc_hgt_u_pft(pfti(l)), z_d_town(l), & z_0_town(l) write (iulog,*) 'clm model is stopping' call endrun() end if ! Magnitude of atmospheric wind ur(l) = max(1.0_r8,sqrt(forc_u(g)*forc_u(g)+forc_v(g)*forc_v(g))) ! Canyon top wind canyontop_wind(fl) = ur(l) * & log( (ht_roof(fl)-z_d_town(l)) / z_0_town(l) ) / & log( (forc_hgt_u_pft(pfti(l))-z_d_town(l)) / z_0_town(l) ) ! U component of canyon wind if (canyon_hwr(fl) < 0.5_r8) then ! isolated roughness flow canyon_u_wind(fl) = canyontop_wind(fl) * exp( -0.5_r8*canyon_hwr(fl)* & (1._r8-(wind_hgt_canyon(fl)/ht_roof(fl))) ) else if (canyon_hwr(fl) < 1.0_r8) then ! wake interference flow canyon_u_wind(fl) = canyontop_wind(fl) * (1._r8+2._r8*(2._r8/rpi - 1._r8)* & (ht_roof(fl)/(ht_roof(fl)/canyon_hwr(fl)) - 0.5_r8)) * & exp(-0.5_r8*canyon_hwr(fl)*(1._r8-(wind_hgt_canyon(fl)/ht_roof(fl)))) else ! skimming flow canyon_u_wind(fl) = canyontop_wind(fl) * (2._r8/rpi) * & exp(-0.5_r8*canyon_hwr(fl)*(1._r8-(wind_hgt_canyon(fl)/ht_roof(fl)))) end if end do ! Compute fluxes - Follows CLM approach for bare soils (Oleson et al 2004) do fl = 1, num_urbanl l = filter_urbanl(fl) g = lgridcell(l) thm_g(l) = forc_t(g) + lapse_rate*forc_hgt_t_pft(pfti(l)) thv_g(l) = forc_th(g)*(1._r8+0.61_r8*forc_q(g)) dth(l) = thm_g(l)-taf(l) dqh(l) = forc_q(g)-qaf(l) dthv = dth(l)*(1._r8+0.61_r8*forc_q(g))+0.61_r8*forc_th(g)*dqh(l) zldis(l) = forc_hgt_u_pft(pfti(l)) - z_d_town(l) ! Initialize Monin-Obukhov length and wind speed including convective velocity call MoninObukIni(ur(l), thv_g(l), dthv, zldis(l), z_0_town(l), um(l), obu(l)) end do ! Initialize conductances wtus_roof(:) = 0._r8 wtus_road_perv(:) = 0._r8 wtus_road_imperv(:) = 0._r8 wtus_sunwall(:) = 0._r8 wtus_shadewall(:) = 0._r8 wtuq_roof(:) = 0._r8 wtuq_road_perv(:) = 0._r8 wtuq_road_imperv(:) = 0._r8 wtuq_sunwall(:) = 0._r8 wtuq_shadewall(:) = 0._r8 ! Make copies so that array sections are not passed in function calls to friction velocity do fl = 1, num_urbanl l = filter_urbanl(fl) z_d_town_loc(l) = z_d_town(l) z_0_town_loc(l) = z_0_town(l) end do ! Start stability iteration do iter = 1,niters ! Get friction velocity, relation for potential ! temperature and humidity profiles of surface boundary layer. if (num_urbanl .gt. 0) then call FrictionVelocity(lbl, ubl, num_urbanl, filter_urbanl, & z_d_town_loc, z_0_town_loc, z_0_town_loc, z_0_town_loc, & obu, iter, ur, um, ustar, & temp1, temp2, temp12m, temp22m, fm, landunit_index=.true.) end if do fl = 1, num_urbanl l = filter_urbanl(fl) g = lgridcell(l) ! Determine aerodynamic resistance to fluxes from urban canopy air to ! atmosphere ramu(l) = 1._r8/(ustar(l)*ustar(l)/um(l)) rahu(l) = 1._r8/(temp1(l)*ustar(l)) rawu(l) = 1._r8/(temp2(l)*ustar(l)) ! Determine magnitude of canyon wind by using horizontal wind determined ! previously and vertical wind from friction velocity (Masson 2000) canyon_wind(fl) = sqrt(canyon_u_wind(fl)**2._r8 + ustar(l)**2._r8) ! Determine canyon_resistance (currently this single resistance determines the ! resistance from urban surfaces (roof, pervious and impervious road, sunlit and ! shaded walls) to urban canopy air, since it is only dependent on wind speed ! Also from Masson 2000. canyon_resistance(fl) = cpair * forc_rho(g) / (11.8_r8 + 4.2_r8*canyon_wind(fl)) end do ! This is the first term in the equation solutions for urban canopy air temperature ! and specific humidity (numerator) and is a landunit quantity do fl = 1, num_urbanl l = filter_urbanl(fl) g = lgridcell(l) taf_numer(l) = thm_g(l)/rahu(l) taf_denom(l) = 1._r8/rahu(l) qaf_numer(l) = forc_q(g)/rawu(l) qaf_denom(l) = 1._r8/rawu(l) ! First term needed for derivative of heat fluxes wtas(l) = 1._r8/rahu(l) wtaq(l) = 1._r8/rawu(l) end do ! Gather other terms for other urban columns for numerator and denominator of ! equations for urban canopy air temperature and specific humidity do pi = 1,maxpatch_urb do fl = 1,num_urbanl l = filter_urbanl(fl) if ( pi <= npfts(l) ) then c = coli(l) + pi - 1 if (ctype(c) == icol_roof) then ! scaled sensible heat conductance wtus(c) = wtlunit_roof(fl)/canyon_resistance(fl) ! unscaled sensible heat conductance wtus_roof(l) = 1._r8/canyon_resistance(fl) if (snowdp(c) > 0._r8) then fwet_roof = min(snowdp(c)/0.05_r8, 1._r8) else fwet_roof = (max(0._r8, h2osoi_liq(c,1)+h2osoi_ice(c,1))/pondmx_urban)**0.666666666666_r8 fwet_roof = min(fwet_roof,1._r8) end if if (qaf(l) > qg(c)) then fwet_roof = 1._r8 end if ! scaled latent heat conductance wtuq(c) = fwet_roof*(wtlunit_roof(fl)/canyon_resistance(fl)) ! unscaled latent heat conductance wtuq_roof(l) = fwet_roof*(1._r8/canyon_resistance(fl)) ! wasteheat from heating/cooling if (trim(urban_hac) == urban_wasteheat_on) then eflx_wasteheat_roof(l) = ac_wasteheat_factor * eflx_urban_ac(c) + & ht_wasteheat_factor * eflx_urban_heat(c) else eflx_wasteheat_roof(l) = 0._r8 end if ! If air conditioning on, always replace heat removed with heat into canyon if (trim(urban_hac) == urban_hac_on .or. trim(urban_hac) == urban_wasteheat_on) then eflx_heat_from_ac_roof(l) = abs(eflx_urban_ac(c)) else eflx_heat_from_ac_roof(l) = 0._r8 end if else if (ctype(c) == icol_road_perv) then ! scaled sensible heat conductance wtus(c) = wtroad_perv(fl)*(1._r8-wtlunit_roof(fl))/canyon_resistance(fl) ! unscaled sensible heat conductance if (wtroad_perv(fl) > 0._r8) then wtus_road_perv(l) = 1._r8/canyon_resistance(fl) else wtus_road_perv(l) = 0._r8 end if ! scaled latent heat conductance wtuq(c) = wtroad_perv(fl)*(1._r8-wtlunit_roof(fl))/canyon_resistance(fl) ! unscaled latent heat conductance if (wtroad_perv(fl) > 0._r8) then wtuq_road_perv(l) = 1._r8/canyon_resistance(fl) else wtuq_road_perv(l) = 0._r8 end if else if (ctype(c) == icol_road_imperv) then ! scaled sensible heat conductance wtus(c) = (1._r8-wtroad_perv(fl))*(1._r8-wtlunit_roof(fl))/canyon_resistance(fl) ! unscaled sensible heat conductance if ((1._r8-wtroad_perv(fl)) > 0._r8) then wtus_road_imperv(l) = 1._r8/canyon_resistance(fl) else wtus_road_imperv(l) = 0._r8 end if if (snowdp(c) > 0._r8) then fwet_road_imperv = min(snowdp(c)/0.05_r8, 1._r8) else fwet_road_imperv = (max(0._r8, h2osoi_liq(c,1)+h2osoi_ice(c,1))/pondmx_urban)**0.666666666666_r8 fwet_road_imperv = min(fwet_road_imperv,1._r8) end if if (qaf(l) > qg(c)) then fwet_road_imperv = 1._r8 end if ! scaled latent heat conductance wtuq(c) = fwet_road_imperv*(1._r8-wtroad_perv(fl))*(1._r8-wtlunit_roof(fl))/canyon_resistance(fl) ! unscaled latent heat conductance if ((1._r8-wtroad_perv(fl)) > 0._r8) then wtuq_road_imperv(l) = fwet_road_imperv*(1._r8/canyon_resistance(fl)) else wtuq_road_imperv(l) = 0._r8 end if else if (ctype(c) == icol_sunwall) then ! scaled sensible heat conductance wtus(c) = canyon_hwr(fl)*(1._r8-wtlunit_roof(fl))/canyon_resistance(fl) ! unscaled sensible heat conductance wtus_sunwall(l) = 1._r8/canyon_resistance(fl) ! scaled latent heat conductance wtuq(c) = 0._r8 ! unscaled latent heat conductance wtuq_sunwall(l) = 0._r8 ! wasteheat from heating/cooling if (trim(urban_hac) == urban_wasteheat_on) then eflx_wasteheat_sunwall(l) = ac_wasteheat_factor * eflx_urban_ac(c) + & ht_wasteheat_factor * eflx_urban_heat(c) else eflx_wasteheat_sunwall(l) = 0._r8 end if ! If air conditioning on, always replace heat removed with heat into canyon if (trim(urban_hac) == urban_hac_on .or. trim(urban_hac) == urban_wasteheat_on) then eflx_heat_from_ac_sunwall(l) = abs(eflx_urban_ac(c)) else eflx_heat_from_ac_sunwall(l) = 0._r8 end if else if (ctype(c) == icol_shadewall) then ! scaled sensible heat conductance wtus(c) = canyon_hwr(fl)*(1._r8-wtlunit_roof(fl))/canyon_resistance(fl) ! unscaled sensible heat conductance wtus_shadewall(l) = 1._r8/canyon_resistance(fl) ! scaled latent heat conductance wtuq(c) = 0._r8 ! unscaled latent heat conductance wtuq_shadewall(l) = 0._r8 ! wasteheat from heating/cooling if (trim(urban_hac) == urban_wasteheat_on) then eflx_wasteheat_shadewall(l) = ac_wasteheat_factor * eflx_urban_ac(c) + & ht_wasteheat_factor * eflx_urban_heat(c) else eflx_wasteheat_shadewall(l) = 0._r8 end if ! If air conditioning on, always replace heat removed with heat into canyon if (trim(urban_hac) == urban_hac_on .or. trim(urban_hac) == urban_wasteheat_on) then eflx_heat_from_ac_shadewall(l) = abs(eflx_urban_ac(c)) else eflx_heat_from_ac_shadewall(l) = 0._r8 end if else write(iulog,*) 'c, ctype, pi = ', c, ctype(c), pi write(iulog,*) 'Column indices for: shadewall, sunwall, road_imperv, road_perv, roof: ' write(iulog,*) icol_shadewall, icol_sunwall, icol_road_imperv, icol_road_perv, icol_roof call endrun( sub//':: ERROR, ctype out of range' ) end if taf_numer(l) = taf_numer(l) + t_grnd(c)*wtus(c) taf_denom(l) = taf_denom(l) + wtus(c) qaf_numer(l) = qaf_numer(l) + qg(c)*wtuq(c) qaf_denom(l) = qaf_denom(l) + wtuq(c) end if end do end do ! Calculate new urban canopy air temperature and specific humidity do fl = 1, num_urbanl l = filter_urbanl(fl) g = lgridcell(l) ! Total waste heat and heat from AC is sum of heat for walls and roofs ! accounting for different surface areas eflx_wasteheat(l) = wtlunit_roof(fl)*eflx_wasteheat_roof(l) + & (1._r8-wtlunit_roof(fl))*(canyon_hwr(fl)*(eflx_wasteheat_sunwall(l) + & eflx_wasteheat_shadewall(l))) ! Limit wasteheat to ensure that we don't get any unrealistically strong ! positive feedbacks due to AC in a warmer climate eflx_wasteheat(l) = min(eflx_wasteheat(l),wasteheat_limit) eflx_heat_from_ac(l) = wtlunit_roof(fl)*eflx_heat_from_ac_roof(l) + & (1._r8-wtlunit_roof(fl))*(canyon_hwr(fl)*(eflx_heat_from_ac_sunwall(l) + & eflx_heat_from_ac_shadewall(l))) ! Calculate traffic heat flux ! Only comes from impervious road eflx_traffic(l) = (1._r8-wtlunit_roof(fl))*(1._r8-wtroad_perv(fl))* & eflx_traffic_factor(l) taf(l) = taf_numer(l)/taf_denom(l) qaf(l) = qaf_numer(l)/qaf_denom(l) wts_sum(l) = wtas(l) + wtus_roof(l) + wtus_road_perv(l) + & wtus_road_imperv(l) + wtus_sunwall(l) + wtus_shadewall(l) wtq_sum(l) = wtaq(l) + wtuq_roof(l) + wtuq_road_perv(l) + & wtuq_road_imperv(l) + wtuq_sunwall(l) + wtuq_shadewall(l) end do ! This section of code is not required if niters = 1 ! Determine stability using new taf and qaf ! TODO: Some of these constants replicate what is in FrictionVelocity and BareGround fluxes should consildate. EBK do fl = 1, num_urbanl l = filter_urbanl(fl) g = lgridcell(l) dth(l) = thm_g(l)-taf(l) dqh(l) = forc_q(g)-qaf(l) tstar = temp1(l)*dth(l) qstar = temp2(l)*dqh(l) thvstar = tstar*(1._r8+0.61_r8*forc_q(g)) + 0.61_r8*forc_th(g)*qstar zeta = zldis(l)*vkc*grav*thvstar/(ustar(l)**2*thv_g(l)) if (zeta >= 0._r8) then !stable zeta = min(2._r8,max(zeta,0.01_r8)) um(l) = max(ur(l),0.1_r8) else !unstable zeta = max(-100._r8,min(zeta,-0.01_r8)) wc = beta(l)*(-grav*ustar(l)*thvstar*zii(l)/thv_g(l))**0.333_r8 um(l) = sqrt(ur(l)*ur(l) + wc*wc) end if obu(l) = zldis(l)/zeta end do end do ! end iteration ! Determine fluxes from canyon surfaces do f = 1, num_urbanp p = filter_urbanp(f) c = pcolumn(p) g = pgridcell(p) l = plandunit(p) ram1(p) = ramu(l) !pass value to global variable ! Upward and downward canopy longwave are zero ulrad(p) = 0._r8 dlrad(p) = 0._r8 ! Derivative of sensible and latent heat fluxes with respect to ! ground temperature if (ctype(c) == icol_roof) then cgrnds(p) = forc_rho(g) * cpair * (wtas(l) + wtus_road_perv(l) + & wtus_road_imperv(l) + wtus_sunwall(l) + wtus_shadewall(l)) * & (wtus_roof(l)/wts_sum(l)) cgrndl(p) = forc_rho(g) * (wtaq(l) + wtuq_road_perv(l) + & wtuq_road_imperv(l) + wtuq_sunwall(l) + wtuq_shadewall(l)) * & (wtuq_roof(l)/wtq_sum(l))*dqgdT(c) else if (ctype(c) == icol_road_perv) then cgrnds(p) = forc_rho(g) * cpair * (wtas(l) + wtus_roof(l) + & wtus_road_imperv(l) + wtus_sunwall(l) + wtus_shadewall(l)) * & (wtus_road_perv(l)/wts_sum(l)) cgrndl(p) = forc_rho(g) * (wtaq(l) + wtuq_roof(l) + & wtuq_road_imperv(l) + wtuq_sunwall(l) + wtuq_shadewall(l)) * & (wtuq_road_perv(l)/wtq_sum(l))*dqgdT(c) else if (ctype(c) == icol_road_imperv) then cgrnds(p) = forc_rho(g) * cpair * (wtas(l) + wtus_roof(l) + & wtus_road_perv(l) + wtus_sunwall(l) + wtus_shadewall(l)) * & (wtus_road_imperv(l)/wts_sum(l)) cgrndl(p) = forc_rho(g) * (wtaq(l) + wtuq_roof(l) + & wtuq_road_perv(l) + wtuq_sunwall(l) + wtuq_shadewall(l)) * & (wtuq_road_imperv(l)/wtq_sum(l))*dqgdT(c) else if (ctype(c) == icol_sunwall) then cgrnds(p) = forc_rho(g) * cpair * (wtas(l) + wtus_roof(l) + & wtus_road_perv(l) + wtus_road_imperv(l) + wtus_shadewall(l)) * & (wtus_sunwall(l)/wts_sum(l)) cgrndl(p) = 0._r8 else if (ctype(c) == icol_shadewall) then cgrnds(p) = forc_rho(g) * cpair * (wtas(l) + wtus_roof(l) + & wtus_road_perv(l) + wtus_road_imperv(l) + wtus_sunwall(l)) * & (wtus_shadewall(l)/wts_sum(l)) cgrndl(p) = 0._r8 end if cgrnd(p) = cgrnds(p) + cgrndl(p)*htvp(c) ! Surface fluxes of momentum, sensible and latent heat taux(p) = -forc_rho(g)*forc_u(g)/ramu(l) tauy(p) = -forc_rho(g)*forc_v(g)/ramu(l) ! Use new canopy air temperature dth(l) = taf(l) - t_grnd(c) if (ctype(c) == icol_roof) then eflx_sh_grnd(p) = -forc_rho(g)*cpair*wtus_roof(l)*dth(l) else if (ctype(c) == icol_road_perv) then eflx_sh_grnd(p) = -forc_rho(g)*cpair*wtus_road_perv(l)*dth(l) else if (ctype(c) == icol_road_imperv) then eflx_sh_grnd(p) = -forc_rho(g)*cpair*wtus_road_imperv(l)*dth(l) else if (ctype(c) == icol_sunwall) then eflx_sh_grnd(p) = -forc_rho(g)*cpair*wtus_sunwall(l)*dth(l) else if (ctype(c) == icol_shadewall) then eflx_sh_grnd(p) = -forc_rho(g)*cpair*wtus_shadewall(l)*dth(l) end if eflx_sh_tot(p) = eflx_sh_grnd(p) eflx_sh_tot_u(p) = eflx_sh_tot(p) dqh(l) = qaf(l) - qg(c) if (ctype(c) == icol_roof) then qflx_evap_soi(p) = -forc_rho(g)*wtuq_roof(l)*dqh(l) else if (ctype(c) == icol_road_perv) then ! Evaporation assigned to soil term if dew or snow ! or if no liquid water available in soil column if (dqh(l) > 0._r8 .or. frac_sno(c) > 0._r8 .or. soilalpha_u(c) .le. 0._r8) then qflx_evap_soi(p) = -forc_rho(g)*wtuq_road_perv(l)*dqh(l) qflx_tran_veg(p) = 0._r8 ! Otherwise, evaporation assigned to transpiration term else qflx_evap_soi(p) = 0._r8 qflx_tran_veg(p) = -forc_rho(g)*wtuq_road_perv(l)*dqh(l) end if qflx_evap_veg(p) = qflx_tran_veg(p) else if (ctype(c) == icol_road_imperv) then qflx_evap_soi(p) = -forc_rho(g)*wtuq_road_imperv(l)*dqh(l) else if (ctype(c) == icol_sunwall) then qflx_evap_soi(p) = 0._r8 else if (ctype(c) == icol_shadewall) then qflx_evap_soi(p) = 0._r8 end if ! SCALED sensible and latent heat flux for error check eflx_sh_grnd_scale(p) = -forc_rho(g)*cpair*wtus(c)*dth(l) qflx_evap_soi_scale(p) = -forc_rho(g)*wtuq(c)*dqh(l) end do ! Check to see that total sensible and latent heat equal the sum of ! the scaled heat fluxes above do fl = 1, num_urbanl l = filter_urbanl(fl) g = lgridcell(l) eflx(l) = -(forc_rho(g)*cpair/rahu(l))*(thm_g(l) - taf(l)) qflx(l) = -(forc_rho(g)/rawu(l))*(forc_q(g) - qaf(l)) eflx_scale(l) = sum(eflx_sh_grnd_scale(pfti(l):pftf(l))) qflx_scale(l) = sum(qflx_evap_soi_scale(pfti(l):pftf(l))) eflx_err(l) = eflx_scale(l) - eflx(l) qflx_err(l) = qflx_scale(l) - qflx(l) end do found = .false. do fl = 1, num_urbanl l = filter_urbanl(fl) if (abs(eflx_err(l)) > 0.01_r8) then found = .true. indexl = l exit end if end do if ( found ) then write(iulog,*)'WARNING: Total sensible heat does not equal sum of scaled heat fluxes for urban columns ',& ' nstep = ',nstep,' indexl= ',indexl,' eflx_err= ',eflx_err(indexl) if (abs(eflx_err(indexl)) > .01_r8) then write(iulog,*)'clm model is stopping - error is greater than .01 W/m**2' write(iulog,*)'eflx_scale = ',eflx_scale(indexl) write(iulog,*)'eflx_sh_grnd_scale: ',eflx_sh_grnd_scale(pfti(indexl):pftf(indexl)) write(iulog,*)'eflx = ',eflx(indexl) call endrun end if end if found = .false. do fl = 1, num_urbanl l = filter_urbanl(fl) ! 4.e-9 kg/m**2/s = 0.01 W/m**2 if (abs(qflx_err(l)) > 4.e-9_r8) then found = .true. indexl = l exit end if end do if ( found ) then write(iulog,*)'WARNING: Total water vapor flux does not equal sum of scaled water vapor fluxes for urban columns ',& ' nstep = ',nstep,' indexl= ',indexl,' qflx_err= ',qflx_err(indexl) if (abs(qflx_err(indexl)) > 4.e-9_r8) then write(iulog,*)'clm model is stopping - error is greater than 4.e-9 kg/m**2/s' write(iulog,*)'qflx_scale = ',qflx_scale(indexl) write(iulog,*)'qflx = ',qflx(indexl) call endrun end if end if ! Gather terms required to determine internal building temperature do pi = 1,maxpatch_urb do fl = 1,num_urbanl l = filter_urbanl(fl) if ( pi <= npfts(l) ) then c = coli(l) + pi - 1 if (ctype(c) == icol_roof) then t_roof_innerl(l) = t_soisno(c,nlevurb) else if (ctype(c) == icol_sunwall) then t_sunwall_innerl(l) = t_soisno(c,nlevurb) else if (ctype(c) == icol_shadewall) then t_shadewall_innerl(l) = t_soisno(c,nlevurb) end if end if end do end do ! Calculate internal building temperature do fl = 1, num_urbanl l = filter_urbanl(fl) lngth_roof = (ht_roof(fl)/canyon_hwr(fl))*wtlunit_roof(fl)/(1._r8-wtlunit_roof(fl)) t_building(l) = (ht_roof(fl)*(t_shadewall_innerl(l) + t_sunwall_innerl(l)) & +lngth_roof*t_roof_innerl(l))/(2._r8*ht_roof(fl)+lngth_roof) end do ! No roots for urban except for pervious road do j = 1, nlevurb do f = 1, num_urbanp p = filter_urbanp(f) c = pcolumn(p) if (ctype(c) == icol_road_perv) then rootr(p,j) = rootr_road_perv(c,j) else rootr(p,j) = 0._r8 end if end do end do do f = 1, num_urbanp p = filter_urbanp(f) c = pcolumn(p) g = pgridcell(p) l = plandunit(p) ! Use urban canopy air temperature and specific humidity to represent ! 2-m temperature and humidity t_ref2m(p) = taf(l) q_ref2m(p) = qaf(l) t_ref2m_u(p) = taf(l) ! 2 m height relative humidity call QSat(t_ref2m(p), forc_pbot(g), e_ref2m, de2mdT, qsat_ref2m, dqsat2mdT) rh_ref2m(p) = min(100._r8, q_ref2m(p) / qsat_ref2m * 100._r8) rh_ref2m_u(p) = rh_ref2m(p) ! Variables needed by history tape t_veg(p) = forc_t(g) ! Add the following to avoid NaN psnsun(p) = 0._r8 psnsha(p) = 0._r8 clm3%g%l%c%p%pps%lncsun(p) = 0._r8 clm3%g%l%c%p%pps%lncsha(p) = 0._r8 clm3%g%l%c%p%pps%vcmxsun(p) = 0._r8 clm3%g%l%c%p%pps%vcmxsha(p) = 0._r8 end do end subroutine UrbanFluxes end module UrbanMod