Context Navigation

source: tags/ORCHIDEE_4_1/ORCHIDEE/src_stomate/stomate_turnover.f90 @ 7852

Last change on this file since 7852 was 7615, checked in by sebastiaan.luyssaert, 2 years ago
Enhanced consistency of variable names: input has been changed in n_input throughout the code and the variable name vegstress introduced in sechiba is now also used in stomate. Enhnaced computational consistency: Pgap_cumul is used in stomate rather than recalculating it before calculating light_tran_to_floor_season. Edited getin_p while checking the code (but no real changes were made) and added several missing stomate and sechiba variables to age_class_distr to improve the 1+1=2 issue when comparing a model run with against a run without age classes. Finally: Write warning 10b in allocation to the history file rather than the out_orchidee file
Property svn:keywords set to `HeadURL Date Author Revision`
File size: 110.2 KB

Line
1	! =================================================================================================================================
2	! MODULE : stomate_turnover.f90
3	!
4	! CONTACT : orchidee-help _at_ listes.ipsl.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	!! \n
22	!_ ================================================================================================================================
23
24	MODULE stomate_turnover
25
26	! modules used:
27	USE xios_orchidee
28	USE ioipsl_para
29	USE stomate_data
30	USE constantes
31	USE grid
32	USE pft_parameters
33	USE sapiens_agriculture, ONLY: crop_harvest
34	USE function_library, ONLY : biomass_to_lai, &
35	check_vegetation_area, check_mass_balance, lai_to_biomass
36	USE time, ONLY : julian_diff
37
38	IMPLICIT NONE
39
40	! private & public routines
41
42	PRIVATE
43	PUBLIC turn_over, turnover_clear, drought_mortality
44
45	LOGICAL, SAVE :: firstcall_turnover = .TRUE. !! first call (true/false)
46	!$OMP THREADPRIVATE(firstcall_turnover)
47	INTEGER(i_std), SAVE :: printlev_loc !! Local level of text output for current module
48	!$OMP THREADPRIVATE(printlev_loc)
49	CONTAINS
50
51
52	!! ================================================================================================================================
53	!! SUBROUTINE : turnover_clear
54	!!
55	!>\BRIEF Set flag ::firstcall_turnover to .TRUE., and therefore activate section 1
56	!! of subroutine turn which writes a message to the output.
57	!!
58	!_ ================================================================================================================================
59
60	SUBROUTINE turnover_clear
61	firstcall_turnover=.TRUE.
62	END SUBROUTINE turnover_clear
63
64
65	!! ================================================================================================================================
66	!! SUBROUTINE : turn_over
67	!!
68	!>\BRIEF Calculate turnover of leaves, roots, fruits and sapwood
69	!! due to aging or climatic induced senescence. Calculate
70	!! herbivory.
71	!!
72	!! DESCRIPTION : This subroutine determines the turnover of leaves and fine
73	!! roots (and stems for grasses) and simulates following processes:
74	!! 1. Mean leaf age is calculated from leaf ages of separate leaf age
75	!! classes. Should actually be recalculated at the end of the routine,
76	!! but it does not change too fast. The mean leaf age is calculated
77	!! using the following equation:
78	!! \latexonly
79	!! \input{turnover_lma_update_eqn1.tex}
80	!! \endlatexonly
81	!! \n
82	!! 2. Meteorological senescence: the detection of the end of the growing
83	!! season and shedding of leaves, fruits and fine roots due to
84	!! unfavourable meteorological conditions. The model distinguishes
85	!! three different types of "climatic" leaf senescence, that do not
86	!! change the age structure: sensitivity to cold temperatures, to lack
87	!! of water, or both. If meteorological conditions are fulfilled, a
88	!! flag ::senescence is set to TRUE. Note that evergreen species do not
89	!! experience climatic senescence. Climatic senescence is triggered by
90	!! sensitivity to cold temperatures where the critical temperature for
91	!! senescence is calculated using the following equation:
92	!! \latexonly
93	!! \input{turnover_temp_crit_eqn2.tex}
94	!! \endlatexonly
95	!! \n
96	!! Climatic senescence is triggered by sensitivity to lack of water
97	!! availability where the moisture availability critical level is
98	!! calculated using the following equation:
99	!! \latexonly
100	!! \input{turnover_moist_crit_eqn3.tex}
101	!! \endlatexonly
102	!! \n
103	!! Climatic senescence is triggered by sensitivity to temperature or to
104	!! lack of water where critical temperature and moisture availability
105	!! are calculated as above.\n
106	!! Trees in climatic senescence lose their fine roots at the same rate
107	!! as they lose their leaves. The rate of biomass loss of both fine
108	!! roots and leaves is presribed through the equation:
109	!! \latexonly
110	!! \input{turnover_clim_senes_biomass_eqn4.tex}
111	!! \endlatexonly
112	!! \n
113	!! with ::leaffall(j) a PFT-dependent time constant which is given in
114	!! ::stomate_constants. In grasses, leaf senescence is extended to
115	!! the whole plant (all carbon pools) except to its carbohydrate reserve.
116	!! 3. Senescence due to aging: the loss of leaves, fruits and biomass due
117	!! to aging At a certain age, leaves fall off, even if the climate
118	!! would allow a green plant all year round. Even if the meteorological
119	!! conditions are favorable for leaf maintenance, plants, and in
120	!! particular, evergreen trees, have to renew their leaves simply because
121	!! the old leaves become inefficient. Roots, fruits (and stems for
122	!! grasses) follow leaves. The ??senescence?? rate varies with leaf
123	!! age. Note that plant is not declared senescent in this case (wchich
124	!! is important for allocation: if the plant loses leaves because of
125	!! their age, it can renew them). The leaf turnover rate due to aging
126	!! of leaves is calculated using the following equation:
127	!! \latexonly
128	!! \input{turnover_age_senes_biomass_eqn5.tex}
129	!! \endlatexonly
130	!! \n
131	!! Drop all leaves if there is a very low leaf mass during senescence.
132	!! After this, the biomass of different carbon pools both for trees
133	!! and grasses is set to zero and the mean leaf age is reset to zero.
134	!! Finally, the leaf fraction and leaf age of the different leaf age
135	!! classes is set to zero. For deciduous trees: next to leaves, also
136	!! fruits and fine roots are dropped.
137	!! For grasses: all aboveground carbon pools, except the carbohydrate
138	!! reserves are affected:
139	!! 4. Update the leaf biomass, leaf age class fraction and the LAI
140	!! Older leaves will fall more frequently than younger leaves and
141	!! therefore the leaf age distribution needs to be recalculated after
142	!! turnover. The fraction of biomass in each leaf class is updated using
143	!! the following equation:
144	!! \latexonly
145	!! \input{turnover_update_LeafAgeDistribution_eqn6.tex}
146	!! \endlatexonly
147	!! \n
148	!! 5. Simulate herbivory activity and update leaf and fruits biomass.
149	!! Herbivore activity affects the biomass of leaves and fruits as well
150	!! as stalks (only for grasses). However, herbivores do not modify leaf
151	!! age structure.
152	!! 6. Calculates fruit turnover for trees. Trees simply lose their fruits
153	!! with a time constant ::longevity_fruit(j), that is set to 90 days for all
154	!! PFTs in ::stomate_constants
155	!! 7. Convert sapwood to heartwood for trees and update heart and softwood
156	!! above and belowground biomass. Sapwood biomass is converted into
157	!! heartwood biomass with a time constant tau ::longevity_sap(j) of 1 year.
158	!! Note that this biomass conversion is not added to "turnover" as the
159	!! biomass is not lost. For the updated heartwood, the sum of new
160	!! heartwood above and new heartwood below after converting sapwood to
161	!! heartwood, is saved as ::hw_new(:). Creation of new heartwood
162	!! decreases the age of the plant ??carbon?? with a factor that is
163	!! determined by: old heartwood ::hw_old(:) divided by the new
164	!! heartwood ::hw_new(:)
165	!!
166	!! RECENT CHANGE(S) : None
167	!!
168	!! MAIN OUTPUT VARIABLES: ::Biomass of leaves, fruits, fine roots and
169	!! sapwood above (latter for grasses only), ::Update leaf age distribution
170	!! with new leaf age class fraction
171	!!
172	!! REFERENCE(S) :
173	!! - Krinner, G., N. Viovy, N. de Noblet-Ducoudre, J. Ogee, J. Polcher, P.
174	!! Friedlingstein, P. Ciais, S. Sitch and I.C. Prentice (2005), A dynamic global
175	!! vegetation model for studies of the coupled atmosphere-biosphere system,
176	!! Global Biogeochemical Cycles, 19, doi:10.1029/2003GB002199.
177	!! - McNaughton, S. J., M. Oesterheld, D. A. Frank and K. J. Williams (1989),
178	!! Ecosystem-level patterns of primary productivity and herbivory in
179	!! terrestrial habitats, Nature, 341, 142-144, 1989.
180	!! - Sitch, S., C. Huntingford, N. Gedney, P. E. Levy, M. Lomas, S. L. Piao,
181	!! Betts, R., Ciais, P., Cox, P., Friedlingstein, P., Jones, C. D.,
182	!! Prentice, I. C. and F. I. Woodward : Evaluation of the terrestrial carbon
183	!! cycle, future plant geography and climate-carbon cycle feedbacks using 5
184	!! dynamic global vegetation models (dgvms), Global Change Biology, 14(9),
185	!! 2015 –2039, 2008.
186	!!
187	!! FLOWCHART :
188	!! \latexonly
189	!! \includegraphics[scale=0.5]{turnover_flowchart_1.png}
190	!! \includegraphics[scale=0.5]{turnover_flowchart_2.png}
191	!! \endlatexonly
192	!! \n
193	!_ ================================================================================================================================
194
195	SUBROUTINE turn_over (npts, dt, PFTpresent, herbivores, &
196	gpp_week, resp_maint_week, &
197	maxvegstress_lastyear, minvegstress_lastyear, vegstress_week, &
198	vegstress_month, t2m_longterm, t2m_month, t2m_week, &
199	veget_max, gdd_from_growthinit, leaf_age, leaf_frac, age, &
200	turnover, plant_status, turnover_time, &
201	circ_class_biomass, circ_class_n, &
202	when_growthinit, longevity_eff_leaf, longevity_eff_sap, longevity_eff_root, &
203	harvest_pool, harvest_type, harvest_cut, harvest_area, &
204	wstress_month, leaf_age_crit, doy_end_gs)
205
206	!! 0. Variable and parameter declaration
207
208	!! 0.1 Input variables
209
210	INTEGER(i_std), INTENT(in) :: npts !! Domain size - number of grid cells
211	!! (unitless)
212	REAL(r_std), INTENT(in) :: dt !! time step (dt_days)
213	LOGICAL, DIMENSION(:,:), INTENT(in) :: PFTpresent !! PFT exists (true/false)
214	REAL(r_std), DIMENSION(:,:), INTENT(in) :: herbivores !! time constant of probability of a leaf to
215	!! be eaten by a herbivore (days)
216	REAL(r_std), DIMENSION(:,:), INTENT(in) :: maxvegstress_lastyear !! last year's maximum moisture availability
217	!! (0-1, unitless)
218	REAL(r_std), DIMENSION(:,:), INTENT(in) :: minvegstress_lastyear !! last year's minimum moisture availability
219	!! (0-1, unitless)
220	REAL(r_std), DIMENSION(:,:), INTENT(in) :: vegstress_week !! "weekly" moisture availability
221	!! (0-1, unitless)
222	REAL(r_std), DIMENSION(:,:), INTENT(in) :: vegstress_month !! "monthly" moisture availability
223	!! (0-1, unitless)
224	REAL(r_std), DIMENSION(:), INTENT(in) :: t2m_longterm !! "longterm" 2-meter temperatures (K)
225	REAL(r_std), DIMENSION(:), INTENT(in) :: t2m_month !! "monthly" 2-meter temperatures (K)
226	REAL(r_std), DIMENSION(:), INTENT(in) :: t2m_week !! "weekly" 2 meter temperatures (K)
227	REAL(r_std), DIMENSION(:,:), INTENT(in) :: veget_max !! "maximal" coverage fraction of a PFT (LAI
228	!! -> infinity) on ground (unitless)
229	REAL(r_std), DIMENSION(:,:), INTENT(in) :: gpp_week !! PFT gross primary productivity
230	REAL(r_std), DIMENSION(:,:), INTENT(in) :: resp_maint_week !! Weekly maintenance respiration
231	REAL(r_std), DIMENSION(:,:), INTENT(in) :: gdd_from_growthinit !! gdd senescence for crop
232	REAL(r_std), DIMENSION(:,:), INTENT(in) :: longevity_eff_root !! Effective root turnover time that accounts
233	!! waterstress (days)
234	REAL(r_std), DIMENSION(:,:), INTENT(in) :: longevity_eff_sap !! Effective sapwood turnover time that accounts
235	!! waterstress (days)
236	REAL(r_std), DIMENSION(:,:), INTENT(in) :: longevity_eff_leaf !! Effective leaf turnover time that accounts
237	!! waterstress (days)
238	REAL(r_std), DIMENSION(:,:), INTENT(in) :: wstress_month !! Water stress factor, based on hum_rel_daily
239	!! (unitless, 0-1)
240	REAL(r_std), DIMENSION(:,:), INTENT(in) :: leaf_age_crit !! critical leaf age (days)
241
242	!! 0.2 Output variables
243	REAL(r_std),DIMENSION(:,:),INTENT(out) :: doy_end_gs !! growing season end day of year (DOY) for
244	!! deciduous PFTs.
245
246	!! 0.3 Modified variables
247
248	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: leaf_age !! age of the leaves (days)
249	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: leaf_frac !! fraction of leaves in leaf age class
250	!! (0-1, unitless)
251	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: harvest_pool !! The wood and biomass that have been
252	!! havested by humans @tex $(gC)$ @endtex
253	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: harvest_type !! Type of management that resulted
254	!! in the harvest (unitless)
255	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: harvest_cut !! Type of cutting that was used for the harvest
256	!! (unitless)
257	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: harvest_area !! Harvested area (m^{2})
258	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: age !! age (years)
259	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: turnover_time !! turnover_time of grasses (days)
260	REAL(r_std), DIMENSION(:,:,:,:), INTENT(out) :: turnover !! Turnover @tex ($gC m^{-2}$) @endtex
261	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(inout) :: circ_class_biomass !! Biomass of the componets of the model
262	!! tree within a circumference
263	!! class @tex $(gC ind^{-1})$ @endtex
264	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: plant_status !! Growth and phenological status of the plant
265	!! Different stati are defined in constantes
266	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: circ_class_n !! Number of individuals in each circ class
267	!! @tex $(m^{-2})$ @endtex
268	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: when_growthinit !! How many days ago was the beginning of
269	!! the growing season (days)
270
271
272	!! 0.4 Local variables
273
274	LOGICAL, DIMENSION(npts) :: shed_rest !! shed the remaining leaves? (true/false)
275	INTEGER(i_std) :: ivm,iele,ilage,ipts !! Index (unitless)
276	INTEGER(i_std) :: ipar,icir,imbc,ij !! Index (unitless)
277	INTEGER(i_std), SAVE :: turn_count !! Counter
278	!$OMP THREADPRIVATE(turn_count)
279	REAL(r_std), DIMENSION(npts,nvm) :: lai !! leaf area index @tex ($m^2 m^{-2}$)
280	REAL(r_std), DIMENSION(npts,nvm) :: leaf_meanage !! mean age of the leaves (days)
281	REAL(r_std), DIMENSION(npts,ncirc) :: dturnover !! Intermediate variable for turnover ??
282	!! @tex ($gC m^{-2}$) @endtex
283	REAL(r_std), DIMENSION(npts) :: vegstress_crit !! critical moisture availability, function
284	!! of last year's moisture availability
285	!! (0-1, unitless)
286	REAL(r_std), DIMENSION(npts) :: tl !! long term annual mean temperature, (C)
287	REAL(r_std), DIMENSION(npts) :: t_crit !! critical senescence temperature, function
288	!! of long term annual temperature (K)
289	REAL(r_std), DIMENSION(npts) :: sapconv !! Sapwood conversion @tex ($gC m^{-2}$)
290	!! @endtex
291	REAL(r_std), DIMENSION(npts) :: hw_old !! old heartwood mass @tex ($gC m^{-2}$)
292	!! @endtex
293	REAL(r_std), DIMENSION(npts) :: hw_new !! new heartwood mass @tex ($gC m^{-2}$)
294	!! @endtex
295	REAL(r_std), DIMENSION(npts) :: lm_old !! old leaf mass @tex ($gC m^{-2}$) @endtex
296	REAL(r_std), DIMENSION(npts) :: init_biomass !! Biomass before turnover. This variable is used
297	!! in IF-statements to ensure that the same
298	!! initial biomass is used for N and C
299	!! @tex ($gC m^{-2}$) @endtex
300	REAL(r_std), DIMENSION(npts,nleafages) :: delta_lm !! leaf mass change for each age class @tex
301	!! ($gC m^{-2}$) @endtex
302	REAL(r_std), DIMENSION(npts) :: turnover_rate !! turnover rate (unitless)
303	REAL(r_std), DIMENSION(npts,nvm) :: new_turnover_time !! instantaneous turnover time (days)
304	REAL(r_std), DIMENSION(npts,nvm) :: harvest_time !! Prescribed harvest time adjusted for the
305	!! longterm temperature at the pixel (days)
306	REAL(r_std) :: branch_turn !! turnover of branches (how much goes
307	!! to litter each day) (cf article CO2fix)
308	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements):: bm_old !! old root mass
309	REAL(r_std), DIMENSION(npts,nvm) :: sapwood_age !! sapwood age (days)
310	REAL(r_std), DIMENSION(npts) :: sw_old !! old sapwood mass
311	!! (gC/(m**2 of nat/agri ground))
312	REAL(r_std), DIMENSION(npts) :: sw_new !! new sapwood mass
313	!! (gC/(m**2 of nat/agri ground)
314	REAL(r_std), DIMENSION(npts,nvm,nmbcomp,nelements)&
315	:: check_intern !! Contains the components of the internal
316	!! mass balance chech for this routine
317	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
318	REAL(r_std), DIMENSION(npts,nvm,nelements) :: closure_intern !! Check closure of internal mass balance
319	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
320	REAL(r_std), DIMENSION(npts,nvm,nelements) :: pool_start, pool_end !! Start and end pool of this routine
321	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
322	REAL(r_std), DIMENSION(npts,nvm) :: veget_max_begin !! temporary storage of veget_max to check
323	!! area conservation (unitless, 0-1)
324	REAL(r_std), DIMENSION(npts,nvm) :: leaf_turn_ageing !! temporary variable for preparing output variable.
325	!! Leaf turnover due to leaf ageing (excl. senescence)
326	!! @tex $(gC m^{-2} dt^{-1})$ @endtex
327	REAL(r_std), DIMENSION(npts,nvm) :: last_plant_status !! Plant status at the start of this routine
328	REAL(r_std), DIMENSION(npts,nvm) :: doy_isenescent !! date of year for isenescence
329	REAL(r_std) :: recycle
330	REAL(r_std),DIMENSION(npts) :: temp_w
331	INTEGER(i_std), DIMENSION(4) :: parts
332	REAL(r_std), DIMENSION(npts,nvm) :: tmp_xios !! temporary variable for send xios
333	REAL(r_std), DIMENSION(npts,nvm) :: ratio_rm_gpp !! temporary variable for send xios
334	LOGICAL :: flag_ratio_sene
335	LOGICAL :: flag_ratio_sene_grass
336	REAL(r_std) :: threshold_ratio
337	!_ ================================================================================================================================
338
339
340	IF (firstcall_turnover) THEN
341	!! Initialize local printlev
342	printlev_loc=get_printlev('turnover')
343	END IF
344
345	IF (printlev_loc>=2) WRITE(numout,*) 'Entering turnover'
346
347	!! Initialize
348	flag_ratio_sene = .TRUE.
349	flag_ratio_sene_grass = .TRUE.
350	threshold_ratio = 1.0
351
352	ratio_rm_gpp(:,:) = zero
353	WHERE (gpp_week .GT. min_stomate)
354	ratio_rm_gpp(:,:) = resp_maint_week(:,:)/gpp_week(:,:)
355	ENDWHERE
356
357	! Debug
358	IF (printlev_loc>=4) THEN
359	DO ipar = 1,nparts
360	DO iele = 1,nelements
361	DO icir =1, ncirc
362	IF(icir == 1 .AND. iele == 1)&
363	WRITE(numout,*) 'Biomass check 01: ',&
364	circ_class_biomass(test_grid,test_pft,1,ipar,iele) * &
365	circ_class_n(test_grid,test_pft,1)
366	ENDDO
367	ENDDO
368	ENDDO
369	ENDIF
370
371	!! 1. first call - output messages
372
373	! Initialize first call
374	IF ( firstcall_turnover ) THEN
375
376	IF(printlev_loc>=4) THEN
377	WRITE(numout,*) 'turnover:'
378	WRITE(numout,*) ' > minimum mean leaf age for senescence ',&
379	'(days) (::min_leaf_age_for_senescence) : ',&
380	min_leaf_age_for_senescence
381	ENDIF
382
383	turn_count = 0
384	firstcall_turnover = .FALSE.
385
386	ENDIF
387	parts(1)=iroot
388	parts(2)=ileaf
389	parts(3)=ifruit
390	parts(4)=isapabove
391
392	!! 2. Initializations
393
394	!! 2.1 set output to zero
395	turnover(:,:,:,:) = zero
396	new_turnover_time(:,:) = zero
397	harvest_time(:,:) = zero
398	doy_end_gs(:,:) = zero
399	doy_isenescent(:,:) = zero
400	leaf_turn_ageing(:,:) = zero
401
402	! Archive plant_status at the start of this routine
403	last_plant_status(:,:) = plant_status(:,:)
404
405	!! 2.2 Initialize check for mass balance closure
406	! The mass balance is calculated at the end of this routine
407	! in section 10.
408	! turnover is always equal to zero, so maybe that term doesn't need
409	! to be included here.
410	IF (err_act.GT.1) THEN
411
412	pool_start = zero
413
414	DO iele = 1,nelements
415	! Biomass pool + turnover
416	DO ipar = 1,nparts
417	DO icir = 1,ncirc
418	! Initial biomass pool
419	pool_start(:,:,iele) = pool_start(:,:,iele) + &
420	circ_class_biomass(:,:,icir,ipar,iele) * &
421	circ_class_n(:,:,icir) * veget_max(:,:)
422	ENDDO
423	! Add turnover to the initial biomass pool
424	pool_start(:,:,iele) = pool_start(:,:,iele) + &
425	(turnover(:,:,ipar,iele) * veget_max(:,:))
426
427	ENDDO
428	ENDDO
429
430	! The biomass harvest pool is expressed in gC pixel-1 So, it
431	! shouldn't be multiplied by veget_max but it should be divided
432	! by area to obtain gC m-2.
433	DO ivm = 1,nvm
434	DO iele = 1,nelements
435	pool_start(:,ivm,iele) = pool_start(:,ivm,iele) + &
436	SUM(harvest_pool(:,ivm,:,iele),2) / &
437	area(:)
438	ENDDO
439	ENDDO
440
441	!! 1.3 Initialize check for area conservation
442	veget_max_begin(:,:) = veget_max(:,:)
443
444	ENDIF ! err_act.GT.1
445
446	!+++CHECK+++
447	! No branch turnover in CAN but it is not clear
448	! why that would be justified. Could be accounted
449	! for somewhere else (seems unlikely).
450	!!$ ! Compute the turnover of branches
451	!!$ ! (2.5% per year : Orchidee standard, Masera 2003, Lehtonen 2004,
452	!!$ ! DeAngelis 1981)
453	!!$ branch_turn = 0.025/(one_year/dt)*ss_branch_turn
454	!++++++++++++
455
456	! Calculate lai
457	DO ivm = 2,nvm
458	DO ipts=1,npts
459	lai(ipts,ivm) = biomass_to_lai(SUM(circ_class_biomass(ipts,ivm,:,ileaf,icarbon)),ivm)
460	ENDDO
461	ENDDO
462	lai(:,1) = zero
463
464	!! 2.3 Initialize check for surface area conservation
465	! Veget_max is a INTENT(in) variable and can therefore
466	! not be changed during the course of this subroutine
467	! No need to check whether the subroutine preserves the
468	! total surface area of the pixel.
469
470	!! 2.4 Recalculate mean leaf age
471	! Mean leaf age is recalculated from leaf ages of separate leaf age
472	! classes. Leaf age is used as a criterion in several of the
473	! following IF/WHERE loops. The mean leaf age is calculated using
474	! the following equation:
475	! \latexonly
476	! \input{turnover_lma_update_eqn1.tex}
477	! \endlatexonly
478	! \n
479	leaf_meanage(:,:) = zero
480
481	DO ilage = 1, nleafages
482	leaf_meanage(:,:) = leaf_meanage(:,:) + &
483	leaf_age(:,:,ilage) * leaf_frac(:,:,ilage)
484	ENDDO
485
486	! Debug
487	IF (printlev_loc>=4) THEN
488	DO ipar = 1,nparts
489	DO iele = 1,nelements
490	DO icir =1, ncirc
491	IF(icir == 1 .AND. iele == 1)&
492	WRITE(numout,*) 'Biomass/turnover check 02: ',&
493	circ_class_biomass(test_grid,test_pft,1,ipar,iele) * &
494	circ_class_n(test_grid,test_pft,1),&
495	turnover(test_grid,test_pft,ipar,iele)
496	ENDDO
497	ENDDO
498	ENDDO
499	ENDIF
500
501	!! 3. Climatic senescence
502
503	! Three different types of "climatic" leaf senescence,
504	! that do not change the age structure.
505
506	DO ivm = 2,nvm ! Loop over # PFTs
507	!! 3.1 Determine if there is climatic senescence.
508	! The climatic senescence can be of three types: sensitivity to
509	! cold temperatures, to lack of water, or both. If meteorological
510	! conditions are fulfilled, a flag senescence is set to TRUE.
511	! Evergreen species do not experience climatic senescence.
512	SELECT CASE ( senescence_type(ivm) )
513
514	CASE ('crop' )
515
516	!+++CHECK+++
517	! This is the original code. CAN took a different approach (see below)
518	! but that different approach need to be carefuly revisited
519	! Crop senescence is based on a GDD criterium as in crop models.
520	! The problem with this approach is the use of a global gdd_senescence. For
521	! C3 crops this value is 2500 but is in the temperate zone not reached by
522	! the end of summer. As such crops keep their leaves through winter and
523	! only exceed this threshold in the next spring. If the harvest finally
524	! takes place, the conditions are good to replant. The fallow season is
525	! therefore reduced to a single day. An alternative approach is now used.
526	!!$ WHERE ( (SUM(circ_class_biomass(:,ivm,:,ileaf,icarbon),2) .GT. zero ) .AND. &
527	!!$ ( leaf_meanage(:,ivm) .GT. min_leaf_age_for_senescence(ivm) )&
528	!!$ .AND.( gdd_from_growthinit(:,ivm) .GT. gdd_senescence(ivm)))
529	!!$
530	!!$ plant_status(:,ivm) = isenescent
531	!!$ doy_isenescent(:,ivm) = julian_diff
532	!!$
533	!!$ ENDWHERE
534
535
536	! This is pure plumbing under time pressure (says the one who
537	! introduced it). A growing season should
538	! last about 4 months where the longterm temperature is 282 K
539	! (crude observation from N. France) before senescence can set
540	! in unless the growing degrees are reached. In colder and warmer
541	! places the growing season is shortened or lengthened. Because
542	! in warmer places the growing season lasts longer, crop can be
543	! planted that required more light. The very long growing season
544	! in warm places does not account for wet/dry seasons but should
545	! be overruled by the gdd_senescence criterium.
546	harvest_time(:,ivm) = 160 + (t2m_longterm(:) - 282) * &
547	ABS(t2m_longterm(:) - 282)
548
549	WHERE ( plant_status(:,ivm) .EQ. icanopy .AND. &
550	leaf_meanage(:,ivm) .GT. min_leaf_age_for_senescence(ivm) .AND. &
551	( gdd_from_growthinit(:,ivm) .GT. gdd_senescence(ivm) .OR. &
552	when_growthinit(:,ivm) .GT. harvest_time(:,ivm) ) )
553
554	plant_status(:,ivm) = isenescent
555	doy_isenescent(:,ivm) = julian_diff
556
557	ENDWHERE
558	!+++++++++++
559
560	CASE ( 'cold' )
561
562	!! 3.1.2 Summergreen species
563	! Climatic senescence is triggered by sensitivity to cold
564	! temperatures as follows:
565	! If biomass is large enough (i.e. when it is greater than zero),
566	! AND (i.e. when leaf mean age is above a certain PFT-dependent
567	! treshold ::min_leaf_age_for_senescence, which is given in
568	! constants), AND the monthly temperature is low enough (i.e.
569	! when monthly temperature ::t2m_month(:) is below a critical
570	! temperature ::t_crit(:), which is calculated in this module),
571	! AND the temperature tendency is negative (i.e. when weekly
572	! temperatures ::t2m_week(:) are lower than monthly
573	! temperatures ::t2m_month(:))
574	! If these conditions are met, senescence is set to TRUE.
575	!
576	! The critical temperature for senescence is calculated using
577	! the following equation:
578	! \latexonly
579	! \input{turnover_temp_crit_eqn2.tex}
580	! \endlatexonly
581	! \n
582	! Critical temperature for senescence may depend on long term
583	! annual mean temperature
584	tl(:) = t2m_longterm(:) - ZeroCelsius
585	t_crit(:) = ZeroCelsius + senescence_temp(ivm,1) + &
586	tl(:) * senescence_temp(ivm,2) + &
587	tl(:)tl(:) senescence_temp(ivm,3)
588
589	!CYmark: ratio_rm_gpp might have complex temporal changes. We keep
590	!min_leaf_age_for_senescence to determine senescent as a safety net.
591	IF (flag_ratio_sene) THEN
592	WHERE ( ( plant_status(:,ivm) .EQ. ipresenescence .OR. &
593	plant_status(:,ivm) .EQ. icanopy ) .AND. &
594	( leaf_meanage(:,ivm).GT.min_leaf_age_for_senescence(ivm)) .AND. &
595	(ratio_rm_gpp(:,ivm) .GT. threshold_ratio ) )
596
597	plant_status(:,ivm) = isenescent
598	doy_isenescent(:,ivm) = julian_diff
599
600	ENDWHERE
601
602	ELSE
603	WHERE ( ( plant_status(:,ivm) .EQ. ipresenescence .OR. &
604	plant_status(:,ivm) .EQ. icanopy) .AND. &
605	( leaf_meanage(:,ivm).GT.min_leaf_age_for_senescence(ivm)).AND. &
606	( t2m_month(:) .LT. t_crit(:) ) .AND. &
607	( t2m_week(:) .LT. t2m_month(:) ) )
608
609	plant_status(:,ivm) = isenescent
610	doy_isenescent(:,ivm) = julian_diff
611
612	ENDWHERE
613	ENDIF
614
615	! Fail safe option. If something goes wrong with t_crit, phenology
616	! will get screwed up, leaf age is getting too high, vcmax is
617	! decreasing, etc. We added an option that should prevent this from
618	! happening. 1.2 is an arbitrary threshold.
619	WHERE ( ( plant_status(:,ivm).EQ.ipresenescence .OR. &
620	plant_status(:,ivm).EQ.icanopy).AND. &
621	when_growthinit(:,ivm).GT.1.2*leaf_age_crit(:,ivm))
622
623	plant_status(:,ivm) = isenescent
624	doy_isenescent(:,ivm) = julian_diff
625
626	ENDWHERE
627
628	CASE ( 'dry' )
629
630	!! 3.1.3 Raingreen species
631	! Climatic senescence is triggered by sensitivity to lack of
632	! water availability as follows:
633	! If biomass is large enough (i.e. when it is greater than zero),
634	! AND (i.e. when leaf mean age is above a certain PFT-dependent
635	! treshold ::min_leaf_age_for_senescence, which is given in
636	! ::stomate_constants), AND the moisture availability drops below
637	! a critical level (i.e. when weekly moisture availability
638	! ::vegstress_week(:,ivm) is below a critical moisture availability
639	! ::vegstress_crit(:), which is calculated in this module),
640	! If these conditions are met, senescence is set to TRUE.
641	!
642	! The moisture availability critical level is calculated using
643	! the following equation:
644	! \latexonly
645	! \input{turnover_moist_crit_eqn3.tex}
646	! \endlatexonly
647	! \n
648	vegstress_crit(:) = &
649	MIN( MAX( minvegstress_lastyear(:,ivm) + hum_frac(ivm) * &
650	( maxvegstress_lastyear(:,ivm) - minvegstress_lastyear(:,ivm) ), &
651	senescence_hum(ivm) ), &
652	nosenescence_hum(ivm) )
653
654	IF (flag_ratio_sene) THEN
655	WHERE ( ( plant_status(:,ivm) .EQ. ipresenescence .OR. &
656	plant_status(:,ivm) .EQ. icanopy) .AND. &
657	(leaf_meanage(:,ivm).GT.min_leaf_age_for_senescence(ivm)).AND.&
658	( ratio_rm_gpp(:,ivm) .GT. threshold_ratio ) )
659
660	plant_status(:,ivm) = isenescent
661	doy_isenescent(:,ivm) = julian_diff
662
663	ENDWHERE
664	ELSE
665	WHERE ( ( plant_status(:,ivm) .EQ. ipresenescence .OR. &
666	plant_status(:,ivm) .EQ. icanopy) .AND. &
667	(leaf_meanage(:,ivm).GT.min_leaf_age_for_senescence(ivm)).AND.&
668	( vegstress_week(:,ivm) .LT. vegstress_crit(:) ) )
669
670	plant_status(:,ivm) = isenescent
671	doy_isenescent(:,ivm) = julian_diff
672
673	ENDWHERE
674	ENDIF
675
676	! Fail safe option. In some wet years (or years with very low LAI) the
677	! vegstress_week never drops below vegstress_crit at some isolated pixels
678	! If that is the case phenology is getting screwed up, leaf age is
679	! getting too high, vcmax is decreaseing, etc. We added an option that
680	! should prevent this to happen. 1.2 is an arbitrary threshold.
681	WHERE ((plant_status(:,ivm).EQ.ipresenescence .OR. &
682	plant_status(:,ivm).EQ.icanopy ) .AND. &
683	when_growthinit(:,ivm).GT.1.2*leaf_age_crit(:,ivm))
684
685	plant_status(:,ivm) = isenescent
686	doy_isenescent(:,ivm) = julian_diff
687
688	ENDWHERE
689
690	CASE ( 'mixed' )
691
692	!! 3.1.4 Mixed criterion: Climatic senescence is triggered
693	!! by sensitivity to temperature or to lack of water
694	! Climatic senescence is triggered by sensitivity to temperature
695	! or to lack of water availability as follows:
696	! If biomass is large enough (i.e. when it is greater than zero),
697	! AND (i.e. when leaf mean age is above a certain PFT-dependent
698	! treshold ::min_leaf_age_for_senescence, which is given in
699	! ::stomate_constants), AND the moisture availability drops below
700	! a critical level (i.e. when weekly moisture availability
701	! ::vegstress_week(:,ivm) is below a critical moisture availability
702	! ::vegstress_crit(:), calculated in this module), OR the monthly
703	! temperature is low enough (i.e. when monthly temperature
704	! ::t2m_month(:) is below a critical temperature ::t_crit(:),
705	! calculated in this module), AND the temperature tendency is
706	! negative (i.e. when weekly temperatures ::t2m_week(:) are
707	! lower than monthly temperatures ::t2m_month(:)).
708	! If these conditions are met, senescence is set to TRUE.
709	vegstress_crit(:) = MIN( MAX( minvegstress_lastyear(:,ivm) + &
710	hum_frac(ivm) * (maxvegstress_lastyear(:,ivm) - &
711	minvegstress_lastyear(:,ivm) ), senescence_hum(ivm) ), &
712	nosenescence_hum(ivm) )
713	tl(:) = t2m_longterm(:) - ZeroCelsius
714	t_crit(:) = ZeroCelsius + senescence_temp(ivm,1) + &
715	tl(:) * senescence_temp(ivm,2) + &
716	tl(:)tl(:) senescence_temp(ivm,3)
717
718	!CYmark: for the moment there is no tree PFT with a 'mixed' mode of
719	!senescence
720	IF ( is_tree(ivm) ) THEN
721
722	! critical temperature for senescence may depend on long
723	! term annual mean temperature
724	WHERE ( ( plant_status(:,ivm) .EQ. ipresenescence .OR. &
725	plant_status(:,ivm) .EQ. icanopy) .AND. &
726	(leaf_meanage(:,ivm).GT.min_leaf_age_for_senescence(ivm)).AND.&
727	( ( vegstress_week(:,ivm) .LT. vegstress_crit(:) ) .OR. &
728	( ( t2m_month(:) .LT. t_crit(:) ) .AND. &
729	( t2m_week(:) .LT. t2m_month(:) ) ) ) )
730
731	plant_status(:,ivm) = isenescent
732	doy_isenescent(:,ivm) = julian_diff
733
734	ENDWHERE
735
736
737	ELSE
738
739	IF (natural(ivm)) THEN
740	! Grasses
741	! Because the DOFOCO increase of LAI is much slower than in
742	! the trunk there is a period at the start of the growing
743	! season where LAI is less than lai_happy. In that period
744	! ::t2_month is still increasing from its winter values but
745	! could be below t_crit. That condition resulted in many
746	! grasslands along the Atlantic and North Sea coast to go
747	! into senescence in May. If someone wants to put it back in
748	! be carefull with the brackets. ((t2m_month(:) .LT. t_crit(:))
749	! .AND. (lai(:,ivm) .GT. lai_happy(ivm)))
750	!!$ WHERE ( ( (plant_status(:,ivm) .EQ. ipresenescence) .AND. &
751	!!$ (leaf_meanage(:,ivm).GT.min_leaf_age_for_senescence(ivm)).AND.&
752	!!$ (t2m_month(:) .LT. ZeroCelsius) .AND. &
753	!!$ (t2m_week(:) .LT. t2m_month(:)) ) )
754	!!$
755	!!$ ! Shed leaves, roots and fruits at a high rate = senescence
756	!!$ turnover_time(:,ivm,ileaf) = leaffall(ivm)
757	!!$ turnover_time(:,ivm,iroot) = leaffall(ivm)
758	!!$ turnover_time(:,ivm,ifruit) = leaffall(ivm)
759	!!$
760	!!$ ! Slowly turnover the structural carbon this allows to store
761	!!$ ! reserves. This structural pool is essential for the
762	!!$ ! functioning of the labile and reserve pools
763	!!$ turnover_time(:,ivm,isapabove) = longevity_eff_sap(:,ivm)
764	!!$ plant_status(:,ivm) = isenescent
765	!!$
766	!!$ ENDWHERE
767
768	!++++CHECK++++
769	! This is still a patch in the trunck (Tag 3.0). The variable
770	! have been changed since ORCHIDEE-CN-CAn so we tried to follow
771	! the idea of the patch rather than its identic code. The idea is
772	! that, if the leaf age gets too high (i.e., too high of fraction
773	! of leaves are in the oldest leaf age class), senescence is
774	! triggered. In the trunk, this is done with lgrassleafage=TRUE.
775	! for sensesence type "mixed", and it triggers the senescence flag.
776	! In ORCHIDEE-CN-CAN, this is the plant_status variable. In the
777	! trunk, the patch increases turnover_time. In CAN, we also increased
778	! the turnover_time, but in a different way.
779	IF(LNVGRASSPATCH)THEN
780	IF(flag_ratio_sene_grass) THEN
781	WHERE ((plant_status(:,ivm) .EQ. ipresenescence .OR. &
782	plant_status(:,ivm) .EQ. icanopy ) .AND. &
783	(ratio_rm_gpp(:,ivm) .GT. threshold_ratio) )
784
785	! Shed leaves, roots and fruits at a high rate = senescence
786	turnover_time(:,ivm,ileaf) = leaffall(ivm)
787	turnover_time(:,ivm,iroot) = leaffall(ivm)
788	turnover_time(:,ivm,ifruit) = leaffall(ivm)
789
790	! Slowly turnover the structural carbon this allows to store
791	! reserves. This structural pool is essential for the
792	! functioning of the labile and reserve pools
793	turnover_time(:,ivm,isapabove) = longevity_eff_sap(:,ivm)
794	plant_status(:,ivm) = isenescent
795	doy_isenescent(:,ivm) = julian_diff
796	ENDWHERE
797	ELSE
798	! This part is the patch.
799	WHERE ((plant_status(:,ivm) .EQ. ipresenescence .OR. &
800	plant_status(:,ivm) .EQ. icanopy ) .AND. &
801	( leaf_frac(:,ivm,nleafages) .GT. 0.6) .AND. &
802	(when_growthinit(:,ivm).GT. 100.))
803
804	! Shed leaves, roots and fruits at a high rate = senescence
805	turnover_time(:,ivm,ileaf) = leaffall(ivm)
806	turnover_time(:,ivm,iroot) = leaffall(ivm)
807	turnover_time(:,ivm,ifruit) = leaffall(ivm)
808
809	! Slowly turnover the structural carbon this allows to store
810	! reserves. This structural pool is essential for the
811	! functioning of the labile and reserve pools
812	turnover_time(:,ivm,isapabove) = longevity_eff_sap(:,ivm)
813	plant_status(:,ivm) = isenescent
814	doy_isenescent(:,ivm) = julian_diff
815	ENDWHERE
816	ENDIF
817	ELSE
818	! These are the lines from CAN above
819	WHERE ( ( (plant_status(:,ivm) .EQ. ipresenescence .OR. &
820	plant_status(:,ivm) .EQ. icanopy) .AND. &
821	(leaf_meanage(:,ivm).GT.min_leaf_age_for_senescence(ivm)).AND.&
822	(t2m_month(:) .LT. ZeroCelsius) .AND. &
823	(t2m_week(:) .LT. t2m_month(:)) ) )
824
825	! Shed leaves, roots and fruits at a high rate = senescence
826	turnover_time(:,ivm,ileaf) = leaffall(ivm)
827	turnover_time(:,ivm,iroot) = leaffall(ivm)
828	turnover_time(:,ivm,ifruit) = leaffall(ivm)
829
830	! Slowly turnover the structural carbon this allows to store
831	! reserves. This structural pool is essential for the
832	! functioning of the labile and reserve pools
833	turnover_time(:,ivm,isapabove) = longevity_eff_sap(:,ivm)
834	plant_status(:,ivm) = isenescent
835	doy_isenescent(:,ivm) = julian_diff
836
837	ENDWHERE
838	ENDIF ! LNVGRASSPATCH
839
840	ELSE
841
842	! Senescence for croplands
843	WHERE ((plant_status(:,ivm) .EQ. ipresenescence .OR. &
844	plant_status(:,ivm) .EQ. icanopy ) .AND. &
845	( leaf_frac(:,ivm,nleafages) .GT. 0.6) .AND. &
846	(when_growthinit(:,ivm) .GT. 100.))
847
848	! Shed leaves, roots and fruits at a high rate = senescence
849	turnover_time(:,ivm,ileaf) = leaffall(ivm)
850	turnover_time(:,ivm,iroot) = leaffall(ivm)
851	turnover_time(:,ivm,ifruit) = leaffall(ivm)
852
853	! Slowly turnover the structural carbon this allows to store
854	! reserves. This structural pool is essential for the
855	! functioning of the labile and reserve pools
856	turnover_time(:,ivm,isapabove) = longevity_eff_sap(:,ivm)
857	plant_status(:,ivm) = isenescent
858	doy_isenescent(:,ivm) = julian_diff
859
860	ENDWHERE
861
862	END IF ! natural
863
864	ENDIF
865
866	CASE ( 'none' )
867
868	!! 3.1.5 Evergreen species
869	! Evergreen species do not experience climatic senescence
870
871	CASE default
872
873	!! 3.1.6 Other cases
874	! In case no climatic senescence type is recognized.
875	WRITE(numout,*) ' turnover: don''t know how to treat this PFT.'
876	WRITE(numout,*) ' number (::ivm) : ',ivm
877	WRITE(numout,*) ' senescence type (::senescence_type(ivm)) : ',&
878	senescence_type(ivm)
879	CALL ipslerr_p (3,'stomate_turnover',&
880	'turnover: don''t know how to treat this PFT','case 1','')
881
882	END SELECT
883
884	! Debug
885	IF (printlev_loc >= 4 .AND. ivm .EQ. test_pft) THEN
886	WRITE(numout,*) 'phenology, plant_status, ',plant_status(:,ivm)
887	ENDIF
888	!-
889
890	!! 3.2 Drop leaves and roots, plus stems and fruits for grasses
891	IF ( is_tree(ivm) ) THEN
892	!! 3.2.1 Trees in climatic senescence lose their fine roots at
893	!! the same rate as they lose their leaves.
894	! The rate of biomass loss of both fine roots and leaves
895	! is presribed through the equation:
896	! \latexonly
897	! \input{turnover_clim_senes_biomass_eqn4.tex}
898	! \endlatexonly
899	! \n
900	! with ::leaffall(ivm) a PFT-dependent time constant which is
901	! given in ::stomate_constants),
902	! Calculate stand-level turnover for
903	! both carbon and nitrogen (gC tree-1)).
904	! For nitrogen, we first calculate the turnover that
905	! is lost from the respective biomass pools, i.e.,
906	! all the nitrogen is lost from the leaf pool. Then
907	! recycle some of the leaf nitrogen into the labile pool
908	! and finaly correct the turnover pool for the recycled
909	! nitrogen.
910
911	DO iele =1, nelements
912	DO ij = 1, 2 ! iroot and ileaf
913	ipar=parts(ij)
914	IF (ipar .EQ. ileaf .AND. iele .EQ. initrogen) THEN
915	recycle=recycle_leaf(ivm)
916	ELSEIF (ipar .EQ. iroot .AND. iele .EQ. initrogen) THEN
917	recycle=recycle_root(ivm)
918	ELSE
919	recycle=zero
920	ENDIF
921	dturnover(:,:) = zero
922	DO icir=1, ncirc
923	WHERE ( plant_status(:,ivm) .EQ. isenescent )
924	dturnover(:,icir) = dt/ leaffall(ivm)* &
925	circ_class_biomass(:,ivm,icir,ipar,iele)
926	ENDWHERE
927	ENDDO
928	circ_class_biomass(:,ivm,:,ilabile,iele)= &
929	circ_class_biomass(:,ivm,:,ilabile,iele) + &
930	dturnover(:,:)*recycle
931	turnover(:,ivm,ipar,iele) = turnover(:,ivm,ipar,iele)+ &
932	SUM(dturnover(:,:) * circ_class_n(:,ivm,:),2) * &
933	( un - recycle )
934	circ_class_biomass(:,ivm,:,ipar,iele) = &
935	circ_class_biomass(:,ivm,:,ipar,iele) - &
936	dturnover(:,:)
937	ENDDO
938	ENDDO
939
940	ELSE
941	!! 3.2.2.1 crops with 'crop' phenological model
942	IF (senescence_type(ivm) .EQ. 'crop') THEN
943	DO ipts = 1,npts
944	IF (plant_status(ipts,ivm) .EQ. isenescent) THEN
945	! Crops are planted every year. So make sure to remove
946	! everything at the end of senescence which for crops
947	! is actually harvest. Next stomate_prescribe should
948	! prescribe a new crop. If sapwood is not removed, some
949	! is left on site and the model will try to grow a crop.
950	! However, there are not enough reserves left so the
951	! crop will die. This results in one year of normal crop
952	! growth (following the prescribe) and then one year of
953	! almost no growth (the year starting from the too low
954	! reserves). This issue has been fixed.
955	CALL crop_harvest(npts, ipts, ivm, dt, &
956	veget_max, turnover, harvest_pool, harvest_type, &
957	harvest_cut, harvest_area, &
958	circ_class_biomass,circ_class_n,leaf_meanage)
959
960	! Update plant_status. Crops were harvested. No living
961	! biomass was left on-site thus we will have to
962	! prescribe a new vegetation during the next growing
963	! season.
964	plant_status(ipts,ivm) = idead
965
966	ENDIF
967	ENDDO
968
969	!! 3.2.2.2 grass based on 'mixed' phenological model
970	ELSEIF (senescence_type(ivm) .EQ. 'mixed') THEN
971	DO iele = 1,nelements ! carbon and nitrogen (gC m-2)
972	DO ij = 1, 4 ! ileaf, iroot, ifruit, isapabove
973	ipar=parts(ij)
974	IF (ipar .EQ. ileaf .AND. iele .EQ. initrogen) THEN
975	recycle=recycle_leaf(ivm)
976	ELSEIF (ipar .EQ. iroot .AND. iele .EQ. initrogen) THEN
977	recycle=recycle_root(ivm)
978	ELSE
979	recycle=zero
980	ENDIF
981	! First calculate the total turnover. Part of the
982	! turnover will be recycled and stored in the labile
983	! pool. The remaining mass is added to the turnover
984	! pool that will eventually become litter. Once the
985	! turnover is known the biomass pool can be updated.
986	dturnover(:,:) = zero
987	WHERE (plant_status(:,ivm) .EQ. isenescent)
988	dturnover(:,1) = circ_class_biomass(:,ivm,1,ipar,iele)*&
989	dt/turnover_time(:,ivm,ipar)
990	ENDWHERE
991	circ_class_biomass(:,ivm,1,ilabile,iele) = &
992	circ_class_biomass(:,ivm,1,ilabile,iele) + &
993	recycle * dturnover(:,1)
994	turnover(:,ivm,ipar,iele) = turnover(:,ivm,ipar,iele)+ &
995	dturnover(:,1) * circ_class_n(:,ivm,1)* (un-recycle)
996	circ_class_biomass(:,ivm,1,ipar,iele) = &
997	circ_class_biomass(:,ivm,1,ipar,iele) - &
998	dturnover(:,1)
999	! Debug
1000	IF (printlev_loc>4) THEN
1001	IF(ivm == test_pft)THEN
1002	WRITE(numout,*) 'total turnover, ', test_grid,test_pft,ipar,iele
1003	WRITE(numout,*) 'dturnover, ', dturnover(test_grid,1)
1004	WRITE(numout,*) 'recycle, ', recycle
1005	WRITE(numout,*) 'turnover, ', turnover(test_grid,test_pft,ipar,iele)
1006	WRITE(numout,*) 'circ_class_biomass, ', &
1007	circ_class_biomass(test_grid,test_pft,1,ipar,iele)
1008	ENDIF
1009	ENDIF
1010
1011	ENDDO
1012	ENDDO
1013
1014	! Debug
1015	IF (printlev_loc>4) THEN
1016	WRITE(numout,*) 'test - recycle, ', ivm, recycle_root(ivm), &
1017	recycle_leaf(ivm)
1018	ENDIF
1019	!-
1020	ELSEIF (senescence_type(ivm) .EQ. 'none') THEN
1021	!! 3.2.2.3 Grasses modeled as an evergreen biome
1022	turnover(:,ivm,:,:) = zero
1023	ELSE
1024	!! 3.2.2.4 Conceptual problem
1025	WRITE(numout,*) 'ERROR: senenescence type not known - 2'
1026	CALL ipslerr_p(3,'stomate_turnover',&
1027	'turnover: senescence type not known','case 2','')
1028	END IF
1029	ENDIF ! tree/grass
1030	ENDDO ! loop over PFTs
1031
1032	! Debug
1033	IF (printlev_loc>=4) THEN
1034	DO ipar = 1,nparts
1035	DO iele = 1,nelements
1036	DO icir =1, ncirc
1037	IF(icir == 1 .AND. iele == 1)&
1038	WRITE(numout,*) 'Biomass/turnover check 03: ',&
1039	circ_class_biomass(test_grid,test_pft,1,ipar,iele) * &
1040	circ_class_n(test_grid,test_pft,1),&
1041	turnover(test_grid,test_pft,ipar,iele)
1042	ENDDO
1043	ENDDO
1044	ENDDO
1045	ENDIF
1046
1047	!! 4. Leaf fall
1048	! At a certain age, leaves fall off, even if the climate would allow a
1049	! green plant all year round. Even if the meteorological conditions are
1050	! favorable for leaf maintenance, plants, and in particular, evergreen
1051	! trees, have to renew their leaves simply because the old leaves become
1052	! inefficient. Another reason for leaf fall may be waterstress. When
1053	! the plant experiences water stress some of the leaves will die.
1054	! Wstress is calculated from the ratio between stressed and unstressed
1055	! GPP. If the wstress is less than 1, this implies that we have to maintain
1056	! a complete canopy (and therefore have the Ra cost) but that only part
1057	! of the canopy will contribute to GPP. Although this is exactely what
1058	! happens on a warm summer day, what we are trying to account for here is
1059	! the impact of a lasting drought (days to weeks). Adaptation to the
1060	! growing conditions are accounted for through KF (see allocation)
1061	! Roots, fruits (and stems) follow leaves. The decay rate varies with
1062	! leaf age. Note that the plant is not declared senescent in this case
1063	! (wchich is important for allocation: if the plant looses leaves because
1064	! of their age, it can renew them).
1065	! Notice that we do not change the reserve pools here (ilabile and
1066	! icarbres). We assume that the tree will move these pools out of old
1067	! leaves and therefore the carbon will not be lost.
1068	! The leaf turnover rate due to aging of leaves is calculated using
1069	! the following equation:
1070	! \latexonly
1071	! \input{turnover_age_senes_biomass_eqn5.tex}
1072	! \endlatexonly
1073	! \n
1074	DO ivm = 2,nvm
1075
1076	!! 4.1 Calculate critical leaf age
1077	! save old leaf mass
1078	lm_old(:) = SUM(circ_class_biomass(:,ivm,:,ileaf,icarbon)* &
1079	circ_class_n(:,ivm,:),2)
1080
1081	!! initialize leaf mass change in age class
1082	delta_lm(:,:) = zero
1083
1084	!! 4.2 Turnover for different leaf age classes
1085	DO ilage = 1,nleafages
1086
1087	turnover_rate(:) = zero
1088
1089	! Age-driven turnover can only occur when there is a canopy
1090	! and if the leaves exceeds a certain age. The plant_status
1091	! of a crop is set for only one day to isenescent. The crop
1092	! is harvested on that day. For the deciudous trees turn-over during senesce
1093	! is accounted for above. For all these reasons, the turnover
1094	! calculated here was limited to the plant_status = icanopy.
1095	WHERE ( (plant_status(:,ivm) .EQ. icanopy .OR. &
1096	plant_status(:,ivm) .EQ. ipresenescence ).AND. &
1097	leaf_age(:,ivm,ilage) .GT. leaf_age_crit(:,ivm) / deux)
1098
1099	! Not clear where this equation comes from or what it is
1100	! based on.
1101	turnover_rate(:) = &
1102	MIN( 0.99_r_std, dt / ( leaf_age_crit(:,ivm) * &
1103	( leaf_age_crit(:,ivm) / leaf_age(:,ivm,ilage) )**quatre ) )
1104
1105	ENDWHERE
1106
1107	IF ( is_tree(ivm) ) THEN
1108
1109	! Stand level turnover (gC m-2) for leaves
1110	! Water stress is increasing (note that ::wstress is thus
1111	! decreasing) so the LAI will be lowered if the decrease is
1112	! small, this will result in a cosmetic change in LAI.
1113	! If water stress is increasing over a the course of a week,
1114	! this decrease may become important. The correction is
1115	! cummulative, for example, ::biomass(ileaf) is decreased
1116	! by 10% in day one and the remaining fraction of
1117	! ::biomass(ileaf) can be decreased by another 10% the next
1118	! day, etc ... The choice to multiply with ::wstress_month
1119	! is arbitrary, we could have well used ::wstress_week. By
1120	! using the monthly value we hope a smoother decrease
1121	! compared to using ::wstress_week or ::wstress_day.
1122	! Update biomass and turnover. During turnover all carbon
1123	! in leaves, roots and fruits will be lost. This is not the
1124	! case for N. Plants are more careful with teir nitrogen
1125	! part of the nitrogen will be resorbed prior to leaf and
1126	! root turnover. The fraction of nitrogen that is resorbed
1127	! is given by ::recycle_leaf. The resorbed nitrogen is
1128	! moved into the labile pool the rest will be lost in
1129	! turnover. Note that no resorprion happens for fruits so
1130	! all nitrogen in a fruit is lost at turnover.
1131	! It is not clear why roots follow the leaves. A turnover
1132	! rate for roots based on longevity_eff_root can be easily calculated
1133	! See crops and grasses where we made use of longevity_eff_root
1134	! Note that this turnover is a bit arbitrairy there is
1135	! no reason why the same fraction of roots and leaves
1136	! should die. This formulation will disturb the allometric
1137	! allocation and will need to be restored by phenological
1138	! growth in the allocation routine.
1139	!+++CHECK+++
1140	!!$ WHERE ( wstress_month(:,ivm) .GE. un-min_stomate)
1141	!!$ temp_w(:) = turnover_rate(:) / &
1142	!!$ ((min_stomate)**(1/tau_hum_month))
1143	!!$ ELSEWHERE( wstress_month(:,ivm) .LT. un-min_stomate)
1144	!!$ temp_w(:) = turnover_rate(:) / &
1145	!!$ ((un-wstress_month(:,ivm))**(1/tau_hum_month))
1146	!!$ END WHERE
1147	!!$
1148	!!$ ! temp_w can be larger than 1. If so there is a chance that we
1149	!!$ ! don have enough leaves to shed. Make sure this does npt happen
1150	!!$ ! and truncate temp_w to 1.
1151	!!$ WHERE ( temp_w(:) .GE. un)
1152	!!$ ! No longer account for wstress
1153	!!$ temp_w(:) = turnover_rate(:)
1154	!!$ ELSE
1155	!!$ temp_w(:) = temp_w(:) * turnover_rate(:)
1156	!!$ ENDWHERE
1157
1158	! This line replaces all lines above. These lines are
1159	! not present in the trunk. Not clear who added them
1160	! and for what reason
1161	temp_w(:) = turnover_rate(:)
1162	!+++++++++++
1163
1164	! Note that delta_lm is in gC m-2 (dturnover is in gC tree-1)
1165	delta_lm(:,ilage) = - temp_w(:) * leaf_frac(:,ivm,ilage) * &
1166	SUM(circ_class_biomass(:,ivm,:,ileaf,icarbon) * &
1167	circ_class_n(:,ivm,:),2)
1168	! Accumulate the leaf turnover due to ageing in gC m-2
1169	leaf_turn_ageing(:,ivm) = leaf_turn_ageing(:,ivm) + delta_lm(:,ilage)
1170
1171	DO iele=1, nelements
1172	! Turnover (gC tree-1): leaves, roots and fruit
1173	! Sapwood turnover is accounted for separatly later
1174	! in this routine
1175	DO ij =1, 3 !iroot, ileaf, ifruit
1176	ipar=parts(ij)
1177	IF (ipar .EQ. ileaf .AND. iele .EQ. initrogen) THEN
1178	recycle=recycle_leaf(ivm)
1179	ELSEIF (ipar .EQ. iroot .AND. iele .EQ. initrogen) THEN
1180	recycle=recycle_root(ivm)
1181	ELSE
1182	recycle= zero
1183	ENDIF
1184	! Note that all of this code is still within an loop over te
1185	! leaf age classes. Avoid double counting by accounting for
1186	! the fractions of the leaf age classes when calculating root,
1187	! sapwood and fruit turnover as well
1188	dturnover(:,:) = zero
1189	DO icir=1, ncirc
1190	IF(ipar .EQ. iroot) THEN
1191	dturnover(:,icir) = leaf_frac(:,ivm,ilage) * dt / &
1192	longevity_eff_root(:,ivm) * &
1193	circ_class_biomass(:,ivm,icir,ipar,iele)
1194	ELSEIF (ipar .EQ. ileaf) THEN
1195	dturnover(:,icir) = leaf_frac(:,ivm,ilage) * temp_w(:) * &
1196	circ_class_biomass(:,ivm,icir,ipar,iele)
1197	ELSEIF(ipar .EQ. ifruit) THEN
1198	dturnover(:,icir) = leaf_frac(:,ivm,ilage) * dt / &
1199	longevity_fruit(ivm) * &
1200	circ_class_biomass(:,ivm,icir,ipar,iele)
1201	ENDIF
1202	ENDDO
1203	! Nitrogen resorption
1204	circ_class_biomass(:,ivm,:,ilabile,iele) = &
1205	circ_class_biomass(:,ivm,:,ilabile,iele) + &
1206	recycle * dturnover(:,:)
1207	!the stand level turnover (gC m-2):
1208	turnover(:,ivm,ipar,iele) = turnover(:,ivm,ipar,iele) + &
1209	SUM(dturnover(:,:)* circ_class_n(:,ivm,:),2) * (un - recycle)
1210	! Update circ_class_biomass
1211	circ_class_biomass(:,ivm,:,ipar,iele) = &
1212	circ_class_biomass(:,ivm,:,ipar,iele) - &
1213	dturnover(:,:)
1214	ENDDO
1215	ENDDO
1216
1217	ELSEIF ( .NOT. is_tree(ivm) .AND. .NOT. natural(ivm) ) THEN
1218
1219	! For crops the sapwood acts as the structure to supports
1220	! leaves and roots and store the reserves and labile.
1221	! Because the growing season is short for crops, the
1222	! turnover of tissues should be kept to a minimum or the GPP
1223	! and NPP will be too low.
1224	! Note that dturnover was redefined so that it can be used
1225	! for both N and C
1226	!+++CHECK+++
1227	!!$ WHERE ( wstress_month(:,ivm) .GE. un-min_stomate)
1228	!!$ temp_w(:) = turnover_rate(:) / &
1229	!!$ ((min_stomate)**(1/tau_hum_month))
1230	!!$ ELSEWHERE( wstress_month(:,ivm) .LT. un-min_stomate)
1231	!!$ temp_w(:) = turnover_rate(:) / &
1232	!!$ ((un-wstress_month(:,ivm))**(1/tau_hum_month))
1233	!!$ END WHERE
1234	!!$ ! Do we have enough biomass to satisfy the turnover?
1235	!!$ WHERE ( temp_w(:) .GE. un)
1236	!!$ ! No longer account for wstress
1237	!!$ temp_w(:) = turnover_rate(:)
1238	!!$ ELSEWHERE
1239	!!$ temp_w(:) = temp_w(:)*turnover_rate(:)
1240	!!$ ENDWHERE
1241	! Most basic approach
1242	temp_w(:) = turnover_rate(:)
1243	!+++++++++++
1244
1245	! Note that delta_lm is in gC m-2 (dturnover is in gC tree-1)
1246	delta_lm(:,ilage) = - temp_w(:) * leaf_frac(:,ivm,ilage) * &
1247	circ_class_biomass(:,ivm,1,ileaf,icarbon)*&
1248	circ_class_n(:,ivm,1)
1249
1250	! Accumulate the leaf turnover due to ageing
1251	leaf_turn_ageing(:,ivm) = leaf_turn_ageing(:,ivm) + delta_lm(:,ilage)
1252
1253	DO iele=1, nelements
1254	! Turnover (gC plant-1):
1255	DO ij =1,4 ! ileaf, iroot, ifruit, isapabove
1256	ipar=parts(ij)
1257	IF (ipar .EQ. ileaf .AND. iele .EQ. initrogen) THEN
1258	recycle=recycle_leaf(ivm)
1259	ELSEIF (ipar .EQ. iroot .AND. iele .EQ. initrogen) THEN
1260	recycle=recycle_root(ivm)
1261	ELSE
1262	recycle= zero
1263	ENDIF
1264	! Note that all of this code is still within an loop over te
1265	! leaf age classes. Avoid double counting by accounting for
1266	! the fractions of the leaf age classes when calculating root,
1267	! sapwood and fruit turnover as well.
1268	dturnover(:,:) = zero
1269	IF(ipar .EQ. iroot) THEN
1270	dturnover(:,1) = circ_class_biomass(:,ivm,1,ipar,iele)*&
1271	dt / longevity_eff_root(:,ivm) * leaf_frac(:,ivm,ilage)
1272	ELSEIF (ipar .EQ. isapabove) THEN
1273	dturnover(:,1) = circ_class_biomass(:,ivm,1,ipar,iele)*&
1274	dt / longevity_eff_sap(:,ivm) * leaf_frac(:,ivm,ilage)
1275	ELSEIF (ipar .EQ. ileaf) THEN
1276	dturnover(:,1) = circ_class_biomass(:,ivm,1,ipar,iele)*&
1277	temp_w(:) * leaf_frac(:,ivm,ilage)
1278	ENDIF
1279	! Resorb some of the nitrogen
1280	circ_class_biomass(:,ivm,1,ilabile,iele) = &
1281	circ_class_biomass(:,ivm,1,ilabile,iele) + &
1282	recycle * dturnover(:,1)
1283	turnover(:,ivm,ipar,iele) = turnover(:,ivm,ipar,iele)+ &
1284	dturnover(:,1)circ_class_n(:,ivm,1)(un - recycle)
1285	! Update circ_class_biomass
1286	circ_class_biomass(:,ivm,1,ipar,iele) = &
1287	circ_class_biomass(:,ivm,1,ipar,iele) - &
1288	dturnover(:,1)
1289	ENDDO
1290	ENDDO
1291
1292	ELSEIF ( .NOT. is_tree(ivm) .AND. natural(ivm) ) THEN
1293
1294	! In CAN grasses were considered managed ecosystem and the turnover
1295	! went into a harvest pool. Here we follow the previous trunk
1296	! definitions and thus assume grasslands are unamanged. The turnover
1297	! should end up in the litter and the within-season turnover should
1298	! be relatively low 10-20% compared to a manage systsem where the
1299	! turnover could be much higher (100%) due to mowing and grazing.
1300	!+++CHECK+++
1301	!!$ WHERE ( wstress_month(:,ivm) .GE. un-min_stomate)
1302	!!$ temp_w(:) = leaf_frac(:,ivm,ilage) * turnover_rate(:) / &
1303	!!$ ((min_stomate)**(1/tau_hum_month))
1304	!!$ ELSEWHERE( wstress_month(:,ivm) .LT. un-min_stomate)
1305	!!$ temp_w(:) = leaf_frac(:,ivm,ilage) * turnover_rate(:) / &
1306	!!$ ((un-wstress_month(:,ivm))**(1/tau_hum_month))
1307	!!$ END WHERE
1308	!!$ ! Do we have enough biomass to satisfy the turnover?
1309	!!$ WHERE ( temp_w(:) .GE. un)
1310	!!$ ! No longer account for wstress
1311	!!$ temp_w(:) = turnover_rate(:)
1312	!!$ ELSEWHERE
1313	!!$ temp_w(:) = temp_w(:)*turnover_rate(:)
1314	!!$ ENDWHERE
1315	! Most basic approach
1316	temp_w(:) = turnover_rate(:)
1317	!+++++++++++
1318
1319	! Note that delta_lm is in gC m-2 (dturnover is in gC tree-1)
1320	delta_lm(:,ilage) = - temp_w(:) * leaf_frac(:,ivm,ilage) * &
1321	circ_class_biomass(:,ivm,1,ileaf,icarbon)*&
1322	circ_class_n(:,ivm,1)
1323
1324	! Accumulate the leaf turnover due to ageing
1325	leaf_turn_ageing(:,ivm) = leaf_turn_ageing(:,ivm) + delta_lm(:,ilage)
1326
1327	DO iele=1, nelements
1328	! Turnover (gC tree-1):
1329	DO ij =1,4 ! ileaf, iroot, ifruit, isapabove
1330	ipar=parts(ij)
1331	IF (ipar .EQ. ileaf .AND. iele .EQ. initrogen) THEN
1332	recycle=recycle_leaf(ivm)
1333	ELSEIF (ipar .EQ. iroot .AND. iele .EQ. initrogen) THEN
1334	recycle=recycle_root(ivm)
1335	ELSE
1336	recycle= zero
1337	ENDIF
1338	! Note that all of this code is still within an loop over te
1339	! leaf age classes. Avoid double counting by accounting for
1340	! the fractions of the leaf age classes when calculating root,
1341	! sapwood and fruit turnover as well.
1342	dturnover(:,:) = zero
1343	IF(ipar .EQ. iroot) THEN
1344	dturnover(:,1) = circ_class_biomass(:,ivm,1,ipar,iele)*&
1345	dt / longevity_eff_root(:,ivm) * leaf_frac(:,ivm,ilage)
1346	ELSEIF (ipar .EQ.isapabove) THEN
1347	dturnover(:,1) =circ_class_biomass(:,ivm,1,ipar,iele)*&
1348	dt / longevity_eff_sap(:,ivm) * leaf_frac(:,ivm,ilage)
1349	ELSEIF (ipar .EQ.ileaf) THEN
1350	dturnover(:,1) = temp_w(:) * leaf_frac(:,ivm,ilage) * &
1351	circ_class_biomass(:,ivm,1,ipar,iele)
1352	ELSEIF (ipar .EQ.ifruit) THEN
1353	dturnover(:,1) = temp_w(:) * leaf_frac(:,ivm,ilage) * &
1354	circ_class_biomass(:,ivm,1,ipar,iele)
1355	ENDIF
1356	! Nitrogen resorption
1357	circ_class_biomass(:,ivm,1,ilabile,iele) = &
1358	circ_class_biomass(:,ivm,1,ilabile,iele) + &
1359	recycle * dturnover(:,1)
1360	! Stand level turnover (gC m-2)
1361	turnover(:,ivm,ipar,iele) = turnover(:,ivm,ipar,iele)+ &
1362	dturnover(:,1)circ_class_n(:,ivm,1) &
1363	(un - recycle)
1364	! Update circ_class_biomass
1365	circ_class_biomass(:,ivm,1,ipar,iele) = &
1366	circ_class_biomass(:,ivm,1,ipar,iele) - &
1367	dturnover(:,1)
1368
1369	ENDDO
1370	ENDDO
1371
1372	ELSE
1373	WRITE(numout,*) 'ERROR: turnover has not been defined - 3'
1374	CALL ipslerr_p(3,'stomate_turnover',&
1375	'turnover: senescence type not known','case 3','')
1376
1377	ENDIF ! is_tree(ivm)
1378
1379	ENDDO ! # leaf ages
1380
1381	!! 4.3 Recalculate the fraction of leaf biomass in each leaf age class.
1382	! Older leaves will be dropped more than younger leaves and therefore
1383	! the leaf age distribution needs to be recalculated after turnover.
1384	! The fraction of biomass in each leaf class is updated using the
1385	! following equation:
1386	! \latexonly
1387	! \input{turnover_update_LeafAgeDistribution_eqn6.tex}
1388	! \endlatexonly
1389	! \n
1390	!
1391	! new fraction = new leaf mass of that fraction / new total leaf mass
1392	! = (old fraction*old total leaf mass ::lm_old(:) + &
1393	! biomass change of that fraction ::delta_lm(:,ipar))/&
1394	! new total leaf mass ::biomass(:,ivm,ileaf)
1395	DO ilage = 1, nleafages
1396	WHERE (SUM(circ_class_biomass(:,ivm,:,ileaf,icarbon)*circ_class_n(:,ivm,:),2) .GT. min_stomate )
1397	leaf_frac(:,ivm,ilage) = (leaf_frac(:,ivm,ilage)*lm_old(:)+ &
1398	delta_lm(:,ilage)) / SUM(circ_class_biomass(:,ivm,:,ileaf,icarbon)* &
1399	circ_class_n(:,ivm,:),2)
1400	ELSEWHERE
1401	leaf_frac(:,ivm,ilage) = zero
1402	ENDWHERE
1403	ENDDO
1404	ENDDO ! loop over PFTs
1405
1406	! Debug
1407	IF (printlev_loc>=4) THEN
1408	DO ipar = 1,nparts
1409	DO iele = 1,nelements
1410	DO icir =1, ncirc
1411	IF(icir == 1 .AND. iele == 1)&
1412	WRITE(numout,*) 'Biomass/turnover check 04: ',&
1413	circ_class_biomass(test_grid,test_pft,1,ipar,iele) * &
1414	circ_class_n(test_grid,test_pft,1),&
1415	turnover(test_grid,test_pft,ipar,iele)
1416	ENDDO
1417	ENDDO
1418	ENDDO
1419	ENDIF
1420
1421	!! 5. Drop all leaves if there is a very low leaf mass during senescence
1422
1423	! Both for deciduous trees and grasses same conditions are checked:
1424	! If biomass is large enough (i.e. when it is greater than zero),
1425	! AND when senescence is set to true
1426	! AND the leaf biomass drops below a critical minimum biomass level
1427	! (i.e. when it is lower than half the minimum initial LAI
1428	! ::lai_initmin(ivm) divided by the specific leaf area ::slainit(ivm),
1429	! ::lai_initmin(ivm) is set to 0.3 in stomate_data.f90 and slainit is a
1430	! PFT-specific constant, If these conditions are met,
1431	! the flag ::shed_rest(:) is set to TRUE.
1432	!
1433	! After this, the biomass of different carbon pools both for trees and
1434	! grasses is set to zero and the mean leaf age is reset to zero. Finally,
1435	! the leaf fraction and leaf age of the different leaf age classes is set
1436	! to zero.
1437	DO ivm = 2,nvm ! Loop over # PFTs
1438	! Debug
1439	IF (ivm .EQ. test_pft .AND. printlev_loc >= 4) THEN
1440	WRITE(numout,*) 'low leaf mass, '
1441	ENDIF
1442	shed_rest(:) = .FALSE.
1443	!! 5.1 For deciduous trees: leaves, fruits and fine roots are dropped
1444	! For deciduous trees: next to leaves, also fruits and fine roots
1445	! are dropped: fruit ::biomass(:,ivm,ifruit) and fine root
1446	! ::biomass(:,ivm,iroot) carbon pools are set to zero.
1447	IF ( is_tree(ivm) .AND. ( senescence_type(ivm) .NE. 'none' ) ) THEN
1448	! Check whether we shed the remaining leaves. The condition
1449	! depends on biomass(ileaf) so first calculate the sheding
1450	! at the tree level. because biomass(ileaf) is not changed
1451	! at the tree level the same statement can then be used at
1452	! the stand level.
1453	! slainit can be used because a threshold is calculated. The
1454	! leaf mass calculated by lai_initmin/2.)/slainit -see below
1455	! is not used in any other calculations.
1456	init_biomass(:) = SUM(circ_class_biomass(:,ivm,:,ileaf,icarbon)*&
1457	circ_class_n(:,ivm,:),2)
1458
1459	DO iele = 1,nelements
1460	DO ij =1,3 !ileaf, iroot,ifruit
1461	ipar=parts(ij)
1462	IF (ipar .EQ. ileaf .AND. iele .EQ. initrogen)THEN
1463	recycle=recycle_leaf(ivm)
1464	ELSEIF (ipar .EQ. iroot .AND. iele .EQ. initrogen ) THEN
1465	recycle=recycle_root(ivm)
1466	ELSE
1467	recycle=zero
1468	ENDIF
1469	dturnover(:,:) = zero
1470	DO icir=1, ncirc
1471	WHERE ( ( init_biomass(:) .GT. zero ) .AND.&
1472	plant_status(:,ivm) .EQ. isenescent .AND. &
1473	( init_biomass(:) .LT. lai_to_biomass(lai_initmin(ivm),ivm)/2))
1474	dturnover(:,icir) = circ_class_biomass(:,ivm,icir,ipar,iele)
1475	circ_class_biomass(:,ivm,icir,ipar,iele)= zero
1476	! Account for resorption at the tree level (gC tree-1)
1477	! The last day of senescence, the recycled nitrogen is moved into
1478	! the reserve pool rather than the labile pool (which is done all
1479	! the other days during senescence). The reason is that the first
1480	! day of dormancy, there are no leaves left and then the labile
1481	! pool is moved into the reserve pool. Deciduous species showed an
1482	! unrealistic one-day spike in the labile pool. By moving the N directly
1483	! into the reserve pool, this spike is avoided.
1484	circ_class_biomass(:,ivm,icir,icarbres,iele) = &
1485	circ_class_biomass(:,ivm,icir,icarbres,iele)+ &
1486	recycle * dturnover(:,icir)
1487	turnover(:,ivm,ipar,iele) = turnover(:,ivm,ipar,iele)+&
1488	dturnover(:,icir)circ_class_n(:,ivm,icir)(un-recycle)
1489	ENDWHERE
1490	ENDDO
1491	ENDDO
1492	ENDDO
1493
1494	! slainit can be used because a threshold is calculated. The
1495	! leaf mass calculated by lai_initmin/2.)/slainit -see below
1496	! is not used in any other calculations.
1497	WHERE(( init_biomass(:) .GT. zero ) .AND. &
1498	plant_status(:,ivm) .EQ. isenescent .AND. &
1499	( init_biomass(:) .LT. &
1500	lai_to_biomass(lai_initmin(ivm),ivm)/2) )
1501	plant_status(:,ivm) = idormant
1502	leaf_meanage(:,ivm) = zero
1503	shed_rest(:) = .TRUE.
1504	ENDWHERE
1505
1506	ELSEIF (.NOT. is_tree(ivm)) THEN
1507
1508	IF (senescence_type(ivm) .EQ. 'crop') THEN
1509
1510	! Nothing should be done because we now
1511	! have a harvest module for crops.
1512
1513	ELSEIF (senescence_type(ivm) .EQ. 'mixed') THEN
1514
1515	!! 6.2 Drop leaves, roots, fruit and sapwood for grasses
1516	! For grasses all aboveground carbon pools, except the
1517	! carbohydrate reserves are affected: fruit
1518	! ::biomass(:,ivm,ifruit,icarbon), fine root
1519	! ::biomass(:,ivm,iroot,icarbon) and sapwood above
1520	! ::biomass(:,ivm,isapabove,icarbon) carbon pools are set
1521	! to zero.
1522	! Shed the remaining leaves if LAI very low.
1523	init_biomass(:) = circ_class_biomass(:,ivm,1,ileaf,icarbon)*&
1524	circ_class_n(:,ivm,1)
1525	DO iele = 1,nelements
1526	! Grasses should live on after senescence. Do not
1527	! shed the sapwood because then the whole
1528	! system collapse because the reserves are
1529	! calculated based on the sapwood. no sapwood = no
1530	! reserves = no phenology = no growth.
1531	!+++ CHECK +++
1532	! The trunk seems to take part of the sapwood for grasses
1533	! in senescence, but we take none here.
1534	DO ij =1,3 !ileaf, iroot, ifruit
1535	ipar=parts(ij)
1536	IF (ipar .EQ. ileaf .AND. iele .EQ. initrogen) THEN
1537	recycle=recycle_leaf(ivm)
1538	ELSEIF (ipar .EQ. iroot .AND. iele .EQ. initrogen) THEN
1539	recycle=recycle_root(ivm)
1540	ELSE
1541	recycle=zero
1542	ENDIF
1543	dturnover(:,:) = zero
1544	! slainit can be used because a threshold is calculated. The
1545	! leaf mass calculated by lai_initmin/2.)/slainit -see below
1546	! is not used in any other calculations.
1547	WHERE ( ( init_biomass(:) .GT. min_stomate ) .AND.&
1548	plant_status(:,ivm) .EQ. isenescent .AND. &
1549	( init_biomass(:) .LT. &
1550	lai_to_biomass(lai_initmin(ivm),ivm)/2) )
1551	! Account for resorption at the tree level (gC tree-1)
1552	! This should be out of the iele loop else labile is
1553	! calculated twice and the mass balance is no longer preserved.
1554	dturnover(:,1) = circ_class_biomass(:,ivm,1,ipar,iele)
1555	circ_class_biomass(:,ivm,1,ipar,iele)= zero
1556	circ_class_biomass(:,ivm,1,ilabile,iele) = &
1557	circ_class_biomass(:,ivm,1,ilabile,iele)+ &
1558	recycle * dturnover(:,1)
1559	turnover(:,ivm,ipar,iele) = turnover(:,ivm,ipar,iele)+ &
1560	dturnover(:,1)circ_class_n(:,ivm,1)(un-recycle)
1561	ENDWHERE
1562	ENDDO
1563	ENDDO
1564	! slainit can be used because a threshold is calculated. The
1565	! leaf mass calculated by lai_initmin/2.)/slainit -see below
1566	! is not used in any other calculations.
1567	WHERE( ( init_biomass(:) .GT. min_stomate ) .AND.&
1568	plant_status(:,ivm) .EQ. isenescent .AND. &
1569	( init_biomass(:) .LT.&
1570	lai_to_biomass(lai_initmin(ivm),ivm)/2) )
1571	shed_rest(:) = .TRUE.
1572	plant_status(:,ivm) = idormant
1573	leaf_meanage(:,ivm) = zero
1574	ENDWHERE
1575
1576	ELSEIF (senescence_type(ivm) .EQ. 'none') THEN
1577	! Nothing should be done
1578
1579	ELSE
1580	WRITE(numout,*) 'ERROR: senescence type is not known - 4'
1581	CALL ipslerr_p(3,'stomate_turnover',&
1582	'turnover: senescence type not known','case 4','')
1583
1584	ENDIF
1585
1586	! Debug
1587	IF (ivm .EQ. test_pft .AND. printlev_loc >= 4) THEN
1588	WRITE(numout,*) 'phenology, shed_rest, ',shed_rest(:)
1589	ENDIF
1590	!-
1591
1592	ENDIF ! is_tree(ivm)
1593
1594
1595	!! Kill PFTs that have no chance to survive their dormancy
1596	! If the plant goes into dormancy the leaves and roots will be
1597	! shed. If the labile and reserve C are empty, the plant won't be
1598	! able to grow again so it should be killed. In most cases this
1599	! will be dealt with in stomate_mark_to_kill but we had a couple of
1600	! pixels that went to idormant at time t and back to ibudavail at
1601	! t+1. Hence the model never passed mortality and crashed in phenology
1602	! because there were no reserve and/or labile carbon pools to
1603	! grow the initial leaves. If a PFT is killed we will set the
1604	! when_growthinit to zero so it is no longer possible to go to
1605	! ibudavail in phenology.
1606	DO ipts = 1,npts
1607
1608	IF(SUM((circ_class_biomass(ipts,ivm,:,ileaf,icarbon) + &
1609	circ_class_biomass(ipts,ivm,:,ilabile,icarbon) + &
1610	circ_class_biomass(ipts,ivm,:,icarbres,icarbon)) * &
1611	circ_class_n(ipts,ivm,:)).LT.min_stomate .AND. &
1612	plant_status(ipts,ivm).EQ.idormant) THEN
1613
1614	! Make sure the plant won't go into budavail or budbreak
1615	! the next day. This PFT should get killed.
1616	when_growthinit(ipts,ivm) = zero
1617
1618	END IF
1619
1620	END DO
1621
1622	!! 5.3 Reset the leaf age structure
1623	! The leaf fraction and leaf age of the different leaf
1624	! age classes is set to zero.
1625	DO ilage = 1, nleafages
1626	WHERE ( shed_rest(:) )
1627	leaf_age(:,ivm,ilage) = zero
1628	leaf_frac(:,ivm,ilage) = zero
1629	ENDWHERE
1630	ENDDO
1631	ENDDO ! loop over PFTs
1632
1633	! Debug
1634	IF (printlev_loc>=4) THEN
1635	DO ipar = 1,nparts
1636	DO iele = 1,nelements
1637	DO icir =1, ncirc
1638	IF(icir == 1 .AND. iele == 1)&
1639	WRITE(numout,*) 'Biomass/turnover check 05: ',&
1640	circ_class_biomass(test_grid,test_pft,1,ipar,iele) * &
1641	circ_class_n(test_grid,test_pft,1),&
1642	turnover(test_grid,test_pft,ipar,iele)
1643	ENDDO
1644	ENDDO
1645	ENDDO
1646	ENDIF
1647
1648	!! 6. Herbivore activity: elephants, cows, gazelles but no lions.
1649
1650	! Herbivore activity affects the biomass of leaves and fruits as well
1651	! as stalks (only for grasses). Herbivore activity does not modify leaf
1652	! age structure. Herbivores ::herbivores(:,ivm) is the time constant of
1653	! probability of a leaf to be eaten by a herbivore, and is calculated in
1654	! ::stomate_season. following Mc Naughton et al. [1989].
1655	IF ( ok_herbivores ) THEN
1656
1657	! If the herbivore activity is allowed (if ::ok_herbivores is true,
1658	! which is set in run.def), remove the amount of biomass consumed
1659	! by herbivory from the leaf biomass ::biomass(:,ivm,ileaf,icarbon) and
1660	! the fruit biomass ::biomass(:,ivm,ifruit,icarbon). The daily amount
1661	! consumed equals the biomass multiplied by 1 day divided by the time
1662	! constant ::herbivores(:,ivm).
1663	DO ivm = 2,nvm ! Loop over # PFTs
1664	IF ( is_tree(ivm) ) THEN
1665	! Debug
1666	IF (ivm .EQ. test_pft .AND. printlev_loc >= 4) THEN
1667	WRITE(numout,*) 'herbivores, ',herbivores
1668	ENDIF
1669	!-
1670	! For trees: only the leaves and fruit carbon pools are affected
1671	! No time to re-allocate nitrogen during browsing
1672	!+++CHECK+++
1673	! Why is root herbivory not accounted for. Could be difficult
1674	! to parameterize but if leaf browsing is accounted for, there
1675	! is no reason to ignore root browsing
1676	init_biomass(:) = SUM(circ_class_biomass(:,ivm,:,ileaf,icarbon)*&
1677	circ_class_n(:,ivm,:),2)
1678	DO iele=1, nelements
1679	! Turnover (gC tree-1)
1680	DO ij =2,3 ! ileaf and ifruit
1681	ipar=parts(ij)
1682	dturnover(:,:) = zero
1683	DO icir=1, ncirc
1684	WHERE (init_biomass(:) .GT. zero .AND.herbivores(:,ivm) .GT. min_stomate)
1685	! Tree Circ level turnover (gC tree-1):
1686	dturnover(:,icir) = circ_class_biomass(:,ivm,icir,ipar,iele) * &
1687	dt/herbivores(:,ivm)
1688	ENDWHERE
1689	ENDDO
1690	!the stand level turnover (gC m-2):
1691	turnover(:,ivm,ipar,iele) = turnover(:,ivm,ipar,iele)+ &
1692	SUM(dturnover(:,:)*circ_class_n(:,ivm,:),2)
1693	! Update circ_class_biomass
1694	circ_class_biomass(:,ivm,:,ipar,iele) = &
1695	circ_class_biomass(:,ivm,:,ipar,iele)-dturnover(:,:)
1696	ENDDO
1697	ENDDO
1698
1699	ELSEIF (.NOT. is_tree(ivm)) THEN
1700
1701	init_biomass(:) = circ_class_biomass(:,ivm,1,ileaf,icarbon)*&
1702	circ_class_n(:,ivm,1)
1703	DO iele=1, nelements
1704	! Turnover (gC tree-1)
1705	DO ij =2,4 !leaf, ifruit, isapwood
1706	ipar=parts(ij)
1707	dturnover(:,:) = zero
1708	WHERE (init_biomass(:) .GT. zero .AND. herbivores(:,ivm) .GT. min_stomate)
1709	! Tree Circ level turnover (gC tree-1):
1710	dturnover(:,1) = circ_class_biomass(:,ivm,1,ipar,iele)* &
1711	dt/herbivores(:,ivm)
1712	ENDWHERE
1713	!the stand level turnover (gC m-2):
1714	turnover(:,ivm,ipar,iele) = turnover(:,ivm,ipar,iele)+&
1715	SUM(dturnover(:,:)*circ_class_n(:,ivm,:),2)
1716	! Update circ_class_biomass
1717	circ_class_biomass(:,ivm,1,ipar,iele) = &
1718	circ_class_biomass(:,ivm,1,ipar,iele)-dturnover(:,1)
1719	ENDDO
1720	ENDDO
1721
1722	ELSE
1723	WRITE(numout,*) 'ERROR: vegetation type is not known'
1724	CALL ipslerr_p(3,'stomate_turnover',&
1725	'vegetation type is not known','','')
1726
1727	ENDIF ! tree/grass?
1728	ENDDO ! loop over PFT
1729	ENDIF ! end herbivores
1730
1731	! Debug
1732	IF (printlev_loc>=4) THEN
1733	DO ipar = 1,nparts
1734	DO iele = 1,nelements
1735	DO icir =1, ncirc
1736	IF(icir == 1 .AND. iele == 1)&
1737	WRITE(numout,*) 'Biomass/turnover check 06: ',&
1738	circ_class_biomass(test_grid,test_pft,1,ipar,iele) * &
1739	circ_class_n(test_grid,test_pft,1),&
1740	turnover(test_grid,test_pft,ipar,iele)
1741	ENDDO
1742	ENDDO
1743	ENDDO
1744	ENDIF
1745
1746	!! 8. Conversion of sapwood to heartwood
1747
1748	! Conversion of sapwood to heartwood both for aboveground and
1749	! belowground sapwood and heartwood. Following LPJ (Sitch et al., 2003),
1750	! sapwood biomass is converted into heartwood biomass with a time
1751	! constant tau ::longevity_sap(ivm) in years. Note that this biomass conversion
1752	! is not added to "turnover" as the biomass is not lost!
1753	DO ivm = 2,nvm ! Loop over # PFTsi
1754	IF ( is_tree(ivm) ) THEN
1755	! For the recalculation of age in 9.2 (in case the vegetation is
1756	! not dynamic ie. ::ok_dgvm is FALSE), the heartwood above and
1757	! below is stored in ::hw_old(:).
1758	IF ( .NOT. ok_dgvm ) THEN
1759	hw_old(:) = SUM(circ_class_biomass(:,ivm,:,iheartabove,icarbon),2) + &
1760	SUM(circ_class_biomass(:,ivm,:,iheartbelow,icarbon),2)
1761	sw_old(:) = SUM(circ_class_biomass(:,ivm,:,isapabove,icarbon),2) + &
1762	SUM(circ_class_biomass(:,ivm,:,isapbelow,icarbon),2)
1763	ENDIF
1764
1765	!! 8.1 Calculate the rate of sapwood to heartwood conversion
1766	! Calculate the rate of sapwood to heartwood conversion with
1767	! the time constant ::longevity_sap(ivm) and update aboveground and
1768	! belowground sapwood ::biomass(:,ivm,isapabove) and
1769	! ::biomass(:,ivm,isapbelow) and heartwood
1770	! ::biomass(:,ivm,iheartabove) and ::biomass(:,ivm,iheartbelow).
1771	! Tree level above and belowground (gC tree-1)
1772	DO iele = 1,nelements
1773	DO icir = 1,ncirc
1774	circ_class_biomass(:,ivm,icir,iheartabove,iele) = &
1775	circ_class_biomass(:,ivm,icir,iheartabove,iele) + &
1776	circ_class_biomass(:,ivm,icir,isapabove,iele) * dt / &
1777	longevity_eff_sap(:,ivm)
1778	circ_class_biomass(:,ivm, icir,isapabove,iele) = &
1779	circ_class_biomass(:,ivm, icir,isapabove,iele) * &
1780	(un - dt / longevity_eff_sap(:,ivm))
1781	circ_class_biomass(:,ivm,icir,iheartbelow,iele) = &
1782	circ_class_biomass(:,ivm, icir,iheartbelow,iele) + &
1783	circ_class_biomass(:,ivm, icir,isapbelow,iele) * dt / &
1784	longevity_eff_sap(:,ivm)
1785	circ_class_biomass(:,ivm, icir,isapbelow,iele) = &
1786	circ_class_biomass(:,ivm, icir,isapbelow,iele) * &
1787	(un - dt / longevity_eff_sap(:,ivm))
1788	ENDDO
1789	ENDDO
1790	!! 8.2 If the vegetation is not dynamic, the age of the plant
1791	!! is decreased.
1792	! The updated heartwood, the sum of new heartwood above and new
1793	! heartwood below after converting sapwood to heartwood, is saved
1794	! as ::hw_new(:). Creation of new heartwood decreases the age of
1795	! the plant with a factor that is determined by: old heartwood
1796	! ::hw_old(:) divided by the new heartwood ::hw_new(:)
1797	IF ( .NOT. ok_dgvm ) THEN
1798	hw_new(:) = SUM(circ_class_biomass(:,ivm,:,iheartabove,icarbon),2) + &
1799	SUM(circ_class_biomass(:,ivm,:,iheartbelow,icarbon),2)
1800	sw_new(:) = SUM(circ_class_biomass(:,ivm,:,isapabove,icarbon),2) + &
1801	SUM(circ_class_biomass(:,ivm,:,isapbelow,icarbon),2)
1802
1803	WHERE ( hw_new(:) .GT. min_stomate )
1804	age(:,ivm) = age(:,ivm) * hw_old(:)/hw_new(:)
1805	ENDWHERE
1806	ENDIF
1807	ENDIF
1808	ENDDO ! loop over PFTs
1809
1810
1811	! Write to output file
1812	CALL xios_orchidee_send_field("HERBIVORES",herbivores)
1813
1814	! the xios operator for leaf_age and leaf_age_crit are maximum so
1815	! the zero values are correctly accounted for in the history file
1816	CALL xios_orchidee_send_field("LEAF_AGE",leaf_meanage)
1817	CALL xios_orchidee_send_field("LEAF_AGE_CRIT",leaf_age_crit)
1818
1819	! By dividing the LEAF_M_MAX_c by LEAF_TURN_AGEING_c one gets
1820	! the annual turnover of leaves during prior to senescence. Based on
1821	! litter traps this value should be 10 to 20%.
1822	CALL xios_orchidee_send_field("LEAF_TURN_AGEING_c",leaf_turn_ageing)
1823
1824	! End of growing season. For phenology ibudbreak indicates the start
1825	! of the growing season. There is no such single-day setting for senescence
1826	! and all stati lasts several days. Use a change in plant_status
1827	! to identify the single-day at which the leaf season has ended.
1828	DO ivm = 2,nvm
1829	DO ipts = 1,npts
1830	! Deciduous trees and grasses should go through isenescent. By
1831	! using .GE.idormant also idead is included.
1832	IF(natural(ivm) .AND. &
1833	last_plant_status(ipts,ivm) .EQ. isenescent .AND. &
1834	plant_status(ipts,ivm) .GE. idormant) THEN
1835	doy_end_gs(ipts,ivm) = julian_diff
1836	ELSEIF (.NOT. natural(ivm) .AND. &
1837	(last_plant_status(ipts,ivm) .EQ. icanopy .OR. &
1838	last_plant_status(ipts,ivm) .EQ. ipresenescence ) .AND. &
1839	plant_status(ipts,ivm) .GE. idormant) THEN
1840	doy_end_gs(ipts,ivm) = julian_diff
1841	ENDIF
1842	ENDDO
1843	ENDDO ! loop of nvm
1844
1845
1846	CALL xios_orchidee_send_field("DOY_END_GS",doy_end_gs)
1847	CALL xios_orchidee_send_field("DOY_ISENE",doy_isenescent)
1848
1849	CALL histwrite_p (hist_id_stomate, 'LEAF_AGE', itime, &
1850	leaf_meanage, npts*nvm, horipft_index)
1851	CALL histwrite_p (hist_id_stomate, 'HERBIVORES', itime, &
1852	herbivores, npts*nvm, horipft_index)
1853
1854
1855
1856	!! 9. Check numerical consistency of this routine
1857
1858	! Debug
1859	IF (printlev_loc>=4) THEN
1860	DO ipar = 1,nparts
1861	DO iele = 1,nelements
1862	DO icir =1, ncirc
1863	IF(icir == 1 .AND. iele == 1)&
1864	WRITE(numout,*) 'Biomass/turnover check 08: ',&
1865	circ_class_biomass(test_grid,test_pft,1,ipar,iele) * &
1866	circ_class_n(test_grid,test_pft,1),&
1867	turnover(test_grid,test_pft,ipar,iele)
1868	ENDDO
1869	ENDDO
1870	ENDDO
1871	ENDIF
1872
1873	IF (err_act.GT.1) THEN
1874
1875	! 9.1 Check surface area
1876	CALL check_vegetation_area("stomate_turnover", npts, veget_max_begin, &
1877	veget_max,'pft')
1878
1879	! 9.2 Calculate final biomass
1880	pool_end(:,:,:) = zero
1881	DO ipar = 1,nparts
1882	DO iele = 1,nelements
1883	DO icir =1, ncirc
1884	pool_end(:,:,iele) = pool_end(:,:,iele) + &
1885	(circ_class_biomass(:,:,icir,ipar,iele) * &
1886	circ_class_n(:,:,icir)* veget_max(:,:))
1887	ENDDO
1888	pool_end(:,:,iele) = pool_end(:,:,iele) + &
1889	(turnover(:,:,ipar,iele) * veget_max(:,:))
1890
1891	ENDDO
1892	ENDDO
1893
1894	! The biomass harvest pool is expressed in gC pixel-1 So, it
1895	! shouldn't be multiplied by veget_max but it should be divided
1896	! by area to obtain gC m-2.
1897	DO ivm = 1,nvm
1898	DO iele = 1,nelements
1899	pool_end(:,ivm,iele) = pool_end(:,ivm,iele) + &
1900	SUM(harvest_pool(:,ivm,:,iele),2) / &
1901	area(:)
1902	ENDDO
1903	ENDDO
1904
1905	!! 9.3 Calculate components of the mass balance
1906	check_intern(:,:,iatm2land,:) = zero
1907	check_intern(:,:,iland2atm,:) = -un * zero
1908	check_intern(:,:,ilat2out,:) = zero
1909	check_intern(:,:,ilat2in,:) = -un * zero
1910	check_intern(:,:,ipoolchange,:) = &
1911	un * (pool_end(:,:,:) - pool_start(:,:,:))
1912	closure_intern(:,:,:) = zero
1913	DO imbc = 1,nmbcomp
1914	DO iele=1,nelements
1915	! Debug
1916	IF (printlev_loc>=4) THEN
1917	WRITE(numout,*) 'check_intern, ivm, imbc, iele, ', imbc, &
1918	iele, check_intern(:,test_pft,imbc,iele)
1919	ENDIF
1920	!-
1921	closure_intern(:,:,iele) = closure_intern(:,:,iele) + &
1922	check_intern(:,:,imbc,iele)
1923	ENDDO
1924	ENDDO
1925
1926	! 9.4 Check mass balance closure
1927	CALL check_mass_balance("stomate_turnover", closure_intern, npts, &
1928	pool_end, pool_start, veget_max,'pft')
1929
1930	ENDIF ! err_act.GT.1
1931
1932	IF(printlev>=3) WRITE(numout,*) 'Leaving turnover'
1933
1934	END SUBROUTINE turn_over
1935
1936
1937	!================================================================================================================================
1938	!! SUBROUTINE : drought_mortality
1939	!!
1940	!>\BRIEF Calculate turnover of sapwood to heartwood induced by drought.
1941	!!
1942	!! DESCRIPTION : This subroutine determines the turnover of sapwood into
1943	!! heartwood induced by drought.
1944	!! RECENT CHANGE(S) : None.
1945	!!
1946	!! MAIN OUTPUT VARIABLES: ::circ_class_biomass
1947	!!
1948	!! REFERENCE(S) :
1949	!!
1950	!! FLOWCHART :
1951	!_
1952	!================================================================================================================================
1953
1954	SUBROUTINE drought_mortality (npts, kill_vessels, vessel_mortality_daily, &
1955	veget_max, biomass_init_drought, bm_to_litter, circ_class_biomass, &
1956	circ_class_n)
1957
1958	!! 0. Variable and parameter declaration
1959
1960	!! 0.1 Input variables
1961	INTEGER(i_std), INTENT(in) :: npts !! Domain size (number of grid cells).
1962	LOGICAL, DIMENSION(:,:), INTENT(in) :: kill_vessels !! Flag to kill vessels at the end of the day following embolism.
1963	REAL(r_std), DIMENSION(:,:), INTENT(in) :: vessel_mortality_daily !! Proportion of daily vessel mortality due to cavitation in the xylem. Equal to
1964	!! a fraction of vl_diff.
1965	REAL(r_std), DIMENSION(:,:), INTENT(in) :: veget_max !! 'Maximal' coverage fraction of a PFT (LAI -> Infinity) on ground (unitless).
1966
1967	!! 0.2 Output variables
1968
1969	!! 0.3 Modified variables
1970	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(inout) :: biomass_init_drought !! Biomass of heartwood or sapwood before onset of drought. Used to compute
1971	!! turnover on same reference biomass in stomate_turnover.f90. Should be the same
1972	!! along entire drought episode.
1973	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: bm_to_litter !! Background mortality of biomass (not senescence-driven).
1974
1975	!! 0.3 Modified variables
1976	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(inout) :: circ_class_biomass !! Biomass of the components of the model tree within a circumference class
1977	!! (gC/ind.).
1978	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: circ_class_n !! Number of individuals in each circ class.
1979
1980	!! 0.4 Local variables
1981	INTEGER(i_std) :: ivm, iele, ilage, ipts !! Index (unitless).
1982	INTEGER(i_std) :: ipar, icir, imbc, ij !! Index (unitless).
1983	REAL(r_std) :: total_sap_biomass !! Total sapwood biomass of model tree.
1984	REAL(r_std), DIMENSION(nparts) :: vessel_turnover !! Turnover for sapwood of model tree.
1985	REAL(r_std), DIMENSION(npts,nvm,ncirc) :: total_cc_biomass !! Total wood biomass of model tree.
1986	REAL(r_std), DIMENSION(npts) :: hw_old !! Old heartwood mass (gC/m2).
1987	REAL(r_std), DIMENSION(npts) :: hw_new !! New heartwood mass (gC/m2).
1988	REAL(r_std), DIMENSION(npts) :: sw_new !! New sapwood mass (gC/m2 of nat./agri. ground).
1989	REAL(r_std), DIMENSION(npts,nvm,nmbcomp,nelements) &
1990	:: check_intern !! Contains the components of the internal mass balance check for this routine
1991	!!(gC/pix./dt).
1992	REAL(r_std), DIMENSION(npts,nvm,nelements) :: closure_intern !! Check closure of internal mass balance (gC/pix./dt).
1993	REAL(r_std), DIMENSION(npts,nvm,nelements) :: pool_start, pool_end !! Start and end pool of this routine gC/pix./dt).
1994	REAL(r_std), DIMENSION(npts,nvm) :: veget_max_begin !! Temporary storage of veget_max to check area vegetation (unitless, [0-1]).
1995
1996	!================================================================================================================================
1997
1998	!! 1.1. Initialize check for mass balance closure
1999
2000	! The mass balance is calculated at the end of this routine
2001	! in section 3.
2002	IF(err_act .GT. 1) THEN
2003	pool_start = zero
2004	DO iele = 1,nelements
2005	! Biomass pool + bm_to_litter.
2006	DO ipar = 1,nparts
2007	DO icir = 1,ncirc
2008	! Initial biomass pool.
2009	pool_start(:,:,iele) = pool_start(:,:,iele) + &
2010	circ_class_biomass(:,:,icir,ipar,iele) * &
2011	circ_class_n(:,:,icir) * veget_max(:,:)
2012	ENDDO
2013	! Add turnover to the initial biomass pool.
2014	pool_start(:,:,iele) = pool_start(:,:,iele) + &
2015	(bm_to_litter(:,:,ipar,iele) * veget_max(:,:))
2016	ENDDO
2017	ENDDO
2018	!! 1.2. Initialize check for area conservation
2019	veget_max_begin(:,:) = veget_max(:,:)
2020
2021	ENDIF ! IF(err_act.GT.1)
2022
2023	! 2. Calculate the rate of sapwood to heartwood conversion with
2024	! the rate of daily mortality ::vessel_mortality_daily(:,ivm).
2025	DO ivm = 2,nvm
2026
2027	IF(is_tree(ivm)) THEN
2028
2029	! Calculation of effect of embolism on heartwood and sapwood
2030	! biomass. Embolism turns sapwood into heartwood. Every time,
2031	! fraction of sapwood biomass is subtracted from sapwood biomass
2032	! and added to heartwood biomass. Distinction between aboveground
2033	! and belowground biomass.
2034	DO ipts = 1,npts
2035	DO iele = 1,nelements
2036	DO icir = 1,ncirc
2037
2038	! Total sapwood in gC tree-1
2039	total_sap_biomass = &
2040	circ_class_biomass(ipts,ivm,icir,isapabove,iele) + &
2041	circ_class_biomass(ipts,ivm,icir,isapbelow,iele)
2042
2043	! Part for drought: If the flag ::kill_vessels(:,ivm) is TRUE,
2044	! and there is sapwood remaining, which means the tree is still
2045	! alive, we compute the turnover induced by drought.
2046	IF(kill_vessels(ipts,ivm) .AND. &
2047	total_sap_biomass .GT. min_stomate) THEN
2048
2049	! First, compute the mass of sapwood that will be
2050	! subtracted from the model tree. Note that vessel_mortality_daily
2051	! contains the increase in mortality (thus not the entire mortality)
2052	! It is multiplied with the biomass at the start of the drought
2053	! to avoid double counting.
2054	vessel_turnover(:) = biomass_init_drought(ipts,ivm,icir,:,iele) * &
2055	vessel_mortality_daily(ipts,ivm)
2056
2057	! Compare ::vessel_turnover(ipar) to the current sapwood
2058	! biomass of the model tree to make sure we do not end up
2059	! with a negative sapwood biomass, which is not realistic. If
2060	! ::vessel_turnover(ipar) is greater than current sapwood
2061	! biomass, we attribute to ::vessel_turnover(ipar) the value
2062	! of the current sapwood biomass, so we will get a null
2063	! sapwood biomass at the end of the turnover calculation.
2064	! Aboveground wood biomass.
2065	IF(vessel_turnover(isapabove) .GE. &
2066	circ_class_biomass(ipts,ivm,icir,isapabove,iele)) THEN
2067	vessel_turnover(isapabove) = &
2068	circ_class_biomass(ipts,ivm,icir,isapabove,iele)
2069	ENDIF
2070
2071	! Use biomass_init_drought in the calculation of
2072	! turnover induced by drought so mortality is not
2073	! overestimated, since ::vessel_mortality_daily(:,ivm) is
2074	! calculated as a proportion of sapwood biomass.
2075	! Biomass of aboveground heartwood. We add dead sapwood to
2076	! heartwood biomass.
2077	circ_class_biomass(ipts,ivm,icir,iheartabove,iele) = &
2078	circ_class_biomass(ipts,ivm,icir,iheartabove,iele) + &
2079	vessel_turnover(isapabove)
2080
2081	! Biomass of aboveground sapwood. We subtract dead sapwood to
2082	! sapwood biomass.
2083	circ_class_biomass(ipts,ivm, icir,isapabove,iele) = &
2084	circ_class_biomass(ipts,ivm,icir,isapabove,iele) - &
2085	vessel_turnover(isapabove)
2086
2087	! Belowground wood biomass.
2088	IF(vessel_turnover(isapbelow) .GE. &
2089	circ_class_biomass(ipts,ivm,icir,isapbelow,iele)) THEN
2090	vessel_turnover(isapbelow) = &
2091	circ_class_biomass(ipts,ivm,icir,isapbelow,iele)
2092	ENDIF
2093
2094	! Biomass of belowground heartwood. We add dead sapwood to
2095	! heartwood biomass.
2096	circ_class_biomass(ipts,ivm,icir,iheartbelow,iele) = &
2097	circ_class_biomass(ipts,ivm, icir,iheartbelow,iele) + &
2098	vessel_turnover(isapbelow)
2099
2100	! Biomass of belowground sapwood. We subtract dead sapwood to
2101	! sapwood biomass.
2102	circ_class_biomass(ipts,ivm, icir,isapbelow,iele) = &
2103	circ_class_biomass(ipts,ivm,icir,isapbelow,iele) - &
2104	vessel_turnover(isapbelow)
2105
2106	ENDIF ! IF(kill_vessels)
2107
2108	! Update total sapwood biomass of model tree after vessel
2109	! moratlity has been dealt with
2110	total_sap_biomass = &
2111	circ_class_biomass(ipts,ivm,icir,isapabove,iele) + &
2112	circ_class_biomass(ipts,ivm,icir,isapbelow,iele)
2113
2114	! Where sapwood biomass becomes zero, model must process
2115	! tree death. The model is still in a DO-loop over iele
2116	! hence the code below will kill the whole tree is either
2117	! icarbon or initrogen is empty.
2118	IF(total_sap_biomass .LE. min_stomate) THEN
2119
2120	! Current biomass of model tree is moved into litter.
2121	! Carbon biomass. We are still in a do-loop over iele
2122	! so we cannot have another do-loop over iele. It was
2123	! written explicitly.
2124	bm_to_litter(ipts,ivm,:,icarbon) = &
2125	bm_to_litter(ipts,ivm,:,icarbon) + &
2126	circ_class_biomass(ipts,ivm,icir,:,icarbon) * &
2127	circ_class_n(ipts,ivm,icir)
2128
2129	! Nitrogen biomass.
2130	bm_to_litter(ipts,ivm,:,initrogen) = &
2131	bm_to_litter(ipts,ivm,:,initrogen) + &
2132	circ_class_biomass(ipts,ivm,icir,:,initrogen) * &
2133	circ_class_n(ipts,ivm,icir)
2134
2135	! To get a closed mass balance, we set biomass and number of
2136	! individuals of dead circumference class to zero since the
2137	! biomass has been moved into litter. We also set the initial
2138	! biomass to zero so that the tree does not revive.
2139	circ_class_biomass(ipts,ivm,icir,:,:) = zero
2140	circ_class_n(ipts,ivm,icir) = zero
2141	biomass_init_drought(ipts,ivm,icir,:,:) = zero
2142
2143	ENDIF ! IF(total_sap_biomass .LE. min_stomate)
2144
2145	ENDDO ! DO icir = 1,ncirc
2146
2147	ENDDO ! DO iele = 1,nelements
2148
2149	ENDDO ! DO ipts = 1,kjpindex
2150
2151	ENDIF ! is_tree
2152
2153	ENDDO ! DO ivm = 2,nvm
2154
2155
2156	!! 3. Check numerical consistency of this routine
2157	IF (err_act.GT.1) THEN
2158
2159	! 3.1. Check surface area
2160	CALL check_vegetation_area("drought_mortality", npts, veget_max_begin, &
2161	veget_max,'pft')
2162
2163	! 3.2. Calculate final biomass
2164	pool_end(:,:,:) = zero
2165	DO ipar = 1,nparts
2166	DO iele = 1,nelements
2167	DO icir =1, ncirc
2168	pool_end(:,:,iele) = pool_end(:,:,iele) + &
2169	(circ_class_biomass(:,:,icir,ipar,iele) * &
2170	circ_class_n(:,:,icir)* veget_max(:,:))
2171	ENDDO
2172	! Added ::bm_to_litter(:,ivm,ipar,iele) here.
2173	pool_end(:,:,iele) = pool_end(:,:,iele) + &
2174	bm_to_litter(:,:,ipar,iele) * veget_max(:,:)
2175	ENDDO
2176	ENDDO
2177
2178	! 3.3. Calculate components of the mass balance
2179	check_intern(:,:,iatm2land,:) = zero
2180	check_intern(:,:,iland2atm,:) = -un * zero
2181	check_intern(:,:,ilat2out,:) = zero
2182	check_intern(:,:,ilat2in,:) = -un * zero
2183	check_intern(:,:,ipoolchange,:) = &
2184	un * (pool_end(:,:,:) - pool_start(:,:,:))
2185	closure_intern(:,:,:) = zero
2186	DO imbc = 1,nmbcomp
2187	DO iele=1,nelements
2188	! Debug
2189	IF (printlev_loc>=4) THEN
2190	WRITE(numout,*) 'check_intern, ivm, imbc, iele, ', imbc, &
2191	iele, check_intern(:,test_pft,imbc,iele)
2192	ENDIF
2193	!-
2194	closure_intern(:,:,iele) = closure_intern(:,:,iele) + &
2195	check_intern(:,:,imbc,iele)
2196	ENDDO
2197	ENDDO
2198
2199	!! 3.4. Check mass balance closure
2200	CALL check_mass_balance("drought_mortality", closure_intern, npts, &
2201	pool_end, pool_start, veget_max,'pft')
2202
2203	ENDIF ! IF(err_act .GT. 1)
2204
2205	IF(printlev>=3) WRITE(numout,*) 'Leaving drought_mortality.'
2206
2207	END SUBROUTINE drought_mortality
2208
2209	END MODULE stomate_turnover

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

Download in other formats: