Context Navigation

function_library.f90 @ 7346

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

Line
1	! =================================================================================================================================
2	! MODULE : function_library
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 Collection of functions that are used throughout the ORCHIDEE code
10	!!
11	!!\n DESCRIPTION: Collection of modules to : (1) convert one variable into another i.e. basal area
12	!! to diameter, diamter to tree height, diameter to crown area, etc. (2) ...
13	!!
14	!! RECENT CHANGE(S): None
15	!!
16	!! REFERENCE(S) :
17	!!
18	!! SVN :
19	!! $HeadURL: svn://forge.ipsl.jussieu.fr/orchidee/branches/ORCHIDEE-DOFOCO/ORCHIDEE/src_stomate/stomate_prescribe.f90 $
20	!! $Date: 2013-01-04 16:50:56 +0100 (Fri, 04 Jan 2013) $
21	!! $Revision: 1126 $
22	!! \n
23	!_ ================================================================================================================================
24
25	MODULE function_library
26
27	! modules used:
28
29	!$ USE ioipsl_para
30	USE pft_parameters
31	USE constantes
32
33	IMPLICIT NONE
34
35	! private & public routines
36
37	PRIVATE
38	PUBLIC calculate_c0_alloc, wood_to_ba_eff_array, wood_to_ba_eff, wood_to_ba, &
39	wood_to_height_eff, wood_to_height, wood_to_qmheight, wood_to_dia_eff, &
40	wood_to_dia, wood_to_qmdia, wood_to_circ, wood_to_cn_array, &
41	wood_to_cn, wood_to_cn_eff, wood_to_cv, wood_to_cv_eff, wood_to_volume, &
42	biomass_to_lai, Nmax, Nmaxyield, Nminyield, distribute_mortality_biomass, &
43	check_biomass_sync, biomass_to_coupled_lai, check_area
44
45	CONTAINS
46
47	!! ================================================================================================================================
48	!! FUNCTION : biomass_to_coupled_lai
49	!!
50	!>\BRIEF Calculate the coupled_LAI based on biomass and veget of each pft
51	!!
52	!! DESCRIPTION : Calculates the lai that coupled with the atmosphere
53	!!
54	!!
55	!!
56	!! RECENT CHANGE(S): None
57	!!
58	!! RETURN VALUE : ::coupled_lai (m2/m2)
59	!!
60	!! REFERENCE(S) :
61	!!
62	!! FLOWCHART : None
63	!! \n
64	!_ ================================================================================================================================
65
66	FUNCTION biomass_to_coupled_lai(biomass_leaf, veget, pft)
67
68	!! 0. Variable and parameter declaration
69
70	!! 0.1 Input variables
71	INTEGER(i_std) :: pft !! PFT number (-)
72	REAL(r_std) :: veget !! 1-Pgap or the vegetation cover
73	REAL(r_std) :: biomass_leaf !! Biomass of an individual tree within a circ
74	!! class @tex $(m^{2} ind^{-1})$ @endtex
75
76	!! 0.2 Output variables
77	REAL(r_std) :: biomass_to_coupled_lai !! The fraction of the LAI that is assumed to interact
78	!! with the atmosphere @tex $(m^{2} m^{-2})$ @endtex
79	!! 0.3 Modified variables
80
81	!! 0.4 Local variables
82	REAL(r_std) :: lai !! Leaf area index
83	!! @tex $(m^{2} m^{-2})$ @endtex
84	!_ ================================================================================================================================
85
86	lai = biomass_to_lai(biomass_leaf, pft)
87
88	!+++CHECK+++
89	! In a closed canopy not all the leaves fully interact with the
90	! atmosphere because leaves can shelter each other. In a more
91	! open canopy most leaves can interact with the atmosphere.
92	! For the moment we are not clear how to calculate the fraction of
93	! the canopy that is coupled to the atmosphere. Also, based on the
94	! simulated evapotranspiration we need the entire canopy to be
95	! coupled to the atmosphere to obtain simulations which are close
96	! to the observations (Jung's upscaled product).
97	IF (veget .GT. min_sechiba) THEN
98
99	biomass_to_coupled_lai = lai
100
101	ELSE
102
103	! There are so many gaps that veget is extremely small. It is fair to
104	! assume that the whole canopy is coupled to the atmosphere.
105	biomass_to_coupled_lai = lai
106
107	ENDIF
108
109	END FUNCTION biomass_to_coupled_lai
110
111
112	!! ================================================================================================================================
113	!! FUNCTION : calculate_c0_alloc
114	!!
115	!>\BRIEF Calculate the baseline root vs sapwood allocation
116	!!
117	!! DESCRIPTION : Calculates the baseline root vs sapwood allocation based on the
118	!! parameters of the pipe model (hydraulic conductivities) and the
119	!! turnover of the different components
120	!!
121	!! RECENT CHANGE(S): None
122	!!
123	!! RETURN VALUE : ::c0_alloc (m)
124	!!
125	!! REFERENCE(S) :
126	!!
127	!! FLOWCHART : None
128	!! \n
129	!_ ================================================================================================================================
130
131	FUNCTION calculate_c0_alloc(pts, pft, tau_eff_root, tau_eff_sap)
132
133	!! 0. Variable and parameter declaration
134
135	!! 0.1 Input variables
136
137	INTEGER(i_std) :: pts !! Pixel number (-)
138	INTEGER(i_std) :: pft !! PFT number (-)
139	REAL(r_std) :: tau_eff_root !! Effective longivety for leaves (days)
140	REAL(r_std) :: tau_eff_sap !! Effective longivety for leaves (days)
141
142	!! 0.2 Output variables
143
144	REAL(r_std) :: calculate_c0_alloc !! quadratic mean height (m)
145
146	!! 0.3 Modified variables
147
148	!! 0.4 Local variables
149	REAL(r_std) :: sapwood_density
150	REAL(r_std) :: qm_dia !! quadratic mean diameter (m)
151
152	!_ ================================================================================================================================
153
154	!! 1. Calculate c0_alloc
155	IF ( is_tree(pft) ) THEN
156
157	sapwood_density = deux * pipe_density(pft) / kilo_to_unit
158	calculate_c0_alloc = sqrt(k_root(pft)/k_sap(pft)tau_eff_sap/tau_eff_rootsapwood_density)
159
160	! Grasses and croplands
161	ELSE
162
163	!+++CHECK+++
164	! Simply copied the same formulation as for trees but note
165	! that the sapwood in trees vs grasses and crops has a very
166	! meaning. In grasses and crops is structural carbon to ensure
167	! that the allocation works. In trees it really is the sapwood
168	sapwood_density = deux * pipe_density(pft) / kilo_to_unit
169	calculate_c0_alloc = sqrt(k_root(pft)/k_sap(pft)tau_eff_sap/tau_eff_rootsapwood_density)
170	!+++++++++++
171
172	ENDIF ! is_tree(j)
173
174	END FUNCTION calculate_c0_alloc
175
176
177
178	!! ================================================================================================================================
179	!! FUNCTION : wood_to_ba_eff_array
180	!!
181	!>\BRIEF Calculate the effective basal area from woody biomass making use of allometric relationships.
182	!!
183	!! DESCRIPTION : Calculate basal area of an individual tree from the woody biomass of that tree making
184	!! use of allometric relationships. Effective basal area accounts for both above and below ground carbon
185	!! and is the basis for the application of the rule of Deleuze and Dhote
186	!!
187	!! RECENT CHANGE(S): None
188	!!
189	!! RETURN VALUE : effective basal area (m2 ind-1)
190	!!
191	!! REFERENCE(S) :
192	!!
193	!! FLOWCHART : None
194	!! \n
195	!_ ================================================================================================================================
196
197	FUNCTION wood_to_ba_eff_array(biomass_temp, npts, pft)
198
199	!! 0. Variable and parameter declaration
200
201	!! 0.1 Input variables
202
203	INTEGER(i_std) :: pft !! PFT number (-)
204	INTEGER(i_std) :: npts !! Number of pixels
205	REAL(r_std), DIMENSION(:,:,:) :: biomass_temp !! Biomass of an individual tree within a circ
206	!! class @tex $(m^{2} ind^{-1})$ @endtex
207
208	!! 0.2 Output variables
209
210	REAL(r_std), DIMENSION(npts,ncirc) :: wood_to_ba_eff_array !! Effective basal area of an individual tree within a circ
211	!! class @tex $(m^{2} ind^{-1})$ @endtex
212
213	!! 0.3 Modified variables
214
215	!! 0.4 Local variables
216	INTEGER(i_std) :: l !! Index
217	REAL(r_std), DIMENSION(npts) :: woodmass_ind !! Woodmass of an individual tree
218	!! @tex $(gC ind{-1})$ @endtex
219
220	!_ ================================================================================================================================
221
222	!! 1. Calculate effective basal area from woodmass
223
224	IF ( is_tree(pft) ) THEN
225
226	DO l = 1,ncirc
227
228	! Woodmass of an individual tree. Note that for the effective basal area
229	! both the above and belowground biomass are used.
230	woodmass_ind(:) = biomass_temp(:,l,isapabove) + biomass_temp(:,l,isapbelow) + &
231	biomass_temp(:,l,iheartabove) + biomass_temp(:,l,iheartbelow)
232
233	! Basal area of that individual (m2 ind-1)
234	wood_to_ba_eff_array(:,l) = (pi/4*(woodmass_ind(:)/&
235	(tree_ff(pft)pipe_density(pft)pipe_tune2(pft))) &
236	(2./pipe_tune3(pft)))(pipe_tune3(pft)/(pipe_tune3(pft)+2))
237
238	ENDDO
239
240	ELSE
241
242	WRITE(numout,*) 'function wood_to_ba_eff_array is not defined for this PFT'
243	WRITE(numout,*) 'pft ',pft
244	CALL ipslerr_p (3,'wood_to_ba_eff_array', &
245	'wood_to_ba_eff_array is not defined for this PFT.', &
246	'See the output file for more details.','')
247
248
249	ENDIF
250
251	END FUNCTION wood_to_ba_eff_array
252
253
254	!! ================================================================================================================================
255	!! FUNCTION : wood_to_ba_eff
256	!!
257	!>\BRIEF Calculate effective basal area from woody biomass making use of allometric relationships
258	!!
259	!! DESCRIPTION : Calculate basal area of an individual tree from the woody biomass of that tree making
260	!! use of allometric relationships. Effective basal area accounts for both above and below ground carbon
261	!! and is the basis for the application of the rule of Deleuze and Dhote.
262	!! (i) woodmass = tree_ff * pipe_densitybaheight
263	!! (ii) height = pipe_tune2 * sqrt(4/piba) * pipe_tune_3
264	!!
265	!! RECENT CHANGE(S): None
266	!!
267	!! RETURN VALUE : effective basal area (m2 ind-1)
268	!!
269	!! REFERENCE(S) :
270	!!
271	!! FLOWCHART : None
272	!! \n
273	!_ ================================================================================================================================
274
275	FUNCTION wood_to_ba_eff(biomass_temp, pft)
276
277	!! 0. Variable and parameter declaration
278
279	!! 0.1 Input variables
280
281	INTEGER(i_std) :: pft !! PFT number (-)
282	REAL(r_std), DIMENSION(:,:) :: biomass_temp !! Biomass of an individual tree within a circ
283	!! class @tex $(m^{2} ind^{-1})$ @endtex
284
285	!! 0.2 Output variables
286
287	REAL(r_std), DIMENSION(ncirc) :: wood_to_ba_eff !! Effective basal area of an individual tree within a circ
288	!! class @tex $(m^{2} ind^{-1})$ @endtex
289
290	!! 0.3 Modified variables
291
292	!! 0.4 Local variables
293	INTEGER(i_std) :: l !! Index
294	REAL(r_std) :: woodmass_ind !! Woodmass of an individual tree
295	!! @tex $(gC ind^{-1})$ @endtex
296
297	!_ ================================================================================================================================
298
299	!! 1. Calculate basal area from woodmass
300
301	IF ( is_tree(pft) ) THEN
302
303	DO l = 1,ncirc
304
305	! Woodmass of an individual tree
306	woodmass_ind = biomass_temp(l,isapabove) + biomass_temp(l,isapbelow) + &
307	biomass_temp(l,iheartabove) + biomass_temp(l,iheartbelow)
308
309	! Basal area of that individual (m2 ind-1)
310	wood_to_ba_eff(l) = (pi/4(woodmass_ind/(tree_ff(pft)pipe_density(pft)*pipe_tune2(pft))) &
311	(2./pipe_tune3(pft)))(pipe_tune3(pft)/(pipe_tune3(pft)+2))
312
313	ENDDO
314
315	ELSE
316
317	WRITE(numout,*) 'pft ',pft
318	CALL ipslerr_p (3,'wood_to_ba_eff', &
319	'wood_to_ba_eff is not defined for this PFT.', &
320	'See the output file for more details.','')
321
322	ENDIF
323
324	END FUNCTION wood_to_ba_eff
325
326
327
328	!! ================================================================================================================================
329	!! FUNCTION : wood_to_ba
330	!!
331	!>\BRIEF Calculate basal area from woody biomass making use of allometric relationships
332	!!
333	!! DESCRIPTION : Calculate basal area of an individual tree from the woody biomass of that tree making
334	!! use of allometric relationships given below. Here basal area is defined in line with its classical
335	!! forestry meaning.
336	!! (i) woodmass = tree_ff * pipe_densitybaheight
337	!! (ii) height = pipe_tune2 * sqrt(4/piba) * pipe_tune_3
338	!!
339	!! RECENT CHANGE(S): None
340	!!
341	!! RETURN VALUE : basal area (m2 ind-1)
342	!!
343	!! REFERENCE(S) :
344	!!
345	!! FLOWCHART : None
346	!! \n
347	!_ ================================================================================================================================
348
349	FUNCTION wood_to_ba(biomass_temp, pft)
350
351	!! 0. Variable and parameter declaration
352
353	!! 0.1 Input variables
354
355	INTEGER(i_std) :: pft !! PFT number (-)
356	REAL(r_std), DIMENSION(:,:) :: biomass_temp !! Biomass of an individual tree within a circ
357	!! class @tex $(m^{2} ind^{-1})$ @endtex
358
359	!! 0.2 Output variables
360
361	REAL(r_std), DIMENSION(ncirc) :: wood_to_ba !! Basal area of an individual tree within a circ
362	!! class @tex $(m^{2} ind^{-1})$ @endtex
363
364	!! 0.3 Modified variables
365
366	!! 0.4 Local variables
367	INTEGER(i_std) :: l !! Index
368	REAL(r_std) :: woodmass_ind !! Woodmass of an individual tree
369	!! @tex $(gC ind^{-1})$ @endtex
370
371	!_ ================================================================================================================================
372
373	!! 1. Calculate basal area from woodmass
374
375	IF ( is_tree(pft) ) THEN
376
377	DO l = 1,ncirc
378
379	! Woodmass of an individual tree
380	woodmass_ind = biomass_temp(l,iheartabove) + biomass_temp(l,isapabove)
381
382
383	! Basal area of that individual (m2 ind-1)
384	wood_to_ba(l) = (pi/4(woodmass_ind/(tree_ff(pft)pipe_density(pft)*pipe_tune2(pft))) &
385	(2./pipe_tune3(pft)))(pipe_tune3(pft)/(pipe_tune3(pft)+2))
386
387	ENDDO
388
389	ELSE
390
391	WRITE(numout,*) 'pft ',pft
392	CALL ipslerr_p (3,'wood_to_ba', &
393	'wood_to_ba is not defined for this PFT.', &
394	'See the output file for more details.','')
395
396	ENDIF
397
398	END FUNCTION wood_to_ba
399
400
401	!! ================================================================================================================================
402	!! FUNCTION : wood_to_height_eff
403	!!
404	!>\BRIEF Calculate the effective tree height from woody biomass making use of allometric relationships
405	!!
406	!! DESCRIPTION : Calculate the effective height of an individual tree from the woody biomass of that tree making
407	!! use of allometric relationships. Effective height makes use of both above and belowground biomass and is
408	!! used in the calculation of the allocation according to deleuze and dhote, hydraulic architecture because also
409	!! the height of the belowground part should be included.
410	!! (i) height(:) = pipe_tune2(j)(4/piba(:))**(pipe_tune3(j)/2) and
411	!! (ii) woodmass_ind = tree_ffpipe_densityba*height
412	!!
413	!! RECENT CHANGE(S): None
414	!!
415	!! RETURN VALUE : height (m)
416	!!
417	!! REFERENCE(S) :
418	!!
419	!! FLOWCHART : None
420	!! \n
421	!_ ================================================================================================================================
422
423	FUNCTION wood_to_height_eff(biomass_temp, pft)
424
425	!! 0. Variable and parameter declaration
426
427	!! 0.1 Input variables
428
429	INTEGER(i_std) :: pft !! PFT number (-)
430	REAL(r_std), DIMENSION(:,:) :: biomass_temp !! Biomass of an individual tree within a circ
431	!! class @tex $(m^{2} ind^{-1})$ @endtex
432
433	!! 0.2 Output variables
434
435	REAL(r_std), DIMENSION(ncirc) :: wood_to_height_eff !! Effective height of an individual tree within a circ
436	!! class (m)
437
438	!! 0.3 Modified variables
439
440	!! 0.4 Local variables
441	INTEGER(i_std) :: l !! Index
442	REAL(r_std) :: woodmass_ind !! Woodmass of an individual tree
443	!! @tex $(gC ind{-1})$ @endtex
444
445	CHARACTER(len=256) :: tmp_var_name !! temporal variable for a text reading (from run.def)
446	REAL(r_std) :: tmp_ratio !! temporal variable for adjusting the height
447	!_ ================================================================================================================================
448
449	!! 1. Calculate height from woodmass
450
451	IF ( is_tree(pft) ) THEN
452
453	DO l = 1,ncirc
454
455	! Woodmass of an individual tree. Both above and belowground biomass are used
456	woodmass_ind = biomass_temp(l,isapabove) + biomass_temp(l,isapbelow) + &
457	biomass_temp(l,iheartabove) + biomass_temp(l,iheartbelow)
458
459	! Height of that individual
460	wood_to_height_eff(l) = (((woodmass_ind/(tree_ff(pft)pipe_density(pft))4/pi)**(pipe_tune3(pft)/2))&
461	pipe_tune2(pft))*(1/(pipe_tune3(pft)/2+1))
462
463	!+++++++ TEMP ++++++++++
464	! This code is only used evaluation of the performance of the multi-layer energy budget.
465	! To reduce the complexity of the tests we want to impose the LAI and its vertical distribution.
466	! The solution is not very elegant but it works. (To start the simulation, a retsart file with height stand is suggested)
467	! This part of code use the CIRC_HEIGHT_RATIO to adjust the original heigth in each circonferance class.
468	! By doing so, the canopy sturcture can be adjust to match the site level observations.
469	IF (ld_fake_height) THEN
470	! WRITE(numout,*)'+++ fake wood to height function +++'
471	!set initial the tmp_ratio reaul to 1.0
472	tmp_ratio=1.0
473	WRITE(tmp_var_name, '(A19,I5.5)') "CIRC_HEIGHT_RATIO__", l
474	CALL getin_p(TRIM(tmp_var_name), tmp_ratio)
475	wood_to_height_eff(l) =tmp_ratio*wood_to_height_eff(l)
476	! WRITE(numout,*) 'SET CIRC_', l, '_HEIGHT = ', wood_to_height_eff(l)
477	ENDIF
478	!++++++++++++++++++++++
479
480	ENDDO
481
482	ELSE
483
484	WRITE(numout,*) 'pft ',pft
485	CALL ipslerr_p (3,'wood_to_height_eff', &
486	'wood_to_height_eff is not defined for this PFT.', &
487	'See the output file for more details.','')
488
489	ENDIF
490
491	END FUNCTION wood_to_height_eff
492
493
494	!! ================================================================================================================================
495	!! FUNCTION : wood_to_height
496	!!
497	!>\BRIEF Calculate tree height from woody biomass making use of allometric relationships
498	!!
499	!! DESCRIPTION : Calculate height of an individual tree from the woody biomass of that tree making
500	!! use of allometric relationships. This is the height used in forestry and for calculating the aerodynamic
501	!! interactions
502	!! (i) height(:) = pipe_tune2(j)(4/piba(:))**(pipe_tune3(j)/2) and
503	!! (ii) woodmass_ind = tree_ffpipe_densityba*height
504	!!
505	!! RECENT CHANGE(S): None
506	!!
507	!! RETURN VALUE : height (m)
508	!!
509	!! REFERENCE(S) :
510	!!
511	!! FLOWCHART : None
512	!! \n
513	!_ ================================================================================================================================
514
515	FUNCTION wood_to_height(biomass_temp, pft)
516
517	!! 0. Variable and parameter declaration
518
519	!! 0.1 Input variables
520
521	INTEGER(i_std) :: pft !! PFT number (-)
522	REAL(r_std), DIMENSION(:,:) :: biomass_temp !! Biomass of an individual tree within a circ
523	!! class @tex $(m^{2} ind^{-1})$ @endtex
524
525	!! 0.2 Output variables
526
527	REAL(r_std), DIMENSION(ncirc) :: wood_to_height !! Height of an individual tree within a circ
528	!! class (m)
529
530	!! 0.3 Modified variables
531
532	!! 0.4 Local variables
533	INTEGER(i_std) :: l !! Index
534	REAL(r_std) :: woodmass_ind !! Woodmass of an individual tree
535	!! @tex $(gC ind{-1})$ @endtex
536
537	!_ ================================================================================================================================
538
539	!! 1. Calculate height from woodmass
540
541	IF ( is_tree(pft) ) THEN
542
543	DO l = 1,ncirc
544
545	! Woodmass of an individual tree (only the aboveground component)
546	woodmass_ind = biomass_temp(l,isapabove) + biomass_temp(l,iheartabove)
547
548	! Height of that individual
549	wood_to_height(l) = (((woodmass_ind/(tree_ff(pft)pipe_density(pft))4/pi)**(pipe_tune3(pft)/2))&
550	pipe_tune2(pft))*(1/(pipe_tune3(pft)/2+1))
551
552	ENDDO
553
554	ELSE
555
556	WRITE(numout,*) 'pft ',pft
557	CALL ipslerr_p (3,'wood_to_height', &
558	'wood_to_height is not defined for this PFT.', &
559	'See the output file for more details.','')
560
561	ENDIF
562
563	END FUNCTION wood_to_height
564
565
566	!! ================================================================================================================================
567	!! FUNCTION : wood_to_qmheight
568	!!
569	!>\BRIEF Calculate the quadratic mean height from the biomass
570	!!
571	!! DESCRIPTION : Calculates the quadratic mean height from the biomass
572	!!
573	!! RECENT CHANGE(S): None
574	!!
575	!! RETURN VALUE : ::qm_height (m)
576	!!
577	!! REFERENCE(S) :
578	!!
579	!! FLOWCHART : None
580	!! \n
581	!_ ================================================================================================================================
582
583	FUNCTION wood_to_qmheight(biomass_temp, ind, pft)
584
585	!! 0. Variable and parameter declaration
586
587	!! 0.1 Input variables
588
589	INTEGER(i_std) :: pft !! PFT number (-)
590	REAL(r_std), DIMENSION(ncirc,nparts) :: biomass_temp !! Biomass of the leaves @tex $(gC m^{-2})$ @endtex
591	REAL(r_std), DIMENSION(ncirc) :: ind !! Number of individuals @tex $(m^{-2})$ @endtex
592
593
594	!! 0.2 Output variables
595
596	REAL(r_std) :: wood_to_qmheight !! quadratic mean height (m)
597
598	!! 0.3 Modified variables
599
600	!! 0.4 Local variables
601	REAL(r_std), DIMENSION(ncirc) :: circ_class_ba !! basal area for each circ_class @tex $(m^{2})$ @endtex
602	REAL(r_std) :: qm_dia !! quadratic mean diameter (m)
603
604	!_ ================================================================================================================================
605
606	!! 1. Calculate qm_height from the biomass
607	IF ( is_tree(pft) ) THEN
608
609	! Basal area at the tree level (m2 tree-1)
610	circ_class_ba(:) = wood_to_ba(biomass_temp(:,:),pft)
611
612	IF (SUM(ind(:)) .NE. zero) THEN
613
614	qm_dia = SQRT( 4/piSUM(circ_class_ba(:)ind(:))/SUM(ind(:)) )
615
616	ELSE
617
618	qm_dia = zero
619
620	ENDIF
621
622	wood_to_qmheight = pipe_tune2(pft)(qm_dia*pipe_tune3(pft))
623
624
625	! Grasses and croplands
626	ELSE
627
628	! Calculate height as a function of the leaf and structural biomass. Use structural
629	! biomass to make sure that the grasslands have a roughness length during the winter
630	! If the biomass increases, vegetation height will increase as well. Divide by
631	! ind(ipts,j) to obtain the height of the individual. biomass(ileaf) is in gC m-2
632	! whereas qm is the height of the individual.
633	IF (SUM(ind(:)) .NE. zero) THEN
634
635	wood_to_qmheight = SUM(biomass_temp(:,ileaf) + biomass_temp(:,isapabove)) / &
636	SUM(ind(:)) * sla(pft) * lai_to_height(pft)
637
638	ELSE
639
640	wood_to_qmheight = zero
641
642	ENDIF
643
644	ENDIF ! is_tree(j)
645
646	END FUNCTION wood_to_qmheight
647
648
649
650	!! ================================================================================================================================
651	!! FUNCTION : wood_to_dia_eff
652	!!
653	!>\BRIEF Calculate effective diameter from woody biomass making use of allometric relationships
654	!!
655	!! DESCRIPTION : Calculate the effective diameter of an individual tree from the woody biomass of that tree making
656	!! use of allometric relationships. Effective diameter accounts for both above and belowground biomass.
657	!! (i) woodmass_ind = tree_ff * pipe_density * height * pi/4dia*2
658	!! (ii) height = pipe_tune2 * dia * pipe_tune3
659	!!
660	!! RECENT CHANGE(S): None
661	!!
662	!! RETURN VALUE : diameter (m)
663	!!
664	!! REFERENCE(S) :
665	!!
666	!! FLOWCHART : None
667	!! \n
668	!_ ================================================================================================================================
669
670	FUNCTION wood_to_dia_eff(biomass_temp, pft)
671
672	!! 0. Variable and parameter declaration
673
674	!! 0.1 Input variables
675
676	INTEGER(i_std) :: pft !! PFT number (-)
677	REAL(r_std), DIMENSION(:,:) :: biomass_temp !! Biomass of an individual tree within a circ
678	!! class @tex $(m^{2} ind^{-1})$ @endtex
679
680	!! 0.2 Output variables
681
682	REAL(r_std), DIMENSION(ncirc) :: wood_to_dia_eff !! Diameter of an individual tree within a circ
683	!! class (m)
684
685	!! 0.3 Modified variables
686
687	!! 0.4 Local variables
688	INTEGER(i_std) :: l !! Index
689	REAL(r_std) :: woodmass_ind !! Woodmass of an individual tree
690	!! @tex $(gC ind^{-1})$ @endtex
691
692	!_ ================================================================================================================================
693
694	!! 1. Calculate basal area from woodmass
695
696	IF ( is_tree(pft) ) THEN
697
698	DO l = 1,ncirc
699
700	! Woodmass of an individual tree
701	woodmass_ind = biomass_temp(l,isapabove) + biomass_temp(l,isapbelow) + &
702	biomass_temp(l,iheartabove) + biomass_temp(l,iheartbelow)
703
704	! Basal area of that individual (m2 ind-1)
705	wood_to_dia_eff(l) = (4/piwoodmass_ind/(tree_ff(pft)pipe_density(pft)pipe_tune2(pft))) * &
706	(1./(2+pipe_tune3(pft)))
707
708	ENDDO
709
710	ELSE
711
712	WRITE(numout,*) 'pft ',pft
713	CALL ipslerr_p (3,'wood_to_dia_eff', &
714	'wood_to_dia_eff is not defined for this PFT.', &
715	'See the output file for more details.','')
716
717	ENDIF
718
719	END FUNCTION wood_to_dia_eff
720
721
722
723	!! ================================================================================================================================
724	!! FUNCTION : wood_to_dia
725	!!
726	!>\BRIEF Calculate diameter from woody biomass making use of allometric relationships
727	!!
728	!! DESCRIPTION : Calculate diameter of an individual tree from the woody biomass of that tree making
729	!! use of allometric relationships. Makes only use of the aboveground biomass and relates to the
730	!! typical forestry diameter (but not normalized at 1.3 m)
731	!! (i) woodmass_ind = tree_ff * pipe_density * height * pi/4dia*2
732	!! (ii) height = pipe_tune2 * dia * pipe_tune3
733	!!
734	!! RECENT CHANGE(S): None
735	!!
736	!! RETURN VALUE : diameter (m)
737	!!
738	!! REFERENCE(S) :
739	!!
740	!! FLOWCHART : None
741	!! \n
742	!_ ================================================================================================================================
743
744	FUNCTION wood_to_dia(biomass_temp, pft)
745
746	!! 0. Variable and parameter declaration
747
748	!! 0.1 Input variables
749
750	INTEGER(i_std) :: pft !! PFT number (-)
751	REAL(r_std), DIMENSION(:,:) :: biomass_temp !! Biomass of an individual tree within a circ
752	!! class @tex $(m^{2} ind^{-1})$ @endtex
753
754	!! 0.2 Output variables
755
756	REAL(r_std), DIMENSION(ncirc) :: wood_to_dia !! Diameter of an individual tree within a circ
757	!! class (m)
758
759	!! 0.3 Modified variables
760
761	!! 0.4 Local variables
762	INTEGER(i_std) :: l !! Index
763	REAL(r_std) :: woodmass_ind !! Woodmass of an individual tree
764	!! @tex $(gC ind^{-1})$ @endtex
765
766	!_ ================================================================================================================================
767
768	!! 1. Calculate basal area from woodmass
769
770	IF ( is_tree(pft) ) THEN
771
772	DO l = 1,ncirc
773
774	! Woodmass of an individual tree
775	woodmass_ind = biomass_temp(l,isapabove) + biomass_temp(l,iheartabove)
776
777	! Basal area of that individual (m2 ind-1)
778	wood_to_dia(l) = (4/piwoodmass_ind/(tree_ff(pft)pipe_density(pft)pipe_tune2(pft))) * &
779	(1./(2+pipe_tune3(pft)))
780
781	ENDDO
782
783	ELSE
784
785	WRITE(numout,*) 'pft ',pft
786	CALL ipslerr_p (3,'wood_to_dia', &
787	'wood_to_dia is not defined for this PFT.', &
788	'See the output file for more details.','')
789
790	ENDIF
791
792	END FUNCTION wood_to_dia
793
794
795	!! ================================================================================================================================
796	!! FUNCTION : wood_to_qmdia
797	!!
798	!>\BRIEF Calculate the quadratic mean diameter from the biomass
799	!!
800	!! DESCRIPTION : Calculates the quadratic mean diameter from the aboveground biomss
801	!!
802	!! RECENT CHANGE(S): None
803	!!
804	!! RETURN VALUE : ::qm_dia (m)
805	!!
806	!! REFERENCE(S) :
807	!!
808	!! FLOWCHART : None
809	!! \n
810	!_ ================================================================================================================================
811
812	FUNCTION wood_to_qmdia(biomass_temp, ind, pft)
813
814	!! 0. Variable and parameter declaration
815
816	!! 0.1 Input variables
817
818	INTEGER(i_std) :: pft !! PFT number (-)
819	REAL(r_std), DIMENSION(ncirc,nparts) :: biomass_temp !! Biomass of the leaves @tex $(gC m^{-2})$ @endtex
820	REAL(r_std), DIMENSION(ncirc) :: ind !! Number of individuals @tex $(m^{-2})$ @endtex
821
822	!! 0.2 Output variables
823
824	REAL(r_std) :: wood_to_qmdia !! quadratic mean diameter (m)
825
826	!! 0.3 Modified variables
827
828	!! 0.4 Local variables
829	REAL(r_std), DIMENSION(ncirc) :: circ_class_ba !! basal area for each circ_class @tex $(m^{2})$ @endtex
830
831	!_ ================================================================================================================================
832
833	!! 1. Calculate qm_dia from the biomass
834	IF ( is_tree(pft) ) THEN
835
836	! Basal area at the tree level (m2 tree-1)
837	circ_class_ba(:) = wood_to_ba(biomass_temp(:,:),pft)
838
839	IF (SUM(ind(:)) .NE. zero) THEN
840
841	wood_to_qmdia = SQRT( 4/piSUM(circ_class_ba(:)ind(:))/SUM(ind(:)) )
842
843	ELSE
844
845	wood_to_qmdia = zero
846
847	ENDIF
848
849
850	! Grasses and croplands
851	ELSE
852
853	wood_to_qmdia = zero
854
855	ENDIF ! is_tree(pft)
856
857	END FUNCTION wood_to_qmdia
858
859
860	!! ================================================================================================================================
861	!! FUNCTION : wood_to_circ
862	!!
863	!>\BRIEF Calculate circumference from woody biomass making use of allometric relationships
864	!!
865	!! DESCRIPTION : All this does it computer the diameter using a different routine, and then
866	!! convert that into a circumference.
867	!!
868	!! RECENT CHANGE(S): None
869	!!
870	!! RETURN VALUE : circumference (m)
871	!!
872	!! REFERENCE(S) :
873	!!
874	!! FLOWCHART : None
875	!! \n
876	!_ ================================================================================================================================
877
878	FUNCTION wood_to_circ(biomass_temp, pft)
879
880	!! 0. Variable and parameter declaration
881
882	!! 0.1 Input variables
883
884	INTEGER(i_std) :: pft !! PFT number (-)
885	REAL(r_std), DIMENSION(:,:) :: biomass_temp !! Biomass of an individual tree within a circ
886	!! class @tex $(m^{2} ind^{-1})$ @endtex
887
888	!! 0.2 Output variables
889
890	REAL(r_std), DIMENSION(ncirc) :: wood_to_circ !! Circumference of an individual tree within a circ
891	!! class (m)
892
893	!! 0.3 Modified variables
894
895	!! 0.4 Local variables
896
897	!_ ================================================================================================================================
898
899	!! 1. Calculate diameter from woodmass
900
901	wood_to_circ(:)=val_exp
902
903	wood_to_circ(:)=wood_to_dia(biomass_temp(:,:),pft)
904
905	! convert to a circumference (m)
906	wood_to_circ(:) = wood_to_circ(:)*pi
907
908	END FUNCTION wood_to_circ
909
910
911
912	!! ================================================================================================================================
913	!! FUNCTION : wood_to_cn_array
914	!!
915	!>\BRIEF Calculate crown area from woody biomass making use of allometric relationships
916	!!
917	!! DESCRIPTION : Calculate crown area of an individual tree from the woody biomass of that tree
918	!! making use of allometric relationship that relates crown area (cn) to diameter (dia) as
919	!! pipe_tune1dia*pipe_tune_exp_coeff where the pipe_tune parameters are pft-specific.
920	!! (i) Basal area is written as a function of wood_mass: woodmass_ind = tree_ffpipe_densityba*height
921	!! (ii) height = pipe_tune2sqrt(4/piba)**pipe_tune3
922	!! (iii) cn = pipe_tune1sqrt(4/piba)**pipe_tune_exp_coeff
923	!!
924	!! RECENT CHANGE(S): None
925	!!
926	!! RETURN VALUE : crown area (m2 ind-1)
927	!!
928	!! REFERENCE(S) :
929	!!
930	!! FLOWCHART : None
931	!! \n
932	!_ ================================================================================================================================
933
934	FUNCTION wood_to_cn_array(biomass_temp, npts, pft)
935
936	!! 0. Variable and parameter declaration
937
938	!! 0.1 Input variables
939
940	INTEGER(i_std) :: pft !! PFT number (-)
941	INTEGER(i_std) :: npts !! Pixel(s), this variable defines the dimensions of ba
942	REAL(r_std), DIMENSION(:,:,:) :: biomass_temp !! Biomass of an individual tree within a circ
943	!! class @tex $(m^{2} ind^{-1})$ @endtex
944
945
946	!! 0.2 Output variables
947
948	REAL(r_std), DIMENSION(npts,ncirc) :: wood_to_cn_array !! Crown area of an individual tree @tex $(m^{2} ind{-1})$ @endtex
949
950	!! 0.3 Modified variables
951
952	!! 0.4 Local variables
953
954	INTEGER(i_std) :: l !! index
955	REAL(r_std), DIMENSION(npts) :: woodmass_ind !! Woodmass of an individual tree @tex $(gC ind{-1})$ @endtex
956
957	!_ ================================================================================================================================
958
959	!! 1. Calculate crown area from basal area
960
961	IF ( is_tree(pft) ) THEN
962
963	DO l = 1,ncirc
964
965	! Woodmass of an individual tree
966	woodmass_ind(:) = biomass_temp(:,l,isapabove) + biomass_temp(:,l,iheartabove)
967
968	! Crown area of that individual
969	! cn = pipe_tune1sqrt(4/piba)**pipe_tune_exp_coeff
970	wood_to_cn_array(:,l) = pipe_tune1(pft) * SQRT( 4/pi(pi/4(woodmass_ind(:)/&
971	(tree_ff(pft)pipe_density(pft)pipe_tune2(pft)))(2./pipe_tune3(pft)))&
972	(pipe_tune3(pft)/(pipe_tune3(pft)+2)) ) ** pipe_tune_exp_coeff(pft)
973
974	ENDDO
975
976	ELSE
977
978	WRITE(numout,*) 'pft ',pft
979	CALL ipslerr_p (3,'wood_to_cn_array', &
980	'wood_to_cn_array is not defined for this PFT.', &
981	'See the output file for more details.','')
982
983	ENDIF
984
985	END FUNCTION wood_to_cn_array
986
987
988	!! ================================================================================================================================
989	!! FUNCTION : wood_to_cn
990	!!
991	!>\BRIEF Calculate crown area from woody biomass making use of allometric relationships
992	!!
993	!! DESCRIPTION : Calculate crown area of an individual tree from the woody biomass of that tree
994	!! making use of allometric relationship that relates crown area (cn) to diameter (dia) as
995	!! pipe_tune1dia*pipe_tune_exp_coeff where the pipe_tune parameters are pft-specific.
996	!! (i) Basal area is written as a function of wood_mass: woodmass_ind = tree_ffpipe_densityba*height
997	!! (ii) height = pipe_tune2sqrt(4/piba)**pipe_tune3
998	!! (iii) cn = pipe_tune1sqrt(4/piba)**pipe_tune_exp_coeff
999	!!
1000	!! RECENT CHANGE(S): None
1001	!!
1002	!! RETURN VALUE : crown area (m2 ind-1)
1003	!!
1004	!! REFERENCE(S) :
1005	!!
1006	!! FLOWCHART : None
1007	!! \n
1008	!_ ================================================================================================================================
1009
1010	FUNCTION wood_to_cn(biomass_temp, pft)
1011
1012	!! 0. Variable and parameter declaration
1013
1014	!! 0.1 Input variables
1015
1016	INTEGER(i_std) :: pft !! PFT number (-)
1017	REAL(r_std), DIMENSION(:,:) :: biomass_temp !! biomass of an individual tree within a circ
1018	!! class @tex $(m^{2} ind^{-1})$ @endtex
1019
1020	!! 0.2 Output variables
1021
1022	REAL(r_std), DIMENSION(ncirc) :: wood_to_cn !! Crown area of an individual tree within a circ
1023	!! class @tex $(m^{2} ind-1)$ @endtex
1024
1025	!! 0.3 Modified variables
1026
1027	!! 0.4 Local variables
1028	INTEGER(i_std) :: l !! Index
1029	REAL(r_std) :: woodmass_ind !! Woodmass of an individual tree
1030	!! @tex $(gC ind{-1})$ @endtex
1031
1032	!_ ================================================================================================================================
1033
1034	!! 1. Calculate crown area from woodmass
1035
1036	IF ( is_tree(pft) ) THEN
1037
1038	DO l = 1,ncirc
1039
1040	! Woodmass of an individual tree
1041	woodmass_ind = biomass_temp(l,isapabove) + biomass_temp(l,iheartabove)
1042
1043	! Crown area of that individual
1044	! cn = pipe_tune1sqrt(4/piba)**pipe_tune_exp_coeff
1045	wood_to_cn(l) = pipe_tune1(pft) * SQRT( 4/pi(pi/4(woodmass_ind/&
1046	(tree_ff(pft)pipe_density(pft)pipe_tune2(pft)))(2./pipe_tune3(pft)))&
1047	(pipe_tune3(pft)&
1048	/(pipe_tune3(pft)+2)) ) ** pipe_tune_exp_coeff(pft)
1049
1050	ENDDO
1051
1052	ELSE
1053
1054	WRITE(numout,*) 'pft ',pft
1055	CALL ipslerr_p (3,'wood_to_cn', &
1056	'wood_to_cn is not defined for this PFT.', &
1057	'See the output file for more details.','')
1058
1059	ENDIF
1060
1061	END FUNCTION wood_to_cn
1062
1063	!! ================================================================================================================================
1064	!! FUNCTION : wood_to_cn_eff
1065	!!
1066	!>\BRIEF Calculate crown area from woody biomass making use of allometric relationships
1067	!!
1068	!! DESCRIPTION : Calculate crown area of an individual tree from the woody biomass of that tree
1069	!! making use of allometric relationship that relates crown area (cn) to diameter (dia) as
1070	!! pipe_tune1dia*pipe_tune_exp_coeff where the pipe_tune parameters are pft-specific.
1071	!! (i) Basal area is written as a function of wood_mass: woodmass_ind = tree_ffpipe_densityba*height
1072	!! (ii) height = pipe_tune2sqrt(4/piba)**pipe_tune3
1073	!! (iii) cn = pipe_tune1sqrt(4/piba)**pipe_tune_exp_coeff
1074	!!
1075	!! NOTE: This is different from wood_to_cn because this uses both the aboveground and
1076	!! belowground wood mass.
1077	!!
1078	!! RECENT CHANGE(S): None
1079	!!
1080	!! RETURN VALUE : crown area (m2 ind-1)
1081	!!
1082	!! REFERENCE(S) :
1083	!!
1084	!! FLOWCHART : None
1085	!! \n
1086	!_ ================================================================================================================================
1087
1088	FUNCTION wood_to_cn_eff(biomass_temp, pft)
1089
1090	!! 0. Variable and parameter declaration
1091
1092	!! 0.1 Input variables
1093
1094	INTEGER(i_std) :: pft !! PFT number (-)
1095	REAL(r_std), DIMENSION(:,:) :: biomass_temp !! biomass of an individual tree within a circ
1096	!! class @tex $(m^{2} ind^{-1})$ @endtex
1097
1098	!! 0.2 Output variables
1099
1100	REAL(r_std), DIMENSION(ncirc) :: wood_to_cn_eff !! Crown area of an individual tree within a circ
1101	!! class @tex $(m^{2} ind-1)$ @endtex
1102
1103	!! 0.3 Modified variables
1104
1105	!! 0.4 Local variables
1106	INTEGER(i_std) :: l !! Index
1107	REAL(r_std) :: woodmass_ind !! Woodmass of an individual tree
1108	!! @tex $(gC ind{-1})$ @endtex
1109
1110	!_ ================================================================================================================================
1111
1112	!! 1. Calculate crown area from woodmass
1113
1114	IF ( is_tree(pft) ) THEN
1115
1116	DO l = 1,ncirc
1117
1118	! Woodmass of an individual tree
1119	woodmass_ind = biomass_temp(l,isapabove) + biomass_temp(l,iheartabove) + &
1120	biomass_temp(l,isapbelow) + biomass_temp(l,iheartbelow)
1121
1122	! Crown area of that individual
1123	! cn = pipe_tune1sqrt(4/piba)**pipe_tune_exp_coeff
1124	wood_to_cn_eff(l) = pipe_tune1(pft) * SQRT( 4/pi(pi/4(woodmass_ind/&
1125	(tree_ff(pft)pipe_density(pft)pipe_tune2(pft)))(2./pipe_tune3(pft)))&
1126	(pipe_tune3(pft)&
1127	/(pipe_tune3(pft)+2)) ) ** pipe_tune_exp_coeff(pft)
1128
1129	ENDDO
1130
1131	ELSE
1132
1133	WRITE(numout,*) 'pft ',pft
1134	CALL ipslerr_p (3,'wood_to_cn_eff', &
1135	'wood_to_cn_eff is not defined for this PFT.', &
1136	'See the output file for more details.','')
1137
1138	ENDIF
1139
1140	END FUNCTION wood_to_cn_eff
1141
1142
1143	!! ================================================================================================================================
1144	!! FUNCTION : wood_to_cv
1145	!!
1146	!>\BRIEF Calculate crown volume from woody biomass making use of allometric relationships
1147	!!
1148	!! DESCRIPTION : Calculate basal areadiameter of an individual tree from the woody biomass of that tree making
1149	!! use of allometric relationships:
1150	!! (i) Basal area is written as a function of wood_mass: woodmass_ind = tree_ffpipe_densityba*height
1151	!! (ii) height = pipe_tune2 * sqrt( 4/piba ) * pipe_tune3
1152	!! (iii) cn = pipe_tune1 * sqrt( 4/piba ) * pipe_tune_exp_coeff
1153	!! (iv) cv = 4.0 / 3.0 * pi * ( SQRT( cn/pi ) ) ** 3
1154	!!
1155	!! RECENT CHANGE(S): None
1156	!!
1157	!! RETURN VALUE : diameter (m3 ind-1)
1158	!!
1159	!! REFERENCE(S) :
1160	!!
1161	!! FLOWCHART : None
1162	!! \n
1163	!_ ================================================================================================================================
1164
1165	FUNCTION wood_to_cv(biomass_temp, pft)
1166
1167	!! 0. Variable and parameter declaration
1168
1169	!! 0.1 Input variables
1170
1171	INTEGER(i_std) :: pft !! PFT number (-)
1172	REAL(r_std), DIMENSION(:,:) :: biomass_temp !! Biomass of an individual tree within a circ
1173	!! class @tex $(m^{2} ind^{-1})$ @endtex
1174
1175	!! 0.2 Output variables
1176
1177	REAL(r_std), DIMENSION(ncirc) :: wood_to_cv !! Crown volume of an individual tree
1178	!! @tex $(m^{2} ind{-1})$ @endtex
1179
1180	!! 0.3 Modified variables
1181
1182	!! 0.4 Local variables
1183	INTEGER(i_std) :: l !! Index
1184	REAL(r_std) :: woodmass_ind !! Woodmass of an individual tree
1185	!! @tex $(gC ind^{-1})$ @endtex
1186	REAL(r_std), DIMENSION(ncirc) :: dia !! Diameter of an individual tree (m)
1187
1188	!_ ================================================================================================================================
1189
1190	!! 1. Calculate crown volume from woodmass
1191
1192	IF ( is_tree(pft) ) THEN
1193
1194	! Crown diameter of the individual tree (m2 ind-1)
1195	dia(:) = (4./piwood_to_cn(biomass_temp(:,:),pft))*(1./2.)
1196	wood_to_cv(:) = pi/6.dia(:)*3.
1197	! WRITE(numout,*) 'wood_to_cv, ', wood_to_cv(:)
1198
1199	ELSE
1200
1201	WRITE(numout,*) 'pft ',pft
1202	CALL ipslerr_p (3,'wood_to_cv', &
1203	'wood_to_cv is not defined for this PFT.', &
1204	'See the output file for more details.','')
1205
1206	ENDIF
1207
1208	END FUNCTION wood_to_cv
1209
1210	!! ================================================================================================================================
1211	!! FUNCTION : wood_to_cv_eff
1212	!!
1213	!>\BRIEF Calculate crown volume from woody biomass making use of allometric relationships
1214	!!
1215	!! DESCRIPTION : Calculate basal areadiameter of an individual tree from the woody biomass of that tree making
1216	!! use of allometric relationships:
1217	!! (i) Basal area is written as a function of wood_mass: woodmass_ind = tree_ffpipe_densityba*height
1218	!! (ii) height = pipe_tune2 * sqrt( 4/piba ) * pipe_tune3
1219	!! (iii) cn = pipe_tune1 * sqrt( 4/piba ) * pipe_tune_exp_coeff
1220	!! (iv) cv = 4.0 / 3.0 * pi * ( SQRT( cn/pi ) ) ** 3
1221	!!
1222	!! NOTE: the difference between this and wood_to_cv is that this include aboveground and belowground
1223	!! biomass in the calculation
1224	!!
1225	!! RECENT CHANGE(S): None
1226	!!
1227	!! RETURN VALUE : diameter (m3 ind-1)
1228	!!
1229	!! REFERENCE(S) :
1230	!!
1231	!! FLOWCHART : None
1232	!! \n
1233	!_ ================================================================================================================================
1234
1235	FUNCTION wood_to_cv_eff(biomass_temp, pft)
1236
1237	!! 0. Variable and parameter declaration
1238
1239	!! 0.1 Input variables
1240
1241	INTEGER(i_std) :: pft !! PFT number (-)
1242	REAL(r_std), DIMENSION(:,:) :: biomass_temp !! Biomass of an individual tree within a circ
1243	!! class @tex $(m^{2} ind^{-1})$ @endtex
1244
1245	!! 0.2 Output variables
1246
1247	REAL(r_std), DIMENSION(ncirc) :: wood_to_cv_eff !! Crown volume of an individual tree
1248	!! @tex $(m^{2} ind{-1})$ @endtex
1249
1250	!! 0.3 Modified variables
1251
1252	!! 0.4 Local variables
1253	INTEGER(i_std) :: l !! Index
1254	REAL(r_std) :: woodmass_ind !! Woodmass of an individual tree
1255	!! @tex $(gC ind^{-1})$ @endtex
1256	REAL(r_std), DIMENSION(ncirc) :: dia !! Diameter of an individual tree (m)
1257
1258	!_ ================================================================================================================================
1259
1260	!! 1. Calculate crown volume from woodmass
1261
1262	IF ( is_tree(pft) ) THEN
1263
1264	! Crown diameter of the individual tree (m2 ind-1)
1265	dia(:) = (4./piwood_to_cn_eff(biomass_temp(:,:),pft))*(1./2.)
1266	wood_to_cv_eff(:) = pi/6.dia(:)*3.
1267	! WRITE(numout,*) 'wood_to_cv_eff, ', wood_to_cv_eff(:)
1268
1269	ELSE
1270
1271	WRITE(numout,*) 'pft ',pft
1272	CALL ipslerr_p (3,'wood_to_cv_eff', &
1273	'wood_to_cv_eff is not defined for this PFT.', &
1274	'See the output file for more details.','')
1275
1276	ENDIF
1277
1278	END FUNCTION wood_to_cv_eff
1279
1280
1281
1282	!! ================================================================================================================================
1283	!! FUNCTION : wood_to_volume
1284	!!
1285	!>\BRIEF This allometric function computes volume as a function of
1286	!! biomass at stand scale. Volume \f$(m^3 m^{-2}) = f(biomass (gC m^{-2}))\f$
1287	!!
1288	!! DESCRIPTION : None
1289	!!
1290	!! RECENT CHANGE(S): None
1291	!!
1292	!! RETURN VALUE : biomass_to_volume
1293	!!
1294	!! REFERENCE(S) : See above, module description.
1295	!!
1296	!! FLOWCHART : None
1297	!! \n
1298	!_ ================================================================================================================================
1299
1300	FUNCTION wood_to_volume(biomass,pft,branch_ratio,inc_branches)
1301
1302	!! 0. Variable and parameter declaration
1303
1304	!! 0.1 Input variables
1305
1306	REAL(r_std), DIMENSION(:) :: biomass !! Stand biomass @tex $(gC m^{-2})$ @endtex
1307	REAL(r_std) :: branch_ratio !! Branch ratio of sap and heartwood biomass
1308	!! unitless
1309	INTEGER(i_std) :: pft !! Plant functional type (unitless)
1310	INTEGER(i_std) :: inc_branches !! Include the branches in the volume calculation?
1311	!! 0: exclude the branches from the volume calculation
1312	!! (thus correct the biomass for the branch ratio)
1313	!! 1: include the branches in the volume calculation
1314	!! (thus use all aboveground biomass)
1315
1316
1317
1318	!! 0.2 Output variables
1319
1320	REAL(r_std) :: wood_to_volume !! The volume of wood per square meter
1321	!! @tex $(m^3 m^{-2})$ @endtex
1322
1323	!! 0.3 Modified variables
1324
1325	!! 0.4 Local variables
1326
1327	REAL(r_std) :: woody_biomass !! Woody biomass at the stand level
1328	!! @tex $(gC m^{-2})$ @endtex
1329
1330	!_ ================================================================================================================================
1331
1332	!! 1. Volume to biomass
1333
1334	! Woody biomass used in the calculation
1335	IF (inc_branches .EQ. 0) THEN
1336
1337	woody_biomass=(biomass(isapabove)+biomass(iheartabove))*(un - branch_ratio)
1338
1339	ELSEIF (inc_branches .EQ. 1) THEN
1340
1341	woody_biomass=(biomass(isapabove)+biomass(iheartabove))
1342
1343	ELSE
1344
1345	ENDIF
1346
1347	! Wood volume expressed in m3 / m2
1348	wood_to_volume = woody_biomass/(pipe_density(pft))
1349
1350	END FUNCTION wood_to_volume
1351
1352
1353
1354	!! ================================================================================================================================
1355	!! FUNCTION : biomass_to_lai
1356	!!
1357	!>\BRIEF Calculate the LAI based on the leaf biomass
1358	!!
1359	!! DESCRIPTION : Calculates the LAI of a PFT/grid square based on the leaf biomass
1360	!!
1361	!! RECENT CHANGE(S): None
1362	!!
1363	!! RETURN VALUE : ::LAI [m2 m{-2}]
1364	!!
1365	!! REFERENCE(S) :
1366	!!
1367	!! FLOWCHART : None
1368	!! \n
1369	!_ ================================================================================================================================
1370
1371	FUNCTION biomass_to_lai(leaf_biomass, pft)
1372
1373	!! 0. Variable and parameter declaration
1374
1375	!! 0.1 Input variables
1376
1377	INTEGER(i_std) :: pft !! PFT number (-)
1378	REAL(r_std) :: leaf_biomass !! Biomass of the leaves
1379	!! @tex $(gC m^{-2})$ @endtex
1380
1381
1382	!! 0.2 Output variables
1383
1384	REAL(r_std) :: biomass_to_lai !! Leaf area index
1385	!! @tex $(m^{2} m^{-2})$ @endtex
1386
1387	!! 0.3 Modified variables
1388
1389	!! 0.4 Local variables
1390	REAL(r_std) :: impose_lai !! LAI read from run.def
1391	!_ ================================================================================================================================
1392
1393	!! 1. Calculate the LAI from the leaf biomass
1394
1395	biomass_to_lai = leaf_biomass * sla(pft)
1396
1397	!!$ !+++++++++ TEMP ++++++++++
1398	!!$ ! This is a perfect place to hack the code to make it run with
1399	!!$ ! constant lai
1400	!!$ WRITE(numout,*) 'WARNING ERROR: Using fake lai values for testing!'
1401	!!$ biomass_to_lai=3.79052
1402	!!$ !+++++++++++++++++++++++++
1403
1404	!+++++++ TEMP ++++++++++
1405	! This code is only used evaluation of the performance of the multi-layer energy budget.
1406	! To reduce the complexity of the tests we want to impose the LAI and its vertical distribution.
1407	! The solution is not very elegant but it works.
1408	! IF (ld_fake_height) THEN
1409	! In order to imposed lai, we read the TOTAL_LAI from run.def
1410	! CALL getin_p('TOTAL_LAI', impose_lai)
1411	! This part of code reset the sla vale to match which alow modeled LAI equal to TOTAL LAI.
1412	! Althought this is ugly way to match the modeled LAI and impose LAI.
1413	! You probably need to go to your ORCHIDEE out put file to find out the suitable SLA value
1414	! and reset it agin in the run.def.
1415	! So, we impose LAI & structure for a quick testing the performance of multilayer energy budget
1416	! without changing the leaf_biomass.
1417	! IF ( leaf_biomass .GT. 0.0) THEN
1418	! sla(pft)=impose_lai/leaf_biomass
1419	! WRITE(numout,'(A,F20.8)') 'USE A FAKE SLA BASED ON imposed LAI/LEAFMASS=', sla(pft)
1420	! ENDIF
1421	! biomass_to_lai=leaf_biomass*sla(pft)
1422	! ENDIF
1423	!++++++++++++++++++++++++
1424
1425	END FUNCTION biomass_to_lai
1426
1427
1428	!! ================================================================================================================================
1429	!! FUNCTION : Nmax
1430	!!
1431	!>\BRIEF This function determines the maximum number of trees per hectare for
1432	!! a given quadratic mean diameter (Dg). It applies the self-thinning principle
1433	!! of Reineke (1933), with Dg instead of mean diameter (Dhote 1999).
1434	!! Parameterization: Dhote (1999) for broad-leaved and Vacchiano (2008)
1435	!! for needle-leaved.
1436	!!
1437	!! DESCRIPTION : None
1438	!!
1439	!! RECENT CHANGE(S): None
1440	!!
1441	!! RETURN VALUE : Nmax
1442	!!
1443	!! REFERENCE(S) : See above, module description.
1444	!!
1445	!! FLOWCHART : None
1446	!! \n
1447	!_ ================================================================================================================================
1448
1449	FUNCTION Nmax(Dg,no_pft)
1450
1451	!! 0. Variable and parameter declaration
1452
1453	!! 0.1 Input variables
1454
1455	REAL(r_std) :: Dg !! Quadratic mean diameter (cm)
1456	INTEGER(i_std) :: no_pft !! Plant functional type (unitless)
1457
1458	!! 0.2 Output variables
1459
1460	REAL(r_std) :: Nmax !! Maximum number of trees according to the self-thinning model
1461
1462	!! 0.3 Modified variables
1463
1464	!! 0.4 Local variables
1465
1466	!_ ================================================================================================================================
1467	!! 1. maximum number of trees per hectare for a given quadratic mean diameter
1468
1469	IF (is_tree(no_pft)) THEN
1470
1471	! thinning curve of the MTC
1472	Nmax = (Dg/alpha_self_thinning(no_pft))**(un/beta_self_thinning(no_pft))
1473
1474	! Truncate the range to avoid huge numbers due the exponental model used to describe self-thinning
1475	Nmax = MIN( Nmax, nmaxtrees(no_pft) )
1476	Nmax = MAX( Nmax, dens_target(no_pft) )
1477
1478	ELSE
1479
1480	WRITE(numout,*) 'Self thinning is not defined for PFT, ', no_pft
1481	CALL ipslerr_p (3,'nmax', &
1482	'Self thinning is not defined for this PFT.', &
1483	'See the output file for more details.','')
1484
1485	ENDIF
1486
1487	END FUNCTION Nmax
1488
1489
1490	!! ================================================================================================================================
1491	!! FUNCTION : Nmaxyield
1492	!!
1493	!>\BRIEF This function determines the maximum number of trees per hectare for
1494	!! a given quadratic mean diameter (Dg). It applies the 75 percentile of the European
1495	!! yield tables
1496	!!
1497	!! DESCRIPTION : None
1498	!!
1499	!! RECENT CHANGE(S): None
1500	!!
1501	!! RETURN VALUE : Nmaxyield
1502	!!
1503	!! REFERENCE(S) : See above, module description.
1504	!!
1505	!! FLOWCHART : None
1506	!! \n
1507	!_ ================================================================================================================================
1508
1509	FUNCTION Nmaxyield(Dg,no_pft)
1510
1511	!! 0. Variable and parameter declaration
1512
1513	!! 0.1 Input variables
1514
1515	REAL(r_std) :: Dg !! Quadratic mean diameter (cm)
1516	INTEGER(i_std) :: no_pft !! Plant functional type (unitless)
1517
1518	!! 0.2 Output variables
1519
1520	REAL(r_std) :: Nmaxyield !! Maximum number of trees according to the self-thinning model
1521
1522	!! 0.3 Modified variables
1523
1524	!! 0.4 Local variables
1525
1526	!_ ================================================================================================================================
1527	!! 1. maximum number of trees per hectare for a given quadratic mean diameter
1528
1529	IF (is_tree(no_pft)) THEN
1530
1531	! thinning curve of the MTC
1532	Nmaxyield = (Dg/alpha_rdi_upper(no_pft))**(un/beta_rdi_upper(no_pft))
1533
1534	! Truncate the range to avoid huge numbers due the exponental model used to describe self-thinning
1535	Nmaxyield = MIN( Nmaxyield, nmaxtrees(no_pft) )
1536	Nmaxyield = MAX( Nmaxyield, dens_target(no_pft) )
1537
1538	ELSE
1539
1540	WRITE(numout,*) 'Self thinning is not defined for PFT, ', no_pft
1541	CALL ipslerr_p (3,'Nmaxyield', &
1542	'Self thinning is not defined for this PFT.', &
1543	'See the output file for more details.','')
1544
1545	ENDIF
1546
1547	END FUNCTION Nmaxyield
1548
1549
1550	!! ================================================================================================================================
1551	!! FUNCTION : Nminyield
1552	!!
1553	!>\BRIEF This function determines the minimum number of trees per hectare for
1554	!! a given quadratic mean diameter (Dg). It applies the 25 percentile of the European
1555	!! yield tables
1556	!!
1557	!! DESCRIPTION : None
1558	!!
1559	!! RECENT CHANGE(S): None
1560	!!
1561	!! RETURN VALUE : Nminyield
1562	!!
1563	!! REFERENCE(S) : See above, module description.
1564	!!
1565	!! FLOWCHART : None
1566	!! \n
1567	!_ ================================================================================================================================
1568
1569	FUNCTION Nminyield(Dg,no_pft)
1570
1571	!! 0. Variable and parameter declaration
1572
1573	!! 0.1 Input variables
1574
1575	REAL(r_std) :: Dg !! Quadratic mean diameter (cm)
1576	INTEGER(i_std) :: no_pft !! Plant functional type (unitless)
1577
1578	!! 0.2 Output variables
1579
1580	REAL(r_std) :: Nminyield !! Minimum number of trees according to the self-thinning model
1581
1582	!! 0.3 Modified variables
1583
1584	!! 0.4 Local variables
1585
1586	!_ ================================================================================================================================
1587	!! 1. minimum number of trees per hectare for a given quadratic mean diameter
1588
1589	IF (is_tree(no_pft)) THEN
1590
1591	! thinning curve of the MTC
1592	Nminyield = (Dg/alpha_rdi_lower(no_pft))**(un/beta_rdi_lower(no_pft))
1593
1594	! Truncate the range to avoid huge numbers due the exponental model used to describe self-thinning
1595
1596	Nminyield = MIN( Nminyield, nmaxtrees(no_pft) )
1597	Nminyield = MAX( Nminyield, dens_target(no_pft) )
1598
1599	ELSE
1600
1601	WRITE(numout,*) 'Self thinning is not defined for PFT, ', no_pft
1602	CALL ipslerr_p (3,'Nminyield', &
1603	'Self thinning is not defined for this PFT.', &
1604	'See the output file for more details.','')
1605
1606
1607	ENDIF
1608
1609	END FUNCTION Nminyield
1610
1611
1612	!! ================================================================================================================================
1613	!! SUBROUTINE : distribute_mortality_biomass
1614	!!
1615	!>\BRIEF Distributes biomass that is going to be killed by natural
1616	!! causes (not self thinning) over circ classes.
1617	!!
1618	!! DESCRIPTION : Mortality is going to kill a certain amount of biomass
1619	!! in forests every day. Since we have circumference classes
1620	!! now for our forests, we need to determine which classes
1621	!! of trees will suffer from this environmental mortality.
1622	!! Right now we are taking an exponential distribution.
1623	!! Notice that this is NOT the same as redistributing biomass
1624	!! after one of the circ classes becomes empty.
1625	!!
1626	!! RECENT CHANGE(S): None
1627	!!
1628	!! MAIN OUTPUT VARIABLE(S): ::circ_class_kill
1629	!!
1630	!! REFERENCE(S) : None
1631	!!
1632	!! FLOWCHART : None
1633	!! \n
1634	!_ ================================================================================================================================
1635	SUBROUTINE distribute_mortality_biomass ( bm_difference, ddf_temp, circ_class_n_temp, &
1636	circ_class_biomass_temp, circ_class_kill_temp )
1637
1638	!! 0. Variable and parameter description
1639
1640	!! 0.1 Input variables
1641	REAL(r_std),INTENT(in) :: bm_difference !! the biomass to distribute
1642	REAL(r_std),INTENT(in) :: ddf_temp !! the death_distribution_factor for this pft
1643	REAL(r_std),DIMENSION(:),INTENT(in) :: circ_class_n_temp !! circ_class_n for this point/PFT
1644	REAL(r_std),DIMENSION(:,:),INTENT(in) :: circ_class_biomass_temp !! circ_class_biomass for
1645	!! this point/PFT
1646	!! 0.2 Output variables
1647
1648	!! 0.3 Modified variables
1649	REAL(r_std),DIMENSION(:),INTENT(inout) :: circ_class_kill_temp !! circ_class_kill for
1650	!! this point/PFT/pool
1651
1652	!! 0.4 Local variables
1653	REAL(r_std), DIMENSION(ncirc) :: biomass_desired !! The biomass that dies naturally
1654	REAL(r_std), DIMENSION(ncirc) :: death_distribution !! The fraction of biomass taken from
1655	!! each circ class for mortality.
1656	LOGICAL :: ldone !! Flag to exit a loop
1657	REAL(r_std) :: scale_factor !!
1658	REAL(r_std) :: sum_total !!
1659	REAL(r_std) :: leftover_bm !! excess biomass that we need to kill
1660	REAL(r_std) :: living_biomass !! summed biomass
1661	REAL(r_std) :: living_trees !! summed biomass
1662	INTEGER :: icir
1663
1664	!_ ================================================================================================================================
1665
1666	IF(ncirc == 1)THEN
1667
1668	! This is the easy case. All of our biomass will be taken from the only
1669	! circumference class that we have.
1670	circ_class_kill_temp(1)=circ_class_kill_temp(1)+&
1671	bm_difference/SUM(circ_class_biomass_temp(1,:))
1672
1673	RETURN
1674	ENDIF
1675
1676	! Here we assume an exponential distribution, arranged so that
1677	! ddf_temp times more biomass is taken from the largest circ class
1678	! compared to the smallest.
1679	biomass_desired(:)=zero
1680	death_distribution(:)=un
1681	scale_factor=ddf_temp**(un/REAL(ncirc-1))
1682	DO icir=2,ncirc
1683	death_distribution(icir)=death_distribution(icir-1)*scale_factor
1684	ENDDO
1685	! Normalize it
1686	sum_total=SUM(death_distribution(:))
1687	death_distribution(:)=death_distribution(:)/sum_total
1688
1689	! Ideally, how much is killed from each class? Be careful to include
1690	! what was killed in self-thinning here! If we don't, we may try to kill
1691	! more biomass than is available in the loop below.
1692	DO icir=1,ncirc
1693	biomass_desired(icir)=death_distribution(icir)*bm_difference+&
1694	circ_class_kill_temp(icir)*SUM(circ_class_biomass_temp(icir,:))
1695	ENDDO
1696
1697	! Right now, we know how much biomass we want to kill in each
1698	! class (biomass_desired). What we will do now is to loop through
1699	! all the circ_classes and see if we have this much biomass
1700	! in each class still alive. The total amount of vegetation still
1701	! alive is circ_class_n_temp(icir)-circ_class_kill_temp(icir). If
1702	! this number of individuals cannot give us the total biomass that
1703	! we need, we keep track of a residual quantity, leftover_bm, which
1704	! we try to take from other circ_classes. It's possible that we
1705	! will have to loop several times, which makes things more complicated.
1706
1707	ldone=.FALSE.
1708	leftover_bm=zero
1709	DO
1710	DO icir=ncirc,1,-1
1711
1712	living_trees=circ_class_n_temp(icir)-circ_class_kill_temp(icir)
1713	living_biomass=&
1714	SUM(circ_class_biomass_temp(icir,:))*living_trees
1715	biomass_desired(icir)=biomass_desired(icir)+leftover_bm
1716
1717	IF(living_biomass .LE. biomass_desired(icir))THEN
1718
1719	! We can't get everything from this class, so we kill whatever is left.
1720	leftover_bm=biomass_desired(icir)-living_biomass
1721	biomass_desired(icir)=zero
1722	circ_class_kill_temp(icir)=circ_class_kill_temp(icir)+&
1723	living_trees
1724
1725	ELSE
1726
1727	! We have enough in this class.
1728	circ_class_kill_temp(icir)=circ_class_kill_temp(icir)+&
1729	biomass_desired(icir)/SUM(circ_class_biomass_temp(icir,:))
1730	biomass_desired(icir)=zero
1731	leftover_bm=zero
1732
1733	ENDIF ! living_biomass .LE. biomass_desired(icir)
1734
1735	ENDDO ! loop over circ classes
1736
1737	IF(leftover_bm .LE. min_stomate) EXIT
1738
1739	! it's possible that we don't have enough biomass left to kill what needs to be
1740	! killed, so everything just dies. I cannot think of a case where this
1741	! would happen, though, since mortality should always be a percentage of the
1742	! total biomass.
1743	ldone=.TRUE.
1744	DO icir=1,ncirc
1745
1746	IF( circ_class_kill_temp(icir) .LT. circ_class_n_temp(icir) ) &
1747	ldone=.FALSE.
1748
1749	ENDDO
1750
1751	IF(ldone)EXIT ! All our biomass is dead, and we still want to kill more!
1752
1753	ENDDO
1754
1755
1756
1757	END SUBROUTINE distribute_mortality_biomass
1758
1759
1760	!! ================================================================================================================================
1761	!! SUBROUTINE : check_biomass_sync
1762	!!
1763	!>\BRIEF
1764	!!
1765	!! DESCRIPTION :
1766	!! RECENT CHANGE(S): None
1767	!!
1768	!! MAIN OUTPUT VARIABLE(S):
1769	!!
1770	!! REFERENCE(S) : None
1771	!!
1772	!! FLOWCHART : None
1773	!! \n
1774	!_ ================================================================================================================================
1775	SUBROUTINE check_biomass_sync ( check_point, npts, biomass, &
1776	circ_class_biomass, circ_class_n , ind, &
1777	lsync, bm_sync)
1778
1779	!! 0. Variable and parameter description
1780
1781	!! 0.1 Input variables
1782	INTEGER(i_std), INTENT(in) :: npts !! Domain size (unitless)
1783	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(in) :: circ_class_biomass !! Biomass of the componets of the model
1784	!! tree within a circumference
1785	!! class @tex $(gC ind^{-1})$ @endtex
1786	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: circ_class_n !! Number of individuals in each circ class
1787	!! @tex $(m^{-2})$ @endtex
1788	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: biomass !! Stand level biomass
1789	!! @tex $(gC m^{-2})$ @endtex
1790	CHARACTER(*),INTENT(in) :: check_point !! A flag to indicate at which
1791	!! point in the code we're doing
1792	!! this check
1793	REAL(r_std), DIMENSION(:,:), INTENT(in) :: ind !! Density of individuals
1794	!! @tex $(m^{-2})$ @endtex
1795
1796	!! 0.2 Output variables
1797	LOGICAL,INTENT(out) :: lsync
1798	REAL(r_std), DIMENSION(:,:,:), INTENT(out) :: bm_sync !! The difference betweeen the
1799	!! biomass in the circ_classes and
1800	!! the total biomass
1801	!! @tex $(gC m^{-2})$ @endtex
1802	!! 0.3 Modified variables
1803
1804	!! 0.4 Local variables
1805	INTEGER :: iele,ipts,ivm,ipar,icir
1806	REAL(r_std) :: total_circ_class_biomass
1807	REAL(r_std),DIMENSION(ncirc) :: tree_size
1808	LOGICAL :: lnegative
1809
1810	!_ ================================================================================================================================
1811
1812	lsync=.TRUE.
1813	lnegative=.FALSE.
1814
1815	bm_sync(:,:,:)=zero
1816
1817	!++++++ TEMP ++++++
1818	! We gain 5-10% speed by skipping this routine
1819
1820	!++++++++++++
1821
1822	! Check to see if the biomass is not equal to the total biomass
1823	! in circ_class_biomass anywhere.
1824	DO ipts=1,npts
1825
1826	DO ivm=1,nvm
1827
1828	! Only woody PFTs have circumference classes therefore
1829	! only woody PFTs need to be syncronized
1830	IF(.NOT. lbypass_cc)THEN
1831	IF(is_tree(ivm)) THEN
1832	tree_size(:)=zero
1833	DO icir=1,ncirc
1834	tree_size(icir)=SUM(circ_class_biomass(ipts,ivm,icir,:,1))
1835	ENDDO
1836	DO icir=2,ncirc
1837	IF(tree_size(icir) .LT. tree_size(icir-1)-min_stomate)THEN
1838	WRITE(numout,*) 'ERROR: stopping in sync'
1839	WRITE(numout,*) check_point
1840	WRITE(numout,*) 'ipts,ivm: ',ipts,ivm
1841	WRITE(numout,*) 'tree_size(icir), tree_size(icir-1), ',&
1842	tree_size(icir), tree_size(icir-1), tree_size(icir) - tree_size(icir-1)
1843	WRITE(numout,*) 'icir, tree_size: ',icir, tree_size(:)
1844	!+++ TEMP +++
1845	!This would not STOP the ORCHIDEE beacause the mass balance is due to imposed LAI
1846	IF(ld_fake_height) THEN
1847	CALL ipslerr_p (2,'check_biomass_sync', &
1848	'The size of the trees in the circ class are not monotonically increasing!',&
1849	'Look in the output file for more details.',&
1850	'')
1851	ELSE
1852	CALL ipslerr_p (3,'check_biomass_sync', &
1853	'The size of the trees in the circ class are not monotonically increasing!',&
1854	'Look in the output file for more details.',&
1855	'')
1856	ENDIF
1857	!++++++++++++
1858	ENDIF
1859	ENDDO
1860	ENDIF
1861	ENDIF
1862
1863	DO iele=1,icarbon
1864
1865	DO ipar=1,nparts
1866
1867	total_circ_class_biomass=zero
1868	DO icir=1,ncirc
1869
1870	total_circ_class_biomass=total_circ_class_biomass+&
1871	circ_class_biomass(ipts,ivm,icir,ipar,iele)*circ_class_n(ipts,ivm,icir)
1872
1873	! Check as well to see if our biomass is ever negative.
1874	! It really should not be.
1875	IF(circ_class_biomass(ipts,ivm,icir,ipar,iele) .LT. -min_stomate)THEN
1876
1877	lnegative=.TRUE.
1878	WRITE(numout,) '!**********************************'
1879	WRITE(numout,*) 'Error: Negative biomass component!'
1880	WRITE(numout,*) 'Check point: ',TRIM(check_point)
1881	WRITE(numout,*) 'circ_class_biomass(ipts,ivm,icir,ipar,iele) ',&
1882	circ_class_biomass(ipts,ivm,icir,ipar,iele)
1883	WRITE(numout,'(A,5I5)') 'ipts,ivm,icir,ipar,iele',ipts,ivm,icir,ipar,iele
1884	WRITE(numout,) '!**********************************'
1885
1886	ENDIF
1887	ENDDO
1888
1889	IF(ABS(biomass(ipts,ivm,ipar,iele) - &
1890	total_circ_class_biomass) .GT. sync_threshold)THEN
1891
1892	WRITE(numout,) '!**********************************'
1893	WRITE(numout,*) 'Biomass and circ_class_biomass are not equal!'
1894	WRITE(numout,*) 'Check point: ',TRIM(check_point)
1895	WRITE(numout,100) 'biomass(ipts,ivm,ipar,iele) ',&
1896	biomass(ipts,ivm,ipar,iele)
1897	WRITE(numout,100) 'total_circ_class_biomass ',&
1898	total_circ_class_biomass
1899	WRITE(numout,100) 'Difference: ',&
1900	ABS(biomass(ipts,ivm,ipar,iele) - total_circ_class_biomass)
1901	WRITE(numout,*) 'ipts,ivm,ipar,iele',ipts,ivm,ipar,iele
1902	WRITE(numout,) '!**********************************'
1903	100 FORMAT(A,E20.10)
1904	lsync=.FALSE.
1905
1906	ENDIF
1907
1908	ENDDO
1909
1910	! we are not going to save the biomass for every component right now,
1911	! just the total
1912	bm_sync(ipts,ivm,iele)=zero
1913
1914	DO ipar=1,nparts
1915
1916
1917	DO icir=1,ncirc
1918
1919	bm_sync(ipts,ivm,iele)=bm_sync(ipts,ivm,iele)+&
1920	circ_class_biomass(ipts,ivm,icir,ipar,iele)*circ_class_n(ipts,ivm,icir)
1921	ENDDO ! ncirc
1922
1923	ENDDO ! nparts
1924
1925	bm_sync(ipts,ivm,iele)=ABS(bm_sync(ipts,ivm,iele)-&
1926	SUM(biomass(ipts,ivm,:,iele)))
1927
1928	ENDDO ! nelements
1929
1930
1931
1932	ENDDO ! loop over PFTs
1933
1934	ENDDO ! loop over points
1935
1936	!---TEMP---
1937	IF(ld_biomass)THEN
1938	WRITE(numout,*) 'Check point: ',TRIM(check_point)
1939	WRITE(numout,*) 'test_pft, test_grid: ',test_pft,test_grid
1940	WRITE(numout,*) 'biomass (ileaf), ', biomass(test_grid,test_pft,ileaf,icarbon)
1941	WRITE(numout,*) 'biomass (iwood), ', biomass(test_grid,test_pft,isapabove,icarbon) + &
1942	biomass(test_grid,test_pft,isapbelow,icarbon) + biomass(test_grid,test_pft,iheartabove,icarbon) + &
1943	biomass(test_grid,test_pft,iheartbelow,icarbon)
1944	WRITE(numout,*) 'biomass (iroot), ', biomass(test_grid,test_pft,iroot,icarbon)
1945	WRITE(numout,'(A,20F14.6)') 'biomassHHH, ',biomass(test_grid,test_pft,:,icarbon)
1946	DO icir=1,ncirc
1947	WRITE(numout,'(A,I1,20F14.6)') 'ccbiomass',icir,circ_class_biomass(test_grid,test_pft,icir,:,icarbon)
1948	ENDDO
1949	WRITE(numout,*) 'circ_class_biomass, ',&
1950	SUM (SUM(circ_class_biomass(test_grid,test_pft,:,:,icarbon),2) * &
1951	circ_class_n(test_grid,test_pft,:))
1952	WRITE(numout,*) 'circ_class_n, ', SUM(circ_class_n(test_grid,test_pft,:))
1953	WRITE(numout,*) 'circ_class_n(:), ', circ_class_n(test_grid,test_pft,:)
1954	WRITE(numout,*) 'ind, ', ind(test_grid,test_pft)
1955	ENDIF
1956
1957	!!$ !----------
1958
1959	IF(.NOT. lsync) THEN
1960	WRITE(numout,*) 'ERROR: stopping in sync #2'
1961	WRITE(numout,*) 'Stopping'
1962	CALL ipslerr_p (3,'check_biomass_sync', &
1963	'circ_class_biomass*circ_class_n is not equal to the total biomass',&
1964	'Look in the output file for more details.',&
1965	'')
1966
1967	ENDIF
1968	IF(lnegative) THEN
1969	WRITE(numout,*) 'ERROR: negative biomass'
1970	WRITE(numout,*) 'Stopping'
1971	CALL ipslerr_p (3,'check_biomass_sync', &
1972	'One of the biomass pools is negative!',&
1973	'Look in the output file for more details.',&
1974	'')
1975	ENDIF
1976
1977	END SUBROUTINE check_biomass_sync
1978
1979	!! ================================================================================================================================
1980	!! SUBROUTINE : check_area
1981	!!
1982	!>\BRIEF: Check whether surface area is preserved after forest
1983	!! management, age_class_distribution and land use change
1984	!!
1985	!! DESCRIPTION : After processes where biomass pools could become zero, e.g., forest management
1986	!! a mass abalnce check alone is not good enough because the check is based on
1987	!! biomass*veget_max. Biomass is zero so veget_max could take any value. Here it is
1988	!! explicitly checked whether veget_max is preserved.
1989	!!
1990	!! RECENT CHANGE(S): None
1991	!!
1992	!! MAIN OUTPUT VARIABLE(S): None
1993	!!
1994	!! REFERENCE(S) : None
1995	!!
1996	!! FLOWCHART : None
1997	!! \n
1998	!_ ================================================================================================================================
1999	SUBROUTINE check_area (check_point, npts, veget_max_begin, veget_max)
2000
2001	!! 0. Variable and parameter description
2002
2003	!! 0.1 Input variables
2004	CHARACTER(*),INTENT(in) :: check_point !! A flag to indicate at which point in the code we're doing
2005	!! this check
2006	INTEGER(i_std), INTENT(in) :: npts !! Number of pixels
2007	REAL(r_std), INTENT(in), DIMENSION(:,:) :: veget_max_begin !! Veget_max at the start of the routine that needs to be checked
2008	REAL(r_std), INTENT(in), DIMENSION(:,:) :: veget_max !! Veget_max at the moment this routine is checked
2009
2010	!! 0.2 Output variables
2011
2012	!! 0.3 Modified variables
2013
2014	!! 0.4 Local variables
2015	INTEGER(i_std) :: ipts,ivm !! Indices
2016	REAL(r_std) :: area_begin !! Temporairy variable to check area conservation
2017	REAL(r_std) :: area_end !! Temporairy variable to check area conservation
2018
2019	!_ ================================================================================================================================
2020
2021	! Which subroutine is being checked?
2022	WRITE(numout,*) 'Check point: ',TRIM(check_point)
2023
2024	! Check for surface area presevation
2025	DO ipts = 1,npts
2026
2027	DO ivm = 1,nvm
2028
2029	! Area balance closure
2030	IF (ivm .EQ. 1) THEN
2031
2032	! Initialize
2033	area_begin = veget_max_begin(ipts,ivm)
2034	area_end = veget_max(ipts,ivm)
2035
2036	ELSEIF (agec_group(ivm-1) .EQ. agec_group(ivm)) THEN
2037
2038	! The current PFT belongs to the same species
2039	! as the previous. Cumulate the areas.
2040	area_begin = area_begin + veget_max_begin(ipts,ivm)
2041	area_end = area_end + veget_max(ipts,ivm)
2042
2043	ELSEIF (agec_group(ivm-1) .NE. agec_group(ivm)) THEN
2044
2045	! The current PFT belongs to a different species
2046	! as the previous. Check the areas of the previous group
2047	! and reset
2048	IF (ABS(area_begin-area_begin) .LE. min_stomate) THEN
2049
2050	! The code preserves the surface area
2051	IF (ld_massbal) WRITE(numout,*) &
2052	'Preserved veget_max in ', check_point, ', igroup, ', &
2053	agec_group(ivm-1)
2054
2055	ELSE
2056
2057	! Surface area is NOT preserved. Write some
2058	! additional information to diagnose the problem
2059	WRITE(numout,*) 'Error: veget_max is not closed in ', check_point
2060	WRITE(numout,*) ' ipts,igroup; ', ipts,agec_group(ivm-1)
2061	WRITE(numout,*) ' Difference is, ', area_end-area_begin
2062	WRITE(numout,*) ' veget_max_begin,veget_max_end: ', area_begin, &
2063	area_end
2064	IF(ld_stop)THEN
2065	CALL ipslerr_p (3,check_point, 'Area error','','')
2066	ENDIF
2067
2068	ENDIF
2069
2070	! Write the surface area to the next
2071	! species group
2072	area_begin = veget_max_begin(ipts,ivm)
2073	area_end = veget_max(ipts,ivm)
2074
2075	ELSE
2076
2077	! Shouldn't be here
2078	CALL ipslerr_p (3,check_point, &
2079	'unexpected case in if-loop','','')
2080	ENDIF
2081
2082	! last pft
2083	IF (ivm .EQ. nvm) THEN
2084
2085	! The code preserves the surface area
2086	IF (ABS(area_begin-area_begin) .LE. min_stomate) THEN
2087
2088	IF (ld_massbal) WRITE(numout,*) &
2089	'Preserved veget_max in ',check_point, ', igroup, ', &
2090	agec_group(ivm)
2091
2092	ELSE
2093
2094	! Surface area is NOT preserved. Write some
2095	! additional information to diagnose the problem
2096	WRITE(numout,*) 'Error: veget_max is not closed in ',check_point
2097	WRITE(numout,*) ' ipts,igroup; ', ipts,agec_group(ivm)
2098	WRITE(numout,*) ' Difference is, ', area_end-area_begin
2099	WRITE(numout,*) ' veget_max_begin,veget_max_end: ', area_begin, &
2100	area_end
2101	IF(ld_stop)THEN
2102	CALL ipslerr_p (3,check_point, 'Area error','','')
2103	ENDIF
2104
2105	ENDIF ! preserve area
2106	ENDIF ! last ivm
2107
2108	ENDDO ! nvm
2109	ENDDO ! npts
2110
2111	END SUBROUTINE check_area
2112
2113
2114	END MODULE function_library

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

Context Navigation

source: branches/publications/ORCHIDEE_CAN_r3069/src_parameters/function_library.f90 @ 7346

Download in other formats: