Context Navigation

source: tags/ORCHIDEE_4_1/ORCHIDEE/src_stomate/sapiens_forestry.f90 @ 7761

Last change on this file since 7761 was 7555, checked in by sebastiaan.luyssaert, 2 years ago
Contributes to ticket #837. Simplified forest management for young stands a bit
File size: 210.6 KB

Line
1	! =================================================================================================================================
2	! MODULE : sapiens_forestry
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 Gathers the main elements for forest management: the "sapiens_forestry_main"
10	!! subroutine, which itself calls a set of subroutines
11	!! sapiens_forestry_clear, clearcut, thinning, harvest, force_load, QsortC, and
12	!! Partition, and a set of functions used in these subroutines.
13	!!
14	!!\n DESCRIPTION: None
15	!!
16	!! RECENT CHANGE(S): None
17	!!
18	!! REFERENCE(S) :
19	!! - Asael, S., 1999. Typologie des peuplements forestiers du massif vosgiens.
20	!! C.R.P.F. Lorraine-Alsace, Nancy, 54 p.
21	!! - Bellassen, V., Le Maire, G., Dhote, J.F., Viovy, N., Ciais, P., 2010.
22	!! Modeling forest management within a global vegetation model â Part 1:
23	!! model structure aSnd general behaviour. Ecological Modelling 221, 2458â2474.
24	!! - Bellassen, V., Le Maire, G., Guin, O., Dhote, J.F., Viovy, N., Ciais, P.,
25	!! 2011a. Modeling forest management within a global vegetation model â Part 2:
26	!! model validation from tree to continental scale. Ecological Modelling 222,
27	!! 57â75.
28	!! - Bellassen, V., Viovy, N., Luyssaert, S., Le Maire, G., Schelhaas, M.J.,
29	!! Ciais, P., 2011b. Reconstruction and attribution of the carbon sink of
30	!! European forests between 1950 and 2000. Global Change Biology 17, 3274-3292.
31	!! - Bottcher, H., Kurz, W.A., Freibauer, A., 2008. Accounting of forest carbon
32	!! sinks and sources under a future climate protocol-factoring out past
33	!! disturbance and management effects on age-class structure. Environmental
34	!! Science & Policy 11, 669-686.
35	!! - Cazin, A., Vallet, P., Dhote, J.F., 2003. Propositions LERFoB pour CARBOFOR
36	!! Atelier modÃ©lisation.
37	!! - Condes, S., Sterba, H., 2005. Derivation of compatible crown width
38	!! equations for some important tree species of Spain. Forest Ecology and
39	!! Management 217, 203-218.
40	!! - Deleuze, C., Pain, O., Dhote, J.F., Herve, J.C., 2004. A flexible radial
41	!! increment model for individual trees in pure even-aged stands. Annals of
42	!! Forest Science 61, 327-335.
43	!! - DhÃŽte, J.-F., HervÃ©, J.-C., 2000. Changements de productivitÃ© dans quatre
44	!! forÃªts de chÃªnes sessiles depuis 1930 : une approche au niveau du peuplement.
45	!! Ann. For. Sci. 57, 651-680.
46	!! - DhÃŽte, J.F., 1999. CompÃ©tition entre classes sociales chez le chÃªne sessile
47	!! et le hÃªtre. Revue ForestiÃšre FranÃ§aise, 309-325.
48	!! - DhÃŽte, J.F., Le MoguÃ©dec, G., 2003. PrÃ©sentation du modÃšle FagacÃ©es.
49	!! INRA, 31 p.
50	!! - IFN, 2008. Raw inventory data - 2005-2007 inventory campaigns, www.ifn.fr.
51	!! - Kolari, P., Pumpanen, J., Rannik, U., Ilvesniemi, H., Hari, P., Berninger,
52	!! F., 2004. Carbon balance of different aged Scots pine forests in Southern
53	!! Finland. Global Change Biology 10, 1106-1119.
54	!! - Lanier, L., 1994. PrÃ©cis de sylviculture. Ecole Nationale du GÃ©nie Rural,
55	!! des Eaux et des ForÃªts (ENGREF), Nancy, 477 p.
56	!! - Law, B.E., Sun, O.J., Campbell, J., Van Tuyl, S., Thornton, P.E., 2003.
57	!! Changes in carbon storage and fluxes in a chronosequence of ponderosa pine.
58	!! Global Change Biology 9, 510-524.
59	!! - Liberloo, M., Calfapietra, C., Lukac, M., Godbold, D., Luos, Z.B., Polle,
60	!! A., Hoosbeek, M.R., Kull, O., Marek, M., Raines, C., Rubino, M., Taylor, G.,
61	!! Scarascia-Mugnozza, G., Ceulemans, R., 2006. Woody biomass production during
62	!! the second rotation of a bio-energy Populus plantation increases in a future
63	!! high CO2 world. Global Change Biology 12, 1094-1106.
64	!! - Liberloo, M., Luyssaert, S., Bellassen, V., Djomo, S.N., Lukac, M.,
65	!! Calfapietra, C., Janssens, I., Hoosbeek, M.R., Viovy, N., Churkina, G.,
66	!! Scarascia-Mugnozza, G., Ceulemans, R., 2010. Bio-energy retains its
67	!! mitigation potential under elevated CO2. PLOS-one 5, e1 1648.
68	!! - Litton, C.M., Raich, J.W., Ryan, M.G., 2007. Carbon allocation in forest
69	!! ecosystems. Global Change Biology 13, 2089-2109.
70	!! - Luyssaert, S., Inglima, I., Jung, M., Richardson, A.D., Reichsteins, M.,
71	!! Papale, D., Piao, S.L., Schulzes, E.D., Wingate, L., Matteucci, G., Aragao,
72	!! L., Aubinet, M., Beers, C., Bernhoffer, C., Black, K.G., Bonal, D.,
73	!! Bonnefond, J.M., Chambers, J., Ciais, P., Cook, B., Davis, K.J., Dolman,
74	!! A.J., Gielen, B., Goulden, M., Grace, J., Granier, A., Grelle, A., Griffis,
75	!! T., Grunwald, T., Guidolotti, G., Hanson, P.J., Harding, R., Hollinger, D.
76	!! Y., Hutyra, L.R., Kolar, P., Kruijt, B., Kutsch, W., Lagergren, F.,
77	!! Laurila, T., Law, B.E., Le Maire, G., Lindroth, A., Loustau, D., Malhi, Y.,
78	!! Mateus, J., Migliavacca, M., Misson, L., Montagnani, L., Moncrieff, J.,
79	!! Moors, E., Munger, J.W., Nikinmaa, E., Ollinger, S.V., Pita, G.,
80	!! Rebmann, C., Roupsard, O., Saigusa, N., Sanz, M.J., Seufert, G., Sierra, C.,
81	!! Smith, M.L., Tang, J., Valentini, R., Vesala, T., Janssens, I.A., 2007. CO2
82	!! balance of boreal, temperate, and tropical forests derived from a global
83	!! database. Global Change Biology 13, 2509-2537.
84	!! - Mokany, K., Raison, R.J., Prokushkin, A.S., 2006. Critical analysis of
85	!! root: shoot ratios in terrestrial biomes. Global Change Biology 12, 84-96.
86	!! - Mund, M., Kummetz, E., Hein, M., Bauer, G.A., Schulze, E.D., 2002. Growth
87	!! and carbon stocks of a spruce forest chronosequence in central Europe.
88	!! Forest Ecology and Management 171, 275-296.
89	!! - Newton, R.F., Amponsah, I.G., 2007. Comparative evaluation of five height-
90	!! diameter models developed for black spruce and jack pine stand-types in
91	!! terms of goodness-of-fit, lack-of-fit and predictive ability. Forest Ecology
92	!! and Management 247, 149-166.
93	!! - Numerical Recipes, 2007. Numerical Recipes: The Art of Scientific Computing.
94	!! Cambridge University Press, 1256 p.
95	!! - Ovington, J.D., Madgwick, H.A.I., 1957. Afforestation and soil reaction.
96	!! Journal of Soil Science 8, 141-149.
97	!! - Pontailler, J.Y., Ceulemans, R., Guittet, J., 1999. Biomass yield of poplar
98	!! after five 2-year coppice rotations. Forestry 72, 157-163.
99	!! - Pretzsch, H., Biber, P., DurskÃœ, J., 2002. The single tree-based stand
100	!! simulator SILVA: construction, application and evaluation. Forest Ecology and
101	!! Management 162, 3.
102	!! - Reineke, L.H., 1933. Perfecting a stand-density index for even-aged forests.
103	!! Journal of Agricultural Research 46, 627-638.
104	!! - Turner, D.P., Ritts, W.D., Cohen, W.B., Maeirsperger, T.K., Gower, S.T.,
105	!! Kirschbaum, A.A., Running, S.W., Zhao, M.S., Wofsy, S.C., Dunn, A.L., Law, B.E.,
106	!! Campbell, J.L., Oechel, W.C., Kwon, H.J., Meyers, T.P., Small, E.E., Kurc, S.A.,
107	!! Gamon, J.A., 2005. Site-level evaluation of satellite-based global terrestrial
108	!! gross primary production and net primary production monitoring. Global Change
109	!! Biology 11, 666-684.
110	!! - Vacchiano, G., Motta, R., Long, J.N., Shaw, J.D., 2008. A density management
111	!! diagram for Scots pine (Pinus sylvestris L.): A tool for assessing the forest's
112	!! protective effect. Forest Ecology and Management 255, 2542-2554.
113	!! - Vieira, I.C.G., de Almeida, A.S., Davidson, E.A., Stone, T.A., de Carvalho,
114	!! C.J.R., Guerrero, J.B., 2003. Classifying successional forests using Landsat
115	!! spectral properties and ecological characteristics in eastern Amazonia. Remote
116	!! Sensing of Environment 87, 470-481.
117	!! - Zianis, D., Mencuccini, M., 2004. On simplifying allometric analyses of forest
118	!! biomass. Forest Ecology and Management 187, 311-332.
119	!!
120	!! SVN :
121	!! $HeadURL: $
122	!! $Date: $
123	!! $Revision: $
124	!! \n
125	!_ ================================================================================================================================
126
127	MODULE sapiens_forestry
128
129	! modules used:
130	#ifdef __NAGFOR
131	USE,INTRINSIC :: IEEE_ARITHMETIC
132	#endif
133
134	USE netcdf
135	USE ioipsl_para
136	USE constantes
137	USE grid
138	USE pft_parameters
139	USE function_library, ONLY: wood_to_volume, wood_to_height, &
140	wood_to_circ, Nmax, &
141	wood_to_dia, sort_circ_class_biomass, &
142	wood_to_qmdia, calculate_rdi, &
143	check_vegetation_area, check_mass_balance
144	USE interpol_help, ONLY: aggregate_p
145	USE stomate_data
146
147	IMPLICIT NONE
148
149	! private & public routines
150
151	PRIVATE
152
153	PUBLIC sapiens_forestry_main, sapiens_forestry_clear, &
154	sapiens_forestry_read_fm, sapiens_forestry_read_species_change, &
155	sapiens_forestry_read_desired_fm, sapiens_forestry_read_litter, &
156	sapiens_forestry_litter_raking, sapiens_forestry_species_change, &
157	sapiens_forestry_flag_species_change, sapiens_forestry_read_spinup_clearcut
158
159	LOGICAL, SAVE :: firstcall_sapiens_forestry = .TRUE. !! first call
160	!$OMP THREADPRIVATE(firstcall_sapiens_forestry)
161	INTEGER(i_std), SAVE :: printlev_loc !! Local level of text output for current module
162	!$OMP THREADPRIVATE(printlev_loc)
163
164	CONTAINS
165
166	!! ================================================================================================================================
167	!! SUBROUTINE : sapiens_forestry_clear
168	!!
169	!>\BRIEF Set the flag ::firstcall_sapiens_forestry to .TRUE. and as such activate section 2 of the subroutine forestry (see below).
170	!!
171	!_ ================================================================================================================================
172
173	SUBROUTINE sapiens_forestry_clear
174	firstcall_sapiens_forestry = .TRUE.
175	END SUBROUTINE sapiens_forestry_clear
176
177
178	!! ================================================================================================================================
179	!! SUBROUTINE : sapiens_forestry_main
180	!!
181	!>\BRIEF This subroutine is the core of the forest management module. It
182	!! has been modified so it only decides if trees are getting cut or not.
183	!!
184	!! DESCRIPTION : This forest management module has been changed quite a bit
185	!! from Valentin's original version (documented in Ballassen et al (2010). The
186	!! biggest reason for these changes is that it is meant to be used with the new
187	!! functional allocation scheme. This scheme already includes circumference
188	!! classes. In each circumference class there are circ_class_n model trees
189	!! which are all identical. Therefore, we don't need to create a picture
190	!! of every tree on the stand since in reality we only have ncirc trees, each
191	!! of which are replicated.
192	!!
193	!! In addition to this, several other processes that were done here were moved
194	!! to other parts of the code. For example, trees are just scheduled to die
195	!! here; the actually killing of the biomass and moving the wood to the harvest
196	!! and litter pools is done in sapiens_kill.f90. All tree establishment after
197	!! clearcutting is done in stomate_prescribe.f90. Self-thinning is a completely
198	!! natural process, and therefore done with environmental mortality in
199	!! stomate_mark_kill.f90.
200
201	!!++++ CHECK REFERENCE ++++
202	!! Liberloo et al. (2010) describes how short
203	!! rotation coppices are simulated and validated.\n
204	!!+++++++++++++++++++++++
205
206	!! The forest management flags have changed as well.
207	!! FM = 1 : No human intervention (ORCHIDEE default)
208	!! 2 : Thinnings based on the RDI, clearcuts based on tree density,
209	!! annual increment, and tree diameter. Thinnings from above and
210	!! from below are determined by the sign of thstrat.
211	!! 3 : Coppices
212	!! 4 : Short rotation coppices
213	!!
214	!! Relative density index (RDI) target of human thinning at the beginning of
215	!! the rotation (unitless). This is the value aimed at by human thinning at
216	!! the beginning of the rotation. It is necessarily lower than 1. The actual
217	!! RDI varies around this target, + or - delta_rdi, when forest_managed = 2.
218	!! When forest_managed = 3, human thinning occurs every year. For FM = 4,
219	!! the thinnings occur as a function of stand age.
220	!!
221	!! RECENT CHANGE(S) : None
222	!!
223	!! MAIN OUTPUT VARIABLE(S): ::circ_class_kill
224	!!
225	!! REFERENCE(S) : See above, module description.
226	!!
227	!! FLOWCHART :
228	!! \n
229	!_ ================================================================================================================================
230
231	SUBROUTINE sapiens_forestry_main (npts, age_stand, last_cut, &
232	circ_class_n, circ_class_kill, forest_managed, spinup_clearcut, &
233	circ_class_biomass, mai, pai, previous_wood_volume, &
234	mai_count, coppice_dens, veget_max, fm_change_map, &
235	species_change_map)
236
237	!! 0. Variable and parameter declaration
238
239	!! 0.1 Input variables
240
241	INTEGER(i_std), INTENT(in) :: npts !! Domain size - number of pixels
242	!! (dimensionless)
243	REAL(r_std), DIMENSION(:,:), INTENT(in) :: veget_max !! "maximal" coverage fraction of a PFT on
244	!! the ground (unitless, 0-1)
245	INTEGER(i_std), DIMENSION(:,:),INTENT(in) :: last_cut !! Years since last thinning (years)
246	INTEGER(i_std), DIMENSION(:,:),INTENT(in) :: age_stand !! Age of stand (years)
247	INTEGER(i_std), DIMENSION(:,:), INTENT(in) :: fm_change_map !! A map which gives the desired FM strategy when
248	!! the PFT will be replanted after a clearcut.
249	!! (1-nvm,unitless)
250	INTEGER(i_std), DIMENSION(:,:), INTENT(in) :: species_change_map !! A map which gives the PFT number that each
251	!! PFT will be replanted as in case of a clearcut.
252	!! (1-nvm,unitless)
253
254	!! 0.2 Output
255
256	!! 0.3 Modified fields
257	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(inout) :: circ_class_biomass !! Biomass components of the model tree
258	!! within a circumference class
259	!! class @tex $(g C ind^{-1})$ @endtex
260	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: circ_class_n !! Number of trees in each circumference
261	!! class (0-2.2+ and poles-large wood, see
262	!! part 12.)
263	INTEGER(i_std), DIMENSION(:,:), INTENT(inout) :: forest_managed !! Forest management flag: 0 = orchidee
264	!! standard, 1= self-thinning only, 2=
265	!! high-stand, 3= high-stand smoothed, 4=
266	!! coppices
267	INTEGER(i_std), DIMENSION (:,:), INTENT(in) :: spinup_clearcut !! Map to indicate clearcut event during spinup
268	!! for a given PFT and pixel (zero = no clearcut; one = clearcut).
269	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(inout) :: circ_class_kill !! Number of trees within a circ that needs
270	!! to be killed @tex $(ind m^{-2})$ @endtex
271	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: mai !! The mean annual increment
272	!! @tex $(m3 / m2 / year)$ @endtex
273	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: pai !! The period annual increment
274	!! @tex $(m3 / m2 / year)$ @endtex
275	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: previous_wood_volume!! The volume of the tree trunks
276	!! in a stand for the previous year.
277	!! @tex $(m3 / m2 )$ @endtex
278	INTEGER(i_std), DIMENSION(:,:),INTENT(inout) :: mai_count !! The number of times we've
279	!! calculated the volume increment
280	!! for a stand
281	REAL(r_std), DIMENSION(:,:),INTENT(inout) :: coppice_dens !! The density of a coppice at the first
282	!! cutting.
283	!! @tex $( 1 / m**2 )$ @endtex
284
285	!! 0.4 Local variables
286
287	INTEGER(i_std) :: icir,ivm,ipts !! Indexes
288	INTEGER(i_std) :: imbc, iele, ipar !! Indexes
289	REAL(r_std) :: trees,total_trees !! a number of trees
290	REAL(r_std), DIMENSION(ncirc) :: diameters_temp !! Trunk diameter calculated from allometry [m]
291	REAL(r_std), DIMENSION(ncirc) :: circ_temp !! Trunk circumference calculated from allometry
292	REAL(r_std), DIMENSION(npts,nvm,ncirc) :: height,dia,cn !! Tree height calculated from allometric
293	!! relationships (m)
294	REAL(r_std) :: ave_tree_height
295	REAL(r_std) :: woody_biomass
296	REAL(r_std) :: current_wood_volume
297	REAL(r_std) :: current_increment
298	REAL(r_std) :: ave_tree_dia
299	INTEGER(i_std) :: nrotations !! The number of rotations that we've done
300	!! for SRC minus one
301	REAL(r_std) :: qm_dia !! quadratic mean diameter of the forest (m)
302	INTEGER(i_std) :: ierr !! error message
303	INTEGER(i_std) :: igroup_new !! Species group to which the species with which
304	!! the stand will be replanted belongs (unitless)
305	INTEGER(i_std) :: igroup_old !! Species group to which the current species
306	!! belongs (unitless)
307	REAL(r_std), DIMENSION(ncirc) :: circ_class_n_new !! variable used to sort circ_class_n
308	REAL(r_std), DIMENSION(ncirc,nparts,nelements) :: circ_class_biomass_new !! variable used to sort circ_class_biomass
309	REAL(r_std), DIMENSION(npts,nvm) :: rdi !! Relative density index (unitless, 0-2)
310	REAL(r_std), DIMENSION(npts,nvm) :: rdi_target_upper !! Upper limit of RDI. When reached the stand
311	!! needs to be thinned (unitless)
312	REAL(r_std), DIMENSION(npts,nvm) :: rdi_target_lower !! Lower limit of RDI. When thinned, thin to this
313	!! stand density (unitless)
314	REAL(r_std), DIMENSION(npts,nvm,nmbcomp,nelements) :: check_intern !! Contains the components of the internal
315	!! mass balance chech for this routine
316	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
317	REAL(r_std), DIMENSION(npts,nvm,nelements) :: closure_intern !! Check closure of internal mass balance
318	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
319	REAL(r_std), DIMENSION(npts,nvm,nelements) :: pool_start !! Start and end pool of this routine
320	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
321	REAL(r_std), DIMENSION(npts,nvm,nelements) :: pool_end !! Start and end pool of this routine
322	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
323	!_ ===============================================================================================================================
324
325	!! 1. Initialize
326
327	IF (firstcall_sapiens_forestry) THEN
328	! Initialize local printlev
329	printlev_loc=get_printlev('sapiens_forestry')
330
331	firstcall_sapiens_forestry=.FALSE.
332	END IF
333
334	IF ( printlev.GE.2) WRITE(numout,*) 'Entering sapiens_forestry_main'
335
336	!! 1.2 Initialize check for mass balance closure
337	! We don't actually have to do this in sapiens_forestry_main, since none of the
338	! pools are changed even touched. The sole purpose of this routine
339	! is to schedule trees for killing which means that only
340	! ::circ_class_kill is changed. Biomass is sorted so we will check
341	! wether nothing goes wrong during sorting
342	IF (err_act.GT.1) THEN
343
344	check_intern(:,:,:,:) = zero
345	pool_start(:,:,:) = zero
346	DO iele = 1,nelements
347
348	DO ipar = 1,nparts
349
350	DO icir = 1,ncirc
351
352	! Initial biomass pool. Use circ_class_biomass
353	! because that variable is the prognostic variable
354	! biomass is syncronized to circ_class_biomass. If
355	! many diameter classes are used, each diameter class
356	! separatly passes all criteria but when combined
357	! a mass balance problem occurs.
358	pool_start(:,:,iele) = pool_start(:,:,iele) + &
359	circ_class_biomass(:,:,icir,ipar,iele) * &
360	circ_class_n(:,:,icir) * veget_max(:,:)
361
362	ENDDO
363
364	ENDDO
365
366	ENDDO
367
368	ENDIF ! err_act.GT.1
369
370	!! 1.3 Initialize check for surface area conservation
371	! Veget_max is a INTENT(in) variable and can therefore
372	! not be changed during the course of this subroutine
373	! No need to check whether the subroutine preserves the
374	! total surface area of the pixel.
375
376	!! 1.4 Sort the diameter classes
377	! A small number of individuals and distributing small labile and
378	! carbres pools from the stand level (at which they are calculated) to
379	! the tree level (at which they are stored) can result in precision
380	! errors. Sometimes these errors are sufficient to overturn the ascending
381	! order of tree biomass that is expected in circ_class_biomass.
382	! sapiens_forestry_main is called after this routine and relies on the
383	! assumption that the biomass is increasing. When needed, sort biomass
384	! to make sure this assumption is satisfied. Sorting is enforced at
385	! several places in the code
386	DO ipts = 1,npts
387	DO ivm = 1,nvm
388	IF (is_tree(ivm)) THEN
389	circ_class_biomass_new = circ_class_biomass(ipts,ivm,:,:,:)
390	circ_class_n_new = circ_class_n(ipts,ivm,:)
391	CALL sort_circ_class_biomass(circ_class_biomass_new, &
392	circ_class_n_new)
393	circ_class_biomass(ipts,ivm,:,:,:) = circ_class_biomass_new
394	circ_class_n(ipts,ivm,:) = circ_class_n_new
395	ENDIF
396	ENDDO
397	ENDDO
398
399
400	!! 2. Check and apply sapiens_forestry_main
401
402	! for every forest stand in the model, we need to decide
403	! if we are thinning or clearcutting.
404	DO ipts=1,npts
405
406	pft:DO ivm=1,nvm
407
408	! Forest management can only be applied to forests
409	IF(.NOT. is_tree(ivm)) CYCLE
410
411	! This PFT is not present, so it can't be managed
412	IF(veget_max(ipts,ivm) == zero) CYCLE
413
414	! There is no biomass, so nothing to manage
415	IF(SUM(SUM(circ_class_biomass(ipts,ivm,:,:,icarbon),1),1) .LE. min_stomate) CYCLE
416
417	! For unmanaged forests, either we apply clearcut during spinup, or
418	! simply do nothing.
419	IF ( forest_managed(ipts,ivm) .EQ. ifm_none ) THEN
420
421	IF ( (spinup_clearcut(ipts,ivm) .EQ. flag_spinup_clearcut) .OR. forced_clear_cut ) THEN
422	CALL clearcut_harvest(circ_class_n(ipts,ivm,:), &
423	circ_class_kill(ipts,ivm,:,ifm_none,icut_clear))
424	ELSE
425	! Otherwise we just get out of the loop
426	CYCLE
427	ENDIF
428
429	ELSE
430
431	! Managed pft, apply to appropriate management
432	IF(printlev_loc>=4)THEN
433	WRITE(numout,*) 'Sapiens_forestry_main: check whether the stand ',&
434	'needs management', ivm
435	CALL flush(numout)
436	ENDIF
437
438	! Convert the aboveground woody biomass into a wood volume,
439	! removing the influence of the branches.
440	current_wood_volume = wood_to_volume(npts,&
441	circ_class_biomass(ipts,ivm,:,:,:),&
442	circ_class_n(ipts,ivm,:),ivm,branch_ratio(ivm),0)
443
444	! The mean annual increment is tricky. Need to calculate it
445	! for the volume, since that is what foresters use and the NPP
446	! would include leaf NPP. However, we cannot calculate a
447	! volume increment immediately after thinning. Instead, the first
448	! year after a thinning or after a clearcut (for consistency)
449	! we just save the current volume.
450	IF(last_cut(ipts,ivm) .EQ. 1)THEN
451
452	! do nothing. We'll save the current wood volume later.
453
454	ELSEIF(last_cut(ipts,ivm) .GT. 1)THEN
455
456	current_increment=&
457	current_wood_volume-previous_wood_volume(ipts,ivm)
458	mai_count(ipts,ivm)=mai_count(ipts,ivm)+1
459	mai(ipts,ivm)=mai(ipts,ivm)*REAL(mai_count(ipts,ivm)-1)/&
460	REAL(mai_count(ipts,ivm))+current_increment/&
461	REAL(mai_count(ipts,ivm))
462
463	! The periodic annual increment is an average over the past few
464	! years. Notice that the first year after a cut that we can
465	! compute an increment for is year 2 (the difference between
466	! year 2 and year 1).
467	IF(last_cut(ipts,ivm) .LT. (n_pai+1))THEN
468
469	pai(ipts,ivm)=pai(ipts,ivm)*REAL(last_cut(ipts,ivm)-2)/&
470	REAL(last_cut(ipts,ivm)-1)+&
471	current_increment/REAL(last_cut(ipts,ivm)-1)
472
473	ELSE
474
475	pai(ipts,ivm)=pai(ipts,ivm)*REAL(n_pai-1)/REAL(n_pai)+&
476	(current_wood_volume-previous_wood_volume(ipts,ivm))&
477	/REAL(n_pai)
478
479	ENDIF
480
481	ELSE
482
483	! Cannot think of a scenario where last_cut is zero.
484	! If it is zero, we thinned or clearcut, but that happens
485	! after this routine and we always increment it before this
486	! routine.
487	WRITE(numout,*) 'ERROR: Sapiens_forestry_main, why is last_cut equal to zero?'
488	WRITE(numout,*) 'ipts, ivm, last_cut(ipts,ivm) ',&
489	ipts, ivm, last_cut(ipts,ivm)
490	IF (err_act.GT.1) CALL ipslerr_p (3, &
491	'sapiens_forestry_main', 'last_cut is zero.',&
492	'This should not happen.',&
493	'Look in the output file for ERROR.')
494	ENDIF
495
496	! We also need to store the wood volume of this year in order to
497	! calculate the PAI for next year.
498	previous_wood_volume(ipts,ivm)=current_wood_volume
499
500	! The following is specific for each forest management strategy
501	IF(forest_managed(ipts,ivm).EQ.ifm_thin)THEN
502
503	! Forests can only be thinned if it is profitable for the
504	! managers to do so. This profitability is based on the
505	! average height of the 100 tallest trees per hectare.
506	! If it's profitable, the stand is thinned to a target
507	! relative density index (RDI), which is the ratio of the
508	! true density to the self thinning density. This target
509	! depends on the management type, but it is always less
510	! than unity. Table 1 in Bellassen et al, Ecological
511	! Modeling (2010) gives values for some of these parameters.
512
513	! First, calculate the average height of the 100 tallest
514	! trees on the hectare. Remember that all our trees are the
515	! same height in a given circ class. Since we don't explicitly
516	! track trees, we just start at the largest circ classs and
517	! add trees to our average until we have 100 (Note that this
518	! number has been externalised as ::ntress_profit)
519	height(ipts,ivm,:)=&
520	wood_to_height(circ_class_biomass(ipts,ivm,:,:,icarbon),&
521	ivm)
522
523	! Debug
524	IF(printlev_loc>=4 .AND. ipts == test_grid .AND. &
525	ivm == test_pft)THEN
526	WRITE(numout,*) 'Sapiens_forestry_main, check for thinning', ivm
527	! just for printing out some information
528	dia(ipts,ivm,:)= &
529	wood_to_dia(circ_class_biomass(ipts,ivm,:,:,icarbon),ivm)
530	DO icir=1,ncirc
531	WRITE(numout,) circ_class_n(ipts,ivm,icir)m2_to_ha,&
532	height(ipts,ivm,icir), &
533	deux*SQRT(cn(ipts,ivm,icir)/pi), dia(ipts,ivm,icir)
534	ENDDO
535	CALL flush(numout)
536	ENDIF
537
538	!-
539	! condition introduced in order to cut the stand at the end
540	! of a spinup run. Put forced_clear_cut=n to avoid
541	! any forced clear-cut
542	IF( forced_clear_cut ) THEN
543	WRITE(numout,*) 'A clear-cut at the begenning of &
544	the simulation is scheduled'
545	CALL clearcut_harvest(circ_class_n(ipts,ivm,:), &
546	circ_class_kill(ipts,ivm,:,ifm_thin,icut_clear))
547	CYCLE pft
548	ENDIF
549
550	! Select the 100 tallest treesb and calculate their
551	! average height
552	total_trees=zero
553	ave_tree_height=zero
554	tree:DO icir=ncirc,1,-1
555	trees=circ_class_n(ipts,ivm,icir)*m2_to_ha
556	IF(trees .GE. ntrees_profit)THEN
557	ave_tree_height=height(ipts,ivm,icir)
558	EXIT tree
559	ELSEIF((trees + total_trees) .GE. ntrees_profit)THEN
560	ave_tree_height=ave_tree_height*total_trees/&
561	REAL(ntrees_profit,r_std)+&
562	(REAL(ntrees_profit,r_std)-total_trees)/&
563	REAL(ntrees_profit,r_std)*&
564	height(ipts,ivm,icir)
565	EXIT tree
566	ELSE
567	ave_tree_height=ave_tree_height*total_trees/&
568	(total_trees+trees)+&
569	trees/(total_trees+trees)*height(ipts,ivm,icir)
570	total_trees=total_trees+trees
571	ENDIF
572
573	ENDDO tree
574
575	! in order to thin, we need to know the relative
576	! density index (RDI) of our stand. This is the
577	! ratio of the current stand density to
578	! the self-thinning density, Nmax.
579	diameters_temp(:)= &
580	wood_to_dia(circ_class_biomass(ipts,ivm,:,:,icarbon),ivm)
581
582	! with these diameters, we can check to see if we
583	! want to clearcut. If the diameters get too big,
584	! harvesting equipment can't handle the trees anymore.
585	! On the other hand, we need to make sure there
586	! are enough big trees that it is worthwhile to cut
587	! them down.
588	total_trees=zero
589	DO icir=1,ncirc
590	IF(diameters_temp(icir) .GT. largest_tree_dia(ivm)) THEN
591	total_trees=total_trees+circ_class_n(ipts,ivm,icir)*&
592	m2_to_ha
593	ENDIF
594	ENDDO
595
596	! Debug
597	IF(printlev_loc>=4 .AND. ipts == test_grid .AND. &
598	ivm == test_pft)THEN
599	WRITE(numout,*) 'Sapiens_forestry_main, testing diameter criterion', ivm
600	WRITE(numout,*) 'total_trees: ', total_trees
601	WRITE(numout,*) 'dens_target: ', dens_target(ivm)
602	WRITE(numout,*) 'largest_tree_dia: ', largest_tree_dia(ivm)
603	WRITE(numout,*) 'diameters_temp: ', diameters_temp(:)
604	WRITE(numout,*) 'circ_class_n (in trees per hectare): ', &
605	circ_class_n(ipts,ivm,:)*m2_to_ha
606	CALL flush(numout)
607	ENDIF
608	!-
609
610	! CONDITION 1: there are less than 100 trees left
611	IF( SUM(circ_class_n(ipts,ivm,:))*m2_to_ha .LT. &
612	dens_target(ivm)) THEN
613
614	! These are conditions which are suitable for switching
615	! from a management to a conservation strategy.
616	! Such a switch is not internalized in the code and
617	! needs to be prescribed. The way to prescribe it is
618	! by using species/management changes and prescribing
619	! the desired management for some of all pixels
620	! as unmanaged.
621	IF(ok_change_species)THEN
622
623	! The criteria for a clearcut are met but check
624	! whether we want to switch to a conservation
625	! strategy on this site.
626	IF(fm_change_map(ipts,ivm) .EQ. ifm_none) THEN
627
628	! Change the management strategy (most likely
629	! high stand management) to conservation
630	forest_managed(ipts,ivm) = ifm_none
631
632	! Debug
633	IF(printlev_loc>=4)THEN
634	WRITE(numout,*) 'species change ',&
635	'- conservation strategy, ',&
636	'density becoming too low'
637	ENDIF
638	!-
639
640	ELSE
641
642	! The species/management change scenario does
643	! not prescribe a change to a conservation strategy
644	! but the criteria for a clearcut are met
645	CALL clearcut_harvest(circ_class_n(ipts,ivm,:), &
646	circ_class_kill(ipts,ivm,:,ifm_thin,icut_clear))
647
648	ENDIF
649
650	ELSE
651
652	! Species/management changes are not considered. The
653	! criteria for a clearcut are met
654	CALL clearcut_harvest(circ_class_n(ipts,ivm,:), &
655	circ_class_kill(ipts,ivm,:,ifm_thin,icut_clear))
656
657	ENDIF
658
659	CYCLE pft
660
661	ENDIF
662
663	! CONDITION 2: There are more than 100 trees for which
664	! the diameter exceeds the threshold, time to cut
665	IF(total_trees .GT. dens_target(ivm))THEN
666
667	! Debug
668	IF(printlev_loc>=4 .AND. ipts == test_grid .AND. &
669	ivm == test_pft)THEN
670	WRITE(numout,*) 'Sapiens_forestry_main, diameter exceeds ',&
671	'treshold, clearcut', ivm
672	WRITE(numout,*) 'total_trees: ', total_trees
673	WRITE(numout,*) 'dens_target: ', dens_target(ivm)
674	WRITE(numout,*) 'largest_tree_dia: ', &
675	largest_tree_dia(ivm)
676	WRITE(numout,*) 'diameters_temp: ', diameters_temp(:)
677	WRITE(numout,*) 'circ_class_n (trees per hectare): ', &
678	circ_class_n(ipts,ivm,:)*m2_to_ha
679	CALL flush(numout)
680	ENDIF
681	!-
682
683	! These are conditions which are suitable for switching
684	! from a management to a conservation strategy. Such a
685	! switch is not internalized in the code and needs to be
686	! prescribed. The way to prescribe it is by using
687	! species/management changes and prescribing the desired
688	! management for some of all pixels as unmanaged.
689	IF(ok_change_species)THEN
690
691	! The criteria for a clearcut are met but check whether
692	! we want to switch to a conservation strategy on this
693	! site.
694	IF( fm_change_map(ipts,ivm) .EQ. ifm_none) THEN
695
696	! Change the management strategy (most likely
697	! high stand management) to conservation
698	forest_managed(ipts,ivm) = ifm_none
699
700	! debug
701	IF(printlev_loc>=4)THEN
702	WRITE(numout,*) 'species change - ',&
703	'conservation strategy, ',&
704	'exceeded largest diameter'
705	ENDIF
706	!-
707
708	ELSE
709
710	! The species/management change scenario does
711	! not prescribe a change to a conservation strategy
712	! but the criteria for a clearcut are met
713	CALL clearcut_harvest(circ_class_n(ipts,ivm,:), &
714	circ_class_kill(ipts,ivm,:,ifm_thin,icut_clear))
715
716	ENDIF
717
718	ELSE
719
720	! Species/management changes are not considered. The
721	! criteria for a clearcut are met
722	CALL clearcut_harvest(circ_class_n(ipts,ivm,:), &
723	circ_class_kill(ipts,ivm,:,ifm_thin,icut_clear))
724
725	ENDIF
726
727	CYCLE pft
728	ENDIF
729
730	! CONTION 3: a clearcut is expected if the trees start
731	! growing too slowly. This is measured by looking at
732	! the mean annual increment, which is simply the trunk
733	! volume divided by the age of the tree. We compare
734	! this to average tree growth over the past five years
735	! (i.e. the tree volume divided by 5), called the
736	! periodic annual increment. If the PAI is smaller than
737	! the MAI, the trees are past their prime and the owner
738	! will cut them and plant anew.
739
740	! We calculate the MAI and the PAI at the same place where
741	! we increment the age of the stand, so here we just have
742	! to compare them. We need to check that we have enough
743	! years from the last cut to calculate a good PAI.
744	IF(PAI(ipts,ivm) .LT. MAI(ipts,ivm) .AND. &
745	last_cut(ipts,ivm) .GT. n_pai)THEN
746
747	! Debug
748	IF(printlev_loc>=4 .AND. ipts == test_grid .AND. &
749	ivm == test_pft)THEN
750	WRITE(numout,*) 'Sapiens_forestry_main, harvesting ',&
751	'increment PAI<MAI'
752	WRITE(numout,*) 'ipts, ivm, mai, pai: ',ipts, &
753	ivm, mai(ipts,ivm), pai(ipts,ivm)
754	ENDIF
755	!-
756
757	IF(ok_change_species)THEN
758
759	! The criteria for a clearcut are met but check
760	! whether we want to switch to a conservation
761	! strategy on this site.
762	IF( fm_change_map(ipts,ivm) .EQ. ifm_none) THEN
763
764	! Change the management strategy (could be
765	! high stand, coppice and short rotation
766	! coppice) to conservation
767	forest_managed(ipts,ivm) = ifm_none
768
769	IF(printlev_loc>=4)THEN
770	WRITE(numout,*) 'species change - ',&
771	'conservation strategy, ', &
772	'stand growth is too low'
773	ENDIF
774
775	ELSE
776
777	! The species/management change scenario does
778	! not prescribe a change to a conservation strategy
779	! but the criteria for a clearcut are met
780	CALL clearcut_harvest(circ_class_n(ipts,ivm,:), &
781	circ_class_kill(ipts,ivm,:,ifm_thin,icut_clear))
782
783	ENDIF
784
785	ELSE
786
787	! Species/management changes are not considered. The
788	! criteria for a clearcut are met
789	CALL clearcut_harvest(circ_class_n(ipts,ivm,:), &
790	circ_class_kill(ipts,ivm,:,ifm_thin,icut_clear))
791
792	ENDIF
793
794	CYCLE pft
795	ENDIF
796
797	! our rdi_target will become a function of both PFT
798	! and diameter/density Calculate rdi and rdi_target
799	! and delta_rdi based on the yield models rdi of the
800	! stand as it is
801	qm_dia = wood_to_qmdia(&
802	circ_class_biomass(ipts,ivm,:,:,icarbon), &
803	circ_class_n(ipts,ivm,:), ivm)
804
805	! The self-thinning and yield relationships are for
806	! diameters expressed in cm
807	rdi(ipts,ivm) = (SUM(circ_class_n(ipts,ivm,:))*m2_to_ha)/ &
808	Nmax(qm_dia*m_to_cm,ivm)
809
810	! rdi_target_upper (based on the yield tables)
811	! rdi depends on the observed self-thinning
812	! relationship and rdi_target_upper
813	rdi_target_upper(ipts,ivm) = MAX(MIN( a_rdi_upper_man(ivm)+ &
814	qm_diam_to_cmb_rdi_upper_man(ivm),c_rdi_upper_man(ivm)),&
815	d_rdi_upper_man(ivm))
816	rdi_target_lower(ipts,ivm) = MAX(MIN(a_rdi_lower_man(ivm)+ &
817	qm_diam_to_cmb_rdi_lower_man(ivm),c_rdi_lower_man(ivm)),&
818	d_rdi_lower_man(ivm))
819
820	! Debug
821	IF(printlev_loc>=4 .AND. ipts == test_grid .AND. &
822	ivm == test_pft)THEN
823	WRITE(numout,*) 'rdi', rdi(ipts,ivm)
824	WRITE(numout,*) 'rdi_target_upper calcul', &
825	rdi_target_upper(ipts,ivm)
826	WRITE(numout,*) 'Sapiens_forestry_main, rdi print ',rdi(ipts,ivm),&
827	rdi_target_upper(ipts,ivm),&
828	rdi_target_lower(ipts,ivm)
829	WRITE(numout,) "Nmax ",Nmax(qm_diam_to_cm,ivm)
830	WRITE(numout,*) "a_rdi_upper_man,b_rdi_upper_man ", &
831	a_rdi_upper_man(ivm),b_rdi_upper_man(ivm)
832	ENDIF
833	!-
834
835	! we only need to thin if the rdi of our stand is outside
836	! the acceptable rdi range given by rdi_target_upper and
837	! rdi_target_lower
838	IF(rdi(ipts,ivm) .GT. rdi_target_upper(ipts,ivm)) THEN
839
840	! we need the circumference of each model tree in meters
841	circ_temp(:)= wood_to_circ(&
842	circ_class_biomass(ipts,ivm,:,:,icarbon),ivm)
843
844	! Debug
845	IF(printlev_loc>=4 .AND. ipts == test_grid .AND. &
846	ivm == test_pft)THEN
847	DO icir=1,ncirc
848	WRITE(numout,*) 'Sapiens_forestry_main thinning now: ',&
849	ipts,ivm,icir
850	WRITE(numout,*) 'Sapiens_forestry_main height, ccn, ',&
851	height(ipts,ivm,icir),&
852	circ_class_n(ipts,ivm,icir)*m2_to_ha
853	ENDDO
854	WRITE(numout,*) "circ_temp ",circ_temp(:)
855	WRITE(numout,*) "NOTE: assumes ncirc = 3"
856	WRITE(numout,*) "circ_class_biomass(1) ",&
857	SUM(circ_class_biomass(ipts,ivm,1,:,icarbon))
858	WRITE(numout,*) "circ_class_biomass(2) ",&
859	SUM(circ_class_biomass(ipts,ivm,2,:,icarbon))
860	WRITE(numout,*) "circ_class_biomass(3) ",&
861	SUM(circ_class_biomass(ipts,ivm,3,:,icarbon))
862	ENDIF
863	!-
864
865	! we always thin down to the lower end of this range,
866	! which then gives our forest a few years to grow before
867	! we need to thin again
868	CALL thinning(taumin(ivm),taumax(ivm),thinstrat(ivm),&
869	diameters_temp(:),circ_temp(:),&
870	circ_class_kill(ipts,ivm,:,ifm_thin,icut_thin),&
871	circ_class_n(ipts,ivm,:),rdi_target_lower(ipts,ivm),&
872	ivm,ipts)
873
874	! Debug
875	IF(printlev_loc>=4 .AND. ipts == test_grid .AND. &
876	ivm == test_pft)THEN
877	WRITE(numout,*) 'Sapiens_forestry_main, circ_class_kill ',&
878	circ_class_kill(ipts,ivm,:,ifm_thin,icut_thin)
879	WRITE(numout,*) 'circ_class_n(ipts,ivm,:) ',&
880	circ_class_n(ipts,ivm,:)
881	DO icir=1,ncirc
882	WRITE(numout,*) 'Sapiens_forestry_main, circ_biomass(iparts)',&
883	icir,circ_class_biomass(ipts,ivm,icir,:,icarbon)
884	ENDDO
885	ENDIF
886	!-
887
888	ENDIF ! rdi .GT. rdi_target_upper
889
890	! Coppicing
891	ELSEIF(forest_managed(ipts,ivm).EQ.ifm_cop)THEN
892
893	! This is for coppicing.
894	! It makes no sense to coppice an evergreen tree species.
895	IF(is_needleleaf(ivm))THEN
896
897	WRITE(numout,*) 'WARNING', is_needleleaf(ivm)
898	WRITE(numout,*) 'Sapiens_forestry_main, Coppicing (FM = 3) should ' &
899	// 'only be done on deciduous trees.'
900	WRITE(numout,*) 'ipts,ivm,forest_managed(ipts,ivm): ',&
901	ipts,ivm,forest_managed(ipts,ivm)
902	IF (err_act.GT.1) CALL ipslerr (2,'sapiens_forestry_main', &
903	'Don t coppice conifers','','')
904
905	ENDIF
906
907	! Depending on the age of the stand we have different options.
908	! The current system is to do nothing until the diameter of the
909	! trees is greater than a certain value. Once this happens, we
910	! coppice, saving the number of trees that exist in this stand.
911	! The trees then regrow until they reach this diameter again.
912	! Some trees will have died due to mortality. At the second
913	! (and subsequent) cuts, we redistribute the below ground root
914	! mass to have the same number of trees as the first cut,
915	! killing the fine roots and harvesting aboveground biomass.
916
917	! Check first to see if the average diameter of the
918	! stand is greater than the prescribed treshold parameter.
919	diameters_temp(:)=wood_to_dia(&
920	circ_class_biomass(ipts,ivm,:,:,icarbon),ivm)
921	ave_tree_dia=zero
922	DO icir=1,ncirc
923	ave_tree_dia=ave_tree_dia+diameters_temp(icir)*&
924	circ_class_n(ipts,ivm,icir)
925	ENDDO
926	ave_tree_dia=ave_tree_dia/SUM(circ_class_n(ipts,ivm,:))
927
928	! Debug
929	IF(printlev_loc>=4)THEN
930	WRITE(numout,*) 'Sapiens_forestry_main, Do we coppice? '
931	WRITE(numout,*) 'ipts,ivm,SUM(circ_class_n(ipts,ivm,:)) ',&
932	ipts,ivm,SUM(circ_class_n(ipts,ivm,:))
933	WRITE(numout,*) 'ave_tree_dia,coppice_diameter, '&
934	//'coppice_dens(ipts,ivm) ',&
935	ave_tree_dia,coppice_diameter(ivm),coppice_dens(ipts,ivm)
936	ENDIF
937	!-
938
939	! We have already coppiced once, so we know our target
940	! tree density.
941	IF(ave_tree_dia .GE. coppice_diameter(ivm) .AND. &
942	coppice_dens(ipts,ivm) .GT. zero)THEN
943
944	! Debug
945	IF(printlev_loc>=4) WRITE(numout,*) 'Sapiens_forestry_main, '&
946	// 'Coppicing second cut. '
947	!-
948
949	IF(ok_change_species)THEN
950
951	! The criteria coppicing are met but check whether
952	! we want to switch to a conservation strategy on
953	! this site.
954	IF(fm_change_map(ipts,ivm) .EQ. ifm_none) THEN
955
956	! Change the management strategy in this case
957	! coppicing to conservation
958	forest_managed(ipts,ivm) = ifm_none
959
960	ELSEIF(fm_change_map(ipts,ivm).EQ.ifm_thin .OR. &
961	fm_change_map(ipts,ivm).EQ.ifm_src) THEN
962
963	! We would like to change the forest management so
964	! this is a good opportunity. However, we can not
965	! simply coppice because then the roots will stay
966	! on the site. We have two options depending on
967	! whether we also want to change the species or
968	! not.
969
970	! Calculate the pft corresponding with the youngest
971	! age class of this species.
972	igroup_new = species_change_map(ipts,ivm)
973	igroup_old = agec_group(ivm)
974
975	! Debug
976	IF(printlev_loc>=4)THEN
977	WRITE(numout,*) 'replant coppice 1+?'
978	WRITE(numout,*) 'igroup_new, igroup_old, ivm, ', &
979	igroup_new, igroup_old, ivm
980	ENDIF
981	!-
982
983	IF(igroup_new.EQ.igroup_old)THEN
984
985	! We have the intention to replant with the same
986	! species and change to high stand management so
987	! it is better not to harvest the trees and let
988	! the coppice develop in standards.
989	forest_managed(ipts,ivm) = ifm_thin
990
991	ELSE
992
993	! We want to replace this coppice by another
994	! species and to change the management strategy.
995	! We have to make sure that the whole stand is
996	! harvested. If so biomass will be zero and
997	! the sapiens_forestry_flag_species_change routine will flag this
998	! PFT to be replanted in species_change the
999	! management will be changed
1000	CALL clearcut_harvest(circ_class_n(ipts,ivm,:), &
1001	circ_class_kill(ipts,ivm,:,ifm_thin,icut_clear))
1002
1003	ENDIF !species_change
1004
1005	ELSE
1006
1007	! The species/management change scenario does
1008	! not prescribe a change to a conservation strategy
1009	! but the criteria for coppicing are met. We have
1010	! already coppiced once, so we know our target
1011	! tree density. We just have to mark all the trees
1012	! for coppicing.
1013	circ_class_kill(ipts,ivm,:,ifm_cop,icut_cop2) = &
1014	circ_class_n(ipts,ivm,:)
1015
1016	ENDIF
1017
1018	ELSE
1019
1020	! Species/management changes are not considered. The
1021	! criteria for a coppice are met. We have already
1022	! coppiced once, so we know our target tree density.
1023	! We just have to mark all the trees for coppicing.
1024	circ_class_kill(ipts,ivm,:,ifm_cop,icut_cop2) = &
1025	circ_class_n(ipts,ivm,:)
1026
1027	ENDIF
1028
1029
1030	! First time coppicing still need to save the coppice density
1031	ELSEIF(ave_tree_dia .GE. coppice_diameter(ivm))THEN
1032
1033	! debug
1034	IF(printlev_loc>=4)THEN
1035	WRITE(numout,*) 'Sapiens_forestry_main, Coppicing first cut. '
1036	WRITE(numout,*) 'coppice_dens(ipts,ivm): ',&
1037	SUM(circ_class_n(ipts,ivm,:))
1038	ENDIF
1039
1040	! save this tree density. We will reallocate to this in
1041	! future cuts.
1042	coppice_dens(ipts,ivm)=SUM(circ_class_n(ipts,ivm,:))
1043
1044	IF(ok_change_species)THEN
1045
1046	! The criteria for coppicing are met but check whether
1047	! we want to switch to a conservation strategy on this
1048	! site.
1049	IF( fm_change_map(ipts,ivm) .EQ. ifm_none) THEN
1050
1051	! Change the management strategy in this case
1052	! coppicing to conservation
1053	forest_managed(ipts,ivm) = ifm_none
1054
1055	ELSEIF(fm_change_map(ipts,ivm).EQ.ifm_thin .OR. &
1056	fm_change_map(ipts,ivm).EQ.ifm_src) THEN
1057
1058	! We would like to change the forest management so
1059	! this is a good opportunity. However, we can not
1060	! simply coppice because then the roots will stay
1061	! on the site. We have two options depending on
1062	! whether we also want to change the species or
1063	! not.
1064
1065	! Calculate the pft corresponding with the youngest
1066	! age class of this species.
1067	igroup_new = species_change_map(ipts,ivm)
1068	igroup_old = agec_group(ivm)
1069
1070	! Debug
1071	IF(printlev_loc>=4)THEN
1072	WRITE(numout,*) 'replant coppice 1?'
1073	WRITE(numout,*) 'igroup_new, igroup_old, ivm, ', &
1074	igroup_new, igroup_old, ivm
1075	ENDIF
1076	!-
1077
1078	IF(igroup_new.EQ.igroup_old)THEN
1079
1080	! We have the intention to replant with the same
1081	! species and change to high stand management so
1082	! it is better not to harvest the trees and let
1083	! the coppice develop in standards. Stop coppicing
1084	! and change the management
1085	forest_managed(ipts,ivm) = ifm_thin
1086
1087	ELSE
1088
1089	! We want to replace this coppice by another
1090	! species and by another management strategy.
1091	! We have to make sure that the whole stand is
1092	! harvested. If so, biomass will be zero and
1093	! the sapiens_forestry_flag_species_change routine will flag this
1094	! PFT to be replanted. In species_change the
1095	! management will be changed and the new PFT
1096	! will be set.
1097	CALL clearcut_harvest(circ_class_n(ipts,ivm,:), &
1098	circ_class_kill(ipts,ivm,:,ifm_thin,icut_clear))
1099
1100	ENDIF !species_change
1101
1102	ELSE
1103
1104	! The species/management change scenario does
1105	! not prescribe a change to a conservation strategy
1106	! but the criteria for a coppicing are met. This is
1107	! the first coppice. Mark all the trees
1108	! for coppicing.
1109	circ_class_kill(ipts,ivm,:,ifm_cop,icut_cop1) = &
1110	circ_class_n(ipts,ivm,:)
1111
1112	ENDIF
1113
1114	ELSE
1115
1116	! Species/management changes are not considered. The
1117	! criteria for a coppice are met. Mark all the trees
1118	! for coppicing.
1119	circ_class_kill(ipts,ivm,:,ifm_cop,icut_cop1) = &
1120	circ_class_n(ipts,ivm,:)
1121
1122	ENDIF
1123
1124	ELSE
1125
1126	! We don't need to coppice anything.
1127
1128	ENDIF
1129
1130	! This is for short rotation coppices.
1131	ELSEIF(forest_managed(ipts,ivm).EQ.ifm_src)THEN
1132
1133	! It makes no sense to coppice an evergreen tree. SRC is limited
1134	! to poplar and willow species.
1135	IF(.NOT. is_deciduous(ivm))THEN
1136
1137	WRITE(numout,*) 'ERROR: Sapiens_forestry_main, Short rotation ' &
1138	// ' coppicing (FM = 4)' &
1139	// 'should only be done on deciduous trees.'
1140	WRITE(numout,*) 'ipts,ivm,forest_managed(ipts,ivm): ',&
1141	ipts,ivm,forest_managed(ipts,ivm)
1142	IF (err_act.GT.1) CALL ipslerr (3,'ERROR: sapiens_forestry_main', &
1143	'SRC should only be done on deciduous trees.',&
1144	'Look in the output file for ERROR.','')
1145
1146	ENDIF
1147
1148	! Debug
1149	IF(printlev_loc>=4)THEN
1150	WRITE(numout,*) 'Sapiens_forestry_main, Do we SRC? '
1151	WRITE(numout,*) 'ipts,ivm,SUM(circ_class_n(ipts,ivm,:)) ',&
1152	ipts,ivm,SUM(circ_class_n(ipts,ivm,:))
1153	WRITE(numout,*) 'age_stand(ipts,ivm),src_rot_length(ivm),' &
1154	// ' coppice_dens(ipts,ivm): ',&
1155	age_stand(ipts,ivm),src_rot_length(ivm), &
1156	coppice_dens(ipts,ivm)
1157	WRITE(numout,*) 'src_nrots(ivm) ',&
1158	src_nrots(ivm)
1159	ENDIF
1160	!-
1161
1162	! NOTE: SRC is not implemented as a historic forest
1163	! management strategy but only used as a future strategy.
1164	! Hence, for the moment we never have the situation where we
1165	! decide we want to convert an SRC into an unmanaged forest.
1166	! This is not a very realist conversion anyway. Contrary to
1167	! ifm_none, ifm_cop and ifm_thin we don not check whether
1168	! the future management strategy is ifm_none.
1169
1170	! The qualifications for doing a short rotation coppice are
1171	! fairly simple. It's based on the age of the stand.
1172	! We harvest all the aboveground biomass. After a certain
1173	! number of rotations, we harvest the aboveground and kill
1174	! the belowground.
1175	! http://www.coppiceresources.co.uk/SRC.asp says that 30
1176	! years between total harvests.
1177	IF(MOD(age_stand(ipts,ivm),src_rot_length(ivm)) == 0)THEN
1178
1179	nrotations=age_stand(ipts,ivm)/src_rot_length(ivm)
1180
1181	IF(nrotations == src_nrots(ivm))THEN
1182
1183	! Debug
1184	IF(printlev_loc>=4) WRITE(numout,*) 'Sapiens_forestry_main, ' &
1185	// 'SRC final harvest. '
1186	!-
1187
1188	! This is the final harvest, where we will kill
1189	! all belowground biomass.
1190	circ_class_kill(ipts,ivm,:,ifm_src,icut_cop3)=&
1191	circ_class_n(ipts,ivm,:)
1192
1193	ELSEIF(nrotations == 1)THEN
1194
1195	! Debug
1196	IF(printlev_loc>=4) WRITE(numout,*) 'Sapiens_forestry_main, '&
1197	// 'SRC first harvest. '
1198	!-
1199
1200	! This is the first harvest, which will lead
1201	! to coppices. Save this tree density. We will
1202	! reallocate to this in future cuts.
1203	coppice_dens(ipts,ivm)=SUM(circ_class_n(ipts,ivm,:))
1204
1205	! Mark all trees for the first coppice cut.
1206	circ_class_kill(ipts,ivm,:,ifm_src,icut_cop1)= &
1207	circ_class_n(ipts,ivm,:)
1208
1209	ELSEIF(nrotations < src_nrots(ivm))THEN
1210
1211	! Debug
1212	IF(printlev_loc>=4) WRITE(numout,*) 'Sapiens_forestry_main, ' &
1213	// 'SRC standard harvest. '
1214	!-
1215
1216	! This is a standard harvest where we take the
1217	! aboveground biomass and leave the belowground.
1218	! Mark all trees for a standard coppice cut.
1219	circ_class_kill(ipts,ivm,:,ifm_src,icut_cop2)=&
1220	circ_class_n(ipts,ivm,:)
1221
1222	ELSE
1223
1224	WRITE(numout,*) 'ERROR: Sapiens_forestry_main, How did we arrive '&
1225	// 'here? We must have '
1226	WRITE(numout,*) 'somehow skipped the harvest in SRC. '
1227	WRITE(numout,*) 'ipts, ivm: ',ipts,ivm
1228	WRITE(numout,*) 'nrotations,src_nrots(ivm): ',&
1229	nrotations,src_nrots(ivm)
1230	IF(err_act.GT.1) CALL ipslerr(3,'ERROR: sapiens_forestry_main', &
1231	'Somehow skipped harvest in SRC.',&
1232	'Look in the output file for ERROR.','')
1233
1234	ENDIF ! check to see how many rotations we've done
1235
1236
1237	ENDIF ! checking if it's a harvest year
1238
1239	ELSE
1240
1241	WRITE(numout,*) 'Sapiens_forestry_main, We do not have any action '&
1242	//'for this management strategy'
1243	WRITE(numout,*) 'forest_managed(ipts,ivm),ipts,ivm ', &
1244	forest_managed(ipts,ivm),ipts,ivm
1245	IF(err_act.GT.1) CALL ipslerr(3,'ERROR: sapiens_forestry_main', &
1246	'No action for this management strategy.',&
1247	'Look in the output file for ERROR.','')
1248
1249	ENDIF ! test of management strategy
1250
1251	END IF
1252
1253	ENDDO pft
1254
1255	ENDDO
1256
1257	! Debug
1258	IF(printlev_loc>=4)THEN
1259	DO ipts=1,npts
1260	DO ivm=1,nvm
1261
1262	DO icir = 1,ncirc
1263	IF(ipts == test_grid .AND. ivm == test_pft)THEN
1264	WRITE(numout,*) 'Sapiens_forestry_main, How many trees do we kill?'
1265	WRITE(numout,*) 'Name icir test_pft circ_class_kill '
1266	WRITE(numout,*) 'Thin ',icir,test_pft,&
1267	circ_class_kill(test_grid,test_pft,:,ifm_thin,icut_thin)
1268	WRITE(numout,*) 'Clearcut ',icir,test_pft,&
1269	circ_class_kill(test_grid,test_pft,:,ifm_thin,icut_clear)
1270	WRITE(numout,*) 'Coppice 1 ',icir,test_pft,&
1271	circ_class_kill(test_grid,test_pft,:,ifm_cop,icut_cop1)
1272	WRITE(numout,*) 'Coppice 2 ',icir,test_pft,&
1273	circ_class_kill(test_grid,test_pft,:,ifm_cop,icut_cop2)
1274	WRITE(numout,*) 'SRC 1 ',icir,test_pft,&
1275	circ_class_kill(test_grid,test_pft,:,ifm_src,icut_cop1)
1276	WRITE(numout,*) 'SRC 2 ',icir,test_pft,&
1277	circ_class_kill(test_grid,test_pft,:,ifm_src,icut_cop2)
1278	WRITE(numout,*) 'SRC 3 ',icir,test_pft,&
1279	circ_class_kill(test_grid,test_pft,:,ifm_src,icut_cop3)
1280	ENDIF
1281	ENDDO
1282
1283	ENDDO
1284	ENDDO
1285	ENDIF
1286	!-
1287
1288
1289	!! 4. Check mass balance closure
1290
1291	IF (err_act.GT.1) THEN
1292
1293	! 3.2 Mass balance closure
1294	! 3.2.1 Calculate final biomass
1295	pool_end(:,:,:) = zero
1296	DO ipar = 1,nparts
1297	DO iele = 1,nelements
1298	DO icir = 1,ncirc
1299	pool_end(:,:,iele) = pool_end(:,:,iele) + &
1300	(circ_class_biomass(:,:,icir,ipar,iele) * &
1301	circ_class_n(:,:,icir) * veget_max(:,:))
1302	ENDDO
1303	ENDDO
1304	ENDDO
1305
1306	! 3.2.2 Calculate mass balance
1307	! Common processes
1308	DO iele=1,nelements
1309	check_intern(:,:,ipoolchange,iele) = -un * (pool_end(:,:,iele) - &
1310	pool_start(:,:,iele))
1311	ENDDO
1312
1313	closure_intern(:,:,:) = zero
1314	DO imbc = 1,nmbcomp
1315	DO iele=1,nelements
1316	! Debug
1317	closure_intern(:,:,iele) = closure_intern(:,:,iele) + &
1318	check_intern(:,:,imbc,iele)
1319	ENDDO
1320	ENDDO
1321
1322	! 3.3 Check mass balance closure
1323	CALL check_mass_balance("sapiens_forestry_main", closure_intern, npts, &
1324	pool_end, pool_start, veget_max, 'pft')
1325
1326	ENDIF ! err_act.GT.1
1327
1328
1329	IF (printlev >= 4) WRITE(numout,*) 'Leaving sapiens_forestry_main'
1330
1331	END SUBROUTINE sapiens_forestry_main
1332
1333
1334
1335	!! ================================================================================================================================
1336	!! SUBROUTINE : clearcut_harvest
1337	!!
1338	!>\BRIEF Marks all the trees at this point and this PFT for clearing
1339	!!
1340	!! DESCRIPTION : This marks all the trees in this stand for killing. The
1341	!! difference between this and clearcut_noharvest is that
1342	!! these trees are moved into the harvest pools in stomate_kill,
1343	!! where only the branches should be left onsite.
1344	!!
1345	!! RECENT CHANGE(S): None
1346	!!
1347	!! MAIN OUTPUT VARIABLE(S): ::circ_class_kill; the number of trees in each circumference class to kill
1348	!!
1349	!! REFERENCE(S) : See above, module description.
1350	!!
1351	!! FLOWCHART :
1352	!! \latexonly
1353	!! \includegraphics[scale=0.5]{clearcutflow.jpg}
1354	!! \endlatexonly
1355	!! \n
1356	!_ ================================================================================================================================
1357
1358	SUBROUTINE clearcut_harvest(circ_class_n_temp, circ_class_kill_temp)
1359
1360	IMPLICIT NONE
1361
1362	!! 0 Variable and parameter declaration
1363
1364	!! 0.1 Input variables
1365	REAL(r_std), DIMENSION(:), INTENT(in) :: circ_class_n_temp !! Number of trees in each circumference
1366	!! class
1367	!! 0.2 Output variables
1368	REAL(r_std), DIMENSION(:), INTENT(out) :: circ_class_kill_temp !! Number of trees within a circ that needs
1369	!! to be killed @tex $(ind m^{-2})$ @endtex
1370
1371	!! 0.3 Modified variables
1372
1373	!! 0.4 Local variables
1374
1375
1376	!_ ================================================================================================================================
1377
1378
1379	! This routine has gotten much simpler. All we want to do here is to
1380	! indicate how many trees we want to kill. The actual killing of
1381	! the trees will be done in lpj_gap. As this is a clearcut, we
1382	! want to kill all of the trees.
1383	circ_class_kill_temp(:)=circ_class_n_temp(:)
1384
1385
1386	END SUBROUTINE clearcut_harvest
1387
1388	!! ================================================================================================================================
1389	!! SUBROUTINE : thinning
1390	!!
1391	!>\BRIEF Thinning of the stand in gridpoint i for PFT j according to the
1392	!! management type (self-thinning, anthropogenic thinning, smoothed anthropogenic
1393	!! thinning or coppicing). Selects which trees are thinned in order to respect
1394	!! self-thinning rules or rdi_objective. Only trees with circumference lower than
1395	!! circ_lim can be thinned.
1396	!!
1397	!! DESCRIPTION : None
1398	!!
1399	!! RECENT CHANGE(S): None
1400	!!
1401	!! MAIN OUTPUT VARIABLE(S): ::circ_ij0; circumference of all individual trees
1402	!! before thinning for 1 ha at gridpoint i and for PFT j (m), ::vol_thinned
1403	!! Volume of thinned trees including waste wood \f$(m^3 ha^{-1})\f$ all other readjusted
1404	!! stand-level variables: ::rdi, ::ind, ::av_height, ...
1405	!!
1406	!! REFERENCE(S) : See above, module description.
1407	!!
1408	!! FLOWCHART :
1409	!! \latexonly
1410	!! \includegraphics[scale=0.5]{thinningflow.jpg}
1411	!! \endlatexonly
1412	!! \n
1413	!_ ================================================================================================================================
1414
1415	SUBROUTINE thinning(tmin,tmax,thstrat,diameters_temp, &
1416	circ_temp,circ_class_kill_temp,circ_class_n_temp,&
1417	rdi_limit,ivm,ipts)
1418
1419	!! 0 Variable and parameter declaration
1420
1421	!! 0.1 Input variables
1422	REAL(r_std), INTENT(in) :: tmin
1423	REAL(r_std), INTENT(in) :: tmax
1424	REAL(r_std), INTENT(in) :: thstrat
1425	REAL(r_std), INTENT(in) :: rdi_limit
1426	REAL(r_std), DIMENSION(:), INTENT(in) :: diameters_temp
1427	REAL(r_std), DIMENSION(:), INTENT(in) :: circ_temp
1428	REAL(r_std), DIMENSION(:), INTENT(in) :: circ_class_n_temp
1429	INTEGER(i_std),INTENT(IN) :: ivm
1430	INTEGER(i_std),INTENT(IN) :: ipts
1431
1432
1433	!! 0.2 Output variables
1434	REAL(r_std), DIMENSION(:), INTENT(out) :: circ_class_kill_temp
1435
1436	!! 0.3 Modified variables
1437
1438	!! 0.4 Local variables
1439	REAL(r_std), DIMENSION(ncirc) :: probability
1440	REAL(r_std), DIMENSION(ncirc) :: circ_class_kill_pot
1441	REAL(r_std) :: circmax
1442	REAL(r_std) :: circmin
1443	REAL(r_std) :: current_rdi
1444	REAL(r_std) :: potential_rdi
1445	REAL(r_std) :: dead_trees
1446	REAL(r_std) :: previous_step
1447	INTEGER :: icir,endcir,startcir,inccir,nsteps_tried
1448
1449
1450	INTEGER :: istep1
1451
1452	!_ ================================================================================================================================
1453
1454	! We are going to assign a thinning probability to every
1455	! circumference class, based on Eq. 12 in Bellassen 2010.
1456	! He does it there for individual trees. The equation
1457	! differs based on the thinning strategy used (from above or below).
1458	! Note that although the biomasses are sorted this does not guarantee
1459	! that the circumferences are sorted as well because circumference
1460	! largely depends on the heartwood mass but biomass also contains
1461	! more volatile pools such as leaves and roots and these pools
1462	! are not always in equilibrium (they may have too little leaves/roots for
1463	! their sapwood). If the model experiences some limitations (for example
1464	! C-limitation when doing allocation) there are some rules that first
1465	! allocate to a certain circ_class. If this happened the total biomass
1466	! of a tree in, e.g. the second circ_class may be higher than total biomass
1467	! of a tree in the third circ_class but the circumference of the tree in the
1468	! third class will still exceed that of the second class. Long story to
1469	! justify why MINVAL and MAXVAL are used instead of circ_temp(1) and
1470	! circ_temp(ncirc)
1471	circ_class_kill_temp(:)=zero
1472	circmin=MINVAL(circ_temp(:))
1473	circmax=MAXVAL(circ_temp(:))
1474	probability(:)=zero
1475
1476	DO icir=1,ncirc
1477
1478	IF(thstrat .GE. zero)THEN
1479
1480	! thinning from below
1481	IF(ncirc == 1)THEN
1482
1483	! In this case, circmin=circmax and the denominator will be
1484	! zero. However, if we have only one circ class it doesn't really
1485	! matter what the probability is set to, because all the biomass
1486	! has to come from this circ class. So we set it equal to the
1487	! minimum probability.
1488	probability(icir)=tmin
1489
1490	ELSE
1491
1492	probability(icir)=tmin+(tmax-tmin)*((circmax-circ_temp(icir))/&
1493	(circmax-circmin))**thstrat
1494
1495	ENDIF ! ncirc ==1
1496
1497	ELSE !thstrat
1498
1499	! thinning from above
1500	IF(ncirc == 1)THEN
1501
1502	! Same issue as above
1503	probability(icir)=tmin
1504
1505	ELSE
1506
1507	probability(icir)=tmin+(tmax-tmin)*&
1508	((circ_temp(icir)-circmin)/(circmax-circmin))**ABS(thstrat)
1509
1510	ENDIF !ncirc == 1
1511
1512	ENDIF ! thstrat
1513
1514	! Debug
1515	IF(printlev_loc>=4 .AND. ivm == test_pft .AND. ipts == test_grid)THEN
1516	WRITE(numout,*) 'probability',probability(icir),icir
1517	WRITE(numout,*) 'circ_temp',circ_temp(icir),icir
1518	WRITE(numout,*) 'tmin,tmax',tmin,tmax
1519	WRITE(numout,*) 'circmin,circmax',circmin,circmax
1520	WRITE(numout,*) 'thstrat',thstrat
1521	ENDIF
1522	CALL flush(numout)
1523	!-
1524
1525	ENDDO
1526
1527	! now we need to kill enough trees to get down to the proper density. If we
1528	! are thinning from above, we start from the highest class, and from
1529	! below, we start from the lowest. One thing to keep in mind here is
1530	! that it is possible that we have not removed enough trees after looping
1531	! through all the circ classes one time, so we need to keep looping through
1532	! until our density is low enough.
1533	IF(thstrat .GE. zero)THEN
1534	startcir=1
1535	endcir=ncirc
1536	inccir=1
1537	ELSE
1538	startcir=ncirc
1539	endcir=1
1540	inccir=-1
1541	ENDIF
1542
1543	istep1=0
1544	thin_out:DO
1545	istep1=istep1+1
1546
1547	IF(istep1 .GT. 1000)THEN
1548
1549	! Something is wrong. Why did it take us so many steps?
1550	WRITE(numout,*) 'WARNING: Taking too many steps in thinning'
1551	WRITE(numout,*) 'ipts,ivm: ',ipts,ivm
1552	WRITE(numout,*) 'current_rdi, rdi_limit: ',current_rdi,rdi_limit
1553	IF(istep1 .GT. 1001)THEN
1554	WRITE(numout,*) 'Exiting the loop'
1555	CALL flush(numout)
1556	EXIT thin_out
1557	ENDIF
1558
1559	ENDIF
1560
1561	DO icir=startcir,endcir,inccir
1562
1563	! check the density to see if we have gotten enough trees.
1564	! Notice we only want to include the trees that will still
1565	! be living at the end of this step
1566	current_rdi=calculate_rdi(diameters_temp(:),circ_class_n_temp(:)-&
1567	circ_class_kill_temp(:),ivm)
1568
1569	! Debug
1570	IF(printlev_loc>=4 .AND. ivm == test_pft .AND. &
1571	ipts == test_grid)THEN
1572	WRITE(numout,*) 'Debugging thinning',ivm
1573	WRITE(numout,*) 'icir,startcir,endcir,inccir',&
1574	icir,startcir,endcir,inccir
1575	WRITE(numout,*) 'rdi_limit',rdi_limit
1576	WRITE(numout,*) 'current_rdi', current_rdi
1577	CALL flush(numout)
1578	ENDIF
1579	!-
1580
1581	IF( current_rdi .LE. rdi_limit)THEN
1582	EXIT thin_out
1583	ENDIF
1584
1585	! if not, we need to kill some or all of the trees in this class
1586	! this is the maximum number of trees we can kill in this circ class
1587	dead_trees=probability(icir)*&
1588	(circ_class_n_temp(icir)-circ_class_kill_temp(icir))
1589
1590	! if we kill all these trees, can we get the RDI as low as we want?
1591	circ_class_kill_pot(:)=circ_class_kill_temp(:)
1592	circ_class_kill_pot(icir)=circ_class_kill_temp(icir)+dead_trees
1593	potential_rdi=calculate_rdi(diameters_temp(:),&
1594	circ_class_n_temp(:)-circ_class_kill_pot(:),ivm)
1595
1596	! Debug
1597	IF(printlev_loc>=4 .AND. ivm == test_pft .AND. ipts == test_grid)THEN
1598	WRITE(numout,*) 'potential_rdi,rdi_limit',potential_rdi,rdi_limit
1599	WRITE(numout,*) 'diameters_tmep',diameters_temp(icir)
1600	WRITE(numout,*) 'circ_class_n_temp',circ_class_n_temp(icir)
1601	WRITE(numout,*) 'circ_class_kill_pot',circ_class_kill_pot(icir)
1602	WRITE(numout,*) 'circ_class_kill_temp',circ_class_kill_temp(icir)
1603	WRITE(numout,*) 'dead_trees',dead_trees
1604	WRITE(numout,*) 'probability',probability(icir)
1605	CALL flush(numout)
1606	ENDIF
1607	!-
1608
1609	IF(circ_class_kill_pot(icir) .LT. -min_stomate)THEN
1610	WRITE(numout,*) 'ERROR: Negative value of circ_class_kill_pot.'
1611	WRITE(numout,*) 'ipts, ivm, icir: ',ipts,ivm,icir
1612	WRITE(numout,*) 'pot ',circ_class_kill_pot(icir)
1613	WRITE(numout,*) 'kill ',circ_class_kill_temp(icir)
1614	WRITE(numout,*) 'n ',circ_class_n_temp(icir)
1615	WRITE(numout,*) 'dead_trees ',dead_trees
1616	WRITE(numout,*) 'probability',probability(icir)
1617	IF (err_act.GT.1) CALL ipslerr_p (3,'forestry', &
1618	'Bad value of circ_class_kill_pot.',&
1619	'Look in the output file for ERROR.',&
1620	'')
1621	ENDIF
1622
1623	IF(potential_rdi .LE. rdi_limit)THEN
1624
1625	! yes we can. How do we know how many trees to kill, though?
1626	! The problem is that the RDI depends on the diameters of the
1627	! remaining trees in a non-linear fashion. Let us do a simple
1628	! bisection search. If we are here, we know that
1629	! circ_class_kill_pot(icir) is too high. So cut it in half and
1630	! try again.
1631	previous_step=dead_trees/deux
1632	circ_class_kill_pot(icir)=circ_class_kill_temp(icir)+previous_step
1633	nsteps_tried=0
1634	inner_loop:DO
1635
1636	potential_rdi=calculate_rdi(diameters_temp(:),&
1637	circ_class_n_temp(:)-circ_class_kill_pot(:),ivm)
1638	nsteps_tried=nsteps_tried+1
1639
1640	!+++CHECK+++
1641	! this value should be externalised
1642	! Our RDI is within our tolerance of the target, so we're good.
1643	IF( ABS(potential_rdi - rdi_limit) .LE. 0.001_r_std)THEN
1644	EXIT inner_loop
1645	ENDIF
1646	!+++++++++++
1647
1648	! Is our potential RDI too high or too low?
1649	previous_step=previous_step/deux
1650	IF(potential_rdi .GT. rdi_limit)THEN
1651
1652	! We need to kill more trees, so we add half of our
1653	! previous step size back.
1654	circ_class_kill_pot(icir)=circ_class_kill_pot(icir)+&
1655	previous_step
1656
1657	ELSE
1658
1659	! We need to kill fewer trees, so we take half of our
1660	! previous step size away.
1661	circ_class_kill_pot(icir)=circ_class_kill_pot(icir)-&
1662	previous_step
1663
1664	ENDIF
1665
1666	! Debug
1667	IF(ivm == test_pft .AND. ipts == test_grid .AND. &
1668	printlev_loc>=4)THEN
1669	WRITE(numout,*) 'nsteps_tried ',nsteps_tried
1670	WRITE(numout,*) 'pot ',circ_class_kill_pot(:)
1671	WRITE(numout,*) 'kill ',circ_class_kill_temp(:)
1672	WRITE(numout,*) 'n ',circ_class_n_temp(:)
1673	CALL flush(numout)
1674	ENDIF
1675	!-
1676
1677	!+++CHECK+++
1678	! this value should be externalised
1679	IF(nsteps_tried .GT. 100)THEN
1680	WRITE(numout,*) 'ERROR:Trying too many steps in '&
1681	//'the bisection search in forestry thinning.'
1682	WRITE(numout,*) 'ivm: ',ivm
1683	WRITE(numout,*) 'pot ',circ_class_kill_pot(:)
1684	WRITE(numout,*) 'kill ',circ_class_kill_temp(:)
1685	WRITE(numout,*) 'n ',circ_class_n_temp(:)
1686	IF(err_act.GT.1) CALL ipslerr_p (3,'ERROR: forestry', &
1687	'Too many steps in the bisection.',&
1688	'Look in the output file for ERROR.','')
1689	ENDIF
1690	!+++++++++++
1691
1692	ENDDO inner_loop
1693
1694	! Schedule the trees for killing.
1695	circ_class_kill_temp(icir)=circ_class_kill_pot(icir)
1696	EXIT thin_out
1697
1698	ELSE
1699
1700	! No we can't, so we have to kill all possible trees and
1701	! move onto the next class
1702	circ_class_kill_temp(icir)=circ_class_kill_temp(icir)+dead_trees
1703
1704	ENDIF
1705
1706
1707	ENDDO
1708
1709	! check to see if we have any trees left. We should!
1710	! If not, there is a problem, and we need to catch it
1711	! to avoid an infinite loop.
1712	IF( (SUM(circ_class_n_temp(:)) - SUM(circ_class_kill_temp(:))) .LE. &
1713	min_stomate)THEN
1714	WRITE(numout,*) 'ERROR: We are thinning, but we have no trees left!'
1715	WRITE(numout,*) 'ivm ',ivm
1716	WRITE(numout,*) 'kill ',circ_class_kill_temp(:)
1717	WRITE(numout,*) 'n ',circ_class_n_temp(:)
1718	IF(err_act.GT.1) CALL ipslerr_p (3,'ERROR: forestry', &
1719	'Thinning, but no trees left.',&
1720	'Look in the output file for ERROR.','')
1721	ENDIF
1722
1723	ENDDO thin_out
1724
1725	END SUBROUTINE thinning
1726
1727
1728	!! ================================================================================================================================
1729	!! SUBROUTINE sapiens_forestry_species_change
1730	!!
1731	!>\BRIEF When a PFT is harvest a new PFT with another species and management strategy
1732	!! is replanted.
1733	!!
1734	!! DESCRIPTION : When a PFT is harvest a new PFT with another species and management strategy
1735	!! is replanted. The routine creates new values for veget_max_new which are then
1736	!! used in sapiens_lcchange to really change the PFTs.
1737	!!
1738	!! The routine sapiens_forestry_species_change is only used when ::ok_change_species is .TRUE. and
1739	!! makes use of the information contained in the species_change map and the
1740	!! desired_fm_map. As an alternative to these maps single values can be specified
1741	!! in the run.def for species_change_force and fm_change_force. Although in
1742	!! species change and management change are more or less independently called in
1743	!! in the code the objective is to call them at the same time. This independence
1744	!! was used a bit for debugging but hasn't been tested.
1745	!!
1746	!!
1747	!! RECENT CHANGE(S) : None
1748	!!
1749	!! MAIN OUTPUT VARIABLE(S): ::forest_managed, ::veget_max_new
1750	!!
1751	!! REFERENCE(S) :
1752	!!
1753	!! FLOWCHART :
1754	!! \n
1755	!_ ================================================================================================================================
1756
1757
1758	SUBROUTINE sapiens_forestry_species_change(npts, lpft_replant, veget_max_new_species, forest_managed, &
1759	species_change_map, veget_max, fm_change_map)
1760
1761	IMPLICIT NONE
1762
1763	!! 0. Variable and parameter declaration
1764
1765	!! 0.1 Input variables
1766	INTEGER(i_std), INTENT(in) :: npts !! Domain size (unitless)
1767	INTEGER(i_std), DIMENSION(:,:), INTENT(in) :: species_change_map !! A map which gives the PFT number that each
1768	!! PFT will be replanted as in case of a clearcut.
1769	!! (1-nvm,unitless)
1770	INTEGER(i_std), DIMENSION(:,:), INTENT(in) :: fm_change_map !! A map which gives the desired FM strategy when
1771	!! the PFT will be replanted after a clearcut.
1772	!! (1-nvm,unitless)
1773
1774	!! 0.2 Output variables
1775
1776
1777	!! 0.3 Modified variables
1778	REAL(r_std), DIMENSION(:,:), INTENT(inout) :: veget_max !! "Maximal" coverage fraction of a PFT (LAI
1779	!! -> infinity) on ground. Includes
1780	!! the nobio fraction, so may sum to
1781	!! less than unity.
1782	REAL(r_std), DIMENSION(:,:),INTENT(out) :: veget_max_new_species !! New "maximal" coverage fraction of a PFT
1783	!! (LAI -> infinity) Includes
1784	!! the nobio fraction, so may sum to
1785	!! less than unity. (unitless)
1786	LOGICAL, DIMENSION(:,:), INTENT(inout) :: lpft_replant !! Set to true if a PFT has been clearcut
1787	!! and needs to be replaced by another species
1788	INTEGER(i_std), DIMENSION (:,:), INTENT(inout) :: forest_managed !! forest management flag (is the forest
1789	!! being managed?)
1790
1791	!! 0.4 Local variables
1792	INTEGER :: ipts,ivm,count !! Indices
1793	INTEGER :: ivm_young_new_species !! Indices
1794	INTEGER :: ivm_young_old_species !! Indices
1795	REAL(r_std) :: diff !! Difference in surface area before and after
1796	!! a species change (should be zero!)
1797	REAL(r_std) :: temp_veget_max
1798	!================================================================================================================================
1799
1800	! Initialize variable(s)
1801	veget_max_new_species(:,:)=zero
1802	count=zero
1803
1804	! Check whether replanting is required and if so update
1805	! the management strategy where needed.
1806	DO ipts = 1,npts
1807
1808	! Search the PFT that needs to be replanted and replant it
1809	DO ivm=1,nvm
1810
1811	IF (is_tree(ivm)) THEN
1812
1813	! Was the PFT harvested this year? If so, there is an
1814	! opportunity to replant
1815	IF(lpft_replant(ipts,ivm))THEN
1816
1817	! Keep track of the number of species changes. In
1818	! case none occur veget_max_new should equal veget_max
1819	! else the mass balance check in LCC will think
1820	! there was a change in PFTs and will generate errors.
1821	count = count + 1
1822
1823	! Plant the new PFT in the youngest age class. Find
1824	! out which PFT is the youngest member of the group
1825	! that we want to regrow.
1826	ivm_young_new_species = zero
1827	DO
1828
1829	! Repeat this loop until the matching age class is found
1830	! or an error occurs
1831	ivm_young_new_species=ivm_young_new_species+1
1832
1833	IF(ivm_young_new_species .GT. nvm)THEN
1834	WRITE(numout,*) "ERROR: Can not find young PFT for this age class."
1835	WRITE(numout,*) "ipts, ivm: ",ipts, ivm
1836	WRITE(numout,*) "species_change_map(ipts,:): ",species_change_map(ipts,:)
1837	CALL ipslerr_p (3,'sapiens_forestry_species_change', &
1838	'Could not find the youngest PFT in this age class.',&
1839	'Something may be wrong with the species_change_map.',&
1840	'Check the out_orchidee file for details.')
1841	ENDIF
1842
1843	! It is confirmed that we found ivm_young_new_species for the pft
1844	! under consideration. This test requires that the age classes are
1845	! defined such that the youngest age class comes first and that the
1846	! age classes of the same species are grouped (see sapiens_lcchange.f90
1847	! where the age classes are defined)
1848	IF(agec_group(ivm_young_new_species) == species_change_map(ipts,ivm))THEN
1849	EXIT
1850	ENDIF
1851
1852	ENDDO ! Searching youngest age class for this species
1853
1854	! The lpft_replant flag is set in sapiens_forestry or stomate_kill.
1855	! Soon after those routines the killing takes place and when
1856	! age classes are used veget_max is updated. This means that
1857	! the veget_max of the oldest age class that was just harvested
1858	! has been moved to the youngest age class of the species that
1859	! was harvested. Possibly there was already something in this youngest
1860	! age class from previous time steps. The veget_max of the youngest
1861	! ages class may thus be larger than the veget_max of the age class
1862	! that was harvested at this time step.
1863	ivm_young_old_species = start_index(agec_group(ivm))
1864
1865	! Is there really a species change?
1866	IF (ivm .EQ. ivm_young_new_species) THEN
1867
1868	! lpft_replant is true but there is no real species change
1869	! Changes from PFTs earlier in the list could already have
1870	! resulted in an increase of veget_max_new_species
1871	!veget_max_new_species(ipts,ivm) = veget_max_new_species(ipts,ivm) + &
1872	! veget_max(ipts,ivm)
1873	veget_max_new_species(ipts,ivm_young_new_species) = &
1874	veget_max_new_species(ipts,ivm_young_new_species) + &
1875	veget_max(ipts,ivm)
1876
1877	!---TEMP---
1878	IF(printlev_loc>=4) THEN
1879	WRITE(numout,*) 'Species change - PFT is the same', ivm
1880	WRITE(numout,*) 'Management type, ', &
1881	forest_managed(ipts,ivm)
1882	ENDIF
1883	!----------
1884
1885	ELSE
1886
1887	! There is a real species change. Add the ground area of the
1888	! old PFT to the ground area of the new PFT.
1889	! +++CHECK+++
1890	! It seems that in relatively few cases the veget_max of the
1891	! killed PFT has not been moved to the youngest age class
1892	! of that PFT - Need to find when this happens and fix the
1893	! problem where it occurs (may be it is related to FM3?)
1894	! but here we added a quick patch. Note that
1895	! this code is also used in the case where the youngest age
1896	! class of the PFT was removed. The code moves the whole PFT
1897	! at once the old PFT can be set to zero.
1898	IF (veget_max(ipts,ivm) .GT. min_stomate) THEN
1899
1900	! Replant with the new species
1901	veget_max_new_species(ipts,ivm_young_new_species) = &
1902	veget_max_new_species(ipts,ivm_young_new_species) + &
1903	veget_max(ipts,ivm)
1904	veget_max_new_species(ipts,ivm) = zero
1905
1906	! Patch the old species - in land cover change
1907	! (adjust_delta_veget) we check whether the changes make
1908	! sense. Those checks will fail unless we make it look
1909	! as if the the youngest PFT was replanted by another
1910	! species
1911	IF (ivm .EQ. ivm_young_old_species) THEN
1912	WRITE(numout,*) &
1913	'1 st class of age cant be replaced by itself' , ivm, ivm_young_old_species
1914	ELSE
1915	!veget_max(ipts,ivm_young_old_species) = &
1916	! veget_max(ipts,ivm_young_old_species) + veget_max(ipts,ivm)
1917	!veget_max(ipts,ivm) = zero
1918	ENDIF
1919	!---TEMP---
1920	IF(printlev_loc>=4) THEN
1921	WRITE(numout,*) 'Species change - patch', ivm
1922	WRITE(numout,*) 'Future PFT',veget_max(ipts,ivm_young_old_species)
1923	WRITE(numout,*) 'current PFT',veget_max(ipts,ivm)
1924	WRITE(numout,*) 'Management type, ', forest_managed(ipts,ivm)
1925	ENDIF
1926	!----------
1927
1928	! Pure case
1929	ELSEIF (veget_max(ipts,ivm_young_old_species) .GT. min_stomate) THEN
1930
1931	! The PFT with the replant label is for example 47. Because the
1932	! PFT has been correctly killed it was replanted in PFT 45 but
1933	! for PFT 45 the replant flag is FALSE so the surface area will
1934	! be copied from veget_max(ipts,45) to
1935	! veget_max_new_species(ipts,45). We will first undo the copy
1936	! from 47 to 45. Say that,for example, PFT 45 (young_old) had to
1937	! become PFT 13 (young_new). We will then replant the correct
1938	! surface are with PFT 13.
1939	veget_max_new_species(ipts,ivm_young_old_species) = &
1940	veget_max_new_species(ipts,ivm_young_old_species) - &
1941	veget_max(ipts,ivm_young_old_species)
1942	veget_max_new_species(ipts,ivm_young_new_species) = &
1943	veget_max_new_species(ipts,ivm_young_new_species) + &
1944	veget_max(ipts,ivm_young_old_species)
1945
1946	!---TEMP---
1947	IF(printlev_loc>=4) THEN
1948	WRITE(numout,*) 'Species change - pure case', ivm
1949	WRITE(numout,*) 'management type, ', &
1950	forest_managed(ipts,ivm_young_new_species)
1951	ENDIF
1952	!----------
1953
1954	! Unexpected case
1955	ELSE
1956
1957	WRITE(numout,*) 'WARNING: species change - unexpected case'
1958	IF(err_act.GT.1)THEN
1959	WRITE(numout,*) "ipts, ivm: ",ipts, ivm
1960	CALL ipslerr_p (3,'sapiens_forestry', &
1961	'Species change','unexpected case in IF-loop','')
1962	ENDIF
1963
1964	ENDIF
1965	! +++++++++++
1966
1967	ENDIF ! Anthropogenic species change?
1968
1969	! Use the new management strategy. The way the desired management
1970	! maps were made follows the species logic. A conifer species will
1971	! never have FM 3 but a deciduous species can. If a conifer is
1972	! changed to a deciduous we need to use the FM of the new deciduous
1973	! species in this example FM3. This is reflected in the code by
1974	! using the index fom ivm_young_new_species in the left and the right
1975	! hand side. Note that we don't care whether the new PFT has already
1976	! some biomass in it or not. Irrespective of its condition we change
1977	! FM to the desired strategy.
1978	forest_managed(ipts,ivm_young_new_species) = &
1979	fm_change_map(ipts,ivm_young_new_species)
1980
1981	!---TEMP---
1982	IF(printlev_loc>=4)THEN
1983	WRITE(numout,*) 'Anthropogenic species change ?'
1984	WRITE(numout,*) 'ipts, ivm, ',ipts, ivm
1985	WRITE(numout,*) 'ivm_young_old_species',ivm_young_old_species
1986	WRITE(numout,*) 'Species we want to change to, ',species_change_map(ipts,ivm)
1987	WRITE(numout,*) 'youngest age class for that species, ', &
1988	ivm_young_new_species
1989	WRITE(numout,*) 'forest_managed, ',forest_managed(ipts,ivm_young_new_species)
1990	ENDIF
1991	!----------
1992
1993	ELSE
1994
1995	! No replanting !
1996	! veget_max_new will be the same as the old_veget_max.
1997	! We add it here just in case another PFT was added to this one
1998	! above.
1999	veget_max_new_species(ipts,ivm) = &
2000	veget_max_new_species(ipts,ivm)+veget_max(ipts,ivm)
2001
2002	! Given that forest_managed is an inout variable nothing should
2003	! be done because we just keep the initial management strategy
2004
2005	ENDIF
2006
2007	ELSE
2008
2009	! Bare soil, grassland and cropland
2010	! veget_max_new will be the same as the old_veget_max.
2011	! We add it here just in case another PFT was added to this one
2012	! above.
2013	veget_max_new_species(ipts,ivm) = &
2014	veget_max_new_species(ipts,ivm) + veget_max(ipts,ivm)
2015
2016	ENDIF
2017
2018	ENDDO !nvm
2019
2020	! Debug
2021	IF(printlev_loc>=4)THEN
2022	WRITE(numout,*) 'End of routine - final values'
2023	WRITE(numout,*) 'pixel,PFT, veget_max, veget_max_new_species'
2024	DO ivm=1,nvm
2025	WRITE(numout,*) ipts, ivm,veget_max(ipts,ivm),&
2026	veget_max_new_species(ipts,ivm)
2027	ENDDO
2028
2029	! Write warning
2030	diff = SUM(veget_max_new_species(ipts,:)-veget_max(ipts,:))
2031	IF( ABS(diff) .GT. min_stomate )THEN
2032	WRITE(numout,*) 'WARNING: species - change surface area is not preserved'
2033	CALL ipslerr_p (2,'sapiens_forestry', &
2034	& 'sapiens_forestry_species_change','surface area is not preserved','')
2035	ELSE
2036	WRITE(numout,*) 'Losses and gains cancel each other out'
2037	ENDIF
2038	diff = veget_max_new_species(ipts,1)-veget_max(ipts,1)
2039	IF( ABS(diff) .GT. min_stomate )THEN
2040	WRITE(numout,*) 'WARNING: species - area of PFT1 is not preserved'
2041	CALL ipslerr_p (2,'sapiens_forestry', &
2042	& 'sapiens_forestry_species_change','PFT1 is not preserved','')
2043	ELSE
2044	WRITE(numout,*) 'No leakage to PFT1'
2045	ENDIF
2046	ENDIF
2047	!-
2048
2049	ENDDO ! npts
2050
2051	! Reset this variable
2052	lpft_replant(:,:)=.FALSE.
2053
2054	! In case no species changes occured veget_max_new
2055	! should equal veget_max else the mass balance
2056	! check in LCC will think there was a change in
2057	! PFTs and will generate errors.
2058	IF (count .LT. min_stomate) THEN
2059	veget_max_new_species(:,:) = veget_max(:,:)
2060	ENDIF
2061
2062	END SUBROUTINE sapiens_forestry_species_change
2063
2064
2065	!! ================================================================================================================================
2066	!! SUBROUTINE : sapiens_forestry_flag_species_change
2067	!!
2068	!>\BRIEF: Check whether there is an opportunity to replant a different species
2069	!!
2070	!! DESCRIPTION : PFTs that died or were harvested should all be empty. This
2071	!! is the information that is used to decide whether there is
2072	!! an opportunity to replant.
2073	!!
2074	!! RECENT CHANGE(S): None
2075	!!
2076	!! MAIN OUTPUT VARIABLE(S): None
2077	!!
2078	!! REFERENCE(S) : None
2079	!!
2080	!! FLOWCHART : None
2081	!! \n
2082	!_ ================================================================================================================================
2083	SUBROUTINE sapiens_forestry_flag_species_change(npts, veget_max, circ_class_biomass, circ_class_n, &
2084	lpft_replant, forest_managed)
2085
2086	!! 0. Variable and parameter description
2087
2088	!! 0.1 Input variables
2089	INTEGER(i_std), INTENT(in) :: npts !! Number of pixels
2090	REAL(r_std), DIMENSION(:,:), INTENT(in) :: veget_max !! Veget_max at the moment this routine is checked
2091	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(in) :: circ_class_biomass !! ciconference level biomass @tex $(gC m^{-2})$ @endtex
2092	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: circ_class_n !! Density of individuals
2093	INTEGER(i_std), DIMENSION (:,:), INTENT(in) :: forest_managed !! forest management flag (is the forest
2094	!! being managed?)
2095
2096	!! 0.2 Output variables
2097
2098	!! 0.3 Modified variables
2099	LOGICAL, DIMENSION(:,:), INTENT(inout) :: lpft_replant !! Set to true if a PFT has been clearcut
2100	!! and needs to be replaced by another species
2101
2102	!! 0.4 Local variables
2103	INTEGER(i_std) :: ipts,ivm !! Indices
2104	INTEGER(i_std) :: error !! Counter for the numbers of errors
2105	REAL(r_std) :: biomass
2106	REAL(r_std) :: ind
2107
2108	!===============================================================================================================================
2109
2110	IF(printlev_loc.GT.4) WRITE(numout,*) 'Entering: change sapiens_forestry_flag_species_change'
2111
2112	!! 1. Initilaze variables
2113	error = zero
2114	!! 2. Check whether there is an opportunity to change the
2115	! species. If so flag the opportunity so we can deal
2116	! with the details later.
2117	DO ivm = 2,nvm
2118
2119	! Only forest PFT can be replanted for the moment
2120	IF(is_tree(ivm))THEN
2121
2122	DO ipts = 1,npts
2123
2124	biomass=SUM(SUM(circ_class_biomass(ipts,ivm,:,:,icarbon),2))
2125	ind=SUM(circ_class_n(ipts,ivm,:))
2126
2127	! If the flag is already true, don't change it
2128	IF(.NOT.lpft_replant(ipts,ivm))THEN
2129
2130	! There is no veget_max. Check whether these PFTs are indeed
2131	! empty. If not we have a serious problem
2132	IF(veget_max(ipts,ivm).LT.min_stomate)THEN
2133
2134	IF(biomass.GT.min_stomate .OR. ind.GT.min_stomate)THEN
2135
2136	WRITE(numout,*) 'ERROR - no veget_max but biomass and/or n',ipts,ivm
2137	WRITE(numout,*) 'biomass, ind, fm, ',biomass,&
2138	ind, forest_managed(ipts,ivm)
2139	error = error+1
2140
2141	ELSEIF(biomass.LT.min_stomate .AND. ind .LT.min_stomate)THEN
2142
2143	! This is the situation we expect. Do nothing
2144
2145	ELSE
2146
2147	! May be we overlooked a condition
2148	WRITE(numout,*) 'ERROR - overlooked a possible case 1'
2149	WRITE(numout,*) 'biomass, n, fm, ',biomass,&
2150	ind, forest_managed(ipts,ivm)
2151	error = error+1
2152
2153	ENDIF
2154
2155	ELSEIF(veget_max(ipts,ivm).GT.min_stomate)THEN
2156
2157	! There is veget_max. This is where the interesting cases should be
2158	IF(biomass.GT.min_stomate .AND. ind.GT.min_stomate)THEN
2159
2160	! This is a simple case: there is still biomass and individuals
2161	! so the stand is alive and there is no opportunity to change
2162	! the species. Note that in sapiens_forestry ::forest_managed
2163	! may already have been changed in case that the desired fm is
2164	! conservation. It is already done in forestry because if
2165	! fm 2 to 4 is changed to conservation we will not harvest but
2166	! leave the biomass on site.
2167	lpft_replant(ipts,ivm) = .FALSE.
2168
2169	ELSEIF( biomass.LT.min_stomate .AND. ind.LT.min_stomate)THEN
2170
2171	! The stand was killed or harvest so if the PFT is under
2172	! human management here is an opportunity to replant.
2173	IF(forest_managed(ipts,ivm).EQ.ifm_none)THEN
2174
2175	! If the stand is under a conservation strategy the species will
2176	! not change. Stands under conservation will switch FM without
2177	! a loss of biomass. Hence, there is no chance to replant.
2178	! If the stand is under conservation and dies from natural causes
2179	! it will be replanted with the same species
2180	lpft_replant(ipts,ivm) = .FALSE.
2181
2182	! Debug
2183	IF(printlev_loc.GT.4) THEN
2184	WRITE(numout,*) 'Replant - no biomass & no ind, conservation',ipts,ivm
2185	WRITE(numout,*) 'biomass, ind, fm, ',biomass,ind, forest_managed(ipts,ivm)
2186	WRITE(numout,*) 'lpft_replant, ',lpft_replant(ipts,ivm)
2187	END IF
2188	! -
2189
2190	ELSEIF(forest_managed(ipts,ivm).NE.ifm_none)THEN
2191
2192	! The stand was managed. This is a chance to change the species.
2193	! May be we wont use the chance but that will only be checked at
2194	! the end of the year in sapiens_forestry_species_change (sapiens_forestry.f90)
2195	lpft_replant(ipts,ivm) = .TRUE.
2196
2197	! Debug
2198	IF(printlev_loc.GT.4) THEN
2199	WRITE(numout,*) 'Replant - no biomass & no ind, management',ipts,ivm
2200	WRITE(numout,*) 'biomass, ind, fm, ',biomass,ind, forest_managed(ipts,ivm)
2201	WRITE(numout,*) 'lpft_replant, ',lpft_replant(ipts,ivm)
2202
2203	END IF
2204	!-
2205
2206	ELSE
2207
2208	WRITE(numout,*) 'ERROR - overlooked a possible option 1', ipts, ivm
2209	WRITE(numout,*) 'biomass, ind, fm, ',biomass, ind, forest_managed(ipts,ivm)
2210	WRITE(numout,*) 'lpft_replant, ',lpft_replant(ipts,ivm)
2211	error = error+1
2212
2213	ENDIF
2214
2215	ELSE
2216
2217	WRITE(numout,*) 'ERROR - overlooked a possible case 2'
2218	WRITE(numout,*) 'biomass, ind, fm, ',biomass,ind, forest_managed(ipts,ivm)
2219	WRITE(numout,*) 'lpft_replant, ',lpft_replant(ipts,ivm)
2220	error = error+1
2221	ENDIF
2222
2223	ELSE
2224
2225	WRITE(numout,*) 'ERROR - overlooked a possible case 3'
2226	WRITE(numout,*) 'biomass, ind, fm, ',biomass,ind, forest_managed(ipts,ivm)
2227	WRITE(numout,*) 'lpft_replant, ',lpft_replant(ipts,ivm)
2228	error = error+1
2229
2230	ENDIF
2231
2232	ENDIF ! lpft_replant
2233
2234	ENDDO ! ivm
2235
2236	ENDIF ! is_tree
2237
2238	ENDDO ! ipts
2239
2240	IF(error.GT.zero)THEN
2241	CALL ipslerr_p(3,'sapiens_forestry_flag_species_change','ERROR: this subroutine did not function as expected','','')
2242	ENDIF
2243
2244	END SUBROUTINE sapiens_forestry_flag_species_change
2245
2246
2247	!! ================================================================================================================================
2248	!! SUBROUTINE : sapiens_forestry_read_fm
2249	!!
2250	!>\BRIEF Read in a map that gives the forest management strategy
2251	!! for each pixel and each PFT.
2252	!!
2253	!! DESCRIPTION :
2254	!!
2255	!!
2256	!! FM = 1 : No human intervention (ORCHIDEE default)
2257	!! 2 : Thinnings based on the RDI, clearcuts based on tree density,
2258	!! annual increment, and tree diameter. Thinnings from above and
2259	!! from below are determined by the sign of thstrat.
2260	!! 3 : Coppices
2261	!! 4 : Short rotation coppices
2262	!!
2263	!! NOTE: This routine was mostly copied from slowproc where the PFTmap is read in.
2264	!! Grid interpolation is used, but only to look at the nearby pixels to see
2265	!! see which management strategy is dominant.
2266	!!
2267	!! RECENT CHANGE(S) : None
2268	!!
2269	!! MAIN OUTPUT VARIABLE(S): ::forest_managed
2270	!!
2271	!! REFERENCE(S) :
2272	!!
2273	!! FLOWCHART :
2274	!! \n
2275	!_ ================================================================================================================================
2276
2277	SUBROUTINE sapiens_forestry_read_fm ( npts, lalo, neighbours, resolution, contfrac, forest_managed )
2278
2279	!! 0. Variable and parameter declaration
2280
2281	!! 0.1 Input variables
2282	INTEGER(i_std), INTENT(in) :: npts !! Domain size - number of pixels
2283	!! (dimensionless)
2284	REAL(r_std), DIMENSION(:,:), INTENT(in) :: lalo !! Vector of latitude and longitudes (beware of the order !)
2285	INTEGER(i_std), DIMENSION(:,:), INTENT(in) :: neighbours !!
2286	REAL(r_std), DIMENSION(:,:), INTENT(in) :: resolution !! The size in m of each grid-box in X and Y
2287	REAL(r_std), DIMENSION(:), INTENT(in) :: contfrac !! Fraction of continent in the grid
2288
2289	!! 0.2 Output
2290	INTEGER(i_std), DIMENSION(:,:), INTENT(out) :: forest_managed !! Forest management flag: 0 = orchidee
2291	!! standard, 1= self-thinning only, 2=
2292	!! high-stand, 3= high-stand smoothed, 4=
2293	!! coppices
2294
2295	!! 0.3 Modified fields
2296
2297
2298	!! 0.4 Local variables
2299	CHARACTER(LEN=80) :: filename !! A string to hold the file name
2300	LOGICAL :: debug=.FALSE. !! A flag to print out debugging information.
2301	INTEGER(i_std) :: fid !! The ID of the NetCDF file.
2302	INTEGER(i_std) :: nb_coord !! The number of coordinates in the NetCDF file
2303	INTEGER(i_std) :: nb_gat !!
2304	INTEGER(i_std) :: nb_var !! The number of variables in the NetCDF file
2305	INTEGER(i_std) :: nb_dim !! The number of dimensions in the NetCDF file
2306	INTEGER(i_std) :: iml, jml, lml,ivm !! indices
2307	INTEGER(i_std) :: ip, inbv, jp !! indices
2308	LOGICAL :: l_ex !! A flag which indicates if a variable
2309	!! exists in the NetCDF file
2310	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_lu, lon_lu !! The latitude and longitude read in from
2311	!! the NetCDF file
2312	INTEGER,DIMENSION(flio_max_var_dims) :: l_d_w !! List of the dimension lengths of the variable
2313	!! in the NetCDF file
2314	INTEGER(i_std) :: ipts !! index
2315	INTEGER(i_std) :: ALLOC_ERR !! A flag tripped if we have an error in allocation
2316	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: fmmap_r !! The map read in from the NetCDF file
2317	INTEGER(i_std),DIMENSION(:,:,:),ALLOCATABLE :: fmmap_i !! The integer form of the map read in
2318	INTEGER(i_std) :: closest_lat !! The index of the closest latitude we found.
2319	INTEGER(i_std) :: closest_lon !! The index of the closest longitude we found.
2320	REAL(r_std) :: distance !! The distance from the current point to the
2321	!! point on the map
2322	REAL(r_std) :: closest_dist !! The distance to the closet point we've found.
2323	INTEGER :: large_int !! A number which indicates that the grid data
2324	!! is not available
2325	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lat_ful, lon_ful
2326	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: mask
2327	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area
2328	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: sub_index
2329	CHARACTER(LEN=30) :: callsign
2330	INTEGER(i_std) :: ibvm, nbvmax,ifm, ivma
2331	LOGICAL :: ok_interpol ! optionnal return of aggregate_2d
2332	REAL(r_std) :: sum_fm,sumf
2333	REAL(r_std),DIMENSION(nfm_types) :: fm_sum
2334	LOGICAL :: ltemp ! temporary logical variable
2335	!_ ================================================================================================================================
2336
2337	IF ( printlev >= 4 ) WRITE(numout,*) 'Entering sapiens_forestry_read_fm'
2338
2339	! This integer has to be large enough that it never shows up on the map, without being
2340	! so large that is causes overflows. Since none of the points on the map should be
2341	! larger than the number of FM strategies we have (nfm_types), this is sufficiently big.
2342	large_int = nvmap*nfm_types+1
2343
2344	!
2345	!Config Key = FM_FILE
2346	!Config Desc = Name of file from which the forest management map is to be read
2347	!Config If = OK_STOMATE
2348	!Config Def = FMmap.nc
2349	!Config Help = The name of the file to be opened to read a forest management
2350	!Config map (including a layer for every PFT) is given here.
2351	!Config Units = [FILE]
2352	filename = 'FMmap.nc'
2353	CALL getin_p('FM_FILE',filename)
2354
2355	IF (is_root_prc) THEN
2356	IF (debug) THEN
2357	WRITE(numout,*) "Entering sapiens_forestry_read_fm. Debug mode."
2358	WRITE (*,'(/," --> fliodmpf")')
2359	CALL fliodmpf (TRIM(filename))
2360	WRITE (*,'(/," --> flioopfd")')
2361	ENDIF
2362	CALL flioopfd (TRIM(filename),fid,nb_dim=nb_coord,nb_var=nb_var,nb_gat=nb_gat)
2363	IF (debug) THEN
2364	WRITE (*,'(" Number of coordinate in the file : ",I2)') nb_coord
2365	WRITE (*,'(" Number of variables in the file : ",I2)') nb_var
2366	WRITE (*,'(" Number of global attributes in the file : ",I2)') nb_gat
2367	ENDIF
2368	ENDIF
2369	CALL bcast(nb_coord)
2370	CALL bcast(nb_var)
2371	CALL bcast(nb_gat)
2372
2373	! This finds the number of longitude points in the file.
2374	IF (is_root_prc) &
2375	CALL flioinqv (fid,v_n="lon",l_ex=l_ex,nb_dims=nb_dim,len_dims=l_d_w)
2376	CALL bcast(l_d_w)
2377	iml=l_d_w(1)
2378	WRITE(numout,*) "FM Map: iml =",iml
2379
2380	! This finds the number of latitude points in the file.
2381	IF (is_root_prc) &
2382	CALL flioinqv (fid,v_n="lat",l_ex=l_ex,nb_dims=nb_dim,len_dims=l_d_w)
2383	CALL bcast(l_d_w)
2384	jml=l_d_w(1)
2385	WRITE(numout,*) "FM Map: jml =",jml
2386
2387	! Now find the number of PFTs in the file. If this is not equal to the number
2388	! of PFTs that we actually have, that's a problem and we'll crash.
2389	IF (is_root_prc) &
2390	CALL flioinqv (fid,v_n="FM_STRAT",l_ex=l_ex,nb_dims=nb_dim,len_dims=l_d_w)
2391	CALL bcast(l_d_w)
2392	lml=l_d_w(3)
2393
2394	IF (lml /= nvmap) THEN
2395	WRITE(numout,*) 'lml = ',lml
2396	WRITE(numout,*) 'nvmap = ',nvmap
2397	WRITE(numout,*) 'Stopping. '
2398	CALL ipslerr_p (3,'sapiens_forestry', &
2399	& 'Problem with forest management strategy map.','lml /= nvmap', &
2400	& '(number of pft must be equal)')
2401	ENDIF
2402	!
2403
2404	! Allocate the map that will be read in
2405	WRITE(numout,*) 'Reading the Forest Management strategy file'
2406	!
2407	ALLOC_ERR=-1
2408	ALLOCATE(fmmap_r(iml,jml,nvmap), STAT=ALLOC_ERR)
2409	IF (ALLOC_ERR/=0) THEN
2410	WRITE(numout,*) "ERROR IN ALLOCATION of fmmap_r : ",ALLOC_ERR
2411
2412	ENDIF
2413	ALLOC_ERR=-1
2414	ALLOCATE(fmmap_i(iml,jml,nvmap), STAT=ALLOC_ERR)
2415	IF (ALLOC_ERR/=0) THEN
2416	WRITE(numout,*) "ERROR IN ALLOCATION of fmmap_i : ",ALLOC_ERR
2417	CALL ipslerr_p (3,'sapiens_forestry_read_fm', 'error in fmmpap_i','','')
2418	ENDIF
2419
2420	! This reads in the map that is in the file
2421	IF (is_root_prc) THEN
2422	fmmap_r(:,:,:)=large_int*2.0
2423	CALL fliogetv (fid,"FM_STRAT", fmmap_r, start=(/ 1, 1, 1 /), count=(/ iml, jml, nvmap /))
2424	! Right now the values are all real, but they should be integers.
2425	! Careful, NINT might not work if the precision of fmmap_r is not single.
2426	! In that case, IDNINT should be used.
2427	DO ip=1,iml
2428	DO jp=1,jml
2429	DO ivma=1,nvmap
2430	! There will be some fill values in here. If we pass
2431	! a large fill value to NINT, it crashes. So let's
2432	! test for it
2433	#ifdef __NAGFOR
2434	ltemp=IEEE_IS_NAN(fmmap_r(ip,jp,ivma))
2435	#else
2436	ltemp=isnan(fmmap_r(ip,jp,ivma))
2437	#endif
2438	IF(ltemp)THEN
2439	fmmap_i(ip,jp,ivma)=large_int
2440	ELSE
2441	IF(fmmap_r(ip,jp,ivma) .GE. 0.0 .AND. fmmap_r(ip,jp,ivma) < large_int)THEN
2442	fmmap_i(ip,jp,ivma)=NINT(fmmap_r(ip,jp,ivma))
2443	ELSE
2444	! This value should be big enough that we don't barl ourselves
2445	! below.
2446	fmmap_i(ip,jp,ivma)=large_int
2447	ENDIF
2448
2449	ENDIF
2450	ENDDO
2451	ENDDO
2452	ENDDO
2453	ENDIF
2454
2455	CALL bcast(fmmap_i)
2456
2457	! Now I need to get the latitude and longitude
2458	! First, get the axes from the map file.
2459	ALLOC_ERR=-1
2460	ALLOCATE(lat_lu(jml), STAT=ALLOC_ERR)
2461	IF (ALLOC_ERR/=0) THEN
2462	WRITE(numout,*) "ERROR IN ALLOCATION of lat_lu : ",ALLOC_ERR
2463	CALL ipslerr_p (3,'sapiens_forestry_read_fm', 'error in lat_lu','','')
2464	ENDIF
2465	ALLOC_ERR=-1
2466	ALLOCATE(lon_lu(iml), STAT=ALLOC_ERR)
2467	IF (ALLOC_ERR/=0) THEN
2468	WRITE(numout,*) "ERROR IN ALLOCATION of lon_lu : ",ALLOC_ERR
2469	CALL ipslerr_p (3,'sapiens_forestry_read_fm', 'error in lon_lu','','')
2470	ENDIF
2471	IF (is_root_prc) THEN
2472	CALL fliogstc (fid, x_axis=lon_lu,y_axis=lat_lu)
2473	ENDIF
2474	CALL bcast(lon_lu)
2475	CALL bcast(lat_lu)
2476
2477	! Now I can interpolate. Remember that we are dealing with integer
2478	! values and management strategies. Therefore, we cannot do a strict
2479	! interpolation, since that might gives values of the strategy to be
2480	! 2.3, or something ridiculous. We cannot just round the number, either,
2481	! since we could have some squares with 1 (no management) and some with
2482	! 3 (coppices); An average of that would be 2 (high stands), which doesn't
2483	! make any sense. So we look to see what the most prevalent type of
2484	! management of nearby pixels is, and then just take that.
2485	ALLOC_ERR=-1
2486	ALLOCATE(lat_ful(iml,jml), STAT=ALLOC_ERR)
2487	IF (ALLOC_ERR/=0) THEN
2488	WRITE(numout,*) "ERROR IN ALLOCATION of lat_ful : ",ALLOC_ERR
2489	CALL ipslerr_p (3,'sapiens_forestry_read_fm', 'error in lat_ful','','')
2490	ENDIF
2491	ALLOC_ERR=-1
2492	ALLOCATE(lon_ful(iml,jml), STAT=ALLOC_ERR)
2493	IF (ALLOC_ERR/=0) THEN
2494	WRITE(numout,*) "ERROR IN ALLOCATION of lon_ful : ",ALLOC_ERR
2495	CALL ipslerr_p (3,'sapiens_forestry_read_fm', 'error in lon_ful','','')
2496	ENDIF
2497	!
2498	DO ip=1,iml
2499	lon_ful(ip,:)=lon_lu(ip)
2500	ENDDO
2501	DO jp=1,jml
2502	lat_ful(:,jp)=lat_lu(jp)
2503	ENDDO
2504	!
2505	! Mask of permitted variables.
2506	!
2507	ALLOC_ERR=-1
2508	ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR)
2509	IF (ALLOC_ERR/=0) THEN
2510	WRITE(numout,*) "ERROR IN ALLOCATION of mask : ",ALLOC_ERR
2511	CALL ipslerr_p (3,'sapiens_forestry_read_fm', 'error in mask','','')
2512	ENDIF
2513	!
2514	mask(:,:) = 0
2515	DO ip=1,iml
2516	DO jp=1,jml
2517	! If at least one PFT has a management strategy, we should interpolate
2518	! this grid point.
2519	sum_fm=SUM(fmmap_i(ip,jp,:))
2520	IF ( sum_fm .GT. min_sechiba .AND. sum_fm .LT. large_int) THEN
2521	mask(ip,jp) = 1
2522	IF (debug) THEN
2523	WRITE(numout,*) "update : SUM(fmmap(",ip,jp,")) = ",sum_fm
2524	ENDIF
2525	ENDIF
2526	ENDDO
2527	ENDDO
2528	!
2529	!
2530	! The number of maximum vegetation map points in the GCM grid should
2531	! also be computed and not imposed here.
2532	!
2533	nbvmax = 200
2534	!
2535	callsign="Forest Management map"
2536	!
2537	ok_interpol = .FALSE.
2538	DO WHILE ( .NOT. ok_interpol )
2539	WRITE(numout,*) "Projection arrays for ",callsign," : "
2540	WRITE(numout,*) "nbvmax = ",nbvmax
2541	!
2542	ALLOC_ERR=-1
2543	ALLOCATE(sub_index(npts, nbvmax,2), STAT=ALLOC_ERR)
2544	IF (ALLOC_ERR/=0) THEN
2545	WRITE(numout,*) "ERROR IN ALLOCATION of sub_index : ",ALLOC_ERR
2546	CALL ipslerr_p (3,'sapiens_forestry_read_fm', 'error in sub_index','','')
2547	ENDIF
2548	sub_index(:,:,:)=0
2549
2550	ALLOC_ERR=-1
2551	ALLOCATE(sub_area(npts, nbvmax), STAT=ALLOC_ERR)
2552	IF (ALLOC_ERR/=0) THEN
2553	WRITE(numout,*) "ERROR IN ALLOCATION of sub_area : ",ALLOC_ERR
2554	CALL ipslerr_p (3,'sapiens_forestry_read_fm', 'error in sub_area','','')
2555	ENDIF
2556	sub_area(:,:)=zero
2557	!
2558	CALL aggregate_p(npts, lalo, neighbours, resolution, contfrac, &
2559	& iml, jml, lon_ful, lat_ful, mask, callsign, &
2560	& nbvmax, sub_index, sub_area, ok_interpol)
2561	!
2562	IF ( .NOT. ok_interpol ) THEN
2563	DEALLOCATE(sub_area)
2564	DEALLOCATE(sub_index)
2565	ENDIF
2566	!
2567	nbvmax = nbvmax * 2
2568	ENDDO
2569
2570	! The FM maps with 28 PFTs for Europe were made on the exact grid that
2571	! is used in the FG4 simulations and therefore avoids the need for
2572	! interpolation. Moreover, the background value (used outside Europe and
2573	! within Europe where the PFT is not present) was no fm (=1). Hence, ifm
2574	! can never by zero (as it should be). Because all difficulties were
2575	! tackled when making the FM maps, reading the maps is straightforward.
2576	! The FM maps with 15 PFTs for the world were made on a 0.25 grid which
2577	! has a much higher resolution than most of the applications (FG1,
2578	! FG2, and FG3), hence, regridding will be necessary. Also to make better
2579	! looking maps land without the specific forest PFT got FM strategy zero.
2580	! Regridding comes with at least two risks: (1) the FM maps use the most
2581	! frequent FM strategy within the search area. This could be "no
2582	! management" because the minority of the pixels in the search area
2583	! contains forest. At the same time the PFT maps will also need to be
2584	! regridded but they use an average value. This way it seems possible that
2585	! we have inconsistencies between the PFT maps (forest present) and the FM
2586	! maps (no management strategy available -> fm = zero). (2) We could also
2587	! have used 1 as the background values. The maps would then look less
2588	! realistic but for the simulations itself it would make little difference
2589	! because the PFT is not present at the pixels were the background values
2590	! are set. When regridding this could result in an edge effect of having
2591	! some unmanaged pixels at the edge of a region with managed forests. The
2592	! code below is suitable for reading both types of FM maps (FG4 vs
2593	! FG1, FG2 and FG3).
2594	DO ipts = 1, npts
2595	! For this point, we need to see what dominant FM types are
2596	! nearby.
2597	sumf=zero
2598	! We need to do this for every PFT
2599	DO ivma=1,nvmap
2600	! If it's not a forest, we don't care.
2601	IF(.NOT. is_tree(start_index(ivma)))THEN
2602	forest_managed(ipts,ivma)=0
2603	ELSE
2604	fm_sum(:)=0.0
2605	DO ibvm=1, nbvmax
2606	! Leave the do loop if all sub areas are treated, sub_area <= 0
2607	IF ( sub_area(ipts,ibvm) <= zero ) EXIT
2608	ip = sub_index(ipts,ibvm,1)
2609	jp = sub_index(ipts,ibvm,2)
2610	ifm=fmmap_i(ip,jp,ivma)
2611	! Only calculate the surface area of the forested pixels.
2612	! This approach should ensure that the regridding of the
2613	! FM maps and the PFT maps will use exactly the same grid.
2614	IF (ifm.NE.zero) THEN
2615	fm_sum(ifm)=fm_sum(ifm)+sub_area(ipts,ibvm)
2616	END IF
2617	ENDDO
2618	! Whichever FM type has the most area, we use that one.
2619	IF (SUM(fm_sum(:)).EQ.zero) THEN
2620	! This is a pixel without forest. Because of the consistency
2621	! between the PFT and the FM maps it is safe to set FM to zero.
2622	! The biggest advantage of doing so is that the variable
2623	! forest_managed should show an intuitive map when being plotted.
2624	forest_managed(ipts,ivma) = zero
2625	ELSE
2626	forest_managed(ipts,ivma)=MAXLOC(fm_sum(:),1)
2627	ENDIF
2628	ENDIF
2629	ENDDO
2630	ENDDO
2631
2632	WRITE(numout,*) 'Done with interpolating the FM map'
2633
2634	!
2635	DEALLOCATE(fmmap_i)
2636	DEALLOCATE(fmmap_r)
2637	DEALLOCATE(lat_lu,lon_lu)
2638	DEALLOCATE(lat_ful,lon_ful)
2639	DEALLOCATE(mask)
2640	IF(ALLOCATED(sub_index)) DEALLOCATE(sub_index)
2641	IF(ALLOCATED(sub_area)) DEALLOCATE(sub_area)
2642
2643	IF ( printlev >= 5 ) WRITE(numout,*) 'Leaving sapiens_forestry_read_fm'
2644
2645	END SUBROUTINE sapiens_forestry_read_fm
2646
2647	!! ================================================================================================================================
2648	!! SUBROUTINE : sapiens_forestry_read_spinup_clearcut
2649	!!
2650	!>\BRIEF Read in a map that gives the state of clearcut during spinup
2651	!! for each pixel and each PFT.
2652	!!
2653	!! DESCRIPTION :
2654	!!
2655	!!
2656	!!
2657	!! NOTE: This routine was mostly copied from slowproc where the PFTmap is read in.
2658	!! Grid interpolation is used, but only to look at the nearby pixels to see
2659	!! see which management strategy is dominant.
2660	!!
2661	!! RECENT CHANGE(S) : None
2662	!!
2663	!! MAIN OUTPUT VARIABLE(S): ::spinup_clearcut
2664	!!
2665	!! REFERENCE(S) :
2666	!!
2667	!! FLOWCHART :
2668	!! \n
2669	!_ ================================================================================================================================
2670
2671	SUBROUTINE sapiens_forestry_read_spinup_clearcut ( npts, lalo, neighbours, resolution, contfrac, spinup_clearcut )
2672
2673	!! 0. Variable and parameter declaration
2674
2675	!! 0.1 Input variables
2676	INTEGER(i_std), INTENT(in) :: npts !! Domain size - number of pixels
2677	!! (dimensionless)
2678	REAL(r_std), DIMENSION(:,:), INTENT(in) :: lalo !! Vector of latitude and longitudes (beware of the order !)
2679	INTEGER(i_std), DIMENSION(:,:), INTENT(in) :: neighbours !!
2680	REAL(r_std), DIMENSION(:,:), INTENT(in) :: resolution !! The size in m of each grid-box in X and Y
2681	REAL(r_std), DIMENSION(:), INTENT(in) :: contfrac !! Fraction of continent in the grid
2682
2683	!! 0.2 Output
2684	INTEGER(i_std), DIMENSION(:,:), INTENT(out) :: spinup_clearcut !! Forest management flag: 0 = orchidee
2685	!! standard, 1= self-thinning only, 2=
2686	!! high-stand, 3= high-stand smoothed, 4=
2687	!! coppices
2688
2689	!! 0.3 Modified fields
2690
2691
2692	!! 0.4 Local variables
2693	CHARACTER(LEN=80) :: filename !! A string to hold the file name
2694	LOGICAL :: debug=.FALSE. !! A flag to print out debugging information.
2695	INTEGER(i_std) :: fid !! The ID of the NetCDF file.
2696	INTEGER(i_std) :: nb_coord !! The number of coordinates in the NetCDF file
2697	INTEGER(i_std) :: nb_gat !!
2698	INTEGER(i_std) :: nb_var !! The number of variables in the NetCDF file
2699	INTEGER(i_std) :: nb_dim !! The number of dimensions in the NetCDF file
2700	INTEGER(i_std) :: iml, jml, lml,ivm !! indices
2701	INTEGER(i_std) :: ip, inbv, jp !! indices
2702	LOGICAL :: l_ex !! A flag which indicates if a variable
2703	!! exists in the NetCDF file
2704	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_lu, lon_lu !! The latitude and longitude read in from
2705	!! the NetCDF file
2706	INTEGER,DIMENSION(flio_max_var_dims) :: l_d_w !! List of the dimension lengths of the variable
2707	!! in the NetCDF file
2708	INTEGER(i_std) :: ipts !! index
2709	INTEGER(i_std) :: ALLOC_ERR !! A flag tripped if we have an error in allocation
2710	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: fmmap_r !! The map read in from the NetCDF file
2711	INTEGER(i_std),DIMENSION(:,:,:),ALLOCATABLE :: fmmap_i !! The integer form of the map read in
2712	INTEGER(i_std) :: closest_lat !! The index of the closest latitude we found.
2713	INTEGER(i_std) :: closest_lon !! The index of the closest longitude we found.
2714	REAL(r_std) :: distance !! The distance from the current point to the
2715	!! point on the map
2716	REAL(r_std) :: closest_dist !! The distance to the closet point we've found.
2717	INTEGER :: large_int !! A number which indicates that the grid data
2718	!! is not available
2719	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lat_ful, lon_ful
2720	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: mask
2721	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area
2722	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: sub_index
2723	CHARACTER(LEN=30) :: callsign
2724	INTEGER(i_std) :: ibvm, nbvmax,ifm, ivma
2725	LOGICAL :: ok_interpol ! optionnal return of aggregate_2d
2726	REAL(r_std) :: sum_fm,sumf
2727	LOGICAL :: ltemp ! temporary logical variable
2728	!_ ================================================================================================================================
2729
2730	IF ( printlev >= 4 ) WRITE(numout,*) 'Entering sapiens_forestry_read_spinup_clearcut'
2731
2732	! This integer has to be large enough that it never shows up on the map, without being
2733	! so large that is causes overflows. Since none of the points on the map should be
2734	! larger than the number of FM strategies we have (nfm_types), this is sufficiently big.
2735	large_int = flag_spinup_clearcut + 1
2736
2737	!
2738	!Config Key = FM_FILE
2739	!Config Desc = Name of file to be read
2740	!Config If = OK_STOMATE
2741	!Config Def = FMmap.nc
2742	!Config Help = The name of the file to be opened to read a forest management
2743	!Config map (including a layer for every PFT) is given here.
2744	!Config Units = [FILE]
2745	filename = 'spinup_clearcut.nc'
2746	CALL getin_p('SPINUP_CLEARCUT_FILE',filename)
2747
2748	IF (is_root_prc) THEN
2749	IF (debug) THEN
2750	WRITE(numout,*) "Entering sapiens_forestry_read_spinup_clearcut. Debug mode."
2751	WRITE (*,'(/," --> fliodmpf")')
2752	CALL fliodmpf (TRIM(filename))
2753	WRITE (*,'(/," --> flioopfd")')
2754	ENDIF
2755	CALL flioopfd (TRIM(filename),fid,nb_dim=nb_coord,nb_var=nb_var,nb_gat=nb_gat)
2756	IF (debug) THEN
2757	WRITE (*,'(" Number of coordinate in the file : ",I2)') nb_coord
2758	WRITE (*,'(" Number of variables in the file : ",I2)') nb_var
2759	WRITE (*,'(" Number of global attributes in the file : ",I2)') nb_gat
2760	ENDIF
2761	ENDIF
2762	CALL bcast(nb_coord)
2763	CALL bcast(nb_var)
2764	CALL bcast(nb_gat)
2765
2766	! This finds the number of longitude points in the file.
2767	IF (is_root_prc) &
2768	CALL flioinqv (fid,v_n="lon",l_ex=l_ex,nb_dims=nb_dim,len_dims=l_d_w)
2769	CALL bcast(l_d_w)
2770	iml=l_d_w(1)
2771	WRITE(numout,*) "spinup_clearcut map: iml =",iml
2772
2773	! This finds the number of latitude points in the file.
2774	IF (is_root_prc) &
2775	CALL flioinqv (fid,v_n="lat",l_ex=l_ex,nb_dims=nb_dim,len_dims=l_d_w)
2776	CALL bcast(l_d_w)
2777	jml=l_d_w(1)
2778	WRITE(numout,*) "spinup_clearcut map: jml =",jml
2779
2780	! Now find the number of PFTs in the file. If this is not equal to the number
2781	! of PFTs that we actually have, that's a problem and we'll crash.
2782	IF (is_root_prc) &
2783	CALL flioinqv (fid,v_n="Clearcut",l_ex=l_ex,nb_dims=nb_dim,len_dims=l_d_w)
2784	CALL bcast(l_d_w)
2785	lml=l_d_w(3)
2786
2787	IF (lml /= nvmap) THEN
2788	WRITE(numout,*) 'lml = ',lml
2789	WRITE(numout,*) 'nvmap = ',nvmap
2790	WRITE(numout,*) 'Stopping. '
2791	CALL ipslerr_p (3,'sapiens_forestry', &
2792	& 'Problem with spinup clearcut map.','lml /= nvmap', &
2793	& '(number of pft must be equal)')
2794	ENDIF
2795	!
2796
2797	! Allocate the map that will be read in
2798	WRITE(numout,*) 'Reading the spinup clearcut strategy file'
2799	!
2800	ALLOC_ERR=-1
2801	ALLOCATE(fmmap_r(iml,jml,nvmap), STAT=ALLOC_ERR)
2802	IF (ALLOC_ERR/=0) THEN
2803	WRITE(numout,*) "ERROR IN ALLOCATION of fmmap_r : ",ALLOC_ERR
2804
2805	ENDIF
2806	ALLOC_ERR=-1
2807	ALLOCATE(fmmap_i(iml,jml,nvmap), STAT=ALLOC_ERR)
2808	IF (ALLOC_ERR/=0) THEN
2809	WRITE(numout,*) "ERROR IN ALLOCATION of fmmap_i : ",ALLOC_ERR
2810	CALL ipslerr_p (3,'sapiens_forestry_read_spinup_clearcut', 'error in fmmpap_i','','')
2811	ENDIF
2812
2813	! This reads in the map that is in the file
2814	IF (is_root_prc) THEN
2815	fmmap_r(:,:,:)=zero
2816	CALL fliogetv (fid,"Clearcut", fmmap_r, start=(/ 1, 1, 1 /), count=(/ iml, jml, nvmap /))
2817	! Right now the values are all real, but they should be integers.
2818	! Careful, NINT might not work if the precision of fmmap_r is not single.
2819	! In that case, IDNINT should be used.
2820	DO ip=1,iml
2821	DO jp=1,jml
2822	DO ivma=1,nvmap
2823	! There will be some fill values in here. If we pass
2824	! a large fill value to NINT, it crashes. So let's
2825	! test for it
2826	#ifdef __NAGFOR
2827	ltemp=IEEE_IS_NAN(fmmap_r(ip,jp,ivma))
2828	#else
2829	ltemp=isnan(fmmap_r(ip,jp,ivma))
2830	#endif
2831	IF(ltemp)THEN
2832	fmmap_i(ip,jp,ivma)=zero
2833	ELSE
2834	IF(fmmap_r(ip,jp,ivma) .GE. 0.0 .AND. fmmap_r(ip,jp,ivma) .LE. flag_spinup_clearcut)THEN
2835	fmmap_i(ip,jp,ivma)=NINT(fmmap_r(ip,jp,ivma))
2836	ELSE
2837	! For the moment there is a strict rule that values of input
2838	! cannot be bigger than 1.
2839	CALL ipslerr_p (3,'sapiens_forestry_read_spinup_clearcut', 'value is bigger than 1','','')
2840	ENDIF
2841
2842	ENDIF
2843	ENDDO
2844	ENDDO
2845	ENDDO
2846	ENDIF
2847
2848	CALL bcast(fmmap_i)
2849
2850	! Now I need to get the latitude and longitude
2851	! First, get the axes from the map file.
2852	ALLOC_ERR=-1
2853	ALLOCATE(lat_lu(jml), STAT=ALLOC_ERR)
2854	IF (ALLOC_ERR/=0) THEN
2855	WRITE(numout,*) "ERROR IN ALLOCATION of lat_lu : ",ALLOC_ERR
2856	CALL ipslerr_p (3,'sapiens_forestry_read_spinup_clearcut', 'error in lat_lu','','')
2857	ENDIF
2858	ALLOC_ERR=-1
2859	ALLOCATE(lon_lu(iml), STAT=ALLOC_ERR)
2860	IF (ALLOC_ERR/=0) THEN
2861	WRITE(numout,*) "ERROR IN ALLOCATION of lon_lu : ",ALLOC_ERR
2862	CALL ipslerr_p (3,'sapiens_forestry_read_spinup_clearcut', 'error in lon_lu','','')
2863	ENDIF
2864	IF (is_root_prc) THEN
2865	CALL fliogstc (fid, x_axis=lon_lu,y_axis=lat_lu)
2866	ENDIF
2867	CALL bcast(lon_lu)
2868	CALL bcast(lat_lu)
2869
2870	! Now I can interpolate. Remember that we are dealing with integer
2871	! values and management strategies. Therefore, we cannot do a strict
2872	! interpolation, since that might gives values of the strategy to be
2873	! 2.3, or something ridiculous. We cannot just round the number, either,
2874	! since we could have some squares with 1 (no management) and some with
2875	! 3 (coppices); An average of that would be 2 (high stands), which doesn't
2876	! make any sense. So we look to see what the most prevalent type of
2877	! management of nearby pixels is, and then just take that.
2878	ALLOC_ERR=-1
2879	ALLOCATE(lat_ful(iml,jml), STAT=ALLOC_ERR)
2880	IF (ALLOC_ERR/=0) THEN
2881	WRITE(numout,*) "ERROR IN ALLOCATION of lat_ful : ",ALLOC_ERR
2882	CALL ipslerr_p (3,'sapiens_forestry_read_spinup_clearcut', 'error in lat_ful','','')
2883	ENDIF
2884	ALLOC_ERR=-1
2885	ALLOCATE(lon_ful(iml,jml), STAT=ALLOC_ERR)
2886	IF (ALLOC_ERR/=0) THEN
2887	WRITE(numout,*) "ERROR IN ALLOCATION of lon_ful : ",ALLOC_ERR
2888	CALL ipslerr_p (3,'sapiens_forestry_read_spinup_clearcut', 'error in lon_ful','','')
2889	ENDIF
2890	!
2891	DO ip=1,iml
2892	lon_ful(ip,:)=lon_lu(ip)
2893	ENDDO
2894	DO jp=1,jml
2895	lat_ful(:,jp)=lat_lu(jp)
2896	ENDDO
2897	!
2898	! Mask of permitted variables.
2899	!
2900	ALLOC_ERR=-1
2901	ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR)
2902	IF (ALLOC_ERR/=0) THEN
2903	WRITE(numout,*) "ERROR IN ALLOCATION of mask : ",ALLOC_ERR
2904	CALL ipslerr_p (3,'sapiens_forestry_read_spinup_clearcut', 'error in mask','','')
2905	ENDIF
2906	!
2907	mask(:,:) = 1
2908	! DO ip=1,iml
2909	! DO jp=1,jml
2910	! ! If at least one PFT has a management strategy, we should interpolate
2911	! ! this grid point.
2912	! sum_fm=SUM(fmmap_i(ip,jp,:))
2913	! IF ( sum_fm .GT. min_sechiba .AND. sum_fm .LT. large_int) THEN
2914	! mask(ip,jp) = 1
2915	! IF (debug) THEN
2916	! WRITE(numout,*) "update : SUM(fmmap(",ip,jp,")) = ",sum_fm
2917	! ENDIF
2918	! ENDIF
2919	! ENDDO
2920	! ENDDO
2921	!
2922	!
2923	! The number of maximum vegetation map points in the GCM grid should
2924	! also be computed and not imposed here.
2925	!
2926	nbvmax = 200
2927	!
2928	callsign="Spinup clearcut map"
2929	!
2930	ok_interpol = .FALSE.
2931	DO WHILE ( .NOT. ok_interpol )
2932	WRITE(numout,*) "Projection arrays for ",callsign," : "
2933	WRITE(numout,*) "nbvmax = ",nbvmax
2934	!
2935	ALLOC_ERR=-1
2936	ALLOCATE(sub_index(npts, nbvmax,2), STAT=ALLOC_ERR)
2937	IF (ALLOC_ERR/=0) THEN
2938	WRITE(numout,*) "ERROR IN ALLOCATION of sub_index : ",ALLOC_ERR
2939	CALL ipslerr_p (3,'sapiens_forestry_read_spinup_clearcut', 'error in sub_index','','')
2940	ENDIF
2941	sub_index(:,:,:)=0
2942
2943	ALLOC_ERR=-1
2944	ALLOCATE(sub_area(npts, nbvmax), STAT=ALLOC_ERR)
2945	IF (ALLOC_ERR/=0) THEN
2946	WRITE(numout,*) "ERROR IN ALLOCATION of sub_area : ",ALLOC_ERR
2947	CALL ipslerr_p (3,'sapiens_forestry_read_spinup_clearcut', 'error in sub_area','','')
2948	ENDIF
2949	sub_area(:,:)=zero
2950	!
2951	CALL aggregate_p(npts, lalo, neighbours, resolution, contfrac, &
2952	& iml, jml, lon_ful, lat_ful, mask, callsign, &
2953	& nbvmax, sub_index, sub_area, ok_interpol)
2954	!
2955	IF ( .NOT. ok_interpol ) THEN
2956	DEALLOCATE(sub_area)
2957	DEALLOCATE(sub_index)
2958	ENDIF
2959	!
2960	nbvmax = nbvmax * 2
2961	ENDDO
2962
2963	! Obtain the majority clearcut flag value after regrdding.
2964	DO ipts = 1, npts
2965	! We need to do this for every PFT
2966	DO ivma=1,nvmap
2967	! For this point, we need to see what dominant FM types are
2968	! nearby.
2969	sumf=zero
2970	! If it's not a forest, we don't care.
2971	IF(.NOT. is_tree(start_index(ivma)))THEN
2972	spinup_clearcut(ipts,ivma)=0
2973	ELSE
2974	DO ibvm=1, nbvmax
2975	! Leave the do loop if all sub areas are treated, sub_area <= 0
2976	IF ( sub_area(ipts,ibvm) <= zero ) EXIT
2977	ip = sub_index(ipts,ibvm,1)
2978	jp = sub_index(ipts,ibvm,2)
2979	ifm=fmmap_i(ip,jp,ivma)
2980	sumf = sumf+ifm
2981	!WRITE(numout,*) 'jp: ',jp,'ip: ',ip,'ifm: ',ifm,'sumf: ',sumf,'ibmv: ',ibvm
2982	ENDDO
2983	! As long as any sub-pixel goes through clearcut, the coarse pixel
2984	! goes through clearcut.
2985	IF (sumf .GT. zero) THEN
2986	spinup_clearcut(ipts,ivma) = flag_spinup_clearcut
2987	ELSE
2988	spinup_clearcut(ipts,ivma) = zero
2989	ENDIF
2990	ENDIF
2991	ENDDO
2992	ENDDO
2993
2994	WRITE(numout,*) 'Done with interpolating the spinup clearcut map'
2995
2996	!
2997	DEALLOCATE(fmmap_i)
2998	DEALLOCATE(fmmap_r)
2999	DEALLOCATE(lat_lu,lon_lu)
3000	DEALLOCATE(lat_ful,lon_ful)
3001	DEALLOCATE(mask)
3002	IF(ALLOCATED(sub_index)) DEALLOCATE(sub_index)
3003	IF(ALLOCATED(sub_area)) DEALLOCATE(sub_area)
3004
3005	IF ( printlev >= 5 ) WRITE(numout,*) 'Leaving sapiens_forestry_read_spinup_clearcut'
3006
3007	END SUBROUTINE sapiens_forestry_read_spinup_clearcut
3008
3009	!! ================================================================================================================================
3010	!! SUBROUTINE : sapiens_forestry_read_litter
3011	!!
3012	!>\BRIEF Read in a map that gives the litter demand for each pixel.
3013	!!
3014	!! DESCRIPTION : Reads in a map for the litter demand of each pixel and interpolates
3015	!! to find the demand for all pixels that we are interested in. The map
3016	!! in this case has units of gC/year. We convert it to be gC/m**2/year, since
3017	!! all the other units in the model are given per square meter. The litter
3018	!! demand maps can be calculated based on, for example, the demand needed
3019	!! per head of livestock for storing them indoors in the winter. This practice
3020	!! reached its peak during the 19th century and rapidly declined afterwards
3021	!! in Europe. If you are not interested in historical simulations or improved
3022	!! estimation of soil carbon pools during a historical spinup, you probably
3023	!! don't need this.
3024	!!
3025	!!
3026	!! RECENT CHANGE(S) : None
3027	!!
3028	!! MAIN OUTPUT VARIABLE(S): ::litter_demand
3029	!!
3030	!! REFERENCE(S) :
3031	!!
3032	!! FLOWCHART :
3033	!! \n
3034	!_ ================================================================================================================================
3035
3036	SUBROUTINE sapiens_forestry_read_litter ( npts, lalo, neighbours, resolution, contfrac, litter_demand )
3037
3038	!! 0. Variable and parameter declaration
3039
3040	!! 0.1 Input variables
3041	INTEGER(i_std), INTENT(in) :: npts !! Domain size - number of pixels
3042	!! (dimensionless)
3043	REAL(r_std), DIMENSION(:,:), INTENT(in) :: lalo !! Vector of latitude and longitudes (beware of the order !)
3044	INTEGER(i_std), DIMENSION(:,:), INTENT(in) :: neighbours !!
3045	REAL(r_std), DIMENSION(:,:), INTENT(in) :: resolution !! The size in m of each grid-box in X and Y
3046	REAL(r_std), DIMENSION(:), INTENT(in) :: contfrac !! Fraction of continent in the grid
3047
3048	!! 0.2 Output
3049	REAL(r_std), DIMENSION(:), INTENT(out) :: litter_demand !! The litter removed from each pixel due
3050	!! to litter raking at the end of the year.
3051	!! @tex $(gC year^{-1})$ @endtex
3052
3053	!! 0.3 Modified fields
3054
3055
3056	!! 0.4 Local variables
3057	!
3058	! PARAMETERS...taken from grid.f90. I don't know why this variable is not in constants.f90.
3059	! default resolution (m)
3060	REAL(r_std), PARAMETER :: default_resolution = 250000.
3061
3062	CHARACTER(LEN=80) :: filename !! A string to hold the file name
3063	LOGICAL :: debug=.FALSE. !! A flag to print out debugging information.
3064	INTEGER(i_std) :: fid !! The ID of the NetCDF file.
3065	INTEGER(i_std) :: nb_coord !! The number of coordinates in the NetCDF file
3066	INTEGER(i_std) :: nb_gat !!
3067	INTEGER(i_std) :: nb_var !! The number of variables in the NetCDF file
3068	INTEGER(i_std) :: nb_dim !! The number of dimensions in the NetCDF file
3069	INTEGER(i_std) :: iml, jml, lml,ivmi,xxx !! indices
3070	INTEGER(i_std) :: ip, inbv, jp !! indices
3071	LOGICAL :: l_ex !! A flag which indicates if a variable
3072	!! exists in the NetCDF file
3073	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_lu, lon_lu !! The latitude and longitude read in from
3074	!! the NetCDF file
3075	INTEGER,DIMENSION(flio_max_var_dims) :: l_d_w !! List of the dimension lengths of the variable
3076	!! in the NetCDF file
3077	INTEGER(i_std) :: ipts !! index
3078	INTEGER(i_std) :: ALLOC_ERR !! A flag tripped if we have an error in allocation
3079	REAL(r_std),DIMENSION(:,:),ALLOCATABLE :: littermap !! The map read in from the NetCDF file
3080	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lat_ful, lon_ful
3081	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: mask
3082	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area
3083	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: sub_index
3084	CHARACTER(LEN=30) :: callsign
3085	INTEGER(i_std) :: ibvm, nbvmax,ifm
3086	LOGICAL :: ok_interpol ! optionnal return of aggregate_2d
3087	REAL(r_std) :: area_sum, coslat
3088	REAL(r_std),DIMENSION(nfm_types) :: fm_sum
3089	REAL(r_std), ALLOCATABLE, DIMENSION(:,:,:) :: temp_resolution
3090	REAL(r_std) :: fillvalue
3091	INTEGER :: i_rc,i_v,no_fill
3092	LOGICAL :: ltemp ! temporary logical variable
3093	!_ ================================================================================================================================
3094
3095	IF ( printlev >= 4 ) WRITE(numout,*) 'Entering sapiens_forestry_read_litter'
3096
3097	litter_demand(:)=zero
3098
3099
3100	!
3101	!Config Key = LITTER_FILE
3102	!Config Desc = Name of file from which the litter raking map is to be read
3103	!Config If = OK_STOMATE
3104	!Config Def = litter_map.nc
3105	!Config Help =
3106	!Config Units = [FILE]
3107	!
3108	filename = 'litter_map.nc'
3109	CALL getin_p('LITTER_FILE',filename)
3110
3111	WRITE(numout,*)'reading the litter_map'
3112
3113	IF (is_root_prc) THEN
3114	IF (debug) THEN
3115	WRITE(numout,*) "Entering sapiens_forestry_read_litter. Debug mode."
3116	WRITE (*,'(/," --> fliodmpf")')
3117	CALL fliodmpf (TRIM(filename))
3118	WRITE (*,'(/," --> flioopfd")')
3119	ENDIF
3120	CALL flioopfd (TRIM(filename),fid,nb_dim=nb_coord,nb_var=nb_var,nb_gat=nb_gat)
3121	IF (debug) THEN
3122	WRITE (*,'(" Number of coordinate in the file : ",I2)') nb_coord
3123	WRITE (*,'(" Number of variables in the file : ",I2)') nb_var
3124	WRITE (*,'(" Number of global attributes in the file : ",I2)') nb_gat
3125	ENDIF
3126	ENDIF
3127	CALL bcast(nb_coord)
3128	CALL bcast(nb_var)
3129	CALL bcast(nb_gat)
3130
3131	! I need the fill value of the map. This doesn't seem to available through
3132	! the IOISPL routines that I can see.
3133	! NF90_INQ_VAR_FILL seems to be only valid for NetCDF-'? So I am going to make
3134	! a HUGE assumption here about the fill value of the maps. I hope that none of the values
3135	! in the maps are above this value.
3136	fillvalue=1e20
3137	WRITE(numout,*) "WARNING: You are using a litter demand map. I cannot determine"
3138	WRITE(numout,*) " the fill value of this map. I am assuming that the fill"
3139	WRITE(numout,*) " value is 1e20 and that all your real values are below that."
3140	WRITE(numout,*) " CONFIRM THIS WITH YOUR LITTER MAP!"
3141
3142
3143	! This finds the number of longitude points in the file.
3144	IF (is_root_prc) &
3145	CALL flioinqv (fid,v_n="lon",l_ex=l_ex,nb_dims=nb_dim,len_dims=l_d_w)
3146	CALL bcast(l_d_w)
3147	iml=l_d_w(1)
3148	WRITE(numout,*) "Litter map: iml =",iml
3149
3150	! This finds the number of latitude points in the file.
3151	IF (is_root_prc) &
3152	CALL flioinqv (fid,v_n="lat",l_ex=l_ex,nb_dims=nb_dim,len_dims=l_d_w)
3153	CALL bcast(l_d_w)
3154	jml=l_d_w(1)
3155	WRITE(numout,*) "Litter map: jml =",jml
3156
3157	! Now find the number of PFTs in the file. If this is not equal to the number
3158	! of PFTs that we actually have, that's a problem and we'll crash.
3159	IF (is_root_prc) &
3160	CALL flioinqv (fid,v_n="LITTER_DEMAND",l_ex=l_ex,nb_dims=nb_dim,len_dims=l_d_w)
3161	CALL bcast(l_d_w)
3162
3163	!
3164
3165	! Allocate the map that will be read in
3166	WRITE(numout,*) 'Reading the litter file',l_d_w
3167	!
3168	ALLOC_ERR=-1
3169	ALLOCATE(littermap(iml,jml), STAT=ALLOC_ERR)
3170	IF (ALLOC_ERR/=0) THEN
3171	WRITE(numout,*) "ERROR IN ALLOCATION of littermap : ",ALLOC_ERR
3172	CALL ipslerr_p (3,'sapiens_forestry_read_litter', 'error in littermap','','')
3173	ENDIF
3174	ALLOC_ERR=-1
3175
3176	! This reads in the map that is in the file
3177	IF (is_root_prc) THEN
3178	CALL fliogetv (fid,"LITTER_DEMAND", littermap, start=(/ 1, 1 /), count=(/ iml, jml /))
3179	ENDIF
3180
3181	CALL bcast(littermap)
3182
3183	! Now I need to get the latitude and longitude
3184	! First, get the axes from the map file.
3185	ALLOC_ERR=-1
3186	ALLOCATE(lat_lu(jml), STAT=ALLOC_ERR)
3187	IF (ALLOC_ERR/=0) THEN
3188	WRITE(numout,*) "ERROR IN ALLOCATION of lat_lu : ",ALLOC_ERR
3189	CALL ipslerr_p (3,'sapiens_forestry_read_litter', 'error in lat_lu','','')
3190	ENDIF
3191	ALLOC_ERR=-1
3192	ALLOCATE(lon_lu(iml), STAT=ALLOC_ERR)
3193	IF (ALLOC_ERR/=0) THEN
3194	WRITE(numout,*) "ERROR IN ALLOCATION of lon_lu : ",ALLOC_ERR
3195	CALL ipslerr_p (3,'sapiens_forestry_read_litter', 'error in lon_lu','','')
3196	ENDIF
3197	IF (is_root_prc) THEN
3198	CALL fliogstc (fid, x_axis=lon_lu,y_axis=lat_lu)
3199	ENDIF
3200	! Do we need to close the file??
3201	!ASLA IF (is_root_prc) THEN
3202	!ASLA CALL flinclo(fid)
3203	!ASLA ENDIF
3204	CALL bcast(lon_lu)
3205	CALL bcast(lat_lu)
3206
3207	! Now I can interpolate. This should be straightforward, since
3208	! we are interpolating real numbers.
3209
3210	ALLOC_ERR=-1
3211	ALLOCATE(lat_ful(iml,jml), STAT=ALLOC_ERR)
3212	IF (ALLOC_ERR/=0) THEN
3213	WRITE(numout,*) "ERROR IN ALLOCATION of lat_ful : ",ALLOC_ERR
3214	CALL ipslerr_p (3,'sapiens_forestry_read_litter', 'error in lat_ful','','')
3215	ENDIF
3216	ALLOC_ERR=-1
3217	ALLOCATE(lon_ful(iml,jml), STAT=ALLOC_ERR)
3218	IF (ALLOC_ERR/=0) THEN
3219	WRITE(numout,*) "ERROR IN ALLOCATION of lon_ful : ",ALLOC_ERR
3220	CALL ipslerr_p (3,'sapiens_forestry_read_litter', 'error in lon_ful','','')
3221	ENDIF
3222	!
3223	DO ip=1,iml
3224	lon_ful(ip,:)=lon_lu(ip)
3225	ENDDO
3226	DO jp=1,jml
3227	lat_ful(:,jp)=lat_lu(jp)
3228	ENDDO
3229
3230	! The map is in absolute numbers. We cannot interpolate absolute numbers, we want
3231	! to add them together (or subtract them, if we are going to finer resolution). Therefore,
3232	! let's convert all the map values to per m**2 first, interpolate that, and then
3233	! convert back to absolute numbers.
3234	! This is ugly. Unfortunately, we had to steal this from grid.f90, since the resolution is
3235	! not passed for every grid square, just those that we are interested in.
3236	ALLOC_ERR=-1
3237	ALLOCATE(temp_resolution(iml,jml,2), STAT=ALLOC_ERR)
3238	IF (ALLOC_ERR/=0) THEN
3239	WRITE(numout,*) "ERROR IN ALLOCATION of temp_resolution : ",ALLOC_ERR
3240	CALL ipslerr_p (3,'sapiens_forestry_read_litter', 'error in temp_resolution','','')
3241	ENDIF
3242	DO ip=1,iml
3243	DO jp=1,jml
3244	!
3245	! 2 resolution
3246	!
3247
3248	!
3249	! 2.1 longitude
3250	!
3251
3252	! prevent infinite resolution at the pole
3253	coslat = MAX( COS( lat_ful(ip,jp) * pi/180. ), mincos )
3254	IF ( iml .GT. 1 ) THEN
3255
3256	IF ( ip .EQ. 1 ) THEN
3257	temp_resolution(ip,jp,1) = &
3258	ABS( lon_ful(ip+1,jp) - lon_ful(ip,jp) ) * &
3259	pi/180. * R_Earth * coslat
3260	ELSEIF ( ip .EQ. iml ) THEN
3261	temp_resolution(ip,jp,1) = &
3262	ABS( lon_ful(ip,jp) - lon_ful(ip-1,jp) ) * &
3263	pi/180. * R_Earth * coslat
3264	ELSE
3265	temp_resolution(ip,jp,1) = &
3266	ABS( lon_ful(ip+1,jp) - lon_ful(ip-1,jp) )/2. *&
3267	pi/180. * R_Earth * coslat
3268	ENDIF
3269
3270	ELSE
3271
3272	temp_resolution(ip,jp,1) = default_resolution
3273
3274	ENDIF
3275
3276	!
3277	! 2.2 latitude
3278	!
3279
3280	IF ( jml .GT. 1 ) THEN
3281
3282	IF ( jp .EQ. 1 ) THEN
3283	temp_resolution(ip,jp,2) = &
3284	ABS( lat_ful(ip,jp) - lat_ful(ip,jp+1) ) * &
3285	pi/180. * R_Earth
3286	ELSEIF ( jp .EQ. jml ) THEN
3287	temp_resolution(ip,jp,2) = &
3288	ABS( lat_ful(ip,jp-1) - lat_ful(ip,jp) ) * &
3289	pi/180. * R_Earth
3290	ELSE
3291	temp_resolution(ip,jp,2) = &
3292	ABS( lat_ful(ip,jp-1) - lat_ful(ip,jp+1) )/2. *&
3293	pi/180. * R_Earth
3294	ENDIF
3295
3296	ELSE
3297
3298	temp_resolution(ip,jp,2) = default_resolution
3299
3300	ENDIF
3301
3302	#ifdef __NAGFOR
3303	ltemp=IEEE_IS_NAN(littermap(ip,jp))
3304	#else
3305	ltemp=isnan(littermap(ip,jp))
3306	#endif
3307	IF( .NOT. ltemp)THEN
3308
3309	IF(littermap(ip,jp) .LT. fillvalue)THEN
3310	! Convert to per square meter for the interpolation if this is a real map value.
3311	littermap(ip,jp)=littermap(ip,jp)/temp_resolution(ip,jp,1)/temp_resolution(ip,jp,2)
3312	ENDIF
3313
3314	ENDIF
3315
3316	ENDDO
3317
3318	ENDDO
3319
3320	!
3321	WRITE(numout,*) 'Reading the LITTER file'
3322	!
3323	!
3324
3325	!
3326	! Mask of permitted variables.
3327	!
3328	ALLOC_ERR=-1
3329	ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR)
3330	IF (ALLOC_ERR/=0) THEN
3331	WRITE(numout,*) "ERROR IN ALLOCATION of mask : ",ALLOC_ERR
3332	CALL ipslerr_p (3,'sapiens_forestry_read_litter', 'error in mask','','')
3333	ENDIF
3334	!
3335	mask(:,:) = 0
3336	DO ip=1,iml
3337	DO jp=1,jml
3338	! If the grid has any litter, we should interpolate
3339	! this grid point.
3340	#ifdef __NAGFOR
3341	ltemp=IEEE_IS_NAN(littermap(ip,jp))
3342	#else
3343	ltemp=isnan(littermap(ip,jp))
3344	#endif
3345	IF( .NOT. ltemp)THEN
3346	IF ( littermap(ip,jp) .GT. min_sechiba .AND. littermap(ip,jp) .LT. fillvalue ) THEN
3347	mask(ip,jp) = 1
3348	IF (debug) THEN
3349	WRITE(numout,*) "update : littermap(",ip,jp,")) = ",littermap(ip,jp)
3350	ENDIF
3351	ENDIF
3352	ENDIF
3353	ENDDO
3354	ENDDO
3355	!
3356	!
3357	! The number of maximum vegetation map points in the GCM grid should
3358	! also be computed and not imposed here.
3359	!
3360	nbvmax = 200
3361	!
3362	callsign="Litter map"
3363	!
3364	ok_interpol = .FALSE.
3365	DO WHILE ( .NOT. ok_interpol )
3366	WRITE(numout,*) "Projection arrays for ",callsign," : "
3367	WRITE(numout,*) "nbvmax = ",nbvmax
3368	!
3369	ALLOC_ERR=-1
3370	ALLOCATE(sub_index(npts, nbvmax,2), STAT=ALLOC_ERR)
3371	IF (ALLOC_ERR/=0) THEN
3372	WRITE(numout,*) "ERROR IN ALLOCATION of sub_index : ",ALLOC_ERR
3373	CALL ipslerr_p (3,'sapiens_forestry_read_litter', 'error in sub_index','','')
3374	ENDIF
3375	sub_index(:,:,:)=0
3376
3377	ALLOC_ERR=-1
3378	ALLOCATE(sub_area(npts, nbvmax), STAT=ALLOC_ERR)
3379	IF (ALLOC_ERR/=0) THEN
3380	WRITE(numout,*) "ERROR IN ALLOCATION of sub_area : ",ALLOC_ERR
3381	CALL ipslerr_p (3,'sapiens_forestry_read_litter', 'error in sub_area','','')
3382	ENDIF
3383	sub_area(:,:)=zero
3384	!
3385	CALL aggregate_p(npts, lalo, neighbours, resolution, contfrac, &
3386	& iml, jml, lon_ful, lat_ful, mask, callsign, &
3387	& nbvmax, sub_index, sub_area, ok_interpol)
3388	!
3389	IF ( .NOT. ok_interpol ) THEN
3390	DEALLOCATE(sub_area)
3391	DEALLOCATE(sub_index)
3392	ENDIF
3393	!
3394	nbvmax = nbvmax * 2
3395	ENDDO
3396	!
3397	DO ipts = 1, npts
3398
3399	litter_demand(ipts)=zero
3400	area_sum=zero
3401
3402	DO ibvm=1, nbvmax
3403	! Leave the do loop if all sub areas are treated, sub_area <= 0
3404	IF ( sub_area(ipts,ibvm) <= zero ) EXIT
3405	ip = sub_index(ipts,ibvm,1)
3406	jp = sub_index(ipts,ibvm,2)
3407	area_sum=area_sum+sub_area(ipts,ibvm)
3408	litter_demand(ipts)=litter_demand(ipts)+sub_area(ipts,ibvm)*littermap(ip,jp)
3409	ENDDO
3410
3411	IF(area_sum .LT. min_stomate)THEN
3412	WRITE(numout,*) 'Missing data for the litter map on this land point!'
3413	WRITE(numout,*) 'ipts,lalo',ipts,lalo(ipts,:)
3414	litter_demand(ipts)=zero
3415	ELSE
3416	litter_demand(ipts)=litter_demand(ipts)/area_sum
3417	! This is in gC/year/m**2 from the interpolation. I need to convert it back
3418	! to absolute numbers.
3419	litter_demand(ipts)=litter_demand(ipts)resolution(ipts,1)resolution(ipts,2)
3420	ENDIF
3421
3422	ENDDO
3423
3424	!
3425	DEALLOCATE(littermap)
3426	DEALLOCATE(temp_resolution)
3427	DEALLOCATE(lat_lu,lon_lu)
3428	DEALLOCATE(lat_ful,lon_ful)
3429	DEALLOCATE(mask)
3430	IF(ALLOCATED(sub_index)) DEALLOCATE(sub_index)
3431	IF(ALLOCATED(sub_area)) DEALLOCATE(sub_area)
3432
3433
3434	IF ( printlev >= 5 ) WRITE(numout,*) 'Leaving sapiens_forestry_read_litter'
3435
3436	END SUBROUTINE sapiens_forestry_read_litter
3437
3438	!! ================================================================================================================================
3439	!! SUBROUTINE : sapiens_forestry_litter_raking
3440	!!
3441	!>\BRIEF Redistribute litter between PFTs due to litter raking
3442	!!
3443	!! DESCRIPTION :
3444	!! (1) The transfer of forest litter to crop lands within each grid cell depends on the litter demands.
3445	!! This is detemined by carbon litter demand maps read in stomate.f90.
3446	!! (2) Determine the CN ratio of litter pools, because the CN
3447	!! ratio should be the same before and after the litter raking.
3448	!! (3) Determine the amount of litter removed from each PFT. We already know
3449	!! the total amount af litter to be transfered from the forest
3450	!! to the crops, but it has to be distibuted amongst the PTF. This is determinde
3451	!! by the fraction of above ground litter of a given PFT to the total above ground
3452	!! litter in the grid cell. (The litter fraction of a PFT can very well be
3453	!! different from veget_max).
3454	!! (4) Remove nitrogen from the forest litter pools, such that a
3455	!! constant CN ratio is kept before and after the litter raking.
3456	!! (5) Move the litter of C and N to the crop PFTs in the grid.
3457	!!
3458	!! Be obs on the units in this subroutine - the changes between
3459	!! mass per grid cell and mass per m2, because the litter
3460	!! demand maps gC per grid cell.
3461	!!
3462	!! RECENT CHANGE(S) : Nitrogen was added to the litter raking subroutine during
3463	!! summer 2019
3464	!!
3465	!! MAIN OUTPUT VARIABLE(S): ::litter pools are modified
3466	!!
3467	!! REFERENCE(S) :
3468	!!
3469	!! FLOWCHART :
3470	!! \n
3471	!_ ================================================================================================================================
3472
3473	SUBROUTINE sapiens_forestry_litter_raking ( npts, veget_max, resolution, litter_demand, litter, lrake_frac )
3474
3475	!! 0. Variable and parameter declaration
3476
3477	!! 0.1 Input variables
3478	INTEGER(i_std), INTENT(in) :: npts !! Domain size - number of pixels
3479	!! (dimensionless)
3480	REAL(r_std), DIMENSION(:,:), INTENT(in) :: veget_max !! "Maximal" coverage fraction of a PFT (LAI
3481	!! -> infinity) on ground
3482	REAL(r_std), DIMENSION(:,:), INTENT(in) :: resolution !! The size in m of each grid-box in X and Y
3483	REAL(r_std), DIMENSION(:), INTENT(in) :: litter_demand !! The litter removed from each pixel due
3484	!! to litter raking at the end of the year.
3485	!! @tex $(gC year^{-1})$ @endtex
3486
3487	!! 0.2 Output
3488	REAL(r_std), DIMENSION(:,:), INTENT(out) :: lrake_frac !! Relative amount of litter that raked (-)
3489
3490	!! 0.3 Modified fields
3491	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(inout) :: litter !! Metabolic and structural litter, above
3492	!! and below ground
3493	!! @tex $(gC m^{-2})$ @endtex
3494
3495	!! 0.4 Local variables
3496	INTEGER(i_std) :: ipts, ivm !! Indices
3497	REAL(r_std) :: total_litter !! Total amount of aboveground forest litter
3498	!! on the pixel (gC)
3499	REAL(r_std) :: ind_litter !! Resulting litter per PFT
3500	REAL(r_std) :: crop_frac !! Fraction of croplands on the pixel (-)
3501	REAL(r_std) :: forest_frac !! Fraction of forests on the pixel (-)
3502	REAL(r_std) :: litter_transfer !! Litter that we would like to move (demand)
3503	!! from the forests to the crops
3504	!! @tex $(gC year^{-1})$ @endtex
3505	REAL(r_std),DIMENSION(nlitt) :: litter_ratio !! Litter in a PFT as a fraction of all the
3506	!! forest litter in the pixel (-)
3507	REAL(r_std),DIMENSION(nvm) :: lfrac !! ??? how does this differ from litter_ratio
3508	REAL(r_std), DIMENSION(npts,nvm,nmbcomp,nelements) :: check_intern !! Contains the components of the internal
3509	!! mass balance chech for this routine
3510	!! @tex $(gC m^{-2} dt^{-1})$ @endtex
3511	REAL(r_std), DIMENSION(npts,nvm,nelements) :: closure_intern !! Check closure of internal mass balance
3512	!! @tex $(gC m^{-2} dt^{-1})$ @endtex
3513	REAL(r_std), DIMENSION(npts,nvm,nelements) :: pool_start !! Start pool of this routine
3514	!! @tex $(gC m^{-2} dt^{-1})$ @endtex
3515	REAL(r_std), DIMENSION(npts,nvm,nelements) :: pool_end !! End pool of this routine
3516	!! @tex $(gC m^{-2} dt^{-1})$ @endtex
3517	REAL(r_std) :: temp_sum !! Temporary storage
3518	INTEGER(i_std) :: iele, ilitt !! Indices
3519	INTEGER(i_std) :: ilevs, imbc !! Indices
3520	REAL(r_std), DIMENSION(npts,nvm) :: pft_area !! The absolute area covered by this PFT
3521	!! @tex $(m^2)$ @endtex
3522	REAL(r_std) :: litter_lost !! Litter removed from the forest PFTs
3523	REAL(r_std) :: litter_gained !! Litter added to the crop PFTs
3524	LOGICAL :: lall_litter !! Flag to catch numerical issues
3525	REAL(r_std), DIMENSION(npts,nvm) :: veget_max_begin !! temporary storage of veget_max to check area conservation
3526	REAL(r_std), DIMENSION(nlitt,nvm) :: CN_litter !! CN ratios for the above ground litter pools
3527	!! before litter raking
3528	REAL(r_std), DIMENSION(nlitt,nvm,nelements) :: litter_lost_pft !! The litter lost due to litter raking for each PFT (gC or gN per grid cell)
3529	REAL(r_std), DIMENSION(npts) :: litter_transfer_N !! The total tranfers of N from the forested PFT to the crops (gN per grid cell).
3530	REAL(r_std) :: ind_litter_N !! Resulting N litter per PFT (gN per m2)
3531	!_ ================================================================================================================================
3532
3533	IF (firstcall_sapiens_forestry) THEN
3534	! Initialize local printlev if it is not already done
3535	printlev_loc=get_printlev('sapiens_forestry')
3536
3537	firstcall_sapiens_forestry=.FALSE.
3538	END IF
3539
3540	IF ( printlev_loc >= 3 ) WRITE(numout,*) 'Entering sapiens_forestry_litter_raking'
3541
3542	!! 1. Initialize biomass at first call
3543
3544	!! 1.2 Initialize check for mass balance closure
3545
3546	! Initial values of litter
3547	IF (err_act.GT.1) THEN
3548
3549	check_intern(:,:,:,:) = zero
3550	pool_start(:,:,:) = zero
3551
3552	DO iele = 1,nelements
3553
3554	! Litter pool (gC m-2) * (m2 m-2)
3555	DO ilitt = 1,nlitt
3556	DO ilevs = 1,nlevs
3557	pool_start(:,:,iele) = pool_start(:,:,iele) + &
3558	litter(:,ilitt,:,ilevs,iele) * veget_max(:,:)
3559	ENDDO
3560	ENDDO
3561
3562	ENDDO ! # nelements
3563
3564	!! 1.3 Initialize check for area conservation
3565	veget_max_begin(:,:) = veget_max(:,:)
3566
3567	ENDIF ! err_act.GT.1
3568
3569	!! 1.4 Initialize
3570	lrake_frac(:,:) = zero
3571	litter_transfer_N(:) = zero
3572
3573	!! Surface areas in m2
3574	DO ipts=1,npts
3575	pft_area(ipts,:)=area(ipts)* veget_max(ipts,:)
3576	ENDDO
3577
3578	! We need to find how much litter we can get from the forests for each
3579	! grid point.
3580	DO ipts=1,npts
3581
3582	litter_lost=zero
3583	litter_gained=zero
3584	total_litter=zero
3585	CN_litter(:,:) = zero
3586	litter_lost_pft(:,:,:) = zero
3587
3588
3589	DO ivm=1,nvm
3590	IF(is_tree(ivm))THEN
3591	! We will take all the above ground litter, since that is what
3592	! a farmer could pick up. Unlikely they dug into the soil to
3593	! collect litter. Notice these are in absolute units, not per
3594	! square meter.
3595	total_litter=total_litter+SUM(litter(ipts,:,ivm,iabove,icarbon))&
3596	*pft_area(ipts,ivm)
3597	ENDIF
3598	ENDDO
3599
3600	! We need to determine the CN ratio of the above grpund litter pools
3601	! such that we keep the CN ratios constant before and after the litter
3602	! raking
3603	DO ivm = 1,nvm
3604	IF(is_tree(ivm)) THEN
3605
3606	IF(SUM(litter(ipts,:,ivm,iabove,icarbon)) .GT. min_stomate &
3607	.AND. SUM(litter(ipts,:,ivm,iabove,initrogen)) .GT. min_stomate)THEN
3608
3609	! In case of empty litter N by component
3610	WHERE(litter(ipts,:,ivm,iabove,initrogen) .GT. zero)
3611	CN_litter(:,ivm) =litter(ipts,:,ivm,iabove,icarbon)/ &
3612	&litter(ipts,:,ivm,iabove,initrogen)
3613	ENDWHERE
3614
3615	IF(ivm .EQ. test_pft .AND. printlev_loc.GE.4) THEN
3616	WRITE(numout,*)'ivm',ivm
3617	WRITE(numout,*)'CN_litter before',CN_litter(:,ivm)
3618	ENDIF
3619
3620	ELSEIF(SUM(litter(ipts,:,ivm,iabove,icarbon)) .LT. min_stomate &
3621	.AND. SUM(litter(ipts,:,ivm,iabove,initrogen)) .LT.min_stomate)THEN
3622
3623	! If there us no litter, we cannot take the CN ratio
3624
3625	ELSE
3626
3627	! There is something wrong with the CN ratio of the litter pools
3628	WRITE(numout,*)'ivm',ivm
3629	WRITE(numout,*)'litter C',litter(ipts,:,ivm,iabove,icarbon)
3630	WRITE(numout,*)'litter N',litter(ipts,:,ivm,iabove,initrogen)
3631	CALL ipslerr_p (3,'ERROR: in litter raking','Something is wrong with the CN ratio of litter &
3632	pools', 'One pools is empty, while the other is not','')
3633	ENDIF
3634	ENDIF
3635	ENDDO
3636
3637	! If we have no forest litter, we cannot rake anything.
3638	IF(total_litter .LT. min_stomate)CYCLE
3639
3640	! The forests are tricky. The litter they produce will not
3641	! necessarily be proportional to the area they take up.
3642	! We want to reduce the amount of litter in each forest in a
3643	! way that doesn't leave any litter pools negative.
3644	lfrac(:)=zero
3645	DO ivm=1,nvm
3646	IF(is_tree(ivm))THEN
3647	lfrac(ivm)=SUM(litter(ipts,:,ivm,iabove,icarbon))&
3648	*pft_area(ipts,ivm)/total_litter
3649	ENDIF
3650	ENDDO
3651
3652	! Now compute the amount of the litter that we need to take to
3653	! satisfy our demand.
3654	litter_transfer=litter_demand(ipts)
3655
3656	! We don't take more litter than we have. Putting in
3657	! a logical flag to prevent us from having numerical
3658	! issues.
3659	lall_litter=.FALSE.
3660	IF(litter_transfer .GT. total_litter) THEN
3661	lall_litter=.TRUE.
3662	litter_transfer=total_litter
3663	ENDIF
3664
3665	! How many crop PFTs are present on this grid square?
3666	crop_frac=zero
3667	DO ivm=1,nvm
3668	IF(.NOT. natural(ivm))THEN
3669	crop_frac=crop_frac+veget_max(ipts,ivm)
3670	ENDIF
3671	ENDDO
3672
3673	! Litter raking is only conducted if we have crops in the grid cell.
3674	IF(crop_frac .GT. min_stomate)THEN
3675	DO ivm=1,nvm
3676
3677	IF(veget_max(ipts,ivm) .LE. min_stomate) CYCLE
3678
3679	! If we have a forest, we remove litter. Keep the ratio of
3680	! the litter pools the same.
3681	IF(is_tree(ivm))THEN
3682	ind_litter=SUM(litter(ipts,:,ivm,iabove,icarbon))*&
3683	pft_area(ipts,ivm)
3684
3685	! I guess there is a chance that a forest has no litter.
3686	! This means there is nothing to take away.
3687	IF(ind_litter .LE. min_stomate)CYCLE
3688
3689	litter_ratio(:)=litter(ipts,:,ivm,iabove,icarbon)*&
3690	pft_area(ipts,ivm)/ind_litter
3691	lrake_frac(ipts,ivm)=lfrac(ivm)*litter_transfer/ind_litter
3692
3693	! We don't want to have negative litter values
3694	IF(lall_litter)THEN
3695	ind_litter=zero
3696	ELSE
3697	ind_litter=ind_litter-lfrac(ivm)*litter_transfer
3698	ENDIF
3699	litter_lost=litter_lost-lfrac(ivm)*litter_transfer
3700
3701	! We need to estimate the lost litter of C for each PFT to
3702	! calculate lost N as well. But we also need to estimate, to
3703	! total N lost. The below is total estimates per grid cell
3704
3705	litter_lost_pft(:,ivm,icarbon)=-lfrac(ivm)litter_transfer&
3706	litter_ratio(:)
3707	WHERE(CN_litter(:,ivm) .GT. zero)
3708	litter_lost_pft(:,ivm,initrogen)=-lfrac(ivm)litter_transfer&
3709	litter_ratio(:)/CN_litter(:,ivm)
3710	ENDWHERE
3711	litter_transfer_N(ipts)=litter_transfer_N(ipts)-&
3712	SUM(litter_lost_pft(:,ivm,initrogen))
3713
3714	! This is a problem. We are taking differences of numbers
3715	! which are 1e9 or so, which means we could have a negative
3716	! value that is insignificant. I put the flag to try to
3717	! prevent that.
3718	IF(.NOT. lall_litter .AND. ind_litter .LT. -min_stomate)THEN
3719	WRITE(numout,*) "Took away too much litter!"
3720	WRITE(numout,*) "ipts,ivm:",ipts,ivm
3721	WRITE(numout,*) "lfrac(ivm),ind_litter,litter_transfer: ",&
3722	lfrac(ivm),SUM(litter(ipts,:,ivm,iabove,icarbon))*&
3723	pft_area(ipts,ivm),litter_transfer
3724	ENDIF
3725	litter(ipts,:,ivm,iabove,icarbon)=litter_ratio(:)*ind_litter&
3726	/pft_area(ipts,ivm)
3727	WHERE(CN_litter(:,ivm) .GT. zero)
3728	litter(ipts,:,ivm,iabove,initrogen)=litter_ratio(:)*ind_litter&
3729	/pft_area(ipts,ivm)/CN_litter(:,ivm)
3730	ENDWHERE
3731	! The resulting litter pool can also be determined be
3732	! subtraction now that the litter lost per pft has been
3733	! calculated (litter_lost_pft) to keep constant CN ratio. It
3734	! might be more intuitive when people are looking into the code.
3735	! It has been verifyed that both methods gives the same results.
3736	! litter(ipts,:,ivm,iabove,initrogen)=litter(ipts,:,ivm,iabove,initrogen)+&
3737	! litter_lost_pft(:,ivm,initrogen)/pft_area(ipts,ivm)
3738
3739	!--- DEBUG ---!
3740	IF (printlev_loc>=4) THEN
3741	IF(ipts .EQ. test_grid .and.ivm .EQ. test_pft) THEN
3742	WRITE(numout,*)'ivm',ivm
3743	WRITE(numout,*)'litter',litter(ipts,:,ivm,iabove,initrogen)
3744	WRITE(numout,*)'CNratio after',litter(ipts,:,ivm,iabove,icarbon)/&
3745	litter(ipts,:,ivm,iabove,initrogen)
3746	WRITE(numout,*)'litter_lost_pft carbon',litter_lost_pft(:,ivm,icarbon)
3747	WRITE(numout,*)'litter_lost_pft nitrogen',litter_lost_pft(:,ivm,initrogen)
3748	WRITE(numout,*)'litter_transfer', litter_transfer
3749	WRITE(numout,*)'lfrac', lfrac(ivm)
3750	WRITE(numout,*)'litter_ratio',litter_ratio
3751	ENDIF
3752	ENDIF
3753
3754	ELSEIF( .NOT. natural(ivm))THEN
3755
3756	! If we have a crop, we add litter.
3757	ind_litter=SUM(litter(ipts,:,ivm,iabove,icarbon))* &
3758	pft_area(ipts,ivm)
3759	ind_litter_N=SUM(litter(ipts,:,ivm,iabove,initrogen))* &
3760	pft_area(ipts,ivm)
3761
3762	! Is it possible that we have no crop litter but we
3763	! have veget_max? Maybe.
3764	IF(ind_litter .LT. min_stomate)THEN
3765	litter_ratio(:)=un/REAL(nlitt)
3766	ELSE
3767	litter_ratio(:)=litter(ipts,:,ivm,iabove,icarbon)*&
3768	pft_area(ipts,ivm)/ind_litter
3769	ENDIF
3770	! This is the total litter we have for this PFT in
3771	! this pixel
3772	ind_litter=ind_litter+veget_max(ipts,ivm)/crop_frac*&
3773	litter_transfer
3774	ind_litter_N=ind_litter_N+veget_max(ipts,ivm)/crop_frac*&
3775	litter_transfer_N(ipts)
3776	! could add a dimension to ind_litter, in stead of making a new
3777	! variable.
3778	litter_gained=litter_gained+veget_max(ipts,ivm)/&
3779	crop_frac*litter_transfer
3780	litter(ipts,:,ivm,iabove,icarbon)=litter_ratio(:)*&
3781	ind_litter/pft_area(ipts,ivm)
3782	litter(ipts,:,ivm,iabove,initrogen)=litter_ratio(:)*&
3783	ind_litter_N/pft_area(ipts,ivm)
3784
3785
3786	ENDIF ! checking to see if this PFT is a crop or forest
3787
3788	ENDDO ! loop over PFTs
3789
3790	ELSE
3791
3792	! We have no agricultural PFTs to move the litter to, so we
3793	! assume that no litter raking was done in this grid square.
3794
3795	ENDIF ! checking to see if we have crops
3796
3797	ENDDO ! loop over pixels
3798
3799	!! 2. Check numerical consistency of this routine
3800
3801	IF (err_act.GT.1) THEN
3802
3803	! 2.2 Check surface area
3804	CALL check_vegetation_area("stomate_prescribe", npts, veget_max_begin, &
3805	veget_max,'pft')
3806
3807	! 2.3 Mass balance closure
3808	! 2.3.1 Calculate final biomass
3809	pool_end = zero
3810	DO iele = 1,nelements
3811	DO ilitt = 1,nlitt
3812	DO ilevs = 1,nlevs
3813	pool_end(:,:,iele) = pool_end(:,:,iele) + &
3814	litter(:,ilitt,:,ilevs,iele) * veget_max(:,:)
3815	ENDDO
3816	ENDDO
3817	ENDDO
3818
3819	!! 2.3.2 Calculate mass balance
3820	! Common processes
3821	DO iele=1,nelements
3822	check_intern(:,:,ipoolchange,iele) = -un * &
3823	(pool_end(:,:,iele) - pool_start(:,:,iele))
3824	ENDDO
3825
3826	closure_intern = zero
3827	DO imbc = 1,nmbcomp
3828	DO iele=1,nelements
3829	! Debug
3830	IF (printlev_loc>=4) WRITE(numout,*) &
3831	'check_intern, ivm, imbc, iele, ', imbc, &
3832	iele, check_intern(:,test_pft,imbc,iele)
3833	!-
3834	closure_intern(:,:,iele) = closure_intern(:,:,iele) + &
3835	check_intern(:,:,imbc,iele)
3836	ENDDO
3837	ENDDO
3838
3839	! 4.3.3 Check mass balance closure
3840	CALL check_mass_balance("sapiens_forestry_litter_raking", closure_intern, npts, &
3841	pool_end, pool_start, veget_max, 'pixel')
3842
3843	ENDIF ! err_act.GT.1
3844
3845	IF ( printlev >= 4 ) WRITE(numout,*) 'Leaving sapiens_forestry_litter_raking'
3846
3847	END SUBROUTINE sapiens_forestry_litter_raking
3848
3849
3850	!! ================================================================================================================================
3851	!! SUBROUTINE : sapiens_forestry_read_species_change
3852	!!
3853	!>\BRIEF Read in a map that gives the PFT that needs to be planted
3854	!! when the current s PFT is harvested.
3855	!!
3856	!! DESCRIPTION :
3857	!!
3858	!! NOTE: This routine was mostly copied from slowproc where the PFTmap is read in.
3859	!! Grid interpolation is used, but only to look at the nearby pixels to see
3860	!! see which management strategy is dominant.
3861	!!
3862	!! RECENT CHANGE(S) : None
3863	!!
3864	!! MAIN OUTPUT VARIABLE(S): ::species_map
3865	!!
3866	!! REFERENCE(S) :
3867	!!
3868	!! FLOWCHART :
3869	!! \n
3870	!_ ================================================================================================================================
3871
3872	SUBROUTINE sapiens_forestry_read_species_change ( npts, lalo, neighbours, resolution, contfrac, species_map )
3873
3874	!! 0. Variable and parameter declaration
3875
3876	!! 0.1 Input variables
3877	INTEGER(i_std), INTENT(in) :: npts !! Domain size - number of pixels
3878	!! (dimensionless)
3879	REAL(r_std), DIMENSION(:,:), INTENT(in) :: lalo !! Vector of latitude and longitudes (beware of the order !)
3880	INTEGER(i_std), DIMENSION(:,:), INTENT(in) :: neighbours !!
3881	REAL(r_std), DIMENSION(:,:), INTENT(in) :: resolution !! The size in m of each grid-box in X and Y
3882	REAL(r_std), DIMENSION(:), INTENT(in) :: contfrac !! Fraction of continent in the grid
3883
3884	!! 0.2 Output
3885	INTEGER(i_std), DIMENSION(:,:), INTENT(out) :: species_map !! The number of a PFT that each PFT will
3886	!! be converted into after death.
3887
3888	!! 0.3 Modified fields
3889
3890
3891	!! 0.4 Local variables
3892	CHARACTER(LEN=80) :: filename !! A string to hold the file name
3893	LOGICAL :: debug=.FALSE. !! A flag to print out debugging information.
3894	INTEGER(i_std) :: fid !! The ID of the NetCDF file.
3895	INTEGER(i_std) :: nb_coord !! The number of coordinates in the NetCDF file
3896	INTEGER(i_std) :: nb_gat !!
3897	INTEGER(i_std) :: nb_var !! The number of variables in the NetCDF file
3898	INTEGER(i_std) :: nb_dim !! The number of dimensions in the NetCDF file
3899	INTEGER(i_std) :: iml, jml, lml,ivm !! indices
3900	INTEGER(i_std) :: ip, inbv, jp !! indices
3901	LOGICAL :: l_ex !! A flag which indicates if a variable
3902	!! exists in the NetCDF file
3903	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_lu, lon_lu !! The latitude and longitude read in from
3904	!! the NetCDF file
3905	INTEGER,DIMENSION(flio_max_var_dims) :: l_d_w !! List of the dimension lengths of the variable
3906	!! in the NetCDF file
3907	INTEGER(i_std) :: ipts !! index
3908	INTEGER(i_std) :: ALLOC_ERR !! A flag tripped if we have an error in allocation
3909	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: scmap_r !! The map read in from the NetCDF file
3910	INTEGER(i_std),DIMENSION(:,:,:),ALLOCATABLE :: scmap_i !! The integer form of the map read in
3911	INTEGER(i_std) :: closest_lat !! The index of the closest latitude we found.
3912	INTEGER(i_std) :: closest_lon !! The index of the closest longitude we found.
3913	REAL(r_std) :: distance !! The distance from the current point to the
3914	!! point on the map
3915	REAL(r_std) :: closest_dist !! The distance to the closet point we've found.
3916	INTEGER :: large_int !! A number which indicates that the grid data
3917	!! is not available
3918	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lat_ful, lon_ful
3919	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: mask
3920	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area
3921	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: sub_index
3922	CHARACTER(LEN=30) :: callsign
3923	INTEGER(i_std) :: ibvm, nbvmax,ifm, ivma
3924	LOGICAL :: ok_interpol ! optionnal return of aggregate_2d
3925	REAL(r_std) :: sum_sc,sumf
3926	REAL(r_std),DIMENSION(nvmap) :: sc_sum
3927	LOGICAL :: ltemp ! temporary logical variable
3928	!_ ================================================================================================================================
3929
3930	!IF ( printlev >= 4 )
3931	WRITE(numout,*) 'Entering sapiens_forestry_read_species_change'
3932
3933	! This integer has to be large enough that it never shows up on the map, without being
3934	! so large that is causes overflows. This needs to be larger than the largest possible
3935	! value multiplied by the number of PFTs. Since each value should never be larger than
3936	! the number of PFTs, this should suffice.
3937	large_int = nvmap*nvmap+1
3938
3939	!
3940	!Config Key = SPECIES_CHANGE_FILE
3941	!Config Desc = Name of file from which the species change map is to be read
3942	!Config If = OK_STOMATE
3943	!Config Def = replant_species.nc
3944	!Config Help = The name of the file to be opened to read the species to replant
3945	!Config for each pixel and PFT.
3946	!Config Units = [FILE]
3947	!
3948	filename = 'replant_species.nc'
3949	CALL getin_p('SPECIES_CHANGE_FILE',filename)
3950
3951	IF (is_root_prc) THEN
3952	IF (debug) THEN
3953	WRITE(numout,*) "Entering sapiens_forestry_read_species_change. Debug mode."
3954	WRITE (*,'(/," --> fliodmpf")')
3955	CALL fliodmpf (TRIM(filename))
3956	WRITE (*,'(/," --> flioopfd")')
3957	ENDIF
3958	CALL flioopfd (TRIM(filename),fid,nb_dim=nb_coord,nb_var=nb_var,nb_gat=nb_gat)
3959	IF (debug) THEN
3960	WRITE (*,'(" Number of coordinate in the file : ",I2)') nb_coord
3961	WRITE (*,'(" Number of variables in the file : ",I2)') nb_var
3962	WRITE (*,'(" Number of global attributes in the file : ",I2)') nb_gat
3963	ENDIF
3964	ENDIF
3965	CALL bcast(nb_coord)
3966	CALL bcast(nb_var)
3967	CALL bcast(nb_gat)
3968
3969	! This finds the number of longitude points in the file.
3970	IF (is_root_prc) &
3971	CALL flioinqv (fid,v_n="lon",l_ex=l_ex,nb_dims=nb_dim,len_dims=l_d_w)
3972	CALL bcast(l_d_w)
3973	iml=l_d_w(1)
3974	WRITE(numout,*) "Species change Map: iml =",iml
3975
3976	! This finds the number of latitude points in the file.
3977	IF (is_root_prc) &
3978	CALL flioinqv (fid,v_n="lat",l_ex=l_ex,nb_dims=nb_dim,len_dims=l_d_w)
3979	CALL bcast(l_d_w)
3980	jml=l_d_w(1)
3981	WRITE(numout,*) "Species change Map: jml =",jml
3982
3983	! Now find the number of PFTs in the file. If this is not equal to the number
3984	! of PFTs that we actually have, that's a problem and we'll crash.
3985	IF (is_root_prc) &
3986	CALL flioinqv (fid,v_n="NEWSPECIES",l_ex=l_ex,nb_dims=nb_dim,len_dims=l_d_w)
3987	CALL bcast(l_d_w)
3988	lml=l_d_w(3)
3989
3990	IF (lml /= nvmap) THEN
3991	WRITE(numout,*) 'lml = ',lml
3992	WRITE(numout,*) 'nvmap = ',nvmap
3993	WRITE(numout,*) 'Stopping. '
3994	CALL ipslerr_p (3,'sapians_forestry', &
3995	& 'Problem with species change map.','lml /= nvmap', &
3996	& '(number of pft must be equal)')
3997	ENDIF
3998	!
3999
4000	! Allocate the map that will be read in
4001	WRITE(numout,*) 'Reading the species change map'
4002	!
4003	ALLOC_ERR=-1
4004	ALLOCATE(scmap_r(iml,jml,nvmap), STAT=ALLOC_ERR)
4005	IF (ALLOC_ERR/=0) THEN
4006	WRITE(numout,*) "ERROR IN ALLOCATION of scmap_r : ",ALLOC_ERR
4007	CALL ipslerr_p (3,'sapiens_forestry_read_species_change', 'error in scmap_r','','')
4008	ENDIF
4009	ALLOC_ERR=-1
4010	ALLOCATE(scmap_i(iml,jml,nvmap), STAT=ALLOC_ERR)
4011	IF (ALLOC_ERR/=0) THEN
4012	WRITE(numout,*) "ERROR IN ALLOCATION of scmap_i : ",ALLOC_ERR
4013	CALL ipslerr_p (3,'sapiens_forestry_read_species_change', 'error in scmap_i','','')
4014	ENDIF
4015
4016	! This reads in the map that is in the file
4017	IF (is_root_prc) THEN
4018	scmap_r(:,:,:)=large_int*2.0
4019	CALL fliogetv (fid,"NEWSPECIES", scmap_r, start=(/ 1, 1, 1 /), count=(/ iml, jml, nvmap /))
4020	! Right now the values are all real, but they should be integers.
4021	! Careful, NINT might not work if the precision of scmap_r is not single.
4022	! In that case, IDNINT should be used.
4023	DO ip=1,iml
4024	DO jp=1,jml
4025	DO ivma=1,nvmap
4026	! There will be some fill values in here. If we pass
4027	! a large fill value to NINT, it crashes. So let's
4028	! test for it
4029	#ifdef __NAGFOR
4030	ltemp=IEEE_IS_NAN(scmap_r(ip,jp,ivma))
4031	#else
4032	ltemp=isnan(scmap_r(ip,jp,ivma))
4033	#endif
4034	IF(ltemp)THEN
4035	scmap_i(ip,jp,ivma)=large_int
4036
4037	ELSE
4038	IF(scmap_r(ip,jp,ivma) .GE. 0.0 .AND. scmap_r(ip,jp,ivma) < large_int)THEN
4039	scmap_i(ip,jp,ivma)=NINT(scmap_r(ip,jp,ivma))
4040	ELSE
4041	! This value should be big enough that we don't barl ourselves
4042	! below.
4043	scmap_i(ip,jp,ivma)=large_int
4044	ENDIF
4045	ENDIF
4046	ENDDO
4047	ENDDO
4048	ENDDO
4049	ENDIF
4050
4051	CALL bcast(scmap_i)
4052
4053	! Now I need to get the latitude and longitude
4054	! First, get the axes from the map file.
4055	ALLOC_ERR=-1
4056	ALLOCATE(lat_lu(jml), STAT=ALLOC_ERR)
4057	IF (ALLOC_ERR/=0) THEN
4058	WRITE(numout,*) "ERROR IN ALLOCATION of lat_lu : ",ALLOC_ERR
4059	CALL ipslerr_p (3,'sapiens_forestry_read_species_change', 'error in lat_lu','','')
4060	ENDIF
4061	ALLOC_ERR=-1
4062	ALLOCATE(lon_lu(iml), STAT=ALLOC_ERR)
4063	IF (ALLOC_ERR/=0) THEN
4064	WRITE(numout,*) "ERROR IN ALLOCATION of lon_lu : ",ALLOC_ERR
4065	CALL ipslerr_p (3,'sapiens_forestry_read_species_change', 'error in lon_lu','','')
4066	ENDIF
4067	IF (is_root_prc) THEN
4068	CALL fliogstc (fid, x_axis=lon_lu,y_axis=lat_lu)
4069	ENDIF
4070	CALL bcast(lon_lu)
4071	CALL bcast(lat_lu)
4072
4073	! Now I can interpolate. Remember that we are dealing with integer
4074	! values and PFTs. Therefore, we cannot do a strict
4075	! interpolation, since that might gives values of the PFT to be
4076	! 2.3, or something ridiculous. We cannot just round the number, either,
4077	! since we could have some squares with poplar (PFT 15) and oak (PFT 13),
4078	! rounding would give a completely different PFT (14). So we look to see what the
4079	! most prevalent type of PFT of nearby pixels is, and then just take that.
4080
4081	ALLOC_ERR=-1
4082	ALLOCATE(lat_ful(iml,jml), STAT=ALLOC_ERR)
4083	IF (ALLOC_ERR/=0) THEN
4084	WRITE(numout,*) "ERROR IN ALLOCATION of lat_ful : ",ALLOC_ERR
4085	CALL ipslerr_p (3,'sapiens_forestry_read_species_change', 'error in lat_ful','','')
4086	ENDIF
4087	ALLOC_ERR=-1
4088	ALLOCATE(lon_ful(iml,jml), STAT=ALLOC_ERR)
4089	IF (ALLOC_ERR/=0) THEN
4090	WRITE(numout,*) "ERROR IN ALLOCATION of lon_ful : ",ALLOC_ERR
4091	CALL ipslerr_p (3,'sapiens_forestry_read_species_change', 'error in lon_ful','','')
4092	ENDIF
4093	!
4094	DO ip=1,iml
4095	lon_ful(ip,:)=lon_lu(ip)
4096	ENDDO
4097	DO jp=1,jml
4098	lat_ful(:,jp)=lat_lu(jp)
4099	ENDDO
4100	!
4101	!
4102	WRITE(numout,*) 'Reading the SPECIES CHANGE file'
4103	!
4104	!
4105
4106	!
4107	! Mask of permitted variables.
4108	!
4109	ALLOC_ERR=-1
4110	ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR)
4111	IF (ALLOC_ERR/=0) THEN
4112	WRITE(numout,*) "ERROR IN ALLOCATION of mask : ",ALLOC_ERR
4113	CALL ipslerr_p (3,'sapiens_forestry_read_species_change', 'error in mask','','')
4114	ENDIF
4115	!
4116	mask(:,:) = 0
4117	DO ip=1,iml
4118	DO jp=1,jml
4119	! If at least one PFT has a management strategy, we should interpolate
4120	! this grid point.
4121	sum_sc=SUM(scmap_i(ip,jp,:))
4122	IF ( sum_sc .GT. min_sechiba .AND. sum_sc .LT. large_int) THEN
4123	mask(ip,jp) = 1
4124	IF (debug) THEN
4125	WRITE(numout,*) "update : SUM(scmap(",ip,jp,")) = ",sum_sc
4126	ENDIF
4127	ENDIF
4128	ENDDO
4129	ENDDO
4130	!
4131	!
4132	! The number of maximum vegetation map points in the GCM grid should
4133	! also be computed and not imposed here.
4134	!
4135	nbvmax = 200
4136	!
4137	callsign="Species change map"
4138	!
4139	ok_interpol = .FALSE.
4140	DO WHILE ( .NOT. ok_interpol )
4141	WRITE(numout,*) "Projection arrays for ",callsign," : "
4142	WRITE(numout,*) "nbvmax = ",nbvmax
4143	!
4144	ALLOC_ERR=-1
4145	ALLOCATE(sub_index(npts, nbvmax,2), STAT=ALLOC_ERR)
4146	IF (ALLOC_ERR/=0) THEN
4147	WRITE(numout,*) "ERROR IN ALLOCATION of sub_index : ",ALLOC_ERR
4148	CALL ipslerr_p (3,'sapiens_forestry_read_species_change', 'error in sub_index','','')
4149	ENDIF
4150	sub_index(:,:,:)=0
4151
4152	ALLOC_ERR=-1
4153	ALLOCATE(sub_area(npts, nbvmax), STAT=ALLOC_ERR)
4154	IF (ALLOC_ERR/=0) THEN
4155	WRITE(numout,*) "ERROR IN ALLOCATION of sub_area : ",ALLOC_ERR
4156	CALL ipslerr_p (3,'sapiens_forestry_read_species_change', 'error in sub_area','','')
4157	ENDIF
4158	sub_area(:,:)=zero
4159	!
4160	CALL aggregate_p(npts, lalo, neighbours, resolution, contfrac, &
4161	& iml, jml, lon_ful, lat_ful, mask, callsign, &
4162	& nbvmax, sub_index, sub_area, ok_interpol)
4163	!
4164	IF ( .NOT. ok_interpol ) THEN
4165	DEALLOCATE(sub_area)
4166	DEALLOCATE(sub_index)
4167	ENDIF
4168	!
4169	nbvmax = nbvmax * 2
4170	ENDDO
4171	!
4172	DO ipts = 1, npts
4173	! For this point, we need to see what dominant new species are
4174	! nearby.
4175	sumf=zero
4176	! We need to do this for every PFT
4177	DO ivma=1,nvmap
4178	! If it's not a forest, we don't care. We will
4179	! always replant it as the same PFT.
4180	IF(.NOT. is_tree(start_index(ivma)))THEN
4181	species_map(ipts,ivma)=ivma
4182	ELSE
4183	sc_sum(:)=0.0
4184	DO ibvm=1, nbvmax
4185	! Leave the do loop if all sub areas are treated, sub_area <= 0
4186	IF ( sub_area(ipts,ibvm) <= zero ) EXIT
4187	ip = sub_index(ipts,ibvm,1)
4188	jp = sub_index(ipts,ibvm,2)
4189	ifm=scmap_i(ip,jp,ivma)
4190	sc_sum(ifm)=sc_sum(ifm)+sub_area(ipts,ibvm)
4191	ENDDO
4192	! Whichever FM type has the most area, we use that one.
4193	species_map(ipts,ivma)=MAXLOC(sc_sum(:),1)
4194	ENDIF
4195	ENDDO
4196	ENDDO
4197	!
4198	DEALLOCATE(scmap_i)
4199	DEALLOCATE(scmap_r)
4200	DEALLOCATE(lat_lu,lon_lu)
4201	DEALLOCATE(lat_ful,lon_ful)
4202	DEALLOCATE(mask)
4203	IF(ALLOCATED(sub_index)) DEALLOCATE(sub_index)
4204	IF(ALLOCATED(sub_area)) DEALLOCATE(sub_area)
4205
4206
4207	IF ( printlev >= 5 ) WRITE(numout,*) 'Leaving sapiens_forestry_read_species_change'
4208
4209	END SUBROUTINE sapiens_forestry_read_species_change
4210
4211
4212	!! ================================================================================================================================
4213	!! SUBROUTINE : sapiens_forestry_read_desired_fm
4214	!!
4215	!>\BRIEF Read in a map that gives the desired forest management strategy
4216	!! after a new PFT was planted.
4217	!!
4218	!! DESCRIPTION :
4219	!!
4220	!!
4221	!! FM = 1 : No human intervention (ORCHIDEE default)
4222	!! 2 : Thinnings based on the RDI, clearcuts based on tree density,
4223	!! annual increment, and tree diameter. Thinnings from above and
4224	!! from below are determined by the sign of thstrat.
4225	!! 3 : Coppices
4226	!! 4 : Short rotation coppices
4227	!!
4228	!! NOTE: This routine was mostly copied from slowproc where the PFTmap is read in.
4229	!! Grid interpolation is used, but only to look at the nearby pixels to see
4230	!! see which management strategy is dominant.
4231	!!
4232	!! RECENT CHANGE(S) : None
4233	!!
4234	!! MAIN OUTPUT VARIABLE(S): ::forest_managed
4235	!!
4236	!! REFERENCE(S) :
4237	!!
4238	!! FLOWCHART :
4239	!! \n
4240	!_ ================================================================================================================================
4241
4242	SUBROUTINE sapiens_forestry_read_desired_fm ( npts, lalo, neighbours, resolution, &
4243	contfrac, desired_managed )
4244
4245	!! 0. Variable and parameter declaration
4246
4247	!! 0.1 Input variables
4248	INTEGER(i_std), INTENT(in) :: npts !! Domain size - number of pixels
4249	!! (dimensionless)
4250	REAL(r_std), DIMENSION(:,:), INTENT(in) :: lalo !! Vector of latitude and longitudes (beware of the order !)
4251	INTEGER(i_std), DIMENSION(:,:), INTENT(in) :: neighbours !!
4252	REAL(r_std), DIMENSION(:,:), INTENT(in) :: resolution !! The size in m of each grid-box in X and Y
4253	REAL(r_std), DIMENSION(:), INTENT(in) :: contfrac !! Fraction of continent in the grid
4254
4255	!! 0.2 Output
4256	INTEGER(i_std), DIMENSION(:,:), INTENT(out) :: desired_managed !! Forest management flag: 0 = orchidee
4257	!! standard, 1= self-thinning only, 2=
4258	!! high-stand, 3= high-stand smoothed, 4=
4259	!! coppices
4260
4261	!! 0.3 Modified fields
4262
4263
4264	!! 0.4 Local variables
4265	CHARACTER(LEN=80) :: filename !! A string to hold the file name
4266	LOGICAL :: debug=.FALSE. !! A flag to print out debugging information.
4267	INTEGER(i_std) :: fid !! The ID of the NetCDF file.
4268	INTEGER(i_std) :: nb_coord !! The number of coordinates in the NetCDF file
4269	INTEGER(i_std) :: nb_gat !!
4270	INTEGER(i_std) :: nb_var !! The number of variables in the NetCDF file
4271	INTEGER(i_std) :: nb_dim !! The number of dimensions in the NetCDF file
4272	INTEGER(i_std) :: iml,jml,lml,ivm !! indices
4273	INTEGER(i_std) :: ip, inbv, jp !! indices
4274	LOGICAL :: l_ex !! A flag which indicates if a variable
4275	!! exists in the NetCDF file
4276	REAL(r_std), ALLOCATABLE, DIMENSION(:) :: lat_lu, lon_lu !! The latitude and longitude read in from
4277	!! the NetCDF file
4278	INTEGER,DIMENSION(flio_max_var_dims) :: l_d_w !! List of the dimension lengths of the variable
4279	!! in the NetCDF file
4280	INTEGER(i_std) :: ipts !! index
4281	INTEGER(i_std) :: ALLOC_ERR !! A flag tripped if we have an error in allocation
4282	REAL(r_std),DIMENSION(:,:,:),ALLOCATABLE :: fmmap_r !! The map read in from the NetCDF file
4283	INTEGER(i_std),DIMENSION(:,:,:),ALLOCATABLE :: fmmap_i !! The integer form of the map read in
4284	INTEGER(i_std) :: closest_lat !! The index of the closest latitude we found.
4285	INTEGER(i_std) :: closest_lon !! The index of the closest longitude we found.
4286	REAL(r_std) :: distance !! The distance from the current point to the
4287	!! point on the map
4288	REAL(r_std) :: closest_dist !! The distance to the closet point we've found.
4289	INTEGER :: large_int !! A number which indicates that the grid data
4290	!! is not available
4291	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: lat_ful, lon_ful
4292	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:) :: mask
4293	REAL(r_std), ALLOCATABLE, DIMENSION(:,:) :: sub_area
4294	INTEGER(i_std), ALLOCATABLE, DIMENSION(:,:,:) :: sub_index
4295	CHARACTER(LEN=30) :: callsign
4296	INTEGER(i_std) :: ibvm, nbvmax,ifm, ivma
4297	LOGICAL :: ok_interpol !! optionnal return of aggregate_2d
4298	REAL(r_std) :: sum_fm,sumf
4299	REAL(r_std),DIMENSION(nfm_types) :: fm_sum
4300	LOGICAL :: ltemp !! temporary logical variable
4301	!_ ================================================================================================================================
4302
4303	IF ( printlev >= 4 ) WRITE(numout,*) 'Entering sapiens_forestry_read_desired_fm'
4304
4305	! This integer has to be large enough that it never shows up on the map, without being
4306	! so large that is causes overflows. Since none of the points on the map should be
4307	! larger than the number of FM strategies we have (nfm_types), this is sufficiently big.
4308	large_int = nvmap*nfm_types+1
4309
4310	!Config Key = FM_FILE
4311	!Config Desc = Name of file from which the forest management map is to be read
4312	!Config If = OK_STOMATE
4313	!Config Def = FMmap.nc
4314	!Config Help = The name of the file to be opened to read a forest management
4315	!Config map (including a layer for every PFT) is given here.
4316	!Config Units = [FILE]
4317	filename = 'FM_desired.nc'
4318	!CALL getin_p('FM_FILE',filename)
4319
4320	IF (is_root_prc) THEN
4321	IF (debug) THEN
4322	WRITE(numout,*) "Entering sapiens_forestry_read_desired_fm. Debug mode."
4323	WRITE (*,'(/," --> fliodmpf")')
4324	CALL fliodmpf (TRIM(filename))
4325	WRITE (*,'(/," --> flioopfd")')
4326	ENDIF
4327	CALL flioopfd (TRIM(filename),fid,nb_dim=nb_coord,nb_var=nb_var,nb_gat=nb_gat)
4328	IF (debug) THEN
4329	WRITE (*,'(" Number of coordinate in the file : ",I2)') nb_coord
4330	WRITE (*,'(" Number of variables in the file : ",I2)') nb_var
4331	WRITE (*,'(" Number of global attributes in the file : ",I2)') nb_gat
4332	ENDIF
4333	ENDIF
4334	CALL bcast(nb_coord)
4335	CALL bcast(nb_var)
4336	CALL bcast(nb_gat)
4337
4338	! This finds the number of longitude points in the file.
4339	IF (is_root_prc) &
4340	CALL flioinqv (fid,v_n="lon",l_ex=l_ex,nb_dims=nb_dim,len_dims=l_d_w)
4341	CALL bcast(l_d_w)
4342	iml=l_d_w(1)
4343	WRITE(numout,*) "FM desired Map: iml =",iml
4344
4345	! This finds the number of latitude points in the file.
4346	IF (is_root_prc) &
4347	CALL flioinqv (fid,v_n="lat",l_ex=l_ex,nb_dims=nb_dim,len_dims=l_d_w)
4348	CALL bcast(l_d_w)
4349	jml=l_d_w(1)
4350	WRITE(numout,*) "FM desired Map: jml =",jml
4351
4352	! Now find the number of PFTs in the file. If this is not equal to the number
4353	! of PFTs that we actually have, that's a problem and we'll crash.
4354	IF (is_root_prc) &
4355	CALL flioinqv (fid,v_n="FM_STRAT",l_ex=l_ex,nb_dims=nb_dim,len_dims=l_d_w)
4356	CALL bcast(l_d_w)
4357	lml=l_d_w(3)
4358
4359	IF (lml /= nvmap) THEN
4360	WRITE(numout,*) 'lml = ',lml
4361	WRITE(numout,*) 'nvmap = ',nvmap
4362	WRITE(numout,*) 'Stopping. '
4363	CALL ipslerr_p (3,'sapiens_forestry', &
4364	& 'Problem with the desired forest management strateg map.',&
4365	& 'lml /= nvmap', '(number of pft must be equal)')
4366	ENDIF
4367
4368	! Allocate the map that will be read in
4369	WRITE(numout,*) 'Reading the Desired Forest Management strategy file'
4370	ALLOC_ERR=-1
4371	ALLOCATE(fmmap_r(iml,jml,nvmap), STAT=ALLOC_ERR)
4372	IF (ALLOC_ERR/=0) THEN
4373	WRITE(numout,*) "ERROR IN ALLOCATION of fmmap_r : ",ALLOC_ERR
4374	CALL ipslerr_p (3,'sapiens_forestry_read_desired_fm', 'error in fmmap_r','','')
4375	ENDIF
4376	ALLOC_ERR=-1
4377	ALLOCATE(fmmap_i(iml,jml,nvmap), STAT=ALLOC_ERR)
4378	IF (ALLOC_ERR/=0) THEN
4379	WRITE(numout,*) "ERROR IN ALLOCATION of fmmap_i : ",ALLOC_ERR
4380	CALL ipslerr_p (3,'sapiens_forestry_read_desired_fm', 'error in fmmpap_i','','')
4381	ENDIF
4382
4383	! This reads in the map that is in the file
4384	IF (is_root_prc) THEN
4385	fmmap_r(:,:,:)=large_int*2.0
4386	CALL fliogetv (fid,"FM_STRAT", fmmap_r, start=(/ 1, 1, 1 /), &
4387	count=(/ iml, jml, nvmap /))
4388
4389	! Right now the values are all real, but they should be integers.
4390	! Careful, NINT might not work if the precision of fmmap_r is not
4391	! single. In that case, IDNINT should be used.
4392	DO ip=1,iml
4393	DO jp=1,jml
4394	DO ivma=1,nvmap
4395
4396	! There will be some fill values in here. If we pass
4397	! a large fill value to NINT, it crashes. So let's
4398	! test for it
4399	#ifdef __NAGFOR
4400	ltemp=IEEE_IS_NAN(fmmap_r(ip,jp,ivma))
4401	#else
4402	ltemp=isnan(fmmap_r(ip,jp,ivma))
4403	#endif
4404	IF(ltemp)THEN
4405	fmmap_i(ip,jp,ivma)=large_int
4406	ELSE
4407	IF(fmmap_r(ip,jp,ivma) .GE. 0.0 .AND. fmmap_r(ip,jp,ivma) < large_int) THEN
4408	fmmap_i(ip,jp,ivma) = NINT(fmmap_r(ip,jp,ivma))
4409	ELSE
4410
4411	! This value should be big enough that we don't barl ourselves
4412	! below.
4413	fmmap_i(ip,jp,ivma)=large_int
4414	ENDIF
4415	ENDIF
4416	ENDDO
4417	ENDDO
4418	ENDDO
4419	ENDIF
4420
4421	CALL bcast(fmmap_i)
4422
4423	! Now I need to get the latitude and longitude
4424	! First, get the axes from the map file.
4425	ALLOC_ERR=-1
4426	ALLOCATE(lat_lu(jml), STAT=ALLOC_ERR)
4427	IF (ALLOC_ERR/=0) THEN
4428	WRITE(numout,*) "ERROR IN ALLOCATION of lat_lu : ",ALLOC_ERR
4429	CALL ipslerr_p (3,'sapiens_forestry_read_desired_fm', 'error in lat_lu','','')
4430	ENDIF
4431	ALLOC_ERR=-1
4432	ALLOCATE(lon_lu(iml), STAT=ALLOC_ERR)
4433	IF (ALLOC_ERR/=0) THEN
4434	WRITE(numout,*) "ERROR IN ALLOCATION of lon_lu : ",ALLOC_ERR
4435	CALL ipslerr_p (3,'sapiens_forestry_read_desired_fm', 'error in lon_lu','','')
4436	ENDIF
4437	IF (is_root_prc) THEN
4438	CALL fliogstc (fid, x_axis=lon_lu,y_axis=lat_lu)
4439	ENDIF
4440	CALL bcast(lon_lu)
4441	CALL bcast(lat_lu)
4442
4443	! Now I can interpolate. Remember that we are dealing with integer
4444	! values and management strategies. Therefore, we cannot do a strict
4445	! interpolation, since that might gives values of the strategy to be
4446	! 2.3, or something ridiculous. We cannot just round the number, either,
4447	! since we could have some squares with 1 (no management) and some with
4448	! 3 (coppices); An average of that would be 2 (high stands), which doesn't
4449	! make any sense. So we look to see what the most prevalent type of
4450	! management of nearby pixels is, and then just take that.
4451	ALLOC_ERR=-1
4452	ALLOCATE(lat_ful(iml,jml), STAT=ALLOC_ERR)
4453	IF (ALLOC_ERR/=0) THEN
4454	WRITE(numout,*) "ERROR IN ALLOCATION of lat_ful : ",ALLOC_ERR
4455	CALL ipslerr_p (3,'sapiens_forestry_read_desired_fm', 'error in lat_ful','','')
4456	ENDIF
4457	ALLOC_ERR=-1
4458	ALLOCATE(lon_ful(iml,jml), STAT=ALLOC_ERR)
4459	IF (ALLOC_ERR/=0) THEN
4460	WRITE(numout,*) "ERROR IN ALLOCATION of lon_ful : ",ALLOC_ERR
4461	CALL ipslerr_p (3,'sapiens_forestry_read_desired_fm', 'error in lon_ful','','')
4462	ENDIF
4463
4464	DO ip=1,iml
4465	lon_ful(ip,:)=lon_lu(ip)
4466	ENDDO
4467	DO jp=1,jml
4468	lat_ful(:,jp)=lat_lu(jp)
4469	ENDDO
4470
4471	! Mask of permitted variables.
4472	ALLOC_ERR=-1
4473	ALLOCATE(mask(iml,jml), STAT=ALLOC_ERR)
4474	IF (ALLOC_ERR/=0) THEN
4475	WRITE(numout,*) "ERROR IN ALLOCATION of mask : ",ALLOC_ERR
4476	CALL ipslerr_p (3,'sapiens_forestry_read_desired_fm', 'error in mask','','')
4477	ENDIF
4478
4479	mask(:,:) = 0
4480	DO ip=1,iml
4481	DO jp=1,jml
4482	! If at least one PFT has a management strategy, we should interpolate
4483	! this grid point.
4484	sum_fm=SUM(fmmap_i(ip,jp,:))
4485	IF ( sum_fm .GT. min_sechiba .AND. sum_fm .LT. large_int) THEN
4486	mask(ip,jp) = 1
4487	IF (debug) THEN
4488	WRITE(numout,*) "update : SUM(fmmap(",ip,jp,")) = ",sum_fm
4489	ENDIF
4490	ENDIF
4491	ENDDO
4492	ENDDO
4493
4494
4495	! The number of maximum vegetation map points in the GCM grid should
4496	! also be computed and not imposed here.
4497	nbvmax = 200
4498
4499	callsign="Desired Forest Management map"
4500	ok_interpol = .FALSE.
4501	DO WHILE ( .NOT. ok_interpol )
4502	WRITE(numout,*) "Projection arrays for ",callsign," : "
4503	WRITE(numout,*) "nbvmax = ",nbvmax
4504
4505	ALLOC_ERR=-1
4506	ALLOCATE(sub_index(npts, nbvmax,2), STAT=ALLOC_ERR)
4507	IF (ALLOC_ERR/=0) THEN
4508	WRITE(numout,*) "ERROR IN ALLOCATION of sub_index : ",ALLOC_ERR
4509	CALL ipslerr_p (3,'sapiens_forestry_read_desired_fm', 'error in sub_index','','')
4510	ENDIF
4511	sub_index(:,:,:)=0
4512
4513	ALLOC_ERR=-1
4514	ALLOCATE(sub_area(npts, nbvmax), STAT=ALLOC_ERR)
4515	IF (ALLOC_ERR/=0) THEN
4516	WRITE(numout,*) "ERROR IN ALLOCATION of sub_area : ",ALLOC_ERR
4517	CALL ipslerr_p (3,'sapiens_forestry_read_desired_fm', 'error in sub_area','','')
4518	ENDIF
4519	sub_area(:,:)=zero
4520
4521	CALL aggregate_p(npts, lalo, neighbours, resolution, contfrac, &
4522	& iml, jml, lon_ful, lat_ful, mask, callsign, &
4523	& nbvmax, sub_index, sub_area, ok_interpol)
4524
4525	IF ( .NOT. ok_interpol ) THEN
4526	DEALLOCATE(sub_area)
4527	DEALLOCATE(sub_index)
4528	ENDIF
4529
4530	nbvmax = nbvmax * 2
4531	ENDDO
4532
4533	DO ipts = 1, npts
4534
4535	! For this point, we need to see what dominant FM types are
4536	! nearby.
4537	sumf=zero
4538
4539	! We need to do this for every PFT
4540	DO ivma=1,nvmap
4541
4542	IF(.NOT. is_tree(start_index(ivma)))THEN
4543
4544	! If it's not a forest, we don't care.
4545	desired_managed(ipts,ivma)=0
4546
4547	ELSE
4548
4549	fm_sum(:)=0.0
4550	DO ibvm=1, nbvmax
4551	! Leave the do loop if all sub areas are treated,
4552	! sub_area <= 0
4553	IF ( sub_area(ipts,ibvm) <= zero ) EXIT
4554	ip = sub_index(ipts,ibvm,1)
4555	jp = sub_index(ipts,ibvm,2)
4556	ifm=fmmap_i(ip,jp,ivma)
4557	fm_sum(ifm)=fm_sum(ifm)+sub_area(ipts,ibvm)
4558	ENDDO
4559
4560	! Whichever FM type has the most area, we use that one.
4561	desired_managed(ipts,ivma)=MAXLOC(fm_sum(:),1)
4562	ENDIF
4563	ENDDO
4564	ENDDO
4565
4566	DEALLOCATE(fmmap_i)
4567	DEALLOCATE(fmmap_r)
4568	DEALLOCATE(lat_lu,lon_lu)
4569	DEALLOCATE(lat_ful,lon_ful)
4570	DEALLOCATE(mask)
4571	IF(ALLOCATED(sub_index)) DEALLOCATE(sub_index)
4572	IF(ALLOCATED(sub_area)) DEALLOCATE(sub_area)
4573
4574	IF ( printlev >= 5 ) WRITE(numout,*) 'Leaving sapiens_forestry_read_desired_fm'
4575
4576	END SUBROUTINE sapiens_forestry_read_desired_fm
4577
4578	END MODULE sapiens_forestry

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

Download in other formats: