subroutine crop_growth(avt, bSnowModule, flu, frar, ida, j, je, n, nn, nveg, pit, ra, tmit, tmn, tx, ylc, yls)
!**** PURPOSE: TO SIMULATE PLANT GROWTH
!**** CALLED IN: CRPMD
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
! PARAMETERS & VARIABLES
!
! >>>>> COMMON PARAMETERS & VARIABLES
! actual = actual evapotranspiration, mm
! alai(j, je) = leaf area index
! ald1(icr) = shape parameter for the LAI developement equation
! ald2(icr) = shape parameter for the LAI developement equation
! be(icr) = biomass-energy ratio for crop, kg m2 MJ-1 ha-1 d-1
! blai(icr) = max LAI for crop
! dlai(icr) = fraction of season, when LAI declines
! dm(j, je) = total biomass, kg/ha
! flu(j) = fraction of subbasin area in the basin
! frar(j, je) = fractional areas of hydrotope in subbasin
! g(j, je) = fraction of heat units to maturity accumulated, 1/1
! hi(icr) = harvest index for crop (from database)
! hia(j, je) = harvest index
! hiad(j, je) = harvest index, adjusted
! huharv(j, je) = harvest index heat unit
! hun(icr) = potential heat units required for maturity of crop
! ialpha = switch parameter for CO2 effect on net photosynth.
! ida = current day
! idayx = ida, to calc ndgro - number of growth days
! idlef = code for the day length effect in crop 0/1
! igro(j, je) = vegetation index, =1 if yes
! ilcc(icr) = land cover category --> readcrp
! ndgro = number of growth days
! nucr(j, je) = crop number (database)
! olai(j, je) = alai(j, je) - leaf area index
!block pit = 58.13
! potentl = potential evapotranspiration, mm
! ra(j) = solar radiation in subbasin j, J/cm^2
! rd(j, je) = root depth, mm
! rdmx(icr) = maximum root depth, mm
! rwt(j, je) = fraction of root weight
! sdn = sum N stress days
! sdp = sum P stress days
! sdt = sum temp stress days
! sdw = sum water stress days
! strsn = N stress factor
! strsp = P stress factor
! swh(j, je) = actual transp. by plants, mm
! swp(j, je) = potent. transp. by plants, mm
! tb(icr) = base temp. for plant growth, degree C
! ts = temp. stress
! tsav(j, je) = temp. stress, accumulated
! tx(j) = average daily temp., degree C
! ws(j, je) = water stress
! wsav(j, je) = water stress, accumulated
! ylc(j) = cos(lat()/clt), lat() - lat, clt=57.296
! yls(j) = sin(lat()/clt), lat() - lat, clt=57.296
! >>>>>
! >>>>> STATIC PARAMETERS
! beadj = adjusted on CO2 biomass-energy ratio
! ch = local par
! dayl = day length
! ddm = delta dm
! delg = delta g
! deltalai = delta LAI
! dlef = local par
! duma = local par
! f = fraction of plant's maximum leaf area index corresponding to
! to a given fraction of potential heat units for plant
! ff = delta f for the day
! h = local par
! heat = local par
! par = local par
! reg = local par
! sd = local par
! tgx = local par
! x1 = local par
! x2 = local par
! x3 = local par
! x4 = local par
! xi = local par
! xinc = local par
! xy = local par
! >>>>>
!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
!**** Include common parameters
use vegetation, only : &
vegetation_adjust_energy_ratio, &
vegetation_s_curve, &
vegetation_temperature_stress, &
vegetation_water_stress, &
ald1, &
wsav, &
sdt, &
ialpha, &
blai, &
olai, &
tsav, &
ald2, &
dm, &
tb, &
g, &
alai, &
ilcc, &
ws, &
to, &
dlai, &
potentl, &
sdw, &
rd, &
ts
use nutrient, only : &
nutrient_phosphorus_uptake, &
nutrient_nitrate_uptake, &
strsp, &
strsn
real(dp), dimension(:), intent(in) :: avt
real(dp), dimension(:), intent(in) :: flu
real(dp), dimension(:, :), intent(in) :: frar
integer, intent(in) :: ida
integer, intent(in) :: nn
integer, dimension(:, :), intent(in) :: nveg
real(dp), intent(in) :: pit
real(dp), dimension(:), intent(in) :: ra
real(dp), dimension(:), intent(in) :: tmn
real(dp), dimension(:), intent(in) :: tx
real(dp), dimension(:), intent(in) :: ylc
real(dp), dimension(:), intent(in) :: yls
logical, intent(in) :: bSnowModule
real(dp), intent(in) :: tmit
integer j, je, n
real(dp) beadj, ch, dayl, ddm, delg, deltalai, dlef, duma, f, ff
real(dp) h, heat, par, reg, sd, tgx, x1, x2, x3, x4, xi, xinc, xy
x1 = 0.05
x2 = 0.9
x3 = 10.05
x4 = 50.90
!ToDo: dlef = 0. ! ???
! dlef is not initialized if g(j, je) >=0.5 .AND. idlef.NOT.0
dlef = 0.
!*********************************************************** START IF (IGRO>=1)
!**** CALC DAILY INCREASE IN HEAT UNITS delg & ACCUMULATED HEAT UNITS g()
if (igro(j, je) .ge. 1) then
if (idayx .ne. ida .and. g(j, je) .gt. 0.) then
ndgro = ndgro + 1
idayx = ida
end if
!###########################
!#### SNOW MODULE ####
!###########################
if (bSnowModule) then
delg = (tmit - tb(nucr(j, je))) / hun(nucr(j, je))
else
delg = (tx(j) - tb(nucr(j, je))) / hun(nucr(j, je))
end if
!###########################
!**** CALC day length (from clgen)
xi = ida
sd = .4102 * sin((xi - 80.25) / pit)
ch = - yls(j) * tan(sd) / ylc(j)
if (ch .lt. 1.) then
if (ch .le. - 1.) then
h = 3.1416
else
h = acos(ch)
end if
else
h = 0.
end if
dayl = 7.72 * h
!**** Correction of delg on day length (Jan Gr�fe)
if (g(j, je) .lt. 0.5) then
dlef = 0.25 + 0.75 * (dayl - 8.) / 8.
if (dlef .gt. 1.) dlef = 1.
if (dlef .lt. 0.25) dlef = 0.25
endif
if (idlef .eq. 0) then
delg = delg * 1.
else
delg = delg * dlef
endif
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 - CALC TEMPERATURE STRESS
!###########################
!#### SNOW MODULE ####
!###########################
if (bSnowModule) then
tgx = tmit - tb(nucr(j, je))
else
tgx = tx(j) - tb(nucr(j, je))
end if
!###########################
if (tgx .le. 0.) then
ts = 0.
else
call vegetation_temperature_stress(tgx, j, je, n, avt, nucr, nveg, tmn, to, tx)
end if
!**** COMPUTE ROOT DEPTH
if (ilcc(nucr(j, je)) .le. 2 .or. ilcc(nucr(j, je)) .eq. 4 &
.or. ilcc(nucr(j, je)) .eq. 5 &
.OR.nucr(j, je) == icc ) then ! cover crop included by S. Liersch
rd(j, je) = 2.5 * g(j, je) * rdmx(nucr(j, je))
if (rd(j, je) .gt. rdmx(nucr(j, je))) rd(j, je) = rdmx(nucr(j, je))
else
rd(j, je) = rdmx(nucr(j, je))
endif
!**** COMPUTE potential increase in biomass - ddm
! STANDARD VERSION: WITHOUT CO2 FERTILIZATION
par = .005 * ra(j) * (1. - exp(- .65 * (alai(j, je) + .05)))
ddm = be(nucr(j, je)) * par
if (ddm .lt. 0.) ddm = 0.
!**** COMPUTE potential increase in biomass - ddm
! VERSION for CLIMATE IMPACT ASSESSMENT: WITH CO2 FERTILIZATION
!#### CALL ADJUSTBE to adjust be if ialpha = 1
if (ialpha .eq. 1) then
call vegetation_adjust_energy_ratio(j, je, beadj, be, nucr, tx)
ddm = beadj * par
if (ddm .lt. 0.) ddm = 0.
endif
!#### CALL NUPTAKE & PUPTAKE - CALC N AND P UPTAKE
call nutrient_nitrate_uptake(j, je, nucr(j, je), dm, flu, frar, g, ida, nn, snup)
call nutrient_phosphorus_uptake(j, je, nucr(j, je), dm, flu, frar, g, ida, nn, spup)
!**** CALC crop growth regulating factor: reg
strsn = 1.
strsp = 1.
reg = amin1(real(ws(j, je), 4), real(ts, 4), real(strsn, 4), real(strsp, 4))
if (reg .lt. 0.) reg = 0.
if (reg .gt. 1.) reg = 1.
tsav(j, je) = tsav(j, je) + ts
wsav(j, je) = wsav(j, je) + ws(j, je)
!**** CALC biomass dm() & root weight fraction rwt()
dm(j, je) = dm(j, je) + ddm * reg
rwt(j, je) = (.4 - .2 * g(j, je))
!**** CALC harvest index under favourable conditions - hia()
! & real(dp) harvest index - hiad()
f = g(j, je) / (g(j, je) + &
exp(ald1(nucr(j, je)) - ald2(nucr(j, je)) * g(j, je)))
ff = f - huharv(j, je)
huharv(j, je) = f
if (igro(j, je) .gt. 0) then
if (g(j, je) .gt. 0.5) then
swh(j, je) = swh(j, je) + actual
swp(j, je) = swp(j, je) + potentl
endif
heat = 100. * (1. - (dlai(nucr(j, je)) - g(j, je)))
!#### CALL SCURVE
call vegetation_s_curve(x1, x2, x3, x4)
hia(j, je) = hi(nucr(j, je)) * 100 * g(j, je) / &
(100 * g(j, je) + exp(11.1 - 10. * g(j, je)))
xinc = 100. * swh(j, je) / (swp(j, je) + 1.e-6)
if (xinc .lt. 0.) xinc = 0.
duma = xinc / (xinc + exp(x1 - x2 * xinc))
hiad(j, je) = hia(j, je) * duma
else
endif
!**** CALC SUM STRESS DAYS
sdw = sdw + (1. - ws(j, je)) * flu(j) * frar(j, je)
sdt = sdt + (1. - ts) * flu(j) * frar(j, je)
sdn = sdn + (1. - strsn) * flu(j) * frar(j, je)
sdp = sdp + (1. - strsp) * flu(j) * frar(j, je)
!**** COMPUTE LEAF AREA INDEX - alai()
if (g(j, je) .le. dlai(nucr(j, je))) then
if (alai(j, je) .gt. blai(nucr(j, je))) &
alai(j, je) = blai(nucr(j, je))
xy = sqrt(reg)
deltalai = ff * blai(nucr(j, je)) * &
(1. - exp(5. * (alai(j, je) - blai(nucr(j, je))))) * xy
alai(j, je) = alai(j, je) + deltalai
if (alai(j, je) .gt. blai(nucr(j, je))) &
alai(j, je) = blai(nucr(j, je))
olai(j, je) = alai(j, je)
else
alai(j, je) = 16. * olai(j, je) * (1. - g(j, je)) ** 2
if (alai(j, je) .gt. olai(j, je)) alai(j, je) = olai(j, je)
end if
end if
!*********************************************************** END IF (G<=1)
end if
!*********************************************************** END IF (IGRO>=1)
return
end subroutine crop_growth