Context Navigation

source: branches/publications/ORCHIDEE_gmd-2018-261/src_parameters/function_library.f90

Last change on this file was 3821, checked in by nicolas.vuichard, 8 years ago
add dynamic SLA
File size: 33.8 KB

Line
1	! =================================================================================================================================
2	! MODULE : function_library
3	!
4	! CONTACT : orchidee-help _at_ ipsl.jussieu.fr
5	!
6	! LICENCE : IPSL (2006)
7	! This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
8	!
9	!>\BRIEF Collection of functions that are used throughout the ORCHIDEE code
10	!!
11	!!\n DESCRIPTION: Collection of modules to : (1) convert one variable into another i.e. basal area
12	!! to diameter, diamter to tree height, diameter to crown area, etc. (2) ...
13	!!
14	!! RECENT CHANGE(S): None
15	!!
16	!! REFERENCE(S) :
17	!!
18	!! SVN :
19	!! $HeadURL: svn://forge.ipsl.jussieu.fr/orchidee/branches/ORCHIDEE-DOFOCO/ORCHIDEE/src_stomate/stomate_prescribe.f90 $
20	!! $Date: 2013-01-04 16:50:56 +0100 (Fri, 04 Jan 2013) $
21	!! $Revision: 1126 $
22	!! \n
23	!_ ================================================================================================================================
24
25	MODULE function_library
26
27	! modules used:
28
29	!$ USE ioipsl_para
30	USE pft_parameters
31	USE constantes
32
33	IMPLICIT NONE
34
35	! private & public routines
36
37	PRIVATE
38	PUBLIC calculate_c0_alloc, wood_to_ba_eff, wood_to_ba, &
39	wood_to_qmheight, &
40	wood_to_qmdia, &
41	wood_to_volume, &
42	biomass_to_lai, &
43	lai_to_biomass, &
44	check_biomass_sync
45
46	INTERFACE biomass_to_lai
47	MODULE PROCEDURE biomass_to_lai_0d, biomass_to_lai_1d
48	END INTERFACE
49
50	INTERFACE lai_to_biomass
51	MODULE PROCEDURE lai_to_biomass_0d, lai_to_biomass_1d
52	END INTERFACE
53
54
55	CONTAINS
56
57
58
59	!! ================================================================================================================================
60	!! FUNCTION : calculate_c0_alloc
61	!!
62	!>\BRIEF Calculate the baseline root vs sapwood allocation
63	!!
64	!! DESCRIPTION : Calculates the baseline root vs sapwood allocation based on the
65	!! parameters of the pipe model (hydraulic conductivities) and the
66	!! turnover of the different components
67	!!
68	!! RECENT CHANGE(S): None
69	!!
70	!! RETURN VALUE : ::c0_alloc (m)
71	!!
72	!! REFERENCE(S) :
73	!!
74	!! FLOWCHART : None
75	!! \n
76	!_ ================================================================================================================================
77
78	FUNCTION calculate_c0_alloc(pts, pft, tau_eff_root, tau_eff_sap)
79
80	!! 0. Variable and parameter declaration
81
82	!! 0.1 Input variables
83
84	INTEGER(i_std) :: pts !! Pixel number (-)
85	INTEGER(i_std) :: pft !! PFT number (-)
86	REAL(r_std) :: tau_eff_root !! Effective longivety for leaves (days)
87	REAL(r_std) :: tau_eff_sap !! Effective longivety for leaves (days)
88
89	!! 0.2 Output variables
90
91	REAL(r_std) :: calculate_c0_alloc !! quadratic mean height (m)
92
93	!! 0.3 Modified variables
94
95	!! 0.4 Local variables
96	REAL(r_std) :: sapwood_density
97	REAL(r_std) :: qm_dia !! quadratic mean diameter (m)
98
99	!_ ================================================================================================================================
100
101	!! 1. Calculate c0_alloc
102	IF ( is_tree(pft) ) THEN
103
104	sapwood_density = deux * pipe_density(pft) / kilo_to_unit
105	calculate_c0_alloc = sqrt(k_root(pft)/k_sap(pft)tau_eff_sap/tau_eff_rootsapwood_density)
106
107	! Grasses and croplands
108	ELSE
109
110	!+++CHECK+++
111	! Simply copied the same formulation as for trees but note
112	! that the sapwood in trees vs grasses and crops has a very
113	! meaning. In grasses and crops is structural carbon to ensure
114	! that the allocation works. In trees it really is the sapwood
115	sapwood_density = deux * pipe_density(pft) / kilo_to_unit
116	calculate_c0_alloc = sqrt(k_root(pft)/k_sap(pft)tau_eff_sap/tau_eff_rootsapwood_density)
117	!+++++++++++
118
119	ENDIF ! is_tree(j)
120
121	END FUNCTION calculate_c0_alloc
122
123
124
125
126
127	!! ================================================================================================================================
128	!! FUNCTION : wood_to_ba_eff
129	!!
130	!>\BRIEF Calculate effective basal area from woody biomass making use of allometric relationships
131	!!
132	!! DESCRIPTION : Calculate basal area of an individual tree from the woody biomass of that tree making
133	!! use of allometric relationships. Effective basal area accounts for both above and below ground carbon
134	!! and is the basis for the application of the rule of Deleuze and Dhote.
135	!! (i) woodmass = tree_ff * pipe_densitybaheight
136	!! (ii) height = pipe_tune2 * sqrt(4/piba) * pipe_tune_3
137	!!
138	!! RECENT CHANGE(S): None
139	!!
140	!! RETURN VALUE : effective basal area (m2 ind-1)
141	!!
142	!! REFERENCE(S) :
143	!!
144	!! FLOWCHART : None
145	!! \n
146	!_ ================================================================================================================================
147
148	FUNCTION wood_to_ba_eff(biomass_temp, pft)
149
150	!! 0. Variable and parameter declaration
151
152	!! 0.1 Input variables
153
154	INTEGER(i_std) :: pft !! PFT number (-)
155	REAL(r_std), DIMENSION(:) :: biomass_temp !! Biomass of an individual tree within a circ
156	!! class @tex $(m^{2} ind^{-1})$ @endtex
157
158	!! 0.2 Output variables
159
160	REAL(r_std), DIMENSION(ncirc) :: wood_to_ba_eff !! Effective basal area of an individual tree within a circ
161	!! class @tex $(m^{2} ind^{-1})$ @endtex
162
163	!! 0.3 Modified variables
164
165	!! 0.4 Local variables
166	INTEGER(i_std) :: l !! Index
167	REAL(r_std) :: woodmass_ind !! Woodmass of an individual tree
168	!! @tex $(gC ind^{-1})$ @endtex
169
170	!_ ================================================================================================================================
171
172	!! 1. Calculate basal area from woodmass
173
174	IF ( is_tree(pft) ) THEN
175
176	DO l = 1,ncirc
177
178	! Woodmass of an individual tree
179	woodmass_ind = biomass_temp(isapabove) + biomass_temp(isapbelow) + &
180	biomass_temp(iheartabove) + biomass_temp(iheartbelow)
181
182	! Basal area of that individual (m2 ind-1)
183	wood_to_ba_eff(l) = (pi/4(woodmass_ind/(tree_ff(pft)pipe_density(pft)*pipe_tune2(pft))) &
184	(2./pipe_tune3(pft)))(pipe_tune3(pft)/(pipe_tune3(pft)+2))
185
186	ENDDO
187
188	ELSE
189
190	WRITE(numout,*) 'pft ',pft
191	CALL ipslerr_p (3,'wood_to_ba_eff', &
192	'wood_to_ba_eff is not defined for this PFT.', &
193	'See the output file for more details.','')
194
195	ENDIF
196
197	END FUNCTION wood_to_ba_eff
198
199
200
201	!! ================================================================================================================================
202	!! FUNCTION : wood_to_ba
203	!!
204	!>\BRIEF Calculate basal area from woody biomass making use of allometric relationships
205	!!
206	!! DESCRIPTION : Calculate basal area of an individual tree from the woody biomass of that tree making
207	!! use of allometric relationships given below. Here basal area is defined in line with its classical
208	!! forestry meaning.
209	!! (i) woodmass = tree_ff * pipe_densitybaheight
210	!! (ii) height = pipe_tune2 * sqrt(4/piba) * pipe_tune_3
211	!!
212	!! RECENT CHANGE(S): None
213	!!
214	!! RETURN VALUE : basal area (m2 ind-1)
215	!!
216	!! REFERENCE(S) :
217	!!
218	!! FLOWCHART : None
219	!! \n
220	!_ ================================================================================================================================
221
222	FUNCTION wood_to_ba(biomass_temp, pft)
223
224	!! 0. Variable and parameter declaration
225
226	!! 0.1 Input variables
227
228	INTEGER(i_std) :: pft !! PFT number (-)
229	REAL(r_std), DIMENSION(:) :: biomass_temp !! Biomass of an individual tree within a circ
230	!! class @tex $(m^{2} ind^{-1})$ @endtex
231
232	!! 0.2 Output variables
233
234	REAL(r_std) :: wood_to_ba !! Basal area of an individual tree within a circ
235	!! class @tex $(m^{2} ind^{-1})$ @endtex
236
237	!! 0.3 Modified variables
238
239	!! 0.4 Local variables
240	REAL(r_std) :: woodmass_ind !! Woodmass of an individual tree
241	!! @tex $(gC ind^{-1})$ @endtex
242
243	!_ ================================================================================================================================
244
245	!! 1. Calculate basal area from woodmass
246
247	IF ( is_tree(pft) ) THEN
248
249
250	! Woodmass of an individual tree
251	woodmass_ind = biomass_temp(iheartabove) + biomass_temp(isapabove)
252
253
254	! Basal area of that individual (m2 ind-1)
255	wood_to_ba = (pi/4(woodmass_ind/(tree_ff(pft)pipe_density(pft)*pipe_tune2(pft))) &
256	(2./pipe_tune3(pft)))(pipe_tune3(pft)/(pipe_tune3(pft)+2))
257
258
259	ELSE
260
261	WRITE(numout,*) 'pft ',pft
262	CALL ipslerr_p (3,'wood_to_ba', &
263	'wood_to_ba is not defined for this PFT.', &
264	'See the output file for more details.','')
265
266	ENDIF
267
268	END FUNCTION wood_to_ba
269
270
271
272
273
274	!! ================================================================================================================================
275	!! FUNCTION : wood_to_qmheight
276	!!
277	!>\BRIEF Calculate the quadratic mean height from the biomass
278	!!
279	!! DESCRIPTION : Calculates the quadratic mean height from the biomass
280	!!
281	!! RECENT CHANGE(S): None
282	!!
283	!! RETURN VALUE : ::qm_height (m)
284	!!
285	!! REFERENCE(S) :
286	!!
287	!! FLOWCHART : None
288	!! \n
289	!_ ================================================================================================================================
290
291	FUNCTION wood_to_qmheight(biomass_temp, ind, pft)
292
293	!! 0. Variable and parameter declaration
294
295	!! 0.1 Input variables
296
297	INTEGER(i_std) :: pft !! PFT number (-)
298	REAL(r_std), DIMENSION(nparts) :: biomass_temp !! Biomass of the leaves @tex $(gC m^{-2})$ @endtex
299	REAL(r_std), DIMENSION(ncirc) :: ind !! Number of individuals @tex $(m^{-2})$ @endtex
300
301
302	!! 0.2 Output variables
303
304	REAL(r_std) :: wood_to_qmheight !! quadratic mean height (m)
305
306	!! 0.3 Modified variables
307
308	!! 0.4 Local variables
309	REAL(r_std), DIMENSION(ncirc) :: circ_class_ba !! basal area for each circ_class @tex $(m^{2})$ @endtex
310	REAL(r_std) :: qm_dia !! quadratic mean diameter (m)
311
312	!_ ================================================================================================================================
313
314	!! 1. Calculate qm_height from the biomass
315	IF ( is_tree(pft) ) THEN
316
317	! Basal area at the tree level (m2 tree-1)
318	circ_class_ba(:) = wood_to_ba(biomass_temp(:),pft)
319
320	IF (SUM(ind(:)) .NE. zero) THEN
321
322	qm_dia = SQRT( 4/piSUM(circ_class_ba(:)ind(:))/SUM(ind(:)) )
323
324	ELSE
325
326	qm_dia = zero
327
328	ENDIF
329
330	wood_to_qmheight = pipe_tune2(pft)(qm_dia*pipe_tune3(pft))
331
332
333	! Grasses and croplands
334	ELSE
335
336	! Calculate height as a function of the leaf and structural biomass. Use structural
337	! biomass to make sure that the grasslands have a roughness length during the winter
338	! If the biomass increases, vegetation height will increase as well. Divide by
339	! ind(ipts,j) to obtain the height of the individual. biomass(ileaf) is in gC m-2
340	! whereas qm is the height of the individual.
341	IF (SUM(ind(:)) .NE. zero) THEN
342
343	wood_to_qmheight = (biomass_temp(ileaf) + biomass_temp(isapabove)) / &
344	SUM(ind(:)) * sla(pft) * lai_to_height(pft)
345
346	ELSE
347
348	wood_to_qmheight = zero
349
350	ENDIF
351
352	ENDIF ! is_tree(j)
353
354	END FUNCTION wood_to_qmheight
355
356
357
358
359
360
361	!! ================================================================================================================================
362	!! FUNCTION : wood_to_qmdia
363	!!
364	!>\BRIEF Calculate the quadratic mean diameter from the biomass
365	!!
366	!! DESCRIPTION : Calculates the quadratic mean diameter from the aboveground biomss
367	!!
368	!! RECENT CHANGE(S): None
369	!!
370	!! RETURN VALUE : ::qm_dia (m)
371	!!
372	!! REFERENCE(S) :
373	!!
374	!! FLOWCHART : None
375	!! \n
376	!_ ================================================================================================================================
377
378	FUNCTION wood_to_qmdia(biomass_temp, ind, pft)
379
380	!! 0. Variable and parameter declaration
381
382	!! 0.1 Input variables
383
384	INTEGER(i_std) :: pft !! PFT number (-)
385	REAL(r_std), DIMENSION(nparts) :: biomass_temp !! Biomass of the leaves @tex $(gC m^{-2})$ @endtex
386	REAL(r_std), DIMENSION(ncirc) :: ind !! Number of individuals @tex $(m^{-2})$ @endtex
387
388	!! 0.2 Output variables
389
390	REAL(r_std) :: wood_to_qmdia !! quadratic mean diameter (m)
391
392	!! 0.3 Modified variables
393
394	!! 0.4 Local variables
395	REAL(r_std), DIMENSION(ncirc) :: circ_class_ba !! basal area for each circ_class @tex $(m^{2})$ @endtex
396
397	!_ ================================================================================================================================
398
399	!! 1. Calculate qm_dia from the biomass
400	IF ( is_tree(pft) ) THEN
401
402	! Basal area at the tree level (m2 tree-1)
403	circ_class_ba(:) = wood_to_ba(biomass_temp(:),pft)
404
405	IF (SUM(ind(:)) .NE. zero) THEN
406
407	wood_to_qmdia = SQRT( 4/piSUM(circ_class_ba(:)ind(:))/SUM(ind(:)) )
408
409	ELSE
410
411	wood_to_qmdia = zero
412
413	ENDIF
414
415
416	! Grasses and croplands
417	ELSE
418
419	wood_to_qmdia = zero
420
421	ENDIF ! is_tree(pft)
422
423	END FUNCTION wood_to_qmdia
424
425
426	!! ================================================================================================================================
427	!! FUNCTION : wood_to_volume
428	!!
429	!>\BRIEF This allometric function computes volume as a function of
430	!! biomass at stand scale. Volume \f$(m^3 m^{-2}) = f(biomass (gC m^{-2}))\f$
431	!!
432	!! DESCRIPTION : None
433	!!
434	!! RECENT CHANGE(S): None
435	!!
436	!! RETURN VALUE : biomass_to_volume
437	!!
438	!! REFERENCE(S) : See above, module description.
439	!!
440	!! FLOWCHART : None
441	!! \n
442	!_ ================================================================================================================================
443
444	FUNCTION wood_to_volume(biomass,pft,branch_ratio,inc_branches)
445
446	!! 0. Variable and parameter declaration
447
448	!! 0.1 Input variables
449
450	REAL(r_std), DIMENSION(:) :: biomass !! Stand biomass @tex $(gC m^{-2})$ @endtex
451	REAL(r_std) :: branch_ratio !! Branch ratio of sap and heartwood biomass
452	!! unitless
453	INTEGER(i_std) :: pft !! Plant functional type (unitless)
454	INTEGER(i_std) :: inc_branches !! Include the branches in the volume calculation?
455	!! 0: exclude the branches from the volume calculation
456	!! (thus correct the biomass for the branch ratio)
457	!! 1: include the branches in the volume calculation
458	!! (thus use all aboveground biomass)
459
460
461
462	!! 0.2 Output variables
463
464	REAL(r_std) :: wood_to_volume !! The volume of wood per square meter
465	!! @tex $(m^3 m^{-2})$ @endtex
466
467	!! 0.3 Modified variables
468
469	!! 0.4 Local variables
470
471	REAL(r_std) :: woody_biomass !! Woody biomass at the stand level
472	!! @tex $(gC m^{-2})$ @endtex
473
474	!_ ================================================================================================================================
475
476	!! 1. Volume to biomass
477
478	! Woody biomass used in the calculation
479	IF (inc_branches .EQ. 0) THEN
480
481	woody_biomass=(biomass(isapabove)+biomass(iheartabove))*(un - branch_ratio)
482
483	ELSEIF (inc_branches .EQ. 1) THEN
484
485	woody_biomass=(biomass(isapabove)+biomass(iheartabove))
486
487	ELSE
488
489	ENDIF
490
491	! Wood volume expressed in m3 / m2
492	wood_to_volume = woody_biomass/(pipe_density(pft))
493
494	END FUNCTION wood_to_volume
495
496
497
498	!! ================================================================================================================================
499	!! FUNCTION : biomass_to_lai_0d
500	!!
501	!>\BRIEF Calculate the LAI based on the leaf biomass
502	!!
503	!! DESCRIPTION : Calculates the LAI of a PFT/grid square based on the leaf biomass
504	!!
505	!! RECENT CHANGE(S): None
506	!!
507	!! RETURN VALUE : ::LAI [m2 m{-2}]
508	!!
509	!! REFERENCE(S) :
510	!!
511	!! FLOWCHART : None
512	!! \n
513	!_ ================================================================================================================================
514
515	FUNCTION biomass_to_lai_0d(leaf_biomass, pft)
516
517	!! 0. Variable and parameter declaration
518
519	!! 0.1 Input variables
520
521	INTEGER(i_std) :: pft !! PFT number (-)
522	REAL(r_std) :: leaf_biomass !! Biomass of the leaves
523	!! @tex $(gC m^{-2})$ @endtex
524
525
526	!! 0.2 Output variables
527
528	REAL(r_std) :: biomass_to_lai_0d !! Leaf area index
529	!! @tex $(m^{2} m^{-2})$ @endtex
530
531	!! 0.3 Modified variables
532
533	!! 0.4 Local variables
534	REAL(r_std) :: impose_lai !! LAI read from run.def
535	!_ ================================================================================================================================
536
537	!! 1. Calculate the LAI from the leaf biomass
538	IF(sla_dyn) THEN
539	biomass_to_lai_0d = log(1.+ ext_coeff_N(pft)* leaf_biomass * slainit(pft))/(ext_coeff_N(pft))
540	ELSE
541	biomass_to_lai_0d = leaf_biomass * sla(pft)
542	ENDIF
543
544	!!$ !+++++++++ TEMP ++++++++++
545	!!$ ! This is a perfect place to hack the code to make it run with
546	!!$ ! constant lai
547	!!$ WRITE(numout,*) 'WARNING ERROR: Using fake lai values for testing!'
548	!!$ biomass_to_lai=3.79052
549	!!$ !+++++++++++++++++++++++++
550
551	!+++++++ TEMP ++++++++++
552	! This code is only used evaluation of the performance of the multi-layer energy budget.
553	! To reduce the complexity of the tests we want to impose the LAI and its vertical distribution.
554	! The solution is not very elegant but it works.
555	! IF (ld_fake_height) THEN
556	! In order to imposed lai, we read the TOTAL_LAI from run.def
557	! CALL getin_p('TOTAL_LAI', impose_lai)
558	! This part of code reset the sla vale to match which alow modeled LAI equal to TOTAL LAI.
559	! Althought this is ugly way to match the modeled LAI and impose LAI.
560	! You probably need to go to your ORCHIDEE out put file to find out the suitable SLA value
561	! and reset it agin in the run.def.
562	! So, we impose LAI & structure for a quick testing the performance of multilayer energy budget
563	! without changing the leaf_biomass.
564	! IF ( leaf_biomass .GT. 0.0) THEN
565	! sla(pft)=impose_lai/leaf_biomass
566	! WRITE(numout,'(A,F20.8)') 'USE A FAKE SLA BASED ON imposed LAI/LEAFMASS=', sla(pft)
567	! ENDIF
568	! biomass_to_lai=leaf_biomass*sla(pft)
569	! ENDIF
570	!++++++++++++++++++++++++
571
572	END FUNCTION biomass_to_lai_0d
573
574	FUNCTION biomass_to_lai_1d(leaf_biomass, npts, pft)
575
576	!! 0. Variable and parameter declaration
577
578	!! 0.1 Input variables
579	INTEGER(i_std) :: npts
580	INTEGER(i_std) :: pft !! PFT number (-)
581	REAL(r_std), DIMENSION(npts) :: leaf_biomass !! Biomass of the leaves
582	!! @tex $(gC m^{-2})$ @endtex
583
584
585	!! 0.2 Output variables
586
587	REAL(r_std), DIMENSION(npts) :: biomass_to_lai_1d !! Leaf area index
588	!! @tex $(m^{2} m^{-2})$ @endtex
589	!_ ================================================================================================================================
590
591	!! 1. Calculate the LAI from the leaf biomass
592	IF(sla_dyn) THEN
593	biomass_to_lai_1d = log(1.+ ext_coeff_N(pft)* leaf_biomass * slainit(pft))/(ext_coeff_N(pft))
594	ELSE
595	biomass_to_lai_1d = leaf_biomass * sla(pft)
596	ENDIF
597
598	END FUNCTION biomass_to_lai_1d
599
600
601	FUNCTION lai_to_biomass_0d(lai, pft)
602
603	!! 0. Variable and parameter declaration
604
605	!! 0.1 Input variables
606	INTEGER(i_std) :: pft !! PFT number (-)
607	REAL(r_std) :: lai !! Leaf Area Index
608	!! @tex $(gC m^{-2})$ @endtex
609
610
611	!! 0.2 Output variables
612
613	REAL(r_std) :: lai_to_biomass_0d !! Leaf area index
614	!! @tex $(m^{2} m^{-2})$ @endtex
615	!_ ================================================================================================================================
616
617	!! 1. Calculate the LAI from the leaf biomass
618	IF(sla_dyn) THEN
619	lai_to_biomass_0d = ( exp(lai*ext_coeff_N(pft)) - 1.) / &
620	(ext_coeff_N(pft) * slainit(pft))
621	ELSE
622	lai_to_biomass_0d = lai / sla(pft)
623	ENDIF
624
625	END FUNCTION lai_to_biomass_0d
626
627
628	FUNCTION lai_to_biomass_1d(lai, npts, pft)
629
630	!! 0. Variable and parameter declaration
631
632	!! 0.1 Input variables
633	INTEGER(i_std) :: npts
634	INTEGER(i_std) :: pft !! PFT number (-)
635	REAL(r_std), DIMENSION(npts) :: lai !! Leaf Area Index
636	!! @tex $(gC m^{-2})$ @endtex
637
638
639	!! 0.2 Output variables
640
641	REAL(r_std), DIMENSION(npts) :: lai_to_biomass_1d !! Leaf area index
642	!! @tex $(m^{2} m^{-2})$ @endtex
643	!_ ================================================================================================================================
644
645	!! 1. Calculate the LAI from the leaf biomass
646	IF(sla_dyn) THEN
647	lai_to_biomass_1d = ( exp(lai*ext_coeff_N(pft)) - 1.) / &
648	(ext_coeff_N(pft) * slainit(pft))
649	ELSE
650	lai_to_biomass_1d = lai / sla(pft)
651	ENDIF
652
653	END FUNCTION lai_to_biomass_1d
654
655
656
657
658	!! ================================================================================================================================
659	!! SUBROUTINE : check_biomass_sync
660	!!
661	!>\BRIEF
662	!!
663	!! DESCRIPTION :
664	!! RECENT CHANGE(S): None
665	!!
666	!! MAIN OUTPUT VARIABLE(S):
667	!!
668	!! REFERENCE(S) : None
669	!!
670	!! FLOWCHART : None
671	!! \n
672	!_ ================================================================================================================================
673	SUBROUTINE check_biomass_sync ( check_point, npts, biomass, &
674	circ_class_biomass, circ_class_n , ind, &
675	lsync, bm_sync)
676
677	!! 0. Variable and parameter description
678
679	!! 0.1 Input variables
680	INTEGER(i_std), INTENT(in) :: npts !! Domain size (unitless)
681	REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(in) :: circ_class_biomass !! Biomass of the componets of the model
682	!! tree within a circumference
683	!! class @tex $(gC ind^{-1})$ @endtex
684	REAL(r_std), DIMENSION(:,:,:), INTENT(in) :: circ_class_n !! Number of individuals in each circ class
685	!! @tex $(m^{-2})$ @endtex
686	REAL(r_std), DIMENSION(:,:,:,:), INTENT(in) :: biomass !! Stand level biomass
687	!! @tex $(gC m^{-2})$ @endtex
688	CHARACTER(*),INTENT(in) :: check_point !! A flag to indicate at which
689	!! point in the code we're doing
690	!! this check
691	REAL(r_std), DIMENSION(:,:), INTENT(in) :: ind !! Density of individuals
692	!! @tex $(m^{-2})$ @endtex
693
694	!! 0.2 Output variables
695	LOGICAL,INTENT(out) :: lsync
696	REAL(r_std), DIMENSION(:,:,:), INTENT(out) :: bm_sync !! The difference betweeen the
697	!! biomass in the circ_classes and
698	!! the total biomass
699	!! @tex $(gC m^{-2})$ @endtex
700	!! 0.3 Modified variables
701
702	!! 0.4 Local variables
703	INTEGER :: iele,ipts,ivm,ipar,icir
704	REAL(r_std) :: total_circ_class_biomass
705	REAL(r_std),DIMENSION(ncirc) :: tree_size
706	LOGICAL :: lnegative
707
708	!_ ================================================================================================================================
709
710	lsync=.TRUE.
711	lnegative=.FALSE.
712
713	bm_sync(:,:,:)=zero
714
715	!++++++ TEMP ++++++
716	! We gain 5-10% speed by skipping this routine
717
718	!++++++++++++
719
720	! Check to see if the biomass is not equal to the total biomass
721	! in circ_class_biomass anywhere.
722	DO ipts=1,npts
723
724	DO ivm=1,nvm
725
726	! Only woody PFTs have circumference classes therefore
727	! only woody PFTs need to be syncronized
728	IF(.NOT. lbypass_cc)THEN
729	IF(is_tree(ivm)) THEN
730	tree_size(:)=zero
731	DO icir=1,ncirc
732	tree_size(icir)=SUM(circ_class_biomass(ipts,ivm,icir,:,1))
733	ENDDO
734	DO icir=2,ncirc
735	IF(tree_size(icir) .LT. tree_size(icir-1)-min_stomate)THEN
736	WRITE(numout,*) 'ERROR: stopping in sync'
737	WRITE(numout,*) check_point
738	WRITE(numout,*) 'ipts,ivm: ',ipts,ivm
739	WRITE(numout,*) 'tree_size(icir), tree_size(icir-1), ',&
740	tree_size(icir), tree_size(icir-1), tree_size(icir) - tree_size(icir-1)
741	WRITE(numout,*) 'icir, tree_size: ',icir, tree_size(:)
742	!+++ TEMP +++
743	!This would not STOP the ORCHIDEE beacause the mass balance is due to imposed LAI
744	IF(ld_fake_height) THEN
745	CALL ipslerr_p (2,'check_biomass_sync', &
746	'The size of the trees in the circ class are not monotonically increasing!',&
747	'Look in the output file for more details.',&
748	'')
749	ELSE
750	CALL ipslerr_p (3,'check_biomass_sync', &
751	'The size of the trees in the circ class are not monotonically increasing!',&
752	'Look in the output file for more details.',&
753	'')
754	ENDIF
755	!++++++++++++
756	ENDIF
757	ENDDO
758	ENDIF
759	ENDIF
760
761	DO iele=1,icarbon
762
763	DO ipar=1,nparts
764
765	total_circ_class_biomass=zero
766	DO icir=1,ncirc
767
768	total_circ_class_biomass=total_circ_class_biomass+&
769	circ_class_biomass(ipts,ivm,icir,ipar,iele)*circ_class_n(ipts,ivm,icir)
770
771	! Check as well to see if our biomass is ever negative.
772	! It really should not be.
773	IF(circ_class_biomass(ipts,ivm,icir,ipar,iele) .LT. -min_stomate)THEN
774
775	lnegative=.TRUE.
776	WRITE(numout,) '!**********************************'
777	WRITE(numout,*) 'Error: Negative biomass component!'
778	WRITE(numout,*) 'Check point: ',TRIM(check_point)
779	WRITE(numout,*) 'circ_class_biomass(ipts,ivm,icir,ipar,iele) ',&
780	circ_class_biomass(ipts,ivm,icir,ipar,iele)
781	WRITE(numout,'(A,5I5)') 'ipts,ivm,icir,ipar,iele',ipts,ivm,icir,ipar,iele
782	WRITE(numout,) '!**********************************'
783
784	ENDIF
785	ENDDO
786
787	IF(ABS(biomass(ipts,ivm,ipar,iele) - &
788	total_circ_class_biomass) .GT. sync_threshold)THEN
789
790	WRITE(numout,) '!**********************************'
791	WRITE(numout,*) 'Biomass and circ_class_biomass are not equal!'
792	WRITE(numout,*) 'Check point: ',TRIM(check_point)
793	WRITE(numout,100) 'biomass(ipts,ivm,ipar,iele) ',&
794	biomass(ipts,ivm,ipar,iele)
795	WRITE(numout,100) 'total_circ_class_biomass ',&
796	total_circ_class_biomass
797	WRITE(numout,100) 'Difference: ',&
798	ABS(biomass(ipts,ivm,ipar,iele) - total_circ_class_biomass)
799	WRITE(numout,*) 'ipts,ivm,ipar,iele',ipts,ivm,ipar,iele
800	WRITE(numout,) '!**********************************'
801	100 FORMAT(A,E20.10)
802	lsync=.FALSE.
803
804	ENDIF
805
806	ENDDO
807
808	! we are not going to save the biomass for every component right now,
809	! just the total
810	bm_sync(ipts,ivm,iele)=zero
811
812	DO ipar=1,nparts
813
814
815	DO icir=1,ncirc
816
817	bm_sync(ipts,ivm,iele)=bm_sync(ipts,ivm,iele)+&
818	circ_class_biomass(ipts,ivm,icir,ipar,iele)*circ_class_n(ipts,ivm,icir)
819	ENDDO ! ncirc
820
821	ENDDO ! nparts
822
823	bm_sync(ipts,ivm,iele)=ABS(bm_sync(ipts,ivm,iele)-&
824	SUM(biomass(ipts,ivm,:,iele)))
825
826	ENDDO ! nelements
827
828
829
830	ENDDO ! loop over PFTs
831
832	ENDDO ! loop over points
833
834	!---TEMP---
835	IF(ld_biomass)THEN
836	WRITE(numout,*) 'Check point: ',TRIM(check_point)
837	WRITE(numout,*) 'test_pft, test_grid: ',test_pft,test_grid
838	WRITE(numout,*) 'biomass (ileaf), ', biomass(test_grid,test_pft,ileaf,icarbon)
839	WRITE(numout,*) 'biomass (iwood), ', biomass(test_grid,test_pft,isapabove,icarbon) + &
840	biomass(test_grid,test_pft,isapbelow,icarbon) + biomass(test_grid,test_pft,iheartabove,icarbon) + &
841	biomass(test_grid,test_pft,iheartbelow,icarbon)
842	WRITE(numout,*) 'biomass (iroot), ', biomass(test_grid,test_pft,iroot,icarbon)
843	WRITE(numout,'(A,20F14.6)') 'biomassHHH, ',biomass(test_grid,test_pft,:,icarbon)
844	DO icir=1,ncirc
845	WRITE(numout,'(A,I1,20F14.6)') 'ccbiomass',icir,circ_class_biomass(test_grid,test_pft,icir,:,icarbon)
846	ENDDO
847	WRITE(numout,*) 'circ_class_biomass, ',&
848	SUM (SUM(circ_class_biomass(test_grid,test_pft,:,:,icarbon),2) * &
849	circ_class_n(test_grid,test_pft,:))
850	WRITE(numout,*) 'circ_class_n, ', SUM(circ_class_n(test_grid,test_pft,:))
851	WRITE(numout,*) 'circ_class_n(:), ', circ_class_n(test_grid,test_pft,:)
852	WRITE(numout,*) 'ind, ', ind(test_grid,test_pft)
853	ENDIF
854
855	!!$ !----------
856
857	IF(.NOT. lsync) THEN
858	WRITE(numout,*) 'ERROR: stopping in sync #2'
859	WRITE(numout,*) 'Stopping'
860	CALL ipslerr_p (3,'check_biomass_sync', &
861	'circ_class_biomass*circ_class_n is not equal to the total biomass',&
862	'Look in the output file for more details.',&
863	'')
864
865	ENDIF
866	IF(lnegative) THEN
867	WRITE(numout,*) 'ERROR: negative biomass'
868	WRITE(numout,*) 'Stopping'
869	CALL ipslerr_p (3,'check_biomass_sync', &
870	'One of the biomass pools is negative!',&
871	'Look in the output file for more details.',&
872	'')
873	ENDIF
874
875	END SUBROUTINE check_biomass_sync
876
877	END MODULE function_library
878

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

Download in other formats: