From 1b3179fbe50fd4e02224908892273955d3f24e26 Mon Sep 17 00:00:00 2001 From: Remko NIJZINK Date: Thu, 19 May 2022 18:52:40 +0200 Subject: [PATCH 1/6] corrected fpar calculations --- VOM_Fortran/VOM-code/transpmodel.f90 | 49 +++++++++++++--------------- 1 file changed, 23 insertions(+), 26 deletions(-) diff --git a/VOM_Fortran/VOM-code/transpmodel.f90 b/VOM_Fortran/VOM-code/transpmodel.f90 index 1f6caf9..1af9d84 100644 --- a/VOM_Fortran/VOM-code/transpmodel.f90 +++ b/VOM_Fortran/VOM-code/transpmodel.f90 @@ -1651,25 +1651,25 @@ subroutine vom_gstom () select case(i_lai_function) case(1) ! no dynamic LAI, fpar is set to 1 - jactt(:,ii) = (1.d0 - p_E ** (-(i_alpha * fpar_lt(ii) * par_h(th_)) & + jactt(:,ii) = (1.d0 - p_E ** (-(i_alpha * par_h(th_)) & & / jmaxt_h(:))) * jmaxt_h(:) ! (3.23), (Out[311]) case(2) ! dynamic LAI, with fpar-calculation - jactt(:,ii) = (1.d0 - p_E ** (-(i_alpha * fpar_lt(ii) * par_h(th_)) & + jactt(:,ii) = (1.d0 - p_E ** (-(i_alpha * par_h(th_)) & & / jmaxt_h(:))) * jmaxt_h(:) ! (3.23), (Out[311]) case(3) ! shaded and sunlit, with diffuse and direct radiation - jactt(:,ii) = ( (1.d0 - p_E ** (-(i_alpha * fpar_lt(ii) * (pardir_h(th_) + pardiff_h(th_)) ) & + jactt(:,ii) = ( (1.d0 - p_E ** (-(i_alpha * (pardir_h(th_) + pardiff_h(th_)) ) & & / jmaxt_h(:))) * jmaxt_h(:) * frac_sunt(ii) ) + & - & ( (1.d0 - p_E ** (-(i_alpha * fpar_lt(ii) * pardiff_h(th_) ) & + & ( (1.d0 - p_E ** (-(i_alpha * pardiff_h(th_) ) & & / jmaxt_h(:))) * jmaxt_h(:) * frac_shadet(ii) ) case(4) ! shaded and sunlit, diffuse and direct radiation, different jact-values shaded-sunlit - jactt(:,ii) = (1.d0 - p_E ** (-(i_alpha * fpar_lt(ii) * (pardir_h(th_) + pardiff_h(th_)) ) & - & / jmaxt_h(:))) * jmaxt_h(:) + jactt(:,ii) = (1.d0 - p_E ** (-(i_alpha * (pardir_h(th_) + pardiff_h(th_)) ) & + & / jmaxt_h(:))) * jmaxt_h(:) * frac_sunt(ii) - jactts(:,ii) = (1.d0 - p_E ** (-(i_alpha * fpar_lt(ii) * pardiff_h(th_) ) & - & / jmaxts_h(:))) * jmaxts_h(:) + jactts(:,ii) = (1.d0 - p_E ** (-(i_alpha * pardiff_h(th_) ) & + & / jmaxts_h(:))) * jmaxts_h(:) * frac_shadet(ii) end select end do @@ -1689,17 +1689,17 @@ subroutine vom_gstom () & / jmaxg_h(:))) * jmaxg_h(:) ! (3.23), (Out[311]) case(3) ! shaded and sunlit, with diffuse and direct radiation - jactg(:,ii) = ( (1.d0 - p_E ** (-(i_alpha * fpar_lg(ii) * (pardir_h(th_) + pardiff_h(th_)) ) & + jactg(:,ii) = ( (1.d0 - p_E ** (-(i_alpha * (pardir_h(th_) + pardiff_h(th_)) ) & & / jmaxg_h(:))) * jmaxg_h(:) * frac_sung(ii) ) + & - & ( (1.d0 - p_E ** (-(i_alpha * fpar_lg(ii) * pardiff_h(th_) ) & + & ( (1.d0 - p_E ** (-(i_alpha * pardiff_h(th_) ) & & / jmaxg_h(:))) * jmaxg_h(:) * frac_shadeg(ii) ) case(4) ! shaded and sunlit, diffuse and direct radiation, different jact-values shaded-sunlit - jactg(:,ii) = (1.d0 - p_E ** (-(i_alpha * fpar_lg(ii) * (pardir_h(th_) + pardiff_h(th_)) ) & - & / jmaxg_h(:))) * jmaxg_h(:) + jactg(:,ii) = (1.d0 - p_E ** (-(i_alpha * (pardir_h(th_) + pardiff_h(th_)) ) & + & / jmaxg_h(:))) * jmaxg_h(:) * frac_sung(ii) - jactgs(:,ii) = (1.d0 - p_E ** (-(i_alpha * fpar_lg(ii) * pardiff_h(th_) ) & - & / jmaxgs_h(:))) * jmaxgs_h(:) + jactgs(:,ii) = (1.d0 - p_E ** (-(i_alpha * pardiff_h(th_) ) & + & / jmaxgs_h(:))) * jmaxgs_h(:) * frac_shadeg(ii) end select @@ -1739,11 +1739,8 @@ subroutine vom_gstom () gstomt = 0.d0 endif transpt = p_a * vd_h(th_) * gstomt ! (3.28) transpiration rate in mol/s - if(i_lai_function == 4) then - etmt__ = frac_sunt(2) * (transpt * 18.d0) / (10.d0 ** 6.d0) ! transpiration rate in m/s - else - etmt__ = (transpt * 18.d0) / (10.d0 ** 6.d0) ! transpiration rate in m/s - end if + etmt__ = (transpt * 18.d0) / (10.d0 ** 6.d0) ! transpiration rate in m/s + ! * calculate stomatal conductance grasses @@ -1774,7 +1771,7 @@ subroutine vom_gstom () transpg(:, :, :) = p_a * vd_h(th_) * gstomg(:, :, :) ! (3.28) transpiration rate in mol/s if(i_lai_function == 4) then do ii = 1,3 - etmg__(:, :, ii) = frac_sung(ii) * (transpg(:, :, ii) * 18.d0) / (10.d0 ** 6.d0) ! transpiration rate in m/s + etmg__(:, :, ii) = (transpg(:, :, ii) * 18.d0) / (10.d0 ** 6.d0) ! transpiration rate in m/s end do else etmg__(:,:,:) = (transpg(:, :, :) * 18.d0) / (10.d0 ** 6.d0) ! transpiration rate in m/s @@ -1813,7 +1810,7 @@ subroutine vom_gstom () endif transpts = p_a * vd_h(th_) * gstomts ! (3.28) transpiration rate in mol/s - etmts__ = frac_shadet(2) * (transpts * 18.d0) / (10.d0 ** 6.d0) ! transpiration rate in m/s + etmts__ = (transpts * 18.d0) / (10.d0 ** 6.d0) ! transpiration rate in m/s else gstomts = 0.d0 transpts = 0.d0 @@ -1851,7 +1848,7 @@ subroutine vom_gstom () transpgs(:,:,:) = p_a * vd_h(th_) * gstomgs(:,:,:) ! (3.28) transpiration rate in mol/s do ii = 1,3 - etmgs__(:,:,ii) = frac_shadeg(ii) * (transpgs(:,:,ii) * 18.d0) / (10.d0 ** 6.d0) ! transpiration rate in m/s + etmgs__(:,:,ii) = (transpgs(:,:,ii) * 18.d0) / (10.d0 ** 6.d0) ! transpiration rate in m/s end do else gstomgs(:,:,:) = 0.d0 @@ -2246,7 +2243,7 @@ subroutine vom_add_hourly () do ii = 1,3 !loop for LAI values if(i_lai_function == 4) then - asst__(:,ii) = frac_sunt(ii) * ( (4.d0 * ca_h(th_) * gstomt + 8.d0 * gammastar & + asst__(:,ii) = ( (4.d0 * ca_h(th_) * gstomt + 8.d0 * gammastar & & * gstomt + jactt(:,ii) - 4.d0 * rlt_h(:) - SQRT((-4.d0 & & * ca_h(th_) * gstomt + 8.d0 * gammastar * gstomt & & + jactt(:,ii) - 4.d0 * rlt_h(:)) ** 2.d0 + 16.d0 & @@ -2254,7 +2251,7 @@ subroutine vom_add_hourly () & + jactt(:,ii) + 8.d0 * rlt_h(:)))) / 8.d0 ) ! (3.22) ; (Out[319]) - assts__(:,ii) = frac_shadet(ii) * ( (4.d0 * ca_h(th_) * gstomts + 8.d0 * gammastar & + assts__(:,ii) = ( (4.d0 * ca_h(th_) * gstomts + 8.d0 * gammastar & & * gstomts + jactts(:,ii) - 4.d0 * rlts_h(:) - SQRT((-4.d0 & & * ca_h(th_) * gstomts + 8.d0 * gammastar * gstomts & & + jactts(:,ii) - 4.d0 * rlts_h(:)) ** 2.d0 + 16.d0 & @@ -2279,7 +2276,7 @@ subroutine vom_add_hourly () do ii = 1,3 !loop for LAI values do jj = 1,3 !loop for CAI_g values if(i_lai_function == 4) then - assg__(jj,:,ii) = frac_sung(ii) * ( (4.d0 * ca_h(th_) * gstomg(jj,:,ii) + 8.d0 * gammastar & + assg__(jj,:,ii) = ( (4.d0 * ca_h(th_) * gstomg(jj,:,ii) + 8.d0 * gammastar & & * gstomg(jj,:,ii) + jactg(:,ii) - 4.d0 * rlg_h(:) & & - SQRT((-4.d0 * ca_h(th_) * gstomg(jj,:,ii) + 8.d0 & & * gammastar * gstomg(jj,:,ii) + jactg(:,ii) - 4.d0 & @@ -2287,7 +2284,7 @@ subroutine vom_add_hourly () & * gstomg(jj,:,ii) * (8.d0 * ca_h(th_) * gstomg(jj,:,ii) & & + jactg(:,ii) + 8.d0 * rlg_h(:)))) / 8.d0) ! (3.22); (Out[319]) - assgs__(jj,:,ii) = frac_shadeg(ii) * ( (4.d0 * ca_h(th_) * gstomgs(jj,:,ii) + 8.d0 * gammastar & + assgs__(jj,:,ii) = ( (4.d0 * ca_h(th_) * gstomgs(jj,:,ii) + 8.d0 * gammastar & & * gstomgs(jj,:,ii) + jactgs(:,ii) - 4.d0 * rlgs_h(:) & & - SQRT((-4.d0 * ca_h(th_) * gstomgs(jj,:,ii) + 8.d0 & & * gammastar * gstomgs(jj,:,ii) + jactgs(:,ii) - 4.d0 & From 1332f51bceb977147c4144f414e7cd393fe6f7af Mon Sep 17 00:00:00 2001 From: Remko NIJZINK Date: Fri, 20 May 2022 10:01:19 +0200 Subject: [PATCH 2/6] changed fpar calc again --- VOM_Fortran/VOM-code/transpmodel.f90 | 45 +++++++++++++++++++++------- 1 file changed, 34 insertions(+), 11 deletions(-) diff --git a/VOM_Fortran/VOM-code/transpmodel.f90 b/VOM_Fortran/VOM-code/transpmodel.f90 index 1af9d84..6a082d1 100644 --- a/VOM_Fortran/VOM-code/transpmodel.f90 +++ b/VOM_Fortran/VOM-code/transpmodel.f90 @@ -1563,7 +1563,13 @@ subroutine vom_gstom () INTEGER:: jj INTEGER:: ilay REAL*8 :: lag_sun(3) - REAL*8 :: lat_sun(3) + REAL*8 :: lat_sun(3) + REAL*8 :: fpar_sunt + REAL*8 :: fpar_sung + REAL*8 :: fpar_shadet + REAL*8 :: fpar_shadeg + + REAL*8 :: kappa select case(i_lai_function) @@ -1659,17 +1665,25 @@ subroutine vom_gstom () & / jmaxt_h(:))) * jmaxt_h(:) ! (3.23), (Out[311]) case(3) ! shaded and sunlit, with diffuse and direct radiation + + fpar_sunt = 1.0d0 - p_E ** (-lat_sun(ii) * kappa * sqrt(i_alpha_abs) ) + fpar_shadet = 1.0d0 - p_E ** (- ( lai_lt(ii) - lat_sun(ii)) * kappa * sqrt(i_alpha_abs) ) + jactt(:,ii) = ( (1.d0 - p_E ** (-(i_alpha * (pardir_h(th_) + pardiff_h(th_)) ) & - & / jmaxt_h(:))) * jmaxt_h(:) * frac_sunt(ii) ) + & + & / jmaxt_h(:))) * jmaxt_h(:) * fpar_sunt) + & & ( (1.d0 - p_E ** (-(i_alpha * pardiff_h(th_) ) & - & / jmaxt_h(:))) * jmaxt_h(:) * frac_shadet(ii) ) + & / jmaxt_h(:))) * jmaxt_h(:) * fpar_shadet ) - case(4) ! shaded and sunlit, diffuse and direct radiation, different jact-values shaded-sunlit + case(4) ! shaded and sunlit, diffuse and direct radiation, different jact-values shaded-sunlit + + fpar_sunt = 1.0d0 - p_E ** (-lat_sun(ii) * kappa * sqrt(i_alpha_abs) ) + fpar_shadet = 1.0d0 - p_E ** (- ( lai_lt(ii) - lat_sun(ii)) * kappa * sqrt(i_alpha_abs) ) + jactt(:,ii) = (1.d0 - p_E ** (-(i_alpha * (pardir_h(th_) + pardiff_h(th_)) ) & - & / jmaxt_h(:))) * jmaxt_h(:) * frac_sunt(ii) + & / jmaxt_h(:))) * jmaxt_h(:) * fpar_sunt jactts(:,ii) = (1.d0 - p_E ** (-(i_alpha * pardiff_h(th_) ) & - & / jmaxts_h(:))) * jmaxts_h(:) * frac_shadet(ii) + & / jmaxts_h(:))) * jmaxts_h(:) * fpar_shadet end select end do @@ -1688,18 +1702,27 @@ subroutine vom_gstom () jactg(:,ii) = (1.d0 - p_E ** (-(i_alpha * fpar_lg(ii) * par_h(th_)) & & / jmaxg_h(:))) * jmaxg_h(:) ! (3.23), (Out[311]) - case(3) ! shaded and sunlit, with diffuse and direct radiation + case(3) ! shaded and sunlit, with diffuse and direct radiation. + + fpar_sung = 1.0d0 - p_E ** (-lag_sun(ii) * kappa * sqrt(i_alpha_abs) ) + fpar_shadeg = 1.0d0 - p_E ** (- ( lai_lg(ii) - lag_sun(ii)) * kappa * sqrt(i_alpha_abs) ) + jactg(:,ii) = ( (1.d0 - p_E ** (-(i_alpha * (pardir_h(th_) + pardiff_h(th_)) ) & - & / jmaxg_h(:))) * jmaxg_h(:) * frac_sung(ii) ) + & + & / jmaxg_h(:))) * jmaxg_h(:) * fpar_sung ) + & & ( (1.d0 - p_E ** (-(i_alpha * pardiff_h(th_) ) & - & / jmaxg_h(:))) * jmaxg_h(:) * frac_shadeg(ii) ) + & / jmaxg_h(:))) * jmaxg_h(:) * fpar_shadeg ) case(4) ! shaded and sunlit, diffuse and direct radiation, different jact-values shaded-sunlit + + fpar_sung = 1.0d0 - p_E ** (-lag_sun(ii) * kappa * sqrt(i_alpha_abs) ) + fpar_shadeg = 1.0d0 - p_E ** (- ( lai_lg(ii) - lag_sun(ii)) * kappa * sqrt(i_alpha_abs) ) + + jactg(:,ii) = (1.d0 - p_E ** (-(i_alpha * (pardir_h(th_) + pardiff_h(th_)) ) & - & / jmaxg_h(:))) * jmaxg_h(:) * frac_sung(ii) + & / jmaxg_h(:))) * jmaxg_h(:) * fpar_sung jactgs(:,ii) = (1.d0 - p_E ** (-(i_alpha * pardiff_h(th_) ) & - & / jmaxgs_h(:))) * jmaxgs_h(:) * frac_shadeg(ii) + & / jmaxgs_h(:))) * jmaxgs_h(:) * fpar_shadeg end select From 5e200eeac8af397a5fcedf93da39911afe6b54e4 Mon Sep 17 00:00:00 2001 From: Remko NIJZINK Date: Fri, 20 May 2022 16:42:45 +0200 Subject: [PATCH 3/6] added new calculation for PARsun and PARshade --- VOM_Fortran/VOM-code/transpmodel.f90 | 152 +++++++++++++++++++++------ 1 file changed, 117 insertions(+), 35 deletions(-) diff --git a/VOM_Fortran/VOM-code/transpmodel.f90 b/VOM_Fortran/VOM-code/transpmodel.f90 index 6a082d1..9feab3c 100644 --- a/VOM_Fortran/VOM-code/transpmodel.f90 +++ b/VOM_Fortran/VOM-code/transpmodel.f90 @@ -1570,7 +1570,13 @@ subroutine vom_gstom () REAL*8 :: fpar_shadeg - REAL*8 :: kappa + REAL*8 :: kappa_t + REAL*8 :: kappa_g + REAL*8 :: PARsun + REAL*8 :: PARshade + REAL*8 :: par_shade_avg + REAL*8 :: par_sct + select case(i_lai_function) case(1) !no dynamic LAI, no shaded/sunlit fractions @@ -1587,15 +1593,15 @@ subroutine vom_gstom () case(2) !dynamic LAI, no shaded/sunlit fractions ! * extinction coefficient (Xiao et al. (2015) eq.6, Campbell and Norman (1998) eq. 15.4) - kappa = sqrt(i_chi_t**2+tan(phi_zenith(th_))**2)/(i_chi_t+1.774*(i_chi_t+1.182)**(-0.733) ) + kappa_t = sqrt(i_chi_t**2+tan(phi_zenith(th_))**2)/(i_chi_t+1.774*(i_chi_t+1.182)**(-0.733) ) ! * fraction of absorbed radiation per crown area (Beer-lambert, Xiao et al. (2015) eq.5, Campbell and Norman (1998) eq. 15.6) - fpar_lt(:) = 1.0d0 - p_E ** (-lai_lt(:) * kappa * sqrt(i_alpha_abs) ) + fpar_lt(:) = 1.0d0 - p_E ** (-lai_lt(:) * kappa_t * sqrt(i_alpha_abs) ) ! * extinction coefficient, Xiao et al. (2015) - kappa = sqrt(i_chi_g**2+tan(phi_zenith(th_))**2)/(i_chi_g+1.774*(i_chi_g+1.182)**(-0.733) ) + kappa_g = sqrt(i_chi_g**2+tan(phi_zenith(th_))**2)/(i_chi_g+1.774*(i_chi_g+1.182)**(-0.733) ) ! * fraction of absorbed radiation per crown area grasses (Beer-lambert) - fpar_lg(:) = 1.0d0 - p_E ** (-lai_lg(:) * kappa * sqrt(i_alpha_abs) ) + fpar_lg(:) = 1.0d0 - p_E ** (-lai_lg(:) * kappa_g * sqrt(i_alpha_abs) ) ! * frac_shadeg(:) = 0.0d0 @@ -1608,22 +1614,23 @@ subroutine vom_gstom () case(3, 4) !dynamic LAI, with shaded/sunlit fractions ! * extinction coefficient (Xiao et al. (2015) eq.6, Campbell and Norman (1998) eq. 15.4) - kappa = sqrt(i_chi_t**2+tan(phi_zenith(th_))**2)/(i_chi_t+1.774*(i_chi_t+1.182)**(-0.733) ) + kappa_t = sqrt(i_chi_t**2+tan(phi_zenith(th_))**2)/(i_chi_t+1.774*(i_chi_t+1.182)**(-0.733) ) ! * fraction of absorbed radiation per crown area (Beer-lambert, Xiao et al. (2015) eq.5, Campbell and Norman (1998) eq. 15.6) - fpar_lt(:) = 1.0d0 - p_E ** (-lai_lt(:) * kappa * sqrt(i_alpha_abs) ) + fpar_lt(:) = 1.0d0 - p_E ** (-lai_lt(:) * kappa_t * sqrt(i_alpha_abs) ) ! * extinction coefficient, Xiao et al. (2015) - kappa = sqrt(i_chi_g**2+tan(phi_zenith(th_))**2)/(i_chi_g+1.774*(i_chi_g+1.182)**(-0.733) ) + kappa_g = sqrt(i_chi_g**2+tan(phi_zenith(th_))**2)/(i_chi_g+1.774*(i_chi_g+1.182)**(-0.733) ) ! * fraction of absorbed radiation per crown area grasses (Beer-lambert) - fpar_lg(:) = 1.0d0 - p_E ** (-lai_lg(:) * kappa * sqrt(i_alpha_abs) ) + fpar_lg(:) = 1.0d0 - p_E ** (-lai_lg(:) * kappa_g * sqrt(i_alpha_abs) ) + lag_sun(:) = 0.0d0 ! * estimate shaded and sunlit fractions do ii = 1,3 !Eq15.23 from Campbell and Norman (1998) - lag_sun(ii) = (1.0d0 - p_E**( -lai_lg(ii) * kappa ) )/kappa + lag_sun(ii) = (1.0d0 - p_E**( -lai_lg(ii) * kappa_g ) )/kappa_g !fractions of shade and sunlit leaf areas frac_sung(ii) = MIN( (lag_sun(ii) / lai_lg(ii)), 1.0d0) @@ -1635,7 +1642,7 @@ subroutine vom_gstom () do ii = 1,3 !Eq15.23 from Campbell and Norman (1998) - lat_sun(ii) = (1.0d0 - p_E**( -lai_lt(ii) * kappa ) )/kappa + lat_sun(ii) = (1.0d0 - p_E**( -lai_lt(ii) * kappa_t ) )/kappa_t !fractions of shade and sunlit leaf areas frac_sunt(ii) = MIN( (lat_sun(ii) / lai_lt(ii) ), 1.0d0) @@ -1666,24 +1673,62 @@ subroutine vom_gstom () case(3) ! shaded and sunlit, with diffuse and direct radiation - fpar_sunt = 1.0d0 - p_E ** (-lat_sun(ii) * kappa * sqrt(i_alpha_abs) ) - fpar_shadet = 1.0d0 - p_E ** (- ( lai_lt(ii) - lat_sun(ii)) * kappa * sqrt(i_alpha_abs) ) + !---------------------------------------------------------- + !determine average PAR on shaded leaves + + !at the top of the canopy, shaded leaves receive PARdiff + !at the bottom of canopy, shaded leaves receive + !take the exponential weighted average of top and bottom (p260, Campbell&Norman): + par_shade_avg = (pardiff_h(th_) * ( 1.d0 - p_E ** ( -sqrt(i_alpha_abs) * kappa_t * lai_lt(ii) ) ) ) / & + ( sqrt(i_alpha_abs) * kappa_t * lai_lt(ii) ) + + !need to add the scattered direct PAR to the shaded PAR: + par_sct = 0.5* ( pardir_h(th_) * p_E ** (-sqrt(i_alpha_abs) * kappa_t * lai_lt(ii) ) - & + pardir_h(th_) * p_E ** (-kappa_t * lai_lt(ii) ) ) + + !final PAR-shade is the sum, multiplied with the absorptivity + PARshade = (par_sct + par_shade_avg) * i_alpha_abs + + !---------------------------------------------------------- + !sunlit leaves + PARsun = i_alpha_abs * (kappa_t * pardir_h(th_) ) + PARshade - jactt(:,ii) = ( (1.d0 - p_E ** (-(i_alpha * (pardir_h(th_) + pardiff_h(th_)) ) & - & / jmaxt_h(:))) * jmaxt_h(:) * fpar_sunt) + & - & ( (1.d0 - p_E ** (-(i_alpha * pardiff_h(th_) ) & - & / jmaxt_h(:))) * jmaxt_h(:) * fpar_shadet ) + jactt(:,ii) = ( (1.d0 - p_E ** (-(i_alpha * PARsun ) & + & / jmaxt_h(:))) * jmaxt_h(:) * lat_sun(ii) ) + & + & ( (1.d0 - p_E ** (-(i_alpha * PARshade ) & + & / jmaxt_h(:))) * jmaxt_h(:) * (lai_lt(ii) - lat_sun(ii) ) ) case(4) ! shaded and sunlit, diffuse and direct radiation, different jact-values shaded-sunlit - fpar_sunt = 1.0d0 - p_E ** (-lat_sun(ii) * kappa * sqrt(i_alpha_abs) ) - fpar_shadet = 1.0d0 - p_E ** (- ( lai_lt(ii) - lat_sun(ii)) * kappa * sqrt(i_alpha_abs) ) + + !---------------------------------------------------------- + !determine average PAR on shaded leaves + + !at the top of the canopy, shaded leaves receive PARdiff + !at the bottom of canopy, shaded leaves receive + !take the exponential weighted average of top and bottom (p260, Campbell&Norman): + par_shade_avg = (pardiff_h(th_) * ( 1.d0 - p_E ** ( -sqrt(i_alpha_abs) * kappa_t * lai_lt(ii) ) ) ) / & + ( sqrt(i_alpha_abs) * kappa_t * lai_lt(ii) ) + + !need to add the scattered direct PAR to the shaded PAR: + par_sct = 0.5* ( pardir_h(th_) * p_E ** (-sqrt(i_alpha_abs) * kappa_t * lai_lt(ii) ) - & + pardir_h(th_) * p_E ** (-kappa_t * lai_lt(ii) ) ) + + !final PAR-shade is the sum, multiplied with the absorptivity + PARshade = (par_sct + par_shade_avg) * i_alpha_abs + + !---------------------------------------------------------- + !sunlit leaves + PARsun = i_alpha_abs * (kappa_t * pardir_h(th_) ) + PARshade + - jactt(:,ii) = (1.d0 - p_E ** (-(i_alpha * (pardir_h(th_) + pardiff_h(th_)) ) & - & / jmaxt_h(:))) * jmaxt_h(:) * fpar_sunt + jactt(:,ii) = (1.d0 - p_E ** (-(i_alpha * PARsun ) & + & / jmaxt_h(:))) * jmaxt_h(:) * lat_sun(ii) + + jactts(:,ii) = (1.d0 - p_E ** (-(i_alpha * PARshade ) & + & / jmaxts_h(:))) * jmaxts_h(:) * (lai_lt(ii) - lat_sun(ii) ) + - jactts(:,ii) = (1.d0 - p_E ** (-(i_alpha * pardiff_h(th_) ) & - & / jmaxts_h(:))) * jmaxts_h(:) * fpar_shadet end select end do @@ -1704,25 +1749,62 @@ subroutine vom_gstom () case(3) ! shaded and sunlit, with diffuse and direct radiation. - fpar_sung = 1.0d0 - p_E ** (-lag_sun(ii) * kappa * sqrt(i_alpha_abs) ) - fpar_shadeg = 1.0d0 - p_E ** (- ( lai_lg(ii) - lag_sun(ii)) * kappa * sqrt(i_alpha_abs) ) + !---------------------------------------------------------- + !determine average PAR on shaded leaves + + !at the top of the canopy, shaded leaves receive PARdiff + !at the bottom of canopy, shaded leaves receive + !take the exponential weighted average of top and bottom (p260, Campbell&Norman): + par_shade_avg = (pardiff_h(th_) * ( 1.d0 - p_E ** ( -sqrt(i_alpha_abs) * kappa_g * lai_lg(ii) ) ) ) / & + ( sqrt(i_alpha_abs) * kappa_g * lai_lg(ii) ) + + !need to add the scattered direct PAR to the shaded PAR: + par_sct = 0.5* ( pardir_h(th_) * p_E ** (-sqrt(i_alpha_abs) * kappa_g * lai_lg(ii) ) - & + pardir_h(th_) * p_E ** (-kappa_g * lai_lg(ii) ) ) + + !final PAR-shade is the sum, multiplied with the absorptivity + PARshade = (par_sct + par_shade_avg) * i_alpha_abs - jactg(:,ii) = ( (1.d0 - p_E ** (-(i_alpha * (pardir_h(th_) + pardiff_h(th_)) ) & - & / jmaxg_h(:))) * jmaxg_h(:) * fpar_sung ) + & - & ( (1.d0 - p_E ** (-(i_alpha * pardiff_h(th_) ) & - & / jmaxg_h(:))) * jmaxg_h(:) * fpar_shadeg ) + !---------------------------------------------------------- + !sunlit leaves + PARsun = i_alpha_abs * (kappa_g * pardir_h(th_) ) + PARshade + + jactg(:,ii) = ( (1.d0 - p_E ** (-(i_alpha * PARsun ) & + & / jmaxg_h(:))) * jmaxg_h(:) * lag_sun(ii) ) + & + & ( (1.d0 - p_E ** (-(i_alpha * PARshade ) & + & / jmaxg_h(:))) * jmaxg_h(:) * ( lai_lg(ii) - lat_sun(ii) ) ) case(4) ! shaded and sunlit, diffuse and direct radiation, different jact-values shaded-sunlit - fpar_sung = 1.0d0 - p_E ** (-lag_sun(ii) * kappa * sqrt(i_alpha_abs) ) - fpar_shadeg = 1.0d0 - p_E ** (- ( lai_lg(ii) - lag_sun(ii)) * kappa * sqrt(i_alpha_abs) ) + !---------------------------------------------------------- + !determine average PAR on shaded leaves + + !at the top of the canopy, shaded leaves receive PARdiff + !at the bottom of canopy, shaded leaves receive + !take the exponential weighted average of top and bottom (p260, Campbell&Norman): + + par_shade_avg = (pardiff_h(th_) * ( 1.d0 - p_E ** ( -sqrt(i_alpha_abs) * kappa_g * lai_lg(ii) ) ) ) / & + ( sqrt(i_alpha_abs) * kappa_g * lai_lg(ii) ) + !need to add the scattered direct PAR to the shaded PAR: + par_sct = 0.5* ( pardir_h(th_) * p_E ** (-sqrt(i_alpha_abs) * kappa_g * lai_lg(ii) ) - & + pardir_h(th_) * p_E ** (-kappa_g * lai_lg(ii) ) ) + + !final PAR-shade is the sum, multiplied with the absorptivity + PARshade = (par_sct + par_shade_avg) * i_alpha_abs + + + + !---------------------------------------------------------- + !sunlit leaves + PARsun = i_alpha_abs * (kappa_g * pardir_h(th_) ) + PARshade + - jactg(:,ii) = (1.d0 - p_E ** (-(i_alpha * (pardir_h(th_) + pardiff_h(th_)) ) & - & / jmaxg_h(:))) * jmaxg_h(:) * fpar_sung + jactg(:,ii) = (1.d0 - p_E ** (-(i_alpha * PARsun ) & + & / jmaxg_h(:))) * jmaxg_h(:) * lag_sun(ii) - jactgs(:,ii) = (1.d0 - p_E ** (-(i_alpha * pardiff_h(th_) ) & - & / jmaxgs_h(:))) * jmaxgs_h(:) * fpar_shadeg + jactgs(:,ii) = (1.d0 - p_E ** (-(i_alpha * PARshade ) & + & / jmaxgs_h(:))) * jmaxgs_h(:) * ( lai_lg(ii) - lat_sun(ii) ) end select From 2e026324bc44aec70aaccf2dc69e8ca5a3444d02 Mon Sep 17 00:00:00 2001 From: Remko NIJZINK Date: Fri, 20 May 2022 21:23:40 +0200 Subject: [PATCH 4/6] corrected difffuse extinction coefficient --- VOM_Fortran/VOM-code/transpmodel.f90 | 19 +++++++++++-------- 1 file changed, 11 insertions(+), 8 deletions(-) diff --git a/VOM_Fortran/VOM-code/transpmodel.f90 b/VOM_Fortran/VOM-code/transpmodel.f90 index 9feab3c..c61e7bb 100644 --- a/VOM_Fortran/VOM-code/transpmodel.f90 +++ b/VOM_Fortran/VOM-code/transpmodel.f90 @@ -1576,6 +1576,7 @@ subroutine vom_gstom () REAL*8 :: PARshade REAL*8 :: par_shade_avg REAL*8 :: par_sct + REAL*8 :: kappa_d select case(i_lai_function) @@ -1613,6 +1614,8 @@ subroutine vom_gstom () case(3, 4) !dynamic LAI, with shaded/sunlit fractions + kappa_d = 0.719d0 + ! * extinction coefficient (Xiao et al. (2015) eq.6, Campbell and Norman (1998) eq. 15.4) kappa_t = sqrt(i_chi_t**2+tan(phi_zenith(th_))**2)/(i_chi_t+1.774*(i_chi_t+1.182)**(-0.733) ) @@ -1679,8 +1682,8 @@ subroutine vom_gstom () !at the top of the canopy, shaded leaves receive PARdiff !at the bottom of canopy, shaded leaves receive !take the exponential weighted average of top and bottom (p260, Campbell&Norman): - par_shade_avg = (pardiff_h(th_) * ( 1.d0 - p_E ** ( -sqrt(i_alpha_abs) * kappa_t * lai_lt(ii) ) ) ) / & - ( sqrt(i_alpha_abs) * kappa_t * lai_lt(ii) ) + par_shade_avg = (pardiff_h(th_) * ( 1.d0 - p_E ** ( -sqrt(i_alpha_abs) * kappa_d * lai_lt(ii) ) ) ) / & + ( sqrt(i_alpha_abs) * kappa_d* lai_lt(ii) ) !need to add the scattered direct PAR to the shaded PAR: par_sct = 0.5* ( pardir_h(th_) * p_E ** (-sqrt(i_alpha_abs) * kappa_t * lai_lt(ii) ) - & @@ -1707,8 +1710,8 @@ subroutine vom_gstom () !at the top of the canopy, shaded leaves receive PARdiff !at the bottom of canopy, shaded leaves receive !take the exponential weighted average of top and bottom (p260, Campbell&Norman): - par_shade_avg = (pardiff_h(th_) * ( 1.d0 - p_E ** ( -sqrt(i_alpha_abs) * kappa_t * lai_lt(ii) ) ) ) / & - ( sqrt(i_alpha_abs) * kappa_t * lai_lt(ii) ) + par_shade_avg = (pardiff_h(th_) * ( 1.d0 - p_E ** ( -sqrt(i_alpha_abs) * kappa_d * lai_lt(ii) ) ) ) / & + ( sqrt(i_alpha_abs) * kappa_d * lai_lt(ii) ) !need to add the scattered direct PAR to the shaded PAR: par_sct = 0.5* ( pardir_h(th_) * p_E ** (-sqrt(i_alpha_abs) * kappa_t * lai_lt(ii) ) - & @@ -1755,8 +1758,8 @@ subroutine vom_gstom () !at the top of the canopy, shaded leaves receive PARdiff !at the bottom of canopy, shaded leaves receive !take the exponential weighted average of top and bottom (p260, Campbell&Norman): - par_shade_avg = (pardiff_h(th_) * ( 1.d0 - p_E ** ( -sqrt(i_alpha_abs) * kappa_g * lai_lg(ii) ) ) ) / & - ( sqrt(i_alpha_abs) * kappa_g * lai_lg(ii) ) + par_shade_avg = (pardiff_h(th_) * ( 1.d0 - p_E ** ( -sqrt(i_alpha_abs) * kappa_d * lai_lg(ii) ) ) ) / & + ( sqrt(i_alpha_abs) * kappa_d * lai_lg(ii) ) !need to add the scattered direct PAR to the shaded PAR: par_sct = 0.5* ( pardir_h(th_) * p_E ** (-sqrt(i_alpha_abs) * kappa_g * lai_lg(ii) ) - & @@ -1783,8 +1786,8 @@ subroutine vom_gstom () !at the bottom of canopy, shaded leaves receive !take the exponential weighted average of top and bottom (p260, Campbell&Norman): - par_shade_avg = (pardiff_h(th_) * ( 1.d0 - p_E ** ( -sqrt(i_alpha_abs) * kappa_g * lai_lg(ii) ) ) ) / & - ( sqrt(i_alpha_abs) * kappa_g * lai_lg(ii) ) + par_shade_avg = (pardiff_h(th_) * ( 1.d0 - p_E ** ( -sqrt(i_alpha_abs) * kappa_d * lai_lg(ii) ) ) ) / & + ( sqrt(i_alpha_abs) * kappa_d * lai_lg(ii) ) !need to add the scattered direct PAR to the shaded PAR: par_sct = 0.5* ( pardir_h(th_) * p_E ** (-sqrt(i_alpha_abs) * kappa_g * lai_lg(ii) ) - & From 9cf8aacaa6a5a0637dc541cbf4b8ae1f1cbb6038 Mon Sep 17 00:00:00 2001 From: Remko NIJZINK Date: Thu, 2 Jun 2022 11:28:26 +0200 Subject: [PATCH 5/6] corrected PARsun calculation (brackets) --- VOM_Fortran/VOM-code/transpmodel.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/VOM_Fortran/VOM-code/transpmodel.f90 b/VOM_Fortran/VOM-code/transpmodel.f90 index c61e7bb..f80de3f 100644 --- a/VOM_Fortran/VOM-code/transpmodel.f90 +++ b/VOM_Fortran/VOM-code/transpmodel.f90 @@ -1694,7 +1694,7 @@ subroutine vom_gstom () !---------------------------------------------------------- !sunlit leaves - PARsun = i_alpha_abs * (kappa_t * pardir_h(th_) ) + PARshade + PARsun = i_alpha_abs * (kappa_t * pardir_h(th_) + PARshade ) jactt(:,ii) = ( (1.d0 - p_E ** (-(i_alpha * PARsun ) & & / jmaxt_h(:))) * jmaxt_h(:) * lat_sun(ii) ) + & @@ -1722,7 +1722,7 @@ subroutine vom_gstom () !---------------------------------------------------------- !sunlit leaves - PARsun = i_alpha_abs * (kappa_t * pardir_h(th_) ) + PARshade + PARsun = i_alpha_abs * (kappa_t * pardir_h(th_) + PARshade ) jactt(:,ii) = (1.d0 - p_E ** (-(i_alpha * PARsun ) & From f6e816ed31e75649c0d176af0f6c50edbf7c1799 Mon Sep 17 00:00:00 2001 From: Remko NIJZINK Date: Thu, 2 Jun 2022 16:24:09 +0200 Subject: [PATCH 6/6] corrected PARshade grasses --- VOM_Fortran/VOM-code/transpmodel.f90 | 7 +++---- 1 file changed, 3 insertions(+), 4 deletions(-) diff --git a/VOM_Fortran/VOM-code/transpmodel.f90 b/VOM_Fortran/VOM-code/transpmodel.f90 index f80de3f..56ade5a 100644 --- a/VOM_Fortran/VOM-code/transpmodel.f90 +++ b/VOM_Fortran/VOM-code/transpmodel.f90 @@ -1770,8 +1770,8 @@ subroutine vom_gstom () !---------------------------------------------------------- !sunlit leaves - PARsun = i_alpha_abs * (kappa_g * pardir_h(th_) ) + PARshade - + PARsun = i_alpha_abs * (kappa_t * pardir_h(th_) + PARshade ) + jactg(:,ii) = ( (1.d0 - p_E ** (-(i_alpha * PARsun ) & & / jmaxg_h(:))) * jmaxg_h(:) * lag_sun(ii) ) + & & ( (1.d0 - p_E ** (-(i_alpha * PARshade ) & @@ -1800,8 +1800,7 @@ subroutine vom_gstom () !---------------------------------------------------------- !sunlit leaves - PARsun = i_alpha_abs * (kappa_g * pardir_h(th_) ) + PARshade - + PARsun = i_alpha_abs * (kappa_t * pardir_h(th_) + PARshade ) jactg(:,ii) = (1.d0 - p_E ** (-(i_alpha * PARsun ) & & / jmaxg_h(:))) * jmaxg_h(:) * lag_sun(ii)