source: branches/publications/ORCHIDEE-CN-P-MIMICS_r7301/ORCHIDEE/src_parameters/function_library.f90 @ 7346

Last change on this file since 7346 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
25MODULE 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_root*sapwood_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_root*sapwood_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_density*ba*height
125!! (ii) height = pipe_tune2 * sqrt(4/pi*ba) ** 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_density*ba*height
199!! (ii) height = pipe_tune2 * sqrt(4/pi*ba) ** 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/pi*SUM(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/pi*SUM(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 m**3 / m**2
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 [m**2 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,*) '!***********************************'
707100                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
783END MODULE function_library
784
Note: See TracBrowser for help on using the repository browser.