source: branches/publications/ORCHIDEE_CAN_NHA/src_parameters/function_library.f90 @ 8066

Last change on this file since 8066 was 7235, checked in by yitong.yao, 4 years ago

update src_parameters
File size: 98.5 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_array, wood_to_ba_eff, wood_to_ba, &
39       wood_to_height_eff, wood_to_height, wood_to_qmheight, wood_to_dia_eff, &
40       wood_to_dia, wood_to_qmdia, wood_to_circ, wood_to_cn_array, &
41       wood_to_cn, wood_to_cn_eff, wood_to_cv, wood_to_cv_eff, wood_to_volume, &
42       biomass_to_lai, Nmax, Nmaxyield, Nminyield, distribute_mortality_biomass, &
43       check_biomass_sync, biomass_to_coupled_lai, stem_volume_f, vxmax_f, vrmax_f,  &
44       root_volume_f, vmax_l_f, canopy_leaf_area_f, grav_potential_f, stomatal_conductance_hydro_f,&
45       kstem_f, kroot_f, kleaf_f, kupper_f,kupper_check,P50_vary_class, distribute_mortality_biomass_yyt,distribute_mortality_biomass_new 
46  CONTAINS
47
48
49
50!! ================================================================================================================================
51! tzjh hydraulic architecture function
52
53! ###########################################################################
54! Define water storage and gravitational potential functions
55! ###########################################################################
56!Stem volume function: Calculate stem volume (m^3) from height and dbh
57FUNCTION stem_volume_f(height,dbh)
58REAL, INTENT(IN) :: height, dbh
59REAL :: stem_volume_f
60   stem_volume_f = height*pi*(dbh/2)**2
61END FUNCTION stem_volume_f
62
63!Stem water storage function: Calculate maximum stem water storage volume in
64!mmol from tree size and water content
65FUNCTION vxmax_f(stem_volume,w_density_stem)
66REAL, INTENT(IN) :: stem_volume,w_density_stem
67REAL :: vxmax_f
68   vxmax_f = stem_volume*w_density_stem
69END FUNCTION vxmax_f
70
71!Root water storage function
72!Calculate maximum root water storage volume in mmol from the root water
73!content, ratio of root to shoot investment and root density     
74FUNCTION vrmax_f(stem_volume, wood_density, root_shoot_ratio, rwc_root)
75REAL, INTENT(IN) :: stem_volume, wood_density, root_shoot_ratio, rwc_root
76REAL :: stem_mass, root_mass,vrmax_f
77   stem_mass = stem_volume * wood_density
78   root_mass = stem_mass *root_shoot_ratio
79   vrmax_f=root_mass*rwc_root
80END FUNCTION vrmax_f
81
82! Root volume function Calculate the volume of the roots in m3 from the stem
83! mass, ratio of root to
84! shoot investment, and root density
85FUNCTION root_volume_f(stem_volume, wood_density, root_shoot_ratio, root_density)
86REAL, INTENT(IN) :: stem_volume, wood_density, root_shoot_ratio, root_density
87REAL :: root_volume_f, stem_mass, root_mass
88        stem_mass = stem_volume*wood_density !g
89        root_mass = stem_mass*root_shoot_ratio !#g
90        root_volume_f = (root_mass/root_density)/(100*100*100) !m3
91END FUNCTION root_volume_f
92
93!Leaf water storage function
94!Calculate the maximum water volume of all of the leaves in mmol from the tree
95!size and leaf dry matter content
96FUNCTION vmax_l_f(dbh, LDMC)
97REAL, INTENT(IN) :: dbh, LDMC
98REAL :: vmax_l_f,leaf_dry_mass, leaf_fresh_mass, water_mass, water_mass_mmol
99   !calculates the leaf dry mass (in kg) from the diameter of the trunk (in
100   !cm), according to the allometry equations in Djomo et al. 2010 Forest
101   !Ecology and Management
102   leaf_dry_mass = -0.1009 + 0.0626*(dbh*100) + 0.0027*(dbh*100)**2
103   !calculates leaf fresh mass (in kg) from the leaf dry matter content (g
104   !dry mass g-1 fresh mass)
105   leaf_fresh_mass = leaf_dry_mass*(1/LDMC)
106   !calculates the mass of water in the leaf (in kg)
107   water_mass = leaf_fresh_mass - leaf_dry_mass
108   ! converts water mass into mmol for consistency with other units
109   water_mass_mmol =  water_mass*1000*1000/18
110   vmax_l_f = water_mass_mmol
111END FUNCTION vmax_l_f
112
113! Canopy leaf area function: Calculate the leaf area of the canopy from the tree
114! size and SLA
115FUNCTION canopy_leaf_area_f(dbh, sla_hydro)
116REAL, INTENT(IN) :: dbh, sla_hydro
117REAL :: leaf_dry_mass, canopy_leaf_area_f
118   leaf_dry_mass = -0.1009 + 0.0626*(dbh*100) + 0.0027*(dbh*100)**2
119   canopy_leaf_area_f = leaf_dry_mass*sla_hydro
120END FUNCTION canopy_leaf_area_f
121
122! Gravitational potential function
123! Calculate how much of a water potential gradient is needed to pull water up
124! along the height of the tree
125! This assumes that the stem capacitor is at half the total tree height, and all
126! of the leaves are at the tree height
127! This also treats the distance the water has to be pulled up from the root to
128! the ground level as negligible- I think quantitatively it should be small
129! compared to the effect of stem height, and I think it would be hard to
130! accurately model this, since it depends on the spread of the roots, the soil
131! depth at which most water uptake occurs, the root length to mass ratio- but we
132! should talk about this
133FUNCTION grav_potential_f(height)
134REAL, INTENT(IN) :: height
135REAL ::  grav_potential_f
136   grav_potential_f = (height*1000*9.81)/1000000
137END FUNCTION grav_potential_f
138
139
140! Define physiology functions
141! #####################################################################################
142! Stomatal conductance function for hydraulics
143!FUNCTION stomatal_conductance_hydro_f(gmax, light, klight, VPD, gpsi, psi_leaf, gpsi_50, gmin)
144!REAL, INTENT(IN) :: gmax, light, klight, VPD, gpsi, psi_leaf, gpsi_50, gmin
145!REAL :: stomatal_conductance_hydro_f
146!   stomatal_conductance_hydro_f =  (gmax*(light/(ligh + klight))*(EXP(-VPD))/(1 + EXP(gpsi*(psi_leaf +gpsi_50)))) + gmin
147!END FUNCTION stomatal_conductance_hydro_f
148
149FUNCTION stomatal_conductance_hydro_f(gmax, light,klight, gpsi, psi_leaf, gpsi_50, gmin)
150REAL, INTENT(IN) :: gmax, light, klight, gpsi, psi_leaf, gpsi_50, gmin
151REAL :: stomatal_conductance_hydro_f
152   stomatal_conductance_hydro_f =  (gmax*(light/(light + klight)))/(1 + EXP(gpsi*(psi_leaf +gpsi_50))) + gmin
153END FUNCTION stomatal_conductance_hydro_f
154
155! Stem conductivity (cavitation) function
156! Calculates the stem conductivity for the next timestep using the psi_stem for
157! the next timestep
158FUNCTION kstem_f(kmax_stem, P50_stem, a_stem, psi_stem)
159REAL, INTENT(IN) :: kmax_stem, P50_stem, a_stem, psi_stem
160REAL :: kstem_f
161   kstem_f = kmax_stem/(1 + exp(a_stem*(psi_stem + P50_stem)))
162END FUNCTION kstem_f
163
164! Root conductivity (cavitation) function
165! Calculates the root conductivity for the next timestep using the psi_root for
166! the next timestep
167FUNCTION kroot_f(kmax_root, P50_root, a_root, psi_root)
168REAL, INTENT(IN) :: kmax_root, P50_root, a_root, psi_root
169REAL :: kroot_f
170   kroot_f = kmax_root/(1 + EXP(a_root*(psi_root + P50_root)))
171END FUNCTION kroot_f
172
173! Leaf conductivity (cavitation) function
174! Calculates the leaf conductivity for the next timestep using the psi_leaf for
175! the next timestep
176FUNCTION kleaf_f(kmax_leaf, P50_leaf, a_leaf, psi_leaf)
177REAL, INTENT(IN) :: kmax_leaf, P50_leaf, a_leaf, psi_leaf
178REAL :: kleaf_f
179   kleaf_f = kmax_leaf/(1 + EXP(a_leaf*(psi_leaf + P50_leaf)))
180END FUNCTION kleaf_f
181
182! Plant conductivity function
183! Optimized with the psi_leaf function to calculate k_leaf and psi_leaf
184! simultaneously, because these variables are functions of each other. 
185FUNCTION kupper_f(kmax_leaf, P50_leaf, a_leaf, psi_leaf, kstem)
186REAL, INTENT(IN) :: kmax_leaf, P50_leaf, a_leaf, psi_leaf, kstem
187REAL :: kleaf, kupper_f
188   kleaf = kmax_leaf/(1 + EXP(a_leaf*(psi_leaf + P50_leaf)))
189   kupper_f = 1/((1/kleaf) + (1/(2*kstem)))
190END FUNCTION kupper_f
191
192FUNCTION kupper_check(kmax_leaf, P50_leaf, a_leaf, psi_leaf, kstem)
193REAL, INTENT(IN) :: kmax_leaf, P50_leaf, a_leaf, psi_leaf, kstem
194REAL :: kleaf, kupper_check
195   kleaf = kmax_leaf/(1 + EXP(a_leaf*(psi_leaf + P50_leaf)))
196   kupper_check = 1/((1/kleaf) + (1/(2*kstem)))
197   print*,'IIIIIIIIIII check', kmax_leaf,a_leaf,psi_leaf,P50_leaf,kleaf,kupper_check
198END FUNCTION kupper_check
199
200
201
202
203
204
205!! ================================================================================================================================
206
207
208!! ================================================================================================================================
209!! FUNCTION    : biomass_to_coupled_lai
210!!
211!>\BRIEF        Calculate the coupled_LAI based on biomass and veget of each pft
212!!
213!! DESCRIPTION : Calculates the lai that coupled with the atmosphere 
214!!
215!!             
216!!
217!! RECENT CHANGE(S): None
218!!
219!! RETURN VALUE : ::coupled_lai (m2/m2)
220!!
221!! REFERENCE(S) :
222!!
223!! FLOWCHART    : None
224!! \n
225!_ ================================================================================================================================
226   
227  FUNCTION biomass_to_coupled_lai(biomass_leaf, veget, pft)
228
229 !! 0. Variable and parameter declaration
230
231    !! 0.1 Input variables
232    INTEGER(i_std)                                       :: pft                    !! PFT number (-)
233    REAL(r_std)                                          :: veget                  !! 1-Pgap or the vegetation cover
234    REAL(r_std)                                          :: biomass_leaf           !! Biomass of an individual tree within a circ
235                                                                                   !! class @tex $(m^{2} ind^{-1})$ @endtex
236
237    !! 0.2 Output variables         
238    REAL(r_std)                                          :: biomass_to_coupled_lai !! The fraction of the LAI that is assumed to interact
239                                                                                   !! with the atmosphere @tex $(m^{2} m^{-2})$ @endtex
240    !! 0.3 Modified variables
241
242    !! 0.4 Local variables
243    REAL(r_std)                                          :: lai                    !! Leaf area index
244                                                                                   !! @tex $(m^{2} m^{-2})$ @endtex               
245!_ ================================================================================================================================
246 
247    lai = biomass_to_lai(biomass_leaf, pft) 
248   
249    !+++CHECK+++
250    ! In a closed canopy not all the leaves fully interact with the
251    ! atmosphere because leaves can shelter each other. In a more
252    ! open canopy most leaves can interact with the atmosphere.
253    ! For the moment we are not clear how to calculate the fraction of
254    ! the canopy that is coupled to the atmosphere. Also, based on the
255    ! simulated evapotranspiration we need the entire canopy to be
256    ! coupled to the atmosphere to obtain simulations which are close
257    ! to the observations (Jung's upscaled product).
258    IF (veget .GT. min_sechiba) THEN
259     
260       biomass_to_coupled_lai = lai   
261
262    ELSE
263   
264       ! There are so many gaps that veget is extremely small. It is fair to
265       ! assume that the whole canopy is coupled to the atmosphere.
266       biomass_to_coupled_lai = lai
267     
268    ENDIF
269
270  END FUNCTION biomass_to_coupled_lai
271 
272 
273!! ================================================================================================================================
274!! FUNCTION     : calculate_c0_alloc
275!!
276!>\BRIEF        Calculate the baseline root vs sapwood allocation
277!!
278!! DESCRIPTION : Calculates the baseline root vs sapwood allocation based on the
279!! parameters of the pipe model (hydraulic conductivities) and the
280!! turnover of the different components             
281!!
282!! RECENT CHANGE(S): None
283!!
284!! RETURN VALUE : ::c0_alloc (m)
285!!
286!! REFERENCE(S) :
287!!
288!! FLOWCHART    : None
289!! \n
290!_ ================================================================================================================================
291   
292  FUNCTION calculate_c0_alloc(pts, pft, tau_eff_root, tau_eff_sap)
293
294 !! 0. Variable and parameter declaration
295
296    !! 0.1 Input variables
297
298    INTEGER(i_std)                             :: pts               !! Pixel number (-)
299    INTEGER(i_std)                             :: pft               !! PFT number (-)
300    REAL(r_std)                                :: tau_eff_root      !! Effective longivety for leaves (days)
301    REAL(r_std)                                :: tau_eff_sap       !! Effective longivety for leaves (days)
302   
303    !! 0.2 Output variables
304               
305    REAL(r_std)                                :: calculate_c0_alloc  !! quadratic mean height (m)
306
307    !! 0.3 Modified variables
308
309    !! 0.4 Local variables
310    REAL(r_std)                                :: sapwood_density
311    REAL(r_std)                                :: qm_dia            !! quadratic mean diameter (m)
312
313!_ ================================================================================================================================
314 
315    !! 1. Calculate c0_alloc
316    IF ( is_tree(pft) ) THEN
317
318       sapwood_density = deux * pipe_density(pft) / kilo_to_unit
319       calculate_c0_alloc = sqrt(k_root(pft)/k_sap(pft)*tau_eff_sap/tau_eff_root*sapwood_density)
320
321    ! Grasses and croplands   
322    ELSE
323
324       !+++CHECK+++
325       ! Simply copied the same formulation as for trees but note
326       ! that the sapwood in trees vs grasses and crops has a very
327       ! meaning. In grasses and crops is structural carbon to ensure
328       ! that the allocation works. In trees it really is the sapwood
329       sapwood_density = deux * pipe_density(pft) / kilo_to_unit
330       calculate_c0_alloc = sqrt(k_root(pft)/k_sap(pft)*tau_eff_sap/tau_eff_root*sapwood_density)
331       !+++++++++++
332
333    ENDIF ! is_tree(j)
334
335 END FUNCTION calculate_c0_alloc
336
337
338
339!! ================================================================================================================================
340!! FUNCTION     : wood_to_ba_eff_array
341!!
342!>\BRIEF        Calculate the effective basal area from woody biomass making use of allometric relationships.
343!!
344!! DESCRIPTION : Calculate basal area of an individual tree from the woody biomass of that tree making
345!! use of allometric relationships. Effective basal area accounts for both above and below ground carbon
346!! and is the basis for the application of the rule of Deleuze and Dhote
347!!
348!! RECENT CHANGE(S): None
349!!
350!! RETURN VALUE : effective basal area (m2 ind-1)
351!!
352!! REFERENCE(S) :
353!!
354!! FLOWCHART    : None
355!! \n
356!_ ================================================================================================================================
357   
358  FUNCTION wood_to_ba_eff_array(biomass_temp, npts, pft)
359
360 !! 0. Variable and parameter declaration
361
362    !! 0.1 Input variables
363
364    INTEGER(i_std)                                          :: pft                  !! PFT number (-)
365    INTEGER(i_std)                                          :: npts                 !! Number of pixels
366    REAL(r_std), DIMENSION(:,:,:)                           :: biomass_temp         !! Biomass of an individual tree within a circ
367                                                                                    !! class @tex $(m^{2} ind^{-1})$ @endtex
368
369    !! 0.2 Output variables
370               
371    REAL(r_std), DIMENSION(npts,ncirc)                      :: wood_to_ba_eff_array !! Effective basal area of an individual tree within a circ
372                                                                                    !! class @tex $(m^{2} ind^{-1})$ @endtex
373
374    !! 0.3 Modified variables
375
376    !! 0.4 Local variables
377    INTEGER(i_std)                                          :: l                    !! Index
378    REAL(r_std), DIMENSION(npts)                            :: woodmass_ind         !! Woodmass of an individual tree
379                                                                                    !! @tex $(gC ind{-1})$ @endtex
380
381!_ ================================================================================================================================
382 
383 !! 1. Calculate effective basal area from woodmass
384
385    IF ( is_tree(pft) ) THEN
386       
387       DO l = 1,ncirc
388         
389          ! Woodmass of an individual tree. Note that for the effective basal area
390          ! both the above and belowground biomass are used.
391          woodmass_ind(:) = biomass_temp(:,l,isapabove) + biomass_temp(:,l,isapbelow) + &
392             biomass_temp(:,l,iheartabove) + biomass_temp(:,l,iheartbelow)
393
394          ! Basal area of that individual (m2 ind-1) 
395          wood_to_ba_eff_array(:,l) = (pi/4*(woodmass_ind(:)/&
396               (tree_ff(pft)*pipe_density(pft)*pipe_tune2(pft))) &
397               **(2./pipe_tune3(pft)))**(pipe_tune3(pft)/(pipe_tune3(pft)+2))
398                 
399       ENDDO
400
401    ELSE
402
403       WRITE(numout,*) 'function wood_to_ba_eff_array is not defined for this PFT'
404       WRITE(numout,*) 'pft ',pft
405       CALL ipslerr_p (3,'wood_to_ba_eff_array', &
406            'wood_to_ba_eff_array is not defined for this PFT.', &
407            'See the output file for more details.','')
408
409
410    ENDIF
411
412  END FUNCTION wood_to_ba_eff_array
413
414
415!! ================================================================================================================================
416!! FUNCTION     : wood_to_ba_eff
417!!
418!>\BRIEF        Calculate effective basal area from woody biomass making use of allometric relationships
419!!
420!! DESCRIPTION :  Calculate basal area of an individual tree from the woody biomass of that tree making
421!! use of allometric relationships. Effective basal area accounts for both above and below ground carbon
422!! and is the basis for the application of the rule of Deleuze and Dhote.
423!! (i) woodmass = tree_ff * pipe_density*ba*height
424!! (ii) height = pipe_tune2 * sqrt(4/pi*ba) ** pipe_tune_3 
425!!
426!! RECENT CHANGE(S): None
427!!
428!! RETURN VALUE : effective basal area (m2 ind-1)
429!!
430!! REFERENCE(S) :
431!!
432!! FLOWCHART    : None
433!! \n
434!_ ================================================================================================================================
435 
436 FUNCTION wood_to_ba_eff(biomass_temp, pft)
437
438 !! 0. Variable and parameter declaration
439
440    !! 0.1 Input variables
441
442    INTEGER(i_std)                                          :: pft                  !! PFT number (-)
443    REAL(r_std), DIMENSION(:,:)                             :: biomass_temp         !! Biomass of an individual tree within a circ
444                                                                                    !! class @tex $(m^{2} ind^{-1})$ @endtex
445
446    !! 0.2 Output variables
447               
448    REAL(r_std), DIMENSION(ncirc)                           :: wood_to_ba_eff       !! Effective basal area of an individual tree within a circ
449                                                                                    !! class @tex $(m^{2} ind^{-1})$ @endtex
450
451    !! 0.3 Modified variables
452
453    !! 0.4 Local variables
454    INTEGER(i_std)                                          :: l                    !! Index
455    REAL(r_std)                                             :: woodmass_ind         !! Woodmass of an individual tree
456                                                                                    !! @tex $(gC ind^{-1})$ @endtex
457
458!_ ================================================================================================================================
459 
460 !! 1. Calculate basal area from woodmass
461
462    IF ( is_tree(pft) ) THEN
463       
464       DO l = 1,ncirc
465         
466          ! Woodmass of an individual tree
467          woodmass_ind = biomass_temp(l,isapabove) + biomass_temp(l,isapbelow) + &
468             biomass_temp(l,iheartabove) + biomass_temp(l,iheartbelow)
469
470          ! Basal area of that individual (m2 ind-1) 
471          wood_to_ba_eff(l) = (pi/4*(woodmass_ind/(tree_ff(pft)*pipe_density(pft)*pipe_tune2(pft))) &
472                  **(2./pipe_tune3(pft)))**(pipe_tune3(pft)/(pipe_tune3(pft)+2))
473                 
474       ENDDO
475
476    ELSE
477
478       WRITE(numout,*) 'pft ',pft
479       CALL ipslerr_p (3,'wood_to_ba_eff', &
480            'wood_to_ba_eff is not defined for this PFT.', &
481            'See the output file for more details.','')
482
483    ENDIF
484
485 END FUNCTION wood_to_ba_eff
486
487
488
489!! ================================================================================================================================
490!! FUNCTION     : wood_to_ba
491!!
492!>\BRIEF        Calculate basal area from woody biomass making use of allometric relationships
493!!
494!! DESCRIPTION : Calculate basal area of an individual tree from the woody biomass of that tree making
495!! use of allometric relationships given below. Here basal area is defined in line with its classical
496!! forestry meaning.
497!! (i) woodmass = tree_ff * pipe_density*ba*height
498!! (ii) height = pipe_tune2 * sqrt(4/pi*ba) ** pipe_tune_3 
499!!
500!! RECENT CHANGE(S): None
501!!
502!! RETURN VALUE : basal area (m2 ind-1)
503!!
504!! REFERENCE(S) :
505!!
506!! FLOWCHART    : None
507!! \n
508!_ ================================================================================================================================
509 
510 FUNCTION wood_to_ba(biomass_temp, pft)
511
512 !! 0. Variable and parameter declaration
513
514    !! 0.1 Input variables
515
516    INTEGER(i_std)                                          :: pft               !! PFT number (-)
517    REAL(r_std), DIMENSION(:,:)                             :: biomass_temp      !! Biomass of an individual tree within a circ
518                                                                                 !! class @tex $(m^{2} ind^{-1})$ @endtex
519
520    !! 0.2 Output variables
521               
522    REAL(r_std), DIMENSION(ncirc)                           :: wood_to_ba        !! Basal area of an individual tree within a circ
523                                                                                 !! class @tex $(m^{2} ind^{-1})$ @endtex
524
525    !! 0.3 Modified variables
526
527    !! 0.4 Local variables
528    INTEGER(i_std)                                          :: l                 !! Index
529    REAL(r_std)                                             :: woodmass_ind      !! Woodmass of an individual tree
530                                                                                 !! @tex $(gC ind^{-1})$ @endtex
531
532!_ ================================================================================================================================
533 
534 !! 1. Calculate basal area from woodmass
535
536    IF ( is_tree(pft) ) THEN
537       
538       DO l = 1,ncirc
539         
540          ! Woodmass of an individual tree
541          woodmass_ind = biomass_temp(l,iheartabove) + biomass_temp(l,isapabove)
542
543
544          ! Basal area of that individual (m2 ind-1) 
545          wood_to_ba(l) = (pi/4*(woodmass_ind/(tree_ff(pft)*pipe_density(pft)*pipe_tune2(pft))) &
546                  **(2./pipe_tune3(pft)))**(pipe_tune3(pft)/(pipe_tune3(pft)+2))
547                 
548       ENDDO
549
550    ELSE
551
552       WRITE(numout,*) 'pft ',pft
553       CALL ipslerr_p (3,'wood_to_ba', &
554            'wood_to_ba is not defined for this PFT.', &
555            'See the output file for more details.','')
556
557    ENDIF
558
559 END FUNCTION wood_to_ba
560
561
562!! ================================================================================================================================
563!! FUNCTION     : wood_to_height_eff
564!!
565!>\BRIEF        Calculate the effective tree height from woody biomass making use of allometric relationships
566!!
567!! DESCRIPTION : Calculate the effective height of an individual tree from the woody biomass of that tree making
568!! use of allometric relationships. Effective height makes use of both above and belowground biomass and is
569!! used in the calculation of the allocation according to deleuze and dhote, hydraulic architecture because also
570!! the height of the belowground part should be included.
571!! (i) height(:) = pipe_tune2(j)*(4/pi*ba(:))**(pipe_tune3(j)/2) and
572!! (ii) woodmass_ind = tree_ff*pipe_density*ba*height   
573!!
574!! RECENT CHANGE(S): None
575!!
576!! RETURN VALUE : height (m)
577!!
578!! REFERENCE(S) :
579!!
580!! FLOWCHART    : None
581!! \n
582!_ ================================================================================================================================
583   
584  FUNCTION wood_to_height_eff(biomass_temp, pft)
585
586 !! 0. Variable and parameter declaration
587
588    !! 0.1 Input variables
589
590    INTEGER(i_std)                                          :: pft                !! PFT number (-)
591    REAL(r_std), DIMENSION(:,:)                             :: biomass_temp       !! Biomass of an individual tree within a circ
592                                                                                  !! class @tex $(m^{2} ind^{-1})$ @endtex
593
594    !! 0.2 Output variables
595               
596    REAL(r_std), DIMENSION(ncirc)                           :: wood_to_height_eff !! Effective height of an individual tree within a circ
597                                                                                  !! class (m)
598
599    !! 0.3 Modified variables
600
601    !! 0.4 Local variables
602    INTEGER(i_std)                                          :: l                  !! Index
603    REAL(r_std)                                             :: woodmass_ind       !! Woodmass of an individual tree
604                                                                                  !! @tex $(gC ind{-1})$ @endtex
605                                                                                 
606    CHARACTER(len=256)  :: tmp_var_name                                           !! temporal variable for a text reading (from run.def)   
607    REAL(r_std)         :: tmp_ratio                                              !! temporal variable for adjusting the height   
608!_ ================================================================================================================================
609 
610 !! 1. Calculate height from woodmass
611
612    IF ( is_tree(pft) ) THEN
613       
614       DO l = 1,ncirc
615         
616          ! Woodmass of an individual tree. Both above and belowground biomass are used
617          woodmass_ind = biomass_temp(l,isapabove) + biomass_temp(l,isapbelow) + &
618             biomass_temp(l,iheartabove) + biomass_temp(l,iheartbelow)
619
620          ! Height of that individual
621          wood_to_height_eff(l) = (((woodmass_ind/(tree_ff(pft)*pipe_density(pft))*4/pi)**(pipe_tune3(pft)/2))&
622               *pipe_tune2(pft))**(1/(pipe_tune3(pft)/2+1))
623         
624          !+++++++ TEMP ++++++++++
625          ! This code is only used evaluation of the performance of the multi-layer energy budget.
626          ! To reduce the complexity of the tests we want to impose the LAI and its vertical distribution.
627          ! The solution is not very elegant but it works.
628          ! This part of code use the CIRC_HEIGHT_RATIO to adjust the original heigth in each circonferance class.
629          ! By doing so, the canopy sturcture can be adjust to match the site level observations.
630               IF (ld_fake_height)  THEN
631          !        WRITE(numout,*)'+++ fake wood to height function +++'
632                   !set initial the tmp_ratio reaul to 1.0
633                   tmp_ratio=1.0
634                   WRITE(tmp_var_name, '(A19,I5.5)') "CIRC_HEIGHT_RATIO__", l
635                   CALL getin_p(TRIM(tmp_var_name),  tmp_ratio)
636                   wood_to_height_eff(l) =tmp_ratio*wood_to_height_eff(l)       
637           !       WRITE(numout,*) 'SET CIRC_', l, '_HEIGHT = ', wood_to_height_eff(l)
638               ENDIF                       
639           !++++++++++++++++++++++               
640 
641       ENDDO
642
643    ELSE
644
645       WRITE(numout,*) 'pft ',pft
646       CALL ipslerr_p (3,'wood_to_height_eff', &
647            'wood_to_height_eff is not defined for this PFT.', &
648            'See the output file for more details.','')
649
650    ENDIF
651
652 END FUNCTION wood_to_height_eff
653
654
655!! ================================================================================================================================
656!! FUNCTION     : wood_to_height
657!!
658!>\BRIEF        Calculate tree height from woody biomass making use of allometric relationships
659!!
660!! DESCRIPTION : Calculate height of an individual tree from the woody biomass of that tree making
661!! use of allometric relationships. This is the height used in forestry and for calculating the aerodynamic
662!! interactions
663!! (i) height(:) = pipe_tune2(j)*(4/pi*ba(:))**(pipe_tune3(j)/2) and
664!! (ii) woodmass_ind = tree_ff*pipe_density*ba*height   
665!!
666!! RECENT CHANGE(S): None
667!!
668!! RETURN VALUE : height (m)
669!!
670!! REFERENCE(S) :
671!!
672!! FLOWCHART    : None
673!! \n
674!_ ================================================================================================================================
675   
676  FUNCTION wood_to_height(biomass_temp, pft)
677
678 !! 0. Variable and parameter declaration
679
680    !! 0.1 Input variables
681
682    INTEGER(i_std)                                          :: pft                !! PFT number (-)
683    REAL(r_std), DIMENSION(:,:)                             :: biomass_temp       !! Biomass of an individual tree within a circ
684                                                                                  !! class @tex $(m^{2} ind^{-1})$ @endtex
685
686    !! 0.2 Output variables
687               
688    REAL(r_std), DIMENSION(ncirc)                           :: wood_to_height     !! Height of an individual tree within a circ
689                                                                                  !! class (m)
690
691    !! 0.3 Modified variables
692
693    !! 0.4 Local variables
694    INTEGER(i_std)                                          :: l                  !! Index
695    REAL(r_std)                                             :: woodmass_ind       !! Woodmass of an individual tree
696                                                              !! @tex $(gC ind{-1})$ @endtex
697
698!_ ================================================================================================================================
699 
700 !! 1. Calculate height from woodmass
701
702    IF ( is_tree(pft) ) THEN
703       
704       DO l = 1,ncirc
705         
706          ! Woodmass of an individual tree (only the aboveground component)
707          woodmass_ind = biomass_temp(l,isapabove) + biomass_temp(l,iheartabove)
708
709          ! Height of that individual
710          wood_to_height(l) = (((woodmass_ind/(tree_ff(pft)*pipe_density(pft))*4/pi)**(pipe_tune3(pft)/2))&
711               *pipe_tune2(pft))**(1/(pipe_tune3(pft)/2+1))
712           
713       ENDDO
714
715    ELSE
716
717       WRITE(numout,*) 'pft ',pft
718       CALL ipslerr_p (3,'wood_to_height', &
719            'wood_to_height is not defined for this PFT.', &
720            'See the output file for more details.','')
721
722    ENDIF
723
724 END FUNCTION wood_to_height
725
726
727!! ================================================================================================================================
728!! FUNCTION     : wood_to_qmheight
729!!
730!>\BRIEF        Calculate the quadratic mean height from the biomass
731!!
732!! DESCRIPTION : Calculates the quadratic mean height from the biomass
733!!
734!! RECENT CHANGE(S): None
735!!
736!! RETURN VALUE : ::qm_height (m)
737!!
738!! REFERENCE(S) :
739!!
740!! FLOWCHART    : None
741!! \n
742!_ ================================================================================================================================
743   
744  FUNCTION wood_to_qmheight(biomass_temp, ind, pft)
745
746 !! 0. Variable and parameter declaration
747
748    !! 0.1 Input variables
749
750    INTEGER(i_std)                             :: pft               !! PFT number (-)
751    REAL(r_std), DIMENSION(ncirc,nparts)       :: biomass_temp      !! Biomass of the leaves @tex $(gC m^{-2})$ @endtex
752    REAL(r_std), DIMENSION(ncirc)              :: ind               !! Number of individuals @tex $(m^{-2})$ @endtex
753
754
755    !! 0.2 Output variables
756               
757    REAL(r_std)                                :: wood_to_qmheight  !! quadratic mean height (m)
758
759    !! 0.3 Modified variables
760
761    !! 0.4 Local variables
762    REAL(r_std), DIMENSION(ncirc)              :: circ_class_ba     !! basal area for each circ_class @tex $(m^{2})$ @endtex
763    REAL(r_std)                                :: qm_dia            !! quadratic mean diameter (m)
764
765!_ ================================================================================================================================
766 
767    !! 1. Calculate qm_height from the biomass
768    IF ( is_tree(pft) ) THEN
769
770       ! Basal area at the tree level (m2 tree-1)
771       circ_class_ba(:) = wood_to_ba(biomass_temp(:,:),pft) 
772       
773       IF (SUM(ind(:)) .NE. zero) THEN
774
775          qm_dia = SQRT( 4/pi*SUM(circ_class_ba(:)*ind(:))/SUM(ind(:)) )
776         
777       ELSE
778         
779          qm_dia = zero
780
781       ENDIF
782       
783       wood_to_qmheight = pipe_tune2(pft)*(qm_dia**pipe_tune3(pft))
784             
785
786    ! Grasses and croplands   
787    ELSE
788
789       ! Calculate height as a function of the leaf and structural biomass. Use structural
790       ! biomass to make sure that the grasslands have a roughness length during the winter
791       ! If the biomass increases, vegetation height will increase as well. Divide by
792       ! ind(ipts,j) to obtain the height of the individual. biomass(ileaf) is in gC m-2
793       ! whereas qm is the height of the individual.
794       IF (SUM(ind(:)) .NE. zero) THEN
795
796          wood_to_qmheight = SUM(biomass_temp(:,ileaf) + biomass_temp(:,isapabove)) / &
797               SUM(ind(:)) * sla(pft) * lai_to_height(pft)
798
799       ELSE
800
801          wood_to_qmheight = zero
802
803       ENDIF
804
805    ENDIF ! is_tree(j)
806
807 END FUNCTION wood_to_qmheight
808
809
810
811!! ================================================================================================================================
812!! FUNCTION     : wood_to_dia_eff
813!!
814!>\BRIEF        Calculate effective diameter from woody biomass making use of allometric relationships
815!!
816!! DESCRIPTION : Calculate the effective diameter of an individual tree from the woody biomass of that tree making
817!! use of allometric relationships. Effective diameter accounts for both above and belowground biomass.
818!! (i) woodmass_ind = tree_ff * pipe_density * height * pi/4*dia**2
819!! (ii) height = pipe_tune2 * dia * pipe_tune3
820!!
821!! RECENT CHANGE(S): None
822!!
823!! RETURN VALUE : diameter (m)
824!!
825!! REFERENCE(S) :
826!!
827!! FLOWCHART    : None
828!! \n
829!_ ================================================================================================================================
830   
831  FUNCTION wood_to_dia_eff(biomass_temp, pft)
832
833 !! 0. Variable and parameter declaration
834
835    !! 0.1 Input variables
836
837    INTEGER(i_std)                                          :: pft               !! PFT number (-)
838    REAL(r_std), DIMENSION(:,:)                             :: biomass_temp      !! Biomass of an individual tree within a circ
839                                                                                 !! class @tex $(m^{2} ind^{-1})$ @endtex
840
841    !! 0.2 Output variables
842               
843    REAL(r_std), DIMENSION(ncirc)                           :: wood_to_dia_eff   !! Diameter of an individual tree within a circ
844                                                                                 !! class (m)
845
846    !! 0.3 Modified variables
847
848    !! 0.4 Local variables
849    INTEGER(i_std)                                          :: l                 !! Index
850    REAL(r_std)                                             :: woodmass_ind      !! Woodmass of an individual tree
851                                                                                 !! @tex $(gC ind^{-1})$ @endtex
852
853!_ ================================================================================================================================
854 
855 !! 1. Calculate basal area from woodmass
856
857    IF ( is_tree(pft) ) THEN
858       
859       DO l = 1,ncirc
860         
861          ! Woodmass of an individual tree
862          woodmass_ind = biomass_temp(l,isapabove) + biomass_temp(l,isapbelow) + &
863             biomass_temp(l,iheartabove) + biomass_temp(l,iheartbelow)
864
865          ! Basal area of that individual (m2 ind-1) 
866          wood_to_dia_eff(l) = (4/pi*woodmass_ind/(tree_ff(pft)*pipe_density(pft)*pipe_tune2(pft))) ** &
867                  (1./(2+pipe_tune3(pft)))
868                 
869       ENDDO
870
871    ELSE
872
873       WRITE(numout,*) 'pft ',pft
874       CALL ipslerr_p (3,'wood_to_dia_eff', &
875            'wood_to_dia_eff is not defined for this PFT.', &
876            'See the output file for more details.','')
877
878    ENDIF
879
880 END FUNCTION wood_to_dia_eff
881
882
883
884!! ================================================================================================================================
885!! FUNCTION     : wood_to_dia
886!!
887!>\BRIEF        Calculate diameter from woody biomass making use of allometric relationships
888!!
889!! DESCRIPTION : Calculate diameter of an individual tree from the woody biomass of that tree making
890!! use of allometric relationships. Makes only use of the aboveground biomass and relates to the
891!! typical forestry diameter (but not normalized at 1.3 m)
892!! (i) woodmass_ind = tree_ff * pipe_density * height * pi/4*dia**2
893!! (ii) height = pipe_tune2 * dia * pipe_tune3
894!!
895!! RECENT CHANGE(S): None
896!!
897!! RETURN VALUE : diameter (m)
898!!
899!! REFERENCE(S) :
900!!
901!! FLOWCHART    : None
902!! \n
903!_ ================================================================================================================================
904   
905  FUNCTION wood_to_dia(biomass_temp, pft)
906
907 !! 0. Variable and parameter declaration
908
909    !! 0.1 Input variables
910
911    INTEGER(i_std)                                          :: pft               !! PFT number (-)
912    REAL(r_std), DIMENSION(:,:)                             :: biomass_temp      !! Biomass of an individual tree within a circ
913                                                                                 !! class @tex $(m^{2} ind^{-1})$ @endtex
914
915    !! 0.2 Output variables
916               
917    REAL(r_std), DIMENSION(ncirc)                           :: wood_to_dia       !! Diameter of an individual tree within a circ
918                                                                                 !! class (m)
919
920    !! 0.3 Modified variables
921
922    !! 0.4 Local variables
923    INTEGER(i_std)                                          :: l                 !! Index
924    REAL(r_std)                                             :: woodmass_ind      !! Woodmass of an individual tree
925                                                                                 !! @tex $(gC ind^{-1})$ @endtex
926
927!_ ================================================================================================================================
928 
929 !! 1. Calculate basal area from woodmass
930
931    IF ( is_tree(pft) ) THEN
932       
933       DO l = 1,ncirc
934         
935          ! Woodmass of an individual tree
936          woodmass_ind = biomass_temp(l,isapabove) + biomass_temp(l,iheartabove)
937
938          ! Basal area of that individual (m2 ind-1) 
939          wood_to_dia(l) = (4/pi*woodmass_ind/(tree_ff(pft)*pipe_density(pft)*pipe_tune2(pft))) ** &
940                  (1./(2+pipe_tune3(pft)))
941                 
942       ENDDO
943
944    ELSE
945
946       WRITE(numout,*) 'pft ',pft
947       CALL ipslerr_p (3,'wood_to_dia', &
948            'wood_to_dia is not defined for this PFT.', &
949            'See the output file for more details.','')
950
951    ENDIF
952
953 END FUNCTION wood_to_dia
954
955
956!! ================================================================================================================================
957!! FUNCTION     : wood_to_qmdia
958!!
959!>\BRIEF        Calculate the quadratic mean diameter from the biomass
960!!
961!! DESCRIPTION : Calculates the quadratic mean diameter from the aboveground biomss
962!!
963!! RECENT CHANGE(S): None
964!!
965!! RETURN VALUE : ::qm_dia (m)
966!!
967!! REFERENCE(S) :
968!!
969!! FLOWCHART    : None
970!! \n
971!_ ================================================================================================================================
972   
973  FUNCTION wood_to_qmdia(biomass_temp, ind, pft)
974
975 !! 0. Variable and parameter declaration
976
977    !! 0.1 Input variables
978
979    INTEGER(i_std)                             :: pft               !! PFT number (-)
980    REAL(r_std), DIMENSION(ncirc,nparts)       :: biomass_temp      !! Biomass of the leaves @tex $(gC m^{-2})$ @endtex
981    REAL(r_std), DIMENSION(ncirc)              :: ind               !! Number of individuals @tex $(m^{-2})$ @endtex
982
983    !! 0.2 Output variables
984               
985    REAL(r_std)                                :: wood_to_qmdia     !! quadratic mean diameter (m)
986
987    !! 0.3 Modified variables
988
989    !! 0.4 Local variables
990    REAL(r_std), DIMENSION(ncirc)              :: circ_class_ba     !! basal area for each circ_class @tex $(m^{2})$ @endtex
991
992!_ ================================================================================================================================
993 
994    !! 1. Calculate qm_dia from the biomass
995    IF ( is_tree(pft) ) THEN
996
997       ! Basal area at the tree level (m2 tree-1)
998       circ_class_ba(:) = wood_to_ba(biomass_temp(:,:),pft) 
999       
1000       IF (SUM(ind(:)) .NE. zero) THEN
1001
1002          wood_to_qmdia = SQRT( 4/pi*SUM(circ_class_ba(:)*ind(:))/SUM(ind(:)) )
1003         
1004       ELSE
1005         
1006          wood_to_qmdia = zero
1007
1008       ENDIF
1009
1010
1011    ! Grasses and croplands   
1012    ELSE
1013
1014       wood_to_qmdia = zero   
1015       
1016    ENDIF ! is_tree(pft)
1017
1018 END FUNCTION wood_to_qmdia
1019
1020
1021!! ================================================================================================================================
1022!! FUNCTION     : wood_to_circ
1023!!
1024!>\BRIEF        Calculate circumference from woody biomass making use of allometric relationships
1025!!
1026!! DESCRIPTION : All this does it computer the diameter using a different routine, and then
1027!!               convert that into a circumference. 
1028!!
1029!! RECENT CHANGE(S): None
1030!!
1031!! RETURN VALUE : circumference (m)
1032!!
1033!! REFERENCE(S) :
1034!!
1035!! FLOWCHART    : None
1036!! \n
1037!_ ================================================================================================================================
1038   
1039  FUNCTION wood_to_circ(biomass_temp, pft)
1040
1041 !! 0. Variable and parameter declaration
1042
1043    !! 0.1 Input variables
1044
1045    INTEGER(i_std)                                          :: pft               !! PFT number (-)
1046    REAL(r_std), DIMENSION(:,:)                             :: biomass_temp      !! Biomass of an individual tree within a circ
1047                                                                                 !! class @tex $(m^{2} ind^{-1})$ @endtex
1048
1049    !! 0.2 Output variables
1050               
1051    REAL(r_std), DIMENSION(ncirc)                           :: wood_to_circ      !! Circumference of an individual tree within a circ
1052                                                                                 !! class (m)
1053
1054    !! 0.3 Modified variables
1055
1056    !! 0.4 Local variables
1057
1058!_ ================================================================================================================================
1059 
1060    !! 1. Calculate diameter from woodmass
1061
1062    wood_to_circ(:)=val_exp
1063
1064    wood_to_circ(:)=wood_to_dia(biomass_temp(:,:),pft)
1065
1066    ! convert to a circumference (m)
1067    wood_to_circ(:) =  wood_to_circ(:)*pi
1068
1069 END FUNCTION wood_to_circ
1070
1071
1072
1073!! ================================================================================================================================
1074!! FUNCTION     : wood_to_cn_array
1075!!
1076!>\BRIEF        Calculate crown area from woody biomass making use of allometric relationships
1077!!
1078!! DESCRIPTION : Calculate crown area of an individual tree from the woody biomass of that tree
1079!! making use of allometric relationship that relates crown area (cn) to diameter (dia) as
1080!! pipe_tune1*dia**pipe_tune_exp_coeff where the pipe_tune parameters are pft-specific.
1081!! (i) Basal area is written as a function of wood_mass: woodmass_ind = tree_ff*pipe_density*ba*height
1082!! (ii) height = pipe_tune2*sqrt(4/pi*ba)**pipe_tune3
1083!! (iii) cn = pipe_tune1*sqrt(4/pi*ba)**pipe_tune_exp_coeff   
1084!!
1085!! RECENT CHANGE(S): None
1086!!
1087!! RETURN VALUE : crown area (m2 ind-1)
1088!!
1089!! REFERENCE(S) :
1090!!
1091!! FLOWCHART    : None
1092!! \n
1093!_ ================================================================================================================================
1094
1095  FUNCTION wood_to_cn_array(biomass_temp, npts, pft)
1096
1097 !! 0. Variable and parameter declaration
1098
1099    !! 0.1 Input variables
1100
1101    INTEGER(i_std)                           :: pft              !! PFT number (-)
1102    INTEGER(i_std)                           :: npts             !! Pixel(s), this variable defines the dimensions of ba
1103    REAL(r_std), DIMENSION(:,:,:)            :: biomass_temp     !! Biomass of an individual tree within a circ
1104                                                                 !! class @tex $(m^{2} ind^{-1})$ @endtex
1105   
1106
1107    !! 0.2 Output variables
1108               
1109    REAL(r_std), DIMENSION(npts,ncirc)       :: wood_to_cn_array !! Crown area of an individual tree @tex $(m^{2} ind{-1})$ @endtex
1110
1111    !! 0.3 Modified variables
1112
1113    !! 0.4 Local variables
1114 
1115    INTEGER(i_std)                           :: l                !! index
1116    REAL(r_std), DIMENSION(npts)             :: woodmass_ind     !! Woodmass of an individual tree @tex $(gC ind{-1})$ @endtex
1117
1118!_ ================================================================================================================================
1119 
1120 !! 1. Calculate crown area from basal area
1121
1122    IF ( is_tree(pft) ) THEN
1123       
1124       DO l = 1,ncirc
1125
1126          ! Woodmass of an individual tree
1127          woodmass_ind(:) = biomass_temp(:,l,isapabove) + biomass_temp(:,l,iheartabove)
1128
1129          ! Crown area of that individual
1130          ! cn = pipe_tune1*sqrt(4/pi*ba)**pipe_tune_exp_coeff
1131          wood_to_cn_array(:,l) = pipe_tune1(pft) * SQRT( 4/pi*(pi/4*(woodmass_ind(:)/&
1132               (tree_ff(pft)*pipe_density(pft)*pipe_tune2(pft)))**(2./pipe_tune3(pft)))**&
1133               (pipe_tune3(pft)/(pipe_tune3(pft)+2)) ) ** pipe_tune_exp_coeff(pft)
1134
1135       ENDDO
1136
1137    ELSE
1138
1139       WRITE(numout,*) 'pft ',pft
1140       CALL ipslerr_p (3,'wood_to_cn_array', &
1141            'wood_to_cn_array is not defined for this PFT.', &
1142            'See the output file for more details.','')
1143
1144    ENDIF
1145
1146  END FUNCTION wood_to_cn_array
1147
1148
1149!! ================================================================================================================================
1150!! FUNCTION     : wood_to_cn
1151!!
1152!>\BRIEF        Calculate crown area from woody biomass making use of allometric relationships
1153!!
1154!! DESCRIPTION : Calculate crown area of an individual tree from the woody biomass of that tree
1155!! making use of allometric relationship that relates crown area (cn) to diameter (dia) as
1156!! pipe_tune1*dia**pipe_tune_exp_coeff where the pipe_tune parameters are pft-specific.
1157!! (i) Basal area is written as a function of wood_mass: woodmass_ind = tree_ff*pipe_density*ba*height
1158!! (ii) height = pipe_tune2*sqrt(4/pi*ba)**pipe_tune3
1159!! (iii) cn = pipe_tune1*sqrt(4/pi*ba)**pipe_tune_exp_coeff   
1160!!
1161!! RECENT CHANGE(S): None
1162!!
1163!! RETURN VALUE : crown area (m2 ind-1)
1164!!
1165!! REFERENCE(S) :
1166!!
1167!! FLOWCHART    : None
1168!! \n
1169!_ ================================================================================================================================
1170   
1171  FUNCTION wood_to_cn(biomass_temp, pft)
1172
1173 !! 0. Variable and parameter declaration
1174
1175    !! 0.1 Input variables
1176
1177    INTEGER(i_std)                                          :: pft               !! PFT number (-)
1178    REAL(r_std), DIMENSION(:,:)                             :: biomass_temp      !! biomass of an individual tree within a circ
1179                                                                                 !! class @tex $(m^{2} ind^{-1})$ @endtex
1180
1181    !! 0.2 Output variables
1182               
1183    REAL(r_std), DIMENSION(ncirc)                           :: wood_to_cn        !! Crown area of an individual tree within a circ
1184                                                                                 !! class @tex $(m^{2} ind-1)$ @endtex
1185
1186    !! 0.3 Modified variables
1187
1188    !! 0.4 Local variables
1189    INTEGER(i_std)                                          :: l                 !! Index
1190    REAL(r_std)                                             :: woodmass_ind      !! Woodmass of an individual tree
1191                                                                                 !! @tex $(gC ind{-1})$ @endtex
1192
1193!_ ================================================================================================================================
1194 
1195 !! 1. Calculate crown area from woodmass
1196
1197    IF ( is_tree(pft) ) THEN
1198       
1199       DO l = 1,ncirc
1200         
1201          ! Woodmass of an individual tree
1202          woodmass_ind = biomass_temp(l,isapabove) + biomass_temp(l,iheartabove)
1203 
1204          ! Crown area of that individual
1205          ! cn = pipe_tune1*sqrt(4/pi*ba)**pipe_tune_exp_coeff
1206          wood_to_cn(l) = pipe_tune1(pft) * SQRT( 4/pi*(pi/4*(woodmass_ind/&
1207               (tree_ff(pft)*pipe_density(pft)*pipe_tune2(pft)))**(2./pipe_tune3(pft)))**&
1208               (pipe_tune3(pft)&
1209               /(pipe_tune3(pft)+2)) ) ** pipe_tune_exp_coeff(pft)
1210               
1211       ENDDO
1212
1213    ELSE
1214
1215       WRITE(numout,*) 'pft ',pft
1216       CALL ipslerr_p (3,'wood_to_cn', &
1217            'wood_to_cn is not defined for this PFT.', &
1218            'See the output file for more details.','')
1219
1220    ENDIF
1221
1222 END FUNCTION wood_to_cn
1223
1224!! ================================================================================================================================
1225!! FUNCTION     : wood_to_cn_eff
1226!!
1227!>\BRIEF        Calculate crown area from woody biomass making use of allometric relationships
1228!!
1229!! DESCRIPTION : Calculate crown area of an individual tree from the woody biomass of that tree
1230!! making use of allometric relationship that relates crown area (cn) to diameter (dia) as
1231!! pipe_tune1*dia**pipe_tune_exp_coeff where the pipe_tune parameters are pft-specific.
1232!! (i) Basal area is written as a function of wood_mass: woodmass_ind = tree_ff*pipe_density*ba*height
1233!! (ii) height = pipe_tune2*sqrt(4/pi*ba)**pipe_tune3
1234!! (iii) cn = pipe_tune1*sqrt(4/pi*ba)**pipe_tune_exp_coeff   
1235!!
1236!! NOTE: This is different from wood_to_cn because this uses both the aboveground and
1237!!       belowground wood mass.
1238!!
1239!! RECENT CHANGE(S): None
1240!!
1241!! RETURN VALUE : crown area (m2 ind-1)
1242!!
1243!! REFERENCE(S) :
1244!!
1245!! FLOWCHART    : None
1246!! \n
1247!_ ================================================================================================================================
1248   
1249  FUNCTION wood_to_cn_eff(biomass_temp, pft)
1250
1251 !! 0. Variable and parameter declaration
1252
1253    !! 0.1 Input variables
1254
1255    INTEGER(i_std)                                          :: pft               !! PFT number (-)
1256    REAL(r_std), DIMENSION(:,:)                             :: biomass_temp      !! biomass of an individual tree within a circ
1257                                                                                 !! class @tex $(m^{2} ind^{-1})$ @endtex
1258
1259    !! 0.2 Output variables
1260               
1261    REAL(r_std), DIMENSION(ncirc)                           :: wood_to_cn_eff    !! Crown area of an individual tree within a circ
1262                                                                                 !! class @tex $(m^{2} ind-1)$ @endtex
1263
1264    !! 0.3 Modified variables
1265
1266    !! 0.4 Local variables
1267    INTEGER(i_std)                                          :: l                 !! Index
1268    REAL(r_std)                                             :: woodmass_ind      !! Woodmass of an individual tree
1269                                                                                 !! @tex $(gC ind{-1})$ @endtex
1270
1271!_ ================================================================================================================================
1272 
1273 !! 1. Calculate crown area from woodmass
1274
1275    IF ( is_tree(pft) ) THEN
1276       
1277       DO l = 1,ncirc
1278         
1279          ! Woodmass of an individual tree
1280          woodmass_ind = biomass_temp(l,isapabove) + biomass_temp(l,iheartabove) + &
1281               biomass_temp(l,isapbelow) + biomass_temp(l,iheartbelow)
1282 
1283          ! Crown area of that individual
1284          ! cn = pipe_tune1*sqrt(4/pi*ba)**pipe_tune_exp_coeff
1285          wood_to_cn_eff(l) = pipe_tune1(pft) * SQRT( 4/pi*(pi/4*(woodmass_ind/&
1286               (tree_ff(pft)*pipe_density(pft)*pipe_tune2(pft)))**(2./pipe_tune3(pft)))**&
1287               (pipe_tune3(pft)&
1288               /(pipe_tune3(pft)+2)) ) ** pipe_tune_exp_coeff(pft)
1289               
1290       ENDDO
1291
1292    ELSE
1293
1294       WRITE(numout,*) 'pft ',pft
1295       CALL ipslerr_p (3,'wood_to_cn_eff', &
1296            'wood_to_cn_eff is not defined for this PFT.', &
1297            'See the output file for more details.','')
1298
1299    ENDIF
1300
1301  END FUNCTION wood_to_cn_eff
1302
1303
1304!! ================================================================================================================================
1305!! FUNCTION     : wood_to_cv
1306!!
1307!>\BRIEF        Calculate crown volume from woody biomass making use of allometric relationships
1308!!
1309!! DESCRIPTION : Calculate basal areadiameter of an individual tree from the woody biomass of that tree making
1310!! use of allometric relationships:
1311!! (i) Basal area is written as a function of wood_mass: woodmass_ind = tree_ff*pipe_density*ba*height
1312!! (ii) height = pipe_tune2 * sqrt( 4/pi*ba ) ** pipe_tune3
1313!! (iii) cn = pipe_tune1 * sqrt( 4/pi*ba ) ** pipe_tune_exp_coeff 
1314!! (iv) cv = 4.0 / 3.0 * pi * ( SQRT( cn/pi ) ) ** 3
1315!!
1316!! RECENT CHANGE(S): None
1317!!
1318!! RETURN VALUE : diameter (m3 ind-1)
1319!!
1320!! REFERENCE(S) :
1321!!
1322!! FLOWCHART    : None
1323!! \n
1324!_ ================================================================================================================================
1325   
1326  FUNCTION wood_to_cv(biomass_temp, pft)
1327
1328 !! 0. Variable and parameter declaration
1329
1330    !! 0.1 Input variables
1331
1332    INTEGER(i_std)                                          :: pft               !! PFT number (-)
1333    REAL(r_std), DIMENSION(:,:)                             :: biomass_temp      !! Biomass of an individual tree within a circ
1334                                                                                 !! class @tex $(m^{2} ind^{-1})$ @endtex
1335
1336    !! 0.2 Output variables
1337               
1338    REAL(r_std), DIMENSION(ncirc)                           :: wood_to_cv        !! Crown volume of an individual tree
1339                                                                                 !! @tex $(m^{2} ind{-1})$ @endtex
1340
1341    !! 0.3 Modified variables
1342
1343    !! 0.4 Local variables
1344    INTEGER(i_std)                                          :: l                 !! Index
1345    REAL(r_std)                                             :: woodmass_ind      !! Woodmass of an individual tree
1346                                                                                 !! @tex $(gC ind^{-1})$ @endtex
1347    REAL(r_std), DIMENSION(ncirc)                           :: dia               !! Diameter of an individual tree (m)
1348
1349!_ ================================================================================================================================
1350 
1351 !! 1. Calculate crown volume from woodmass
1352
1353    IF ( is_tree(pft) ) THEN
1354
1355       ! Crown diameter of the individual tree (m2 ind-1) 
1356       dia(:) = (4./pi*wood_to_cn(biomass_temp(:,:),pft))**(1./2.)
1357       wood_to_cv(:) = pi/6.*dia(:)**3. 
1358       ! WRITE(numout,*) 'wood_to_cv, ', wood_to_cv(:)
1359
1360    ELSE
1361
1362       WRITE(numout,*) 'pft ',pft
1363       CALL ipslerr_p (3,'wood_to_cv', &
1364            'wood_to_cv is not defined for this PFT.', &
1365            'See the output file for more details.','')
1366
1367    ENDIF
1368
1369 END FUNCTION wood_to_cv 
1370
1371   !! ================================================================================================================================
1372!! FUNCTION     : wood_to_cv_eff
1373!!
1374!>\BRIEF        Calculate crown volume from woody biomass making use of allometric relationships
1375!!
1376!! DESCRIPTION : Calculate basal areadiameter of an individual tree from the woody biomass of that tree making
1377!! use of allometric relationships:
1378!! (i) Basal area is written as a function of wood_mass: woodmass_ind = tree_ff*pipe_density*ba*height
1379!! (ii) height = pipe_tune2 * sqrt( 4/pi*ba ) ** pipe_tune3
1380!! (iii) cn = pipe_tune1 * sqrt( 4/pi*ba ) ** pipe_tune_exp_coeff 
1381!! (iv) cv = 4.0 / 3.0 * pi * ( SQRT( cn/pi ) ) ** 3
1382!!
1383!! NOTE: the difference between this and wood_to_cv is that this include aboveground and belowground
1384!!       biomass in the calculation
1385!!
1386!! RECENT CHANGE(S): None
1387!!
1388!! RETURN VALUE : diameter (m3 ind-1)
1389!!
1390!! REFERENCE(S) :
1391!!
1392!! FLOWCHART    : None
1393!! \n
1394!_ ================================================================================================================================
1395   
1396  FUNCTION wood_to_cv_eff(biomass_temp, pft)
1397
1398 !! 0. Variable and parameter declaration
1399
1400    !! 0.1 Input variables
1401
1402    INTEGER(i_std)                                          :: pft               !! PFT number (-)
1403    REAL(r_std), DIMENSION(:,:)                             :: biomass_temp      !! Biomass of an individual tree within a circ
1404                                                                                 !! class @tex $(m^{2} ind^{-1})$ @endtex
1405
1406    !! 0.2 Output variables
1407               
1408    REAL(r_std), DIMENSION(ncirc)                           :: wood_to_cv_eff    !! Crown volume of an individual tree
1409                                                                                 !! @tex $(m^{2} ind{-1})$ @endtex
1410
1411    !! 0.3 Modified variables
1412
1413    !! 0.4 Local variables
1414    INTEGER(i_std)                                          :: l                 !! Index
1415    REAL(r_std)                                             :: woodmass_ind      !! Woodmass of an individual tree
1416                                                                                 !! @tex $(gC ind^{-1})$ @endtex
1417    REAL(r_std), DIMENSION(ncirc)                           :: dia               !! Diameter of an individual tree (m)
1418
1419!_ ================================================================================================================================
1420 
1421 !! 1. Calculate crown volume from woodmass
1422
1423    IF ( is_tree(pft) ) THEN
1424
1425       ! Crown diameter of the individual tree (m2 ind-1) 
1426       dia(:) = (4./pi*wood_to_cn_eff(biomass_temp(:,:),pft))**(1./2.)
1427       wood_to_cv_eff(:) = pi/6.*dia(:)**3. 
1428       ! WRITE(numout,*) 'wood_to_cv_eff, ', wood_to_cv_eff(:)
1429
1430    ELSE
1431
1432       WRITE(numout,*) 'pft ',pft
1433       CALL ipslerr_p (3,'wood_to_cv_eff', &
1434            'wood_to_cv_eff is not defined for this PFT.', &
1435            'See the output file for more details.','')
1436
1437    ENDIF
1438
1439  END FUNCTION wood_to_cv_eff
1440
1441
1442
1443!! ================================================================================================================================
1444!! FUNCTION     : wood_to_volume
1445!!
1446!>\BRIEF        This allometric function computes volume as a function of
1447!! biomass at stand scale. Volume \f$(m^3 m^{-2}) = f(biomass (gC m^{-2}))\f$
1448!!
1449!! DESCRIPTION : None
1450!!
1451!! RECENT CHANGE(S): None
1452!!
1453!! RETURN VALUE : biomass_to_volume
1454!!
1455!! REFERENCE(S) : See above, module description.
1456!!
1457!! FLOWCHART    : None
1458!! \n
1459!_ ================================================================================================================================
1460
1461  FUNCTION wood_to_volume(biomass,pft,branch_ratio,inc_branches)
1462
1463 !! 0. Variable and parameter declaration
1464
1465    !! 0.1 Input variables
1466
1467    REAL(r_std), DIMENSION(:)   :: biomass              !! Stand biomass @tex $(gC m^{-2})$ @endtex 
1468    REAL(r_std)                 :: branch_ratio         !! Branch ratio of sap and heartwood biomass
1469                                                        !! unitless
1470    INTEGER(i_std)              :: pft                  !! Plant functional type (unitless)
1471    INTEGER(i_std)              :: inc_branches         !! Include the branches in the volume calculation?
1472                                                        !! 0: exclude the branches from the volume calculation
1473                                                        !! (thus correct the biomass for the branch ratio)
1474                                                        !! 1: include the branches in the volume calculation
1475                                                        !! (thus use all aboveground biomass)
1476                                   
1477   
1478
1479    !! 0.2 Output variables
1480
1481    REAL(r_std)                 :: wood_to_volume       !! The volume of wood per square meter
1482                                                        !!  @tex $(m^3 m^{-2})$ @endtex
1483
1484    !! 0.3 Modified variables
1485
1486    !! 0.4 Local variables
1487   
1488    REAL(r_std)                 :: woody_biomass        !! Woody biomass at the stand level
1489                                                        !! @tex $(gC m^{-2})$ @endtex
1490
1491!_ ================================================================================================================================
1492
1493 !! 1. Volume to biomass
1494
1495    ! Woody biomass used in the calculation
1496    IF (inc_branches .EQ. 0) THEN
1497
1498       woody_biomass=(biomass(isapabove)+biomass(iheartabove))*(un - branch_ratio)
1499
1500    ELSEIF (inc_branches .EQ. 1) THEN
1501
1502       woody_biomass=(biomass(isapabove)+biomass(iheartabove))
1503
1504    ELSE
1505
1506    ENDIF
1507
1508    ! Wood volume expressed in m**3 / m**2
1509    wood_to_volume = woody_biomass/(pipe_density(pft))
1510
1511  END FUNCTION wood_to_volume
1512
1513
1514
1515!! ================================================================================================================================
1516!! FUNCTION     : biomass_to_lai
1517!!
1518!>\BRIEF        Calculate the LAI based on the leaf biomass
1519!!
1520!! DESCRIPTION : Calculates the LAI of a PFT/grid square based on the leaf biomass
1521!!
1522!! RECENT CHANGE(S): None
1523!!
1524!! RETURN VALUE : ::LAI [m**2 m**{-2}]
1525!!
1526!! REFERENCE(S) :
1527!!
1528!! FLOWCHART    : None
1529!! \n
1530!_ ================================================================================================================================
1531   
1532  FUNCTION biomass_to_lai(leaf_biomass, pft)
1533
1534 !! 0. Variable and parameter declaration
1535
1536    !! 0.1 Input variables
1537
1538    INTEGER(i_std)                                          :: pft               !! PFT number (-)
1539    REAL(r_std)                                             :: leaf_biomass      !! Biomass of the leaves
1540                                                                                 !! @tex $(gC m^{-2})$ @endtex
1541
1542
1543    !! 0.2 Output variables
1544               
1545    REAL(r_std)                                             :: biomass_to_lai    !! Leaf area index
1546                                                                                 !! @tex $(m^{2} m^{-2})$ @endtex
1547
1548    !! 0.3 Modified variables
1549
1550    !! 0.4 Local variables
1551    REAL(r_std)                                             :: impose_lai         !! LAI read from run.def
1552!_ ================================================================================================================================
1553 
1554    !! 1. Calculate the LAI from the leaf biomass
1555
1556    biomass_to_lai = leaf_biomass * sla(pft)
1557     
1558!!$    !+++++++++ TEMP ++++++++++
1559!!$    ! This is a perfect place to hack the code to make it run with
1560!!$    ! constant lai
1561!!$    WRITE(numout,*) 'WARNING ERROR: Using fake lai values for testing!'
1562!!$    biomass_to_lai=3.79052
1563!!$    !+++++++++++++++++++++++++
1564
1565            !+++++++ TEMP ++++++++++
1566            ! This code is only used evaluation of the performance of the multi-layer energy budget.
1567            ! To reduce the complexity of the tests we want to impose the LAI and its vertical distribution.
1568            ! The solution is not very elegant but it works.
1569              IF (ld_fake_height) THEN
1570            ! In order to imposed lai, we read the TOTAL_LAI from run.def
1571              CALL getin_p('TOTAL_LAI', impose_lai)
1572            ! This part of code reset the sla vale to match which alow modeled LAI equal to TOTAL LAI.
1573            ! Althought this is ugly way to match the modeled LAI and impose LAI.
1574            ! You probably need to go to your ORCHIDEE out put file to find out the suitable SLA value
1575            ! and reset it agin in the run.def.
1576            ! So, we impose LAI & structure for a quick testing the performance of multilayer energy budget
1577            ! without changing the leaf_biomass.   
1578                 IF ( leaf_biomass .GT. 0.0) THEN
1579                    sla(pft)=impose_lai/leaf_biomass
1580                      WRITE(numout,'(A,F20.8)') 'USE A FAKE SLA BASED ON imposed LAI/LEAFMASS=', sla(pft)
1581                 ENDIF
1582                    biomass_to_lai=leaf_biomass*sla(pft) 
1583              ENDIF 
1584           !++++++++++++++++++++++++
1585
1586     END FUNCTION biomass_to_lai
1587
1588
1589!! ================================================================================================================================
1590!! FUNCTION     : Nmax
1591!!
1592!>\BRIEF        This function determines the maximum number of trees per hectare for
1593!! a given quadratic mean diameter (Dg). It applies the self-thinning principle
1594!! of Reineke (1933), with Dg instead of mean diameter (Dhote 1999).
1595!! Parameterization: Dhote (1999) for broad-leaved and Vacchiano (2008)
1596!! for needle-leaved.
1597!!
1598!! DESCRIPTION : None
1599!!
1600!! RECENT CHANGE(S): None
1601!!
1602!! RETURN VALUE : Nmax
1603!!
1604!! REFERENCE(S)   : See above, module description.
1605!!
1606!! FLOWCHART   : None
1607!! \n
1608!_ ================================================================================================================================
1609   
1610  FUNCTION Nmax(Dg,no_pft)
1611
1612 !! 0. Variable and parameter declaration
1613
1614    !! 0.1 Input variables
1615
1616    REAL(r_std)      :: Dg              !! Quadratic mean diameter (cm)
1617    INTEGER(i_std)   :: no_pft          !! Plant functional type (unitless)
1618
1619    !! 0.2 Output variables
1620
1621    REAL(r_std)      :: Nmax            !! Maximum number of trees according to the self-thinning model
1622
1623    !! 0.3 Modified variables
1624
1625    !! 0.4 Local variables
1626
1627!_ ================================================================================================================================
1628 !! 1. maximum number of trees per hectare for a given quadratic mean diameter
1629
1630    IF (is_tree(no_pft)) THEN
1631
1632       ! thinning curve of the MTC
1633       Nmax = (Dg/alpha_self_thinning(no_pft))**(un/beta_self_thinning(no_pft))
1634
1635       ! Truncate the range to avoid huge numbers due the exponental model used to describe self-thinning
1636       Nmax = MIN( Nmax, nmaxtrees(no_pft) )
1637       Nmax = MAX( Nmax, dens_target(no_pft) )
1638
1639    ELSE
1640
1641       WRITE(numout,*) 'Self thinning is not defined for PFT, ', no_pft
1642       CALL ipslerr_p (3,'nmax', &
1643            'Self thinning is not defined for this PFT.', &
1644            'See the output file for more details.','')
1645
1646    ENDIF
1647
1648  END FUNCTION Nmax
1649
1650
1651!! ================================================================================================================================
1652!! FUNCTION     : Nmaxyield
1653!!
1654!>\BRIEF        This function determines the maximum number of trees per hectare for
1655!! a given quadratic mean diameter (Dg). It applies the 75 percentile  of the European
1656!! yield tables
1657!!
1658!! DESCRIPTION : None
1659!!
1660!! RECENT CHANGE(S): None
1661!!
1662!! RETURN VALUE : Nmaxyield
1663!!
1664!! REFERENCE(S)   : See above, module description.
1665!!
1666!! FLOWCHART   : None
1667!! \n
1668!_ ================================================================================================================================
1669   
1670  FUNCTION Nmaxyield(Dg,no_pft)
1671
1672 !! 0. Variable and parameter declaration
1673
1674    !! 0.1 Input variables
1675
1676    REAL(r_std)      :: Dg              !! Quadratic mean diameter (cm)
1677    INTEGER(i_std)   :: no_pft          !! Plant functional type (unitless)
1678
1679    !! 0.2 Output variables
1680
1681    REAL(r_std)      :: Nmaxyield       !! Maximum number of trees according to the self-thinning model
1682
1683    !! 0.3 Modified variables
1684
1685    !! 0.4 Local variables
1686
1687!_ ================================================================================================================================
1688 !! 1. maximum number of trees per hectare for a given quadratic mean diameter
1689
1690    IF (is_tree(no_pft)) THEN
1691
1692       ! thinning curve of the MTC
1693       Nmaxyield = (Dg/alpha_rdi_upper(no_pft))**(un/beta_rdi_upper(no_pft))
1694
1695       ! Truncate the range to avoid huge numbers due the exponental model used to describe self-thinning
1696       Nmaxyield = MIN( Nmaxyield, nmaxtrees(no_pft) )
1697       Nmaxyield = MAX( Nmaxyield, dens_target(no_pft) )
1698
1699    ELSE
1700
1701       WRITE(numout,*) 'Self thinning is not defined for PFT, ', no_pft
1702       CALL ipslerr_p (3,'Nmaxyield', &
1703            'Self thinning is not defined for this PFT.', &
1704            'See the output file for more details.','')
1705
1706    ENDIF
1707
1708  END FUNCTION Nmaxyield
1709
1710
1711!! ================================================================================================================================
1712!! FUNCTION     : Nminyield
1713!!
1714!>\BRIEF        This function determines the minimum number of trees per hectare for
1715!! a given quadratic mean diameter (Dg). It applies the 25 percentile of the European
1716!! yield tables
1717!!
1718!! DESCRIPTION : None
1719!!
1720!! RECENT CHANGE(S): None
1721!!
1722!! RETURN VALUE : Nminyield
1723!!
1724!! REFERENCE(S)   : See above, module description.
1725!!
1726!! FLOWCHART   : None
1727!! \n
1728!_ ================================================================================================================================
1729   
1730  FUNCTION Nminyield(Dg,no_pft)
1731
1732 !! 0. Variable and parameter declaration
1733
1734    !! 0.1 Input variables
1735
1736    REAL(r_std)      :: Dg              !! Quadratic mean diameter (cm)
1737    INTEGER(i_std)   :: no_pft          !! Plant functional type (unitless)
1738
1739    !! 0.2 Output variables
1740
1741    REAL(r_std)      :: Nminyield       !! Minimum number of trees according to the self-thinning model
1742
1743    !! 0.3 Modified variables
1744
1745    !! 0.4 Local variables
1746
1747!_ ================================================================================================================================
1748 !! 1. minimum number of trees per hectare for a given quadratic mean diameter
1749
1750    IF (is_tree(no_pft)) THEN
1751
1752       ! thinning curve of the MTC
1753       Nminyield = (Dg/alpha_rdi_lower(no_pft))**(un/beta_rdi_lower(no_pft))
1754
1755       ! Truncate the range to avoid huge numbers due the exponental model used to describe self-thinning
1756       
1757       Nminyield = MIN( Nminyield, nmaxtrees(no_pft) )
1758       Nminyield = MAX( Nminyield, dens_target(no_pft) )
1759
1760    ELSE
1761
1762       WRITE(numout,*) 'Self thinning is not defined for PFT, ', no_pft
1763       CALL ipslerr_p (3,'Nminyield', &
1764            'Self thinning is not defined for this PFT.', &
1765            'See the output file for more details.','')
1766
1767
1768    ENDIF
1769
1770  END FUNCTION Nminyield
1771
1772
1773!! ================================================================================================================================
1774!! SUBROUTINE  : distribute_mortality_biomass
1775!!
1776!>\BRIEF       Distributes biomass that is going to be killed by natural
1777!!             causes (not self thinning) over circ classes.
1778!!
1779!! DESCRIPTION : Mortality is going to kill a certain amount of biomass
1780!!               in forests every day.  Since we have circumference classes
1781!!               now for our forests, we need to determine which classes
1782!!               of trees will suffer from this environmental mortality.
1783!!               Right now we are taking an exponential distribution.
1784!!               Notice that this is NOT the same as redistributing biomass
1785!!               after one of the circ classes becomes empty.
1786!!
1787!! RECENT CHANGE(S): None
1788!!
1789!! MAIN OUTPUT VARIABLE(S): ::circ_class_kill
1790!!
1791!! REFERENCE(S) : None
1792!!
1793!! FLOWCHART    : None
1794!! \n
1795!_ ================================================================================================================================
1796  SUBROUTINE distribute_mortality_biomass ( bm_difference, ddf_temp, circ_class_n_temp, &
1797       circ_class_biomass_temp, circ_class_kill_temp )
1798
1799 !! 0. Variable and parameter description
1800
1801    !! 0.1 Input variables
1802    REAL(r_std),INTENT(in)                                    :: bm_difference !! the biomass to distribute
1803    REAL(r_std),INTENT(in)                                    :: ddf_temp !! the death_distribution_factor for this pft
1804    REAL(r_std),DIMENSION(:),INTENT(in)                       :: circ_class_n_temp !! circ_class_n for this point/PFT
1805    REAL(r_std),DIMENSION(:,:),INTENT(in)                     :: circ_class_biomass_temp !! circ_class_biomass for
1806                                                                          !! this point/PFT
1807    !! 0.2 Output variables
1808
1809    !! 0.3 Modified variables
1810    REAL(r_std),DIMENSION(:),INTENT(inout)                    :: circ_class_kill_temp !! circ_class_kill for
1811                                                                                      !! this point/PFT/pool
1812    !! 0.4 Local variables
1813    REAL(r_std), DIMENSION(ncirc)                             :: biomass_desired    !! The biomass that dies naturally
1814    REAL(r_std), DIMENSION(ncirc)                             :: death_distribution !! The fraction of biomass taken from
1815                                                                                    !! each circ class for mortality.
1816    LOGICAL                                                   :: ldone              !! Flag to exit a loop
1817    REAL(r_std)                                               :: scale_factor       !!
1818    REAL(r_std)                                               :: sum_total          !!
1819    REAL(r_std)                                               :: leftover_bm        !! excess biomass that we need to kill
1820    REAL(r_std)                                               :: living_biomass      !! summed biomass
1821    REAL(r_std)                                               :: living_trees       !! summed biomass
1822    INTEGER                                                   :: icir
1823
1824!_ ================================================================================================================================
1825
1826    IF(ncirc == 1)THEN
1827
1828       ! This is the easy case.  All of our biomass will be taken from the only
1829       ! circumference class that we have.
1830       circ_class_kill_temp(1)=circ_class_kill_temp(1)+&
1831            bm_difference/SUM(circ_class_biomass_temp(1,:))
1832
1833       RETURN
1834    ENDIF
1835
1836    ! Here we assume an exponential distribution, arranged so that
1837    ! ddf_temp times more biomass is taken from the largest circ class
1838    ! compared to the smallest.
1839    biomass_desired(:)=zero
1840    death_distribution(:)=un
1841    scale_factor=ddf_temp**(un/REAL(ncirc-1))
1842    DO icir=2,ncirc
1843       death_distribution(icir)=death_distribution(icir-1)*scale_factor
1844    ENDDO
1845    ! Normalize it
1846    ! added by yitong yao Feb 24 2020 15:33 
1847    ! only kill circ meet the condition
1848    ! added by yitong yao Feb 24 2020 15:40
1849
1850    sum_total=SUM(death_distribution(:))
1851    death_distribution(:)=death_distribution(:)/sum_total
1852
1853    ! Ideally, how much is killed from each class?  Be careful to include
1854    ! what was killed in self-thinning here!  If we don't, we may try to kill
1855    ! more biomass than is available in the loop below.
1856    DO icir=1,ncirc
1857       biomass_desired(icir)=death_distribution(icir)*bm_difference !+&
1858    !        circ_class_kill_temp(icir)*SUM(circ_class_biomass_temp(icir,:))
1859    ENDDO
1860
1861    ! Right now, we know how much biomass we want to kill in each
1862    ! class (biomass_desired).  What we will do now is to loop through
1863    ! all the circ_classes and see if we have this much biomass
1864    ! in each class still alive.  The total amount of vegetation still
1865    ! alive is circ_class_n_temp(icir)-circ_class_kill_temp(icir).  If
1866    ! this number of individuals cannot give us the total biomass that
1867    ! we need, we keep track of a residual quantity, leftover_bm, which
1868    ! we try to take from other circ_classes.  It's possible that we
1869    ! will have to loop several times, which makes things more complicated.
1870
1871    ldone=.FALSE.
1872    leftover_bm=zero
1873    DO
1874       DO icir=ncirc,1,-1
1875
1876          living_trees=circ_class_n_temp(icir)-circ_class_kill_temp(icir)
1877          living_biomass=&
1878               SUM(circ_class_biomass_temp(icir,:))*living_trees
1879          biomass_desired(icir)=biomass_desired(icir)+leftover_bm
1880
1881          IF(living_biomass .LE. biomass_desired(icir))THEN
1882
1883             ! We can't get everything from this class, so we kill whatever is left.
1884             leftover_bm=biomass_desired(icir)-living_biomass
1885             biomass_desired(icir)=zero
1886             circ_class_kill_temp(icir)=circ_class_kill_temp(icir)+&
1887                  living_trees
1888
1889          ELSE
1890
1891             ! We have enough in this class.
1892             circ_class_kill_temp(icir)=circ_class_kill_temp(icir)+&
1893                  biomass_desired(icir)/SUM(circ_class_biomass_temp(icir,:))
1894             biomass_desired(icir)=zero
1895             leftover_bm=zero
1896
1897          ENDIF ! living_biomass .LE. biomass_desired(icir)
1898
1899       ENDDO ! loop over circ classes
1900
1901       IF(leftover_bm .LE. min_stomate) EXIT
1902
1903       ! it's possible that we don't have enough biomass left to kill what needs to be
1904       ! killed, so everything just dies.  I cannot think of a case where this
1905       ! would happen, though, since mortality should always be a percentage of the
1906       ! total biomass.
1907       ldone=.TRUE.
1908       DO icir=1,ncirc
1909
1910          IF( circ_class_kill_temp(icir) .LT. circ_class_n_temp(icir) ) &
1911               ldone=.FALSE. 
1912
1913       ENDDO
1914
1915       IF(ldone)EXIT ! All our biomass is dead, and we still want to kill more!
1916
1917    ENDDO
1918
1919   
1920
1921  END SUBROUTINE distribute_mortality_biomass
1922 
1923  ! added by yitong yao 07 Jan 10:05
1924       
1925
1926  ! added by yitong yao Feb 24 2020 15:47
1927  ! new distribute routine
1928
1929  SUBROUTINE distribute_mortality_biomass_yyt ( bm_difference, ddf_temp, circ_class_n_temp, &
1930       circ_class_biomass_temp, circ_class_kill_temp, circ_class_mor_temp )
1931
1932 !! 0. Variable and parameter description
1933
1934    !! 0.1 Input variables
1935    REAL(r_std),INTENT(in)                                    :: bm_difference !! the biomass to distribute
1936    REAL(r_std),INTENT(in)                                    :: ddf_temp !! the death_distribution_factor for this pft
1937    REAL(r_std),DIMENSION(:),INTENT(in)                       :: circ_class_n_temp !! circ_class_n for this point/PFT
1938    REAL(r_std),DIMENSION(:,:),INTENT(in)                     :: circ_class_biomass_temp !! circ_class_biomass for
1939                                                                          !! this point/PFT
1940    !! 0.2 Output variables
1941
1942    !! 0.3 Modified variables
1943    REAL(r_std),DIMENSION(:),INTENT(inout)                    :: circ_class_kill_temp !! circ_class_kill for
1944                                                                                      !! this point/PFT/pool
1945    REAL(r_std),DIMENSION(:),INTENT(in)                       :: circ_class_mor_temp 
1946    !! 0.4 Local variables
1947    REAL(r_std), DIMENSION(ncirc)                             :: biomass_desired    !! The biomass that dies naturally
1948    REAL(r_std), DIMENSION(ncirc)                             :: death_distribution !! The fraction of biomass taken from
1949                                                                                    !! each circ class for mortality.
1950    LOGICAL                                                   :: ldone              !! Flag to exit a loop
1951    REAL(r_std)                                               :: scale_factor       !!
1952    REAL(r_std)                                               :: sum_total          !!
1953    REAL(r_std)                                               :: leftover_bm        !! excess biomass that we need to kill
1954    REAL(r_std)                                               :: living_biomass      !! summed biomass
1955    REAL(r_std)                                               :: living_trees       !! summed biomass
1956    INTEGER                                                   :: icir
1957
1958!_ ================================================================================================================================
1959
1960    IF(ncirc == 1)THEN
1961
1962       ! This is the easy case.  All of our biomass will be taken from the only
1963       ! circumference class that we have.
1964       circ_class_kill_temp(1)=circ_class_kill_temp(1)+&
1965            bm_difference/SUM(circ_class_biomass_temp(1,:))
1966
1967       RETURN
1968    ENDIF
1969
1970    ! Here we assume an exponential distribution, arranged so that
1971    ! ddf_temp times more biomass is taken from the largest circ class
1972    ! compared to the smallest.
1973    biomass_desired(:)=zero
1974    death_distribution(:)=un
1975    scale_factor=ddf_temp**(un/REAL(ncirc-1))
1976    DO icir=2,ncirc
1977       death_distribution(icir)=death_distribution(icir-1)*scale_factor
1978    ENDDO
1979    ! Normalize it
1980    ! added by yitong yao Feb 24 2020 15:33 
1981    ! only kill circ meet the condition
1982    DO icir=1,ncirc
1983        IF (circ_class_mor_temp(icir) ==  0) THEN
1984            death_distribution(icir)=0
1985        ENDIF
1986    ENDDO 
1987    ! added by yitong yao Feb 24 2020 15:40
1988
1989    sum_total=SUM(death_distribution(:))
1990    death_distribution(:)=death_distribution(:)/sum_total
1991
1992    ! Ideally, how much is killed from each class?  Be careful to include
1993    ! what was killed in self-thinning here!  If we don't, we may try to kill
1994    ! more biomass than is available in the loop below.
1995    DO icir=1,ncirc
1996       biomass_desired(icir)=death_distribution(icir)*bm_difference !+&
1997    !        circ_class_kill_temp(icir)*SUM(circ_class_biomass_temp(icir,:))
1998    ENDDO
1999
2000    ! Right now, we know how much biomass we want to kill in each
2001    ! class (biomass_desired).  What we will do now is to loop through
2002    ! all the circ_classes and see if we have this much biomass
2003    ! in each class still alive.  The total amount of vegetation still
2004    ! alive is circ_class_n_temp(icir)-circ_class_kill_temp(icir).  If
2005    ! this number of individuals cannot give us the total biomass that
2006    ! we need, we keep track of a residual quantity, leftover_bm, which
2007    ! we try to take from other circ_classes.  It's possible that we
2008    ! will have to loop several times, which makes things more complicated.
2009
2010    ldone=.FALSE.
2011    leftover_bm=zero
2012    DO
2013       DO icir=ncirc,1,-1
2014
2015          living_trees=circ_class_n_temp(icir)-circ_class_kill_temp(icir)
2016          living_biomass=&
2017               SUM(circ_class_biomass_temp(icir,:))*living_trees
2018          biomass_desired(icir)=biomass_desired(icir)+leftover_bm
2019
2020          IF(living_biomass .LE. biomass_desired(icir))THEN
2021
2022             ! We can't get everything from this class, so we kill whatever is left.
2023             leftover_bm=biomass_desired(icir)-living_biomass
2024             biomass_desired(icir)=zero
2025             circ_class_kill_temp(icir)=circ_class_kill_temp(icir)+&
2026                  living_trees
2027
2028          ELSE
2029
2030             ! We have enough in this class.
2031             circ_class_kill_temp(icir)=circ_class_kill_temp(icir)+&
2032                  biomass_desired(icir)/SUM(circ_class_biomass_temp(icir,:))
2033             biomass_desired(icir)=zero
2034             leftover_bm=zero
2035
2036          ENDIF ! living_biomass .LE. biomass_desired(icir)
2037
2038       ENDDO ! loop over circ classes
2039
2040       IF(leftover_bm .LE. min_stomate) EXIT
2041
2042       ! it's possible that we don't have enough biomass left to kill what needs to be
2043       ! killed, so everything just dies.  I cannot think of a case where this
2044       ! would happen, though, since mortality should always be a percentage of the
2045       ! total biomass.
2046       ldone=.TRUE.
2047       DO icir=1,ncirc
2048
2049          IF( circ_class_kill_temp(icir) .LT. circ_class_n_temp(icir) ) &
2050               ldone=.FALSE. 
2051
2052       ENDDO
2053
2054       IF(ldone)EXIT ! All our biomass is dead, and we still want to kill more!
2055
2056    ENDDO
2057
2058   
2059
2060  END SUBROUTINE distribute_mortality_biomass_yyt
2061 
2062
2063 
2064
2065  ! added by yitong yao March 18 2020 20:42
2066  ! new distribute routine
2067
2068  SUBROUTINE distribute_mortality_biomass_new ( bm_difference, ddf_temp, circ_class_n_temp, &
2069       circ_class_biomass_temp, circ_class_kill_temp, circ_class_mor_temp )
2070
2071 !! 0. Variable and parameter description
2072
2073    !! 0.1 Input variables
2074    REAL(r_std),INTENT(in)                                    :: bm_difference !! the biomass to distribute
2075    REAL(r_std),INTENT(in)                                    :: ddf_temp !! the death_distribution_factor for this pft
2076    REAL(r_std),DIMENSION(:),INTENT(in)                       :: circ_class_n_temp !! circ_class_n for this point/PFT
2077    REAL(r_std),DIMENSION(:,:),INTENT(in)                     :: circ_class_biomass_temp !! circ_class_biomass for
2078                                                                          !! this point/PFT
2079    !! 0.2 Output variables
2080
2081    !! 0.3 Modified variables
2082    REAL(r_std),DIMENSION(:),INTENT(inout)                    :: circ_class_kill_temp !! circ_class_kill for
2083                                                                                      !! this point/PFT/pool
2084    REAL(r_std),DIMENSION(:),INTENT(in)                       :: circ_class_mor_temp 
2085    !! 0.4 Local variables
2086    REAL(r_std), DIMENSION(ncirc)                             :: biomass_desired    !! The biomass that dies naturally
2087    REAL(r_std), DIMENSION(ncirc)                             :: death_distribution !! The fraction of biomass taken from
2088                                                                                    !! each circ class for mortality.
2089    LOGICAL                                                   :: ldone              !! Flag to exit a loop
2090    REAL(r_std)                                               :: scale_factor       !!
2091    REAL(r_std)                                               :: sum_total          !!
2092    REAL(r_std)                                               :: leftover_bm        !! excess biomass that we need to kill
2093    REAL(r_std)                                               :: living_biomass      !! summed biomass
2094    REAL(r_std)                                               :: living_trees       !! summed biomass
2095    INTEGER                                                   :: icir
2096
2097!_ ================================================================================================================================
2098
2099    IF(ncirc == 1)THEN
2100
2101       ! This is the easy case.  All of our biomass will be taken from the only
2102       ! circumference class that we have.
2103       circ_class_kill_temp(1)=circ_class_kill_temp(1)+&
2104            bm_difference/SUM(circ_class_biomass_temp(1,:))
2105
2106       RETURN
2107    ENDIF
2108
2109    ! Here we assume an exponential distribution, arranged so that
2110    ! ddf_temp times more biomass is taken from the largest circ class
2111    ! compared to the smallest.
2112    biomass_desired(:)=zero
2113    death_distribution(:)=un
2114    scale_factor=ddf_temp**(un/REAL(ncirc-1))
2115    DO icir=2,ncirc
2116       death_distribution(icir)=death_distribution(icir-1)*scale_factor
2117    ENDDO
2118    ! Normalize it
2119    ! added by yitong yao Feb 24 2020 15:33 
2120    ! only kill circ meet the condition
2121    DO icir=1,ncirc
2122        IF (circ_class_mor_temp(icir) ==  0) THEN
2123            death_distribution(icir)=0
2124        ENDIF
2125    ENDDO 
2126    ! added by yitong yao Feb 24 2020 15:40
2127
2128    sum_total=SUM(death_distribution(:))
2129    death_distribution(:)=death_distribution(:)/sum_total
2130
2131    ! Ideally, how much is killed from each class?  Be careful to include
2132    ! what was killed in self-thinning here!  If we don't, we may try to kill
2133    ! more biomass than is available in the loop below.
2134    DO icir=1,ncirc
2135       biomass_desired(icir)=death_distribution(icir)*bm_difference !+&
2136    !        circ_class_kill_temp(icir)*SUM(circ_class_biomass_temp(icir,:))
2137    ENDDO
2138
2139    ! Right now, we know how much biomass we want to kill in each
2140    ! class (biomass_desired).  What we will do now is to loop through
2141    ! all the circ_classes and see if we have this much biomass
2142    ! in each class still alive.  The total amount of vegetation still
2143    ! alive is circ_class_n_temp(icir)-circ_class_kill_temp(icir).  If
2144    ! this number of individuals cannot give us the total biomass that
2145    ! we need, we keep track of a residual quantity, leftover_bm, which
2146    ! we try to take from other circ_classes.  It's possible that we
2147    ! will have to loop several times, which makes things more complicated.
2148
2149    ldone=.FALSE.
2150    leftover_bm=zero
2151    DO
2152       DO icir=ncirc,1,-1
2153
2154          living_trees=circ_class_n_temp(icir)-circ_class_kill_temp(icir)
2155          living_biomass=&
2156               SUM(circ_class_biomass_temp(icir,:))*living_trees
2157          biomass_desired(icir)=biomass_desired(icir)+leftover_bm
2158
2159          IF(living_biomass .LE. biomass_desired(icir))THEN
2160
2161             ! We can't get everything from this class, so we kill whatever is left.
2162             leftover_bm=biomass_desired(icir)-living_biomass
2163             biomass_desired(icir)=zero
2164             circ_class_kill_temp(icir)=circ_class_kill_temp(icir)+&
2165                  living_trees
2166
2167          ELSE
2168
2169             ! We have enough in this class.
2170             circ_class_kill_temp(icir)=circ_class_kill_temp(icir)+&
2171                  biomass_desired(icir)/SUM(circ_class_biomass_temp(icir,:))
2172             biomass_desired(icir)=zero
2173             leftover_bm=zero
2174
2175          ENDIF ! living_biomass .LE. biomass_desired(icir)
2176
2177       ENDDO ! loop over circ classes
2178
2179       IF(leftover_bm .LE. min_stomate) EXIT
2180
2181       ! it's possible that we don't have enough biomass left to kill what needs to be
2182       ! killed, so everything just dies.  I cannot think of a case where this
2183       ! would happen, though, since mortality should always be a percentage of the
2184       ! total biomass.
2185       ldone=.TRUE.
2186       DO icir=1,ncirc
2187
2188          IF( circ_class_kill_temp(icir) .LT. circ_class_n_temp(icir) ) &
2189               ldone=.FALSE. 
2190
2191       ENDDO
2192
2193       IF(ldone)EXIT ! All our biomass is dead, and we still want to kill more!
2194
2195    ENDDO
2196
2197   
2198
2199  END SUBROUTINE distribute_mortality_biomass_new
2200 
2201
2202
2203
2204  ! added by yitong yao March 18 2020 20:42
2205
2206
2207
2208
2209  ! added by yitong yao Feb 24 2020 15:47
2210
2211
2212
2213 
2214  SUBROUTINE P50_vary_class (circ_class_dia_temp,circ_class_P50_temp)
2215  !! 0. variable and parameter description
2216  REAL(r_std),DIMENSION(ncirc),INTENT(out)       :: circ_class_P50_temp
2217  REAL(r_std),DIMENSION(ncirc),INTENT(in)        :: circ_class_dia_temp
2218  INTEGER                                       :: icir
2219 
2220  DO icir=1,ncirc 
2221     IF (circ_class_dia_temp(icir) .GT. 0.50) THEN
2222        circ_class_P50_temp(icir)= -1.0 / circ_class_dia_temp(icir)
2223     ELSE
2224        circ_class_P50_temp(icir)=circ_class_dia_temp(icir) * 3.125 - 3.5625
2225     ENDIF
2226  ENDDO
2227 ! dia * 4.75 - 3.675
2228  END SUBROUtINE P50_vary_class
2229  !!
2230  ! added by yitong yao 07 Jan 10:06 above line
2231
2232!! ================================================================================================================================
2233!! SUBROUTINE  : check_biomass_sync
2234!!
2235!>\BRIEF       
2236!!
2237!! DESCRIPTION :
2238!! RECENT CHANGE(S): None
2239!!
2240!! MAIN OUTPUT VARIABLE(S):
2241!!
2242!! REFERENCE(S) : None
2243!!
2244!! FLOWCHART    : None
2245!! \n
2246!_ ================================================================================================================================
2247  SUBROUTINE check_biomass_sync ( check_point, npts, biomass, &
2248       circ_class_biomass, circ_class_n , ind, &
2249       lsync, bm_sync)
2250
2251 !! 0. Variable and parameter description
2252
2253    !! 0.1 Input variables
2254    INTEGER(i_std), INTENT(in)                                 :: npts                 !! Domain size (unitless)
2255    REAL(r_std), DIMENSION(:,:,:,:,:), INTENT(in)              :: circ_class_biomass   !! Biomass of the componets of the model 
2256                                                                                       !! tree within a circumference
2257                                                                                       !! class @tex $(gC ind^{-1})$ @endtex 
2258    REAL(r_std), DIMENSION(:,:,:), INTENT(in)                  :: circ_class_n         !! Number of individuals in each circ class
2259                                                                                       !! @tex $(m^{-2})$ @endtex
2260    REAL(r_std), DIMENSION(:,:,:,:), INTENT(in)                :: biomass              !! Stand level biomass
2261                                                                                       !! @tex $(gC m^{-2})$ @endtex
2262    CHARACTER(*),INTENT(in)                                    :: check_point          !! A flag to indicate at which
2263                                                                                       !! point in the code we're doing
2264                                                                                       !! this check
2265    REAL(r_std), DIMENSION(:,:), INTENT(in)                    :: ind                  !! Density of individuals
2266                                                                                       !! @tex $(m^{-2})$ @endtex
2267
2268    !! 0.2 Output variables
2269    LOGICAL,INTENT(out)                                        :: lsync
2270    REAL(r_std), DIMENSION(:,:,:), INTENT(out)                 :: bm_sync              !! The difference betweeen the
2271                                                                                       !! biomass in the circ_classes and
2272                                                                                       !! the total biomass
2273                                                                                       !! @tex $(gC m^{-2})$ @endtex
2274    !! 0.3 Modified variables
2275
2276    !! 0.4 Local variables
2277    INTEGER                                                    :: iele,ipts,ivm,ipar,icir
2278    REAL(r_std)                                                :: total_circ_class_biomass
2279    REAL(r_std),DIMENSION(ncirc)                               :: tree_size
2280    LOGICAL                                                    :: lnegative
2281   
2282!_ ================================================================================================================================
2283
2284    lsync=.TRUE.
2285    lnegative=.FALSE.
2286
2287    bm_sync(:,:,:)=zero
2288
2289    !++++++ TEMP ++++++
2290    ! We gain 5-10% speed by skipping this routine
2291       
2292    !++++++++++++
2293
2294    ! Check to see if the biomass is not equal to the total biomass
2295    ! in circ_class_biomass anywhere.
2296    DO ipts=1,npts
2297
2298       DO ivm=1,nvm
2299
2300          ! Only woody PFTs have circumference classes therefore
2301          ! only woody PFTs need to be syncronized
2302          IF(.NOT. lbypass_cc)THEN
2303             IF(is_tree(ivm)) THEN
2304                tree_size(:)=zero
2305                DO icir=1,ncirc
2306                   tree_size(icir)=SUM(circ_class_biomass(ipts,ivm,icir,:,1))
2307                ENDDO
2308                DO icir=2,ncirc
2309                   IF(tree_size(icir) .LT. tree_size(icir-1)-min_stomate)THEN
2310                      WRITE(numout,*) 'ERROR: stopping in sync'
2311                      WRITE(numout,*) check_point
2312                      WRITE(numout,*) 'ipts,ivm: ',ipts,ivm
2313                      WRITE(numout,*) 'tree_size(icir), tree_size(icir-1), ',&
2314                           tree_size(icir), tree_size(icir-1), tree_size(icir) - tree_size(icir-1) 
2315                      WRITE(numout,'(A,10F20.10)') 'icir, tree_size: ',icir, tree_size(:)
2316                      CALL ipslerr_p (3,'check_biomass_sync', &
2317                           'The size of the trees in the circ class are not monotonically increasing!',&
2318                           'Look in the output file for more details.',&
2319                           '')
2320                   ENDIF
2321                ENDDO
2322             ENDIF
2323          ENDIF
2324         
2325          DO iele=1,icarbon
2326               
2327             DO ipar=1,nparts
2328               
2329                total_circ_class_biomass=zero
2330                DO icir=1,ncirc
2331                   
2332                   total_circ_class_biomass=total_circ_class_biomass+&
2333                        circ_class_biomass(ipts,ivm,icir,ipar,iele)*circ_class_n(ipts,ivm,icir)
2334
2335                   ! Check as well to see if our biomass is ever negative.
2336                   ! It really should not be.
2337                   IF(circ_class_biomass(ipts,ivm,icir,ipar,iele) .LT. -min_stomate)THEN
2338
2339                      lnegative=.TRUE.
2340                      WRITE(numout,*) '!***********************************'
2341                      WRITE(numout,*) 'Error: Negative biomass component!'
2342                      WRITE(numout,*) 'Check point: ',TRIM(check_point)
2343                      WRITE(numout,*) 'circ_class_biomass(ipts,ivm,icir,ipar,iele) ',&
2344                           circ_class_biomass(ipts,ivm,icir,ipar,iele)
2345                      WRITE(numout,'(A,5I5)') 'ipts,ivm,icir,ipar,iele',ipts,ivm,icir,ipar,iele
2346                      WRITE(numout,*) '!***********************************'
2347
2348                   ENDIF
2349                ENDDO
2350
2351                IF(ABS(biomass(ipts,ivm,ipar,iele) -  &
2352                     total_circ_class_biomass) .GT. sync_threshold)THEN
2353
2354                   WRITE(numout,*) '!***********************************'
2355                   WRITE(numout,*) 'Biomass and circ_class_biomass are not equal!'
2356                   WRITE(numout,*) 'Check point: ',TRIM(check_point)
2357                   WRITE(numout,100) 'biomass(ipts,ivm,ipar,iele) ',&
2358                        biomass(ipts,ivm,ipar,iele)
2359                   WRITE(numout,100) 'total_circ_class_biomass ',&
2360                        total_circ_class_biomass
2361                   WRITE(numout,100) 'Difference: ',&
2362                        ABS(biomass(ipts,ivm,ipar,iele) - total_circ_class_biomass)
2363                   WRITE(numout,*) 'ipts,ivm,ipar,iele',ipts,ivm,ipar,iele
2364                   WRITE(numout,*) '!***********************************'
2365100                FORMAT(A,E20.10)
2366                   lsync=.FALSE.
2367
2368                ENDIF
2369
2370             ENDDO
2371
2372             ! we are not going to save the biomass for every component right now,
2373             ! just the total
2374             bm_sync(ipts,ivm,iele)=zero
2375
2376             DO ipar=1,nparts
2377                   
2378
2379                DO icir=1,ncirc
2380
2381                   bm_sync(ipts,ivm,iele)=bm_sync(ipts,ivm,iele)+&
2382                        circ_class_biomass(ipts,ivm,icir,ipar,iele)*circ_class_n(ipts,ivm,icir)
2383                ENDDO ! ncirc
2384
2385             ENDDO ! nparts
2386
2387             bm_sync(ipts,ivm,iele)=ABS(bm_sync(ipts,ivm,iele)-&
2388                  SUM(biomass(ipts,ivm,:,iele)))
2389
2390          ENDDO ! nelements
2391
2392
2393
2394       ENDDO ! loop over PFTs
2395
2396    ENDDO ! loop over points
2397 
2398    !---TEMP---
2399    IF(ld_biomass)THEN
2400       WRITE(numout,*) 'Check point: ',TRIM(check_point)
2401       WRITE(numout,*) 'test_pft, test_grid: ',test_pft,test_grid
2402       WRITE(numout,*) 'biomass (ileaf), ', biomass(test_grid,test_pft,ileaf,icarbon)
2403       WRITE(numout,*) 'biomass (iwood), ', biomass(test_grid,test_pft,isapabove,icarbon) + &
2404           biomass(test_grid,test_pft,isapbelow,icarbon) + biomass(test_grid,test_pft,iheartabove,icarbon) + &
2405           biomass(test_grid,test_pft,iheartbelow,icarbon)
2406       WRITE(numout,*) 'biomass (iroot), ', biomass(test_grid,test_pft,iroot,icarbon)
2407       WRITE(numout,'(A,20F14.6)') 'biomassHHH, ',biomass(test_grid,test_pft,:,icarbon)
2408       DO icir=1,ncirc
2409          WRITE(numout,'(A,I1,20F14.6)') 'ccbiomass',icir,circ_class_biomass(test_grid,test_pft,icir,:,icarbon)
2410       ENDDO
2411       WRITE(numout,*) 'circ_class_biomass, ',&
2412            SUM (SUM(circ_class_biomass(test_grid,test_pft,:,:,icarbon),2) * &
2413            circ_class_n(test_grid,test_pft,:))
2414       WRITE(numout,*) 'circ_class_n, ', SUM(circ_class_n(test_grid,test_pft,:))
2415       WRITE(numout,*) 'circ_class_n(:), ', circ_class_n(test_grid,test_pft,:)
2416       WRITE(numout,*) 'ind, ', ind(test_grid,test_pft)
2417    ENDIF
2418
2419!!$    !----------
2420
2421    IF(.NOT. lsync) THEN
2422       WRITE(numout,*) 'ERROR: stopping in sync #2'
2423       WRITE(numout,*) 'Stopping'
2424       CALL ipslerr_p (3,'check_biomass_sync', &
2425            'circ_class_biomass*circ_class_n is not equal to the total biomass',&
2426            'Look in the output file for more details.',&
2427            '')
2428
2429    ENDIF
2430    IF(lnegative) THEN
2431       WRITE(numout,*) 'ERROR: negative biomass'
2432       WRITE(numout,*) 'Stopping'
2433       CALL ipslerr_p (3,'check_biomass_sync', &
2434            'One of the biomass pools is negative!',&
2435            'Look in the output file for more details.',&
2436            '')
2437    ENDIF
2438
2439  END SUBROUTINE check_biomass_sync
2440
2441END MODULE function_library
2442
Note: See TracBrowser for help on using the repository browser.