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/EDBtranMod.F90 b/biogeophys/EDBtranMod.F90 index 9c5475059c..6966257453 100644 --- a/biogeophys/EDBtranMod.F90 +++ b/biogeophys/EDBtranMod.F90 @@ -1,352 +1,343 @@ module EDBtranMod - - !------------------------------------------------------------------------------ - ! !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 - ! - type(ed_cohort_type), pointer :: currentCohort ! current cohort - type(ed_patch_type) , pointer :: currentPatch ! current patch - !------------------------------------------------------------------------------ + + !------------------------------------------------------------------------------------- + ! Description: + ! + ! ------------------------------------------------------------------------------------ + + use pftconMod , only : pftcon + use clm_varcon , only : tfrz + 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 + public :: get_active_suction_layers contains - - !------------------------------------------------------------------------------ - subroutine btran_ed( bounds, p, sites, nsites, hsites, & - soilstate_inst, waterstate_inst, temperature_inst, energyflux_inst) - ! - ! !DESCRIPTION: - ! - ! !USES: - 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 - ! + + ! ==================================================================================== + + logical function check_layer_water(h2o_liq_vol, tempk) + + implicit none + ! Arguments + 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 + + ! ===================================================================================== + + 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 :: iv !leaf layer - integer :: s !site - integer :: c !column - integer :: j !soil layer - integer :: ft ! plant functional type index - !---------------------------------------------------------------------- - - ! 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 + integer :: s ! site + integer :: j ! soil layer !------------------------------------------------------------------------------ - - 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) - - 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 - - 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) - - 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) - ) - - if (patch%is_veg(p)) then - - c = patch%column(p) - s = hsites(c) - - currentPatch => map_clmpatch_to_edpatch(sites(s), p) - - 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 - - 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 - currentPatch%rootr_ft(FT,j) = currentPatch%rootr_ft(FT,j)/currentPatch%btran_ft(FT) - else - currentPatch%rootr_ft(FT,j) = 0._r8 - end if + + associate( & + numlevgrnd => ctrl_parms%numlevgrnd ) + + do s = 1,nsites + 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 - - 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 - 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) - 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 - + else + bc_out(s)%active_suction_gl(:) = .false. + end if + end do + + end associate + end subroutine get_active_suction_layers + + ! ===================================================================================== + + 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_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 + ! --------------------------------------------------------------------------------- + + ! 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 + type(ed_cohort_type),pointer :: ccohort ! Current cohort 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 + + bc_out(s)%rootr_pagl(:,:) = 0._r8 + + 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 ( 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) + + 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 + 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 + ! 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_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_pagl(ifp,j) = bc_out(s)%rootr_pagl(ifp,j) + & + 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 + 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_pagl(ifp,:)) + 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) + 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 + 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/EDPhotosynthesisMod.F90 b/biogeophys/EDPhotosynthesisMod.F90 index c5691b08f6..7aeef5eb2a 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..7aedb32eb3 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,100 @@ 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 + + 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(inout) :: 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 + + 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 + + ! 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 + + 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 + ! ** 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 +1065,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_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) = & - forc_par_d*cpatch%fabd_sun_z(CL,ft,iv) + & - forc_par_i*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) = & - forc_par_d*cpatch%fabd_sha_z(CL,ft,iv) + & - forc_par_i*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) @@ -1066,9 +1087,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/EDTypesMod.F90 b/main/EDTypesMod.F90 index 5b9df778f0..d69ebf6dda 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, public :: 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 e62a9d75ba..79ae9d6009 100644 --- a/main/FatesInterfaceMod.F90 +++ b/main/FatesInterfaceMod.F90 @@ -15,34 +15,118 @@ module FatesInterfaceMod ! PUBLIC API!!!! ! ------------------------------------------------------------------------------------ - use EDtypesMod , only : ed_site_type + use EDtypesMod , only : ed_site_type, & + numPatchesPerCol, & + ctrl_parms + + use shr_kind_mod , only : r8 => shr_kind_r8 ! INTERF-TODO: REMOVE THIS + use clm_varpar , only : nlevgrnd + + ! ------------------------------------------------------------------------------------ + ! 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. + ! Naming conventions: _gl means ground layer dimensions + ! _pa means patch dimensions + ! _rb means radiation band + ! ------------------------------------------------------------------------------------ + type, public :: bc_in_type + + ! The actual number of FATES' ED patches + integer :: npatches + + ! Downwelling direct beam radiation (patch,radiation-band) [W/m2] + real(r8), allocatable :: solad_parb(:,:) + + ! 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_gl(:) + + ! Effective porosity = porosity - vol_ic, of layers in each site [-] + real(r8), allocatable :: eff_porosity_gl(:) + + ! volumetric soil water at saturation (porosity) + 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(:) + + ! Site level filter for uptake response functions + logical :: filter_btran + + + 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(:) + + ! 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(:,:) + + ! 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 + + 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 :: zero_bcs end type fates_interface_type + + public :: set_fates_ctrlparms + + contains - ! ------------------------------------------------------------------------------------ + ! ==================================================================================== ! INTERF-TODO: THIS IS A PLACE-HOLDER ROUTINE, NOT CALLED YET... subroutine fates_clean(this) @@ -62,4 +146,172 @@ subroutine fates_clean(this) return end subroutine fates_clean + + ! ==================================================================================== + + + subroutine allocate_bcin(bc_in) + + ! --------------------------------------------------------------------------------- + ! Allocate and Initialze the FATES boundary condition vectors + ! --------------------------------------------------------------------------------- + + implicit none + type(bc_in_type), intent(inout) :: bc_in + + ! Allocate input boundaries + + ! Radiation + allocate(bc_in%solad_parb(numPatchesPerCol,ctrl_parms%numSWBands)) + allocate(bc_in%solai_parb(numPatchesPerCol,ctrl_parms%numSWBands)) + + ! Hydrology + 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 + + 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(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)) + + + return + end subroutine allocate_bcout + + ! ==================================================================================== + + 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_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)%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 + + 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