Context Navigation

source: branches/publications/ORCHIDEE_CAN_r3069/src_stomate/stomate_kill.f90 @ 7475

Last change on this file since 7475 was 3069, checked in by sebastiaan.luyssaert, 9 years ago
PROD: tested gloabally for the zoomed European grid for 20 years. Do not use the previous commit as there were problems with svn it does not contain the described changes. This commit contains the bug fixes to species change and added functionality to change forest management following mortality or a harvest.
File size: 84.5 KB

Line
1	! =================================================================================================================================
2	! MODULE : stomate_kill
3	!
4	! CONTACT : orchidee-help _at_ ipsl.jussieu.fr
5	!
6	! LICENCE : IPSL (2006)
7	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
8	!
9	!>\BRIEF Simulate mortality of individuals and update biomass, litter and
10	!! stand density of the PFT
11	!!
12	!!\n DESCRIPTION : Simulate mortality of individuals and update biomass, litter and
13	!! stand density of the PFT. This module replaces lpj_gap/kill. Biomass actually
14	!! dies here and is moved to the litter pools. This does not work with the DGVM
15	!! at the moment.
16	!!
17	!! RECENT CHANGE(S): None
18	!!
19	!! REFERENCE(S) :
20	!! - Sitch, S., B. Smith, et al. (2003), Evaluation of ecosystem dynamics,
21	!! plant geography and terrestrial carbon cycling in the LPJ dynamic
22	!! global vegetation model, Global Change Biology, 9, 161-185.\n
23	!! - Waring, R. H. (1983). "Estimating forest growth and efficiency in relation
24	!! to canopy leaf area." Advances in Ecological Research 13: 327-354.\n
25	!!
26	!! SVN :
27	!! $HeadURL: svn://forge.ipsl.jussieu.fr/orchidee/branches/ORCHIDEE-DOFOCO/ORCHIDEE/src_stomate/stomate_kill.f90 $
28	!! $Date: 2013-09-27 09:59:24 +0200 (Fri, 27 Sep 2013) $
29	!! $Revision: 1485 $
30	!! \n
31	!_ ================================================================================================================================
32
33	MODULE stomate_kill
34
35	! modules used:
36
37	USE stomate_data
38	USE pft_parameters
39	USE ioipsl_para
40	USE function_library, ONLY : distribute_mortality_biomass,&
41	nmax, wood_to_dia, check_area
42	USE constantes
43	USE stomate_prescribe
44	USE sapiens_lcchange, ONLY : merge_biomass_pfts
45
46	IMPLICIT NONE
47
48	! private & public routines
49
50	PRIVATE
51	PUBLIC gap_prognostic, gap_clear, gap_clean
52
53	! Variable declaration
54
55	LOGICAL, SAVE :: firstcall = .TRUE. !! first call flag
56	!$OMP THREADPRIVATE(firstcall)
57
58	CONTAINS
59
60
61	!! ================================================================================================================================
62	!! SUBROUTINE : gap_clear
63	!!
64	!>\BRIEF Set the firstcall flag back to .TRUE. to prepare for the next simulation.
65	!_ ================================================================================================================================
66
67	SUBROUTINE gap_clear
68	firstcall = .TRUE.
69	END SUBROUTINE gap_clear
70
71
72	!! ================================================================================================================================
73	!! SUBROUTINE : gap_prognostic
74	!!
75	!>\BRIEF Transfer dead biomass to litter and update stand density for trees
76	!! which die of natural causes.
77	!!
78	!! DESCRIPTION : This routine has a simple purpose: kill individuals by natural causes.
79	!! The variable circ_class_kill indicates how many individuals need to be killed by the various
80	!! processes in the code, and is calculated in other modules. gap_prognostic takes
81	!! this variable and decides on the correct order of which to actually kill the
82	!! trees, making sure that we don't double count.
83	!!
84	!! RECENT CHANGE(S):
85	!!
86	!! MAIN OUTPUT VARIABLE(S): ::circ_class_biomass; biomass, ::ind density of individuals,
87	!! ::mortality mortality (fraction of
88	!! trees that is dying per time step)
89	!!
90	!! REFERENCE(S) :
91	!! - Sitch, S., B. Smith, et al. (2003), Evaluation of ecosystem dynamics,
92	!! plant geography and terrestrial carbon cycling in the LPJ dynamic
93	!! global vegetation model, Global Change Biology, 9, 161-185.
94	!! - Waring, R. H. (1983). "Estimating forest growth and efficiency in relation
95	!! to canopy leaf area." Advances in Ecological Research 13: 327-354.
96	!!
97	!! FLOWCHART : None
98	!!\n
99	!_ ================================================================================================================================
100
101	SUBROUTINE gap_prognostic (npts, biomass, ind, bm_to_litter, &
102	circ_class_biomass, circ_class_kill, circ_class_n, &
103	veget_max)
104
105	!! 0. Variable and parameter declaration
106
107	!! 0.1 Input variables
108	INTEGER(i_std), INTENT(in) :: npts !! Domain size (-)
109	REAL(r_std),DIMENSION(npts,nvm),INTENT(in) :: veget_max !! Maximum fraction of vegetation type
110
111	!! 0.2 Output variables
112
113	!! 0.3 Modified variables
114	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), &
115	INTENT(inout) :: biomass !! Biomass @tex $(gC m^{-2}) $@endtex
116	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: ind !! Stand level number of individuals
117	!! @tex $(m^{-2})$ @endtex
118	REAL(r_std), DIMENSION(npts,nvm,ncirc,nparts,nelements), &
119	INTENT(inout) :: circ_class_biomass !! Biomass of the components of the model
120	!! tree within a circumference
121	!! class @tex $(gC ind^{-1})$ @endtex
122	REAL(r_std), DIMENSION(npts,nvm,ncirc), INTENT(inout) :: circ_class_n !! Number of individuals in each circ class
123	!! @tex $(m^{-2})$ @endtex
124	REAL(r_std), DIMENSION(npts,nvm,ncirc,nfm_types,ncut_times),&
125	INTENT(inout) :: circ_class_kill !! Number of trees within a circ that needs
126	!! to be killed @tex $(ind m^{-2})$ @endtex
127	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements), &
128	INTENT(inout) :: bm_to_litter !! Biomass transfer to litter
129	!! @tex $(gC m^{-2})$ @endtex
130
131	!! 0.4 Local variables
132	INTEGER(i_std) :: ipts, ivm, ipar !! Indices
133	INTEGER(i_std) :: iele, icir, imbc !! Indices
134	INTEGER(i_std) :: ifm,icut !! Indices
135	INTEGER,DIMENSION(nvm) :: nkilled !! the number of grid points at which a given
136	!! given PFT's biomass has been reduced to zero
137	REAL(r_std), DIMENSION(npts,nvm,nmbcomp,nelements) :: check_intern !! Contains the components of the internal
138	!! mass balance chech for this routine
139	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
140	REAL(r_std), DIMENSION(npts,nvm,nelements) :: closure_intern !! Check closure of internal mass balance
141	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
142	REAL(r_std), DIMENSION(npts,nvm,nelements) :: pool_start !! Start pool of this routine
143	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
144	REAL(r_std), DIMENSION(npts,nvm,nelements) :: pool_end !! End pool of this routine
145	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
146
147	!_ ================================================================================================================================
148
149	!! 1. Set firstcall flag
150
151	IF ( firstcall ) THEN
152
153	firstcall = .FALSE.
154
155	ENDIF
156
157	IF (bavard.GE.2) WRITE(numout,*) 'Entering gap. Use constant mortality = ', &
158	control%ok_constant_mortality
159
160	!! 2. Initialize variables
161	nkilled(:)=0
162
163	! 2.1 Initialize check for mass balance closure
164	! The mass balance is calculated at the end of this routine
165	! in section .
166	pool_start = zero
167	DO ipar = 1,nparts
168	DO iele = 1,nelements
169	! Initial litter and biomass
170	pool_start(:,:,iele) = pool_start(:,:,iele) + &
171	(biomass(:,:,ipar,iele) * veget_max(:,:)) + &
172	(bm_to_litter(:,:,ipar,iele) * veget_max(:,:))
173	ENDDO
174	ENDDO
175
176	!! 2.2 Initialize check for surface area conservation
177	! Veget_max is a INTENT(in) variable and can therefore
178	! not be changed during the course of this subroutine
179	! No need to check whether the subroutine preserves the
180	! total surface area of the pixel.
181
182
183	!! 3. Kill plants
184
185	DO ipts = 1,npts ! loop over land_points
186
187	DO ivm = 2,nvm ! loop over #PFT
188
189	IF(veget_max(ipts,ivm) == zero)THEN
190	! this vegetation type is not present, so no reason to do the
191	! calculation
192	CYCLE
193	ENDIF
194
195	!! First, all the plants that are supposed to be killed by fire are killed.
196	!! We do not yet know how to handle this, so this will have to be taken
197	!! into account when SPITFIRE is coupled.
198	DO iele = 1,nelements
199
200	!! All of the natural death is grouped
201	!! together in one pool since all of the biomass will be left on
202	!! site (moved to the litter pools).
203	DO ifm=1,nfm_types
204
205	DO icut=1,ncut_times
206
207	! Since we gave circ_class_kill the same dimensions as the
208	! harvest pool, we need to explicitly say which combinations
209	! of ifm and icut are killed in which ways. Currently, I am
210	! putting all natural death into ifm_none, even if it happens
211	! in another management type (for example, self-thinning will
212	! happen with ifm_thin, but it's really a natural death).
213	IF(ifm .NE. ifm_none)CYCLE
214	IF((icut .NE. icut_thin) .AND. (icut .NE. icut_clear))CYCLE
215
216	! Killing is done a little differently for trees and non-trees, since
217	! circ_classes are not defined for non-trees.
218	IF(is_tree(ivm))THEN
219
220	DO icir=1,ncirc
221
222	IF(test_pft == ivm .AND. test_grid == ipts .AND. ld_kill)&
223	WRITE(numout,*) 'killing before: ',ipts,ivm,icir,ifm,icut,&
224	circ_class_kill(ipts,ivm,icir,ifm,icut),&
225	circ_class_n(ipts,ivm,icir)
226
227	! move the dead biomass to the respective litter pool
228	bm_to_litter(ipts,ivm,:,iele) = bm_to_litter(ipts,ivm,:,iele) + &
229	circ_class_biomass(ipts,ivm,icir,:,iele)*&
230	circ_class_kill(ipts,ivm,icir,ifm,icut)
231
232	! remove the number of individuals that died
233	circ_class_n(ipts,ivm,icir) = circ_class_n(ipts,ivm,icir) - &
234	circ_class_kill(ipts,ivm,icir,ifm,icut)
235
236
237	! remove the dead biomass
238	biomass(ipts,ivm,:,iele) = biomass(ipts,ivm,:,iele) - &
239	circ_class_biomass(ipts,ivm,icir,:,iele)*&
240	circ_class_kill(ipts,ivm,icir,ifm,icut)
241
242	! adjust the number of individuals
243	ind(ipts,ivm)=ind(ipts,ivm)-&
244	circ_class_kill(ipts,ivm,icir,ifm,icut)
245
246	! Here, it's possible that the number of individuals left
247	! will be a very, very small amount. If it's very low,
248	! we will just move all that biomass to the dead pool
249	! and set the number of individuals equal to zero,
250	! effectively killing the circ class.
251	IF(circ_class_n(ipts,ivm,icir) .LE. min_stomate)THEN
252
253	bm_to_litter(ipts,ivm,:,iele) = bm_to_litter(ipts,ivm,:,iele) + &
254	circ_class_biomass(ipts,ivm,icir,:,iele)*&
255	circ_class_n(ipts,ivm,icir)
256	biomass(ipts,ivm,:,iele) = biomass(ipts,ivm,:,iele) - &
257	circ_class_biomass(ipts,ivm,icir,:,iele)*&
258	circ_class_n(ipts,ivm,icir)
259	ind(ipts,ivm)=ind(ipts,ivm)-&
260	circ_class_n(ipts,ivm,icir)
261
262	! zero the number of individuals in this circ class
263	circ_class_n(ipts,ivm,icir) = circ_class_n(ipts,ivm,icir) - &
264	circ_class_n(ipts,ivm,icir)
265
266	ENDIF
267
268	!!$ IF(test_pft == ivm .AND. test_grid == ipts)&
269	!!$ WRITE(numout,*) 'killing after: ',ipts,ivm,icir,ifm,icut,&
270	!!$ circ_class_kill(ipts,ivm,icir,ifm,icut),&
271	!!$ circ_class_n(ipts,ivm,icir)
272	ENDDO ! loop over circ classes
273
274
275	ELSE
276
277	!! WARNING !!
278	! This is not very clean. We don't use circ classes for grasses
279	! and crops. All mortality is done based on biomass. However,
280	! it's cleaner here to just pass circ_class_kill from lpj_kill,
281	! and it is in line with what is done for the forests. So we take
282	! the total grass/crop biomass that is scheduled for killing in
283	! lpj_kill and define circ_class_kill(:,:,inatural,1) and
284	! circ_class_biomass(:,:,1,:,:) so that the product of these two
285	! matches the amount of biomass scheduled for killing. We do NOT
286	! change the number of individuals after killing, either.
287
288	! Let us compute the fraction of plants which die of natural
289	! causes in this step.
290	! Individuals are not killed for grasses and crops, so we cannot
291	! calculate the mortality as a function of the individuals
292
293	! We renormalize by circ_class_n below. If it's zero,
294	! we have an error. Of course, if it's zero there is no
295	! reason to be here in the first place.
296	IF(circ_class_n(ipts,ivm,1) .LE. min_stomate)CYCLE
297
298	! move the dead biomass to the respective litter pool
299	bm_to_litter(ipts,ivm,:,iele) = bm_to_litter(ipts,ivm,:,iele) + &
300	circ_class_biomass(ipts,ivm,1,:,iele)*&
301	circ_class_kill(ipts,ivm,1,ifm,icut)
302
303	! remove the dead biomass
304	biomass(ipts,ivm,:,iele) = biomass(ipts,ivm,:,iele) - &
305	circ_class_biomass(ipts,ivm,1,:,iele)*&
306	circ_class_kill(ipts,ivm,1,ifm,icut)
307
308	! circ_class_biomass has to change such that the biomass
309	! and circ_class_biomass*circ_class_n are in sync. So we
310	! will just sync it here. This has to be done this way
311	! since circ_class_n never changes for grass/crops.
312	circ_class_biomass(ipts,ivm,1,:,iele) = biomass(ipts,ivm,:,iele)/&
313	circ_class_n(ipts,ivm,1)
314
315	! In order for something to be prescribed, ind must be zero.
316	! We don't change ind and circ_class_n if we just kill some
317	! of the biomass, but if we kill all of it we must set
318	! it to zero in order for it to be prescribed.
319	! IF(ivm == test_pft) WRITE(numout,*) 'Killing? ', &
320	! SUM(biomass(ipts,ivm,:,iele))
321	IF(SUM(biomass(ipts,ivm,:,iele)) .LE. min_stomate)THEN
322	circ_class_n(ipts,ivm,:)=zero
323	ind(ipts,ivm)=zero
324	ENDIF
325
326
327	ENDIF ! checking for a tree
328
329	! This stuff should never been killed again.
330	! To make sure of that, reset all the
331	! variables to zero.
332	circ_class_kill(ipts,ivm,:,ifm,icut)=zero
333	ENDDO
334	ENDDO
335
336	ENDDO ! loop over elements
337
338	ENDDO ! loop over pfts
339
340	END DO ! loop over land points
341
342	!! 4. Check mass balance closure
343
344	! Consider this whole routine as a black box with incoming and outgoing
345	! fluxes and a change in the mass of the box. Express in absolute
346	! units gC (or gN), hence, multiply with dt and veget_max. In most routines
347	! veget_max does not change and could be omitted but a general approach
348	! was prefered.
349
350	!! 4.1 Calculate pools at the end of the routine
351	pool_end = zero
352	DO ipar = 1,nparts
353	DO iele = 1,nelements
354	pool_end(:,:,iele) = pool_end(:,:,iele) + &
355	(biomass(:,:,ipar,iele) * veget_max(:,:)) + &
356	(bm_to_litter(:,:,ipar,iele) * veget_max(:,:))
357	ENDDO
358	ENDDO
359
360	!! 4.2 Calculate components of the mass balance
361	check_intern(:,:,iatm2land,icarbon) = zero
362	check_intern(:,:,iland2atm,icarbon) = -un * zero
363	check_intern(:,:,ilat2out,icarbon) = zero
364	check_intern(:,:,ilat2in,icarbon) = -un * zero
365	check_intern(:,:,ipoolchange,icarbon) = -un * &
366	(pool_end(:,:,icarbon) - pool_start(:,:,icarbon))
367	closure_intern = zero
368	DO imbc = 1,nmbcomp
369	closure_intern(:,:,icarbon) = closure_intern(:,:,icarbon) + &
370	check_intern(:,:,imbc,icarbon)
371	ENDDO
372
373	! Write outcome
374	DO ipts=1,npts
375	DO ivm=1,nvm
376	IF(ABS(closure_intern(ipts,ivm,icarbon)) .LE. min_stomate)THEN
377	IF (ld_massbal) WRITE(numout,*) 'Mass balance closure in stomate_kill.f90'
378	ELSE
379	WRITE(numout,*) 'Error: mass balance is not closed in stomate_kill.f90'
380	WRITE(numout,*) ' ipts,ivm; ', ipts,ivm
381	WRITE(numout,*) ' Difference is, ', closure_intern(ipts,ivm,icarbon)
382	WRITE(numout,*) ' pool_end,pool_start: ', pool_end(ipts,ivm,icarbon), &
383	pool_start(ipts,ivm,icarbon)
384	IF(ld_stop)THEN
385	CALL ipslerr_p (3,'gap_prognostic', 'Mass balance error.','','')
386	ENDIF
387	ENDIF
388	ENDDO
389	ENDDO
390
391	IF (bavard.GE.2) WRITE(numout,*) 'Leaving gap_prognostic'
392
393	END SUBROUTINE gap_prognostic
394
395
396	!! ================================================================================================================================
397	!! SUBROUTINE : gap_clean
398	!!
399	!>\BRIEF After biomass has been killed, there are some clean-up operations
400	!! to do.
401	!!
402	!! DESCRIPTION : If all the biomass has been killed for a PFT, we want to reset
403	!! some counters. There is also the chance that we will have
404	!! killed all the biomass in a circ class of trees. If this is
405	!! the case, we want to redistribute the biomass in the remaining
406	!! classes so that all circ classes have some biomass in them.
407	!! A circ class with no biomass causes allocation to crash.
408	!!
409	!! RECENT CHANGE(S):
410	!!
411	!! MAIN OUTPUT VARIABLE(S): ::circ_class_biomass
412	!!
413	!! REFERENCE(S) :
414	!!
415	!! FLOWCHART : None
416	!!\n
417	!_ ================================================================================================================================
418
419	SUBROUTINE gap_clean (npts, biomass, circ_class_biomass, &
420	circ_class_n, age, everywhere, leaf_age, &
421	leaf_frac, senescence, when_growthinit, circ_class_dist, &
422	age_stand, PFTpresent, last_cut, mai_count, &
423	veget_max, npp_longterm, KF, co2_to_bm, ind, &
424	wstress_season, lm_lastyearmax, lignin_struc, lignin_wood, &
425	carbon, litter, vir_wstress_fac, veget, dt_days, &
426	forest_managed, bm_to_litter, lab_fac, gdd_from_growthinit,&
427	gdd_midwinter, time_hum_min, gdd_m5_dormance, ncd_dormance,&
428	moiavail_month, moiavail_week, gpp_week, gpp_daily, resp_maint,&
429	resp_growth, npp_daily, use_reserve, mai, pai, ngd_minus5,&
430	rue_longterm, previous_wood_volume,&
431	tau_eff_leaf, tau_eff_sap, tau_eff_root, moiavail_growingseason,&
432	orphan_flux, k_latosa_adapt, lpft_replant)
433
434	!! 0. Variable and parameter declaration
435
436	!! 0.1 Input variables
437	INTEGER(i_std), INTENT(in) :: npts !! Domain size (-)
438	REAL(r_std), DIMENSION(ncirc), INTENT(in) :: circ_class_dist !! The probability distribution of trees
439	!! in a circ class in case of a
440	!! redistribution (unitless).
441	REAL(r_std), DIMENSION(:,:), INTENT(in) :: vir_wstress_fac !! Water stress factor, based on hum_rel_daily
442	!! (unitless, 0-1)
443	REAL(r_std), DIMENSION(:,:), INTENT(in) :: veget !! Fraction of vegetation type including
444	!! non-biological fraction (unitless)
445	REAL(r_std), DIMENSION(:,:), INTENT(in) :: tau_eff_root !! Effective root turnover time that accounts
446	!! waterstress (days)
447	REAL(r_std), DIMENSION(:,:), INTENT(in) :: tau_eff_sap !! Effective sapwood turnover time that accounts
448	!! waterstress (days)
449	REAL(r_std), DIMENSION(:,:), INTENT(in) :: tau_eff_leaf !! Effective leaf turnover time that accounts
450	!! waterstress (days)
451	REAL(r_std), INTENT(in) :: dt_days !! Time step of vegetation dynamics for stomate (days)
452	INTEGER(i_std), DIMENSION (:,:), INTENT(in) :: forest_managed !! forest management flag
453	LOGICAL, DIMENSION(npts,nvm), INTENT(in) :: lpft_replant !! Set to true if a PFT has been clearcut
454	!! and needs to be replaced by another species
455
456	!! 0.2 Output variables
457
458	!! 0.3 Modified variables
459	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: biomass !! Biomass @tex $(gC m^{-2}) $@endtex
460	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(inout) :: circ_class_biomass !! Biomass of the components of the model
461	!! tree within a circumference
462	!! class @tex $(gC ind^{-1})$ @endtex
463	LOGICAL, DIMENSION(:,:), INTENT(inout) :: senescence !! Is the plant senescent? (only for deciduous
464	!! trees - carbohydrate reserve) (true/false)
465	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: age !! Mean age (years)
466	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: when_growthinit !! How many days ago was the beginning of the
467	!! growing season (days)
468	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: everywhere !! Is the PFT everywhere in the grid box or very
469	!! localized (after its introduction)
470	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: leaf_age !! Leaf age (days)
471	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: leaf_frac !! Fraction of leaves in leaf age class
472	!! (unitless;0-1)
473	INTEGER(i_std), DIMENSION(:,:), INTENT(inout) :: age_stand !! Age of the forest stand (years)
474	INTEGER(i_std), DIMENSION(:,:), INTENT(inout) :: last_cut !! Years since last thinning (years)
475	INTEGER(i_std), DIMENSION(:,:), INTENT(inout) :: mai_count !! The number of times we've
476	!! calculated the volume increment
477	!! for a stand
478	LOGICAL, DIMENSION(:,:), INTENT(inout) :: PFTpresent !! PFT present (0 or 1)
479	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: KF !! Scaling factor to convert sapwood mass into leaf
480	!! mass (m)
481	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: k_latosa_adapt !! Leaf to sapwood area adapted for long
482	!! term water stress (m)
483	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: veget_max !! "maximal" coverage fraction of a PFT on the ground
484	!! (unitless, 0-1)
485	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: co2_to_bm !! CO2 taken from the atmosphere to get C to create
486	!! the seedlings @tex (gC.m^{-2}dt^{-1})$ @endtex
487	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ind !! Density of individuals at the stand level
488	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: wstress_season !! Water stress factor, based on hum_rel_daily
489	!! (unitless, 0-1)
490
491	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: npp_longterm !! "long term" net primary productivity
492	!! @tex ($gC m^{-2} year^{-1}$) @endtex
493	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: lm_lastyearmax !! last year's maximum leaf mass for each PFT
494	!! @tex ($gC m^{-2}$) @endtex
495	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: lignin_struc !! ratio Lignine/Carbon in structural litter,
496	!! above and below ground
497	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: lignin_wood !! ratio Lignine/Carbon in woody litter,
498	!! above and below ground
499	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: carbon !! carbon pool: active, slow, or passive
500	!! @tex ($gC m^{-2}$) @endtex
501	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(inout) :: litter !! metabolic and structural litter, above and
502	!! below ground @tex ($gC m^{-2}$) @endtex
503	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: circ_class_n !! Number of individuals in each circ class
504	!! @tex $(ind m^{-2})$ @endtex
505	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: bm_to_litter !! Transfer of biomass to litter
506	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
507	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: lab_fac !! Activity of labile pool factor (??units??)
508	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_from_growthinit !! growing degree days, since growthinit
509	!! for crops
510	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_midwinter !! Growing degree days (K), since midwinter
511	!! (for phenology) - this is written to the
512	!! history files
513	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: time_hum_min !! Time elapsed since strongest moisture
514	!! availability (days)
515	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gdd_m5_dormance !! Growing degree days (K), threshold -5 deg
516	!! C (for phenology)
517	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ncd_dormance !! Number of chilling days (days), since
518	!! leaves were lost (for phenology)
519	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: moiavail_month !! "Monthly" moisture availability (0 to 1,
520	!! unitless)
521	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: moiavail_week !! "Weekly" moisture availability
522	!! (0 to 1, unitless)
523	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gpp_week !! Mean weekly gross primary productivity
524	!! @tex $(gC m^{-2} day^{-1})$ @endtex
525	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: gpp_daily !! Daily gross primary productivity
526	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
527	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_maint !! Maintenance respiration
528	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
529	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: resp_growth !! Growth respiration
530	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
531	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: npp_daily !! Net primary productivity
532	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
533	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: use_reserve !! Mass taken from carbohydrate reserve
534	!! @tex $(gC m^{-2})$ @endtex
535	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: mai !! The mean annual increment
536	!! @tex $(m3 / m2 / year)$ @endtex
537	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: pai !! The period annual increment
538	!! @tex $(m3 / m2 / year)$ @endtex
539	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: ngd_minus5 !! Number of growing days (days), threshold
540	!! -5 deg C (for phenology)
541	!!$ REAL(r_std), DIMENSION(:,:), INTENT(inout) :: t_photo_stress !! Temperature stress for CEXCHANGE
542	!!$ !! photosynthesis (0-1)
543	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: rue_longterm !! Longterm radiation use efficiency
544	!! (??units??)
545	REAL(r_std), DIMENSION(npts,nvm), INTENT(inout) :: previous_wood_volume !! The volume of the tree trunks
546	!! in a stand for the previous year.
547	!! @tex $(m3 / m2 )$ @endtex
548	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: moiavail_growingseason !! Mean growingseason moisture
549	!! availability (0 to 1, unitless)
550	REAL(r_std), DIMENSION(:,:,:,:), INTENT(inout) :: orphan_flux !! Storage for fluxes of PFTs that no longer exist.
551	!! Following the total destruction of a PFT in LCC
552	!! (veget_max_new = 0), a flux from before LCC needs a
553	!! storage @tex $(gC m^{-2} dtslow^{-1})$ @endtex
554
555	!! 0.4 Local variables
556	REAL(r_std) :: share_young !! Share of the veget_max of the youngest age class
557	!! (unitless, 0-1)
558	INTEGER(i_std) :: ipts,ivm,ilage !! Indices
559	INTEGER(i_std) :: icir,jcir,ilev !! Indices
560	INTEGER(i_std) :: ipar,iele,imbc !! Indices
561	INTEGER(i_std) :: ilit,icarb !! Indices
562	INTEGER(i_std) :: inc !! Indices
563	LOGICAL :: lredistribute !! Flag if we need to redistribute individuals
564	!! among the circ classes
565	INTEGER,DIMENSION(nvm) :: nkilled !! the number of grid points at which a given
566	!! given PFT's biomass has been reduced to
567	!! zero
568	REAL(r_std), DIMENSION(npts,nvm,nmbcomp,nelements)&
569	:: check_intern !! Contains the components of the internal
570	!! mass balance chech for this routine
571	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
572	REAL(r_std), DIMENSION(npts,nvm,nelements) :: closure_intern !! Check closure of internal mass balance
573	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
574	REAL(r_std), DIMENSION(npts,nvm,nelements) :: pool_start !! Start pool of this routine
575	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
576	REAL(r_std), DIMENSION(npts,nvm,nelements) :: pool_end !! End pool of this routine
577	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
578	REAL(r_std), DIMENSION(ncirc) :: old_trees_left
579	REAL(r_std), DIMENSION(ncirc) :: circ_class_n_new
580	REAL(r_std), DIMENSION(ncirc,nparts) :: circ_class_biomass_new
581	REAL(r_std) :: trees_needed
582	LOGICAL :: lyoungest
583	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements):: new_biomass !! Temporary variable for biomass
584	!! @tex $(gC.m^{-2})$ @endtex
585	REAL(r_std), DIMENSION(npts,nvm,nleafages) :: new_leaf_frac !! Temporary variable for leaf_frac (unitless;0-1)
586	REAL(r_std), DIMENSION(npts,nvm) :: new_ind !! Temporary variable for ind @tex $(m^{-2})$ @endtex
587	REAL(r_std), DIMENSION(ncirc) :: est_circ_class_n !! Temporary variable for circ_class_n of the
588	!! established vegetation @tex $(ind m^{-2})$ @endtex
589	REAL(r_std), DIMENSION(ncirc) :: tmp_circ_class_n !! Temporary variable for circ_class_n of the
590	!! established vegetation @tex $(ind m^{-2})$ @endtex
591	REAL(r_std), DIMENSION(npts,nvm,ncirc) :: new_circ_class_n !! Temporary variable for circ_class_n of the
592	!! new vegetation
593	!! @tex $(ind m^{-2})$ @endtex
594	REAL(r_std), DIMENSION(ncirc,nparts,nelements) :: est_circ_class_biomass !! Temporary variable for circ_class_biomass of the
595	!! established vegetation @tex $(g C ind^{-1})$ @endtex
596	REAL(r_std), DIMENSION(ncirc,nparts,nelements) :: tmp_circ_class_biomass !! Temporary variable for circ_class_biomass of the
597	!! established vegetation @tex $(g C ind^{-1})$ @endtex
598	REAL(r_std), DIMENSION(npts,nvm,ncirc,nparts,nelements) &
599	:: new_circ_class_biomass !! Temporary variable for circ_class_biomass of the
600	!! new vegetation @tex $(g C ind^{-1})$ @endtex
601	REAL(r_std), DIMENSION(ncirc) :: est_circ_class_dia !! Variable to store temporary values for
602	!! circ_class (m)
603	REAL(r_std), DIMENSION(ncirc) :: tmp_circ_class_dia !! Variable to store temporary values for
604	!! circ_class (m)
605	REAL(r_std), DIMENSION(npts,nvm,ncirc) :: new_circ !! Temporary variable for circ (m)
606	REAL(r_std), DIMENSION(npts,nvm) :: new_co2_to_bm !! Temporary variable for co2_to_bm
607	!! @tex (gC.m^{-2}dt^{-1})$ @endtex
608	REAL(r_std), DIMENSION(npts,nvm) :: new_age !! Temporary variable for age (years)
609	REAL(r_std), DIMENSION(npts,nvm) :: new_lm_lastyearmax !! Temporary variable for lm_last_year
610	!! @tex ($gC m^{-2}$) @endtex
611	REAL(r_std), DIMENSION(npts,nvm) :: new_everywhere !! Temporary variable for everywhere (unitless, 0-1)
612	REAL(r_std), DIMENSION(npts,nvm) :: dum_when_growthinit !! Dummy for when_growthinit (days)
613	REAL(r_std), DIMENSION(npts,nvm) :: dum_npp_longterm !! Dummy for npp_long_term
614	!! @tex ($gC m^{-2} year^{-1}$) @endtex
615	INTEGER(i_std) :: ivma
616	LOGICAL, DIMENSION(npts,nvm) :: dum_PFTpresent !! Dummy for PFTpresent (0 or 1)
617	LOGICAL, DIMENSION(npts,nvm) :: dum_senescence !! Dummy for senescence (true/false)
618	REAL(r_std) :: moi_bank !! bank to store moiavail_growingseason of the cleared
619	!! PFT's (unitless; 0-1)
620	REAL(r_std), DIMENSION(npts,nvm) :: loss_gain !! The same as delta_veget but distributed of all
621	!! age classes and thus taking the age-classes into
622	!! account (unitless, 0-1)
623	REAL(r_std) :: total_losses !! Sum of losses only in delta_veget (unitless, 0-1)
624	REAL(r_std), DIMENSION(nlitt,nlevs,nelements) :: litter_bank !! Bank to store the litter that becomes available when
625	!! part of a PFT is removed @tex ($gC m^{-2}$) @endtex
626	REAL(r_std), DIMENSION(ncarb) :: soil_bank !! Bank to store soil carbon that becomes available when
627	!! part of a PFT is removed @tex ($gC m^{-2}$) @endtex
628	REAL(r_std), DIMENSION(nlevs) :: struct_ltr_bank !! Bank to store the lignin in structural litter that
629	!! becomes available when part of a PFT is removed
630	!! (0-1,unitless)
631	REAL(r_std),DIMENSION(nlevs) :: woody_ltr_bank !! Bank to store the lignin in woody litter that
632	!! becomes available when part of a PFT is removed
633	!! (0-1,unitless)
634	REAL(r_std),DIMENSION(nvm,nparts,nelements) :: fresh_litter !! Pool of fresh litter that is being released during
635	!! LCC @tex ($gC m^{-2}$) @endtex
636	REAL(r_std),DIMENSION(nparts,nelements) :: fresh_ltr_bank !! Bank to store all the fresh litter from site
637	!! clearing during LCC. Weighted value of fresh_litter
638	!! @tex ($gC m^{-2}$) @endtex
639	INTEGER(i_std) :: iyoung
640	REAL(r_std), DIMENSION(npts,nvm,norphans,nelements) &
641	:: orphan_flux_local !! Storage for fluxes of PFTs that no longer exist.
642	!! This is only for co2bm at the moment, since that
643	!! is the only one used in this routine. Following the
644	!! total destruction of a PFT by death
645	!! (veget_max_new = 0), a flux from before death needs
646	!! storage @tex $(gC m^{-2} dtslow^{-1})$ @endtex
647	REAL(r_std), DIMENSION(nlitt,nlevs) :: litter_weight_young !! The fraction of litter on the young
648	!! PFT.
649	!! @tex $-$ @endtex
650	REAL(r_std), DIMENSION(nparts) :: v1 !! temporay variable for sorting circ_class_biomass
651	REAL(r_std) :: v2 !! temporay variables for sorting circ_class_n
652	REAL(r_std) :: sort !! flag triggering sorting routine (0 or 1)
653	REAL(r_std), DIMENSION(ncirc) :: tree_size !! temporay variables for checking the biomass
654	!! for each circ class @tex ($gC tree^{-1}$) @endtex
655	REAL(r_std), DIMENSION(npts,nvm) :: veget_max_begin !! Temporairy variable to check area conservation
656
657	!_ ================================================================================================================================
658
659	IF (bavard.GE.2) WRITE(numout,*) 'Entering gap_clean.'
660
661	!! 1. Initialize check for mass balance closure
662
663	! The mass balance is calculated at the end of this routine
664	! in section .
665	pool_start = zero
666	DO iele = 1,nelements
667	DO ipar = 1,nparts
668	! Initial biomass
669	pool_start(:,:,iele) = pool_start(:,:,iele) + &
670	(biomass(:,:,ipar,iele) * veget_max(:,:))
671	ENDDO
672
673	! The litter and soil carbon pools can be moved around
674	! if age classes are changed.
675
676	! Litter pool (gC m-2) * (m2 m-2)
677	DO ilit = 1,nlitt
678	DO ilev = 1,nlevs
679	pool_start(:,:,iele) = pool_start(:,:,iele) + &
680	litter(:,ilit,:,ilev,iele) * veget_max(:,:)
681	ENDDO
682	ENDDO
683
684	! Soil carbon (gC m-2) * (m2 m-2)
685	DO icarb = 1,ncarb
686	pool_start(:,:,iele) = pool_start(:,:,iele) + &
687	carbon(:,icarb,:) * veget_max(:,:)
688	ENDDO
689
690	ENDDO
691
692	! It is possible that the orphan fluxes are not zero here. For example,
693	! if we prescribed biomass at the beginning of this timestep,
694	! ico2bm will be some value. We are only interested in what happens in this
695	! step for the moment, so we use a local variable for this flux.
696	orphan_flux_local(:,:,:,:) = zero
697
698	! Initialize area check
699	veget_max_begin(:,:) = veget_max(:,:)
700
701	!! 2. Redistribute biomass
702
703	DO ipts = 1,npts ! loop over land_points
704
705	DO ivm = 2,nvm ! loop over #PFT
706
707	IF(veget_max(ipts,ivm) == zero)THEN
708	! this vegetation type is not present, so no reason to do the
709	! calculation.
710	CYCLE
711	ENDIF
712
713	! First we check to see if one of the circ classes is empty. We only need
714	! to do this if there is some biomass, since we don't care if all
715	! of the circ classes are empty.
716	IF(SUM(biomass(ipts,ivm,:,icarbon)) .GT. min_stomate)THEN
717
718	IF(is_tree(ivm))THEN
719
720	! By setting this to .TRUE. redistribution
721	! will be taken care of every day. Initially this flag was
722	! only true at the end of the year but that resuted in
723	! spikes and pulses.
724	lredistribute=.TRUE.
725
726	IF(lredistribute) THEN
727
728	!---TEMP---
729	IF(ld_biomass .AND. ivm==test_pft)THEN
730	WRITE(numout,*) 'Start of redistributing biomass in gap_clean'
731	WRITE(numout,*) 'ipts,ivm',ipts,ivm
732	WRITE(numout,*) 'circ_class_n, ',circ_class_n(ipts,ivm,:)
733	WRITE(numout,*) 'circ_class_biomass, ', &
734	SUM(circ_class_biomass(ipts,ivm,:,:,icarbon),2)
735	ENDIF
736	!-----------
737
738	! What we want to do is recreate the circumference class distribution,
739	! since if we have hit this point we have at least one class that
740	! is empty and this will cause the allocation to crash.
741
742	! The first step is to decide of the distribution of individuals
743	! that we want in each class, for example, a uniform or exponential
744	! distribution. This was made this an input parameter. Note
745	! that this is a decisison with very high consequences. The diamter
746	! range may vary but the distribution is always the same.
747
748	! Now we need to populate the new distribution of trees among
749	! the circumference classes, and move the heartwood and sap masses
750	! to the new distributions. We also need to move the other pools.
751	! This is tricky because the allometric relations for the new
752	! tree sizes will NOT give the same amount of biomass for
753	! the non-woody pools. Let us first try it just distributing
754	! everything equally, and then if that doesn't work we can
755	! redistribute the non-woody pools more cleverly, even if that
756	! means we change the total amount of biomass we have in this
757	! stand. I want to try it this way because we conserve biomass
758	! this way, and I hope that the stress on the trees will be
759	! small enough that it doesn't cause problems. This
760	! redistribution should only happen rarely, so the trees
761	! should have a chance to equilibrate.
762	circ_class_n_new(:)=circ_class_dist(:)*SUM(circ_class_n(ipts,ivm,:))
763	circ_class_biomass_new(:,:)=zero
764
765	! now we populate the new classes
766	old_trees_left(:)=circ_class_n(ipts,ivm,:)
767
768	! The subsequent calculation nicely closes the carbon
769	! balance but within the precision 10-16. If there are
770	! enough trees (first case in the IF-loop), we introduce
771	! a precision error in the number of trees. This implies
772	! that when the number of trees is very small, the
773	! precision issue can become a numerical issue. This is
774	! especially true for the largest circumference classes
775	! because they contains the fewest individuals.
776	! We observed the largest circumference classes
777	! containing 10-4 to 10-5 gC less than the previous
778	! class. The rest of the code shouldn't care too much
779	! whether the circumferences classes are in order or not
780	! but to reduce the chances this problem occurs
781	! we inverted the DO-loops. That way the precision error
782	! occurs in carried to the smaller diameter classes
783	! which have more individuals and so the precision error
784	! stays a precision issue. Note that this routine is
785	! called every day so the problem is likely to be
786	! corrected the next day.
787	DO icir=ncirc,1,-1
788
789	trees_needed=circ_class_n_new(icir)
790
791	DO jcir=ncirc,1,-1
792
793	IF(trees_needed .LE. zero)EXIT
794
795	!+++TEMP+++
796	IF(ld_biomass .AND. ivm==test_pft)THEN
797	WRITE(numout,*) 'start of loop'
798	WRITE(numout,*) 'circ_class_biomass_new, ',&
799	icir,SUM(circ_class_biomass_new(icir,:))
800	WRITE(numout,*) 'old_trees_left, ', &
801	jcir, old_trees_left(jcir)
802	WRITE(numout,*) 'trees needed, ',trees_needed
803	ENDIF
804	!++++++++++
805
806	IF(old_trees_left(jcir) .GE. trees_needed )THEN
807
808	! we can get all the trees we need from this class
809	! and don't have to continue searching
810	circ_class_biomass_new(icir,:)=circ_class_biomass_new(icir,:)+&
811	circ_class_biomass(ipts,ivm,jcir,:,icarbon)*trees_needed
812	old_trees_left(jcir)=old_trees_left(jcir)-trees_needed
813	trees_needed = zero
814
815	!+++TEMP+++
816	IF(ld_biomass .AND. ivm==test_pft)THEN
817	WRITE(numout,*) 'enough trees'
818	ENDIF
819	!++++++++++
820
821	EXIT
822
823	ELSE
824
825	! the trees in this class are not sufficient, so we
826	! need to move all of them to the new class.
827	circ_class_biomass_new(icir,:)=circ_class_biomass_new(icir,:)+&
828	circ_class_biomass(ipts,ivm,jcir,:,icarbon)*&
829	old_trees_left(jcir)
830	trees_needed = trees_needed-old_trees_left(jcir)
831	old_trees_left(jcir)=zero
832
833	!+++TEMP+++
834	IF(ld_biomass .AND. ivm==test_pft)THEN
835	WRITE(numout,*) 'not enough trees'
836	ENDIF
837	!++++++++++
838
839	ENDIF
840
841	!+++TEMP+++
842	IF(ld_biomass .AND. ivm==test_pft)THEN
843	WRITE(numout,*) 'start of loop'
844	WRITE(numout,*) 'circ_class_biomass_new, ',&
845	icir,SUM(circ_class_biomass_new(icir,:))
846	WRITE(numout,*) 'old_trees_left, ', &
847	jcir, old_trees_left(jcir)
848	WRITE(numout,*) 'trees needed, ',trees_needed
849	ENDIF
850	!++++++++++
851
852	ENDDO ! jcir=1,ncirc
853
854	ENDDO ! icir=1,ncirc
855
856	! right now, circ_class_biomass_new gives the total biomass
857	! in each class, but it should only be for a model tree.
858	! So let's normalize it.
859	DO icir=1,ncirc
860	circ_class_biomass_new(icir,:)=circ_class_biomass_new(icir,:)/&
861	circ_class_n_new(icir)
862	ENDDO
863
864	! Check whether the order was preserved. We expect that the
865	! biomass of an individual tree increases with an increasing
866	! circ class. This is probably not important for the correct
867	! functioning of the model. We try to preserve the order as
868	! an additional mean to control the quality of the simulations
869	! If the order was not preserved we will sort the biomasses.
870	tree_size(:)=zero
871	DO icir=1,ncirc
872
873	tree_size(icir)=SUM(circ_class_biomass_new(icir,:))
874
875	ENDDO
876	sort = zero
877	DO icir=2,ncirc
878
879	IF(tree_size(icir) .LT. tree_size(icir-1))THEN
880
881	! The order is not preserved
882	sort = un
883	EXIT
884
885	ENDIF
886
887	ENDDO
888
889	! Sort the biomasses. Note that at the same time we also
890	! have to sort circ_class_n. We made use of the Shell sort
891	! method
892	IF (sort==1) THEN
893
894	!+++TEMP+++
895	IF(ld_biomass)THEN
896	WRITE(numout,*) 'gap_clean: sort begin', &
897	SUM(circ_class_biomass_new(:,:),2),&
898	circ_class_n_new(:)
899	ENDIF
900	!++++++++++
901
902	inc = 1
903	DO
904	inc=3*inc+1
905	IF (inc .GT. ncirc) EXIT
906	ENDDO
907	DO
908	inc = inc/3
909	DO icir=inc+1,ncirc
910	v1(:) = circ_class_biomass_new(icir,:)
911	v2 = circ_class_n_new(icir)
912	jcir=icir
913	DO
914	IF (SUM(circ_class_biomass_new(jcir-inc,:)) .LE. SUM(v1(:))) EXIT
915	circ_class_biomass_new(jcir,:) = &
916	circ_class_biomass_new(jcir-inc,:)
917	circ_class_n_new(jcir) = circ_class_n_new(jcir-inc)
918	jcir = jcir-inc
919	IF (jcir .LE. inc) EXIT
920	ENDDO
921	circ_class_biomass_new(jcir,:)=v1(:)
922	circ_class_n_new(jcir)=v2
923	ENDDO
924	IF (inc .LE. 1) EXIT
925	ENDDO
926
927	!+++TEMP+++
928	IF(ld_biomass)THEN
929	WRITE(numout,*) 'gap_clean: sort end', &
930	SUM(circ_class_biomass_new(:,:),2),&
931	circ_class_n_new(:)
932	ENDIF
933	!++++++++++
934
935	ENDIF
936
937	! now update the variables that pass this information around
938	DO icir=1,ncirc
939	circ_class_biomass(ipts,ivm,icir,:,icarbon)=circ_class_biomass_new(icir,:)
940	circ_class_n(ipts,ivm,icir)=circ_class_n_new(icir)
941	ENDDO
942
943
944	ENDIF ! redistribute
945
946	!---TEMP---
947	IF(ld_biomass .AND. ivm==test_pft)THEN
948	WRITE(numout,*) 'End of redistributing biomass in gap_clean'
949	WRITE(numout,*) 'ipts,ivm',ipts,ivm
950	WRITE(numout,*) 'circ_class_n, ',circ_class_n(ipts,ivm,:)
951	WRITE(numout,*) 'circ_class_biomass, ', &
952	SUM(circ_class_biomass(ipts,ivm,:,:,icarbon),2)
953	ENDIF
954	!----------
955
956	ELSE
957
958	! we have biomass here but it's not a tree. Don't need to do anything.
959
960	ENDIF
961
962	ELSE
963
964	! All the biomass was killed for this site. Some flags to
965	! be reset in this case. It is OK to kill the PFT even if
966	! we are using species change but we can not replant yet
967	! because we want to replant with a different species this is
968	! basically the same as a land cover change and so it is best
969	! taken care of with the sapiens_lcchange code. To end up
970	! here we need a veget_max for this PFT but the biomass is
971	! has to be total. Another check is through using ::PFTpresent.
972	IF(PFTpresent(ipts,ivm))THEN
973
974	nkilled(ivm)=nkilled(ivm)+1 ! purely an informational variable
975
976	!! 1.5 Reinitialize vegetation characteristics in STOMATE
977	senescence(ipts,ivm) = .FALSE.
978	age(ipts,ivm) = zero
979	when_growthinit(ipts,ivm) = large_value
980
981	! Specific variables for forest stands
982	IF(is_tree(ivm))THEN
983
984	last_cut(ipts,ivm) = zero
985	mai_count(ipts,ivm) = zero
986	age_stand(ipts,ivm) = zero
987
988	ENDIF
989
990	!! 1.6 Update leaf ages
991	DO ilage = 1, nleafages
992
993	leaf_age(ipts,ivm,ilage) = zero
994	leaf_frac(ipts,ivm,ilage) = zero
995
996	ENDDO
997
998	! This used to be done in lpj_kill. If we have grasses,
999	! we reset the long term GPP. This value as 500 in the
1000	! old code. We externalized the variable, but we are
1001	! leaving the default at 500. Many PFTs should have
1002	! a long term value of less than 500, so someone could
1003	! do some work on this in the future.
1004	IF (.NOT. is_tree(ivm) .AND. .NOT. control%ok_constant_mortality) THEN
1005
1006	npp_longterm(ipts,ivm) = npp_reset_value(ivm)
1007
1008	ENDIF
1009
1010	ENDIF ! is PFT present?
1011
1012	ENDIF ! check if biomass is equal to zero
1013
1014	! If we have no vegetation left, this stand has been completely cleared.
1015	! We have to reset the age, regardless if it was for natural or human reasons.
1016	IF(is_tree(ivm))THEN
1017
1018	! This is essentially the same code that is found in land cover change,
1019	! since the same process happens there. Notice that here we do not
1020	! need to use the _bank variables, since we assume that if a stand dies
1021	! due to natural causes it will always be replaced by the same species.
1022	! We also assume that if a stand is clearcut, the model will replace it
1023	! by the same species. To do so otherwise would require something more
1024	! like a DGVM. That is what is being done in the species change code.
1025	! If the species change code is used, the PFT will be replanted at the
1026	! end of the year. Not in the middle of the year as being done here.
1027
1028	! If we are using age classes, we need to reset a lot of counters
1029	! and change veget_max, moving veget_max from this age class to
1030	! the lowest age class. If veget_max is greater than zero and there
1031	! is no biomass, prescribe will attempt to grow trees here in the
1032	! next timestep. This is fine if we have no age classes, but it doesn't
1033	! make any sense to prescribe trees that are 50 years old. We should
1034	! only prescribe trees which are 0 years old.
1035	IF( SUM(biomass(ipts,ivm,:,icarbon)) .LT. min_stomate .AND. &
1036	.NOT.lpft_replant(ipts,ivm))THEN
1037
1038	!---TEMP---
1039	IF(ld_forestry .OR. ld_kill)THEN
1040	WRITE(numout,*) 'Getting reset in stomate_kill'
1041	WRITE(numout,*) 'ivm,ipts',ivm,ipts
1042	ENDIF
1043	!----------
1044
1045	IF(nagec .GT. 1)THEN
1046
1047	! First, we need to find the youngest age class of this PFT.
1048	iyoung=start_index(agec_group(ivm))
1049	IF(ivm == iyoung)THEN
1050	lyoungest=.TRUE.
1051	ELSE
1052	lyoungest=.FALSE.
1053	ENDIF
1054
1055	IF(lyoungest)THEN
1056
1057	! We don't need to do anything. Things will be handled properly with
1058	! prescribe in the next step.
1059
1060	ELSE
1061
1062	! All of our counters for this stand are zero. We need to merge that
1063	! with the youngest age class. There are a couple possibilities.
1064	! One is that nothing exists right now in the youngest age class, in
1065	! which case the veget_max of the old age class becomes that of the
1066	! youngest and we prescribe. Another possibility is that something
1067	! does exist in the youngest age class, in which case we prescribe
1068	! to the current age class and then merge biomass with the youngest.
1069	! A third is that both the youngest age class and the current age
1070	! class have died. In that case we can do the same thing as if there
1071	! is no vegetation in the young age class, we just have to make sure
1072	! to add the veget_max and not replace it.
1073	! Veget_max after PFT expansion. As ::veget_max is used in the
1074	! IF-statements we won't update it until the end of this routine. If
1075	! we would updatebefore that, the same PFT may be subject to
1076	! conflicting actions.
1077	share_young = veget_max(ipts,iyoung) / &
1078	( veget_max(ipts,iyoung) + veget_max(ipts,ivm) )
1079
1080	! We also need a scaling factor which includes the litter
1081	DO ilev=1,nlevs
1082	DO ilit=1,nlitt
1083	IF(litter(ipts,ilit,iyoung,ilev,icarbon) .GT. min_stomate)THEN
1084
1085	litter_weight_young(ilit,ilev)=&
1086	litter(ipts,ilit,iyoung,ilev,icarbon)*&
1087	veget_max(ipts,iyoung)/ &
1088	(litter(ipts,ilit,ivm,ilev,icarbon)*veget_max(ipts,ivm) + &
1089	litter(ipts,ilit,iyoung,ilev,icarbon)*&
1090	veget_max(ipts,iyoung))
1091
1092	ELSE
1093
1094	litter_weight_young(ilit,ilev)=zero
1095
1096	ENDIF
1097	END DO
1098	ENDDO
1099
1100	IF( SUM(biomass(ipts,iyoung,:,icarbon)) .LT. min_stomate)THEN
1101
1102	IF(ld_agec)THEN
1103	WRITE(numout,*) 'Merging our biomass to the youngest age class.'
1104	WRITE(numout,*) 'No current biomass in the youngest age class.'
1105	WRITE(numout,*) 'ipts,iyoung,ivm: ',ipts,iyoung,ivm
1106	WRITE(numout,*) 'share_young: ',share_young
1107	WRITE(numout,*) 'veget_max(young and current): ',&
1108	veget_max(ipts,iyoung),veget_max(ipts,ivm)
1109	ENDIF
1110
1111	! Is there anything in the smallest age class? It
1112	! does not matter if the youngest age class died in this
1113	! step or not.
1114	veget_max(ipts,iyoung)=veget_max(ipts,iyoung)+veget_max(ipts,ivm)
1115
1116	moiavail_growingseason(ipts,iyoung) = &
1117	share_young * moiavail_growingseason(ipts,iyoung) + &
1118	moiavail_growingseason(ipts,ivm) * &
1119	(un - share_young)
1120
1121	! I am not sure why we merge wstress here. wstress is recalculated
1122	! everyday from moiavail_growingseason.
1123	wstress_season(ipts,iyoung) = &
1124	share_young * wstress_season(ipts,iyoung) + &
1125	wstress_season(ipts,ivm) * (un - share_young)
1126
1127	! The litter variables also need to be merged, since these will not
1128	! get updated in prescribe in the next step.
1129	litter(ipts,:,iyoung,:,:) = &
1130	share_young * litter(ipts,:,iyoung,:,:) + &
1131	litter(ipts,:,ivm,:,:) * &
1132	(un - share_young)
1133	carbon(ipts,:,iyoung) = &
1134	share_young * carbon(ipts,:,iyoung) + &
1135	carbon(ipts,:,ivm) * &
1136	(un - share_young)
1137	DO ilev=1,nlevs
1138	lignin_struc(ipts,iyoung,ilev) = &
1139	litter_weight_young(istructural,ilev) * &
1140	lignin_struc(ipts,iyoung,ilev) + &
1141	(un - litter_weight_young(istructural,ilev)) * &
1142	lignin_struc(ipts,ivm,ilev)
1143	lignin_wood(ipts,iyoung,ilev) = &
1144	litter_weight_young(iwoody,ilev) * &
1145	lignin_wood(ipts,iyoung,ilev) + &
1146	(un - litter_weight_young(iwoody,ilev) ) * &
1147	lignin_wood(ipts,ivm,ilev)
1148	ENDDO
1149
1150	bm_to_litter(ipts,iyoung,:,:) = &
1151	share_young * bm_to_litter(ipts,iyoung,:,:) + &
1152	bm_to_litter(ipts,ivm,:,:) * &
1153	(un - share_young)
1154
1155
1156	ELSE
1157
1158	! There are two important things to do here. First,
1159	! we need to prescribe biomass to this PFT, since it
1160	! is currently empty. Then we need to merge this biomass
1161	! with what is already in the youngest age class of this PFT.
1162
1163	! We want to make use of prescribe_prognostic but it should
1164	! be noted that the youngest PFT already contains biomass.
1165	! Therefore we will create a temporary PFT without biomass which
1166	! can be used to establish the new vegetation within the
1167	! youngest. For most of the INTENT(inout) variables of
1168	! prescribe_prognostic we receive temporary variables back this
1169	! helps us to calculate the weighted mean of the newly established
1170	! vegetation and the vegetation that is already there. For some
1171	! other variables we receive dummies as we are simply not
1172	! going to use these values but will continue using the values of the
1173	! vegetation that is already there.
1174
1175	IF(ld_agec)THEN
1176	WRITE(numout,*) 'Merging our biomass to the youngest age class.'
1177	WRITE(numout,*) 'Biomass already in the youngest age class.'
1178	WRITE(numout,*) 'ipts,iyoung,ivm: ',ipts,iyoung,ivm
1179	WRITE(numout,*) 'share_young: ',share_young
1180	WRITE(numout,*) 'veget_max(young and current): ',&
1181	veget_max(ipts,iyoung),veget_max(ipts,ivm)
1182	ENDIF
1183
1184	!! 5.2.3.1 Initialize new and dummy variables
1185	new_biomass(:,:,:,:) = zero
1186	new_circ_class_biomass(:,:,:,:,:) = zero
1187	new_ind(:,:) = zero
1188	new_circ_class_n(:,:,:) = zero
1189	new_lm_lastyearmax(:,:) = zero
1190	new_age(:,:) = zero
1191	new_leaf_frac(:,:,:) = zero
1192	new_co2_to_bm(:,:) = zero
1193	new_everywhere(:,:) = zero
1194	dum_when_growthinit(:,:) = large_value
1195	dum_npp_longterm(:,:) = zero
1196	dum_PFTpresent(:,:) = .TRUE.
1197	dum_senescence(:,:) = .FALSE.
1198
1199	!! 5.2.3.2 Merge the new and already available vegetation
1200
1201	! Assign value to moiavail_growingseason
1202	! We do this differently from LCC, because we are replanting
1203	! the same stand that we just cut.
1204	moiavail_growingseason(ipts,iyoung) = &
1205	moiavail_growingseason(ipts,iyoung) * &
1206	share_young + moiavail_growingseason(ipts,ivm) * &
1207	(un - share_young)
1208	wstress_season(ipts,iyoung) = wstress_season(ipts,iyoung) * &
1209	share_young + wstress_season(ipts,ivm) * (un - share_young)
1210
1211	!! 5.2.3.3 Initialize the newly established vegetation:
1212	! Initialize the newly established vegetation: density of
1213	! individuals, biomass, allocation factors, leaf age distribution,
1214	! etc to some reasonable value ::veget_max is not updated yet,
1215	! therefore loss_gain is passed in the argument list.
1216	! For loss_gain, we only want to replant the current PFT.
1217	! If lpft_replat is TRUE we should never be here in the
1218	! first place but to be safe lpft_replant is passed into
1219	! prescribe. This is an optional variable that will prevent
1220	! prescribing biomass if species change is used.
1221
1222	!---TEMP----
1223	IF(ld_species)THEN
1224	IF(lpft_replant(ipts,ivm))THEN
1225	WRITE(numout,*) 'ERROR - gap_clean should never be here'
1226	CALL ipslerr_p(3,'gap_clean in stomate_kill.f90',&
1227	'species change was activated, replanting needed',&
1228	'do not replant in gap_clean, wait for lcchange','')
1229	ENDIF
1230	ENDIF
1231	!-----------
1232
1233	loss_gain(:,:)=zero
1234	loss_gain(ipts,ivm)=veget_max(ipts,ivm)
1235	CALL prescribe_prognostic (npts,loss_gain,veget, dt_days,dum_PFTpresent, &
1236	new_everywhere, dum_when_growthinit, new_biomass, new_leaf_frac, &
1237	new_ind, new_circ_class_n, new_circ_class_biomass, new_co2_to_bm, &
1238	forest_managed, KF, &
1239	dum_senescence, new_age, dum_npp_longterm, new_lm_lastyearmax, &
1240	tau_eff_leaf, tau_eff_sap, tau_eff_root, k_latosa_adapt, &
1241	lpft_replant)
1242
1243	!! 5.2.3.3.1 Merge biomass of new and already available trees
1244	! Unlike LCC, we are only at this point if we are dealing
1245	! with forests, so we don't have to check for that here.
1246
1247	! Copy the number of individuals and the biomass of the established
1248	! vegetation to a temporary variable. In the subroutine ::circ_class_n
1249	! and ::circ_class_biomass will be overwritten with the characteristics
1250	! of the merged vegetation
1251	! The term "established" here refers to the youngest age class, while
1252	! "tmp" is the prescribed vegetation that we just created.
1253	est_circ_class_n(:) = circ_class_n(ipts,iyoung,:)
1254	est_circ_class_biomass(:,:,:) = circ_class_biomass(ipts,iyoung,:,:,:)
1255	tmp_circ_class_n(:) = new_circ_class_n(ipts,ivm,:)
1256	tmp_circ_class_biomass(:,:,:) = new_circ_class_biomass(ipts,ivm,:,:,:)
1257
1258	! Store circumference (m) in a temporary variable that can be
1259	! sorted from small to large - note that the dimension of these
1260	! variables are 2*ncirc
1261	est_circ_class_dia(:) = &
1262	wood_to_dia(est_circ_class_biomass(:,:,icarbon),ivm)
1263	tmp_circ_class_dia(:) = &
1264	wood_to_dia(tmp_circ_class_biomass(:,:,icarbon),ivm)
1265
1266	!!$ !---TEMP---
1267	!!$ IF(ld_species)THEN
1268	!!$ WRITE(numout,*) 'stomate_kill share_young, ',ipts, ivm, share_young
1269	!!$ WRITE(numout,*) 'rest, ', SUM(circ_class_n(ipts,ivm,:)), &
1270	!!$ SUM(est_circ_class_n), SUM(tmp_circ_class_n), &
1271	!!$ SUM(SUM(circ_class_biomass(ipts,ivm,:,:,icarbon),2),1), &
1272	!!$ SUM(est_circ_class_biomass), &
1273	!!$ SUM(tmp_circ_class_biomass), &
1274	!!$ ipts, iyoung, SUM(est_circ_class_dia), SUM(tmp_circ_class_dia), &
1275	!!$ lpft_replant(ipts,ivm)
1276	!!$ ENDIF
1277	!!$ !----------
1278
1279	! Merge the biomass
1280	CALL merge_biomass_pfts(npts, share_young, circ_class_n, &
1281	est_circ_class_n, tmp_circ_class_n, circ_class_biomass, &
1282	est_circ_class_biomass, &
1283	tmp_circ_class_biomass, ind, biomass, circ_class_dist, &
1284	ipts, iyoung, est_circ_class_dia, tmp_circ_class_dia)
1285
1286	!! 5.2.3.4 Calculate the PFT characteristics of the merged PFT
1287	! Take the weighted mean of the existing vegetation and the new
1288	! vegetation in this PFT. Note that co2_to_bm is in gC. m-2 dt-1
1289	! so we should also take the weighted mean (rather than sum if
1290	! this where absolute values).
1291	lm_lastyearmax(ipts,iyoung) = share_young * &
1292	lm_lastyearmax(ipts,iyoung) + &
1293	(un - share_young) * new_lm_lastyearmax(ipts,ivm)
1294	age(ipts,iyoung) = share_young * age(ipts,iyoung) + &
1295	(un - share_young) * new_age(ipts,ivm)
1296	leaf_frac(ipts,iyoung,:) = share_young * leaf_frac(ipts,iyoung,:) + &
1297	(un - share_young) * new_leaf_frac(ipts,ivm,:)
1298	co2_to_bm(ipts,iyoung) = share_young * co2_to_bm(ipts,iyoung) + &
1299	(un - share_young) * new_co2_to_bm(ipts,ivm)
1300
1301	! Update the orphan flux for the CO2 turned into biomass by prescribe,
1302	! as well as the area of vegetation regrown.
1303	orphan_flux_local(ipts,ivm,ico2bm,icarbon) = new_co2_to_bm(ipts,ivm)
1304	orphan_flux_local(ipts,ivm,ivegnew,icarbon) = veget_max(ipts,ivm)
1305
1306	! Everywhere deals with the migration of vegetation. Copy the
1307	! status of the most migrated vegetation for the whole PFT
1308	everywhere(ipts,iyoung) = MAX(everywhere(ipts,iyoung), &
1309	new_everywhere(ipts,ivm))
1310
1311	! The new soil&litter pools are the weighted mean of the newly
1312	! established vegetation for that PFT and the vegetation that
1313	! already exists in that PFT. Notice that we do not use the
1314	! "bank" concept present in LCC because we are replanting
1315	! the same PFT which already existed, just in a younger class.
1316	litter(ipts,:,iyoung,:,:) = share_young * litter(ipts,:,iyoung,:,:) + &
1317	(un - share_young) * litter(ipts,:,ivm,:,:)
1318	carbon(ipts,:,iyoung) = share_young * carbon(ipts,:,iyoung) + &
1319	(un - share_young) * carbon(ipts,:,ivm)
1320	DO ilev=1,nlevs
1321	lignin_struc(ipts,iyoung,ilev) = &
1322	litter_weight_young(istructural,ilev) * &
1323	lignin_struc(ipts,iyoung,ilev) + &
1324	(un - litter_weight_young(istructural,ilev)) * &
1325	lignin_struc(ipts,ivm,ilev)
1326	lignin_wood(ipts,iyoung,ilev) = &
1327	litter_weight_young(iwoody,ilev) * &
1328	lignin_wood(ipts,iyoung,ilev) + &
1329	(un - litter_weight_young(iwoody,ilev) ) * &
1330	lignin_wood(ipts,ivm,ilev)
1331	ENDDO
1332	bm_to_litter(ipts,iyoung,:,:) = share_young * &
1333	bm_to_litter(ipts,iyoung,:,:) + &
1334	(un - share_young) * bm_to_litter(ipts,ivm,:,:)
1335
1336	! Now move the veget_max from the current class to the young one.
1337	veget_max(ipts,iyoung)=veget_max(ipts,iyoung)+veget_max(ipts,ivm)
1338
1339	ENDIF
1340
1341	! Reset all these variables, since the PFT is dead.
1342	PFTpresent(ipts,ivm) = .FALSE.
1343	veget_max(ipts,ivm) = zero
1344	lm_lastyearmax(ipts,ivm) = zero
1345	age(ipts,ivm) = zero
1346	leaf_frac(ipts,ivm,:) = zero
1347	co2_to_bm(ipts,ivm) = zero
1348	everywhere(ipts,ivm) = zero
1349	litter(ipts,:,ivm,:,:) = zero
1350	carbon(ipts,:,ivm) = zero
1351	lignin_struc(ipts,ivm,:) = zero
1352	lignin_wood(ipts,ivm,:) = zero
1353	bm_to_litter(ipts,ivm,:,:) = zero
1354	biomass(ipts,ivm,:,:) = zero
1355	ind(ipts,ivm) = zero
1356
1357	! Keep ::biomass and circ_class_biomass synchronized
1358	circ_class_biomass(ipts,ivm,:,:,:) = zero
1359	circ_class_n(ipts,ivm,:) = zero
1360
1361	ENDIF ! If this is the youngest age class
1362
1363	ENDIF ! If we are using age classes
1364
1365	ENDIF ! All biomass in ivm and plant the same species
1366
1367	ENDIF ! is_tree
1368
1369	ENDDO ! loop over pfts
1370
1371	END DO ! loop over land points
1372
1373	!! 3. Check mass balance closure
1374
1375	! Consider this whole routine as a black box with incoming and outgoing
1376	! fluxes and a change in the mass of the box. Express in absolute
1377	! units gC (or gN), hence, multiply with dt and veget_max. In most routines
1378	! veget_max does not change and could be omitted but a general approach
1379	! was prefered.
1380
1381	!! 3.1 Calculate pools at the end of the routine
1382	pool_end = zero
1383	DO iele = 1,nelements
1384	DO ipar = 1,nparts
1385	pool_end(:,:,iele) = pool_end(:,:,iele) + &
1386	(biomass(:,:,ipar,iele) * veget_max(:,:))
1387	ENDDO
1388
1389	! Litter pool (gC m-2) * (m2 m-2)
1390	DO ilit = 1,nlitt
1391	DO ilev = 1,nlevs
1392	pool_end(:,:,iele) = pool_end(:,:,iele) + &
1393	litter(:,ilit,:,ilev,iele) * veget_max(:,:)
1394	ENDDO
1395	ENDDO
1396
1397	! Soil carbon (gC m-2) * (m2 m-2)
1398	DO icarb = 1,ncarb
1399	pool_end(:,:,iele) = pool_end(:,:,iele) + &
1400	carbon(:,icarb,:) * veget_max(:,:)
1401	ENDDO
1402	ENDDO
1403
1404	!! 4.2 Calculate components of the mass balance
1405	check_intern(:,:,iatm2land,icarbon) = orphan_flux_local(:,:,ico2bm,icarbon) * &
1406	orphan_flux_local(:,:,ivegnew,icarbon)
1407	check_intern(:,:,iland2atm,icarbon) = -un * zero
1408	check_intern(:,:,ilat2out,icarbon) = zero
1409	check_intern(:,:,ilat2in,icarbon) = -un * zero
1410	check_intern(:,:,ipoolchange,icarbon) = -un * &
1411	(pool_end(:,:,icarbon) - pool_start(:,:,icarbon))
1412	closure_intern = zero
1413	DO imbc = 1,nmbcomp
1414	closure_intern(:,:,icarbon) = closure_intern(:,:,icarbon) + &
1415	check_intern(:,:,imbc,icarbon)
1416	ENDDO
1417
1418	! Write outcome. Need to some over all PFTs here because vegetation can
1419	! move between them when an age class dies
1420	DO ipts=1,npts
1421	IF(ABS(SUM(closure_intern(ipts,:,icarbon))) .LE. min_stomate)THEN
1422	IF (ld_massbal) WRITE(numout,*) 'Mass balance closure in gap_clean'
1423	ELSE
1424	WRITE(numout,*) 'Error: mass balance is not closed in gap_clean',&
1425	SUM(closure_intern(ipts,:,icarbon))
1426	DO ivm=1,nvm
1427	WRITE(numout,*) ' ipts,ivm; ', ipts,ivm
1428	WRITE(numout,*) ' Difference is, ', closure_intern(ipts,ivm,icarbon)
1429	WRITE(numout,*) ' pool_end,pool_start: ', &
1430	pool_end(ipts,ivm,icarbon), pool_start(ipts,ivm,icarbon)
1431	WRITE(numout,*) ' orphan_fluxes: ', &
1432	orphan_flux_local(ipts,ivm,ico2bm,icarbon),&
1433	orphan_flux_local(ipts,ivm,ivegnew,icarbon)
1434	WRITE(numout,*) ' veget_max: ', veget_max(ipts,ivm)
1435	ENDDO
1436	IF(ld_stop)THEN
1437	CALL ipslerr_p (3,'gap_clean', 'Mass balance error.','','')
1438	ENDIF
1439	ENDIF
1440	ENDDO
1441
1442	! Check preservation of surface area
1443	! Needs to be further tested. The implementation was started
1444	! because it looked like we had a bug here but half way it was
1445	! realised that the bug had to be somewhere else (it was fixed).
1446	! CALL check_area('gap_clean', npts, veget_max_begin, veget_max)
1447
1448	IF (bavard.GE.2) WRITE(numout,*) 'Leaving gap_clean'
1449
1450	END SUBROUTINE gap_clean
1451
1452	END MODULE stomate_kill

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

Download in other formats: