! FILE vegfun.f95 ! ! SUBROUTINES IN THIS FILE CALLED FROM ! subroutine wstress(j, je, k, n) crpmd, vegmd ! subroutine tstress(tgx, j, je, n) crpmd, vegmd ! subroutine scurve(x1, x2, x3, x4) crpmd, adjustbe ! subroutine ascrv(x1, x2, x3, x4, x5, x6) readcrp ! subroutine adjustbe(j, je, xtem, beadj) crpmd module vegetation use utilities, only : dp implicit none ! initial Leaf Area Index for land use type real(dp), save, dimension(:), allocatable :: alai0 integer, parameter :: nih = 7 ! water use rate - depth parameter, used in wstress real(dp) :: ub = 3.065 ! alpha factor for vegetation real(dp), save :: alfa ! plant transpiration, mm real(dp), save :: ep ! temperature stress factor real(dp), save :: ts ! plant evaporation, potential, mm real(dp), save :: potentl real(dp), save :: sdt = 0. ! sum temp stress days real(dp), save :: sdw = 0. ! sum water stress days ! |hour |time threshold used to define dormant , dimension(subtot) real(dp), save, dimension(:), allocatable :: dormhr ! dormancy threshold for day length logical, save :: bDormancy = .false. real(dp), save :: epco = 1. ! limit plant water demand compensation from lower layers. ! pap calc in readbas from ub, ub=3.065 blockdata real(dp), save :: uob ! switch parameter: to calc CO2 effect on net photosynthesis? integer, save :: ialpha ! switch parameter: to calc CO2 effect on transpiration? integer, save :: ibeta ! switch parameter: 3 / 4 - C3 or C4 crop? real(dp), save :: ic3c4 ! atm CO2 in the reference period, ppm real(dp), save :: co2ref ! atm CO2 in the scenario period, ppm real(dp), save :: co2sce real(dp), save, dimension(:), allocatable :: gwqLastday ! maximum additional additional uptake rate in wetlands, mm real(dp), save :: rzmaxup = 0. ! min day length, h, calc. in readwet real(dp), save, dimension(:), allocatable :: daylmn ! day length in subbasin, h, calc in readcli real(dp), save, dimension(:), allocatable :: daylen ! index for dormant period integer, save, dimension(:, :), allocatable :: idorm ! crop residue in two upper soil layers, kg / ha real(dp), save, dimension(:, :, :), allocatable :: rsd ! total biomass, kg / ha real(dp), save, dimension(:, :), allocatable :: dm ! fraction heat units accumulated, 1 / 1 real(dp), save, dimension(:, :), allocatable :: g ! alai(j, je) - leaf area index real(dp), save, dimension(:, :), allocatable :: olai ! leaf area index real(dp), save, dimension(:, :), allocatable :: alai ! root depth, mm real(dp), save, dimension(:, :), allocatable :: rd ! water stress factor real(dp), save, dimension(:, :), allocatable :: ws ! water stress, sum for the growth period real(dp), save, dimension(:, :), allocatable :: wsav ! temp stress, sum for the growth period real(dp), save, dimension(:, :), allocatable :: tsav ! base temp. for plant growth, degree C real(dp), save, dimension(:), allocatable :: tb ! max LAI for crop real(dp), save, dimension(:), allocatable :: blai ! fraction of growing season, when LAI declines real(dp), save, dimension(:), allocatable :: dlai ! minimum Leaf Area Index (for forest and natural vegetation) real(dp), save, dimension(:), allocatable :: almn ! specific leaf area, m2 / kg, LAI / SLA in kg / m2 real(dp), save, dimension(:), allocatable :: sla ! land cover category for vegetation iv integer, save, dimension(:), allocatable :: ilcc ! shape parameter for the LAI developement equation real(dp), save, dimension(:), allocatable :: ald1 ! shape parameter for the LAI developement equation real(dp), save, dimension(:), allocatable :: ald2 ! optimal temperature for plant growth, degrees C real(dp), save, dimension(:), allocatable :: to ! output ids integer :: & transpiration_output_id = 0, & vegetation_water_stress_oid = 0, & vegetation_temp_stress_oid = 0, & heat_unit_fract_output_id = 0, & biomass_total_output_id = 0, & leaf_area_index_output_id = 0, & root_depth_output_id = 0 namelist / VEGETATION_PARAMETERS / & co2ref, & co2sce, & epco, & ialpha, & ibeta, & ic3c4, & rzmaxup, & ub, & bDormancy contains !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& subroutine vegetation_initialise(mb, nlut, mcrdb, meap, lat) use input, only : get_config_fid use output, only : output_register_hydrotope_var integer, intent(in) :: mb, nlut, mcrdb, meap real(dp), intent(in) :: lat(:) read(get_config_fid(), VEGETATION_PARAMETERS) call vegetation_allocate(mb, nlut, mcrdb, meap) call vegetation_prepare_parameters(lat) ! output variables transpiration_output_id = output_register_hydrotope_var("transpiration", .false.) vegetation_water_stress_oid = & output_register_hydrotope_var("vegetation_water_stress") vegetation_temp_stress_oid = & output_register_hydrotope_var("vegetation_temperature_stress") heat_unit_fract_output_id = output_register_hydrotope_var("heat_unit_fraction") biomass_total_output_id = output_register_hydrotope_var("biomass_total") leaf_area_index_output_id = output_register_hydrotope_var("leaf_area_index") root_depth_output_id = output_register_hydrotope_var("root_depth") end subroutine vegetation_initialise subroutine vegetation_prepare_parameters(lat) real(dp), intent(in) :: lat(:) real(dp) :: ch(size(lat)), hrs(size(lat)) ch = .439 * abs(tan(lat / 57.296)) hrs = merge(0., acos(ch), ch >= 1.) daylmn = 7.72 * hrs if (bDormancy) then dormhr = merge(1., (abs(lat) - 20.) / 20., abs(lat) > 40.) dormhr = merge(-0.1, dormhr, abs(lat) < 20.) else dormhr = - 1. end if end subroutine vegetation_prepare_parameters subroutine vegetation_allocate(mb, nlut, mcrdb, meap) integer, intent(in) :: mb, nlut, mcrdb, meap allocate(alai(mb, meap)) allocate(alai0(nlut)) allocate(ald1(mcrdb)) allocate(ald2(mcrdb)) allocate(almn(mcrdb)) allocate(blai(mcrdb)) allocate(daylen(mb)) allocate(daylmn(mb)) allocate(dlai(mcrdb)) allocate(dm(mb, meap)) allocate(dormhr(mb)) allocate(g(mb, meap)) allocate(gwqLastday(mb)) allocate(idorm(mb, meap)) allocate(ilcc(mcrdb)) allocate(olai(mb, meap)) allocate(rd(mb, meap)) allocate(rsd(mb, meap, 2)) allocate(sla(mcrdb)) allocate(tb(mcrdb)) allocate(tsav(mb, meap)) allocate(ws(mb, meap)) allocate(wsav(mb, meap)) allocate(to(mcrdb)) alai(:, :) = 0. alai0 = 0. ald1 = 0. ald2 = 0. almn = 0. blai = 0. daylen = 0. daylmn = 0. dlai = 0. dm = 0. dormhr = 0. g = 0. gwqLastday = 0. idorm = 0 ilcc = 0 olai(: , :) = 0. rd = 0. rsd = 0. sla = 0. tb = 0. tsav = 0. uob = 1. - exp(- ub) ws = 0. wsav = 0. to = 0. end subroutine vegetation_allocate subroutine dealloc_vegetation deallocate(alai) deallocate(alai0) deallocate(ald1) deallocate(ald2) deallocate(almn) deallocate(blai) deallocate(daylen) deallocate(daylmn) deallocate(dlai) deallocate(dm) deallocate(g) deallocate(gwqLastday) deallocate(idorm) deallocate(ilcc) deallocate(olai) deallocate(rd) deallocate(rsd) deallocate(sla) deallocate(tb) deallocate(tsav) deallocate(ws) deallocate(wsav) deallocate(to) end subroutine dealloc_vegetation subroutine vegetation_water_stress(j, je, k, n, wet, disLastDay, additionalGwUptake, bWAM_Module, dart, daycounter, fc, frar, humi, icc, ida, iy, iyr, mstruc, nbyr, nn, nucr, nveg, rdmx, sbar, sep, ste, z) !**** PURPOSE: THIS SUBROUTINE ESTIMATES WATER STRESS FACTOR. ! IT DISTRIBUTES POTENTIAL PLANT EVAPORATION THROUGH ! THE ROOT ZONE AND CALCULATES ACTUAL PLANT WATER USE ! BASED ON SOIL WATER AVAILABILITY !**** CALLED IN: CRPMD, VEGMD !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! PARAMETERS & VARIABLES ! ! >>>>> COMMON PARAMETERS & VARIABLES ! alfa = calculated in adjustbe ! co2ref = CO2 in reference period, ppm ! co2sce = CO2 in simulated period, ppm ! disLastDay = discharge of last day, m3/s ! additionalGwUptake = total additional groundwater uptake, m3/s ! ep = plant evaporation, mm ! g(j, je) = fraction heat units accumulated, 1/1 ! humi(j) = air humidity, % ! ialpha = switch parameter: to calc CO2 effect on net photosynthesis? ! ibeta = switch parameter: to calc CO2 effect on transpiration? ! ida = current day ! ilcc(iv) = land cover category for vegetation iv ! iwshd = number of hydrotope to print from wstress(), if iwstr = 1 ! iwssb = number of subbasin to print from wstress(), if iwstr = 1 ! iwstr = switch code to print from wstress() ! nn = number of soil layers ! nucr(j, je) = crop number (database) ! nveg(j, je) = vegetation number (database) ! potentl = plant evaporation, potential, mm ! precip = precipitation, mm ! rd(j, je) = root depth, mm ! rdmx(iv) = max root depth, mm ! sep = water percolation ! ste(j, je, l) = water storage, recalc here, mm !block ub = 3.065 ! ul(l, k) = upper limit water content in layer, mm ! uob = pap calc in readbas from ub, ub=3.065 blockdata ! ws(j, je) = water stress factor ! z(l, k) = depth of layer, mm ! >>>>> ! >>>>> STATIC PARAMETERS ! beta = beta factor for plant transpiration under higher CO2 ! ccfi = local par ! epin = local par ! gx = local par ! ir = local par ! l = local par ! num = local par ! sum = local par ! u = local par |mm H2O |water uptake by plants in each soil layer ! ul4 = local par ! xx = local par |mm H2O |water uptake by plants from all layers ! xx1 = local par ! xx2 = local par ! xx3 = local par ! >>>>> !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !**** Include common parameters use management, only : & TSubbasin, TWaterUser, management_is_transfer_subbasin, management_subbasin_pointer, & management_user_pointer, management_is_active_period, wam_correct_irrigationdemand, & wam_plantdemand_summary, wam_plantdemand_summary, wam_irrigationdemand_summary use landuse, only : landuse_is_cropland real(dp), dimension(:), intent(inout) :: additionalGwUptake real(dp), dimension(:), intent(inout) :: disLastDay logical, intent(in) :: bWAM_Module real(dp), dimension(:), intent(in) :: dart integer, intent(in) :: daycounter real(dp), dimension(:, :), intent(in) :: fc real(dp), dimension(:, :), intent(in) :: frar real(dp), dimension(:), intent(in) :: humi integer, intent(in) :: icc integer, intent(in) :: ida integer, intent(in) :: iy integer, intent(in) :: iyr integer, dimension(:, :, :), intent(in) :: mstruc integer, intent(in) :: nbyr integer, intent(in) :: nn integer, dimension(:, :), intent(in) :: nucr integer, dimension(:, :), intent(in) :: nveg real(dp), dimension(:), intent(in) :: rdmx real(dp), dimension(:), intent(in) :: sbar real(dp), intent(inout) :: sep real(dp), dimension(:, :, :), intent(inout) :: ste real(dp), dimension(:, :), intent(in) :: z integer j, je, n, k, ir, l, num, wet real(dp) beta, ccfi, epin, gx, sum, ul4, xx, xx1, xx2, xx3 real(dp) areahyd, additionalUptake, co2sce_interpolated real(dp), dimension(10) :: u ! epco as a global calibration parameter in * .bsn real(dp) zu, sump real(dp) rdmx_temp !#### WATER MANAGEMENT MODULE #### real(dp) :: area ! hydrotop area, m2 !#### WATER MANAGEMENT MODULE #### TYPE (TSubbasin), POINTER :: pS ! pointer on subbasins !#### WATER MANAGEMENT MODULE #### TYPE (TWaterUser), POINTER :: pWU ! pointer on actual TWU u = 0. potentl = ep ! The factor "epco" is introduced to limit plant water demand ! compensation from lower layers. ! Without epco (or epco=1) water demand which could not be met ! by an individual layer was moved to the next lower layer. ! epco=0 means no compensation. epco=0.5 means 50% compensation ! from next lower layer. ! It is basically a calibration parameter and reflects ! the plants flexibility to respond to water limitations ! with increased root growth. -- Pia Gottschalk !**** SL !epco=0.25 ! epco is now a calibration paramter in *.bsn !**** SL !**** Including factor beta for plant transpiration under higher CO2 ! as a constant beta coupled to alfa if (ialpha .eq. 1 .and. ibeta .eq. 1) then ! Calculation of transient CO2 increase = simple interpolation ! from start year CO2 concentration (co2ref) ! to end year CO2 concentration (co2sce): co2sce_interpolated = ((co2sce - co2ref) / (nbyr - 1)) * (iy - 1) + co2ref ! beta = alfa * co2sce / co2ref ! Correction of calculation of beta, Pia Gottschalk, 23/08/2011 beta = alfa * co2ref / co2sce_interpolated ep = ep * beta end if epin = ep xx = 0. if (landuse_is_cropland(n) ) then num = nucr(j, je) else num = nveg(j, je) end if !**** CALC plant transpiration ep and water stress ws() ! Correction from F. Badeck: ul4 = ul()*0.6, 60% instead of 25% before ! Correction from P. Gottschalk, 19/07/11: ! ul4 = ul()*0.5, 50% instead of 25% (SWIM original) before ! Correction from P. Gottschalk, 19/07/11, is based on ! "Leitfaden zur Beregnung landwirtschaftlicher Kulturen", ! LVLF, September 2005, Tabelle 3 if (ep .le. 0.) then ws(j, je) = 1. else l = 1 ir = 0 sum = 0. sump = 0. ! select case (ilcc(num)) ! case (1, 2, 4, 5) ! rd(j, je) = 2.5 * g(j, je) * rdmx(num) ! if ( rd(j, je) > rdmx(num) ) rd(j, je)= rdmx(num) ! !if (sol_rd < 10.) sol_rd = 10. ! case default ! rd(j, je) = rdmx(num) ! end select ! Pia: rdmx (max. rooting depth) is adjusted to the depth ! of the arable layer or max rooting depth (from crop.dat) ! if the depth of the arable layer is greater than the ! crop specific max root depth. ! To adjust hydrotop specific water uptake the variable ! "rdmx_temp" is introduced here locally. ! This prevents that the loop over the arable layer to calculate ! water uptake generates a water stress artefact. ! Before, water demand was calculated along the whole root profile ! but the loop only ran over the arable layer. ! Note, this effectively reduces plant water stress ! because max. root depth for winter wheat is set to 2 m ! whereas the arable layer is usually set to a shallower layer. rdmx_temp = min(rdmx(num), z(nn, k)) if (ilcc(num) .le. 2 .or. ilcc(num) .eq. 4 .or. ilcc(num) .eq. 5 & .OR. nucr(j, je) == icc ) then ! why not include cover crop? rd(j, je) = 2.5 * g(j, je) * rdmx_temp if (rd(j, je) .gt. rdmx_temp) then rd(j, je) = rdmx_temp end if else rd(j, je) = rdmx_temp end if xx = 0. zu = 0. do l = 1, nn if (ir .gt. 0) exit if (rd(j, je) .le. z(l, k)) then gx = rd(j, je) ir = l else gx = z(l, k) end if if (rd(j, je) .le. 0.) then sum = ep / uob else if (num .eq. 69) then sum = ep * (1. - exp(- ub)) / uob else sum = ep * (1. - exp(- ub * gx / rd(j, je))) / uob end if end if sump = sum ! u(l) = sum - xx u(l) = sum - sump + (sump - xx) * epco ! Pia: Commented out as the threshold from which plant water ! uptake is reduced due to increasing matrix potential (here 0.5) ! refers to the plant available water (=fc) and not to saturation (=ul) ! Note, this reduces plant water stress. ! ul4 = ul(l, k) * 0.5 ul4 = fc(l, k) * 0.75 ! Pia: Commented out: linear decline of plant water uptake below fc*0.5 ! u(l) = u(l) * ste(j, je, l) / ul4 ! Pia: Now using SWAT approach which features a exponential decline ! of plant water uptake below fc*0.5 if (ste(j, je, l) .lt. ul4 ) u(l) = u(l) * exp(5 * (ste(j, je, l) / ul4 - 1)) if (ste(j, je, l) .lt. u(l) ) u(l) = ste(j, je, l) ste(j, je, l) = ste(j, je, l) - u(l) xx = xx + u(l) zu = z(l, k) end do ! l = 1, nn ws(j, je) = xx / ep !################################# !#### WATER MANAGEMENT MODULE #### !################################# ! calculate plant water demand for irrigated crops if (bWAM_Module) then if (management_is_transfer_subbasin(j) ) then if (landuse_is_cropland(n) .AND. mstruc(j, je, 7) > 0 & ! if HRU is cropland .AND. if value in column IRR in * .str file is > 0 .AND. nucr(j, je) /= icc ) then ! if hydrotope is irrigated cropland, no cover crop (icc) pS => management_subbasin_pointer(j) ! pointer on current subbasin (could be source or destination) pWU => management_user_pointer(pS % pos_irr) ! pointer on current water user area = frar(j, je) * sbar(j) ! HRU area in m^2 !----------------------------------------------------------------- ! Plant and total irrigation requirements are calculated here. ! If irr_opt == 1, these requirements are written to the output files, but ! the amount of water used for irrigation will be overwritten by input time series in: wam_withdraw_Transfer_Out() !!! !----------------------------------------------------------------- if (pWU % wu_opt == 4 .AND. management_is_active_period(iyr, ida, pWU) ) then ! if source is not external but another subbasin ! calculate plant demand but only for output purposes because under this option, ! water is provided (input time series) independently from requirements. ! xx = water uptake by plants from all layers [mm] ! ep = plant transpiration [mm] pWU % plantDemand(daycounter) = real(pWU % plantDemand(daycounter) + (xx - ep) / 1000. * area / 86400. * (- 1.)) ! mm - > m3 / s ! calculate actual irrigation demand including losses depending on irrigation practices pWU % irrigationDemand(daycounter) = wam_correct_irrigationDemand(pWU, pWU % plantDemand(daycounter)) ! summarise total plant demand and irrigation demand for summary output file wam_plantDemand_summary(daycounter) = & wam_plantDemand_summary(daycounter) + pWU % plantDemand(daycounter) wam_irrigationDemand_summary(daycounter) = & wam_irrigationDemand_summary(daycounter) + pWU % irrigationDemand(daycounter) end if ! ( pWU % wu_opt == 4 ) end if ! ( landuse_is_cropland(n) .AND. mstrucx(j, je, 1) >= 1 ) end if ! ( wam_is_transfer_subbasin(j) ) end if ! ( wam_bTransfer ) !################################# ep = xx end if ! if (ep .le. 0.) !**** Correction on air humidity if (ida .gt. 100 .and. humi(j) .gt. 75.) then xx2 = potentl - ep ! ccfi = (humi(j) - 75.)/10. ! ccfi = 2.*(humi(j) - 75.)/10. ccfi = 1. if (ccfi .gt. 1.) ccfi = 1. if (xx2 .gt. 0.) then ! xx3=min(sep, ccfi*xx2) xx3 = ccfi * xx2 if (xx3 .gt. sep) xx3 = sep ep = ep + xx3 sep = sep - xx3 if (sep .lt. 0.) sep = 0. ws(j, je) = ep / potentl end if end if !!!! new Riperian zone Fred, uptake limited by river inflow into subbasin ! (disLastDay) if (wet.ge.1.and.disLastDay(j)>0.) then !if(j.eq.1) print*, "+++++++", disLastDay(j) ! calculate transpiration deficite xx1 xx1 = potentl - ep if (xx1.gt.0.) then ! uptake is limitited to rzmaxup if (xx1>rzmaxup) xx1=rzmaxup ! calc hydrotop area areahyd=frar(j,je)*dart(j) ! calculate additional uptake in m3 additionalUptake = xx1*0.5*areahyd*1000./86400. ! and avoid negative flow if ((disLastDay(j).gt.additionalUptake)) then ep = ep + 0.5 * xx1 sep = sep - 0.5 * xx1 ws(j,je) = ep / potentl additionalGwUptake(j) = additionalGwUptake(j) + additionalUptake disLastDay(j) = disLastDay(j) - additionalUptake else additionalUptake = 0. end if end if !if(j.eq.1) print*, "-------", disLastDay(j) if(disLastDay(j).lt.0.) print*, "!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!", j, disLastDay(j) end if return end subroutine vegetation_water_stress !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& subroutine vegetation_temperature_stress(tgx, j, je, n, avt, nucr, nveg, tmn, to, tx) !**** PURPOSE: COMPUTES TEMPERATURE STRESS FOR CROP/VEGETATION GROWTH !**** CALLED IN: GROWTH, VEGMD !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! PARAMETERS & VARIABLES ! ! >>>>> COMMON PARAMETERS & VARIABLES ! avt(j) = average annual daily temerature., degree C ! nucr(j, je) = crop number ! nveg(j, je) = vegetation number ! tb(iv) = base temp for plant growth, degree C ! tmn(j) = daily min temp, from readcli, degree C ! to(iv) = opt temp for growth, degree C ! ts = temp stress factor ! tx(j) = daily aver temp, degree C ! >>>>> ! >>>>> STATIC PARAMETERS ! num = local par ! rto = local par ! rto1 = local par ! rto2 = local par ! >>>>> !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !**** Include common parameters use landuse, only : landuse_is_cropland use snow, only : bSnowModule, tmin, tx_tmp, tmit, tmin_tmp real(dp), dimension(:), intent(in) :: avt integer, dimension(:, :), intent(in) :: nucr integer, dimension(:, :), intent(in) :: nveg real(dp), dimension(:), intent(in) :: tmn real(dp), dimension(:), intent(in) :: to real(dp), dimension(:), intent(in) :: tx integer j, je, n, num real(dp) tgx, rto, rto1, rto2 !########################### !#### SNOW MODULE #### !########################### if (bSnowModule) then tx_tmp = tmit tmin_tmp = tmin else tx_tmp = tx(j) tmin_tmp = tmn(j) end if !########################### if (landuse_is_cropland(n) ) then num = nucr(j, je) else num = nveg(j, je) endif rto = 0. rto1 = 0. rto2 = 0. !**** CALC temp stress factor ts if (tx_tmp .le. to(num)) then tgx = (to(num) + tb(num)) / (to(num) - tb(num)) rto1 = tgx * (to(num) - tx_tmp) ** 2 rto2 = rto1 / (tx_tmp + 1.e-6) if (rto2 .le. 200.) then ts = exp(- 0.1054 * rto2) else ts = 0. end if else tgx = 2. * to(num) - tb(num) - tx_tmp rto = ((to(num)-tx_tmp) / (tgx + 1.e-6)) ** 2 if (rto .le. 200.) then ts = exp(- 0.1054 * rto) else ts = 0. end if endif if (tmin_tmp .le. (avt(j) - 15.)) ts = 0. return end subroutine vegetation_temperature_stress !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& subroutine vegetation_s_curve(x1, x2, x3, x4) !**** PURPOSE: THIS SUBROUTINE COMPUTES S CURVE PARAMETERS x1 & x2 ! GIVEN 2 (X, Y) POINTS !**** CALLED IN: CRPMD !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ real(dp) x1, x2, x3, x4, xx xx =log(x3 / x1 - x3) x2 = (xx -log(x4 / x2 - x4)) / (x4 - x3) x1 = xx + x3 * x2 return end subroutine vegetation_s_curve !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& subroutine vegetation_s_curve_parameters(x1, x2, x3, x4, x5, x6) !**** PURPOSE: THIS SUBROUTINE COMPUTES S CURVE PARAMETERS x5 & x6 ! GIVEN 2 (X, Y) POINTS !**** CALLED IN: READCRP !~~~~~~~~~~~ ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ real(dp) xx, x1, x2, x3, x4, x5, x6 xx = 0.0 x5 = 0.0 x6 = 0.0 xx = Log(x3 / x1 - x3) x6 = (xx - Log(x4 / x2 - x4)) / (x4 - x3) x5 = xx + (x3 * x6) return end subroutine vegetation_s_curve_parameters !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& subroutine vegetation_adjust_energy_ratio(j, je, beadj, be, nucr, tx) !**** PURPOSE: this subroutine adjusts biomass-energy ratio for crop ! for higher atmosphere CO2 concentrations ! F. Wechsung method (temperature-dependent) !**** CALLED IN: CRPMD !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! PARAMETERS & VARIABLES ! ! >>>>> COMMON PARAMETERS & VARIABLES ! alfa = alpha factor for vegetation ! be(icr) = biomass-energy ratio for crop ! co2ref = atm CO2 in the reference period, ppm ! co2sce = atm CO2 in the scenario period, ppm ! ic3c4 = switch parameter: 3/4 - C3 or C4 crop? ! nucr(j, je) = crop number (database) ! tx(j) = average daily temp., degree C ! >>>>> ! >>>>> STATIC PARAMETERS ! aa = coef ! bb = coef ! c1 = co2ref ! c2 = co2sce ! cc = coef ! ci1 = CO2 internal corr co2ref ! ci2 = CO2 internal corr co2sce ! deg = local par ! del = local par ! delq = local par ! xlnalfa = local par ! >>>>> !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !**** Include common parameters use snow, only : bSnowModule, tx_tmp, tmit real(dp), dimension(:), intent(in) :: be integer, dimension(:, :), intent(in) :: nucr real(dp), dimension(:), intent(in) :: tx integer j, je real(dp) beadj, aa, bb, c1, c2, cc, ci1, ci2, deg, del, delq, xlnalfa !########################### !#### SNOW MODULE #### !########################### if (bSnowModule) then tx_tmp = tmit else tx_tmp = tx(j) end if !########################### !**** CALC xlnalfa = f(tx) aa = 0.3898E-2 bb = 0.3769E-5 cc = 0.3697E-4 c1 = co2ref c2 = co2sce ci1 = c1 * 0.7 ci2 = c2 * 0.7 del = ci2 - ci1 delq = ci2 * ci2 - ci1 * ci1 xlnalfa = aa * del - bb * delq + cc * del * tx_tmp !**** CALC degree - different for C3 and C4 crops if (ic3c4 .eq. 4) then deg = 0.36 else deg = 0.6 endif !**** CALC alfa and beadj alfa = exp(xlnalfa) ** deg if (alfa .lt. 1) alfa = 1. beadj = be(nucr(j, je)) * alfa return end subroutine vegetation_adjust_energy_ratio !&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&&& subroutine vegetation_process(j, je, k, n, wet, disLastDay, additionalGwUptake, avt, bWAM_Module, be, cva, dart, daycounter, fc, flu, frar, huharv, humi, hun, icc, ida, iy, iyr, mstruc, nbyr, nn, nucr, nveg, ra, rdmx, rwt, sbar, sep, snup, spup, ste, sumfc, swe, tmn, to, tx, uap, z, alai, olai) !**** PURPOSE: CALC daily growth of plant biomass for natural vegetation !**** CALLED IN: HYDROTOP !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! PARAMETERS & VARIABLES ! ! >>>>> COMMON PARAMETERS & VARIABLES ! alai(j, je) = leaf area index ! ald1(iv) = shape parameter for the LAI developement equation ! for veget iv ! ald2(iv) = shape parameter for the LAI developement equation ! for veget iv ! almn(iv) = LAI minimum (for forest and natural vegetation) ! be(iv) = biomass-energy ratio for crop, kg m2 MJ-1 ha-1 d-1 ! blai(iv) = max LAI for crop ! cva(j, je) = vegetation cover, kg/ha ! daylen(j) = day length in subbasin, h, calc in readcli ! daylmn(j) = min day length, h, calc. in readwet ! dlai(iv) = fraction of season, when LAI declines ! dm(j, je) = total biomass, kg/ha ! flu(j) = fraction of subbasin area in the basin ! g(j, je) = fraction of heat units to maturity accumulated ! huharv(j, je) = harvest index heat unit ! hun(iv) = potential heat units required for maturity of crop ! ida = current day ! idorm(j, je) = index for dormant period ! ih1, ih2, ih3 = hydrotopes for HYDROTOPE PRINTOUT ! isb1, isb2, isb3 = subbasins for HYDROTOPE PRINTOUT ! iy = current year as counter (1, ..., nbyr) ! iyr = current year ! nveg(j, je) = number of vegetation from crop database ! olai(j, je) = alai(j, je) - leaf area index ! ra(j) = solar radiation in subbasin j, J/cm^2 ! rd(j, je) = root depth, mm ! rdmx(iv) = maximum root depth, mm ! rsd(j, je, 2) = residue, kg/ha ! rwt(j, je) = fraction of root weight ! sdt = sum temp stress days ! sdw = sum water stress days ! sla(iv) = specific leaf area, m2/kg, LAI/SLA in kg/m2 ! tb(iv) = base temperature for plant growth, degrees C ! ts = temp. stress ! tsav(j, je) = temp. stress, accumulated ! tx(j) = daily aver temp in the subbasin, degrees C ! uap = P uptake, kg/ha ! ws(j, je) = water stress ! wsav(j, je) = water stress, accumulated ! >>>>> ! >>>>> STATIC PARAMETERS ! ddm = delta dm ! delg = delta g ! deltalai = delta LAI ! f = fraction of plant's maximum leaf area index corresponding ! to a given fraction of potential heat units for plant ! ff = delta f for the day ! par = local par ! reg = local par ! resnew = new residue ! tgx = local par ! xx = local par ! yy = local par ! zz = local par ! >>>>> !~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ !**** Include common parameters use snow, only: bSnowModule, tmit, tx_tmp use nutrient, only : nutrient_phosphorus_uptake, nutrient_nitrate_uptake real(dp), dimension(:, :), intent(inout) :: alai real(dp), dimension(:, :), intent(inout) :: olai real(dp), dimension(:), intent(inout) :: additionalGwUptake real(dp), dimension(:), intent(inout) :: disLastDay real(dp), dimension(:), intent(in) :: avt logical, intent(in) :: bWAM_Module real(dp), dimension(:), intent(in) :: be real(dp), dimension(:, :), intent(inout) :: cva real(dp), dimension(:), intent(in) :: dart integer, intent(in) :: daycounter real(dp), dimension(:, :), intent(in) :: fc real(dp), dimension(:), intent(in) :: flu real(dp), dimension(:, :), intent(in) :: frar real(dp), dimension(:, :), intent(inout) :: huharv real(dp), dimension(:), intent(in) :: humi integer, dimension(:), intent(in) :: hun integer, intent(in) :: icc integer, intent(in) :: ida integer, intent(in) :: iy integer, intent(in) :: iyr integer, dimension(:, :, :), intent(in) :: mstruc integer, intent(in) :: nbyr integer, intent(in) :: nn integer, dimension(:, :), intent(in) :: nucr integer, dimension(:, :), intent(in) :: nveg real(dp), dimension(:), intent(in) :: ra real(dp), dimension(:), intent(in) :: rdmx real(dp), dimension(:, :), intent(inout) :: rwt real(dp), dimension(:), intent(in) :: sbar real(dp), intent(inout) :: sep real(dp), dimension(:, :), intent(inout) :: snup real(dp), dimension(:, :), intent(inout) :: spup real(dp), dimension(:, :, :), intent(inout) :: ste real(dp), dimension(:), intent(in) :: sumfc real(dp), dimension(:, :), intent(in) :: swe real(dp), dimension(:), intent(in) :: tmn real(dp), dimension(:), intent(in) :: to real(dp), dimension(:), intent(in) :: tx real(dp), intent(out) :: uap real(dp), dimension(:, :), intent(in) :: z integer j, je, k, n, wet real(dp) ddm, delg, deltalai, f, ff, par, reg, resnew, tgx, xx, yy, zz integer :: veg_end real(dp) :: veg_wstress, veg_dieoff, dieoff ! ToDo: veg_end (end of vegetation period to restart growth) ! this parameter is to be defined somewhere as input parameter, not fixed like here! veg_end = 365 veg_wstress = .1 veg_dieoff = .99 dieoff = 0. !########################### !#### SNOW MODULE #### !########################### if (bSnowModule) then tx_tmp = tmit else tx_tmp = tx(j) end if !########################### uap = 0. ts = 0. !**** CALC land cover for natural vegetation: cva(j, je) = 1000. * alai(j, je) + rsd(j, je, 1) !**** CHECK if start of dormant period ! CALC rsd(), dm(), alai(), g() ! Residue allocation formula is changed: to check! ! if (idorm(j, je).eq.0) then ! if (daylen(j)-1..lt.daylmn(j)) then ! if (iy.ne.1.or.ida.ge.180) then if (bDormancy) then if (idorm(j, je) == 0 ) then if (daylen(j) + dormhr(j) < daylmn(j) ) then if (iy > 1.OR.ida >= veg_end) then idorm(j, je) = 1 resnew = (olai(j, je) - alai(j, je)) * 10000. / sla(nveg(j, je)) ! sla = specific leaf area, m2 / kg, LAI / SLA in kg / m2 rsd(j, je, 1) = rsd(j, je, 1) + resnew * 0.5 rsd(j, je, 2) = rsd(j, je, 2) + resnew * 0.5 dm(j, je) = dm(j, je) - resnew alai(j, je) = almn(nveg(j, je)) g(j, je) = 0. else idorm(j, je) = 1 alai(j, je) = almn(nveg(j, je)) g(j, je) = 0. end if end if end if else ! ( bDormancy ) if (idorm(j, je) .eq. 0) then if (daylen(j) - 1. .lt. daylmn(j)) then if (iy .ne. 1 .or. ida .ge. 180) then idorm(j, je) = 1 resnew = (olai(j, je) - alai(j, je)) * 10000. / sla(nveg(j, je)) ! sla = specific leaf area, m2 / kg, LAI / SLA in kg / m2 rsd(j, je, 1) = rsd(j, je, 1) + resnew * 0.5 rsd(j, je, 2) = rsd(j, je, 2) + resnew * 0.5 dm(j, je) = dm(j, je) - resnew alai(j, je) = almn(nveg(j, je)) g(j, je) = 0. else idorm(j, je) = 1 alai(j, je) = almn(nveg(j, je)) g(j, je) = 0. end if end if end if end if ! ( bDormancy ) !**** check if end of dormant period if (bDormancy) then if ((idorm(j, je) >= 1).AND.(daylen(j) + dormhr(j) >= daylmn(j)).AND.ida < veg_end ) then idorm(j, je) = 0 end if else if (idorm(j, je) .ge. 1) then if (daylen(j) - 1. .ge. daylmn(j)) then idorm(j, je) = 0 end if end if end if ! ( bDormancy ) !**** Assuming that the root depth of natural vegetation is not allowed to decrease rd(j, je) = rdmx(nveg(j, je)) !#### CALL WSTRESS call vegetation_water_stress(j, je, k, n, wet, disLastDay, additionalGwUptake, bWAM_Module, dart, daycounter, fc, frar, humi, icc, ida, iy, iyr, mstruc, nbyr, nn, nucr, nveg, rdmx, sbar, sep, ste, z) !sl begin !########################################################### !**** This code snippet has been included to better account for vegetation growth ! controlled by water stress rather than temperature or daylength. ! Vegetation growth is limited to water availability using the soil water index. ! This part is only active if parameter nat_veg is 1 (*.bsn) if (bDormancy) then ! swe(j, jea)/sumfc(k) = soil water index (0-1) ! swe = soil water content [mm] ! sumfc = sum of field capacity in soil [mm] if (swe(j, je) / sumfc(k) <= veg_wstress ) then idorm(j, je) = 1 ! vegetation stops growing (is dormant due to water stress) olai(j, je) = alai(j, je) ! olai = alai of previous day alai(j, je) = alai(j, je) * veg_dieoff ! assuming dying vegetation up to a certain amount per time step alai(j, je) = max(alai(j, je), almn(nveg(j, je))) ! calcuate residue from difference between acutal lai and lai of previous day resnew = (olai(j, je) - alai(j, je)) * 10000. / sla(nveg(j, je)) ! sla = specific leaf area, m2 / kg, LAI / SLA in kg / m2 resnew = max(0.0_dp, resnew) rsd(j, je, 1) = rsd(j, je, 1) + resnew * .5 ! half of this residue is contributing to residue of upper storage rsd(j, je, 2) = rsd(j, je, 2) + resnew * .5 ! the other half contributes to second storage dm(j, je) = dm(j, je) - resnew ! reduce biomass by value of residue dm = max(0.0_dp, dm) else idorm(j, je) = 0 end if ! fraction of heat units and alai need to be initialized after growing season ! ToDo: the definition of the growing season needs re-thinking! if (idorm(j, je) == 0 ) then if (g(j, je) >= 1. ) then g(j, je) = 0. !alai(j, je) = almn(nveg(j, je)) end if end if if (ida == veg_end) then g(j, je) = 0. end if end if ! if (bDormancy) !########################################################### !sl end ! if (idorm(j, je).ge.1) goto 10 if (idorm(j, je) == 0) then !**** COMPUTE DAILY INCREASE IN HEAT UNITS delg delg = (tx_tmp - tb(nveg(j, je))) / hun(nveg(j, je)) if (delg .lt. 0.) delg = 0. g(j, je) = g(j, je) + delg if (g(j, je) .gt. 1.) g(j, je) = 1. !*********************************************************** START IF (G<=1) !**** GROWTH SEASON if (g(j, je) .le. 1) then !#### CALL TSTRESS: COMPUTE TEMPERATURE STRESS tgx = tx_tmp - tb(nveg(j, je)) if (tgx .le. 0.) then ts = 0. else call vegetation_temperature_stress(tgx, j, je, n, avt, nucr, nveg, tmn, to, tx) end if !**** CALC daily biomass increase: ddm par = .005 * ra(j) * (1. - exp(- .65 * (alai(j, je) + .05))) ddm = be(nveg(j, je)) * par if (ddm .lt. 0.) ddm = 0. !#### CALL NUPTAKE, PUPTAKE: CALCULATE N AND P UPTAKE call nutrient_nitrate_uptake(j, je, nveg(j, je), dm, flu, frar, g, ida, nn, snup) call nutrient_phosphorus_uptake(j, je, nveg(j, je), dm, flu, frar, g, ida, nn, spup) !**** CALC BIOMASS dm(), ROOT WEIGHT rwt(9 xx = dm(j, je) + ddm reg = amin1(real(ws(j, je), 4), real(ts, 4)) if (reg .lt. 0.) reg = 0. if (reg .gt. 1.) reg = 1. dm(j, je) = dm(j, je) + ddm * reg rwt(j, je) = (.4 - .2 * g(j, je)) tsav(j, je) = tsav(j, je) + ts wsav(j, je) = wsav(j, je) + ws(j, je) !**** CALC f, ff, huharv() f = g(j, je) / (g(j, je) + exp(ald1(nveg(j, je)) - ald2(nveg(j, je)) * g(j, je))) ff = f - huharv(j, je) huharv(j, je) = f !**** CALC SUM STRESS DAYS sdw = sdw + (1. - ws(j, je)) * flu(j) sdt = sdt + (1. - ts) * flu(j) !**** CALC LAI and adjust for lower limit of LAI for forest alnm() if (g(j, je) .le. dlai(nveg(j, je)) ) then if (alai(j, je) .gt. blai(nveg(j, je)) ) alai(j, je) = blai(nveg(j, je)) deltalai = ff * blai(nveg(j, je)) * (1. - exp(5. * (alai(j, je) - blai(nveg(j, je))))) * sqrt(reg) alai(j, je) = alai(j, je) + deltalai if (alai(j, je) .gt. blai(nveg(j, je)) ) alai(j, je) = blai(nveg(j, je)) olai(j, je) = alai(j, je) if (alai(j, je) .lt. almn(nveg(j, je)) ) alai(j, je) = almn(nveg(j, je)) else yy = sqrt(1. - g(j, je)) zz = 1. / sqrt(1. - dlai(nveg(j, je))) alai(j, je) = zz * olai(j, je) * yy if (alai(j, je) .lt. almn(nveg(j, je)) ) alai(j, je) = almn(nveg(j, je)) end if ! ( g(j, je) .le. dlai(nveg(j, je)) ) end if ! (g(j, je) .le. 1) !*********************************************************** END IF (G<=1) end if ! ( idorm(j, je) == 0) end subroutine vegetation_process subroutine vegetation_store_output(j, je) use output, only : output_store_hydrotope_value integer, intent(in) :: j, je call output_store_hydrotope_value(vegetation_water_stress_oid, j, je, ws(j, je)) call output_store_hydrotope_value(vegetation_temp_stress_oid, j, je, ts) call output_store_hydrotope_value(heat_unit_fract_output_id, j, je, g(j, je)) call output_store_hydrotope_value(biomass_total_output_id, j, je, dm(j, je)) call output_store_hydrotope_value(leaf_area_index_output_id, j, je, alai(j, je)) call output_store_hydrotope_value(root_depth_output_id, j, je, rd(j, je)) end subroutine vegetation_store_output end module vegetation