From a4d23b512a9173af78cd00d182bbcedd332f7f38 Mon Sep 17 00:00:00 2001 From: Keith Oleson Date: Mon, 4 Nov 2024 13:29:54 -0700 Subject: [PATCH 1/3] Make ALBGRD and ALBGRI active by default --- src/biogeophys/SurfaceAlbedoType.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/biogeophys/SurfaceAlbedoType.F90 b/src/biogeophys/SurfaceAlbedoType.F90 index e90caeecbb..10819a47b6 100644 --- a/src/biogeophys/SurfaceAlbedoType.F90 +++ b/src/biogeophys/SurfaceAlbedoType.F90 @@ -242,12 +242,12 @@ subroutine InitHistory(this, bounds) this%albgrd_col(begc:endc,:) = spval call hist_addfld2d (fname='ALBGRD', units='proportion', type2d='numrad', & avgflag='A', long_name='ground albedo (direct)', & - ptr_col=this%albgrd_col, default='inactive') + ptr_col=this%albgrd_col, default='active') this%albgri_col(begc:endc,:) = spval call hist_addfld2d (fname='ALBGRI', units='proportion', type2d='numrad', & avgflag='A', long_name='ground albedo (indirect)', & - ptr_col=this%albgri_col, default='inactive') + ptr_col=this%albgri_col, default='active') if (use_SSRE) then this%albdSF_patch(begp:endp,:) = spval From e72e947773ce7ca6035c7ba37fc58072c080ae17 Mon Sep 17 00:00:00 2001 From: Keith Oleson Date: Wed, 6 Nov 2024 12:37:51 -0700 Subject: [PATCH 2/3] Implement urbanl and urbanc coszen filters --- src/biogeophys/UrbanAlbedoMod.F90 | 219 +++++++++++++++++------------- 1 file changed, 122 insertions(+), 97 deletions(-) diff --git a/src/biogeophys/UrbanAlbedoMod.F90 b/src/biogeophys/UrbanAlbedoMod.F90 index 2e854a0497..25555b2da8 100644 --- a/src/biogeophys/UrbanAlbedoMod.F90 +++ b/src/biogeophys/UrbanAlbedoMod.F90 @@ -78,9 +78,14 @@ subroutine UrbanAlbedo (bounds, num_urbanl, filter_urbanl, & type(surfalb_type) , intent(inout) :: surfalb_inst ! ! !LOCAL VARIABLES: - integer :: fl,fp,fc,g,l,p,c,ib ! indices + integer :: fl,fp,fc,g,l,p,c,ib,f,fn ! indices and counters integer :: ic ! 0=unit incoming direct; 1=unit incoming diffuse - integer :: num_solar ! counter + integer, allocatable :: filter_urbanl_coszen_gt0(:) ! filter for urban landunits with coszen > 0 + integer :: num_urbanl_coszen_gt0 ! number of urban landunits with coszen > 0 + integer, allocatable :: filter_nourbanl_coszen_gt0(:) ! filter for urban landunits with coszen <= 0 + integer :: num_nourbanl_coszen_gt0 ! number of urban landunits with coszen <= 0 + integer, allocatable :: filter_urbanc_coszen_gt0(:) ! filter for urban columns with coszen > 0 + integer :: num_urbanc_coszen_gt0 ! number of urban columns with coszen > 0 real(r8) :: coszen (bounds%begl:bounds%endl) ! cosine solar zenith angle for next time step (landunit) real(r8) :: zen (bounds%begl:bounds%endl) ! solar zenith angle (radians) real(r8) :: sdir (bounds%begl:bounds%endl, numrad) ! direct beam solar radiation on horizontal surface @@ -165,9 +170,19 @@ subroutine UrbanAlbedo (bounds, num_urbanl, filter_urbanl, & begl => bounds%begl , & vf_sr => urbanparams_inst%vf_sr , & ! Input: [real(r8) (:) ] view factor of sky for road vf_sw => urbanparams_inst%vf_sw , & ! Input: [real(r8) (:) ] view factor of sky for one wall - endl => bounds%endl & + endl => bounds%endl , & + begc => bounds%begc , & + endc => bounds%endc & ) + !TODO KO - Indentation needs to be done properly throughout this module. Will do later + ! to make the changes easier to review now. + + ! Allocate urbanl and urbanc coszen filters + allocate(filter_urbanl_coszen_gt0(bounds%endl-bounds%begl+1)) + allocate(filter_nourbanl_coszen_gt0(bounds%endl-bounds%begl+1)) + allocate(filter_urbanc_coszen_gt0(bounds%endc-bounds%begc+1)) + ! ---------------------------------------------------------------------------- ! Solar declination and cosine solar zenith angle and zenith angle for ! next time step @@ -228,15 +243,36 @@ subroutine UrbanAlbedo (bounds, num_urbanl, filter_urbanl, & end do end do - ! ---------------------------------------------------------------------------- - ! Urban Code - ! ---------------------------------------------------------------------------- - - num_solar = 0 + ! Populate urbanl and urbanc coszen filters + f = 0 + fn = 0 do fl = 1,num_urbanl l = filter_urbanl(fl) - if (coszen(l) > 0._r8) num_solar = num_solar + 1 + if (coszen(l) > 0._r8) then + f = f + 1 + filter_urbanl_coszen_gt0(f) = l + else + fn = fn + 1 + filter_nourbanl_coszen_gt0(fn) = l + end if end do + num_urbanl_coszen_gt0 = f + num_nourbanl_coszen_gt0 = fn + + f = 0 + do fc = 1,num_urbanc + c = filter_urbanc(fc) + l = col%landunit(c) + if (coszen(l) > 0._r8) then + f = f + 1 + filter_urbanc_coszen_gt0(f) = c + end if + end do + num_urbanc_coszen_gt0 = f + + ! ---------------------------------------------------------------------------- + ! Urban Code + ! ---------------------------------------------------------------------------- do ib = 1,numrad do fl = 1,num_urbanl @@ -267,16 +303,15 @@ subroutine UrbanAlbedo (bounds, num_urbanl, filter_urbanl, & end do ! ---------------------------------------------------------------------------- - ! Only do the rest if all coszen are positive + ! Use the urban coszen filters to only calculate quantities when coszen > 0 ! ---------------------------------------------------------------------------- - if (num_solar > 0)then ! Set constants - solar fluxes are per unit incoming flux do ib = 1,numrad - do fl = 1,num_urbanl - l = filter_urbanl(fl) + do fl = 1,num_urbanl_coszen_gt0 + l = filter_urbanl_coszen_gt0(fl) sdir(l,ib) = 1._r8 sdif(l,ib) = 1._r8 end do @@ -287,9 +322,11 @@ subroutine UrbanAlbedo (bounds, num_urbanl, filter_urbanl, & if (num_urbanl > 0) then call incident_direct (bounds, & - num_urbanl, filter_urbanl, & + num_urbanl_coszen_gt0, & + filter_urbanl_coszen_gt0, & + num_nourbanl_coszen_gt0, & + filter_nourbanl_coszen_gt0, & canyon_hwr(begl:endl), & - coszen(begl:endl), & zen(begl:endl), & sdir(begl:endl, :), & sdir_road(begl:endl, :), & @@ -300,9 +337,10 @@ subroutine UrbanAlbedo (bounds, num_urbanl, filter_urbanl, & ! Incident diffuse radiation for ! (a) roof and (b) road and both walls in urban canyon. - if (num_urbanl > 0) then + if (num_urbanl_coszen_gt0 > 0) then call incident_diffuse (bounds, & - num_urbanl, filter_urbanl, & + num_urbanl_coszen_gt0, & + filter_urbanl_coszen_gt0, & canyon_hwr(begl:endl), & sdif(begl:endl, :), & sdif_road(begl:endl, :), & @@ -316,7 +354,8 @@ subroutine UrbanAlbedo (bounds, num_urbanl, filter_urbanl, & ic = 0 call SnowAlbedo(bounds, & num_urbanc, filter_urbanc, & - coszen(begl:endl), & + num_urbanc_coszen_gt0, & + filter_urbanc_coszen_gt0, & ic, & albsnd_roof(begl:endl, :), & albsnd_improad(begl:endl, :), & @@ -326,7 +365,8 @@ subroutine UrbanAlbedo (bounds, num_urbanl, filter_urbanl, & ic = 1 call SnowAlbedo(bounds, & num_urbanc, filter_urbanc, & - coszen(begl:endl), & + num_urbanc_coszen_gt0, & + filter_urbanc_coszen_gt0, & ic, & albsni_roof(begl:endl, :), & albsni_improad(begl:endl, :), & @@ -336,8 +376,8 @@ subroutine UrbanAlbedo (bounds, num_urbanl, filter_urbanl, & ! Combine snow-free and snow albedos do ib = 1,numrad - do fc = 1,num_urbanc - c = filter_urbanc(fc) + do fc = 1,num_urbanc_coszen_gt0 + c = filter_urbanc_coszen_gt0(fc) l = col%landunit(c) if (ctype(c) == icol_roof) then alb_roof_dir_s(l,ib) = alb_roof_dir(l,ib)*(1._r8-frac_sno(c)) & @@ -362,10 +402,10 @@ subroutine UrbanAlbedo (bounds, num_urbanl, filter_urbanl, & ! 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 > 0) then + if (num_urbanl_coszen_gt0 > 0) then call net_solar (bounds, & - num_urbanl, filter_urbanl, & - coszen (begl:endl), & + num_urbanl_coszen_gt0, & + filter_urbanl_coszen_gt0, & canyon_hwr (begl:endl), & wtroad_perv (begl:endl), & sdir (begl:endl, :), & @@ -423,6 +463,7 @@ subroutine UrbanAlbedo (bounds, num_urbanl, filter_urbanl, & albgrd(c,ib) = sref_improad_dir(l,ib) albgri(c,ib) = sref_improad_dif(l,ib) endif + !TODO KO - Clean up these SNICAR comments (not needed?) ! add new snicar albedo variables for history fields if (coszen(l) > 0._r8) then albgrd_hst(c,ib) = albgrd(c,ib) @@ -444,7 +485,6 @@ subroutine UrbanAlbedo (bounds, num_urbanl, filter_urbanl, & ! end add new snicar end do end do - end if end associate @@ -452,7 +492,7 @@ end subroutine UrbanAlbedo !----------------------------------------------------------------------- subroutine SnowAlbedo (bounds , & - num_urbanc, filter_urbanc, coszen, ind , & + num_urbanc, filter_urbanc, num_urbanc_coszen_gt0, filter_urbanc_coszen_gt0, ind , & albsn_roof, albsn_improad, albsn_perroad, & waterstatebulk_inst) ! @@ -466,8 +506,9 @@ subroutine SnowAlbedo (bounds , & type(bounds_type), intent(in) :: bounds integer , intent(in) :: num_urbanc ! number of urban columns in clump integer , intent(in) :: filter_urbanc(:) ! urban column filter + integer , intent(in) :: num_urbanc_coszen_gt0 ! number of urban columns with coszen > 0 + integer , intent(in) :: filter_urbanc_coszen_gt0(:) ! filter for urban columns with coszen > 0 integer , intent(in) :: ind ! 0=direct beam, 1=diffuse radiation - real(r8), intent(in) :: coszen ( bounds%begl: ) ! cosine solar zenith angle [landunit] real(r8), intent(out):: albsn_roof ( bounds%begl: , 1: ) ! roof snow albedo by waveband [landunit, numrad] real(r8), intent(out):: albsn_improad ( bounds%begl: , 1: ) ! impervious road snow albedo by waveband [landunit, numrad] real(r8), intent(out):: albsn_perroad ( bounds%begl: , 1: ) ! pervious road snow albedo by waveband [landunit, numrad] @@ -488,7 +529,6 @@ subroutine SnowAlbedo (bounds , & SHR_ASSERT_ALL_FL(numrad == 2, sourcefile, __LINE__) ! Enforce expected array sizes - SHR_ASSERT_ALL_FL((ubound(coszen) == (/bounds%endl/)), sourcefile, __LINE__) SHR_ASSERT_ALL_FL((ubound(albsn_roof) == (/bounds%endl, numrad/)), sourcefile, __LINE__) SHR_ASSERT_ALL_FL((ubound(albsn_improad) == (/bounds%endl, numrad/)), sourcefile, __LINE__) SHR_ASSERT_ALL_FL((ubound(albsn_perroad) == (/bounds%endl, numrad/)), sourcefile, __LINE__) @@ -497,10 +537,10 @@ subroutine SnowAlbedo (bounds , & caller = 'UrbanAlbedoMod:SnowAlbedo', & h2osno_total = h2osno_total(bounds%begc:bounds%endc)) - do fc = 1,num_urbanc - c = filter_urbanc(fc) + do fc = 1,num_urbanc_coszen_gt0 + c = filter_urbanc_coszen_gt0(fc) l = col%landunit(c) - if (coszen(l) > 0._r8 .and. h2osno_total(c) > 0._r8) then + if (h2osno_total(c) > 0._r8) then if (col%itype(c) == icol_roof) then albsn_roof(l,1) = snal0 albsn_roof(l,2) = snal1 @@ -528,8 +568,9 @@ subroutine SnowAlbedo (bounds , & end subroutine SnowAlbedo !----------------------------------------------------------------------- - subroutine incident_direct (bounds , & - num_urbanl, filter_urbanl, canyon_hwr, coszen, zen , & + subroutine incident_direct (bounds, & + num_urbanl_coszen_gt0, filter_urbanl_coszen_gt0, & + num_nourbanl_coszen_gt0, filter_nourbanl_coszen_gt0, canyon_hwr, zen, & sdir, sdir_road, sdir_sunwall, sdir_shadewall) ! ! !DESCRIPTION: @@ -567,10 +608,11 @@ subroutine incident_direct (bounds , & ! ! !ARGUMENTS: type(bounds_type), intent(in) :: bounds - integer , intent(in) :: num_urbanl ! number of urban landunits - integer , intent(in) :: filter_urbanl(:) ! urban landunit filter + integer , intent(in) :: num_urbanl_coszen_gt0 ! number of urban landunits with coszen > 0 + integer , intent(in) :: filter_urbanl_coszen_gt0(:) ! filter for urban landunits with coszen > 0 + integer , intent(in) :: num_nourbanl_coszen_gt0 ! number of urban landunits with coszen <= 0 + integer , intent(in) :: filter_nourbanl_coszen_gt0(:) ! filter for urban landunits with coszen <= 0 real(r8), intent(in) :: canyon_hwr( bounds%begl: ) ! ratio of building height to street width [landunit] - real(r8), intent(in) :: coszen( bounds%begl: ) ! cosine solar zenith angle [landunit] real(r8), intent(in) :: zen( bounds%begl: ) ! solar zenith angle (radians) [landunit] real(r8), intent(in) :: sdir( bounds%begl: , 1: ) ! direct beam solar radiation incident on horizontal surface [landunit, numrad] real(r8), intent(out) :: sdir_road( bounds%begl: , 1: ) ! direct beam solar radiation incident on road per unit incident flux [landunit, numrad] @@ -595,26 +637,22 @@ subroutine incident_direct (bounds , & ! Enforce expected array sizes SHR_ASSERT_ALL_FL((ubound(canyon_hwr) == (/bounds%endl/)), sourcefile, __LINE__) - SHR_ASSERT_ALL_FL((ubound(coszen) == (/bounds%endl/)), sourcefile, __LINE__) SHR_ASSERT_ALL_FL((ubound(zen) == (/bounds%endl/)), sourcefile, __LINE__) SHR_ASSERT_ALL_FL((ubound(sdir) == (/bounds%endl, numrad/)), sourcefile, __LINE__) SHR_ASSERT_ALL_FL((ubound(sdir_road) == (/bounds%endl, numrad/)), sourcefile, __LINE__) SHR_ASSERT_ALL_FL((ubound(sdir_sunwall) == (/bounds%endl, numrad/)), sourcefile, __LINE__) SHR_ASSERT_ALL_FL((ubound(sdir_shadewall) == (/bounds%endl, numrad/)), sourcefile, __LINE__) - do fl = 1,num_urbanl - l = filter_urbanl(fl) - 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 + do fl = 1,num_urbanl_coszen_gt0 + l = filter_urbanl_coszen_gt0(fl) + theta0(l) = asin(min( (1._r8/(canyon_hwr(l)*tan(max(zen(l),0.000001_r8)))), 1._r8 )) + tanzen(l) = tan(zen(l)) end do do ib = 1,numrad - do fl = 1,num_urbanl - l = filter_urbanl(fl) - if (coszen(l) > 0._r8) then + do fl = 1,num_urbanl_coszen_gt0 + l = filter_urbanl_coszen_gt0(fl) sdir_shadewall(l,ib) = 0._r8 ! incident solar radiation on wall and road integrated over all canyon orientations (0 <= theta <= pi/2) @@ -628,21 +666,20 @@ subroutine incident_direct (bounds , & 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 fl = 1,num_nourbanl_coszen_gt0 + l = filter_nourbanl_coszen_gt0(fl) + sdir_road(l,ib) = 0._r8 + sdir_sunwall(l,ib) = 0._r8 + sdir_shadewall(l,ib) = 0._r8 end do - do fl = 1,num_urbanl - l = filter_urbanl(fl) - 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(subgrid_index=l, subgrid_level=subgrid_level_landunit, msg=errmsg(sourcefile, __LINE__)) - endif + do fl = 1,num_urbanl_coszen_gt0 + l = filter_urbanl_coszen_gt0(fl) + 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(subgrid_index=l, subgrid_level=subgrid_level_landunit, msg=errmsg(sourcefile, __LINE__)) endif end do @@ -650,9 +687,8 @@ subroutine incident_direct (bounds , & ! sum sroad and swall over all canyon orientations (0 <= theta <= pi/2) if (numchk) then - do fl = 1,num_urbanl - l = filter_urbanl(fl) - if (coszen(l) > 0._r8) then + do fl = 1,num_urbanl_coszen_gt0 + l = filter_urbanl_coszen_gt0(fl) sumr = 0._r8 sumw = 0._r8 num = 0._r8 @@ -670,11 +706,9 @@ subroutine incident_direct (bounds , & end do err2(l) = sumr/num - sdir_road(l,ib) err3(l) = sumw/num - sdir_sunwall(l,ib) - endif end do - do fl = 1,num_urbanl - l = filter_urbanl(fl) - if (coszen(l) > 0._r8) then + do fl = 1,num_urbanl_coszen_gt0 + l = filter_urbanl_coszen_gt0(fl) 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' @@ -685,7 +719,6 @@ subroutine incident_direct (bounds , & write (iulog,*) 'clm model is stopping' call endrun(subgrid_index=l, subgrid_level=subgrid_level_landunit, msg=errmsg(sourcefile, __LINE__)) end if - end if end do end if @@ -695,7 +728,7 @@ end subroutine incident_direct !----------------------------------------------------------------------- subroutine incident_diffuse (bounds, & - num_urbanl, filter_urbanl, canyon_hwr, & + num_urbanl_coszen_gt0, filter_urbanl_coszen_gt0, canyon_hwr, & sdif, sdif_road, sdif_sunwall, sdif_shadewall, & urbanparams_inst) ! @@ -707,8 +740,8 @@ subroutine incident_diffuse (bounds, & ! ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds - integer , intent(in) :: num_urbanl ! number of urban landunits - integer , intent(in) :: filter_urbanl(:) ! urban landunit filter + integer , intent(in) :: num_urbanl_coszen_gt0 ! number of urban landunits with coszen > 0 + integer , intent(in) :: filter_urbanl_coszen_gt0(:) ! filter for urban landunits with coszen > 0 real(r8) , intent(in) :: canyon_hwr ( bounds%begl: ) ! ratio of building height to street width [landunit] real(r8) , intent(in) :: sdif ( bounds%begl: , 1: ) ! diffuse solar radiation incident on horizontal surface [landunit, numrad] real(r8) , intent(out) :: sdif_road ( bounds%begl: , 1: ) ! diffuse solar radiation incident on road [landunit, numrad] @@ -738,8 +771,8 @@ subroutine incident_diffuse (bounds, & ! diffuse solar and conservation check. need to convert wall fluxes to ground area - do fl = 1,num_urbanl - l = filter_urbanl(fl) + do fl = 1,num_urbanl_coszen_gt0 + l = filter_urbanl_coszen_gt0(fl) sdif_road(l,ib) = sdif(l,ib) * vf_sr(l) sdif_sunwall(l,ib) = sdif(l,ib) * vf_sw(l) sdif_shadewall(l,ib) = sdif(l,ib) * vf_sw(l) @@ -750,8 +783,8 @@ subroutine incident_diffuse (bounds, & ! error check - do fl = 1, num_urbanl - l = filter_urbanl(fl) + do fl = 1, num_urbanl_coszen_gt0 + l = filter_urbanl_coszen_gt0(fl) if (abs(err(l)) > 0.001_r8) then write (iulog,*) 'urban diffuse solar radiation balance error',err(l) write (iulog,*) 'clm model is stopping' @@ -767,7 +800,7 @@ end subroutine incident_diffuse !----------------------------------------------------------------------- subroutine net_solar (bounds , & - num_urbanl, filter_urbanl, coszen, canyon_hwr, wtroad_perv, sdir, sdif , & + num_urbanl_coszen_gt0, filter_urbanl_coszen_gt0, 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, & @@ -782,9 +815,8 @@ subroutine net_solar (bounds ! ! !ARGUMENTS: type (bounds_type), intent(in) :: bounds - integer , intent(in) :: num_urbanl ! number of urban landunits - integer , intent(in) :: filter_urbanl(:) ! urban landunit filter - real(r8), intent(in) :: coszen ( bounds%begl: ) ! cosine solar zenith angle [landunit] + integer, intent(in) :: num_urbanl_coszen_gt0 ! number of urban landunits with coszen > 0 + integer, intent(in) :: filter_urbanl_coszen_gt0(:) ! filter for urban landunits with coszen > 0 real(r8), intent(in) :: canyon_hwr ( bounds%begl: ) ! ratio of building height to street width [landunit] real(r8), intent(in) :: wtroad_perv ( bounds%begl: ) ! weight of pervious road wrt total road [landunit] real(r8), intent(in) :: sdir ( bounds%begl: , 1: ) ! direct beam solar radiation incident on horizontal surface [landunit, numrad] @@ -897,7 +929,6 @@ subroutine net_solar (bounds !----------------------------------------------------------------------- ! Enforce expected array sizes - SHR_ASSERT_ALL_FL((ubound(coszen) == (/bounds%endl/)), sourcefile, __LINE__) SHR_ASSERT_ALL_FL((ubound(canyon_hwr) == (/bounds%endl/)), sourcefile, __LINE__) SHR_ASSERT_ALL_FL((ubound(wtroad_perv) == (/bounds%endl/)), sourcefile, __LINE__) SHR_ASSERT_ALL_FL((ubound(sdir) == (/bounds%endl, numrad/)), sourcefile, __LINE__) @@ -948,15 +979,14 @@ subroutine net_solar (bounds ! Calculate impervious road - do fl = 1,num_urbanl - l = filter_urbanl(fl) + do fl = 1,num_urbanl_coszen_gt0 + l = filter_urbanl_coszen_gt0(fl) wtroad_imperv(l) = 1._r8 - wtroad_perv(l) end do do ib = 1,numrad - do fl = 1,num_urbanl - l = filter_urbanl(fl) - if (coszen(l) > 0._r8) then + do fl = 1,num_urbanl_coszen_gt0 + l = filter_urbanl_coszen_gt0(fl) ! initial absorption and reflection for road and both walls. ! distribute reflected radiation to sky, road, and walls @@ -1057,7 +1087,6 @@ subroutine net_solar (bounds sref_perroad_dif(l,ib) = perroad_r_sky_dif(l) sref_sunwall_dif(l,ib) = sunwall_r_sky_dif(l) sref_shadewall_dif(l,ib) = shadewall_r_sky_dif(l) - endif end do @@ -1080,9 +1109,8 @@ subroutine net_solar (bounds ! ! do separately for direct beam and diffuse - do fl = 1,num_urbanl - l = filter_urbanl(fl) - if (coszen(l) > 0._r8) then + do fl = 1,num_urbanl_coszen_gt0 + l = filter_urbanl_coszen_gt0(fl) ! reflected direct beam @@ -1279,20 +1307,17 @@ subroutine net_solar (bounds canyon_alb_dif(l) = sref_canyon_dif(l) / max(stot_dif(l), 1.e-06_r8) canyon_alb_dir(l) = sref_canyon_dir(l) / max(stot_dir(l), 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 - l = filter_urbanl(fl) - if (coszen(l) > 0._r8) then - sref_roof_dir(l,ib) = alb_roof_dir(l,ib) * sdir(l,ib) - sref_roof_dif(l,ib) = alb_roof_dif(l,ib) * sdif(l,ib) - sabs_roof_dir(l,ib) = sdir(l,ib) - sref_roof_dir(l,ib) - sabs_roof_dif(l,ib) = sdif(l,ib) - sref_roof_dif(l,ib) - end if + do fl = 1,num_urbanl_coszen_gt0 + l = filter_urbanl_coszen_gt0(fl) + sref_roof_dir(l,ib) = alb_roof_dir(l,ib) * sdir(l,ib) + sref_roof_dif(l,ib) = alb_roof_dif(l,ib) * sdif(l,ib) + sabs_roof_dir(l,ib) = sdir(l,ib) - sref_roof_dir(l,ib) + sabs_roof_dif(l,ib) = sdif(l,ib) - sref_roof_dif(l,ib) end do end do ! end of radiation band loop From 58ae949b74d82f52d1c2e93ef211b669ae870f80 Mon Sep 17 00:00:00 2001 From: Keith Oleson Date: Thu, 7 Nov 2024 17:52:24 -0700 Subject: [PATCH 3/3] Fix indentation --- src/biogeophys/UrbanAlbedoMod.F90 | 960 +++++++++++++++--------------- 1 file changed, 476 insertions(+), 484 deletions(-) diff --git a/src/biogeophys/UrbanAlbedoMod.F90 b/src/biogeophys/UrbanAlbedoMod.F90 index 25555b2da8..3e84f73176 100644 --- a/src/biogeophys/UrbanAlbedoMod.F90 +++ b/src/biogeophys/UrbanAlbedoMod.F90 @@ -175,9 +175,6 @@ subroutine UrbanAlbedo (bounds, num_urbanl, filter_urbanl, & endc => bounds%endc & ) - !TODO KO - Indentation needs to be done properly throughout this module. Will do later - ! to make the changes easier to review now. - ! Allocate urbanl and urbanc coszen filters allocate(filter_urbanl_coszen_gt0(bounds%endl-bounds%begl+1)) allocate(filter_nourbanl_coszen_gt0(bounds%endl-bounds%begl+1)) @@ -307,184 +304,179 @@ subroutine UrbanAlbedo (bounds, num_urbanl, filter_urbanl, & ! ---------------------------------------------------------------------------- - ! Set constants - solar fluxes are per unit incoming flux + ! Set constants - solar fluxes are per unit incoming flux - do ib = 1,numrad - do fl = 1,num_urbanl_coszen_gt0 - l = filter_urbanl_coszen_gt0(fl) - sdir(l,ib) = 1._r8 - sdif(l,ib) = 1._r8 - end do + do ib = 1,numrad + do fl = 1,num_urbanl_coszen_gt0 + l = filter_urbanl_coszen_gt0(fl) + sdir(l,ib) = 1._r8 + sdif(l,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 > 0) then - call incident_direct (bounds, & - num_urbanl_coszen_gt0, & - filter_urbanl_coszen_gt0, & - num_nourbanl_coszen_gt0, & - filter_nourbanl_coszen_gt0, & - canyon_hwr(begl:endl), & - zen(begl:endl), & - sdir(begl:endl, :), & - sdir_road(begl:endl, :), & - sdir_sunwall(begl:endl, :), & - sdir_shadewall(begl:endl, :)) - end if - - ! Incident diffuse radiation for - ! (a) roof and (b) road and both walls in urban canyon. - - if (num_urbanl_coszen_gt0 > 0) then - call incident_diffuse (bounds, & - num_urbanl_coszen_gt0, & - filter_urbanl_coszen_gt0, & - canyon_hwr(begl:endl), & - sdif(begl:endl, :), & - sdif_road(begl:endl, :), & - sdif_sunwall(begl:endl, :), & - sdif_shadewall(begl:endl, :), & - urbanparams_inst) - end if - - ! Get snow albedos for roof and impervious and pervious road - if (num_urbanl > 0) then - ic = 0 - call SnowAlbedo(bounds, & - num_urbanc, filter_urbanc, & - num_urbanc_coszen_gt0, & - filter_urbanc_coszen_gt0, & - ic, & - albsnd_roof(begl:endl, :), & - albsnd_improad(begl:endl, :), & - albsnd_perroad(begl:endl, :), & - waterstatebulk_inst) - - ic = 1 - call SnowAlbedo(bounds, & - num_urbanc, filter_urbanc, & - num_urbanc_coszen_gt0, & - filter_urbanc_coszen_gt0, & - ic, & - albsni_roof(begl:endl, :), & - albsni_improad(begl:endl, :), & - albsni_perroad(begl:endl, :), & - waterstatebulk_inst) - end if - - ! Combine snow-free and snow albedos - do ib = 1,numrad - do fc = 1,num_urbanc_coszen_gt0 - c = filter_urbanc_coszen_gt0(fc) - l = col%landunit(c) - if (ctype(c) == icol_roof) then - alb_roof_dir_s(l,ib) = alb_roof_dir(l,ib)*(1._r8-frac_sno(c)) & - + albsnd_roof(l,ib)*frac_sno(c) - alb_roof_dif_s(l,ib) = alb_roof_dif(l,ib)*(1._r8-frac_sno(c)) & - + albsni_roof(l,ib)*frac_sno(c) - else if (ctype(c) == icol_road_imperv) then - alb_improad_dir_s(l,ib) = alb_improad_dir(l,ib)*(1._r8-frac_sno(c)) & - + albsnd_improad(l,ib)*frac_sno(c) - alb_improad_dif_s(l,ib) = alb_improad_dif(l,ib)*(1._r8-frac_sno(c)) & - + albsni_improad(l,ib)*frac_sno(c) - else if (ctype(c) == icol_road_perv) then - alb_perroad_dir_s(l,ib) = alb_perroad_dir(l,ib)*(1._r8-frac_sno(c)) & - + albsnd_perroad(l,ib)*frac_sno(c) - alb_perroad_dif_s(l,ib) = alb_perroad_dif(l,ib)*(1._r8-frac_sno(c)) & - + albsni_perroad(l,ib)*frac_sno(c) - end if - end do + ! Incident direct beam radiation for + ! (a) roof and (b) road and both walls in urban canyon + + if (num_urbanl > 0) then + call incident_direct (bounds, & + num_urbanl_coszen_gt0, & + filter_urbanl_coszen_gt0, & + num_nourbanl_coszen_gt0, & + filter_nourbanl_coszen_gt0, & + canyon_hwr(begl:endl), & + zen(begl:endl), & + sdir(begl:endl, :), & + sdir_road(begl:endl, :), & + sdir_sunwall(begl:endl, :), & + sdir_shadewall(begl:endl, :)) + end if + + ! Incident diffuse radiation for + ! (a) roof and (b) road and both walls in urban canyon. + + if (num_urbanl_coszen_gt0 > 0) then + call incident_diffuse (bounds, & + num_urbanl_coszen_gt0, & + filter_urbanl_coszen_gt0, & + canyon_hwr(begl:endl), & + sdif(begl:endl, :), & + sdif_road(begl:endl, :), & + sdif_sunwall(begl:endl, :), & + sdif_shadewall(begl:endl, :), & + urbanparams_inst) + end if + + ! Get snow albedos for roof and impervious and pervious road + if (num_urbanl > 0) then + ic = 0 + call SnowAlbedo(bounds, & + num_urbanc, filter_urbanc, & + num_urbanc_coszen_gt0, & + filter_urbanc_coszen_gt0, & + ic, & + albsnd_roof(begl:endl, :), & + albsnd_improad(begl:endl, :), & + albsnd_perroad(begl:endl, :), & + waterstatebulk_inst) + + ic = 1 + call SnowAlbedo(bounds, & + num_urbanc, filter_urbanc, & + num_urbanc_coszen_gt0, & + filter_urbanc_coszen_gt0, & + ic, & + albsni_roof(begl:endl, :), & + albsni_improad(begl:endl, :), & + albsni_perroad(begl:endl, :), & + waterstatebulk_inst) + end if + + ! Combine snow-free and snow albedos + do ib = 1,numrad + do fc = 1,num_urbanc_coszen_gt0 + c = filter_urbanc_coszen_gt0(fc) + l = col%landunit(c) + if (ctype(c) == icol_roof) then + alb_roof_dir_s(l,ib) = alb_roof_dir(l,ib)*(1._r8-frac_sno(c)) & + + albsnd_roof(l,ib)*frac_sno(c) + alb_roof_dif_s(l,ib) = alb_roof_dif(l,ib)*(1._r8-frac_sno(c)) & + + albsni_roof(l,ib)*frac_sno(c) + else if (ctype(c) == icol_road_imperv) then + alb_improad_dir_s(l,ib) = alb_improad_dir(l,ib)*(1._r8-frac_sno(c)) & + + albsnd_improad(l,ib)*frac_sno(c) + alb_improad_dif_s(l,ib) = alb_improad_dif(l,ib)*(1._r8-frac_sno(c)) & + + albsni_improad(l,ib)*frac_sno(c) + else if (ctype(c) == icol_road_perv) then + alb_perroad_dir_s(l,ib) = alb_perroad_dir(l,ib)*(1._r8-frac_sno(c)) & + + albsnd_perroad(l,ib)*frac_sno(c) + alb_perroad_dif_s(l,ib) = alb_perroad_dif(l,ib)*(1._r8-frac_sno(c)) & + + albsni_perroad(l,ib)*frac_sno(c) + end if 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 + ! 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_coszen_gt0 > 0) then + call net_solar (bounds, & + num_urbanl_coszen_gt0, & + filter_urbanl_coszen_gt0, & + canyon_hwr (begl:endl), & + wtroad_perv (begl:endl), & + sdir (begl:endl, :), & + sdif (begl:endl, :), & + alb_improad_dir_s (begl:endl, :), & + alb_perroad_dir_s (begl:endl, :), & + alb_wall_dir (begl:endl, :), & + alb_roof_dir_s (begl:endl, :), & + alb_improad_dif_s (begl:endl, :), & + alb_perroad_dif_s (begl:endl, :), & + alb_wall_dif (begl:endl, :), & + alb_roof_dif_s (begl:endl, :), & + sdir_road (begl:endl, :), & + sdir_sunwall (begl:endl, :), & + sdir_shadewall (begl:endl, :), & + sdif_road (begl:endl, :), & + sdif_sunwall (begl:endl, :), & + sdif_shadewall (begl:endl, :), & + sref_improad_dir (begl:endl, :), & + sref_perroad_dir (begl:endl, :), & + sref_sunwall_dir (begl:endl, :), & + sref_shadewall_dir (begl:endl, :), & + sref_roof_dir (begl:endl, :), & + sref_improad_dif (begl:endl, :), & + sref_perroad_dif (begl:endl, :), & + sref_sunwall_dif (begl:endl, :), & + sref_shadewall_dif (begl:endl, :), & + sref_roof_dif (begl:endl, :), & + urbanparams_inst, solarabs_inst) + end if - if (num_urbanl_coszen_gt0 > 0) then - call net_solar (bounds, & - num_urbanl_coszen_gt0, & - filter_urbanl_coszen_gt0, & - canyon_hwr (begl:endl), & - wtroad_perv (begl:endl), & - sdir (begl:endl, :), & - sdif (begl:endl, :), & - alb_improad_dir_s (begl:endl, :), & - alb_perroad_dir_s (begl:endl, :), & - alb_wall_dir (begl:endl, :), & - alb_roof_dir_s (begl:endl, :), & - alb_improad_dif_s (begl:endl, :), & - alb_perroad_dif_s (begl:endl, :), & - alb_wall_dif (begl:endl, :), & - alb_roof_dif_s (begl:endl, :), & - sdir_road (begl:endl, :), & - sdir_sunwall (begl:endl, :), & - sdir_shadewall (begl:endl, :), & - sdif_road (begl:endl, :), & - sdif_sunwall (begl:endl, :), & - sdif_shadewall (begl:endl, :), & - sref_improad_dir (begl:endl, :), & - sref_perroad_dir (begl:endl, :), & - sref_sunwall_dir (begl:endl, :), & - sref_shadewall_dir (begl:endl, :), & - sref_roof_dir (begl:endl, :), & - sref_improad_dif (begl:endl, :), & - sref_perroad_dif (begl:endl, :), & - sref_sunwall_dif (begl:endl, :), & - sref_shadewall_dif (begl:endl, :), & - sref_roof_dif (begl:endl, :), & - urbanparams_inst, solarabs_inst) - end if + ! ---------------------------------------------------------------------------- + ! Map urban output to surfalb_inst components + ! ---------------------------------------------------------------------------- - ! ---------------------------------------------------------------------------- - ! Map urban output to surfalb_inst components - ! ---------------------------------------------------------------------------- - - ! Set albgrd and albgri (ground albedos) and albd and albi (surface albedos) - - do ib = 1,numrad - do fc = 1,num_urbanc - c = filter_urbanc(fc) - l = col%landunit(c) - if (ctype(c) == icol_roof) then - albgrd(c,ib) = sref_roof_dir(l,ib) - albgri(c,ib) = sref_roof_dif(l,ib) - else if (ctype(c) == icol_sunwall) then - albgrd(c,ib) = sref_sunwall_dir(l,ib) - albgri(c,ib) = sref_sunwall_dif(l,ib) - else if (ctype(c) == icol_shadewall) then - albgrd(c,ib) = sref_shadewall_dir(l,ib) - albgri(c,ib) = sref_shadewall_dif(l,ib) - else if (ctype(c) == icol_road_perv) then - albgrd(c,ib) = sref_perroad_dir(l,ib) - albgri(c,ib) = sref_perroad_dif(l,ib) - else if (ctype(c) == icol_road_imperv) then - albgrd(c,ib) = sref_improad_dir(l,ib) - albgri(c,ib) = sref_improad_dif(l,ib) - endif - !TODO KO - Clean up these SNICAR comments (not needed?) -! add new snicar albedo variables for history fields - if (coszen(l) > 0._r8) then - albgrd_hst(c,ib) = albgrd(c,ib) - albgri_hst(c,ib) = albgri(c,ib) - end if -! end add new snicar - end do - do fp = 1,num_urbanp - p = filter_urbanp(fp) - c = patch%column(p) - l = patch%landunit(p) - albd(p,ib) = albgrd(c,ib) - albi(p,ib) = albgri(c,ib) -! add new snicar albedo variables for history fields - if (coszen(l) > 0._r8) then - albd_hst(p,ib) = albd(p,ib) - albi_hst(p,ib) = albi(p,ib) - end if -! end add new snicar - end do + ! Set albgrd and albgri (ground albedos) and albd and albi (surface albedos) + + do ib = 1,numrad + do fc = 1,num_urbanc + c = filter_urbanc(fc) + l = col%landunit(c) + if (ctype(c) == icol_roof) then + albgrd(c,ib) = sref_roof_dir(l,ib) + albgri(c,ib) = sref_roof_dif(l,ib) + else if (ctype(c) == icol_sunwall) then + albgrd(c,ib) = sref_sunwall_dir(l,ib) + albgri(c,ib) = sref_sunwall_dif(l,ib) + else if (ctype(c) == icol_shadewall) then + albgrd(c,ib) = sref_shadewall_dir(l,ib) + albgri(c,ib) = sref_shadewall_dif(l,ib) + else if (ctype(c) == icol_road_perv) then + albgrd(c,ib) = sref_perroad_dir(l,ib) + albgri(c,ib) = sref_perroad_dif(l,ib) + else if (ctype(c) == icol_road_imperv) then + albgrd(c,ib) = sref_improad_dir(l,ib) + albgri(c,ib) = sref_improad_dif(l,ib) + endif + if (coszen(l) > 0._r8) then + albgrd_hst(c,ib) = albgrd(c,ib) + albgri_hst(c,ib) = albgri(c,ib) + end if + end do + do fp = 1,num_urbanp + p = filter_urbanp(fp) + c = patch%column(p) + l = patch%landunit(p) + albd(p,ib) = albgrd(c,ib) + albi(p,ib) = albgri(c,ib) + if (coszen(l) > 0._r8) then + albd_hst(p,ib) = albd(p,ib) + albi_hst(p,ib) = albi(p,ib) + end if end do + end do end associate @@ -653,19 +645,19 @@ subroutine incident_direct (bounds, & do fl = 1,num_urbanl_coszen_gt0 l = filter_urbanl_coszen_gt0(fl) - sdir_shadewall(l,ib) = 0._r8 + sdir_shadewall(l,ib) = 0._r8 - ! incident solar radiation on wall and road integrated over all canyon orientations (0 <= theta <= pi/2) + ! 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)))) + 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 + ! 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) + swall_projected = (sdir_shadewall(l,ib) + sdir_sunwall(l,ib)) * canyon_hwr(l) + err1(l) = sdir(l,ib) - (sdir_road(l,ib) + swall_projected) end do do fl = 1,num_nourbanl_coszen_gt0 l = filter_nourbanl_coszen_gt0(fl) @@ -689,36 +681,36 @@ subroutine incident_direct (bounds, & if (numchk) then do fl = 1,num_urbanl_coszen_gt0 l = filter_urbanl_coszen_gt0(fl) - 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) + 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) end do do fl = 1,num_urbanl_coszen_gt0 l = filter_urbanl_coszen_gt0(fl) - 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(subgrid_index=l, subgrid_level=subgrid_level_landunit, msg=errmsg(sourcefile, __LINE__)) - 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(subgrid_index=l, subgrid_level=subgrid_level_landunit, msg=errmsg(sourcefile, __LINE__)) - end if + 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(subgrid_index=l, subgrid_level=subgrid_level_landunit, msg=errmsg(sourcefile, __LINE__)) + 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(subgrid_index=l, subgrid_level=subgrid_level_landunit, msg=errmsg(sourcefile, __LINE__)) + end if end do end if @@ -803,8 +795,8 @@ subroutine net_solar (bounds num_urbanl_coszen_gt0, filter_urbanl_coszen_gt0, 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, & + 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 , & urbanparams_inst, solarabs_inst) @@ -988,329 +980,329 @@ subroutine net_solar (bounds do fl = 1,num_urbanl_coszen_gt0 l = filter_urbanl_coszen_gt0(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. + ! 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(l) = 0.0_r8 + road_r_dir(l) = 0.0_r8 + improad_a_dir(l) = (1._r8-alb_improad_dir(l,ib)) * sdir_road(l,ib) + improad_r_dir(l) = alb_improad_dir(l,ib) * sdir_road(l,ib) + improad_r_sky_dir(l) = improad_r_dir(l) * vf_sr(l) + improad_r_sunwall_dir(l) = improad_r_dir(l) * vf_wr(l) + improad_r_shadewall_dir(l) = improad_r_dir(l) * vf_wr(l) + road_a_dir(l) = road_a_dir(l) + improad_a_dir(l)*wtroad_imperv(l) + road_r_dir(l) = road_r_dir(l) + improad_r_dir(l)*wtroad_imperv(l) + + perroad_a_dir(l) = (1._r8-alb_perroad_dir(l,ib)) * sdir_road(l,ib) + perroad_r_dir(l) = alb_perroad_dir(l,ib) * sdir_road(l,ib) + perroad_r_sky_dir(l) = perroad_r_dir(l) * vf_sr(l) + perroad_r_sunwall_dir(l) = perroad_r_dir(l) * vf_wr(l) + perroad_r_shadewall_dir(l) = perroad_r_dir(l) * vf_wr(l) + road_a_dir(l) = road_a_dir(l) + perroad_a_dir(l)*wtroad_perv(l) + road_r_dir(l) = road_r_dir(l) + perroad_r_dir(l)*wtroad_perv(l) + + road_r_sky_dir(l) = road_r_dir(l) * vf_sr(l) + road_r_sunwall_dir(l) = road_r_dir(l) * vf_wr(l) + road_r_shadewall_dir(l) = road_r_dir(l) * vf_wr(l) + + sunwall_a_dir(l) = (1._r8-alb_wall_dir(l,ib)) * sdir_sunwall(l,ib) + sunwall_r_dir(l) = alb_wall_dir(l,ib) * sdir_sunwall(l,ib) + sunwall_r_sky_dir(l) = sunwall_r_dir(l) * vf_sw(l) + sunwall_r_road_dir(l) = sunwall_r_dir(l) * vf_rw(l) + sunwall_r_shadewall_dir(l) = sunwall_r_dir(l) * vf_ww(l) + + shadewall_a_dir(l) = (1._r8-alb_wall_dir(l,ib)) * sdir_shadewall(l,ib) + shadewall_r_dir(l) = alb_wall_dir(l,ib) * sdir_shadewall(l,ib) + shadewall_r_sky_dir(l) = shadewall_r_dir(l) * vf_sw(l) + shadewall_r_road_dir(l) = shadewall_r_dir(l) * vf_rw(l) + shadewall_r_sunwall_dir(l) = shadewall_r_dir(l) * vf_ww(l) + + ! diffuse + + road_a_dif(l) = 0.0_r8 + road_r_dif(l) = 0.0_r8 + improad_a_dif(l) = (1._r8-alb_improad_dif(l,ib)) * sdif_road(l,ib) + improad_r_dif(l) = alb_improad_dif(l,ib) * sdif_road(l,ib) + improad_r_sky_dif(l) = improad_r_dif(l) * vf_sr(l) + improad_r_sunwall_dif(l) = improad_r_dif(l) * vf_wr(l) + improad_r_shadewall_dif(l) = improad_r_dif(l) * vf_wr(l) + road_a_dif(l) = road_a_dif(l) + improad_a_dif(l)*wtroad_imperv(l) + road_r_dif(l) = road_r_dif(l) + improad_r_dif(l)*wtroad_imperv(l) + + perroad_a_dif(l) = (1._r8-alb_perroad_dif(l,ib)) * sdif_road(l,ib) + perroad_r_dif(l) = alb_perroad_dif(l,ib) * sdif_road(l,ib) + perroad_r_sky_dif(l) = perroad_r_dif(l) * vf_sr(l) + perroad_r_sunwall_dif(l) = perroad_r_dif(l) * vf_wr(l) + perroad_r_shadewall_dif(l) = perroad_r_dif(l) * vf_wr(l) + road_a_dif(l) = road_a_dif(l) + perroad_a_dif(l)*wtroad_perv(l) + road_r_dif(l) = road_r_dif(l) + perroad_r_dif(l)*wtroad_perv(l) + + road_r_sky_dif(l) = road_r_dif(l) * vf_sr(l) + road_r_sunwall_dif(l) = road_r_dif(l) * vf_wr(l) + road_r_shadewall_dif(l) = road_r_dif(l) * vf_wr(l) + + sunwall_a_dif(l) = (1._r8-alb_wall_dif(l,ib)) * sdif_sunwall(l,ib) + sunwall_r_dif(l) = alb_wall_dif(l,ib) * sdif_sunwall(l,ib) + sunwall_r_sky_dif(l) = sunwall_r_dif(l) * vf_sw(l) + sunwall_r_road_dif(l) = sunwall_r_dif(l) * vf_rw(l) + sunwall_r_shadewall_dif(l) = sunwall_r_dif(l) * vf_ww(l) + + shadewall_a_dif(l) = (1._r8-alb_wall_dif(l,ib)) * sdif_shadewall(l,ib) + shadewall_r_dif(l) = alb_wall_dif(l,ib) * sdif_shadewall(l,ib) + shadewall_r_sky_dif(l) = shadewall_r_dif(l) * vf_sw(l) + shadewall_r_road_dif(l) = shadewall_r_dif(l) * vf_rw(l) + shadewall_r_sunwall_dif(l) = shadewall_r_dif(l) * vf_ww(l) + + ! initialize sum of direct and diffuse solar absorption and reflection for road and both walls + + sabs_improad_dir(l,ib) = improad_a_dir(l) + sabs_perroad_dir(l,ib) = perroad_a_dir(l) + sabs_sunwall_dir(l,ib) = sunwall_a_dir(l) + sabs_shadewall_dir(l,ib) = shadewall_a_dir(l) + + sabs_improad_dif(l,ib) = improad_a_dif(l) + sabs_perroad_dif(l,ib) = perroad_a_dif(l) + sabs_sunwall_dif(l,ib) = sunwall_a_dif(l) + sabs_shadewall_dif(l,ib) = shadewall_a_dif(l) + + sref_improad_dir(l,ib) = improad_r_sky_dir(l) + sref_perroad_dir(l,ib) = perroad_r_sky_dir(l) + sref_sunwall_dir(l,ib) = sunwall_r_sky_dir(l) + sref_shadewall_dir(l,ib) = shadewall_r_sky_dir(l) + + sref_improad_dif(l,ib) = improad_r_sky_dif(l) + sref_perroad_dif(l,ib) = perroad_r_sky_dif(l) + sref_sunwall_dif(l,ib) = sunwall_r_sky_dif(l) + sref_shadewall_dif(l,ib) = shadewall_r_sky_dif(l) + + 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_coszen_gt0 + l = filter_urbanl_coszen_gt0(fl) + + ! reflected direct beam + + do iter_dir = 1, n + ! step (1) + + stot(l) = (sunwall_r_road_dir(l) + shadewall_r_road_dir(l))*canyon_hwr(l) + + road_a_dir(l) = 0.0_r8 + road_r_dir(l) = 0.0_r8 + improad_a_dir(l) = (1._r8-alb_improad_dir(l,ib)) * stot(l) + improad_r_dir(l) = alb_improad_dir(l,ib) * stot(l) + road_a_dir(l) = road_a_dir(l) + improad_a_dir(l)*wtroad_imperv(l) + road_r_dir(l) = road_r_dir(l) + improad_r_dir(l)*wtroad_imperv(l) + perroad_a_dir(l) = (1._r8-alb_perroad_dir(l,ib)) * stot(l) + perroad_r_dir(l) = alb_perroad_dir(l,ib) * stot(l) + road_a_dir(l) = road_a_dir(l) + perroad_a_dir(l)*wtroad_perv(l) + road_r_dir(l) = road_r_dir(l) + perroad_r_dir(l)*wtroad_perv(l) + + stot(l) = road_r_sunwall_dir(l)/canyon_hwr(l) + shadewall_r_sunwall_dir(l) + sunwall_a_dir(l) = (1._r8-alb_wall_dir(l,ib)) * stot(l) + sunwall_r_dir(l) = alb_wall_dir(l,ib) * stot(l) + + stot(l) = road_r_shadewall_dir(l)/canyon_hwr(l) + sunwall_r_shadewall_dir(l) + shadewall_a_dir(l) = (1._r8-alb_wall_dir(l,ib)) * stot(l) + shadewall_r_dir(l) = alb_wall_dir(l,ib) * stot(l) + + ! step (2) - ! direct beam + sabs_improad_dir(l,ib) = sabs_improad_dir(l,ib) + improad_a_dir(l) + sabs_perroad_dir(l,ib) = sabs_perroad_dir(l,ib) + perroad_a_dir(l) + sabs_sunwall_dir(l,ib) = sabs_sunwall_dir(l,ib) + sunwall_a_dir(l) + sabs_shadewall_dir(l,ib) = sabs_shadewall_dir(l,ib) + shadewall_a_dir(l) + + ! step (3) - road_a_dir(l) = 0.0_r8 - road_r_dir(l) = 0.0_r8 - improad_a_dir(l) = (1._r8-alb_improad_dir(l,ib)) * sdir_road(l,ib) - improad_r_dir(l) = alb_improad_dir(l,ib) * sdir_road(l,ib) improad_r_sky_dir(l) = improad_r_dir(l) * vf_sr(l) improad_r_sunwall_dir(l) = improad_r_dir(l) * vf_wr(l) improad_r_shadewall_dir(l) = improad_r_dir(l) * vf_wr(l) - road_a_dir(l) = road_a_dir(l) + improad_a_dir(l)*wtroad_imperv(l) - road_r_dir(l) = road_r_dir(l) + improad_r_dir(l)*wtroad_imperv(l) - perroad_a_dir(l) = (1._r8-alb_perroad_dir(l,ib)) * sdir_road(l,ib) - perroad_r_dir(l) = alb_perroad_dir(l,ib) * sdir_road(l,ib) perroad_r_sky_dir(l) = perroad_r_dir(l) * vf_sr(l) perroad_r_sunwall_dir(l) = perroad_r_dir(l) * vf_wr(l) perroad_r_shadewall_dir(l) = perroad_r_dir(l) * vf_wr(l) - road_a_dir(l) = road_a_dir(l) + perroad_a_dir(l)*wtroad_perv(l) - road_r_dir(l) = road_r_dir(l) + perroad_r_dir(l)*wtroad_perv(l) road_r_sky_dir(l) = road_r_dir(l) * vf_sr(l) road_r_sunwall_dir(l) = road_r_dir(l) * vf_wr(l) road_r_shadewall_dir(l) = road_r_dir(l) * vf_wr(l) - sunwall_a_dir(l) = (1._r8-alb_wall_dir(l,ib)) * sdir_sunwall(l,ib) - sunwall_r_dir(l) = alb_wall_dir(l,ib) * sdir_sunwall(l,ib) sunwall_r_sky_dir(l) = sunwall_r_dir(l) * vf_sw(l) sunwall_r_road_dir(l) = sunwall_r_dir(l) * vf_rw(l) sunwall_r_shadewall_dir(l) = sunwall_r_dir(l) * vf_ww(l) - shadewall_a_dir(l) = (1._r8-alb_wall_dir(l,ib)) * sdir_shadewall(l,ib) - shadewall_r_dir(l) = alb_wall_dir(l,ib) * sdir_shadewall(l,ib) shadewall_r_sky_dir(l) = shadewall_r_dir(l) * vf_sw(l) shadewall_r_road_dir(l) = shadewall_r_dir(l) * vf_rw(l) shadewall_r_sunwall_dir(l) = shadewall_r_dir(l) * vf_ww(l) - ! diffuse + ! step (4) + + sref_improad_dir(l,ib) = sref_improad_dir(l,ib) + improad_r_sky_dir(l) + sref_perroad_dir(l,ib) = sref_perroad_dir(l,ib) + perroad_r_sky_dir(l) + sref_sunwall_dir(l,ib) = sref_sunwall_dir(l,ib) + sunwall_r_sky_dir(l) + sref_shadewall_dir(l,ib) = sref_shadewall_dir(l,ib) + shadewall_r_sky_dir(l) + + ! step (5) + + crit = max(road_a_dir(l), sunwall_a_dir(l), shadewall_a_dir(l)) + 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(subgrid_index=l, subgrid_level=subgrid_level_landunit, msg=errmsg(sourcefile, __LINE__)) + endif + + ! reflected diffuse + + do iter_dif = 1, n + ! step (1) + + stot(l) = (sunwall_r_road_dif(l) + shadewall_r_road_dif(l))*canyon_hwr(l) + road_a_dif(l) = 0.0_r8 + road_r_dif(l) = 0.0_r8 + improad_a_dif(l) = (1._r8-alb_improad_dif(l,ib)) * stot(l) + improad_r_dif(l) = alb_improad_dif(l,ib) * stot(l) + road_a_dif(l) = road_a_dif(l) + improad_a_dif(l)*wtroad_imperv(l) + road_r_dif(l) = road_r_dif(l) + improad_r_dif(l)*wtroad_imperv(l) + perroad_a_dif(l) = (1._r8-alb_perroad_dif(l,ib)) * stot(l) + perroad_r_dif(l) = alb_perroad_dif(l,ib) * stot(l) + road_a_dif(l) = road_a_dif(l) + perroad_a_dif(l)*wtroad_perv(l) + road_r_dif(l) = road_r_dif(l) + perroad_r_dif(l)*wtroad_perv(l) + + stot(l) = road_r_sunwall_dif(l)/canyon_hwr(l) + shadewall_r_sunwall_dif(l) + sunwall_a_dif(l) = (1._r8-alb_wall_dif(l,ib)) * stot(l) + sunwall_r_dif(l) = alb_wall_dif(l,ib) * stot(l) + + stot(l) = road_r_shadewall_dif(l)/canyon_hwr(l) + sunwall_r_shadewall_dif(l) + shadewall_a_dif(l) = (1._r8-alb_wall_dif(l,ib)) * stot(l) + shadewall_r_dif(l) = alb_wall_dif(l,ib) * stot(l) + + ! step (2) + + sabs_improad_dif(l,ib) = sabs_improad_dif(l,ib) + improad_a_dif(l) + sabs_perroad_dif(l,ib) = sabs_perroad_dif(l,ib) + perroad_a_dif(l) + sabs_sunwall_dif(l,ib) = sabs_sunwall_dif(l,ib) + sunwall_a_dif(l) + sabs_shadewall_dif(l,ib) = sabs_shadewall_dif(l,ib) + shadewall_a_dif(l) + + ! step (3) - road_a_dif(l) = 0.0_r8 - road_r_dif(l) = 0.0_r8 - improad_a_dif(l) = (1._r8-alb_improad_dif(l,ib)) * sdif_road(l,ib) - improad_r_dif(l) = alb_improad_dif(l,ib) * sdif_road(l,ib) improad_r_sky_dif(l) = improad_r_dif(l) * vf_sr(l) improad_r_sunwall_dif(l) = improad_r_dif(l) * vf_wr(l) improad_r_shadewall_dif(l) = improad_r_dif(l) * vf_wr(l) - road_a_dif(l) = road_a_dif(l) + improad_a_dif(l)*wtroad_imperv(l) - road_r_dif(l) = road_r_dif(l) + improad_r_dif(l)*wtroad_imperv(l) - perroad_a_dif(l) = (1._r8-alb_perroad_dif(l,ib)) * sdif_road(l,ib) - perroad_r_dif(l) = alb_perroad_dif(l,ib) * sdif_road(l,ib) perroad_r_sky_dif(l) = perroad_r_dif(l) * vf_sr(l) perroad_r_sunwall_dif(l) = perroad_r_dif(l) * vf_wr(l) perroad_r_shadewall_dif(l) = perroad_r_dif(l) * vf_wr(l) - road_a_dif(l) = road_a_dif(l) + perroad_a_dif(l)*wtroad_perv(l) - road_r_dif(l) = road_r_dif(l) + perroad_r_dif(l)*wtroad_perv(l) road_r_sky_dif(l) = road_r_dif(l) * vf_sr(l) road_r_sunwall_dif(l) = road_r_dif(l) * vf_wr(l) road_r_shadewall_dif(l) = road_r_dif(l) * vf_wr(l) - sunwall_a_dif(l) = (1._r8-alb_wall_dif(l,ib)) * sdif_sunwall(l,ib) - sunwall_r_dif(l) = alb_wall_dif(l,ib) * sdif_sunwall(l,ib) sunwall_r_sky_dif(l) = sunwall_r_dif(l) * vf_sw(l) sunwall_r_road_dif(l) = sunwall_r_dif(l) * vf_rw(l) sunwall_r_shadewall_dif(l) = sunwall_r_dif(l) * vf_ww(l) - shadewall_a_dif(l) = (1._r8-alb_wall_dif(l,ib)) * sdif_shadewall(l,ib) - shadewall_r_dif(l) = alb_wall_dif(l,ib) * sdif_shadewall(l,ib) shadewall_r_sky_dif(l) = shadewall_r_dif(l) * vf_sw(l) - shadewall_r_road_dif(l) = shadewall_r_dif(l) * vf_rw(l) - shadewall_r_sunwall_dif(l) = shadewall_r_dif(l) * vf_ww(l) - - ! initialize sum of direct and diffuse solar absorption and reflection for road and both walls - - sabs_improad_dir(l,ib) = improad_a_dir(l) - sabs_perroad_dir(l,ib) = perroad_a_dir(l) - sabs_sunwall_dir(l,ib) = sunwall_a_dir(l) - sabs_shadewall_dir(l,ib) = shadewall_a_dir(l) - - sabs_improad_dif(l,ib) = improad_a_dif(l) - sabs_perroad_dif(l,ib) = perroad_a_dif(l) - sabs_sunwall_dif(l,ib) = sunwall_a_dif(l) - sabs_shadewall_dif(l,ib) = shadewall_a_dif(l) - - sref_improad_dir(l,ib) = improad_r_sky_dir(l) - sref_perroad_dir(l,ib) = perroad_r_sky_dir(l) - sref_sunwall_dir(l,ib) = sunwall_r_sky_dir(l) - sref_shadewall_dir(l,ib) = shadewall_r_sky_dir(l) - - sref_improad_dif(l,ib) = improad_r_sky_dif(l) - sref_perroad_dif(l,ib) = perroad_r_sky_dif(l) - sref_sunwall_dif(l,ib) = sunwall_r_sky_dif(l) - sref_shadewall_dif(l,ib) = shadewall_r_sky_dif(l) - - 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_coszen_gt0 - l = filter_urbanl_coszen_gt0(fl) - - ! reflected direct beam - - do iter_dir = 1, n - ! step (1) - - stot(l) = (sunwall_r_road_dir(l) + shadewall_r_road_dir(l))*canyon_hwr(l) + shadewall_r_road_dif(l) = shadewall_r_dif(l) * vf_rw(l) + shadewall_r_sunwall_dif(l) = shadewall_r_dif(l) * vf_ww(l) - road_a_dir(l) = 0.0_r8 - road_r_dir(l) = 0.0_r8 - improad_a_dir(l) = (1._r8-alb_improad_dir(l,ib)) * stot(l) - improad_r_dir(l) = alb_improad_dir(l,ib) * stot(l) - road_a_dir(l) = road_a_dir(l) + improad_a_dir(l)*wtroad_imperv(l) - road_r_dir(l) = road_r_dir(l) + improad_r_dir(l)*wtroad_imperv(l) - perroad_a_dir(l) = (1._r8-alb_perroad_dir(l,ib)) * stot(l) - perroad_r_dir(l) = alb_perroad_dir(l,ib) * stot(l) - road_a_dir(l) = road_a_dir(l) + perroad_a_dir(l)*wtroad_perv(l) - road_r_dir(l) = road_r_dir(l) + perroad_r_dir(l)*wtroad_perv(l) + ! step (4) - stot(l) = road_r_sunwall_dir(l)/canyon_hwr(l) + shadewall_r_sunwall_dir(l) - sunwall_a_dir(l) = (1._r8-alb_wall_dir(l,ib)) * stot(l) - sunwall_r_dir(l) = alb_wall_dir(l,ib) * stot(l) + sref_improad_dif(l,ib) = sref_improad_dif(l,ib) + improad_r_sky_dif(l) + sref_perroad_dif(l,ib) = sref_perroad_dif(l,ib) + perroad_r_sky_dif(l) + sref_sunwall_dif(l,ib) = sref_sunwall_dif(l,ib) + sunwall_r_sky_dif(l) + sref_shadewall_dif(l,ib) = sref_shadewall_dif(l,ib) + shadewall_r_sky_dif(l) - stot(l) = road_r_shadewall_dir(l)/canyon_hwr(l) + sunwall_r_shadewall_dir(l) - shadewall_a_dir(l) = (1._r8-alb_wall_dir(l,ib)) * stot(l) - shadewall_r_dir(l) = alb_wall_dir(l,ib) * stot(l) + ! step (5) - ! step (2) - - sabs_improad_dir(l,ib) = sabs_improad_dir(l,ib) + improad_a_dir(l) - sabs_perroad_dir(l,ib) = sabs_perroad_dir(l,ib) + perroad_a_dir(l) - sabs_sunwall_dir(l,ib) = sabs_sunwall_dir(l,ib) + sunwall_a_dir(l) - sabs_shadewall_dir(l,ib) = sabs_shadewall_dir(l,ib) + shadewall_a_dir(l) + crit = max(road_a_dif(l), sunwall_a_dif(l), shadewall_a_dif(l)) + 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(subgrid_index=l, subgrid_level=subgrid_level_landunit, msg=errmsg(sourcefile, __LINE__)) + endif - ! step (3) + ! total reflected by canyon - sum of solar reflection to sky from canyon. + ! project wall fluxes to horizontal surface + + sref_canyon_dir(l) = 0.0_r8 + sref_canyon_dif(l) = 0.0_r8 + sref_canyon_dir(l) = sref_canyon_dir(l) + sref_improad_dir(l,ib)*wtroad_imperv(l) + sref_canyon_dif(l) = sref_canyon_dif(l) + sref_improad_dif(l,ib)*wtroad_imperv(l) + sref_canyon_dir(l) = sref_canyon_dir(l) + sref_perroad_dir(l,ib)*wtroad_perv(l) + sref_canyon_dif(l) = sref_canyon_dif(l) + sref_perroad_dif(l,ib)*wtroad_perv(l) + sref_canyon_dir(l) = sref_canyon_dir(l) + (sref_sunwall_dir(l,ib) + sref_shadewall_dir(l,ib))*canyon_hwr(l) + sref_canyon_dif(l) = sref_canyon_dif(l) + (sref_sunwall_dif(l,ib) + sref_shadewall_dif(l,ib))*canyon_hwr(l) + + ! total absorbed by canyon. project wall fluxes to horizontal surface + + sabs_canyon_dir(l) = 0.0_r8 + sabs_canyon_dif(l) = 0.0_r8 + sabs_canyon_dir(l) = sabs_canyon_dir(l) + sabs_improad_dir(l,ib)*wtroad_imperv(l) + sabs_canyon_dif(l) = sabs_canyon_dif(l) + sabs_improad_dif(l,ib)*wtroad_imperv(l) + sabs_canyon_dir(l) = sabs_canyon_dir(l) + sabs_perroad_dir(l,ib)*wtroad_perv(l) + sabs_canyon_dif(l) = sabs_canyon_dif(l) + sabs_perroad_dif(l,ib)*wtroad_perv(l) + sabs_canyon_dir(l) = sabs_canyon_dir(l) + (sabs_sunwall_dir(l,ib) + sabs_shadewall_dir(l,ib))*canyon_hwr(l) + sabs_canyon_dif(l) = sabs_canyon_dif(l) + (sabs_sunwall_dif(l,ib) + sabs_shadewall_dif(l,ib))*canyon_hwr(l) + + ! 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(l) = sdir_road(l,ib) + (sdir_sunwall(l,ib) + sdir_shadewall(l,ib))*canyon_hwr(l) + stot_dif(l) = sdif_road(l,ib) + (sdif_sunwall(l,ib) + sdif_shadewall(l,ib))*canyon_hwr(l) + + err = stot_dir(l) + stot_dif(l) & + - (sabs_canyon_dir(l) + sabs_canyon_dif(l) + sref_canyon_dir(l) + sref_canyon_dif(l)) + 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(l) + write(iulog,*)' stot_dif = ',stot_dif(l) + write(iulog,*)' sabs_canyon_dir = ',sabs_canyon_dir(l) + write(iulog,*)' sabs_canyon_dif = ',sabs_canyon_dif(l) + write(iulog,*)' sref_canyon_dir = ',sref_canyon_dir(l) + write(iulog,*)' sref_canyon_dif = ',sref_canyon_dir(l) + write(iulog,*) 'clm model is stopping' + call endrun(subgrid_index=l, subgrid_level=subgrid_level_landunit, msg=errmsg(sourcefile, __LINE__)) + endif - improad_r_sky_dir(l) = improad_r_dir(l) * vf_sr(l) - improad_r_sunwall_dir(l) = improad_r_dir(l) * vf_wr(l) - improad_r_shadewall_dir(l) = improad_r_dir(l) * vf_wr(l) + ! canyon albedo - perroad_r_sky_dir(l) = perroad_r_dir(l) * vf_sr(l) - perroad_r_sunwall_dir(l) = perroad_r_dir(l) * vf_wr(l) - perroad_r_shadewall_dir(l) = perroad_r_dir(l) * vf_wr(l) - - road_r_sky_dir(l) = road_r_dir(l) * vf_sr(l) - road_r_sunwall_dir(l) = road_r_dir(l) * vf_wr(l) - road_r_shadewall_dir(l) = road_r_dir(l) * vf_wr(l) - - sunwall_r_sky_dir(l) = sunwall_r_dir(l) * vf_sw(l) - sunwall_r_road_dir(l) = sunwall_r_dir(l) * vf_rw(l) - sunwall_r_shadewall_dir(l) = sunwall_r_dir(l) * vf_ww(l) - - shadewall_r_sky_dir(l) = shadewall_r_dir(l) * vf_sw(l) - shadewall_r_road_dir(l) = shadewall_r_dir(l) * vf_rw(l) - shadewall_r_sunwall_dir(l) = shadewall_r_dir(l) * vf_ww(l) - - ! step (4) - - sref_improad_dir(l,ib) = sref_improad_dir(l,ib) + improad_r_sky_dir(l) - sref_perroad_dir(l,ib) = sref_perroad_dir(l,ib) + perroad_r_sky_dir(l) - sref_sunwall_dir(l,ib) = sref_sunwall_dir(l,ib) + sunwall_r_sky_dir(l) - sref_shadewall_dir(l,ib) = sref_shadewall_dir(l,ib) + shadewall_r_sky_dir(l) - - ! step (5) - - crit = max(road_a_dir(l), sunwall_a_dir(l), shadewall_a_dir(l)) - 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(subgrid_index=l, subgrid_level=subgrid_level_landunit, msg=errmsg(sourcefile, __LINE__)) - endif - - ! reflected diffuse - - do iter_dif = 1, n - ! step (1) - - stot(l) = (sunwall_r_road_dif(l) + shadewall_r_road_dif(l))*canyon_hwr(l) - road_a_dif(l) = 0.0_r8 - road_r_dif(l) = 0.0_r8 - improad_a_dif(l) = (1._r8-alb_improad_dif(l,ib)) * stot(l) - improad_r_dif(l) = alb_improad_dif(l,ib) * stot(l) - road_a_dif(l) = road_a_dif(l) + improad_a_dif(l)*wtroad_imperv(l) - road_r_dif(l) = road_r_dif(l) + improad_r_dif(l)*wtroad_imperv(l) - perroad_a_dif(l) = (1._r8-alb_perroad_dif(l,ib)) * stot(l) - perroad_r_dif(l) = alb_perroad_dif(l,ib) * stot(l) - road_a_dif(l) = road_a_dif(l) + perroad_a_dif(l)*wtroad_perv(l) - road_r_dif(l) = road_r_dif(l) + perroad_r_dif(l)*wtroad_perv(l) - - stot(l) = road_r_sunwall_dif(l)/canyon_hwr(l) + shadewall_r_sunwall_dif(l) - sunwall_a_dif(l) = (1._r8-alb_wall_dif(l,ib)) * stot(l) - sunwall_r_dif(l) = alb_wall_dif(l,ib) * stot(l) - - stot(l) = road_r_shadewall_dif(l)/canyon_hwr(l) + sunwall_r_shadewall_dif(l) - shadewall_a_dif(l) = (1._r8-alb_wall_dif(l,ib)) * stot(l) - shadewall_r_dif(l) = alb_wall_dif(l,ib) * stot(l) - - ! step (2) - - sabs_improad_dif(l,ib) = sabs_improad_dif(l,ib) + improad_a_dif(l) - sabs_perroad_dif(l,ib) = sabs_perroad_dif(l,ib) + perroad_a_dif(l) - sabs_sunwall_dif(l,ib) = sabs_sunwall_dif(l,ib) + sunwall_a_dif(l) - sabs_shadewall_dif(l,ib) = sabs_shadewall_dif(l,ib) + shadewall_a_dif(l) - - ! step (3) - - improad_r_sky_dif(l) = improad_r_dif(l) * vf_sr(l) - improad_r_sunwall_dif(l) = improad_r_dif(l) * vf_wr(l) - improad_r_shadewall_dif(l) = improad_r_dif(l) * vf_wr(l) - - perroad_r_sky_dif(l) = perroad_r_dif(l) * vf_sr(l) - perroad_r_sunwall_dif(l) = perroad_r_dif(l) * vf_wr(l) - perroad_r_shadewall_dif(l) = perroad_r_dif(l) * vf_wr(l) - - road_r_sky_dif(l) = road_r_dif(l) * vf_sr(l) - road_r_sunwall_dif(l) = road_r_dif(l) * vf_wr(l) - road_r_shadewall_dif(l) = road_r_dif(l) * vf_wr(l) - - sunwall_r_sky_dif(l) = sunwall_r_dif(l) * vf_sw(l) - sunwall_r_road_dif(l) = sunwall_r_dif(l) * vf_rw(l) - sunwall_r_shadewall_dif(l) = sunwall_r_dif(l) * vf_ww(l) - - shadewall_r_sky_dif(l) = shadewall_r_dif(l) * vf_sw(l) - shadewall_r_road_dif(l) = shadewall_r_dif(l) * vf_rw(l) - shadewall_r_sunwall_dif(l) = shadewall_r_dif(l) * vf_ww(l) - - ! step (4) - - sref_improad_dif(l,ib) = sref_improad_dif(l,ib) + improad_r_sky_dif(l) - sref_perroad_dif(l,ib) = sref_perroad_dif(l,ib) + perroad_r_sky_dif(l) - sref_sunwall_dif(l,ib) = sref_sunwall_dif(l,ib) + sunwall_r_sky_dif(l) - sref_shadewall_dif(l,ib) = sref_shadewall_dif(l,ib) + shadewall_r_sky_dif(l) - - ! step (5) - - crit = max(road_a_dif(l), sunwall_a_dif(l), shadewall_a_dif(l)) - 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(subgrid_index=l, subgrid_level=subgrid_level_landunit, msg=errmsg(sourcefile, __LINE__)) - endif - - ! total reflected by canyon - sum of solar reflection to sky from canyon. - ! project wall fluxes to horizontal surface - - sref_canyon_dir(l) = 0.0_r8 - sref_canyon_dif(l) = 0.0_r8 - sref_canyon_dir(l) = sref_canyon_dir(l) + sref_improad_dir(l,ib)*wtroad_imperv(l) - sref_canyon_dif(l) = sref_canyon_dif(l) + sref_improad_dif(l,ib)*wtroad_imperv(l) - sref_canyon_dir(l) = sref_canyon_dir(l) + sref_perroad_dir(l,ib)*wtroad_perv(l) - sref_canyon_dif(l) = sref_canyon_dif(l) + sref_perroad_dif(l,ib)*wtroad_perv(l) - sref_canyon_dir(l) = sref_canyon_dir(l) + (sref_sunwall_dir(l,ib) + sref_shadewall_dir(l,ib))*canyon_hwr(l) - sref_canyon_dif(l) = sref_canyon_dif(l) + (sref_sunwall_dif(l,ib) + sref_shadewall_dif(l,ib))*canyon_hwr(l) - - ! total absorbed by canyon. project wall fluxes to horizontal surface - - sabs_canyon_dir(l) = 0.0_r8 - sabs_canyon_dif(l) = 0.0_r8 - sabs_canyon_dir(l) = sabs_canyon_dir(l) + sabs_improad_dir(l,ib)*wtroad_imperv(l) - sabs_canyon_dif(l) = sabs_canyon_dif(l) + sabs_improad_dif(l,ib)*wtroad_imperv(l) - sabs_canyon_dir(l) = sabs_canyon_dir(l) + sabs_perroad_dir(l,ib)*wtroad_perv(l) - sabs_canyon_dif(l) = sabs_canyon_dif(l) + sabs_perroad_dif(l,ib)*wtroad_perv(l) - sabs_canyon_dir(l) = sabs_canyon_dir(l) + (sabs_sunwall_dir(l,ib) + sabs_shadewall_dir(l,ib))*canyon_hwr(l) - sabs_canyon_dif(l) = sabs_canyon_dif(l) + (sabs_sunwall_dif(l,ib) + sabs_shadewall_dif(l,ib))*canyon_hwr(l) - - ! 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(l) = sdir_road(l,ib) + (sdir_sunwall(l,ib) + sdir_shadewall(l,ib))*canyon_hwr(l) - stot_dif(l) = sdif_road(l,ib) + (sdif_sunwall(l,ib) + sdif_shadewall(l,ib))*canyon_hwr(l) - - err = stot_dir(l) + stot_dif(l) & - - (sabs_canyon_dir(l) + sabs_canyon_dif(l) + sref_canyon_dir(l) + sref_canyon_dif(l)) - 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(l) - write(iulog,*)' stot_dif = ',stot_dif(l) - write(iulog,*)' sabs_canyon_dir = ',sabs_canyon_dir(l) - write(iulog,*)' sabs_canyon_dif = ',sabs_canyon_dif(l) - write(iulog,*)' sref_canyon_dir = ',sref_canyon_dir(l) - write(iulog,*)' sref_canyon_dif = ',sref_canyon_dir(l) - write(iulog,*) 'clm model is stopping' - call endrun(subgrid_index=l, subgrid_level=subgrid_level_landunit, msg=errmsg(sourcefile, __LINE__)) - endif - - ! canyon albedo - - canyon_alb_dif(l) = sref_canyon_dif(l) / max(stot_dif(l), 1.e-06_r8) - canyon_alb_dir(l) = sref_canyon_dir(l) / max(stot_dir(l), 1.e-06_r8) + canyon_alb_dif(l) = sref_canyon_dif(l) / max(stot_dif(l), 1.e-06_r8) + canyon_alb_dir(l) = sref_canyon_dir(l) / max(stot_dir(l), 1.e-06_r8) end do ! end of landunit loop - ! Refected and absorbed solar radiation per unit incident radiation for roof + ! Reflected and absorbed solar radiation per unit incident radiation for roof do fl = 1,num_urbanl_coszen_gt0 l = filter_urbanl_coszen_gt0(fl)