diff --git a/cime_config/testdefs/ExpectedTestFails.xml b/cime_config/testdefs/ExpectedTestFails.xml index 9735c02f4c..f5a9e3dce4 100644 --- a/cime_config/testdefs/ExpectedTestFails.xml +++ b/cime_config/testdefs/ExpectedTestFails.xml @@ -308,20 +308,6 @@ - - - FAIL - #2325 - - - - - - FAIL - #2325 - - - FAIL diff --git a/src/biogeophys/BalanceCheckMod.F90 b/src/biogeophys/BalanceCheckMod.F90 index b3efe6e525..65d2a1f81a 100644 --- a/src/biogeophys/BalanceCheckMod.F90 +++ b/src/biogeophys/BalanceCheckMod.F90 @@ -956,6 +956,7 @@ subroutine BalanceCheck( bounds, & indexp = maxloc( abs(errsol(bounds%begp:bounds%endp)), 1 , mask = (errsol(bounds%begp:bounds%endp) /= spval) ) + bounds%begp -1 indexg = patch%gridcell(indexp) + indexc = patch%column(indexp) write(iulog,*)'WARNING:: BalanceCheck, solar radiation balance error (W/m2)' write(iulog,*)'nstep = ',nstep write(iulog,*)'errsol = ',errsol(indexp) @@ -969,7 +970,18 @@ subroutine BalanceCheck( bounds, & write(iulog,*)'forc_solai(1) = ',forc_solai(indexg,1) write(iulog,*)'forc_solai(2) = ',forc_solai(indexg,2) write(iulog,*)'forc_tot = ',forc_solad(indexg,1)+forc_solad(indexg,2) & - +forc_solai(indexg,1)+forc_solai(indexg,2) + +forc_solai(indexg,1)+forc_solai(indexg,2) + + write(iulog,*)'coszen_col:',surfalb_inst%coszen_col(indexc) + write(iulog,*)'fabd:',fabd(indexp,:) + write(iulog,*)'fabi:',fabi(indexp,:) + write(iulog,*)'albd:',albd(indexp,:) + write(iulog,*)'albi:',albi(indexp,:) + write(iulog,*)'ftdd:',ftdd(indexp,:) + write(iulog,*)'ftid:',ftid(indexp,:) + write(iulog,*)'ftii:',ftii(indexp,:) + + write(iulog,*)'CTSM is stopping' call endrun(subgrid_index=indexp, subgrid_level=subgrid_level_patch, msg=errmsg(sourcefile, __LINE__)) end if diff --git a/src/biogeophys/SurfaceAlbedoMod.F90 b/src/biogeophys/SurfaceAlbedoMod.F90 index 6628f0fa4d..e8216aafdf 100644 --- a/src/biogeophys/SurfaceAlbedoMod.F90 +++ b/src/biogeophys/SurfaceAlbedoMod.F90 @@ -35,6 +35,8 @@ module SurfaceAlbedoMod public :: SurfaceAlbedo_readnl public :: SurfaceAlbedoInitTimeConst public :: SurfaceAlbedo ! Surface albedo and two-stream fluxes + public :: UpdateZenithAngles + ! ! !PRIVATE MEMBER FUNCTIONS: private :: SoilAlbedo ! Determine ground surface albedo @@ -224,7 +226,57 @@ subroutine SurfaceAlbedoInitTimeConst(bounds) end subroutine SurfaceAlbedoInitTimeConst + ! ----------------------------------------------------------------------- + + subroutine UpdateZenithAngles(bounds, surfalb_inst, nextsw_cday, declinp1) + + ! Incorporate surface slopes to generate column level zenith angles + + use clm_varctl , only : downscale_hillslope_meteorology + use shr_orb_mod + + type(bounds_type) , intent(in) :: bounds ! bounds + type(surfalb_type) , intent(inout) :: surfalb_inst + real(r8) , intent(in) :: nextsw_cday ! calendar day at Greenwich (1.00, ..., days/year) + real(r8) , intent(in) :: declinp1 ! declination angle (radians) for next time step + + integer :: c,g ! indices + + associate( & + azsun_grc => surfalb_inst%azsun_grc , & ! Output: [real(r8) (:) ] cosine of solar zenith angle + coszen_grc => surfalb_inst%coszen_grc , & ! Output: [real(r8) (:) ] cosine of solar zenith angle + coszen_col => surfalb_inst%coszen_col) ! Output: [real(r8) (:) ] cosine of solar zenith angle + + + ! Cosine solar zenith angle for next time step + ! First calculate grid-scale zenith based on time and location, no topographic effects + do g = bounds%begg,bounds%endg + coszen_grc(g) = shr_orb_cosz (nextsw_cday, grc%lat(g), grc%lon(g), declinp1) + azsun_grc(g) = shr_orb_azimuth(nextsw_cday, grc%lat(g), grc%lon(g), declinp1, acos(coszen_grc(g))) + end do + + ! calculate local incidence angle based on column slope and aspect + do c = bounds%begc,bounds%endc + g = col%gridcell(c) + if (col%is_hillslope_column(c) .and. downscale_hillslope_meteorology) then + + ! hill_slope is [m/m], convert to radians with atan function + coszen_col(c) = shr_orb_cosinc(acos(coszen_grc(g)),azsun_grc(g),atan(col%hill_slope(c)),col%hill_aspect(c)) + + if(coszen_grc(g) > 0._r8 .and. coszen_col(c) < 0._r8) coszen_col(c) = 0._r8 + else + coszen_col(c) = coszen_grc(g) + endif + + end do + + end associate + + return + end subroutine UpdateZenithAngles + !----------------------------------------------------------------------- + subroutine SurfaceAlbedo(bounds,nc, & num_nourbanc, filter_nourbanc, & num_nourbanp, filter_nourbanp, & @@ -256,13 +308,13 @@ subroutine SurfaceAlbedo(bounds,nc, & ! only computed over active points. ! ! !USES: - use shr_orb_mod + !use shr_orb_mod use clm_time_manager , only : get_nstep use abortutils , only : endrun use clm_varctl , only : use_subgrid_fluxes, use_snicar_frc, use_fates use CLMFatesInterfaceMod, only : hlm_fates_interface_type use landunit_varcon , only : istsoil - use clm_varctl , only : downscale_hillslope_meteorology + ! !ARGUMENTS: type(bounds_type) , intent(in) :: bounds ! bounds @@ -427,32 +479,15 @@ subroutine SurfaceAlbedo(bounds,nc, & fabi_sha_z => surfalb_inst%fabi_sha_z_patch & ! Output: [real(r8) (:,:) ] absorbed shaded leaf diffuse PAR (per unit lai+sai) for each canopy layer ) - ! Cosine solar zenith angle for next time step - - do g = bounds%begg,bounds%endg - coszen_grc(g) = shr_orb_cosz (nextsw_cday, grc%lat(g), grc%lon(g), declinp1) - end do - - do c = bounds%begc,bounds%endc - g = col%gridcell(c) - if (col%is_hillslope_column(c) .and. downscale_hillslope_meteorology) then - ! calculate local incidence angle based on column slope and aspect - zenith_angle = acos(coszen_grc(g)) - - azsun_grc(g) = shr_orb_azimuth(nextsw_cday, grc%lat(g), grc%lon(g), declinp1, zenith_angle) - ! hill_slope is [m/m], convert to radians - coszen_col(c) = shr_orb_cosinc(zenith_angle,azsun_grc(g),atan(col%hill_slope(c)),col%hill_aspect(c)) - - if(coszen_grc(g) > 0._r8 .and. coszen_col(c) < 0._r8) coszen_col(c) = 0._r8 - - else - coszen_col(c) = coszen_grc(g) - endif - end do + ! Apply column level zenith angles to the patch level do fp = 1,num_nourbanp p = filter_nourbanp(fp) c = patch%column(p) - coszen_patch(p) = coszen_col(c) + if(patch%is_fates(p))then + coszen_patch(p) = spval + else + coszen_patch(p) = coszen_col(c) + end if end do ! Initialize output because solar radiation only done if coszen > 0 @@ -1045,10 +1080,7 @@ subroutine SurfaceAlbedo(bounds,nc, & if (use_fates) then - call clm_fates%wrap_canopy_radiation(bounds, nc, & - num_vegsol, filter_vegsol, & - coszen_patch(bounds%begp:bounds%endp), & - fcansno(bounds%begp:bounds%endp), surfalb_inst) + call clm_fates%wrap_canopy_radiation(bounds, nc, fcansno(bounds%begp:bounds%endp), surfalb_inst) else diff --git a/src/main/clm_driver.F90 b/src/main/clm_driver.F90 index f93143d9e3..38ff291f6c 100644 --- a/src/main/clm_driver.F90 +++ b/src/main/clm_driver.F90 @@ -52,6 +52,7 @@ module clm_driver use AerosolMod , only : AerosolMasses use SnowSnicarMod , only : SnowAge_grain use SurfaceAlbedoMod , only : SurfaceAlbedo + use SurfaceAlbedoMod , only : UpdateZenithAngles use UrbanAlbedoMod , only : UrbanAlbedo ! use SurfaceRadiationMod , only : SurfaceRadiation, CanopySunShadeFracs @@ -1224,7 +1225,15 @@ subroutine clm_drv(doalb, nextsw_cday, declinp1, declin, rstwr, nlend, rdate, ro ! Determine albedos for next time step ! ============================================================================ - if (doalb) then + if (doalb .or. use_fates) then + ! During branch runs and non continue_run restarts, the doalb flag + ! does not trigger correctly for fates runs (and non-fates?), and thus + ! the zenith angles are not calculated and ready when radiation scattering + ! needs to occur. + call UpdateZenithAngles(bounds_clump,surfalb_inst, nextsw_cday, declinp1) + end if + + if (doalb ) then ! Albedos for non-urban columns call t_startf('surfalb') diff --git a/src/utils/clmfates_interfaceMod.F90 b/src/utils/clmfates_interfaceMod.F90 index 269189d1b7..4663366ffe 100644 --- a/src/utils/clmfates_interfaceMod.F90 +++ b/src/utils/clmfates_interfaceMod.F90 @@ -2752,28 +2752,31 @@ end subroutine wrap_WoodProducts ! ====================================================================================== - subroutine wrap_canopy_radiation(this, bounds_clump, nc, & - num_vegsol, filter_vegsol, coszen, fcansno, surfalb_inst) + subroutine wrap_canopy_radiation(this, bounds_clump, nc, fcansno, surfalb_inst) + + ! Pass boundary conditions (zenith angle, ground albedo and snow-cover) + ! to FATES, and have fates call normalized radiation scattering. + ! These methods are normalized because it is the zenith of the next + ! timestep, and because of atmospheric model coupling, we need to provide + ! an albedo. This normalized solution will be scaled by the actual + ! downwelling radiation during the wrap sunshade fraction call + ! Arguments class(hlm_fates_interface_type), intent(inout) :: this type(bounds_type), intent(in) :: bounds_clump - ! filter for vegetated pfts with coszen>0 integer , intent(in) :: nc ! clump index - integer , intent(in) :: num_vegsol - integer , intent(in) :: filter_vegsol(num_vegsol) - ! cosine solar zenith angle for next time step - real(r8) , intent(in) :: coszen( bounds_clump%begp: ) real(r8) , intent(in) :: fcansno( bounds_clump%begp: ) type(surfalb_type) , intent(inout) :: surfalb_inst ! locals - integer :: s,c,p,ifp,icp + integer :: s,c,p,ifp,g call t_startf('fates_wrapcanopyradiation') associate(& + coszen_col => surfalb_inst%coszen_col , & !in albgrd_col => surfalb_inst%albgrd_col , & !in albgri_col => surfalb_inst%albgri_col , & !in albd => surfalb_inst%albd_patch , & !out @@ -2787,46 +2790,50 @@ subroutine wrap_canopy_radiation(this, bounds_clump, nc, & do s = 1, this%fates(nc)%nsites c = this%f2hmap(nc)%fcolumn(s) + g = col%gridcell(c) - do ifp = 1,this%fates(nc)%sites(s)%youngest_patch%patchno - p = ifp+col%patchi(c) + this%fates(nc)%bc_in(s)%coszen = coszen_col(c) - if( any(filter_vegsol==p) )then + ! FATES only does short-wave radiation scattering when + ! the zenith angle is positive. In the future, FATES + ! may perform calculations on diffuse radiation during + ! dawn and dusk when the sun is below the horizon, but not yet - this%fates(nc)%bc_in(s)%filter_vegzen_pa(ifp) = .true. - this%fates(nc)%bc_in(s)%coszen_pa(ifp) = coszen(p) - this%fates(nc)%bc_in(s)%fcansno_pa(ifp) = fcansno(p) - this%fates(nc)%bc_in(s)%albgr_dir_rb(:) = albgrd_col(c,:) - this%fates(nc)%bc_in(s)%albgr_dif_rb(:) = albgri_col(c,:) + do ifp = 1,this%fates(nc)%sites(s)%youngest_patch%patchno + p = ifp+col%patchi(c) + this%fates(nc)%bc_in(s)%fcansno_pa(ifp) = fcansno(p) + end do - else + + if(coszen_col(c) > 0._r8) then - this%fates(nc)%bc_in(s)%filter_vegzen_pa(ifp) = .false. + this%fates(nc)%bc_in(s)%albgr_dir_rb(:) = albgrd_col(c,:) + this%fates(nc)%bc_in(s)%albgr_dif_rb(:) = albgri_col(c,:) + else - end if + ! This will ensure a crash in FATES if it tries + this%fates(nc)%bc_in(s)%albgr_dir_rb(:) = spval + this%fates(nc)%bc_in(s)%albgr_dif_rb(:) = spval + + end if - end do end do - call FatesNormalizedCanopyRadiation(this%fates(nc)%nsites, & - this%fates(nc)%sites, & - this%fates(nc)%bc_in, & - this%fates(nc)%bc_out) + call FatesNormalizedCanopyRadiation( & + this%fates(nc)%sites, & + this%fates(nc)%bc_in, & + this%fates(nc)%bc_out) ! Pass FATES BC's back to HLM ! ----------------------------------------------------------------------------------- - do icp = 1,num_vegsol - p = filter_vegsol(icp) - c = patch%column(p) - s = this%f2hmap(nc)%hsites(c) - ! do if structure here and only pass natveg columns - ifp = p-col%patchi(c) + do s = 1, this%fates(nc)%nsites + + c = this%f2hmap(nc)%fcolumn(s) + + do ifp = 1, this%fates(nc)%sites(s)%youngest_patch%patchno + + p = ifp+col%patchi(c) - if(.not.this%fates(nc)%bc_in(s)%filter_vegzen_pa(ifp) )then - write(iulog,*) 's,p,ifp',s,p,ifp - write(iulog,*) 'Not all patches on the natveg column were passed to canrad',patch%sp_pftorder_index(p) -! call endrun(msg=errMsg(sourcefile, __LINE__)) - else albd(p,:) = this%fates(nc)%bc_out(s)%albd_parb(ifp,:) albi(p,:) = this%fates(nc)%bc_out(s)%albi_parb(ifp,:) fabd(p,:) = this%fates(nc)%bc_out(s)%fabd_parb(ifp,:) @@ -2834,9 +2841,10 @@ subroutine wrap_canopy_radiation(this, bounds_clump, nc, & ftdd(p,:) = this%fates(nc)%bc_out(s)%ftdd_parb(ifp,:) ftid(p,:) = this%fates(nc)%bc_out(s)%ftid_parb(ifp,:) ftii(p,:) = this%fates(nc)%bc_out(s)%ftii_parb(ifp,:) - end if + + end do end do - + end associate call t_stopf('fates_wrapcanopyradiation')