From 2a0ff4dc554a7518275e74b55484bac8b4f2fe3a Mon Sep 17 00:00:00 2001 From: Marcos Longo Date: Tue, 9 Jan 2024 11:06:35 -0800 Subject: [PATCH] Fixed multiple bugs in the implementation, and added a logical variable (use_huestis_eff_cosz) that allows switching between the use of the "effective" cosine of zenith angle and the actual value (ED-2.2 default). --- ED/build/make/include.mk.macosx | 28 ++-- ED/src/driver/ed_met_driver.f90 | 2 +- ED/src/dynamics/farq_katul.f90 | 12 -- ED/src/dynamics/radiate_driver.f90 | 17 ++- ED/src/init/ed_params.f90 | 26 ++-- ED/src/io/ed_xml_config.f90 | 42 +++--- ED/src/memory/canopy_radiation_coms.f90 | 169 ++++++++++++++---------- ED/src/utils/radiate_utils.f90 | 8 +- 8 files changed, 168 insertions(+), 136 deletions(-) diff --git a/ED/build/make/include.mk.macosx b/ED/build/make/include.mk.macosx index 69ed3365c..48bd49726 100644 --- a/ED/build/make/include.mk.macosx +++ b/ED/build/make/include.mk.macosx @@ -5,7 +5,7 @@ #----- Define make (gnu make works best). -------------------------------------------------# -MAKE=/usr/bin/make +MAKE=/usr/local/bin/make #------------------------------------------------------------------------------------------# #----- Libraries. -------------------------------------------------------------------------# @@ -20,8 +20,8 @@ BASE=$(ED_ROOT)/build/ # libraries compiled with the same compiler you set for F_COMP and C_COMP. You may still # # be able to compile without HDF5 but it will not run. # #------------------------------------------------------------------------------------------# -ZLIB_PATH=/Users/mlongo/Util/ED2_Libs/zlib-1.2.11 -HDF5_PATH=/Users/mlongo/Util/ED2_Libs/hdf5-1.10.2 +ZLIB_PATH=/usr/local +HDF5_PATH=/usr/local HDF5_INCS=-I$(HDF5_PATH)/include HDF5C_INCS=$(HDF5_INCS) HDF5_LIBS=-L$(ZLIB_PATH)/lib -lz -L$(HDF5_PATH)/lib -lhdf5_fortran -lhdf5 -lhdf5_hl @@ -51,10 +51,9 @@ USE_INTERF=1 #################################### COMPILER SETTINGS ##################################### CMACH=MAC_OS_X FC_TYPE=GNU -F_COMP=gfortran -C_COMP=gcc -LOADER=gfortran -C_LOADER=gcc +F_COMP=/usr/local/bin/gfortran-13 +C_COMP=/usr/local/bin/gcc-13 +LOADER=/usr/local/bin/gfortran-13 LIBS= MOD_EXT=mod ############################################################################################ @@ -88,20 +87,17 @@ endif ifeq ($(KIND_COMP),A) F_OPTS= -O0 -ffree-line-length-none -g -fimplicit-none -Wall -finit-real=snan \ -finit-integer=-2147483648 -ffpe-trap=invalid,zero,overflow,underflow \ - -fcheck=all -frecursive -fsignaling-nans -Werror -mmacosx-version-min=10.12 \ - -fopenmp -static - C_OPTS= -O0 -DLITTLE -g -static + -fcheck=all -frecursive -fsignaling-nans -Werror -fopenmp -fbacktrace -static + C_OPTS= -O0 -DLITTLE -g -fbacktrace -static LOADER_OPTS= -O0 -ffree-line-length-none -g -fimplicit-none -Wall -finit-real=snan \ -finit-integer=-2147483648 -ffpe-trap=invalid,zero,overflow,underflow \ - -fcheck=all -frecursive -fsignaling-nans -Werror \ - -mmacosx-version-min=10.12 -fopenmp + -fcheck=all -frecursive -fsignaling-nans -Werror -fopenmp -fbacktrace #---------------------------------------------------------------------------------------# endif ifeq ($(KIND_COMP),E) - F_OPTS= -O3 -ffree-line-length-none -frecursive -mmacosx-version-min=10.12 -fopenmp \ - -static - C_OPTS= -O0 -DLITTLE -g -static - LOADER_OPTS= -O3 -ffree-line-length-none -frecursive -mmacosx-version-min=10.12 -fopenmp + F_OPTS= -O3 -ffree-line-length-none -frecursive -fopenmp -fbacktrace -static + C_OPTS= -O3 -DLITTLE -g -fbacktrace -static + LOADER_OPTS= -O3 -ffree-line-length-none -frecursive -fopenmp -fbacktrace #---------------------------------------------------------------------------------------# endif #------------------------------------------------------------------------------------------# diff --git a/ED/src/driver/ed_met_driver.f90 b/ED/src/driver/ed_met_driver.f90 index 32d419586..01ec03b7a 100644 --- a/ED/src/driver/ed_met_driver.f90 +++ b/ED/src/driver/ed_met_driver.f90 @@ -1387,7 +1387,7 @@ subroutine update_met_drivers(cgrid) case default chapman_prev = & mean_chapman(cgrid%lon(ipy),cgrid%lat(ipy),prevmet_timea & - ,dt_radinterp,met_frq(iformat,iv)) + ,dt_radinterp,met_frq(iformat,iv),.true.) end select !---------------------------------------------------------------! diff --git a/ED/src/dynamics/farq_katul.f90 b/ED/src/dynamics/farq_katul.f90 index 31ac9b186..529d34a04 100644 --- a/ED/src/dynamics/farq_katul.f90 +++ b/ED/src/dynamics/farq_katul.f90 @@ -409,9 +409,6 @@ subroutine optimization_solver8(ib,ipft,lambda, real(kind=8) :: test_dcidg real(kind=8) :: test_dfcdg real(kind=8) :: test_dfedg - real(kind=8) :: test_fc_light - real(kind=8) :: test_fc_rubp - real(kind=8) :: test_fc_3rd real(kind=8) :: opt_ci_light real(kind=8) :: opt_ci_rubp real(kind=8) :: opt_ci_3rd @@ -852,15 +849,6 @@ subroutine photosynthesis_stomata_solver8(ib,gsc,limit_case, real(kind=8) :: a,b,c !! Coefficients of the quadratic equation to solve ci real(kind=8) :: rad !! sqrt(b2-4ac) real(kind=8) :: dbdg,dcdg !! derivatives of b,c wrt. gsc - real(kind=8) :: ci_rubp !! ci for rubp-limited scenario - real(kind=8) :: dcidg_rubp !! derivative of ci wrt. gsc for rubp-limited scenario - real(kind=8) :: dfcdg_rubp !! derivative of fc wrt. gsc for rubp-limited scenario - real(kind=8) :: ci_light !! ci for light-limited scenario - real(kind=8) :: dcidg_light !! derivative of ci wrt. gsc for light-limited scenario - real(kind=8) :: dfcdg_light !! derivative of fc wrt. gsc for light-limited scenario - real(kind=8) :: ci_3rd !! ci for TPU/CO2-limited scenario - real(kind=8) :: dcidg_3rd !! derivative of ci wrt. gsc for TPU/CO2-limited scenario - real(kind=8) :: dfcdg_3rd !! derivative of fc wrt. gsc for TPU/CO2-limited scenario !------------------------------------------------------------------------------------! diff --git a/ED/src/dynamics/radiate_driver.f90 b/ED/src/dynamics/radiate_driver.f90 index bcf2868b6..debf12cbc 100644 --- a/ED/src/dynamics/radiate_driver.f90 +++ b/ED/src/dynamics/radiate_driver.f90 @@ -16,6 +16,7 @@ subroutine canopy_radiation(cgrid) , patchtype ! ! structure use ed_para_coms , only : nthreads ! ! intent(in) use canopy_radiation_coms , only : cosz_min & ! intent(in) + , use_huestis_eff_cosz & ! intent(in) , rshort_twilight_min & ! intent(in) , find_eff_cosz ! ! intent(in) use consts_coms , only : pio180 ! ! intent(in) @@ -91,12 +92,16 @@ subroutine canopy_radiation(cgrid) ! Find the two logicals, that will tell whether the time exceeds the ! ! (lower) thresholds for twilight and daytime. ! !---------------------------------------------------------------------------! - 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 + if (use_huestis_eff_cosz) then + daytime = cgrid%cosz (ipy) > cosz_min .and. & + cpoly%met(isi)%rshort > rshort_twilight_min + twilight = cpoly%cosaoi (isi) > cosz_min .and. & + cpoly%met(isi)%rshort > rshort_twilight_min + else + daytime = cpoly%cosaoi (isi) > cosz_min .and. & + cpoly%met(isi)%rshort > rshort_twilight_min + twilight = cpoly%met(isi)%rshort > rshort_twilight_min + end if !---------------------------------------------------------------------------! diff --git a/ED/src/init/ed_params.f90 b/ED/src/init/ed_params.f90 index d60dd8238..b26e66e91 100644 --- a/ED/src/init/ed_params.f90 +++ b/ED/src/init/ed_params.f90 @@ -7305,13 +7305,13 @@ subroutine init_can_rad_params() , wood_emiss_tir & ! intent(out) , fvis_beam_def & ! intent(out) , fvis_diff_def & ! intent(out) - , fnir_beam_def & ! intent(out) - , fnir_diff_def & ! intent(out) , snow_albedo_vis & ! intent(out) , snow_albedo_nir & ! intent(out) , snow_emiss_tir & ! intent(out) , rshort_twilight_min & ! intent(out) - , cosz_min ! ! intent(out) + , cosz_min & ! intent(out) + , use_huestis_eff_cosz & ! intent(out) + , dzen_ref ! ! intent(out) use pft_coms , only : is_grass & ! intent(in) , is_tropical & ! intent(in) , is_conifer ! ! intent(in) @@ -7496,18 +7496,26 @@ subroutine init_can_rad_params() ! 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) + cosz_min = cos(89.95*pio180) ! cos(89.5*pio180) rshort_twilight_min = 0.5 !---------------------------------------------------------------------------------------! + !---------------------------------------------------------------------------------------! + ! Flag to determine whether or not to use the modified Chapman function to find the ! + ! "effective" cosine of zenith angle that accounts for the Earth's curvature (and so ! + ! it allows solving twilight irradiance). To fall back to ED-2.2 default, set this ! + ! flag to .false. ! + !---------------------------------------------------------------------------------------! + use_huestis_eff_cosz = .true. + !---------------------------------------------------------------------------------------! !---------------------------------------------------------------------------------------! ! This variable is the resolution (in degrees) for deriving the look-up table for ! ! the modified Chapman's function. ! !---------------------------------------------------------------------------------------! - dzen_ref = 5.0d-2 + dzen_ref = 1.0d-2 !---------------------------------------------------------------------------------------! @@ -8303,7 +8311,7 @@ subroutine init_derived_params_after_xml() , fnir_diff_def & ! intent(out) , cosz_min8 & ! intent(out) , nzen_ref & ! intent(out) - , zend_ref & ! intent(out) + , zen_ref & ! intent(out) , huestis_ref & ! intent(out) , leaf_scatter_vis & ! intent(out) , leaf_backscatter_vis & ! intent(out) @@ -9024,10 +9032,10 @@ subroutine init_derived_params_after_xml() ! 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)) + nzen_ref = 1 + ceiling( 1.80d2 / dzen_ref) + allocate(zen_ref (nzen_ref)) allocate(huestis_ref(nzen_ref)) - call set_huestis_lut(nzen_ref,dzen_ref,zend_ref,huestis_ref) + call set_huestis_lut() !---------------------------------------------------------------------------------------! diff --git a/ED/src/io/ed_xml_config.f90 b/ED/src/io/ed_xml_config.f90 index 478337cdb..63d7028f9 100644 --- a/ED/src/io/ed_xml_config.f90 +++ b/ED/src/io/ed_xml_config.f90 @@ -220,22 +220,24 @@ recursive subroutine read_ed_xml_config(filename) 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 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 getConfigINT ('use_huestis_eff_cosz','can_rad',i,ival,texist) + if (texist) use_huestis_eff_cosz = ival == 1 + call getConfigREAL ('dzen_ref' ,'can_rad',i,rval,texist) + if (texist) dzen_ref = rval call libxml2f90__ll_selecttag('UP','config',1) !move back up to top level end do @@ -1906,7 +1908,13 @@ subroutine write_ed_xml_config 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 ) + if (use_huestis_eff_cosz) then + ival = 1 + else + ival = 0 + end if + call putConfigINT ("use_huestis_eff_cosz", ival ) + call putConfigREAL8("dzen_ref" , dzen_ref ) call libxml2f90_ll_closetag("can_read") diff --git a/ED/src/memory/canopy_radiation_coms.f90 b/ED/src/memory/canopy_radiation_coms.f90 index 89dd384c9..d3a1f5242 100644 --- a/ED/src/memory/canopy_radiation_coms.f90 +++ b/ED/src/memory/canopy_radiation_coms.f90 @@ -310,12 +310,14 @@ module canopy_radiation_coms !---------------------------------------------------------------------------------------! ! Define variables for computing the modified Chapman function. ! !---------------------------------------------------------------------------------------! + !----- Use the modified Chapman (.true. means yes, .false. sets eff_cosz = cosz). ------! + logical :: use_huestis_eff_cosz !----- 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 + real(kind=8), dimension(:), allocatable :: zen_ref !----- Reference values for the modified Chapman function (Huestis et al. 2001). -------! real(kind=8), dimension(:), allocatable :: huestis_ref !---------------------------------------------------------------------------------------! @@ -478,7 +480,7 @@ end subroutine nullify_radscratch ! 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) + subroutine set_huestis_lut() use consts_coms, only : erad & ! intent(in) , ehgt & ! intent(in) , pio1808 & ! intent(in) @@ -487,22 +489,17 @@ subroutine set_huestis_lut(nzen,dzen,zend,huestis) , 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_zen ! Sine of zen 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 + real(kind=8) :: integ_kern ! Integral of the kernel from 0 to the current angle. + integer :: i ! Row counter integer :: j ! Column counter !------------------------------------------------------------------------------------! @@ -517,63 +514,80 @@ subroutine set_huestis_lut(nzen,dzen,zend,huestis) !------------------------------------------------------------------------------------! ! Define the reference zenith angles. ! !------------------------------------------------------------------------------------! - do i=1,nzen - zend (i) = min(1.80d2, dzen * dble(i-1) ) + do i=1,nzen_ref + zen_ref(i) = min(1.80d2, dzen_ref * dble(i-1) ) end do !------------------------------------------------------------------------------------! !------------------------------------------------------------------------------------! - ! Initialise integrator and the modified Chapman function for the first ! - ! element, which should be always 1. ! + ! Initialise the modified Chapman function for the first element, which should ! + ! be always 1. ! !------------------------------------------------------------------------------------! - integ_kern = 0.d0 - huestis(1) = 1.d0 + integ_kern = 0.d0 + huestis_ref(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. ! + ! Integrate kernel. 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 ) - !---------------------------------------------------------------------------------! - - + i_loop: do i=2,nzen_ref !---------------------------------------------------------------------------------! - ! 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). ! + ! Find the zenith of the reference angle, and initialise integrator. ! !---------------------------------------------------------------------------------! - 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 + integ_kern = 0.d0 + sin_zen = sin( zen_ref(i) * pio1808 ) !---------------------------------------------------------------------------------! + !---------------------------------------------------------------------------------! + ! Integrate kernel. ! + !---------------------------------------------------------------------------------! + j_loop: do j=2,i + !----- Find the mid points for integrand, and the trigonometric functions. ----! + lambda = 5.d-1 * ( zen_ref(j-1) + zen_ref(j) ) + dlambda = zen_ref(j) - zen_ref(j-1) + 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_zen / 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. + !------------------------------------------------------------------------------! + integ_kern = integ_kern + kernel * dlambda * pio1808 + !------------------------------------------------------------------------------! + end do j_loop !---------------------------------------------------------------------------------! - ! Update the kernel integral and find the modified Chapman function. ! + ! Find the modified Chapman function at zen = zen_ref(i). ! !---------------------------------------------------------------------------------! - integ_kern = integ_kern + kernel * dlambda_ref * pio1808 - huestis(i) = min( exp(lnexp_max8), 1.d0 + xcurve*sin(zend*pio1808) * integ_kern) + huestis_ref(i) = min( exp(lnexp_max8), 1.d0 + xcurve * sin_zen * integ_kern) !---------------------------------------------------------------------------------! - end do + end do i_loop !------------------------------------------------------------------------------------! return @@ -593,8 +607,8 @@ end subroutine set_huestis_lut ! cosine of the zenith angle. ! !---------------------------------------------------------------------------------------! real(kind=4) function find_eff_cosz(cosz) - use consts_coms, only : pio1808 & ! intent(in) - , tiny_offset ! ! intent(in) + use consts_coms, only : pio1808 & ! intent(in) + , tiny_num ! ! intent(in) implicit none !----- Arguments. -------------------------------------------------------------------! real(kind=4), intent(in ) :: cosz @@ -610,34 +624,47 @@ real(kind=4) function find_eff_cosz(cosz) !------------------------------------------------------------------------------------! - - !----- 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 )) !------------------------------------------------------------------------------------! + ! Check whether to actually find the effective cosine of the zenith angle or just ! + ! use the actual one. ! + !------------------------------------------------------------------------------------! + if (use_huestis_eff_cosz) then + !----- 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 - !------------------------------------------------------------------------------------! + !---------------------------------------------------------------------------------! + ! 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_ref(iprev) == huestis_ref(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) + + !---------------------------------------------------------------------------------! + ! The effective cosine of the zenith angle is the inverse of the Chapman ! + ! function. ! + !---------------------------------------------------------------------------------! + find_eff_cosz = sngloff( 1.d0 / chapman, tiny_num) + !---------------------------------------------------------------------------------! + else + !---------------------------------------------------------------------------------! + ! Do not use the effective cosine of the zenith angle. ! + !---------------------------------------------------------------------------------! + find_eff_cosz = cosz + !---------------------------------------------------------------------------------! + end if !------------------------------------------------------------------------------------! return diff --git a/ED/src/utils/radiate_utils.f90 b/ED/src/utils/radiate_utils.f90 index b14dea861..2f488b347 100644 --- a/ED/src/utils/radiate_utils.f90 +++ b/ED/src/utils/radiate_utils.f90 @@ -328,7 +328,7 @@ subroutine short_bdown_weissnorman(rshort_full,atm_prss,cosz,eff_cosz,par_beam,p ! 4, 5, and 10 of WN85. ! !---------------------------------------------------------------------------------! nir_beam_pot = ( nir_beam_top & - * exp ( nir_beam_expext * (atm_prss/prefsea) * chapman) - w10 ) + * 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 @@ -404,8 +404,8 @@ end subroutine short_bdown_weissnorman ! of two competing plant species. Int. J. Biometeorol. 54, 283-295. ! ! doi:10.1007/s00484-009-0279-3 (BX10). ! !---------------------------------------------------------------------------------------! - subroutine short_bdown_clearidx(rshort_full,cosz,par_beam,par_diff,nir_beam,nir_diff & - ,rshort_diff) + subroutine short_bdown_clearidx(rshort_full,cosz,eff_cosz,par_beam,par_diff,nir_beam & + ,nir_diff,rshort_diff) use consts_coms , only : solar & ! intent(in) , lnexp_min & ! intent(in) , lnexp_max ! ! intent(in) @@ -836,7 +836,7 @@ end subroutine dump_radinfo ! This subroutine calculates angle of incidence based on local slope and aspect. ! !---------------------------------------------------------------------------------------! subroutine angle_of_incid(aoi,eff_cosz,solar_hour_aspect,slope,terrain_aspect) - use consts_coms, only : tiny_offset ! ! intent(in) + use rk4_coms, only : tiny_offset ! ! intent(in) implicit none !----- Arguments. -------------------------------------------------------------------! real, intent(in) :: eff_cosz ! Effective cosine of zenithal angle