Context Navigation

source: branches/publications/ORCHIDEE_CAN_r2290/src_stomate/stomate_turnover.f90 @ 5242

Last change on this file since 5242 was 2275, checked in by sebastiaan.luyssaert, 10 years ago
DEV: works on a single pixel. Experimenting with adaptation to waterstress. Committed to migrate the code to Curie. Do not use for regional simulations. Updates will follow.
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 122.8 KB

Line
1	! =================================================================================================================================
2	! MODULE : stomate_turnover.f90
3	!
4	! CONTACT : orchidee-help _at_ ipsl.jussieu.fr
5	!
6	! LICENCE : IPSL (2006)
7	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
8	!
9	!>\BRIEF This module manages the end of the growing season and calculates herbivory
10	! and turnover of leaves, fruits, fine roots.
11	!!
12	!!\n DESCRIPTION: This subroutine calculates leaf senescence due to climatic conditions or
13	!! as a function of leaf age and new LAI, and subsequent turnover of the different plant
14	!! biomass compartments (sections 1 to 6), herbivory (section 7), fruit turnover for trees
15	!! (section 8) and sapwood conversion (section 9).
16	!!
17	!! RECENT CHANGE(S): None
18	!!
19	!! SVN :
20	!! $HeadURL$
21	!! $Date$
22	!! $Revision$
23	!! \n
24	!_ ================================================================================================================================
25
26	MODULE stomate_turnover
27
28	! modules used:gg
29	USE xios_orchidee
30	USE ioipsl_para
31	USE stomate_data
32	USE constantes
33	USE pft_parameters
34	USE sapiens_agriculture, ONLY: crop_harvest
35	USE function_library, ONLY : biomass_to_lai
36
37	IMPLICIT NONE
38
39	! private & public routines
40
41	PRIVATE
42	PUBLIC turnover_diagnostic, turnover_prognostic, turnover_clear
43
44	LOGICAL, SAVE :: firstcall = .TRUE. !! first call (true/false)
45	!$OMP THREADPRIVATE(firstcall)
46
47	CONTAINS
48
49
50	!! ================================================================================================================================
51	!! SUBROUTINE : turnover_clear
52	!!
53	!>\BRIEF Set flag ::firstcall to .TRUE., and therefore activate section 1
54	!! of subroutine turn which writes a message to the output.
55	!!
56	!_ ================================================================================================================================
57
58	SUBROUTINE turnover_clear
59	firstcall=.TRUE.
60	END SUBROUTINE turnover_clear
61
62
63	!! ================================================================================================================================
64	!! SUBROUTINE : turn
65	!!
66	!>\BRIEF Calculate turnover of leaves, roots, fruits and sapwood due to aging or climatic
67	!! induced senescence. Calculate herbivory.
68	!!
69	!! DESCRIPTION : This subroutine determines the turnover of leaves and fine roots (and stems for grasses)
70	!! and simulates following processes:
71	!! 1. Mean leaf age is calculated from leaf ages of separate leaf age classes. Should actually
72	!! be recalculated at the end of the routine, but it does not change too fast. The mean leaf
73	!! age is calculated using the following equation:
74	!! \latexonly
75	!! \input{turnover_lma_update_eqn1.tex}
76	!! \endlatexonly
77	!! \n
78	!! 2. Meteorological senescence: the detection of the end of the growing season and shedding
79	!! of leaves, fruits and fine roots due to unfavourable meteorological conditions.
80	!! The model distinguishes three different types of "climatic" leaf senescence, that do not
81	!! change the age structure: sensitivity to cold temperatures, to lack of water, or both.
82	!! If meteorological conditions are fulfilled, a flag ::senescence is set to TRUE. Note
83	!! that evergreen species do not experience climatic senescence.
84	!! Climatic senescence is triggered by sensitivity to cold temperatures where the critical
85	!! temperature for senescence is calculated using the following equation:
86	!! \latexonly
87	!! \input{turnover_temp_crit_eqn2.tex}
88	!! \endlatexonly
89	!! \n
90	!! Climatic senescence is triggered by sensitivity to lack of water availability where the
91	!! moisture availability critical level is calculated using the following equation:
92	!! \latexonly
93	!! \input{turnover_moist_crit_eqn3.tex}
94	!! \endlatexonly
95	!! \n
96	!! Climatic senescence is triggered by sensitivity to temperature or to lack of water where
97	!! critical temperature and moisture availability are calculated as above.\n
98	!! Trees in climatic senescence lose their fine roots at the same rate as they lose their leaves.
99	!! The rate of biomass loss of both fine roots and leaves is presribed through the equation:
100	!! \latexonly
101	!! \input{turnover_clim_senes_biomass_eqn4.tex}
102	!! \endlatexonly
103	!! \n
104	!! with ::leaffall(j) a PFT-dependent time constant which is given in
105	!! ::stomate_constants. In grasses, leaf senescence is extended to the whole plant
106	!! (all carbon pools) except to its carbohydrate reserve.
107	!! 3. Senescence due to aging: the loss of leaves, fruits and biomass due to aging
108	!! At a certain age, leaves fall off, even if the climate would allow a green plant
109	!! all year round. Even if the meteorological conditions are favorable for leaf maintenance,
110	!! plants, and in particular, evergreen trees, have to renew their leaves simply because the
111	!! old leaves become inefficient. Roots, fruits (and stems for grasses) follow leaves.
112	!! The ??senescence?? rate varies with leaf age. Note that plant is not declared senescent
113	!! in this case (wchich is important for allocation: if the plant loses leaves because of
114	!! their age, it can renew them). The leaf turnover rate due to aging of leaves is calculated
115	!! using the following equation:
116	!! \latexonly
117	!! \input{turnover_age_senes_biomass_eqn5.tex}
118	!! \endlatexonly
119	!! \n
120	!! Drop all leaves if there is a very low leaf mass during senescence. After this, the biomass
121	!! of different carbon pools both for trees and grasses is set to zero and the mean leaf age
122	!! is reset to zero. Finally, the leaf fraction and leaf age of the different leaf age classes
123	!! is set to zero. For deciduous trees: next to leaves, also fruits and fine roots are dropped.
124	!! For grasses: all aboveground carbon pools, except the carbohydrate reserves are affected:
125	!! 4. Update the leaf biomass, leaf age class fraction and the LAI
126	!! Older leaves will fall more frequently than younger leaves and therefore the leaf age
127	!! distribution needs to be recalculated after turnover. The fraction of biomass in each
128	!! leaf class is updated using the following equation:
129	!! \latexonly
130	!! \input{turnover_update_LeafAgeDistribution_eqn6.tex}
131	!! \endlatexonly
132	!! \n
133	!! 5. Simulate herbivory activity and update leaf and fruits biomass. Herbivore activity
134	!! affects the biomass of leaves and fruits as well as stalks (only for grasses).
135	!! However, herbivores do not modify leaf age structure.
136	!! 6. Calculates fruit turnover for trees. Trees simply lose their fruits with a time
137	!! constant ::tau_fruit(j), that is set to 90 days for all PFTs in ::stomate_constants
138	!! 7. Convert sapwood to heartwood for trees and update heart and softwood above and
139	!! belowground biomass. Sapwood biomass is converted into heartwood biomass
140	!! with a time constant tau ::tau_sap(j) of 1 year. Note that this biomass conversion
141	!! is not added to "turnover" as the biomass is not lost. For the updated heartwood,
142	!! the sum of new heartwood above and new heartwood below after converting sapwood to
143	!! heartwood, is saved as ::hw_new(:). Creation of new heartwood decreases the age of
144	!! the plant ??carbon?? with a factor that is determined by: old heartwood ::hw_old(:)
145	!! divided by the new heartwood ::hw_new(:)
146	!!
147	!! RECENT CHANGE(S) : None
148	!!
149	!! MAIN OUTPUT VARIABLES: ::Biomass of leaves, fruits, fine roots and sapwood above (latter for grasses only),
150	!! ::Update LAI, ::Update leaf age distribution with new leaf age class fraction
151	!!
152	!! REFERENCE(S) :
153	!! - Krinner, G., N. Viovy, N. de Noblet-Ducoudre, J. Ogee, J. Polcher, P.
154	!! Friedlingstein, P. Ciais, S. Sitch and I.C. Prentice (2005), A dynamic global
155	!! vegetation model for studies of the coupled atmosphere-biosphere system, Global
156	!! Biogeochemical Cycles, 19, doi:10.1029/2003GB002199.
157	!! - McNaughton, S. J., M. Oesterheld, D. A. Frank and K. J. Williams (1989),
158	!! Ecosystem-level patterns of primary productivity and herbivory in terrestrial habitats,
159	!! Nature, 341, 142-144, 1989.
160	!! - Sitch, S., C. Huntingford, N. Gedney, P. E. Levy, M. Lomas, S. L. Piao, , Betts, R., Ciais, P., Cox, P.,
161	!! Friedlingstein, P., Jones, C. D., Prentice, I. C. and F. I. Woodward : Evaluation of the terrestrial carbon
162	!! cycle, future plant geography and climate-carbon cycle feedbacks using 5 dynamic global vegetation
163	!! models (dgvms), Global Change Biology, 14(9), 2015â2039, 2008.
164	!!
165	!! FLOWCHART :
166	!! \latexonly
167	!! \includegraphics[scale=0.5]{turnover_flowchart_1.png}
168	!! \includegraphics[scale=0.5]{turnover_flowchart_2.png}
169	!! \endlatexonly
170	!! \n
171	!_ ================================================================================================================================
172
173	SUBROUTINE turnover_diagnostic (npts, dt, PFTpresent, herbivores, maxmoiavail_lastyear, &
174	minmoiavail_lastyear, moiavail_week, moiavail_month, tlong_ref, t2m_month, &
175	t2m_week, veget_max, gdd_from_growthinit, leaf_age, leaf_frac, &
176	age, lai, biomass, turnover, senescence, &
177	turnover_time, forest_managed)
178
179	!! 0. Variable and parameter declaration
180
181	!! 0.1 Input variables
182
183	INTEGER(i_std), INTENT(in) :: npts !! Domain size - number of grid cells
184	!! (unitless)
185	INTEGER(i_std), DIMENSION (npts,nvm), INTENT(in) :: forest_managed !! forest management flag
186	REAL(r_std), INTENT(in) :: dt !! time step (dt_days)
187	LOGICAL, DIMENSION(npts,nvm), INTENT(in) :: PFTpresent !! PFT exists (true/false)
188	REAL(r_std), DIMENSION(npts), INTENT(in) :: tlong_ref !! "long term" 2 meter reference
189	!! temperatures (K)
190	REAL(r_std), DIMENSION(npts), INTENT(in) :: t2m_month !! "monthly" 2-meter temperatures (K)
191	REAL(r_std), DIMENSION(npts), INTENT(in) :: t2m_week !! "weekly" 2 meter temperatures (K)
192	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: herbivores !! time constant of probability of a leaf to
193	!! be eaten by a herbivore (days)
194	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: maxmoiavail_lastyear !! last year's maximum moisture availability
195	!! (0-1, unitless)
196	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: minmoiavail_lastyear !! last year's minimum moisture availability
197	!! (0-1, unitless)
198	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: moiavail_week !! "weekly" moisture availability
199	!! (0-1, unitless)
200	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: moiavail_month !! "monthly" moisture availability
201	!! (0-1, unitless)
202	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: veget_max !! "maximal" coverage fraction of a PFT (LAI
203	!! -> infinity) on ground (unitless)
204	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: gdd_from_growthinit !! gdd senescence for crop
205
206	!! 0.2 Output variables
207
208	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), INTENT(out) :: turnover !! Turnover @tex ($gC m^{-2}$) @endtex
209	LOGICAL, DIMENSION(npts,nvm), INTENT(out) :: senescence !! is the plant senescent? (true/false)
210	!! (interesting only for deciduous trees:
211	!! carbohydrate reserve)
212	!! 0.3 Modified variables
213
214	REAL(r_std), DIMENSION(npts,nvm,nleafages), INTENT(inout) :: leaf_age !! age of the leaves (days)
215	REAL(r_std), DIMENSION(npts,nvm,nleafages), INTENT(inout) :: leaf_frac !! fraction of leaves in leaf age class
216	!! (0-1, unitless)
217	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: age !! age (years)
218	REAL(r_std), DIMENSION(npts,nvm), INTENT(in) :: lai !! leaf area index @tex ($m^2 m^{-2}$)
219	!! @endtex
220	REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(inout) :: turnover_time !! turnover_time of grasses (days)
221	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), &
222	INTENT(inout) :: biomass !! biomass @tex ($gC m^{-2}$) @endtex
223
224	!! 0.4 Local variables
225
226	REAL(r_std), DIMENSION(npts,nvm) :: leaf_meanage !! mean age of the leaves (days)
227	REAL(r_std), DIMENSION(npts) :: dturnover !! Intermediate variable for turnover ??
228	!! @tex ($gC m^{-2}$) @endtex
229	REAL(r_std), DIMENSION(npts) :: moiavail_crit !! critical moisture availability, function
230	!! of last year's moisture availability
231	!! (0-1, unitless)
232	REAL(r_std), DIMENSION(npts) :: tl !! long term annual mean temperature, (C)
233	REAL(r_std), DIMENSION(npts) :: t_crit !! critical senescence temperature, function
234	!! of long term annual temperature (K)
235	LOGICAL, DIMENSION(npts) :: shed_rest !! shed the remaining leaves? (true/false)
236	REAL(r_std), DIMENSION(npts) :: sapconv !! Sapwood conversion @tex ($gC m^{-2}$)
237	!! @endtex
238	REAL(r_std), DIMENSION(npts) :: hw_old !! old heartwood mass @tex ($gC m^{-2}$)
239	!! @endtex
240	REAL(r_std), DIMENSION(npts) :: hw_new !! new heartwood mass @tex ($gC m^{-2}$)
241	!! @endtex
242	REAL(r_std), DIMENSION(npts) :: lm_old !! old leaf mass @tex ($gC m^{-2}$) @endtex
243	REAL(r_std), DIMENSION(npts,nleafages) :: delta_lm !! leaf mass change for each age class @tex
244	!! ($gC m^{-2}$) @endtex
245	REAL(r_std), DIMENSION(npts) :: turnover_rate !! turnover rate (unitless)
246	REAL(r_std), DIMENSION(npts,nvm) :: leaf_age_crit !! critical leaf age (days)
247	REAL(r_std), DIMENSION(npts,nvm) :: new_turnover_time !! instantaneous turnover time (days)
248	INTEGER(i_std) :: j,m,k !! Index (unitless)
249	REAL(r_std) :: branch_turn !! turnover of branches (how much goes
250	!! to litter each day) (cf article CO2fix)
251	REAL(r_std), SAVE :: ss_branch_turn !! Variations for sensitivity analysis
252	!$OMP THREADPRIVATE(ss_branch_turn)
253	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements) :: bm_old !! old root mass
254	REAL(r_std), DIMENSION(npts,nvm) :: sapwood_age !! sapwood age (days)
255	REAL(r_std), DIMENSION(npts) :: sw_old !! old sapwood mass
256	!! (gC/(m**2 of nat/agri ground))
257	REAL(r_std), DIMENSION(npts) :: sw_new !! new sapwood mass
258	!! (gC/(m**2 of nat/agri ground))
259	INTEGER(i_std), SAVE :: turn_count !!
260	!$OMP THREADPRIVATE(turn_count)
261	!_ ================================================================================================================================
262
263	IF (bavard.GE.3) WRITE(numout,*) 'Entering turnover'
264
265	!! 1. first call - output messages
266
267	IF ( firstcall ) THEN
268
269	WRITE(numout,*) 'turnover:'
270
271	WRITE(numout,*) ' > minimum mean leaf age for senescence (days) (::min_leaf_age_for_senescence) : '&
272	,min_leaf_age_for_senescence
273
274	turn_count = 0
275
276	!Config Key = SS_BRANCH_TURN
277	!Config Desc =
278	!Config If =
279	!Config Def =
280	!Config Help =
281	!Config Units =
282	!!!!!!!! needs a default value
283	ss_branch_turn=un
284	CALL getin_p ('SS_BRANCH_TURN',ss_branch_turn)
285
286	firstcall = .FALSE.
287
288
289	ENDIF
290
291	!! 2. Initializations
292
293	!! 2.1 set output to zero
294	turnover(:,:,:,:) = zero
295	new_turnover_time(:,:) = zero
296	senescence(:,:) = .FALSE.
297	sapwood_age(:,:)=un
298
299	! save old biomass
300	bm_old(:,:,:,:) = biomass(:,:,:,:)
301	! Compute the turnover of branches
302	! (2.5% per year : Orchidee standard, Masera 2003, Lehtonen 2004, DeAngelis 1981)
303	branch_turn = 0.025/(one_year/dt)*ss_branch_turn
304
305	!! 2.2 Recalculate mean leaf age
306	! Mean leaf age is recalculated from leaf ages of separate leaf age classes. Should actually be recalculated at the
307	! end of this routine, but it does not change too fast.
308	! The mean leaf age is calculated using the following equation:
309	! \latexonly
310	! \input{turnover_lma_update_eqn1.tex}
311	! \endlatexonly
312	! \n
313	leaf_meanage(:,:) = zero
314
315	DO m = 1, nleafages
316	leaf_meanage(:,:) = leaf_meanage(:,:) + leaf_age(:,:,m) * leaf_frac(:,:,m)
317	ENDDO
318
319
320	IF (ld_trnov) THEN
321	WRITE(numout,*) 'leaf_meanage in turnover', leaf_meanage(:,test_pft)
322	WRITE(numout,*) 'leaf_age in turnover', leaf_age(:,test_pft,:)
323	WRITE(numout,*) 'leaf_frac in turnover', leaf_frac(:,test_pft,:)
324	ENDIF
325
326	CALL histwrite_p (hist_id_stomate, 'LEAFAGE', itime, &
327	leaf_meanage, npts*nvm, horipft_index)
328
329
330	!! 3. Climatic senescence
331
332	! Three different types of "climatic" leaf senescence,
333	! that do not change the age structure.
334	DO j = 2,nvm ! Loop over # PFTs
335
336	!! 3.1 Determine if there is climatic senescence.
337	! The climatic senescence can be of three types:
338	! sensitivity to cold temperatures, to lack of water, or both. If meteorological conditions are
339	! fulfilled, a flag senescence is set to TRUE.
340	! Evergreen species do not experience climatic senescence.
341
342	SELECT CASE ( senescence_type(j) )
343
344
345	CASE ('crop' )!for crop senescence is based on a GDD criterium as in crop models
346	WHERE ( ( biomass(:,j,ileaf,icarbon) .GT. zero ) .AND. &
347	( leaf_meanage(:,j) .GT. min_leaf_age_for_senescence(j) ) .AND. &
348	( gdd_from_growthinit(:,j) .GT. gdd_senescence(j)))
349	senescence(:,j) = .TRUE.
350	ENDWHERE
351
352	CASE ( 'cold' )
353
354	!! 3.1.1 Summergreen species: Climatic senescence is triggered by sensitivity to cold temperatures
355	! Climatic senescence is triggered by sensitivity to cold temperatures as follows:
356	! If biomass is large enough (i.e. when it is greater than zero),
357	! AND (i.e. when leaf mean age is above a certain PFT-dependent treshold ::min_leaf_age_for_senescence,
358	! which is given in ::stomate_constants),
359	! AND the monthly temperature is low enough (i.e. when monthly temperature ::t2m_month(:) is below a critical
360	! temperature ::t_crit(:), which is calculated in this module),
361	! AND the temperature tendency is negative (i.e. when weekly temperatures ::t2m_week(:) are lower than monthly
362	! temperatures ::t2m_month(:))
363	! If these conditions are met, senescence is set to TRUE.
364	!
365	! The critical temperature for senescence is calculated using the following equation:
366	! \latexonly
367	! \input{turnover_temp_crit_eqn2.tex}
368	! \endlatexonly
369	! \n
370	!
371	! Critical temperature for senescence may depend on long term annual mean temperature
372	tl(:) = tlong_ref(:) - ZeroCelsius
373	t_crit(:) = ZeroCelsius + senescence_temp(j,1) + &
374	tl(:) * senescence_temp(j,2) + &
375	tl(:)tl(:) senescence_temp(j,3)
376
377	WHERE ( ( biomass(:,j,ileaf,icarbon) .GT. zero ) .AND. &
378	( leaf_meanage(:,j) .GT. min_leaf_age_for_senescence(j) ) .AND. &
379	( t2m_month(:) .LT. t_crit(:) ) .AND. ( t2m_week(:) .LT. t2m_month(:) ) )
380
381	senescence(:,j) = .TRUE.
382
383	ENDWHERE
384
385	CASE ( 'dry' )
386
387	!! 3.1.2 Raingreen species: Climatic senescence is triggered by sensitivity to lack of water availability
388	! Climatic senescence is triggered by sensitivity to lack of water availability as follows:
389	! If biomass is large enough (i.e. when it is greater than zero),
390	! AND (i.e. when leaf mean age is above a certain PFT-dependent treshold ::min_leaf_age_for_senescence,
391	! which is given in ::stomate_constants),
392	! AND the moisture availability drops below a critical level (i.e. when weekly moisture availability
393	! ::moiavail_week(:,j) is below a critical moisture availability ::moiavail_crit(:),
394	! which is calculated in this module),
395	! If these conditions are met, senescence is set to TRUE.
396	!
397	! The moisture availability critical level is calculated using the following equation:
398	! \latexonly
399	! \input{turnover_moist_crit_eqn3.tex}
400	! \endlatexonly
401	! \n
402	moiavail_crit(:) = &
403	MIN( MAX( minmoiavail_lastyear(:,j) + hum_frac(j) * &
404	( maxmoiavail_lastyear(:,j) - minmoiavail_lastyear(:,j) ), &
405	senescence_hum(j) ), &
406	nosenescence_hum(j) )
407
408	WHERE ( ( biomass(:,j,ileaf,icarbon) .GT. zero ) .AND. &
409	( leaf_meanage(:,j) .GT. min_leaf_age_for_senescence(j) ) .AND. &
410	( moiavail_week(:,j) .LT. moiavail_crit(:) ) )
411
412	senescence(:,j) = .TRUE.
413
414	ENDWHERE
415
416	CASE ( 'mixed' )
417
418	!! 3.1.3 Mixed criterion: Climatic senescence is triggered by sensitivity to temperature or to lack of water
419	! Climatic senescence is triggered by sensitivity to temperature or to lack of water availability as follows:
420	! If biomass is large enough (i.e. when it is greater than zero),
421	! AND (i.e. when leaf mean age is above a certain PFT-dependent treshold ::min_leaf_age_for_senescence,
422	! which is given in ::stomate_constants),
423	! AND the moisture availability drops below a critical level (i.e. when weekly moisture availability
424	! ::moiavail_week(:,j) is below a critical moisture availability ::moiavail_crit(:), calculated in this module),
425	! OR
426	! the monthly temperature is low enough (i.e. when monthly temperature ::t2m_month(:) is below a critical
427	! temperature ::t_crit(:), calculated in this module),
428	! AND the temperature tendency is negative (i.e. when weekly temperatures ::t2m_week(:) are lower than
429	! monthly temperatures ::t2m_month(:)).
430	! If these conditions are met, senescence is set to TRUE.
431	moiavail_crit(:) = &
432	MIN( MAX( minmoiavail_lastyear(:,j) + hum_frac(j) * &
433	(maxmoiavail_lastyear(:,j) - minmoiavail_lastyear(:,j) ), &
434	senescence_hum(j) ), &
435	nosenescence_hum(j) )
436
437	tl(:) = tlong_ref(:) - ZeroCelsius
438	t_crit(:) = ZeroCelsius + senescence_temp(j,1) + &
439	tl(:) * senescence_temp(j,2) + &
440	tl(:)tl(:) senescence_temp(j,3)
441
442	IF ( is_tree(j) ) THEN
443	! critical temperature for senescence may depend on long term annual mean temperature
444	WHERE ( ( biomass(:,j,ileaf,icarbon) .GT. zero ) .AND. &
445	( leaf_meanage(:,j) .GT. min_leaf_age_for_senescence(j) ) .AND. &
446	( ( moiavail_week(:,j) .LT. moiavail_crit(:) ) .OR. &
447	( ( t2m_month(:) .LT. t_crit(:) ) .AND. ( t2m_week(:) .LT. t2m_month(:) ) ) ) )
448	senescence(:,j) = .TRUE.
449	ENDWHERE
450
451	ELSE
452
453	! The nparts dimension was introduced for the prognstic part. Use dimension 1 so that
454	! no new variable needs to be defined
455	turnover_time(:,j,1) = max_turnover_time(j)
456
457	WHERE ( ( biomass(:,j,ileaf,icarbon) .GT. zero ) .AND. &
458	( leaf_meanage(:,j) .GT. min_leaf_age_for_senescence(j) ) .AND. &
459	( ( moiavail_week(:,j) .LT. moiavail_crit(:) )))
460	turnover_time(:,j,1) = max_turnover_time(j) * &
461	(1.- (1.- (moiavail_week(:,j)/ moiavail_crit(:)))**2)
462	ENDWHERE
463
464	WHERE ( turnover_time(:,j,1) .LT. min_turnover_time(j) )
465	turnover_time(:,j,1) = min_turnover_time(j)
466	ENDWHERE
467
468	WHERE ((( biomass(:,j,ileaf,icarbon) .GT. zero ) .AND. &
469	( leaf_meanage(:,j) .GT. min_leaf_age_for_senescence(j) ) .AND. &
470	((t2m_month(:) .LT. t_crit(:)) .AND. (lai(:,j) .GT. lai_max(j)/4.) .OR. &
471	(t2m_month(:) .LT. ZeroCelsius)) .AND. ( t2m_week(:) .LT. t2m_month(:) )))
472	turnover_time(:,j,1)= leaffall(j)
473	ENDWHERE
474
475	ENDIF
476
477
478	!! Evergreen species do not experience climatic senescence
479	CASE ( 'none' )
480
481
482	!! In case no climatic senescence type is recognized.
483	CASE default
484
485	WRITE(numout,*) ' turnover: don''t know how to treat this PFT.'
486	WRITE(numout,*) ' number (::j) : ',j
487	WRITE(numout,*) ' senescence type (::senescence_type(j)) : ',senescence_type(j)
488
489	STOP
490
491	END SELECT
492
493	!! 3.2 Drop leaves and roots, plus stems and fruits for grasses
494
495	IF ( is_tree(j) ) THEN
496
497	!! 3.2.1 Trees in climatic senescence lose their fine roots at the same rate as they lose their leaves.
498	! The rate of biomass loss of both fine roots and leaves is presribed through the equation:
499	! \latexonly
500	! \input{turnover_clim_senes_biomass_eqn4.tex}
501	! \endlatexonly
502	! \n
503	! with ::leaffall(j) a PFT-dependent time constant which is given in ::stomate_constants),
504	WHERE ( senescence(:,j) )
505
506	turnover(:,j,ileaf,icarbon) = biomass(:,j,ileaf,icarbon) * dt / leaffall(j)
507	turnover(:,j,iroot,icarbon) = biomass(:,j,iroot,icarbon) * dt / leaffall(j)
508
509	ENDWHERE
510
511	ELSE
512
513	!! 3.2.2 In grasses, leaf senescence is extended to the whole plant
514	! In grasses, leaf senescence is extended to the whole plant (all carbon pools) except to its
515	! carbohydrate reserve.
516
517	IF (senescence_type(j) .EQ. 'crop') THEN
518	! 3.2.2.1 crops with 'crop' phenological model
519	WHERE ( senescence(:,j) )
520	turnover(:,j,ileaf,icarbon) = biomass(:,j,ileaf,icarbon) * dt / leaffall(j)
521	turnover(:,j,iroot,icarbon) = biomass(:,j,iroot,icarbon) * dt / leaffall(j)
522	turnover(:,j,isapabove,icarbon) = biomass(:,j,isapabove,icarbon) * dt / leaffall(j)
523	turnover(:,j,ifruit,icarbon) = biomass(:,j,ifruit,icarbon) * dt /leaffall(j)
524	ENDWHERE
525	ELSE
526
527	! 3.2.2.2 grass or crops based on 'mixed' phenological model
528	WHERE (turnover_time(:,j,1) .LT. max_turnover_time(j))
529	turnover(:,j,ileaf,icarbon) = biomass(:,j,ileaf,icarbon) * dt / turnover_time(:,j,1)
530	turnover(:,j,isapabove,icarbon) = biomass(:,j,isapabove,icarbon) * dt / turnover_time(:,j,1)
531	turnover(:,j,iroot,icarbon) = biomass(:,j,iroot,icarbon) * dt / turnover_time(:,j,1)
532	turnover(:,j,ifruit,icarbon) = biomass(:,j,ifruit,icarbon) * dt / turnover_time(:,j,1)
533	ENDWHERE
534	ENDIF
535	ENDIF ! tree/grass
536
537	biomass(:,j,ileaf,icarbon) = biomass(:,j,ileaf,icarbon) - turnover(:,j,ileaf,icarbon)
538	biomass(:,j,isapabove,icarbon) = biomass(:,j,isapabove,icarbon) - turnover(:,j,isapabove,icarbon)
539	biomass(:,j,iroot,icarbon) = biomass(:,j,iroot,icarbon) - turnover(:,j,iroot,icarbon)
540	biomass(:,j,ifruit,icarbon) = biomass(:,j,ifruit,icarbon) - turnover(:,j,ifruit,icarbon)
541
542	ENDDO ! loop over PFTs
543
544	!! 4. Leaf fall
545	! At a certain age, leaves fall off, even if the climate would allow a green plant
546	! all year round. Even if the meteorological conditions are favorable for leaf maintenance,
547	! plants, and in particular, evergreen trees, have to renew their leaves simply because the
548	! old leaves become inefficient.
549	! Roots, fruits (and stems) follow leaves. The decay rate varies with leaf age.
550	! Note that plant is not declared senescent in this case (wchich is important for allocation:
551	! if the plant loses leaves because of their age, it can renew them).
552	!
553	! The leaf turnover rate due to aging of leaves is calculated using the following equation:
554	! \latexonly
555	! \input{turnover_age_senes_biomass_eqn5.tex}
556	! \endlatexonly
557	! \n
558	DO j = 2,nvm ! Loop over # PFTs
559
560	!! save old leaf mass
561	lm_old(:) = biomass(:,j,ileaf,icarbon)
562
563	!! initialize leaf mass change in age class
564	delta_lm(:,:) = zero
565
566	IF ( is_tree(j) .OR. (.NOT. natural(j)) ) THEN
567
568	!! 4.1 Trees: leaves, roots, fruits roots and fruits follow leaves.
569
570	!! 4.1.1 Critical age: prescribed for trees
571	leaf_age_crit(:,j) = tau_leaf(j)
572
573	ELSE
574
575	!! 4.2 Grasses: leaves, roots, fruits, sap follow leaves.
576
577	!! 4.2.1 Critical leaf age depends on long-term temperature
578	! Critical leaf age depends on long-term temperature
579	! generally, lower turnover in cooler climates.
580	leaf_age_crit(:,j) = &
581	MIN( tau_leaf(j) * leaf_age_crit_coeff(1) , &
582	MAX( tau_leaf(j) * leaf_age_crit_coeff(2) , &
583	tau_leaf(j) - leaf_age_crit_coeff(3) * &
584	( tlong_ref(:)-ZeroCelsius - leaf_age_crit_tref ) ) )
585
586	END IF
587
588	! 4.2.2 Loop over leaf age classes
589	DO m = 1, nleafages
590
591	turnover_rate(:) = zero
592
593	WHERE ( leaf_age(:,j,m) .GT. leaf_age_crit(:,j)/2. )
594
595	turnover_rate(:) = &
596	MIN( 0.99_r_std, dt / ( leaf_age_crit(:,j) * &
597	( leaf_age_crit(:,j) / leaf_age(:,j,m) )**quatre ) )
598
599	ENDWHERE
600
601	dturnover(:) = biomass(:,j,ileaf,icarbon) * leaf_frac(:,j,m) * turnover_rate(:)
602	turnover(:,j,ileaf,icarbon) = turnover(:,j,ileaf,icarbon) + dturnover(:)
603	biomass(:,j,ileaf,icarbon) = biomass(:,j,ileaf,icarbon) - dturnover(:)
604
605	! save leaf mass change
606	delta_lm(:,m) = - dturnover(:)
607
608	dturnover(:) = biomass(:,j,iroot,icarbon) * leaf_frac(:,j,m) * turnover_rate(:)
609	turnover(:,j,iroot,icarbon) = turnover(:,j,iroot,icarbon) + dturnover(:)
610	biomass(:,j,iroot,icarbon) = biomass(:,j,iroot,icarbon) - dturnover(:)
611
612	dturnover(:) = biomass(:,j,ifruit,icarbon) * leaf_frac(:,j,m) * turnover_rate(:)
613	turnover(:,j,ifruit,icarbon) = turnover(:,j,ifruit,icarbon) + dturnover(:)
614	biomass(:,j,ifruit,icarbon) = biomass(:,j,ifruit,icarbon) - dturnover(:)
615
616	IF (.NOT. is_tree(j)) THEN
617	dturnover(:) = biomass(:,j,isapabove,icarbon) * leaf_frac(:,j,m) * turnover_rate(:)
618	turnover(:,j,isapabove,icarbon) = turnover(:,j,isapabove,icarbon) + dturnover(:)
619	biomass(:,j,isapabove,icarbon) = biomass(:,j,isapabove,icarbon) - dturnover(:)
620	ENDIF
621
622	ENDDO
623
624
625
626	!! 4.3 Recalculate the fraction of leaf biomass in each leaf age class.
627	! Older leaves will fall more fast than younger leaves and therefore
628	! the leaf age distribution needs to be recalculated after turnover.
629	! The fraction of biomass in each leaf class is updated using the following equation:
630	! \latexonly
631	! \input{turnover_update_LeafAgeDistribution_eqn6.tex}
632	! \endlatexonly
633	! \n
634	!
635	! new fraction = new leaf mass of that fraction / new total leaf mass
636	! = (old fraction*old total leaf mass ::lm_old(:) + biomass change of that fraction ::delta_lm(:,m) ) /
637	! new total leaf mass ::biomass(:,j,ileaf
638	DO m = 1, nleafages
639
640	WHERE ( biomass(:,j,ileaf,icarbon) .GT. zero )
641
642	leaf_frac(:,j,m) = ( leaf_frac(:,j,m)*lm_old(:) + delta_lm(:,m) ) / biomass(:,j,ileaf,icarbon)
643
644	ELSEWHERE
645
646	leaf_frac(:,j,m) = zero
647
648	ENDWHERE
649
650	ENDDO
651
652	ENDDO ! loop over PFTs
653
654	!! 5. New (provisional) LAI
655	! ::lai(:,j) is determined from the leaf biomass ::biomass(:,j,ileaf,icarbon) and the
656	! specific leaf surface :: sla(j) (m^2 gC^{-1})
657	! The leaf area index is updated using the following equation:
658	! \latexonly
659	! \input{turnover_update_LAI_eqn7.tex}
660	! \endlatexonly
661	! \n
662
663	! lai(:,ibare_sechiba) = zero
664	! DO j = 2, nvm ! Loop over # PFTs
665	! lai(:,j) = biomass(:,j,ileaf,icarbon) * sla(j)
666	! ENDDO
667
668	!! 6. Definitely drop all leaves if there is a very low leaf mass during senescence.
669
670	! Both for deciduous trees and grasses same conditions are checked:
671	! If biomass is large enough (i.e. when it is greater than zero),
672	! AND when senescence is set to true
673	! AND the leaf biomass drops below a critical minimum biomass level (i.e. when it is lower than half
674	! the minimum initial LAI ::lai_initmin(j) divided by the specific leaf area ::sla(j),
675	! ::lai_initmin(j) is set to 0.3 in stomate_data.f90 and sla is a constant that is set to 0.015366 m2/gC),
676	! If these conditions are met, the flag ::shed_rest(:) is set to TRUE.
677	!
678	! After this, the biomass of different carbon pools both for trees and grasses is set to zero
679	! and the mean leaf age is reset to zero.
680	! Finally, the leaf fraction and leaf age of the different leaf age classes is set to zero.
681	DO j = 2,nvm ! Loop over # PFTs
682
683	shed_rest(:) = .FALSE.
684
685	!! 6.1 For deciduous trees: next to leaves, also fruits and fine roots are dropped
686	! For deciduous trees: next to leaves, also fruits and fine roots are dropped: fruit ::biomass(:,j,ifruit)
687	! and fine root ::biomass(:,j,iroot) carbon pools are set to zero.
688	IF ( is_tree(j) .AND. ( senescence_type(j) .NE. 'none' ) ) THEN
689
690	! check whether we shed the remaining leaves
691	WHERE ( ( biomass(:,j,ileaf,icarbon) .GT. zero ) .AND. senescence(:,j) .AND. &
692	( biomass(:,j,ileaf,icarbon) .LT. (lai_initmin(j) / 2.)/sla(j) ) )
693
694	shed_rest(:) = .TRUE.
695
696	turnover(:,j,ileaf,icarbon) = turnover(:,j,ileaf,icarbon) + biomass(:,j,ileaf,icarbon)
697	turnover(:,j,iroot,icarbon) = turnover(:,j,iroot,icarbon) + biomass(:,j,iroot,icarbon)
698	turnover(:,j,ifruit,icarbon) = turnover(:,j,ifruit,icarbon) + biomass(:,j,ifruit,icarbon)
699
700	biomass(:,j,ileaf,icarbon) = zero
701	biomass(:,j,iroot,icarbon) = zero
702	biomass(:,j,ifruit,icarbon) = zero
703
704	! reset leaf age
705	leaf_meanage(:,j) = zero
706
707	ENDWHERE
708
709	ENDIF
710
711	!! 6.2 For grasses: all aboveground carbon pools, except the carbohydrate reserves are affected:
712	! For grasses: all aboveground carbon pools, except the carbohydrate reserves are affected:
713	! fruit ::biomass(:,j,ifruit,icarbon), fine root ::biomass(:,j,iroot,icarbon) and sapwood above
714	! ::biomass(:,j,isapabove,icarbon) carbon pools are set to zero.
715	IF ( .NOT. is_tree(j) ) THEN
716
717	! Shed the remaining leaves if LAI very low.
718	WHERE ( ( biomass(:,j,ileaf,icarbon) .GT. zero ) .AND. senescence(:,j) .AND. &
719	( biomass(:,j,ileaf,icarbon) .LT. (lai_initmin(j) / 2.)/sla(j) ))
720
721	shed_rest(:) = .TRUE.
722
723	turnover(:,j,ileaf,icarbon) = turnover(:,j,ileaf,icarbon) + biomass(:,j,ileaf,icarbon)
724	turnover(:,j,isapabove,icarbon) = turnover(:,j,isapabove,icarbon) + biomass(:,j,isapabove,icarbon)
725	turnover(:,j,iroot,icarbon) = turnover(:,j,iroot,icarbon) + biomass(:,j,iroot,icarbon)
726	turnover(:,j,ifruit,icarbon) = turnover(:,j,ifruit,icarbon) + biomass(:,j,ifruit,icarbon)
727
728	biomass(:,j,ileaf,icarbon) = zero
729	biomass(:,j,isapabove,icarbon) = zero
730	biomass(:,j,iroot,icarbon) = zero
731	biomass(:,j,ifruit,icarbon) = zero
732
733	! reset leaf age
734	leaf_meanage(:,j) = zero
735
736	ENDWHERE
737
738	ENDIF
739
740	!! 6.3 Reset the leaf age structure: the leaf fraction and leaf age of the different leaf age classes is set to zero.
741
742	DO m = 1, nleafages
743
744	WHERE ( shed_rest(:) )
745
746	leaf_age(:,j,m) = zero
747	leaf_frac(:,j,m) = zero
748
749	ENDWHERE
750
751	ENDDO
752
753
754
755	!! 6.4 drop part of the branches and mortality of smaller trees until average height reaches 10 m
756	!! (herbivory, water shortage, ...). This "youth mortality" does not affect tree density as young
757	!! dead trees are immediately replaced.
758	!
759	IF ( control%forest_management ) THEN
760
761	IF ( is_tree(j) ) THEN
762
763	WHERE (forest_managed(:,j) > 0)
764
765	! Drop part of the branches
766	turnover(:,j,isapabove,icarbon) = branch_ratio(j)branch_turnbiomass(:,j,isapabove,icarbon)
767	turnover(:,j,iheartabove,icarbon) = branch_ratio(j)branch_turnbiomass(:,j,iheartabove,icarbon)
768	turnover(:,j,isapbelow,icarbon) = branch_ratio(j)branch_turnbiomass(:,j,isapbelow,icarbon)
769	turnover(:,j,iheartbelow,icarbon) = branch_ratio(j)branch_turnbiomass(:,j,iheartbelow,icarbon)
770
771	biomass(:,j,isapabove,icarbon) = biomass(:,j,isapabove,icarbon) - turnover(:,j,isapabove,icarbon)
772	biomass(:,j,iheartabove,icarbon) = biomass(:,j,iheartabove,icarbon) - turnover(:,j,iheartabove,icarbon)
773	biomass(:,j,isapbelow,icarbon) = biomass(:,j,isapbelow,icarbon) - turnover(:,j,isapbelow,icarbon)
774	biomass(:,j,iheartbelow,icarbon) = biomass(:,j,iheartbelow,icarbon) - turnover(:,j,iheartbelow,icarbon)
775
776	END WHERE
777
778	END IF
779
780	END IF ! (control%forest_management)
781
782	ENDDO ! loop over PFTs
783
784	!! 7. Herbivore activity: elephants, cows, gazelles but no lions.
785
786	! Herbivore activity affects the biomass of leaves and fruits as well
787	! as stalks (only for grasses). Herbivore activity does not modify leaf
788	! age structure. Herbivores ::herbivores(:,j) is the time constant of
789	! probability of a leaf to be eaten by a herbivore, and is calculated in
790	! ::stomate_season. following Mc Naughton et al. [1989].
791
792	IF ( control%ok_herbivory ) THEN
793
794	! If the herbivore activity is allowed (if ::control%ok_herbivory is true, which is set in run.def),
795	! remove the amount of biomass consumed by herbivory from the leaf biomass ::biomass(:,j,ileaf,icarbon) and
796	! the fruit biomass ::biomass(:,j,ifruit,icarbon).
797	! The daily amount consumed equals the biomass multiplied by 1 day divided by the time constant ::herbivores(:,j).
798	DO j = 2,nvm ! Loop over # PFTs
799
800	IF ( is_tree(j) ) THEN
801
802	!! For trees: only the leaves and fruit carbon pools are affected
803
804	WHERE (biomass(:,j,ileaf,icarbon) .GT. zero)
805
806	! added by shilong
807	WHERE (herbivores(:,j).GT. zero)
808
809	dturnover(:) = biomass(:,j,ileaf,icarbon) * dt / herbivores(:,j)
810	turnover(:,j,ileaf,icarbon) = turnover(:,j,ileaf,icarbon) + dturnover(:)
811	biomass(:,j,ileaf,icarbon) = biomass(:,j,ileaf,icarbon) - dturnover(:)
812
813	dturnover(:) = biomass(:,j,ifruit,icarbon) * dt / herbivores(:,j)
814	turnover(:,j,ifruit,icarbon) = turnover(:,j,ifruit,icarbon) + dturnover(:)
815	biomass(:,j,ifruit,icarbon) = biomass(:,j,ifruit,icarbon) - dturnover(:)
816
817	ENDWHERE
818
819	ENDWHERE
820
821	ELSE
822
823	! For grasses: all aboveground carbon pools are affected: leaves, fruits and sapwood above
824	WHERE ( biomass(:,j,ileaf,icarbon) .GT. zero )
825
826	! added by shilong
827	WHERE (herbivores(:,j) .GT. zero)
828
829	dturnover(:) = biomass(:,j,ileaf,icarbon) * dt / herbivores(:,j)
830	turnover(:,j,ileaf,icarbon) = turnover(:,j,ileaf,icarbon) + dturnover(:)
831	biomass(:,j,ileaf,icarbon) = biomass(:,j,ileaf,icarbon) - dturnover(:)
832
833	dturnover(:) = biomass(:,j,isapabove,icarbon) * dt / herbivores(:,j)
834	turnover(:,j,isapabove,icarbon) = turnover(:,j,isapabove,icarbon) + dturnover(:)
835	biomass(:,j,isapabove,icarbon) = biomass(:,j,isapabove,icarbon) - dturnover(:)
836
837	dturnover(:) = biomass(:,j,ifruit,icarbon) * dt / herbivores(:,j)
838	turnover(:,j,ifruit,icarbon) = turnover(:,j,ifruit,icarbon) + dturnover(:)
839	biomass(:,j,ifruit,icarbon) = biomass(:,j,ifruit,icarbon) - dturnover(:)
840
841	ENDWHERE
842
843	ENDWHERE
844
845	ENDIF ! tree/grass?
846
847	ENDDO ! loop over PFTs
848
849	ENDIF ! end herbivores
850
851	!! 8. Fruit turnover for trees
852
853	! Fruit turnover for trees: trees simply lose their fruits with a time constant ::tau_fruit(j),
854	! that is set to 90 days for all PFTs in ::stomate_constants
855
856	DO k = 1,nelements
857
858	DO j = 2,nvm ! Loop over # PFTs
859
860	IF ( is_tree(j) ) THEN
861
862	dturnover(:) = biomass(:,j,ifruit,k) * dt / tau_fruit(j)
863	turnover(:,j,ifruit,k) = turnover(:,j,ifruit,k) + dturnover(:)
864	biomass(:,j,ifruit,k) = biomass(:,j,ifruit,k) - dturnover(:)
865
866	ENDIF
867
868	ENDDO ! loop over PFTs
869
870	END DO
871
872	!! 9 Conversion of sapwood to heartwood both for aboveground and belowground sapwood and heartwood.
873
874	! Following LPJ (Sitch et al., 2003), sapwood biomass is converted into heartwood biomass
875	! with a time constant tau ::tau_sap(j) of 1 year.
876	! Note that this biomass conversion is not added to "turnover" as the biomass is not lost!
877	DO j = 2,nvm ! Loop over # PFTs
878
879	IF ( is_tree(j) ) THEN
880
881	!! For the recalculation of age in 9.2 (in case the vegetation is not dynamic ie. ::control%ok_dgvm is FALSE),
882	!! the heartwood above and below is stored in ::hw_old(:).
883	IF ( .NOT. control%ok_dgvm ) THEN
884
885	hw_old(:) = biomass(:,j,iheartabove,icarbon) + biomass(:,j,iheartbelow,icarbon)
886	sw_old(:) = biomass(:,j,isapabove,icarbon) + biomass(:,j,isapbelow,icarbon)
887
888	ENDIF
889
890	!! 9.1 Calculate the rate of sapwood to heartwood conversion
891	! Calculate the rate of sapwood to heartwood conversion with the time constant ::tau_sap(j)
892	! and update aboveground and belowground sapwood ::biomass(:,j,isapabove) and ::biomass(:,j,isapbelow)
893	! and heartwood ::biomass(:,j,iheartabove) and ::biomass(:,j,iheartbelow).
894
895	DO k = 1,nelements
896
897	! Above the ground
898	sapconv(:) = biomass(:,j,isapabove,k) * dt / tau_sap(j)
899	biomass(:,j,isapabove,k) = biomass(:,j,isapabove,k) - sapconv(:)
900	biomass(:,j,iheartabove,k) = biomass(:,j,iheartabove,k) + sapconv(:)
901
902	! Below the ground
903	sapconv(:) = biomass(:,j,isapbelow,k) * dt / tau_sap(j)
904	biomass(:,j,isapbelow,k) = biomass(:,j,isapbelow,k) - sapconv(:)
905	biomass(:,j,iheartbelow,k) = biomass(:,j,iheartbelow,k) + sapconv(:)
906
907	END DO
908
909	!! 9.2 If the vegetation is not dynamic, the age of the plant is decreased.
910	! The updated heartwood, the sum of new heartwood above and new heartwood below after
911	! converting sapwood to heartwood, is saved as ::hw_new(:) .
912	! Creation of new heartwood decreases the age of the plant with a factor that is determined by:
913	! old heartwood ::hw_old(:) divided by the new heartwood ::hw_new(:)
914	IF ( .NOT. control%ok_dgvm ) THEN
915
916	hw_new(:) = biomass(:,j,iheartabove,icarbon) + biomass(:,j,iheartbelow,icarbon)
917	sw_new(:) = biomass(:,j,isapabove,icarbon) + biomass(:,j,isapbelow,icarbon)
918
919	WHERE ( hw_new(:) .GT. zero )
920
921	age(:,j) = age(:,j) * hw_old(:)/hw_new(:)
922	sapwood_age(:,j) = sapwood_age(:,j) * sw_old(:)/sw_new(:)
923
924	ENDWHERE
925
926	ENDIF
927
928	ENDIF
929
930	ENDDO ! loop over PFTs
931
932
933	CALL xios_orchidee_send_field("HERBIVORES",herbivores)
934	CALL xios_orchidee_send_field("LEAF_AGE",leaf_meanage)
935
936
937	! Write mean leaf age and time constant of probability of a leaf to be eaten by a herbivore
938	! to the stomate output file.
939	!CALL histwrite_p (hist_id_stomate, 'LEAFAGE', itime, &
940	! leaf_meanage, npts*nvm, horipft_index)
941	CALL histwrite_p (hist_id_stomate, 'HERBIVORES', itime, &
942	herbivores, npts*nvm, horipft_index)
943
944	IF (bavard.GE.4) WRITE(numout,*) 'Leaving turnover'
945
946	END SUBROUTINE turnover_diagnostic
947
948
949
950	!! ================================================================================================================================
951	!! SUBROUTINE : turnover_prognostic
952	!!
953	!>\BRIEF Calculate turnover of leaves, roots, fruits and sapwood due to aging or climatic
954	!! induced senescence. Calculate herbivory.
955	!!
956	!! DESCRIPTION : This subroutine determines the turnover of leaves and fine roots (and stems for grasses)
957	!! and simulates following processes:
958	!! 1. Mean leaf age is calculated from leaf ages of separate leaf age classes. Should actually
959	!! be recalculated at the end of the routine, but it does not change too fast. The mean leaf
960	!! age is calculated using the following equation:
961	!! \latexonly
962	!! \input{turnover_lma_update_eqn1.tex}
963	!! \endlatexonly
964	!! \n
965	!! 2. Meteorological senescence: the detection of the end of the growing season and shedding
966	!! of leaves, fruits and fine roots due to unfavourable meteorological conditions.
967	!! The model distinguishes three different types of "climatic" leaf senescence, that do not
968	!! change the age structure: sensitivity to cold temperatures, to lack of water, or both.
969	!! If meteorological conditions are fulfilled, a flag ::senescence is set to TRUE. Note
970	!! that evergreen species do not experience climatic senescence.
971	!! Climatic senescence is triggered by sensitivity to cold temperatures where the critical
972	!! temperature for senescence is calculated using the following equation:
973	!! \latexonly
974	!! \input{turnover_temp_crit_eqn2.tex}
975	!! \endlatexonly
976	!! \n
977	!! Climatic senescence is triggered by sensitivity to lack of water availability where the
978	!! moisture availability critical level is calculated using the following equation:
979	!! \latexonly
980	!! \input{turnover_moist_crit_eqn3.tex}
981	!! \endlatexonly
982	!! \n
983	!! Climatic senescence is triggered by sensitivity to temperature or to lack of water where
984	!! critical temperature and moisture availability are calculated as above.\n
985	!! Trees in climatic senescence lose their fine roots at the same rate as they lose their leaves.
986	!! The rate of biomass loss of both fine roots and leaves is presribed through the equation:
987	!! \latexonly
988	!! \input{turnover_clim_senes_biomass_eqn4.tex}
989	!! \endlatexonly
990	!! \n
991	!! with ::leaffall(j) a PFT-dependent time constant which is given in
992	!! ::stomate_constants. In grasses, leaf senescence is extended to the whole plant
993	!! (all carbon pools) except to its carbohydrate reserve.
994	!! 3. Senescence due to aging: the loss of leaves, fruits and biomass due to aging
995	!! At a certain age, leaves fall off, even if the climate would allow a green plant
996	!! all year round. Even if the meteorological conditions are favorable for leaf maintenance,
997	!! plants, and in particular, evergreen trees, have to renew their leaves simply because the
998	!! old leaves become inefficient. Roots, fruits (and stems for grasses) follow leaves.
999	!! The ??senescence?? rate varies with leaf age. Note that plant is not declared senescent
1000	!! in this case (wchich is important for allocation: if the plant loses leaves because of
1001	!! their age, it can renew them). The leaf turnover rate due to aging of leaves is calculated
1002	!! using the following equation:
1003	!! \latexonly
1004	!! \input{turnover_age_senes_biomass_eqn5.tex}
1005	!! \endlatexonly
1006	!! \n
1007	!! Drop all leaves if there is a very low leaf mass during senescence. After this, the biomass
1008	!! of different carbon pools both for trees and grasses is set to zero and the mean leaf age
1009	!! is reset to zero. Finally, the leaf fraction and leaf age of the different leaf age classes
1010	!! is set to zero. For deciduous trees: next to leaves, also fruits and fine roots are dropped.
1011	!! For grasses: all aboveground carbon pools, except the carbohydrate reserves are affected:
1012	!! 4. Update the leaf biomass, leaf age class fraction and the LAI
1013	!! Older leaves will fall more frequently than younger leaves and therefore the leaf age
1014	!! distribution needs to be recalculated after turnover. The fraction of biomass in each
1015	!! leaf class is updated using the following equation:
1016	!! \latexonly
1017	!! \input{turnover_update_LeafAgeDistribution_eqn6.tex}
1018	!! \endlatexonly
1019	!! \n
1020	!! 5. Simulate herbivory activity and update leaf and fruits biomass. Herbivore activity
1021	!! affects the biomass of leaves and fruits as well as stalks (only for grasses).
1022	!! However, herbivores do not modify leaf age structure.
1023	!! 6. Calculates fruit turnover for trees. Trees simply lose their fruits with a time
1024	!! constant ::tau_fruit(j), that is set to 90 days for all PFTs in ::stomate_constants
1025	!! 7. Convert sapwood to heartwood for trees and update heart and softwood above and
1026	!! belowground biomass. Sapwood biomass is converted into heartwood biomass
1027	!! with a time constant tau ::tau_sap(j) of 1 year. Note that this biomass conversion
1028	!! is not added to "turnover" as the biomass is not lost. For the updated heartwood,
1029	!! the sum of new heartwood above and new heartwood below after converting sapwood to
1030	!! heartwood, is saved as ::hw_new(:). Creation of new heartwood decreases the age of
1031	!! the plant ??carbon?? with a factor that is determined by: old heartwood ::hw_old(:)
1032	!! divided by the new heartwood ::hw_new(:)
1033	!!
1034	!! RECENT CHANGE(S) : None
1035	!!
1036	!! MAIN OUTPUT VARIABLES: ::Biomass of leaves, fruits, fine roots and sapwood above (latter for grasses only),
1037	!! ::Update LAI, ::Update leaf age distribution with new leaf age class fraction
1038	!!
1039	!! REFERENCE(S) :
1040	!! - Krinner, G., N. Viovy, N. de Noblet-Ducoudre, J. Ogee, J. Polcher, P.
1041	!! Friedlingstein, P. Ciais, S. Sitch and I.C. Prentice (2005), A dynamic global
1042	!! vegetation model for studies of the coupled atmosphere-biosphere system, Global
1043	!! Biogeochemical Cycles, 19, doi:10.1029/2003GB002199.
1044	!! - McNaughton, S. J., M. Oesterheld, D. A. Frank and K. J. Williams (1989),
1045	!! Ecosystem-level patterns of primary productivity and herbivory in terrestrial habitats,
1046	!! Nature, 341, 142-144, 1989.
1047	!! - Sitch, S., C. Huntingford, N. Gedney, P. E. Levy, M. Lomas, S. L. Piao, , Betts, R., Ciais, P., Cox, P.,
1048	!! Friedlingstein, P., Jones, C. D., Prentice, I. C. and F. I. Woodward : Evaluation of the terrestrial carbon
1049	!! cycle, future plant geography and climate-carbon cycle feedbacks using 5 dynamic global vegetation
1050	!! models (dgvms), Global Change Biology, 14(9), 2015â2039, 2008.
1051	!!
1052	!! FLOWCHART :
1053	!! \latexonly
1054	!! \includegraphics[scale=0.5]{turnover_flowchart_1.png}
1055	!! \includegraphics[scale=0.5]{turnover_flowchart_2.png}
1056	!! \endlatexonly
1057	!! \n
1058	!_ ================================================================================================================================
1059
1060	SUBROUTINE turnover_prognostic (npts, dt, PFTpresent, herbivores, &
1061	maxmoiavail_lastyear, minmoiavail_lastyear, moiavail_week, moiavail_month, &
1062	tlong_ref, t2m_longterm, t2m_month, t2m_week, veget_max, &
1063	gdd_from_growthinit, leaf_age, leaf_frac, age, &
1064	biomass, turnover, senescence, turnover_time, &
1065	circ_class_biomass, circ_class_n, ind, allow_initpheno, &
1066	when_growthinit, tau_eff_leaf, tau_eff_sap, tau_eff_root, &
1067	harvest_pool, harvest_type, harvest_cut, harvest_area, &
1068	resolution, wstress_month)
1069
1070	!! 0. Variable and parameter declaration
1071
1072	!! 0.1 Input variables
1073
1074	INTEGER(i_std), INTENT(in) :: npts !! Domain size - number of grid cells
1075	!! (unitless)
1076	REAL(r_std), INTENT(in) :: dt !! time step (dt_days)
1077	LOGICAL, DIMENSION(:,:), INTENT(in) :: PFTpresent !! PFT exists (true/false)
1078	REAL(r_std),DIMENSION(:,:), INTENT(in) :: resolution !! The length in m of each grid-box in X and Y
1079	REAL(r_std), DIMENSION(:,:), INTENT(in) :: herbivores !! time constant of probability of a leaf to
1080	!! be eaten by a herbivore (days)
1081	REAL(r_std), DIMENSION(:,:), INTENT(in) :: maxmoiavail_lastyear !! last year's maximum moisture availability
1082	!! (0-1, unitless)
1083	REAL(r_std), DIMENSION(:,:), INTENT(in) :: minmoiavail_lastyear !! last year's minimum moisture availability
1084	!! (0-1, unitless)
1085	REAL(r_std), DIMENSION(:,:), INTENT(in) :: moiavail_week !! "weekly" moisture availability
1086	!! (0-1, unitless)
1087	REAL(r_std), DIMENSION(:,:), INTENT(in) :: moiavail_month !! "monthly" moisture availability
1088	!! (0-1, unitless)
1089	REAL(r_std), DIMENSION(:), INTENT(in) :: tlong_ref !! "long term" 2 meter reference
1090	!! temperatures (K)
1091	REAL(r_std), DIMENSION(:), INTENT(in) :: t2m_longterm !! "longterm" 2-meter temperatures (K)
1092	REAL(r_std), DIMENSION(:), INTENT(in) :: t2m_month !! "monthly" 2-meter temperatures (K)
1093	REAL(r_std), DIMENSION(:), INTENT(in) :: t2m_week !! "weekly" 2 meter temperatures (K)
1094	REAL(r_std), DIMENSION(:,:), INTENT(in) :: veget_max !! "maximal" coverage fraction of a PFT (LAI
1095	!! -> infinity) on ground (unitless)
1096	REAL(r_std), DIMENSION(:,:), INTENT(in) :: gdd_from_growthinit !! gdd senescence for crop
1097	REAL(r_std), DIMENSION(:,:), INTENT(in) :: when_growthinit !! How many days ago was the beginning of
1098	!! the growing season (days)
1099	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: circ_class_n !! Number of individuals in each circ class
1100	!! @tex $(m^{-2})$ @endtex
1101	LOGICAL, DIMENSION(:,:), INTENT(in) :: allow_initpheno !! Are we allowed to declare the
1102	!! beginning of the growing season?
1103	!! and the end of senescence (true/false)
1104	REAL(r_std), DIMENSION(:,:), INTENT(in) :: tau_eff_root !! Effective root turnover time that accounts
1105	!! waterstress (days)
1106	REAL(r_std), DIMENSION(:,:), INTENT(in) :: tau_eff_sap !! Effective sapwood turnover time that accounts
1107	!! waterstress (days)
1108	REAL(r_std), DIMENSION(:,:), INTENT(in) :: tau_eff_leaf !! Effective leaf turnover time that accounts
1109	!! waterstress (days)
1110	REAL(r_std), DIMENSION(:,:), INTENT(in) :: wstress_month !! Water stress factor, based on hum_rel_daily
1111	!! (unitless, 0-1)
1112
1113
1114	!! 0.2 Output variables
1115
1116
1117	!! 0.3 Modified variables
1118
1119	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: leaf_age !! age of the leaves (days)
1120	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: leaf_frac !! fraction of leaves in leaf age class
1121	!! (0-1, unitless)
1122	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: harvest_pool !! The wood and biomass that have been
1123	!! havested by humans @tex $(gC)$ @endtex
1124	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: harvest_type !! Type of management that resulted
1125	!! in the harvest (unitless)
1126	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: harvest_cut !! Type of cutting that was used for the harvest
1127	!! (unitless)
1128	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: harvest_area !! Harvested area (m^{2})
1129	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: age !! age (years)
1130	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ind !! Density of individuals
1131	!! @tex $(m^{-2})$ @endtex
1132	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: turnover_time !! turnover_time of grasses (days)
1133	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: biomass !! biomass @tex ($gC m^{-2}$) @endtex
1134	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: turnover !! Turnover @tex ($gC m^{-2}$) @endtex
1135	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(inout) :: circ_class_biomass !! Biomass of the componets of the model
1136	!! tree within a circumference
1137	!! class @tex $(gC ind^{-1})$ @endtex
1138	LOGICAL, DIMENSION(:,:), INTENT(inout) :: senescence !! is the plant senescent? (true/false)
1139	!! (interesting only for deciduous trees:
1140	!! carbohydrate reserve)
1141
1142	!! 0.4 Local variables
1143
1144	LOGICAL, DIMENSION(npts) :: shed_rest !! shed the remaining leaves? (true/false)
1145	INTEGER(i_std) :: ivm,iele,ilage,ipts !! Index (unitless)
1146	INTEGER(i_std) :: ipar,icirc,imbc !! Index (unitless)
1147	INTEGER(i_std), SAVE :: turn_count !!
1148	!$OMP THREADPRIVATE(turn_count)
1149	REAL(r_std), DIMENSION(npts,nvm) :: leaf_meanage !! mean age of the leaves (days)
1150	REAL(r_std), DIMENSION(npts,nvm) :: lai !! leaf area index @tex ($m^2 m^{-2}$)
1151	!! @endtex
1152	REAL(r_std), DIMENSION(npts) :: dturnover !! Intermediate variable for turnover ??
1153	!! @tex ($gC m^{-2}$) @endtex
1154	REAL(r_std), DIMENSION(npts) :: moiavail_crit !! critical moisture availability, function
1155	!! of last year's moisture availability
1156	!! (0-1, unitless)
1157	REAL(r_std), DIMENSION(npts) :: tl !! long term annual mean temperature, (C)
1158	REAL(r_std), DIMENSION(npts) :: t_crit !! critical senescence temperature, function
1159	!! of long term annual temperature (K)
1160	REAL(r_std), DIMENSION(npts) :: sapconv !! Sapwood conversion @tex ($gC m^{-2}$)
1161	!! @endtex
1162	REAL(r_std), DIMENSION(npts) :: hw_old !! old heartwood mass @tex ($gC m^{-2}$)
1163	!! @endtex
1164	REAL(r_std), DIMENSION(npts) :: hw_new !! new heartwood mass @tex ($gC m^{-2}$)
1165	!! @endtex
1166	REAL(r_std), DIMENSION(npts) :: lm_old !! old leaf mass @tex ($gC m^{-2}$) @endtex
1167	REAL(r_std), DIMENSION(npts,nleafages) :: delta_lm !! leaf mass change for each age class @tex
1168	!! ($gC m^{-2}$) @endtex
1169	REAL(r_std), DIMENSION(npts) :: turnover_rate !! turnover rate (unitless)
1170	REAL(r_std), DIMENSION(npts,nvm) :: leaf_age_crit !! critical leaf age (days)
1171	REAL(r_std), DIMENSION(npts,nvm) :: new_turnover_time !! instantaneous turnover time (days)
1172	REAL(r_std), DIMENSION(npts,nvm) :: harvest_time !! Prescribed harvest time adjusted for the
1173	!! longterm temperature at the pixel (days)
1174	REAL(r_std) :: branch_turn !! turnover of branches (how much goes
1175	!! to litter each day) (cf article CO2fix)
1176	REAL(r_std), SAVE :: ss_branch_turn !! Variations for sensitivity analysis
1177	!$OMP THREADPRIVATE(ss_branch_turn)
1178	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements):: bm_old !! old root mass
1179	REAL(r_std), DIMENSION(npts,nvm) :: sapwood_age !! sapwood age (days)
1180	REAL(r_std), DIMENSION(npts) :: sw_old !! old sapwood mass
1181	!! (gC/(m**2 of nat/agri ground))
1182	REAL(r_std), DIMENSION(npts) :: sw_new !! new sapwood mass
1183	!! (gC/(m**2 of nat/agri ground)
1184	REAL(r_std), DIMENSION(npts) :: pixel_area !! surface area of the pixel @tex ($m^{2}$) @endtex
1185	REAL(r_std), DIMENSION(npts,nvm,nmbcomp,nelements)&
1186	:: check_intern !! Contains the components of the internal
1187	!! mass balance chech for this routine
1188	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
1189	REAL(r_std), DIMENSION(npts,nvm,nelements) :: closure_intern !! Check closure of internal mass balance
1190	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
1191	REAL(r_std), DIMENSION(npts,nvm,nelements) :: pool_start, pool_end !! Start and end pool of this routine
1192	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
1193
1194	!_ ================================================================================================================================
1195
1196	IF (bavard.GE.3) WRITE(numout,*) 'Entering prognostic turnover'
1197
1198	!! 1. first call - output messages
1199
1200	IF ( firstcall ) THEN
1201
1202	WRITE(numout,*) 'turnover:'
1203
1204	WRITE(numout,*) ' > minimum mean leaf age for senescence (days) (::min_leaf_age_for_senescence) : '&
1205	,min_leaf_age_for_senescence
1206
1207	turn_count = 0
1208
1209	!Config Key = SS_BRANCH_TURN
1210	!Config Desc =
1211	!Config If =
1212	!Config Def =
1213	!Config Help =
1214	!Config Units =
1215
1216	! this needs to be initialized, else it crashes in the getin_p call
1217	ss_branch_turn=un
1218	CALL getin_p ('SS_BRANCH_TURN',ss_branch_turn)
1219
1220	firstcall = .FALSE.
1221
1222	ENDIF
1223
1224
1225	!! 2. Initializations
1226
1227	!! 2.1 set output to zero
1228	turnover(:,:,:,:) = zero
1229	new_turnover_time(:,:) = zero
1230	harvest_time(:,:) = zero
1231
1232	!! 2.2 Initialize check for mass balance closure
1233	! The mass balance is calculated at the end of this routine
1234	! in section 10.
1235	! turnover is always equal to zero, so maybe that term doesn't need
1236	! to be included here.
1237	pool_start = zero
1238	DO ipar = 1,nparts
1239	DO iele = 1,nelements
1240	! Initial biomass pool
1241	pool_start(:,:,iele) = pool_start(:,:,iele) + &
1242	(biomass(:,:,ipar,iele) * veget_max(:,:)) + &
1243	(turnover(:,:,ipar,iele) * veget_max(:,:))
1244	ENDDO
1245	ENDDO
1246
1247	! Surface area (m2) of the pixel
1248	pixel_area(:) = resolution(:,1) * resolution(:,2)
1249
1250	! The biomass harvest pool is expressed in gC pixel-1 So, it
1251	! shouldn't be multiplied by veget_max but it should be divided
1252	! by pixel_area to obtain gC m-2.
1253	DO ivm = 1,nvm
1254	pool_start(:,ivm,icarbon) = pool_start(:,ivm,icarbon) + &
1255	SUM(harvest_pool(:,ivm,:,icarbon),2) / &
1256	pixel_area(:)
1257	ENDDO
1258
1259	! save old biomass
1260	bm_old(:,:,:,:) = biomass(:,:,:,:)
1261
1262	! Compute the turnover of branches
1263	! (2.5% per year : Orchidee standard, Masera 2003, Lehtonen 2004, DeAngelis 1981)
1264	branch_turn = 0.025/(one_year/dt)*ss_branch_turn
1265
1266	! Calculate lai
1267	DO ivm = 2,nvm
1268	DO ipts=1,npts
1269	lai(ipts,ivm) = biomass_to_lai(biomass(ipts,ivm,ileaf,icarbon),ivm)
1270	ENDDO
1271	ENDDO
1272	lai(:,1) = zero
1273
1274	!! 2.2 Recalculate mean leaf age
1275	! Mean leaf age is recalculated from leaf ages of separate leaf age classes.
1276	! Leaf age is used as a criterion in several of the following IF/WHERE loops
1277	! The mean leaf age is calculated using the following equation:
1278	! \latexonly
1279	! \input{turnover_lma_update_eqn1.tex}
1280	! \endlatexonly
1281	! \n
1282	leaf_meanage(:,:) = zero
1283
1284	DO ilage = 1, nleafages
1285
1286	leaf_meanage(:,:) = leaf_meanage(:,:) + leaf_age(:,:,ilage) * leaf_frac(:,:,ilage)
1287
1288	ENDDO
1289
1290	!! 3. Climatic senescence
1291
1292	! Three different types of "climatic" leaf senescence,
1293	! that do not change the age structure.
1294	DO ivm = 2,nvm ! Loop over # PFTs
1295
1296	!! 3.1 Determine if there is climatic senescence.
1297	! The climatic senescence can be of three types: sensitivity to cold temperatures,
1298	! to lack of water, or both. If meteorological conditions are fulfilled, a flag
1299	! senescence is set to TRUE. Evergreen species do not experience climatic senescence.
1300	SELECT CASE ( senescence_type(ivm) )
1301
1302	!! 3.1.1 Crops
1303	! Crop senescence is based on a GDD criterium as in crop models.
1304	CASE ('crop' )
1305
1306	!+++CHECK+++
1307	! This is pure plumbing under time pressure (says the one who introduced it).
1308	! Switch senescence to true and it will stay true until begin_leaves becomes true
1309	! A growing season should last about 4 months where the longterm temperature
1310	! is 282 K (crude observation from N. France) before senescence can set in unless
1311	! the growing degrees are reached. In colder and warmer places the growing season
1312	! is shortened or lengthened. Because in warmer places the growing season lasts
1313	! longer, crop can be planted that required more light.
1314	harvest_time(:,ivm) = 130 + 3 * ( t2m_longterm(:) - 282)
1315
1316	WHERE ( ( biomass(:,ivm,ileaf,icarbon) .GT. zero ) .AND. &
1317	( leaf_meanage(:,ivm) .GT. min_leaf_age_for_senescence(ivm) ) .AND. &
1318	( ( gdd_from_growthinit(:,ivm) .GT. gdd_senescence(ivm) ) .OR. &
1319	( ( when_growthinit(:,ivm) .GT. harvest_time(:,ivm) ) ) ) )
1320
1321	senescence(:,ivm) = .TRUE.
1322
1323	ENDWHERE
1324	!+++++++++++
1325
1326
1327	!! 3.1.2 Summergreen species
1328	! Climatic senescence is triggered by sensitivity to cold temperatures as follows:
1329	! If biomass is large enough (i.e. when it is greater than zero),
1330	! AND (i.e. when leaf mean age is above a certain PFT-dependent treshold ::min_leaf_age_for_senescence,
1331	! which is given in constants),
1332	! AND the monthly temperature is low enough (i.e. when monthly temperature ::t2m_month(:) is below a critical
1333	! temperature ::t_crit(:), which is calculated in this module),
1334	! AND the temperature tendency is negative (i.e. when weekly temperatures ::t2m_week(:) are lower than monthly
1335	! temperatures ::t2m_month(:))
1336	! If these conditions are met, senescence is set to TRUE.
1337	!
1338	! The critical temperature for senescence is calculated using the following equation:
1339	! \latexonly
1340	! \input{turnover_temp_crit_eqn2.tex}
1341	! \endlatexonly
1342	! \n
1343	! Critical temperature for senescence may depend on long term annual mean temperature
1344	CASE ( 'cold' )
1345
1346	tl(:) = tlong_ref(:) - ZeroCelsius
1347	t_crit(:) = ZeroCelsius + senescence_temp(ivm,1) + &
1348	tl(:) * senescence_temp(ivm,2) + &
1349	tl(:)tl(:) senescence_temp(ivm,3)
1350
1351	WHERE ( ( biomass(:,ivm,ileaf,icarbon) .GT. zero ) .AND. &
1352	( leaf_meanage(:,ivm) .GT. min_leaf_age_for_senescence(ivm) ) .AND. &
1353	( t2m_month(:) .LT. t_crit(:) ) .AND. &
1354	( t2m_week(:) .LT. t2m_month(:) ) )
1355
1356	senescence(:,ivm) = .TRUE.
1357
1358	ENDWHERE
1359
1360	!! 3.1.3 Raingreen species
1361	! Climatic senescence is triggered by sensitivity to lack of water availability as follows:
1362	! If biomass is large enough (i.e. when it is greater than zero),
1363	! AND (i.e. when leaf mean age is above a certain PFT-dependent treshold ::min_leaf_age_for_senescence,
1364	! which is given in ::stomate_constants),
1365	! AND the moisture availability drops below a critical level (i.e. when weekly moisture availability
1366	! ::moiavail_week(:,ivm) is below a critical moisture availability ::moiavail_crit(:),
1367	! which is calculated in this module),
1368	! If these conditions are met, senescence is set to TRUE.
1369	!
1370	! The moisture availability critical level is calculated using the following equation:
1371	! \latexonly
1372	! \input{turnover_moist_crit_eqn3.tex}
1373	! \endlatexonly
1374	! \n
1375	CASE ( 'dry' )
1376
1377	moiavail_crit(:) = &
1378	MIN( MAX( minmoiavail_lastyear(:,ivm) + hum_frac(ivm) * &
1379	( maxmoiavail_lastyear(:,ivm) - minmoiavail_lastyear(:,ivm) ), &
1380	senescence_hum(ivm) ), &
1381	nosenescence_hum(ivm) )
1382
1383	WHERE ( ( biomass(:,ivm,ileaf,icarbon) .GT. zero ) .AND. &
1384	( leaf_meanage(:,ivm) .GT. min_leaf_age_for_senescence(ivm) ) .AND. &
1385	( moiavail_week(:,ivm) .LT. moiavail_crit(:) ) )
1386
1387	senescence(:,ivm) = .TRUE.
1388
1389	ENDWHERE
1390
1391	!! 3.1.4 Mixed criterion: Climatic senescence is triggered by sensitivity to temperature or to lack of water
1392	! Climatic senescence is triggered by sensitivity to temperature or to lack of water availability as follows:
1393	! If biomass is large enough (i.e. when it is greater than zero),
1394	! AND (i.e. when leaf mean age is above a certain PFT-dependent treshold ::min_leaf_age_for_senescence,
1395	! which is given in ::stomate_constants),
1396	! AND the moisture availability drops below a critical level (i.e. when weekly moisture availability
1397	! ::moiavail_week(:,ivm) is below a critical moisture availability ::moiavail_crit(:), calculated in this module),
1398	! OR
1399	! the monthly temperature is low enough (i.e. when monthly temperature ::t2m_month(:) is below a critical
1400	! temperature ::t_crit(:), calculated in this module),
1401	! AND the temperature tendency is negative (i.e. when weekly temperatures ::t2m_week(:) are lower than
1402	! monthly temperatures ::t2m_month(:)).
1403	! If these conditions are met, senescence is set to TRUE.
1404	CASE ( 'mixed' )
1405
1406	moiavail_crit(:) = MIN( MAX( minmoiavail_lastyear(:,ivm) + hum_frac(ivm) * &
1407	(maxmoiavail_lastyear(:,ivm) - minmoiavail_lastyear(:,ivm) ), &
1408	senescence_hum(ivm) ), nosenescence_hum(ivm) )
1409	tl(:) = tlong_ref(:) - ZeroCelsius
1410	t_crit(:) = ZeroCelsius + senescence_temp(ivm,1) + &
1411	tl(:) * senescence_temp(ivm,2) + &
1412	tl(:)tl(:) senescence_temp(ivm,3)
1413
1414	IF ( is_tree(ivm) ) THEN
1415
1416	! critical temperature for senescence may depend on long term annual mean temperature
1417	WHERE ( ( biomass(:,ivm,ileaf,icarbon) .GT. zero ) .AND. &
1418	( leaf_meanage(:,ivm) .GT. min_leaf_age_for_senescence(ivm) ) .AND. &
1419	( ( moiavail_week(:,ivm) .LT. moiavail_crit(:) ) .OR. &
1420	( ( t2m_month(:) .LT. t_crit(:) ) .AND. ( t2m_week(:) .LT. t2m_month(:) ) ) ) )
1421
1422	senescence(:,ivm) = .TRUE.
1423
1424	ENDWHERE
1425
1426	! Grasses
1427	ELSE
1428
1429	!!$ !+++CHECK+++
1430	!!$ ! There is no longer this kind of turnover outside the senescence period.
1431	!!$ ! The physiological basis of this turnover is unclear because the turn
1432	!!$ ! over due to ageing has been accounted for in Â§4. Not sure whether the
1433	!!$ ! values for ::turnover_time calculated here are ever used because turnover
1434	!!$ ! from ageing in Â§4 recalculates the turnover times.
1435	!!$ WHERE ( ( biomass(:,ivm,ileaf,icarbon) .GT. zero ) .AND. &
1436	!!$ ( leaf_meanage(:,ivm) .GT. min_leaf_age_for_senescence(ivm) ) .AND. &
1437	!!$ ( ( moiavail_week(:,ivm) .LT. moiavail_crit(:) )))
1438	!!$
1439	!!$ turnover_time(:,ivm,ileaf) = tau_eff_leaf(:,ivm) * &
1440	!!$ (un - (un - (moiavail_week(:,ivm)/ moiavail_crit(:)))**deux)
1441	!!$ turnover_time(:,ivm,iroot) = tau_eff_root(:,ivm) * &
1442	!!$ (un - (un - (moiavail_week(:,ivm)/ moiavail_crit(:)))**deux)
1443	!!$ turnover_time(:,ivm,isapabove) = tau_eff_sap(:,ivm) * &
1444	!!$ (un - (un - (moiavail_week(:,ivm)/ moiavail_crit(:)))**deux)
1445	!!$ turnover_time(:,ivm,ifruit) = tau_fruit(ivm) * &
1446	!!$ (un - (un - (moiavail_week(:,ivm)/ moiavail_crit(:)))**deux)
1447	!!$
1448	!!$ ENDWHERE
1449	!!$ !+++++++++++
1450
1451	! Because the DOFOCO increase of LAI is much slower than in the trunk there
1452	! is a period at the start of the growing season where LAI is less than lai_happy
1453	! In that period t2_month is still increasing from its winter values but could be
1454	! below t_crit. That condition resulted in many grasslands along the Atlantic and
1455	! North Sea coast to go into senescence in May. If someone wants to put it back in
1456	! be carefull with the brackets.((t2m_month(:) .LT. t_crit(:)) .AND. (lai(:,ivm)
1457	! .GT. lai_happy(ivm)))
1458	WHERE ( ( (biomass(:,ivm,ileaf,icarbon) .GT. zero) .AND. &
1459	(leaf_meanage(:,ivm) .GT. min_leaf_age_for_senescence(ivm)) .AND. &
1460	(t2m_month(:) .LT. ZeroCelsius) .AND. (t2m_week(:) .LT. t2m_month(:)) ) )
1461
1462	! Shed leaves, roots and fruits at a high rate = senescence
1463	turnover_time(:,ivm,ileaf) = leaffall(ivm)
1464	turnover_time(:,ivm,iroot) = leaffall(ivm)
1465	turnover_time(:,ivm,ifruit) = leaffall(ivm)
1466
1467	! Slowly turnover the structural carbon this allows to store
1468	! reserves. This structural pool is essential for the functioning
1469	! of the labile and reserve pools
1470	turnover_time(:,ivm,isapabove) = tau_eff_sap(:,ivm)
1471	senescence(:,ivm) = .TRUE.
1472
1473	ENDWHERE
1474
1475	ENDIF
1476
1477	!! 3.1.5 Evergreen species
1478	! Evergreen species do not experience climatic senescence
1479	CASE ( 'none' )
1480
1481	!! 3.1.6 Other cases
1482	! In case no climatic senescence type is recognized.
1483	CASE default
1484
1485	WRITE(numout,*) ' turnover: don''t know how to treat this PFT.'
1486	WRITE(numout,*) ' number (::ivm) : ',ivm
1487	WRITE(numout,*) ' senescence type (::senescence_type(ivm)) : ',senescence_type(ivm)
1488
1489	STOP
1490
1491	END SELECT
1492
1493	!---TEMP---
1494	IF (ivm .EQ. test_pft .AND. ld_pheno) THEN
1495	WRITE(numout,*) 'phenology, senescence, ',senescence(:,ivm)
1496	ENDIF
1497	!----------
1498
1499	!! 3.2 Drop leaves and roots, plus stems and fruits for grasses
1500	IF ( is_tree(ivm) ) THEN
1501
1502	!! 3.2.1 Trees in climatic senescence
1503	! Trees in climate senescence lose their fine roots at the same rate as they lose their leaves.
1504	! The rate of biomass loss of both fine roots and leaves is presribed through the equation:
1505	! \latexonly
1506	! \input{turnover_clim_senes_biomass_eqn4.tex}
1507	! \endlatexonly
1508	! \n
1509	! with ::leaffall(ivm) a PFT-dependent time constant which is given in ::stomate_constants).
1510
1511	! Notice that we do not change the carbohydrate (icarbres) and labile (ilabile) pools here. We
1512	! assume that trees can anticipate senescence and move these pools
1513	! to more permanent storage areas so that are not lost, unlike
1514	! in the case of harvesting.
1515	WHERE ( senescence(:,ivm) )
1516
1517	! Turnover at the stand level (gC m-2)
1518	turnover(:,ivm,ileaf,icarbon) = biomass(:,ivm,ileaf,icarbon) * dt / leaffall(ivm)
1519	turnover(:,ivm,iroot,icarbon) = biomass(:,ivm,iroot,icarbon) * dt / leaffall(ivm)
1520
1521	! Update biomass at the stand level (gC m-2)
1522	biomass(:,ivm,ileaf,icarbon) = biomass(:,ivm,ileaf,icarbon) - turnover(:,ivm,ileaf,icarbon)
1523	biomass(:,ivm,iroot,icarbon) = biomass(:,ivm,iroot,icarbon) - turnover(:,ivm,iroot,icarbon)
1524
1525	ENDWHERE
1526
1527	! Turnover at tree level (gC tree-1)
1528	DO icirc = 1,ncirc
1529
1530	WHERE ( senescence(:,ivm) )
1531
1532	circ_class_biomass(:,ivm,icirc,ileaf,icarbon) = circ_class_biomass(:,ivm,icirc,ileaf,icarbon) * &
1533	(un - dt / leaffall(ivm))
1534	circ_class_biomass(:,ivm,icirc,iroot,icarbon) = circ_class_biomass(:,ivm,icirc,iroot,icarbon) * &
1535	(un - dt / leaffall(ivm))
1536
1537	ENDWHERE
1538
1539	ENDDO
1540
1541	! Grasses and crops
1542	ELSE
1543
1544	!! 3.2.2.1 crops with 'crop' phenological model
1545	IF (senescence_type(ivm) .EQ. 'crop') THEN
1546
1547	DO ipts = 1,npts
1548
1549	IF (senescence(ipts,ivm)) THEN
1550
1551	! Crops are planted every year. So make sure to remove everything
1552	! at the end of senescence which for crops is actually harvest.
1553	! Next stomate_prescribe should prescribe a new crop. If sapwood
1554	! is not removed, some is left on site and the model will try to
1555	! grow a crop. However, there are not enough reserves left so the
1556	! crop will die. This results in one year of normal crop growth
1557	! (following the prescribe) and then one year of almost no growth
1558	! (the year starting from the too low reserves). This issue has
1559	! been fixed.
1560	CALL crop_harvest(npts, ipts, ivm, dt, &
1561	veget_max, turnover, harvest_pool, harvest_type, &
1562	harvest_cut, harvest_area, resolution, biomass, &
1563	ind, leaf_meanage)
1564
1565	ENDIF
1566
1567	ENDDO
1568
1569	!! 3.2.2.2 grass or crops based on 'mixed' phenological model
1570	! Phenological turnover
1571	ELSEIF (senescence_type(ivm) .EQ. 'mixed') THEN
1572
1573	WHERE (senescence(:,ivm))
1574
1575	turnover(:,ivm,ileaf,icarbon) = biomass(:,ivm,ileaf,icarbon) * dt / turnover_time(:,ivm,ileaf)
1576	turnover(:,ivm,isapabove,icarbon) = biomass(:,ivm,isapabove,icarbon) * dt / turnover_time(:,ivm,isapabove)
1577	turnover(:,ivm,iroot,icarbon) = biomass(:,ivm,iroot,icarbon) * dt / turnover_time(:,ivm,iroot)
1578	turnover(:,ivm,ifruit,icarbon) = biomass(:,ivm,ifruit,icarbon) * dt / turnover_time(:,ivm,ifruit)
1579
1580	ENDWHERE
1581
1582	! Update the biomass pools for grasses and crops
1583	biomass(:,ivm,ileaf,icarbon) = biomass(:,ivm,ileaf,icarbon) - turnover(:,ivm,ileaf,icarbon)
1584	biomass(:,ivm,isapabove,icarbon) = biomass(:,ivm,isapabove,icarbon) - turnover(:,ivm,isapabove,icarbon)
1585	biomass(:,ivm,iroot,icarbon) = biomass(:,ivm,iroot,icarbon) - turnover(:,ivm,iroot,icarbon)
1586	biomass(:,ivm,ifruit,icarbon) = biomass(:,ivm,ifruit,icarbon) - turnover(:,ivm,ifruit,icarbon)
1587
1588	!! 3.2.2.3 Grasses modeled as an evergreen biome
1589	ELSEIF (senescence_type(ivm) .EQ. 'none') THEN
1590
1591	turnover(:,ivm,:,icarbon) = zero
1592
1593	!! 3.2.2.4 Conceptual problem
1594	ELSE
1595
1596	WRITE(numout,*) 'ERROR: senenescence type not known'
1597
1598	END IF
1599
1600	! Synchronize circ_class_biomass and biomass
1601	DO ipar = 1,nparts
1602
1603	WHERE (ind(:,ivm) .GT. zero)
1604
1605	circ_class_biomass(:,ivm,1,ipar,icarbon) = &
1606	biomass(:,ivm,ipar,icarbon) / ind(:,ivm)
1607
1608	ELSEWHERE
1609
1610	circ_class_biomass(:,ivm,1,ipar,icarbon) = zero
1611
1612	ENDWHERE
1613
1614	ENDDO ! nparts
1615
1616	ENDIF ! tree/grass
1617
1618	!---TEMP---
1619	IF (ivm .EQ. test_pft .AND. ld_trnov) THEN
1620	WRITE(numout,*) 'Check turnover'
1621	! DO ipar = 1,nparts
1622	! WRITE(numout,*) 'biomass vs circ_class_biomass (gC m-2), ', biomass(1,ivm,ipar,icarbon), &
1623	! SUM(circ_class_biomass(1,ivm,:,ipar,icarbon)*circ_class_n(1,ivm,:))
1624	! ENDDO
1625	WRITE(numout,*) 'TOTAL biomass vs circ_class_biomass, ', SUM(biomass(1,ivm,:,icarbon)), &
1626	SUM( SUM(circ_class_biomass(1,ivm,:,:,icarbon),2) * circ_class_n(1,ivm,:),1)
1627	ENDIF
1628	!----------
1629
1630	ENDDO ! loop over PFTs
1631
1632
1633	!! 4. Leaf fall
1634
1635	! At a certain age, leaves fall off, even if the climate would allow a
1636	! green plant all year round. Even if the meteorological conditions are
1637	! favorable for leaf maintenance, plants, and in particular, evergreen
1638	! trees, have to renew their leaves simply because the old leaves become
1639	! inefficient. Another reason for leaf fall may be waterstress. When
1640	! the plant experiences water stress some of the leaves will die.
1641	! Wstress is calculated from the ratio between stressed and unstressed
1642	! GPP. If the wstress is less than 1, this implies that we have to maintain
1643	! a complete canopy (and therefore have the Ra cost) but that only part
1644	! of the canopy will contribute to GPP. Although this is exactely what
1645	! happens on a warm summer day, what we are trying to account for here is
1646	! the impact of a lasting drought (days to weeks). Adaptation to the
1647	! growing conditions are accounted for through KF (see allocation).
1648	!
1649	! Roots, fruits (and stems) follow leaves. The decay rate varies with
1650	! leaf age. Note that plant is not declared senescent in this case
1651	! (wchich is important for allocation: if the plant loses leaves because
1652	! of their age, it can renew them).
1653	!
1654	! Notice that we do not change the reserve pools here (ilabile and
1655	! icarbres). We assume that the tree will move these pools out of old
1656	! leaves and therefore the carbon will not be lost. Trees are smart.
1657	!
1658	! The leaf turnover rate due to aging of leaves is calculated using
1659	! the following equation:
1660	! \latexonly
1661	! \input{turnover_age_senes_biomass_eqn5.tex}
1662	! \endlatexonly
1663	! \n
1664	DO ivm = 2,nvm ! Loop over # PFTs
1665
1666	!! 4.1 Calculate critical leaf age
1667	! save old leaf mass
1668	lm_old(:) = biomass(:,ivm,ileaf,icarbon)
1669
1670	!! initialize leaf mass change in age class
1671	delta_lm(:,:) = zero
1672
1673	! Trees and crops
1674	! leaves, roots, fruits roots and fruits follow leaves.
1675	IF ( is_tree(ivm) .OR. (.NOT. natural(ivm)) ) THEN
1676
1677	! Critical age: prescribed for trees
1678	leaf_age_crit(:,ivm) = tau_eff_leaf(:,ivm)
1679
1680	! Grasses
1681	! leaves, roots, fruits, sap follow leaves.
1682	ELSE
1683
1684	! Critical leaf age depends on long-term temperature
1685	! generally, lower turnover in cooler climates.
1686	leaf_age_crit(:,ivm) = MIN( tau_eff_leaf(:,ivm) * leaf_age_crit_coeff(1) , &
1687	MAX( tau_eff_leaf(:,ivm) * leaf_age_crit_coeff(2) , &
1688	tau_eff_leaf(:,ivm) - leaf_age_crit_coeff(3) * &
1689	( tlong_ref(:)-ZeroCelsius - leaf_age_crit_tref ) ) )
1690
1691	END IF
1692
1693	!! 4.2 Turnover for different leaf age classes
1694	DO ilage = 1, nleafages
1695
1696	turnover_rate(:) = zero
1697
1698	WHERE ( leaf_age(:,ivm,ilage) .GT. leaf_age_crit(:,ivm) / deux)
1699
1700	turnover_rate(:) = MIN( 0.99_r_std, dt / ( leaf_age_crit(:,ivm) * &
1701	( leaf_age_crit(:,ivm) / leaf_age(:,ivm,ilage) )**quatre ) )
1702
1703	ENDWHERE
1704
1705	IF (is_tree(ivm)) THEN
1706
1707	! Stand level turnover (gC m-2)
1708	! Leaves
1709	! Water stress is increasing (note that ::wstress is thus
1710	! decreasing)so the LAI will be lowered if the decrease is
1711	! small, this will result in a cosmetic change in LAI.
1712	! If water stress is increasing over a the course of a week,
1713	! this decrease may become important. The correction is
1714	! cummulative, for example, ::biomass(ileaf) is decreased
1715	! by 10% in day one and the remaining fraction of
1716	! ::biomass(ileaf) can be decreased by another 10% the next
1717	! day, etc ... The choice to multiply with ::wstress_month
1718	! is arbitrary, we could have well used ::wstress_week. By
1719	! using the monthly value we hope a smoother decrease
1720	! compared to using ::wstress_week or ::wstress_day.
1721
1722	!+++CHECK+++
1723	! Switch off the drought feedback
1724	! Do we have enough biomass to satisfy the turnover?
1725	!!$ WHERE ( (wstress_month(:,ivm)**(1/tau_hum_month)) .LT. un )
1726	!!$
1727	!!$ ! Account for water stress
1728	!!$ turnover_rate(:) = turnover_rate(:) / &
1729	!!$ (wstress_month(:,ivm)**(1/tau_hum_month))
1730	!!$
1731	!!$ ENDWHERE
1732	!+++++++++++
1733
1734	! Update biomass and turnover pools
1735	dturnover(:) = biomass(:,ivm,ileaf,icarbon) * leaf_frac(:,ivm,ilage) * turnover_rate(:)
1736	biomass(:,ivm,ileaf,icarbon) = biomass(:,ivm,ileaf,icarbon) - dturnover(:)
1737	turnover(:,ivm,ileaf,icarbon) = turnover(:,ivm,ileaf,icarbon) + dturnover(:)
1738
1739	! save leaf mass change
1740	delta_lm(:,ilage) = - dturnover(:)
1741
1742	!+++CHECK+++
1743	! It is not clear why roots follow the leaves. A turnover rate for roots based on
1744	! tau_eff_root can be easily calculated
1745	! Roots
1746	dturnover(:) = biomass(:,ivm,iroot,icarbon) * leaf_frac(:,ivm,ilage) * turnover_rate(:)
1747	biomass(:,ivm,iroot,icarbon) = biomass(:,ivm,iroot,icarbon) - dturnover(:)
1748	turnover(:,ivm,iroot,icarbon) = turnover(:,ivm,iroot,icarbon) + dturnover(:)
1749	!+++++++++++
1750
1751	! Fruit
1752	dturnover(:) = biomass(:,ivm,ifruit,icarbon) * leaf_frac(:,ivm,ilage) * turnover_rate(:)
1753	biomass(:,ivm,ifruit,icarbon) = biomass(:,ivm,ifruit,icarbon) - dturnover(:)
1754	turnover(:,ivm,ifruit,icarbon) = turnover(:,ivm,ifruit,icarbon) + dturnover(:)
1755
1756	! Tree level turnover (gC tree-1): leaves, roots and fruit
1757	DO icirc = 1,ncirc
1758
1759	circ_class_biomass(:,ivm,icirc,ileaf,icarbon) = circ_class_biomass(:,ivm,icirc,ileaf,icarbon) * &
1760	(un - leaf_frac(:,ivm,ilage) * turnover_rate(:))
1761	circ_class_biomass(:,ivm,icirc,iroot,icarbon) = circ_class_biomass(:,ivm,icirc,iroot,icarbon) * &
1762	(un - leaf_frac(:,ivm,ilage) * turnover_rate(:))
1763	circ_class_biomass(:,ivm,icirc,ifruit,icarbon) = circ_class_biomass(:,ivm,icirc,ifruit,icarbon) * &
1764	(un - leaf_frac(:,ivm,ilage) * turnover_rate(:))
1765
1766	ENDDO
1767
1768	! Crops
1769	ELSEIF (.NOT. natural(ivm)) THEN
1770
1771	! The sapwood acts as the structure to supports leaves and roots and store the
1772	! reserves and labile. Because the growing season is short for crops, the
1773	! turnover of tissues should be kept to a minimum or the GPP and NPP will be
1774	! too low.
1775	dturnover(:) = biomass(:,ivm,ileaf,icarbon) * leaf_frac(:,ivm,ilage) * &
1776	turnover_rate(:) / (wstress_month(:,ivm)**(1/tau_hum_month))
1777
1778	! Do we have enough biomass to satisfy the turnover?
1779	WHERE ( dturnover(:) .GE. biomass(:,ivm,ileaf,icarbon) )
1780
1781	! No longer account for wstress
1782	dturnover(:) = biomass(:,ivm,ileaf,icarbon) * &
1783	leaf_frac(:,ivm,ilage) * turnover_rate(:)
1784
1785	ENDWHERE
1786
1787	! Update biomass and turnover pools
1788	turnover(:,ivm,ileaf,icarbon) = turnover(:,ivm,ileaf,icarbon) + dturnover(:)
1789	biomass(:,ivm,ileaf,icarbon) = biomass(:,ivm,ileaf,icarbon) - dturnover(:)
1790
1791	! save leaf mass change
1792	delta_lm(:,ilage) = - dturnover(:)
1793
1794	dturnover(:) = biomass(:,ivm,iroot,icarbon) * dt / tau_eff_root(:,ivm)
1795	turnover(:,ivm,iroot,icarbon) = turnover(:,ivm,iroot,icarbon) + dturnover(:)
1796	biomass(:,ivm,iroot,icarbon) = biomass(:,ivm,iroot,icarbon) - dturnover(:)
1797
1798	dturnover(:) = biomass(:,ivm,ifruit,icarbon) * leaf_frac(:,ivm,ilage) * turnover_rate(:)
1799	turnover(:,ivm,ifruit,icarbon) = turnover(:,ivm,ifruit,icarbon) + dturnover(:)
1800	biomass(:,ivm,ifruit,icarbon) = biomass(:,ivm,ifruit,icarbon) - dturnover(:)
1801
1802	! Synchronize circ_class_biomass and biomass
1803	DO ipar = 1,nparts
1804
1805	WHERE (ind(:,ivm) .GT. zero)
1806
1807	circ_class_biomass(:,ivm,1,ipar,icarbon) = &
1808	biomass(:,ivm,ipar,icarbon) / ind(:,ivm)
1809
1810	ELSEWHERE
1811
1812	circ_class_biomass(:,ivm,1,ipar,icarbon) = zero
1813
1814	ENDWHERE
1815
1816	ENDDO ! nparts
1817
1818	! Grasses
1819	ELSEIF ( ( .NOT. is_tree(ivm) ) .AND. natural(ivm) ) THEN
1820
1821	! Grasses are now modeled as an evergreen biome, turnover is thus needed.
1822	! This section was implemented much later then the previous blocks and we
1823	! therefore followed our own advice and used tau_sap and tau_root rather
1824	! then having all biomass pool to follow leaf turnover.
1825	dturnover(:) = biomass(:,ivm,ileaf,icarbon) * leaf_frac(:,ivm,ilage) * &
1826	turnover_rate(:) / (wstress_month(:,ivm)**(1/tau_hum_month))
1827
1828	! Do we have enough biomass to satisfy the turnover?
1829	WHERE ( dturnover(:) .GE. biomass(:,ivm,ileaf,icarbon) )
1830
1831	! No longer account for wstress
1832	dturnover(:) = biomass(:,ivm,ileaf,icarbon) * &
1833	leaf_frac(:,ivm,ilage) * turnover_rate(:)
1834
1835	ENDWHERE
1836
1837	! Grasses are no longer considered natural and so it is assumed that the
1838	! aboveground turnover ends up in the harvest pool (mowing or grazing).
1839	! No idea whether this is anything close to reality but it is better to
1840	! assuming that there is no lateral transport. After testing over Europe
1841	! the annual harvest_pool for grasses was 128 Tg C year compared to the
1842	! estimated/observed 161 Tg C year. It is left to grassland people to further
1843	! improve this. Note that strictly speaking there is still no harvest, all
1844	! we did is moving the natural turnover (due to longevity) into the harvest
1845	! pool.
1846	harvest_pool(:,ivm,1,icarbon) = harvest_pool(:,ivm,1,icarbon) + &
1847	dturnover(:) * veget_max(:,ivm) * pixel_area(:)
1848	biomass(:,ivm,ileaf,icarbon) = biomass(:,ivm,ileaf,icarbon) - dturnover(:)
1849
1850	! save leaf mass change
1851	delta_lm(:,ilage) = - dturnover(:)
1852
1853	! Seeds are eaten
1854	dturnover(:) = biomass(:,ivm,ifruit,icarbon) * dt / tau_eff_leaf(:,ivm)
1855	harvest_pool(:,ivm,1,icarbon) = harvest_pool(:,ivm,1,icarbon) + &
1856	dturnover(:) * veget_max(:,ivm) * pixel_area(:)
1857	biomass(:,ivm,ifruit,icarbon) = biomass(:,ivm,ifruit,icarbon) - dturnover(:)
1858
1859	! Half the sapwood is eaten, half is added to the litter
1860	dturnover(:) = biomass(:,ivm,isapabove,icarbon) * dt / tau_eff_sap(:,ivm)
1861	harvest_pool(:,ivm,1,icarbon) = harvest_pool(:,ivm,1,icarbon) + &
1862	0.5 * dturnover(:) * veget_max(:,ivm) * pixel_area(:)
1863	turnover(:,ivm,isapabove,icarbon) = turnover(:,ivm,isapabove,icarbon) + &
1864	0.5 * dturnover(:)
1865	biomass(:,ivm,isapabove,icarbon) = biomass(:,ivm,isapabove,icarbon) - dturnover(:)
1866
1867	! Fill associated harvest variables
1868	harvest_type(:,ivm) = ifm_grass
1869	harvest_cut(:,ivm) = icut_grass
1870	harvest_area(:,ivm) = veget_max(:,ivm) * pixel_area(:)
1871
1872	! Make sure the harvest does not enter the turnover pool
1873	turnover(:,ivm,ileaf,icarbon) = zero
1874	turnover(:,ivm,ifruit,icarbon) = zero
1875
1876	! Root turnover at the stand level (gC m-2). Roots are added to the litter pool
1877	dturnover(:) = biomass(:,ivm,iroot,icarbon) * dt / tau_eff_root(:,ivm)
1878	turnover(:,ivm,iroot,icarbon) = turnover(:,ivm,iroot,icarbon) + dturnover(:)
1879	biomass(:,ivm,iroot,icarbon) = biomass(:,ivm,iroot,icarbon) - dturnover(:)
1880
1881	! Synchronize circ_class_biomass and biomass
1882	DO ipar = 1,nparts
1883
1884	WHERE (ind(:,ivm) .GT. zero)
1885
1886	circ_class_biomass(:,ivm,1,ipar,icarbon) = &
1887	biomass(:,ivm,ipar,icarbon) / ind(:,ivm)
1888
1889	ELSEWHERE
1890
1891	circ_class_biomass(:,ivm,1,ipar,icarbon) = zero
1892
1893	ENDWHERE
1894
1895	ENDDO ! nparts
1896
1897	ELSE
1898
1899	WRITE(numout,*) 'ERROR: turnover has not been defined'
1900
1901	ENDIF ! is_tree(ivm)
1902
1903	!---TEMP---
1904	IF (ivm .EQ. test_pft .AND. ld_trnov) THEN
1905	WRITE(numout,*) 'Check turnover - senesence'
1906	! DO ipar = 1,nparts
1907	! WRITE(numout,*) 'biomass vs circ_class_biomass (gC m-2), ', biomass(test_grid,ivm,ipar,icarbon), &
1908	! SUM(circ_class_biomass(test_grid,ivm,:,ipar,icarbon)*circ_class_n(test_grid,ivm,:))
1909	! ENDDO
1910	WRITE(numout,*) 'TOTAL biomass vs circ_class_biomass, ', SUM(biomass(test_grid,ivm,:,icarbon)), &
1911	SUM( SUM(circ_class_biomass(test_grid,ivm,:,:,icarbon),2) * circ_class_n(test_grid,ivm,:),1)
1912	ENDIF
1913	!----------
1914
1915	ENDDO
1916
1917	!! 4.3 Recalculate the fraction of leaf biomass in each leaf age class.
1918	! Older leaves will fall faster than younger leaves and therefore
1919	! the leaf age distribution needs to be recalculated after turnover.
1920	! The fraction of biomass in each leaf class is updated using the following equation:
1921	! \latexonly
1922	! \input{turnover_update_LeafAgeDistribution_eqn6.tex}
1923	! \endlatexonly
1924	! \n
1925	! new fraction = new leaf mass of that fraction / new total leaf mass
1926	! = (old fraction*old total leaf mass ::lm_old(:) + biomass change of that fraction ::delta_lm(:,m) ) /
1927	! new total leaf mass ::biomass(:,ivm,ileaf
1928	DO ilage = 1, nleafages
1929
1930	WHERE ( biomass(:,ivm,ileaf,icarbon) .GT. zero )
1931
1932	leaf_frac(:,ivm,ilage) = ( leaf_frac(:,ivm,ilage)*lm_old(:) + delta_lm(:,ilage) ) / biomass(:,ivm,ileaf,icarbon)
1933
1934	ELSEWHERE
1935
1936	leaf_frac(:,ivm,ilage) = zero
1937
1938	ENDWHERE
1939
1940	ENDDO
1941
1942	ENDDO ! loop over PFTs
1943
1944
1945	!! 5. Definitely drop all leaves if there is a very low leaf mass during senescence.
1946
1947	! Both for deciduous trees and grasses same conditions are checked:
1948	! If biomass is large enough (i.e. when it is greater than zero),
1949	! AND when senescence is set to true
1950	! AND the leaf biomass drops below a critical minimum biomass level (i.e. when it is lower than half
1951	! the minimum initial LAI ::lai_initmin(ivm) divided by the specific leaf area ::sla(ivm),
1952	! ::lai_initmin(ivm) is set to 0.3 in stomate_data.f90 and sla is a constant that is set to 0.015366 m2/gC),
1953	! If these conditions are met, the flag ::shed_rest(:) is set to TRUE.
1954	!
1955	! After this, the biomass of different carbon pools both for trees and grasses is set to zero
1956	! and the mean leaf age is reset to zero.
1957	! Finally, the leaf fraction and leaf age of the different leaf age classes is set to zero.
1958	DO ivm = 2,nvm ! Loop over # PFTs
1959
1960	shed_rest(:) = .FALSE.
1961
1962	!! 5.1 For deciduous trees: next to leaves, also fruits and fine roots are dropped
1963	! For deciduous trees: next to leaves, also fruits and fine roots are dropped: fruit ::biomass(:,ivm,ifruit)
1964	! and fine root ::biomass(:,ivm,iroot) carbon pools are set to zero.
1965	IF ( is_tree(ivm) .AND. ( senescence_type(ivm) .NE. 'none' ) ) THEN
1966
1967	! Check whether we shed the remaining leaves. The condition depends on biomass(ileaf)
1968	! so first calculate the sheding at the tree level. because biomass(ileaf) is not
1969	! changed at the tree level the same statement can then be used at the stand level.
1970	DO icirc = 1,ncirc
1971
1972	! check whether we shed the remaining leaves
1973	WHERE ( ( biomass(:,ivm,ileaf,icarbon) .GT. zero ) .AND. senescence(:,ivm) .AND. &
1974	( biomass(:,ivm,ileaf,icarbon) .LT. (lai_initmin(ivm) / 2.)/sla(ivm) ) )
1975
1976	! Tree level (gC tree-1)
1977	circ_class_biomass(:,ivm,icirc,ileaf,icarbon) = zero
1978	circ_class_biomass(:,ivm,icirc,iroot,icarbon) = zero
1979	circ_class_biomass(:,ivm,icirc,ifruit,icarbon) = zero
1980
1981	ENDWHERE
1982
1983	ENDDO !ncirc
1984
1985	! check whether we shed the remaining leaves
1986	WHERE ( ( biomass(:,ivm,ileaf,icarbon) .GT. zero ) .AND. senescence(:,ivm) .AND. &
1987	( biomass(:,ivm,ileaf,icarbon) .LT. (lai_initmin(ivm) / 2.)/sla(ivm) ) )
1988
1989	shed_rest(:) = .TRUE.
1990
1991	! Stand level (gc m-2)
1992	turnover(:,ivm,ileaf,icarbon) = turnover(:,ivm,ileaf,icarbon) + biomass(:,ivm,ileaf,icarbon)
1993	turnover(:,ivm,iroot,icarbon) = turnover(:,ivm,iroot,icarbon) + biomass(:,ivm,iroot,icarbon)
1994	turnover(:,ivm,ifruit,icarbon) = turnover(:,ivm,ifruit,icarbon) + biomass(:,ivm,ifruit,icarbon)
1995	biomass(:,ivm,ileaf,icarbon) = zero
1996	biomass(:,ivm,iroot,icarbon) = zero
1997	biomass(:,ivm,ifruit,icarbon) = zero
1998
1999	! reset leaf age
2000	leaf_meanage(:,ivm) = zero
2001
2002	ENDWHERE
2003
2004	ENDIF
2005
2006	!! 5.2 For grasses: all aboveground carbon pools, except the carbohydrate reserves are affected:
2007	! For grasses: all aboveground carbon pools, except the carbohydrate reserves are affected:
2008	! fruit ::biomass(:,ivm,ifruit,icarbon), fine root ::biomass(:,ivm,iroot,icarbon) and sapwood above
2009	! ::biomass(:,ivm,isapabove,icarbon) carbon pools are set to zero.
2010	IF ( .NOT. is_tree(ivm) ) THEN
2011
2012	IF (senescence_type(ivm) .EQ. 'crop') THEN
2013
2014	!!$ ! Shed the remaining leaves if LAI very low.
2015	!!$ WHERE ( ( biomass(:,ivm,ileaf,icarbon) .GT. zero ) .AND. senescence(:,ivm) .AND. &
2016	!!$ ( biomass(:,ivm,ileaf,icarbon) .LT. (lai_initmin(ivm) / deux) / sla(ivm) ) )
2017	!!$
2018	!!$ shed_rest(:) = .TRUE.
2019	!!$
2020	!!$ ! Crops are planted every year. So make sure to remove everything at the end of senescence
2021	!!$ ! which for crops is actually harvest. Next stomate_prescribe should prescribe a new crop.
2022	!!$ ! If sapwood is not removed, some is left on site and the model will try to grow a crop.
2023	!!$ ! However, there are not enough reserves left so the crop will die. This results in one year
2024	!!$ ! of normal crop growth (following the prescribe) and then one year of almost no growth (the
2025	!!$ ! year starting from the too low reserves).
2026	!!$ turnover(:,ivm,ileaf,icarbon) = turnover(:,ivm,ileaf,icarbon) + biomass(:,ivm,ileaf,icarbon) + &
2027	!!$ biomass(:,ivm,icarbres,icarbon) + biomass(:,ivm,ilabile,icarbon)
2028	!!$ turnover(:,ivm,iroot,icarbon) = turnover(:,ivm,iroot,icarbon) + biomass(:,ivm,iroot,icarbon)
2029	!!$ turnover(:,ivm,ifruit,icarbon) = turnover(:,ivm,ifruit,icarbon) + biomass(:,ivm,ifruit,icarbon)
2030	!!$ turnover(:,ivm,isapabove,icarbon) = turnover(:,ivm,isapabove,icarbon) + biomass(:,ivm,isapabove,icarbon) + &
2031	!!$ biomass(:,ivm,isapbelow,icarbon)
2032	!!$ turnover(:,ivm,isapabove,icarbon) = turnover(:,ivm,isapabove,icarbon) + biomass(:,ivm,iheartabove,icarbon) + &
2033	!!$ biomass(:,ivm,iheartbelow,icarbon)
2034	!!$ biomass(:,ivm,isapabove,icarbon) = zero
2035	!!$ biomass(:,ivm,ileaf,icarbon) = zero
2036	!!$ biomass(:,ivm,iroot,icarbon) = zero
2037	!!$ biomass(:,ivm,ifruit,icarbon) = zero
2038	!!$ biomass(:,ivm,icarbres,icarbon) = zero
2039	!!$ biomass(:,ivm,ilabile,icarbon) = zero
2040	!!$
2041	!!$ ! To be sure put overlooked residuals also in the turnover pool
2042	!!$ biomass(:,ivm,iheartbelow,icarbon) = zero
2043	!!$ biomass(:,ivm,iheartabove,icarbon) = zero
2044	!!$ biomass(:,ivm,isapbelow,icarbon) = zero
2045	!!$
2046	!!$ ! No individuals left
2047	!!$ ind(:,ivm) = zero
2048	!!$
2049	!!$ ! reset leaf age
2050	!!$ leaf_meanage(:,ivm) = zero
2051	!!$
2052	!!$ ENDWHERE
2053
2054	ELSEIF (senescence_type(ivm) .EQ. 'mixed') THEN
2055
2056	! Shed the remaining leaves if LAI very low.
2057	WHERE ( ( biomass(:,ivm,ileaf,icarbon) .GT. zero ) .AND. senescence(:,ivm) .AND. &
2058	( biomass(:,ivm,ileaf,icarbon) .LT. (lai_initmin(ivm) / deux) / sla(ivm) ) )
2059
2060	shed_rest(:) = .TRUE.
2061
2062	! Grasses should live on after senescence. Do not shed the sapwood because then the whole
2063	! system collapse because the reserves are calculated based on the sapwood. no sapwood = no
2064	! reserves = no phenology = no growth.
2065	turnover(:,ivm,ileaf,icarbon) = turnover(:,ivm,ileaf,icarbon) + biomass(:,ivm,ileaf,icarbon) + &
2066	biomass(:,ivm,icarbres,icarbon) + biomass(:,ivm,ilabile,icarbon)
2067	turnover(:,ivm,iroot,icarbon) = turnover(:,ivm,iroot,icarbon) + biomass(:,ivm,iroot,icarbon)
2068	turnover(:,ivm,ifruit,icarbon) = turnover(:,ivm,ifruit,icarbon) + biomass(:,ivm,ifruit,icarbon)
2069	biomass(:,ivm,ileaf,icarbon) = zero
2070	biomass(:,ivm,iroot,icarbon) = zero
2071	biomass(:,ivm,ifruit,icarbon) = zero
2072	biomass(:,ivm,icarbres,icarbon) = zero
2073	biomass(:,ivm,ilabile,icarbon) = zero
2074
2075	! Note that grasses stay alive so do not change the number of individuals!
2076
2077	! reset leaf age
2078	leaf_meanage(:,ivm) = zero
2079
2080	ENDWHERE
2081
2082	ELSEIF (senescence_type(ivm) .EQ. 'none') THEN
2083
2084	! Nothing should be done
2085
2086	ELSE
2087
2088	WRITE(numout,*) 'ERROR: senescence type is not known'
2089
2090	ENDIF
2091
2092	!---TEMP---
2093	IF (ivm .EQ. test_pft .AND. ld_pheno) THEN
2094	WRITE(numout,*) 'phenology, shed_rest, ',shed_rest(:)
2095	ENDIF
2096	!----------
2097
2098	! Synchronize circ_class_biomass and biomass
2099	DO ipar = 1,nparts
2100
2101	WHERE (ind(:,ivm) .GT. zero)
2102
2103	circ_class_biomass(:,ivm,1,ipar,icarbon) = &
2104	biomass(:,ivm,ipar,icarbon) / ind(:,ivm)
2105
2106	ELSEWHERE
2107
2108	circ_class_biomass(:,ivm,1,ipar,icarbon) = zero
2109
2110	ENDWHERE
2111
2112	ENDDO ! nparts
2113
2114	ENDIF ! is_tree(ivm)
2115
2116	!---TEMP---
2117	IF (ivm .EQ. test_pft .AND. ld_trnov) THEN
2118	WRITE(numout,*) 'Check turnover - leaf sheding'
2119	WRITE(numout,*) 'ifruit biomass vs circ_class_biomass (gC m-2), ', biomass(test_grid,ivm,ifruit,icarbon), &
2120	SUM(circ_class_biomass(test_grid,ivm,:,ifruit,icarbon)*circ_class_n(test_grid,ivm,:))
2121
2122	WRITE(numout,*) 'TOTAL biomass vs circ_class_biomass, ', SUM(biomass(test_grid,ivm,:,icarbon)), &
2123	SUM( SUM(circ_class_biomass(test_grid,ivm,:,:,icarbon),2) * circ_class_n(test_grid,ivm,:),1)
2124	ENDIF
2125	!----------
2126
2127	!! 5.3 Reset the leaf age structure: the leaf fraction and leaf age
2128	! of the different leaf age classes is set to zero.
2129	DO ilage = 1, nleafages
2130
2131	WHERE ( shed_rest(:) )
2132
2133	leaf_age(:,ivm,ilage) = zero
2134	leaf_frac(:,ivm,ilage) = zero
2135
2136	ENDWHERE
2137
2138	ENDDO
2139
2140	ENDDO ! loop over PFTs
2141
2142
2143	!! 6. Herbivore activity: elephants, cows, gazelles but no lions.
2144
2145	! Herbivore activity affects the biomass of leaves and fruits as well
2146	! as stalks (only for grasses). Herbivore activity does not modify leaf
2147	! age structure. Herbivores ::herbivores(:,ivm) is the time constant of
2148	! probability of a leaf to be eaten by a herbivore, and is calculated in
2149	! ::stomate_season. following Mc Naughton et al. [1989].
2150	IF ( control%ok_herbivory ) THEN
2151
2152	! If the herbivore activity is allowed (if ::control%ok_herbivory is true, which is set in run.def),
2153	! remove the amount of biomass consumed by herbivory from the leaf biomass ::biomass(:,ivm,ileaf,icarbon) and
2154	! the fruit biomass ::biomass(:,ivm,ifruit,icarbon).
2155	! The daily amount consumed equals the biomass multiplied by 1 day divided by the time constant ::herbivores(:,ivm).
2156	DO ivm = 2,nvm ! Loop over # PFTs
2157
2158	IF ( is_tree(ivm) ) THEN
2159
2160	!---TEMP---
2161	IF (ivm .EQ. test_pft .AND. ld_trnov) THEN
2162	WRITE(numout,*) 'herbivores, ',herbivores
2163	ENDIF
2164	!---------
2165
2166	!! For trees: only the leaves and fruit carbon pools are affected
2167	!+++CHECK+++
2168	! Why is root herbivory not accounted for. Could be difficult to parameterize but if
2169	! leaf browsing is accounted for, there is no reason to ignore root browsing
2170	WHERE (biomass(:,ivm,ileaf,icarbon) .GT. zero)
2171
2172	WHERE (herbivores(:,ivm).GT. zero)
2173
2174	! Stand level (gC m-2)
2175	dturnover(:) = biomass(:,ivm,ileaf,icarbon) * dt / herbivores(:,ivm)
2176	biomass(:,ivm,ileaf,icarbon) = biomass(:,ivm,ileaf,icarbon) - dturnover(:)
2177	turnover(:,ivm,ileaf,icarbon) = turnover(:,ivm,ileaf,icarbon) + dturnover(:)
2178
2179	dturnover(:) = biomass(:,ivm,ifruit,icarbon) * dt / herbivores(:,ivm)
2180	biomass(:,ivm,ifruit,icarbon) = biomass(:,ivm,ifruit,icarbon) - dturnover(:)
2181	turnover(:,ivm,ifruit,icarbon) = turnover(:,ivm,ifruit,icarbon) + dturnover(:)
2182
2183	ENDWHERE
2184
2185	ENDWHERE
2186
2187	! Tree level herbivory
2188	DO icirc = 1,ncirc
2189
2190	!! For trees: only the leaves and fruit carbon pools are affected
2191	WHERE (biomass(:,ivm,ileaf,icarbon) .GT. zero)
2192
2193	WHERE (herbivores(:,ivm).GT. zero)
2194	! Tree level (gC tree-1)
2195	circ_class_biomass(:,ivm,icirc,ileaf,icarbon) = circ_class_biomass(:,ivm,icirc,ileaf,icarbon) * &
2196	(un - dt / herbivores(:,ivm))
2197	circ_class_biomass(:,ivm,icirc,ifruit,icarbon) = circ_class_biomass(:,ivm,icirc,ifruit,icarbon) * &
2198	(un - dt / herbivores(:,ivm))
2199	ENDWHERE
2200
2201	ENDWHERE
2202
2203	ENDDO ! ncirc
2204	!+++++++++
2205
2206	! grasses and croplands
2207	ELSE
2208
2209	! For grasses: all aboveground carbon pools are affected: leaves, fruits and sapwood above
2210	WHERE ( biomass(:,ivm,ileaf,icarbon) .GT. zero )
2211
2212	WHERE (herbivores(:,ivm) .GT. zero)
2213
2214	dturnover(:) = biomass(:,ivm,ileaf,icarbon) * dt / herbivores(:,ivm)
2215	turnover(:,ivm,ileaf,icarbon) = turnover(:,ivm,ileaf,icarbon) + dturnover(:)
2216	biomass(:,ivm,ileaf,icarbon) = biomass(:,ivm,ileaf,icarbon) - dturnover(:)
2217
2218	dturnover(:) = biomass(:,ivm,isapabove,icarbon) * dt / herbivores(:,ivm)
2219	turnover(:,ivm,isapabove,icarbon) = turnover(:,ivm,isapabove,icarbon) + dturnover(:)
2220	biomass(:,ivm,isapabove,icarbon) = biomass(:,ivm,isapabove,icarbon) - dturnover(:)
2221
2222	dturnover(:) = biomass(:,ivm,ifruit,icarbon) * dt / herbivores(:,ivm)
2223	turnover(:,ivm,ifruit,icarbon) = turnover(:,ivm,ifruit,icarbon) + dturnover(:)
2224	biomass(:,ivm,ifruit,icarbon) = biomass(:,ivm,ifruit,icarbon) - dturnover(:)
2225
2226	ENDWHERE
2227
2228	ENDWHERE
2229
2230	! Synchronize circ_class_biomass and biomass
2231	DO ipar = 1,nparts
2232
2233	WHERE (ind(:,ivm) .GT. zero)
2234
2235	circ_class_biomass(:,ivm,1,ipar,icarbon) = &
2236	biomass(:,ivm,ipar,icarbon) / ind(:,ivm)
2237
2238	ELSEWHERE
2239
2240	circ_class_biomass(:,ivm,1,ipar,icarbon) = zero
2241
2242	ENDWHERE
2243
2244	ENDDO ! nparts
2245
2246	ENDIF ! tree/grass?
2247
2248	!---TEMP---
2249	IF (ivm .EQ. test_pft .AND. ld_trnov) THEN
2250	WRITE(numout,*) 'Check turnover - herbivory '
2251	WRITE(numout,*) 'TOTAL biomass vs circ_class_biomass, ', (biomass(test_grid,ivm,:,icarbon)), &
2252	( SUM(circ_class_biomass(test_grid,ivm,:,:,icarbon),2) * circ_class_n(test_grid,ivm,:))
2253	ENDIF
2254	!----------
2255
2256	ENDDO ! loop over PFTs
2257
2258	ENDIF ! end herbivores
2259
2260
2261	!! 7. Fruit turnover for trees
2262
2263	! Fruit turnover for trees: trees simply lose their fruits
2264	! with a time constant ::tau_fruit(ivm),
2265	! that is set to 90 days for all PFTs in ::constants_mtc
2266	! Again, we assume that the trees do not lose either of the
2267	! reserve pools (ilabile or icarbres).
2268	DO iele = 1,nelements
2269
2270	DO ivm = 2,nvm ! Loop over # PFTs
2271
2272	IF ( is_tree(ivm) ) THEN
2273
2274	! Stand level (gC m-2)
2275	dturnover(:) = biomass(:,ivm,ifruit,iele) * dt / tau_fruit(ivm)
2276	biomass(:,ivm,ifruit,iele) = biomass(:,ivm,ifruit,iele) - dturnover(:)
2277	turnover(:,ivm,ifruit,iele) = turnover(:,ivm,ifruit,iele) + dturnover(:)
2278
2279	! Tree level (gC tree-1)
2280	circ_class_biomass(:,ivm,:,ifruit,iele) = circ_class_biomass(:,ivm,:,ifruit,iele) * &
2281	(un - dt / tau_fruit(ivm))
2282
2283	!---TEMP---
2284	IF (ivm .EQ. test_pft .AND. ld_trnov) THEN
2285	WRITE(numout,*) 'Check turnover - fruits'
2286	WRITE(numout,*) 'biomass(ifruit), circ_class_biomass(ifruit), ',biomass(test_grid,ivm,ifruit,icarbon), &
2287	SUM(circ_class_biomass(test_grid,ivm,:,ifruit,icarbon)*circ_class_n(test_grid,ivm,:))
2288	WRITE(numout,*) 'turnover, ', dt/tau_fruit(ivm)
2289	ENDIF
2290	!----------
2291
2292	ENDIF ! is_tree
2293
2294	!---TEMP---
2295	IF (ivm .EQ. test_pft .AND. ld_trnov) THEN
2296	WRITE(numout,*) 'biomass vs circ_class_biomass (gC m-2), ', biomass(test_grid,ivm,ifruit,icarbon), &
2297	SUM(circ_class_biomass(test_grid,ivm,:,ifruit,icarbon)*circ_class_n(test_grid,ivm,:))
2298	WRITE(numout,*) 'TOTAL biomass vs circ_class_biomass, ', SUM(biomass(test_grid,ivm,:,icarbon)), &
2299	SUM( SUM(circ_class_biomass(test_grid,ivm,:,:,icarbon),2) * circ_class_n(test_grid,ivm,:),1)
2300	ENDIF
2301	!----------
2302
2303	ENDDO ! loop over PFTs
2304
2305	END DO ! nelements
2306
2307
2308	!! 8. Conversion of sapwood to heartwood both for aboveground and belowground sapwood and heartwood.
2309
2310	! Following LPJ (Sitch et al., 2003), sapwood biomass is converted into heartwood biomass
2311	! with a time constant tau ::tau_sap(j) of 1 year.
2312	! Note that this biomass conversion is not added to "turnover" as the biomass is not lost!
2313	DO ivm = 2,nvm ! Loop over # PFTs
2314
2315	IF ( is_tree(ivm) ) THEN
2316
2317	! For the recalculation of age in 9.2 (in case the vegetation is not dynamic ie. ::control%ok_dgvm is FALSE),
2318	! the heartwood above and below is stored in ::hw_old(:).
2319	IF ( .NOT. control%ok_dgvm ) THEN
2320
2321	hw_old(:) = biomass(:,ivm,iheartabove,icarbon) + biomass(:,ivm,iheartbelow,icarbon)
2322	sw_old(:) = biomass(:,ivm,isapabove,icarbon) + biomass(:,ivm,isapbelow,icarbon)
2323
2324	ENDIF
2325
2326	!! 8.1 Calculate the rate of sapwood to heartwood conversion
2327	! Calculate the rate of sapwood to heartwood conversion with the
2328	! time constant ::tau_sap(ivm) and update aboveground and
2329	! belowground sapwood ::biomass(:,ivm,isapabove) and ::biomass(:,ivm,isapbelow)
2330	! and heartwood ::biomass(:,ivm,iheartabove) and ::biomass(:,ivm,iheartbelow).
2331	DO iele = 1,nelements
2332
2333	! Stand level above the ground (gC m-2)
2334	sapconv(:) = biomass(:,ivm,isapabove,iele) * dt / tau_eff_sap(:,ivm)
2335	biomass(:,ivm,iheartabove,iele) = biomass(:,ivm,iheartabove,iele) + sapconv(:)
2336	biomass(:,ivm,isapabove,iele) = biomass(:,ivm,isapabove,iele) - sapconv(:)
2337
2338	! Stand level below the ground (gC m-2)
2339	sapconv(:) = biomass(:,ivm,isapbelow,iele) * dt / tau_eff_sap(:,ivm)
2340	biomass(:,ivm,iheartbelow,iele) = biomass(:,ivm,iheartbelow,iele) + sapconv(:)
2341	biomass(:,ivm,isapbelow,iele) = biomass(:,ivm,isapbelow,iele) - sapconv(:)
2342
2343	DO icirc = 1,ncirc
2344	! Tree level above and belowground (gC tree-1)
2345	circ_class_biomass(:,ivm,icirc,iheartabove,iele) = &
2346	circ_class_biomass(:,ivm,icirc,iheartabove,iele) + &
2347	circ_class_biomass(:,ivm,icirc,isapabove,iele) * dt / &
2348	tau_eff_sap(:,ivm)
2349	circ_class_biomass(:,ivm,icirc,isapabove,iele) = &
2350	circ_class_biomass(:,ivm,icirc,isapabove,iele) * &
2351	(un - dt / tau_eff_sap(:,ivm))
2352	circ_class_biomass(:,ivm,icirc,iheartbelow,iele) = &
2353	circ_class_biomass(:,ivm,icirc,iheartbelow,iele) + &
2354	circ_class_biomass(:,ivm,icirc,isapbelow,iele) * dt / &
2355	tau_eff_sap(:,ivm)
2356	circ_class_biomass(:,ivm,icirc,isapbelow,iele) = &
2357	circ_class_biomass(:,ivm,icirc,isapbelow,iele) * &
2358	(un - dt / tau_eff_sap(:,ivm))
2359
2360	ENDDO
2361
2362	ENDDO
2363
2364	!---TEMP---
2365	IF (ivm .EQ. test_pft .AND. ld_trnov) THEN
2366	WRITE(numout,*) 'Check turnover - sap & heartwood'
2367	! DO ipar = 1, nparts
2368	! WRITE(numout,*) 'biomass vs circ_class_biomass (gC m-2), ', biomass(test_grid,ivm,ipar,icarbon), &
2369	! SUM(circ_class_biomass(test_grid,ivm,:,ipar,icarbon)*circ_class_n(test_grid,ivm,:))
2370	! ENDDO
2371	WRITE(numout,*) 'TOTAL biomass vs circ_class_biomass, ', SUM(biomass(test_grid,ivm,:,icarbon)), &
2372	SUM( SUM(circ_class_biomass(test_grid,ivm,:,:,icarbon),2) * circ_class_n(test_grid,ivm,:),1)
2373	ENDIF
2374	!----------
2375
2376	!! 8.2 If the vegetation is not dynamic, the age of the plant is decreased.
2377	! The updated heartwood, the sum of new heartwood above and new heartwood below after
2378	! converting sapwood to heartwood, is saved as ::hw_new(:) .
2379	! Creation of new heartwood decreases the age of the plant with a factor that is determined by:
2380	! old heartwood ::hw_old(:) divided by the new heartwood ::hw_new(:)
2381	IF ( .NOT. control%ok_dgvm ) THEN
2382
2383	hw_new(:) = biomass(:,ivm,iheartabove,icarbon) + biomass(:,ivm,iheartbelow,icarbon)
2384	sw_new(:) = biomass(:,ivm,isapabove,icarbon) + biomass(:,ivm,isapbelow,icarbon)
2385
2386	WHERE ( hw_new(:) .GT. zero )!
2387
2388	age(:,ivm) = age(:,ivm) * hw_old(:)/hw_new(:)
2389
2390	! Sapwood age doesn't seem to be used, but it causes problems because it is not initialised in this
2391	! routine and it is a local variable
2392	!!$ sapwood_age(:,ivm) = sapwood_age(:,ivm) * sw_old(:)/sw_new(:)
2393
2394	ENDWHERE
2395
2396	ENDIF
2397
2398	ENDIF
2399
2400	ENDDO ! loop over PFTs
2401
2402	!---TEMP---
2403	IF (ld_trnov) THEN
2404	WRITE(numout,*) 'leaf_meanage in turnover', leaf_meanage(:,test_pft)
2405	WRITE(numout,*) 'leaf_age in turnover', leaf_age(:,test_pft,:)
2406	WRITE(numout,*) 'leaf_frac in turnover', leaf_frac(:,test_pft,:)
2407	ENDIF
2408	!----------
2409
2410	! Write mean leaf age and time constant of probability of a leaf to be eaten by a herbivore
2411	! to the stomate output file.
2412	CALL histwrite_p (hist_id_stomate, 'LEAFAGE', itime, &
2413	leaf_meanage, npts*nvm, horipft_index)
2414	CALL histwrite_p (hist_id_stomate, 'HERBIVORES', itime, &
2415	herbivores, npts*nvm, horipft_index)
2416
2417
2418	!! 10. Check mass balance closure
2419
2420	!! 10.1 Calculate pools at the end of the routine
2421	pool_end = zero
2422	DO ipar = 1,nparts
2423	DO iele = 1,nelements
2424	pool_end(:,:,iele) = pool_end(:,:,iele) + &
2425	(biomass(:,:,ipar,iele) * veget_max(:,:)) + &
2426	(turnover(:,:,ipar,iele) * veget_max(:,:))
2427	ENDDO
2428	ENDDO
2429
2430	! The biomass harvest pool is expressed in gC pixel-1 So, it
2431	! shouldn't be multiplied by veget_max but it should be divided
2432	! by pixel_area to obtain gC m-2.
2433	DO ivm = 1,nvm
2434	pool_end(:,ivm,icarbon) = pool_end(:,ivm,icarbon) + &
2435	SUM(harvest_pool(:,ivm,:,icarbon),2) / &
2436	pixel_area(:)
2437	ENDDO
2438
2439	!! 10.2 Calculate components of the mass balance
2440	check_intern(:,:,iatm2land,icarbon) = zero
2441	check_intern(:,:,iland2atm,icarbon) = -un * zero
2442	check_intern(:,:,ilat2out,icarbon) = zero
2443	check_intern(:,:,ilat2in,icarbon) = -un * zero
2444	check_intern(:,:,ipoolchange,icarbon) = -un * (pool_end(:,:,icarbon) - pool_start(:,:,icarbon))
2445	closure_intern = zero
2446	DO imbc = 1,nmbcomp
2447	closure_intern(:,:,icarbon) = closure_intern(:,:,icarbon) + check_intern(:,:,imbc,icarbon)
2448	ENDDO
2449
2450	! Write outcome
2451	DO ipts = 1,npts
2452	DO ivm = 1,nvm
2453	IF (closure_intern(ipts,ivm,icarbon) .LT. min_stomate .AND. &
2454	closure_intern(ipts,ivm,icarbon) .GT. -min_stomate) THEN
2455	IF (ld_massbal) WRITE(numout,*) 'Mass balance closure: stomate_turnover.f90', ipts, ivm
2456	ELSE
2457	WRITE(numout,*) 'Error: mass balance is not closed in stomate_turnover.f90'
2458	WRITE(numout,*) ' ipts,ivm; ', ipts,ivm
2459	WRITE(numout,*) ' Difference is, ', closure_intern(ipts,ivm,icarbon)
2460	WRITE(numout,*) ' pool_end,pool_start: ', pool_end(ipts,ivm,icarbon), &
2461	pool_start(ipts,ivm,icarbon)
2462	ENDIF
2463	ENDDO
2464	ENDDO
2465
2466	IF (bavard.GE.4) WRITE(numout,*) 'Leaving prognostic turnover'
2467
2468	END SUBROUTINE turnover_prognostic
2469
2470	END MODULE stomate_turnover

Note: See TracBrowser for help on using the repository browser.

Download in other formats: