subroutine vegetation_process(j, je, k, n, wet, 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(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, 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