Context Navigation

source: branches/publications/ORCHIDEE_CAN_r2290/src_parameters/function_library.f90 @ 7844

Last change on this file since 7844 was 2275, checked in by sebastiaan.luyssaert, 10 years ago
DEV: works on a single pixel. Experimenting with adaptation to waterstress. Committed to migrate the code to Curie. Do not use for regional simulations. Updates will follow.
File size: 72.7 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
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	STOP
245
246	ENDIF
247
248	END FUNCTION wood_to_ba_eff_array
249
250
251	!! ================================================================================================================================
252	!! FUNCTION : wood_to_ba_eff
253	!!
254	!>\BRIEF Calculate effective basal area from woody biomass making use of allometric relationships
255	!!
256	!! DESCRIPTION : Calculate basal area of an individual tree from the woody biomass of that tree making
257	!! use of allometric relationships. Effective basal area accounts for both above and below ground carbon
258	!! and is the basis for the application of the rule of Deleuze and Dhote.
259	!! (i) woodmass = tree_ff * pipe_densitybaheight
260	!! (ii) height = pipe_tune2 * sqrt(4/piba) * pipe_tune_3
261	!!
262	!! RECENT CHANGE(S): None
263	!!
264	!! RETURN VALUE : effective basal area (m2 ind-1)
265	!!
266	!! REFERENCE(S) :
267	!!
268	!! FLOWCHART : None
269	!! \n
270	!_ ================================================================================================================================
271
272	FUNCTION wood_to_ba_eff(biomass_temp, pft)
273
274	!! 0. Variable and parameter declaration
275
276	!! 0.1 Input variables
277
278	INTEGER(i_std) :: pft !! PFT number (-)
279	REAL(r_std), DIMENSION(:,:) :: biomass_temp !! Biomass of an individual tree within a circ
280	!! class @tex $(m^{2} ind^{-1})$ @endtex
281
282	!! 0.2 Output variables
283
284	REAL(r_std), DIMENSION(ncirc) :: wood_to_ba_eff !! Effective basal area of an individual tree within a circ
285	!! class @tex $(m^{2} ind^{-1})$ @endtex
286
287	!! 0.3 Modified variables
288
289	!! 0.4 Local variables
290	INTEGER(i_std) :: l !! Index
291	REAL(r_std) :: woodmass_ind !! Woodmass of an individual tree
292	!! @tex $(gC ind^{-1})$ @endtex
293
294	!_ ================================================================================================================================
295
296	!! 1. Calculate basal area from woodmass
297
298	IF ( is_tree(pft) ) THEN
299
300	DO l = 1,ncirc
301
302	! Woodmass of an individual tree
303	woodmass_ind = biomass_temp(l,isapabove) + biomass_temp(l,isapbelow) + &
304	biomass_temp(l,iheartabove) + biomass_temp(l,iheartbelow)
305
306	! Basal area of that individual (m2 ind-1)
307	wood_to_ba_eff(l) = (pi/4(woodmass_ind/(tree_ff(pft)pipe_density(pft)*pipe_tune2(pft))) &
308	(2./pipe_tune3(pft)))(pipe_tune3(pft)/(pipe_tune3(pft)+2))
309
310	ENDDO
311
312	ELSE
313
314	WRITE(numout,*) 'pft ',pft
315	STOP 'function wood_to_ba_eff is not defined for this PFT'
316
317	ENDIF
318
319	END FUNCTION wood_to_ba_eff
320
321
322
323	!! ================================================================================================================================
324	!! FUNCTION : wood_to_ba
325	!!
326	!>\BRIEF Calculate basal area from woody biomass making use of allometric relationships
327	!!
328	!! DESCRIPTION : Calculate basal area of an individual tree from the woody biomass of that tree making
329	!! use of allometric relationships given below. Here basal area is defined in line with its classical
330	!! forestry meaning.
331	!! (i) woodmass = tree_ff * pipe_densitybaheight
332	!! (ii) height = pipe_tune2 * sqrt(4/piba) * pipe_tune_3
333	!!
334	!! RECENT CHANGE(S): None
335	!!
336	!! RETURN VALUE : basal area (m2 ind-1)
337	!!
338	!! REFERENCE(S) :
339	!!
340	!! FLOWCHART : None
341	!! \n
342	!_ ================================================================================================================================
343
344	FUNCTION wood_to_ba(biomass_temp, pft)
345
346	!! 0. Variable and parameter declaration
347
348	!! 0.1 Input variables
349
350	INTEGER(i_std) :: pft !! PFT number (-)
351	REAL(r_std), DIMENSION(:,:) :: biomass_temp !! Biomass of an individual tree within a circ
352	!! class @tex $(m^{2} ind^{-1})$ @endtex
353
354	!! 0.2 Output variables
355
356	REAL(r_std), DIMENSION(ncirc) :: wood_to_ba !! Basal area of an individual tree within a circ
357	!! class @tex $(m^{2} ind^{-1})$ @endtex
358
359	!! 0.3 Modified variables
360
361	!! 0.4 Local variables
362	INTEGER(i_std) :: l !! Index
363	REAL(r_std) :: woodmass_ind !! Woodmass of an individual tree
364	!! @tex $(gC ind^{-1})$ @endtex
365
366	!_ ================================================================================================================================
367
368	!! 1. Calculate basal area from woodmass
369
370	IF ( is_tree(pft) ) THEN
371
372	DO l = 1,ncirc
373
374	! Woodmass of an individual tree
375	woodmass_ind = biomass_temp(l,iheartabove) + biomass_temp(l,isapabove)
376
377
378	! Basal area of that individual (m2 ind-1)
379	wood_to_ba(l) = (pi/4(woodmass_ind/(tree_ff(pft)pipe_density(pft)*pipe_tune2(pft))) &
380	(2./pipe_tune3(pft)))(pipe_tune3(pft)/(pipe_tune3(pft)+2))
381
382	ENDDO
383
384	ELSE
385
386	WRITE(numout,*) 'pft ',pft
387	STOP 'function wood_to_ba is not defined for this PFT'
388
389	ENDIF
390
391	END FUNCTION wood_to_ba
392
393
394	!! ================================================================================================================================
395	!! FUNCTION : wood_to_height_eff
396	!!
397	!>\BRIEF Calculate the effective tree height from woody biomass making use of allometric relationships
398	!!
399	!! DESCRIPTION : Calculate the effective height of an individual tree from the woody biomass of that tree making
400	!! use of allometric relationships. Effective height makes use of both above and belowground biomass and is
401	!! used in the calculation of the allocation according to deleuze and dhote, hydraulic architecture because also
402	!! the height of the belowground part should be included.
403	!! (i) height(:) = pipe_tune2(j)(4/piba(:))**(pipe_tune3(j)/2) and
404	!! (ii) woodmass_ind = tree_ffpipe_densityba*height
405	!!
406	!! RECENT CHANGE(S): None
407	!!
408	!! RETURN VALUE : height (m)
409	!!
410	!! REFERENCE(S) :
411	!!
412	!! FLOWCHART : None
413	!! \n
414	!_ ================================================================================================================================
415
416	FUNCTION wood_to_height_eff(biomass_temp, pft)
417
418	!! 0. Variable and parameter declaration
419
420	!! 0.1 Input variables
421
422	INTEGER(i_std) :: pft !! PFT number (-)
423	REAL(r_std), DIMENSION(:,:) :: biomass_temp !! Biomass of an individual tree within a circ
424	!! class @tex $(m^{2} ind^{-1})$ @endtex
425
426	!! 0.2 Output variables
427
428	REAL(r_std), DIMENSION(ncirc) :: wood_to_height_eff !! Effective height of an individual tree within a circ
429	!! class (m)
430
431	!! 0.3 Modified variables
432
433	!! 0.4 Local variables
434	INTEGER(i_std) :: l !! Index
435	REAL(r_std) :: woodmass_ind !! Woodmass of an individual tree
436	!! @tex $(gC ind{-1})$ @endtex
437
438	!_ ================================================================================================================================
439
440	!! 1. Calculate height from woodmass
441
442	IF ( is_tree(pft) ) THEN
443
444	DO l = 1,ncirc
445
446	! Woodmass of an individual tree. Both above and belowground biomass are used
447	woodmass_ind = biomass_temp(l,isapabove) + biomass_temp(l,isapbelow) + &
448	biomass_temp(l,iheartabove) + biomass_temp(l,iheartbelow)
449
450	! Height of that individual
451	wood_to_height_eff(l) = (((woodmass_ind/(tree_ff(pft)pipe_density(pft))4/pi)**(pipe_tune3(pft)/2))&
452	pipe_tune2(pft))*(1/(pipe_tune3(pft)/2+1))
453
454	ENDDO
455
456	ELSE
457
458	WRITE(numout,*) 'pft ',pft
459	STOP 'function wood_to_height_eff is not defined for this PFT'
460
461	ENDIF
462
463	END FUNCTION wood_to_height_eff
464
465
466	!! ================================================================================================================================
467	!! FUNCTION : wood_to_height
468	!!
469	!>\BRIEF Calculate tree height from woody biomass making use of allometric relationships
470	!!
471	!! DESCRIPTION : Calculate height of an individual tree from the woody biomass of that tree making
472	!! use of allometric relationships. This is the height used in forestry and for calculating the aerodynamic
473	!! interactions
474	!! (i) height(:) = pipe_tune2(j)(4/piba(:))**(pipe_tune3(j)/2) and
475	!! (ii) woodmass_ind = tree_ffpipe_densityba*height
476	!!
477	!! RECENT CHANGE(S): None
478	!!
479	!! RETURN VALUE : height (m)
480	!!
481	!! REFERENCE(S) :
482	!!
483	!! FLOWCHART : None
484	!! \n
485	!_ ================================================================================================================================
486
487	FUNCTION wood_to_height(biomass_temp, pft)
488
489	!! 0. Variable and parameter declaration
490
491	!! 0.1 Input variables
492
493	INTEGER(i_std) :: pft !! PFT number (-)
494	REAL(r_std), DIMENSION(:,:) :: biomass_temp !! Biomass of an individual tree within a circ
495	!! class @tex $(m^{2} ind^{-1})$ @endtex
496
497	!! 0.2 Output variables
498
499	REAL(r_std), DIMENSION(ncirc) :: wood_to_height !! Height of an individual tree within a circ
500	!! class (m)
501
502	!! 0.3 Modified variables
503
504	!! 0.4 Local variables
505	INTEGER(i_std) :: l !! Index
506	REAL(r_std) :: woodmass_ind !! Woodmass of an individual tree
507	!! @tex $(gC ind{-1})$ @endtex
508
509	!_ ================================================================================================================================
510
511	!! 1. Calculate height from woodmass
512
513	IF ( is_tree(pft) ) THEN
514
515	DO l = 1,ncirc
516
517	! Woodmass of an individual tree (only the aboveground component)
518	woodmass_ind = biomass_temp(l,isapabove) + biomass_temp(l,iheartabove)
519
520	! Height of that individual
521	wood_to_height(l) = (((woodmass_ind/(tree_ff(pft)pipe_density(pft))4/pi)**(pipe_tune3(pft)/2))&
522	pipe_tune2(pft))*(1/(pipe_tune3(pft)/2+1))
523
524	ENDDO
525
526	ELSE
527
528	WRITE(numout,*) 'pft ',pft
529	STOP 'function wood_to_height is not defined for this PFT'
530
531	ENDIF
532
533	END FUNCTION wood_to_height
534
535
536	!! ================================================================================================================================
537	!! FUNCTION : wood_to_qmheight
538	!!
539	!>\BRIEF Calculate the quadratic mean height from the biomass
540	!!
541	!! DESCRIPTION : Calculates the quadratic mean height from the biomass
542	!!
543	!! RECENT CHANGE(S): None
544	!!
545	!! RETURN VALUE : ::qm_height (m)
546	!!
547	!! REFERENCE(S) :
548	!!
549	!! FLOWCHART : None
550	!! \n
551	!_ ================================================================================================================================
552
553	FUNCTION wood_to_qmheight(biomass_temp, ind, pft)
554
555	!! 0. Variable and parameter declaration
556
557	!! 0.1 Input variables
558
559	INTEGER(i_std) :: pft !! PFT number (-)
560	REAL(r_std), DIMENSION(ncirc,nparts) :: biomass_temp !! Biomass of the leaves @tex $(gC m^{-2})$ @endtex
561	REAL(r_std), DIMENSION(ncirc) :: ind !! Number of individuals @tex $(m^{-2})$ @endtex
562
563
564	!! 0.2 Output variables
565
566	REAL(r_std) :: wood_to_qmheight !! quadratic mean height (m)
567
568	!! 0.3 Modified variables
569
570	!! 0.4 Local variables
571	REAL(r_std), DIMENSION(ncirc) :: circ_class_ba !! basal area for each circ_class @tex $(m^{2})$ @endtex
572	REAL(r_std) :: qm_dia !! quadratic mean diameter (m)
573
574	!_ ================================================================================================================================
575
576	!! 1. Calculate qm_height from the biomass
577	IF ( is_tree(pft) ) THEN
578
579	! Basal area at the tree level (m2 tree-1)
580	circ_class_ba(:) = wood_to_ba(biomass_temp(:,:),pft)
581
582	IF (SUM(ind(:)) .NE. zero) THEN
583
584	qm_dia = SQRT( 4/piSUM(circ_class_ba(:)ind(:))/SUM(ind(:)) )
585
586	ELSE
587
588	qm_dia = zero
589
590	ENDIF
591
592	wood_to_qmheight = pipe_tune2(pft)(qm_dia*pipe_tune3(pft))
593
594
595	! Grasses and croplands
596	ELSE
597
598	! Calculate height as a function of the leaf and structural biomass. Use structural
599	! biomass to make sure that the grasslands have a roughness length during the winter
600	! If the biomass increases, vegetation height will increase as well. Divide by
601	! ind(ipts,j) to obtain the height of the individual. biomass(ileaf) is in gC m-2
602	! whereas qm is the height of the individual.
603	IF (SUM(ind(:)) .NE. zero) THEN
604
605	wood_to_qmheight = SUM(biomass_temp(:,ileaf) + biomass_temp(:,isapabove)) / &
606	SUM(ind(:)) * sla(pft) * lai_to_height(pft)
607
608	ELSE
609
610	wood_to_qmheight = zero
611
612	ENDIF
613
614	ENDIF ! is_tree(j)
615
616	END FUNCTION wood_to_qmheight
617
618
619
620	!! ================================================================================================================================
621	!! FUNCTION : wood_to_dia_eff
622	!!
623	!>\BRIEF Calculate effective diameter from woody biomass making use of allometric relationships
624	!!
625	!! DESCRIPTION : Calculate the effective diameter of an individual tree from the woody biomass of that tree making
626	!! use of allometric relationships. Effective diameter accounts for both above and belowground biomass.
627	!! (i) woodmass_ind = tree_ff * pipe_density * height * pi/4dia*2
628	!! (ii) height = pipe_tune2 * dia * pipe_tune3
629	!!
630	!! RECENT CHANGE(S): None
631	!!
632	!! RETURN VALUE : diameter (m)
633	!!
634	!! REFERENCE(S) :
635	!!
636	!! FLOWCHART : None
637	!! \n
638	!_ ================================================================================================================================
639
640	FUNCTION wood_to_dia_eff(biomass_temp, pft)
641
642	!! 0. Variable and parameter declaration
643
644	!! 0.1 Input variables
645
646	INTEGER(i_std) :: pft !! PFT number (-)
647	REAL(r_std), DIMENSION(:,:) :: biomass_temp !! Biomass of an individual tree within a circ
648	!! class @tex $(m^{2} ind^{-1})$ @endtex
649
650	!! 0.2 Output variables
651
652	REAL(r_std), DIMENSION(ncirc) :: wood_to_dia_eff !! Diameter of an individual tree within a circ
653	!! class (m)
654
655	!! 0.3 Modified variables
656
657	!! 0.4 Local variables
658	INTEGER(i_std) :: l !! Index
659	REAL(r_std) :: woodmass_ind !! Woodmass of an individual tree
660	!! @tex $(gC ind^{-1})$ @endtex
661
662	!_ ================================================================================================================================
663
664	!! 1. Calculate basal area from woodmass
665
666	IF ( is_tree(pft) ) THEN
667
668	DO l = 1,ncirc
669
670	! Woodmass of an individual tree
671	woodmass_ind = biomass_temp(l,isapabove) + biomass_temp(l,isapbelow) + &
672	biomass_temp(l,iheartabove) + biomass_temp(l,iheartbelow)
673
674	! Basal area of that individual (m2 ind-1)
675	wood_to_dia_eff(l) = (4/piwoodmass_ind/(tree_ff(pft)pipe_density(pft)pipe_tune2(pft))) * &
676	(1./(2+pipe_tune3(pft)))
677
678	ENDDO
679
680	ELSE
681
682	WRITE(numout,*) 'pft ',pft
683	STOP 'function wood_to_dia_eff is not defined for this PFT'
684
685	ENDIF
686
687	END FUNCTION wood_to_dia_eff
688
689
690
691	!! ================================================================================================================================
692	!! FUNCTION : wood_to_dia
693	!!
694	!>\BRIEF Calculate diameter from woody biomass making use of allometric relationships
695	!!
696	!! DESCRIPTION : Calculate diameter of an individual tree from the woody biomass of that tree making
697	!! use of allometric relationships. Makes only use of the aboveground biomass and relates to the
698	!! typical forestry diameter (but not normalized at 1.3 m)
699	!! (i) woodmass_ind = tree_ff * pipe_density * height * pi/4dia*2
700	!! (ii) height = pipe_tune2 * dia * pipe_tune3
701	!!
702	!! RECENT CHANGE(S): None
703	!!
704	!! RETURN VALUE : diameter (m)
705	!!
706	!! REFERENCE(S) :
707	!!
708	!! FLOWCHART : None
709	!! \n
710	!_ ================================================================================================================================
711
712	FUNCTION wood_to_dia(biomass_temp, pft)
713
714	!! 0. Variable and parameter declaration
715
716	!! 0.1 Input variables
717
718	INTEGER(i_std) :: pft !! PFT number (-)
719	REAL(r_std), DIMENSION(:,:) :: biomass_temp !! Biomass of an individual tree within a circ
720	!! class @tex $(m^{2} ind^{-1})$ @endtex
721
722	!! 0.2 Output variables
723
724	REAL(r_std), DIMENSION(ncirc) :: wood_to_dia !! Diameter of an individual tree within a circ
725	!! class (m)
726
727	!! 0.3 Modified variables
728
729	!! 0.4 Local variables
730	INTEGER(i_std) :: l !! Index
731	REAL(r_std) :: woodmass_ind !! Woodmass of an individual tree
732	!! @tex $(gC ind^{-1})$ @endtex
733
734	!_ ================================================================================================================================
735
736	!! 1. Calculate basal area from woodmass
737
738	IF ( is_tree(pft) ) THEN
739
740	DO l = 1,ncirc
741
742	! Woodmass of an individual tree
743	woodmass_ind = biomass_temp(l,isapabove) + biomass_temp(l,iheartabove)
744
745	! Basal area of that individual (m2 ind-1)
746	wood_to_dia(l) = (4/piwoodmass_ind/(tree_ff(pft)pipe_density(pft)pipe_tune2(pft))) * &
747	(1./(2+pipe_tune3(pft)))
748
749	ENDDO
750
751	ELSE
752
753	WRITE(numout,*) 'pft ',pft
754	STOP 'function wood_to_dia is not defined for this PFT'
755
756	ENDIF
757
758	END FUNCTION wood_to_dia
759
760
761	!! ================================================================================================================================
762	!! FUNCTION : wood_to_qmdia
763	!!
764	!>\BRIEF Calculate the quadratic mean diameter from the biomass
765	!!
766	!! DESCRIPTION : Calculates the quadratic mean diameter from the aboveground biomss
767	!!
768	!! RECENT CHANGE(S): None
769	!!
770	!! RETURN VALUE : ::qm_dia (m)
771	!!
772	!! REFERENCE(S) :
773	!!
774	!! FLOWCHART : None
775	!! \n
776	!_ ================================================================================================================================
777
778	FUNCTION wood_to_qmdia(biomass_temp, ind, pft)
779
780	!! 0. Variable and parameter declaration
781
782	!! 0.1 Input variables
783
784	INTEGER(i_std) :: pft !! PFT number (-)
785	REAL(r_std), DIMENSION(ncirc,nparts) :: biomass_temp !! Biomass of the leaves @tex $(gC m^{-2})$ @endtex
786	REAL(r_std), DIMENSION(ncirc) :: ind !! Number of individuals @tex $(m^{-2})$ @endtex
787
788	!! 0.2 Output variables
789
790	REAL(r_std) :: wood_to_qmdia !! quadratic mean diameter (m)
791
792	!! 0.3 Modified variables
793
794	!! 0.4 Local variables
795	REAL(r_std), DIMENSION(ncirc) :: circ_class_ba !! basal area for each circ_class @tex $(m^{2})$ @endtex
796
797	!_ ================================================================================================================================
798
799	!! 1. Calculate qm_dia from the biomass
800	IF ( is_tree(pft) ) THEN
801
802	! Basal area at the tree level (m2 tree-1)
803	circ_class_ba(:) = wood_to_ba(biomass_temp(:,:),pft)
804
805	IF (SUM(ind(:)) .NE. zero) THEN
806
807	wood_to_qmdia = SQRT( 4/piSUM(circ_class_ba(:)ind(:))/SUM(ind(:)) )
808
809	ELSE
810
811	wood_to_qmdia = zero
812
813	ENDIF
814
815
816	! Grasses and croplands
817	ELSE
818
819	wood_to_qmdia = zero
820
821	ENDIF ! is_tree(pft)
822
823	END FUNCTION wood_to_qmdia
824
825
826	!! ================================================================================================================================
827	!! FUNCTION : wood_to_circ
828	!!
829	!>\BRIEF Calculate circumference from woody biomass making use of allometric relationships
830	!!
831	!! DESCRIPTION : All this does it computer the diameter using a different routine, and then
832	!! convert that into a circumference.
833	!!
834	!! RECENT CHANGE(S): None
835	!!
836	!! RETURN VALUE : circumference (m)
837	!!
838	!! REFERENCE(S) :
839	!!
840	!! FLOWCHART : None
841	!! \n
842	!_ ================================================================================================================================
843
844	FUNCTION wood_to_circ(biomass_temp, pft)
845
846	!! 0. Variable and parameter declaration
847
848	!! 0.1 Input variables
849
850	INTEGER(i_std) :: pft !! PFT number (-)
851	REAL(r_std), DIMENSION(:,:) :: biomass_temp !! Biomass of an individual tree within a circ
852	!! class @tex $(m^{2} ind^{-1})$ @endtex
853
854	!! 0.2 Output variables
855
856	REAL(r_std), DIMENSION(ncirc) :: wood_to_circ !! Circumference of an individual tree within a circ
857	!! class (m)
858
859	!! 0.3 Modified variables
860
861	!! 0.4 Local variables
862
863	!_ ================================================================================================================================
864
865	!! 1. Calculate diameter from woodmass
866
867	wood_to_circ(:)=val_exp
868
869	wood_to_circ(:)=wood_to_dia(biomass_temp(:,:),pft)
870
871	! convert to a circumference (m)
872	wood_to_circ(:) = wood_to_circ(:)*pi
873
874	END FUNCTION wood_to_circ
875
876
877
878	!! ================================================================================================================================
879	!! FUNCTION : wood_to_cn_array
880	!!
881	!>\BRIEF Calculate crown area from woody biomass making use of allometric relationships
882	!!
883	!! DESCRIPTION : Calculate crown area of an individual tree from the woody biomass of that tree
884	!! making use of allometric relationship that relates crown area (cn) to diameter (dia) as
885	!! pipe_tune1dia*pipe_tune_exp_coeff where the pipe_tune parameters are pft-specific.
886	!! (i) Basal area is written as a function of wood_mass: woodmass_ind = tree_ffpipe_densityba*height
887	!! (ii) height = pipe_tune2sqrt(4/piba)**pipe_tune3
888	!! (iii) cn = pipe_tune1sqrt(4/piba)**pipe_tune_exp_coeff
889	!!
890	!! RECENT CHANGE(S): None
891	!!
892	!! RETURN VALUE : crown area (m2 ind-1)
893	!!
894	!! REFERENCE(S) :
895	!!
896	!! FLOWCHART : None
897	!! \n
898	!_ ================================================================================================================================
899
900	FUNCTION wood_to_cn_array(biomass_temp, npts, pft)
901
902	!! 0. Variable and parameter declaration
903
904	!! 0.1 Input variables
905
906	INTEGER(i_std) :: pft !! PFT number (-)
907	INTEGER(i_std) :: npts !! Pixel(s), this variable defines the dimensions of ba
908	REAL(r_std), DIMENSION(:,:,:) :: biomass_temp !! Biomass of an individual tree within a circ
909	!! class @tex $(m^{2} ind^{-1})$ @endtex
910
911
912	!! 0.2 Output variables
913
914	REAL(r_std), DIMENSION(npts,ncirc) :: wood_to_cn_array !! Crown area of an individual tree @tex $(m^{2} ind{-1})$ @endtex
915
916	!! 0.3 Modified variables
917
918	!! 0.4 Local variables
919
920	INTEGER(i_std) :: l !! index
921	REAL(r_std), DIMENSION(npts) :: woodmass_ind !! Woodmass of an individual tree @tex $(gC ind{-1})$ @endtex
922
923	!_ ================================================================================================================================
924
925	!! 1. Calculate crown area from basal area
926
927	IF ( is_tree(pft) ) THEN
928
929	DO l = 1,ncirc
930
931	! Woodmass of an individual tree
932	woodmass_ind(:) = biomass_temp(:,l,isapabove) + biomass_temp(:,l,iheartabove)
933
934	! Crown area of that individual
935	! cn = pipe_tune1sqrt(4/piba)**pipe_tune_exp_coeff
936	wood_to_cn_array(:,l) = pipe_tune1(pft) * SQRT( 4/pi(pi/4(woodmass_ind(:)/&
937	(tree_ff(pft)pipe_density(pft)pipe_tune2(pft)))(2./pipe_tune3(pft)))&
938	(pipe_tune3(pft)/(pipe_tune3(pft)+2)) ) ** pipe_tune_exp_coeff(pft)
939
940	ENDDO
941
942	ELSE
943
944	WRITE(numout,*) 'pft ',pft
945	STOP 'function wood_to_cn_array is not defined for this PFT'
946
947	ENDIF
948
949	END FUNCTION wood_to_cn_array
950
951
952	!! ================================================================================================================================
953	!! FUNCTION : wood_to_cn
954	!!
955	!>\BRIEF Calculate crown area from woody biomass making use of allometric relationships
956	!!
957	!! DESCRIPTION : Calculate crown area of an individual tree from the woody biomass of that tree
958	!! making use of allometric relationship that relates crown area (cn) to diameter (dia) as
959	!! pipe_tune1dia*pipe_tune_exp_coeff where the pipe_tune parameters are pft-specific.
960	!! (i) Basal area is written as a function of wood_mass: woodmass_ind = tree_ffpipe_densityba*height
961	!! (ii) height = pipe_tune2sqrt(4/piba)**pipe_tune3
962	!! (iii) cn = pipe_tune1sqrt(4/piba)**pipe_tune_exp_coeff
963	!!
964	!! RECENT CHANGE(S): None
965	!!
966	!! RETURN VALUE : crown area (m2 ind-1)
967	!!
968	!! REFERENCE(S) :
969	!!
970	!! FLOWCHART : None
971	!! \n
972	!_ ================================================================================================================================
973
974	FUNCTION wood_to_cn(biomass_temp, pft)
975
976	!! 0. Variable and parameter declaration
977
978	!! 0.1 Input variables
979
980	INTEGER(i_std) :: pft !! PFT number (-)
981	REAL(r_std), DIMENSION(:,:) :: biomass_temp !! biomass of an individual tree within a circ
982	!! class @tex $(m^{2} ind^{-1})$ @endtex
983
984	!! 0.2 Output variables
985
986	REAL(r_std), DIMENSION(ncirc) :: wood_to_cn !! Crown area of an individual tree within a circ
987	!! class @tex $(m^{2} ind-1)$ @endtex
988
989	!! 0.3 Modified variables
990
991	!! 0.4 Local variables
992	INTEGER(i_std) :: l !! Index
993	REAL(r_std) :: woodmass_ind !! Woodmass of an individual tree
994	!! @tex $(gC ind{-1})$ @endtex
995
996	!_ ================================================================================================================================
997
998	!! 1. Calculate crown area from woodmass
999
1000	IF ( is_tree(pft) ) THEN
1001
1002	DO l = 1,ncirc
1003
1004	! Woodmass of an individual tree
1005	woodmass_ind = biomass_temp(l,isapabove) + biomass_temp(l,iheartabove)
1006
1007	! Crown area of that individual
1008	! cn = pipe_tune1sqrt(4/piba)**pipe_tune_exp_coeff
1009	wood_to_cn(l) = pipe_tune1(pft) * SQRT( 4/pi(pi/4(woodmass_ind/&
1010	(tree_ff(pft)pipe_density(pft)pipe_tune2(pft)))(2./pipe_tune3(pft)))&
1011	(pipe_tune3(pft)&
1012	/(pipe_tune3(pft)+2)) ) ** pipe_tune_exp_coeff(pft)
1013
1014	ENDDO
1015
1016	ELSE
1017
1018	WRITE(numout,*) 'pft ',pft
1019	STOP 'function wood_to_cn is not defined for this PFT'
1020
1021	ENDIF
1022
1023	END FUNCTION wood_to_cn
1024
1025	!! ================================================================================================================================
1026	!! FUNCTION : wood_to_cn_eff
1027	!!
1028	!>\BRIEF Calculate crown area from woody biomass making use of allometric relationships
1029	!!
1030	!! DESCRIPTION : Calculate crown area of an individual tree from the woody biomass of that tree
1031	!! making use of allometric relationship that relates crown area (cn) to diameter (dia) as
1032	!! pipe_tune1dia*pipe_tune_exp_coeff where the pipe_tune parameters are pft-specific.
1033	!! (i) Basal area is written as a function of wood_mass: woodmass_ind = tree_ffpipe_densityba*height
1034	!! (ii) height = pipe_tune2sqrt(4/piba)**pipe_tune3
1035	!! (iii) cn = pipe_tune1sqrt(4/piba)**pipe_tune_exp_coeff
1036	!!
1037	!! NOTE: This is different from wood_to_cn because this uses both the aboveground and
1038	!! belowground wood mass.
1039	!!
1040	!! RECENT CHANGE(S): None
1041	!!
1042	!! RETURN VALUE : crown area (m2 ind-1)
1043	!!
1044	!! REFERENCE(S) :
1045	!!
1046	!! FLOWCHART : None
1047	!! \n
1048	!_ ================================================================================================================================
1049
1050	FUNCTION wood_to_cn_eff(biomass_temp, pft)
1051
1052	!! 0. Variable and parameter declaration
1053
1054	!! 0.1 Input variables
1055
1056	INTEGER(i_std) :: pft !! PFT number (-)
1057	REAL(r_std), DIMENSION(:,:) :: biomass_temp !! biomass of an individual tree within a circ
1058	!! class @tex $(m^{2} ind^{-1})$ @endtex
1059
1060	!! 0.2 Output variables
1061
1062	REAL(r_std), DIMENSION(ncirc) :: wood_to_cn_eff !! Crown area of an individual tree within a circ
1063	!! class @tex $(m^{2} ind-1)$ @endtex
1064
1065	!! 0.3 Modified variables
1066
1067	!! 0.4 Local variables
1068	INTEGER(i_std) :: l !! Index
1069	REAL(r_std) :: woodmass_ind !! Woodmass of an individual tree
1070	!! @tex $(gC ind{-1})$ @endtex
1071
1072	!_ ================================================================================================================================
1073
1074	!! 1. Calculate crown area from woodmass
1075
1076	IF ( is_tree(pft) ) THEN
1077
1078	DO l = 1,ncirc
1079
1080	! Woodmass of an individual tree
1081	woodmass_ind = biomass_temp(l,isapabove) + biomass_temp(l,iheartabove) + &
1082	biomass_temp(l,isapbelow) + biomass_temp(l,iheartbelow)
1083
1084	! Crown area of that individual
1085	! cn = pipe_tune1sqrt(4/piba)**pipe_tune_exp_coeff
1086	wood_to_cn_eff(l) = pipe_tune1(pft) * SQRT( 4/pi(pi/4(woodmass_ind/&
1087	(tree_ff(pft)pipe_density(pft)pipe_tune2(pft)))(2./pipe_tune3(pft)))&
1088	(pipe_tune3(pft)&
1089	/(pipe_tune3(pft)+2)) ) ** pipe_tune_exp_coeff(pft)
1090
1091	ENDDO
1092
1093	ELSE
1094
1095	WRITE(numout,*) 'pft ',pft
1096	STOP 'function wood_to_cn_eff is not defined for this PFT'
1097
1098	ENDIF
1099
1100	END FUNCTION wood_to_cn_eff
1101
1102
1103	!! ================================================================================================================================
1104	!! FUNCTION : wood_to_cv
1105	!!
1106	!>\BRIEF Calculate crown volume from woody biomass making use of allometric relationships
1107	!!
1108	!! DESCRIPTION : Calculate basal areadiameter of an individual tree from the woody biomass of that tree making
1109	!! use of allometric relationships:
1110	!! (i) Basal area is written as a function of wood_mass: woodmass_ind = tree_ffpipe_densityba*height
1111	!! (ii) height = pipe_tune2 * sqrt( 4/piba ) * pipe_tune3
1112	!! (iii) cn = pipe_tune1 * sqrt( 4/piba ) * pipe_tune_exp_coeff
1113	!! (iv) cv = 4.0 / 3.0 * pi * ( SQRT( cn/pi ) ) ** 3
1114	!!
1115	!! RECENT CHANGE(S): None
1116	!!
1117	!! RETURN VALUE : diameter (m3 ind-1)
1118	!!
1119	!! REFERENCE(S) :
1120	!!
1121	!! FLOWCHART : None
1122	!! \n
1123	!_ ================================================================================================================================
1124
1125	FUNCTION wood_to_cv(biomass_temp, pft)
1126
1127	!! 0. Variable and parameter declaration
1128
1129	!! 0.1 Input variables
1130
1131	INTEGER(i_std) :: pft !! PFT number (-)
1132	REAL(r_std), DIMENSION(:,:) :: biomass_temp !! Biomass of an individual tree within a circ
1133	!! class @tex $(m^{2} ind^{-1})$ @endtex
1134
1135	!! 0.2 Output variables
1136
1137	REAL(r_std), DIMENSION(ncirc) :: wood_to_cv !! Crown volume of an individual tree
1138	!! @tex $(m^{2} ind{-1})$ @endtex
1139
1140	!! 0.3 Modified variables
1141
1142	!! 0.4 Local variables
1143	INTEGER(i_std) :: l !! Index
1144	REAL(r_std) :: woodmass_ind !! Woodmass of an individual tree
1145	!! @tex $(gC ind^{-1})$ @endtex
1146	REAL(r_std), DIMENSION(ncirc) :: dia !! Diameter of an individual tree (m)
1147
1148	!_ ================================================================================================================================
1149
1150	!! 1. Calculate crown volume from woodmass
1151
1152	IF ( is_tree(pft) ) THEN
1153
1154	! Crown diameter of the individual tree (m2 ind-1)
1155	dia(:) = (4./piwood_to_cn(biomass_temp(:,:),pft))*(1./2.)
1156	wood_to_cv(:) = pi/6.dia(:)*3.
1157	! WRITE(numout,*) 'wood_to_cv, ', wood_to_cv(:)
1158
1159	ELSE
1160
1161	WRITE(numout,*) 'pft ',pft
1162	STOP 'function wood_to_cv is not defined for this PFT'
1163
1164	ENDIF
1165
1166	END FUNCTION wood_to_cv
1167
1168	!! ================================================================================================================================
1169	!! FUNCTION : wood_to_cv_eff
1170	!!
1171	!>\BRIEF Calculate crown volume from woody biomass making use of allometric relationships
1172	!!
1173	!! DESCRIPTION : Calculate basal areadiameter of an individual tree from the woody biomass of that tree making
1174	!! use of allometric relationships:
1175	!! (i) Basal area is written as a function of wood_mass: woodmass_ind = tree_ffpipe_densityba*height
1176	!! (ii) height = pipe_tune2 * sqrt( 4/piba ) * pipe_tune3
1177	!! (iii) cn = pipe_tune1 * sqrt( 4/piba ) * pipe_tune_exp_coeff
1178	!! (iv) cv = 4.0 / 3.0 * pi * ( SQRT( cn/pi ) ) ** 3
1179	!!
1180	!! NOTE: the difference between this and wood_to_cv is that this include aboveground and belowground
1181	!! biomass in the calculation
1182	!!
1183	!! RECENT CHANGE(S): None
1184	!!
1185	!! RETURN VALUE : diameter (m3 ind-1)
1186	!!
1187	!! REFERENCE(S) :
1188	!!
1189	!! FLOWCHART : None
1190	!! \n
1191	!_ ================================================================================================================================
1192
1193	FUNCTION wood_to_cv_eff(biomass_temp, pft)
1194
1195	!! 0. Variable and parameter declaration
1196
1197	!! 0.1 Input variables
1198
1199	INTEGER(i_std) :: pft !! PFT number (-)
1200	REAL(r_std), DIMENSION(:,:) :: biomass_temp !! Biomass of an individual tree within a circ
1201	!! class @tex $(m^{2} ind^{-1})$ @endtex
1202
1203	!! 0.2 Output variables
1204
1205	REAL(r_std), DIMENSION(ncirc) :: wood_to_cv_eff !! Crown volume of an individual tree
1206	!! @tex $(m^{2} ind{-1})$ @endtex
1207
1208	!! 0.3 Modified variables
1209
1210	!! 0.4 Local variables
1211	INTEGER(i_std) :: l !! Index
1212	REAL(r_std) :: woodmass_ind !! Woodmass of an individual tree
1213	!! @tex $(gC ind^{-1})$ @endtex
1214	REAL(r_std), DIMENSION(ncirc) :: dia !! Diameter of an individual tree (m)
1215
1216	!_ ================================================================================================================================
1217
1218	!! 1. Calculate crown volume from woodmass
1219
1220	IF ( is_tree(pft) ) THEN
1221
1222	! Crown diameter of the individual tree (m2 ind-1)
1223	dia(:) = (4./piwood_to_cn_eff(biomass_temp(:,:),pft))*(1./2.)
1224	wood_to_cv_eff(:) = pi/6.dia(:)*3.
1225	! WRITE(numout,*) 'wood_to_cv_eff, ', wood_to_cv_eff(:)
1226
1227	ELSE
1228
1229	WRITE(numout,*) 'pft ',pft
1230	STOP 'function wood_to_cv_eff is not defined for this PFT'
1231
1232	ENDIF
1233
1234	END FUNCTION wood_to_cv_eff
1235
1236
1237
1238	!! ================================================================================================================================
1239	!! FUNCTION : wood_to_volume
1240	!!
1241	!>\BRIEF This allometric function computes volume as a function of
1242	!! biomass at stand scale. Volume \f$(m^3 m^{-2}) = f(biomass (gC m^{-2}))\f$
1243	!!
1244	!! DESCRIPTION : None
1245	!!
1246	!! RECENT CHANGE(S): None
1247	!!
1248	!! RETURN VALUE : biomass_to_volume
1249	!!
1250	!! REFERENCE(S) : See above, module description.
1251	!!
1252	!! FLOWCHART : None
1253	!! \n
1254	!_ ================================================================================================================================
1255
1256	FUNCTION wood_to_volume(biomass,pft,branch_ratio,inc_branches)
1257
1258	!! 0. Variable and parameter declaration
1259
1260	!! 0.1 Input variables
1261
1262	REAL(r_std), DIMENSION(:) :: biomass !! Stand biomass @tex $(gC m^{-2})$ @endtex
1263	REAL(r_std) :: branch_ratio !! Branch ratio of sap and heartwood biomass
1264	!! unitless
1265	INTEGER(i_std) :: pft !! Plant functional type (unitless)
1266	INTEGER(i_std) :: inc_branches !! Include the branches in the volume calculation?
1267	!! 0: exclude the branches from the volume calculation
1268	!! (thus correct the biomass for the branch ratio)
1269	!! 1: include the branches in the volume calculation
1270	!! (thus use all aboveground biomass)
1271
1272
1273
1274	!! 0.2 Output variables
1275
1276	REAL(r_std) :: wood_to_volume !! The volume of wood per square meter
1277	!! @tex $(m^3 m^{-2})$ @endtex
1278
1279	!! 0.3 Modified variables
1280
1281	!! 0.4 Local variables
1282
1283	REAL(r_std) :: woody_biomass !! Woody biomass at the stand level
1284	!! @tex $(gC m^{-2})$ @endtex
1285
1286	!_ ================================================================================================================================
1287
1288	!! 1. Volume to biomass
1289
1290	! Woody biomass used in the calculation
1291	IF (inc_branches .EQ. 0) THEN
1292
1293	woody_biomass=(biomass(isapabove)+biomass(iheartabove))*(un - branch_ratio)
1294
1295	ELSEIF (inc_branches .EQ. 1) THEN
1296
1297	woody_biomass=(biomass(isapabove)+biomass(iheartabove))
1298
1299	ELSE
1300
1301	ENDIF
1302
1303	! Wood volume expressed in m3 / m2
1304	wood_to_volume = woody_biomass/(pipe_density(pft))
1305
1306	END FUNCTION wood_to_volume
1307
1308
1309
1310	!! ================================================================================================================================
1311	!! FUNCTION : biomass_to_lai
1312	!!
1313	!>\BRIEF Calculate the LAI based on the leaf biomass
1314	!!
1315	!! DESCRIPTION : Calculates the LAI of a PFT/grid square based on the leaf biomass
1316	!!
1317	!! RECENT CHANGE(S): None
1318	!!
1319	!! RETURN VALUE : ::LAI [m2 m{-2}]
1320	!!
1321	!! REFERENCE(S) :
1322	!!
1323	!! FLOWCHART : None
1324	!! \n
1325	!_ ================================================================================================================================
1326
1327	FUNCTION biomass_to_lai(leaf_biomass, pft)
1328
1329	!! 0. Variable and parameter declaration
1330
1331	!! 0.1 Input variables
1332
1333	INTEGER(i_std) :: pft !! PFT number (-)
1334	REAL(r_std) :: leaf_biomass !! Biomass of the leaves
1335	!! @tex $(gC m^{-2})$ @endtex
1336
1337
1338	!! 0.2 Output variables
1339
1340	REAL(r_std) :: biomass_to_lai !! Leaf area index
1341	!! @tex $(m^{2} m^{-2})$ @endtex
1342
1343	!! 0.3 Modified variables
1344
1345	!! 0.4 Local variables
1346
1347	!_ ================================================================================================================================
1348
1349	!! 1. Calculate the LAI from the leaf biomass
1350
1351	biomass_to_lai = leaf_biomass * sla(pft)
1352
1353	!!$ !+++++++++ TEMP ++++++++++
1354	!!$ ! This is a perfect place to hack the code to make it run with
1355	!!$ ! constant lai
1356	!!$ WRITE(numout,*) 'WARNING ERROR: Using fake lai values for testing!'
1357	!!$ biomass_to_lai=3.79052
1358	!!$ !+++++++++++++++++++++++++
1359
1360	END FUNCTION biomass_to_lai
1361
1362
1363
1364
1365
1366
1367
1368
1369
1370
1371
1372
1373
1374
1375
1376
1377
1378
1379
1380
1381	!! ================================================================================================================================
1382	!! FUNCTION : Nmax
1383	!!
1384	!>\BRIEF This function determines the maximum number of trees per hectare for
1385	!! a given quadratic mean diameter (Dg). It applies the self-thinning principle
1386	!! of Reineke (1933), with Dg instead of mean diameter (Dhote 1999).
1387	!! Parameterization: Dhote (1999) for broad-leaved and Vacchiano (2008)
1388	!! for needle-leaved.
1389	!!
1390	!! DESCRIPTION : None
1391	!!
1392	!! RECENT CHANGE(S): None
1393	!!
1394	!! RETURN VALUE : Nmax
1395	!!
1396	!! REFERENCE(S) : See above, module description.
1397	!!
1398	!! FLOWCHART : None
1399	!! \n
1400	!_ ================================================================================================================================
1401
1402	FUNCTION Nmax(Dg,no_pft)
1403
1404	!! 0. Variable and parameter declaration
1405
1406	!! 0.1 Input variables
1407
1408	REAL(r_std) :: Dg !! Quadratic mean diameter (cm)
1409	INTEGER(i_std) :: no_pft !! Plant functional type (unitless)
1410
1411	!! 0.2 Output variables
1412
1413	REAL(r_std) :: Nmax !! Maximum number of trees according to the self-thinning model
1414
1415	!! 0.3 Modified variables
1416
1417	!! 0.4 Local variables
1418
1419	!_ ================================================================================================================================
1420	!! 1. maximum number of trees per hectare for a given quadratic mean diameter
1421
1422	IF (is_tree(no_pft)) THEN
1423
1424	! thinning curve of the MTC
1425	Nmax = (Dg/alpha_self_thinning(no_pft))**(un/beta_self_thinning(no_pft))
1426
1427	! Truncate the range to avoid huge numbers due the exponental model used to describe self-thinning
1428	Nmax = MIN( Nmax, nmaxtrees(no_pft) )
1429	Nmax = MAX( Nmax, dens_target(no_pft) )
1430
1431	ELSE
1432
1433	WRITE(numout,*) 'Self thinning is not defined for PFT, ', no_pft
1434	STOP
1435
1436	ENDIF
1437
1438	END FUNCTION Nmax
1439
1440
1441	!! ================================================================================================================================
1442	!! FUNCTION : Nmaxyield
1443	!!
1444	!>\BRIEF This function determines the maximum number of trees per hectare for
1445	!! a given quadratic mean diameter (Dg). It applies the 75 percentile of the European
1446	!! yield tables
1447	!!
1448	!! DESCRIPTION : None
1449	!!
1450	!! RECENT CHANGE(S): None
1451	!!
1452	!! RETURN VALUE : Nmaxyield
1453	!!
1454	!! REFERENCE(S) : See above, module description.
1455	!!
1456	!! FLOWCHART : None
1457	!! \n
1458	!_ ================================================================================================================================
1459
1460	FUNCTION Nmaxyield(Dg,no_pft)
1461
1462	!! 0. Variable and parameter declaration
1463
1464	!! 0.1 Input variables
1465
1466	REAL(r_std) :: Dg !! Quadratic mean diameter (cm)
1467	INTEGER(i_std) :: no_pft !! Plant functional type (unitless)
1468
1469	!! 0.2 Output variables
1470
1471	REAL(r_std) :: Nmaxyield !! Maximum number of trees according to the self-thinning model
1472
1473	!! 0.3 Modified variables
1474
1475	!! 0.4 Local variables
1476
1477	!_ ================================================================================================================================
1478	!! 1. maximum number of trees per hectare for a given quadratic mean diameter
1479
1480	IF (is_tree(no_pft)) THEN
1481
1482	! thinning curve of the MTC
1483	Nmaxyield = (Dg/alpha_rdi_upper(no_pft))**(un/beta_rdi_upper(no_pft))
1484
1485	! Truncate the range to avoid huge numbers due the exponental model used to describe self-thinning
1486	Nmaxyield = MIN( Nmaxyield, nmaxtrees(no_pft) )
1487	Nmaxyield = MAX( Nmaxyield, dens_target(no_pft) )
1488
1489	ELSE
1490
1491	WRITE(numout,*) 'Self thinning is not defined for PFT, ', no_pft
1492	STOP
1493
1494	ENDIF
1495
1496	END FUNCTION Nmaxyield
1497
1498
1499	!! ================================================================================================================================
1500	!! FUNCTION : Nminyield
1501	!!
1502	!>\BRIEF This function determines the minimum number of trees per hectare for
1503	!! a given quadratic mean diameter (Dg). It applies the 25 percentile of the European
1504	!! yield tables
1505	!!
1506	!! DESCRIPTION : None
1507	!!
1508	!! RECENT CHANGE(S): None
1509	!!
1510	!! RETURN VALUE : Nminyield
1511	!!
1512	!! REFERENCE(S) : See above, module description.
1513	!!
1514	!! FLOWCHART : None
1515	!! \n
1516	!_ ================================================================================================================================
1517
1518	FUNCTION Nminyield(Dg,no_pft)
1519
1520	!! 0. Variable and parameter declaration
1521
1522	!! 0.1 Input variables
1523
1524	REAL(r_std) :: Dg !! Quadratic mean diameter (cm)
1525	INTEGER(i_std) :: no_pft !! Plant functional type (unitless)
1526
1527	!! 0.2 Output variables
1528
1529	REAL(r_std) :: Nminyield !! Minimum number of trees according to the self-thinning model
1530
1531	!! 0.3 Modified variables
1532
1533	!! 0.4 Local variables
1534
1535	!_ ================================================================================================================================
1536	!! 1. minimum number of trees per hectare for a given quadratic mean diameter
1537
1538	IF (is_tree(no_pft)) THEN
1539
1540	! thinning curve of the MTC
1541	Nminyield = (Dg/alpha_rdi_lower(no_pft))**(un/beta_rdi_lower(no_pft))
1542
1543	! Truncate the range to avoid huge numbers due the exponental model used to describe self-thinning
1544
1545	Nminyield = MIN( Nminyield, nmaxtrees(no_pft) )
1546	Nminyield = MAX( Nminyield, dens_target(no_pft) )
1547
1548	ELSE
1549
1550	WRITE(numout,*) 'Self thinning is not defined for PFT, ', no_pft
1551	STOP
1552
1553	ENDIF
1554
1555	END FUNCTION Nminyield
1556
1557
1558	!! ================================================================================================================================
1559	!! SUBROUTINE : distribute_mortality_biomass
1560	!!
1561	!>\BRIEF Distributes biomass that is going to be killed by natural
1562	!! causes (not self thinning) over circ classes.
1563	!!
1564	!! DESCRIPTION : Mortality is going to kill a certain amount of biomass
1565	!! in forests every day. Since we have circumference classes
1566	!! now for our forests, we need to determine which classes
1567	!! of trees will suffer from this environmental mortality.
1568	!! Right now we are taking an exponential distribution.
1569	!! Notice that this is NOT the same as redistributing biomass
1570	!! after one of the circ classes becomes empty.
1571	!!
1572	!! RECENT CHANGE(S): None
1573	!!
1574	!! MAIN OUTPUT VARIABLE(S): ::circ_class_kill
1575	!!
1576	!! REFERENCE(S) : None
1577	!!
1578	!! FLOWCHART : None
1579	!! \n
1580	!_ ================================================================================================================================
1581	SUBROUTINE distribute_mortality_biomass ( bm_difference, ddf_temp, circ_class_n_temp, &
1582	circ_class_biomass_temp, circ_class_kill_temp )
1583
1584	!! 0. Variable and parameter description
1585
1586	!! 0.1 Input variables
1587	REAL(r_std),INTENT(in) :: bm_difference !! the biomass to distribute
1588	REAL(r_std),INTENT(in) :: ddf_temp !! the death_distribution_factor for this pft
1589	REAL(r_std),DIMENSION(:),INTENT(in) :: circ_class_n_temp !! circ_class_n for this point/PFT
1590	REAL(r_std),DIMENSION(:,:),INTENT(in) :: circ_class_biomass_temp !! circ_class_biomass for
1591	!! this point/PFT
1592
1593	!! 0.2 Output variables
1594
1595
1596	!! 0.3 Modified variables
1597	REAL(r_std),DIMENSION(:),INTENT(inout) :: circ_class_kill_temp !! circ_class_kill for
1598	!! this point/PFT/pool
1599
1600
1601	!! 0.4 Local variables
1602	REAL(r_std), DIMENSION(ncirc) :: biomass_desired !! The biomass that dies naturally
1603	REAL(r_std), DIMENSION(ncirc) :: death_distribution !! The fraction of biomass taken from
1604	!! each circ class for mortality.
1605	LOGICAL :: ldone !! Flag to exit a loop
1606	REAL(r_std) :: scale_factor !!
1607	REAL(r_std) :: sum_total !!
1608	REAL(r_std) :: leftover_bm !! excess biomass that we need to kill
1609	REAL(r_std) :: living_biomass !! summed biomass
1610	REAL(r_std) :: living_trees !! summed biomass
1611	INTEGER :: icir
1612
1613	!_ ================================================================================================================================
1614
1615	IF(ncirc == 1)THEN
1616
1617	! This is the easy case. All of our biomass will be taken from the only
1618	! circumference class that we have.
1619	circ_class_kill_temp(1)=circ_class_kill_temp(1)+&
1620	bm_difference/SUM(circ_class_biomass_temp(1,:))
1621
1622	RETURN
1623	ENDIF
1624
1625	! Here we assume an exponential distribution, arranged so that
1626	! ddf_temp times more biomass is taken from the largest circ class
1627	! compared to the smallest.
1628	biomass_desired(:)=zero
1629	death_distribution(:)=un
1630	scale_factor=ddf_temp**(un/REAL(ncirc-1))
1631	DO icir=2,ncirc
1632	death_distribution(icir)=death_distribution(icir-1)*scale_factor
1633	ENDDO
1634	! Normalize it
1635	sum_total=SUM(death_distribution(:))
1636	death_distribution(:)=death_distribution(:)/sum_total
1637
1638	! Ideally, how much is killed from each class? Be careful to include
1639	! what was killed in self-thinning here! If we don't, we may try to kill
1640	! more biomass than is available in the loop below.
1641	DO icir=1,ncirc
1642	biomass_desired(icir)=death_distribution(icir)*bm_difference+&
1643	circ_class_kill_temp(icir)*SUM(circ_class_biomass_temp(icir,:))
1644	ENDDO
1645
1646	! Right now, we know how much biomass we want to kill in each
1647	! class (biomass_desired). What we will do now is to loop through
1648	! all the circ_classes and see if we have this much biomass
1649	! in each class still alive. The total amount of vegetation still
1650	! alive is circ_class_n_temp(icir)-circ_class_kill_temp(icir). If
1651	! this number of individuals cannot give us the total biomass that
1652	! we need, we keep track of a residual quantity, leftover_bm, which
1653	! we try to take from other circ_classes. It's possible that we
1654	! will have to loop several times, which makes things more complicated.
1655
1656	ldone=.FALSE.
1657	leftover_bm=zero
1658	DO
1659	DO icir=ncirc,1,-1
1660
1661	living_trees=circ_class_n_temp(icir)-circ_class_kill_temp(icir)
1662	living_biomass=&
1663	SUM(circ_class_biomass_temp(icir,:))*living_trees
1664	biomass_desired(icir)=biomass_desired(icir)+leftover_bm
1665
1666	IF(living_biomass .LE. biomass_desired(icir))THEN
1667
1668	! We can't get everything from this class, so we kill whatever is left.
1669	leftover_bm=biomass_desired(icir)-living_biomass
1670	biomass_desired(icir)=zero
1671	circ_class_kill_temp(icir)=circ_class_kill_temp(icir)+&
1672	living_trees
1673
1674	ELSE
1675
1676	! We have enough in this class.
1677	circ_class_kill_temp(icir)=circ_class_kill_temp(icir)+&
1678	biomass_desired(icir)/SUM(circ_class_biomass_temp(icir,:))
1679	biomass_desired(icir)=zero
1680	leftover_bm=zero
1681
1682	ENDIF ! living_biomass .LE. biomass_desired(icir)
1683
1684	ENDDO ! loop over circ classes
1685
1686	IF(leftover_bm .LE. min_stomate) EXIT
1687
1688	! it's possible that we don't have enough biomass left to kill what needs to be
1689	! killed, so everything just dies. I cannot think of a case where this
1690	! would happen, though, since mortality should always be a percentage of the
1691	! total biomass.
1692	ldone=.TRUE.
1693	DO icir=1,ncirc
1694
1695	IF( circ_class_kill_temp(icir) .LT. circ_class_n_temp(icir) ) &
1696	ldone=.FALSE.
1697
1698	ENDDO
1699
1700	IF(ldone)EXIT ! All our biomass is dead, and we still want to kill more!
1701
1702	ENDDO
1703
1704
1705
1706	END SUBROUTINE distribute_mortality_biomass
1707
1708
1709	!! ================================================================================================================================
1710	!! SUBROUTINE : check_biomass_sync
1711	!!
1712	!>\BRIEF
1713	!!
1714	!! DESCRIPTION :
1715	!! RECENT CHANGE(S): None
1716	!!
1717	!! MAIN OUTPUT VARIABLE(S):
1718	!!
1719	!! REFERENCE(S) : None
1720	!!
1721	!! FLOWCHART : None
1722	!! \n
1723	!_ ================================================================================================================================
1724	SUBROUTINE check_biomass_sync ( check_point, npts, biomass, &
1725	circ_class_biomass, circ_class_n , ind, &
1726	lsync, bm_sync)
1727
1728	!! 0. Variable and parameter description
1729
1730	!! 0.1 Input variables
1731	INTEGER(i_std), INTENT(in) :: npts !! Domain size (unitless)
1732	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(in) :: circ_class_biomass !! Biomass of the componets of the model
1733	!! tree within a circumference
1734	!! class @tex $(gC ind^{-1})$ @endtex
1735	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: circ_class_n !! Number of individuals in each circ class
1736	!! @tex $(m^{-2})$ @endtex
1737	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: biomass !! Stand level biomass
1738	!! @tex $(gC m^{-2})$ @endtex
1739	CHARACTER(*),INTENT(in) :: check_point !! A flag to indicate at which
1740	!! point in the code we're doing
1741	!! this check
1742	REAL(r_std), DIMENSION(:,:), INTENT(in) :: ind !! Density of individuals
1743	!! @tex $(m^{-2})$ @endtex
1744
1745	!! 0.2 Output variables
1746	LOGICAL,INTENT(out) :: lsync
1747	REAL(r_std), DIMENSION(:,:,:), INTENT(out) :: bm_sync !! The difference betweeen the
1748	!! biomass in the circ_classes and
1749	!! the total biomass
1750	!! @tex $(gC m^{-2})$ @endtex
1751	!! 0.3 Modified variables
1752
1753	!! 0.4 Local variables
1754	INTEGER :: iele,ipts,ivm,ipar,icir
1755	REAL(r_std) :: total_circ_class_biomass
1756	LOGICAL :: lnegative
1757	!_ ================================================================================================================================
1758
1759	lsync=.TRUE.
1760	lnegative=.FALSE.
1761
1762	bm_sync(:,:,:)=zero
1763
1764	!++++++ TEMP ++++++
1765	! We gain 5-10% speed by skipping this routine
1766	! RETURN
1767	!++++++++++++
1768
1769	! Check to see if the biomass is not equal to the total biomass
1770	! in circ_class_biomass anywhere.
1771	DO ipts=1,npts
1772
1773	DO ivm=1,nvm
1774
1775	! Only woody PFTs have circumference classes therefore
1776	! only woody PFTs need to be syncronized
1777	!!$ IF (is_tree(ivm)) THEN
1778
1779	DO iele=1,icarbon
1780
1781	DO ipar=1,nparts
1782
1783	total_circ_class_biomass=zero
1784	DO icir=1,ncirc
1785
1786	total_circ_class_biomass=total_circ_class_biomass+&
1787	circ_class_biomass(ipts,ivm,icir,ipar,iele)*circ_class_n(ipts,ivm,icir)
1788
1789	! Check as well to see if our biomass is ever negative.
1790	! It really should not be.
1791	IF(circ_class_biomass(ipts,ivm,icir,ipar,iele) .LT. -min_stomate)THEN
1792
1793	lnegative=.TRUE.
1794	WRITE(numout,) '!**********************************'
1795	WRITE(numout,*) 'Error: Negative biomass component!'
1796	WRITE(numout,*) 'Check point: ',TRIM(check_point)
1797	WRITE(numout,*) 'circ_class_biomass(ipts,ivm,icir,ipar,iele) ',&
1798	circ_class_biomass(ipts,ivm,icir,ipar,iele)
1799	WRITE(numout,'(A,5I5)') 'ipts,ivm,icir,ipar,iele',ipts,ivm,icir,ipar,iele
1800	WRITE(numout,) '!**********************************'
1801
1802	ENDIF
1803
1804	ENDDO
1805
1806	IF(ABS(biomass(ipts,ivm,ipar,iele) - &
1807	total_circ_class_biomass) .GT. sync_threshold)THEN
1808
1809	WRITE(numout,) '!**********************************'
1810	WRITE(numout,*) 'Biomass and circ_class_biomass are not equal!'
1811	WRITE(numout,*) 'Check point: ',TRIM(check_point)
1812	WRITE(numout,100) 'biomass(ipts,ivm,ipar,iele) ',&
1813	biomass(ipts,ivm,ipar,iele)
1814	WRITE(numout,100) 'total_circ_class_biomass ',&
1815	total_circ_class_biomass
1816	WRITE(numout,100) 'Difference: ',&
1817	ABS(biomass(ipts,ivm,ipar,iele) - total_circ_class_biomass)
1818	WRITE(numout,*) 'ipts,ivm,ipar,iele',ipts,ivm,ipar,iele
1819	WRITE(numout,) '!**********************************'
1820	100 FORMAT(A,E20.10)
1821	lsync=.FALSE.
1822
1823	ENDIF
1824
1825	ENDDO
1826
1827	! we are not going to save the biomass for every component right now,
1828	! just the total
1829	bm_sync(ipts,ivm,iele)=zero
1830
1831	DO ipar=1,nparts
1832
1833	DO icir=1,ncirc
1834
1835	bm_sync(ipts,ivm,iele)=bm_sync(ipts,ivm,iele)+&
1836	circ_class_biomass(ipts,ivm,icir,ipar,iele)*circ_class_n(ipts,ivm,icir)
1837	ENDDO ! ncirc
1838
1839	ENDDO ! nparts
1840
1841	bm_sync(ipts,ivm,iele)=ABS(bm_sync(ipts,ivm,iele)-&
1842	SUM(biomass(ipts,ivm,:,iele)))
1843
1844	ENDDO ! nelements
1845
1846	!!$ ENDIF ! is_tree
1847
1848	ENDDO ! loop over PFTs
1849
1850	ENDDO ! loop over points
1851
1852	!---TEMP---
1853	IF(ld_biomass)THEN
1854	WRITE(numout,*) 'Check point: ',TRIM(check_point)
1855	WRITE(numout,*) 'test_pft, test_grid: ',test_pft,test_grid
1856	!!$ WRITE(numout,*) 'biomass (ileaf), ', biomass(test_grid,test_pft,ileaf,icarbon)
1857	!!$ WRITE(numout,*) 'biomass (iwood), ', biomass(test_grid,test_pft,isapabove,icarbon) + &
1858	!!$ biomass(test_grid,test_pft,isapbelow,icarbon) + biomass(test_grid,test_pft,iheartabove,icarbon) + &
1859	!!$ biomass(test_grid,test_pft,iheartbelow,icarbon)
1860	!!$ WRITE(numout,*) 'biomass (iroot), ', biomass(test_grid,test_pft,iroot,icarbon)
1861	WRITE(numout,'(A,20F14.6)') 'biomassHHH, ',biomass(test_grid,test_pft,:,icarbon)
1862	DO icir=1,ncirc
1863	WRITE(numout,'(A,I1,20F14.6)') 'ccbiomass',icir,circ_class_biomass(test_grid,test_pft,icir,:,icarbon)
1864	ENDDO
1865	WRITE(numout,*) 'circ_class_biomass, ',&
1866	SUM (SUM(circ_class_biomass(test_grid,test_pft,:,:,icarbon),2) * &
1867	circ_class_n(test_grid,test_pft,:))
1868	WRITE(numout,*) 'circ_class_n, ', SUM(circ_class_n(test_grid,test_pft,:))
1869	WRITE(numout,*) 'circ_class_n(:), ', circ_class_n(test_grid,test_pft,:)
1870	WRITE(numout,*) 'ind, ', ind(test_grid,test_pft)
1871	ENDIF
1872	!!$ !----------
1873
1874	IF(.NOT. lsync) THEN
1875	WRITE(numout,*) 'ERROR: stopping in sync'
1876	WRITE(numout,*) 'Stopping'
1877	STOP
1878	ENDIF
1879	IF(lnegative) THEN
1880	WRITE(numout,*) 'ERROR: negative biomass'
1881	WRITE(numout,*) 'Stopping'
1882	STOP
1883	ENDIF
1884
1885	END SUBROUTINE check_biomass_sync
1886
1887	END MODULE function_library
1888

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

Download in other formats: