From 230e3a8159759a9f47bf5f7324006d4d74a365ad Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Thu, 16 Jun 2016 14:58:52 -0700 Subject: [PATCH 1/9] pruned unnecessary variables from EDBtran, removed some confusing redundant calls to zeroing and initializing energyflux_inst%btran_patch in EDBtran --- biogeophys/EDBtranMod.F90 | 7 ------- 1 file changed, 7 deletions(-) diff --git a/biogeophys/EDBtranMod.F90 b/biogeophys/EDBtranMod.F90 index 9c5475059c..030eb4f2b2 100644 --- a/biogeophys/EDBtranMod.F90 +++ b/biogeophys/EDBtranMod.F90 @@ -121,10 +121,7 @@ subroutine btran_ed( bounds, p, sites, nsites, hsites, & sucsat => soilstate_inst%sucsat_col , & ! Input: [real(r8) (:,:) ] minimum soil suction (mm) watsat => soilstate_inst%watsat_col , & ! Input: [real(r8) (:,:) ] volumetric soil water at saturation (porosity) - watdry => soilstate_inst%watdry_col , & ! Input: [real(r8) (:,:) ] btran parameter for btran=0 - watopt => soilstate_inst%watopt_col , & ! Input: [real(r8) (:,:) ] btran parameter for btran = 1 bsw => soilstate_inst%bsw_col , & ! Input: [real(r8) (:,:) ] Clapp and Hornberger "b" - soilbeta => soilstate_inst%soilbeta_col , & ! Input: [real(r8) (:) ] soil wetness relative to field capacity sand => soilstate_inst%sandfrac_patch , & ! Input: [real(r8) (:) ] % sand of soil rootr => soilstate_inst%rootr_patch , & ! Output: [real(r8) (:,:) ] Fraction of water uptake in each layer @@ -135,7 +132,6 @@ subroutine btran_ed( bounds, p, sites, nsites, hsites, & t_soisno => temperature_inst%t_soisno_col , & ! Input: [real(r8) (:,:) ] soil temperature (Kelvin) btran => energyflux_inst%btran_patch , & ! Output: [real(r8) (:) ] transpiration wetness factor (0 to 1) - btran2 => energyflux_inst%btran2_patch , & ! Output: [real(r8) (:) ] rresis => energyflux_inst%rresis_patch & ! Output: [real(r8) (:,:) ] root resistance by layer (0-1) (nlevgrnd) ) @@ -172,8 +168,6 @@ subroutine btran_ed( bounds, p, sites, nsites, hsites, & end if end do !j - btran(p) = currentPatch%btran_ft(1) !FIX(RF,032414) for TRF where is this used? - ! Normalize root resistances to get layer contribution to ET do j = 1,nlevgrnd if (currentPatch%btran_ft(FT) > 0.0_r8) then @@ -196,7 +190,6 @@ subroutine btran_ed( bounds, p, sites, nsites, hsites, & do j = 1,nlevgrnd rootr(p,j) = 0._r8 - btran(p) = 0.0_r8 do FT = 1,numpft_ed if(sum(pftgs) > 0._r8)then !prevent problem with the first timestep - might fail !bit-retart test as a result? FIX(RF,032414) From 86e4bd38fb612973e556592b501157a75384152c Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Mon, 20 Jun 2016 17:32:53 -0700 Subject: [PATCH 2/9] First pass of adding the bc_in and bc_out memory structures. Prototyped on sun/shade fraction call. Compiled, untested. --- biogeochem/EDPhysiologyMod.F90 | 17 +-- biogeophys/EDAccumulateFluxesMod.F90 | 5 + biogeophys/EDPhotosynthesisMod.F90 | 2 +- biogeophys/EDSurfaceAlbedoMod.F90 | 199 +++++++++++++++------------ main/EDMainMod.F90 | 9 +- main/FatesInterfaceMod.F90 | 143 +++++++++++++++++-- 6 files changed, 256 insertions(+), 119 deletions(-) diff --git a/biogeochem/EDPhysiologyMod.F90 b/biogeochem/EDPhysiologyMod.F90 index c5aff001f1..4f00de3b8b 100755 --- a/biogeochem/EDPhysiologyMod.F90 +++ b/biogeochem/EDPhysiologyMod.F90 @@ -70,7 +70,7 @@ subroutine canopy_derivs( currentPatch ) end subroutine canopy_derivs ! ============================================================================ - subroutine non_canopy_derivs( currentPatch, temperature_inst, soilstate_inst, waterstate_inst) + subroutine non_canopy_derivs( currentPatch, temperature_inst ) ! ! !DESCRIPTION: ! Returns time differentials of the state vector @@ -80,8 +80,6 @@ subroutine non_canopy_derivs( currentPatch, temperature_inst, soilstate_inst, wa ! !ARGUMENTS type(ed_patch_type) , intent(inout) :: currentPatch type(temperature_type) , intent(in) :: temperature_inst - type(soilstate_type) , intent(in) :: soilstate_inst - type(waterstate_type) , intent(in) :: waterstate_inst ! ! !LOCAL VARIABLES: integer c,p @@ -108,7 +106,7 @@ subroutine non_canopy_derivs( currentPatch, temperature_inst, soilstate_inst, wa ! update fragmenting pool fluxes call cwd_input(currentPatch) - call cwd_out( currentPatch, temperature_inst, soilstate_inst, waterstate_inst) + call cwd_out( currentPatch, temperature_inst) do p = 1,numpft_ed currentPatch%dseed_dt(p) = currentPatch%seeds_in(p) - currentPatch%seed_decay(p) - currentPatch%seed_germination(p) @@ -244,7 +242,7 @@ subroutine phenology( currentSite, temperature_inst, waterstate_inst) ! ! !USES: use clm_varcon, only : tfrz - use clm_time_manager, only : get_days_per_year, get_curr_date + use clm_time_manager, only : get_curr_date use clm_time_manager, only : get_ref_date, timemgr_datediff use EDTypesMod, only : udata use PatchType , only : patch @@ -1145,7 +1143,7 @@ subroutine fragmentation_scaler( currentPatch, temperature_inst ) ! !USES: use shr_const_mod , only : SHR_CONST_PI, SHR_CONST_TKFRZ use EDSharedParamsMod , only : EDParamsShareInst - use PatchType , only : patch + ! ! !ARGUMENTS type(ed_patch_type) , intent(inout) :: currentPatch @@ -1154,7 +1152,7 @@ subroutine fragmentation_scaler( currentPatch, temperature_inst ) ! !LOCAL VARIABLES: logical :: use_century_tfunc = .false. type(ed_site_type), pointer :: currentSite - integer :: c,p,j + integer :: p,j real(r8) :: t_scalar real(r8) :: w_scalar real(r8) :: catanf ! hyperbolic temperature function from CENTURY @@ -1173,7 +1171,6 @@ subroutine fragmentation_scaler( currentPatch, temperature_inst ) ! c = currentPatch%siteptr%clmcolumn p = currentPatch%clm_pno - c = patch%column(p) ! set "froz_q10" parameter froz_q10 = EDParamsShareInst%froz_q10 @@ -1204,7 +1201,7 @@ subroutine fragmentation_scaler( currentPatch, temperature_inst ) end subroutine fragmentation_scaler ! ============================================================================ - subroutine cwd_out( currentPatch, temperature_inst, soilstate_inst, waterstate_inst) + subroutine cwd_out( currentPatch, temperature_inst ) ! ! !DESCRIPTION: ! Simple CWD fragmentation Model @@ -1217,8 +1214,6 @@ subroutine cwd_out( currentPatch, temperature_inst, soilstate_inst, waterstate_i ! !ARGUMENTS type(ed_patch_type) , intent(inout), target :: currentPatch type(temperature_type) , intent(in) :: temperature_inst - type(soilstate_type) , intent(in) :: soilstate_inst - type(waterstate_type) , intent(in) :: waterstate_inst ! ! !LOCAL VARIABLES: type(ed_site_type), pointer :: currentSite diff --git a/biogeophys/EDAccumulateFluxesMod.F90 b/biogeophys/EDAccumulateFluxesMod.F90 index ad32348ca0..006a03537e 100644 --- a/biogeophys/EDAccumulateFluxesMod.F90 +++ b/biogeophys/EDAccumulateFluxesMod.F90 @@ -55,6 +55,11 @@ subroutine AccumulateFluxes_ED(bounds, p, sites, nsites, hsites , photosyns_inst psncanopy => photosyns_inst%psncanopy_patch & ! Output: [real(r8) (:,:)] canopy scale photosynthesis umol CO2 /m**2/ s ) + + ! INTERF-TODO: WHY IS THIS BEING UPDATED? + ! IT IS JUST GOING TO BE ZEROED A THE END OF THE FUNCTION + ! THAT CALLS THIS SUBROUTINE (CANOPYFLUXES), AND IT WON'T + ! BE USED BETWEEN NOW AND THEN fpsn(p) = psncanopy(p) if (patch%is_veg(p)) then diff --git a/biogeophys/EDPhotosynthesisMod.F90 b/biogeophys/EDPhotosynthesisMod.F90 index 058cbc749c..f4f9b6442c 100644 --- a/biogeophys/EDPhotosynthesisMod.F90 +++ b/biogeophys/EDPhotosynthesisMod.F90 @@ -59,7 +59,7 @@ subroutine Photosynthesis_ED (bounds, fn, filterp, esat_tv, eair, oair, cair, & real(r8) , intent(in) :: eair( bounds%begp: ) ! vapor pressure of canopy air (Pa) real(r8) , intent(in) :: oair( bounds%begp: ) ! Atmospheric O2 partial pressure (Pa) real(r8) , intent(in) :: cair( bounds%begp: ) ! Atmospheric CO2 partial pressure (Pa) - real(r8) , intent(inout) :: rb( bounds%begp: ) ! boundary layer resistance (s/m) + real(r8) , intent(in) :: rb( bounds%begp: ) ! boundary layer resistance (s/m) real(r8) , intent(in) :: dayl_factor( bounds%begp: ) ! scalar (0-1) for daylength type(ed_site_type) , intent(inout), target :: sites(nsites) integer , intent(in) :: nsites diff --git a/biogeophys/EDSurfaceAlbedoMod.F90 b/biogeophys/EDSurfaceAlbedoMod.F90 index 3110bd63a4..d9b84bd503 100644 --- a/biogeophys/EDSurfaceAlbedoMod.F90 +++ b/biogeophys/EDSurfaceAlbedoMod.F90 @@ -17,7 +17,9 @@ module EDSurfaceRadiationMod use shr_log_mod , only : errMsg => shr_log_errMsg use decompMod , only : bounds_type use clm_varpar , only : numrad, nclmax - + use FatesInterfaceMod , only : bc_in_type, & + bc_out_type + implicit none private @@ -28,7 +30,11 @@ module EDSurfaceRadiationMod real(r8), public :: albice(numrad) = & ! albedo land ice by waveband (1=vis, 2=nir) (/ 0.80_r8, 0.55_r8 /) - + + ! INTERF-TODO: THIS NEEDS SOME CONSISTENCY AND SHOULD BE SET IN THE INTERFACE + ! WITH OTHER DIMENSIONS + integer, parameter :: ipar = 1 ! The band index for PAR + contains subroutine ED_Norman_Radiation (bounds, & @@ -950,85 +956,93 @@ subroutine ED_Norman_Radiation (bounds, & end associate end subroutine ED_Norman_Radiation - -subroutine ED_SunShadeFracs(cpatch,forc_par_d,forc_par_i,fsun) - - - use clm_varctl , only : iulog - - ! Arguments In - - real(r8),intent(in) :: forc_par_d - real(r8),intent(in) :: forc_par_i - - ! Arguments inout - type (ed_patch_type),intent(inout), target :: cpatch ! c"urrent" patch - - - ! Arguments Out - real(r8),intent(out) :: fsun - - ! locals - real(r8) :: sunlai - real(r8) :: shalai - integer :: CL - integer :: FT - integer :: iv - - - ! zero out various datas - cpatch%ed_parsun_z(:,:,:) = 0._r8 - cpatch%ed_parsha_z(:,:,:) = 0._r8 - cpatch%ed_laisun_z(:,:,:) = 0._r8 - cpatch%ed_laisha_z(:,:,:) = 0._r8 - - fsun = 0._r8 - sunlai = 0._r8 - shalai = 0._r8 - - ! Loop over patches to calculate laisun_z and laisha_z for each layer. - ! Derive canopy laisun, laisha, and fsun from layer sums. - ! If sun/shade big leaf code, nrad=1 and fsun_z(p,1) and tlai_z(p,1) from - ! SurfaceAlbedo is canopy integrated so that layer value equals canopy value. - - ! cpatch%f_sun is calculated in the surface_albedo routine... - - do CL = 1, cpatch%NCL_p - do FT = 1,numpft_ed - do iv = 1, cpatch%nrad(CL,ft) !NORMAL CASE. - - ! FIX(SPM,040114) - existing comment - ! ** Should this be elai or tlai? Surely we only do radiation for elai? - - cpatch%ed_laisun_z(CL,ft,iv) = cpatch%elai_profile(CL,ft,iv) * & - cpatch%f_sun(CL,ft,iv) - - if ( DEBUG ) write(iulog,*) 'edsurfRad 570 ',cpatch%elai_profile(CL,ft,iv) - if ( DEBUG ) write(iulog,*) 'edsurfRad 571 ',cpatch%f_sun(CL,ft,iv) - - cpatch%ed_laisha_z(CL,ft,iv) = cpatch%elai_profile(CL,ft,iv) * & - (1._r8 - cpatch%f_sun(CL,ft,iv)) - - end do - - !needed for the VOC emissions, etc. - sunlai = sunlai + sum(cpatch%ed_laisun_z(CL,ft,1:cpatch%nrad(CL,ft))) - shalai = shalai + sum(cpatch%ed_laisha_z(CL,ft,1:cpatch%nrad(CL,ft))) - - end do - end do - - if(sunlai+shalai > 0._r8)then - fsun = sunlai / (sunlai+shalai) - else - fsun = 0._r8 - endif - - if(fsun > 1._r8)then - write(iulog,*) 'too much leaf area in profile', fsun, & - cpatch%lai,sunlai,shalai - endif - + ! ====================================================================================== + + subroutine ED_SunShadeFracs(sites,nsites,bc_in,bc_out) + + use clm_varctl , only : iulog + + ! Arguments + + type(ed_site_type),intent(inout),target :: sites(nsites) + integer,intent(in) :: nsites + type(bc_in_type),intent(in) :: bc_in(nsites) + type(bc_out_type),intent(out) :: bc_out(nsites) + + + ! locals + + type (ed_patch_type),pointer :: cpatch ! c"urrent" patch + real(r8) :: sunlai + real(r8) :: shalai + integer :: CL + integer :: FT + integer :: iv + integer :: s + integer :: ifp + + + do s = 1,nsites + + ifp = 0 + cpatch => sites(s)%oldest_patch + do while (associated(cpatch)) + + ifp=ifp+1 + + ! zero out various datas + cpatch%ed_parsun_z(:,:,:) = 0._r8 + cpatch%ed_parsha_z(:,:,:) = 0._r8 + cpatch%ed_laisun_z(:,:,:) = 0._r8 + cpatch%ed_laisha_z(:,:,:) = 0._r8 + + bc_out(s)%fsun_pa(ifp) = 0._r8 + sunlai = 0._r8 + shalai = 0._r8 + + ! Loop over patches to calculate laisun_z and laisha_z for each layer. + ! Derive canopy laisun, laisha, and fsun from layer sums. + ! If sun/shade big leaf code, nrad=1 and fsun_z(p,1) and tlai_z(p,1) from + ! SurfaceAlbedo is canopy integrated so that layer value equals canopy value. + + ! cpatch%f_sun is calculated in the surface_albedo routine... + + do CL = 1, cpatch%NCL_p + do FT = 1,numpft_ed + do iv = 1, cpatch%nrad(CL,ft) !NORMAL CASE. + + ! FIX(SPM,040114) - existing comment + ! ** Should this be elai or tlai? Surely we only do radiation for elai? + + cpatch%ed_laisun_z(CL,ft,iv) = cpatch%elai_profile(CL,ft,iv) * & + cpatch%f_sun(CL,ft,iv) + + if ( DEBUG ) write(iulog,*) 'edsurfRad 570 ',cpatch%elai_profile(CL,ft,iv) + if ( DEBUG ) write(iulog,*) 'edsurfRad 571 ',cpatch%f_sun(CL,ft,iv) + + cpatch%ed_laisha_z(CL,ft,iv) = cpatch%elai_profile(CL,ft,iv) * & + (1._r8 - cpatch%f_sun(CL,ft,iv)) + + end do + + !needed for the VOC emissions, etc. + sunlai = sunlai + sum(cpatch%ed_laisun_z(CL,ft,1:cpatch%nrad(CL,ft))) + shalai = shalai + sum(cpatch%ed_laisha_z(CL,ft,1:cpatch%nrad(CL,ft))) + + end do + end do + + if(sunlai+shalai > 0._r8)then + bc_out(s)%fsun_pa(ifp) = sunlai / (sunlai+shalai) + else + bc_out(s)%fsun_pa(ifp) = 0._r8 + endif + + if(bc_out(s)%fsun_pa(ifp) > 1._r8)then + write(iulog,*) 'too much leaf area in profile', bc_out(s)%fsun_pa(ifp), & + cpatch%lai,sunlai,shalai + endif + ! Absorbed PAR profile through canopy ! If sun/shade big leaf code, nrad=1 and fluxes from SurfaceAlbedo ! are canopy integrated so that layer values equal big leaf values. @@ -1044,21 +1058,21 @@ subroutine ED_SunShadeFracs(cpatch,forc_par_d,forc_par_i,fsun) if ( DEBUG ) then write(iulog,*) 'edsurfRad 653 ', cpatch%ed_parsun_z(CL,ft,iv) - write(iulog,*) 'edsurfRad 654 ', forc_par_d - write(iulog,*) 'edsurfRad 655 ', forc_par_i + write(iulog,*) 'edsurfRad 654 ', bc_in(s)%solad_pa(ifp,ipar) + write(iulog,*) 'edsurfRad 655 ', bc_in(s)%solai_pa(ifp,ipar) write(iulog,*) 'edsurfRad 656 ', cpatch%fabd_sun_z(CL,ft,iv) write(iulog,*) 'edsurfRad 657 ', cpatch%fabi_sun_z(CL,ft,iv) endif cpatch%ed_parsun_z(CL,ft,iv) = & - forc_par_d*cpatch%fabd_sun_z(CL,ft,iv) + & - forc_par_i*cpatch%fabi_sun_z(CL,ft,iv) + bc_in(s)%solad_pa(ifp,ipar)*cpatch%fabd_sun_z(CL,ft,iv) + & + bc_in(s)%solai_pa(ifp,ipar)*cpatch%fabi_sun_z(CL,ft,iv) if ( DEBUG )write(iulog,*) 'edsurfRad 663 ', cpatch%ed_parsun_z(CL,ft,iv) cpatch%ed_parsha_z(CL,ft,iv) = & - forc_par_d*cpatch%fabd_sha_z(CL,ft,iv) + & - forc_par_i*cpatch%fabi_sha_z(CL,ft,iv) + bc_in(s)%solad_pa(ifp,ipar)*cpatch%fabd_sha_z(CL,ft,iv) + & + bc_in(s)%solai_pa(ifp,ipar)*cpatch%fabi_sha_z(CL,ft,iv) if ( DEBUG ) write(iulog,*) 'edsurfRad 669 ', cpatch%ed_parsha_z(CL,ft,iv) @@ -1066,9 +1080,14 @@ subroutine ED_SunShadeFracs(cpatch,forc_par_d,forc_par_i,fsun) end do !FT end do !CL - return - - end subroutine ED_SunShadeFracs + cpatch => cpatch%younger + enddo + + + enddo + return + +end subroutine ED_SunShadeFracs ! ! MOVE TO THE INTERFACE diff --git a/main/EDMainMod.F90 b/main/EDMainMod.F90 index 29a97ec261..f920f39709 100755 --- a/main/EDMainMod.F90 +++ b/main/EDMainMod.F90 @@ -6,7 +6,6 @@ module EDMainMod use shr_kind_mod , only : r8 => shr_kind_r8 - use decompMod , only : bounds_type use clm_varctl , only : iulog use atm2lndType , only : atm2lnd_type use SoilStateType , only : soilstate_type @@ -77,7 +76,7 @@ subroutine ed_ecosystem_dynamics(currentSite, & call disturbance_rates(currentSite) ! Integrate state variables from annual rates to daily timestep - call ed_integrate_state_variables(currentSite, soilstate_inst, temperature_inst, waterstate_inst) + call ed_integrate_state_variables(currentSite, temperature_inst ) !****************************************************************************** ! Reproduction, Recruitment and Cohort Dynamics : controls cohort organisation @@ -134,7 +133,7 @@ subroutine ed_ecosystem_dynamics(currentSite, & end subroutine ed_ecosystem_dynamics !-------------------------------------------------------------------------------! - subroutine ed_integrate_state_variables(currentSite, soilstate_inst, temperature_inst, waterstate_inst) + subroutine ed_integrate_state_variables(currentSite, temperature_inst ) ! ! !DESCRIPTION: ! FIX(SPM,032414) refactor so everything goes through interface @@ -143,9 +142,7 @@ subroutine ed_integrate_state_variables(currentSite, soilstate_inst, temperature ! ! !ARGUMENTS: type(ed_site_type) , intent(in) :: currentSite - type(soilstate_type) , intent(in) :: soilstate_inst type(temperature_type) , intent(in) :: temperature_inst - type(waterstate_type) , intent(in) :: waterstate_inst ! ! !LOCAL VARIABLES: type(ed_patch_type) , pointer :: currentPatch @@ -217,7 +214,7 @@ subroutine ed_integrate_state_variables(currentSite, soilstate_inst, temperature write(6,*)'DEBUG18: calling non_canopy_derivs with pno= ',currentPatch%clm_pno endif - call non_canopy_derivs( currentPatch, temperature_inst, soilstate_inst, waterstate_inst ) + call non_canopy_derivs( currentPatch, temperature_inst ) !update state variables simultaneously according to derivatives for this time period. do p = 1,numpft_ed diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index e62a9d75ba..deb69ef1c9 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -15,34 +15,107 @@ module FatesInterfaceMod ! PUBLIC API!!!! ! ------------------------------------------------------------------------------------ - use EDtypesMod , only : ed_site_type + use EDtypesMod , only : ed_site_type, & + numPatchesPerCol + use shr_kind_mod , only : r8 => shr_kind_r8 ! INTERF-TODO: REMOVE THIS + ! ------------------------------------------------------------------------------------ + ! Certain dimension information used by FATES is dictated by the the driver + ! or the host model, or perhaps may be some compromise between what FATES will want + ! in a best-case scenario and what space the driver/host will allow based on its + ! memory constraints (most-likely due to IO) + ! ------------------------------------------------------------------------------------ + + type, private :: fates_dims_type + + integer :: numSWBands ! Maximum number of broad-bands in the short-wave radiation + ! specturm to track + ! (typically 2 as a default, VIS/NIR, in ED variants <2016) + + + end type fates_dims_type + + + + ! ------------------------------------------------------------------------------------ + ! Notes on types + ! For floating point arrays, it is sometimes the convention to define the arrays as + ! POINTER instead of ALLOCATABLE. This usually achieves the same result with subtle + ! differences. POINTER arrays can point to scalar values, discontinuous array slices + ! or alias other variables, ALLOCATABLES cannnot. According to S. Lionel + ! (Intel-Forum Post), ALLOCATABLES are better perfomance wise as long as they point + ! to contiguous memory spaces and do not alias other variables, the case here. + ! ------------------------------------------------------------------------------------ + + type, public :: bc_in_type + + ! The actual number of FATES' ED patches + integer :: npatches + + ! Downwelling direct beam radiation (patch,broad-band) [W/m2?] + real(r8),allocatable :: solad_pa(:,:) + + ! Downwelling diffuse (I-ndirect) radiation (patch,broad-band) [W/m2] + real(r8),allocatable :: solai_pa(:,:) + + + end type bc_in_type + + + type, public :: bc_out_type + + ! Sunlit fraction of the canopy for this patch [0-1] + real(r8),allocatable :: fsun_pa(:) + + end type bc_out_type + + type, public :: fates_interface_type ! This is the root of the ED/FATES hierarchy of instantaneous state variables - ! ie the root of the linked lists. Each path list is currently associated - ! with a grid-cell, this is intended to be migrated to columns - ! prev: type(ed_site_type)::ed_allsites_inst + ! ie the root of the linked lists. Each path list is currently associated with a + ! grid-cell, this is intended to be migrated to columns integer :: nsites type(ed_site_type), allocatable :: sites(:) + + ! These are boundary conditions that the FATES models are required to be filled. + ! These values are filled by the driver or HLM. Once filled, these have an + ! intent(in) status. Each site has a derived type structure, which may include + ! a scalar for site level data, a patch vector, potentially cohort vectors (but + ! not yet atm) and other dimensions such as soil-depth or pft. These vectors + ! are initialized by maximums, and the allocations are static in time to avoid + ! having to allocate/de-allocate memory + + type(bc_in_type), allocatable :: bc_in(:) + + ! These are the boundary conditions that the FATES model returns to its HLM or + ! driver. It has the same allocation strategy and similar vector types. - ! INTERF-TODO ADD THE DLM->FATES BOUNDARY CONDITION CLASS - ! These are boundary condition variables populated by the DLM - ! type(fates_bc_type) :: fatesbc - + type(bc_out_type), allocatable :: bc_out(:) + contains - ! Procedures for initializing FATES threaded memory and communicators -! procedure, public :: fates_clean + procedure, public :: allocate_bcs + procedure, public :: zero_bcs end type fates_interface_type + ! ------------------------------------------------------------------------------------ + ! Dimension information is independent of which clump it is on + ! and since these are typically read only variables and never updated on the clump + ! we need not attach these to the threaded instances of the fates_interface_type + ! ------------------------------------------------------------------------------------ + + type(fates_dims_type) :: fates_dims + + + contains - ! ------------------------------------------------------------------------------------ + ! ==================================================================================== ! INTERF-TODO: THIS IS A PLACE-HOLDER ROUTINE, NOT CALLED YET... subroutine fates_clean(this) @@ -62,4 +135,52 @@ subroutine fates_clean(this) return end subroutine fates_clean + + ! ==================================================================================== + + + subroutine allocate_bcs(this,s) + + ! --------------------------------------------------------------------------------- + ! Allocate and Initialze the FATES boundary condition vectors + ! --------------------------------------------------------------------------------- + + implicit none + class(fates_interface_type), intent(inout) :: this + integer, intent(in) :: s + + ! Allocate input boundaries + + allocate(this%bc_in(s)%solad_pa(numPatchesPerCol,fates_dims%numSWBands)) + allocate(this%bc_in(s)%solai_pa(numPatchesPerCol,fates_dims%numSWBands)) + + ! Allocate output boundaries + + allocate(this%bc_out(s)%fsun_pa(numPatchesPerCol)) + + return + end subroutine allocate_bcs + + ! ==================================================================================== + + subroutine zero_bcs(this,s) + + implicit none + class(fates_interface_type), intent(inout) :: this + integer, intent(in) :: s + + ! Input boundaries + + this%bc_in(s)%solad_pa(:,:) = 0.0_r8 + this%bc_in(s)%solai_pa(:,:) = 0.0_r8 + + ! Output boundaries + + this%bc_out(s)%fsun_pa(:) = 0.0_r8 + + + return + end subroutine zero_bcs + + end module FatesInterfaceMod From 91ff6e51191e2508b654363aaba7fc784c0e2825 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 21 Jun 2016 18:43:55 -0700 Subject: [PATCH 3/9] incremental changes to bc passing, mostly related to btran --- biogeophys/EDBtranMod.F90 | 172 +++++++++++++++++++++++-------------- main/FatesInterfaceMod.F90 | 57 +++++++++--- 2 files changed, 154 insertions(+), 75 deletions(-) diff --git a/biogeophys/EDBtranMod.F90 b/biogeophys/EDBtranMod.F90 index 030eb4f2b2..e24ce418dc 100644 --- a/biogeophys/EDBtranMod.F90 +++ b/biogeophys/EDBtranMod.F90 @@ -1,22 +1,18 @@ module EDBtranMod + + !------------------------------------------------------------------------------------- + ! Description: + ! + ! ------------------------------------------------------------------------------------ - !------------------------------------------------------------------------------ - ! !DESCRIPTION: - ! This routine accumulates NPP, GPP and respiration of each cohort over the course of each 24 hour period. - ! The fluxes are stored per cohort, and the npp_clm (etc) fluxes are calcualted in EDPhotosynthesis - ! This routine cannot be in EDPhotosynthesis because EDPhotosynthesis is a loop and therefore would - ! erroneously add these things up multiple times. - ! Rosie Fisher. March 2014. - ! - ! !USES: use pftconMod , only : pftcon use EDTypesMod , only : ed_patch_type, ed_cohort_type, numpft_ed use EDEcophysContype , only : EDecophyscon ! implicit none private - ! - public :: BTRAN_ED + + public :: btran_ed ! type(ed_cohort_type), pointer :: currentCohort ! current cohort type(ed_patch_type) , pointer :: currentPatch ! current patch @@ -25,37 +21,21 @@ module EDBtranMod contains !------------------------------------------------------------------------------ - subroutine btran_ed( bounds, p, sites, nsites, hsites, & - soilstate_inst, waterstate_inst, temperature_inst, energyflux_inst) - ! - ! !DESCRIPTION: - ! - ! !USES: + subroutine btran_ed( sites, nsites, bc_in, bc_out) + + use shr_kind_mod , only : r8 => shr_kind_r8 - use shr_const_mod , only : shr_const_pi - use decompMod , only : bounds_type - use clm_varpar , only : nlevgrnd - use clm_varctl , only : iulog - use clm_varcon , only : tfrz, denice, denh2o - use SoilStateType , only : soilstate_type - use WaterStateType , only : waterstate_type - use TemperatureType , only : temperature_type - use EnergyFluxType , only : energyflux_type - use GridcellType , only : grc - use ColumnType , only : col - use PatchType , only : patch - use EDTypesMod , only : ed_site_type, map_clmpatch_to_edpatch - ! - ! !ARGUMENTS - type(bounds_type) , intent(in) :: bounds ! clump bounds - integer , intent(in) :: p ! patch/'p' - type(ed_site_type) , intent(inout), target :: sites(nsites) - integer , intent(in) :: nsites - integer , intent(in) :: hsites(bounds%begc:bounds%endc) - type(soilstate_type) , intent(inout) :: soilstate_inst - type(waterstate_type) , intent(in) :: waterstate_inst - type(temperature_type) , intent(in) :: temperature_inst - type(energyflux_type) , intent(inout) :: energyflux_inst + use EDtypesMod , only : ed_patch_type, ed_site_type + use FatesInterfaceMod , only : bc_in_type, & + bc_out_type + + ! Arguments + + type(ed_site_type),intent(inout),target :: sites(nsites) + integer,intent(in) :: nsites + type(bc_in_type),intent(in) :: bc_in(nsites) + type(bc_out_type),intent(out) :: bc_out(nsites) + ! ! !LOCAL VARIABLES: integer :: iv !leaf layer @@ -63,6 +43,11 @@ subroutine btran_ed( bounds, p, sites, nsites, hsites, & integer :: c !column integer :: j !soil layer integer :: ft ! plant functional type index + + + ! The root effective root fraction for different pfts and soil layers + real(r8) :: rootr(numpft_ed,nlevgrnd) + !---------------------------------------------------------------------- ! Inputs to model from CLM. To be read in through an input file for the purposes of this program. @@ -113,38 +98,94 @@ subroutine btran_ed( bounds, p, sites, nsites, hsites, & real(r8) :: temprootr !------------------------------------------------------------------------------ - associate(& - dz => col%dz , & ! Input: [real(r8) (:,:) ] layer depth (m) - smpso => pftcon%smpso , & ! Input: soil water potential at full stomatal opening (mm) - smpsc => pftcon%smpsc , & ! Input: soil water potential at full stomatal closure (mm) + do s = 1,nsites + + ifp = 0 + cpatch => sites(s)%oldest_patch + do while (associated(cpatch)) + ifp=ifp+1 + + do FT = 1,numpft_ed + currentPatch%btran_ft(FT) = 0.0_r8 + do j = 1,nlevgrnd + + if(bc_in(s)%smp_sl(j) > -900.0_r8 ) then ! The flag for frozen or no water is -999 + + smp_node = max(smpsc(FT), bc_in(s)%smp_sl(j)) + + rresis_ft = min( (bc_in(s)%eff_porosity_sl(j)/bc_in(s)%watsat(j))* & + (smp_node - smpsc(FT)) / (smpso(FT) - smpsc(FT)), 1._r8) + + + currentPatch%rootr_ft(FT,j) = currentPatch%rootfr_ft(FT,j)*rresis_ft + + ! root water uptake is not linearly proportional to root density, + ! to allow proper deep root funciton. Replace with equations from SPA/Newman. FIX(RF,032414) + ! currentPatch%rootr_ft(FT,j) = currentPatch%rootfr_ft(FT,j)**0.3*rresis_FT(FT,j)/ & + ! sum(currentPatch%rootfr_ft(FT,1:nlevgrnd)**0.3) + currentPatch%btran_ft(FT) = currentPatch%btran_ft(FT) + currentPatch%rootr_ft(FT,j) + + else + currentPatch%rootr_ft(FT,j) = 0._r8 + end if - sucsat => soilstate_inst%sucsat_col , & ! Input: [real(r8) (:,:) ] minimum soil suction (mm) - watsat => soilstate_inst%watsat_col , & ! Input: [real(r8) (:,:) ] volumetric soil water at saturation (porosity) - bsw => soilstate_inst%bsw_col , & ! Input: [real(r8) (:,:) ] Clapp and Hornberger "b" - sand => soilstate_inst%sandfrac_patch , & ! Input: [real(r8) (:) ] % sand of soil - rootr => soilstate_inst%rootr_patch , & ! Output: [real(r8) (:,:) ] Fraction of water uptake in each layer + end do !j - h2osoi_ice => waterstate_inst%h2osoi_ice_col , & ! Input: [real(r8) (:,:) ] ice lens (kg/m2) - h2osoi_vol => waterstate_inst%h2osoi_vol_col , & ! Input: [real(r8) (:,:) ] volumetric soil water (0<=h2osoi_vol<=watsat) [m3/m3] - h2osoi_liq => waterstate_inst%h2osoi_liq_col , & ! Input: [real(r8) (:,:) ] liquid water (kg/m2) + ! Normalize root resistances to get layer contribution to ET + do j = 1,nlevgrnd + if (currentPatch%btran_ft(FT) > 0.0_r8) then + currentPatch%rootr_ft(FT,j) = currentPatch%rootr_ft(FT,j)/currentPatch%btran_ft(FT) + else + currentPatch%rootr_ft(FT,j) = 0._r8 + end if + end do - t_soisno => temperature_inst%t_soisno_col , & ! Input: [real(r8) (:,:) ] soil temperature (Kelvin) + end do !PFT - btran => energyflux_inst%btran_patch , & ! Output: [real(r8) (:) ] transpiration wetness factor (0 to 1) - rresis => energyflux_inst%rresis_patch & ! Output: [real(r8) (:,:) ] root resistance by layer (0-1) (nlevgrnd) - ) - - if (patch%is_veg(p)) then + ! PFT-averaged point level root fraction for extraction purposese. + ! This probably needs to be weighted by actual transpiration from each pft. FIX(RF,032414). + pftgs(:) = 0._r8 + currentCohort => currentPatch%tallest + do while(associated(currentCohort)) + pftgs(currentCohort%pft) = pftgs(currentCohort%pft) + currentCohort%gscan * currentCohort%n + currentCohort => currentCohort%shorter + enddo + + + ! Process the boundary output, this is necessary for calculating the soil-moisture + ! sink term across the different layers in driver/host. Photosynthesis will + ! pass the host a total transpiration for the patch. This needs rootr to be + ! distributed over the soil layers. - c = patch%column(p) - s = hsites(c) + do j = 1,nlevgrnd + bc_out(s)%rootr_pa(ifp,j) = 0._r8 + do FT = 1,numpft_ed + if(sum(pftgs) > 0._r8)then !prevent problem with the first timestep - might fail + !bit-retart test as a result? FIX(RF,032414) + bc_out(s)%rootr_pa(ifp,j) = bc_out(s)%rootr_pa(ifp,j) + & + currentPatch%rootr_ft(FT,j) * pftgs(ft)/sum(pftgs) + else + bc_out(s)%rootr_pa(ifp,j) = bc_out(s)%rootr_pa(ifp,j) + & + currentPatch%rootr_ft(FT,j) * 1./numpft_ed + end if + enddo + enddo - currentPatch => map_clmpatch_to_edpatch(sites(s), p) - do FT = 1,numpft_ed - currentPatch%btran_ft(FT) = 0.0_r8 - do j = 1,nlevgrnd + + cpatch => cpatch%younger + end do + + end do + + + + + + do FT = 1,numpft_ed + currentPatch%btran_ft(FT) = 0.0_r8 + do j = 1,nlevgrnd !Root resistance factors vol_ice = min(watsat(c,j), h2osoi_ice(c,j)/(dz(c,j)*denice)) @@ -199,6 +240,7 @@ subroutine btran_ed( bounds, p, sites, nsites, hsites, & end if enddo enddo + !--------------------------------------------------------------------------------------- diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index deb69ef1c9..200271e876 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -18,7 +18,7 @@ module FatesInterfaceMod use EDtypesMod , only : ed_site_type, & numPatchesPerCol use shr_kind_mod , only : r8 => shr_kind_r8 ! INTERF-TODO: REMOVE THIS - + use clm_varpar , only : nlevgrnd ! ------------------------------------------------------------------------------------ ! Certain dimension information used by FATES is dictated by the the driver @@ -32,8 +32,7 @@ module FatesInterfaceMod integer :: numSWBands ! Maximum number of broad-bands in the short-wave radiation ! specturm to track ! (typically 2 as a default, VIS/NIR, in ED variants <2016) - - + end type fates_dims_type @@ -53,12 +52,20 @@ module FatesInterfaceMod ! The actual number of FATES' ED patches integer :: npatches - ! Downwelling direct beam radiation (patch,broad-band) [W/m2?] - real(r8),allocatable :: solad_pa(:,:) + ! Downwelling direct beam radiation (patch,broad-band) [W/m2] + real(r8), allocatable :: solad_pa(:,:) ! Downwelling diffuse (I-ndirect) radiation (patch,broad-band) [W/m2] - real(r8),allocatable :: solai_pa(:,:) + real(r8), allocatable :: solai_pa(:,:) + ! Soil suction potential of layers in each site, negative, [mm] + real(r8), allocatable :: smp_sl(:) + + ! Effective porosity = porosity - vol_ic, of layers in each site [-] + real(r8), allocatable :: eff_porosity_sl(:) + + ! volumetric soil water at saturation (porosity) + real(r8), allocatable :: watsat_sl(:) end type bc_in_type @@ -68,6 +75,18 @@ module FatesInterfaceMod ! Sunlit fraction of the canopy for this patch [0-1] real(r8),allocatable :: fsun_pa(:) + ! Root soil water stress (resistance) by layer + ! (diagnostic, should not be used by HLM) + real(r8),allocatable :: rresis_pa(:,:) + + ! Effective fraction of roots in each soil layer + ! (diagnostic, should not be used by HLM) + real(r8), allocatable :: rootr_pa(:,:) + + ! Integrated (vertically) transpiration wetness factor (0 to 1) + ! (diagnostic, should not be used by HLM) + real(r8), allocatable :: btran_pa(:) + end type bc_out_type @@ -151,12 +170,25 @@ subroutine allocate_bcs(this,s) ! Allocate input boundaries + ! Radiation allocate(this%bc_in(s)%solad_pa(numPatchesPerCol,fates_dims%numSWBands)) allocate(this%bc_in(s)%solai_pa(numPatchesPerCol,fates_dims%numSWBands)) + + ! Hydrology + allocate(this%bc_in(s)%smp_sl(nlevgrnd)) + allocate(this%bc_in(s)%eff_porosity_sl(nlevgrnd)) + allocate(this%bc_in(s)%watsat_sl(nlevgrnd)) ! Allocate output boundaries - + + ! Radiation allocate(this%bc_out(s)%fsun_pa(numPatchesPerCol)) + + ! Hydrology + allocate(this%bc_out(s)%rresis_pa(numPatchesPerCol,nlevgrnd)) + allocate(this%bc_out(s)%rootr_pa(numPatchesPerCol,nlevgrnd)) + allocate(this%bc_out(s)%btran_pa(numPatchesPerCol)) + return end subroutine allocate_bcs @@ -171,13 +203,18 @@ subroutine zero_bcs(this,s) ! Input boundaries - this%bc_in(s)%solad_pa(:,:) = 0.0_r8 - this%bc_in(s)%solai_pa(:,:) = 0.0_r8 + this%bc_in(s)%solad_pa(:,:) = 0.0_r8 + this%bc_in(s)%solai_pa(:,:) = 0.0_r8 + this%bc_in(s)%smp_sl(:) = 0.0_r8 + this%bc_in(s)%eff_porosity_sl(:) = 0.0_r8 + this%bc_in(s)%watsat_sl(:) = 0.0_r8 ! Output boundaries this%bc_out(s)%fsun_pa(:) = 0.0_r8 - + this%bc_out(s)%rresis_pa(:) = 0.0_r8 + this%bc_out(s)%rootr_pa(:) = 0.0_r8 + this%bc_out(s)%btran_pa(:) = 0.0_r8 return end subroutine zero_bcs From ff2d1ff7f9e1406da20b2f77222695317b98de9f Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Tue, 28 Jun 2016 08:52:22 -0700 Subject: [PATCH 4/9] boundary condition work on Fates: prototype working for 1) allocation of bc structures, 2) use of bc structures in sunshade_fracs and btran, and 3) use of control parameters sent from host through a dictionary to allocate bc structures --- biogeophys/EDBtranMod.F90 | 625 +++++++++++++----------------- biogeophys/EDSurfaceAlbedoMod.F90 | 19 +- main/EDTypesMod.F90 | 22 ++ main/FatesInterfaceMod.F90 | 191 ++++++--- 4 files changed, 438 insertions(+), 419 deletions(-) diff --git a/biogeophys/EDBtranMod.F90 b/biogeophys/EDBtranMod.F90 index e24ce418dc..d9f5b464f8 100644 --- a/biogeophys/EDBtranMod.F90 +++ b/biogeophys/EDBtranMod.F90 @@ -4,384 +4,291 @@ module EDBtranMod ! Description: ! ! ------------------------------------------------------------------------------------ - - use pftconMod , only : pftcon - use EDTypesMod , only : ed_patch_type, ed_cohort_type, numpft_ed - use EDEcophysContype , only : EDecophyscon - ! - implicit none - private - - public :: btran_ed - ! - type(ed_cohort_type), pointer :: currentCohort ! current cohort - type(ed_patch_type) , pointer :: currentPatch ! current patch - !------------------------------------------------------------------------------ + + use pftconMod , only : pftcon + use EDTypesMod , only : ed_site_type, & + ed_patch_type, & + ed_cohort_type, & + numpft_ed, & + ctrl_parms + use shr_kind_mod , only : r8 => shr_kind_r8 + use FatesInterfaceMod , only : bc_in_type, & + bc_out_type + use clm_varctl , only : iulog !INTERF-TODO: THIS SHOULD BE MOVED + + ! + implicit none + private + + public :: btran_ed + ! + type(ed_cohort_type), pointer :: currentCohort ! current cohort + type(ed_patch_type) , pointer :: currentPatch ! current patch contains - - !------------------------------------------------------------------------------ - subroutine btran_ed( sites, nsites, bc_in, bc_out) - - - use shr_kind_mod , only : r8 => shr_kind_r8 - use EDtypesMod , only : ed_patch_type, ed_site_type - use FatesInterfaceMod , only : bc_in_type, & - bc_out_type - - ! Arguments - - type(ed_site_type),intent(inout),target :: sites(nsites) - integer,intent(in) :: nsites - type(bc_in_type),intent(in) :: bc_in(nsites) - type(bc_out_type),intent(out) :: bc_out(nsites) - - ! - ! !LOCAL VARIABLES: - integer :: iv !leaf layer - integer :: s !site - integer :: c !column - integer :: j !soil layer - integer :: ft ! plant functional type index - - - ! The root effective root fraction for different pfts and soil layers - real(r8) :: rootr(numpft_ed,nlevgrnd) - - !---------------------------------------------------------------------- - - ! Inputs to model from CLM. To be read in through an input file for the purposes of this program. - integer, parameter :: nv = 5 ! Number of canopy layers - real(r8) :: xksat ! maximum hydraulic conductivity of soil [mm/s] - real(r8) :: s1 ! HC intermediate - real(r8) :: swp_mpa(nlevgrnd) ! matrix potential - MPa - real(r8) :: hk(nlevgrnd) ! hydraulic conductivity [mm h2o/s] - real(r8) :: rootxsecarea ! root X-sectional area (m2) - real(r8) :: rootmass(nlevgrnd) ! root mass in each layer (g) - real(r8) :: rootlength(nlevgrnd) ! root length in each layer (m) - real(r8) :: soilr1(nlevgrnd) ! soil-to-root resistance in each layer (MPa s m2 mmol-1) - real(r8) :: soilr2(nlevgrnd) ! internal root resistance in each layer (MPa s m2 mmol-1) - real(r8) :: rs ! intermediate variable - real(r8) :: soilr_z(nlevgrnd) ! soil-to-xylem resistance in each layer (MPa s m2 mmol-1) - real(r8) :: lsoil(nlevgrnd) ! hydraulic conductivity in each soil layer - - real(r8) :: estevap(nlevgrnd) ! potential suction from each soil layer (mmol m-2 s-1) - real(r8) :: totestevap ! potential suction from each soil layer (mmol m-2 s-1) - real(r8) :: fraction_uptake(nlevgrnd) ! Uptake of water from each soil layer (-) - real(r8) :: maxevap(nlevgrnd) ! potential suction from each soil layer (mmol m-2 s-1) - real(r8) :: totmaxevap ! potential suction from each soil layer (mmol m-2 s-1) - real(r8) :: fleaf ! fraction of leaves in each canopy layer - - ! Model parameters - real(r8) :: head = 0.009807_r8 ! head of pressure (MPa/m) - real(r8) :: rootdens = 0.5e6_r8 ! root density, g biomass m-3 root - real(r8) :: pi = shr_const_pi - real(r8) :: vol_ice ! partial volume of ice lens in layer - real(r8) :: eff_porosity ! effective porosity in layer - real(r8) :: vol_liq ! partial volume of liquid water in layer - real(r8) :: s_node ! vol_liq/eff_porosity - real(r8) :: smp_node ! matrix potential - - ! To be read in from pft file ultimately. - real(r8) :: minlwp = -2.5_r8 ! minimum leaf water potential in MPa - real(r8) :: rootrad = 0.001_r8 ! root radius in metres - - ! Outputs to CLM_SPA - real(r8) :: weighted_SWP ! weighted apparent soil water potential: MPa. - real(r8) :: canopy_soil_resistance(nv) ! Resistance experienced by each canopy layer: MPa s m2 mmol-1 - - ! SPA Pointers from CLM type. - logical, parameter :: SPA_soil=.false. ! Is the BTRAN model SPA or CLM? FIX(SPM,032414) ed - make this a namelist var - - real(r8) :: rresis_ft(numpft_ed,nlevgrnd) ! resistance to water uptake per pft and soil layer. - real(r8) :: pftgs(numpft_ed) ! pft weighted stomatal conductance s/m - real(r8) :: temprootr - !------------------------------------------------------------------------------ - - - do s = 1,nsites - - ifp = 0 - cpatch => sites(s)%oldest_patch - do while (associated(cpatch)) - ifp=ifp+1 - - do FT = 1,numpft_ed - currentPatch%btran_ft(FT) = 0.0_r8 - do j = 1,nlevgrnd - - if(bc_in(s)%smp_sl(j) > -900.0_r8 ) then ! The flag for frozen or no water is -999 - - smp_node = max(smpsc(FT), bc_in(s)%smp_sl(j)) - - rresis_ft = min( (bc_in(s)%eff_porosity_sl(j)/bc_in(s)%watsat(j))* & - (smp_node - smpsc(FT)) / (smpso(FT) - smpsc(FT)), 1._r8) - - - currentPatch%rootr_ft(FT,j) = currentPatch%rootfr_ft(FT,j)*rresis_ft - - ! root water uptake is not linearly proportional to root density, - ! to allow proper deep root funciton. Replace with equations from SPA/Newman. FIX(RF,032414) - ! currentPatch%rootr_ft(FT,j) = currentPatch%rootfr_ft(FT,j)**0.3*rresis_FT(FT,j)/ & - ! sum(currentPatch%rootfr_ft(FT,1:nlevgrnd)**0.3) - currentPatch%btran_ft(FT) = currentPatch%btran_ft(FT) + currentPatch%rootr_ft(FT,j) - - else - currentPatch%rootr_ft(FT,j) = 0._r8 - end if - - end do !j - - ! Normalize root resistances to get layer contribution to ET - do j = 1,nlevgrnd - if (currentPatch%btran_ft(FT) > 0.0_r8) then - currentPatch%rootr_ft(FT,j) = currentPatch%rootr_ft(FT,j)/currentPatch%btran_ft(FT) - else - currentPatch%rootr_ft(FT,j) = 0._r8 - end if - end do - - end do !PFT - - ! PFT-averaged point level root fraction for extraction purposese. - ! This probably needs to be weighted by actual transpiration from each pft. FIX(RF,032414). - pftgs(:) = 0._r8 - currentCohort => currentPatch%tallest - do while(associated(currentCohort)) - pftgs(currentCohort%pft) = pftgs(currentCohort%pft) + currentCohort%gscan * currentCohort%n - currentCohort => currentCohort%shorter - enddo - - - ! Process the boundary output, this is necessary for calculating the soil-moisture - ! sink term across the different layers in driver/host. Photosynthesis will - ! pass the host a total transpiration for the patch. This needs rootr to be - ! distributed over the soil layers. - - do j = 1,nlevgrnd - bc_out(s)%rootr_pa(ifp,j) = 0._r8 - do FT = 1,numpft_ed - if(sum(pftgs) > 0._r8)then !prevent problem with the first timestep - might fail - !bit-retart test as a result? FIX(RF,032414) - bc_out(s)%rootr_pa(ifp,j) = bc_out(s)%rootr_pa(ifp,j) + & - currentPatch%rootr_ft(FT,j) * pftgs(ft)/sum(pftgs) - else - bc_out(s)%rootr_pa(ifp,j) = bc_out(s)%rootr_pa(ifp,j) + & - currentPatch%rootr_ft(FT,j) * 1./numpft_ed - end if - enddo - enddo - - - - cpatch => cpatch%younger - end do - - end do - - - - - do FT = 1,numpft_ed - currentPatch%btran_ft(FT) = 0.0_r8 - do j = 1,nlevgrnd - - !Root resistance factors - vol_ice = min(watsat(c,j), h2osoi_ice(c,j)/(dz(c,j)*denice)) - eff_porosity = watsat(c,j)-vol_ice - vol_liq = min(eff_porosity, h2osoi_liq(c,j)/(dz(c,j)*denh2o)) - if (vol_liq <= 0._r8 .or. t_soisno(c,j) <= tfrz-2._r8) then - currentPatch%rootr_ft(FT,j) = 0._r8 - else - s_node = max(vol_liq/eff_porosity,0.01_r8) - smp_node = max(smpsc(FT), -sucsat(c,j)*s_node**(-bsw(c,j))) - !FIX(RF,032414) for junipers - rresis_ft(FT,j) = min( (eff_porosity/watsat(c,j))* & - (smp_node - smpsc(FT)) / (smpso(FT) - smpsc(FT)), 1._r8) - - currentPatch%rootr_ft(FT,j) = currentPatch%rootfr_ft(FT,j)*rresis_FT(FT,j) - ! root water uptake is not linearly proportional to root density, - ! to allow proper deep root funciton. Replace with equations from SPA/Newman. FIX(RF,032414) - ! currentPatch%rootr_ft(FT,j) = currentPatch%rootfr_ft(FT,j)**0.3*rresis_FT(FT,j)/ & - ! sum(currentPatch%rootfr_ft(FT,1:nlevgrnd)**0.3) - currentPatch%btran_ft(FT) = currentPatch%btran_ft(FT) + currentPatch%rootr_ft(FT,j) - end if - end do !j - - ! Normalize root resistances to get layer contribution to ET - do j = 1,nlevgrnd - if (currentPatch%btran_ft(FT) > 0.0_r8) then - currentPatch%rootr_ft(FT,j) = currentPatch%rootr_ft(FT,j)/currentPatch%btran_ft(FT) - else - currentPatch%rootr_ft(FT,j) = 0._r8 - end if - end do - - end do !PFT - - ! PFT-averaged point level root fraction for extraction purposese. - ! This probably needs to be weighted by actual transpiration from each pft. FIX(RF,032414). - pftgs(:) = 0._r8 - currentCohort => currentPatch%tallest - do while(associated(currentCohort)) - pftgs(currentCohort%pft) = pftgs(currentCohort%pft) + currentCohort%gscan * currentCohort%n - currentCohort => currentCohort%shorter - enddo - - do j = 1,nlevgrnd - rootr(p,j) = 0._r8 - do FT = 1,numpft_ed - if(sum(pftgs) > 0._r8)then !prevent problem with the first timestep - might fail - !bit-retart test as a result? FIX(RF,032414) - rootr(p,j) = rootr(p,j) + currentPatch%rootr_ft(FT,j) * pftgs(ft)/sum(pftgs) - else - rootr(p,j) = rootr(p,j) + currentPatch%rootr_ft(FT,j) * 1./numpft_ed - end if - enddo - enddo - - + ! ==================================================================================== + + subroutine btran_ed( sites, nsites, bc_in, bc_out) + + ! --------------------------------------------------------------------------------- + ! Calculate the transpiration wetness function (BTRAN) and the root uptake + ! distribution (ROOTR). + ! Boundary conditions in: bc_in(s)%eff_porosity_sl(j) unfrozen porosity + ! bc_in(s)%watsat_sl(j) porosity + ! bc_in(s)%active_uptake_sl(j) frozen/not frozen + ! bc_in(s)%smp_sl(j) suction + ! Boundary conditions out: bc_out(s)%rootr_pa root uptake distribution + ! bc_out(s)%btran_pa wetness factor + ! --------------------------------------------------------------------------------- + + + + ! Arguments + + type(ed_site_type),intent(inout),target :: sites(nsites) + integer,intent(in) :: nsites + type(bc_in_type),intent(in) :: bc_in(nsites) + type(bc_out_type),intent(inout) :: bc_out(nsites) + + ! + ! !LOCAL VARIABLES: + type(ed_patch_type),pointer :: cpatch ! Current Patch Pointer + integer :: s ! site + integer :: j ! soil layer + integer :: ifp ! patch vector index for the site + integer :: ft ! plant functional type index + real(r8) :: smp_node ! matrix potential + real(r8) :: rresis ! suction limitation to transpiration independent + ! of root density + real(r8) :: pftgs(numpft_ed) ! pft weighted stomatal conductance s/m + real(r8) :: temprootr + !------------------------------------------------------------------------------ + + associate( & + numlevgrnd => ctrl_parms%numlevgrnd , & + smpsc => pftcon%smpsc , & ! INTERF-TODO: THESE SHOULD BE FATES PARAMETERS + smpso => pftcon%smpso & ! INTERF-TODO: THESE SHOULD BE FATES PARAMETERS + ) + + do s = 1,nsites + + ifp = 0 + cpatch => sites(s)%oldest_patch + do while (associated(cpatch)) + ifp=ifp+1 + + ! THIS SHOULD REALLY BE A COHORT LOOP ONCE WE HAVE rootfr_ft FOR COHORTS (RGK) + + do ft = 1,numpft_ed + cpatch%btran_ft(ft) = 0.0_r8 + do j = 1,numlevgrnd + + ! Calculations are only relevant where liquid water exists + ! see clm_fates%wrap_btran for calculation with CLM/ALM + + if( bc_in(s)%active_uptake_sl(j) ) then + + smp_node = max(smpsc(ft), bc_in(s)%smp_sl(j)) + + rresis = min( (bc_in(s)%eff_porosity_sl(j)/bc_in(s)%watsat_sl(j))* & + (smp_node - smpsc(ft)) / (smpso(ft) - smpsc(ft)), 1._r8) + + cpatch%rootr_ft(ft,j) = cpatch%rootfr_ft(ft,j)*rresis + + ! root water uptake is not linearly proportional to root density, + ! to allow proper deep root funciton. Replace with equations from SPA/Newman. FIX(RF,032414) + ! cpatch%rootr_ft(ft,j) = cpatch%rootfr_ft(ft,j)**0.3*rresis_ft(ft,j)/ & + ! sum(cpatch%rootfr_ft(ft,1:nlevgrnd)**0.3) + cpatch%btran_ft(ft) = cpatch%btran_ft(ft) + cpatch%rootr_ft(ft,j) + + else + cpatch%rootr_ft(ft,j) = 0._r8 + end if + + end do !j + + ! Normalize root resistances to get layer contribution to ET + do j = 1,numlevgrnd + if (cpatch%btran_ft(ft) > 0.0_r8) then + cpatch%rootr_ft(ft,j) = cpatch%rootr_ft(ft,j)/cpatch%btran_ft(ft) + else + cpatch%rootr_ft(ft,j) = 0._r8 + end if + end do + + end do !PFT + + ! PFT-averaged point level root fraction for extraction purposese. + ! This probably needs to be weighted by actual transpiration from each pft. FIX(RF,032414). + pftgs(:) = 0._r8 + currentCohort => cpatch%tallest + do while(associated(currentCohort)) + pftgs(currentCohort%pft) = pftgs(currentCohort%pft) + currentCohort%gscan * currentCohort%n + currentCohort => currentCohort%shorter + enddo + + ! Process the boundary output, this is necessary for calculating the soil-moisture + ! sink term across the different layers in driver/host. Photosynthesis will + ! pass the host a total transpiration for the patch. This needs rootr to be + ! distributed over the soil layers. + + do j = 1,numlevgrnd + bc_out(s)%rootr_pa(ifp,j) = 0._r8 + do ft = 1,numpft_ed + if(sum(pftgs) > 0._r8)then !prevent problem with the first timestep - might fail + !bit-retart test as a result? FIX(RF,032414) + bc_out(s)%rootr_pa(ifp,j) = bc_out(s)%rootr_pa(ifp,j) + & + cpatch%rootr_ft(ft,j) * pftgs(ft)/sum(pftgs) + else + bc_out(s)%rootr_pa(ifp,j) = bc_out(s)%rootr_pa(ifp,j) + & + cpatch%rootr_ft(ft,j) * 1./numpft_ed + end if + enddo + enddo + + !weight patch level output BTRAN for the + bc_out(s)%btran_pa(ifp) = 0.0_r8 + do ft = 1,numpft_ed + if(sum(pftgs) > 0._r8)then !prevent problem with the first timestep - might fail + !bit-retart test as a result? FIX(RF,032414) + bc_out(s)%btran_pa(ifp) = bc_out(s)%btran_pa(ifp) + cpatch%btran_ft(ft) * pftgs(ft)/sum(pftgs) + else + bc_out(s)%btran_pa(ifp) = bc_out(s)%btran_pa(ifp) + cpatch%btran_ft(ft) * 1./numpft_ed + end if + enddo + + + + temprootr = sum(bc_out(s)%rootr_pa(ifp,:)) + if(temprootr /= 1.0_r8)then + write(iulog,*) 'error with rootr in canopy fluxes',temprootr + if(temprootr > 0._r8)then + do j = 1,numlevgrnd + bc_out(s)%rootr_pa(ifp,j) = bc_out(s)%rootr_pa(ifp,j)/temprootr + enddo + end if + end if + + cpatch => cpatch%younger + end do + end do + + end associate + + end subroutine btran_ed + + ! ========================================================================================= !--------------------------------------------------------------------------------------- ! SPA based recalculation of BTRAN and water uptake. !--------------------------------------------------------------------------------------- - if (SPA_soil) then ! normal case don't run this. - rootr(p,:) = 0._r8 - do FT = 1,numpft_ed - - ! Soil Physics - do j = 1,nlevgrnd - ! CLM water retention curve. Clapp and Hornberger equation. - s1 = max(h2osoi_vol(c,j)/watsat(c,j), 0.01_r8) - s1 = min(1.0_r8,s1) - smp_node = -sucsat(c,j)*s1**(-bsw(c,j)) - swp_mpa(j) = smp_node *10.0_r8/1000000.0_r8 !convert from mm to Mpa - - ! CLM hydraulic conductivity curve. - ! As opposed to the Richard's equation solution in SoilHydrology.Mod - ! the conductivity here is defined in the middle of the layer in question, not at the edge... - xksat = 0.0070556_r8 * (10._r8**(-0.884_r8+0.0153_r8*sand(p)) ) - hk(j) = xksat*s1**(2._r8*bsw(c,j)+2._r8) !removed the ice from here to avoid 1st ts crashing - enddo - - ! Root resistance - rootxsecarea=3.14159*rootrad**2 - do j = 1,nlevgrnd - rootmass(j) = EDecophyscon%soilbeta(FT) * currentPatch%rootfr_ft(FT,j) - rootlength(j) = rootmass(j)/(rootdens*rootxsecarea) !m m-3 soil - Lsoil(j) = hk(j)/1000/head !converts from mms-1 to ms-1 and then to m2 s-1 MPa-1 - if(Lsoil(j) < 1e-35_r8.or.currentPatch%rootfr_ft(ft,j) <= 0.0_r8)then !prevent floating point error - soilr_z(j) = 1e35_r8 - soilr2(j) = 1e35_r8 - else - ! Soil-to-root water uptake from Newman (1969). - rs = sqrt (1._r8 / (rootlength(j) * pi)) - soilr1(j) = log(rs/rootrad) / (2.0_r8 * pi * rootlength(j) * Lsoil(j) * dz(c,j)) - ! convert from MPa s m2 m-3 to MPa s m2 mmol-1 - soilr1(j) = soilr1(j) * 1E-6_r8 * 18_r8 * 0.001_r8 - ! second component of below ground resistance is related to root hydraulics - soilr2(j) = EDecophyscon%rootresist(FT)/(rootmass(j)*dz(c,j)) - soilr_z(j) = soilr1(j)+soilr2(j) - end if - enddo +! if (SPA_soil) then ! normal case don't run this. +! rootr(p,:) = 0._r8 +! do FT = 1,numpft_ed + +! ! Soil Physics +! do j = 1,nlevgrnd +! ! CLM water retention curve. Clapp and Hornberger equation. +! s1 = max(h2osoi_vol(c,j)/watsat(c,j), 0.01_r8) +! s1 = min(1.0_r8,s1) +! smp_node = -sucsat(c,j)*s1**(-bsw(c,j)) +! swp_mpa(j) = smp_node *10.0_r8/1000000.0_r8 !convert from mm to Mpa + +! ! CLM hydraulic conductivity curve. +! ! As opposed to the Richard's equation solution in SoilHydrology.Mod +! ! the conductivity here is defined in the middle of the layer in question, not at the edge... +! xksat = 0.0070556_r8 * (10._r8**(-0.884_r8+0.0153_r8*sand(p)) ) +! hk(j) = xksat*s1**(2._r8*bsw(c,j)+2._r8) !removed the ice from here to avoid 1st ts crashing +! enddo + +! ! Root resistance +! rootxsecarea=3.14159*rootrad**2 +! do j = 1,nlevgrnd +! rootmass(j) = EDecophyscon%soilbeta(FT) * cpatch%rootfr_ft(FT,j) +! rootlength(j) = rootmass(j)/(rootdens*rootxsecarea) !m m-3 soil +! Lsoil(j) = hk(j)/1000/head !converts from mms-1 to ms-1 and then to m2 s-1 MPa-1 +! if(Lsoil(j) < 1e-35_r8.or.cpatch%rootfr_ft(ft,j) <= 0.0_r8)then !prevent floating point error +! soilr_z(j) = 1e35_r8 +! soilr2(j) = 1e35_r8 +! else +! ! Soil-to-root water uptake from Newman (1969). +! rs = sqrt (1._r8 / (rootlength(j) * pi)) +! soilr1(j) = log(rs/rootrad) / (2.0_r8 * pi * rootlength(j) * Lsoil(j) * dz(c,j)) +! ! convert from MPa s m2 m-3 to MPa s m2 mmol-1 +! soilr1(j) = soilr1(j) * 1E-6_r8 * 18_r8 * 0.001_r8 +! ! second component of below ground resistance is related to root hydraulics +! soilr2(j) = EDecophyscon%rootresist(FT)/(rootmass(j)*dz(c,j)) +! soilr_z(j) = soilr1(j)+soilr2(j) +! end if +! enddo ! Aggregate soil layers - totestevap=0._r8 - weighted_SWP=0._r8 - estevap=0._r8 - fraction_uptake=0._r8 - canopy_soil_resistance=0._r8 !Reset Counters - totmaxevap = 0._r8 +! totestevap=0._r8 +! weighted_SWP=0._r8 +! estevap=0._r8 +! fraction_uptake=0._r8 +! canopy_soil_resistance=0._r8 !Reset Counters +! totmaxevap = 0._r8 ! Estimated max transpiration from LWP gradient / soil resistance - do j = 1,nlevgrnd - estevap(j) = (swp_mpa(j) - minlwp)/(soilr_z(j)) - estevap(j) = max(0._r8,estevap(j)) ! no negative uptake - maxevap(j) = (0.0_r8 - minlwp)/(soilr2(j)) - enddo - totestevap = sum(estevap) - totmaxevap = sum(maxevap) +! do j = 1,nlevgrnd +! estevap(j) = (swp_mpa(j) - minlwp)/(soilr_z(j)) +! estevap(j) = max(0._r8,estevap(j)) ! no negative uptake +! maxevap(j) = (0.0_r8 - minlwp)/(soilr2(j)) +! enddo +! totestevap = sum(estevap) +! totmaxevap = sum(maxevap) ! Weighted soil water potential - do j = 1,nlevgrnd - if(totestevap > 0._r8)then - fraction_uptake(j) = estevap(j)/totestevap !Fraction of total ET taken from this soil layer - else - estevap(j) = 0._r8 - fraction_uptake(j)=1._r8/nlevgrnd - end if - weighted_SWP = weighted_SWP + swp_mpa(j) * estevap(j) - enddo - - - if(totestevap > 0._r8)then - weighted_swp = weighted_swp/totestevap - ! weight SWP for the total evaporation - else - write(iulog,*) 'empty soil', totestevap - ! error check - weighted_swp = minlwp - end if +! do j = 1,nlevgrnd +! if(totestevap > 0._r8)then +! fraction_uptake(j) = estevap(j)/totestevap !Fraction of total ET taken from this soil layer +! else +! estevap(j) = 0._r8 +! fraction_uptake(j)=1._r8/nlevgrnd +! end if +! weighted_SWP = weighted_SWP + swp_mpa(j) * estevap(j) +! enddo + +! if(totestevap > 0._r8)then +! weighted_swp = weighted_swp/totestevap +! ! weight SWP for the total evaporation +! else +! write(iulog,*) 'empty soil', totestevap +! ! error check +! weighted_swp = minlwp +! end if ! Weighted soil-root resistance. Aggregate the conductances (1/soilR) for each soil layer - do iv = 1,nv !leaf layers - fleaf = 1.0_r8/nv - do j = 1,nlevgrnd !root layers - ! Soil resistance for each canopy layer is related to leaf area - ! The conductance of the root system to the - ! whole canopy is reduced by the fraction of leaves in this layer... - canopy_soil_resistance(iv) = canopy_soil_resistance(iv)+fleaf * 1.0_r8/(soilr_z(j)) - enddo - ! Turn aggregated conductance back into resistance. mmol MPa-1 s-1 m-2 to MPa s m2 mmol-1 - canopy_soil_resistance(iv) = 1./canopy_soil_resistance(iv) - enddo - - currentPatch%btran_ft(FT) = totestevap/totmaxevap - do j = 1,nlevgrnd - if(sum(pftgs) > 0._r8)then !prevent problem with the first timestep - might fail - !bit-retart test as a result? FIX(RF,032414) - rootr(p,j) = rootr(p,j) + fraction_uptake(j) * pftgs(ft)/sum(pftgs) - else - rootr(p,j) = rootr(p,j) + fraction_uptake(j) * 1./numpft_ed - end if - enddo - - enddo !pft loop - - end if ! +! do iv = 1,nv !leaf layers +! fleaf = 1.0_r8/nv +! do j = 1,nlevgrnd !root layers +! ! Soil resistance for each canopy layer is related to leaf area +! ! The conductance of the root system to the +! ! whole canopy is reduced by the fraction of leaves in this layer... +! canopy_soil_resistance(iv) = canopy_soil_resistance(iv)+fleaf * 1.0_r8/(soilr_z(j)) +! enddo +! ! Turn aggregated conductance back into resistance. mmol MPa-1 s-1 m-2 to MPa s m2 mmol-1 +! canopy_soil_resistance(iv) = 1./canopy_soil_resistance(iv) +! enddo +! +! cpatch%btran_ft(FT) = totestevap/totmaxevap +! do j = 1,nlevgrnd +! if(sum(pftgs) > 0._r8)then !prevent problem with the first timestep - might fail +! !bit-retart test as a result? FIX(RF,032414) +! rootr(p,j) = rootr(p,j) + fraction_uptake(j) * pftgs(ft)/sum(pftgs) +! else +! rootr(p,j) = rootr(p,j) + fraction_uptake(j) * 1./numpft_ed +! end if +! enddo +! enddo !pft loop +! end if ! !--------------------------------------------------------------------------------------- ! end of SPA based recalculation of BTRAN and water uptake. !--------------------------------------------------------------------------------------- - !weight patch level output BTRAN for the - btran(p) = 0.0_r8 - do FT = 1,numpft_ed - if(sum(pftgs) > 0._r8)then !prevent problem with the first timestep - might fail - !bit-retart test as a result? FIX(RF,032414) - btran(p) = btran(p) + currentPatch%btran_ft(FT) * pftgs(ft)/sum(pftgs) - else - btran(p) = btran(p) + currentPatch%btran_ft(FT) * 1./numpft_ed - end if - enddo - - temprootr = sum(rootr(p,:)) - if(temprootr /= 1.0_r8)then - !write(iulog,*) 'error with rootr in canopy fluxes',sum(rootr(p,:)) - if(temprootr > 0._r8)then - do j = 1,nlevgrnd - rootr(p,j) = rootr(p,j) / temprootr - enddo - end if - end if - - else ! edpatch - currentPatch%btran_ft(1:numpft_ed) = 1._r8 - end if ! edpatch - - end associate - - end subroutine btran_ed + end module EDBtranMod diff --git a/biogeophys/EDSurfaceAlbedoMod.F90 b/biogeophys/EDSurfaceAlbedoMod.F90 index d9b84bd503..a8df5e600c 100644 --- a/biogeophys/EDSurfaceAlbedoMod.F90 +++ b/biogeophys/EDSurfaceAlbedoMod.F90 @@ -962,16 +962,16 @@ subroutine ED_SunShadeFracs(sites,nsites,bc_in,bc_out) use clm_varctl , only : iulog + implicit none + ! Arguments - type(ed_site_type),intent(inout),target :: sites(nsites) integer,intent(in) :: nsites type(bc_in_type),intent(in) :: bc_in(nsites) - type(bc_out_type),intent(out) :: bc_out(nsites) + type(bc_out_type),intent(inout) :: bc_out(nsites) ! locals - type (ed_patch_type),pointer :: cpatch ! c"urrent" patch real(r8) :: sunlai real(r8) :: shalai @@ -983,20 +983,24 @@ subroutine ED_SunShadeFracs(sites,nsites,bc_in,bc_out) do s = 1,nsites - + ifp = 0 cpatch => sites(s)%oldest_patch + do while (associated(cpatch)) ifp=ifp+1 - + + if( DEBUG ) write(iulog,*) 'edsurfRad_5600',ifp,s,cpatch%NCL_p,numpft_ed + ! zero out various datas cpatch%ed_parsun_z(:,:,:) = 0._r8 cpatch%ed_parsha_z(:,:,:) = 0._r8 cpatch%ed_laisun_z(:,:,:) = 0._r8 cpatch%ed_laisha_z(:,:,:) = 0._r8 - + bc_out(s)%fsun_pa(ifp) = 0._r8 + sunlai = 0._r8 shalai = 0._r8 @@ -1009,6 +1013,9 @@ subroutine ED_SunShadeFracs(sites,nsites,bc_in,bc_out) do CL = 1, cpatch%NCL_p do FT = 1,numpft_ed + + if( DEBUG ) write(iulog,*) 'edsurfRad_5601',CL,FT,cpatch%nrad(CL,ft) + do iv = 1, cpatch%nrad(CL,ft) !NORMAL CASE. ! FIX(SPM,040114) - existing comment diff --git a/main/EDTypesMod.F90 b/main/EDTypesMod.F90 index e167744af9..5dd71ff2cd 100755 --- a/main/EDTypesMod.F90 +++ b/main/EDTypesMod.F90 @@ -106,6 +106,28 @@ module EDTypesMod integer , allocatable :: pft_levscpf_ed(:) integer , allocatable :: scls_levscpf_ed(:) + + + type, private :: ctrl_parms_type + + + ! These parameters are dictated by FATES internals + + + + + ! These parameters are dictated by the host model or driver + + integer :: numSWBands ! Maximum number of broad-bands in the short-wave radiation + ! specturm to track + ! (typically 2 as a default, VIS/NIR, in ED variants <2016) + + integer :: numlevgrnd ! Number of soil layers + + end type ctrl_parms_type + + type(ctrl_parms_type), public :: ctrl_parms + !************************************ !** COHORT type structure ** diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 200271e876..ef10b388df 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -16,26 +16,12 @@ module FatesInterfaceMod ! ------------------------------------------------------------------------------------ use EDtypesMod , only : ed_site_type, & - numPatchesPerCol + numPatchesPerCol, & + ctrl_parms + use shr_kind_mod , only : r8 => shr_kind_r8 ! INTERF-TODO: REMOVE THIS use clm_varpar , only : nlevgrnd - ! ------------------------------------------------------------------------------------ - ! Certain dimension information used by FATES is dictated by the the driver - ! or the host model, or perhaps may be some compromise between what FATES will want - ! in a best-case scenario and what space the driver/host will allow based on its - ! memory constraints (most-likely due to IO) - ! ------------------------------------------------------------------------------------ - - type, private :: fates_dims_type - - integer :: numSWBands ! Maximum number of broad-bands in the short-wave radiation - ! specturm to track - ! (typically 2 as a default, VIS/NIR, in ED variants <2016) - - end type fates_dims_type - - ! ------------------------------------------------------------------------------------ ! Notes on types @@ -67,6 +53,11 @@ module FatesInterfaceMod ! volumetric soil water at saturation (porosity) real(r8), allocatable :: watsat_sl(:) + ! If there is no liquid volume of water in a soil layer, or + ! if the layer is 2 degrees below freezing, the layer will not + ! be deemed active for water uptake via transpiration and photosynthesis + logical, allocatable :: active_uptake_sl(:) + end type bc_in_type @@ -77,7 +68,8 @@ module FatesInterfaceMod ! Root soil water stress (resistance) by layer ! (diagnostic, should not be used by HLM) - real(r8),allocatable :: rresis_pa(:,:) +! real(r8),allocatable :: rresis_pa(:,:) ! not used by host, not calculated + ! yet by FATES ! Effective fraction of roots in each soil layer ! (diagnostic, should not be used by HLM) @@ -117,19 +109,12 @@ module FatesInterfaceMod contains - procedure, public :: allocate_bcs procedure, public :: zero_bcs end type fates_interface_type - ! ------------------------------------------------------------------------------------ - ! Dimension information is independent of which clump it is on - ! and since these are typically read only variables and never updated on the clump - ! we need not attach these to the threaded instances of the fates_interface_type - ! ------------------------------------------------------------------------------------ - - type(fates_dims_type) :: fates_dims - + + public :: set_fates_ctrlparms contains @@ -158,40 +143,49 @@ end subroutine fates_clean ! ==================================================================================== - subroutine allocate_bcs(this,s) + subroutine allocate_bcin(bc_in) ! --------------------------------------------------------------------------------- ! Allocate and Initialze the FATES boundary condition vectors ! --------------------------------------------------------------------------------- - + implicit none - class(fates_interface_type), intent(inout) :: this - integer, intent(in) :: s + type(bc_in_type), intent(inout) :: bc_in ! Allocate input boundaries ! Radiation - allocate(this%bc_in(s)%solad_pa(numPatchesPerCol,fates_dims%numSWBands)) - allocate(this%bc_in(s)%solai_pa(numPatchesPerCol,fates_dims%numSWBands)) - + allocate(bc_in%solad_pa(numPatchesPerCol,ctrl_parms%numSWBands)) + allocate(bc_in%solai_pa(numPatchesPerCol,ctrl_parms%numSWBands)) + ! Hydrology - allocate(this%bc_in(s)%smp_sl(nlevgrnd)) - allocate(this%bc_in(s)%eff_porosity_sl(nlevgrnd)) - allocate(this%bc_in(s)%watsat_sl(nlevgrnd)) + allocate(bc_in%smp_sl(ctrl_parms%numlevgrnd)) + allocate(bc_in%eff_porosity_sl(ctrl_parms%numlevgrnd)) + allocate(bc_in%watsat_sl(ctrl_parms%numlevgrnd)) + allocate(bc_in%active_uptake_sl(ctrl_parms%numlevgrnd)) - ! Allocate output boundaries + return + end subroutine allocate_bcin + + subroutine allocate_bcout(bc_out) + ! --------------------------------------------------------------------------------- + ! Allocate and Initialze the FATES boundary condition vectors + ! --------------------------------------------------------------------------------- + + implicit none + type(bc_out_type), intent(inout) :: bc_out + + ! Radiation - allocate(this%bc_out(s)%fsun_pa(numPatchesPerCol)) + allocate(bc_out%fsun_pa(numPatchesPerCol)) ! Hydrology - allocate(this%bc_out(s)%rresis_pa(numPatchesPerCol,nlevgrnd)) - allocate(this%bc_out(s)%rootr_pa(numPatchesPerCol,nlevgrnd)) - allocate(this%bc_out(s)%btran_pa(numPatchesPerCol)) - - + allocate(bc_out%rootr_pa(numPatchesPerCol,ctrl_parms%numlevgrnd)) + allocate(bc_out%btran_pa(numPatchesPerCol)) + return - end subroutine allocate_bcs + end subroutine allocate_bcout ! ==================================================================================== @@ -203,21 +197,110 @@ subroutine zero_bcs(this,s) ! Input boundaries - this%bc_in(s)%solad_pa(:,:) = 0.0_r8 - this%bc_in(s)%solai_pa(:,:) = 0.0_r8 - this%bc_in(s)%smp_sl(:) = 0.0_r8 - this%bc_in(s)%eff_porosity_sl(:) = 0.0_r8 - this%bc_in(s)%watsat_sl(:) = 0.0_r8 + this%bc_in(s)%solad_pa(:,:) = 0.0_r8 + this%bc_in(s)%solai_pa(:,:) = 0.0_r8 + this%bc_in(s)%smp_sl(:) = 0.0_r8 + this%bc_in(s)%eff_porosity_sl(:) = 0.0_r8 + this%bc_in(s)%watsat_sl(:) = 0.0_r8 + this%bc_in(s)%active_uptake_sl(:) = .false. ! Output boundaries - this%bc_out(s)%fsun_pa(:) = 0.0_r8 - this%bc_out(s)%rresis_pa(:) = 0.0_r8 - this%bc_out(s)%rootr_pa(:) = 0.0_r8 - this%bc_out(s)%btran_pa(:) = 0.0_r8 + this%bc_out(s)%fsun_pa(:) = 0.0_r8 + this%bc_out(s)%rootr_pa(:,:) = 0.0_r8 + this%bc_out(s)%btran_pa(:) = 0.0_r8 return end subroutine zero_bcs + + ! ==================================================================================== + + subroutine set_fates_ctrlparms(tag,dimval) + + ! --------------------------------------------------------------------------------- + ! INTERF-TODO: NEED ALLOWANCES FOR REAL AND CHARACTER ARGS.. + ! Certain model control parameters and dimensions used by FATES are dictated by + ! the the driver or the host mode. To see which parameters should be filled here + ! please also look at the ctrl_parms_type in FATESTYpeMod, in the section listing + ! components dictated by the host model. + ! + ! Some important points: + ! 1. Calls to this function are likely from the clm_fates module in the HLM. + ! 2. The calls should be preceeded by a flush function. + ! 3. All values in ctrl_parm (FATESTypesMod.F90) that are classified as + ! 'dictated by the HLM' must be listed in this subroutine + ! 4. Should look like this: + ! + ! call set_fates_ctrlparms('flush_to_unset') + ! call set_fates_ctrlparms('num_sw_bbands',numrad) ! or other variable + ! ... + ! call set_fates_ctrlparms('num_lev_ground',nlevgrnd) ! or other variable + ! call set_fates_ctrlparms('check_allset') + ! + ! RGK-2016 + ! --------------------------------------------------------------------------------- + + ! Arguments + integer, optional, intent(in) :: dimval + character(len=*),intent(in) :: tag + + ! local variables + logical :: all_set + integer, parameter :: unset_int = -999 + real(r8), parameter :: unset_double = -999.9 + + select case (trim(tag)) + case('flush_to_unset') + + write(*,*) 'Flushing FATES control parameters prior to transfer from host' + ctrl_parms%numSwBands = unset_int + ctrl_parms%numlevgrnd = unset_int + + case('check_allset') + + if(ctrl_parms%numSWBands .eq. unset_int) then + write(*,*) 'FATES dimension/parameter unset: num_sw_rad_bbands' + ! INTERF-TODO: FATES NEEDS INTERNAL end_run + ! end_run('MESSAGE') + end if + + if(ctrl_parms%numlevgrnd .eq. unset_int) then + write(*,*) 'FATES dimension/parameter unset: numlevground' + ! INTERF-TODO: FATES NEEDS INTERNAL end_run + ! end_run('MESSAGE') + end if + + write(*,*) 'Checked. All control parameters sent to FATES.' + + case default + + if(present(dimval))then + select case (trim(tag)) + + case('num_sw_bbands') + + ctrl_parms%numSwBands = dimval + write(*,*) 'Transfering num_sw_bbands = ',dimval,' to FATES' + + case('num_lev_ground') + + ctrl_parms%numlevgrnd = dimval + write(*,*) 'Transfering num_lev_ground = ',dimval,' to FATES' + + case default + write(*,*) 'tag not recognized:',trim(tag) + ! end_run + end select + else + write(*,*) 'no value was provided for the tag' + end if + + end select + + + return + end subroutine set_fates_ctrlparms + end module FatesInterfaceMod From 3ee5771a3225b49250689ebfe5d9a1f90242c2a8 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Fri, 1 Jul 2016 09:38:47 -0700 Subject: [PATCH 5/9] 1) Changed naming convention of some boundary condition variables per suggestion by D Lawrence. 2) moved a logical step related to btran to be inside the fates btran calculation, to do this I moved the soil volume and temperature variables to boundary conditions and removed the active flag logical --- biogeophys/EDBtranMod.F90 | 29 ++++++------- biogeophys/EDSurfaceAlbedoMod.F90 | 12 +++--- main/FatesInterfaceMod.F90 | 70 ++++++++++++++++--------------- 3 files changed, 57 insertions(+), 54 deletions(-) diff --git a/biogeophys/EDBtranMod.F90 b/biogeophys/EDBtranMod.F90 index d9f5b464f8..c0fc917956 100644 --- a/biogeophys/EDBtranMod.F90 +++ b/biogeophys/EDBtranMod.F90 @@ -6,6 +6,7 @@ module EDBtranMod ! ------------------------------------------------------------------------------------ use pftconMod , only : pftcon + use clm_varcon , only : tfrz use EDTypesMod , only : ed_site_type, & ed_patch_type, & ed_cohort_type, & @@ -34,11 +35,11 @@ subroutine btran_ed( sites, nsites, bc_in, bc_out) ! --------------------------------------------------------------------------------- ! Calculate the transpiration wetness function (BTRAN) and the root uptake ! distribution (ROOTR). - ! Boundary conditions in: bc_in(s)%eff_porosity_sl(j) unfrozen porosity - ! bc_in(s)%watsat_sl(j) porosity - ! bc_in(s)%active_uptake_sl(j) frozen/not frozen - ! bc_in(s)%smp_sl(j) suction - ! Boundary conditions out: bc_out(s)%rootr_pa root uptake distribution + ! Boundary conditions in: bc_in(s)%eff_porosity_gl(j) unfrozen porosity + ! bc_in(s)%watsat_gl(j) porosity + ! bc_in(s)%active_uptake_gl(j) frozen/not frozen + ! bc_in(s)%smp_gl(j) suction + ! Boundary conditions out: bc_out(s)%rootr_pagl root uptake distribution ! bc_out(s)%btran_pa wetness factor ! --------------------------------------------------------------------------------- @@ -86,12 +87,12 @@ subroutine btran_ed( sites, nsites, bc_in, bc_out) ! Calculations are only relevant where liquid water exists ! see clm_fates%wrap_btran for calculation with CLM/ALM - - if( bc_in(s)%active_uptake_sl(j) ) then + + if ( bc_in(s)%h2o_liqvol_gl(j) .gt. 0._r8 .and. bc_in(s)%tempk_gl(j) .gt. tfrz-2._r8) then - smp_node = max(smpsc(ft), bc_in(s)%smp_sl(j)) + smp_node = max(smpsc(ft), bc_in(s)%smp_gl(j)) - rresis = min( (bc_in(s)%eff_porosity_sl(j)/bc_in(s)%watsat_sl(j))* & + rresis = min( (bc_in(s)%eff_porosity_gl(j)/bc_in(s)%watsat_gl(j))* & (smp_node - smpsc(ft)) / (smpso(ft) - smpsc(ft)), 1._r8) cpatch%rootr_ft(ft,j) = cpatch%rootfr_ft(ft,j)*rresis @@ -134,14 +135,14 @@ subroutine btran_ed( sites, nsites, bc_in, bc_out) ! distributed over the soil layers. do j = 1,numlevgrnd - bc_out(s)%rootr_pa(ifp,j) = 0._r8 + bc_out(s)%rootr_pagl(ifp,j) = 0._r8 do ft = 1,numpft_ed if(sum(pftgs) > 0._r8)then !prevent problem with the first timestep - might fail !bit-retart test as a result? FIX(RF,032414) - bc_out(s)%rootr_pa(ifp,j) = bc_out(s)%rootr_pa(ifp,j) + & + bc_out(s)%rootr_pagl(ifp,j) = bc_out(s)%rootr_pagl(ifp,j) + & cpatch%rootr_ft(ft,j) * pftgs(ft)/sum(pftgs) else - bc_out(s)%rootr_pa(ifp,j) = bc_out(s)%rootr_pa(ifp,j) + & + bc_out(s)%rootr_pagl(ifp,j) = bc_out(s)%rootr_pagl(ifp,j) + & cpatch%rootr_ft(ft,j) * 1./numpft_ed end if enddo @@ -160,12 +161,12 @@ subroutine btran_ed( sites, nsites, bc_in, bc_out) - temprootr = sum(bc_out(s)%rootr_pa(ifp,:)) + temprootr = sum(bc_out(s)%rootr_pagl(ifp,:)) if(temprootr /= 1.0_r8)then write(iulog,*) 'error with rootr in canopy fluxes',temprootr if(temprootr > 0._r8)then do j = 1,numlevgrnd - bc_out(s)%rootr_pa(ifp,j) = bc_out(s)%rootr_pa(ifp,j)/temprootr + bc_out(s)%rootr_pagl(ifp,j) = bc_out(s)%rootr_pagl(ifp,j)/temprootr enddo end if end if diff --git a/biogeophys/EDSurfaceAlbedoMod.F90 b/biogeophys/EDSurfaceAlbedoMod.F90 index a8df5e600c..7aedb32eb3 100644 --- a/biogeophys/EDSurfaceAlbedoMod.F90 +++ b/biogeophys/EDSurfaceAlbedoMod.F90 @@ -1065,21 +1065,21 @@ subroutine ED_SunShadeFracs(sites,nsites,bc_in,bc_out) if ( DEBUG ) then write(iulog,*) 'edsurfRad 653 ', cpatch%ed_parsun_z(CL,ft,iv) - write(iulog,*) 'edsurfRad 654 ', bc_in(s)%solad_pa(ifp,ipar) - write(iulog,*) 'edsurfRad 655 ', bc_in(s)%solai_pa(ifp,ipar) + write(iulog,*) 'edsurfRad 654 ', bc_in(s)%solad_parb(ifp,ipar) + write(iulog,*) 'edsurfRad 655 ', bc_in(s)%solai_parb(ifp,ipar) write(iulog,*) 'edsurfRad 656 ', cpatch%fabd_sun_z(CL,ft,iv) write(iulog,*) 'edsurfRad 657 ', cpatch%fabi_sun_z(CL,ft,iv) endif cpatch%ed_parsun_z(CL,ft,iv) = & - bc_in(s)%solad_pa(ifp,ipar)*cpatch%fabd_sun_z(CL,ft,iv) + & - bc_in(s)%solai_pa(ifp,ipar)*cpatch%fabi_sun_z(CL,ft,iv) + bc_in(s)%solad_parb(ifp,ipar)*cpatch%fabd_sun_z(CL,ft,iv) + & + bc_in(s)%solai_parb(ifp,ipar)*cpatch%fabi_sun_z(CL,ft,iv) if ( DEBUG )write(iulog,*) 'edsurfRad 663 ', cpatch%ed_parsun_z(CL,ft,iv) cpatch%ed_parsha_z(CL,ft,iv) = & - bc_in(s)%solad_pa(ifp,ipar)*cpatch%fabd_sha_z(CL,ft,iv) + & - bc_in(s)%solai_pa(ifp,ipar)*cpatch%fabi_sha_z(CL,ft,iv) + bc_in(s)%solad_parb(ifp,ipar)*cpatch%fabd_sha_z(CL,ft,iv) + & + bc_in(s)%solai_parb(ifp,ipar)*cpatch%fabi_sha_z(CL,ft,iv) if ( DEBUG ) write(iulog,*) 'edsurfRad 669 ', cpatch%ed_parsha_z(CL,ft,iv) diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index ef10b388df..1394145a97 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -31,6 +31,9 @@ module FatesInterfaceMod ! or alias other variables, ALLOCATABLES cannnot. According to S. Lionel ! (Intel-Forum Post), ALLOCATABLES are better perfomance wise as long as they point ! to contiguous memory spaces and do not alias other variables, the case here. + ! Naming conventions: _gl means ground layer dimensions + ! _pa means patch dimensions + ! _rb means radiation band ! ------------------------------------------------------------------------------------ type, public :: bc_in_type @@ -38,25 +41,27 @@ module FatesInterfaceMod ! The actual number of FATES' ED patches integer :: npatches - ! Downwelling direct beam radiation (patch,broad-band) [W/m2] - real(r8), allocatable :: solad_pa(:,:) + ! Downwelling direct beam radiation (patch,radiation-band) [W/m2] + real(r8), allocatable :: solad_parb(:,:) - ! Downwelling diffuse (I-ndirect) radiation (patch,broad-band) [W/m2] - real(r8), allocatable :: solai_pa(:,:) + ! Downwelling diffuse (I-ndirect) radiation (patch,radiation-band) [W/m2] + real(r8), allocatable :: solai_parb(:,:) ! Soil suction potential of layers in each site, negative, [mm] - real(r8), allocatable :: smp_sl(:) + real(r8), allocatable :: smp_gl(:) ! Effective porosity = porosity - vol_ic, of layers in each site [-] - real(r8), allocatable :: eff_porosity_sl(:) + real(r8), allocatable :: eff_porosity_gl(:) ! volumetric soil water at saturation (porosity) - real(r8), allocatable :: watsat_sl(:) + real(r8), allocatable :: watsat_gl(:) + + ! Temperature of ground layers [K] + real(r8), allocatable :: tempk_gl(:) + + ! Liquid volume in ground layer + real(r8), allocatable :: h2o_liqvol_gl(:) - ! If there is no liquid volume of water in a soil layer, or - ! if the layer is 2 degrees below freezing, the layer will not - ! be deemed active for water uptake via transpiration and photosynthesis - logical, allocatable :: active_uptake_sl(:) end type bc_in_type @@ -66,14 +71,8 @@ module FatesInterfaceMod ! Sunlit fraction of the canopy for this patch [0-1] real(r8),allocatable :: fsun_pa(:) - ! Root soil water stress (resistance) by layer - ! (diagnostic, should not be used by HLM) -! real(r8),allocatable :: rresis_pa(:,:) ! not used by host, not calculated - ! yet by FATES - ! Effective fraction of roots in each soil layer - ! (diagnostic, should not be used by HLM) - real(r8), allocatable :: rootr_pa(:,:) + real(r8), allocatable :: rootr_pagl(:,:) ! Integrated (vertically) transpiration wetness factor (0 to 1) ! (diagnostic, should not be used by HLM) @@ -155,14 +154,15 @@ subroutine allocate_bcin(bc_in) ! Allocate input boundaries ! Radiation - allocate(bc_in%solad_pa(numPatchesPerCol,ctrl_parms%numSWBands)) - allocate(bc_in%solai_pa(numPatchesPerCol,ctrl_parms%numSWBands)) + allocate(bc_in%solad_parb(numPatchesPerCol,ctrl_parms%numSWBands)) + allocate(bc_in%solai_parb(numPatchesPerCol,ctrl_parms%numSWBands)) ! Hydrology - allocate(bc_in%smp_sl(ctrl_parms%numlevgrnd)) - allocate(bc_in%eff_porosity_sl(ctrl_parms%numlevgrnd)) - allocate(bc_in%watsat_sl(ctrl_parms%numlevgrnd)) - allocate(bc_in%active_uptake_sl(ctrl_parms%numlevgrnd)) + allocate(bc_in%smp_gl(ctrl_parms%numlevgrnd)) + allocate(bc_in%eff_porosity_gl(ctrl_parms%numlevgrnd)) + allocate(bc_in%watsat_gl(ctrl_parms%numlevgrnd)) + allocate(bc_in%tempk_gl(ctrl_parms%numlevgrnd)) + allocate(bc_in%h2o_liqvol_gl(ctrl_parms%numlevgrnd)) return end subroutine allocate_bcin @@ -181,9 +181,10 @@ subroutine allocate_bcout(bc_out) allocate(bc_out%fsun_pa(numPatchesPerCol)) ! Hydrology - allocate(bc_out%rootr_pa(numPatchesPerCol,ctrl_parms%numlevgrnd)) + allocate(bc_out%rootr_pagl(numPatchesPerCol,ctrl_parms%numlevgrnd)) allocate(bc_out%btran_pa(numPatchesPerCol)) + return end subroutine allocate_bcout @@ -197,18 +198,19 @@ subroutine zero_bcs(this,s) ! Input boundaries - this%bc_in(s)%solad_pa(:,:) = 0.0_r8 - this%bc_in(s)%solai_pa(:,:) = 0.0_r8 - this%bc_in(s)%smp_sl(:) = 0.0_r8 - this%bc_in(s)%eff_porosity_sl(:) = 0.0_r8 - this%bc_in(s)%watsat_sl(:) = 0.0_r8 - this%bc_in(s)%active_uptake_sl(:) = .false. + this%bc_in(s)%solad_parb(:,:) = 0.0_r8 + this%bc_in(s)%solai_parb(:,:) = 0.0_r8 + this%bc_in(s)%smp_gl(:) = 0.0_r8 + this%bc_in(s)%eff_porosity_gl(:) = 0.0_r8 + this%bc_in(s)%watsat_gl(:) = 0.0_r8 + this%bc_in(s)%tempk_gl(:) = 0.0_r8 + this%bc_in(s)%h2o_liqvol_gl(:) = 0.0_r8 ! Output boundaries - this%bc_out(s)%fsun_pa(:) = 0.0_r8 - this%bc_out(s)%rootr_pa(:,:) = 0.0_r8 - this%bc_out(s)%btran_pa(:) = 0.0_r8 + this%bc_out(s)%fsun_pa(:) = 0.0_r8 + this%bc_out(s)%rootr_pagl(:,:) = 0.0_r8 + this%bc_out(s)%btran_pa(:) = 0.0_r8 return end subroutine zero_bcs From 34822f97fdef3d29438348bded03a748a904b2bb Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Sun, 3 Jul 2016 18:10:03 -0400 Subject: [PATCH 6/9] Added a call to determine active layers with plant water uptake. Could not embed this logic inside the btran because the active layer filters are also needed for determining which layers suction can be calculated for, but this is outside of FATES perview, so a call before btran was made. --- biogeophys/EDBtranMod.F90 | 66 +++++++++++++++++++++++++++++++------- main/FatesInterfaceMod.F90 | 8 ++++- 2 files changed, 61 insertions(+), 13 deletions(-) diff --git a/biogeophys/EDBtranMod.F90 b/biogeophys/EDBtranMod.F90 index c0fc917956..b227055950 100644 --- a/biogeophys/EDBtranMod.F90 +++ b/biogeophys/EDBtranMod.F90 @@ -22,15 +22,58 @@ module EDBtranMod private public :: btran_ed - ! - type(ed_cohort_type), pointer :: currentCohort ! current cohort - type(ed_patch_type) , pointer :: currentPatch ! current patch + public :: get_active_suction_layers contains - ! ==================================================================================== + ! ==================================================================================== + + logical function check_layer_water(h2o_liq_vol, tempk) + + implicit none + ! Arguments + real(r8) :: h2o_liq_vol + real(r8) :: tempk + + if ( h2o_liq_vol .gt. 0._r8 .and. tempk .gt. tfrz-2._r8) then + check_layer_water = .true. + else + check_layer_water = .false. + end if + return + end function check_layer_water + + ! ===================================================================================== + + subroutine get_active_suction_layers(sites,nsites,bc_in,bc_out) + + ! Arguments + + type(ed_site_type),intent(inout),target :: sites(nsites) + integer,intent(in) :: nsites + type(bc_in_type),intent(in) :: bc_in(nsites) + type(bc_out_type),intent(inout) :: bc_out(nsites) + + ! !LOCAL VARIABLES: + integer :: s ! site + integer :: j ! soil layer + !------------------------------------------------------------------------------ + + associate( & + numlevgrnd => ctrl_parms%numlevgrnd ) + + do s = 1,nsites + do j = 1,numlevgrnd + bc_out(s)%active_suction_gl(j) = check_layer_water(bc_in(s)%h2o_liqvol_gl(j),bc_in(s)%tempk_gl(j) ) + end do + end do + + end associate + end subroutine get_active_suction_layers + + ! ===================================================================================== - subroutine btran_ed( sites, nsites, bc_in, bc_out) + subroutine btran_ed( sites, nsites, bc_in, bc_out) ! --------------------------------------------------------------------------------- ! Calculate the transpiration wetness function (BTRAN) and the root uptake @@ -42,8 +85,6 @@ subroutine btran_ed( sites, nsites, bc_in, bc_out) ! Boundary conditions out: bc_out(s)%rootr_pagl root uptake distribution ! bc_out(s)%btran_pa wetness factor ! --------------------------------------------------------------------------------- - - ! Arguments @@ -55,6 +96,7 @@ subroutine btran_ed( sites, nsites, bc_in, bc_out) ! ! !LOCAL VARIABLES: type(ed_patch_type),pointer :: cpatch ! Current Patch Pointer + type(ed_cohort_type),pointer :: ccohort ! Current cohort pointer integer :: s ! site integer :: j ! soil layer integer :: ifp ! patch vector index for the site @@ -88,7 +130,7 @@ subroutine btran_ed( sites, nsites, bc_in, bc_out) ! Calculations are only relevant where liquid water exists ! see clm_fates%wrap_btran for calculation with CLM/ALM - if ( bc_in(s)%h2o_liqvol_gl(j) .gt. 0._r8 .and. bc_in(s)%tempk_gl(j) .gt. tfrz-2._r8) then + if ( check_layer_water(bc_in(s)%h2o_liqvol_gl(j),bc_in(s)%tempk_gl(j)) ) then smp_node = max(smpsc(ft), bc_in(s)%smp_gl(j)) @@ -123,10 +165,10 @@ subroutine btran_ed( sites, nsites, bc_in, bc_out) ! PFT-averaged point level root fraction for extraction purposese. ! This probably needs to be weighted by actual transpiration from each pft. FIX(RF,032414). pftgs(:) = 0._r8 - currentCohort => cpatch%tallest - do while(associated(currentCohort)) - pftgs(currentCohort%pft) = pftgs(currentCohort%pft) + currentCohort%gscan * currentCohort%n - currentCohort => currentCohort%shorter + ccohort => cpatch%tallest + do while(associated(ccohort)) + pftgs(ccohort%pft) = pftgs(ccohort%pft) + ccohort%gscan * ccohort%n + ccohort => ccohort%shorter enddo ! Process the boundary output, this is necessary for calculating the soil-moisture diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 1394145a97..7266f1d889 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -71,6 +71,11 @@ module FatesInterfaceMod ! Sunlit fraction of the canopy for this patch [0-1] real(r8),allocatable :: fsun_pa(:) + ! Logical stating whether a ground layer can have water uptake by plants + ! The only condition right now is that liquid water exists + ! The name (suction) is used to indicate that soil suction should be calculated + logical, allocatable :: active_suction_gl(:) + ! Effective fraction of roots in each soil layer real(r8), allocatable :: rootr_pagl(:,:) @@ -181,6 +186,7 @@ subroutine allocate_bcout(bc_out) allocate(bc_out%fsun_pa(numPatchesPerCol)) ! Hydrology + allocate(bc_out%active_suction_gl(ctrl_parms%numlevgrnd)) allocate(bc_out%rootr_pagl(numPatchesPerCol,ctrl_parms%numlevgrnd)) allocate(bc_out%btran_pa(numPatchesPerCol)) @@ -207,7 +213,7 @@ subroutine zero_bcs(this,s) this%bc_in(s)%h2o_liqvol_gl(:) = 0.0_r8 ! Output boundaries - + this%bc_out(s)%active_suction_gl(:) = .false. this%bc_out(s)%fsun_pa(:) = 0.0_r8 this%bc_out(s)%rootr_pagl(:,:) = 0.0_r8 this%bc_out(s)%btran_pa(:) = 0.0_r8 From 36fbb223e01d1967355b097f04d4d404e3c2f660 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Sun, 3 Jul 2016 20:12:05 -0400 Subject: [PATCH 7/9] optimized the logical function that tests if liquid soil water is present for uptake. --- biogeophys/EDBtranMod.F90 | 17 +++++++++-------- 1 file changed, 9 insertions(+), 8 deletions(-) diff --git a/biogeophys/EDBtranMod.F90 b/biogeophys/EDBtranMod.F90 index b227055950..a35374498b 100644 --- a/biogeophys/EDBtranMod.F90 +++ b/biogeophys/EDBtranMod.F90 @@ -32,13 +32,14 @@ logical function check_layer_water(h2o_liq_vol, tempk) implicit none ! Arguments - real(r8) :: h2o_liq_vol - real(r8) :: tempk - - if ( h2o_liq_vol .gt. 0._r8 .and. tempk .gt. tfrz-2._r8) then - check_layer_water = .true. - else - check_layer_water = .false. + real(r8),intent(in) :: h2o_liq_vol + real(r8),intent(in) :: tempk + + check_layer_water = .false. + if ( h2o_liq_vol .gt. 0._r8 ) then + if ( tempk .gt. tfrz-2._r8) then + check_layer_water = .true. + end if end if return end function check_layer_water @@ -64,7 +65,7 @@ subroutine get_active_suction_layers(sites,nsites,bc_in,bc_out) do s = 1,nsites do j = 1,numlevgrnd - bc_out(s)%active_suction_gl(j) = check_layer_water(bc_in(s)%h2o_liqvol_gl(j),bc_in(s)%tempk_gl(j) ) + bc_out(s)%active_suction_gl(j) = check_layer_water( bc_in(s)%h2o_liqvol_gl(j),bc_in(s)%tempk_gl(j) ) end do end do From 678e3df4834f81eb2d5f1836fda21ea10d6cf83f Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Mon, 18 Jul 2016 14:56:27 -0400 Subject: [PATCH 8/9] Added an input filter to the wrap_btran that operates on the exposed veg filter. --- biogeophys/EDBtranMod.F90 | 33 ++++++++++++++++++++------------- main/FatesInterfaceMod.F90 | 3 +++ 2 files changed, 23 insertions(+), 13 deletions(-) diff --git a/biogeophys/EDBtranMod.F90 b/biogeophys/EDBtranMod.F90 index a35374498b..ebb41293cd 100644 --- a/biogeophys/EDBtranMod.F90 +++ b/biogeophys/EDBtranMod.F90 @@ -36,6 +36,7 @@ logical function check_layer_water(h2o_liq_vol, tempk) real(r8),intent(in) :: tempk check_layer_water = .false. + if ( h2o_liq_vol .gt. 0._r8 ) then if ( tempk .gt. tfrz-2._r8) then check_layer_water = .true. @@ -64,9 +65,13 @@ subroutine get_active_suction_layers(sites,nsites,bc_in,bc_out) numlevgrnd => ctrl_parms%numlevgrnd ) do s = 1,nsites - do j = 1,numlevgrnd - bc_out(s)%active_suction_gl(j) = check_layer_water( bc_in(s)%h2o_liqvol_gl(j),bc_in(s)%tempk_gl(j) ) - end do + if (bc_in(s)%filter_btran) then + do j = 1,numlevgrnd + bc_out(s)%active_suction_gl(j) = check_layer_water( bc_in(s)%h2o_liqvol_gl(j),bc_in(s)%tempk_gl(j) ) + end do + else + bc_out(s)%active_suction_gl(:) = .false. + end if end do end associate @@ -116,7 +121,9 @@ subroutine btran_ed( sites, nsites, bc_in, bc_out) ) do s = 1,nsites - + + bc_out(s)%rootr_pagl(:,:) = 0._r8 + ifp = 0 cpatch => sites(s)%oldest_patch do while (associated(cpatch)) @@ -130,13 +137,13 @@ subroutine btran_ed( sites, nsites, bc_in, bc_out) ! Calculations are only relevant where liquid water exists ! see clm_fates%wrap_btran for calculation with CLM/ALM - + if ( check_layer_water(bc_in(s)%h2o_liqvol_gl(j),bc_in(s)%tempk_gl(j)) ) then smp_node = max(smpsc(ft), bc_in(s)%smp_gl(j)) rresis = min( (bc_in(s)%eff_porosity_gl(j)/bc_in(s)%watsat_gl(j))* & - (smp_node - smpsc(ft)) / (smpso(ft) - smpsc(ft)), 1._r8) + (smp_node - smpsc(ft)) / (smpso(ft) - smpsc(ft)), 1._r8) cpatch%rootr_ft(ft,j) = cpatch%rootfr_ft(ft,j)*rresis @@ -183,10 +190,10 @@ subroutine btran_ed( sites, nsites, bc_in, bc_out) if(sum(pftgs) > 0._r8)then !prevent problem with the first timestep - might fail !bit-retart test as a result? FIX(RF,032414) bc_out(s)%rootr_pagl(ifp,j) = bc_out(s)%rootr_pagl(ifp,j) + & - cpatch%rootr_ft(ft,j) * pftgs(ft)/sum(pftgs) + cpatch%rootr_ft(ft,j) * pftgs(ft)/sum(pftgs) else bc_out(s)%rootr_pagl(ifp,j) = bc_out(s)%rootr_pagl(ifp,j) + & - cpatch%rootr_ft(ft,j) * 1./numpft_ed + cpatch%rootr_ft(ft,j) * 1./numpft_ed end if enddo enddo @@ -202,11 +209,9 @@ subroutine btran_ed( sites, nsites, bc_in, bc_out) end if enddo - - temprootr = sum(bc_out(s)%rootr_pagl(ifp,:)) - if(temprootr /= 1.0_r8)then - write(iulog,*) 'error with rootr in canopy fluxes',temprootr + if(abs(1.0_r8-temprootr) > 1.0e-9_r8)then + write(iulog,*) 'error with rootr in canopy fluxes',temprootr,sum(pftgs),sum(cpatch%rootr_ft(1:2,:),dim=2) if(temprootr > 0._r8)then do j = 1,numlevgrnd bc_out(s)%rootr_pagl(ifp,j) = bc_out(s)%rootr_pagl(ifp,j)/temprootr @@ -216,11 +221,13 @@ subroutine btran_ed( sites, nsites, bc_in, bc_out) cpatch => cpatch%younger end do + + end do end associate - end subroutine btran_ed + end subroutine btran_ed ! ========================================================================================= diff --git a/main/FatesInterfaceMod.F90 b/main/FatesInterfaceMod.F90 index 7266f1d889..79ae9d6009 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -62,6 +62,9 @@ module FatesInterfaceMod ! Liquid volume in ground layer real(r8), allocatable :: h2o_liqvol_gl(:) + ! Site level filter for uptake response functions + logical :: filter_btran + end type bc_in_type From 0e6f6ee0a45fdac7226a2f7185572d16d2ac5162 Mon Sep 17 00:00:00 2001 From: Ryan Knox Date: Wed, 20 Jul 2016 02:39:33 -0700 Subject: [PATCH 9/9] modified some print statements --- biogeophys/EDBtranMod.F90 | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/biogeophys/EDBtranMod.F90 b/biogeophys/EDBtranMod.F90 index ebb41293cd..6966257453 100644 --- a/biogeophys/EDBtranMod.F90 +++ b/biogeophys/EDBtranMod.F90 @@ -210,13 +210,11 @@ subroutine btran_ed( sites, nsites, bc_in, bc_out) enddo temprootr = sum(bc_out(s)%rootr_pagl(ifp,:)) - if(abs(1.0_r8-temprootr) > 1.0e-9_r8)then + if(abs(1.0_r8-temprootr) > 1.0e-10_r8 .and. temprootr > 1.0e-10_r8)then write(iulog,*) 'error with rootr in canopy fluxes',temprootr,sum(pftgs),sum(cpatch%rootr_ft(1:2,:),dim=2) - if(temprootr > 0._r8)then - do j = 1,numlevgrnd - bc_out(s)%rootr_pagl(ifp,j) = bc_out(s)%rootr_pagl(ifp,j)/temprootr - enddo - end if + do j = 1,numlevgrnd + bc_out(s)%rootr_pagl(ifp,j) = bc_out(s)%rootr_pagl(ifp,j)/temprootr + enddo end if cpatch => cpatch%younger