| blob | aa02408401d54f8ede4293f1ef4e272df28e90f5 |
1 !***********************************************************************
2 !* GECROS for early and late covering crops *
3 !* Genotype-by-Environment interaction on CROp growth Simulator *
4 !* *
5 !* Author: Xinyou YIN *
6 !* Crop and Weed Ecology Group *
7 !* Wageningen University & Research Centre *
8 !* PO Box 430, 6700 AK Wageningen, Netherlands *
9 !* *
10 !* Modified and extended for winter crops *
11 !* by Joachim INGWERSEN *
12 !* Biogeophysics Group *
13 !* University of Hohenheim *
14 !* *
15 !***********************************************************************
17 MODULE module_sf_gecros
19 implicit none
21 !Runtime constant variables defined within this subroutine
22 !Crop-specific parameters (former *.tbl file)
23 REAL, PARAMETER :: EG = 1. !Efficiency of germination
24 REAL, PARAMETER :: CFV = 0.48 !Carbon fraction in vegetative organs
25 REAL, PARAMETER :: YGV = 0.81 !Grwoth efficiency for vegetative growth
26 REAL, PARAMETER :: FFAT = 0.02 !Fraction of fat in storage organs
27 REAL, PARAMETER :: FLIG = 0.06 !Fraction of lignin in storage organs
28 REAL, PARAMETER :: FOAC = 0.02 !Fraction of organic acids in storage organs
29 REAL, PARAMETER :: FMIN = 0.02 !Fraction of minerals in storage organs
30 REAL, PARAMETER :: LNCI = 0.03586 !Initial value of LNC (g N g-1)
31 REAL, PARAMETER :: TBD = 0.0 !Base temperature of phenology (C)
32 REAL, PARAMETER :: TOD = 22.5 !Optimum temperature of phenology (C)
33 REAL, PARAMETER :: TCD = 37.0 !Ceiling temperature of phenology (C)
34 REAL, PARAMETER :: TSEN = 1.0 !Curvature of temperature response (1)
35 REAL, PARAMETER :: SPSP = 0.2 !DS for start of photoperiod-sensitive phase
36 REAL, PARAMETER :: EPSP = 0.7 !DS for end of photoperiod-sensitive phase
37 REAL, PARAMETER :: CDMHT = 492.6 !Stem dry weight per unit plant height (g m-2 m2)
38 REAL, PARAMETER :: PMEH = 0.6468 !Fraction of sigmoid curve inflexion in plant height growth period (1)
39 REAL, PARAMETER :: ESDI = 1.35 !ESD for indeterminate crops
40 REAL, PARAMETER :: PMES = 0.50 !Fraction of sigmoid curve inflexion in seed growth growth period (1)
41 REAL, PARAMETER :: TBDV = -1.3 !Base temperature of vernalization (C)
42 REAL, PARAMETER :: TODV = 4.9 !Optimum temperature of vernalization (C)
43 REAL, PARAMETER :: TCDV = 15.7 !Ceiling temperature of vernalization (C)
44 REAL, PARAMETER :: NUPTX = 0.5/24./3600. !Maximum crop nitrogen uptake (g N m-2 s-1)
45 REAL, PARAMETER :: SLA0 = 0.0237 !Specific leaf area constant (m2 leaf g-1)
46 REAL, PARAMETER :: SLNMIN = 0.35 !Base SLN for photosynthesis (g N m-2 leaf)
47 REAL, PARAMETER :: RNCMIN = 0.005 !Minimum N concentration in roots (g N g-1)
48 REAL, PARAMETER :: STEMNC = 0.01 !Minimum N concentration in stem (g N g-1)
49 REAL, PARAMETER :: RLVDS = 0.0904 !Rate of C turnover from dead leaves to litter (d-1)
50 REAL, PARAMETER :: DSCRIT = 0.225 !DS at which winter dormancy ends (1)
51 REAL, PARAMETER :: EAJMAX = 48270.!Energy of activation for JMAX (J mol-1)
52 REAL, PARAMETER :: XVN = 24.96 !Slope of linear relationship between JMAX and leaf N (mu-mol CO2 s-1 g-1 N)
53 REAL, PARAMETER :: XJN = 49.92 !Slope of linear relationship between VCMX and leaf N (mu-mol e- s-1 g-1 N)
54 REAL, PARAMETER :: NPL = 390.0 !Number of plants per m2
55 REAL, PARAMETER :: SEEDW = 0.0475 !Seed weight (g seed-1)
56 REAL, PARAMETER :: SEEDNC = 0.020 !Seed N concentration (g N g-1)
57 REAL, PARAMETER :: BLD = 25.58 !Leaf angle (from horizontal) (degree)
58 REAL, PARAMETER :: HTMX = 1.1 !Maximum plant height (m)
59 REAL, PARAMETER :: MTDV = 46.12 !Minimum thermal days for vegetative growth phase (d)
60 REAL, PARAMETER :: MTDR = 40.22 !Minimum thermal days for reproductive phase (d)
61 REAL, PARAMETER :: PSEN = -.104 !Photoperiod sensitivity of phenological development (h-1)
62 REAL, PARAMETER :: C3C4 = 1. !C3C4 = 1. for C3 crops; C3C4 = -1. for C4 crops.
64 !Runtime-constant general parameters
65 REAL, PARAMETER :: LEGUME = -1. !LEGUME = 1. for leguminous crops; = -1. for non-leguminous crops.
66 REAL, PARAMETER :: DETER = 1. !DETER = 1. for determinate crops; = -1. for indeterminate crops.
67 REAL, PARAMETER :: SLP = -1. !SLP = 1. for short-day crops; = -1. for long-day crops.
68 REAL, PARAMETER :: NSUP1 = 0. !NH4-N supply in g/m2/s
69 REAL, PARAMETER :: NSUP2 = 1./172800. !NO3-N supply in g/m2/s
70 REAL, PARAMETER :: CO2A=350. !Atmospherric CO2 concentration (ppm)
71 REAL, PARAMETER :: FCRSH=0.5 !Initial fraction of C in shoot (-)
72 REAL, PARAMETER :: FNRSH=0.63 !Initial fraction of N in shoot (-)
73 REAL, PARAMETER :: PNPRE=0.7 !Proportion of seed N that comes from non-structural stems during seed fill
74 REAL, PARAMETER :: CB=0.75 !Factor for initial N concentration of seed fill
75 REAL, PARAMETER :: CX=1.00 !Factor for final N concentration of seed fill
76 REAL, PARAMETER :: TM=1.5 !DS when transition from CB to CX is fastest
77 REAL, PARAMETER :: RSS=100. !Soil resistance to evaporation (s/m)
78 REAL, PARAMETER :: LS=0. !Lodging sverity (value between zero and unity)
79 REAL, PARAMETER :: WRB=0.25 !Critical root weight density (g/m2/cm depth)
80 REAL, PARAMETER :: THETA=0.7 !Convexity for light response curve of electron transport (J2) in photosynthesis
81 REAL, PARAMETER :: PNLS=1. !Fraction of dead leaf N incorporated into soil litter N
82 REAL, PARAMETER :: INSP=-2. !Inclination of sun angle (degree)
83 REAL, PARAMETER :: CCFIX=6. !Carbon cost of symbiotic N fixation (g C/g N)
84 REAL, PARAMETER :: RDMX=130. !Maximum rooting depth (m)
85 REAL, PARAMETER :: TCP=86400. !Conversion factor from day into sec
86 REAL, PARAMETER :: Z1244=0.272727272727 !12./44.
87 REAL, PARAMETER :: LOG005=-2.995732274 !LOG(0.05)
89 !Runtime constant variables computed in driver.f
90 REAL :: YGO !Growth efficiency of storage organs (g C/g C)
91 REAL :: CFO !Carbon fraction in storage organs (g C/g)
92 REAL :: LNCMIN !Minimum N concentration in leaves (g N/g)
93 REAL :: CLVI !Initial value of CLV
94 REAL :: CRTI !Initial value of CRT
95 REAL :: NLVI !Initial value of NLV
96 REAL :: NRTI !Initial value of NRT
97 REAL :: HTI !Initial value of HT
98 REAL :: RDI !Initial value of RD
100 !*** Debugging on/off
101 LOGICAL :: debugging=.false.
103 CONTAINS
105 SUBROUTINE gecros (DOY, DT, CROP, RB, RT, RTS, FB, SNOWH, & !I
106 WN, SFCTMP, EAIR, RSD, RLD, PRCP, WUL, WLL, WCMIN, LWIDTH, & !I & !I
107 STATE_GECROS, & !H
108 ATRJC, ATRJS, FSR, FRSU, ARSWSU, ARSWSH) !O
110 IMPLICIT NONE
112 ! character(len=19), INTENT(IN) :: nowdate ! string of current date e.g. 2012-30-10_12:00:00
113 INTEGER, INTENT(IN) :: CROP !CROP=1 -> early-covering crop, CROP=2 -> later covering crop
114 REAL, INTENT(IN) :: DOY !Julian day of the current time step
115 REAL, INTENT(IN) :: DT !Integration time step (s)
116 REAL, INTENT(IN) :: RB !Leaf boundary layer resistance (s/m)
117 REAL, INTENT(IN) :: RT !Aerodynamic canopy resistance (s/m)
118 REAL, INTENT(IN) :: RTS !Aerodynamic ground resistance (s/m)
119 REAL, INTENT(IN) :: FB !Fraction of vegetation cover by snow
120 REAL, INTENT(IN) :: SNOWH !Snow height (cm)
121 REAL, INTENT(IN) :: WN !Wind speed (m/s)
122 REAL, INTENT(IN) :: SFCTMP !Air temperature (K)
123 REAL, INTENT(IN) :: EAIR !Vapour pressure (Pa)
124 REAL, INTENT(IN) :: RSD !Downwelling shortwave radiation (W/m2)
125 REAL, INTENT(IN) :: RLD !Downwelling longwave raidation (W/m2)
126 REAL, INTENT(IN) :: PRCP !Precipiration (kg/m2/s)
127 REAL, INTENT(IN) :: WUL !Soil water in the upper layer (L/m2 or mm)
128 REAL, INTENT(IN) :: WLL !Soil water in the lower layer (L/m2 or mm)
129 REAL, INTENT(IN) :: WCMIN !Minimum soil watter content (m3/m3)
130 REAL, INTENT(IN) :: LWIDTH !Leaf width (m) read in from MPTABLE.TBL
131 REAL, INTENT(OUT) :: ATRJC !Absorbed global radiation by canopy (W/m2)
132 REAL, INTENT(OUT) :: ATRJS !Absorbed global radiation by soil (W/m2)
133 REAL, INTENT(OUT) :: FSR !Reflected shortwave radiation (W/m2)
134 REAL, INTENT(OUT) :: FRSU !Fraction of sunlit leaves in the canopy
135 REAL, INTENT(OUT) :: ARSWSU !Actual stomatal resistance of sunlit leaves (s/m)
136 REAL, INTENT(OUT) :: ARSWSH !Actual stomatal resistance of shaded leaves (s/m)
137 REAL, DIMENSION(1:60), INTENT(INOUT) :: STATE_GECROS !Array with Gecros state variables
139 !*** Gecros state variables
140 REAL :: DS !Development stage
141 REAL :: LAI !Green leaf area index
142 REAL :: TLAI !Total leaf area index
143 REAL :: CTDU !Cumulative thermal day unit (d)
144 REAL :: CVDU !Cumulative thermal day unit of vernalization (d)
145 REAL :: CLV !Carbon in living leaves (g C/m2)
146 REAL :: CLVD !Carbon in dead leaves (g C/m2)
147 REAL :: CSST !Carbon in structural stems (g C/m2)
148 REAL :: CSO !Carbon in storage organs (g C/m2)
149 REAL :: CSRT !Carbon in structural roots (g C/m2)
150 REAL :: CRTD !Carbon in dead roots (g C/m2)
151 REAL :: CLVDS !Carbon in dead leaves that have become litter in soil (g C/m2)
152 REAL :: NRT !Nitrogen in living roots (g N/m2)
153 REAL :: NST !Nitrogen in stems (g N/m2)
154 REAL :: NLV !Nitrogen in living leaves (g N/m2)
155 REAL :: NSO !Nitrogen in storage organs (g N/m2)
156 REAL :: TNLV !Total leaf nitrogen (including N in senesced leaves) (g N/m2)
157 REAL :: NLVD !Nitrogen in dead leaves (g N/m2)
158 REAL :: NRTD !Nitrogen in dead roots (g N/m2)
159 REAL :: CRVS !Carbon in stem reserves (g C/m2)
160 REAL :: CRVR !Carbon in root reserves (g C/m2)
161 REAL :: NREOE !NRES accumulated till the end of seed number determining period (g N/m2)
162 REAL :: NREOF !NRES accumulated till the moment at which seed fill starts (g N/m2)
163 REAL :: DCDSR !Short fall of carbon demand for seed fill in previous time steps (g C/m2)
164 REAL :: DCDTR !Short fall of carbon demand for structural stems in previous time steps (g C/m2)
165 REAL :: SLNB !SLN in bottom leaves of canopy (g N/m2)
166 REAL :: LAIC !Carbon determined LAI (m2/m2)
167 REAL :: RMUL !Total of RMUN+RMUA+RMUS+RMLD
168 REAL :: NDEMP !Crop nitrogen demand of the previous time step (g N/m2/s)
169 REAL :: NSUPP !Nitrogen supply of the previous time step (g N/m2/s)
170 REAL :: NFIXT !Total symbiotically fixed nitrogen during growth (g N/m2)
171 REAL :: NFIXR !Reserve pool of symbiotically fixed nitrogen (g N/m2)
172 REAL :: DCDTP !Carbon demand for structural stem growth of the previous time step (g C/m2/s)
173 REAL :: HT !Canopy height (m)
174 REAL :: TPCAN !Cumulative canopy photosynthesis over growth period (g CO2/m2)
175 REAL :: TRESP !Total crop respiratory cost during growth (g CO2/m2)
176 REAL :: TNUPT !Total crop nitrogen uptake during growth (g N/m2)
177 REAL :: LITNT !Total litter nitrogen entering soil during growth (g N/m2)
178 REAL :: WSO !Yield (g/m2)
179 REAL :: WSTRAW !Straw (g/m2)
180 REAL :: GrainNC !Nitrogen in grain (kg N/ha)
181 REAL :: StrawNC !Nitrogen in straw (kg N/ha)
182 REAL :: APCAN !Actual gross canopy photosynthesis (g CO2/m2/d)
184 character(Len=12) :: inputstring
185 INTEGER :: nowday, minutes
186 REAL :: SC !Solar constant (W/m2)
187 REAL :: SINLD !Seasonal offset of sine of solar height (-)
188 REAL :: COSLD !Amplitude of sine of solar height (-)
189 REAL :: DAYL !Day length (h)
190 REAL :: DDLP !Day length for photoperiodism (h)
191 REAL :: DSINBE !DSINB to correct for lower atmospheric transmission at lower solar elevation (s/d)
192 REAL :: DVP !Vapour pressure (kPa)
193 REAL :: WNM !Wind speed (m/s)
194 REAL :: TAIR !Air temperature (C)
195 REAL :: ROOTD !Rooting depth (m)
196 REAL :: WLV !Dry weight of living leaves (g/m2)
197 REAL :: WST !Dry weight of stems (g/m2)
198 REAL :: WRT !Dry weight of roots (g/m2)
199 REAL :: WSH !Dry weight of shoot (g/m2)
200 REAL :: WLVD !Dry weight of dead leaves (g/m2)
201 REAL :: WRTD !Dry weight of dead roots (g/m2)
202 REAL :: CRT !Carbon in living roots (g C/m2)
203 REAL :: CSH !Carbon in living shoot (g C/m2)
204 REAL :: NSH !Nitrogen in living shoot (g N/m2)
205 REAL :: DLAI !LAI of dead leaves still hanging on stem (m2/m2)
206 REAL :: ESD !DS for end of seed-number determing phase (-)
207 REAL :: KCRN !EXtinction coefficient of root nitrogen (m2/g C)
208 REAL :: NFIXD !Crop demand-determined NFIX (g N/2/s)
209 REAL :: KR !Extinction coefficient of root weight density over soil septh (1/cm)
210 REAL :: HNCCR !Critical shoot N concentration ( g N/g)
211 REAL :: FVPD !Slope of linear effect of VPD of intercelluar to ambient CO2 ratio (1/kPa)
212 REAL :: NRETS !Total crop-residue N returned to soil (g N/m2)
213 REAL :: WCUL !Soil water content in uppler layer (m3/m3)
214 REAL :: DWSUP !Water supply for evapotranspiration (mm/s)
215 REAL :: CSRTN !N determined CSRT (g C/m2)
216 REAL :: NRES !Estimated vegetative-organ N remobilizable for seed fill (g N/m2)
217 REAL :: ONC !N concentration in storage organ (g N/g)
218 REAL :: RNC !N concentration in roots (g N/g)
219 REAL :: LNC !N concentration in living leaves (g N/g)
220 REAL :: KL !Extinction coefficient diffuse component for PAE (m2/m2 leaf)
221 REAL :: CTOT !Total C in living shoots and roots (g C/m2)
222 REAL :: NSHH !N in shoots (excluding dead leaves incorporated into soil litter (g N/m2)
223 REAL :: NTOT !Total N in living shoots and roots (g N/m2)
224 REAL :: WSHH !Dry weight of shoots organs (excluding shadded leaves) (g/m2)
225 REAL :: TSN !Total seed number (seeds/m2)
226 REAL :: HNC !Actural N concentration in living shoot (g N/g)
227 REAL :: PSO !Protein content in storage organs (g protein/m2)
228 REAL :: KLN !Intermediate variable to compute KN (g N/m2 leaf)
229 REAL :: NBK !Intermediate variable to compute KN (g N/m2 leaf)
230 REAL :: KW !Wind speed extinction coefficient in canopy (m2/m2 leaf)
231 REAL :: WTOT !Dry weight of total living plant parts (g/m2)
232 REAL :: CCHKIN !C in crop accumulated since start of simulation
233 REAL :: NCHKIN !N in crop accumulated since start of simulation
234 REAL :: LCRT !Rate of C loss in rooots because of senescence (g C/m2/s)
235 REAL :: TSW !Thousand seed weight (g)
236 REAL :: PNC !N concentration in living plant material (g N/m2)
237 REAL :: NDEMD !Defiency-driven N demand (g N/m2/s)
238 REAL :: FCRVR !Fraction of new root C partitioned to root reserves (g C/g C)
239 REAL :: KN !Leaf N extinction coefficient in the canopy (m2/m2 leaf)
240 REAL :: CCHK !Diff. between C added to the crop since start and net total C fluxes rel. to CCHKIN (%)
241 REAL :: HI !Harvest index (g/g)
242 REAL :: LWRT !Rate of loss of root biomass because of senescence (g/m2/s)
243 REAL :: NCHK !As CCHK but for N
244 REAL :: LAIN !N-determined LAI (m2/m2)
245 REAL :: LNRT !Rate of loss of N because of senescence (g N/m2/s)
246 REAL :: LWLVM !Intermediate variable to compute LWLV (g/m2/s)
247 REAL :: SLNNT !Value of SLNT with small plant-N increment (g N/m2 leaf)
248 REAL :: SLNBC !SLNB calculated from exponential N profile in canopy (g N/m2 leaf)
249 REAL :: SLN !Specific leaf N content (g N/m2 leaf)
250 REAL :: SLA !Specific leaf area (m2 leaf/g)
251 REAL :: LWLV !Rate of loss of leaf weight because of leaf senescence (g/m2/s)
252 REAL :: AESOIL !Actual soil evaporation (mm/s)
253 REAL :: APCANN !APCANS with small plant-N increment (g CO2/m2/s)
254 REAL :: APCANS !Actual standing canopy CO2 assimilation (g CO2/m2/s)
255 REAL :: ATCAN !Actual canopy transpiration (mm/s)
256 REAL :: DAPAR !PAR absorbed by canopy (J/m2/s)
257 REAL :: DIFS !Daytime average soil-air temperature difference (C)
258 REAL :: DIFSH !Daytime average shaded leaf-air temperature difference (C)
259 REAL :: DIFSU !Daytime average sunlit leaf-air temperature difference (C)
260 REAL :: FRSH !Fraction of shaded leaves in the canopy
261 REAL :: HOD !Hour of the day
262 REAL :: PESOIL !Potential soil evaporation (mm/s)
263 REAL :: PPCAN !Potential canopy assimilation (g CO2/m2/s)
264 REAL :: PTCAN !Potential canopy transpiration (mm/s)
265 REAL :: RCAN !Canopy resistance (s/m)
266 REAL :: RSLNB !Rate of change in SLNB (g N/m2 leaf/s)
267 REAL :: LNLV !Rate of loss of leaf N because of senescence (g N/m2/s)
268 REAL :: LCLV !Rate of loss of leaf C because of senescence (g C/m2/s)
269 REAL :: RMRE !Residual maintenance respiration (g CO2/m2/s)
270 REAL :: TAVSS !Soil surface temperature (C)
271 REAL :: RMN !RM calculated with a small increment in plant N (g CO2/m2/s)
272 REAL :: VDU !Rate of Vernalization day unit increase (d/s)
273 REAL :: TDU !Rate of thermal day unit increase (d/s)
274 REAL :: RM !Non-growth components of respiration, excluding the cost of N fixation (g CO2/m2/s)
275 REAL :: LVDS !Rate of transfer of C from dead leaves to litter (g C/m2/s)
276 REAL :: NFIXE !Available energy-determined NFIX (g N/m2/s)
277 REAL :: DVR !Phenological development rate (1/s)
278 REAL :: RX !Respiratory cost of N fixation (g CO2/m2/s)
279 REAL :: RNSUPP !Rate of change in NSUPP (g N/m2/s)
280 REAL :: NSUP !N supply to crop (g N/m2/s)
281 REAL :: NFIX !Symbioticall fixed N (g N/m2/s)
282 REAL :: LITN !Litter N entering soil (g N/m2/s)
283 REAL :: LITC !Litter C entering soil (g C/m2/s)
284 REAL :: FDH !Expected relative growth rate of plant height (1/s)
285 REAL :: FDS !Expected relative growth rate of storage organs (1/s)
286 REAL :: SHSAN !SHSA calculated with a small increment in plant-N (g C/g C/s)
287 REAL :: SHSA !Relative shoot activity (g C/g C/s)
288 REAL :: RMUS !Respiratory cost of mineral uptake (g CO2/m2/s)
289 REAL :: NDEMAD !Intermediate variable related to NDEM (g N/m2/s)
290 REAL :: NDEMA !Activity-driven NDEM (g N/m2/s)
291 REAL :: NCR !Intermediate variable (g N/g C)
292 REAL :: DERI !First order-derivative of SHSA with respect to crop N/C ratio (g C/g N/s)
293 REAL :: ASSA !Assimilates available from current photosynthesis for growth (g CO2/m2/s)
294 REAL :: RNFIXR !Rate of change in NFIXR (g N/m2/s)
295 REAL :: RNDEMP !Rate of change in NDEMP (g N/m2/s)
296 REAL :: RMLD !Respiration due to phloem loading of C assimilates to root (g CO2/m2/s)
297 REAL :: RCSRT !Rate of change in CSRT (g C/m2/s)
298 REAL :: NUPTN !Nitrate-N uptake by the crop (g N/m2/s)
299 REAL :: NUPTA !Ammonium-N uptake by the crop (g N/m2/s)
300 REAL :: NDEM !Crop N demand (g N/m2/s)SUBROUTINE ENERGY
301 REAL :: FNSH !Fraction of newly absorbed N partitioned to shoot (g N/g N)
302 REAL :: FCSH !Fraction of new V partitioned to shoot (g C/g C)
303 REAL :: DCSR !C supply from current photosynthesis for root growth (g C/m2/s)
304 REAL :: SLNT !SLN for top leaves in canopy (g N/m2 leaf)
305 REAL :: RMUN !Respiratory cost of nitrate-N uptake (g CO2/m2/s)
306 REAL :: RMUA !Respiratory cost of ammonium-N uptake (g CO2/m2/s)
307 REAL :: NUPT !Crop N uptake (g N/m2/s)
308 REAL :: DCDSC !C demand for seed filling at current time step (g C/m2/s)
309 REAL :: DCDS !C demand for filling of storage organgs at current time step (g C/m2/s)
310 REAL :: FLWCS !Flow of current assimilated C to storage organs (g C/m2/s)
311 REAL :: DCSS !C supply from current photosynthesis for shoot growth (g C/m2/s)
312 REAL :: DCST !C suppyl from current photosynthesis for structural stem growth (g C/m2/s)
313 REAL :: FCSO !Fraction of new C partitioned to storage organs (g C/m2/s)
314 REAL :: RRMUL !Rate of change in RMUL (g CO2/m2/s)
315 REAL :: RHT !Rate of change in HT (m/s)
316 REAL :: IFSH !Integral factor of stresses on plant height growth (-)
317 REAL :: GAP !Gap betweeb C supply and C demand for seed growth (g C/m2/s)
318 REAL :: CREMSI !Intermediate variable to compute CREMS (g C/m2/s)
319 REAL :: CREMS !C remobilized from stem reserves to storage organs (g C/m2/s)
320 REAL :: CREMRI !Intermediate variable to compute CREMR (g C/m2/s)
321 REAL :: CREMR !C remobilized from root reserves to storage organs (g C/m2/s)
322 REAL :: RWSO !Rate of change in storage organs(g/m2/s)
323 REAL :: RWRT !Rate of change in roots (g/m2/s)
324 REAL :: RRD !Rate of change in RD (cm/s)
325 REAL :: RDCDTP !Rate of change in DCDTP (g C/m2/s)
326 REAL :: RDCDSR !Rate of change in DCDSR (g C/m2/s)
327 REAL :: RCSST !Rate of change in CSST (g C/m2/s)
328 REAL :: RCSO !Rate of change in CSO (g C/m2/s)
329 REAL :: RCRVR !Rate of change in CRVR (g C/m2/s)
330 REAL :: FLWCT !Flow of assimilated C to structural stems (g C/m2/s)
331 REAL :: FCSST !Fraction of new shoot C partitioned to structural stems (g C/m2/s)
332 REAL :: FCLV !Fraction of new shoot C partitioned to leaves (g C/m2/s)
333 REAL :: DCDTC !C demand of structural stem growth at current time step (g C/m2/s)
334 REAL :: DCDT !C demand for the growth of structural stems (g C/m2/s)
335 REAL :: FCRVS !Fraction of new C paritioned to stem reserves (g C/m2/s)
336 REAL :: RCLV !Rate of change in CLV (g C/m2/s)
337 REAL :: RCRVS !Rate of change in CRVS (g C/m2/s)
338 REAL :: RDCDTR !Rate of change in DCDTR (g C/m2/s)
339 REAL :: RESTOT !Total respiratory cost (g CO2/m2/s)
340 REAL :: RG !Growth respiration (g CO2/m2/s)
341 REAL :: RLAI !Rate of change in LAI (m2 leaf/m2/s)
342 REAL :: RNLV !Rate of change in NLV (g N/m2/s)
343 REAL :: RNREOE !Rate of change in NREOE (g N/m2/s)
344 REAL :: RNREOF !Rate of change in NREOF (g N/m2/s)
345 REAL :: RNRES !Rate of change in NRES (g N/m2/s)
346 REAL :: RNRT !Rate of change in NRT (g N/m2/s)
347 REAL :: RNSO !Rate of change in NSO (g N/m2/s)
348 REAL :: RNST !Rate of change in NSZ (g N/m2/s)
349 REAL :: RRP !Respiration/photosynthesis ratio (-)
350 REAL :: RTNLV !Rate of change in TNLV (g N/m2/s)
351 REAL :: RWLV !Rate of change in WLV (g/m2/s)
352 REAL :: RWST !Rate of change in WST (g/m2/s)
353 REAL :: GLAT !Latitude (degree)
354 REAL :: SD1 !Thickness of upper evaporative soil layer (cm) (equals RDI)
356 !*** write STATE_GECROS array into Gecros variables
357 DS = STATE_GECROS(1)
358 CTDU = STATE_GECROS(2)
359 CVDU = STATE_GECROS(3)
360 CLV = STATE_GECROS(4)
361 CLVD = STATE_GECROS(5)
362 CSST = STATE_GECROS(6)
363 CSO = STATE_GECROS(7)
364 CSRT = STATE_GECROS(8)
365 CRTD = STATE_GECROS(9)
366 CLVDS = STATE_GECROS(10)
367 NRT = STATE_GECROS(11)
368 NST = STATE_GECROS(12)
369 NLV = STATE_GECROS(13)
370 NSO = STATE_GECROS(14)
371 TNLV = STATE_GECROS(15)
372 NLVD = STATE_GECROS(16)
373 NRTD = STATE_GECROS(17)
374 CRVS = STATE_GECROS(18)
375 CRVR = STATE_GECROS(19)
376 NREOE = STATE_GECROS(20)
377 NREOF = STATE_GECROS(21)
378 DCDSR = STATE_GECROS(22)
379 DCDTR = STATE_GECROS(23)
380 SLNB = STATE_GECROS(24)
381 LAIC = STATE_GECROS(25)
382 RMUL = STATE_GECROS(26)
383 NDEMP = STATE_GECROS(27)
384 NSUPP = STATE_GECROS(28)
385 NFIXT = STATE_GECROS(29)
386 NFIXR = STATE_GECROS(30)
387 DCDTP = STATE_GECROS(31)
388 HT = STATE_GECROS(32)
389 ROOTD = STATE_GECROS(33)
390 TPCAN = STATE_GECROS(34)
391 TRESP = STATE_GECROS(35)
392 TNUPT = STATE_GECROS(36)
393 LITNT = STATE_GECROS(37)
394 GLAT = STATE_GECROS(44)
395 WSO = STATE_GECROS(45)
396 WSTRAW = STATE_GECROS(46)
397 GrainNC = STATE_GECROS(47)
398 StrawNC = STATE_GECROS(48)
399 LAI = STATE_GECROS(49)
400 TLAI = STATE_GECROS(50)
401 SD1 = STATE_GECROS(52)
403 ! Used for debugging
404 ! if (nowdate(9:12).eq.'2500') then
405 ! write(*,*) nowdate, ' 1 ', DS
406 ! write(*,*) nowdate, ' 2 ', CTDU
407 ! write(*,*) nowdate, ' 3 ', CVDU
408 ! write(*,*) nowdate, ' 4 ', CLV
409 ! write(*,*) nowdate, ' 5 ', CLVD
410 ! write(*,*) nowdate, ' 6 ', CSST
411 ! write(*,*) nowdate, ' 7 ', CSO
412 ! write(*,*) nowdate, ' 8 ', CSRT
413 ! write(*,*) nowdate, ' 9 ', CRTD
414 ! write(*,*) nowdate, ' 10 ', CLVDS
415 ! write(*,*) nowdate, ' 11 ', NRT
416 ! write(*,*) nowdate, ' 12 ', NST
417 ! write(*,*) nowdate, ' 13 ', NLV
418 ! write(*,*) nowdate, ' 14 ', NSO
419 ! write(*,*) nowdate, ' 15 ', TNLV
420 ! write(*,*) nowdate, ' 16 ', NLVD
421 ! write(*,*) nowdate, ' 17 ', NRTD
422 ! write(*,*) nowdate, ' 18 ', CRVS
423 ! write(*,*) nowdate, ' 19 ', CRVR
424 ! write(*,*) nowdate, ' 20 ', NREOE
425 ! write(*,*) nowdate, ' 21 ', NREOF
426 ! write(*,*) nowdate, ' 22 ', DCDSR
427 ! write(*,*) nowdate, ' 23 ', DCDTR
428 ! write(*,*) nowdate, ' 24 ', SLNB
429 ! write(*,*) nowdate, ' 25 ', LAIC
430 ! write(*,*) nowdate, ' 26 ', RMUL
431 ! write(*,*) nowdate, ' 27 ', NDEMP
432 ! write(*,*) nowdate, ' 28 ', NSUPP
433 ! write(*,*) nowdate, ' 29 ', NFIXT
434 ! write(*,*) nowdate, ' 30 ', NFIXR
435 ! write(*,*) nowdate, ' 31 ', DCDTP
436 ! write(*,*) nowdate, ' 32 ', HT
437 ! write(*,*) nowdate, ' 33 ', ROOTD
438 ! write(*,*) nowdate, ' 34 ', TPCAN
439 ! write(*,*) nowdate, ' 35 ', TRESP
440 ! write(*,*) nowdate, ' 36 ', TNUPT
441 ! write(*,*) nowdate, ' 37 ', LITNT
442 ! read(*,*)
443 ! endif
445 PPCAN=0.
446 APCANS=0.
447 APCANN=0.
448 APCAN=0.
449 PTCAN=0.
450 ATCAN=0.
451 PESOIL=0.
452 AESOIL=0.
453 DIFS=0.
454 DIFSU=0.
455 DIFSH=0.
456 DAPAR=0.
457 RCAN=0.
458 DVR=0.
460 nowday = INT(DOY)
461 HOD = float(nint((DOY-int(DOY))*86400.))/3600.
463 ! Conversion from K to C
464 TAIR = SFCTMP - 273.15
466 ! Conversion of rel. humidity into VP (kPa)
467 DVP = EAIR*0.001 !Converts EAIR from Pa to kPa
468 WNM = MAX (0.1, WN)
470 ! Photoperiod, solar constant and daily extraterrestrial radiation
471 CALL ASTRO(aint(DOY),GLAT,INSP,SC,SINLD,COSLD,DAYL,DDLP,DSINBE)
473 ! Plant weights (g/m2)
474 WLV = CLV / CFV
475 WST = CSST / CFV + CRVS/0.444
476 WSO = CSO / CFO
477 WRT = CSRT / CFV + CRVR/0.444
478 WLVD = CLVD / CFV
479 WRTD = CRTD / CFV
481 ! Carbon in shoot and root (g C/m2)
482 CSH = CLV + CSST + CRVS + CSO
483 CRT = CSRT + CRVR
485 ! Nitrogen in shoot (g C/m2)
486 NSH = NST + NLV + NSO
488 ! Extinction coefficient of root nitrogen (m2/g C)
489 KCRN = -LOG005/6.3424/CFV/WRB/RDMX
491 ! DS for end of seed number determining period
492 ESD = INSW(DETER, ESDI, 1.)
494 ! Dead leaves still hanging on the plant (m2/m2
495 DLAI = (CLVD-CLVDS)/CFV*SLA0
497 ! Total leaf area index (dead plus living leaves)
498 TLAI = MAX(0.01,LAI + DLAI)
500 ! Crop demand-determined NFIX
501 NFIXD = MAX(0., NDEMP - NSUPP)
503 ! Critical shoot nitrogen concentration g N g-1
504 HNCCR = LNCI*EXP(-.4*DS)
506 !Extinction coefficient of root weight density over the soil depth cm-1
507 KR = -LOG005/RDMX !cm-1
509 ! Slope of linear effect of VPD on intercelluar to ambient CO2 ratio (1/kPa)
510 FVPD = INSW (C3C4, 0.195127, 0.116214)
512 ! Total crop-residue nitrogen returned to soil (g N/m2)
513 NRETS = LITNT+INSW(DS-2.,0.,NLV+NST+NRT+NFIXR+(CLVD-CLVDS)/ &
514 CFV*LNCMIN*(1.+PNLS)/2.)
516 !if (nowdate(1:12).eq.'201110061930') then
517 !write(*,*) 'Bis hier her bin ich gekommen', ROOTD
518 !endif
520 ! Soil water content of the upper and lower layer (m3/m3)
521 WCUL = (WUL+WCMIN*10.*ROOTD)/10./ROOTD
523 ! Daily water supply for evapotranspiration (mm/s)
524 DWSUP = MAX(.1/TCP,WUL/TCP+.1/TCP)
526 ! Nitrogen-determined CSRT (g C/m2)
527 CSRTN = 1./KCRN*LOG(1.+KCRN*MAX(0.,(NRT*CFV-CRVR*RNCMIN))/RNCMIN)
529 ! Straw weight (g/m2)
530 WSTRAW = WLV + WST + (WLVD - CLVDS/CFV)
532 ! Shoot weight (g/m2)
533 WSH = WLV + WST + WSO
535 ! Estimated vegetative-organ N remobilizable for seed growth (g N/m2)
536 NRES = NREOF + (NREOE-NREOF)*(ESD-1.)/NOTNUL(MIN(DS,ESD)-1.)
538 ! Nitrogen concentration in living leaves LAI(g N/m2)
539 LNC = NLV / NOTNUL(WLV)
541 ! Nitrogen concentration in roots (g N/m2)
542 RNC = NRT / NOTNUL(WRT)
544 ! Nitrogen concentration seeds (g N/m2)
545 ONC = INSW(-WSO, NSO/NOTNUL(WSO), 0.)
547 ! Nitrogen in grain and straw in kg N per ha
548 GrainNC = NSO*10.
549 StrawNC = (NLV+NST)*10.
551 ! Extinction coefficient of nitrogen and wind
552 CALL KDIFF (TLAI,BLD*3.141592654/180.,0.2, KL)
554 CTOT = CSH + CRT
555 NSHH = NSH +(WLVD-CLVDS/CFV)*LNCMIN
556 NTOT = NSH + NRT
557 WSHH = WSH + (WLVD-CLVDS/CFV)
558 TSN = NRES/PNPRE/SEEDNC/SEEDW
559 HNC = NSH / NOTNUL(WSH)
561 ! Amount of seed protein
562 PSO = 6.25*WSO*ONC
564 KLN = KL*(TNLV-SLNMIN*TLAI)
566 NBK = SLNMIN*(1.-EXP(-KL*TLAI))
567 KW = KL
568 WTOT = WSH + WRT
570 ! Crop carbon balance check
571 CCHKIN = CTOT + CLVD + CRTD -CLVI-CRTI
573 ! Crop nitrogen balance check
574 NCHKIN = NTOT + NLVD + NRTD -NLVI-NRTI
576 LCRT = MAX(MIN(CSRT-1.E-4,CSRT-MIN(CSRTN,CSRT)),0.)/TCP !Eq. 43, p. 38, DELT ok, gC m-2 s-1
578 FCRVR = INSW(CSRTN-CSRT, 1., 0.)
579 PNC = NTOT / NOTNUL(WTOT)
581 KN = 1./TLAI*LOG(MAX(1.001,(KLN+NBK)/(KLN*EXP(-KL*TLAI)+NBK)))
583 TSW = WSO/NOTNUL(TSN)*1000.
585 NDEMD = INSW(DS-1., WSH*(HNCCR-HNC)*(1.+NRT/MAX(1E-2,NSH))/TCP, 0.) !Eq. 20, p.25, DELT ok, g N m-2 s-1
587 ! Biomass formation
588 HI = WSO / NOTNUL(WSHH)
589 CCHK = (CCHKIN-(TPCAN-TRESP)*Z1244)/NOTNUL(CCHKIN)*100.
590 NCHK = (NCHKIN-TNUPT)/NOTNUL(TNUPT)*100.
591 LWRT = LCRT/CFV !g m-2 s-1
593 ! Leaf area development
594 LAIN = LOG(1.+ KN*MAX(0.,NLV)/SLNMIN)/KN
595 LNRT = LWRT*RNCMIN !g N m-2 s-1
597 ! Green leaves still hanging on the plant (m2/m2)
598 LAI = MAX(0.01,MIN(LAIN, LAIC))
600 ! Leaf senescence
601 ! The equation differs somewhat from Yin. The right hand side is not divided by deltaT
602 ! This is done by computing LWLV
603 ! LWLVM = (LAIC-MIN(LAIC,LAIN))/SLA0/DELT !Eq. 42, p.36, DELT ok, g m-2 s-1
605 LWLVM = (LAIC-MIN(LAIC,LAIN))/SLA0 !Eq. 42, p.36, DELT ok, g m-2 s-1
607 ! Specific leaf nitrogen and its profile in the canopy
608 SLN = NLV/LAI
609 SLNT = NLV*KN/(1.-EXP(-KN*LAI))
610 SLNBC = NLV*KN*EXP(-KN*LAI)/(1.-EXP(-KN*LAI))
611 SLNNT = (NLV+0.001*NLV)*KN /(1.-EXP(-KN*LAI))
612 SLA = LAI/NOTNUL(WLV)
614 ! ji, 1402, Stay-green maize -> no senescense
615 LWLV = MIN((WLV-1.E-5)/TCP, (LWLVM+REANOR(ESD-DS,LWLVM)*0.03*WLV)/TCP) !g m-2 s-1, gecheckt ok
617 !write(*,*) nowdate, DWSUP,CO2A,LS,EAJMAX,XVN,XJN,THETA,WCUL,FVPD,RB,RT,RTS
618 !read(*,*)
619 ! Call TOTPT: Computes daily total photosynthesis and stomatal resistance
620 CALL TOTPT(HOD,DS,SC,SINLD,COSLD,DAYL,DSINBE,RSD,TAIR,DVP,WNM,C3C4,LAI, &
621 TLAI,HT,LWIDTH,ROOTD,SD1,RSS,BLD,KN,KW,SLN,SLNT,SLNNT,SLNMIN,FB,&
622 DWSUP,CO2A,LS,EAJMAX,XVN,XJN,THETA,WCUL,FVPD,RB,RT,RTS, &
623 PPCAN,APCANS,APCANN,APCAN,PTCAN,ATCAN,PESOIL,AESOIL,DIFS,DIFSU, &
624 DIFSH,DAPAR,RCAN,ATRJC,ATRJS,FSR,FRSU,FRSH,ARSWSU,ARSWSH)
626 RSLNB = (SLNBC-SLNB)/TCP !!Eq. 38, p. 35, DELT ok, g N m-2 leaf s-1
628 LNLV = MIN(LWLV,LWLVM)*LNCMIN + (LWLV-MIN(LWLV,LWLVM))*LNC !g m-2 s-1
629 LCLV = LWLV*CFV !g m-2 s-1
631 ! Residual maintenance respiration g m-2 s-1
632 RMRE = MAX(MIN(44./12.*0.218*(NTOT-WSH*LNCMIN-WRT*RNCMIN)/TCP &
633 ,APCAN-(1.E-5+RMUL)/TCP), 0.) !0.218 has the unit g C g-1 N d-1!
635 TAVSS = TAIR + DIFS
637 ! Call TUNIT: Computes cumulative thermal units (CTDU)
638 CALL TUNIT (1.*CROP, DS,TAIR,MAX(0.,DIFS),DAYL,TBD,TOD,TCD,TBDV,TODV,TCDV,TSEN,TDU,VDU)
640 RMN = MAX(0., MIN(APCAN-1.E-5/TCP,RMUL/TCP) + MAX(MIN(44./12.*0.218* &
641 (1.001*NTOT-WSH*LNCMIN-WRT*RNCMIN)/TCP,APCAN-(1.E-5+RMUL)/TCP), 0.))
643 RM = MAX(0., MIN(APCAN-1.E-5/TCP,RMUL/TCP) + RMRE) !gC m-2 s-1
645 ! Daily and total C and N returns from crop to soil
646 ! ji: RLVDS eingefuehrt. Standardmaessig war der Wert von RLVDS auf 0.1 hard-gecoded
647 ! CLVD: Amount of carbon in dead leaves
648 ! CLVDS: Amount of carbon in dead leaves that have become litter in soil
650 IF(DS.lt.0.25) THEN
651 LVDS = (CLVD-CLVDS)/TCP !gC m-2 s-1
652 ELSE
653 LVDS = RLVDS*(CLVD-CLVDS)*(TAVSS-TBD)/(TOD-TBD)/TCP !gC m-2 s-1
654 ENDIF
656 CALL PHENO (1.*CROP,DS,SLP,DDLP,SPSP,EPSP,PSEN,MTDV,MTDR,TDU,CVDU,DVR)
658 NFIXE = MAX(0., APCAN-(1.E-5+RM)/TCP)/CCFIX*Z1244
660 ! Daily carbon flow for seed filling
661 CALL BETAF(DVR,1.,PMES,LIMIT(1.,2.,DS)-1., FDS)
663 LITC = LCRT + LVDS !gC m-2 s-1
664 LITN = LNRT + LVDS/CFV *LNCMIN*PNLS !gN m-2 s-1
666 CALL BETAF(DVR,(1.+ESD)/2.,PMEH*(1.+ESD)/2.,MIN((1.+ESD)/2.,DS), FDH)
668 NSUP = NSUP1 + NSUP2 !gN m-2 s-1
670 NFIX = INSW (LEGUME, 0., MIN(NFIXE, NFIXD)) !N fixation
672 RNSUPP = (NSUP - NSUPP)/DT !RNSUPP g N m-2 s-2
673 RX = 44./12.*(CCFIX*NFIX)
675 !*** Current photo-assimilates (g CO2 m-2 s-1) for growth, and R/P ratio
676 ASSA = APCAN - RM - RX
678 !*** Crop nitrogen demand and uptake (g N m-2 s-1)
680 SHSA = Z1244 * YGV*MAX(1E-16,APCAN -RM -RX)/ MAX(0.1,CSH)
681 SHSAN = Z1244 * YGV*(APCANN-RMN-RX)/ MAX(0.1,CSH)
682 DERI = MAX(0.,(SHSAN - SHSA)/(0.001*MAX(0.01,NTOT)/MAX(0.1,CTOT)))
684 RMUS = 0.06*0.05/0.454*YGV*ASSA
685 NDEMA = CRT * SHSA**2/NOTNUL(DERI)
686 NCR = INSW(SLNT-SLNMIN,0.,MIN(NUPTX,NDEMA))/(YGV*MAX(1E-16,APCANS-RM-RX)*Z1244)
687 NDEMAD = INSW(LNC-1.5*LNCI, MAX(NDEMA, NDEMD), 0.)
689 !*** Nitrogen partitioning between shoots and roots
690 FNSH = 1./(1.+NCR*DERI/SHSA*CSH/MAX(1E-2,CRT)*NRT/MAX(1E-2,NSH))
692 !*** Carbon partitioning among organs and reserve pools
693 FCSH = 1./(1.+NCR*DERI/SHSA)
695 !*** ji 10.02.14, winter dormancy
696 !*** ji 16.06.15, ecc/lcc switch
697 IF (CROP==1) THEN
698 NDEM = INSW(DS-DSCRIT,INSW(SLNMIN-SLN+1.E-5, MIN(NUPTX,.01*NDEMAD), 0.), &
699 INSW(SLNMIN-SLN+1.E-5, MIN(NUPTX,NDEMAD), 0.))
700 ELSE
701 NDEM = INSW(SLNMIN-SLN+1.E-5, MIN(NUPTX,NDEMAD), 0.)
702 ENDIF
704 DCSR = Z1244*(1.-FCSH)*ASSA
705 NUPTN = MIN(NSUP2, NSUP2/NOTNUL(NSUP)*MAX(0.,NDEM-NFIXR/TCP))
706 RCSRT = Z1244*ASSA*(1.-FCSH)*(1.-FCRVR)*YGV - LCRT
707 RNFIXR = NFIX - MIN(NDEM,NFIXR/TCP)
708 RMLD = 0.06*(1.-FCSH)*ASSA
709 RNDEMP = (NDEM - NDEMP)/DT !as for RNSUPP, DELT ok
710 NUPTA = MIN(NSUP1, NSUP1/NOTNUL(NSUP)*MAX(0.,NDEM-NFIXR/TCP))
712 ! Carbon supply from current photo-assimilates for shoot & root growth
713 DCSS = Z1244* FCSH *ASSA
714 NUPT = MAX(0., NUPTA + NUPTN + MIN(NDEM, NFIXR/TCP))
715 RMUA = 44./12.*0.17*NUPTA
717 CALL SINKG(DS,1.,TSN*SEEDW*CFO,YGO,FDS,DCDSR,DCSS,DT,&
718 DCDSC,DCDS,FLWCS)
720 ! Maintenance and total respiration (g CO2 m-2 d-1)
721 RMUN = 44./12.*2.05*NUPTN
723 ! Daily carbon flow for structural stem growth
724 DCST = DCSS - FLWCS
725 FCSO = FLWCS/DCSS
727 RRMUL = (RMUN+RMUA+RMUS+RMLD-RMUL)/DT
729 GAP = MAX(0., DCDS-DCSS)
731 !*** ji, 16.06.15, ecc/lcc switch
732 IF (CROP==1) THEN
733 IFSH = INSW(DCST-1E-11, 1., LIMIT(0.,1.,DCST/NOTNUL(DCDTP)))
734 ELSE
735 IFSH = LIMIT(0.,1.,DCST/NOTNUL(DCDTP))
736 ENDIF
738 CREMSI = MIN(0.94*CRVS, CRVS/NOTNUL(CRVS+CRVR)*GAP)/0.94
739 CREMRI = MIN(0.94*CRVR, CRVR/NOTNUL(CRVS+CRVR)*GAP)/0.94
740 CREMS = INSW(DCDS-DCSS, 0., CREMSI)
741 CREMR = INSW(DCDS-DCSS, 0., CREMRI)
743 IF (CROP==1) THEN
744 RHT = MIN(HTMX-HT, FDH*HTMX*IFSH)
745 ELSE
746 RHT = MIN(HTMX-HT, FDH*HTMX*INSW(DCST-1E-4, 1., LIMIT(0.,1.,DCST/NOTNUL(DCDTP))))
747 ENDIF
749 IF (CROP==1) THEN
750 CALL SINKG(DS,.1,CDMHT*HTMX*CFV,YGV,FDH*IFSH,DCDTR,DCST,DT, &
751 DCDTC,DCDT,FLWCT)
752 ELSE
753 CALL SINKG(DS,.0,CDMHT*HTMX*CFV,YGV,FDH*IFSH,DCDTR,DCST,DT, &
754 DCDTC,DCDT,FLWCT)
755 ENDIF
757 RCRVR = FCRVR*DCSR - CREMR
759 IF (CROP==1) THEN
760 RDCDTP = (DCDTC-DCDTP)/DT
761 ELSE
762 DCDTP = DCDTC
763 RDCDTP = 0.
764 ENDIF
766 !ji, the MIN function avoids that a situation occurs during which
767 !no assimilates are partioned to the leaves. FCSST max. 85%
768 IF (CROP==1) THEN
769 FCSST = MIN(.85,INSW(DS-(ESD+0.2), FLWCT/MAX(1E-16,DCSS), 0.))
770 !FCSST = INSW(DS-(ESD+0.2), FLWCT/DCSS, 0.) original version
771 ELSE
772 FCSST = MIN(.85,INSW(DS-(ESD+0.2), FLWCT/MAX(1E-16,DCSS), 0.))
773 ENDIF
775 RCSO = Z1244*ASSA*FCSH*FCSO*YGO + 0.94*(CREMS+CREMR)*YGO
776 RWSO = RCSO / CFO
777 RDCDSR = MAX(0., (DCDSC-RCSO/YGO))-(FLWCS-MIN(DCDSC,DCSS))
778 RWRT = RCSRT/CFV + RCRVR/0.444
779 RCSST = Z1244*ASSA* FCSH * FCSST *YGV
781 FCLV = REAAND(LAIN-LAIC,ESD-DS)*(1.-FCSO-FCSST)
783 RRD = INSW(ROOTD-RDMX, MIN((RDMX-ROOTD)/TCP,(RWRT+LWRT)/(WRB+KR* &
784 (WRT+WRTD))), 0.)
786 RCLV = Z1244*ASSA* FCSH * FCLV * YGV - LCLV
788 FCRVS = 1. - FCLV - FCSO - FCSST
790 RDCDTR = MAX(0., (DCDTC-RCSST/YGV))-(FLWCT-MIN(DCDTC,DCST))
791 RCRVS = FCRVS*DCSS - CREMS
792 RWLV = RCLV / CFV
794 RNRES = NUPT-(LNCMIN*(RCLV+LCLV)+RNCMIN*(RCSRT+LCRT)+STEMNC* &
795 RCSST)/CFV
797 RG = 44./12.*((1.-YGV)/YGV*(RCLV+RCSST+RCSRT+LCLV+LCRT)+ &
798 (1.-YGO)/YGO* RCSO)
800 RWST = RCSST/ CFV + RCRVS/0.444
801 RNREOF = INSW (DS-1.0, RNRES, 0.)
802 RNREOE = INSW (DS-ESD, RNRES, 0.)
804 RESTOT = RM+RX+RG + 44./12.*0.06*(CREMS+CREMR)
806 CALL RNACC (FNSH,NUPT,RWST,STEMNC,LNCMIN,RNCMIN,LNC,RNC,NLV,NRT,WLV,WRT, &
807 DT,CB,CX,TM,DS,SEEDNC,RWSO,LNLV,LNRT, RNRT,RNST,RNLV,RTNLV,RNSO)
809 RRP = RESTOT / APCAN
811 CALL RLAIC(1.*CROP,DS,SLA0,RWLV,LAIC,KN,NLV,RNLV,SLNB,RSLNB, RLAI)
813 ! Used for debugging
814 ! if (nowdate(9:12).eq.'2500') then
815 ! write(*,*) nowdate, ' 1 ', DVR
816 ! write(*,*) nowdate, ' 2 ', TDU
817 ! write(*,*) nowdate, ' 3 ', VDU
818 ! write(*,*) nowdate, ' 4 '
819 ! write(*,*) nowdate, ' 5 ', RCLV
820 ! write(*,*) nowdate, ' 6 ', LCLV
821 ! write(*,*) nowdate, ' 7 ', RCSST
822 ! write(*,*) nowdate, ' 8 ', RCSO
823 ! write(*,*) nowdate, ' 9 ', RCSRT
824 ! write(*,*) nowdate, ' 10 ',LCRT
825 ! write(*,*) nowdate, ' 11 ',LVDS
826 ! write(*,*) nowdate, ' 12 ',RNRT
827 ! write(*,*) nowdate, ' 13 ',RNST
828 ! write(*,*) nowdate, ' 14 ',RNLV
829 ! write(*,*) nowdate, ' 15 ',RNSO
830 ! write(*,*) nowdate, ' 16 ',RTNLV
831 ! write(*,*) nowdate, ' 17 ',LNLV
832 ! write(*,*) nowdate, ' 18 ',LNRT
833 ! write(*,*) nowdate, ' 19 ',RCRVS
834 ! write(*,*) nowdate, ' 20 ',RCRVR
835 ! write(*,*) nowdate, ' 21 ',RNREOE
836 ! write(*,*) nowdate, ' 22 ',RNREOF
837 ! write(*,*) nowdate, ' 23 ',RDCDSR
838 ! write(*,*) nowdate, ' 24 ',RDCDTR
839 ! write(*,*) nowdate, ' 25 ',RSLNB
840 ! write(*,*) nowdate, ' 26 ',RLAI
841 ! write(*,*) nowdate, ' 27 ',RRMUL
842 ! write(*,*) nowdate, ' 27 ',RRMUL
843 ! write(*,*) nowdate, ' 28 ',RNDEMP
844 ! write(*,*) nowdate, ' 29 ',RNSUPP
845 ! write(*,*) nowdate, ' 30 ',NFIX
846 ! write(*,*) nowdate, ' 31 ',RNFIXR
847 ! write(*,*) nowdate, ' 32 ',RDCDTP
848 ! write(*,*) nowdate, ' 33 ',RHT
849 ! write(*,*) nowdate, ' 34 ',RRD
850 ! write(*,*) nowdate, ' 35 ',APCAN
851 ! write(*,*) nowdate, ' 36 ',RESTOT
852 ! write(*,*) nowdate, ' 37 ',NUPT
853 ! write(*,*) nowdate, ' 38 ',LITN
854 ! read(*,*)
855 ! endif
857 ! Integration of ODEs
858 DS = MAX(0., INTGRL(DS, DVR, DT))
859 CTDU = MAX(0., INTGRL(CTDU, TDU, DT))
860 CVDU = MAX(0., INTGRL(CVDU, VDU, DT))
861 CLV = MAX(0., INTGRL (CLV, RCLV, DT))
862 CLVD = MAX(0., INTGRL (CLVD, LCLV, DT))
863 CSST = MAX(0., INTGRL (CSST, RCSST, DT))
864 CSO = MAX(0., INTGRL (CSO, RCSO, DT))
865 CSRT = MAX(0., INTGRL (CSRT, RCSRT, DT))
866 CRTD = MAX(0., INTGRL (CRTD, LCRT, DT))
867 CLVDS = MAX(0.,INTGRL (CLVDS, LVDS, DT))
868 NRT = MAX(0.,INTGRL (NRT, RNRT, DT))
869 NST = MAX(0.,INTGRL (NST, RNST, DT))
870 NLV = MAX(0.,INTGRL (NLV, RNLV, DT))
871 NSO = MAX(0.,INTGRL (NSO, RNSO, DT))
872 TNLV = MAX(0.,INTGRL (TNLV, RTNLV, DT))
873 NLVD = MAX(0.,INTGRL (NLVD, LNLV, DT))
874 NRTD = MAX(0.,INTGRL (NRTD, LNRT, DT))
875 CRVS = MAX(0.,INTGRL (CRVS, RCRVS, DT))
876 CRVR = MAX(0.,INTGRL (CRVR, RCRVR, DT))
877 NREOE = MAX(0.,INTGRL(NREOE, RNREOE, DT))
878 NREOF = MAX(0.,INTGRL(NREOF, RNREOF, DT))
879 DCDSR = MAX(0.,INTGRL(DCDSR, RDCDSR, DT))
880 DCDTR = MAX(0.,INTGRL(DCDTR, RDCDTR, DT))
881 SLNB = MAX(0.,INTGRL(SLNB, RSLNB, DT))
882 LAIC = MAX(0.,INTGRL(LAIC, RLAI, DT))
883 RMUL = MAX(0.,INTGRL(RMUL, RRMUL, DT))
884 NDEMP = MAX(0.,INTGRL(NDEMP, RNDEMP, DT))
885 NSUPP = MAX(0.,INTGRL(NSUPP, RNSUPP, DT))
886 NFIXT = MAX(0.,INTGRL(NFIXT, NFIX, DT))
887 NFIXR = MAX(0.,INTGRL(NFIXR, RNFIXR, DT))
888 DCDTP = MAX(0.,INTGRL(DCDTP, RDCDTP, DT))
889 HT = MAX(0.,INTGRL(HT, RHT, DT))
890 ROOTD = MAX(1.,INTGRL(ROOTD, RRD, DT))
891 TPCAN = MAX(0.,INTGRL(TPCAN, APCAN, DT))
892 TRESP = MAX(0.,INTGRL(TRESP, RESTOT, DT))
893 TNUPT = MAX(0.,INTGRL(TNUPT, NUPT, DT))
894 LITNT = MAX(1e-3,INTGRL(LITNT, LITN, DT))
896 ! Write updated Gecros variables into STATE_GECROS array
897 STATE_GECROS(1) = DS
898 STATE_GECROS(2) = CTDU
899 STATE_GECROS(3) = CVDU
900 STATE_GECROS(4) = CLV
901 STATE_GECROS(5) = CLVD
902 STATE_GECROS(6) = CSST
903 STATE_GECROS(7) = CSO
904 STATE_GECROS(8) = CSRT
905 STATE_GECROS(9) = CRTD
906 STATE_GECROS(10) = CLVDS
907 STATE_GECROS(11) = NRT
908 STATE_GECROS(12) = NST
909 STATE_GECROS(13) = NLV
910 STATE_GECROS(14) = NSO
911 STATE_GECROS(15) = TNLV
912 STATE_GECROS(16) = NLVD
913 STATE_GECROS(17) = NRTD
914 STATE_GECROS(18) = CRVS
915 STATE_GECROS(19) = CRVR
916 STATE_GECROS(20) = NREOE
917 STATE_GECROS(21) = NREOF
918 STATE_GECROS(22) = DCDSR
919 STATE_GECROS(23) = DCDTR
920 STATE_GECROS(24) = SLNB
921 STATE_GECROS(25) = LAIC
922 STATE_GECROS(26) = RMUL
923 STATE_GECROS(27) = NDEMP
924 STATE_GECROS(28) = NSUPP
925 STATE_GECROS(29) = NFIXT
926 STATE_GECROS(30) = NFIXR
927 STATE_GECROS(31) = DCDTP
928 STATE_GECROS(32) = HT
929 STATE_GECROS(33) = ROOTD
930 STATE_GECROS(34) = TPCAN
931 STATE_GECROS(35) = TRESP
932 STATE_GECROS(36) = TNUPT
933 STATE_GECROS(37) = LITNT
934 STATE_GECROS(45) = WSO
935 STATE_GECROS(46) = WSTRAW
936 STATE_GECROS(47) = GrainNC
937 STATE_GECROS(48) = StrawNC
938 STATE_GECROS(49) = LAI
939 STATE_GECROS(50) = TLAI
941 ! Used for debugging
942 ! if(debugging) then
943 ! if (nowdate(9:12).eq.'1200') then
944 ! write(*,*) nowdate, RSD
945 ! write(*,*) "DS: ", DS
946 ! write(*,*) "CTDU: ",CTDU
947 ! write(*,*) "CVDU: ",CVDU
948 ! write(*,*) "Carbon in living leaves CLV: ",CLV
949 ! write(*,*) "Carbon in dead leaves CLVD: ",CLVD
950 ! write(*,*) "Carbon in structural stems CSST: ",CSST
951 ! write(*,*) "Carbon in storage organs CSO: ",CSO
952 ! write(*,*) "Carbon in structural roots CSRT: ",CSRT
953 ! write(*,*) "Carbon in living roots CRT: ",CRT
954 ! write(*,*) "Carbon in dead roots CRTD: ",CRTD
955 ! write(*,*) "Rooting depth ROOTD: ",ROOTD
956 ! write(*,*) "Carbon litter on soil CLVDS: ",CLVDS
957 ! write(*,*) "Carbon in stem reserves CRVS: ",CRVS
958 ! write(*,*) "Carbon in root reserves CRVR: ",CRVR
959 ! write(*,*) "Total leaf area index TLAI: ",TLAI
960 ! write(*,*) "Dead leaf area index DLAI: ",DLAI
961 ! write(*,*) "Living leaf area index GLAI: ",LAI, LAIC, LAIN
962 ! write(*,*) "Specific leaf area index SLA: ",SLA
963 ! write(*,*) "Carbon balance check CCHK: ",CCHK
964 ! write(*,*) "Nitrogen balance check NCHK: ",NCHK
965 ! read(*,*)
966 ! endif
967 ! endif
969 ! ----------------------------------------------------------------------
970 END SUBROUTINE gecros
971 ! ----------------------------------------------------------------------
973 !******************** SUBROUTINES FOR CROP SIMULATION *******************
974 !*----------------------------------------------------------------------*
975 !* SUBROUTINE TUNIT *
976 !* Purpose: This subroutine calculates the daily amount of thermal day *
977 !* *
978 !* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) *
979 !* *
980 !* name type meaning units class *
981 !* ---- ---- ------- ----- ----- *
982 !* DS R4 Development stage - I *
983 !* TMAX R4 Daily maximum temperature oC I *
984 !* TMIN R4 Daily minimum temperature oC I *
985 !* DIF R4 Daytime plant-air temperature differential oC I *
986 !* DAYL R4 Astronomic daylength (base = 0 degrees) h I *
987 !* TBD R4 Base temperature for phenology oC I *
988 !* TOD R4 Optimum temperature for phenology oC I *
989 !* TCD R4 Ceiling temperature for phenology oC I *
990 !* TSEN R4 Curvature for temperature response - I *
991 !* TDU R4 Rate of thermal days increase d s-1 O *
992 !* VDU R4 Rate of vernalization days increase d s-1 O *
993 !*----------------------------------------------------------------------*
994 SUBROUTINE TUNIT(CROP,DS,TAIR,DIF,DAYL,TBD,TOD,TCD,TBDV,TODV,TCDV,TSEN,TDU,VDU)
995 IMPLICIT REAL (A-Z)
996 INTEGER I
998 !*---assuming development rate at supra-optimum temperatures during
999 !* the reproductive phase equals that at the optimum temperature
1000 IF (DS.GT.1.) THEN
1001 TAIR = MIN (TAIR,TOD)
1002 ELSE
1003 TAIR = TAIR
1004 ENDIF
1006 !*---vernalization response (Lenz 2007, p. 26)
1007 !*---Ingwersen 29.6.2010
1008 IF (TAIR.LT.TBDV .OR. TAIR.GT.TCDV) THEN
1009 TUV = 0.
1010 ELSE
1011 TUV = (((TCDV-TAIR)/(TCDV-TODV))*((TAIR-TBDV)/(TODV-TBDV))**((TODV-TBDV)/(TCDV-TODV)))**TSEN
1012 ENDIF
1014 !*---instantaneous thermal unit based on bell-shaped temperature response
1015 IF (TAIR.LT.TBD .OR. TAIR.GT.TCD) THEN
1016 TU = 0.
1017 ELSE
1018 TU = (((TCD-TAIR)/(TCD-TOD))*((TAIR-TBD)/(TOD-TBD))**((TOD-TBD)/(TCD-TOD)))**TSEN
1019 ENDIF
1021 !*---daily thermal unit
1022 !*** ji, 16.6.15, ecc/lcc switch
1023 IF (CROP==1) THEN
1024 TDU = TU/TCP
1025 ELSE
1026 TDU = INSW(DS-2.,TU/TCP,0.)
1027 ENDIF
1028 VDU = TUV/TCP
1030 RETURN
1031 END SUBROUTINE TUNIT
1034 !*----------------------------------------------------------------------*
1035 !* SUBROUTINE PHENO *
1036 !* Purpose: This subroutine calculates phenological development rate. *
1037 !* *
1038 !* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) *
1039 !* *
1040 !* name type meaning units class *
1041 !* ---- ---- ------- ----- ----- *
1042 !* DS R4 Development stage - I *
1043 !* SLP R4 Crop type(1. for short-day,-1. for long-day) - I *
1044 !* DDLP R4 Daylength for photoperiodism h I *
1045 !* SPSP R4 DS for start of photoperiod-sensitive phase - I *
1046 !* EPSP R4 DS for end of photoperiod-sensitive phase - I *
1047 !* PSEN R4 Photoperiod sensitivity (+ for SD, - for LD) h-1 I *
1048 !* MTDV R4 Minimum thermal days for vegetative phase d I *
1049 !* MTDR R4 Minimum thermal days for reproductive phase d I *
1050 !* TDU R4 Thermal unit - I *
1051 !* CVDU R4 Cumulative vernalization unit (ji, 29.6) - I *
1052 !* DVR R4 Development rate s-1 O *
1053 !*----------------------------------------------------------------------*
1054 SUBROUTINE PHENO (CROP,DS,SLP,DDLP,SPSP,EPSP,PSEN,MTDV,MTDR,TDU,CVDU,DVR)
1055 IMPLICIT REAL (A-Z)
1057 !*---determining if it is for short-day or long-day crop
1058 IF (SLP.LT.0.) THEN
1059 MOP = 18. !minimum optimum photoperiod for long-day crop
1060 DLP = MIN(MOP,DDLP)
1061 ELSE
1062 MOP = 11. !maximum optimum photoperiod for short-day crop
1063 DLP = MAX(MOP,DDLP)
1064 ENDIF
1066 !*---effect of photoperiod on development rate
1067 IF (DS.LT.SPSP .OR. DS.GT.EPSP) THEN
1068 EFP = 1.
1069 ELSE
1070 EFP = MAX(0., 1.-PSEN*(DLP-MOP))
1071 ENDIF
1073 !*---effect of vernalization (ji, 21.6.10)
1074 !*** ji, 16.6.15, ecc/lcc switch
1075 IF (CROP==1) THEN
1076 EFV = CVDU**5./(22.5**5. + CVDU**5.)
1077 ELSE
1078 EFV = 1.0
1079 ENDIF
1081 !*---development rate of vegetative and reproductive phases
1082 !*---extended for vernalization according to Lenz (2007); ji: 21.6.2010
1083 IF (DS.LE.0.4) THEN
1084 DVR = 1./MTDV*TDU*EFP*EFV
1085 ENDIF
1087 IF (DS.GT.0.4 .AND. DS.LE.1.0) THEN
1088 DVR = 1./MTDV*TDU*EFP
1089 ENDIF
1091 IF (DS.GT.1.0) THEN
1092 DVR = 1./MTDR*TDU
1093 ENDIF
1095 RETURN
1096 END SUBROUTINE PHENO
1099 !*----------------------------------------------------------------------*
1100 !* SUBROUTINE RNACC *
1101 !* Purpose: This subroutine calculates rate of N accumulation in organs*
1102 !* *
1103 !* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) *
1104 !* *
1105 !* name type meaning units class *
1106 !* ---- ---- ------- ----- ----- *
1107 !* FNSH R4 Fraction of new N partitioned to shoot - I *
1108 !* NUPT R4 Nitrogen uptake at a time step gN/m2/s I *
1109 !* RWST R4 Rate of stem weight g/m2/s I *
1110 !* STEMNC R4 Nitrogen concentration in stem gN/g I *
1111 !* LNCMIN R4 Minimum N concentration in leaf gN/g I *
1112 !* RNCMIN R4 Minimum N concentration in root gN/g I *
1113 !* LNC R4 Nitrogen concentration in leaf gN/g I *
1114 !* RNC R4 Nitrogen concentration in root gN/g I *
1115 !* NLV R4 Canopy (green)leaf N content gN/m2 I *
1116 !* NRT R4 (living)root N content gN/m2 I *
1117 !* WLV R4 Canopy (green)leaf weight g/m2 I *
1118 !* WRT R4 (living)Root weight g/m2 I *
1119 !* DELT R4 Time step of simulation s I *
1120 !* CB R4 Factor for initial N concent. of seed-fill - I *
1121 !* CX R4 Factor for final N concent. of seed-fill - I *
1122 !* TM R4 DS when transition from CB to CX is fastest - I *
1123 !* DS R4 Development stage - I *
1124 !* SEEDNC R4 Standard seed N concentration gN/g I *
1125 !* RWSO R4 growth rate of seed g/m2/s I *
1126 !* LNLV R4 Loss rate of NLV due to senescence gN/m2/s I *
1127 !* LNRT R4 Loss rate of NRT due to senescence gN/m2/s I *
1128 !* RNRT R4 rate of N accumulation in root gN/m2/s O *
1129 !* RNST R4 rate of N accumulation in stem gN/m2/s O *
1130 !* RNLV R4 rate of N accumulation in leaf gN/m2/s O *
1131 !* RTNLV R4 Positive value of RNLV gN/m2/s O *
1132 !* RNSO R4 rate of N accumulation in seed(storage organ)gN/m2/s O *
1133 !*----------------------------------------------------------------------*
1134 SUBROUTINE RNACC (FNSH,NUPT,RWST,STEMNC,LNCMIN,RNCMIN,LNC,RNC, &
1135 NLV,NRT,WLV,WRT,DELT,CB,CX,TM,DS,SEEDNC, &
1136 RWSO,LNLV,LNRT, RNRT,RNST,RNLV,RTNLV,RNSO)
1137 IMPLICIT REAL (A-Z)
1139 !*---amount of N partitioned to shoot
1140 NSHN = FNSH * NUPT
1142 !*---leaf N (NLVA) or root N (NRTA) available for remobilization within a day
1143 NLVA = INSW(LNCMIN-LNC, NLV-WLV*LNCMIN, 0.)/TCP
1144 NRTA = INSW(RNCMIN-RNC, NRT-WRT*RNCMIN, 0.)/TCP
1146 NTA = NLVA + NRTA
1148 !*---rate of N accumulation in stem
1149 RNST = RWST * INSW(-NTA,STEMNC,0.)
1151 !*---expected N dynamics during seed(storage organ) filling
1152 CDS = CB+(CX-CB)*(4.-TM-DS)/(2.-TM)*(DS-1.)**(1./(2.-TM))
1153 ENSNC = LIMIT(CB,CX,CDS) * SEEDNC
1155 !*---rate of N accumulation in seed
1156 NGS = NSHN - RNST - ENSNC*RWSO
1157 NONC = MAX(0.,INSW(NTA+NGS,(NTA+NSHN-RNST)/NOTNUL(RWSO),ENSNC))
1158 RNSO = RWSO*NONC
1160 !*---rate of N accumulation in leaf!
1162 NLVN = INSW(NTA+NGS,-NLVA-LNLV,-NLVA/NOTNUL(NTA)*(-NGS)-LNLV)
1163 GNLV = INSW(NGS, NLVN, NSHN-RNST-RNSO-LNLV)
1164 RNLV = MAX (-NLV+1.E-7, GNLV)
1165 RTNLV = MAX(0., RNLV)
1167 !*---rate of N accumulation in root
1168 NRTN = INSW(NTA+NGS, NUPT-NSHN-NRTA-LNRT, &
1169 NUPT-NSHN-NRTA/NOTNUL(NTA)*(-NGS)-LNRT)
1170 GNRT = INSW(NGS, NRTN, NUPT-NSHN-LNRT)
1171 RNRT = MAX (-NRT+5.E-8, GNRT)
1173 RETURN
1174 END SUBROUTINE RNACC
1177 !*----------------------------------------------------------------------*
1178 !* SUBROUTINE RLAIC *
1179 !* Purpose: This subroutine calculates the daily increase of leaf *
1180 ! area index (m2 leaf/m2 ground/day). *
1181 ! *
1182 !* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) *
1183 !* *
1184 !* name type meaning units class *
1185 !* ---- ---- ------- ----- ----- *
1186 !* DS R4 Development stage - I *
1187 !* SLA0 R4 Specific leaf area constant m2 g-1 I *
1188 !* RWLV R4 Rate of increment in leaf weight g m-2 s-1 I *
1189 !* LAI R4 Leaf area index m2 m-2 I *
1190 !* KN R4 Leaf nitrogen extinction coefficient m2 m-2 I *
1191 !* NLV R4 Total leaf nitrogen content in a canopy g m-2 I *
1192 !* RNLV R4 Rate of increment in NLV g m-2 s-1 I *
1193 !* SLNB R4 Nitrogen content of bottom leaves g m-2 I *
1194 !* RSLNB R4 Rate of increment in SLNB g m-2 s-1 I *
1195 !* RLAI R4 Rate of increment in leaf area index m2 m-2s-1 O *
1196 !*----------------------------------------------------------------------*
1197 SUBROUTINE RLAIC(CROP,DS,SLA0,RWLV,LAI,KN,NLV,RNLV,SLNB,RSLNB, RLAI)
1198 IMPLICIT REAL (A-Z)
1200 SLNB = MAX(1E-2, SLNB)
1201 !*---rate of LAI driven by carbon supply
1202 RLAI = INSW(RWLV, MAX(-LAI+1.E-5,SLA0*RWLV), SLA0*RWLV)
1204 !*---rate of LAI driven by nitrogen during juvenile phase
1205 !*** ji, 16.6.15, ecc/lcc switch
1206 IF ((CROP==2) .AND. (LAI.LT.1.5) .AND. (DS.LT.0.75)) THEN
1207 RLAI = MAX(0.,(SLNB*RNLV-NLV*RSLNB)/SLNB/(SLNB+KN*NLV))
1208 ENDIF
1211 RETURN
1212 END SUBROUTINE RLAIC
1215 !*----------------------------------------------------------------------*
1216 !* SUBROUTINE BETAF *
1217 !* Purpose: This subroutine calculates the dynamics of expected growth *
1218 !* of sinks, based on the beta sigmoid growth equation *
1219 !* *
1220 !* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) *
1221 !* *
1222 !* name type meaning units class *
1223 !* ---- ---- ------- ----- ----- *
1224 !* DVR R4 Development rate s-1 I *
1225 !* TE R4 Stage at which sink growth stops - I *
1226 !* TX R4 Stage at which sink growth rate is maximal - I *
1227 !* TI R4 Stage of a day considered - I *
1228 !* FD R4 Relative expected growth of a sink at a day s-1 O *
1229 !*----------------------------------------------------------------------*
1230 SUBROUTINE BETAF(DVRX,TE,TX,TI, FD)
1232 REAL, INTENT(IN) :: DVRX, TE, TX, TI
1233 REAL, INTENT(OUT) :: FD
1235 FD = DVRX*(2.*TE-TX)*(TE-TI)/TE/(TE-TX)**2*(TI/TE)**(TX/(TE-TX))
1236 !FD = DVRX*(TE-TI)/(TE-TX)*(TI/TX)**(TX/(TE-TX)) !Eq. a 1a Yin et al 2003
1238 END SUBROUTINE BETAF
1241 !*----------------------------------------------------------------------*
1242 !* SUBROUTINE SINKG *
1243 !* Purpose: This subroutine calculates carbon demand for sink growth. *
1244 !* *
1245 !* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) *
1246 !* *
1247 !* name type meaning units class *
1248 !* ---- ---- ------- ----- ----- *
1249 !* DS R4 Development stage - I *
1250 !* SSG R4 Stage at which sink growth starts - I *
1251 !* TOTC R4 Total carbon in a sink at end of its growth g C/m2 I *
1252 !* YG R4 Growth efficiency g C/g C I *
1253 !* FD R4 Relative expected growth of a sink at a day s-1 I *
1254 !* DCDR R4 Shortfall of C demand in previous days g C/m2 I *
1255 !* DCS R4 Daily C supply for sink growth g C/m2/s I *
1256 !* DELT R4 Time step of integration s I *
1257 !* DCDC R4 C demand of the current day g C/m2/s O *
1258 !* DCD R4 Daily C demand for sink growth g C/m2/s O *
1259 !* FLWC R4 Flow of current assimilated C to sink g C/m2/s O *
1260 !*----------------------------------------------------------------------*
1261 SUBROUTINE SINKG(DS,SSG,TOTC,YG,FD,DCDR,DCS,DELT,DCDC,DCD,FLWC)
1262 IMPLICIT REAL (A-Z)
1264 !*---expected demand for C of the current time step
1265 DCDC = INSW (DS-SSG, 0., TOTC/YG*FD)
1267 !*---total demand for C at the time step considered
1268 DCD = DCDC + MAX(0.,DCDR)/DELT
1270 !*---flow of current assimilated carbon to sink
1271 FLWC = MIN(DCD, DCS)
1273 RETURN
1274 END SUBROUTINE SINKG
1277 !*----------------------------------------------------------------------*
1278 !* SUBROUTINE ASTRO (from the SUCROS model) *
1279 !* Purpose: This subroutine calculates astronomic daylength, *
1280 !* diurnal radiation characteristics such as the daily *
1281 !* integral of sine of solar elevation and solar constant. *
1282 !* *
1283 !* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) *
1284 !* name type meaning units class *
1285 !* ---- ---- ------- ----- ----- *
1286 !* DOY R4 Daynumber (Jan 1st = 1) - I *
1287 !* LAT R4 Latitude of the site degree I *
1288 !* INSP R4 Inclination of sun angle for computing DDLP degree I *
1289 !* SC R4 Solar constant J m-2 s-1 O *
1290 !* SINLD R4 Seasonal offset of sine of solar height - O *
1291 !* COSLD R4 Amplitude of sine of solar height - O *
1292 !* DAYL R4 Astronomic daylength (base = 0 degrees) h O *
1293 !* DDLP R4 Photoperiodic daylength h O *
1294 !* DSINBE R4 Daily total of effective solar height s d-1 O *
1295 !* *
1296 !* FATAL ERROR CHECKS (execution terminated, message) *
1297 !* condition: LAT > 67, LAT < -67 *
1298 !* *
1299 !* FILE usage : none *
1300 !*----------------------------------------------------------------------*
1301 SUBROUTINE ASTRO (DOY,LAT,INSP,SC,SINLD,COSLD,DAYL,DDLP,DSINBE)
1302 IMPLICIT REAL (A-Z)
1304 !*---PI and conversion factor from degrees to radians
1305 PI = 3.141592654
1306 RAD = PI/180.
1307 !*---check on input range of parameters
1308 IF (LAT.GT.67.) STOP 'ERROR IN ASTRO: LAT> 67'
1309 IF (LAT.LT.-67.) STOP 'ERROR IN ASTRO: LAT>-67'
1311 !*---declination of the sun as function of daynumber (DOY)
1312 DEC = -ASIN (SIN (23.45*RAD)*COS (2.*PI*(DOY+10.)/365.))
1314 !*---SINLD, COSLD and AOB are intermediate variables
1315 SINLD = SIN (RAD*LAT)*SIN (DEC)
1316 COSLD = COS (RAD*LAT)*COS (DEC)
1317 AOB = SINLD/COSLD
1319 !*---daylength (DAYL)
1320 DAYL = 12.0*(1.+2.*ASIN (AOB)/PI)
1321 DDLP = 12.0*(1.+2.*ASIN((-SIN(INSP*RAD)+SINLD)/COSLD)/PI)
1323 DSINB = 3600.*(DAYL*SINLD+24.*COSLD*SQRT (1.-AOB*AOB)/PI)
1324 DSINBE = 3600.*(DAYL*(SINLD+0.4*(SINLD*SINLD+COSLD*COSLD*0.5))+ &
1325 12.0*COSLD*(2.0+3.0*0.4*SINLD)*SQRT (1.-AOB*AOB)/PI)
1327 !*---solar constant (SC)
1328 SC = 1367.*(1.+0.033*COS(2.*PI*(DOY-10.)/365.))
1330 RETURN
1331 END SUBROUTINE ASTRO
1334 !*----------------------------------------------------------------------*
1335 !* SUBROUTINE TOTPT *
1336 !* Purpose: This subroutine calculates daily total gross photosynthesis*
1337 !* and transpiration by performing a Gaussian integration *
1338 !* over time. At five different times of the day, temperature *
1339 !* and radiation are computed to determine assimilation *
1340 !* and transpiration whereafter integration takes place. *
1341 !* *
1342 !* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) *
1343 !* *
1344 !* name type meaning units class *
1345 !* ---- ---- ------- ----- ----- *
1346 !* HOD R4 Hour of the day h I *
1347 !* SC R4 Solar constant J m-2 s-1 I *
1348 !* SINLD R4 Seasonal offset of sine of solar height - I *
1349 !* COSLD R4 Amplitude of sine of solar height - I *
1350 !* DAYL R4 Astronomic daylength (base = 0 degrees) h I *
1351 !* DSINBE R4 Daily total of effective solar height s d-1 I *
1352 !* RSD R4 Global radiation J m-2 s-1 I *
1353 !* TAIR R4 Air temperature oC I *
1354 !* DVP R4 Vapour pressure kPa I *
1355 !* WNM R4 daily average wind speed (>=0.1 m/s) m s-1 I *
1356 !* C3C4 R4 Crop type (=1 for C3, -1 for C4 crops) - I *
1357 !* LAI R4 (green)Leaf area index m2 m-2 I *
1358 !* TLAI R4 Total Leaf area index m2 m-2 I *
1359 !* HT R4 Plant height m I *
1360 !* LWIDTH R4 Leaf width m I *
1361 !* RD R4 Rooting depth cm I *
1362 !* SD1 R4 Depth of evaporative upper soil layer cm I *
1363 !* RSS R4 Soil resistance,equivalent to leaf stomata s m-1 I *
1364 !* BLD R4 Leaf angle from horizontal degree I *
1365 !* KN R4 Leaf nitrogen extinction coefficient m2 m-2 I *
1366 !* KW R4 Windspeed extinction coefficient in canopy m2 m-2 I *
1367 !* SLN R4 Average leaf nitrogen content in canopy g m-2 I *
1368 !* SLNT R4 Top-leaf nitrogen content g m-2 I *
1369 !* SLNN R4 Value of SLNT with small plant-N increment g m-2 I *
1370 !* SLNMIN R4 Minimum or base SLNT for photosynthesis g m-2 I *
1371 !* DWSUP R4 Daily water supply for evapotranspiration mm s-1 I *
1372 !* CO2A R4 Ambient CO2 concentration ml m-3 I *
1373 !* LS R4 Lodging severity - I *
1374 !* EAJMAX R4 Energy of activation for Jmax J mol-1 I *
1375 !* XVN R4 Slope of linearity between Vcmax & leaf N umol/g/s I *
1376 !* XJN R4 Slope of linearity between Jmax & leaf N umol/g/s I *
1377 !* THETA R4 Convexity for light response of e-transport - I *
1378 !* WCUL R4 Water content of the upper soil layer m3 m-3 I *
1379 !* FVPD R4 Slope for linear effect of VPD on Ci/Ca (kPa)-1 I *
1380 !* PPCAN R4 Potential canopy CO2 assimilation g m-2 s-1 O *
1381 !* APCANS R4 Actual standing-canopy CO2 assimilation g m-2 s-1 O *
1382 !* APCANN R4 APCANS with small plant-N increment g m-2 s-1 O *
1383 !* APCAN R4 Actual canopy CO2 assimilation g m-2 s-1 O *
1384 !* PTCAN R4 Potential canopy transpiration mm s-1 O *
1385 !* ATCAN R4 Actual canopy transpiration mm s-1 O *
1386 !* PESOIL R4 Potential soil evaporation mm s-1 O *
1387 !* AESOIL R4 Actual soil evaporation mm s-1 O *
1388 !* DIFS R4 Soil-air temp. difference oC O *
1389 !* DIFSU R4 Sunlit leaf-air temp. diff. oC O *
1390 !* DIFSH R4 Shaded leaf-air temp. diff. oC O *
1391 !* DAPAR R4 PAR absorbed by crop canopy J m-2 s-1 O *
1392 !* RCAN R4 Canopy resistance s/m 0 *
1393 !*----------------------------------------------------------------------*
1394 SUBROUTINE TOTPT(HOD,DS, SC,SINLD,COSLD,DAYL,DSINBE,RSD,TAIR, DVP, &
1395 WNM,C3C4,LAI,TLAI,HT,LWIDTH,RD,SD1,RSS,BLD,KN,KW, &
1396 SLN,SLNT,SLNN,SLNMIN,FB,DWSUP,CO2A,LS,EAJMAX, &
1397 XVN,XJN,THETA,WCUL,FVPD,RB,RT,RTS,PPCAN,APCANS, &
1398 APCANN,APCAN,PTCAN,ATCAN,PESOIL,AESOIL,DIFS,DIFSU,DIFSH,DAPAR, &
1399 RCAN, ATRJC, ATRJS, FSR, FRSU, FRSH, ARSWSU, ARSWSH)
1401 REAL, INTENT(IN) :: HOD,DS,SC,SINLD,COSLD,DAYL,DSINBE,RSD,TAIR, DVP, &
1402 WNM,C3C4,LAI,TLAI,HT,LWIDTH,RD,SD1,RSS,BLD,KN,KW, &
1403 SLN,SLNT,SLNN,SLNMIN,FB,DWSUP,CO2A,LS,EAJMAX, &
1404 XVN,XJN,THETA,WCUL,FVPD,RT,RTS !RB
1406 REAL, INTENT(INOUT) :: PPCAN,APCANS,APCANN,APCAN,PTCAN,ATCAN,PESOIL, &
1407 AESOIL,DIFS,DIFSU,DIFSH,DAPAR,RCAN,ATRJC,ATRJS,FSR,FRSU,FRSH,ARSWSU,ARSWSH
1408 REAL :: ACO2I,ADIFS,ADIFSH,ADIFSU,ANIRSH,ANIRSU,ANRAD,APAR,APARSH
1409 REAL :: APARSU,ASVP,ATMTR,ATRJSH,ATRJSU,ATSH,ATSU,AV_Albedo
1410 REAL :: AV_RSWSH,AV_RSWSU,BL,CUMRSD,DATRJC,DATRJS,DNRAD,FRDF
1411 REAL :: GBHC,GBHLF,GBHSH,GBHSU,IAE,IAP,IAPL,IAPN,IAPNN,IAPS,IAT,IPE
1412 REAL :: IPH,IPHSOIL,IPP,IPPL,IPT,IRDL,KB,KBPNIR,KBPPAR,KDPNIR,KDPPAR
1413 REAL :: NIR,NIRDF,NIRDR,NPSH,NPSHN,NPSU,NPSUN,NRADS,NRADSH,PANSH,PANSU
1414 REAL :: PAR,NRADSU,PARDF,PARDR,PASSH,PASSU,PCBNIR,PCBPAR,PCDNIR,PCDPAR
1415 REAL :: PHCAN,PHSH,PHSOIL,PHSU,PI,PLFSH,PLFSU,PSNIR,PSPAR,PT1,PTSH,PTSU
1416 REAL :: RBHS,RBHSH,RBHSU,RBWS,RBWSH,RBWSU,RSWSH,RSWSU,SCPNIR
1417 REAL :: SCPPAR,SINB,SLOPSH,SLOPSU,WND,WSUP,WSUP1,Albedo,RB, DSsw
1418 REAL :: TLEAFSH, TLEAFSU
1420 PI = 3.141592654
1422 !*---output-variables set to 0. and five different times of a day(HOUR)
1423 PPCAN = 0.
1424 APCANS = 0.
1425 APCANN = 0.
1426 APCAN = 0.
1427 PTCAN = 0.
1428 PHCAN = 0.
1429 ATCAN = 0.
1430 PESOIL = 0.
1431 AESOIL = 0.
1432 PHSOIL = 0.
1433 DIFS = 0.
1434 DIFSU = 0.
1435 DIFSH = 0.
1436 DAPAR = 0.
1437 DNRAD = 0.
1438 DATRJC = 0.
1439 DATRJS = 0.
1440 AV_RSWSU = 0.
1441 AV_RSWSH = 0.
1442 AV_Albedo = 0.
1443 CUMRSD = 0.
1445 !*---sine of solar elevation
1446 SINB = MAX (.01, SINLD+COSLD*COS(2.*PI*(HOD-12.)/24.))
1448 !*---daytime course of water supply
1449 WSUP = DWSUP
1450 WSUP1 = WSUP*SD1/RD
1452 !*---daytime course of wind speed
1453 WND = WNM !m s-1
1455 !*---total incoming PAR and NIR
1456 PAR = (1.-FB)*0.5*RSD !J m-2 s-1
1457 NIR = (1.-FB)*0.5*RSD !J m-2 s-1
1459 !*---diffuse light fraction (FRDF) from atmospheric transmission (ATMTR)
1460 ATMTR = PAR/(0.5*SC*SINB) ! unitless, fraction
1462 IF (ATMTR.LE.0.22) THEN
1463 FRDF = 1.
1464 ELSE IF (ATMTR.GT.0.22 .AND. ATMTR.LE.0.35) THEN
1465 FRDF = 1.-6.4*(ATMTR-0.22)**2
1466 ELSE
1467 FRDF = 1.47-1.66*ATMTR
1468 ENDIF
1470 FRDF = MAX (FRDF, 0.15+0.85*(1.-EXP (-0.1/SINB)))
1472 !*---incoming diffuse PAR (PARDF) and direct PAR (PARDR)
1473 PARDF = PAR * FRDF !J m-2 s-1
1474 PARDR = PAR - PARDF !J m-2 s-1
1476 !*---incoming diffuse NIR (NIRDF) and direct NIR (NIRDR)
1477 NIRDF = NIR * FRDF !J m-2 s-1
1478 NIRDR = NIR - NIRDF !J m-2 s-1
1480 !*---extinction and reflection coefficients
1481 BL = BLD*PI/180. !leaf angle, conversion to radians
1482 CALL KBEAM (SINB,BL,KB)
1484 SCPPAR = 0.2 !leaf scattering coefficient for PAR
1485 SCPNIR = 0.8 !leaf scattering coefficient for NIR
1486 CALL KDIFF (TLAI,BL,SCPPAR, KDPPAR)
1487 CALL KDIFF (TLAI,BL,SCPNIR, KDPNIR)
1489 CALL REFL (SCPPAR,KB, KBPPAR,PCBPAR)
1490 CALL REFL (SCPNIR,KB, KBPNIR,PCBNIR)
1492 PCDPAR = 0.057 !canopy diffuse PAR reflection coefficient
1493 PCDNIR = 0.389 !canopy diffuse NIR reflection coefficient
1495 !*---fraction of sunlit and shaded components in canopy
1496 !* ji 4.7.11 LAI -> TLAI
1497 FRSU = 1./KB/TLAI*(1.-EXP(-KB*TLAI))
1498 FRSH = 1.-FRSU
1500 !*---leaf boundary layer conductance for canopy, sunlit and shaded leaves !m s-1
1501 GBHLF = 0.01*SQRT(WND/LWIDTH)
1502 RBHSU = RB/(FRSU*LAI) !boundary layer resistance to heat,sunlit part !s m-1
1503 RBWSU = RB/(FRSU*LAI) !boundary layer resistance to H2O, sunlit part
1504 RBHSH = RB/(FRSH*LAI) !boundary layer resistance to heat,shaded part
1505 RBWSH = RB/(FRSH*LAI) !boundary layer resistance to H2O, shaded part
1506 RCAN = WNM*10.
1508 !*---boundary layer resistance for soil !s m-1
1509 RBHS = 172.*SQRT(0.05/MAX(0.1,WND*EXP(-KW*TLAI)))
1510 RBWS = 0.93*RBHS
1512 !*---photosynthetically active nitrogen for sunlit and shaded leaves
1513 CALL PAN (SLNT,SLNMIN,LAI,KN,KB, NPSU,NPSH)
1514 CALL PAN (SLNN,SLNMIN,LAI,KN,KB, NPSUN,NPSHN)
1516 !*---absorbed PAR and NIR by sunlit leaves and shaded leaves
1517 !* ji 4.7.11 LAI->TLAI
1518 !* Ansonsten fuehrt es dazu, dass nach dem Abreifen die Albedo auf >0.7 ansteigt und vom Bestand viel
1519 !* zu wenig Energie mehr absorbiert wird
1520 CALL LIGAB (SCPPAR,KB,KBPPAR,KDPPAR,PCBPAR,PCDPAR,PARDR,PARDF,TLAI,APARSU,APARSH)
1521 CALL LIGAB (SCPNIR,KB,KBPNIR,KDPNIR,PCBNIR,PCDNIR,NIRDR,NIRDF,TLAI,ANIRSU,ANIRSH)
1522 APAR = APARSU+APARSH !J m-2 s-1
1524 !*---absorbed total radiation (PAR+NIR) by sunlit and shaded leaves
1525 ATRJSU = APARSU+ANIRSU !J m-2 s-1
1526 ATRJSH = APARSH+ANIRSH !J m-2 s-1
1527 ATRJC = ATRJSH + ATRJSU !J m-2 s-1
1529 !*---absorbed total radiation (PAR+NIR) by soil
1530 PSPAR = 0.1 !soil PAR reflection
1531 PSNIR = INSW(WCUL-0.5, 0.52-0.68*WCUL, 0.18) !soil NIR reflection
1532 ATRJS=(1.-PSPAR)*(PARDR*EXP(-KBPPAR*TLAI)+PARDF*EXP(-KDPPAR*TLAI)) &
1533 +(1.-PSNIR)*(NIRDR*EXP(-KBPNIR*TLAI)+NIRDF*EXP(-KDPNIR*TLAI))
1535 !* ji 4.7.11 Berechnung der Albedo eingefuegt. Nachts wird die Albedo auf einen Durchschnittswert gesetzt (.2)
1536 FSR = RSD-ATRJC-ATRJS
1538 if (RSD.gt.0) then
1539 Albedo = (RSD-ATRJC-ATRJS)/RSD
1540 else
1541 Albedo = .2
1542 endif
1544 !*---instantaneous potential photosynthesis and transpiration
1545 CALL PPHTR(FRSU,TAIR,DVP,CO2A,C3C4,FVPD,APARSU,NPSU,RBWSU,RBHSU, &
1546 RT*FRSU,ATRJSU,ATMTR,EAJMAX,XVN,XJN,THETA,PLFSU, &
1547 PTSU,PHSU,RSWSU,NRADSU,SLOPSU)
1549 CALL PPHTR(FRSH,TAIR,DVP,CO2A,C3C4,FVPD,APARSH,NPSH,RBWSH,RBHSH, &
1550 RT*FRSH,ATRJSH,ATMTR,EAJMAX,XVN,XJN,THETA,PLFSH, &
1551 PTSH,PHSH,RSWSH,NRADSH,SLOPSH)
1553 IPP = PLFSU + PLFSH !gCO2/m2/s
1554 IPT = PTSU + PTSH !mm s-1
1555 IPH = PHSU + PHSH !J m-2 s-1
1556 ANRAD = NRADSU + NRADSH !J m-2 s-1
1558 !*--- PT1: Potential transpiration using water from the upper evaporative layer (mm s-1)
1559 !*--- SD1: thickness of upper evaporative layer (cm); default value: 5 cm
1560 !*--- RD: Rooting depth (cm)
1561 PT1 = IPT * SD1/RD
1563 !*---instantaneous potential soil evaporation
1564 CALL PEVAP (TAIR,DVP,RSS,RTS,RBWS,RBHS,ATRJS,ATMTR, &
1565 PT1,WSUP1,IPE,IPHSOIL,NRADS)
1567 !*---instantaneous actual soil evaporation, actual canopy
1568 !* transpiration and photosynthesis
1569 IAE = MIN (IPE,IPE/(PT1+IPE)*WSUP1) !actual soil evaporation mm s-1
1570 IAT = MIN (IPT,PT1/(PT1+IPE)*WSUP1+WSUP-WSUP1) !actual transpiration mm s-1
1571 ATSU = PTSU/IPT*IAT !actual transpiration of sunlit leaves mm s-1
1572 ATSH = PTSH/IPT*IAT !actual transpiration of shaded leaves mm s-1
1574 CALL DIFLA (NRADS,IAE,RBHS,RTS, ADIFS)
1576 CALL APHTR (TAIR,APARSU,DVP,CO2A,C3C4,FVPD,NRADSU,ATSU,PTSU, &
1577 RT*FRSU,RBHSU,RBWSU,RSWSU,SLOPSU,NPSU,NPSUN, &
1578 EAJMAX,XVN,XJN,THETA,PASSU,PANSU,ADIFSU,ARSWSU)
1580 CALL APHTR (TAIR,APARSH,DVP,CO2A,C3C4,FVPD,NRADSH,ATSH,PTSH, &
1581 RT*FRSH,RBHSH,RBWSH,RSWSH,SLOPSH,NPSH,NPSHN, &
1582 EAJMAX,XVN,XJN,THETA,PASSH,PANSH,ADIFSH,ARSWSH)
1584 IAPS = PASSU + PASSH !actual canopy photosynthesis gCO2 m-2 s-1
1585 IAPN = PANSU + PANSH !actual canopy photosynthesis with a small N increment gCO2 m-2 s-1
1586 IAP = IAPS
1587 IAPNN = IAPN
1589 !*---integration of assimilation and transpiration to a daily total
1590 PPCAN = IPP
1591 APCANS = IAPS
1592 APCANN = IAPNN
1593 APCAN = IAP
1594 PTCAN = IPT
1595 PHCAN = IPH
1596 ATCAN = IAT
1597 PESOIL = IPE
1598 AESOIL = IAE
1599 PHSOIL = IPHSOIL
1600 DIFS = ADIFS
1601 DIFSU = ADIFSU
1602 DIFSH = ADIFSH
1603 DAPAR = APAR
1604 DNRAD = ANRAD !net absorbed radiation by canopy
1605 DATRJC = ATRJC !absorbed radiation by canopy
1606 DATRJS = ATRJS !absorbed radiation by soil
1607 CUMRSD = RSD !incoming solar radiation
1608 AV_RSWSU = RSWSU
1609 AV_RSWSH = RSWSH
1610 AV_Albedo = Albedo
1611 TLEAFSU = TAIR + DIFSU
1612 TLEAFSH = TAIR + DIFSH
1614 RETURN
1615 END SUBROUTINE TOTPT
1618 !*----------------------------------------------------------------------*
1619 !* SUBROUTINE PPHTR *
1620 !* Purpose: This subroutine calculates potential leaf photosynthesis *
1621 !* and transpiration. *
1622 !* *
1623 !* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) *
1624 !* name type meaning units class *
1625 !* ---- ---- ------- ----- ----- *
1626 !* FRAC R4 Fraction of leaf classes (sunlit vs shaded) - I *
1627 !* TAIR R4 Air temperature oC I *
1628 !* DVP R4 Vapour pressure kPa I *
1629 !* CO2A R4 Ambient CO2 concentration ml m-3 I *
1630 !* C3C4 R4 Crop type (=1. for C3, -1 for C4 crops) - I *
1631 !* FVPD R4 Slope for linear effect of VPD on Ci/Ca (kPa)-1 I *
1632 !* PAR R4 Absorbed photosynth. active radiation J m-2 s-1 I *
1633 !* NP R4 Photosynthetically active N content g m-2 I *
1634 !* RBW R4 Leaf boundary layer resistance to water s m-1 I *
1635 !* RBH R4 Leaf boundary layer resistance to heat s m-1 I *
1636 !* RT R4 Turbulence resistance s m-1 I *
1637 !* ATRJ R4 Absorbed global radiation J m-2 s-1 I *
1638 !* ATMTR R4 Atmospheric transmissivity - I *
1639 !* EAJMAX R4 Energy of activation for Jmax J mol-1 I *
1640 !* XVN R4 Slope of linearity between Vcmax & leaf N umol/g/s I *
1641 !* XJN R4 Slope of linearity between Jmax & leaf N umol/g/s I *
1642 !* THETA R4 Convexity for light response of e-transport - I *
1643 !* PLF R4 Potential leaf photosynthesis gCO2/m2/s O *
1644 !* PT R4 Potential leaf transpiration mm s-1 O *
1645 !* PH R4 Potential leaf sensible heat flux J m-2 s-1 O *
1646 !* RSW R4 Potential stomatal resistance to water s m-1 O *
1647 !* NRADC R4 Net leaf absorbed radiation J m-2 s-1 O *
1648 !* SLOPEL R4 Slope of saturated vapour pressure curve kPa oC-1 O *
1649 !*----------------------------------------------------------------------*
1650 SUBROUTINE PPHTR(FRAC,TAIR,DVP,CO2A,C3C4,FVPD,PAR,NP,RBW,RBH,RT, &
1651 ATRJ,ATMTR,EAJMAX,XVN,XJN,THETA,PLF,PT,PH,RSW,NRADC,SLOPEL)
1652 IMPLICIT REAL (A-Z)
1654 !*---first-round calculation to determine leaf temperature
1655 CALL ICO2 (TAIR,DVP,FVPD,CO2A,C3C4, SVP,FCO2I)
1657 CALL PHOTO(C3C4,PAR,TAIR,FCO2I,NP,EAJMAX,XVN,XJN,THETA,FPLF, &
1658 FLRD)
1660 VPD = MAX (0., SVP- DVP)
1661 SLOPE = 4158.6 * SVP/(TAIR + 239.)**2
1662 CALL GCRSW(FPLF,FLRD,TAIR,CO2A,FCO2I,RBW,RT, FRSW)
1664 CALL PTRAN(FRSW,RT,RBW,RBH,ATRJ,ATMTR,FRAC,TAIR,DVP, &
1665 SLOPE, VPD, FPT, FPH, FNRADC)
1667 CALL DIFLA (FNRADC,FPT,RBH,RT, FDIF)
1669 TLEAF = TAIR + FDIF
1671 !*---second-round calculation to determine potential photosynthesis
1672 !* and transpiration
1673 CALL ICO2 (TLEAF,DVP,FVPD,CO2A,C3C4, SVPL,CO2I)
1674 CALL PHOTO (C3C4,PAR,TLEAF,CO2I,NP,EAJMAX,XVN,XJN,THETA,PLF,LRD)
1676 SLOPEL = (SVPL-SVP)/NOTNUL(TLEAF-TAIR)
1678 CALL GCRSW (PLF,LRD,TLEAF,CO2A,CO2I,RBW,RT, RSW)
1679 CALL PTRAN (RSW,RT,RBW,RBH,ATRJ,ATMTR,FRAC,TLEAF,DVP, &
1680 SLOPEL,VPD, PT, PH, NRADC)
1682 CALL DIFLA (FNRADC,FPT,RBH,RT, FDIF)
1684 TLEAF = TAIR + FDIF
1686 !*---third-round calculation to determine potential photosynthesis
1687 !* and transpiration
1688 CALL ICO2 (TLEAF,DVP,FVPD,CO2A,C3C4, SVPL,CO2I)
1689 CALL PHOTO (C3C4,PAR,TLEAF,CO2I,NP,EAJMAX,XVN,XJN,THETA,PLF,LRD)
1691 SLOPEL = (SVPL-SVP)/NOTNUL(TLEAF-TAIR)
1693 CALL GCRSW (PLF,LRD,TLEAF,CO2A,CO2I,RBW,RT, RSW)
1694 CALL PTRAN (RSW,RT,RBW,RBH,ATRJ,ATMTR,FRAC,TLEAF,DVP, &
1695 SLOPEL,VPD, PT, PH, NRADC)
1697 CALL DIFLA (FNRADC,FPT,RBH,RT, FDIF)
1699 TLEAF = TAIR + FDIF
1701 !*---fourth-round calculation to determine potential photosynthesis
1702 !* and transpiration
1703 CALL ICO2 (TLEAF,DVP,FVPD,CO2A,C3C4, SVPL,CO2I)
1704 CALL PHOTO (C3C4,PAR,TLEAF,CO2I,NP,EAJMAX,XVN,XJN,THETA,PLF,LRD)
1706 SLOPEL = (SVPL-SVP)/NOTNUL(TLEAF-TAIR)
1708 CALL GCRSW (PLF,LRD,TLEAF,CO2A,CO2I,RBW,RT, RSW)
1709 CALL PTRAN (RSW,RT,RBW,RBH,ATRJ,ATMTR,FRAC,TLEAF,DVP, &
1710 SLOPEL,VPD, PT, PH, NRADC)
1712 RETURN
1713 END SUBROUTINE PPHTR
1716 !*----------------------------------------------------------------------*
1717 !* SUBROUTINE PEVAP *
1718 !* Purpose: This subroutine calculates potential soil evaporation. *
1719 !* *
1720 !* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) *
1721 !* name type meaning units class *
1722 !* ---- ---- ------- ----- ----- *
1723 !* TAIR R4 Air temperature oC I *
1724 !* DVP R4 Vapour pressure kPa I *
1725 !* RSS R4 Soil resistance,equivalent to leaf stomata s m-1 I *
1726 !* RTS R4 Turbulence resistance for soil s m-1 I *
1727 !* RBWS R4 Soil boundary layer resistance to water s m-1 I *
1728 !* RBHS R4 Soil boundary layer resistance to heat s m-1 I *
1729 !* ATRJS R4 Absorbed global radiation by soil J m-2 s-1 I *
1730 !* ATMTR R4 Atmospheric transmissivity - I *
1731 !* PT1 R4 Potential leaf transpiration using water mm s-1 I *
1732 !* from upper evaporative soil layer *
1733 !* WSUP1 R4 Water supply from upper evaporative soil mm s-1 I *
1734 !* layer for evapotranspiration *
1735 !* PESOIL R4 Potential soil evaporation mm s-1 O *
1736 !* PHSOIL R4 Potential soil sensible heat flux J m-2 s-1 O *
1737 !* NRADS R4 Net soil absorbed radiation J m-2 s-1 O *
1738 !*----------------------------------------------------------------------*
1739 SUBROUTINE PEVAP (TAIR,DVP,RSS,RTS,RBWS,RBHS,ATRJS,ATMTR, &
1740 PT1,WSUP1,PESOIL,PHSOIL,NRADS)
1741 IMPLICIT REAL (A-Z)
1743 !*--- first-round calculation to estimate soil surface temperature (TAVS)
1744 SVP = 0.611*EXP(17.4*TAIR/(TAIR+239.))
1745 VPD = MAX (0., SVP-DVP)
1746 SLOPE = 4158.6 * SVP/(TAIR + 239.)**2
1747 CALL PTRAN(RSS,RTS,RBWS,RBHS,ATRJS,ATMTR,1.,TAIR,DVP, &
1748 SLOPE,VPD, FPE, FPH, FNRADS)
1749 FPESOL = MAX(0., FPE)
1750 FAESOL = MIN(FPESOL,FPESOL/(PT1+FPESOL)*WSUP1)
1751 CALL DIFLA (FNRADS,FAESOL,RBHS,RTS, FDIFS)
1752 TAVS = TAIR + FDIFS
1754 !*---second-round calculation to estimate potential soil evaporation
1755 SVPS = 0.611*EXP(17.4*TAVS/(TAVS+239.))
1756 SLOPES = (SVPS-SVP)/NOTNUL(FDIFS)
1758 CALL PTRAN(RSS,RTS,RBWS,RBHS,ATRJS,ATMTR,1.,TAVS,DVP, &
1759 SLOPES,VPD, PE, PH, NRADS)
1760 PESOIL = MAX(0., PE)
1761 PHSOIL = PH
1763 RETURN
1764 END SUBROUTINE PEVAP
1766 !*----------------------------------------------------------------------*
1767 !* SUBROUTINE APHTR *
1768 !* Purpose: This subroutine calculates actual leaf photosynthesis when *
1769 !* water stress occurs. *
1770 !* *
1771 !* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) *
1772 !* name type meaning units class *
1773 !* ---- ---- ------- ----- ----- *
1774 !* TAIR R4 Air temperature oC I *
1775 !* PAR R4 Absorbed photosynth. active radiation J m-2 s-1 I *
1776 !* DVP R4 Vapour pressure kPa I *
1777 !* CO2A R4 Ambient CO2 concentration ml m-3 I *
1778 !* C3C4 R4 Crop type (=1. for C3, -1 for C4 crops) - I *
1779 !* FVPD R4 Slope for linear effect of VPD on Ci/Ca (kPa)-1 I *
1780 !* NRADC R4 Net leaf absorbed radiation J m-2 s-1 I *
1781 !* AT R4 Actual leaf transpiration mm s-1 I *
1782 !* PT R4 Potential leaf transpiration mm s-1 I *
1783 !* RT R4 Turbulence resistance s m-1 I *
1784 !* RBH R4 Leaf boundary layer resistance to heat s m-1 I *
1785 !* RBW R4 Leaf boundary layer resistance to water s m-1 I *
1786 !* RSW R4 Potential stomatal resistance to water s m-1 I *
1787 !* SLOPEL R4 Slope of saturated vapour pressure curve kPa oC-1 I *
1788 !* NP R4 Photosynthet. active leaf N content g m-2 I *
1789 !* NPN R4 NP with small plant-N increment g m-2 I *
1790 !* EAJMAX R4 Energy of activation for Jmax J mol-1 I *
1791 !* XVN R4 Slope of linearity between Vcmax & leaf N umol/g/s I *
1792 !* XJN R4 Slope of linearity between Jmax & leaf N umol/g/s I *
1793 !* THETA R4 Convexity for light response of e-transport - I *
1794 !* PLFAS R4 Actual leaf photosynthesis gCO2/m2/s O *
1795 !* PLFAN R4 PLFAS with small plant-N increment gCO2/m2/s O *
1796 !* ADIF R4 Actual leaf-air temperature difference oC O *
1797 !* ARSW R4 Actual stomatal resistance to water s m-1 O *
1798 !*----------------------------------------------------------------------*
1799 SUBROUTINE APHTR(TAIR,PAR,DVP,CO2A,C3C4,FVPD,NRADC,AT,PT,RT,RBH, &
1800 RBW,RSW,SLOPEL,NP,NPN,EAJMAX,XVN,XJN,THETA, &
1801 PLFAS,PLFAN,ADIF,ARSW)
1802 IMPLICIT REAL (A-Z)
1804 PSYCH = 0.067 !psychrometric constant (kPa/oC)
1806 !*---leaf temperature if water stress occurs
1807 CALL DIFLA (NRADC,AT,RBH,RT, ADIF)
1808 ATLEAF = TAIR + ADIF
1810 !*---stomatal resistance to water if water stress occurs
1811 ARSW = (PT-AT)*(SLOPEL*(RBH+RT)+PSYCH*(RBW+RT))/AT/PSYCH+PT/AT*RSW
1813 !*---potential photosynthesis at the new leaf temperature
1814 CALL ICO2 (ATLEAF,DVP,FVPD,CO2A,C3C4, SVPA,ACO2I)
1815 CALL PHOTO(C3C4,PAR,ATLEAF,ACO2I,NPN,EAJMAX,XVN,XJN,THETA,APLFN,ARDN)
1816 CALL PHOTO(C3C4,PAR,ATLEAF,ACO2I,NP,EAJMAX,XVN,XJN,THETA,APLF,ARD)
1818 !*---actual photosynthesis under water stress condition
1819 PLFAS = (1.6*RSW+1.3*RBW+RT)/(1.6*ARSW+1.3*RBW+RT)*(APLF-ARD)+ARD
1820 PLFAN = (1.6*RSW+1.3*RBW+RT)/(1.6*ARSW+1.3*RBW+RT)*(APLFN-ARDN)+ARDN
1822 RETURN
1823 END SUBROUTINE APHTR
1826 !*----------------------------------------------------------------------*
1827 !* SUBROUTINE PTRAN *
1828 !* Purpose: This subroutine calculates leaf transpiration, using the *
1829 !* Penman-Monteith equation *
1830 !* *
1831 !* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) *
1832 !* name type meaning units class *
1833 !* ---- ---- ------- ----- ----- *
1834 !* RSW R4 Potential stomatal resistance to water s m-1 I *
1835 !* RT R4 Turbulence resistance s m-1 I *
1836 !* RBW R4 Leaf boundary layer resistance to water s m-1 I *
1837 !* RBH R4 Leaf boundary layer resistance to heat s m-1 I *
1838 !* ATRJ R4 Absorbed global radiation J m-2 s-1 I *
1839 !* ATMTR R4 Atmospheric transmissivity - I *
1840 !* FRAC R4 Fraction of leaf classes (sunlit vs shaded)- I *
1841 !* TLEAF R4 Leaf temperature oC I *
1842 !* DVP R4 Vapour pressure kPa I *
1843 !* SLOPE R4 Slope of saturated vapour pressure curve kPa oC-1 I *
1844 !* VPD R4 Saturation vapour pressure deficit of air kPa I *
1845 !* PT R4 Potential leaf transpiration mm s-1 O *
1846 !* PH R4 Potential leaf sensible heat flux J m-2 s-1 O *
1847 !* NRADC R4 Net leaf absorbed radiation J m-2 s-1 O *
1848 !*----------------------------------------------------------------------*
1849 SUBROUTINE PTRAN (RSW,RT,RBW,RBH,ATRJ, &
1850 ATMTR,FRAC,TLEAF,DVP,SLOPE,VPD,PT,PH,NRADC)
1851 IMPLICIT REAL (A-Z)
1853 !*---some physical constants
1854 BOLTZM = 5.668E-8 !Stefan-Boltzmann constant(J/m2/s/K4)
1855 LHVAP = 2.4E6 !latent heat of water vaporization(J/kg)
1856 VHCA = 1200. !volumetric heat capacity (J/m3/oC)
1857 PSYCH = 0.067 !psychrometric constant (kPa/oC)
1859 !*---net absorbed radiation
1860 CLEAR = MAX(0., MIN(1., (ATMTR-0.25)/0.45)) !sky clearness
1861 BBRAD = BOLTZM*(TLEAF +273.)**4
1862 RLWN = BBRAD*(0.56-0.079*SQRT(DVP*10.))*(0.1+0.9*CLEAR)*FRAC
1863 NRADC = ATRJ - RLWN
1865 !*---intermediate variable related to resistances
1866 PSR = PSYCH*(RBW+RT+RSW)/(RBH+RT)
1868 !*---Compute PT
1869 !*---radiation-determined term
1870 PTR = NRADC*SLOPE/(SLOPE+PSR)/LHVAP
1872 !*---vapour pressure-determined term
1873 PTD = (VHCA*VPD/(RBH+RT))/(SLOPE+PSR)/LHVAP
1875 !*---potential evaporation or transpiration
1876 PT = MAX(1.E-10,PTR+PTD)
1878 !*---ji
1879 !*---Compute PH in W/m2 (see Bolan p. 202)
1880 !*---radiation-determined term
1881 PHR = NRADC*PSR/(SLOPE+PSR)
1883 !*---vapour pressure-determined term
1884 PHD = (VHCA*VPD/(RBH+RT))/(SLOPE+PSR)
1886 !*---potential evaporation or transpiration
1887 PH = PHR-PHD
1889 RETURN
1890 END SUBROUTINE PTRAN
1893 !*----------------------------------------------------------------------*
1894 !* SUBROUTINE DIFLA *
1895 !* Purpose: This subroutine calculates leaf(canopy)-air temperature *
1896 !* differential. *
1897 !* *
1898 !* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) *
1899 !* name type meaning units class *
1900 !* ---- ---- ------- ----- ----- *
1901 !* NRADC R4 Net leaf absorbed radiation J m-2 s-1 I *
1902 !* PT R4 Potential leaf transpiration mm s-1 I *
1903 !* RBH R4 Leaf boundary layer resistance to heat s m-1 I *
1904 !* RT R4 Turbulence resistance s m-1 I *
1905 !* DIF R4 Leaf-air temperature difference oC O *
1906 !*----------------------------------------------------------------------*
1907 SUBROUTINE DIFLA (NRADC,PT,RBH,RT, DIF)
1908 IMPLICIT REAL (A-Z)
1910 LHVAP = 2.4E6 !latent heat of water vaporization(J/kg)
1911 VHCA = 1200. !volumetric heat capacity (J/m3/oC)
1913 DIF = LIMIT (-25., 25., (NRADC-LHVAP*PT)*(RBH+RT)/VHCA)
1915 RETURN
1916 END SUBROUTINE DIFLA
1917 !*----------------------------------------------------------------------*
1918 !* SUBROUTINE ICO2 *
1919 !* Purpose: This subroutine calculates the internal CO2 concentration *
1920 !* as affected by vapour pressure deficit. *
1921 !* *
1922 !* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) *
1923 !* name type meaning units class *
1924 !* ---- ---- ------- ----- ----- *
1925 !* TLEAF R4 Leaf temperature oC I *
1926 !* DVP R4 Vapour pressure kPa I *
1927 !* FVPD R4 Slope for linear effect of VPDL on Ci/Ca (kPa)-1 I *
1928 !* CO2A R4 Ambient CO2 concentration ml m-3 I *
1929 !* C3C4 R4 Crop type (=1. for C3, -1 for C4 crops) - I *
1930 !* SVPL R4 Saturated vapour pressure of leaf kPa O *
1931 !* CO2I R4 intercellular CO2 concentration ml m-3 O *
1932 !*----------------------------------------------------------------------*
1933 SUBROUTINE ICO2 (TLEAF,DVP,FVPD,CO2A,C3C4, SVPL,CO2I)
1934 IMPLICIT REAL (A-Z)
1936 !*---air-to-leaf vapour pressure deficit
1937 SVPL = 0.611 * EXP(17.4 * TLEAF / (TLEAF + 239.))
1938 VPDL = MAX (0., SVPL - DVP)
1940 !*---Michaelis-Menten const. for CO2 at 25oC (umol/mol)
1941 KMC25 = INSW(C3C4, 650., 404.9) !greater KMC25 for C4 than C3
1943 !*---Michaelis-Menten const. for O2 at 25oC (mmol/mol)
1944 KMO25 = INSW(C3C4, 450., 278.4) !greater KMO25 for C4 than C3
1946 !*---CO2 compensation point in absence of dark respiration (GAMMAX)
1947 O2 = 210. !oxygen concentration(mmol/mol)
1948 EAVCMX = 65330. !energy of activation for Vcmx(J/mol)
1949 EAKMC = 79430. !energy of activation for KMC (J/mol)
1950 EAKMO = 36380. !energy of activation for KMO (J/mol)
1951 EARD = 46390. !energy of activation for dark respiration(J/mol)
1952 RDVX25 = 0.0089 !ratio of dark respiration to Vcmax at 25oC
1953 TO = 298.15
1955 KMC = KMC25*EXP((1./TO-1./(TLEAF+273.))*EAKMC/8.314)
1956 KMO = KMO25*EXP((1./TO-1./(TLEAF+273.))*EAKMO/8.314)
1957 GAMMAX = 0.5*EXP(-3.3801+5220./(TLEAF+273.)/8.314)*O2*KMC/KMO
1959 !*---CO2 compensation point (GAMMA)
1960 RDVCX = RDVX25*EXP((1./TO-1./(TLEAF+273.))*(EARD-EAVCMX)/8.314)
1961 GAMMA0 = (GAMMAX+RDVCX*KMC*(1.+O2/KMO))/(1.-RDVCX)
1962 GAMMA_ = INSW (C3C4, GAMMA0/10., GAMMA0)
1964 !*---internal/ambient CO2 ratio, based on data of Morison & Gifford (1983)
1965 RCICA = 1.-(1.-GAMMA_/CO2A)*(0.14+FVPD*VPDL)
1967 !*---intercellular CO2 concentration
1968 CO2I = RCICA * CO2A
1970 RETURN
1971 END SUBROUTINE ICO2
1974 !*----------------------------------------------------------------------*
1975 !* SUBROUTINE GCRSW *
1976 !* Purpose: This subroutine calculates overall leaf conductance *
1977 !* for CO2 (GC) and the stomatal resistance to water (RSW). *
1978 !* *
1979 !* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) *
1980 !* name type meaning units class *
1981 !* ---- ---- ------- ----- ----- *
1982 !* PLEAF R4 Gross leaf photosynthesis gCO2/m2/s I *
1983 !* RDLEAF R4 Leaf dark respiration gCO2/m2/s I *
1984 !* TLEAF R4 Leaf temperature oC I *
1985 !* CO2A R4 Ambient CO2 concentration ml m-3 I *
1986 !* CO2I R4 Internal CO2 concentration ml m-3 I *
1987 !* RT R4 Turbulence resistance s m-1 I *
1988 !* RBW R4 Leaf boundary layer resistance to water s m-1 I *
1989 !* RSW R4 Potential stomatal resistance to water s m-1 O *
1990 !*----------------------------------------------------------------------*
1991 SUBROUTINE GCRSW (PLEAF,RDLEAF,TLEAF,CO2A,CO2I,RBW,RT, RSW)
1992 IMPLICIT REAL (A-Z)
1994 !*---potential conductance for CO2
1995 !*ji 13.7.11 MAX routine faengt negative Leitfaehigkeiten ab
1996 GC = MAX(1E-6,(PLEAF-RDLEAF)*(273.+TLEAF)/0.53717/(CO2A-CO2I))
1998 !*---potential stomatal resistance to water
1999 RSW = MAX(1E-30, 1./GC - RBW*1.3 - RT)/1.6
2001 RETURN
2002 END SUBROUTINE GCRSW
2005 !*----------------------------------------------------------------------*
2006 !* SUBROUTINE PAN *
2007 !* Purpose: This subroutine calculates photosynthetically active *
2008 !* nitrogen content for sunlit and shaded parts of canopy. *
2009 !* *
2010 !* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) *
2011 !* name type meaning units class *
2012 !* ---- ---- ------- ----- ----- *
2013 !* SLNT R4 Top-leaf nitrogen content g m-2 I *
2014 !* SLNMIN R4 Minimum or base SLNT for photosynthesis g m-2 I *
2015 !* LAI R4 (green)Leaf area index m2 m-2 I *
2016 !* KN R4 Leaf nitrogen extinction coefficient m2 m-2 I *
2017 !* KB R4 Direct beam radiation extinction coeff. m2 m-2 I *
2018 !* NPSU R4 Photosynthet. active N for sunlit leaves g m-2 O *
2019 !* NPSH R4 Photosynthet. active N for shaded leaves g m-2 O *
2020 !*----------------------------------------------------------------------*
2021 SUBROUTINE PAN(SLNT,SLNMIN,LAI,KN,KB, NPSU,NPSH)
2022 IMPLICIT REAL (A-Z)
2024 !*---total photosynthetic nitrogen in canopy
2025 NPC = SLNT*(1.-EXP(-KN*LAI))/KN-SLNMIN*LAI
2027 !*---photosynthetic nitrogen for sunlit and shaded parts of canopy
2028 NPSU = SLNT*(1.-EXP(-(KN+KB)*LAI))/(KN+KB)-SLNMIN*(1.-EXP(-KB*LAI))/KB
2029 NPSH = NPC-NPSU
2031 RETURN
2032 END SUBROUTINE PAN
2035 !*----------------------------------------------------------------------*
2036 !* SUBROUTINE PHOTO *
2037 !* Purpose: This subroutine calculates leaf photosynthesis and dark *
2038 !* respiration, based on a renewed Farquhar biochemistry *
2039 !* (cf Yin et al.2004. Plant, Cell & Environment 27:1211-1222)*
2040 !* *
2041 !* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) *
2042 !* name type meaning units class *
2043 !* ---- ---- ------- ----- ----- *
2044 !* C3C4 R4 Crop type (=1. for C3, -1. for C4 crops) - I *
2045 !* PAR R4 Leaf absorbed photosynth. active radiance J m-2 s-1 I *
2046 !* TLEAF R4 Leaf temperature oC I *
2047 !* CO2I R4 Intercellular CO2 concentration ml m-3 I *
2048 !* NP R4 Photosynthetically active leaf N content g m-2 I *
2049 !* EAJMAX R4 Energy of activation for Jmax J mol-1 I *
2050 !* XVN R4 Slope of linearity between Vcmax & leaf N umol/g/s I *
2051 !* XJN R4 Slope of linearity between Jmax & leaf N umol/g/s I *
2052 !* THETA R4 Convexity for light response of e-transport - I *
2053 !* PLEAF R4 Gross leaf photosynthesis gCO2/m2/s O *
2054 !* RDLEAF R4 Leaf dark respiration gCO2/m2/s O *
2055 !*----------------------------------------------------------------------*
2056 SUBROUTINE PHOTO(C3C4,PAR,TLEAF,CO2I,NP,EAJMAX,XVN,XJN, &
2057 THETA,PLEAF,RDLEAF)
2058 IMPLICIT REAL (A-Z)
2060 !*---Michaelis-Menten constants for CO2 and O2 at 25oC
2061 IF (C3C4.LT.0.) THEN
2062 KMC25 = 650. !greater KMC25 for C4 than C3; unit:(umol/mol)
2063 KMO25 = 450. !greater KMO25 for C4 than C3; unit:(mmol/mol)
2064 ELSE
2065 KMC25 = 404.9 !unit:(umol/mol)
2066 KMO25 = 278.4 !unit:(mmol/mol)
2067 ENDIF
2069 !*---other constants related to the Farquhar-type photosynthesis model
2070 O2 = 210. !oxygen concentration(mmol/mol)
2071 EAVCMX = 65330. !energy of activation for Vcmx(J/mol)
2072 EAKMC = 79430. !energy of activation for KMC (J/mol)
2073 EAKMO = 36380. !energy of activation for KMO (J/mol)
2074 EARD = 46390. !energy of activation for dark respiration(J/mol)
2075 DEJMAX = 200000. !energy of deactivation for JMAX (J/mol)
2076 SJ = 650. !entropy term in JT equation (J/mol/K)
2077 PHI2M = 0.85 !maximum electron transport efficiency of PS II
2078 HH = 3. !number of protons required to synthesise 1 ATP
2079 KTMP = 1.0 !Factor for reducing photosynthesis in case of T<5C
2080 JTMAX = 3.12
2081 TO = 298.15
2083 !*---PAR photon flux in umol/m2/s absorbed by leaf photo-systems
2084 UPAR = 4.56*PAR !4.56 conversion factor in umol/J
2086 !*---Michaelis-Menten constants for CO2 and O2 respectively
2087 KMC = KMC25*EXP((1./TO-1./(TLEAF+273.))*EAKMC/8.314)
2088 KMO = KMO25*EXP((1./TO-1./(TLEAF+273.))*EAKMO/8.314)
2090 !*---CO2 compensation point in the absence of dark respiration
2091 GAMMAX = 0.5*EXP(-3.3801+5220./(TLEAF+273.)/8.314)*O2*KMC/KMO
2093 !*---Arrhenius function for the effect of temperature on carboxylation
2094 VCT = EXP((1./TO-1./(TLEAF+273.))*EAVCMX/8.314)
2096 !*---function for the effect of temperature on electron transport
2097 JT = EXP((1./TO-1./(TLEAF+273.))*EAJMAX/8.314)* &
2098 (1.+EXP(SJ/8.314-DEJMAX/TO/8.314))/ &
2099 (1.+EXP(SJ/8.314-1./(TLEAF+273.) *DEJMAX/8.314))
2101 VCMX = XVN*VCT*NP
2102 JMAX = XJN*JT *NP
2104 !*---CO2 concentration at carboxylation site & electron pathways and
2105 !* their stoichiometries
2106 FPSEUD = 0. !assuming no pseudocyclic e- transport
2107 IF (C3C4.LT.0.) THEN
2108 ZZ = 0.2 !CO2 leakage from bundle-sheath to mesophyll
2109 CC = 10.*CO2I !to mimic C4 CO2 concentrating mechanism
2110 SF = 2.*(CC-GAMMAX)/(1.-ZZ)
2111 FQ = 1.- FPSEUD- 2.*(4.*CC+8.*GAMMAX)/HH/(SF+3.*CC+7.*GAMMAX)
2112 FCYC = FQ
2113 ELSE
2114 CC = CO2I
2115 SF = 0.
2116 FQ = 0.
2117 FCYC = 1.-(FPSEUD*HH*(SF+3.*CC+7.*GAMMAX)/(4.*CC+8.*GAMMAX)+1.)/ &
2118 (HH*(SF+3.*CC+7.*GAMMAX)/(4.*CC+8.*GAMMAX)-1.)
2119 ENDIF
2121 !*--- electron transport rate in dependence on PAR photon flux
2122 ALPHA2 = (1.-FCYC)/(1.+(1.-FCYC)/PHI2M)
2123 X = ALPHA2*UPAR/MAX(1.E-10,JMAX)
2124 J2 = JMAX*(1.+X-((1.+X)**2.-4.*X*THETA)**0.5)/2./THETA
2126 !*---rates of carboxylation limited by Rubisco and electron transport
2127 VC = VCMX * CC/(CC + KMC*(O2/KMO+1.))
2128 VJ = J2 * CC*(2.+FQ-FCYC)/HH/(SF+3.*CC+7.*GAMMAX)/(1.-FCYC)
2130 !*---gross rate of leaf photosynthesis
2131 ALF = (1.-GAMMAX/CC)*MIN(VC,VJ)
2132 PLEAF = MAX(1.E-10, (1.E-6)*44.*ALF)
2134 !*---rate of leaf dark respiration
2135 RDVX25 = 0.0089 !ratio of dark respiration to Vcmax at 25oC
2136 RDT = EXP((1./TO-1./(TLEAF+273.))*EARD/8.314)
2137 RDLEAF = (1.E-6)*44. *RDVX25*(XVN*NP) * RDT
2139 RETURN
2140 END SUBROUTINE PHOTO
2143 !*----------------------------------------------------------------------*
2144 !* SUBROUTINE REFL *
2145 !* Purpose: This subroutine calculates reflection coefficients. *
2146 !* *
2147 !* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) *
2148 !* name type meaning units class *
2149 !* ---- ---- ------- ----- ----- *
2150 !* SCP R4 Leaf scattering coefficient - I *
2151 !* KB R4 Direct beam radiation extinction coeff. m2 m-2 I *
2152 !* KBP R4 Scattered beam radiation extinction coeff. m2 m-2 O *
2153 !* PCB R4 Canopy beam radiation reflection coeff. - O *
2154 !*----------------------------------------------------------------------*
2155 SUBROUTINE REFL (SCP,KB, KBP,PCB)
2156 IMPLICIT REAL (A-Z)
2158 !*--- scattered beam radiation extinction coefficient
2159 KBP = KB*SQRT(1.-SCP)
2161 !*---canopy reflection coefficient for horizontal leaves
2162 PH = (1.-SQRT(1.-SCP))/(1.+SQRT(1.-SCP))
2164 !*---Canopy beam radiation reflection coefficient
2165 PCB = 1.-EXP(-2.*PH*KB/(1.+KB))
2167 RETURN
2168 END SUBROUTINE REFL
2171 !*----------------------------------------------------------------------*
2172 !* SUBROUTINE LIGAB *
2173 !* Purpose: This subroutine calculates absorbed light for sunlit and *
2174 !* shaded leaves. *
2175 !* *
2176 !* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) *
2177 !* name type meaning units class *
2178 !* ---- ---- ------- ----- ----- *
2179 !* SCP R4 Leaf scattering coefficient - I *
2180 !* KB R4 Direct beam radiation extinction coeff. m2 m-2 I *
2181 !* KBP R4 Scattered beam radiation extinction coeff. m2 m-2 I *
2182 !* KDP R4 Diffuse radiation extinction coefficient m2 m-2 I *
2183 !* PCB R4 Canopy beam radiation reflection coeff. - I *
2184 !* PCD R4 Canopy diffuse radiation reflection coeff. - I *
2185 !* IB0 R4 Incident direct-beam radiation J m-2 s-1 I *
2186 !* ID0 R4 Incident diffuse radiation J m-2 s-1 I *
2187 !* LAI R4 (green)Leaf area index m2 m-2 I *
2188 !* ISU R4 Absorbed radiation by sunlit leaves J m-2 s-1 O *
2189 !* ISH R4 Absorbed radiation by shaded leaves J m-2 s-1 O *
2190 !*----------------------------------------------------------------------*
2191 SUBROUTINE LIGAB (SCP,KB,KBP,KDP,PCB,PCD,IB0,ID0,LAI, ISU,ISH)
2192 IMPLICIT REAL (A-Z)
2194 !*---total absorbed light by canopy
2195 IC = (1.-PCB)*MAX(1e-30,IB0)*(1.-EXP(-KBP*LAI))+ &
2196 (1.-PCD)*MAX(1e-30,ID0)*(1.-EXP(-KDP*LAI))
2198 !*---absorbed light by sunlit and shaded fractions of canopy
2199 ISU = (1.-SCP)*MAX(1e-30,IB0)*(1.-EXP(-KB *LAI))+(1.-PCD)*MAX(1e-30,ID0)/(KDP+KB)* &
2200 KDP*(1.-EXP(-(KDP+KB)*LAI))+MAX(1e-30,IB0)*((1.-PCB)/(KBP+KB)*KBP* &
2201 (1.-EXP(-(KBP+KB)*LAI))-(1.-SCP)*(1.-EXP(-2.*KB*LAI))/2.)
2203 ISH = IC-ISU
2205 RETURN
2206 END SUBROUTINE LIGAB
2209 !*----------------------------------------------------------------------*
2210 !* SUBROUTINE KBEAM *
2211 !* Purpose: This subroutine calculates extinction coefficient for *
2212 !* direct beam radiation. *
2213 !* *
2214 !* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) *
2215 !* name type meaning units class *
2216 !* ---- ---- ------- ----- ----- *
2217 !* SINB R4 Sine of solar elevation - I *
2218 !* BL R4 Leaf angle (from horizontal) radians I *
2219 !* KB R4 Direct beam radiation extinction coeff. m2 m-2 O *
2220 !*----------------------------------------------------------------------*
2221 SUBROUTINE KBEAM (SINB,BL, KB)
2222 IMPLICIT REAL (A-Z)
2224 !*---solar elevation in radians
2225 B = ASIN(SINB)
2227 !*---average projection of leaves in the direction of a solar beam
2228 IF (SINB.GE.SIN(BL)) THEN
2229 OAV = SINB*COS(BL)
2230 ELSE
2231 OAV = 2./3.141592654*(SINB*COS(BL)*ASIN(TAN(B)/TAN(BL)) &
2232 +((SIN(BL))**2-SINB**2)**0.5)
2233 ENDIF
2235 !*---beam radiation extinction coefficient
2236 KB = OAV/SINB
2238 RETURN
2239 END SUBROUTINE KBEAM
2242 !*----------------------------------------------------------------------*
2243 !* SUBROUTINE KDIFF *
2244 !* Purpose: This subroutine calculates extinction coefficient for *
2245 !* diffuse radiation. *
2246 !* *
2247 !* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) *
2248 !* name type meaning units class *
2249 !* ---- ---- ------- ----- ----- *
2250 !* LAI R4 Total leaf area index m2 m-2 I *
2251 !* BL R4 Leaf angle (from horizontal) radians I *
2252 !* SCP R4 Leaf scattering coefficient - I *
2253 !* KDP R4 Diffuse radiation extinction coefficient m2 m-2 O *
2254 !*----------------------------------------------------------------------*
2255 SUBROUTINE KDIFF (LAI,BL,SCP, KDP)
2256 IMPLICIT REAL (A-Z)
2258 PI = 3.141592654
2260 !*---extinction coefficient of beam lights from 15, 45 and 75 elevations
2261 CALL KBEAM (SIN(15.*PI/180.),BL, KB15)
2262 CALL KBEAM (SIN(45.*PI/180.),BL, KB45)
2263 CALL KBEAM (SIN(75.*PI/180.),BL, KB75)
2265 !*---diffuse light extinction coefficient
2266 KDP = -1./LAI*LOG(0.178*EXP(-KB15*(1.-SCP)**0.5*LAI) &
2267 +0.514*EXP(-KB45*(1.-SCP)**0.5*LAI) &
2268 +0.308*EXP(-KB75*(1.-SCP)**0.5*LAI))
2270 RETURN
2271 END SUBROUTINE KDIFF
2273 !*----------------------------------------------------------------------*
2274 !* FUNCTION INTGRL *
2275 !* Purpose: This function integrates a differential equation *
2276 !* using the Euler method. Substitutes the intrinsic FST function *
2277 !* FORMAL PARAMETERS: (I=input,O=output,C=control,IN=init,T=time) *
2278 !* name type meaning units class *
2279 !* ---- ---- ------- ----- ----- *
2280 !* POOL R4 Currenmt stock size g X m-2 I *
2281 !* RATE R4 Rate of change 1/d I *
2282 !* DT R4 Time step d I *
2283 !* INTGRL R4 Returns the stock size in next time step g X m-2 0 *
2284 FUNCTION INTGRL(POOL, RATE, DT)
2285 IMPLICIT NONE
2286 real:: INTGRL, POOL, RATE, DT
2288 INTGRL = POOL + RATE*DT
2290 RETURN
2291 END FUNCTION INTGRL
2293 !*----------------------------------------------------------------------*
2294 !* FUNCTION NOTNUL *
2295 !* Y = NOTNUL(X) *
2296 !* Y is equal to X but 1.0 in case of x=0.0. Note that X is *
2297 !* evaluated without any tolerance interval *
2298 !* name type meaning units class *
2299 !* ---- ---- ------- ----- ----- *
2300 !* X R4 Variable to be checked I *
2301 !*----------------------------------------------------------------------*
2302 FUNCTION NOTNUL(X)
2303 IMPLICIT NONE
2304 real:: NOTNUL, X
2306 if (X.ne.0) then
2307 NOTNUL = X
2308 else
2309 NOTNUL = 1
2310 endif
2312 RETURN
2313 END FUNCTION NOTNUL
2315 !*----------------------------------------------------------------------*
2316 !* FUNCTION LIMIT *
2317 !* Y = LIMIT(XL, XH, X) *
2318 !* Y is equal to X but limited between XL and XH *
2319 !* Y - Returned as X bounded on [Xl,XH] *
2320 !* XL- Lower bound of X *
2321 !* XH- Upper bound of X *
2322 !* ---- ---- ------- *
2323 !* X R4 Variable to be checken I *
2324 !*----------------------------------------------------------------------*
2325 FUNCTION LIMIT(XL, XH, X)
2326 IMPLICIT NONE
2327 real:: LIMIT, XL, XH, X
2329 if (X.ge.XL.and.X.le.XH) then
2330 LIMIT = X
2331 else
2332 if (X.gt.XH) then
2333 LIMIT = XH
2334 else
2335 LIMIT = XL
2336 endif
2337 endif
2339 RETURN
2340 END FUNCTION LIMIT
2342 !*----------------------------------------------------------------------*
2343 !* FUNCTION INSW *
2344 !* Y = INSW(X, Y1, Y2) *
2345 !* Input switch. Y is set equal to Y1 orY2 depending on the value of X *
2346 !* Y - Returned as either Y1 or Y2 *
2347 !* X - Control variable *
2348 !* Y1- Returned value of Y if X<0 *
2349 !* Y2- Returned value of Y if X>=0 *
2350 !* ---- ---- ------- *
2351 !* X R4 Variable to be checken I *
2352 !*----------------------------------------------------------------------*
2353 FUNCTION INSW(X, Y1, Y2)
2354 implicit none
2355 real :: INSW, X, Y1, Y2
2357 if (X.lt.0) then
2358 INSW = Y1
2359 else
2360 INSW = Y2
2361 endif
2363 RETURN
2364 END FUNCTION INSW
2367 !*----------------------------------------------------------------------*
2368 !* FUNCTION REAAND *
2369 !* Y = REAAND(X1, X2) *
2370 !* Returns 1.0 if both input variables are positive, otherwise Y=0.0 *
2371 !* ---- ---- ------- *
2372 !* X1 R4 1. variable to be checked I *
2373 !* X2 R4 2. variable to be checked I *
2374 !*----------------------------------------------------------------------*
2375 FUNCTION REAAND(X1, X2)
2376 implicit none
2377 real :: REAAND, X1, X2
2379 if (X1.gt.0.and.X2.gt.0) then
2380 REAAND = 1.
2381 else
2382 REAAND = 0.
2383 endif
2385 RETURN
2386 END FUNCTION REAAND
2389 !*----------------------------------------------------------------------*
2390 !* FUNCTION REAnOR *
2391 !* Y = REANOR(X1, X2) *
2392 !* Returns 1.0 if both input variables are less than or equal to 0., *
2393 !* otherwise Y=0.
2394 !* ---- ---- ------- *
2395 !* X1 R4 1. variable to be checked I *
2396 !* X2 R4 2. variable to be checked I *
2397 !*----------------------------------------------------------------------*
2398 FUNCTION REANOR(X1, X2)
2399 implicit none
2400 real :: REANOR, X1, X2
2402 if (X1.le.0.and.X2.le.0) then
2403 REANOR = 1.
2404 else
2405 REANOR = 0.
2406 endif
2408 RETURN
2409 END FUNCTION REANOR
2411 END MODULE module_sf_gecros