diff --git a/ED/src/driver/ed_met_driver.f90 b/ED/src/driver/ed_met_driver.f90 index e6fb7449a..32d419586 100644 --- a/ED/src/driver/ed_met_driver.f90 +++ b/ED/src/driver/ed_met_driver.f90 @@ -861,7 +861,8 @@ subroutine update_met_drivers(cgrid) , current_time ! ! intent(in) use canopy_air_coms , only : ubmin & ! intent(in) , ustmin ! ! intent(in) - use canopy_radiation_coms, only : cosz_min ! ! intent(in) + use canopy_radiation_coms, only : cosz_min & ! intent(in) + , find_eff_cosz ! ! intent(in) use consts_coms , only : day_sec & ! intent(in) , t00 & ! intent(in) , t3ple & ! intent(in) @@ -882,7 +883,7 @@ subroutine update_met_drivers(cgrid) use lapse , only : calc_met_lapse ! ! sub-routine use update_derived_utils , only : update_model_time_dm ! ! sub-routine use radiate_utils , only : ed_zen & ! function - , mean_daysecz & ! function + , mean_chapman & ! function , solar_radiation_breakdown & ! sub-routine , dump_radinfo ! ! sub-routine @@ -915,8 +916,8 @@ subroutine update_met_drivers(cgrid) real(kind=4) :: snden ! snow density (kg/m3) real(kind=4) :: fice ! Ice fraction precipication real(kind=4) :: fliq ! Liquid fraction precipitation - real(kind=4) :: secz_prev ! Mean of sec(zenith angle) - previous - real(kind=4) :: secz_next ! Mean of sec(zenith angle) - next + real(kind=4) :: chapman_prev ! Mean of Chapman function - previous + real(kind=4) :: chapman_next ! Mean of Chapman function - next real(kind=4) :: fperp_prev ! Perpendicular flux - previous real(kind=4) :: fperp_next ! Perpendicular flux - next logical :: night_next @@ -932,7 +933,8 @@ subroutine update_met_drivers(cgrid) !------------------------------------------------------------------------------------! do ipy=1,cgrid%npolygons cgrid%met(ipy)%vels = 0.0 - cgrid%cosz(ipy) = ed_zen(cgrid%lon(ipy),cgrid%lat(ipy),current_time) + cgrid%cosz (ipy) = ed_zen(cgrid%lon(ipy),cgrid%lat(ipy),current_time) + cgrid%eff_cosz(ipy) = find_eff_cosz(cgrid%cosz(ipy)) end do !------------------------------------------------------------------------------------! @@ -956,9 +958,9 @@ subroutine update_met_drivers(cgrid) !------------------------------------------------------------------------! mprev = 1 prevmet_timea = current_time + !------------------------------------------------------------------------! case (0,3) - !------------------------------------------------------------------------! ! Time dependent variables, but with no time interpolation. ! !------------------------------------------------------------------------! @@ -967,6 +969,7 @@ subroutine update_met_drivers(cgrid) mprev = int(float(nseconds_elapsed) / met_frq(iformat,iv)) + 1 dtprev = floor(real(nseconds_elapsed)/met_frq(iformat,iv)) & * met_frq(iformat,iv) - real(nseconds_elapsed) + !------------------------------------------------------------------------! @@ -1023,10 +1026,12 @@ subroutine update_met_drivers(cgrid) else do ipy = 1,cgrid%npolygons !------------------------------------------------------------------! - ! Check whether this is day time or night time. We should ! - ! only do the full interpolation thing during the day, the night ! - ! should have no incoming radiation until we incorporate fireflies ! - ! (or a good twilight scheme) in the model. ! + ! Check whether this is daytime or night/twilight time. We ! + ! should carry out the sun-angle dependent interpolation only when ! + ! there is direct light. Although we check for daytime using the ! + ! actual cosine of the zenith angle, we interpolate irradiance ! + ! using the effective cosine of zenith angle, which prevents near- ! + ! singularities close to sunset times. ! !------------------------------------------------------------------! if (cgrid%cosz(ipy) > cosz_min) then !---------------------------------------------------------------! @@ -1035,12 +1040,13 @@ subroutine update_met_drivers(cgrid) !---------------------------------------------------------------! select case (imetavg) case (0) - secz_prev = mean_daysecz(cgrid%lon(ipy),cgrid%lat(ipy) & - ,prevmet_timea,dt_radinterp,0.) + chapman_prev = & + mean_chapman(cgrid%lon(ipy),cgrid%lat(ipy),prevmet_timea & + ,dt_radinterp,0.,.false.) case default - secz_prev = mean_daysecz(cgrid%lon(ipy),cgrid%lat(ipy) & - ,prevmet_timea,dt_radinterp & - ,met_frq(iformat,iv)) + chapman_prev = & + mean_chapman(cgrid%lon(ipy),cgrid%lat(ipy),prevmet_timea & + ,dt_radinterp,met_frq(iformat,iv),.true.) end select !---------------------------------------------------------------! @@ -1050,7 +1056,7 @@ subroutine update_met_drivers(cgrid) ! Decide whether we want to use the previous or the next ! ! data based on both the zenith angle. ! !---------------------------------------------------------------! - night_prev = secz_prev == 0. + night_prev = chapman_prev == 0. !---------------------------------------------------------------! @@ -1072,15 +1078,15 @@ subroutine update_met_drivers(cgrid) ! secant and scaling back with the current cosine of zenith ! ! angle. ! !------------------------------------------------------------! - fperp_prev = cgrid%metinput(ipy)%nbdsf(mprev) * secz_prev - cgrid%met(ipy)%nir_beam = fperp_prev * cgrid%cosz(ipy) + fperp_prev = cgrid%metinput(ipy)%nbdsf(mprev) * chapman_prev + cgrid%met(ipy)%nir_beam = fperp_prev * cgrid%eff_cosz(ipy) end if !---------------------------------------------------------------! !----- Assume all future stuff to be zero, we won't use them. --! - secz_next = 0. + chapman_next = 0. night_next = .true. fperp_next = 0. wprev = 1.0 @@ -1090,8 +1096,8 @@ subroutine update_met_drivers(cgrid) else !----- Night time, assign it zero. -----------------------------! - secz_prev = 0. - secz_next = 0. + chapman_prev = 0. + chapman_next = 0. night_prev = .true. night_next = .true. fperp_prev = 0. @@ -1112,8 +1118,8 @@ subroutine update_met_drivers(cgrid) if (print_radinterp) then call dump_radinfo(cgrid,ipy,mprev,mnext,prevmet_timea & ,nextmet_timea,met_frq(iformat,iv),wprev,wnext & - ,secz_prev,secz_next,fperp_prev,fperp_next & - ,night_prev,night_next,'nbdsf') + ,chapman_prev,chapman_next,fperp_prev & + ,fperp_next,night_prev,night_next,'nbdsf') end if !------------------------------------------------------------------! end do @@ -1133,26 +1139,27 @@ subroutine update_met_drivers(cgrid) end do else do ipy = 1,cgrid%npolygons - !------------------------------------------------------------------! - ! Check whether this is day time or night time. We should ! - ! only do the full interpolation thing during the day, the night ! - ! should have no incoming radiation until we incorporate fireflies ! - ! (or a good twilight scheme) in the model. ! + ! Check whether this is day/twilight time or night time. We ! + ! should carry out the sun-angle dependent interpolation only when ! + ! there is diffuse light. Interpolation is done with the ! + ! effective cosine of zenith angle, which should be safe for ! + ! interpolations as long as cosz_min is not way too close to zero. ! !------------------------------------------------------------------! - if (cgrid%cosz(ipy) > cosz_min) then + if (cgrid%eff_cosz(ipy) > cosz_min) then !---------------------------------------------------------------! ! Define the normalisation factors for the previous and the ! ! next time. ! !---------------------------------------------------------------! select case (imetavg) case (0) - secz_prev = mean_daysecz(cgrid%lon(ipy),cgrid%lat(ipy) & - ,prevmet_timea,dt_radinterp,0.) + chapman_prev = & + mean_chapman(cgrid%lon(ipy),cgrid%lat(ipy),prevmet_timea & + ,dt_radinterp,0.,.true.) case default - secz_prev = mean_daysecz(cgrid%lon(ipy),cgrid%lat(ipy) & - ,prevmet_timea,dt_radinterp & - ,met_frq(iformat,iv)) + chapman_prev = & + mean_chapman(cgrid%lon(ipy),cgrid%lat(ipy),prevmet_timea & + ,dt_radinterp,met_frq(iformat,iv),.true.) end select !---------------------------------------------------------------! @@ -1162,7 +1169,7 @@ subroutine update_met_drivers(cgrid) ! Decide whether we want to use the previous or the next ! ! data based on both the zenith angle. ! !---------------------------------------------------------------! - night_prev = secz_prev == 0. + night_prev = chapman_prev == 0. !---------------------------------------------------------------! @@ -1184,15 +1191,15 @@ subroutine update_met_drivers(cgrid) ! secant and scaling back with the current cosine of zenith ! ! angle. ! !------------------------------------------------------------! - fperp_prev = cgrid%metinput(ipy)%nddsf(mprev) * secz_prev - cgrid%met(ipy)%nir_diffuse = fperp_prev * cgrid%cosz(ipy) + fperp_prev = cgrid%metinput(ipy)%nddsf(mprev) * chapman_prev + cgrid%met(ipy)%nir_diffuse = fperp_prev * cgrid%eff_cosz(ipy) end if !---------------------------------------------------------------! !----- Assume all future stuff to be zero, we won't use them. --! - secz_next = 0. + chapman_next = 0. night_next = .true. fperp_next = 0. wprev = 1.0 @@ -1202,8 +1209,8 @@ subroutine update_met_drivers(cgrid) else !----- Night time, assign it zero. -----------------------------! - secz_prev = 0. - secz_next = 0. + chapman_prev = 0. + chapman_next = 0. night_prev = .true. night_next = .true. fperp_prev = 0. @@ -1224,8 +1231,8 @@ subroutine update_met_drivers(cgrid) if (print_radinterp) then call dump_radinfo(cgrid,ipy,mprev,mnext,prevmet_timea & ,nextmet_timea,met_frq(iformat,iv),wprev,wnext & - ,secz_prev,secz_next,fperp_prev,fperp_next & - ,night_prev,night_next,'nddsf') + ,chapman_prev,chapman_next,fperp_prev & + ,fperp_next,night_prev,night_next,'nddsf') end if !------------------------------------------------------------------! end do @@ -1246,10 +1253,12 @@ subroutine update_met_drivers(cgrid) do ipy = 1,cgrid%npolygons !------------------------------------------------------------------! - ! Check whether this is day time or night time. We should ! - ! only do the full interpolation thing during the day, the night ! - ! should have no incoming radiation until we incorporate fireflies ! - ! (or a good twilight scheme) in the model. ! + ! Check whether this is daytime or night/twilight time. We ! + ! should carry out the sun-angle dependent interpolation only when ! + ! there is direct light. Although we check for daytime using the ! + ! actual cosine of the zenith angle, we interpolate irradiance ! + ! using the effective cosine of zenith angle, which prevents near- ! + ! singularities close to sunset times. ! !------------------------------------------------------------------! if (cgrid%cosz(ipy) > cosz_min) then !---------------------------------------------------------------! @@ -1258,12 +1267,13 @@ subroutine update_met_drivers(cgrid) !---------------------------------------------------------------! select case (imetavg) case (0) - secz_prev = mean_daysecz(cgrid%lon(ipy),cgrid%lat(ipy) & - ,prevmet_timea,dt_radinterp,0.) + chapman_prev = & + mean_chapman(cgrid%lon(ipy),cgrid%lat(ipy),prevmet_timea & + ,dt_radinterp,0.,.false.) case default - secz_prev = mean_daysecz(cgrid%lon(ipy),cgrid%lat(ipy) & - ,prevmet_timea,dt_radinterp & - ,met_frq(iformat,iv)) + chapman_prev = & + mean_chapman(cgrid%lon(ipy),cgrid%lat(ipy),prevmet_timea & + ,dt_radinterp,met_frq(iformat,iv),.false.) end select !---------------------------------------------------------------! @@ -1273,7 +1283,7 @@ subroutine update_met_drivers(cgrid) ! Decide whether we want to use the previous or the next ! ! data based on both the zenith angle. ! !---------------------------------------------------------------! - night_prev = secz_prev == 0. + night_prev = chapman_prev == 0. !---------------------------------------------------------------! @@ -1295,15 +1305,15 @@ subroutine update_met_drivers(cgrid) ! secant and scaling back with the current cosine of zenith ! ! angle. ! !------------------------------------------------------------! - fperp_prev = cgrid%metinput(ipy)%vbdsf(mprev) * secz_prev - cgrid%met(ipy)%par_beam = fperp_prev * cgrid%cosz(ipy) + fperp_prev = cgrid%metinput(ipy)%vbdsf(mprev) * chapman_prev + cgrid%met(ipy)%par_beam = fperp_prev * cgrid%eff_cosz(ipy) end if !---------------------------------------------------------------! !----- Assume all future stuff to be zero, we won't use them. --! - secz_next = 0. + chapman_next = 0. night_next = .true. fperp_next = 0. wprev = 1.0 @@ -1313,8 +1323,8 @@ subroutine update_met_drivers(cgrid) else !----- Night time, assign it zero. -----------------------------! - secz_prev = 0. - secz_next = 0. + chapman_prev = 0. + chapman_next = 0. night_prev = .true. night_next = .true. fperp_prev = 0. @@ -1335,8 +1345,8 @@ subroutine update_met_drivers(cgrid) if (print_radinterp) then call dump_radinfo(cgrid,ipy,mprev,mnext,prevmet_timea & ,nextmet_timea,met_frq(iformat,iv),wprev,wnext & - ,secz_prev,secz_next,fperp_prev,fperp_next & - ,night_prev,night_next,'vbdsf') + ,chapman_prev,chapman_next,fperp_prev & + ,fperp_next,night_prev,night_next,'vbdsf') end if !------------------------------------------------------------------! end do @@ -1358,24 +1368,26 @@ subroutine update_met_drivers(cgrid) do ipy = 1,cgrid%npolygons !------------------------------------------------------------------! - ! Check whether this is day time or night time. We should ! - ! only do the full interpolation thing during the day, the night ! - ! should have no incoming radiation until we incorporate fireflies ! - ! (or a good twilight scheme) in the model. ! + ! Check whether this is day/twilight time or night time. We ! + ! should carry out the sun-angle dependent interpolation only when ! + ! there is diffuse light. Interpolation is done with the ! + ! effective cosine of zenith angle, which should be safe for ! + ! interpolations as long as cosz_min is not way too close to zero. ! !------------------------------------------------------------------! - if (cgrid%cosz(ipy) > cosz_min) then + if (cgrid%eff_cosz(ipy) > cosz_min) then !---------------------------------------------------------------! ! Define the normalisation factors for the previous and the ! ! next time. ! !---------------------------------------------------------------! select case (imetavg) case (0) - secz_prev = mean_daysecz(cgrid%lon(ipy),cgrid%lat(ipy) & - ,prevmet_timea,dt_radinterp,0.) + chapman_prev = & + mean_chapman(cgrid%lon(ipy),cgrid%lat(ipy),prevmet_timea & + ,dt_radinterp,0.,.true.) case default - secz_prev = mean_daysecz(cgrid%lon(ipy),cgrid%lat(ipy) & - ,prevmet_timea,dt_radinterp & - ,met_frq(iformat,iv)) + chapman_prev = & + mean_chapman(cgrid%lon(ipy),cgrid%lat(ipy),prevmet_timea & + ,dt_radinterp,met_frq(iformat,iv)) end select !---------------------------------------------------------------! @@ -1385,7 +1397,7 @@ subroutine update_met_drivers(cgrid) ! Decide whether we want to use the previous or the next ! ! data based on both the zenith angle. ! !---------------------------------------------------------------! - night_prev = secz_prev == 0. + night_prev = chapman_prev == 0. !---------------------------------------------------------------! @@ -1407,15 +1419,15 @@ subroutine update_met_drivers(cgrid) ! secant and scaling back with the current cosine of zenith ! ! angle. ! !------------------------------------------------------------! - fperp_prev = cgrid%metinput(ipy)%vddsf(mprev) * secz_prev - cgrid%met(ipy)%par_diffuse = fperp_prev * cgrid%cosz(ipy) + fperp_prev = cgrid%metinput(ipy)%vddsf(mprev) * chapman_prev + cgrid%met(ipy)%par_diffuse = fperp_prev * cgrid%eff_cosz(ipy) end if !---------------------------------------------------------------! !----- Assume all future stuff to be zero, we won't use them. --! - secz_next = 0. + chapman_next = 0. night_next = .true. fperp_next = 0. wprev = 1.0 @@ -1425,8 +1437,8 @@ subroutine update_met_drivers(cgrid) else !----- Night time, assign it zero. -----------------------------! - secz_prev = 0. - secz_next = 0. + chapman_prev = 0. + chapman_next = 0. night_prev = .true. night_next = .true. fperp_prev = 0. @@ -1447,8 +1459,8 @@ subroutine update_met_drivers(cgrid) if (print_radinterp) then call dump_radinfo(cgrid,ipy,mprev,mnext,prevmet_timea & ,nextmet_timea,met_frq(iformat,iv),wprev,wnext & - ,secz_prev,secz_next,fperp_prev,fperp_next & - ,night_prev,night_next,'vddsf') + ,chapman_prev,chapman_next,fperp_prev & + ,fperp_next,night_prev,night_next,'vddsf') end if !------------------------------------------------------------------! end do @@ -1779,10 +1791,12 @@ subroutine update_met_drivers(cgrid) do ipy= 1,cgrid%npolygons !------------------------------------------------------------------! - ! Check whether this is day time or night time. We should ! - ! only do the full interpolation thing during the day, the night ! - ! should have no incoming radiation until we incorporate fireflies ! - ! (or a good twilight scheme) in the model. ! + ! Check whether this is daytime or night/twilight time. We ! + ! should carry out the sun-angle dependent interpolation only when ! + ! there is direct light. Although we check for daytime using the ! + ! actual cosine of the zenith angle, we interpolate irradiance ! + ! using the effective cosine of zenith angle, which prevents near- ! + ! singularities close to sunset times. ! !------------------------------------------------------------------! if (cgrid%cosz(ipy) > cosz_min) then @@ -1792,17 +1806,19 @@ subroutine update_met_drivers(cgrid) !---------------------------------------------------------------! select case (imetavg) case (0) - secz_prev = mean_daysecz(cgrid%lon(ipy),cgrid%lat(ipy) & - ,prevmet_timea,dt_radinterp,0.) - secz_next = mean_daysecz(cgrid%lon(ipy),cgrid%lat(ipy) & - ,nextmet_timea,dt_radinterp,0.) + chapman_prev = & + mean_chapman(cgrid%lon(ipy),cgrid%lat(ipy),prevmet_timea & + ,dt_radinterp,0.,.false.) + chapman_next = & + mean_chapman(cgrid%lon(ipy),cgrid%lat(ipy),nextmet_timea & + ,dt_radinterp,0.,.false.) case default - secz_prev = mean_daysecz(cgrid%lon(ipy),cgrid%lat(ipy) & - ,prevmet_timea,dt_radinterp & - ,met_frq(iformat,iv)) - secz_next = mean_daysecz(cgrid%lon(ipy),cgrid%lat(ipy) & - ,nextmet_timea,dt_radinterp & - ,met_frq(iformat,iv)) + chapman_prev = & + mean_chapman(cgrid%lon(ipy),cgrid%lat(ipy),prevmet_timea & + ,dt_radinterp,met_frq(iformat,iv),.false.) + chapman_next = & + mean_chapman(cgrid%lon(ipy),cgrid%lat(ipy),nextmet_timea & + ,dt_radinterp,met_frq(iformat,iv),.false.) end select !---------------------------------------------------------------! @@ -1812,8 +1828,8 @@ subroutine update_met_drivers(cgrid) ! Decide whether we want to use the previous or the next ! ! data based on the zenith angle. ! !---------------------------------------------------------------! - night_prev = secz_prev == 0. - night_next = secz_next == 0. + night_prev = chapman_prev == 0. + night_next = chapman_next == 0. !---------------------------------------------------------------! @@ -1837,8 +1853,8 @@ subroutine update_met_drivers(cgrid) ! current cosine of zenith angle. ! !------------------------------------------------------------! fperp_prev = 0. - fperp_next = cgrid%metinput(ipy)%nbdsf(mnext) * secz_next - cgrid%met(ipy)%nir_beam = fperp_next * cgrid%cosz(ipy) + fperp_next = cgrid%metinput(ipy)%nbdsf(mnext) * chapman_next + cgrid%met(ipy)%nir_beam = fperp_next * cgrid%eff_cosz(ipy) elseif (night_next) then !-----------------------------------------------------------! ! Dusk time, the next time is zero so we only have ! @@ -1846,22 +1862,22 @@ subroutine update_met_drivers(cgrid) ! the next time using the previous secant and scaling back ! ! with the current cosine of zenith angle. ! !-----------------------------------------------------------! - fperp_prev = cgrid%metinput(ipy)%nbdsf(mprev) * secz_prev + fperp_prev = cgrid%metinput(ipy)%nbdsf(mprev) * chapman_prev fperp_next = 0. - cgrid%met(ipy)%nir_beam = fperp_prev * cgrid%cosz(ipy) + cgrid%met(ipy)%nir_beam = fperp_prev * cgrid%eff_cosz(ipy) else !----- Daytime, use both previous and next values. ---------! - fperp_next = cgrid%metinput(ipy)%nbdsf(mnext) * secz_next - fperp_prev = cgrid%metinput(ipy)%nbdsf(mprev) * secz_prev - cgrid%met(ipy)%nir_beam = cgrid%cosz(ipy) & + fperp_next = cgrid%metinput(ipy)%nbdsf(mnext) * chapman_next + fperp_prev = cgrid%metinput(ipy)%nbdsf(mprev) * chapman_prev + cgrid%met(ipy)%nir_beam = cgrid%eff_cosz(ipy) & * ( fperp_next * wnext & + fperp_prev * wprev ) end if !---------------------------------------------------------------! else !----- Night time, assign it zero. -----------------------------! - secz_prev = 0. - secz_next = 0. + chapman_prev = 0. + chapman_next = 0. night_prev = .true. night_next = .true. fperp_prev = 0. @@ -1880,8 +1896,8 @@ subroutine update_met_drivers(cgrid) if (print_radinterp) then call dump_radinfo(cgrid,ipy,mprev,mnext,prevmet_timea & ,nextmet_timea,met_frq(iformat,iv),wprev,wnext & - ,secz_prev,secz_next,fperp_prev,fperp_next & - ,night_prev,night_next,'nbdsf') + ,chapman_prev,chapman_next,fperp_prev & + ,fperp_next,night_prev,night_next,'nbdsf') end if !------------------------------------------------------------------! end do @@ -1904,12 +1920,13 @@ subroutine update_met_drivers(cgrid) do ipy= 1,cgrid%npolygons !------------------------------------------------------------------! - ! Check whether this is day time or night time. We should ! - ! only do the full interpolation thing during the day, the night ! - ! should have no incoming radiation until we incorporate fireflies ! - ! (or a good twilight scheme) in the model. ! + ! Check whether this is day/twilight time or night time. We ! + ! should carry out the sun-angle dependent interpolation only when ! + ! there is diffuse light. Interpolation is done with the ! + ! effective cosine of zenith angle, which should be safe for ! + ! interpolations as long as cosz_min is not way too close to zero. ! !------------------------------------------------------------------! - if (cgrid%cosz(ipy) > cosz_min) then + if (cgrid%eff_cosz(ipy) > cosz_min) then !---------------------------------------------------------------! ! Define the normalisation factors for the previous and the ! @@ -1917,17 +1934,19 @@ subroutine update_met_drivers(cgrid) !---------------------------------------------------------------! select case (imetavg) case (0) - secz_prev = mean_daysecz(cgrid%lon(ipy),cgrid%lat(ipy) & - ,prevmet_timea,dt_radinterp,0.) - secz_next = mean_daysecz(cgrid%lon(ipy),cgrid%lat(ipy) & - ,nextmet_timea,dt_radinterp,0.) + chapman_prev = & + mean_chapman(cgrid%lon(ipy),cgrid%lat(ipy),prevmet_timea & + ,dt_radinterp,0.,.true.) + chapman_next = & + mean_chapman(cgrid%lon(ipy),cgrid%lat(ipy),nextmet_timea & + ,dt_radinterp,0.,.true.) case default - secz_prev = mean_daysecz(cgrid%lon(ipy),cgrid%lat(ipy) & - ,prevmet_timea,dt_radinterp & - ,met_frq(iformat,iv)) - secz_next = mean_daysecz(cgrid%lon(ipy),cgrid%lat(ipy) & - ,nextmet_timea,dt_radinterp & - ,met_frq(iformat,iv)) + chapman_prev = & + mean_chapman(cgrid%lon(ipy),cgrid%lat(ipy),prevmet_timea & + ,dt_radinterp,met_frq(iformat,iv),.true.) + chapman_next = & + mean_chapman(cgrid%lon(ipy),cgrid%lat(ipy),nextmet_timea & + ,dt_radinterp,met_frq(iformat,iv),.true.) end select !---------------------------------------------------------------! @@ -1937,8 +1956,8 @@ subroutine update_met_drivers(cgrid) ! Decide whether we want to use the previous or the next ! ! data based on the zenith angle. ! !---------------------------------------------------------------! - night_prev = secz_prev == 0. - night_next = secz_next == 0. + night_prev = chapman_prev == 0. + night_next = chapman_next == 0. !---------------------------------------------------------------! @@ -1962,37 +1981,37 @@ subroutine update_met_drivers(cgrid) ! current cosine of zenith angle. ! !------------------------------------------------------------! fperp_prev = 0. - fperp_next = cgrid%metinput(ipy)%nddsf(mnext) * secz_next - cgrid%met(ipy)%nir_diffuse = fperp_next * cgrid%cosz(ipy) + fperp_next = cgrid%metinput(ipy)%nddsf(mnext) * chapman_next + cgrid%met(ipy)%nir_diffuse = fperp_next * cgrid%eff_cosz(ipy) elseif (night_next) then !------------------------------------------------------------! ! Dusk time, the next time is zero so we only have ! ! meaningful information for the previous time: we scale the ! ! next time using the previous secant and scaling back with ! - ! the current cosine of zenith angle. ! + ! the current effective cosine of zenith angle. ! !------------------------------------------------------------! - fperp_prev = cgrid%metinput(ipy)%nddsf(mprev) * secz_prev + fperp_prev = cgrid%metinput(ipy)%nddsf(mprev) * chapman_prev fperp_next = 0. - cgrid%met(ipy)%nir_diffuse = fperp_prev * cgrid%cosz(ipy) + cgrid%met(ipy)%nir_diffuse = fperp_prev * cgrid%eff_cosz(ipy) + !------------------------------------------------------------! else !----- Daytime, use both previous and next values. ----------! - fperp_next = cgrid%metinput(ipy)%nddsf(mnext) * secz_next - fperp_prev = cgrid%metinput(ipy)%nddsf(mprev) * secz_prev - cgrid%met(ipy)%nir_diffuse = cgrid%cosz(ipy) & + fperp_next = cgrid%metinput(ipy)%nddsf(mnext) * chapman_next + fperp_prev = cgrid%metinput(ipy)%nddsf(mprev) * chapman_prev + cgrid%met(ipy)%nir_diffuse = cgrid%eff_cosz(ipy) & * ( fperp_next * wnext & + fperp_prev * wprev ) end if !---------------------------------------------------------------! else !----- Night time, assign it zero. -----------------------------! - secz_prev = 0. - secz_next = 0. + chapman_prev = 0. + chapman_next = 0. night_prev = .true. night_next = .true. fperp_prev = 0. fperp_next = 0. cgrid%met(ipy)%nir_diffuse = 0.0 - !---------------------------------------------------------------! end if !------------------------------------------------------------------! @@ -2005,8 +2024,8 @@ subroutine update_met_drivers(cgrid) if (print_radinterp) then call dump_radinfo(cgrid,ipy,mprev,mnext,prevmet_timea & ,nextmet_timea,met_frq(iformat,iv),wprev,wnext & - ,secz_prev,secz_next,fperp_prev,fperp_next & - ,night_prev,night_next,'nddsf') + ,chapman_prev,chapman_next,fperp_prev & + ,fperp_next,night_prev,night_next,'nddsf') end if !------------------------------------------------------------------! end do @@ -2027,10 +2046,12 @@ subroutine update_met_drivers(cgrid) do ipy= 1,cgrid%npolygons !------------------------------------------------------------------! - ! Check whether this is day time or night time. We should ! - ! only do the full interpolation thing during the day, the night ! - ! should have no incoming radiation until we incorporate fireflies ! - ! (or a good twilight scheme) in the model. ! + ! Check whether this is daytime or night/twilight time. We ! + ! should carry out the sun-angle dependent interpolation only when ! + ! there is direct light. Although we check for daytime using the ! + ! actual cosine of the zenith angle, we interpolate irradiance ! + ! using the effective cosine of zenith angle, which prevents near- ! + ! singularities close to sunset times. ! !------------------------------------------------------------------! if (cgrid%cosz(ipy) > cosz_min) then @@ -2040,17 +2061,19 @@ subroutine update_met_drivers(cgrid) !---------------------------------------------------------------! select case (imetavg) case (0) - secz_prev = mean_daysecz(cgrid%lon(ipy),cgrid%lat(ipy) & - ,prevmet_timea,dt_radinterp,0.) - secz_next = mean_daysecz(cgrid%lon(ipy),cgrid%lat(ipy) & - ,nextmet_timea,dt_radinterp,0.) + chapman_prev = & + mean_chapman(cgrid%lon(ipy),cgrid%lat(ipy),prevmet_timea & + ,dt_radinterp,0.,.true.) + chapman_next = & + mean_chapman(cgrid%lon(ipy),cgrid%lat(ipy),nextmet_timea & + ,dt_radinterp,0.,.true.) case default - secz_prev = mean_daysecz(cgrid%lon(ipy),cgrid%lat(ipy) & - ,prevmet_timea,dt_radinterp & - ,met_frq(iformat,iv)) - secz_next = mean_daysecz(cgrid%lon(ipy),cgrid%lat(ipy) & - ,nextmet_timea,dt_radinterp & - ,met_frq(iformat,iv)) + chapman_prev = & + mean_chapman(cgrid%lon(ipy),cgrid%lat(ipy),prevmet_timea & + ,dt_radinterp,met_frq(iformat,iv),.true.) + chapman_next = & + mean_chapman(cgrid%lon(ipy),cgrid%lat(ipy),nextmet_timea & + ,dt_radinterp,met_frq(iformat,iv),.true.) end select !---------------------------------------------------------------! @@ -2060,8 +2083,8 @@ subroutine update_met_drivers(cgrid) ! Decide whether we want to use the previous or the next ! ! data based on the zenith angle. ! !---------------------------------------------------------------! - night_prev = secz_prev == 0. - night_next = secz_next == 0. + night_prev = chapman_prev == 0. + night_next = chapman_next == 0. !---------------------------------------------------------------! @@ -2085,8 +2108,8 @@ subroutine update_met_drivers(cgrid) ! current cosine of zenith angle. ! !------------------------------------------------------------! fperp_prev = 0. - fperp_next = cgrid%metinput(ipy)%vbdsf(mnext) * secz_next - cgrid%met(ipy)%par_beam = fperp_next * cgrid%cosz(ipy) + fperp_next = cgrid%metinput(ipy)%vbdsf(mnext) * chapman_next + cgrid%met(ipy)%par_beam = fperp_next * cgrid%eff_cosz(ipy) elseif (night_next) then !------------------------------------------------------------! ! Dusk time, the next time is zero so we only have ! @@ -2094,22 +2117,22 @@ subroutine update_met_drivers(cgrid) ! next time using the previous secant and scaling back with ! ! the current cosine of zenith angle. ! !------------------------------------------------------------! - fperp_prev = cgrid%metinput(ipy)%vbdsf(mprev) * secz_prev + fperp_prev = cgrid%metinput(ipy)%vbdsf(mprev) * chapman_prev fperp_next = 0. - cgrid%met(ipy)%par_beam = fperp_prev * cgrid%cosz(ipy) + cgrid%met(ipy)%par_beam = fperp_prev * cgrid%eff_cosz(ipy) else !----- Daytime, use both previous and next values. ----------! - fperp_next = cgrid%metinput(ipy)%vbdsf(mnext) * secz_next - fperp_prev = cgrid%metinput(ipy)%vbdsf(mprev) * secz_prev - cgrid%met(ipy)%par_beam = cgrid%cosz(ipy) & + fperp_next = cgrid%metinput(ipy)%vbdsf(mnext) * chapman_next + fperp_prev = cgrid%metinput(ipy)%vbdsf(mprev) * chapman_prev + cgrid%met(ipy)%par_beam = cgrid%eff_cosz(ipy) & * ( fperp_next * wnext & - + fperp_prev * wprev ) + + fperp_prev * wprev ) end if !---------------------------------------------------------------! else !----- Night time, assign it zero. -----------------------------! - secz_prev = 0. - secz_next = 0. + chapman_prev = 0. + chapman_next = 0. night_prev = .true. night_next = .true. fperp_prev = 0. @@ -2128,8 +2151,8 @@ subroutine update_met_drivers(cgrid) if (print_radinterp) then call dump_radinfo(cgrid,ipy,mprev,mnext,prevmet_timea & ,nextmet_timea,met_frq(iformat,iv),wprev,wnext & - ,secz_prev,secz_next,fperp_prev,fperp_next & - ,night_prev,night_next,'vbdsf') + ,chapman_prev,chapman_next,fperp_prev & + ,fperp_next,night_prev,night_next,'vbdsf') end if !------------------------------------------------------------------! end do @@ -2151,12 +2174,13 @@ subroutine update_met_drivers(cgrid) do ipy= 1,cgrid%npolygons !------------------------------------------------------------------! - ! Check whether this is day time or night time. We should ! - ! only do the full interpolation thing during the day, the night ! - ! should have no incoming radiation until we incorporate fireflies ! - ! (or a good twilight scheme) in the model. ! + ! Check whether this is day/twilight time or night time. We ! + ! should carry out the sun-angle dependent interpolation only when ! + ! there is diffuse light. Interpolation is done with the ! + ! effective cosine of zenith angle, which should be safe for ! + ! interpolations as long as cosz_min is not way too close to zero. ! !------------------------------------------------------------------! - if (cgrid%cosz(ipy) > cosz_min) then + if (cgrid%eff_cosz(ipy) > cosz_min) then !---------------------------------------------------------------! ! Define the normalisation factors for the previous and the ! @@ -2164,17 +2188,19 @@ subroutine update_met_drivers(cgrid) !---------------------------------------------------------------! select case (imetavg) case (0) - secz_prev = mean_daysecz(cgrid%lon(ipy),cgrid%lat(ipy) & - ,prevmet_timea,dt_radinterp,0.) - secz_next = mean_daysecz(cgrid%lon(ipy),cgrid%lat(ipy) & - ,nextmet_timea,dt_radinterp,0.) + chapman_prev = & + mean_chapman(cgrid%lon(ipy),cgrid%lat(ipy),prevmet_timea & + ,dt_radinterp,0.,.true.) + chapman_next = & + mean_chapman(cgrid%lon(ipy),cgrid%lat(ipy),nextmet_timea & + ,dt_radinterp,0.,.true.) case default - secz_prev = mean_daysecz(cgrid%lon(ipy),cgrid%lat(ipy) & - ,prevmet_timea,dt_radinterp & - ,met_frq(iformat,iv)) - secz_next = mean_daysecz(cgrid%lon(ipy),cgrid%lat(ipy) & - ,nextmet_timea,dt_radinterp & - ,met_frq(iformat,iv)) + chapman_prev = & + mean_chapman(cgrid%lon(ipy),cgrid%lat(ipy),prevmet_timea & + ,dt_radinterp,met_frq(iformat,iv),.true.) + chapman_next = & + mean_chapman(cgrid%lon(ipy),cgrid%lat(ipy),nextmet_timea & + ,dt_radinterp,met_frq(iformat,iv),.true.) end select !---------------------------------------------------------------! @@ -2184,8 +2210,8 @@ subroutine update_met_drivers(cgrid) ! Decide whether we want to use the previous or the next ! ! data based on the zenith angle. ! !---------------------------------------------------------------! - night_prev = secz_prev == 0. - night_next = secz_next == 0. + night_prev = chapman_prev == 0. + night_next = chapman_next == 0. !---------------------------------------------------------------! @@ -2209,8 +2235,8 @@ subroutine update_met_drivers(cgrid) ! current cosine of zenith angle. ! !------------------------------------------------------------! fperp_prev = 0. - fperp_next = cgrid%metinput(ipy)%vddsf(mnext) * secz_next - cgrid%met(ipy)%par_diffuse = fperp_next * cgrid%cosz(ipy) + fperp_next = cgrid%metinput(ipy)%vddsf(mnext) * chapman_next + cgrid%met(ipy)%par_diffuse = fperp_next * cgrid%eff_cosz(ipy) elseif (night_next) then !------------------------------------------------------------! ! Dusk time, the next time is zero so we only have ! @@ -2218,22 +2244,22 @@ subroutine update_met_drivers(cgrid) ! next time using the previous secant and scaling back with ! ! the current cosine of zenith angle. ! !------------------------------------------------------------! - fperp_prev = cgrid%metinput(ipy)%vddsf(mprev) * secz_prev + fperp_prev = cgrid%metinput(ipy)%vddsf(mprev) * chapman_prev fperp_next = 0. - cgrid%met(ipy)%par_diffuse = fperp_prev * cgrid%cosz(ipy) + cgrid%met(ipy)%par_diffuse = fperp_prev * cgrid%eff_cosz(ipy) else !----- Daytime, use both previous and next values. ----------! - fperp_next = cgrid%metinput(ipy)%vddsf(mnext) * secz_next - fperp_prev = cgrid%metinput(ipy)%vddsf(mprev) * secz_prev - cgrid%met(ipy)%par_diffuse = cgrid%cosz(ipy) & + fperp_next = cgrid%metinput(ipy)%vddsf(mnext) * chapman_next + fperp_prev = cgrid%metinput(ipy)%vddsf(mprev) * chapman_prev + cgrid%met(ipy)%par_diffuse = cgrid%eff_cosz(ipy) & * ( fperp_next * wnext & + fperp_prev * wprev ) end if !---------------------------------------------------------------! else !----- Night time, assign it zero. -----------------------------! - secz_prev = 0. - secz_next = 0. + chapman_prev = 0. + chapman_next = 0. night_prev = .true. night_next = .true. fperp_prev = 0. @@ -2252,8 +2278,8 @@ subroutine update_met_drivers(cgrid) if (print_radinterp) then call dump_radinfo(cgrid,ipy,mprev,mnext,prevmet_timea & ,nextmet_timea,met_frq(iformat,iv),wprev,wnext & - ,secz_prev,secz_next,fperp_prev,fperp_next & - ,night_prev,night_next,'vddsf') + ,chapman_prev,chapman_next,fperp_prev & + ,fperp_next,night_prev,night_next,'vddsf') end if !------------------------------------------------------------------! end do diff --git a/ED/src/dynamics/euler_driver.f90 b/ED/src/dynamics/euler_driver.f90 index 416104399..d965b977e 100644 --- a/ED/src/dynamics/euler_driver.f90 +++ b/ED/src/dynamics/euler_driver.f90 @@ -112,7 +112,7 @@ subroutine euler_timestep(cgrid) ,cmet%par_diffuse,cmet%nir_beam,cmet%nir_diffuse & ,cmet%geoht,cpoly%lsl(isi),cpoly%ntext_soil(:,isi) & ,cpoly%green_leaf_factor(:,isi),cgrid%lon(ipy) & - ,cgrid%lat(ipy),cgrid%cosz(ipy)) + ,cgrid%lat(ipy),cgrid%cosz(ipy),cgrid%eff_cosz(ipy)) !------------------------------------------------------------------------------! diff --git a/ED/src/dynamics/heun_driver.f90 b/ED/src/dynamics/heun_driver.f90 index c3c9631b6..1e8e29e86 100644 --- a/ED/src/dynamics/heun_driver.f90 +++ b/ED/src/dynamics/heun_driver.f90 @@ -112,7 +112,7 @@ subroutine heun_timestep(cgrid) ,cmet%par_diffuse,cmet%nir_beam,cmet%nir_diffuse & ,cmet%geoht,cpoly%lsl(isi),cpoly%ntext_soil(:,isi) & ,cpoly%green_leaf_factor(:,isi),cgrid%lon(ipy) & - ,cgrid%lat(ipy),cgrid%cosz(ipy)) + ,cgrid%lat(ipy),cgrid%cosz(ipy),cgrid%eff_cosz(ipy)) !------------------------------------------------------------------------------! diff --git a/ED/src/dynamics/hybrid_driver.f90 b/ED/src/dynamics/hybrid_driver.f90 index 7aaf1d59f..8e489d682 100644 --- a/ED/src/dynamics/hybrid_driver.f90 +++ b/ED/src/dynamics/hybrid_driver.f90 @@ -122,7 +122,7 @@ subroutine hybrid_timestep(cgrid) ,cpoly%lsl(isi),cpoly%ntext_soil(:,isi) & ,cpoly%green_leaf_factor(:,isi) & ,cgrid%lon(ipy),cgrid%lat(ipy) & - ,cgrid%cosz(ipy)) + ,cgrid%cosz(ipy),cgrid%eff_cosz(ipy)) diff --git a/ED/src/dynamics/radiate_driver.f90 b/ED/src/dynamics/radiate_driver.f90 index bea0ea3e1..bcf2868b6 100644 --- a/ED/src/dynamics/radiate_driver.f90 +++ b/ED/src/dynamics/radiate_driver.f90 @@ -16,7 +16,8 @@ subroutine canopy_radiation(cgrid) , patchtype ! ! structure use ed_para_coms , only : nthreads ! ! intent(in) use canopy_radiation_coms , only : cosz_min & ! intent(in) - , rshort_twilight_min ! ! intent(in) + , rshort_twilight_min & ! intent(in) + , find_eff_cosz ! ! intent(in) use consts_coms , only : pio180 ! ! intent(in) use grid_coms , only : nzg & ! intent(in) , nzs ! ! intent(in) @@ -59,7 +60,8 @@ subroutine canopy_radiation(cgrid) polyloop: do ipy = 1,cgrid%npolygons !----- Find the solar zenith angle [cosz] -------------------------------------! - cgrid%cosz(ipy) = ed_zen(cgrid%lon(ipy),cgrid%lat(ipy),current_time) + cgrid%cosz (ipy) = ed_zen(cgrid%lon(ipy),cgrid%lat(ipy),current_time) + cgrid%eff_cosz(ipy) = find_eff_cosz(cgrid%cosz(ipy)) !------------------------------------------------------------------------------! cpoly => cgrid%polygon(ipy) @@ -80,18 +82,21 @@ subroutine canopy_radiation(cgrid) * (mod(current_time%time + cgrid%lon(ipy) / 15. + 24., 24.) - 12.) sloperad = cpoly%slope(isi) * pio180 aspectrad = cpoly%aspect(isi) * pio180 - call angle_of_incid(cpoly%cosaoi(isi),cgrid%cosz(ipy),hrangl & + call angle_of_incid(cpoly%cosaoi(isi),cgrid%eff_cosz(ipy),hrangl & ,sloperad,aspectrad) !---------------------------------------------------------------------------! !---------------------------------------------------------------------------! - ! Find the two logicals, that will tell which part of the day we are at ! - ! least. ! + ! Find the two logicals, that will tell whether the time exceeds the ! + ! (lower) thresholds for twilight and daytime. ! !---------------------------------------------------------------------------! - daytime = cpoly%cosaoi(isi) > cosz_min .and. & + daytime = cgrid%cosz (ipy) > cosz_min .and. & + cpoly%cosaoi (isi) > cosz_min .and. & + cpoly%met(isi)%rshort > rshort_twilight_min + twilight = cgrid%eff_cosz(ipy) > cosz_min .and. & + cpoly%cosaoi (isi) > cosz_min .and. & cpoly%met(isi)%rshort > rshort_twilight_min - twilight = cpoly%met(isi)%rshort > rshort_twilight_min !---------------------------------------------------------------------------! @@ -110,8 +115,8 @@ subroutine canopy_radiation(cgrid) !---------------------------------------------------------------------------! maxcohort = 1 do ipa = 1,csite%npatches - cpatch=>csite%patch(ipa) - if ( cpatch%ncohorts>maxcohort ) maxcohort = cpatch%ncohorts + cpatch => csite%patch(ipa) + if ( cpatch%ncohorts > maxcohort ) maxcohort = cpatch%ncohorts end do !---------------------------------------------------------------------------! diff --git a/ED/src/dynamics/rk4_driver.F90 b/ED/src/dynamics/rk4_driver.F90 index ae758614c..7b1123bf9 100644 --- a/ED/src/dynamics/rk4_driver.F90 +++ b/ED/src/dynamics/rk4_driver.F90 @@ -128,7 +128,7 @@ subroutine rk4_timestep(cgrid) ,cmet%par_diffuse,cmet%nir_beam,cmet%nir_diffuse & ,cmet%geoht,cpoly%lsl(isi),cpoly%ntext_soil(:,isi) & ,cpoly%green_leaf_factor(:,isi),cgrid%lon(ipy) & - ,cgrid%lat(ipy),cgrid%cosz(ipy)) + ,cgrid%lat(ipy),cgrid%cosz(ipy),cgrid%eff_cosz(ipy)) !------------------------------------------------------------------------------! !------------------------------------------------------------------------------! diff --git a/ED/src/dynamics/rk4_integ_utils.f90 b/ED/src/dynamics/rk4_integ_utils.f90 index 46db5c36a..e830417f8 100644 --- a/ED/src/dynamics/rk4_integ_utils.f90 +++ b/ED/src/dynamics/rk4_integ_utils.f90 @@ -255,7 +255,7 @@ end subroutine odeint subroutine copy_met_2_rk4site(mzg,atm_ustar,atm_theiv,atm_vpdef,atm_theta,atm_tmp & ,atm_shv,atm_co2,zoff,exner,pcpg,qpcpg,dpcpg,prss,rshort & ,rlong,par_beam,par_diffuse,nir_beam,nir_diffuse,geoht & - ,lsl,ntext_soil,green_leaf_factor,lon,lat,cosz) + ,lsl,ntext_soil,green_leaf_factor,lon,lat,cosz,eff_cosz) use ed_max_dims , only : n_pft ! ! intent(in) use rk4_coms , only : rk4site ! ! structure use canopy_air_coms, only : ustmin8 ! ! intent(in) @@ -295,6 +295,7 @@ subroutine copy_met_2_rk4site(mzg,atm_ustar,atm_theiv,atm_vpdef,atm_theta,atm_tm real , intent(in) :: lon real , intent(in) :: lat real , intent(in) :: cosz + real , intent(in) :: eff_cosz !------------------------------------------------------------------------------------! @@ -326,6 +327,7 @@ subroutine copy_met_2_rk4site(mzg,atm_ustar,atm_theiv,atm_vpdef,atm_theta,atm_tm rk4site%lon = dble(lon ) rk4site%lat = dble(lat ) rk4site%cosz = dble(cosz ) + rk4site%eff_cosz = dble(eff_cosz ) rk4site%green_leaf_factor(:) = dble(green_leaf_factor(:)) !------------------------------------------------------------------------------------! diff --git a/ED/src/dynamics/rk4_misc.f90 b/ED/src/dynamics/rk4_misc.f90 index f847831ca..5c0b44ba6 100644 --- a/ED/src/dynamics/rk4_misc.f90 +++ b/ED/src/dynamics/rk4_misc.f90 @@ -3984,7 +3984,9 @@ subroutine print_rk4patch(y,csite,ipa) write (unit=*,fmt='(a,1x,es12.4)') ' Downward SW radiation : ',rk4site%rshort write (unit=*,fmt='(a,1x,es12.4)') ' Downward LW radiation : ',rk4site%rlong write (unit=*,fmt='(a,1x,es12.4)') ' Zenith angle (deg) : ' & - ,acos(rk4site%cosz) / pio1808 + ,acos(rk4site%cosz ) / pio1808 + write (unit=*,fmt='(a,1x,es12.4)') ' Effect. zenith angle (deg) : ' & + ,acos(rk4site%eff_cosz) / pio1808 write (unit=*,fmt='(80a)') ('=',k=1,80) write (unit=*,fmt='(a)' ) 'Leaf information (only those resolvable are shown): ' diff --git a/ED/src/init/ed_params.f90 b/ED/src/init/ed_params.f90 index 98b14d90b..d60dd8238 100644 --- a/ED/src/init/ed_params.f90 +++ b/ED/src/init/ed_params.f90 @@ -7311,8 +7311,7 @@ subroutine init_can_rad_params() , snow_albedo_nir & ! intent(out) , snow_emiss_tir & ! intent(out) , rshort_twilight_min & ! intent(out) - , cosz_min & ! intent(out) - , cosz_min8 ! ! intent(out) + , cosz_min ! ! intent(out) use pft_coms , only : is_grass & ! intent(in) , is_tropical & ! intent(in) , is_conifer ! ! intent(in) @@ -7327,18 +7326,17 @@ subroutine init_can_rad_params() !---------------------------------------------------------------------------------------! ! The following parameters are used to split the shortwave radiation into visible ! - ! and near-infrared radiation. ! + ! and near-infrared radiation. The NIR counterparts are defined in sub-routine ! + ! init_derived_params_after_xml() ! !---------------------------------------------------------------------------------------! - fvis_beam_def = 0.43 - fnir_beam_def = 1.0 - fvis_beam_def - fvis_diff_def = 0.52 - fnir_diff_def = 1.0 - fvis_diff_def + fvis_beam_def = 0.43 ! The remainder will be defined as fnir_beam_def + fvis_diff_def = 0.52 ! The remainder will be defined as fnir_diff_def !---------------------------------------------------------------------------------------! !---------------------------------------------------------------------------------------! - ! Clumping factor. This factor indicates the degree of clumpiness of leaves. !a + ! Clumping factor. This factor indicates the degree of clumpiness of leaves. ! ! 0 -- black hole ! ! 1 -- homogeneous, no clumping. ! !---------------------------------------------------------------------------------------! @@ -7491,12 +7489,25 @@ subroutine init_can_rad_params() !---------------------------------------------------------------------------------------! - ! These variables are the thresholds for things that should be computed during the ! - ! day time hours only. ! + ! Double-precision version of the minimum cosine of zenith angle for which we ! + ! assume day time conditions (i.e., when there is any chance for direct solar ! + ! irradiance) . If applied to effective cosine of zenith angle, this defines when ! + ! twilight conditions (i.e., there is any solar irradiance) are met or exceeded. ! + ! Also rshort_twilight_min provides an extra threshold for when twilight conditions are ! + ! met or exceeded. ! !---------------------------------------------------------------------------------------! + cosz_min = cos(89.9*pio180) ! cos(89.5*pio180) rshort_twilight_min = 0.5 - cosz_min = cos(89.*pio180) !cos(89.5*pio180) - cosz_min8 = dble(cosz_min) + !---------------------------------------------------------------------------------------! + + + + + !---------------------------------------------------------------------------------------! + ! This variable is the resolution (in degrees) for deriving the look-up table for ! + ! the modified Chapman's function. ! + !---------------------------------------------------------------------------------------! + dzen_ref = 5.0d-2 !---------------------------------------------------------------------------------------! @@ -8273,6 +8284,10 @@ subroutine init_derived_params_after_xml() use ed_therm_lib , only : calc_veg_hcap ! ! function use canopy_radiation_coms, only : ihrzrad & ! intent(in) , cci_hmax & ! intent(in) + , fvis_beam_def & ! intent(in) + , fvis_diff_def & ! intent(in) + , cosz_min & ! intent(in) + , dzen_ref & ! intent(in) , leaf_trans_vis & ! intent(in) , leaf_reflect_vis & ! intent(in) , wood_trans_vis & ! intent(in) @@ -8284,6 +8299,12 @@ subroutine init_derived_params_after_xml() , leaf_emiss_tir & ! intent(in) , wood_emiss_tir & ! intent(in) , orient_factor & ! intent(in) + , fnir_beam_def & ! intent(out) + , fnir_diff_def & ! intent(out) + , cosz_min8 & ! intent(out) + , nzen_ref & ! intent(out) + , zend_ref & ! intent(out) + , huestis_ref & ! intent(out) , leaf_scatter_vis & ! intent(out) , leaf_backscatter_vis & ! intent(out) , wood_scatter_vis & ! intent(out) @@ -8296,7 +8317,8 @@ subroutine init_derived_params_after_xml() , wood_backscatter_tir & ! intent(out) , phi1 & ! intent(out) , phi2 & ! intent(out) - , mu_bar ! ! intent(out) + , mu_bar & ! intent(out) + , set_huestis_lut ! ! sub-routine use rk4_coms , only : effarea_heat & ! intent(in) , effarea_evap & ! intent(in) , effarea_transp ! ! intent(in) @@ -8971,6 +8993,45 @@ subroutine init_derived_params_after_xml() end select !---------------------------------------------------------------------------------------! + + + + !---------------------------------------------------------------------------------------! + ! The following parameters are used to split the shortwave radiation into visible ! + ! and near-infrared radiation. The VIS counterparts have been defined either at sub- ! + ! routine init_can_rad_params or through the XML. ! + !---------------------------------------------------------------------------------------! + fnir_beam_def = 1.0 - fvis_beam_def + fnir_diff_def = 1.0 - fvis_diff_def + !---------------------------------------------------------------------------------------! + + + + + !---------------------------------------------------------------------------------------! + ! Double-precision version of the minimum cosine of zenith angle for which we ! + ! assume day time conditions (i.e., when there is any chance for direct solar ! + ! irradiance) . If applied to effective cosine of zenith angle, this defines when ! + ! twilight conditions (i.e., there is any solar irradiance) are met or exceeded. ! + !---------------------------------------------------------------------------------------! + cosz_min8 = dble(cosz_min ) + !---------------------------------------------------------------------------------------! + + + + + !---------------------------------------------------------------------------------------! + ! Assign the dimensions for the modified Chapman function look-up table, and ! + ! set the reference zenith angles and modified Chapman function values. ! + !---------------------------------------------------------------------------------------! + nzen_ref = ceiling( 1.80d2 / dzen_ref) + allocate(zend_ref (nzen_ref)) + allocate(huestis_ref(nzen_ref)) + call set_huestis_lut(nzen_ref,dzen_ref,zend_ref,huestis_ref) + !---------------------------------------------------------------------------------------! + + + !---------------------------------------------------------------------------------------! ! Scattering coefficients. Contrary to ED-2.1, these values are based on the ! ! description by by Sellers (1985) and the CLM technical manual, which includes the ! diff --git a/ED/src/io/ed_init_history.f90 b/ED/src/io/ed_init_history.f90 index f4d91cc49..c3f49f076 100644 --- a/ED/src/io/ed_init_history.f90 +++ b/ED/src/io/ed_init_history.f90 @@ -610,6 +610,8 @@ subroutine fill_history_grid_p11(cgrid,ipy,py_index) ,'TOTAL_BASAL_AREA_RECRUIT ' ,dsetrank,iparallel,.false.,foundvar) call hdf_getslab_r(cgrid%cosz (ipy:ipy) & ,'COSZ ' ,dsetrank,iparallel,.true. ,foundvar) + call hdf_getslab_r(cgrid%eff_cosz (ipy:ipy) & + ,'EFF_COSZ ' ,dsetrank,iparallel,.true. ,foundvar) call hdf_getslab_r(cgrid%cbudget_initialstorage (ipy:ipy) & ,'CBUDGET_INITIALSTORAGE_PY ',dsetrank,iparallel,.true. ,foundvar) call hdf_getslab_r(cgrid%cbudget_nep (ipy:ipy) & diff --git a/ED/src/io/ed_xml_config.f90 b/ED/src/io/ed_xml_config.f90 index 22e741ecd..478337cdb 100644 --- a/ED/src/io/ed_xml_config.f90 +++ b/ED/src/io/ed_xml_config.f90 @@ -194,7 +194,7 @@ recursive subroutine read_ed_xml_config(filename) end if - !******* MET PARAMS + !******* CANOPY AIR PARAMS call libxml2f90__ll_selectlist(TRIM(FILENAME)) call libxml2f90__ll_selecttag('ACT','config',1) !select upper level tag call libxml2f90__ll_exist('DOWN','can_air',ntag) !get number of met tags @@ -211,6 +211,37 @@ recursive subroutine read_ed_xml_config(filename) end if + !******* CANOPY RADIATION PARAMS + call libxml2f90__ll_selectlist(TRIM(FILENAME)) + call libxml2f90__ll_selecttag('ACT','config',1) !select upper level tag + call libxml2f90__ll_exist('DOWN','can_rad',ntag) !get number of met tags + print*,"CAN_RAD READ FROM FILE ::",ntag + if (ntag >= 1) then + do i=1,ntag + call libxml2f90__ll_selecttag('DOWN','can_rad',i) + + call getConfigREAL ('fvis_beam_def' ,'can_rad',i,rval,texist) + if (texist) fvis_beam_def = sngloff(rval,tiny_offset) + call getConfigREAL ('fvis_diff_def' ,'can_rad',i,rval,texist) + if (texist) fvis_diff_def = sngloff(rval,tiny_offset) + call getConfigREAL ('snow_albedo_vis' ,'can_rad',i,rval,texist) + if (texist) snow_albedo_vis = sngloff(rval,tiny_offset) + call getConfigREAL ('snow_albedo_nir' ,'can_rad',i,rval,texist) + if (texist) snow_albedo_nir = sngloff(rval,tiny_offset) + call getConfigREAL ('snow_emiss_tir' ,'can_rad',i,rval,texist) + if (texist) snow_emiss_tir = sngloff(rval,tiny_offset) + call getConfigREAL ('cosz_min' ,'can_rad',i,rval,texist) + if (texist) cosz_min = sngloff(rval,tiny_offset) + call getConfigREAL ('rshort_twilight_min','can_rad',i,rval,texist) + if (texist) rshort_twilight_min = sngloff(rval,tiny_offset) + call getConfigREAL ('chap_dzen' ,'can_rad',i,rval,texist) + if (texist) chap_dzen = rval + call libxml2f90__ll_selecttag('UP','config',1) !move back up to top level + + end do + end if + + !******* MISC call libxml2f90__ll_selectlist(TRIM(FILENAME)) @@ -1804,6 +1835,7 @@ subroutine write_ed_xml_config implicit none ! integer :: ival integer(4) :: i,ival + character(512) :: xfilout ! integer :: i ! character*(*) :: filename @@ -1864,6 +1896,21 @@ subroutine write_ed_xml_config + + !************ CANOPY RADIATION ***************** + call libxml2f90_ll_opentag("can_rad") + call putConfigREAL ("fvis_beam_def" , fvis_beam_def ) + call putConfigREAL ("fvis_diff_def" , fvis_diff_def ) + call putConfigREAL ("snow_albedo_vis" , snow_albedo_vis ) + call putConfigREAL ("snow_albedo_nir" , snow_albedo_nir ) + call putConfigREAL ("snow_emiss_tir" , snow_emiss_tir ) + call putConfigREAL ("cosz_min" , cosz_min ) + call putConfigREAL ("rshort_twilight_min" , rshort_twilight_min ) + call putConfigREAL8("chap_dzen" , chap_dzen ) + call libxml2f90_ll_closetag("can_read") + + + !************ MISC ***************** call libxml2f90_ll_opentag("misc") call putConfigINT("restart_mode",ied_init_mode) diff --git a/ED/src/memory/canopy_radiation_coms.f90 b/ED/src/memory/canopy_radiation_coms.f90 index fc4d65398..89dd384c9 100644 --- a/ED/src/memory/canopy_radiation_coms.f90 +++ b/ED/src/memory/canopy_radiation_coms.f90 @@ -307,6 +307,20 @@ module canopy_radiation_coms !---------------------------------------------------------------------------------------! + !---------------------------------------------------------------------------------------! + ! Define variables for computing the modified Chapman function. ! + !---------------------------------------------------------------------------------------! + !----- Zenith angle resolution of the look-up table for the Chapman function. ----------! + real(kind=8) :: dzen_ref + !----- Dimension of the look-up table bins. --------------------------------------------! + integer :: nzen_ref + !----- Reference zenith angles for the Chapman function. -------------------------------! + real(kind=8), dimension(:), allocatable :: zend_ref + !----- Reference values for the modified Chapman function (Huestis et al. 2001). -------! + real(kind=8), dimension(:), allocatable :: huestis_ref + !---------------------------------------------------------------------------------------! + + !=======================================================================================! !=======================================================================================! @@ -443,6 +457,193 @@ subroutine nullify_radscratch(cradscr) end subroutine nullify_radscratch !=======================================================================================! !=======================================================================================! + + + + + + + !=======================================================================================! + !=======================================================================================! + ! This function calculates the effect of sun angle increasing the optical depth in a ! + ! sphere. This allows accounting for the diffuse irradiance at twilight, using the ! + ! modified Chapman function following H01. The computation of the Chapman function is ! + ! computationally demanding, therefore, we pre-calculate the values for a broad range ! + ! of angles. We do not account for the terrain altitude as this effect is minimal and ! + ! would require substantial changes in the code. ! + ! ! + ! Reference: ! + ! ! + ! Huestis DL. 2001. Accurate evaluation of the chapman function for atmospheric ! + ! attenuation. J. Quant. Spectrosc. Radiat. Transf., 69: 709-721. ! + ! doi:10.1016/S0022-4073(00)00107-2 ! + !---------------------------------------------------------------------------------------! + subroutine set_huestis_lut(nzen,dzen,zend,huestis) + use consts_coms, only : erad & ! intent(in) + , ehgt & ! intent(in) + , pio1808 & ! intent(in) + , tiny_num8 & ! intent(in) + , lnexp_min8 & ! intent(in) + , lnexp_max8 ! ! intent(in) + + implicit none + !----- Arguments. -------------------------------------------------------------------! + integer , intent(in) :: nzen ! Number of bins + real(kind=8) , intent(in) :: dzen ! Bin width + real(kind=8), dimension(nzen), intent(out) :: zend ! Reference zenith angle + real(kind=8), dimension(nzen), intent(out) :: huestis ! Modified Chapman function + !----- Local variables. -------------------------------------------------------------! + real(kind=8) :: xcurve ! Curvature ratio + real(kind=8) :: lambda ! Integrating element for zenith angles + real(kind=8) :: dlambda ! Integrating width for zenith angles + real(kind=8) :: sin_zend ! Sine of zend + real(kind=8) :: sin_lambda ! Sine of lambda + real(kind=8) :: cos_lambda ! Cosine of lambda + real(kind=8) :: ln_kernel ! Natural logarithm of the kernel + real(kind=8) :: kernel ! Kernel + real(kind=8) :: integ_kernel ! Integral of the kernel from 0 to the current angle. + integer :: i ! Row counter + integer :: j ! Column counter + !------------------------------------------------------------------------------------! + + + !------------------------------------------------------------------------------------! + ! Dimensionless curvature ratio. ! + !------------------------------------------------------------------------------------! + xcurve = dble( erad / ehgt ) + !------------------------------------------------------------------------------------! + + + !------------------------------------------------------------------------------------! + ! Define the reference zenith angles. ! + !------------------------------------------------------------------------------------! + do i=1,nzen + zend (i) = min(1.80d2, dzen * dble(i-1) ) + end do + !------------------------------------------------------------------------------------! + + + !------------------------------------------------------------------------------------! + ! Initialise integrator and the modified Chapman function for the first ! + ! element, which should be always 1. ! + !------------------------------------------------------------------------------------! + integ_kern = 0.d0 + huestis(1) = 1.d0 + !------------------------------------------------------------------------------------! + + + !------------------------------------------------------------------------------------! + ! Integrate kernel for the first bin. We use a staggered approach for computing ! + ! the kernels and deriving the modified Chapman function evaluation. ! + !------------------------------------------------------------------------------------! + do i=2,nzen + !----- Find the mid points for integrand, and the trigonometric functions. -------! + lambda = 5.d-1 * ( zend(i-1) + zend(i) ) + dlambda = zend(i) - zend(i-1) + sin_zend = sin( zend(i) * pio1808 ) + sin_lambda = sin( lambda * pio1808 ) + cos_lambda = cos( lambda * pio1808 ) + !---------------------------------------------------------------------------------! + + + + !---------------------------------------------------------------------------------! + ! Find kernel. When the sine of lambda approaches zero, the kernel function ! + ! becomes undefined, so we use the limit values instead, as our goal is to ! + ! integrate the function. Also, we cap the natural logarithm of the kernel to ! + ! avoid floating point exceptions (though it should be safe with double ! + ! precision). ! + !---------------------------------------------------------------------------------! + if ( abs(sin_lambda) < tiny_num8 ) then + kernel = 0.d0 + elseif (abs(1.d0+cos_lambda) < tiny_num8) then + kernel = exp(lnexp_max8) + else + ln_kernel = xcurve * ( 1.d0 - sin_zend / sin_lambda ) + ln_kernel = max(lnexp_min8,min(lnexp_max8,ln_kernel)) + kernel = exp(ln_kernel) / ( 1.d0 + cos_lambda ) + end if + !---------------------------------------------------------------------------------! + + + + !---------------------------------------------------------------------------------! + ! Update the kernel integral and find the modified Chapman function. ! + !---------------------------------------------------------------------------------! + integ_kern = integ_kern + kernel * dlambda_ref * pio1808 + huestis(i) = min( exp(lnexp_max8), 1.d0 + xcurve*sin(zend*pio1808) * integ_kern) + !---------------------------------------------------------------------------------! + end do + !------------------------------------------------------------------------------------! + + return + end subroutine set_huestis_lut + !=======================================================================================! + !=======================================================================================! + + + + + + + !=======================================================================================! + !=======================================================================================! + ! This sub-routine interpolates the modified Chapman function (after H01) by inter- ! + ! polating the pre-calculated values from the lookup table, then finds the effective ! + ! cosine of the zenith angle. ! + !---------------------------------------------------------------------------------------! + real(kind=4) function find_eff_cosz(cosz) + use consts_coms, only : pio1808 & ! intent(in) + , tiny_offset ! ! intent(in) + implicit none + !----- Arguments. -------------------------------------------------------------------! + real(kind=4), intent(in ) :: cosz + !----- Local variables. -------------------------------------------------------------! + real(kind=8) :: zen + integer :: iprev + integer :: inext + real(kind=8) :: pwr_prev + real(kind=8) :: pwr_next + real(kind=8) :: chapman + !---- External functions. -----------------------------------------------------------! + real(kind=4), external :: sngloff + !------------------------------------------------------------------------------------! + + + + !----- Find the zenith angle and the nearest look-up table points. ------------------! + zen = acos( dble(cosz) ) / pio1808 + iprev = max(1 ,floor ( zen / dzen_ref )) + inext = min(nzen_ref,ceiling( zen / dzen_ref )) + !------------------------------------------------------------------------------------! + + + !------------------------------------------------------------------------------------! + ! Retrieve either the exact value from the look-up table (when the zenith ! + ! angle matches a value in the look-up table or when both values of the modified ! + ! Chapman function are identical) or log-linearly interpolate the values. ! + !------------------------------------------------------------------------------------! + if (huestis(iprev) == huestis(inext)) then + chapman = huestis_ref(iprev) + else + pwr_next = ( zen - zen_ref(iprev) ) / ( zen_ref(inext) - zen_ref(iprev)) + pwr_prev = 1.d0 - pwr_next + chapman = huestis_ref(iprev) ** pwr_prev * huestis_ref(inext) ** pwr_next + end if + !------------------------------------------------------------------------------------! + + + !------------------------------------------------------------------------------------! + ! The effective cosine of the zenith angle is the inverse of the Chapman ! + ! function. ! + !------------------------------------------------------------------------------------! + find_eff_cosz = sngloff( 1.d0 / chapman, tiny_num) + !------------------------------------------------------------------------------------! + + return + end function find_eff_cosz + !=======================================================================================! + !=======================================================================================! end module canopy_radiation_coms !==========================================================================================! !==========================================================================================! diff --git a/ED/src/memory/consts_coms.F90 b/ED/src/memory/consts_coms.F90 index 2345913b0..4ea3552b1 100644 --- a/ED/src/memory/consts_coms.F90 +++ b/ED/src/memory/consts_coms.F90 @@ -68,6 +68,7 @@ Module consts_coms , b_gcgs => gcgs & ! intent(in) , b_gg => gg & ! intent(in) , b_erad => erad & ! intent(in) + , b_ehgt => ehgt & ! intent(in) , b_spcon => spcon & ! intent(in) , b_spconkm => spconkm & ! intent(in) , b_eradi => eradi & ! intent(in) @@ -243,6 +244,7 @@ Module consts_coms real , parameter :: gcgs = b_gcgs real , parameter :: gg = b_gg real , parameter :: erad = b_erad + real , parameter :: ehgt = b_ehgt real , parameter :: spcon = b_spcon real , parameter :: spconkm = b_spconkm real , parameter :: eradi = b_eradi @@ -394,11 +396,11 @@ Module consts_coms !---------------------------------------------------------------------------------------! ! Universal constants ! !---------------------------------------------------------------------------------------! - real, parameter :: stefan = 5.6696e-8 ! Stefan-Boltzmann constant [ W/m�/K^4] - real, parameter :: boltzmann = 1.3806503e-23 ! Boltzmann constant [m�kg/s�/K] - real, parameter :: t00 = 273.15 ! 0�C [ �C] + real, parameter :: stefan = 5.6696e-8 ! Stefan-Boltzmann constant [ W/m2/K4] + real, parameter :: boltzmann = 1.3806503e-23 ! Boltzmann constant [m2kg/s2/K] + real, parameter :: t00 = 273.15 ! 0�C [ degC] real, parameter :: rmol = 8.314510 ! Molar gas constant [ J/mol/K] - real, parameter :: volmol = 0.022710980 ! Molar volume at STP [ m�] + real, parameter :: volmol = 0.022710980 ! Molar volume at STP [ m3] real, parameter :: volmoll = volmol*1e3 ! Molar volume at STP [ L] !---------------------------------------------------------------------------------------! @@ -439,10 +441,11 @@ Module consts_coms ! General Earth properties ! !---------------------------------------------------------------------------------------! real, parameter :: vonk = 0.40 ! Von K�rm�n constant [ ---] - real, parameter :: grav = 9.80665 ! Gravity acceleration [ m/s�] + real, parameter :: grav = 9.80665 ! Gravity acceleration [ m/s2] real, parameter :: erad = 6370997. ! Earth radius [ m] real, parameter :: erad2 = 2.*erad ! Earth diameter [ m] - real, parameter :: solar = 1.3533e3 ! Solar constant [ W/m�] + real, parameter :: ehgt = 8500. ! Earth's scale height (~ RT/g) [ m] + real, parameter :: solar = 1.3533e3 ! Solar constant [ W/m2] real, parameter :: p00 = 1.e5 ! Reference pressure [ Pa] real, parameter :: prefsea = 101325. ! Reference sea level pressure [ Pa] real, parameter :: p00i = 1. / p00 ! 1/p00 [ 1/Pa] @@ -466,9 +469,9 @@ Module consts_coms ! 10.11 (MU08). ! ! These terms could be easily made function of temperature in the future if needed be. ! !---------------------------------------------------------------------------------------! - real, parameter :: th_diff0 = 1.89e-5 ! Air thermal diffusivity [ m�/s] + real, parameter :: th_diff0 = 1.89e-5 ! Air thermal diffusivity [ m2/s] real, parameter :: dth_diff = 0.007 ! Temperature dependency slope [ 1/K] - real, parameter :: kin_visc0 = 1.33e-5 ! Kinematic viscosity [ m�/s] + real, parameter :: kin_visc0 = 1.33e-5 ! Kinematic viscosity [ m2/s] real, parameter :: dkin_visc = 0.007 ! Temperature dependency slope [ 1/K] !---------------------------------------------------------------------------------------! @@ -516,8 +519,8 @@ Module consts_coms !---------------------------------------------------------------------------------------! ! Liquid water properties ! !---------------------------------------------------------------------------------------! - real, parameter :: wdns = 1.000e3 ! Liquid water density [ kg/m�] - real, parameter :: wdnsi = 1./wdns ! Inverse of liquid water density [ m�/kg] + real, parameter :: wdns = 1.000e3 ! Liquid water density [ kg/m3] + real, parameter :: wdnsi = 1./wdns ! Inverse of liquid water density [ m3/kg] real, parameter :: cliq = 4.186e3 ! Liquid water specific heat (Cl) [ J/kg/K] real, parameter :: cliqi = 1./cliq ! Inverse of water heat capacity [ kg K/J] !---------------------------------------------------------------------------------------! @@ -527,12 +530,12 @@ Module consts_coms !---------------------------------------------------------------------------------------! ! Ice properties ! !---------------------------------------------------------------------------------------! - real, parameter :: idns = 9.167e2 ! "Hard" ice density [ kg/m�] - real, parameter :: idnsi = 1./idns ! Inverse of ice density [ m�/kg] - real, parameter :: fdns = 2.000e2 ! Frost density [ kg/m�] - real, parameter :: fdnsi = 1./fdns ! Inverse of frost density [ m�/kg] - real, parameter :: fsdns = 1.000e2 ! Fresh snow density [ kg/m�] - real, parameter :: fsdnsi = 1./fsdns ! Inverse of liquid water density [ m�/kg] + real, parameter :: idns = 9.167e2 ! "Hard" ice density [ kg/m3] + real, parameter :: idnsi = 1./idns ! Inverse of ice density [ m3/kg] + real, parameter :: fdns = 2.000e2 ! Frost density [ kg/m3] + real, parameter :: fdnsi = 1./fdns ! Inverse of frost density [ m3/kg] + real, parameter :: fsdns = 1.000e2 ! Fresh snow density [ kg/m3] + real, parameter :: fsdnsi = 1./fsdns ! Inverse of liquid water density [ m3/kg] real, parameter :: cice = 2.093e3 ! Ice specific heat (Ci) [ J/kg/K] real, parameter :: cicei = 1. / cice ! Inverse of ice heat capacity [ kg K/J] !---------------------------------------------------------------------------------------! diff --git a/ED/src/memory/ed_state_vars.F90 b/ED/src/memory/ed_state_vars.F90 index d72654bec..2d0d2f8a7 100644 --- a/ED/src/memory/ed_state_vars.F90 +++ b/ED/src/memory/ed_state_vars.F90 @@ -2740,6 +2740,10 @@ module ed_state_vars ! cosz_min) then - cloud = min(1.,max(0.,(c5 * cosz - rshort_tot) / (c4 * cosz))) - difrat = min(1.,max(0.,0.0604 / ( cosz -0.0223 ) + 0.0683)) + cloud = min(1.,max(0.,(c5 * eff_cosz - rshort_tot) / (c4 * eff_cosz))) + difrat = min(1.,max(0.,0.0604 / ( eff_cosz -0.0223 ) + 0.0683)) difrat = difrat + ( 1. - difrat ) * cloud vnrat = ( c1 - cloud*c2 ) / ( ( c1 - cloud*c3 ) + ( c1 - cloud*c2 )) @@ -214,8 +217,8 @@ end subroutine short2diff_sib ! Weiss, A., J. M. Norman, 1985: Partitioning solar radiation into direct and diffuse, ! ! visible and near-infrared components. Agric. For. Meteorol., 34, 205-213. (WN85) ! !---------------------------------------------------------------------------------------! - subroutine short_bdown_weissnorman(rshort_full,atm_prss,cosz,par_beam,par_diff,nir_beam & - ,nir_diff,rshort_diff) + subroutine short_bdown_weissnorman(rshort_full,atm_prss,cosz,eff_cosz,par_beam,par_diff & + ,nir_beam,nir_diff,rshort_diff) use consts_coms , only : solar & ! intent(in) , prefsea & ! intent(in) , twothirds ! ! intent(in) @@ -229,6 +232,7 @@ subroutine short_bdown_weissnorman(rshort_full,atm_prss,cosz,par_beam,par_diff,n real, intent(in) :: rshort_full ! Incident SW radiation (total) [ W/m2] real, intent(in) :: atm_prss ! Atmospheric pressure [ Pa] real, intent(in) :: cosz ! cos(zenith distance) [ ---] + real, intent(in) :: eff_cosz ! effective cos(zenith distance) [ ---] real, intent(out) :: par_beam ! Incident PAR (direct ) [ W/m2] real, intent(out) :: par_diff ! Incident PAR (diffuse) [ W/m2] real, intent(out) :: nir_beam ! Incident near-infrared (direct ) [ W/m2] @@ -252,8 +256,8 @@ subroutine short_bdown_weissnorman(rshort_full,atm_prss,cosz,par_beam,par_diff,n real :: ratio ! Ratio between obs. and expected [ ---] real :: aux_par ! Auxiliary variable [ ---] real :: aux_nir ! Auxiliary variable [ ---] - real :: secz ! sec(zenith distance) [ ---] - real :: log10secz ! log10[sec(zenith distance)] [ ---] + real :: chapman ! Chapman function [ ---] + real :: log10chapman ! log10[Chapman function] [ ---] real :: w10 ! Minimum Water absorption [ W/m2] !------------------------------------------------------------------------------------! ! Local constants. ! @@ -284,9 +288,12 @@ subroutine short_bdown_weissnorman(rshort_full,atm_prss,cosz,par_beam,par_diff,n par_diff = fvis_diff_def * rshort_diff nir_diff = fnir_diff_def * rshort_diff else - !----- Save 1/cos(zen), which is the secant. We will use this several times. ----! - secz = 1. / cosz - log10secz = log10(secz) + !---------------------------------------------------------------------------------! + ! Save the Chapman function, which is defined as 1/eff_cos(zen). We will use ! + ! this several times. ! + !---------------------------------------------------------------------------------! + chapman = 1. / eff_cosz + log10chapman = log10(chapman) !---------------------------------------------------------------------------------! @@ -300,8 +307,8 @@ subroutine short_bdown_weissnorman(rshort_full,atm_prss,cosz,par_beam,par_diff,n ! 3, and 9 of WN85. ! !---------------------------------------------------------------------------------! par_beam_pot = par_beam_top & - * exp ( par_beam_expext * (atm_prss / prefsea) * secz) * cosz - par_diff_pot = par2diff_sun * (par_beam_top - par_beam_pot) * cosz + * exp ( par_beam_expext * (atm_prss / prefsea) * chapman) * eff_cosz + par_diff_pot = par2diff_sun * (par_beam_top - par_beam_pot) * eff_cosz par_full_pot = par_beam_pot + par_diff_pot !---------------------------------------------------------------------------------! @@ -310,7 +317,8 @@ subroutine short_bdown_weissnorman(rshort_full,atm_prss,cosz,par_beam,par_diff,n !---------------------------------------------------------------------------------! ! Find the NIR absorption of 10 mm of precipitable water, using WN85 eqn. 6. ! !---------------------------------------------------------------------------------! - w10 = solar * 10 ** (wn85_06(1) + log10secz * (wn85_06(2) + wn85_06(3)*log10secz)) + w10 = solar & + * 10. ** (wn85_06(1) + log10chapman * (wn85_06(2) + wn85_06(3)*log10chapman)) !---------------------------------------------------------------------------------! @@ -320,8 +328,9 @@ subroutine short_bdown_weissnorman(rshort_full,atm_prss,cosz,par_beam,par_diff,n ! 4, 5, and 10 of WN85. ! !---------------------------------------------------------------------------------! nir_beam_pot = ( nir_beam_top & - * exp ( nir_beam_expext * (atm_prss/prefsea) * secz) - w10 ) * cosz - nir_diff_pot = nir2diff_sun * ( nir_beam_top - nir_beam_pot - w10 ) * cosz + * exp ( nir_beam_expext * (atm_prss/prefsea) * chapman) - w10 ) + * eff_cosz + nir_diff_pot = nir2diff_sun * ( nir_beam_top - nir_beam_pot - w10 ) * eff_cosz nir_full_pot = nir_beam_pot + nir_diff_pot !---------------------------------------------------------------------------------! @@ -405,6 +414,7 @@ subroutine short_bdown_clearidx(rshort_full,cosz,par_beam,par_diff,nir_beam,nir_ !----- Arguments. -------------------------------------------------------------------! real, intent(in) :: rshort_full ! Incident SW radiation (total) [ W/m2] real, intent(in) :: cosz ! cos(zenith distance) [ ---] + real, intent(in) :: eff_cosz ! effective cos(zenith distance) [ ---] real, intent(out) :: par_beam ! Incident PAR (direct ) [ W/m2] real, intent(out) :: par_diff ! Incident PAR (diffuse) [ W/m2] real, intent(out) :: nir_beam ! Incident near-infrared (direct ) [ W/m2] @@ -445,7 +455,7 @@ subroutine short_bdown_clearidx(rshort_full,cosz,par_beam,par_diff,nir_beam,nir_ nir_diff = rshort_diff - par_diff else !----- Total radiation at the top of the atmosphere [ W/m2], using ED defaults. -! - rshort_toa = solar * cosz + rshort_toa = solar * eff_cosz !---------------------------------------------------------------------------------! @@ -490,7 +500,7 @@ subroutine short_bdown_clearidx(rshort_full,cosz,par_beam,par_diff,nir_beam,nir_ ! include this term. ! !---------------------------------------------------------------------------------! if (apply_bx10_corr) then - fdiff = fdiff_1st / ( (1.0-fdiff_1st) * cosz + fdiff_1st) + fdiff = fdiff_1st / ( (1.0-fdiff_1st) * eff_cosz + fdiff_1st) else fdiff = fdiff_1st end if @@ -568,17 +578,17 @@ end subroutine update_rad_avg ! output. ! !---------------------------------------------------------------------------------------! subroutine dump_radinfo(cgrid,ipy,mprev,mnext,prevmet_timea,nextmet_timea,met_frq & - ,wprev,wnext,secz_prev,secz_next,fperp_prev,fperp_next & + ,wprev,wnext,chapman_prev,chapman_next,fperp_prev,fperp_next & ,night_prev,night_next,variable) - use ed_state_vars , only : edtype ! ! structure - use ed_max_dims , only : str_len ! ! intent(in) - use ed_misc_coms , only : simtime & ! intent(in) - , current_time ! ! intent(in) - use met_driver_coms , only : nbdsf_file & ! intent(in) - , nddsf_file & ! intent(in) - , vbdsf_file & ! intent(in) - , vddsf_file ! ! intent(in) - use canopy_radiation_coms, only : cosz_min ! ! intent(in) + use ed_state_vars , only : edtype ! ! structure + use ed_max_dims , only : str_len ! ! intent(in) + use ed_misc_coms , only : simtime & ! intent(in) + , current_time ! ! intent(in) + use met_driver_coms , only : nbdsf_file & ! intent(in) + , nddsf_file & ! intent(in) + , vbdsf_file & ! intent(in) + , vddsf_file ! ! intent(in) + use canopy_radiation_coms, only : cosz_min ! ! intent(in) implicit none !----- Arguments. -------------------------------------------------------------------! type(edtype) , target :: cgrid @@ -590,8 +600,8 @@ subroutine dump_radinfo(cgrid,ipy,mprev,mnext,prevmet_timea,nextmet_timea,met_fr real , intent(in) :: met_frq real , intent(in) :: wprev real , intent(in) :: wnext - real , intent(in) :: secz_prev - real , intent(in) :: secz_next + real , intent(in) :: chapman_prev + real , intent(in) :: chapman_next real , intent(in) :: fperp_prev real , intent(in) :: fperp_next logical , intent(in) :: night_prev @@ -609,8 +619,8 @@ subroutine dump_radinfo(cgrid,ipy,mprev,mnext,prevmet_timea,nextmet_timea,met_fr logical :: night_now logical :: itsthere !----- Local constants. -------------------------------------------------------------! - character(len=18) , parameter :: fmth = '(3(a,3x),15(a,1x))' - character(len=35) , parameter :: fmtd = '(3(a,3x),3(11x,l1,1x),12(f12.3,1x))' + character(len=18) , parameter :: fmth = '(3(a,3x),16(a,1x))' + character(len=35) , parameter :: fmtd = '(3(a,3x),3(11x,l1,1x),13(f12.3,1x))' character(len=29) , parameter :: fmtt = '(2(i2.2,a),i4.4,1x,3(i2.2,a))' !----- Locally saved variables, to delete and re-create scratch files. --------------! logical , save :: first_nbdsf = .true. @@ -652,11 +662,11 @@ subroutine dump_radinfo(cgrid,ipy,mprev,mnext,prevmet_timea,nextmet_timea,met_fr write (unit=87,fmt=fmth) & ' WHEN_PREV',' WHEN_NOW' & ,' WHEN_NEXT',' NIGHT_PREV',' NIGHT_NOW' & - ,' NIGHT_NEXT',' MET_FRQ',' WPREV' & - ,' WNEXT',' FLUX_PREV',' SECZ_PREV' & - ,' FPERP_PREV',' FLUX_NEXT',' SECZ_NEXT' & - ,' FPERP_NEXT',' FLUX_NOW',' COSZ_NOW' & - ,' FPERP_NOW' + ,' NIGHT_NEXT',' MET_FRQ',' WPREV' & + ,' WNEXT',' FLUX_PREV','CHAPMAN_PREV' & + ,' FPERP_PREV',' FLUX_NEXT','CHAPMAN_NEXT' & + ,' FPERP_NEXT',' FLUX_NOW',' COSZ_NOW' & + ,'EFF_COSZ_NOW',' FPERP_NOW' close (unit=87,status='keep') !------------------------------------------------------------------------------! @@ -691,11 +701,11 @@ subroutine dump_radinfo(cgrid,ipy,mprev,mnext,prevmet_timea,nextmet_timea,met_fr write (unit=87,fmt=fmth) & ' WHEN_PREV',' WHEN_NOW' & ,' WHEN_NEXT',' NIGHT_PREV',' NIGHT_NOW' & - ,' NIGHT_NEXT',' MET_FRQ',' WPREV' & - ,' WNEXT',' FLUX_PREV',' SECZ_PREV' & - ,' FPERP_PREV',' FLUX_NEXT',' SECZ_NEXT' & - ,' FPERP_NEXT',' FLUX_NOW',' COSZ_NOW' & - ,' FPERP_NOW' + ,' NIGHT_NEXT',' MET_FRQ',' WPREV' & + ,' WNEXT',' FLUX_PREV','CHAPMAN_PREV' & + ,' FPERP_PREV',' FLUX_NEXT','CHAPMAN_NEXT' & + ,' FPERP_NEXT',' FLUX_NOW',' COSZ_NOW' & + ,'EFF_COSZ_NOW',' FPERP_NOW' close (unit=87,status='keep') !------------------------------------------------------------------------------! @@ -729,11 +739,11 @@ subroutine dump_radinfo(cgrid,ipy,mprev,mnext,prevmet_timea,nextmet_timea,met_fr write (unit=87,fmt=fmth) & ' WHEN_PREV',' WHEN_NOW' & ,' WHEN_NEXT',' NIGHT_PREV',' NIGHT_NOW' & - ,' NIGHT_NEXT',' MET_FRQ',' WPREV' & - ,' WNEXT',' FLUX_PREV',' SECZ_PREV' & - ,' FPERP_PREV',' FLUX_NEXT',' SECZ_NEXT' & - ,' FPERP_NEXT',' FLUX_NOW',' COSZ_NOW' & - ,' FPERP_NOW' + ,' NIGHT_NEXT',' MET_FRQ',' WPREV' & + ,' WNEXT',' FLUX_PREV','CHAPMAN_PREV' & + ,' FPERP_PREV',' FLUX_NEXT','CHAPMAN_NEXT' & + ,' FPERP_NEXT',' FLUX_NOW',' COSZ_NOW' & + ,'EFF_COSZ_NOW',' FPERP_NOW' close (unit=87,status='keep') !------------------------------------------------------------------------------! @@ -768,11 +778,11 @@ subroutine dump_radinfo(cgrid,ipy,mprev,mnext,prevmet_timea,nextmet_timea,met_fr write (unit=87,fmt=fmth) & ' WHEN_PREV',' WHEN_NOW' & ,' WHEN_NEXT',' NIGHT_PREV',' NIGHT_NOW' & - ,' NIGHT_NEXT',' MET_FRQ',' WPREV' & - ,' WNEXT',' FLUX_PREV',' SECZ_PREV' & - ,' FPERP_PREV',' FLUX_NEXT',' SECZ_NEXT' & - ,' FPERP_NEXT',' FLUX_NOW',' COSZ_NOW' & - ,' FPERP_NOW' + ,' NIGHT_NEXT',' MET_FRQ',' WPREV' & + ,' WNEXT',' FLUX_PREV','CHAPMAN_PREV' & + ,' FPERP_PREV',' FLUX_NEXT','CHAPMAN_NEXT' & + ,' FPERP_NEXT',' FLUX_NOW',' COSZ_NOW' & + ,'EFF_COSZ_NOW',' FPERP_NOW' close (unit=87,status='keep') !------------------------------------------------------------------------------! @@ -797,7 +807,7 @@ subroutine dump_radinfo(cgrid,ipy,mprev,mnext,prevmet_timea,nextmet_timea,met_fr ,nextmet_timea%year,nextmet_timea%hour,'.',nextmet_timea%min & ,'.',nextmet_timea%sec,'UTC' fperp_now = fperp_prev * wprev + fperp_next * wnext - night_now = cgrid%cosz(ipy) <= cosz_min + night_now = cgrid%eff_cosz(ipy) <= cosz_min !------------------------------------------------------------------------------------! @@ -806,8 +816,9 @@ subroutine dump_radinfo(cgrid,ipy,mprev,mnext,prevmet_timea,nextmet_timea,met_fr !------------------------------------------------------------------------------------! open(unit=87,file=trim(thisfile),status='old',action='write',position='append') write(unit=87,fmt=fmtd) when_prev,when_now,when_next,night_prev,night_now,night_next & - ,met_frq,wprev,wnext,flux_prev,secz_prev,fperp_prev,flux_next & - ,secz_next,fperp_next,flux_now,cgrid%cosz(ipy),fperp_now + ,met_frq,wprev,wnext,flux_prev,chapman_prev,fperp_prev & + ,flux_next,chapman_next,fperp_next,flux_now,cgrid%cosz(ipy) & + ,cgrid%eff_cosz(ipy),fperp_now close(unit=87,status='keep') !------------------------------------------------------------------------------------! return @@ -824,37 +835,38 @@ end subroutine dump_radinfo !=======================================================================================! ! This subroutine calculates angle of incidence based on local slope and aspect. ! !---------------------------------------------------------------------------------------! - subroutine angle_of_incid(aoi,cosz,solar_hour_aspect,slope,terrain_aspect) + subroutine angle_of_incid(aoi,eff_cosz,solar_hour_aspect,slope,terrain_aspect) + use consts_coms, only : tiny_offset ! ! intent(in) implicit none !----- Arguments. -------------------------------------------------------------------! - real, intent(in) :: cosz ! Cosine of zenithal angle + real, intent(in) :: eff_cosz ! Effective cosine of zenithal angle real, intent(in) :: slope ! Terrain slope real, intent(in) :: solar_hour_aspect ! horizontal location of the sun defined with ! the same reference as terrain aspect. real, intent(in) :: terrain_aspect ! Terrain aspect real, intent(out) :: aoi ! Angle of incidence !----- Local variables. -------------------------------------------------------------! - real(kind=8) :: cosz8 ! Double prec. counterpart of cosz - real(kind=8) :: sinz8 ! Sine of zenithal angle - real(kind=8) :: slope8 ! Double prec. counterpart of slope - real(kind=8) :: sh_asp8 ! Double prec. counterpart of solar_hour_aspect - real(kind=8) :: terr_asp8 ! Double prec. counterpart of terrain_aspect - real(kind=8) :: aoi8 ! Double prec. counterpart of aoi - !----- Local parameters. ------------------------------------------------------------! - real(kind=8), parameter :: tiny_offset=1.d-20 + real(kind=8) :: eff_cosz8 ! Double prec. counterpart of eff_cosz + real(kind=8) :: eff_sinz8 ! Effective sine of zenithal angle + real(kind=8) :: slope8 ! Double prec. counterpart of slope + real(kind=8) :: sh_asp8 ! Double prec. counterpart of solar_hour_aspect + real(kind=8) :: terr_asp8 ! Double prec. counterpart of terrain_aspect + real(kind=8) :: aoi8 ! Double prec. counterpart of aoi !----- External functions. ----------------------------------------------------------! real , external :: sngloff !------------------------------------------------------------------------------------! - cosz8 = dble(cosz) - sinz8 = sqrt(1.d0-cosz8*cosz8) - slope8 = dble(slope) - sh_asp8 = dble(solar_hour_aspect) - terr_asp8 = dble(terrain_aspect) - if (cosz8 < 0.d0) then + eff_cosz8 = dble(eff_cosz) + eff_sinz8 = sqrt(1.d0-eff_cosz8*eff_cosz8) + slope8 = dble(slope) + sh_asp8 = dble(solar_hour_aspect) + terr_asp8 = dble(terrain_aspect) + if (eff_cosz8 <= 0.d0) then aoi8 = 0.d0 else - aoi8 = max(0.d0, cosz8*dcos(slope8) + sinz8*dsin(slope8)*dcos(sh_asp8-terr_asp8)) + aoi8 = eff_cosz8 * dcos(slope8) & + + eff_sinz8 * dsin(slope8) * dcos(sh_asp8-terr_asp8) + aoi8 = max(0.d0, aoi8) end if aoi = sngloff(aoi8,tiny_offset) @@ -961,31 +973,35 @@ end function ed_zen !=======================================================================================! !=======================================================================================! - ! This function computes the average secant of the daytime zenith angle. In case ! - ! the period of integration is that accounts for the zenith angle is 0. or less than ! - ! one time step, then the average is the actual value. Night-time periods are ignored ! - ! and if there is no daytime value, then we set it to 0. ! + ! This function computes the average Chapman function (secant of the daytime ! + ! effective zenith angle). When the period of integration is less than one time step, ! + ! then the average is the actual value. Night-time periods are ignored and if there is ! + ! no daytime value, then we set it to 0. Twilight may or may not be accounted for, ! + ! depending on the add_twilight flag. ! !---------------------------------------------------------------------------------------! - real function mean_daysecz(plon,plat,whena,dt,tmax) + real function mean_chapman(plon,plat,whena,dt,tmax,add_twilight) use update_derived_utils , only : update_model_time_dm ! ! sub-routine use ed_misc_coms , only : simtime ! ! structure - use canopy_radiation_coms, only : cosz_min ! ! intent(in) + use canopy_radiation_coms, only : cosz_min & ! intent(in) + , find_eff_cosz ! ! sub-routine implicit none !------ Arguments. ------------------------------------------------------------------! - real(kind=4) , intent(in) :: plon - real(kind=4) , intent(in) :: plat - type(simtime), intent(in) :: whena - real(kind=4) , intent(in) :: dt - real(kind=4) , intent(in) :: tmax + real(kind=4) , intent(in) :: plon ! Longitude + real(kind=4) , intent(in) :: plat ! Latitude + type(simtime), intent(in) :: whena ! Current time + real(kind=4) , intent(in) :: dt ! Time step + real(kind=4) , intent(in) :: tmax ! Maximum time + logical , intent(in) :: add_twilight ! Add twilight times to the average? !------ Local variables. ------------------------------------------------------------! - type(simtime) :: now ! Current time - integer :: is ! Step counter - integer :: nsteps ! Number of steps to perform the average - real :: dtfit ! Delta-t that nicely fits within tmax - real :: dtnow ! Delta-t for this time - real :: cosz ! Declination - real :: daytot ! Total time that was daytime - real :: mean_daycosz ! Average cosine of zenith angle + type(simtime) :: now ! Current time + integer :: is ! Step counter + integer :: nsteps ! # of steps to perform the average + real :: dtfit ! Delta-t that nicely fits within tmax + real :: dtnow ! Delta-t for this time + real :: cosz ! Cosine of zenith angle + real :: eff_cosz ! Effective cosine of zenith angle + real :: daytot ! Total time that was daytime + real :: mean_eff_cosz ! Average eff. cosine of zenith angle !------------------------------------------------------------------------------------! @@ -995,12 +1011,13 @@ real function mean_daysecz(plon,plat,whena,dt,tmax) !------------------------------------------------------------------------------------! if (dt >= tmax) then !----- Less than one time, no average necessary. ---------------------------------! - cosz = ed_zen(plon,plat,whena) - if (cosz > cosz_min) then - mean_daysecz = 1.0 / cosz + cosz = ed_zen(plon,plat,whena) + eff_cosz = find_eff_cosz(cosz) + if ( ( cosz > cosz_min) .or. (add_twilight .and. eff_cosz > cosz_min) ) then + mean_chapman = 1.0 / eff_cosz else !----- Night-time, set the mean to zero. --------------------------------------! - mean_daysecz = 0.0 + mean_chapman = 0.0 end if else @@ -1012,8 +1029,8 @@ real function mean_daysecz(plon,plat,whena,dt,tmax) dtfit = tmax / real(nsteps) !---------------------------------------------------------------------------------! - mean_daycosz = 0.0 - daytot = 0.0 + mean_eff_cosz = 0.0 + daytot = 0.0 do is=1,nsteps !----- Get the current time. --------------------------------------------------! now = whena @@ -1021,12 +1038,13 @@ real function mean_daysecz(plon,plat,whena,dt,tmax) call update_model_time_dm(now,dtnow) !----- Get the cosine of the zenith angle. ------------------------------------! - cosz = ed_zen(plon,plat,now) + cosz = ed_zen(plon,plat,now) + eff_cosz = find_eff_cosz(cosz) !----- Add to the integral only if it this value is valid. --------------------! - if (cosz > cosz_min) then - mean_daycosz = mean_daycosz + dtfit * cosz - daytot = daytot + dtfit + if ( ( cosz > cosz_min) .or. (add_twilight .and. eff_cosz > cosz_min) ) then + mean_eff_cosz = mean_eff_cosz + dtfit * eff_cosz + daytot = daytot + dtfit end if !------------------------------------------------------------------------------! end do @@ -1037,17 +1055,17 @@ real function mean_daysecz(plon,plat,whena,dt,tmax) !---------------------------------------------------------------------------------! ! Find the normalisation factor. ! !---------------------------------------------------------------------------------! - if (daytot > 0.0 .and. mean_daycosz > 0.0) then - mean_daycosz = mean_daycosz / daytot - mean_daysecz = 1.0 / mean_daycosz + if (daytot > 0.0 .and. mean_eff_cosz > 0.0) then + mean_eff_cosz = mean_eff_cosz / daytot + mean_chapman = 1.0 / mean_eff_cosz else - mean_daysecz = 0.0 + mean_chapman = 0.0 end if !---------------------------------------------------------------------------------! end if return - end function mean_daysecz + end function mean_chapman !=======================================================================================! !=======================================================================================! end module radiate_utils