Context Navigation

function_library.f90

Last change on this file was 4734, checked in by daniel.goll, 7 years ago

I out commented code related prescribing lai and height
for which the switch was not properly implemented and caused
issue on obelix

DSG

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

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

Context Navigation

source: branches/publications/ORCHIDEE_CN-P_v1.2_r5986/ORCHIDEE/src_parameters/function_library.f90

Download in other formats: