From 1b3179fbe50fd4e02224908892273955d3f24e26 Mon Sep 17 00:00:00 2001 From: Remko NIJZINK Date: Thu, 19 May 2022 18:52:40 +0200 Subject: [PATCH 1/7] 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/7] 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/7] 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/7] 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/7] 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/7] 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) From b08c27ef833f15c0dc8c1384769034e6f0ef5d21 Mon Sep 17 00:00:00 2001 From: Remko NIJZINK Date: Fri, 15 Jul 2022 10:06:53 +0200 Subject: [PATCH 7/7] added shell scripts and edited docs --- VOM_Fortran/VOM-docu/source/quickstart.rst | 48 ++++++++++++++ VOM_Fortran/VOM-scripts/run_vom.sh | 73 ++++++++++++++++++++++ VOM_Fortran/VOM-scripts/vom_parallel.job | 69 ++++++++++++++++++++ 3 files changed, 190 insertions(+) create mode 100755 VOM_Fortran/VOM-scripts/run_vom.sh create mode 100755 VOM_Fortran/VOM-scripts/vom_parallel.job diff --git a/VOM_Fortran/VOM-docu/source/quickstart.rst b/VOM_Fortran/VOM-docu/source/quickstart.rst index a2828bb..886fe70 100755 --- a/VOM_Fortran/VOM-docu/source/quickstart.rst +++ b/VOM_Fortran/VOM-docu/source/quickstart.rst @@ -63,6 +63,54 @@ This can be changed in vom_namelist, or on the command line: -n The VOM_namelist (filename can be different) +To run successfully, these input-files are needed: + +**vom_namelist** + Contains all settings to run the VOM. + +**pars.txt** + Contains the (optimized) vegetation parameters (only needed for single run). + +**dailyweather.prn** + Contains the meteorological forcing. + +By default, the executable looks for the vom_namelist in the current workdirectory. The default directory for the other files is /input, relative to the workdirectory. +This can be changed in vom_namelist, or on the command line: + +-i Inputpath to directory with dailyweather.prn, and optionally pars.txt. + +-o Outputpath for all outputfiles. + +-n The VOM_namelist (filename can be different) + + + +Scripts for running the model +----------------------------- +In the folder VOM-scripts, two shell-scripts can be found for running the VOM: + +.. code-block:: bash + + vom_parallel.job + + + run_vom.sh + + +The model is compiled and executed by the script run_vom.sh. It can be used in the following way: + +.. code-block:: bash + + bash src_sh/run_vom.sh src/VOM/VOM_Fortran/ + + +vom_parallel.job is a job script for the LIST HPC using SLURM. Note that the script needs to be modified for own applications of the VOM on an HPC and serves here as an example. In this example, submitting a job with SLURM functions as follows: + +.. code-block:: bash + + sbatch -o -e -n -J vom_parallel.job + + Model modes ----------------- diff --git a/VOM_Fortran/VOM-scripts/run_vom.sh b/VOM_Fortran/VOM-scripts/run_vom.sh new file mode 100755 index 0000000..efed1f0 --- /dev/null +++ b/VOM_Fortran/VOM-scripts/run_vom.sh @@ -0,0 +1,73 @@ +#!/bin/sh +# coding: utf-8 +#*********************************************************************** +# run_vom.sh +# Compiles and runs the Vegetation Optimality Model (VOM). +# +# +#----------------------------------------------------------------------- +# Authors: Remko Nijzink +# Now at: LIST (Luxembourg Institute of Science and Technology) +#----------------------------------------------------------------------- +# +# Copyright (C) 2022, Luxembourg Institute of Science and Technology, all rights reserved. +# +# Code licensed under GPL version 2 or higher +# SPDX-License-Identifier: GPL-2.0-or-later +# +#*********************************************************************** + + +exe_dir=$1 +inputdir=$2 +input_weather=$3 +input_soil=$4 +nml_input=$5 +outputdir=$6 +restart_dir=$7 + + +date + +#compile code +C_INCLUDE_PATH=$C_INCLUDE_PATH:/usr/include/ +LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/usr/lib/ + +#check if paths are empty, use default otherwise +if [ ! -d "$NC_INCLUDE" ]; then + NC_INCLUDE=/usr/include + echo $NC_INCLUDE +fi + +if [ ! -d "$NC_LIB" ]; then + NC_LIB=/usr/lib + echo $NC_LIB +fi + + + +make --directory $exe_dir FC=gfortran NC_INCLUDE=$NC_INCLUDE NC_LIB=$NC_LIB + +#check if the outputdir exists and else make one +if [ ! -d "$outputdir" ]; then + mkdir $outputdir +fi + +if [ -f "$restart_dir/sce_lastloop.txt" ]; then + cp $restart_dir/* $outputdir +fi + +if [ ! -d "$inputdir" ]; then +mkdir $inputdir +cp $input_weather $inputdir +cp $input_soil $inputdir + +fi + +#run the model +$exe_dir/model.x -i $inputdir -o $outputdir -n $nml_input + +#clean again +make clean --directory $exe_dir + +date diff --git a/VOM_Fortran/VOM-scripts/vom_parallel.job b/VOM_Fortran/VOM-scripts/vom_parallel.job new file mode 100755 index 0000000..abfd8d2 --- /dev/null +++ b/VOM_Fortran/VOM-scripts/vom_parallel.job @@ -0,0 +1,69 @@ +#!/bin/sh +# coding: utf-8 +#*********************************************************************** +# vom_parallel.sh +# Jobscript for the LIST HPC, to run run_vom.sh. +# +# +#----------------------------------------------------------------------- +# Authors: Remko Nijzink +# Now at: LIST (Luxembourg Institute of Science and Technology) +#----------------------------------------------------------------------- +# +# Copyright (C) 2022, Luxembourg Institute of Science and Technology, all rights reserved. +# +# Code licensed under GPL version 2 or higher +# SPDX-License-Identifier: GPL-2.0-or-later +# +#*********************************************************************** +# +# Your job name +#SBATCH --job-name=VOMexp +# +# Use current working directory +# +# +#SBATCH -e /home/int/eva/nijzink/logfiles/err +#SBATCH -o /home/int/eva/nijzink/logfiles/out +#SBATCH --account=wave +# pe (Parallel environment) request. Set your number of processors here. +#SBATCH -n 4 +#SBATCH --time=47:00:00 +# +# The following output will show in the output file. Used for debugging. +echo "Got $SLURM_NNODES processors." +# +# load modules +module load userspace/tr17.10 +module load userspace/spack +module load gcc-8.2.0-gcc-7.2.0-la27vub +module load netcdf-fortran-4.5.2-gcc-8.2.0-y4qfiu6 +# +# to avoid an error message in renku +#export LC_ALL='en_US.UTF-8' +#export LANG='en_US.UTF-8' +#module load null shared use.own conda +# +#files to be tracked with renku +exe_src=$2 +dailyweather_input=$1 + +NC_INCLUDE=/trinity/shared/apps/spack/opt/spack/linux-centos7-x86_64/gcc-8.2.0/netcdf-fortran-4.5.2-y4qfiu6nzx3voevihng2ihu4eunwfc3x/include +NC_LIB=/trinity/shared/apps/spack/opt/spack/linux-centos7-x86_64/gcc-8.2.0/netcdf-fortran-4.5.2-y4qfiu6nzx3voevihng2ihu4eunwfc3x/lib/ + +export NC_INCLUDE +export NC_LIB + +date +source /home/int/eva/nijzink/miniconda3/etc/profile.d/conda.sh +conda activate renku + +#inputs: +#workdir +#dailyweather.prn +#vom_namelist +#VOM-code + +renku run bash src_sh/run_vom.sh $@ + +date