From 270144cb47103839bc2a4c2873a8fd270871f80a Mon Sep 17 00:00:00 2001 From: Daniel Kennedy Date: Mon, 11 Jan 2021 17:14:31 -0700 Subject: [PATCH 01/12] implement iwue with fpsn/gs --- src/biogeophys/PhotosynthesisMod.F90 | 15 +++++++++++++-- 1 file changed, 13 insertions(+), 2 deletions(-) diff --git a/src/biogeophys/PhotosynthesisMod.F90 b/src/biogeophys/PhotosynthesisMod.F90 index 44e28a681c..8886bcd74a 100644 --- a/src/biogeophys/PhotosynthesisMod.F90 +++ b/src/biogeophys/PhotosynthesisMod.F90 @@ -172,6 +172,7 @@ module PhotosynthesisMod real(r8), pointer, private :: psnsun_wp_patch (:) ! patch product-limited sunlit leaf photosynthesis (umol CO2/m**2/s) real(r8), pointer, private :: psnsha_wp_patch (:) ! patch product-limited shaded leaf photosynthesis (umol CO2/m**2/s) + real(r8), pointer, public :: iwue_patch (:) ! patch intrinsic water use efficiency (umol CO2/mol H2O) real(r8), pointer, public :: fpsn_patch (:) ! patch photosynthesis (umol CO2/m**2 ground/s) real(r8), pointer, private :: fpsn_wc_patch (:) ! patch Rubisco-limited photosynthesis (umol CO2/m**2 ground/s) real(r8), pointer, private :: fpsn_wj_patch (:) ! patch RuBP-limited photosynthesis (umol CO2/m**2 ground/s) @@ -313,7 +314,7 @@ subroutine InitAllocate(this, bounds) allocate(this%fpsn_wc_patch (begp:endp)) ; this%fpsn_wc_patch (:) = nan allocate(this%fpsn_wj_patch (begp:endp)) ; this%fpsn_wj_patch (:) = nan allocate(this%fpsn_wp_patch (begp:endp)) ; this%fpsn_wp_patch (:) = nan - + allocate(this%iwue_patch (begp:endp)) ; this%iwue_patch (:) = nan allocate(this%lnca_patch (begp:endp)) ; this%lnca_patch (:) = nan allocate(this%lmrsun_z_patch (begp:endp,1:nlevcan)) ; this%lmrsun_z_patch (:,:) = nan @@ -396,6 +397,11 @@ subroutine InitHistory(this, bounds) avgflag='A', long_name='photosynthesis', & ptr_patch=this%fpsn_patch, set_lake=0._r8, set_urb=0._r8) + this%fpsn_patch(begp:endp) = spval + call hist_addfld1d (fname='IWUE', units='umol CO2 / mol H2O', & + avgflag='A', long_name='intrinsic water usef efficiency', & + ptr_patch=this%iwue_patch, set_lake=0._r8, set_urb=0._r8) + ! Don't by default output this rate limiting step as only makes sense if you are outputing ! the others each time-step this%fpsn_wc_patch(begp:endp) = spval @@ -1849,6 +1855,8 @@ subroutine PhotosynthesisTotal (fn, filterp, & laisun => canopystate_inst%laisun_patch , & ! Input: [real(r8) (:) ] sunlit leaf area laisha => canopystate_inst%laisha_patch , & ! Input: [real(r8) (:) ] shaded leaf area + gs_mol_sun => photosyns_inst%gs_mol_sun_patch , & ! Input: [real(r8) (:) ] sunlit leaf stomatal conductance (umol H2O/m**2/s) + gs_mol_sha => photosyns_inst%gs_mol_sha_patch , & ! Input: [real(r8) (:) ] shaded leaf stomatal conductance (umol H2O/m**2/s) psnsun => photosyns_inst%psnsun_patch , & ! Input: [real(r8) (:) ] sunlit leaf photosynthesis (umol CO2 /m**2/ s) psnsha => photosyns_inst%psnsha_patch , & ! Input: [real(r8) (:) ] shaded leaf photosynthesis (umol CO2 /m**2/ s) rc13_canair => photosyns_inst%rc13_canair_patch , & ! Output: [real(r8) (:) ] C13O2/C12O2 in canopy air @@ -1869,7 +1877,8 @@ subroutine PhotosynthesisTotal (fn, filterp, & fpsn => photosyns_inst%fpsn_patch , & ! Output: [real(r8) (:) ] photosynthesis (umol CO2 /m**2 /s) fpsn_wc => photosyns_inst%fpsn_wc_patch , & ! Output: [real(r8) (:) ] Rubisco-limited photosynthesis (umol CO2 /m**2 /s) fpsn_wj => photosyns_inst%fpsn_wj_patch , & ! Output: [real(r8) (:) ] RuBP-limited photosynthesis (umol CO2 /m**2 /s) - fpsn_wp => photosyns_inst%fpsn_wp_patch & ! Output: [real(r8) (:) ] product-limited photosynthesis (umol CO2 /m**2 /s) + fpsn_wp => photosyns_inst%fpsn_wp_patch , & ! Output: [real(r8) (:) ] product-limited photosynthesis (umol CO2 /m**2 /s) + iwue => photosyns_inst%iwue_patch & ! Output: [real(r8) (:) ] intrinsic water use efficience (umol CO2 / mol H2o) ) if ( use_c14 ) then @@ -1895,6 +1904,8 @@ subroutine PhotosynthesisTotal (fn, filterp, & fpsn_wc(p) = psnsun_wc(p)*laisun(p) + psnsha_wc(p)*laisha(p) fpsn_wj(p) = psnsun_wj(p)*laisun(p) + psnsha_wj(p)*laisha(p) fpsn_wp(p) = psnsun_wp(p)*laisun(p) + psnsha_wp(p)*laisha(p) + ! calculate inherent wue, 1e6 converts gsmol from umol/m2/s to mol/m2/s + iwue(p) = 1.e06_r8*fpsn(p)/(laisun(p)*gs_mol_sun(p)+laisha(p)*gs_mol_sha(p)) end if if (use_cn) then From a4f6423684c3ce5dcd7f2671a305644d9c39ced4 Mon Sep 17 00:00:00 2001 From: Daniel Kennedy Date: Mon, 11 Jan 2021 19:41:40 -0700 Subject: [PATCH 02/12] only compute iwue daytime --- src/biogeophys/CanopyFluxesMod.F90 | 2 +- src/biogeophys/PhotosynthesisMod.F90 | 19 ++++++++++++++----- 2 files changed, 15 insertions(+), 6 deletions(-) diff --git a/src/biogeophys/CanopyFluxesMod.F90 b/src/biogeophys/CanopyFluxesMod.F90 index b96b43187c..bf18016c13 100644 --- a/src/biogeophys/CanopyFluxesMod.F90 +++ b/src/biogeophys/CanopyFluxesMod.F90 @@ -1351,7 +1351,7 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, ! Determine total photosynthesis call PhotosynthesisTotal(fn, filterp, & - atm2lnd_inst, canopystate_inst, photosyns_inst) + atm2lnd_inst, canopystate_inst, photosyns_inst, solarabs_inst) ! Calculate ozone stress. This needs to be done after rssun and rsshade are ! computed by the Photosynthesis routine. However, Photosynthesis also uses the diff --git a/src/biogeophys/PhotosynthesisMod.F90 b/src/biogeophys/PhotosynthesisMod.F90 index 8886bcd74a..ade998a084 100644 --- a/src/biogeophys/PhotosynthesisMod.F90 +++ b/src/biogeophys/PhotosynthesisMod.F90 @@ -399,7 +399,7 @@ subroutine InitHistory(this, bounds) this%fpsn_patch(begp:endp) = spval call hist_addfld1d (fname='IWUE', units='umol CO2 / mol H2O', & - avgflag='A', long_name='intrinsic water usef efficiency', & + avgflag='A', long_name='intrinsic water use efficiency', & ptr_patch=this%iwue_patch, set_lake=0._r8, set_urb=0._r8) ! Don't by default output this rate limiting step as only makes sense if you are outputing @@ -1829,7 +1829,7 @@ end subroutine Photosynthesis !------------------------------------------------------------------------------ subroutine PhotosynthesisTotal (fn, filterp, & - atm2lnd_inst, canopystate_inst, photosyns_inst) + atm2lnd_inst, canopystate_inst, photosyns_inst,solarabs_inst) ! ! Determine total photosynthesis ! @@ -1839,6 +1839,7 @@ subroutine PhotosynthesisTotal (fn, filterp, & type(atm2lnd_type) , intent(in) :: atm2lnd_inst type(canopystate_type) , intent(in) :: canopystate_inst type(photosyns_type) , intent(inout) :: photosyns_inst + type(solarabs_type) , intent(in) :: solarabs_inst ! ! !LOCAL VARIABLES: integer :: f,fp,p,l,g ! indices @@ -1855,8 +1856,10 @@ subroutine PhotosynthesisTotal (fn, filterp, & laisun => canopystate_inst%laisun_patch , & ! Input: [real(r8) (:) ] sunlit leaf area laisha => canopystate_inst%laisha_patch , & ! Input: [real(r8) (:) ] shaded leaf area - gs_mol_sun => photosyns_inst%gs_mol_sun_patch , & ! Input: [real(r8) (:) ] sunlit leaf stomatal conductance (umol H2O/m**2/s) - gs_mol_sha => photosyns_inst%gs_mol_sha_patch , & ! Input: [real(r8) (:) ] shaded leaf stomatal conductance (umol H2O/m**2/s) + par_z => solarabs_inst%parsun_z_patch , & ! Input: [real(r8) (:,:) ] par absorbed per unit lai for canopy layer (w/m**2) + + gs_mol_sun => photosyns_inst%gs_mol_sun_patch , & ! Input: [real(r8) (:,:) ] sunlit leaf stomatal conductance (umol H2O/m**2/s) + gs_mol_sha => photosyns_inst%gs_mol_sha_patch , & ! Input: [real(r8) (:,:) ] shaded leaf stomatal conductance (umol H2O/m**2/s) psnsun => photosyns_inst%psnsun_patch , & ! Input: [real(r8) (:) ] sunlit leaf photosynthesis (umol CO2 /m**2/ s) psnsha => photosyns_inst%psnsha_patch , & ! Input: [real(r8) (:) ] shaded leaf photosynthesis (umol CO2 /m**2/ s) rc13_canair => photosyns_inst%rc13_canair_patch , & ! Output: [real(r8) (:) ] C13O2/C12O2 in canopy air @@ -1905,7 +1908,13 @@ subroutine PhotosynthesisTotal (fn, filterp, & fpsn_wj(p) = psnsun_wj(p)*laisun(p) + psnsha_wj(p)*laisha(p) fpsn_wp(p) = psnsun_wp(p)*laisun(p) + psnsha_wp(p)*laisha(p) ! calculate inherent wue, 1e6 converts gsmol from umol/m2/s to mol/m2/s - iwue(p) = 1.e06_r8*fpsn(p)/(laisun(p)*gs_mol_sun(p)+laisha(p)*gs_mol_sha(p)) + ! iwue calculation not compatible with multi-layer canopy + if (par_z(p,1) > 0._r8) then !daytime + iwue(p) = 1.e06_r8*fpsn(p)/(laisun(p)*gs_mol_sun(p,1)+laisha(p)*gs_mol_sha(p,1)) + else + iwue(p) = spval + write(iulog,*) 'zqz iwue' + end if end if if (use_cn) then From c1160214ece6568103049ce7a2998ab22e7e5533 Mon Sep 17 00:00:00 2001 From: Daniel Kennedy Date: Thu, 14 Jan 2021 10:34:29 -0700 Subject: [PATCH 03/12] working vpd2m --- src/biogeophys/CanopyFluxesMod.F90 | 8 +++++++- src/biogeophys/WaterDiagnosticBulkType.F90 | 12 +++++++++++- 2 files changed, 18 insertions(+), 2 deletions(-) diff --git a/src/biogeophys/CanopyFluxesMod.F90 b/src/biogeophys/CanopyFluxesMod.F90 index bf18016c13..bfb1ea245d 100644 --- a/src/biogeophys/CanopyFluxesMod.F90 +++ b/src/biogeophys/CanopyFluxesMod.F90 @@ -527,6 +527,7 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, q_ref2m => waterdiagnosticbulk_inst%q_ref2m_patch , & ! Output: [real(r8) (:) ] 2 m height surface specific humidity (kg/kg) rh_ref2m_r => waterdiagnosticbulk_inst%rh_ref2m_r_patch , & ! Output: [real(r8) (:) ] Rural 2 m height surface relative humidity (%) rh_ref2m => waterdiagnosticbulk_inst%rh_ref2m_patch , & ! Output: [real(r8) (:) ] 2 m height surface relative humidity (%) + vpd_ref2m => waterdiagnosticbulk_inst%vpd_ref2m_patch , & ! Output: [real(r8) (:) ] 2 m height surface vapor pressure deficit (Pa) rhaf => waterdiagnosticbulk_inst%rh_af_patch , & ! Output: [real(r8) (:) ] fractional humidity of canopy air [dimensionless] qflx_tran_veg => waterfluxbulk_inst%qflx_tran_veg_patch , & ! Output: [real(r8) (:) ] vegetation transpiration (mm H2O/s) (+ = to atm) @@ -1256,6 +1257,12 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, rh_ref2m(p) = min(100._r8, q_ref2m(p) / qsat_ref2m * 100._r8) rh_ref2m_r(p) = rh_ref2m(p) + ! 2 m vapor pressure deficit + vpd_ref2m(p) = e_ref2m*(1-rh_ref2m(p)/100._r8) + + + + ! Human Heat Stress if ( all_human_stress_indices .or. fast_human_stress_indices ) then call KtoC(t_ref2m(p), tc_ref2m(p)) @@ -1434,7 +1441,6 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, end associate - end subroutine CanopyFluxes end module CanopyFluxesMod diff --git a/src/biogeophys/WaterDiagnosticBulkType.F90 b/src/biogeophys/WaterDiagnosticBulkType.F90 index 21cc9d283b..a5d739b7b4 100644 --- a/src/biogeophys/WaterDiagnosticBulkType.F90 +++ b/src/biogeophys/WaterDiagnosticBulkType.F90 @@ -53,6 +53,7 @@ module WaterDiagnosticBulkType real(r8), pointer :: h2osno_top_col (:) ! col top-layer mass of snow [kg] real(r8), pointer :: sno_liq_top_col (:) ! col snow liquid water fraction (mass), top layer [fraction] + real(r8), pointer :: vpd_ref2m_patch (:) ! patch 2 m height surface vapor pressue deficit (Pa) real(r8), pointer :: rh_ref2m_patch (:) ! patch 2 m height surface relative humidity (%) real(r8), pointer :: rh_ref2m_r_patch (:) ! patch 2 m height surface relative humidity - rural (%) real(r8), pointer :: rh_ref2m_u_patch (:) ! patch 2 m height surface relative humidity - urban (%) @@ -190,7 +191,8 @@ subroutine InitBulkAllocate(this, bounds) allocate(this%h2osno_top_col (begc:endc)) ; this%h2osno_top_col (:) = nan allocate(this%sno_liq_top_col (begc:endc)) ; this%sno_liq_top_col (:) = nan - allocate(this%dqgdT_col (begc:endc)) ; this%dqgdT_col (:) = nan + allocate(this%dqgdT_col (begc:endc)) ; this%dqgdT_col (:) = nan + allocate(this%vpd_ref2m_patch (begp:endp)) ; this%vpd_ref2m_patch (:) = nan allocate(this%rh_ref2m_patch (begp:endp)) ; this%rh_ref2m_patch (:) = nan allocate(this%rh_ref2m_u_patch (begp:endp)) ; this%rh_ref2m_u_patch (:) = nan allocate(this%rh_ref2m_r_patch (begp:endp)) ; this%rh_ref2m_r_patch (:) = nan @@ -268,6 +270,14 @@ subroutine InitBulkHistory(this, bounds) long_name=this%info%lname('vertically summed soil cie (veg landunits only)'), & ptr_col=this%h2osoi_ice_tot_col, l2g_scale_type='veg') + this%vpd_ref2m_patch(begp:endp) = spval + call hist_addfld1d ( & + fname=this%info%fname('VPD2M'), & + units='Pa', & + avgflag='A', & + long_name=this%info%lname('2m vapor pressure deficit'), & + ptr_patch=this%vpd_ref2m_patch) + this%rh_ref2m_patch(begp:endp) = spval call hist_addfld1d ( & fname=this%info%fname('RH2M'), & From db253e38081af16e75a8adb1951f970b3bff20e6 Mon Sep 17 00:00:00 2001 From: Daniel Kennedy Date: Tue, 19 Jan 2021 08:17:36 -0700 Subject: [PATCH 04/12] starting over --- src/biogeophys/CanopyFluxesMod.F90 | 10 ++------- src/biogeophys/PhotosynthesisMod.F90 | 26 +++------------------- src/biogeophys/WaterDiagnosticBulkType.F90 | 12 +--------- 3 files changed, 6 insertions(+), 42 deletions(-) diff --git a/src/biogeophys/CanopyFluxesMod.F90 b/src/biogeophys/CanopyFluxesMod.F90 index bfb1ea245d..b96b43187c 100644 --- a/src/biogeophys/CanopyFluxesMod.F90 +++ b/src/biogeophys/CanopyFluxesMod.F90 @@ -527,7 +527,6 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, q_ref2m => waterdiagnosticbulk_inst%q_ref2m_patch , & ! Output: [real(r8) (:) ] 2 m height surface specific humidity (kg/kg) rh_ref2m_r => waterdiagnosticbulk_inst%rh_ref2m_r_patch , & ! Output: [real(r8) (:) ] Rural 2 m height surface relative humidity (%) rh_ref2m => waterdiagnosticbulk_inst%rh_ref2m_patch , & ! Output: [real(r8) (:) ] 2 m height surface relative humidity (%) - vpd_ref2m => waterdiagnosticbulk_inst%vpd_ref2m_patch , & ! Output: [real(r8) (:) ] 2 m height surface vapor pressure deficit (Pa) rhaf => waterdiagnosticbulk_inst%rh_af_patch , & ! Output: [real(r8) (:) ] fractional humidity of canopy air [dimensionless] qflx_tran_veg => waterfluxbulk_inst%qflx_tran_veg_patch , & ! Output: [real(r8) (:) ] vegetation transpiration (mm H2O/s) (+ = to atm) @@ -1257,12 +1256,6 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, rh_ref2m(p) = min(100._r8, q_ref2m(p) / qsat_ref2m * 100._r8) rh_ref2m_r(p) = rh_ref2m(p) - ! 2 m vapor pressure deficit - vpd_ref2m(p) = e_ref2m*(1-rh_ref2m(p)/100._r8) - - - - ! Human Heat Stress if ( all_human_stress_indices .or. fast_human_stress_indices ) then call KtoC(t_ref2m(p), tc_ref2m(p)) @@ -1358,7 +1351,7 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, ! Determine total photosynthesis call PhotosynthesisTotal(fn, filterp, & - atm2lnd_inst, canopystate_inst, photosyns_inst, solarabs_inst) + atm2lnd_inst, canopystate_inst, photosyns_inst) ! Calculate ozone stress. This needs to be done after rssun and rsshade are ! computed by the Photosynthesis routine. However, Photosynthesis also uses the @@ -1441,6 +1434,7 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, end associate + end subroutine CanopyFluxes end module CanopyFluxesMod diff --git a/src/biogeophys/PhotosynthesisMod.F90 b/src/biogeophys/PhotosynthesisMod.F90 index ade998a084..44e28a681c 100644 --- a/src/biogeophys/PhotosynthesisMod.F90 +++ b/src/biogeophys/PhotosynthesisMod.F90 @@ -172,7 +172,6 @@ module PhotosynthesisMod real(r8), pointer, private :: psnsun_wp_patch (:) ! patch product-limited sunlit leaf photosynthesis (umol CO2/m**2/s) real(r8), pointer, private :: psnsha_wp_patch (:) ! patch product-limited shaded leaf photosynthesis (umol CO2/m**2/s) - real(r8), pointer, public :: iwue_patch (:) ! patch intrinsic water use efficiency (umol CO2/mol H2O) real(r8), pointer, public :: fpsn_patch (:) ! patch photosynthesis (umol CO2/m**2 ground/s) real(r8), pointer, private :: fpsn_wc_patch (:) ! patch Rubisco-limited photosynthesis (umol CO2/m**2 ground/s) real(r8), pointer, private :: fpsn_wj_patch (:) ! patch RuBP-limited photosynthesis (umol CO2/m**2 ground/s) @@ -314,7 +313,7 @@ subroutine InitAllocate(this, bounds) allocate(this%fpsn_wc_patch (begp:endp)) ; this%fpsn_wc_patch (:) = nan allocate(this%fpsn_wj_patch (begp:endp)) ; this%fpsn_wj_patch (:) = nan allocate(this%fpsn_wp_patch (begp:endp)) ; this%fpsn_wp_patch (:) = nan - allocate(this%iwue_patch (begp:endp)) ; this%iwue_patch (:) = nan + allocate(this%lnca_patch (begp:endp)) ; this%lnca_patch (:) = nan allocate(this%lmrsun_z_patch (begp:endp,1:nlevcan)) ; this%lmrsun_z_patch (:,:) = nan @@ -397,11 +396,6 @@ subroutine InitHistory(this, bounds) avgflag='A', long_name='photosynthesis', & ptr_patch=this%fpsn_patch, set_lake=0._r8, set_urb=0._r8) - this%fpsn_patch(begp:endp) = spval - call hist_addfld1d (fname='IWUE', units='umol CO2 / mol H2O', & - avgflag='A', long_name='intrinsic water use efficiency', & - ptr_patch=this%iwue_patch, set_lake=0._r8, set_urb=0._r8) - ! Don't by default output this rate limiting step as only makes sense if you are outputing ! the others each time-step this%fpsn_wc_patch(begp:endp) = spval @@ -1829,7 +1823,7 @@ end subroutine Photosynthesis !------------------------------------------------------------------------------ subroutine PhotosynthesisTotal (fn, filterp, & - atm2lnd_inst, canopystate_inst, photosyns_inst,solarabs_inst) + atm2lnd_inst, canopystate_inst, photosyns_inst) ! ! Determine total photosynthesis ! @@ -1839,7 +1833,6 @@ subroutine PhotosynthesisTotal (fn, filterp, & type(atm2lnd_type) , intent(in) :: atm2lnd_inst type(canopystate_type) , intent(in) :: canopystate_inst type(photosyns_type) , intent(inout) :: photosyns_inst - type(solarabs_type) , intent(in) :: solarabs_inst ! ! !LOCAL VARIABLES: integer :: f,fp,p,l,g ! indices @@ -1856,10 +1849,6 @@ subroutine PhotosynthesisTotal (fn, filterp, & laisun => canopystate_inst%laisun_patch , & ! Input: [real(r8) (:) ] sunlit leaf area laisha => canopystate_inst%laisha_patch , & ! Input: [real(r8) (:) ] shaded leaf area - par_z => solarabs_inst%parsun_z_patch , & ! Input: [real(r8) (:,:) ] par absorbed per unit lai for canopy layer (w/m**2) - - gs_mol_sun => photosyns_inst%gs_mol_sun_patch , & ! Input: [real(r8) (:,:) ] sunlit leaf stomatal conductance (umol H2O/m**2/s) - gs_mol_sha => photosyns_inst%gs_mol_sha_patch , & ! Input: [real(r8) (:,:) ] shaded leaf stomatal conductance (umol H2O/m**2/s) psnsun => photosyns_inst%psnsun_patch , & ! Input: [real(r8) (:) ] sunlit leaf photosynthesis (umol CO2 /m**2/ s) psnsha => photosyns_inst%psnsha_patch , & ! Input: [real(r8) (:) ] shaded leaf photosynthesis (umol CO2 /m**2/ s) rc13_canair => photosyns_inst%rc13_canair_patch , & ! Output: [real(r8) (:) ] C13O2/C12O2 in canopy air @@ -1880,8 +1869,7 @@ subroutine PhotosynthesisTotal (fn, filterp, & fpsn => photosyns_inst%fpsn_patch , & ! Output: [real(r8) (:) ] photosynthesis (umol CO2 /m**2 /s) fpsn_wc => photosyns_inst%fpsn_wc_patch , & ! Output: [real(r8) (:) ] Rubisco-limited photosynthesis (umol CO2 /m**2 /s) fpsn_wj => photosyns_inst%fpsn_wj_patch , & ! Output: [real(r8) (:) ] RuBP-limited photosynthesis (umol CO2 /m**2 /s) - fpsn_wp => photosyns_inst%fpsn_wp_patch , & ! Output: [real(r8) (:) ] product-limited photosynthesis (umol CO2 /m**2 /s) - iwue => photosyns_inst%iwue_patch & ! Output: [real(r8) (:) ] intrinsic water use efficience (umol CO2 / mol H2o) + fpsn_wp => photosyns_inst%fpsn_wp_patch & ! Output: [real(r8) (:) ] product-limited photosynthesis (umol CO2 /m**2 /s) ) if ( use_c14 ) then @@ -1907,14 +1895,6 @@ subroutine PhotosynthesisTotal (fn, filterp, & fpsn_wc(p) = psnsun_wc(p)*laisun(p) + psnsha_wc(p)*laisha(p) fpsn_wj(p) = psnsun_wj(p)*laisun(p) + psnsha_wj(p)*laisha(p) fpsn_wp(p) = psnsun_wp(p)*laisun(p) + psnsha_wp(p)*laisha(p) - ! calculate inherent wue, 1e6 converts gsmol from umol/m2/s to mol/m2/s - ! iwue calculation not compatible with multi-layer canopy - if (par_z(p,1) > 0._r8) then !daytime - iwue(p) = 1.e06_r8*fpsn(p)/(laisun(p)*gs_mol_sun(p,1)+laisha(p)*gs_mol_sha(p,1)) - else - iwue(p) = spval - write(iulog,*) 'zqz iwue' - end if end if if (use_cn) then diff --git a/src/biogeophys/WaterDiagnosticBulkType.F90 b/src/biogeophys/WaterDiagnosticBulkType.F90 index a5d739b7b4..21cc9d283b 100644 --- a/src/biogeophys/WaterDiagnosticBulkType.F90 +++ b/src/biogeophys/WaterDiagnosticBulkType.F90 @@ -53,7 +53,6 @@ module WaterDiagnosticBulkType real(r8), pointer :: h2osno_top_col (:) ! col top-layer mass of snow [kg] real(r8), pointer :: sno_liq_top_col (:) ! col snow liquid water fraction (mass), top layer [fraction] - real(r8), pointer :: vpd_ref2m_patch (:) ! patch 2 m height surface vapor pressue deficit (Pa) real(r8), pointer :: rh_ref2m_patch (:) ! patch 2 m height surface relative humidity (%) real(r8), pointer :: rh_ref2m_r_patch (:) ! patch 2 m height surface relative humidity - rural (%) real(r8), pointer :: rh_ref2m_u_patch (:) ! patch 2 m height surface relative humidity - urban (%) @@ -191,8 +190,7 @@ subroutine InitBulkAllocate(this, bounds) allocate(this%h2osno_top_col (begc:endc)) ; this%h2osno_top_col (:) = nan allocate(this%sno_liq_top_col (begc:endc)) ; this%sno_liq_top_col (:) = nan - allocate(this%dqgdT_col (begc:endc)) ; this%dqgdT_col (:) = nan - allocate(this%vpd_ref2m_patch (begp:endp)) ; this%vpd_ref2m_patch (:) = nan + allocate(this%dqgdT_col (begc:endc)) ; this%dqgdT_col (:) = nan allocate(this%rh_ref2m_patch (begp:endp)) ; this%rh_ref2m_patch (:) = nan allocate(this%rh_ref2m_u_patch (begp:endp)) ; this%rh_ref2m_u_patch (:) = nan allocate(this%rh_ref2m_r_patch (begp:endp)) ; this%rh_ref2m_r_patch (:) = nan @@ -270,14 +268,6 @@ subroutine InitBulkHistory(this, bounds) long_name=this%info%lname('vertically summed soil cie (veg landunits only)'), & ptr_col=this%h2osoi_ice_tot_col, l2g_scale_type='veg') - this%vpd_ref2m_patch(begp:endp) = spval - call hist_addfld1d ( & - fname=this%info%fname('VPD2M'), & - units='Pa', & - avgflag='A', & - long_name=this%info%lname('2m vapor pressure deficit'), & - ptr_patch=this%vpd_ref2m_patch) - this%rh_ref2m_patch(begp:endp) = spval call hist_addfld1d ( & fname=this%info%fname('RH2M'), & From 6be7f00fc03da8e1a72f00ef5a4d9787ac591b05 Mon Sep 17 00:00:00 2001 From: Daniel Kennedy Date: Tue, 19 Jan 2021 08:39:48 -0700 Subject: [PATCH 05/12] working vpd2m --- src/biogeophys/CanopyFluxesMod.F90 | 5 +++++ src/biogeophys/WaterDiagnosticBulkType.F90 | 12 +++++++++++- 2 files changed, 16 insertions(+), 1 deletion(-) diff --git a/src/biogeophys/CanopyFluxesMod.F90 b/src/biogeophys/CanopyFluxesMod.F90 index b96b43187c..a3116076d4 100644 --- a/src/biogeophys/CanopyFluxesMod.F90 +++ b/src/biogeophys/CanopyFluxesMod.F90 @@ -527,6 +527,7 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, q_ref2m => waterdiagnosticbulk_inst%q_ref2m_patch , & ! Output: [real(r8) (:) ] 2 m height surface specific humidity (kg/kg) rh_ref2m_r => waterdiagnosticbulk_inst%rh_ref2m_r_patch , & ! Output: [real(r8) (:) ] Rural 2 m height surface relative humidity (%) rh_ref2m => waterdiagnosticbulk_inst%rh_ref2m_patch , & ! Output: [real(r8) (:) ] 2 m height surface relative humidity (%) + vpd_ref2m => waterdiagnosticbulk_inst%vpd_ref2m_patch , & ! Output: [real(r8) (:) ] 2 m height surface vapor pressure deficit (Pa) rhaf => waterdiagnosticbulk_inst%rh_af_patch , & ! Output: [real(r8) (:) ] fractional humidity of canopy air [dimensionless] qflx_tran_veg => waterfluxbulk_inst%qflx_tran_veg_patch , & ! Output: [real(r8) (:) ] vegetation transpiration (mm H2O/s) (+ = to atm) @@ -1256,6 +1257,10 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, rh_ref2m(p) = min(100._r8, q_ref2m(p) / qsat_ref2m * 100._r8) rh_ref2m_r(p) = rh_ref2m(p) + ! 2 m height vapor pressure deficit + vpd_ref2m(p) = e_ref2m*(1-rh_ref2m(p)/100._r8) + + ! Human Heat Stress if ( all_human_stress_indices .or. fast_human_stress_indices ) then call KtoC(t_ref2m(p), tc_ref2m(p)) diff --git a/src/biogeophys/WaterDiagnosticBulkType.F90 b/src/biogeophys/WaterDiagnosticBulkType.F90 index 21cc9d283b..030e9773b8 100644 --- a/src/biogeophys/WaterDiagnosticBulkType.F90 +++ b/src/biogeophys/WaterDiagnosticBulkType.F90 @@ -53,6 +53,7 @@ module WaterDiagnosticBulkType real(r8), pointer :: h2osno_top_col (:) ! col top-layer mass of snow [kg] real(r8), pointer :: sno_liq_top_col (:) ! col snow liquid water fraction (mass), top layer [fraction] + real(r8), pointer :: vpd_ref2m_patch (:) ! patch 2 m height surface vapor pressure deficit (Pa) real(r8), pointer :: rh_ref2m_patch (:) ! patch 2 m height surface relative humidity (%) real(r8), pointer :: rh_ref2m_r_patch (:) ! patch 2 m height surface relative humidity - rural (%) real(r8), pointer :: rh_ref2m_u_patch (:) ! patch 2 m height surface relative humidity - urban (%) @@ -190,7 +191,8 @@ subroutine InitBulkAllocate(this, bounds) allocate(this%h2osno_top_col (begc:endc)) ; this%h2osno_top_col (:) = nan allocate(this%sno_liq_top_col (begc:endc)) ; this%sno_liq_top_col (:) = nan - allocate(this%dqgdT_col (begc:endc)) ; this%dqgdT_col (:) = nan + allocate(this%dqgdT_col (begc:endc)) ; this%dqgdT_col (:) = nan + allocate(this%vpd_ref2m_patch (begp:endp)) ; this%vpd_ref2m_patch (:) = nan allocate(this%rh_ref2m_patch (begp:endp)) ; this%rh_ref2m_patch (:) = nan allocate(this%rh_ref2m_u_patch (begp:endp)) ; this%rh_ref2m_u_patch (:) = nan allocate(this%rh_ref2m_r_patch (begp:endp)) ; this%rh_ref2m_r_patch (:) = nan @@ -268,6 +270,14 @@ subroutine InitBulkHistory(this, bounds) long_name=this%info%lname('vertically summed soil cie (veg landunits only)'), & ptr_col=this%h2osoi_ice_tot_col, l2g_scale_type='veg') + this%vpd_ref2m_patch(begp:endp) = spval + call hist_addfld1d ( & + fname=this%info%fname('VPD2M'), & + units='Pa', & + avgflag='A', & + long_name=this%info%lname('2m vapor pressure deficit'), & + ptr_patch=this%vpd_ref2m_patch) + this%rh_ref2m_patch(begp:endp) = spval call hist_addfld1d ( & fname=this%info%fname('RH2M'), & From 137c6a64223a97fc1a944438b5e9aa5f92dc4e34 Mon Sep 17 00:00:00 2001 From: Daniel Kennedy Date: Tue, 19 Jan 2021 10:10:16 -0700 Subject: [PATCH 06/12] revert to move calculation to canopy fluxes --- src/biogeophys/CanopyFluxesMod.F90 | 5 ----- 1 file changed, 5 deletions(-) diff --git a/src/biogeophys/CanopyFluxesMod.F90 b/src/biogeophys/CanopyFluxesMod.F90 index a3116076d4..b96b43187c 100644 --- a/src/biogeophys/CanopyFluxesMod.F90 +++ b/src/biogeophys/CanopyFluxesMod.F90 @@ -527,7 +527,6 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, q_ref2m => waterdiagnosticbulk_inst%q_ref2m_patch , & ! Output: [real(r8) (:) ] 2 m height surface specific humidity (kg/kg) rh_ref2m_r => waterdiagnosticbulk_inst%rh_ref2m_r_patch , & ! Output: [real(r8) (:) ] Rural 2 m height surface relative humidity (%) rh_ref2m => waterdiagnosticbulk_inst%rh_ref2m_patch , & ! Output: [real(r8) (:) ] 2 m height surface relative humidity (%) - vpd_ref2m => waterdiagnosticbulk_inst%vpd_ref2m_patch , & ! Output: [real(r8) (:) ] 2 m height surface vapor pressure deficit (Pa) rhaf => waterdiagnosticbulk_inst%rh_af_patch , & ! Output: [real(r8) (:) ] fractional humidity of canopy air [dimensionless] qflx_tran_veg => waterfluxbulk_inst%qflx_tran_veg_patch , & ! Output: [real(r8) (:) ] vegetation transpiration (mm H2O/s) (+ = to atm) @@ -1257,10 +1256,6 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, rh_ref2m(p) = min(100._r8, q_ref2m(p) / qsat_ref2m * 100._r8) rh_ref2m_r(p) = rh_ref2m(p) - ! 2 m height vapor pressure deficit - vpd_ref2m(p) = e_ref2m*(1-rh_ref2m(p)/100._r8) - - ! Human Heat Stress if ( all_human_stress_indices .or. fast_human_stress_indices ) then call KtoC(t_ref2m(p), tc_ref2m(p)) From b945d93d925402b2239a311f01feb5c23d30c38d Mon Sep 17 00:00:00 2001 From: Daniel Kennedy Date: Tue, 19 Jan 2021 11:01:55 -0700 Subject: [PATCH 07/12] not handling low transpiration well --- src/biogeophys/CanopyFluxesMod.F90 | 28 +++++++++++++++++++--- src/biogeophys/WaterDiagnosticBulkType.F90 | 12 +++++++++- 2 files changed, 36 insertions(+), 4 deletions(-) diff --git a/src/biogeophys/CanopyFluxesMod.F90 b/src/biogeophys/CanopyFluxesMod.F90 index b96b43187c..bdb770a831 100644 --- a/src/biogeophys/CanopyFluxesMod.F90 +++ b/src/biogeophys/CanopyFluxesMod.F90 @@ -16,7 +16,7 @@ module CanopyFluxesMod use clm_varctl , only : iulog, use_cn, use_lch4, use_c13, use_c14, use_cndv, use_fates, & use_luna, use_hydrstress use clm_varpar , only : nlevgrnd, nlevsno - use clm_varcon , only : namep + use clm_varcon , only : namep, spval use pftconMod , only : pftcon use decompMod , only : bounds_type use ActiveLayerMod , only : active_layer_type @@ -276,7 +276,6 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, !added by K.Sakaguchi for stability formulation real(r8), parameter :: ria = 0.5_r8 ! free parameter for stable formulation (currently = 0.5, "gamma" in Sakaguchi&Zeng,2008) - real(r8) :: dtime ! land model time step (sec) real(r8) :: zldis(bounds%begp:bounds%endp) ! reference height "minus" zero displacement height [m] real(r8) :: zeta ! dimensionless height used in Monin-Obukhov theory @@ -469,7 +468,7 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, swmp80_ref2m_r => humanindex_inst%swmp80_ref2m_r_patch , & ! Output: [real(r8) (:) ] Rural 2 m Swamp Cooler temperature 80% effi (C) sabv => solarabs_inst%sabv_patch , & ! Input: [real(r8) (:) ] solar radiation absorbed by vegetation (W/m**2) - + par_z_sun => solarabs_inst%parsun_z_patch , & ! Input: [real(r8) (:,:) ] par absorbed per unit lai for canopy layer (w/m**2) frac_veg_nosno => canopystate_inst%frac_veg_nosno_patch , & ! Input: [integer (:) ] fraction of vegetation not covered by snow (0 OR 1) [-] elai => canopystate_inst%elai_patch , & ! Input: [real(r8) (:) ] one-sided leaf area index with burying by snow esai => canopystate_inst%esai_patch , & ! Input: [real(r8) (:) ] one-sided stem area index with burying by snow @@ -517,6 +516,7 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, qg_h2osfc => waterdiagnosticbulk_inst%qg_h2osfc_col , & ! Input: [real(r8) (:) ] specific humidity at h2osfc surface [kg/kg] qg => waterdiagnosticbulk_inst%qg_col , & ! Input: [real(r8) (:) ] specific humidity at ground surface [kg/kg] dqgdT => waterdiagnosticbulk_inst%dqgdT_col , & ! Input: [real(r8) (:) ] temperature derivative of "qg" + h2osoi_ice => waterstatebulk_inst%h2osoi_ice_col , & ! Input: [real(r8) (:,:) ] ice lens (kg/m2) h2osoi_vol => waterstatebulk_inst%h2osoi_vol_col , & ! Input: [real(r8) (:,:) ] volumetric soil water (0<=h2osoi_vol<=watsat) [m3/m3] by F. Li and S. Levis h2osoi_liq => waterstatebulk_inst%h2osoi_liq_col , & ! Input: [real(r8) (:,:) ] liquid water (kg/m2) @@ -528,6 +528,8 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, rh_ref2m_r => waterdiagnosticbulk_inst%rh_ref2m_r_patch , & ! Output: [real(r8) (:) ] Rural 2 m height surface relative humidity (%) rh_ref2m => waterdiagnosticbulk_inst%rh_ref2m_patch , & ! Output: [real(r8) (:) ] 2 m height surface relative humidity (%) rhaf => waterdiagnosticbulk_inst%rh_af_patch , & ! Output: [real(r8) (:) ] fractional humidity of canopy air [dimensionless] + vpd_ref2m => waterdiagnosticbulk_inst%vpd_ref2m_patch , & ! Output: [real(r8) (:) ] 2 m height surface vapor pressure deficit (Pa) + wue_ei => waterdiagnosticbulk_inst%wue_ei_patch , & ! Output: [real(r8) (:) ] ecosystem-scale inherent water use efficiency (gC kgH2O-1 hPa) qflx_tran_veg => waterfluxbulk_inst%qflx_tran_veg_patch , & ! Output: [real(r8) (:) ] vegetation transpiration (mm H2O/s) (+ = to atm) qflx_evap_veg => waterfluxbulk_inst%qflx_evap_veg_patch , & ! Output: [real(r8) (:) ] vegetation evaporation (mm H2O/s) (+ = to atm) @@ -538,6 +540,7 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, rssun => photosyns_inst%rssun_patch , & ! Output: [real(r8) (:) ] leaf sunlit stomatal resistance (s/m) (output from Photosynthesis) rssha => photosyns_inst%rssha_patch , & ! Output: [real(r8) (:) ] leaf shaded stomatal resistance (s/m) (output from Photosynthesis) + fpsn => photosyns_inst%fpsn_patch , & ! Input: [real(r8) (:) ] photosynthesis (umol CO2 /m**2 /s) grnd_ch4_cond => ch4_inst%grnd_ch4_cond_patch , & ! Output: [real(r8) (:) ] tracer conductance for boundary layer [m/s] @@ -1256,6 +1259,9 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, rh_ref2m(p) = min(100._r8, q_ref2m(p) / qsat_ref2m * 100._r8) rh_ref2m_r(p) = rh_ref2m(p) + ! 2m vapor pressure deficit + vpd_ref2m(p) = e_ref2m*(1._r8-rh_ref2m(p)/100._r8) + ! Human Heat Stress if ( all_human_stress_indices .or. fast_human_stress_indices ) then call KtoC(t_ref2m(p), tc_ref2m(p)) @@ -1353,6 +1359,22 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, call PhotosynthesisTotal(fn, filterp, & atm2lnd_inst, canopystate_inst, photosyns_inst) + ! Calculate water use efficiency + do f = 1, fn + p = filterp(f) + c = patch%column(p) + g = patch%gridcell(p) + + if ( par_z_sun(p,1) > 0._r8 .and. qflx_tran_veg(p)>0._r8) then + ! 12e-6 converts umolCO2->gC + ! 1e-2 converts Pa->hPa + wue_ei(p) = 12e-6_r8*fpsn(p) * 1e-2_r8*vpd_ref2m(p) / qflx_tran_veg(p) + else + wue_ei(p) = spval + end if + + end do + ! Calculate ozone stress. This needs to be done after rssun and rsshade are ! computed by the Photosynthesis routine. However, Photosynthesis also uses the ! ozone stress computed here. Thus, the ozone stress computed in timestep i is diff --git a/src/biogeophys/WaterDiagnosticBulkType.F90 b/src/biogeophys/WaterDiagnosticBulkType.F90 index 030e9773b8..edf013fc1c 100644 --- a/src/biogeophys/WaterDiagnosticBulkType.F90 +++ b/src/biogeophys/WaterDiagnosticBulkType.F90 @@ -53,6 +53,7 @@ module WaterDiagnosticBulkType real(r8), pointer :: h2osno_top_col (:) ! col top-layer mass of snow [kg] real(r8), pointer :: sno_liq_top_col (:) ! col snow liquid water fraction (mass), top layer [fraction] + real(r8), pointer :: wue_ei_patch (:) ! patch ecosystem-scale inherent water use efficiency (gC kgH2O-1 hPa) real(r8), pointer :: vpd_ref2m_patch (:) ! patch 2 m height surface vapor pressure deficit (Pa) real(r8), pointer :: rh_ref2m_patch (:) ! patch 2 m height surface relative humidity (%) real(r8), pointer :: rh_ref2m_r_patch (:) ! patch 2 m height surface relative humidity - rural (%) @@ -192,7 +193,8 @@ subroutine InitBulkAllocate(this, bounds) allocate(this%sno_liq_top_col (begc:endc)) ; this%sno_liq_top_col (:) = nan allocate(this%dqgdT_col (begc:endc)) ; this%dqgdT_col (:) = nan - allocate(this%vpd_ref2m_patch (begp:endp)) ; this%vpd_ref2m_patch (:) = nan + allocate(this%wue_ei_patch (begp:endp)) ; this%wue_ei_patch (:) = nan + allocate(this%vpd_ref2m_patch (begp:endp)) ; this%vpd_ref2m_patch (:) = nan allocate(this%rh_ref2m_patch (begp:endp)) ; this%rh_ref2m_patch (:) = nan allocate(this%rh_ref2m_u_patch (begp:endp)) ; this%rh_ref2m_u_patch (:) = nan allocate(this%rh_ref2m_r_patch (begp:endp)) ; this%rh_ref2m_r_patch (:) = nan @@ -270,6 +272,14 @@ subroutine InitBulkHistory(this, bounds) long_name=this%info%lname('vertically summed soil cie (veg landunits only)'), & ptr_col=this%h2osoi_ice_tot_col, l2g_scale_type='veg') + this%wue_ei_patch(begp:endp) = spval + call hist_addfld1d ( & + fname=this%info%fname('WUE_EI'), & + units='gC kgH2O-1 hPa', & + avgflag='A', & + long_name=this%info%lname('ecosystem-scale inherenet water use efficiency'), & + ptr_patch=this%wue_ei_patch) + this%vpd_ref2m_patch(begp:endp) = spval call hist_addfld1d ( & fname=this%info%fname('VPD2M'), & From 3e69ec4d048ae577842b5ab3b12f969ffc25f548 Mon Sep 17 00:00:00 2001 From: Daniel Kennedy Date: Tue, 19 Jan 2021 13:22:06 -0700 Subject: [PATCH 08/12] require transpiration tolerance --- src/biogeophys/CanopyFluxesMod.F90 | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/biogeophys/CanopyFluxesMod.F90 b/src/biogeophys/CanopyFluxesMod.F90 index bdb770a831..6b4c0399dd 100644 --- a/src/biogeophys/CanopyFluxesMod.F90 +++ b/src/biogeophys/CanopyFluxesMod.F90 @@ -1365,7 +1365,8 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, c = patch%column(p) g = patch%gridcell(p) - if ( par_z_sun(p,1) > 0._r8 .and. qflx_tran_veg(p)>0._r8) then + if ( par_z_sun(p,1) > 0._r8 .and. qflx_tran_veg(p)>4.e-9_r8) then + ! 4e-9 sets a transpiration tolerance of 0.01 W/m2 ! 12e-6 converts umolCO2->gC ! 1e-2 converts Pa->hPa wue_ei(p) = 12e-6_r8*fpsn(p) * 1e-2_r8*vpd_ref2m(p) / qflx_tran_veg(p) From badf7585213feea3b1b33f8f534ae1521e2ec17c Mon Sep 17 00:00:00 2001 From: Daniel Kennedy Date: Wed, 20 Jan 2021 09:26:53 -0700 Subject: [PATCH 09/12] working iwue and wue_ei, but the latter looks pretty unstable --- src/biogeophys/CanopyFluxesMod.F90 | 13 ++++++++++++- src/biogeophys/PhotosynthesisMod.F90 | 4 ++-- src/biogeophys/WaterDiagnosticBulkType.F90 | 10 ++++++++++ 3 files changed, 24 insertions(+), 3 deletions(-) diff --git a/src/biogeophys/CanopyFluxesMod.F90 b/src/biogeophys/CanopyFluxesMod.F90 index 6b4c0399dd..4251592eaa 100644 --- a/src/biogeophys/CanopyFluxesMod.F90 +++ b/src/biogeophys/CanopyFluxesMod.F90 @@ -324,6 +324,7 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, real(r8) :: qsatldT(bounds%begp:bounds%endp) ! derivative of "qsatl" on "t_veg" real(r8) :: e_ref2m ! 2 m height surface saturated vapor pressure [Pa] real(r8) :: qsat_ref2m ! 2 m height surface saturated specific humidity [kg/kg] + real(r8) :: gs ! canopy conductance for iwue cal [molH2O/m2ground/s] real(r8) :: air(bounds%begp:bounds%endp) ! atmos. radiation temporay set real(r8) :: bir(bounds%begp:bounds%endp) ! atmos. radiation temporay set real(r8) :: cir(bounds%begp:bounds%endp) ! atmos. radiation temporay set @@ -530,6 +531,7 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, rhaf => waterdiagnosticbulk_inst%rh_af_patch , & ! Output: [real(r8) (:) ] fractional humidity of canopy air [dimensionless] vpd_ref2m => waterdiagnosticbulk_inst%vpd_ref2m_patch , & ! Output: [real(r8) (:) ] 2 m height surface vapor pressure deficit (Pa) wue_ei => waterdiagnosticbulk_inst%wue_ei_patch , & ! Output: [real(r8) (:) ] ecosystem-scale inherent water use efficiency (gC kgH2O-1 hPa) + iwue => waterdiagnosticbulk_inst%iwue_patch , & ! Output: [real(r8) (:) ] ecosystem-scale inherent water use efficiency (gC kgH2O-1 hPa) qflx_tran_veg => waterfluxbulk_inst%qflx_tran_veg_patch , & ! Output: [real(r8) (:) ] vegetation transpiration (mm H2O/s) (+ = to atm) qflx_evap_veg => waterfluxbulk_inst%qflx_evap_veg_patch , & ! Output: [real(r8) (:) ] vegetation evaporation (mm H2O/s) (+ = to atm) @@ -537,7 +539,8 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, qflx_ev_snow => waterfluxbulk_inst%qflx_ev_snow_patch , & ! Output: [real(r8) (:) ] evaporation flux from snow (mm H2O/s) [+ to atm] qflx_ev_soil => waterfluxbulk_inst%qflx_ev_soil_patch , & ! Output: [real(r8) (:) ] evaporation flux from soil (mm H2O/s) [+ to atm] qflx_ev_h2osfc => waterfluxbulk_inst%qflx_ev_h2osfc_patch , & ! Output: [real(r8) (:) ] evaporation flux from h2osfc (mm H2O/s) [+ to atm] - + gs_mol_sun => photosyns_inst%gs_mol_sun_patch , & ! Input: [real(r8) (:) ] patch sunlit leaf stomatal conductance (umol H2O/m**2/s) + gs_mol_sha => photosyns_inst%gs_mol_sha_patch , & ! Input: [real(r8) (:) ] patch shaded leaf stomatal conductance (umol H2O/m**2/s) rssun => photosyns_inst%rssun_patch , & ! Output: [real(r8) (:) ] leaf sunlit stomatal resistance (s/m) (output from Photosynthesis) rssha => photosyns_inst%rssha_patch , & ! Output: [real(r8) (:) ] leaf shaded stomatal resistance (s/m) (output from Photosynthesis) fpsn => photosyns_inst%fpsn_patch , & ! Input: [real(r8) (:) ] photosynthesis (umol CO2 /m**2 /s) @@ -1369,9 +1372,17 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, ! 4e-9 sets a transpiration tolerance of 0.01 W/m2 ! 12e-6 converts umolCO2->gC ! 1e-2 converts Pa->hPa + ! 1e-6 converts umolH2O->molH2O wue_ei(p) = 12e-6_r8*fpsn(p) * 1e-2_r8*vpd_ref2m(p) / qflx_tran_veg(p) + gs = 1.e-6_r8*(laisun(p)*gs_mol_sun(p,1)+laisha(p)*gs_mol_sha(p,1)) + if ( gs>0._r8 ) then + iwue(p) = fpsn(p)/gs + else + iwue(p) = spval + end if else wue_ei(p) = spval + iwue(p) = spval end if end do diff --git a/src/biogeophys/PhotosynthesisMod.F90 b/src/biogeophys/PhotosynthesisMod.F90 index 44e28a681c..50b11883df 100644 --- a/src/biogeophys/PhotosynthesisMod.F90 +++ b/src/biogeophys/PhotosynthesisMod.F90 @@ -127,8 +127,8 @@ module PhotosynthesisMod real(r8), pointer, private :: vcmax_z_phs_patch (:,:,:) ! patch maximum rate of carboxylation (umol co2/m**2/s) real(r8), pointer, private :: kp_z_phs_patch (:,:,:) ! patch initial slope of CO2 response curve (C4 plants) real(r8), pointer, private :: tpu_z_phs_patch (:,:,:) ! patch triose phosphate utilization rate (umol CO2/m**2/s) - real(r8), pointer, private :: gs_mol_sun_patch (:,:) ! patch sunlit leaf stomatal conductance (umol H2O/m**2/s) - real(r8), pointer, private :: gs_mol_sha_patch (:,:) ! patch shaded leaf stomatal conductance (umol H2O/m**2/s) + real(r8), pointer, public :: gs_mol_sun_patch (:,:) ! patch sunlit leaf stomatal conductance (umol H2O/m**2/s) + real(r8), pointer, public :: gs_mol_sha_patch (:,:) ! patch shaded leaf stomatal conductance (umol H2O/m**2/s) real(r8), pointer, private :: gs_mol_sun_ln_patch (:,:) ! patch sunlit leaf stomatal conductance averaged over 1 hour before to 1 hour after local noon (umol H2O/m**2/s) real(r8), pointer, private :: gs_mol_sha_ln_patch (:,:) ! patch shaded leaf stomatal conductance averaged over 1 hour before to 1 hour after local noon (umol H2O/m**2/s) real(r8), pointer, private :: ac_patch (:,:) ! patch Rubisco-limited gross photosynthesis (umol CO2/m**2/s) diff --git a/src/biogeophys/WaterDiagnosticBulkType.F90 b/src/biogeophys/WaterDiagnosticBulkType.F90 index edf013fc1c..116d2577b6 100644 --- a/src/biogeophys/WaterDiagnosticBulkType.F90 +++ b/src/biogeophys/WaterDiagnosticBulkType.F90 @@ -53,6 +53,7 @@ module WaterDiagnosticBulkType real(r8), pointer :: h2osno_top_col (:) ! col top-layer mass of snow [kg] real(r8), pointer :: sno_liq_top_col (:) ! col snow liquid water fraction (mass), top layer [fraction] + real(r8), pointer :: iwue_patch (:) ! patch intrinsic water use efficiency (umolCO2/molH2O) real(r8), pointer :: wue_ei_patch (:) ! patch ecosystem-scale inherent water use efficiency (gC kgH2O-1 hPa) real(r8), pointer :: vpd_ref2m_patch (:) ! patch 2 m height surface vapor pressure deficit (Pa) real(r8), pointer :: rh_ref2m_patch (:) ! patch 2 m height surface relative humidity (%) @@ -193,6 +194,7 @@ subroutine InitBulkAllocate(this, bounds) allocate(this%sno_liq_top_col (begc:endc)) ; this%sno_liq_top_col (:) = nan allocate(this%dqgdT_col (begc:endc)) ; this%dqgdT_col (:) = nan + allocate(this%iwue_patch (begp:endp)) ; this%iwue_patch (:) = nan allocate(this%wue_ei_patch (begp:endp)) ; this%wue_ei_patch (:) = nan allocate(this%vpd_ref2m_patch (begp:endp)) ; this%vpd_ref2m_patch (:) = nan allocate(this%rh_ref2m_patch (begp:endp)) ; this%rh_ref2m_patch (:) = nan @@ -272,6 +274,14 @@ subroutine InitBulkHistory(this, bounds) long_name=this%info%lname('vertically summed soil cie (veg landunits only)'), & ptr_col=this%h2osoi_ice_tot_col, l2g_scale_type='veg') + this%iwue_patch(begp:endp) = spval + call hist_addfld1d ( & + fname=this%info%fname('IWUE'), & + units='umolCO2/molH2O', & + avgflag='A', & + long_name=this%info%lname('intrinsic water use efficiency'), & + ptr_patch=this%iwue_patch) + this%wue_ei_patch(begp:endp) = spval call hist_addfld1d ( & fname=this%info%fname('WUE_EI'), & From b3245d4cd1b4dc917a336eb370799ca76fd83f43 Mon Sep 17 00:00:00 2001 From: Daniel Kennedy Date: Mon, 25 Jan 2021 17:26:15 -0700 Subject: [PATCH 10/12] working iWUE at local noon --- src/biogeophys/CanopyFluxesMod.F90 | 49 ++++++++++------------ src/biogeophys/WaterDiagnosticBulkType.F90 | 26 ++++-------- 2 files changed, 30 insertions(+), 45 deletions(-) diff --git a/src/biogeophys/CanopyFluxesMod.F90 b/src/biogeophys/CanopyFluxesMod.F90 index 4251592eaa..759e78951a 100644 --- a/src/biogeophys/CanopyFluxesMod.F90 +++ b/src/biogeophys/CanopyFluxesMod.F90 @@ -15,7 +15,7 @@ module CanopyFluxesMod use abortutils , only : endrun use clm_varctl , only : iulog, use_cn, use_lch4, use_c13, use_c14, use_cndv, use_fates, & use_luna, use_hydrstress - use clm_varpar , only : nlevgrnd, nlevsno + use clm_varpar , only : nlevgrnd, nlevsno, nlevcan use clm_varcon , only : namep, spval use pftconMod , only : pftcon use decompMod , only : bounds_type @@ -221,7 +221,7 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, ! ! !USES: use shr_const_mod , only : SHR_CONST_RGAS - use clm_time_manager , only : get_step_size_real, get_prev_date,is_end_curr_day + use clm_time_manager , only : get_step_size_real, get_prev_date,is_end_curr_day, is_near_local_noon use clm_varcon , only : sb, cpair, hvap, vkc, grav, denice use clm_varcon , only : denh2o, tfrz, tlsai_crit, alpha_aero use clm_varcon , only : c14ratio @@ -400,7 +400,7 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, integer :: ft ! plant functional type index real(r8) :: h2ocan ! total canopy water (mm H2O) real(r8) :: dt_veg_temp(bounds%begp:bounds%endp) - integer :: iv + integer, parameter :: iv=1 ! index for first canopy layer (iwue calculation) logical :: is_end_day ! is end of current day real(r8) :: tl_ini ! leaf temperature from beginning of time step [K] real(r8) :: ts_ini ! stem temperature from beginning of time step [K] @@ -470,6 +470,7 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, sabv => solarabs_inst%sabv_patch , & ! Input: [real(r8) (:) ] solar radiation absorbed by vegetation (W/m**2) par_z_sun => solarabs_inst%parsun_z_patch , & ! Input: [real(r8) (:,:) ] par absorbed per unit lai for canopy layer (w/m**2) + frac_veg_nosno => canopystate_inst%frac_veg_nosno_patch , & ! Input: [integer (:) ] fraction of vegetation not covered by snow (0 OR 1) [-] elai => canopystate_inst%elai_patch , & ! Input: [real(r8) (:) ] one-sided leaf area index with burying by snow esai => canopystate_inst%esai_patch , & ! Input: [real(r8) (:) ] one-sided stem area index with burying by snow @@ -530,8 +531,7 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, rh_ref2m => waterdiagnosticbulk_inst%rh_ref2m_patch , & ! Output: [real(r8) (:) ] 2 m height surface relative humidity (%) rhaf => waterdiagnosticbulk_inst%rh_af_patch , & ! Output: [real(r8) (:) ] fractional humidity of canopy air [dimensionless] vpd_ref2m => waterdiagnosticbulk_inst%vpd_ref2m_patch , & ! Output: [real(r8) (:) ] 2 m height surface vapor pressure deficit (Pa) - wue_ei => waterdiagnosticbulk_inst%wue_ei_patch , & ! Output: [real(r8) (:) ] ecosystem-scale inherent water use efficiency (gC kgH2O-1 hPa) - iwue => waterdiagnosticbulk_inst%iwue_patch , & ! Output: [real(r8) (:) ] ecosystem-scale inherent water use efficiency (gC kgH2O-1 hPa) + iwue_ln => waterdiagnosticbulk_inst%iwue_ln_patch , & ! Output: [real(r8) (:) ] local noon ecosystem-scale inherent water use efficiency (gC kgH2O-1 hPa) qflx_tran_veg => waterfluxbulk_inst%qflx_tran_veg_patch , & ! Output: [real(r8) (:) ] vegetation transpiration (mm H2O/s) (+ = to atm) qflx_evap_veg => waterfluxbulk_inst%qflx_evap_veg_patch , & ! Output: [real(r8) (:) ] vegetation evaporation (mm H2O/s) (+ = to atm) @@ -1363,30 +1363,25 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, atm2lnd_inst, canopystate_inst, photosyns_inst) ! Calculate water use efficiency - do f = 1, fn - p = filterp(f) - c = patch%column(p) - g = patch%gridcell(p) - - if ( par_z_sun(p,1) > 0._r8 .and. qflx_tran_veg(p)>4.e-9_r8) then - ! 4e-9 sets a transpiration tolerance of 0.01 W/m2 - ! 12e-6 converts umolCO2->gC - ! 1e-2 converts Pa->hPa - ! 1e-6 converts umolH2O->molH2O - wue_ei(p) = 12e-6_r8*fpsn(p) * 1e-2_r8*vpd_ref2m(p) / qflx_tran_veg(p) - gs = 1.e-6_r8*(laisun(p)*gs_mol_sun(p,1)+laisha(p)*gs_mol_sha(p,1)) - if ( gs>0._r8 ) then - iwue(p) = fpsn(p)/gs + ! does not support multi-layer canopy + if (nlevcan == 1) then + do f = 1, fn + p = filterp(f) + c = patch%column(p) + g = patch%gridcell(p) + + if ( is_near_local_noon( grc%londeg(g), deltasec=3600 ) .and. fpsn(p)>0._r8 )then + gs = 1.e-6_r8*(laisun(p)*gs_mol_sun(p,iv)+laisha(p)*gs_mol_sha(p,iv)) ! 1e-6 converts umolH2O->molH2O + if ( gs>0._r8 ) then + iwue_ln(p) = fpsn(p)/gs + else + iwue_ln(p) = spval + end if else - iwue(p) = spval + iwue_ln(p) = spval end if - else - wue_ei(p) = spval - iwue(p) = spval - end if - - end do - + end do + end if ! Calculate ozone stress. This needs to be done after rssun and rsshade are ! computed by the Photosynthesis routine. However, Photosynthesis also uses the ! ozone stress computed here. Thus, the ozone stress computed in timestep i is diff --git a/src/biogeophys/WaterDiagnosticBulkType.F90 b/src/biogeophys/WaterDiagnosticBulkType.F90 index 116d2577b6..7bbbd94fdf 100644 --- a/src/biogeophys/WaterDiagnosticBulkType.F90 +++ b/src/biogeophys/WaterDiagnosticBulkType.F90 @@ -17,7 +17,7 @@ module WaterDiagnosticBulkType use decompMod , only : bounds_type use abortutils , only : endrun use clm_varctl , only : use_cn, iulog, use_luna - use clm_varpar , only : nlevgrnd, nlevsno + use clm_varpar , only : nlevgrnd, nlevsno, nlevcan use clm_varcon , only : spval use LandunitType , only : lun use ColumnType , only : col @@ -53,8 +53,7 @@ module WaterDiagnosticBulkType real(r8), pointer :: h2osno_top_col (:) ! col top-layer mass of snow [kg] real(r8), pointer :: sno_liq_top_col (:) ! col snow liquid water fraction (mass), top layer [fraction] - real(r8), pointer :: iwue_patch (:) ! patch intrinsic water use efficiency (umolCO2/molH2O) - real(r8), pointer :: wue_ei_patch (:) ! patch ecosystem-scale inherent water use efficiency (gC kgH2O-1 hPa) + real(r8), pointer :: iwue_ln_patch (:) ! patch intrinsic water use efficiency near local noon (umolCO2/molH2O) real(r8), pointer :: vpd_ref2m_patch (:) ! patch 2 m height surface vapor pressure deficit (Pa) real(r8), pointer :: rh_ref2m_patch (:) ! patch 2 m height surface relative humidity (%) real(r8), pointer :: rh_ref2m_r_patch (:) ! patch 2 m height surface relative humidity - rural (%) @@ -194,8 +193,7 @@ subroutine InitBulkAllocate(this, bounds) allocate(this%sno_liq_top_col (begc:endc)) ; this%sno_liq_top_col (:) = nan allocate(this%dqgdT_col (begc:endc)) ; this%dqgdT_col (:) = nan - allocate(this%iwue_patch (begp:endp)) ; this%iwue_patch (:) = nan - allocate(this%wue_ei_patch (begp:endp)) ; this%wue_ei_patch (:) = nan + allocate(this%iwue_ln_patch (begp:endp)) ; this%iwue_ln_patch (:) = nan allocate(this%vpd_ref2m_patch (begp:endp)) ; this%vpd_ref2m_patch (:) = nan allocate(this%rh_ref2m_patch (begp:endp)) ; this%rh_ref2m_patch (:) = nan allocate(this%rh_ref2m_u_patch (begp:endp)) ; this%rh_ref2m_u_patch (:) = nan @@ -274,21 +272,13 @@ subroutine InitBulkHistory(this, bounds) long_name=this%info%lname('vertically summed soil cie (veg landunits only)'), & ptr_col=this%h2osoi_ice_tot_col, l2g_scale_type='veg') - this%iwue_patch(begp:endp) = spval + this%iwue_ln_patch(begp:endp) = spval call hist_addfld1d ( & - fname=this%info%fname('IWUE'), & + fname=this%info%fname('IWUELN'), & units='umolCO2/molH2O', & - avgflag='A', & - long_name=this%info%lname('intrinsic water use efficiency'), & - ptr_patch=this%iwue_patch) - - this%wue_ei_patch(begp:endp) = spval - call hist_addfld1d ( & - fname=this%info%fname('WUE_EI'), & - units='gC kgH2O-1 hPa', & - avgflag='A', & - long_name=this%info%lname('ecosystem-scale inherenet water use efficiency'), & - ptr_patch=this%wue_ei_patch) + avgflag='A', & + long_name=this%info%lname('local noon intrinsic water use efficiency'), & + ptr_patch=this%iwue_ln_patch, set_lake=spval, set_urb=spval) this%vpd_ref2m_patch(begp:endp) = spval call hist_addfld1d ( & From 2563cc803dec3e9cd9e75856dd29ba8cf527528e Mon Sep 17 00:00:00 2001 From: Daniel Kennedy Date: Fri, 29 Jan 2021 16:50:31 -0700 Subject: [PATCH 11/12] endrun if nlevcan>1 --- src/biogeophys/CanopyFluxesMod.F90 | 3 +++ 1 file changed, 3 insertions(+) diff --git a/src/biogeophys/CanopyFluxesMod.F90 b/src/biogeophys/CanopyFluxesMod.F90 index 759e78951a..a1b60779e0 100644 --- a/src/biogeophys/CanopyFluxesMod.F90 +++ b/src/biogeophys/CanopyFluxesMod.F90 @@ -1381,6 +1381,9 @@ subroutine CanopyFluxes(bounds, num_exposedvegp, filter_exposedvegp, iwue_ln(p) = spval end if end do + else + call endrun(msg=' ERROR: IWUELN calculation not compatible with nlevcan>1 ' // & + errMsg(sourcefile, __LINE__)) end if ! Calculate ozone stress. This needs to be done after rssun and rsshade are ! computed by the Photosynthesis routine. However, Photosynthesis also uses the From 4b2d2572c6b1ab5fbe982d7f926ee28d576bc462 Mon Sep 17 00:00:00 2001 From: Erik Kluzek Date: Mon, 1 Feb 2021 16:20:04 -0700 Subject: [PATCH 12/12] Add new diagnostics to the restart file, so restarts can be exact --- src/biogeophys/WaterDiagnosticBulkType.F90 | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) diff --git a/src/biogeophys/WaterDiagnosticBulkType.F90 b/src/biogeophys/WaterDiagnosticBulkType.F90 index 7bbbd94fdf..b235bd3b43 100644 --- a/src/biogeophys/WaterDiagnosticBulkType.F90 +++ b/src/biogeophys/WaterDiagnosticBulkType.F90 @@ -702,6 +702,22 @@ subroutine RestartBulk(this, bounds, ncid, flag, writing_finidat_interp_dest_fil writing_finidat_interp_dest_file = writing_finidat_interp_dest_file, & waterstatebulk_inst = waterstatebulk_inst) + call restartvar(ncid=ncid, flag=flag, & + varname=this%info%fname('IWUELN'), & + xtype=ncd_double, & + dim1name='pft', & + long_name=this%info%lname('local noon intrinsic water use efficiency'), & + units='umolCO2/molH2O', & + interpinic_flag='interp', readvar=readvar, data=this%iwue_ln_patch) + + call restartvar(ncid=ncid, flag=flag, & + varname=this%info%fname('VPD2M'), & + xtype=ncd_double, & + dim1name='pft', & + long_name=this%info%lname('2m vapor pressure deficit'), & + units='Pa', & + interpinic_flag='interp', readvar=readvar, data=this%vpd_ref2m_patch) + call restartvar(ncid=ncid, flag=flag, & varname=this%info%fname('FWET'), & xtype=ncd_double, &