Skip to content

Commit

Permalink
First implementation of the curvature effects on the "effective zenit…
Browse files Browse the repository at this point in the history 
…h angle". This is

based on the modified Chapman function proposed by Huestis et al. (2001). This function
allows diffuse irradiance to occur even when the sun is (slightly) below horizon, and
thus allowing twilight to be correctly accounted for.

This method should also reduce crashes during sunrise/sunset, although this does not
change that it is a very strongly encouraged to set IMETAVG correctly (a future pull
request will migrate IMETAVG from ED2IN to the met driver header.

This still needs to be tested, more commits with bug fixes likely coming soon.
  • Loading branch information
@mpaiao
mpaiao committed Jan 6, 2024
1 parent 4809be5 commit 83f6c70
Show file tree
Hide file tree
Showing 16 changed files with 723 additions and 342 deletions.
406 changes: 216 additions & 190 deletions ED/src/driver/ed_met_driver.f90
View file

Large diffs are not rendered by default.

2 changes: 1 addition & 1 deletion ED/src/dynamics/euler_driver.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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))
!------------------------------------------------------------------------------!


Expand Down
2 changes: 1 addition & 1 deletion ED/src/dynamics/heun_driver.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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))
!------------------------------------------------------------------------------!


Expand Down
2 changes: 1 addition & 1 deletion ED/src/dynamics/hybrid_driver.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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))



Expand Down
23 changes: 14 additions & 9 deletions ED/src/dynamics/radiate_driver.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand All @@ -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
!---------------------------------------------------------------------------!


Expand All @@ -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
!---------------------------------------------------------------------------!

Expand Down
2 changes: 1 addition & 1 deletion ED/src/dynamics/rk4_driver.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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))
!------------------------------------------------------------------------------!

!------------------------------------------------------------------------------!
Expand Down
4 changes: 3 additions & 1 deletion ED/src/dynamics/rk4_integ_utils.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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
!------------------------------------------------------------------------------------!


Expand Down Expand Up @@ -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(:))
!------------------------------------------------------------------------------------!

Expand Down
4 changes: 3 additions & 1 deletion ED/src/dynamics/rk4_misc.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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): '
Expand Down
87 changes: 74 additions & 13 deletions ED/src/init/ed_params.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand All @@ -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. !
!---------------------------------------------------------------------------------------!
Expand Down Expand Up @@ -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
!---------------------------------------------------------------------------------------!


Expand Down Expand Up @@ -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)
Expand All @@ -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)
Expand All @@ -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)
Expand Down Expand Up @@ -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 !
Expand Down
2 changes: 2 additions & 0 deletions ED/src/io/ed_init_history.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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) &
Expand Down
Loading

0 comments on commit 83f6c70

Please sign in to comment.