Context Navigation

source: tags/ORCHIDEE_4_1/ORCHIDEE/src_parameters/function_library.f90 @ 7704

Last change on this file since 7704 was 7490, checked in by sebastiaan.luyssaert, 2 years ago
Contributes to ticket 795. Mostly a rewrite of the changes committed in r7450. Contains an extra mass balance fix compared to r7450
File size: 150.1 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 pft_parameters
30	USE constantes
31
32	IMPLICIT NONE
33
34	! private & public routines
35
36	PRIVATE
37	PUBLIC calculate_c0_alloc, &
38	wood_to_ba_eff, &
39	wood_to_ba, &
40	wood_to_height, &
41	wood_to_qmheight, &
42	wood_to_height_eff, &
43	wood_to_volume, &
44	wood_to_dia, &
45	wood_to_qmdia, &
46	wood_to_circ, &
47	wood_to_cn_dia, &
48	wood_to_cv, &
49	cc_to_biomass, &
50	biomass_to_cc, &
51	cc_to_lai, &
52	biomass_to_lai, &
53	lai_to_biomass, &
54	Nmax, &
55	calculate_rdi, &
56	check_vegetation_area, &
57	check_mass_balance, &
58	intermediate_mass_balance_check, &
59	check_pixel_area, &
60	check_area_change, &
61	check_variable_change, &
62	check_variable_swap, &
63	check_biomass_change, &
64	biomass_to_coupled_lai, &
65	sort_circ_class_biomass, &
66	Arrhenius, &
67	Arrhenius_modified, &
68	Arrhenius_modified_nodim
69
70	INTERFACE biomass_to_lai
71	MODULE PROCEDURE biomass_to_lai_0d, biomass_to_lai_1d
72	END INTERFACE
73
74	INTERFACE lai_to_biomass
75	MODULE PROCEDURE lai_to_biomass_0d, lai_to_biomass_1d
76	END INTERFACE
77
78	INTERFACE cc_to_biomass
79	MODULE PROCEDURE cc_to_biomass_2d, cc_to_biomass_4d
80	END INTERFACE cc_to_biomass
81
82	INTERFACE biomass_to_cc
83	MODULE PROCEDURE biomass_to_cc_1d, biomass_to_cc_3d
84	END INTERFACE biomass_to_cc
85
86	INTERFACE cc_to_lai
87	MODULE PROCEDURE cc_to_lai_0d, cc_to_lai_1d
88	END INTERFACE
89
90	INTERFACE Arrhenius
91	MODULE PROCEDURE Arrhenius_0d, Arrhenius_1d
92	END INTERFACE
93
94	INTERFACE Arrhenius_modified
95	MODULE PROCEDURE Arrhenius_modified_0d, Arrhenius_modified_1d
96	END INTERFACE
97
98	CONTAINS
99
100
101	!! ================================================================================================================================
102	!! FUNCTION : biomass_to_coupled_lai
103	!!
104	!>\BRIEF Calculate the coupled_LAI based on biomass and veget of each pft
105	!!
106	!! DESCRIPTION : Calculates the lai that is coupled with the atmosphere
107	!!
108	!!
109	!!
110	!! RECENT CHANGE(S): None
111	!!
112	!! RETURN VALUE : ::coupled_lai (m2/m2)
113	!!
114	!! REFERENCE(S) :
115	!!
116	!! FLOWCHART : None
117	!! \n
118	!_ ================================================================================================================================
119
120	FUNCTION biomass_to_coupled_lai(biomass_leaf, veget, pft)
121
122	!! 0. Variable and parameter declaration
123
124	!! 0.1 Input variables
125	INTEGER(i_std) :: pft !! PFT number (-)
126	REAL(r_std) :: veget !! 1-Pgap or the vegetation cover
127	REAL(r_std) :: biomass_leaf !! Biomass of an individual tree within a circ
128	!! class @tex $(m^{2} ind^{-1})$ @endtex
129
130	!! 0.2 Output variables
131	REAL(r_std) :: biomass_to_coupled_lai !! What we get out of this function, depends on the on which
132	!! leaf biomass is passed to the function;
133	!! is it at the tree or the stand level. At the tree level
134	!! it corresponds to the leaf area per tree (m2). At the
135	!! stand levels it is the fraction of the LAI that is
136	!! assumed to interact with the atmosphere
137	!! @tex $(m^{2} m^{-2})$ @endtex
138	!! 0.3 Modified variables
139
140	!! 0.4 Local variables
141	REAL(r_std) :: lai !! Leaf area index
142	!! @tex $(m^{2} m^{-2})$ @endtex
143	!_ ================================================================================================================================
144
145	lai = biomass_to_lai(biomass_leaf, pft)
146
147	!+++CHECK+++
148	! In a closed canopy not all the leaves fully interact with the
149	! atmosphere because leaves can shelter each other. In a more
150	! open canopy most leaves can interact with the atmosphere.
151	! For the moment we are not clear how to calculate the fraction of
152	! the canopy that is coupled to the atmosphere. Also, based on the
153	! simulated evapotranspiration we need the entire canopy to be
154	! coupled to the atmosphere to obtain simulations which are close
155	! to the observations (Jung's upscaled product).
156	IF (veget .GT. min_sechiba) THEN
157
158	biomass_to_coupled_lai = lai
159
160	ELSE
161
162	! There are so many gaps that veget is extremely small. It is fair to
163	! assume that the whole canopy is coupled to the atmosphere.
164	biomass_to_coupled_lai = lai
165
166	ENDIF
167
168	END FUNCTION biomass_to_coupled_lai
169
170
171	!! ================================================================================================================================
172	!! FUNCTION : biomass_to_lai_0d
173	!!
174	!>\BRIEF Calculate the LAI based on the leaf biomass
175	!!
176	!! DESCRIPTION : Calculates the LAI of a PFT/grid square based on the leaf biomass
177	!!
178	!! RECENT CHANGE(S): None
179	!!
180	!! RETURN VALUE : ::LAI [m2 m{-2}]
181	!!
182	!! REFERENCE(S) :
183	!!
184	!! FLOWCHART : None
185	!! \n
186	!_ ================================================================================================================================
187
188	FUNCTION biomass_to_lai_0d(leaf_biomass, pft)
189
190	!! 0. Variable and parameter declaration
191
192	!! 0.1 Input variables
193
194	INTEGER(i_std) :: pft !! PFT number (-)
195	REAL(r_std) :: leaf_biomass !! Biomass of the leaves
196	!! @tex $(gC m^{-2})$ @endtex
197
198
199	!! 0.2 Output variables
200
201	REAL(r_std) :: biomass_to_lai_0d !! Leaf area index
202	!! @tex $(m^{2} m^{-2})$ @endtex
203
204	!! 0.3 Modified variables
205
206	!! 0.4 Local variables
207	REAL(r_std) :: impose_lai !! LAI read from run.def
208	!_ ================================================================================================================================
209
210	!! 1. Calculate the LAI from the leaf biomass
211	IF(sla_dyn) THEN
212	biomass_to_lai_0d = log(1.+ ext_coeff_N(pft) * leaf_biomass * slainit(pft))/(ext_coeff_N(pft))
213	ELSE
214	biomass_to_lai_0d = leaf_biomass * sla(pft)
215	ENDIF
216
217	!!$ !+++HACK+++
218	!!$ ! This is the perfect place to hack the code to make it run with
219	!!$ ! constant lai
220	!!$ WRITE(numout,*) 'WARNING: Using fake lai values for testing!'
221	!!$ biomass_to_lai_1d(:)=3.79052
222	!!$ !++++++++++
223
224	END FUNCTION biomass_to_lai_0d
225
226
227	!! ================================================================================================================================
228	!! FUNCTION : biomass_to_lai_1d
229	!!
230	!>\BRIEF Calculate the LAI based on the leaf biomass
231	!!
232	!! DESCRIPTION : Calculates the LAI of a PFT/grid square based on the leaf biomass
233	!!
234	!! RECENT CHANGE(S): None
235	!!
236	!! RETURN VALUE : ::LAI [m2 m{-2}]
237	!!
238	!! REFERENCE(S) :
239	!!
240	!! FLOWCHART : None
241	!! \n
242	!_ ================================================================================================================================
243
244	FUNCTION biomass_to_lai_1d(leaf_biomass, npts, pft)
245
246	!! 0. Variable and parameter declaration
247
248	!! 0.1 Input variables
249	INTEGER(i_std) :: npts
250	INTEGER(i_std) :: pft !! PFT number (-)
251	REAL(r_std), DIMENSION(npts) :: leaf_biomass !! Biomass of the leaves
252	!! @tex $(gC m^{-2})$ @endtex
253
254	!! 0.2 Output variables
255
256	REAL(r_std), DIMENSION(npts) :: biomass_to_lai_1d !! Leaf area index
257	!! @tex $(m^{2} m^{-2})$ @endtex
258	!_ ================================================================================================================================
259
260	!! 1. Calculate the LAI from the leaf biomass
261	IF(sla_dyn) THEN
262	biomass_to_lai_1d = log(1.+ ext_coeff_N(pft) * leaf_biomass * slainit(pft))/(ext_coeff_N(pft))
263	ELSE
264	biomass_to_lai_1d = leaf_biomass * sla(pft)
265	ENDIF
266
267	!!$ !+++HACK+++
268	!!$ ! This is the perfect place to hack the code to make it run with
269	!!$ ! constant lai
270	!!$ WRITE(numout,*) 'WARNING: Using fake lai values for testing!'
271	!!$ biomass_to_lai_1d(:)=3.79052
272	!!$ !++++++++++
273
274	END FUNCTION biomass_to_lai_1d
275
276
277	!! ================================================================================================================================
278	!! FUNCTION : lai_to_biomass_0d
279	!!
280	!>\BRIEF Calculate leaf biomass from lai
281	!!
282	!! DESCRIPTION : Calculates leaf biomass (m2 m-2) from leaf biomass (gC m-2)
283	!!
284	!! RECENT CHANGE(S): None
285	!!
286	!! RETURN VALUE : ::Leaf biomass [gC m^{-2}]
287	!!
288	!! REFERENCE(S) :
289	!!
290	!! FLOWCHART : None
291	!! \n
292	!_ ================================================================================================================================
293
294	FUNCTION lai_to_biomass_0d(lai, pft)
295
296	!! 0. Variable and parameter declaration
297
298	!! 0.1 Input variables
299	INTEGER(i_std) :: pft !! PFT number (-)
300	REAL(r_std) :: lai !! Leaf Area Index
301	!! @tex $(m^{2} m^{-2})$ @endtex
302
303
304	!! 0.2 Output variables
305
306	REAL(r_std) :: lai_to_biomass_0d !! Leaf area index
307	!! @tex $(gC m^{-2})$ @endtex
308
309	!_ ================================================================================================================================
310
311	!! 1. Calculate the LAI from the leaf biomass
312	IF(sla_dyn) THEN
313	lai_to_biomass_0d = ( exp(lai*ext_coeff_N(pft)) - 1.) / &
314	(ext_coeff_N(pft) * slainit(pft))
315	ELSE
316	lai_to_biomass_0d = lai / sla(pft)
317	ENDIF
318
319	END FUNCTION lai_to_biomass_0d
320
321
322
323	!! ================================================================================================================================
324	!! FUNCTION : lai_to_biomass_1d
325	!!
326	!>\BRIEF Calculate leaf biomass from lai
327	!!
328	!! DESCRIPTION : Calculates leaf biomass (m2 m-2) from leaf biomass (gC m-2)
329	!!
330	!! RECENT CHANGE(S): None
331	!!
332	!! RETURN VALUE : ::Leaf biomass [gC m^{-2}]
333	!!
334	!! REFERENCE(S) :
335	!!
336	!! FLOWCHART : None
337	!! \n
338	!_ ================================================================================================================================
339
340	FUNCTION lai_to_biomass_1d(lai, npts, pft)
341
342	!! 0. Variable and parameter declaration
343
344	!! 0.1 Input variables
345	INTEGER(i_std) :: npts
346	INTEGER(i_std) :: pft !! PFT number (-)
347	REAL(r_std), DIMENSION(npts) :: lai !! Leaf Area Index
348	!! @tex $(m^{2} m^{-2})$ @endtex
349
350
351	!! 0.2 Output variables
352
353	REAL(r_std), DIMENSION(npts) :: lai_to_biomass_1d !! Leaf area index
354	!! @tex $(gC m^{-2})$ @endtex
355
356	!_ ================================================================================================================================
357
358	!! 1. Calculate the LAI from the leaf biomass
359	IF(sla_dyn) THEN
360	lai_to_biomass_1d = ( exp(lai*ext_coeff_N(pft)) - 1.) / &
361	(ext_coeff_N(pft) * slainit(pft))
362	ELSE
363	lai_to_biomass_1d = lai / sla(pft)
364	ENDIF
365
366	END FUNCTION lai_to_biomass_1d
367
368
369	!! ================================================================================================================================
370	!! FUNCTION : calculate_c0_alloc
371	!!
372	!>\BRIEF Calculate the baseline root vs sapwood allocation
373	!!
374	!! DESCRIPTION : Calculates the baseline root vs sapwood allocation based on the
375	!! parameters of the pipe model (hydraulic conductivities) and the
376	!! turnover of the different components
377	!!
378	!! RECENT CHANGE(S): None
379	!!
380	!! RETURN VALUE : ::c0_alloc (m)
381	!!
382	!! REFERENCE(S) :
383	!!
384	!! FLOWCHART : None
385	!! \n
386	!_ ================================================================================================================================
387
388	FUNCTION calculate_c0_alloc(pts, pft, longevity_eff_root, longevity_eff_sap)
389
390	!! 0. Variable and parameter declaration
391
392	!! 0.1 Input variables
393
394	INTEGER(i_std) :: pts !! Pixel number (-)
395	INTEGER(i_std) :: pft !! PFT number (-)
396	REAL(r_std) :: longevity_eff_root !! Effective longivety for leaves (days)
397	REAL(r_std) :: longevity_eff_sap !! Effective longivety for leaves (days)
398
399	!! 0.2 Output variables
400
401	REAL(r_std) :: calculate_c0_alloc !! quadratic mean height (m)
402
403	!! 0.3 Modified variables
404
405	!! 0.4 Local variables
406	REAL(r_std) :: sapwood_density
407	REAL(r_std) :: qm_dia !! quadratic mean diameter (m)
408
409	!_ ================================================================================================================================
410
411	!! 1. Calculate c0_alloc
412	IF ( is_tree(pft) ) THEN
413
414	sapwood_density = deux * pipe_density(pft) / kilo_to_unit
415	calculate_c0_alloc = sqrt(k_belowground(pft)/k_sap(pft)longevity_eff_sap/longevity_eff_rootsapwood_density)
416
417	! Grasses and croplands
418	ELSE
419
420	!+++CHECK+++
421	! Simply copied the same formulation as for trees but note
422	! that the sapwood in trees vs grasses and crops has a very
423	! meaning. In grasses and crops is structural carbon to ensure
424	! that the allocation works. In trees it really is the sapwood
425	sapwood_density = deux * pipe_density(pft) / kilo_to_unit
426	calculate_c0_alloc = sqrt(k_belowground(pft)/k_sap(pft)longevity_eff_sap/longevity_eff_rootsapwood_density)
427	!+++++++++++
428
429	ENDIF ! is_tree(j)
430
431	END FUNCTION calculate_c0_alloc
432
433
434	!! ================================================================================================================================
435	!! FUNCTION : wood_to_ba_eff
436	!!
437	!>\BRIEF Calculate effective basal area from woody biomass making use of allometric relationships
438	!!
439	!! DESCRIPTION : Calculate basal area of an individual tree from the woody biomass of that tree making
440	!! use of allometric relationships. Effective basal area accounts for both above and below ground carbon
441	!! and is the basis for the application of the rule of Deleuze and Dhote.
442	!! (i) woodmass = tree_ff * pipe_densitybaheight
443	!! (ii) height = pipe_tune2 * sqrt(4/piba) * pipe_tune_3
444	!!
445	!! RECENT CHANGE(S): None
446	!!
447	!! RETURN VALUE : effective basal area (m2 ind-1)
448	!!
449	!! REFERENCE(S) :
450	!!
451	!! FLOWCHART : None
452	!! \n
453	!_ ================================================================================================================================
454
455	FUNCTION wood_to_ba_eff(biomass_temp, pft)
456
457	!! 0. Variable and parameter declaration
458
459	!! 0.1 Input variables
460
461	INTEGER(i_std) :: pft !! PFT number (-)
462	REAL(r_std), DIMENSION(:,:) :: biomass_temp !! Biomass of an individual tree within a circ
463	!! class @tex $(m^{2} ind^{-1})$ @endtex
464
465	!! 0.2 Output variables
466
467	REAL(r_std), DIMENSION(ncirc) :: wood_to_ba_eff !! Effective basal area of an individual tree within a circ
468	!! class @tex $(m^{2} ind^{-1})$ @endtex
469
470	!! 0.3 Modified variables
471
472	!! 0.4 Local variables
473	INTEGER(i_std) :: l !! Index
474	REAL(r_std) :: woodmass_ind !! Woodmass of an individual tree
475	!! @tex $(gC ind^{-1})$ @endtex
476	REAL(r_std) :: temp !! temporary variable
477
478	!_ ================================================================================================================================
479
480	!! 1. Calculate basal area from woodmass
481
482	IF ( is_tree(pft) ) THEN
483
484	DO l = 1,ncirc
485
486	! Woodmass of an individual tree
487	woodmass_ind = biomass_temp(l,isapabove) + biomass_temp(l,isapbelow) + &
488	biomass_temp(l,iheartabove) + biomass_temp(l,iheartbelow)
489
490	temp = woodmass_ind/(tree_ff(pft)pipe_density(pft)pipe_tune2(pft))
491
492	IF (temp.LT.EPSILON(un)) THEN
493	! Expect an overflow when taking an exponent of such a small number
494	wood_to_ba_eff(l) = min_stomate
495	ELSE
496	! Basal area of that individual (m2 ind-1)
497	wood_to_ba_eff(l) =(pi/4(woodmass_ind/(tree_ff(pft)pipe_density(pft)*pipe_tune2(pft))) &
498	(2./pipe_tune3(pft)))(pipe_tune3(pft)/(pipe_tune3(pft)+2))
499
500	ENDIF
501
502	ENDDO
503
504	ELSE
505
506	WRITE(numout,*) 'pft ',pft
507	CALL ipslerr_p (3,'wood_to_ba_eff', &
508	'wood_to_ba_eff is not defined for this PFT.', &
509	'See the output file for more details.','')
510
511	ENDIF
512
513	END FUNCTION wood_to_ba_eff
514
515
516	!! ================================================================================================================================
517	!! FUNCTION : wood_to_ba
518	!!
519	!>\BRIEF Calculate basal area from woody biomass making use of allometric relationships
520	!!
521	!! DESCRIPTION : Calculate basal area of an individual tree from the woody biomass of that tree making
522	!! use of allometric relationships given below. Here basal area is defined in line with its classical
523	!! forestry meaning.
524	!! (i) woodmass = tree_ff * pipe_densitybaheight
525	!! (ii) height = pipe_tune2 * sqrt(4/piba) * pipe_tune_3
526	!!
527	!! RECENT CHANGE(S): None
528	!!
529	!! RETURN VALUE : basal area (m2 ind-1)
530	!!
531	!! REFERENCE(S) :
532	!!
533	!! FLOWCHART : None
534	!! \n
535	!_ ================================================================================================================================
536
537	FUNCTION wood_to_ba(biomass_temp, pft)
538
539	!! 0. Variable and parameter declaration
540
541	!! 0.1 Input variables
542
543	INTEGER(i_std) :: pft !! PFT number (-)
544	REAL(r_std), DIMENSION(:,:) :: biomass_temp !! Biomass of an individual tree within a circ
545	!! class @tex $(m^{2} ind^{-1})$ @endtex
546
547	!! 0.2 Output variables
548
549	REAL(r_std), DIMENSION(ncirc) :: wood_to_ba !! Basal area of an individual tree within a circ
550	!! class @tex $(m^{2} ind^{-1})$ @endtex
551
552	!! 0.3 Modified variables
553
554	!! 0.4 Local variables
555	INTEGER(i_std) :: l !! Index
556	REAL(r_std) :: woodmass_ind !! Woodmass of an individual tree
557	!! @tex $(gC ind^{-1})$ @endtex
558	REAL(r_std) :: temp !! Temporary variable
559
560	!_ ================================================================================================================================
561
562	!! 1. Calculate basal area from woodmass
563
564	IF ( is_tree(pft) ) THEN
565
566	DO l = 1,ncirc
567
568	! Woodmass of an individual tree
569	woodmass_ind = biomass_temp(l,iheartabove) + biomass_temp(l,isapabove)
570
571	temp = woodmass_ind/(tree_ff(pft)pipe_density(pft)pipe_tune2(pft))
572
573	IF (temp.LT.EPSILON(un)) THEN
574	! Expect an overflow when taking an exponent of such a small number
575	wood_to_ba(l) = min_stomate
576	ELSE
577	! Basal area of that individual (m2 ind-1)
578	wood_to_ba(l) = (pi/4(woodmass_ind/(tree_ff(pft)pipe_density(pft)*pipe_tune2(pft))) &
579	(2./pipe_tune3(pft)))(pipe_tune3(pft)/(pipe_tune3(pft)+2))
580	ENDIF
581
582	ENDDO
583
584	ELSE
585
586	WRITE(numout,*) 'pft ',pft
587	CALL ipslerr_p (3,'wood_to_ba', &
588	'wood_to_ba is not defined for this PFT.', &
589	'See the output file for more details.','')
590
591	ENDIF
592
593	END FUNCTION wood_to_ba
594
595
596	!! ================================================================================================================================
597	!! FUNCTION : ccbio_to_biomass
598	!!
599	!>\BRIEF Calculate total biomass from circ_class_biomass
600	!!
601	!! DESCRIPTION : circ_class_biomass is expressed in gC tree-1 and thus needs to
602	!! be multiplied by the number of trees per circumference class
603	!! circ_class_n expressed in tree m-2 to calculate the total
604	!! biomass in gC m-2.
605	!!
606	!! RECENT CHANGE(S): None
607	!!
608	!! RETURN VALUE : total biomass (gC(N) m-2)
609	!!
610	!! REFERENCE(S) :
611	!!
612	!! FLOWCHART : None
613	!! \n
614	!_ ================================================================================================================================
615
616	FUNCTION cc_to_biomass_2d(npts,ivm,ccbio,ccind)
617
618	!! 0. Variable and parameter declaration
619
620	!! 0.1 Input variables
621
622	INTEGER(i_std), INTENT(in) :: npts !! Domain size (unitless)
623	INTEGER(i_std), INTENT(in) :: ivm !! number of the specific PFT (unitless)
624	REAL(r_std), DIMENSION(:,:,:) :: ccbio !! circ_class_biomass. Biomass of an individual
625	!! tree @tex $(gC(N) ind^{-1})$ @endtex
626	REAL(r_std), DIMENSION(:) :: ccind !! circ_class_n. Number of trees per circumference
627	!! class
628
629	!! 0.2 Output variables
630
631	REAL(r_std), DIMENSION(nparts,nelements) :: cc_to_biomass_2d !! biomass @tex $(gC(N) m^{2})$ @endtex
632
633	!! 0.3 Modified variables
634
635	!! 0.4 Local variables
636	INTEGER(i_std) :: iele,icir,ipar !! Index
637
638	!_ ================================================================================================================================
639
640	cc_to_biomass_2d(:,:) = zero
641
642	DO iele = 1,nelements
643
644	! Note that when used in 2d ivm is a single PFT
645	! Use MAX to prevent precision (-10e-16) issues.
646	IF ( is_tree(ivm) ) THEN
647
648	DO ipar = 1,nparts
649
650	DO icir = 1,ncirc
651
652	cc_to_biomass_2d(ipar,iele) = MAX(zero,cc_to_biomass_2d(ipar,iele) + &
653	ccbio(icir,ipar,iele) * ccind(icir))
654
655	ENDDO
656
657	ENDDO
658
659	ELSE
660
661	! Grasses and cropland only have one circumference class
662	DO ipar = 1,nparts
663
664	cc_to_biomass_2d(ipar,iele) = MAX(zero,cc_to_biomass_2d(ipar,iele) + &
665	ccbio(1,ipar,iele) * ccind(1))
666
667	ENDDO
668
669	ENDIF
670
671	ENDDO
672
673	END FUNCTION cc_to_biomass_2d
674
675
676	!! ================================================================================================================================
677
678	FUNCTION cc_to_biomass_4d(npts,nvm,ccbio,ccind)
679
680	!! 0. Variable and parameter declaration
681
682	!! 0.1 Input variables
683
684	INTEGER(i_std), INTENT(in) :: npts !! Domain size (unitless)
685	INTEGER(i_std), INTENT(in) :: nvm !! Total number of PFTs (unitless)
686	REAL(r_std), DIMENSION(:,:,:,:,:) :: ccbio !! circ_class_biomass. Biomass of an individual
687	!! tree @tex $(gC(N) ind^{-1})$ @endtex
688	REAL(r_std), DIMENSION(:,:,:) :: ccind !! circ_class_n. Number of trees per circumference
689	!! class
690
691	!! 0.2 Output variables
692
693	REAL(r_std), DIMENSION(npts,nvm,nparts,nelements) :: cc_to_biomass_4d !! biomass @tex $(gC(N) m^{2})$ @endtex
694
695	!! 0.3 Modified variables
696
697	!! 0.4 Local variables
698	INTEGER(i_std) :: iele,icir,ipar,ivm !! Index
699
700	!_ ================================================================================================================================
701
702	cc_to_biomass_4d(:,:,:,:) = zero
703
704	DO iele = 1,nelements
705
706	! Note that when used in 4d, nvm rather than ivm
707	! should be passed. Now the function will loop over all
708	! PFTs instead of just 1 PFT for the 2d version of this
709	! function. Use MAX to prevent precision (-10e-16) issues.
710	DO ivm = 2,nvm
711
712	IF ( is_tree(ivm) ) THEN
713
714	DO ipar = 1,nparts
715
716	DO icir = 1,ncirc
717
718	cc_to_biomass_4d(:,ivm,ipar,iele) = MAX(zero,cc_to_biomass_4d(:,ivm,ipar,iele) + &
719	ccbio(:,ivm,icir,ipar,iele) * ccind(:,ivm,icir))
720
721	ENDDO
722
723	ENDDO
724
725	ELSE
726
727	! Grasses and cropland only have one circumference class
728	DO ipar = 1,nparts
729
730	cc_to_biomass_4d(:,ivm,ipar,iele) = MAX(zero,cc_to_biomass_4d(:,ivm,ipar,iele) + &
731	ccbio(:,ivm,1,ipar,iele) * ccind(:,ivm,1))
732
733	ENDDO
734
735	ENDIF
736
737	ENDDO
738
739	ENDDO
740
741	END FUNCTION cc_to_biomass_4d
742
743
744	!! ================================================================================================================================
745	!! FUNCTION : biomass_to_cc
746	!!
747	!>\BRIEF Calculate circ_class_biomass from biomass
748	!!
749	!! DESCRIPTION : The labile and carbres pools are often dealt with at the stand
750	!! level for simplicity but then need to be stored at the plant
751	!! level. This is needed because labile and carbres C is lost when
752	!! plants are dying. Also, all plant C and N should be stored in a
753	!! single variable.
754	!!
755	!! RECENT CHANGE(S): None
756	!!
757	!! RETURN VALUE : total biomass (gC(N) m-2)
758	!!
759	!! REFERENCE(S) :
760	!!
761	!! FLOWCHART : None
762	!! \n
763	!_ ================================================================================================================================
764
765	FUNCTION biomass_to_cc_1d(bio,ccbio,ccind)
766
767	!! 0. Variable and parameter declaration
768
769	!! 0.1 Input variables
770
771	REAL(r_std) :: bio !! biomass. Biomass of a specific pool of
772	!! a whole stand
773	!! @tex $(gC(N) m^{-2})$ @endtex
774	REAL(r_std), DIMENSION(:) :: ccbio !! circ_class_biomass. Biomass of a specif pool of
775	!! the individual trees in each circumference class
776	!! @tex $(gC(N) ind^{-1})$ @endtex
777	REAL(r_std), DIMENSION(:) :: ccind !! circ_class_n. Number of trees per circumference
778	!! class
779
780	!! 0.2 Output variables
781
782	REAL(r_std), DIMENSION(ncirc) :: biomass_to_cc_1d !! biomass @tex $(gC(N) m^{2})$ @endtex
783
784	!! 0.3 Modified variables
785
786	!! 0.4 Local variables
787
788	INTEGER(i_std) :: icir !! Index
789	REAL(r_std), DIMENSION(ncirc) :: share_ncirc !! ratio of specific pool
790	!! across circumference classes
791
792	!_ ================================================================================================================================
793
794	biomass_to_cc_1d(:) = zero
795
796	! Convert the pools
797	IF (SUM(ccbio(:)*ccind(:)) .GT. zero) THEN
798
799	! Proportional to the present pools
800	share_ncirc(:) = (ccbio(:)ccind(:))/SUM(ccbio(:)ccind(:))
801
802	ELSEIF (SUM(ccind(:)) .GT. zero) THEN
803
804	! There are no pools at present. Make it proportional
805	! to the number of trees
806	share_ncirc(:) = ccind(:)/SUM(ccind(:))
807
808	ELSE
809
810	! There is no biomass in circ_class_biomass or individuals
811	share_ncirc(:) = zero
812
813
814	ENDIF
815
816	share_ncirc(1) = zero
817	share_ncirc(1) = 1 - SUM(share_ncirc(:))
818
819	! Distribute bio over the ncirc circumference classes of biomass_to_cc
820	! biomass_to_cc will eventually be assigned to circ_class_biomass
821	DO icir = 1,ncirc
822
823	IF (ccind(icir) .GT. zero) THEN
824
825	biomass_to_cc_1d(icir) = bio / ccind(icir) * share_ncirc(icir)
826
827	ELSE
828
829	biomass_to_cc_1d(icir) = zero
830
831	ENDIF
832
833	ENDDO
834
835	END FUNCTION biomass_to_cc_1d
836
837	!! ================================================================================================================================
838	FUNCTION biomass_to_cc_3d(bio,ccbio,ccind,npts,nvm)
839
840	!! 0. Variable and parameter declaration
841
842	!! 0.1 Input variables
843	INTEGER(i_std), INTENT(in) :: npts !! Domain size (unitless)
844	INTEGER(i_std), INTENT(in) :: nvm !! Total number of PFTs (unitless)
845	REAL(r_std),DIMENSION(:,:) :: bio !! biomass. Biomass of a specific pool of
846	!! a whole stand. dimensions are npts and nvm
847	!! @tex $(gC(N) m^{-2})$ @endtex
848	REAL(r_std), DIMENSION(:,:,:) :: ccbio !! circ_class_biomass. Biomass of a specif pool of
849	!! the individual trees in each circumference class
850	!! @tex $(gC(N) ind^{-1})$ @endtex
851	REAL(r_std), DIMENSION(:,:,:) :: ccind !! circ_class_n. Number of trees per circumference
852	!! class
853
854	!! 0.2 Output variables
855
856	REAL(r_std), DIMENSION(npts,nvm,ncirc) :: biomass_to_cc_3d !! biomass @tex $(gC(N) m^{2})$ @endtex
857
858	!! 0.3 Modified variables
859
860	!! 0.4 Local variables
861
862	INTEGER(i_std) :: icir,ipts,ivm !! Index
863	REAL(r_std), DIMENSION(ncirc) :: share_ncirc !! ratio of specific pool
864	!! across circumference classes
865
866	!_ ================================================================================================================================
867
868	biomass_to_cc_3d(:,:,:) = zero
869
870	DO ipts = 1,npts
871
872	DO ivm = 1,nvm
873
874	! Convert the pools
875	IF (is_tree(ivm)) THEN
876
877	! Forests have ncirc circ classes
878	IF (SUM(ccbio(ipts,ivm,:)*ccind(ipts,ivm,:)) .GT. zero) THEN
879
880	share_ncirc(:) = (ccbio(ipts,ivm,:)*ccind(ipts,ivm,:)) / &
881	SUM(ccbio(ipts,ivm,:)*ccind(ipts,ivm,:))
882
883	ELSE
884
885	share_ncirc(:) = zero
886
887	ENDIF
888
889	ELSE
890
891	! Grass and crops
892	share_ncirc(:) = zero
893
894	ENDIF
895
896	! Make sure that the shares add up to one
897	! and that grass and croplands are properly
898	! dealt with
899	share_ncirc(1) = zero
900	share_ncirc(1) = 1 - SUM(share_ncirc(:))
901
902	DO icir = 1,ncirc
903
904	IF (ccind(ipts,ivm,icir) .GT. zero) THEN
905
906	biomass_to_cc_3d(ipts,ivm,icir) = bio(ipts,ivm) / &
907	ccind(ipts,ivm,icir) * share_ncirc(icir)
908
909	ELSE
910
911	biomass_to_cc_3d(ipts,ivm,icir) = zero
912
913	ENDIF ! plants are present in this circ_class
914
915	ENDDO ! icir
916
917	ENDDO ! ivm
918
919	ENDDO ! ipts
920
921	END FUNCTION biomass_to_cc_3d
922
923	!! ================================================================================================================================
924
925
926
927
928
929
930
931	!! ================================================================================================================================
932	!! FUNCTION : cc_to_lai_0d
933	!!
934	!>\BRIEF Calculate the LAI based on the circ_class_biomass and circ_class_n
935	!!
936	!! DESCRIPTION : Calculates the LAI of a PFT/grid square based on the leaf biomass
937	!!
938	!! RECENT CHANGE(S): None
939	!!
940	!! RETURN VALUE : ::LAI [m2 m{-2}]
941	!!
942	!! REFERENCE(S) :
943	!!
944	!! FLOWCHART : None
945	!! \n
946	!_ ================================================================================================================================
947
948	FUNCTION cc_to_lai_0d(ccbio,ccind,pft)
949
950	!! 0. Variable and parameter declaration
951
952	!! 0.1 Input variables
953	INTEGER(i_std), INTENT(in) :: pft !! Number of the PFT under study (unitless)
954	REAL(r_std), DIMENSION(:) :: ccbio !! circ_class_biomass. Biomass of an individual
955	!! tree @tex $(gC(N) ind^{-1})$ @endtex
956	REAL(r_std), DIMENSION(:) :: ccind !! circ_class_n. Number of trees per circumference
957	!! class
958
959	!! 0.2 Output variables
960	REAL(r_std) :: cc_to_lai_0d !! leaf area index @tex $(m^{2} m^{-2})$ @endtex
961
962	!! 0.3 Modified variables
963
964	!! 0.4 Local variables
965	INTEGER(i_std) :: icir !! Index
966	REAL(r_std) :: leaf_biomass !! Leaf mass for carbon calculated from
967	!! circ_class_biomass and circ_class_n.
968	!! @tex $(gC m^{-2})$ @endtex
969	!_ ================================================================================================================================
970
971	!! 1. Calculate leaf biomass
972	cc_to_lai_0d = zero
973	leaf_biomass = zero
974
975	IF ( is_tree(pft) ) THEN
976
977	DO icir = 1,ncirc
978
979	leaf_biomass = leaf_biomass + ccbio(icir) * ccind(icir)
980
981	ENDDO
982
983	ELSE
984
985	! Grasses and cropland only have one circumference class
986	leaf_biomass = ccbio(1) * ccind(1)
987
988	ENDIF
989
990	!! 2. Calculate the LAI from the leaf biomass
991	IF(sla_dyn) THEN
992
993	cc_to_lai_0d = log(1.+ ext_coeff_N(pft) * leaf_biomass * &
994	slainit(pft))/(ext_coeff_N(pft))
995
996	ELSE
997
998	cc_to_lai_0d = leaf_biomass * sla(pft)
999
1000	ENDIF
1001
1002	END FUNCTION cc_to_lai_0d
1003
1004
1005	!! ================================================================================================================================
1006	!! FUNCTION : cc_to_lai_1d
1007	!!
1008	!>\BRIEF Calculate the LAI based on the circ_class_biomass and circ_class_n
1009	!!
1010	!! DESCRIPTION : Calculates the LAI of a PFT/grid square based on the leaf biomass
1011	!!
1012	!! RECENT CHANGE(S): None
1013	!!
1014	!! RETURN VALUE : ::LAI [m2 m{-2}]
1015	!!
1016	!! REFERENCE(S) :
1017	!!
1018	!! FLOWCHART : None
1019	!! \n
1020	!_ ================================================================================================================================
1021
1022	FUNCTION cc_to_lai_1d(npts,ccbio,ccind,pft)
1023
1024	!! 0. Variable and parameter declaration
1025
1026	!! 0.1 Input variables
1027	INTEGER(i_std), INTENT(in) :: npts !! Domain size (unitless)
1028	INTEGER(i_std), INTENT(in) :: pft !! number of PFT under study (unitless)
1029	REAL(r_std), DIMENSION(:,:,:,:) :: ccbio !! circ_class_biomass. Biomass of an individual
1030	!! tree @tex $(gC(N) ind^{-1})$ @endtex
1031	REAL(r_std), DIMENSION(:,:) :: ccind !! circ_class_n. Number of trees per circumference
1032	!! class
1033
1034	!! 0.2 Output variables
1035
1036	REAL(r_std), DIMENSION(npts) :: cc_to_lai_1d !! leaf area index @tex $(m^{2} m^{-2})$ @endtex
1037
1038	!! 0.3 Modified variables
1039
1040	!! 0.4 Local variables
1041	INTEGER(i_std) :: icir !! Index
1042	REAL(r_std), DIMENSION(npts) :: leaf_biomass !! Leaf mass for carbon calculated from
1043	!! circ_class_biomass and circ_class_n.
1044	!! @tex $(gC m^{-2})$ @endtex
1045	!_ ================================================================================================================================
1046
1047	!! 1. Calculate leaf biomass for a specific PFT
1048	cc_to_lai_1d(:) = zero
1049	leaf_biomass(:) = zero
1050
1051	! Note that when used in 2d ivm is a single PFT
1052	IF ( is_tree(pft) ) THEN
1053
1054	DO icir = 1,ncirc
1055
1056	leaf_biomass(:) = leaf_biomass(:) + &
1057	ccbio(:,icir,ileaf,icarbon) * ccind(:,icir)
1058
1059	ENDDO
1060
1061	ELSE
1062
1063	! Grasses and cropland only have one circumference class
1064
1065	leaf_biomass(:) = ccbio(:,1,ileaf,icarbon) * ccind(:,1)
1066
1067	ENDIF
1068
1069	!! 2. Calculate the LAI from the leaf biomass
1070	IF(sla_dyn) THEN
1071
1072	cc_to_lai_1d(:) = log(1.+ ext_coeff_N(pft)* leaf_biomass(:) * &
1073	slainit(pft))/(ext_coeff_N(pft))
1074
1075	ELSE
1076
1077	cc_to_lai_1d(:) = leaf_biomass(:) * sla(pft)
1078
1079	ENDIF
1080
1081	END FUNCTION cc_to_lai_1d
1082
1083
1084	!! ================================================================================================================================
1085	!! FUNCTION : wood_to_qmheight
1086	!!
1087	!>\BRIEF Calculate the quadratic mean height from the biomass
1088	!!
1089	!! DESCRIPTION : Calculates the quadratic mean height from the biomass
1090	!!
1091	!! RECENT CHANGE(S): None
1092	!!
1093	!! RETURN VALUE : ::qm_height (m)
1094	!!
1095	!! REFERENCE(S) :
1096	!!
1097	!! FLOWCHART : None
1098	!! \n
1099	!_ ================================================================================================================================
1100
1101	FUNCTION wood_to_qmheight(biomass_temp, ind, pft)
1102
1103	!! 0. Variable and parameter declaration
1104
1105	!! 0.1 Input variables
1106
1107	INTEGER(i_std) :: pft !! PFT number (-)
1108	REAL(r_std), DIMENSION(ncirc,nparts) :: biomass_temp !! Biomass of the leaves @tex $(gC m^{-2})$ @endtex
1109	REAL(r_std), DIMENSION(ncirc) :: ind !! Number of individuals @tex $(m^{-2})$ @endtex
1110
1111
1112	!! 0.2 Output variables
1113
1114	REAL(r_std) :: wood_to_qmheight !! quadratic mean height (m)
1115
1116	!! 0.3 Modified variables
1117
1118	!! 0.4 Local variables
1119	REAL(r_std), DIMENSION(ncirc) :: circ_class_ba !! basal area for each circ_class @tex $(m^{2})$ @endtex
1120	REAL(r_std) :: qm_dia !! quadratic mean diameter (m)
1121
1122	!_ ================================================================================================================================
1123
1124	!! 1. Calculate qm_height from the biomass
1125	IF ( is_tree(pft) ) THEN
1126
1127	! Basal area at the tree level (m2 tree-1)
1128	circ_class_ba(:) = wood_to_ba(biomass_temp(:,:),pft)
1129
1130	IF (SUM(ind(:)) .NE. zero) THEN
1131
1132	qm_dia = SQRT( 4/piSUM(circ_class_ba(:)ind(:))/SUM(ind(:)) )
1133
1134	ELSE
1135
1136	qm_dia = zero
1137
1138	ENDIF
1139
1140	wood_to_qmheight = pipe_tune2(pft)(qm_dia*pipe_tune3(pft))
1141
1142	ELSE
1143
1144	! Grasses and croplands
1145	! Calculate height as a function of the leaf biomass. Ensure that the
1146	! grasslands have a roughness length during the winter. Therefore use a
1147	! lower threshold.
1148	IF (SUM(ind(:)) .NE. zero) THEN
1149
1150	wood_to_qmheight = MAX(0.2,biomass_temp(1,ileaf)ind(1)sla(pft)*lai_to_height(pft))
1151	ELSE
1152
1153	wood_to_qmheight = zero
1154
1155	ENDIF
1156
1157	ENDIF ! is_tree(j)
1158
1159	END FUNCTION wood_to_qmheight
1160
1161
1162	!! ================================================================================================================================
1163	!! FUNCTION : wood_to_qmdia
1164	!!
1165	!>\BRIEF Calculate the quadratic mean diameter from the biomass
1166	!!
1167	!! DESCRIPTION : Calculates the quadratic mean diameter from the aboveground biomass
1168	!!
1169	!! RECENT CHANGE(S): None
1170	!!
1171	!! RETURN VALUE : ::qm_dia (m)
1172	!!
1173	!! REFERENCE(S) :
1174	!!
1175	!! FLOWCHART : None
1176	!! \n
1177	!_ ================================================================================================================================
1178
1179	FUNCTION wood_to_qmdia(biomass_temp, ind, pft)
1180
1181	!! 0. Variable and parameter declaration
1182
1183	!! 0.1 Input variables
1184
1185	INTEGER(i_std) :: pft !! PFT number (-)
1186	REAL(r_std), DIMENSION(ncirc,nparts) :: biomass_temp !! Biomass of the leaves @tex $(gC m^{-2})$ @endtex
1187	REAL(r_std), DIMENSION(ncirc) :: ind !! Number of individuals @tex $(m^{-2})$ @endtex
1188
1189	!! 0.2 Output variables
1190
1191	REAL(r_std) :: wood_to_qmdia !! quadratic mean diameter (m)
1192
1193	!! 0.3 Modified variables
1194
1195	!! 0.4 Local variables
1196	REAL(r_std), DIMENSION(ncirc) :: circ_class_ba !! basal area for each circ_class @tex $(m^{2})$ @endtex
1197
1198	!_ ================================================================================================================================
1199
1200	!! 1. Calculate qm_dia from the biomass
1201	IF ( is_tree(pft) ) THEN
1202
1203	! Basal area at the tree level (m2 tree-1)
1204	circ_class_ba(:) = wood_to_ba(biomass_temp(:,:),pft)
1205
1206	IF (SUM(ind(:)) .NE. zero) THEN
1207
1208	wood_to_qmdia = SQRT( 4/piSUM(circ_class_ba(:)ind(:))/SUM(ind(:)) )
1209
1210	ELSE
1211
1212	wood_to_qmdia = zero
1213
1214	ENDIF
1215
1216
1217	! Grasses and croplands
1218	ELSE
1219
1220	wood_to_qmdia = zero
1221
1222	ENDIF ! is_tree(pft)
1223
1224	END FUNCTION wood_to_qmdia
1225
1226
1227	!! ================================================================================================================================
1228	!! FUNCTION : wood_to_volume
1229	!!
1230	!>\BRIEF This allometric function computes volume as a function of
1231	!! biomass at stand scale. Volume \f$(m^3 m^{-2}) = f(biomass (gC m^{-2}))\f$
1232	!!
1233	!! DESCRIPTION : None
1234	!!
1235	!! RECENT CHANGE(S): None
1236	!!
1237	!! RETURN VALUE : biomass_to_volume
1238	!!
1239	!! REFERENCE(S) : See above, module description.
1240	!!
1241	!! FLOWCHART : None
1242	!! \n
1243	!_ ================================================================================================================================
1244
1245	FUNCTION wood_to_volume(npts,ccbio,ccind,pft,branch_ratio,inc_branches)
1246
1247	!! 0. Variable and parameter declaration
1248
1249	!! 0.1 Input variables
1250	INTEGER(i_std), INTENT(in) :: npts !! Domain size (unitless)
1251	REAL(r_std), INTENT(in), DIMENSION(:,:,:) :: ccbio !! circ_class_biomass or biomass of an
1252	!! individual tree in each circumference
1253	!! class @tex $(gC tree^{-1})$ @endtex
1254	REAL(r_std), INTENT(in), DIMENSION(:) :: ccind !! Number of trees per circumference class
1255	REAL(r_std), INTENT(in) :: branch_ratio !! Branch ratio of sap and heartwood biomass
1256	!! unitless
1257	INTEGER(i_std), INTENT(in) :: pft !! Plant functional type (unitless)
1258	INTEGER(i_std), INTENT(in) :: inc_branches !! Include the branches in the volume calculation?
1259	!! 0: exclude the branches from the volume calculation
1260	!! (thus correct the biomass for the branch ratio)
1261	!! 1: include the branches in the volume calculation
1262	!! (thus use all aboveground biomass)
1263
1264
1265
1266	!! 0.2 Output variables
1267
1268	REAL(r_std) :: wood_to_volume !! The volume of wood per square meter
1269	!! @tex $(m^3 m^{-2})$ @endtex
1270
1271	!! 0.3 Modified variables
1272
1273	!! 0.4 Local variables
1274	REAL(r_std), DIMENSION(nparts,nelements) :: biomass !! Total biomass at the stand level
1275	!! @tex $(gC m^{-2})$ @endtex
1276	REAL(r_std) :: woody_biomass !! Woody biomass at the stand level
1277	!! @tex $(gC m^{-2})$ @endtex
1278
1279	!_ ================================================================================================================================
1280
1281	!! 1. Volume to biomass
1282
1283	! Calculate stand biomass from circ_class_biomass. In cc_to_biomass_2d
1284	! which is used here, npts is not used but we need to pass it to have
1285	! a consistent call structure for cc_to_biomass_2d and cc_to_biomass_4d
1286	biomass = cc_to_biomass(npts,pft,ccbio,ccind)
1287
1288	! Woody biomass used in the calculation
1289	IF (inc_branches .EQ. 0) THEN
1290
1291	woody_biomass = (biomass(isapabove,icarbon)+biomass(iheartabove,icarbon)) * &
1292	(un - branch_ratio)
1293
1294	ELSEIF (inc_branches .EQ. 1) THEN
1295
1296	woody_biomass = (biomass(isapabove,icarbon)+biomass(iheartabove,icarbon))
1297
1298	ELSE
1299
1300	CALL ipslerr_p(3,'ERROR in wood_to_volume',&
1301	'Specify how the branches should be accounted for','','')
1302
1303	ENDIF
1304
1305	! Wood volume expressed in m3 / m2
1306	wood_to_volume = woody_biomass/(pipe_density(pft))
1307
1308	END FUNCTION wood_to_volume
1309
1310
1311	!! ================================================================================================================================
1312	!! FUNCTION : wood_to_height
1313	!!
1314	!>\BRIEF Calculate tree height from woody biomass making use of allometric relationships
1315	!!
1316	!! DESCRIPTION : Calculate height of an individual tree from the woody biomass of that tree making
1317	!! use of allometric relationships. This is the height used in forestry and for calculating the aerodynamic
1318	!! interactions
1319	!! (i) height(:) = pipe_tune2(j)(4/piba(:))**(pipe_tune3(j)/2) and
1320	!! (ii) woodmass_ind = tree_ffpipe_densityba*height
1321	!!
1322	!! RECENT CHANGE(S): None
1323	!!
1324	!! RETURN VALUE : height (m)
1325	!!
1326	!! REFERENCE(S) :
1327	!!
1328	!! FLOWCHART : None
1329	!! \n
1330	!_ ================================================================================================================================
1331
1332	FUNCTION wood_to_height(biomass_temp, pft)
1333
1334	!! 0. Variable and parameter declaration
1335
1336	!! 0.1 Input variables
1337
1338	INTEGER(i_std) :: pft !! PFT number (-)
1339	REAL(r_std), DIMENSION(:,:) :: biomass_temp !! Biomass of an individual tree within a circ
1340	!! class @tex $(m^{2} ind^{-1})$ @endtex
1341
1342	!! 0.2 Output variables
1343
1344	REAL(r_std), DIMENSION(ncirc) :: wood_to_height !! Height of an individual tree within a circ
1345	!! class (m)
1346
1347	!! 0.3 Modified variables
1348
1349	!! 0.4 Local variables
1350	INTEGER(i_std) :: l !! Index
1351	REAL(r_std) :: woodmass_ind !! Woodmass of an individual tree
1352	!! @tex $(gC ind{-1})$ @endtex
1353	REAL(r_std) :: temp !! Temporary variable
1354	!_ ================================================================================================================================
1355
1356	!! 1. Calculate height from woodmass
1357
1358	IF ( is_tree(pft) ) THEN
1359
1360	DO l = 1,ncirc
1361
1362	! Woodmass of an individual tree (only the aboveground component)
1363	woodmass_ind = biomass_temp(l,isapabove) + biomass_temp(l,iheartabove)
1364
1365	temp = (woodmass_ind/(tree_ff(pft)pipe_density(pft))4/pi)
1366
1367	IF (temp.LT.EPSILON(un)) THEN
1368	! Expect an overflow when taking an exponent of such a small number
1369	wood_to_height(l) = min_stomate
1370	ELSE
1371	! Height of that individual
1372	wood_to_height(l) = (((woodmass_ind/(tree_ff(pft)*&
1373	pipe_density(pft))4/pi)(pipe_tune3(pft)/2))&
1374	pipe_tune2(pft))**(1/(pipe_tune3(pft)/2+1))
1375	ENDIF
1376
1377	ENDDO
1378
1379	ELSE
1380
1381	WRITE(numout,*) 'ERROR: wood_to_height, pft ',pft
1382	CALL ipslerr_p (3,'wood_to_height', &
1383	'wood_to_height is not defined for this PFT.', &
1384	'See the output file for more details.','')
1385
1386	ENDIF
1387
1388	END FUNCTION wood_to_height
1389
1390
1391	!! ================================================================================================================================
1392	!! FUNCTION : wood_to_height_eff
1393	!!
1394	!>\BRIEF Calculate the effective tree height from woody biomass making use of allometric relationships
1395	!!
1396	!! DESCRIPTION : Calculate the effective height of an individual tree from the woody biomass of that tree making
1397	!! use of allometric relationships. Effective height makes use of both above and belowground biomass and is
1398	!! used in the calculation of the allocation according to deleuze and dhote, hydraulic architecture because also
1399	!! the height of the belowground part should be included.
1400	!! (i) height(:) = pipe_tune2(j)(4/piba(:))**(pipe_tune3(j)/2) and
1401	!! (ii) woodmass_ind = tree_ffpipe_densityba*height
1402	!!
1403	!! RECENT CHANGE(S): None
1404	!!
1405	!! RETURN VALUE : height (m)
1406	!!
1407	!! REFERENCE(S) :
1408	!!
1409	!! FLOWCHART : None
1410	!! \n
1411	!_ ================================================================================================================================
1412
1413	FUNCTION wood_to_height_eff(biomass_temp, pft)
1414
1415	!! 0. Variable and parameter declaration
1416
1417	!! 0.1 Input variables
1418
1419	INTEGER(i_std) :: pft !! PFT number (-)
1420	REAL(r_std), DIMENSION(:,:) :: biomass_temp !! Biomass of an individual tree within a circ
1421	!! class @tex $(m^{2} ind^{-1})$ @endtex
1422
1423	!! 0.2 Output variables
1424
1425	REAL(r_std), DIMENSION(ncirc) :: wood_to_height_eff !! Effective height of an individual tree within a circ
1426	!! class (m)
1427
1428	!! 0.3 Modified variables
1429
1430	!! 0.4 Local variables
1431	INTEGER(i_std) :: l !! Index
1432	REAL(r_std) :: woodmass_ind !! Woodmass of an individual tree
1433	!! @tex $(gC ind{-1})$ @endtex
1434	REAL(r_std) :: temp !! Temporary variable
1435	!_ ================================================================================================================================
1436
1437	!! 1. Calculate height from woodmass
1438
1439	IF ( is_tree(pft) ) THEN
1440
1441	DO l = 1,ncirc
1442
1443	! Woodmass of an individual tree. Both above and belowground biomass are used
1444	woodmass_ind = biomass_temp(l,isapabove) + biomass_temp(l,isapbelow) + &
1445	biomass_temp(l,iheartabove) + biomass_temp(l,iheartbelow)
1446
1447	temp = (woodmass_ind/(tree_ff(pft)pipe_density(pft))4/pi)
1448
1449	IF (temp.LT.EPSILON(un)) THEN
1450	! Expect an overflow when taking an exponent of such a small number
1451	wood_to_height_eff(l) = min_stomate
1452	ELSE
1453	! Height of that individual
1454	wood_to_height_eff(l) = (((woodmass_ind/(tree_ff(pft)pipe_density(pft))4/pi)**(pipe_tune3(pft)/2))&
1455	pipe_tune2(pft))*(1/(pipe_tune3(pft)/2+1))
1456	ENDIF
1457
1458	ENDDO
1459
1460	ELSE
1461
1462	WRITE(numout,*) 'pft ',pft
1463	CALL ipslerr_p (3,'wood_to_height_eff', &
1464	'wood_to_height_eff is not defined for this PFT.', &
1465	'See the output file for more details.','')
1466
1467	ENDIF
1468
1469	END FUNCTION wood_to_height_eff
1470
1471
1472	!! ================================================================================================================================
1473	!! FUNCTION : wood_to_dia
1474	!!
1475	!>\BRIEF Calculate diameter from woody biomass making use of allometric relationships
1476	!!
1477	!! DESCRIPTION : Calculate diameter of an individual tree from the woody biomass of that tree making
1478	!! use of allometric relationships. Makes only use of the aboveground biomass and relates to the
1479	!! typical forestry diameter (but not normalized at 1.3 m)
1480	!! (i) woodmass_ind = tree_ff * pipe_density * height * pi/4dia*2
1481	!! (ii) height = pipe_tune2 * dia * pipe_tune3
1482	!!
1483	!! RECENT CHANGE(S): None
1484	!!
1485	!! RETURN VALUE : diameter (m)
1486	!!
1487	!! REFERENCE(S) :
1488	!!
1489	!! FLOWCHART : None
1490	!! \n
1491	!_ ================================================================================================================================
1492
1493	FUNCTION wood_to_dia(biomass_temp, pft)
1494
1495
1496	!! 0. Variable and parameter declaration
1497
1498
1499	!! 0.1 Input variables
1500
1501	INTEGER(i_std) :: pft !! PFT number (-)
1502	REAL(r_std), DIMENSION(:,:) :: biomass_temp !! Biomass of an individual tree within a circ
1503	!! class @tex $(m^{2} ind^{-1})$ @endtex
1504	!! 0.2 Output variables
1505
1506	REAL(r_std), DIMENSION(ncirc) :: wood_to_dia !! Diameter of an individual tree within a circ
1507	!! class (m)
1508
1509	!! 0.3 Modified variables
1510
1511	!! 0.4 Local variables
1512	INTEGER(i_std) :: l !! Index
1513	REAL(r_std) :: woodmass_ind !! Woodmass of an individual tree
1514	!! @tex $(gC ind^{-1})$ @endtex
1515	REAL(r_std) :: temp !! temporary variable
1516	!_ ================================================================================================================================
1517
1518	!! 1. Calculate basal area from woodmass
1519
1520	IF ( is_tree(pft) ) THEN
1521
1522	DO l = 1,ncirc
1523
1524	! Woodmass of an individual tree
1525	woodmass_ind = biomass_temp(l,isapabove) + biomass_temp(l,iheartabove)
1526
1527	temp = 4/piwoodmass_ind/(tree_ff(pft)pipe_density(pft)*pipe_tune2(pft))
1528
1529	IF (temp.LT.EPSILON(un)) THEN
1530	! Expect an overflow when taking an exponent of such a small number
1531	wood_to_dia(l) = min_stomate
1532	ELSE
1533	! Basal area of that individual (m2 ind-1)
1534	wood_to_dia(l) = (4/piwoodmass_ind/(tree_ff(pft)pipe_density(pft)pipe_tune2(pft))) * &
1535	(1./(2+pipe_tune3(pft)))
1536	ENDIF
1537	ENDDO
1538
1539	ELSE
1540
1541	WRITE(numout,*) 'pft ',pft
1542	CALL ipslerr_p (3,'wood_to_dia', &
1543	'wood_to_dia is not defined for this PFT.', &
1544	'See the output file for more details.','')
1545
1546	ENDIF
1547
1548	END FUNCTION wood_to_dia
1549
1550
1551	!! ================================================================================================================================
1552	!! FUNCTION : wood_to_circ
1553	!!
1554	!>\BRIEF Calculate circumference from woody biomass making use of allometric relationships
1555	!!
1556	!! DESCRIPTION : All this does it computer the diameter using a different routine, and then
1557	!! convert that into a circumference.
1558	!!
1559	!! RECENT CHANGE(S): None
1560	!!
1561	!! RETURN VALUE : circumference (m)
1562	!!
1563	!! REFERENCE(S) :
1564	!!
1565	!! FLOWCHART : None
1566	!! \n
1567	!_ ================================================================================================================================
1568
1569	FUNCTION wood_to_circ(biomass_temp, pft)
1570
1571	!! 0. Variable and parameter declaration
1572
1573	!! 0.1 Input variables
1574
1575	INTEGER(i_std) :: pft !! PFT number (-)
1576	REAL(r_std), DIMENSION(:,:) :: biomass_temp !! Biomass of an individual tree within a circ
1577	!! class @tex $(m^{2} ind^{-1})$ @endtex
1578
1579	!! 0.2 Output variables
1580
1581	REAL(r_std), DIMENSION(ncirc) :: wood_to_circ !! Circumference of an individual tree within a circ
1582	!! class (m)
1583
1584	!! 0.3 Modified variables
1585
1586	!! 0.4 Local variables
1587
1588	!_ ================================================================================================================================
1589
1590	!! 1. Calculate diameter from woodmass
1591
1592	wood_to_circ(:)=val_exp
1593
1594	wood_to_circ(:)=wood_to_dia(biomass_temp(:,:),pft)
1595
1596	! convert to a circumference (m)
1597	wood_to_circ(:) = wood_to_circ(:)*pi
1598
1599	END FUNCTION wood_to_circ
1600
1601
1602	!! ================================================================================================================================
1603	! SUBROUTINE : wood_to_cn_dia
1604	!!
1605	!>\BRIEF Calculate horizontal crown diameter from woody biomass making use of allometric relationships
1606	!!
1607	!! DESCRIPTION : Calculate tree height and apply the reduction coefficient crown_to_height and crown_vertohor_dia
1608	!!
1609	!! RECENT CHANGE(S): None
1610	!!
1611	!! RETURN VALUE : horizontal crown diameter (m)
1612	!!
1613	!! REFERENCE(S) :
1614	!!
1615	!! FLOWCHART : None
1616	!! \n
1617	!_ ================================================================================================================================
1618	SUBROUTINE wood_to_cn_dia(biomass_temp, ind, pft, dia_hor, dia_ver)
1619
1620	!! 0. Variable and parameter declaration
1621
1622	!! 0.1 Input variables
1623
1624	INTEGER(i_std) :: pft !! PFT number (-)
1625	REAL(r_std), DIMENSION(:,:) :: biomass_temp !! Biomass of an individual tree within a circ
1626	!! class @tex $(m^{2} ind^{-1})$ @endtex
1627	REAL(r_std), DIMENSION(:) :: ind !! Number of individuals within a circ class
1628	!! @tex $(ind m^{2})$ @endtex
1629
1630	!! 0.2 Output variables
1631
1632	REAL(r_std), DIMENSION(ncirc) :: dia_hor !! Horizontal crown diameter
1633	!! for each individual tree within a circ
1634	!! class (m)
1635	REAL(r_std), DIMENSION(ncirc) :: dia_ver !! Vertical crown diameter (in othe words, crown
1636	!! depth) for each individual tree within a circ
1637	!! class (m)
1638
1639	!! 0.3 Modified variables
1640
1641	!! 0.4 Local variables
1642	INTEGER(i_std) :: icir !! Index
1643	REAL(r_std) :: woodmass_ind !! Woodmass of an individual tree
1644	!! @tex $(gC ind{-1})$ @endtex
1645	REAL(r_std) :: temp !! Temporary variable
1646	REAL(r_std) :: cn_vol !! Crown volume of a individual tree within a circ class (m)
1647	REAL(r_std) :: total_cn_vol !! Total crown volume of the stand across all circ classes (m)
1648	REAL(r_std) :: max_height !! Volume of the entire canopy space (surface area * tree height
1649	!! but normalized for surface area) (m)
1650	REAL(r_std) :: adjusted_crown_vertohor_dia !! calculated ratio between vertical and horizontal
1651	!! crown diameter. The calculated value guarantees that the
1652	!! canopys space is filled for 74% (given as ::crown_packing)
1653	REAL(r_std) :: adjusted_crown_to_height !! Calculate ratio between tree height and vertical
1654	!! crown diameter. The calculated value guarantees that the
1655	!! canopys space is filled for 74% (given as ::crown_packing)
1656
1657	!_ ================================================================================================================================
1658
1659	!! 1. Calculate horizontal crown diameter from woodmass
1660
1661	IF ( is_tree(pft) ) THEN
1662
1663	! Initialize
1664	total_cn_vol=zero
1665
1666	! Use wood mass to calculate the prescribed vertical and horizontal crown
1667	! diameters
1668	DO icir = 1,ncirc
1669
1670	! Woodmass of an individual tree (only the aboveground component)
1671	woodmass_ind = biomass_temp(icir,isapabove) + biomass_temp(icir,iheartabove)
1672
1673	temp = (woodmass_ind/(tree_ff(pft)pipe_density(pft))4/pi)
1674
1675	IF (temp.LT.EPSILON(un)) THEN
1676
1677	! Expect an overflow when taking an exponent of such a small number
1678	dia_hor(icir) = min_stomate
1679	dia_ver(icir) = min_stomate
1680
1681	ELSE
1682
1683	! Horizonatal crown diameter of an individual based on the height
1684	! and prescribed ratios.
1685	dia_hor(icir) = crown_vertohor_dia(pft) * crown_to_height(pft) * &
1686	(((woodmass_ind/(tree_ff(pft)*&
1687	pipe_density(pft))4/pi)(pipe_tune3(pft)/2))&
1688	pipe_tune2(pft))**(1/(pipe_tune3(pft)/2+1))
1689	! Vertical crown diameter derived from the height of that individual
1690	dia_ver(icir) = crown_to_height(pft)(((woodmass_ind/(tree_ff(pft)&
1691	pipe_density(pft))4/pi)(pipe_tune3(pft)/2))&
1692	pipe_tune2(pft))**(1/(pipe_tune3(pft)/2+1))
1693
1694	ENDIF
1695
1696	! Calculate the crown volume making use of the prescribed crown diameters
1697	! Notice that we don't have to take the surface area of the stand
1698	! into account. total_bm will be in units of gC / m**2, while
1699	! the canopy volume should be in units of m3 / m2. The surface
1700	! area would cancel out when we take the ratio. The unit of total_cn_vol
1701	! is thus similar to how precipitation is being expressed.
1702	! Right now, total_cn_vol is in m*3/tree tree/m**2 thus m
1703	cn_vol = pi/6dia_hor(icir)dia_hor(icir)*dia_ver(icir)
1704	total_cn_vol=total_cn_vol+cn_vol*ind(icir)
1705	ENDDO
1706
1707	! Space available to be filled with the crown volumes. The same
1708	! unit convention as above is used (m)
1709	max_height = MAXVAL(wood_to_height(biomass_temp(:,:),pft))
1710
1711
1712	IF (total_cn_vol.GT.max_height) THEN
1713
1714	! The prescribed crown diameters take too much space. Adjust the
1715	! horizontal diameters such that the space is more correctly filled
1716	! Given that we have to deal with different tree sizes and that the
1717	! canopies are assumed to be ellipsoids this is a very complex problem
1718	! to calculate exactly. We will make a couple of shortcuts.
1719	! First we will assume that by packing ellipsoids we have the same
1720	! efficiency as close-packing of equal spheres which is around 0.74.
1721	! total_cn_vol = 1/6pidia_hor^2*dia_ver
1722	! dia_hor = crown_vertohor_dia*dia_ver
1723	! <=> total_cn_vol = 1/6picrown_vertohor_dia^2*dia_ver^3
1724	! Reduce total_cn_vol to max_height * 0.74 (close packing) and calculate
1725	! the adjusted crown_vertohor_dia.
1726	! As an alternative scale both the vertyical and horizontal crown dimensions
1727	adjusted_crown_to_height = &
1728	((max_heightcrown_packing6) / &
1729	(pi(crown_vertohor_dia(pft)2)SUM(ind(:)dia_ver(:)3)))*(1./3.)
1730	dia_ver(:) = adjusted_crown_to_height * crown_to_height(pft) * &
1731	wood_to_height(biomass_temp(:,:),pft)
1732	dia_hor(:) = crown_vertohor_dia(pft) * dia_ver(:)
1733
1734	! Error checking
1735	total_cn_vol = zero
1736	DO icir = 1,ncirc
1737	total_cn_vol = total_cn_vol + &
1738	pi/6dia_hor(icir)dia_hor(icir)dia_ver(icir)ind(icir)
1739	END DO
1740	IF (ABS(total_cn_vol-crown_packing*max_height).GT.min_stomate) THEN
1741	WRITE(numout,*) 'ind, ', ind(:)
1742	WRITE(numout,*) 'dia_ver, ', dia_ver(:)
1743	WRITE(numout,*) 'dia_hor, ', dia_hor(:)
1744	WRITE(numout,*) 'adjusted_crown_vertohor_dia, ',adjusted_crown_vertohor_dia
1745	WRITE(numout,*) 'adjusted_crown_to_height, ',adjusted_crown_to_height
1746	WRITE(numout,*) 'recalculated total_cn_vol, ',total_cn_vol
1747	WRITE(numout,) 'crown_packingmax_height, ', crown_packing*max_height
1748	CALL ipslerr_p(3,'problem with crown packing',&
1749	'adjusted horizonatal diameter is wrong',&
1750	'check the calculations in the code','')
1751	END IF
1752
1753	END IF
1754
1755	ELSE
1756
1757	WRITE(numout,*) 'ERROR: wood_to_height, pft ',pft
1758	CALL ipslerr_p (3,'wood_to_height', &
1759	'wood_to_height is not defined for this PFT.', &
1760	'See the output file for more details.','')
1761
1762	ENDIF
1763
1764	END SUBROUTINE wood_to_cn_dia
1765
1766
1767	!! ================================================================================================================================
1768	!! FUNCTION : wood_to_cv
1769	!!
1770	!>\BRIEF Calculate crown volume from woody biomass making use of the vertical and horizontal crown diameters
1771	!!
1772	!! DESCRIPTION : Calculate the volume of an elipsoid. V = 1/6 * pi * dia_hor*2 dia_ver
1773	!!
1774	!! RECENT CHANGE(S): None
1775	!!
1776	!! RETURN VALUE : crown volume (m3 ind-1)
1777	!!
1778	!! REFERENCE(S) :
1779	!!
1780	!! FLOWCHART : None
1781	!! \n
1782	!_ ================================================================================================================================
1783
1784	FUNCTION wood_to_cv(biomass_temp, ind, pft)
1785
1786	!! 0. Variable and parameter declaration
1787
1788	!! 0.1 Input variables
1789
1790	INTEGER(i_std) :: pft !! PFT number (-)
1791	REAL(r_std), DIMENSION(:,:) :: biomass_temp !! Biomass of an individual tree within a circ
1792	!! class @tex $(m^{2} ind^{-1})$ @endtex
1793	REAL(r_std), DIMENSION(: ) :: ind !! Number of individuals within a circ class
1794	!! @tex $(ind m^{2})$ @endtex
1795
1796	!! 0.2 Output variables
1797
1798	REAL(r_std), DIMENSION(ncirc) :: wood_to_cv !! Crown volume of an individual tree
1799	!! @tex $(m^{2} ind{-1})$ @endtex
1800
1801	!! 0.3 Modified variables
1802
1803	!! 0.4 Local variables
1804	INTEGER(i_std) :: l !! Index
1805	REAL(r_std) :: woodmass_ind !! Woodmass of an individual tree
1806	!! @tex $(gC ind^{-1})$ @endtex
1807	REAL(r_std), DIMENSION(ncirc) :: dia_ver !! vertical crown diameter of an individual tree (m)
1808	REAL(r_std), DIMENSION(ncirc) :: dia_hor !! horizontal crown diameter of an individual tree (m)
1809
1810	!_ ================================================================================================================================
1811
1812	!! 1. Calculate crown volume from woodmass
1813
1814	IF ( is_tree(pft) ) THEN
1815
1816	! Crown diameter of the individual tree (m2 ind-1)
1817	CALL wood_to_cn_dia(biomass_temp, ind, pft, dia_hor, dia_ver)
1818	wood_to_cv(:) = pi/6.dia_ver(:)dia_hor(:)*dia_hor(:)
1819	! WRITE(numout,*) 'wood_to_cv_eff, ', wood_to_cv_eff(:)
1820
1821	ELSE
1822
1823	WRITE(numout,*) 'pft ',pft
1824	CALL ipslerr_p (3,'wood_to_cv_eff', &
1825	'wood_to_cv_eff is not defined for this PFT.', &
1826	'See the output file for more details.','')
1827
1828	ENDIF
1829
1830	END FUNCTION wood_to_cv
1831
1832
1833	!! ================================================================================================================================
1834	!! FUNCTION : Nmax
1835	!!
1836	!>\BRIEF This function determines the maximum number of trees per hectare for
1837	!! a given quadratic mean diameter (Dg). It applies the self-thinning principle
1838	!! of Reineke (1933), with Dg instead of mean diameter (Dhote 1999).
1839	!! Parameterization: Dhote (1999) for broad-leaved and Vacchiano (2008)
1840	!! for needle-leaved.
1841	!!
1842	!! DESCRIPTION : None
1843	!!
1844	!! RECENT CHANGE(S): None
1845	!!
1846	!! RETURN VALUE : Nmax
1847	!!
1848	!! REFERENCE(S) : See above, module description.
1849	!!
1850	!! FLOWCHART : None
1851	!! \n
1852	!_ ================================================================================================================================
1853
1854	FUNCTION Nmax(Dg,no_pft)
1855
1856	!! 0. Variable and parameter declaration
1857
1858	!! 0.1 Input variables
1859
1860	REAL(r_std) :: Dg !! Quadratic mean diameter (cm)
1861	INTEGER(i_std) :: no_pft !! Plant functional type (unitless)
1862
1863	!! 0.2 Output variables
1864
1865	REAL(r_std) :: Nmax !! Maximum number of trees according to the self-thinning model in (trees ha-1)
1866
1867	!! 0.3 Modified variables
1868
1869	!! 0.4 Local variables
1870
1871	!_ ================================================================================================================================
1872	!! 1. maximum number of trees per hectare for a given quadratic mean diameter
1873
1874	IF (is_tree(no_pft)) THEN
1875
1876	! thinning curve of the MTC
1877	Nmax = (Dg/alpha_self_thinning(no_pft))**(un/beta_self_thinning(no_pft))
1878
1879	! Truncate the range to avoid huge numbers due the exponental model used to describe self-thinning
1880	! Don't set to ntreesmax else stomate_mark_kill will assume we are at self-thinning and
1881	! won't apply background mortality.
1882	! Nmax = MIN( Nmax, ntreesmax(no_pft))
1883	Nmax = MAX( Nmax, dens_target(no_pft) )
1884
1885	ELSE
1886
1887	WRITE(numout,*) 'Self thinning is not defined for PFT, ', no_pft
1888	CALL ipslerr_p (3,'nmax', &
1889	'Self thinning is not defined for this PFT.', &
1890	'See the output file for more details.','')
1891
1892	ENDIF
1893
1894	END FUNCTION Nmax
1895
1896	!! ================================================================================================================================
1897	!! FUNCTION : calculate_rdi
1898	!!
1899	!>\BRIEF Calculate the relative density index of a stand
1900	!!
1901	!! DESCRIPTION : None
1902	!!
1903	!! RECENT CHANGE(S): None
1904	!!
1905	!! RETURN VALUE : calculate_rdi
1906	!!
1907	!! REFERENCE(S) : Bellassen, V., Le Maire, G., Dhote, J.F., Viovy, N., Ciais, P., 2010.
1908	!! Modeling forest management within a global vegetation model â Part 1:
1909	!! model structure and general behaviour. Ecological Modelling 221, 2458â2474.
1910	!!
1911	!! FLOWCHART : None
1912	!! \n
1913	!_ ================================================================================================================================
1914
1915	FUNCTION calculate_rdi(diameters_temp,circ_class_n_temp,ivm)
1916
1917	!! 0. Variable and parameter declaration
1918
1919	IMPLICIT NONE
1920
1921	!! 0.1 Input variables
1922	REAL(r_std),DIMENSION(ncirc),INTENT(in) :: circ_class_n_temp !! The distribution of individuals
1923	!! for a given grid square and PFT
1924	REAL(r_std),DIMENSION(ncirc),INTENT(in) :: diameters_temp !! The diameters of the trunks
1925	!! for each circ class (m)
1926	INTEGER(i_std),INTENT(IN) :: ivm !! The number of the PFT
1927	!! (unitless)
1928
1929	!! 0.2 Output variables
1930	REAL(r_std) :: calculate_rdi !! The RDI of the stand
1931	!! (unitless, always positive)
1932
1933	!! 0.3 Modified variables
1934
1935	!! 0.4 Local variables
1936	REAL(r_std) :: Dg !! The quadratic mean diameter
1937	REAL(r_std) :: ind_loc
1938	INTEGER(i_std) :: icir !! Indices
1939
1940	!_ ================================================================================================================================
1941
1942	!! 1. Calculate the relative density index for a stand, based on
1943	!! the number of trees in each circumference class and the
1944	!! the diameter of a model tree in each class.
1945	ind_loc=SUM(circ_class_n_temp(:))
1946
1947	Dg=zero
1948
1949	DO icir=1,ncirc
1950	Dg=Dg+circ_class_n_temp(icir)m2_to_ha(diameters_temp(icir)m_to_cm)*2
1951	ENDDO
1952	Dg=SQRT(Dg/(ind_loc*m2_to_ha))
1953
1954	calculate_rdi=(ind_loc*m2_to_ha)/Nmax(Dg,ivm)
1955
1956	END FUNCTION calculate_rdi
1957
1958
1959	!! ================================================================================================================================
1960	!! SUBROUTINE : check_vegetation_area
1961	!!
1962	!>\BRIEF
1963	!!
1964	!! DESCRIPTION :
1965	!! RECENT CHANGE(S): None
1966	!!
1967	!! MAIN OUTPUT VARIABLE(S):
1968	!!
1969	!! REFERENCE(S) : None
1970	!!
1971	!! FLOWCHART : None
1972	!! \n
1973	!_ ================================================================================================================================
1974	SUBROUTINE check_vegetation_area(check_point, npts, veget_max_begin, &
1975	veget_max, test_type, change_nobio)
1976
1977	!! 0. Variable and parameter description
1978
1979	!! 0.1 Input variables
1980	INTEGER(i_std), INTENT(in) :: npts !! Domain size (unitless)
1981	CHARACTER(*),INTENT(in) :: check_point !! A flag to indicate at which
1982	!! point in the code we're doing
1983	!! this check
1984	CHARACTER(*), INTENT(in) :: test_type !! Which test do we want to run on veget_max
1985	REAL(r_std), DIMENSION(:,:), INTENT(in) :: veget_max !! "maximal" coverage fraction of a PFT
1986	!! (LAI -> infinity) on ground. May sum to
1987	!! less than unity if the pixel has
1988	!! nobio area. (unitless; 0-1)
1989	REAL(r_std), DIMENSION(:,:), INTENT(in) :: veget_max_begin !! temporary storage of veget_max to check area conservation
1990	REAL(r_std), DIMENSION(:), INTENT(in), OPTIONAL :: change_nobio !! Change in the nobio fraction due to urbanisation or possibly glacier retreat
1991
1992	!! 0.2 Output variables
1993
1994	!! 0.3 Modified variables
1995
1996	!! 0.4 Local variables
1997	INTEGER(i_std) :: iele, ipts, ivm !! Indices
1998	INTEGER(i_std) :: ipar, icir !! Indices
1999	INTEGER(i_std) :: ivma,iagec,ipft !! Counter
2000	REAL(r_std) :: temp_begin, temp_end !! temporary variables to accumulate veget_max of
2001	!! different age classes of the same species
2002
2003	!_ ================================================================================================================================
2004
2005	SELECT CASE (test_type)
2006
2007	CASE ('pft')
2008
2009	! In all subroutines where veget_max should be preserved at the
2010	! PFT level, the most strict test can be performed
2011	DO ipts=1,npts
2012
2013	DO ivm = 1,nvm
2014
2015	IF (veget_max_begin(ipts,ivm) - veget_max(ipts,ivm) .GT. &
2016	min_stomate ) THEN
2017	WRITE(numout,*) 'Surface area is not preserved in ', &
2018	TRIM(check_point)
2019	WRITE(numout,*) 'Area has changed for PFT, ',ipts, ivm
2020	WRITE(numout,*) 'Change from, ',veget_max_begin(ipts,ivm), &
2021	' to ', veget_max(ipts,ivm)
2022	IF(err_act.GT.1) CALL ipslerr_p (plev, TRIM(check_point), &
2023	'Surface area not preserved','','')
2024	ENDIF
2025
2026	ENDDO !nvm
2027
2028	! One message per pixel
2029	!!$ IF (err_act.GT.1) WRITE(numout,*) 'Surface area preserved in ', &
2030	!!$ TRIM(check_point), ipts
2031
2032	ENDDO ! npts
2033
2034	CASE ('ageclass')
2035
2036	! In subroutines where veget_max can move from one pft to
2037	! another within the same species. This is the test to use
2038	DO ipts=1,npts
2039
2040	DO ivma=1,nvmap
2041
2042	! get pft number
2043	ivm=start_index(ivma)
2044
2045	! If we only have a single age class for this
2046	! PFT, we can simplify the test
2047	IF(nagec_pft(ivma) .EQ. 1) THEN
2048
2049	IF (veget_max_begin(ipts,ivm) - veget_max(ipts,ivm) .GT. &
2050	min_stomate ) THEN
2051	WRITE(numout,*) 'Surface area is not preserved in ', &
2052	TRIM(check_point)
2053	WRITE(numout,*) 'Area has changed for PFT, ',ipts, ivm
2054	WRITE(numout,*) 'Change from, ',veget_max_begin(ipts,ivm), &
2055	' to ', veget_max(ipts,ivm)
2056	IF(err_act.GT.1) CALL ipslerr_p (plev, TRIM(check_point), &
2057	'Surface area not preserved','','')
2058	ENDIF
2059
2060	ELSE ! more than 1 age class
2061
2062	temp_begin = 0
2063	temp_end = 0
2064	DO iagec = nagec_pft(ivma),1,-1
2065
2066	! Accumulate veget_max for all pfts/age classes
2067	! that belong to the same species
2068	ipft = ivm+iagec-1
2069	temp_begin = temp_begin + veget_max_begin(ipts,ipft)
2070	temp_end = temp_end + veget_max(ipts,ipft)
2071
2072	ENDDO
2073
2074	IF (temp_begin - temp_end .GT. min_stomate ) THEN
2075	WRITE(numout,*) 'Surface area is not preserved in ', &
2076	TRIM(check_point)
2077	WRITE(numout,*) 'Area has changed for PFT, ',ipts, ivm
2078	WRITE(numout,*) 'Change from, ',veget_max_begin(ipts,ivm), &
2079	' to ', veget_max(ipts,ivm)
2080	IF(err_act.GT.1) CALL ipslerr_p (plev, TRIM(check_point), &
2081	'Surface area not preserved','','')
2082	ENDIF
2083
2084	ENDIF ! number of age classes
2085
2086	ENDDO ! nvmap
2087
2088	! One message per pixel
2089	!IF (err_act.GT.1) WRITE(numout,*) 'Surface area preserved in ', &
2090	! TRIM(check_point), ipts
2091
2092	ENDDO ! ipts
2093
2094	CASE('pixel')
2095
2096	! Use this test in subroutines where veget_max can be moved
2097	! between species. This is the least strict test and only
2098	! test wheter the total veget_max was preserved.
2099	DO ipts=1,npts
2100
2101	! Change_nobio is the change in the frac_nobio that may occur in
2102	! LCC. This change should be accounted for in the mass balance check.
2103	! ::change_nobio is defined as an optional element and is only present
2104	! when the surface area is checked in sapiens_lcchange
2105	IF(PRESENT(change_nobio))THEN
2106	IF (ABS(SUM(veget_max_begin(ipts,:)) - SUM(veget_max(ipts,:)) - change_nobio(ipts)) .GT. &
2107	min_stomate ) THEN
2108	WRITE(numout,*) 'Surface area is not preserved in ', &
2109	TRIM(check_point)
2110	WRITE(numout,*) 'Area in pixel ', ipts,' has changed form, ', &
2111	SUM(veget_max_begin(ipts,:)), 'to,', SUM(veget_max(ipts,:))
2112	IF(err_act.GT.1) CALL ipslerr_p (plev, TRIM(check_point), &
2113	'Surface area not preserved','','')
2114	ENDIF
2115	ELSE
2116	IF (ABS(SUM(veget_max_begin(ipts,:)) - SUM(veget_max(ipts,:))) .GT. &
2117	min_stomate ) THEN
2118	WRITE(numout,*) 'Surface area is not preserved in ', &
2119	TRIM(check_point)
2120	WRITE(numout,*) 'Area in pixel ', ipts,' has changed form, ', &
2121	SUM(veget_max_begin(ipts,:)), 'to,', SUM(veget_max(ipts,:))
2122	IF(err_act.GT.1) CALL ipslerr_p (plev, TRIM(check_point), &
2123	'Surface area not preserved','','')
2124	ENDIF
2125	ENDIF
2126
2127	ENDDO ! npts
2128
2129	CASE DEFAULT
2130
2131	WRITE(numout,*) 'Error: an unknown type/resolution was requested'
2132	WRITE(numout,*) ' for the surface area check. Check the cases'
2133	WRITE(numout,*) ' defined in src_parameters/function_library or'
2134	WRITE(numout,*) ' the spelling of the type in the CALL argument',&
2135	test_type
2136	CALL ipslerr_p (3,'surface_area_check in function_library',&
2137	'unknown type', 'Check *_out_orchidee for more details','')
2138
2139	END SELECT
2140
2141	END SUBROUTINE check_vegetation_area
2142
2143	!! ================================================================================================================================
2144	!! SUBROUTINE : check_mass balance
2145	!!
2146	!>\BRIEF
2147	!!
2148	!! DESCRIPTION :
2149	!! RECENT CHANGE(S): None
2150	!!
2151	!! MAIN OUTPUT VARIABLE(S):
2152	!!
2153	!! REFERENCE(S) : None
2154	!!
2155	!! FLOWCHART : None
2156	!! \n
2157	!_ ================================================================================================================================
2158	SUBROUTINE check_mass_balance(check_point, closure_intern, npts, pool_end, &
2159	pool_start, veget_max, test_type, pft)
2160
2161	!! 0. Variable and parameter description
2162
2163	!! 0.1 Input variables
2164	INTEGER(i_std), INTENT(in) :: npts !! Domain size (unitless)
2165	CHARACTER(*),INTENT(in) :: check_point !! A flag to indicate at which
2166	!! point in the code we're doing
2167	!! this check
2168	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: closure_intern !! Check closure of internal mass balance
2169	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
2170	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: pool_start !! Start and end pool of this routine
2171	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
2172	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: pool_end !! Start and end pool of this routine
2173	CHARACTER(*), INTENT(in) :: test_type !! Which test do we want to run on veget_max
2174	REAL(r_std), DIMENSION(:,:), INTENT(in) :: veget_max !! "maximal" coverage fraction of a PFT
2175	!! (LAI -> infinity) on ground. May sum to
2176	!! less than unity if the pixel has
2177	!! nobio area. (unitless; 0-1)
2178	INTEGER(i_std), OPTIONAL, INTENT(in) :: pft !! Number of the PFT for which the mass balance has to be checked
2179
2180	!! 0.2 Output variables
2181
2182	!! 0.3 Modified variables
2183
2184	!! 0.4 Local variables
2185	INTEGER :: iele, ipts, ivm !! Indices
2186	INTEGER :: ipar, icir, err !! Indices
2187	INTEGER(i_std) :: ivma, iagec, ipft !! Counter
2188	REAL(r_std) :: temp_begin, temp_end !! temporary variables to accumulate veget_max of
2189
2190	CHARACTER(LEN=10) :: var_name !! Name of element
2191
2192	!_ ================================================================================================================================
2193
2194	SELECT CASE (test_type)
2195
2196	CASE ('ipts')
2197
2198	! Close the mass balance for one specific pixel and pft.
2199	! This :test is used for subroutines that are called within
2200	! a DO ipts = 1,npts and DO ivm = 1,nvm loop. The mass
2201	! balance is checked for a single pft at a single pixel.
2202
2203	DO iele=1,nelements
2204
2205	IF (iele.EQ.icarbon) var_name = 'carbon'
2206	IF (iele.EQ.initrogen) var_name = 'nitrogen'
2207	err = zero
2208
2209	! Mass balance closure
2210	DO ipts=1,npts
2211	IF(ABS(SUM(closure_intern(ipts,:,iele))) .GT. min_stomate)THEN
2212
2213	! Count the number of errors
2214	err = err + 1
2215
2216	! Give extra details of where the mass balance
2217	! problem occurs. Stop the model if the user
2218	! asked for it
2219	IF(err_act.GT.1) THEN
2220	WRITE(numout,*) 'Error: mass balance is not closed in ', &
2221	TRIM(check_point), ' for ', var_name
2222	WRITE(numout,*) ' Problem occurs in pixel, ',ipts
2223	WRITE(numout,*) ' Difference is, ', &
2224	SUM(closure_intern(ipts,:,iele))
2225	WRITE(numout,*) ' pool_end,pool_start: ', &
2226	SUM(pool_end(ipts,:,iele)), &
2227	SUM(pool_start(ipts,:,iele))
2228
2229	! Exit or write warning depening on plev
2230	CALL ipslerr_p(plev, TRIM(check_point), &
2231	'Mass balance error','','')
2232	ENDIF
2233
2234	ENDIF
2235	END DO
2236
2237	!+++CHECK+++
2238	! A lot of output - only write a message if their is a problem
2239	!!$ ! One message per pixel. The err condition is needed because under
2240	!!$ ! err_act = 2 the model will write an error to the out_orchidee file
2241	!!$ ! but will not be stopped. Therefore it should be checked whether we
2242	!!$ ! have an error or not
2243	!!$ IF (err_act.GT.1 .AND. err.EQ.0) THEN
2244	!!$ WRITE(numout,*) 'Mass balance closure in ', &
2245	!!$ TRIM(check_point), ' for ',var_name
2246	!!$ ENDIF
2247	!++++++++++
2248
2249	ENDDO ! nelements
2250
2251	CASE ('pft')
2252
2253	! Close mass balance at the pixel and pft level.
2254	DO iele=1,nelements
2255
2256	IF (iele.EQ.icarbon) var_name = 'carbon'
2257	IF (iele.EQ.initrogen) var_name = 'nitrogen'
2258
2259	DO ipts=1,npts
2260
2261	err = zero
2262
2263	DO ivm=1,nvm
2264
2265	!check for NaN
2266	IF (isnan(pool_start(ipts,ivm,iele))) THEN
2267	WRITE(numout,*) 'Error: found an NaN in ', &
2268	TRIM(check_point), ' for ', var_name
2269	WRITE(numout,*)'At pixel',ipts,' for element',iele,&
2270	' and for PFT',ivm,' pool_start is NaN'
2271	WRITE(numout,*) ' pool_start: ', pool_start(ipts,:,iele)
2272	CALL ipslerr_p(plev,TRIM(check_point),'NaN problem','','')
2273	ENDIF
2274	IF (isnan(pool_end(ipts,ivm,iele))) THEN
2275	WRITE(numout,*) 'Error: found an NaN in ', &
2276	TRIM(check_point), ' for ', var_name
2277	WRITE(numout,*) 'At the pixel',ipts,' for element',iele,&
2278	' and for PFT',ivm,' pool_end is NaN'
2279	WRITE(numout,*) ' pool_end: ', pool_end(ipts,:,iele)
2280	CALL ipslerr_p(plev,TRIM(check_point),'NaN problem','','')
2281	ENDIF
2282	IF (isnan(closure_intern(ipts,ivm,iele))) THEN
2283	WRITE(numout,*) 'Error: found an NaN in ', &
2284	TRIM(check_point), ' for ', var_name
2285	WRITE(numout,*) 'At the pixel',ipts,' for element',iele,&
2286	' and for PFT',ivm,' closure_intern is NaN'
2287	WRITE(numout,*) ' closure_intern: ', closure_intern(ipts,:,iele)
2288	CALL ipslerr_p(plev,TRIM(check_point), 'NaN problem','','')
2289	ENDIF
2290
2291	! Mass balance closure
2292	IF(ABS(closure_intern(ipts,ivm,iele)) .GT. min_stomate)THEN
2293
2294	err = err+1
2295
2296	! If the model is under development or being tested we will
2297	! output extra information to facilitate debugging.
2298	IF(err_act.GT.1)THEN
2299	WRITE(numout,*) 'Error: mass balance is not closed in ', &
2300	TRIM(check_point), ' for ', var_name
2301	WRITE(numout,*) ' ipts, ivm; ', ipts,ivm
2302	WRITE(numout,*) ' Difference is, ', &
2303	closure_intern(ipts,ivm,iele)
2304	WRITE(numout,*) ' pool_end,pool_start: ', &
2305	pool_end(ipts,ivm,iele), &
2306	pool_start(ipts,ivm,iele)
2307
2308	! Exit or write warning depening on plev
2309	CALL ipslerr_p(plev, TRIM(check_point), &
2310	'Mass balance error','','')
2311	ENDIF
2312	ENDIF
2313	ENDDO ! nvm
2314
2315	!+++CHECK+++
2316	! A lot of output - only write a message if their is a problem
2317	!!$ ! One message per pixel
2318	!!$ IF (err_act.GT.1 .AND. err.EQ.zero) THEN
2319	!!$ WRITE(numout,*) 'Mass balance closure in ', &
2320	!!$ TRIM(check_point), ' for ',var_name
2321	!!$ ENDIF
2322	!++++++++++
2323
2324	ENDDO ! npts
2325
2326	ENDDO ! nelements
2327
2328	CASE ('ipft')
2329
2330	! Close mass balance at the pixel and pft level. Check for a
2331	! specific pft. This option is called from within a DO-loop.
2332	DO iele=1,nelements
2333
2334	IF (iele.EQ.icarbon) var_name = 'carbon'
2335	IF (iele.EQ.initrogen) var_name = 'nitrogen'
2336
2337	DO ipts=1,npts
2338
2339	err = zero
2340
2341	!check for NaN
2342	IF (isnan(pool_start(ipts,pft,iele))) THEN
2343	WRITE(numout,*) 'Error: found an NaN in ', &
2344	TRIM(check_point), ' for ', var_name
2345	WRITE(numout,*)'At the pixel',ipts,' for element',iele,&
2346	' and for PFT',pft,' pool_start is NaN'
2347	WRITE(numout,*) ' pool_start: ', pool_start(ipts,:,iele)
2348	CALL ipslerr_p(plev,TRIM(check_point),'NaN problem','','')
2349	ENDIF
2350	IF (isnan(pool_end(ipts,pft,iele))) THEN
2351	WRITE(numout,*) 'Error: found an NaN in ', &
2352	TRIM(check_point), ' for ', var_name
2353	WRITE(numout,*) 'At the pixel',ipts,' for element',iele,&
2354	' and for PFT',pft,' pool_end is NaN'
2355	WRITE(numout,*) ' pool_end: ', pool_end(ipts,:,iele)
2356	CALL ipslerr_p(plev,TRIM(check_point),'NaN problem','','')
2357	ENDIF
2358	IF (isnan(closure_intern(ipts,pft,iele))) THEN
2359	WRITE(numout,*) 'Error: found an NaN in ', &
2360	TRIM(check_point), ' for ', var_name
2361	WRITE(numout,*) 'At the pixel',ipts,' for element',iele,&
2362	' and for PFT',pft,' closure_intern is NaN'
2363	WRITE(numout,*) ' closure_intern: ', closure_intern(ipts,:,iele)
2364	CALL ipslerr_p(plev,TRIM(check_point), 'NaN problem','','')
2365	ENDIF
2366
2367
2368	! Mass balance closure
2369	IF(ABS(closure_intern(ipts,pft,iele)) .GT. min_stomate)THEN
2370
2371	err = err+1
2372
2373	! If the model is under development or being tested we will
2374	! output extra information to facilitate debugging.
2375	IF(err_act.GT.1)THEN
2376	WRITE(numout,*) 'Error: mass balance is not closed in ', &
2377	TRIM(check_point), ' for ', var_name
2378	WRITE(numout,*) ' ipts, pft; ', ipts,pft
2379	WRITE(numout,*) ' Difference is, ', &
2380	closure_intern(ipts,pft,iele)
2381	WRITE(numout,*) ' pool_end,pool_start: ', &
2382	pool_end(ipts,pft,iele), &
2383	pool_start(ipts,pft,iele)
2384
2385	! Exit or write warning depening on plev
2386	CALL ipslerr_p(plev, TRIM(check_point), &
2387	'Mass balance error','','')
2388	ENDIF
2389	ENDIF
2390
2391	ENDDO ! npts
2392
2393	ENDDO ! nelements
2394
2395	CASE ('ageclass')
2396
2397	! This option is only used when err_act.GE.2
2398
2399	! In the subroutine under test, biomass can move
2400	! from one age class to another age class in the
2401	! same species. Test across all age classes for
2402	! the same species.
2403	DO iele=1,nelements
2404
2405	IF (iele.EQ.icarbon) var_name = 'carbon'
2406	IF (iele.EQ.initrogen) var_name = 'nitrogen'
2407
2408	DO ipts=1,npts
2409
2410	err = zero
2411
2412	DO ivma=1,nvmap
2413
2414	! get pft number
2415	ivm=start_index(ivma)
2416
2417	!check for NaN
2418	IF (isnan(pool_start(ipts,ivm,iele))) THEN
2419	WRITE(numout,*) 'Error: found an NaN in ', &
2420	TRIM(check_point), ' for ', var_name
2421	WRITE(numout,*)'At the pixel',ipts,' for element',iele,&
2422	' and for PFT',ivm,' pool_start is NaN'
2423	WRITE(numout,*) ' pool_start: ', pool_start(ipts,:,iele)
2424	CALL ipslerr_p(plev,TRIM(check_point),'NaN problem','','')
2425	ENDIF
2426	IF (isnan(pool_end(ipts,ivm,iele))) THEN
2427	WRITE(numout,*) 'Error: found an NaN in ', &
2428	TRIM(check_point), ' for ', var_name
2429	WRITE(numout,*) 'At the pixel',ipts,' for element',iele,&
2430	' and for PFT',ivm,' pool_end is NaN'
2431	WRITE(numout,*) ' pool_end: ', pool_end(ipts,:,iele)
2432	CALL ipslerr_p(plev,TRIM(check_point),'NaN problem','','')
2433	ENDIF
2434	IF (isnan(closure_intern(ipts,ivm,iele))) THEN
2435	WRITE(numout,*) 'Error: found an NaN in ', &
2436	TRIM(check_point), ' for ', var_name
2437	WRITE(numout,*) 'At the pixel',ipts,' for element',iele,&
2438	' and for PFT',ivm,' closure_intern is NaN'
2439	WRITE(numout,*) ' closure_intern: ', closure_intern(ipts,:,iele)
2440	CALL ipslerr_p(plev,TRIM(check_point), 'NaN problem','','')
2441	ENDIF
2442
2443	! If we only have a single age class for this
2444	! PFT, we can simplify the test
2445	IF(nagec_pft(ivma) .EQ. 1) THEN
2446
2447	! Mass balance closure
2448	IF(ABS(closure_intern(ipts,ivm,iele)) .GT. min_stomate)THEN
2449	err=err+1
2450	WRITE(numout,*) 'Error: mass balance is not closed in ', &
2451	TRIM(check_point), ' for ', var_name
2452	WRITE(numout,*) ' ipts, ivm; ', ipts,ivm
2453	WRITE(numout,*) ' Difference is, ', &
2454	closure_intern(ipts,ivm,iele)
2455	WRITE(numout,*) ' pool_end,pool_start: ', &
2456	pool_end(ipts,ivm,iele), &
2457	pool_start(ipts,ivm,iele)
2458	IF (err_act.GT.1) CALL ipslerr_p (plev, TRIM(check_point), &
2459	'Mass balance error','','')
2460
2461	ENDIF ! test for closure
2462
2463	ELSE ! more than 1 age class
2464
2465	temp_begin = 0
2466	DO iagec = nagec_pft(ivma),1,-1
2467
2468	! Accumulate check_intern for all pfts/age classes
2469	! that belong to the same species
2470	ipft = ivm+iagec-1
2471	temp_begin = temp_begin + closure_intern(ipts,ipft,iele)
2472
2473	ENDDO ! nagec_pft
2474
2475	IF (temp_begin .GT. min_stomate) THEN
2476
2477	! mass balance
2478	err=err+1
2479	WRITE(numout,*) 'Error: mass balance is not closed in ', &
2480	TRIM(check_point), ' for ', var_name
2481	WRITE(numout,*) ' ipts, ivma; ', ipts,ivma
2482	WRITE(numout,*) ' first and last pft, ', &
2483	start_index(ivma), ipft
2484	WRITE(numout,*) ' Difference is, ', &
2485	temp_begin
2486	IF (err_act.GT.1) CALL ipslerr_p (plev, TRIM(check_point), &
2487	'Mass balance error','','')
2488
2489	ENDIF ! test for closure
2490
2491	ENDIF ! number of age classes
2492
2493	ENDDO ! ivmap
2494
2495
2496	!+++CHECK+++
2497	! A lot of output - only write a message if their is a problem
2498	!!$ ! One message per pixel
2499	!!$ IF (err_act.GT.1 .AND. err.EQ.zero) THEN
2500	!!$ WRITE(numout,*) 'Mass balance closure in ', &
2501	!!$ TRIM(check_point), ' for ',var_name
2502	!!$ ENDIF
2503	!++++++++++
2504
2505	ENDDO ! npts
2506
2507	ENDDO ! iele
2508
2509	CASE ('pixel')
2510
2511	! For testing subroutines where biomass can also move
2512	! from one species to another.
2513	DO iele = 1, nelements
2514
2515	IF (iele.EQ.icarbon) var_name = 'carbon'
2516	IF (iele.EQ.initrogen) var_name = 'nitrogen'
2517
2518	DO ipts=1,npts
2519
2520	err = zero
2521
2522	!check for NaN
2523	IF (isnan(SUM(pool_start(ipts,:,iele)))) THEN
2524	WRITE(numout,*) 'Error: found an NaN in ', &
2525	TRIM(check_point), ' for ', var_name
2526	WRITE(numout,*)'At the pixel',ipts,' for element',iele,&
2527	' pool_start is NaN'
2528	WRITE(numout,*) ' pool_start: ', pool_start(ipts,:,iele)
2529	CALL ipslerr_p(plev,TRIM(check_point),'NaN problem','','')
2530	ENDIF
2531	IF (isnan(SUM(pool_end(ipts,:,iele)))) THEN
2532	WRITE(numout,*) 'Error: found an NaN in ', &
2533	TRIM(check_point), ' for ', var_name
2534	WRITE(numout,*) 'At the pixel',ipts,' for element',iele,&
2535	' pool_end is NaN'
2536	WRITE(numout,*) ' pool_end: ', pool_end(ipts,:,iele)
2537	CALL ipslerr_p(plev,TRIM(check_point),'NaN problem','','')
2538	ENDIF
2539	IF (isnan(SUM(closure_intern(ipts,:,iele)))) THEN
2540	WRITE(numout,*) 'Error: found an NaN in ', &
2541	TRIM(check_point), ' for ', var_name
2542	WRITE(numout,*) 'At the pixel',ipts,' for element',iele,&
2543	' closure_intern is NaN'
2544	WRITE(numout,*) ' closure_intern: ', closure_intern(ipts,:,iele)
2545	CALL ipslerr_p(plev,TRIM(check_point), 'NaN problem','','')
2546	ENDIF
2547
2548	! Mass balance closure
2549	IF(ABS(SUM(closure_intern(ipts,:,iele))) .GT. min_stomate)THEN
2550
2551	err = err+1
2552
2553	IF(err_act.GT.1)THEN
2554	WRITE(numout,*) 'Error: mass balance is not closed in ', &
2555	TRIM(check_point), ' for ', var_name
2556	WRITE(numout,*) ' ipts, ', ipts
2557	WRITE(numout,*) ' Difference is, ', &
2558	SUM(closure_intern(ipts,:,iele))
2559	WRITE(numout,*) ' pool_end,pool_start: ', &
2560	SUM(pool_end(ipts,:,iele)), &
2561	SUM(pool_start(ipts,:,iele))
2562	! Exit or write warning depening on plev
2563	CALL ipslerr_p (plev, TRIM(check_point), &
2564	'Mass balance error','','')
2565	ENDIF
2566
2567	ENDIF
2568
2569	!+++CHECK+++
2570	! A lot of output - only write a message if their is a problem
2571	!!$ ! One message per pixel
2572	!!$ IF (err_act.GT.1 .AND. err.EQ.zero) THEN
2573	!!$ WRITE(numout,*) 'Mass balance closure in ', &
2574	!!$ TRIM(check_point), ' for ',var_name
2575	!!$ ENDIF
2576	!++++++++++
2577
2578	ENDDO ! npts
2579
2580	ENDDO ! iele
2581
2582	CASE ('products')
2583
2584	! This option is only used when err_act.GE.2
2585
2586	! testing mass balance for the product pool. Veget_max is not used
2587	! ivm is not defined for closure_intern
2588	DO iele=1,nelements
2589
2590	IF (iele.EQ.icarbon) var_name = 'carbon'
2591	IF (iele.EQ.initrogen) var_name = 'nitrogen'
2592
2593	DO ipts=1,npts
2594
2595	err = zero
2596
2597	IF(ABS(closure_intern(ipts,1,iele)) .GT. min_stomate)THEN
2598
2599	err = err+1
2600	WRITE(numout,*) 'Error: mass balance is not closed'//&
2601	'in basic_product_use', var_name
2602	WRITE(numout,*) 'ipts,ivm,iele ', ipts
2603	WRITE(numout,*) 'Difference is, ', &
2604	closure_intern(ipts,1,iele)
2605	WRITE(numout,*) 'pool_end, pool_start: ', &
2606	pool_end(ipts,1,iele), pool_start(ipts,1,iele)
2607	IF(err_act.GT.1) CALL ipslerr_p (plev,'basic_product_use', &
2608	'Mass balance error.','','')
2609	ENDIF
2610
2611	!+++CHECK+++
2612	! A lot of output - only write a message if their is a problem
2613	!!$ ! One message per pixel
2614	!!$ IF (err_act.GT.1 .AND. err.EQ.zero) THEN
2615	!!$ WRITE(numout,*) 'Mass balance closure in ', &
2616	!!$ TRIM(check_point), ' for ',var_name
2617	!!$ ENDIF
2618	!++++++++++
2619
2620	ENDDO
2621	ENDDO
2622
2623	CASE DEFAULT
2624
2625	WRITE(numout,*) 'Error: an unknown type/resolution was requested'
2626	WRITE(numout,*) ' for the mass balance check. Check the cases'
2627	WRITE(numout,*) ' defined in src_parameters/function_library or'
2628	WRITE(numout,*) ' the spelling of the type in the CALL argument',&
2629	test_type
2630	CALL ipslerr_p (3,'mass_balance_check in function_library',&
2631	'unknown type', 'Check *_out_orchidee for more details','')
2632
2633	END SELECT
2634
2635	END SUBROUTINE check_mass_balance
2636
2637
2638	!! ================================================================================================================================
2639	!! FUNCTION : intermediate_mass_balance_check
2640	!!
2641	!>\BRIEF
2642	!!
2643	!! DESCRIPTION : check whether the mass balance is closed in "the middle" of a stomate_growth_fun_all.
2644	!! Might be difficult to use the exact same code in other subroutines
2645	!!
2646	!! RECENT CHANGE(S): None
2647	!!
2648	!! RETURN VALUE : None
2649	!!
2650	!! REFERENCE(S) :
2651	!!
2652	!! FLOWCHART : None
2653	!! \n
2654	!_ ================================================================================================================================
2655
2656	SUBROUTINE intermediate_mass_balance_check(pool_start, pool_end, circ_class_biomass, &
2657	circ_class_n, veget_max, bm_alloc_tot, gpp_daily, atm_to_bm, dt, npts, &
2658	resp_maint, resp_growth, check_intern_init, ipts, j, label, test_type)
2659
2660	!! 0.1 Input variables
2661	INTEGER(i_std), INTENT(in) :: npts !! Domain size (unitless)
2662	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(in) :: circ_class_biomass !! Biomass components of the model tree
2663	!! within a circumference class
2664	!! class @tex $(g C ind^{-1})$ @endtex
2665	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: circ_class_n !! Number of individuals in each circ class
2666	!! @tex $(m^{-2})$ @endtex
2667	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: pool_start !! Start of this routine
2668	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
2669	REAL(r_std), DIMENSION(:,:), INTENT(in) :: veget_max !! PFT "Maximal" coverage fraction of a PFT
2670	!! @tex $(m^2 m^{-2})$ @endtex
2671	REAL(r_std), DIMENSION(:,:), INTENT(in) :: bm_alloc_tot !! Allocatable biomass for the whole plant
2672	!! @tex $(gC.m^{-2})$ @endtex
2673	REAL(r_std), DIMENSION(:,:), INTENT(in) :: gpp_daily !! PFT gross primary productivity
2674	!! @tex $(gC.m^{-2}dt^{-1})$ @endtex
2675	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: atm_to_bm !! Nitrogen and carbon which is added to the ecosystem to
2676	!! support vegetation growth (gC or gN/m2/day)
2677	REAL(r_std), INTENT(in) :: dt !! Time step of the simulations for stomate (days)
2678	REAL(r_std), DIMENSION(:,:), INTENT(in) :: resp_maint !! PFT maintenance respiration
2679	!! @tex $(gC.m^{-2}dt^{-1})$ @endtex
2680	REAL(r_std), DIMENSION(:,:), INTENT(in) :: resp_growth !! PFT growth respiration
2681	!! @tex $(gC.m^{-2}dt^{-1})$ @endtex
2682	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: check_intern_init !! Contains the components of the internal
2683	!! mass balance chech for this routine
2684	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
2685	INTEGER(i_std), INTENT(in) :: ipts, j !! indices
2686	CHARACTER(len=*), INTENT(in) :: label !! label to identify the mass balance check
2687
2688	!! 0.2 Output variables
2689
2690	!! 0.3 Modified variables
2691	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: pool_end !! Start and end pool of this routine
2692	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
2693	CHARACTER(*), INTENT(in) :: test_type !! Which test do we want to run on veget_max
2694
2695	!! 0.4 Local variables
2696	INTEGER(i_std) :: ipar, iele, icir !! Indices
2697	INTEGER(i_std) :: imbc !! Indices
2698	CHARACTER(len=80) :: message !! Label for error message
2699	REAL(r_std), DIMENSION(npts,nvm,nmbcomp,nelements):: check_intern !! Contains the components of the internal
2700	!! mass balance chech for this routine
2701	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
2702	REAL(r_std), DIMENSION(npts,nvm,nelements) :: closure_intern !! Check closure of internal mass balance
2703	!! @tex $(gC pixel^{-1} dt^{-1})$ @endtex
2704
2705	!_ ================================================================================================================================
2706
2707	SELECT CASE(test_type)
2708	CASE('pft')
2709	! Calculate pool_end
2710	DO ipar = 1,nparts
2711	DO iele = 1,nelements
2712	DO icir = 1,ncirc
2713	pool_end(:,:,iele) = pool_end(:,:,iele) + &
2714	(circ_class_biomass(:,:,icir,ipar,iele) * &
2715	circ_class_n(:,:,icir) * veget_max(:,:))
2716	ENDDO
2717	ENDDO
2718	ENDDO
2719
2720	! Calculate mass balance
2721	! Specific processes for carbon
2722	check_intern(:,:,:,:) = zero
2723	check_intern(:,:,iatm2land,icarbon) = &
2724	check_intern_init(:,:,iatm2land,icarbon) + &
2725	(gpp_daily(:,:) + atm_to_bm(:,:,icarbon)) * &
2726	dt * veget_max(:,:)
2727	check_intern(:,:,iatm2land,initrogen) = &
2728	check_intern_init(:,:,iatm2land,initrogen) + &
2729	atm_to_bm(:,:,initrogen) * dt * veget_max(:,:)
2730	check_intern(:,:,iland2atm,icarbon) = -un * &
2731	(resp_maint(:,:) + resp_growth(:,:)) * &
2732	veget_max(:,:)
2733
2734	! Common processes for icarbon and initrogen
2735	DO iele=1,nelements
2736	check_intern(:,:,ipoolchange,iele) = -un * (pool_end(:,:,iele) - &
2737	pool_start(:,:,iele))
2738	ENDDO
2739
2740	closure_intern = zero
2741	DO imbc = 1,nmbcomp
2742	DO iele=1,nelements
2743	closure_intern(:,:,iele) = closure_intern(:,:,iele) + &
2744	check_intern(:,:,imbc,iele)
2745	ENDDO
2746	ENDDO
2747
2748	CASE ('ipft')
2749
2750	! Calculate pool_end
2751	DO ipar = 1,nparts
2752	DO iele = 1,nelements
2753	DO icir = 1,ncirc
2754	pool_end(ipts,j,iele) = pool_end(ipts,j,iele) + &
2755	(circ_class_biomass(ipts,j,icir,ipar,iele) * &
2756	circ_class_n(ipts,j,icir) * veget_max(ipts,j))
2757	ENDDO
2758	ENDDO
2759	ENDDO
2760
2761	! Calculate mass balance
2762	! Specific processes for carbon
2763	check_intern(:,:,:,:) = zero
2764	check_intern(ipts,j,iatm2land,icarbon) = &
2765	check_intern_init(ipts,j,iatm2land,icarbon) + &
2766	(gpp_daily(ipts,j) + atm_to_bm(ipts,j,icarbon)) * &
2767	dt * veget_max(ipts,j)
2768	check_intern(ipts,j,iatm2land,initrogen) = &
2769	check_intern_init(ipts,j,iatm2land,initrogen) + &
2770	atm_to_bm(ipts,j,initrogen) * dt * veget_max(ipts,j)
2771	check_intern(ipts,j,iland2atm,icarbon) = -un * &
2772	(resp_maint(ipts,j) + resp_growth(ipts,j)) * &
2773	veget_max(ipts,j)
2774
2775	! Common processes for icarbon and initrogen
2776	DO iele=1,nelements
2777	check_intern(ipts,j,ipoolchange,iele) = -un * (pool_end(ipts,j,iele) - &
2778	pool_start(ipts,j,iele))
2779	ENDDO
2780
2781	closure_intern = zero
2782	DO imbc = 1,nmbcomp
2783	DO iele=1,nelements
2784	closure_intern(ipts,j,iele) = closure_intern(ipts,j,iele) + &
2785	check_intern(ipts,j,imbc,iele)
2786	ENDDO
2787	ENDDO
2788
2789	CASE DEFAULT
2790
2791	CALL ipslerr_p(3,'intermediate mass balance check', &
2792	'test_type not correcty defined','case not specified','')
2793
2794	END SELECT
2795
2796	message = 'Intermediate mass balance check '//TRIM(label)
2797	WRITE(numout,*) message, ' for pft, ', j
2798	CALL check_mass_balance(message, closure_intern, npts, pool_end, &
2799	pool_start, veget_max, test_type, j)
2800
2801	END SUBROUTINE intermediate_mass_balance_check
2802
2803
2804	!! ================================================================================================================================
2805	!! FUNCTION : check_area_change
2806	!!
2807	!>\BRIEF
2808	!!
2809	!! DESCRIPTION : check whether the change in vegetation is cancelled out by the
2810	!! change in non-biological land covers
2811	!!
2812	!! RECENT CHANGE(S): None
2813	!!
2814	!! RETURN VALUE : None
2815	!!
2816	!! REFERENCE(S) :
2817	!!
2818	!! FLOWCHART : None
2819	!! \n
2820	!_ ================================================================================================================================
2821
2822	SUBROUTINE check_area_change(check_point, npts, frac_nobio, frac_nobio_new, loss_gain)
2823
2824	!! 0.1 Input variables
2825	INTEGER, INTENT(in) :: npts !! Domain size - number of pixels (unitless)
2826	REAL(r_std), DIMENSION(:,:),INTENT(in) :: frac_nobio !! Fraction of grid cell covered by lakes, land
2827	!! ice, cities, ... before land cover change (unitless, 0-1)
2828	REAL(r_std),DIMENSION(:,:),INTENT(in) :: frac_nobio_new !! Fraction of grid cell covered by lakes, land
2829	!! ice, cities, ... after land cover change (unitless, 0-1)
2830	REAL(r_std), DIMENSION(:,:), INTENT(in) :: loss_gain !! losses and gains due to LCC distributed over all
2831	!! age classes and thus taking the age-classes into
2832	!! account (unitless, 0-1)
2833	CHARACTER(*),INTENT(in) :: check_point !! A flag to indicate at which
2834	!! point in the code we're doing
2835	!! this check
2836
2837	!! 0.2 Output variables
2838
2839	!! 0.3 Modified variables
2840
2841	!! 0.4 Local variables
2842	INTEGER(i_std) :: ipts !! Index
2843	INTEGER(i_std), DIMENSION(npts) :: count !! counter
2844	REAL(r_std) :: total !! temporary variable for total area of a pixel (unitless, 0-1)
2845	REAL(r_std), DIMENSION(npts) :: change_nobio !! Change in the non biological fraction in a pixel (0-1, unitless)
2846
2847	!_ ================================================================================================================================
2848
2849	! Quality check. It is now expected that the changes cancel out each other
2850	! and are exactley zero.
2851	count(:) = zero
2852	change_nobio(:) = SUM(frac_nobio_new(:,:) - frac_nobio(:,:),2)
2853	WHERE(ABS(SUM(loss_gain(:,:),2)+change_nobio(:)).GT.10*EPSILON(un))
2854	count(:) = un
2855	END WHERE
2856
2857	! Error checking
2858	IF(SUM(count(:)).GT.zero)THEN
2859
2860	! Refine the error message
2861	DO ipts = 1,npts
2862	IF (ABS(SUM(loss_gain(ipts,:))+change_nobio(ipts)).GT.10*EPSILON(un)) THEN
2863	total = ABS(SUM(loss_gain(ipts,:))+change_nobio(ipts))
2864	WRITE(numout,*) 'ipts, change_nobio, loss_gain, mismatch, ', &
2865	ipts, change_nobio(ipts), SUM(loss_gain(ipts,:)), total
2866	WRITE(numout,*) 'ipts, frac_nobio_new, frac_nobio. ', &
2867	ipts, SUM(frac_nobio_new(ipts,:)), SUM(frac_nobio(ipts,:))
2868	END IF
2869	END DO
2870	WRITE(numout,*) 'Found ,',SUM(count(:)), &
2871	' pixels for which the total change differs from zero (=0)'
2872	CALL ipslerr_p(3,TRIM(check_point), &
2873	'The change in surface areas do not cancel out to zero', '','')
2874
2875	END IF
2876
2877	END SUBROUTINE check_area_change
2878
2879
2880	!! ================================================================================================================================
2881	!! FUNCTION : check_pixel_area
2882	!!
2883	!>\BRIEF
2884	!!
2885	!! DESCRIPTION : check whether the vegetation and non-biological fractions
2886	!! still add up to one (=1).
2887	!!
2888	!! RECENT CHANGE(S): None
2889	!!
2890	!! RETURN VALUE : None
2891	!!
2892	!! REFERENCE(S) :
2893	!!
2894	!! FLOWCHART : None
2895	!! \n
2896	!_ ================================================================================================================================
2897
2898	SUBROUTINE check_pixel_area(check_point, npts, veget_max, frac_nobio)
2899
2900	!! 0.1 Input variables
2901	INTEGER, INTENT(in) :: npts !! Domain size - number of pixels (unitless)
2902	REAL(r_std), DIMENSION(:,:),INTENT(in) :: veget_max !! Cover fraction of the vegetation (unitless, 0-1)
2903	REAL(r_std),DIMENSION(:,:),INTENT(in) :: frac_nobio !! Fraction of grid cell covered by lakes, land
2904	!! ice, cities, ... (unitless, 0-1)
2905	CHARACTER(*),INTENT(in) :: check_point !! A flag to indicate at which
2906	!! point in the code we're doing
2907	!! this check
2908
2909	!! 0.2 Output variables
2910
2911	!! 0.3 Modified variables
2912
2913	!! 0.4 Local variables
2914	INTEGER(i_std) :: ipts !! Index
2915	INTEGER(i_std), DIMENSION(npts) :: count !! counter
2916
2917	!_ ================================================================================================================================
2918
2919	! Quality check. It is still expected that the different vegetation
2920	! fractions in each pixel sums up to exactly one.
2921	count(:) = zero
2922	WHERE(ABS(un-SUM(veget_max(:,:),2)-SUM(frac_nobio(:,:),2)).GT.10*EPSILON(un))
2923	count(:) = un
2924	END WHERE
2925
2926	! Error checking
2927	IF(SUM(count(:)).GT.zero)THEN
2928
2929	! Refine the error message
2930	DO ipts = 1,npts
2931	IF (ABS(un-SUM(veget_max(ipts,:))-SUM(frac_nobio(ipts,:))).GT.10*EPSILON(un)) THEN
2932	WRITE(numout,*) 'ipts, frac_nobio, veget_max_new, mismatch, ', ipts, &
2933	SUM(frac_nobio(ipts,:)), SUM(veget_max(ipts,:)), &
2934	un-SUM(veget_max(ipts,:))-SUM(frac_nobio(ipts,:))
2935	END IF
2936	END DO
2937	WRITE(numout,*) 'Found ,',SUM(count(:)), ' pixels for which the total surface area &
2938	is not one'
2939	CALL ipslerr_p(3,TRIM(check_point), &
2940	'The sum of the vegetation and non-biological fractions',&
2941	'differs from one','')
2942
2943	END IF
2944
2945	END SUBROUTINE check_pixel_area
2946
2947
2948	!! ================================================================================================================================
2949	!! FUNCTION : check_variable_change
2950	!!
2951	!>\BRIEF
2952	!!
2953	!! DESCRIPTION : check whether the relationship between veget_max, veget_max_new and
2954	!! loss_gain holds
2955	!!
2956	!! RECENT CHANGE(S): None
2957	!!
2958	!! RETURN VALUE : None
2959	!!
2960	!! REFERENCE(S) :
2961	!!
2962	!! FLOWCHART : None
2963	!! \n
2964	!_ ================================================================================================================================
2965
2966	SUBROUTINE check_variable_change(check_point, npts, veget_max, veget_max_new, loss_gain)
2967
2968	!! 0.1 Input variables
2969	INTEGER, INTENT(in) :: npts !! Domain size - number of pixels (unitless)
2970	REAL(r_std), DIMENSION(:,:),INTENT(in) :: veget_max !! Cover fraction of the vegetation before land cover change (unitless, 0-1)
2971	REAL(r_std),DIMENSION(:,:),INTENT(in) :: veget_max_new !! Cover fraction of the vegetation after land cover change (unitless, 0-1)
2972	REAL(r_std), DIMENSION(:,:), INTENT(in) :: loss_gain !! losses and gains due to LCC distributed over all
2973	!! age classes and thus taking the age-classes into
2974	!! account (unitless, 0-1)
2975	CHARACTER(*),INTENT(in) :: check_point !! A flag to indicate at which
2976	!! point in the code we're doing
2977	!! this check
2978
2979	!! 0.2 Output variables
2980
2981	!! 0.3 Modified variables
2982
2983	!! 0.4 Local variables
2984	INTEGER(i_std) :: ipts !! Index
2985	INTEGER(i_std), DIMENSION(npts) :: count !! counter
2986	REAL(r_std) :: total !! temporary variable for total area of a pixel (unitless, 0-1)
2987	!_ ================================================================================================================================
2988
2989	! Quality check. In slowproc veget_max_new was calculated as veget_max
2990	! plus loss_gain. If all went well veget_max, loss_gain and
2991	! veget_max_new are not changed in between slowproc and
2992	! sapiens_lcchange. One of the reason this assumption is valid is
2993	! because LCC is calculated at the first day of the year whereas
2994	! the other processes that can change veget_max are calculated
2995	! the last day of the year. If the above is not true we have to
2996	! carefully check what is going on withe veget_max, loss_gain
2997	! and/ot veget_max_new.
2998	count(:) = zero
2999	WHERE(ABS(SUM(veget_max_new(:,:)-loss_gain(:,:)-veget_max(:,:),2)).GT.10*EPSILON(un))
3000	count(:) = un
3001	END WHERE
3002
3003	! Error checking
3004	IF(SUM(count(:)).GT.zero)THEN
3005
3006	! Refine the error message
3007	DO ipts = 1,npts
3008	IF (ABS(SUM(veget_max_new(ipts,:)-loss_gain(ipts,:) - &
3009	veget_max(ipts,:))).GT.10*EPSILON(un)) THEN
3010	total = ABS(SUM(veget_max_new(ipts,:)-loss_gain(ipts,:) - &
3011	veget_max(ipts,:)))
3012	WRITE(numout,*) 'ipts, veget_max_new, veget_max, loss_gain, mismatch, ', &
3013	ipts, SUM(veget_max_new(ipts,:)), SUM(veget_max(ipts,:)), &
3014	SUM(loss_gain(ipts,:)), total
3015	END IF
3016	END DO
3017	WRITE(numout,*) 'Found ,',SUM(count(:)), ' pixels for which the total surface area &
3018	is not one'
3019	CALL ipslerr_p(3,TRIM(check_point),&
3020	'Too big changes in surface area when updating veget_max_new','','')
3021
3022	END IF
3023
3024	END SUBROUTINE check_variable_change
3025
3026
3027	!! ================================================================================================================================
3028	!! FUNCTION : check_variable_swap
3029	!!
3030	!>\BRIEF
3031	!!
3032	!! DESCRIPTION : at one point in sapiens_lcchange veget_max_new should be
3033	!! equal to veget_max. Check whether this is the case
3034	!!
3035	!! RECENT CHANGE(S): None
3036	!!
3037	!! RETURN VALUE : None
3038	!!
3039	!! REFERENCE(S) :
3040	!!
3041	!! FLOWCHART : None
3042	!! \n
3043	!_ ================================================================================================================================
3044
3045	SUBROUTINE check_variable_swap(check_point, npts, veget_max, veget_max_new)
3046
3047	!! 0.1 Input variables
3048	INTEGER, INTENT(in) :: npts !! Domain size - number of pixels (unitless)
3049	REAL(r_std), DIMENSION(:,:),INTENT(in) :: veget_max !! Cover fraction of the vegetation before land cover change (unitless, 0-1)
3050	REAL(r_std),DIMENSION(:,:),INTENT(in) :: veget_max_new !! Cover fraction of the vegetation after land cover change (unitless, 0-1)
3051	CHARACTER(*),INTENT(in) :: check_point !! A flag to indicate at which
3052	!! point in the code we're doing
3053	!! this check
3054
3055	!! 0.2 Output variables
3056
3057	!! 0.3 Modified variables
3058
3059	!! 0.4 Local variables
3060	INTEGER(i_std) :: ipts !! Index
3061	INTEGER(i_std), DIMENSION(npts) :: count !! counter
3062	REAL(r_std) :: total !! temporary variable for total area of a pixel (unitless, 0-1)
3063	!_ ================================================================================================================================
3064
3065	count(:) = zero
3066	WHERE(ABS(SUM(veget_max(:,:) - veget_max_new(:,:),2)).GT.10*EPSILON(un))
3067	count(:) = un
3068	END WHERE
3069
3070	! Error checking
3071	IF(SUM(count(:)).GT.zero)THEN
3072
3073	! Refine the error message
3074	DO ipts = 1,npts
3075	IF (ABS(SUM(veget_max(ipts,:) - veget_max_new(ipts,:))).GT.10*EPSILON(un)) THEN
3076	total = ABS(SUM(veget_max(ipts,:) - veget_max_new(ipts,:)))
3077	WRITE(numout,*) 'ipts, veget_max, veget_max_new, mismatch, ', &
3078	ipts, SUM(veget_max(ipts,:)), SUM(veget_max_new(ipts,:)), total
3079	END IF
3080	END DO
3081	WRITE(numout,*) 'Found ,',SUM(count(:)), &
3082	' pixels for which the updated veget_max differs from veget_max_new'
3083	CALL ipslerr_p(3,TRIM(check_point),&
3084	'Updated veget_max but the update differs from veget_max_new', &
3085	'improve slowproc_adjust_delta_vegetmax','')
3086
3087	END IF
3088
3089	END SUBROUTINE check_variable_swap
3090
3091
3092	! ================================================================================================================================
3093	!! SUBROUTINE : check_biomass_change
3094	!!
3095	!>\BRIEF Check mass conservation for the biomass pool in sapiens_lcchange
3096	!!
3097	!! DESCRIPTION : Within sapiens_lcchange, biomass of pfts may change but in a
3098	!! predictive way depending on whether there was a loss or gain in veget_max.
3099	!!
3100	!! RECENT CHANGE(S) :
3101	!!
3102	!! MAIN OUTPUT VARIABLE(S) : None
3103	!!
3104	!! REFERENCES : None
3105	!!
3106	!! FLOWCHART :
3107	!! \n
3108	!_ ================================================================================================================================
3109
3110	SUBROUTINE check_biomass_change(npts, dt_days, area, circ_class_biomass, &
3111	circ_class_n, veget_max, loss_gain, atm_to_bm_old, atm_to_bm, &
3112	harvest_pool, harvest_pool_old, fresh_litter, bm_to_litter_old, &
3113	biomass_start, biomass_end, n_call)
3114
3115	!! 0. Variable declaration
3116
3117	!! 0.1 Input variables
3118	INTEGER, INTENT(in) :: npts !! Domain size - number of pixels (unitless)
3119	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(in) :: circ_class_biomass !! Biomass components of the model tree within a
3120	!! circumference class @tex $(g C ind^{-1})$ @endtex
3121	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: circ_class_n !! Number of individuals in each circ class
3122	!! @tex $(ind m^{-2})$ @endtex
3123	REAL(r_std), DIMENSION(:,:), INTENT(in) :: veget_max !! "maximal" coverage fraction of a PFT on the ground
3124	!! May sum to
3125	!! less than unity if the pixel has
3126	!! nobio area. (unitless, 0-1)
3127	REAL(r_std), DIMENSION(:,:), INTENT(in) :: loss_gain !! losses and gains due to LCC distributed over all
3128	!! age classes and thus taking the age-classes into
3129	!! account (unitless, 0-1)
3130	INTEGER(i_std), INTENT(in) :: n_call !! Number indicating whether it is the first or second call
3131	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: atm_to_bm_old !! Temporary variable to store atm_to_bm at the start of
3132	!! the routine @tex (gC.m^{-2}dt^{-1})$ @endtex
3133	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: atm_to_bm !! C and N taken from the atmosphere to get C and N to
3134	!! create the seedlings @tex (gC.m^{-2}dt^{-1})$ @endtex
3135	REAL(r_std), INTENT(in) :: dt_days !! Time step of vegetation dynamics for stomate (days)
3136	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: harvest_pool !! The wood and biomass that have been
3137	!! havested by humans @tex $(gC)$ @endtex
3138	REAL(r_std),DIMENSION(:,:,:,:), INTENT(in) :: fresh_litter !! Pool of fresh litter that is being released during @tex $(gC)$ @endtex
3139	REAL(r_std), DIMENSION(:), INTENT(in) :: area !! surface area of the pixel @tex ($m^{2}$) @endtex
3140	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: harvest_pool_old !! Temporary variable to check mass conservation
3141	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: bm_to_litter_old !! Temporary transfer of biomass to litter
3142	!! @tex $(gC m^{-2} dtslow^{-1})$ @endtex
3143
3144
3145	!! 0.2 Output variables
3146
3147	!! 0.3 Modified variables
3148	REAL(r_std), DIMENSION(:,:,:),INTENT(inout) :: biomass_start !! biomass before LCC in gC or gN.
3149	REAL(r_std), DIMENSION(:,:,:),INTENT(inout) :: biomass_end !! biomass after LCC in gC or gN.
3150
3151	!! 0.4 Local variables
3152	INTEGER(i_std) :: iele,icir,ivm,ipts !! Indices
3153	REAL(r_std), DIMENSION(nelements) :: check !! Som of pools and fluxes gC or gN pixel-1.
3154	!_ ================================================================================================================================
3155
3156	IF (n_call.EQ.1) THEN
3157	! Calculate the total biomass at the start of the routine
3158	DO ipts = 1,npts
3159	DO ivm = 1,nvm
3160	DO iele = 1, nelements
3161	DO icir = 1,ncirc
3162	! Multiply with area to make the units easier to understand.
3163	! The unit of biomass_start is gC pixel-1
3164	biomass_start(ipts,ivm,iele) = biomass_start(ipts,ivm,iele) + &
3165	SUM(circ_class_biomass(ipts,ivm,icir,:,iele)) * &
3166	circ_class_n(ipts,ivm,icir) * veget_max(ipts,ivm) * area(ipts)
3167	END DO
3168	END DO
3169	END DO
3170	END DO
3171	END IF ! n_call.eq.1
3172
3173	IF (n_call.EQ.2) THEN
3174	! First calculate the total biomass at the end of the routine
3175	DO ipts = 1,npts
3176	DO ivm = 1,nvm
3177	DO iele = 1, nelements
3178	DO icir = 1,ncirc
3179	! Multiply with area to make the units easier to understand.
3180	! The unit of biomass_start is gC pixel-1
3181	biomass_end(ipts,ivm,iele) = biomass_end(ipts,ivm,iele) + &
3182	SUM(circ_class_biomass(ipts,ivm,icir,:,iele)) * &
3183	circ_class_n(ipts,ivm,icir) * veget_max(ipts,ivm) *area(ipts)
3184	END DO
3185	END DO
3186
3187	! Check for mass conservation
3188	check(:) = zero
3189	IF (loss_gain(ipts,ivm).GE.zero) THEN
3190
3191	! Gain in surface area. Biomass_end includes the newly established
3192	! vegetation. The C and N in this vegetation was taken from the air
3193	! and stored in atm_to_bm. Subtracting atm_to_bm from biomass_end should
3194	! give biomass_start. biopmass_start and biomass_end are in
3195	! gG or gN pixel-1. Atm_to_bm is in gC m-2 dt-1.
3196	check(:) = (biomass_start(ipts,ivm,:) - biomass_end(ipts,ivm,:) + &
3197	(atm_to_bm(ipts,ivm,:) - atm_to_bm_old(ipts,ivm,:)) * &
3198	veget_max(ipts,ivm) * area(ipts) * dt_days)/area(ipts)
3199
3200	ELSE
3201
3202	! Loss in surface area. When surface area is lost, the biomass was
3203	! moved into the harvest_pool and fresh_litter. If all these pools are
3204	! accounted for in biomass_end, it should be possible to reconstruct
3205	! biomass_start. biomass_end, biomass_start ar in gC or gN pixel-1.
3206	! fresh_litter is in gC or gN m-2. harvest_pool is in gC or gN pixel-1.
3207	! Note that in case of a species change there will already be some
3208	! biomass in the harvest_pool. In case of a classic land cover change
3209	! this is extremely unlikley for the moment (but this could be the
3210	! case if salavage logging following wind, fires or pests happens in
3211	! the same time step). Fresh_litter contains part of biomass that was
3212	! removed as well as bm_to_litter. To check for mass balance conservation
3213	! of the biomass we need to account for the biomass that went into the
3214	! fresh_litter but not for the bm_to_litter. The latter comes from a
3215	! previous mortality event.
3216	check(:) = (biomass_start(ipts,ivm,:) - biomass_end(ipts,ivm,:) - &
3217	SUM(fresh_litter(ipts,ivm,:,:)-bm_to_litter_old(ipts,ivm,:,:),1) * &
3218	ABS(loss_gain(ipts,ivm)) * area(ipts) - &
3219	SUM(harvest_pool(ipts,ivm,:,:) - harvest_pool_old(ipts,ivm,:,:),1)) / &
3220	area(ipts)
3221
3222	END IF
3223
3224	! Check for errors
3225	DO iele = 1,nelements
3226	IF (ABS(check(iele)).GT.min_stomate) THEN
3227	WRITE(numout,*) 'ipts, ivm, iele, ', ipts, ivm, iele
3228	WRITE(numout,*) 'loss_gain, check, ', loss_gain(ipts,ivm), check(iele)
3229	WRITE(numout,*) 'biomass_start, ', biomass_start(ipts,ivm,iele)
3230	WRITE(numout,*) 'biomass_end, ', biomass_end(ipts,ivm,iele)
3231	WRITE(numout,*) 'diff, ', biomass_end(ipts,ivm,iele)-biomass_start(ipts,ivm,iele)
3232	WRITE(numout,*) 'harvest_pool, ', SUM(harvest_pool(ipts,ivm,:,iele))
3233	WRITE(numout,*) 'fresh_litter from biomass, ', &
3234	SUM(fresh_litter(ipts,ivm,:,iele)-bm_to_litter_old(ipts,ivm,:,iele),1) * &
3235	ABS(loss_gain(ipts,ivm)) * area(ipts)
3236	WRITE(numout,) 'atm_to_bm, ', (atm_to_bm(ipts,ivm,iele) - atm_to_bm_old(ipts,ivm,iele)) &
3237	veget_max(ipts,ivm) * area(ipts) * dt_days
3238	CALL ipslerr_p(3,'sapiens_lcchange','Mass balance error for circ_class_biomass,',&
3239	'harvest_pool, fresh_litter and/or atm_to_bm','')
3240	ENDIF
3241	END DO
3242
3243	END DO
3244
3245	END DO
3246
3247	ENDIF ! n_call.eq.2
3248
3249	END SUBROUTINE check_biomass_change
3250
3251
3252	!! ================================================================================================================================
3253	!! FUNCTION : sort_circ_class_biomass
3254	!!
3255	!>\BRIEF
3256	!!
3257	!! DESCRIPTION : sort the diameters in circ_class_biomass from low to high and
3258	!! change the order in circ_class_n accordingly
3259	!!
3260	!! RECENT CHANGE(S): None
3261	!!
3262	!! RETURN VALUE : ::circ_class_biomass (gC tree-1),circ_class_n (tree m-2)
3263	!!
3264	!! REFERENCE(S) :
3265	!!
3266	!! FLOWCHART : None
3267	!! \n
3268	!_ ================================================================================================================================
3269
3270	SUBROUTINE sort_circ_class_biomass(circ_class_biomass_new,circ_class_n_new)
3271
3272	!! 0. Variable and parameter declaration
3273
3274	!! 0.1 Input variables
3275
3276	!! 0.2 Output variables
3277
3278	!! 0.3 Modified variables
3279	REAL(r_std), DIMENSION(:), INTENT(inout) :: circ_class_n_new
3280	REAL(r_std), DIMENSION(:,:,:), INTENT(inout) :: circ_class_biomass_new
3281
3282	!! 0.4 Local variables
3283	INTEGER(i_std) :: inc,icir,jcir !! Indices
3284	REAL(r_std), DIMENSION(nparts,nelements) :: v1 !! temporay variable for sorting circ_class_biomass
3285	REAL(r_std) :: v2 !! temporay variables for sorting circ_class_n
3286	REAL(r_std) :: sort !! flag triggering sorting routine (0 or 1)
3287	REAL(r_std), DIMENSION(ncirc) :: tree_size !! temporay variables for checking the biomass
3288	!! for each circ class @tex ($gC tree^{-1}$) @endtex
3289
3290	!_ ================================================================================================================================
3291
3292	! Calculate the total biomass in each circumference class
3293	tree_size(:)=zero
3294	DO icir=1,ncirc
3295
3296	! Only check the order of the carbon pools
3297	tree_size(icir)=SUM(circ_class_biomass_new(icir,:,icarbon))
3298
3299	ENDDO
3300
3301	! Check whether the order was preserved
3302	sort = zero
3303	DO icir=2,ncirc
3304
3305	IF(tree_size(icir) .LT. tree_size(icir-1))THEN
3306
3307	! The order is not preserved
3308	sort = un
3309	EXIT
3310
3311	ENDIF
3312
3313	ENDDO
3314
3315	! Sort the biomasses. Note that at the same time we also
3316	! have to sort circ_class_n. We made use of the Shell sort
3317	! method
3318	IF (sort==1) THEN
3319
3320	! Based on nummerical recipies in Fortran
3321	inc = 1
3322
3323	DO
3324	inc=3*inc+1
3325	IF (inc .GT. ncirc) EXIT
3326
3327	ENDDO
3328
3329	DO
3330	inc = inc/3
3331
3332	DO icir=inc+1,ncirc
3333
3334	v1(:,:) = circ_class_biomass_new(icir,:,:)
3335	v2 = circ_class_n_new(icir)
3336	jcir=icir
3337
3338	DO
3339
3340	! Use the carbon pool to sort the circumference classes
3341	! but when moving circ classes around, also move the
3342	! other elements, i.e. nitrogen.
3343	IF (SUM(circ_class_biomass_new(jcir-inc,:,icarbon)) .LE. &
3344	SUM(v1(:,icarbon))) EXIT
3345	circ_class_biomass_new(jcir,:,:) = &
3346	circ_class_biomass_new(jcir-inc,:,:)
3347	circ_class_n_new(jcir) = circ_class_n_new(jcir-inc)
3348	jcir = jcir-inc
3349	IF (jcir .LE. inc) EXIT
3350
3351	ENDDO
3352
3353	circ_class_biomass_new(jcir,:,:)=v1(:,:)
3354	circ_class_n_new(jcir)=v2
3355
3356	ENDDO
3357
3358	IF (inc .LE. 1) EXIT
3359
3360	ENDDO
3361
3362	!!$ ! Debug
3363	!!$ WRITE(numout,*) 'sort end - carbon', &
3364	!!$ SUM(circ_class_biomass_new(:,:,icarbon),2),&
3365	!!$ circ_class_n_new(:)
3366	!!$ WRITE(numout,*) 'sort end - nitrogen', &
3367	!!$ SUM(circ_class_biomass_new(:,:,initrogen),2)
3368	!!$ !-
3369
3370	ENDIF
3371
3372	END SUBROUTINE sort_circ_class_biomass
3373
3374	!! ================================================================================================================================
3375	!! SUBROUTINE : Arrhenius_1d
3376	!!
3377	!>\BRIEF Calculate temperature dependency based on Arrhenius
3378	!!
3379	!! DESCRIPTION :
3380	!!
3381	!! RECENT CHANGE(S): None
3382	!!
3383	!! MAIN OUTPUT VARIABLE(S): :: val_arrhenius
3384	!!
3385	!! REFERENCE(S) :
3386	!!
3387	!! FLOWCHART : None
3388	!_ ================================================================================================================================
3389
3390	FUNCTION Arrhenius_1d (kjpindex,temp,ref_temp,energy_act) RESULT ( val_arrhenius )
3391	!! 0.1 Input variables
3392
3393	INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size (-)
3394	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: temp !! Temperature (K)
3395	REAL(r_std), INTENT(in) :: ref_temp !! Temperature of reference (K)
3396	REAL(r_std),INTENT(in) :: energy_act !! Activation Energy (J mol-1)
3397
3398	!! 0.2 Result
3399
3400	REAL(r_std), DIMENSION(kjpindex) :: val_arrhenius !! Temperature dependance based on
3401	!! a Arrhenius function (-)
3402	!_ ================================================================================================================================
3403
3404	val_arrhenius(:)=EXP(((temp(:)-ref_temp)energy_act)/(ref_tempRR*(temp(:))))
3405	END FUNCTION Arrhenius_1d
3406
3407	!!
3408	!================================================================================================================================
3409	!! SUBROUTINE : Arrhenius_0d
3410	!!
3411	!>\BRIEF Calculate temperature dependency based on Arrhenius for a
3412	!! single pixel
3413	!!
3414	!! DESCRIPTION :
3415	!!
3416	!! RECENT CHANGE(S): None
3417	!!
3418	!! MAIN OUTPUT VARIABLE(S): :: val_arrhenius
3419	!!
3420	!! REFERENCE(S) :
3421	!!
3422	!! FLOWCHART : None
3423	!_
3424	!================================================================================================================================
3425
3426	FUNCTION Arrhenius_0d (temp,ref_temp,energy_act) RESULT ( val_arrhenius )
3427	!! 0.1 Input variables
3428
3429	REAL(r_std),INTENT(in) :: temp !! Temperature (K)
3430	REAL(r_std), INTENT(in) :: ref_temp !! Temperature of reference (K)
3431	REAL(r_std),INTENT(in) :: energy_act !! Activation Energy (J mol-1)
3432
3433	!! 0.2 Result
3434
3435	REAL(r_std) :: val_arrhenius !! Temperature dependance based on
3436	!! a Arrhenius function (-)
3437	!_
3438	!================================================================================================================================
3439
3440	val_arrhenius=EXP(((temp-ref_temp)energy_act)/(ref_tempRR*(temp)))
3441	END FUNCTION Arrhenius_0d
3442
3443
3444	!! ================================================================================================================================
3445	!! SUBROUTINE : Arrhenius_modified_1d
3446	!!
3447	!>\BRIEF Calculate temperature dependency based on Arrhenius
3448	!!
3449	!! DESCRIPTION :
3450	!!
3451	!! RECENT CHANGE(S): None
3452	!!
3453	!! MAIN OUTPUT VARIABLE(S): :: val_arrhenius
3454	!!
3455	!! REFERENCE(S) :
3456	!!
3457	!! FLOWCHART : None
3458	!_ ================================================================================================================================
3459
3460	FUNCTION Arrhenius_modified_1d (kjpindex,temp,ref_temp,energy_act,energy_deact,entropy) RESULT ( val_arrhenius )
3461	!! 0.1 Input variables
3462
3463	INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size (-)
3464	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: temp !! Temperature (K)
3465	REAL(r_std), INTENT(in) :: ref_temp !! Temperature of reference (K)
3466	REAL(r_std),INTENT(in) :: energy_act !! Activation Energy (J mol-1)
3467	REAL(r_std),INTENT(in) :: energy_deact !! Deactivation Energy (J mol-1)
3468	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: entropy !! Entropy term (J K-1 mol-1)
3469
3470	!! 0.2 Result
3471
3472	REAL(r_std), DIMENSION(kjpindex) :: val_arrhenius !! Temperature dependance based on
3473	!! a Arrhenius function (-)
3474	!_ ================================================================================================================================
3475
3476	val_arrhenius(:)=EXP(((temp(:)-ref_temp)energy_act)/(ref_tempRR*(temp(:)))) &
3477	* (1. + EXP( (ref_temp * entropy(:) - energy_deact) / (ref_temp * RR ))) &
3478	/ (1. + EXP( (temp(:) * entropy(:) - energy_deact) / ( RR*temp(:))))
3479
3480	END FUNCTION Arrhenius_modified_1d
3481
3482	!! ================================================================================================================================
3483	!! SUBROUTINE : Arrhenius_modified_0d
3484	!!
3485	!>\BRIEF Calculate temperature dependency based on Arrhenius
3486	!!
3487	!! DESCRIPTION :
3488	!!
3489	!! RECENT CHANGE(S): None
3490	!!
3491	!! MAIN OUTPUT VARIABLE(S): :: val_arrhenius
3492	!!
3493	!! REFERENCE(S) :
3494	!!
3495	!! FLOWCHART : None
3496	!_ ================================================================================================================================
3497
3498	FUNCTION Arrhenius_modified_0d (kjpindex,temp,ref_temp,energy_act,energy_deact,entropy) RESULT ( val_arrhenius )
3499	!! 0.1 Input variables
3500
3501	INTEGER(i_std),INTENT(in) :: kjpindex !! Domain size (-)
3502	REAL(r_std),DIMENSION(kjpindex),INTENT(in) :: temp !! Temperature (K)
3503	REAL(r_std), INTENT(in) :: ref_temp !! Temperature of reference (K)
3504	REAL(r_std),INTENT(in) :: energy_act !! Activation Energy (J mol-1)
3505	REAL(r_std),INTENT(in) :: energy_deact !! Deactivation Energy (J mol-1)
3506	REAL(r_std),INTENT(in) :: entropy !! Entropy term (J K-1 mol-1)
3507
3508	!! 0.2 Result
3509
3510	REAL(r_std), DIMENSION(kjpindex) :: val_arrhenius !! Temperature dependance based on
3511	!! a Arrhenius function (-)
3512	!_ ================================================================================================================================
3513
3514	val_arrhenius(:)=EXP(((temp(:)-ref_temp)energy_act)/(ref_tempRR*(temp(:)))) &
3515	* (1. + EXP( (ref_temp * entropy - energy_deact) / (ref_temp * RR ))) &
3516	/ (1. + EXP( (temp(:) * entropy - energy_deact) / ( RR*temp(:))))
3517
3518	END FUNCTION Arrhenius_modified_0d
3519
3520	!!
3521	!================================================================================================================================
3522	!! SUBROUTINE : Arrhenius_modified_nodim
3523	!!
3524	!>\BRIEF Calculate temperature dependency based on Arrhenius.
3525	!! This function is for a single pixel output.
3526	!!
3527	!! DESCRIPTION :
3528	!!
3529	!! RECENT CHANGE(S): None
3530	!!
3531	!! MAIN OUTPUT VARIABLE(S): :: val_arrhenius
3532	!!
3533	!! REFERENCE(S) :
3534	!!
3535	!! FLOWCHART : None
3536	!_
3537	!================================================================================================================================
3538
3539	FUNCTION Arrhenius_modified_nodim (temp,ref_temp,energy_act,energy_deact,entropy) RESULT (val_arrhenius)
3540	!! 0.1 Input variables
3541
3542	REAL(r_std),INTENT(in) :: temp !! Temperature (K)
3543	REAL(r_std),INTENT(in) :: ref_temp !! Temperature of reference (K)
3544	REAL(r_std),INTENT(in) :: energy_act !! Activation Energy (J mol-1)
3545	REAL(r_std),INTENT(in) :: energy_deact !! Deactivation Energy (J mol-1)
3546	REAL(r_std),INTENT(in) :: entropy !! Entropy term (J K-1 mol-1)
3547
3548	!! 0.2 Result
3549
3550	REAL(r_std) :: val_arrhenius !! Temperature dependance based on
3551	!! a Arrhenius function (-)
3552	!_
3553	!================================================================================================================================
3554
3555	val_arrhenius=EXP(((temp-ref_temp)energy_act)/(ref_tempRR*(temp))) &
3556	* (1. + EXP( (ref_temp * entropy - energy_deact) / (ref_temp * RR ))) &
3557	/ (1. + EXP( (temp * entropy - energy_deact) / ( RR*temp)))
3558
3559	END FUNCTION Arrhenius_modified_nodim
3560
3561	END MODULE function_library

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

Download in other formats: