source: tags/ORCHIDEE_1_9_5/ORCHIDEE/src_stomate/stomate_lpj.f90 @ 8

Last change on this file since 8 was 8, checked in by orchidee, 14 years ago

import first tag equivalent to CVS orchidee_1_9_5 + OOL_1_9_5
File size: 44.0 KB
Line 
1! Stomate: phenology, allocation, etc.
2!
3! authors: A. Botta, P. Friedlingstein, C. Morphopoulos, N. Viovy, et al.
4!
5! bits and pieces put together by G. Krinner
6!
7! version 0.0: August 1998
8!
9! $Header: /home/ssipsl/CVSREP/ORCHIDEE/src_stomate/stomate_lpj.f90,v 1.26 2010/08/05 16:02:37 ssipsl Exp $
10! IPSL (2006)
11!  This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
12!
13MODULE stomate_lpj
14
15  ! modules used:
16
17  USE ioipsl
18  USE grid
19  USE stomate_constants
20  USE lpj_constraints
21  USE lpj_pftinout
22  USE lpj_kill
23  USE lpj_crown
24  USE lpj_fire
25  USE lpj_gap
26  USE lpj_light
27  USE lpj_establish
28  USE lpj_cover
29  USE stomate_prescribe
30  USE stomate_phenology
31  USE stomate_alloc
32  USE stomate_npp
33  USE stomate_turnover
34  USE stomate_litter
35  USE stomate_soilcarbon
36  USE stomate_vmax
37  USE stomate_assimtemp
38  USE constantes_veg
39  USE stomate_lcchange
40  !  USE Write_Field_p
41
42  IMPLICIT NONE
43
44  ! private & public routines
45
46  PRIVATE
47  PUBLIC StomateLpj,StomateLpj_clear
48
49  ! first call
50  LOGICAL, SAVE                                           :: firstcall = .TRUE.
51
52CONTAINS
53
54  SUBROUTINE StomateLpj_clear
55
56    CALL prescribe_clear
57    CALL phenology_clear
58    CALL npp_calc_clear
59    CALL turn_clear
60    CALL soilcarbon_clear
61    CALL constraints_clear
62    CALL establish_clear
63    CALL fire_clear
64    CALL gap_clear
65    CALL light_clear
66    CALL pftinout_clear
67    CALL alloc_clear
68  END SUBROUTINE StomateLpj_clear
69
70  SUBROUTINE StomateLpj (npts, dt_days, EndOfYear, &
71       neighbours, resolution, &
72       clay, herbivores, &
73       tsurf_daily, tsoil_daily, t2m_daily, t2m_min_daily, &
74       litterhum_daily, soilhum_daily, &
75       maxmoiavail_lastyear, minmoiavail_lastyear, &
76       gdd0_lastyear, precip_lastyear, &
77       moiavail_month, moiavail_week, tlong_ref, t2m_month, t2m_week, &
78       tsoil_month, soilhum_month, &
79       gdd_m5_dormance, gdd_midwinter, ncd_dormance, ngd_minus5, &
80       turnover_longterm, gpp_daily, time_lowgpp, &
81       time_hum_min, maxfpc_lastyear, resp_maint_part, &
82       PFTpresent, age, fireindex, firelitter, &
83       leaf_age, leaf_frac, biomass, ind, adapted, regenerate, &
84       senescence, when_growthinit, &
85       litterpart, litter, dead_leaves, carbon, black_carbon, lignin_struc, &
86       veget_max, veget, npp_longterm, lm_lastyearmax, veget_lastlight, &
87       everywhere, need_adjacent, RIP_time, &
88       lai, rprof,npp_daily, turnover_daily, turnover_time,&
89       control_moist, control_temp, soilcarbon_input, &
90       co2_to_bm, co2_fire, resp_hetero, resp_maint, resp_growth, &
91       height, deadleaf_cover, vcmax, vjmax, &
92       t_photo_min, t_photo_opt, t_photo_max,bm_to_litter, &
93       prod10,prod100,flux10, flux100, veget_max_new, &
94       convflux,cflux_prod10,cflux_prod100, harvest_above, lcchange)
95
96    !
97    ! 0 declarations
98    !
99
100    ! 0.1 input
101
102    ! Domain size
103    INTEGER(i_std), INTENT(in)                                           :: npts
104    ! time step of Stomate in days
105    REAL(r_std), INTENT(in)                                        :: dt_days
106    ! indices of the 8 neighbours of each grid point (1=N, 2=NE, 3=E, 4=SE, 5=S, 6=SW, 7=W, 8=NW)
107    INTEGER(i_std), DIMENSION(npts,8), INTENT(in)                 :: neighbours
108    ! resolution at each grid point in m (1=E-W, 2=N-S)
109    REAL(r_std), DIMENSION(npts,2), INTENT(in)                     :: resolution
110    ! clay fraction
111    REAL(r_std), DIMENSION(npts), INTENT(in)                       :: clay
112    ! time constant of probability of a leaf to be eaten by a herbivore (days)
113    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)                  :: herbivores
114    ! daily surface temperatures (K)
115    REAL(r_std), DIMENSION(npts), INTENT(in)                       :: tsurf_daily
116    ! daily soil temperatures (K)
117    REAL(r_std), DIMENSION(npts,nbdl), INTENT(in)                  :: tsoil_daily
118    ! daily 2 meter temperatures (K)
119    REAL(r_std), DIMENSION(npts), INTENT(in)                       :: t2m_daily
120    ! daily minimum 2 meter temperatures (K)
121    REAL(r_std), DIMENSION(npts), INTENT(in)                       :: t2m_min_daily
122    ! daily litter humidity
123    REAL(r_std), DIMENSION(npts), INTENT(in)                       :: litterhum_daily
124    ! daily soil humidity
125    REAL(r_std), DIMENSION(npts,nbdl), INTENT(in)                  :: soilhum_daily
126    ! last year's maximum moisture availability
127    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)                  :: maxmoiavail_lastyear
128    ! last year's minimum moisture availability
129    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)                  :: minmoiavail_lastyear
130    ! last year's GDD0
131    REAL(r_std), DIMENSION(npts), INTENT(in)                       :: gdd0_lastyear
132    ! lastyear's precipitation (mm/year)
133    REAL(r_std), DIMENSION(npts), INTENT(in)                       :: precip_lastyear
134    ! "monthly" moisture availability
135    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)                  :: moiavail_month
136    ! "weekly" moisture availability
137    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)                  :: moiavail_week
138    ! "long term" 2 meter reference temperatures (K)
139    REAL(r_std), DIMENSION(npts), INTENT(in)                       :: tlong_ref
140    ! "monthly" 2-meter temperatures (K)
141    REAL(r_std), DIMENSION(npts), INTENT(in)                       :: t2m_month
142    ! "weekly" 2-meter temperatures (K)
143    REAL(r_std), DIMENSION(npts), INTENT(in)                       :: t2m_week
144    ! "monthly" soil temperatures (K)
145    REAL(r_std), DIMENSION(npts,nbdl), INTENT(in)                  :: tsoil_month
146    ! "monthly" soil humidity
147    REAL(r_std), DIMENSION(npts,nbdl), INTENT(in)                  :: soilhum_month
148    ! growing degree days, threshold -5 deg C (for phenology)
149    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)                   :: gdd_m5_dormance
150    ! growing degree days, since midwinter (for phenology)
151    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)                   :: gdd_midwinter
152    ! number of chilling days, since leaves were lost (for phenology)
153    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)                  :: ncd_dormance
154    ! number of growing days, threshold -5 deg C (for phenology)
155    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)                  :: ngd_minus5
156    ! "long term" turnover rate (gC/(m**2 of ground)/year)
157    REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(in)           :: turnover_longterm
158    ! daily gross primary productivity (gC/(m**2 of ground)/day)
159    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)                  :: gpp_daily
160    ! duration of dormance (d)
161    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)                  :: time_lowgpp
162    ! time elapsed since strongest moisture availability (d)
163    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)                  :: time_hum_min
164    ! last year's maximum fpc for each natural PFT, on ground
165    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)                  :: maxfpc_lastyear
166    ! maintenance respiration of different plant parts (gC/day/m**2 of ground)
167    REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(in)             :: resp_maint_part
168
169    ! 0.2 modified fields
170
171    ! PFT exists
172    LOGICAL, DIMENSION(npts,nvm), INTENT(inout)                  :: PFTpresent
173    ! age (years)
174    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)               :: age
175    ! Probability of fire
176    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)          :: fireindex
177    ! Longer term litter above the ground, gC/m**2 of ground
178    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)          :: firelitter
179    ! leaf age (days)
180    REAL(r_std), DIMENSION(npts,nvm,nleafages), INTENT(inout)     :: leaf_age
181    ! fraction of leaves in leaf age class
182    REAL(r_std), DIMENSION(npts,nvm,nleafages), INTENT(inout)     :: leaf_frac
183    ! biomass (gC/(m**2 of ground))
184    REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(inout)        :: biomass
185    ! density of individuals (1/(m**2 of ground))
186    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)               :: ind
187    ! Winter too cold? between 0 and 1
188    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)               :: adapted
189    ! Winter sufficiently cold? between 0 and 1
190    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)               :: regenerate
191    ! is the plant senescent? (only for deciduous trees - carbohydrate reserve)
192    LOGICAL, DIMENSION(npts,nvm), INTENT(inout)                  :: senescence
193    ! how many days ago was the beginning of the growing season
194    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)               :: when_growthinit
195    ! fraction of litter above the ground belonging to different PFTs
196    REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(inout)         :: litterpart
197    ! metabolic and structural litter, above and below ground (gC/(m**2 of ground))
198    REAL(r_std), DIMENSION(npts,nlitt,nvm,nlevs), INTENT(inout)  :: litter
199    ! dead leaves on ground, per PFT, metabolic and structural,
200    !   in gC/(m**2 of ground)
201    REAL(r_std), DIMENSION(npts,nvm,nlitt), INTENT(inout)         :: dead_leaves
202    ! carbon pool: active, slow, or passive,(gC/(m**2 of ground))
203    REAL(r_std), DIMENSION(npts,ncarb,nvm), INTENT(inout)    :: carbon
204    ! black carbon on the ground (gC/(m**2 of total ground))
205    REAL(r_std), DIMENSION(npts), INTENT(inout)                    :: black_carbon
206    ! ratio Lignine/Carbon in structural litter, above and below ground,
207    ! (gC/(m**2 of ground))
208    REAL(r_std), DIMENSION(npts,nvm,nlevs), INTENT(inout)    :: lignin_struc
209    ! "maximal" coverage fraction of a PFT (LAI -> infinity) on ground
210    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)               :: veget_max
211    ! fractional coverage on ground, taking into account LAI (=grid-scale fpc)
212    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)               :: veget
213    ! "long term" net primary productivity (gC/(m**2 of ground)/year)
214    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)               :: npp_longterm
215    ! last year's maximum leaf mass, for each PFT (gC/(m**2 of ground))
216    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)               :: lm_lastyearmax
217    ! vegetation fractions (on ground) after last light competition
218    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)               :: veget_lastlight
219    ! is the PFT everywhere in the grid box or very localized (after its introduction)
220    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)               :: everywhere
221    ! in order for this PFT to be introduced, does it have to be present in an
222    !   adjacent grid box?
223    LOGICAL, DIMENSION(npts,nvm), INTENT(inout)                  :: need_adjacent
224    ! How much time ago was the PFT eliminated for the last time (y)
225    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)               :: RIP_time
226    ! Turnover_time of leaves for grasses (d)
227    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)          :: turnover_time
228
229    ! 0.3 output
230
231    ! leaf area index OF AN INDIVIDUAL PLANT
232    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)                 :: lai
233    ! root depth. This will, one day, be a prognostic variable. It will be calculated by
234    ! STOMATE (save in restart file & give to hydrology module!), probably somewhere
235    ! in the allocation routine. For the moment, it is prescribed.
236    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)                 :: rprof
237    ! net primary productivity (gC/day/(m**2 of ground))
238    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)                 :: npp_daily
239    ! Turnover rates (gC/(m**2 of ground)/day)
240    REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(out)          :: turnover_daily
241    ! moisture control of heterotrophic respiration
242    REAL(r_std), DIMENSION(npts,nlevs), INTENT(inout)                :: control_moist
243    ! temperature control of heterotrophic respiration, above and below
244    REAL(r_std), DIMENSION(npts,nlevs), INTENT(inout)                :: control_temp
245    ! quantity of carbon going into carbon pools from litter decomposition
246    !   (gC/(m**2 of ground)/day)
247    REAL(r_std), DIMENSION(npts,ncarb,nvm), INTENT(inout)      :: soilcarbon_input
248    ! co2 taken up (gC/(m**2 of total ground)/day)
249    !NV devient 2D
250    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)                      :: co2_to_bm
251    ! carbon emitted into the atmosphere by fire (living and dead biomass)
252    ! (in gC/m**2 of total ground/day)
253    !NV devient 2D
254    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)                      :: co2_fire
255    ! heterotrophic respiration (gC/day/m**2 of total ground)
256    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: resp_hetero
257    ! maintenance respiration (gC/day/(m**2 of total ground))
258    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)                 :: resp_maint
259    ! growth respiration (gC/day/(m**2 of total ground))
260    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)                 :: resp_growth
261    ! height of vegetation (m)
262    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)                 :: height
263    ! fraction of soil covered by dead leaves
264    REAL(r_std), DIMENSION(npts), INTENT(out)                      :: deadleaf_cover
265    ! Maximum rate of carboxylation
266    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)                 :: vcmax
267    ! Maximum rate of RUbp regeneration
268    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)                 :: vjmax
269    ! Minimum temperature for photosynthesis (deg C)
270    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)                 :: t_photo_min
271    ! Optimum temperature for photosynthesis (deg C)
272    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)                 :: t_photo_opt
273    ! Maximum temperature for photosynthesis (deg C)
274    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)                 :: t_photo_max
275    ! conversion of biomass to litter (gC/(m**2 of ground)) / day
276    REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(out)          :: bm_to_litter
277    !
278    ! new "maximal" coverage fraction of a PFT (LAI -> infinity)
279    REAL(r_std), DIMENSION(npts,nvm),INTENT(inout)                 :: veget_max_new 
280
281    ! products remaining in the 10/100 year-turnover pool after the annual release for each compartment
282    ! (10 or 100 + 1 : input from year of land cover change)
283    REAL(r_std),DIMENSION(npts,0:10), INTENT(inout)                        :: prod10
284    REAL(r_std),DIMENSION(npts,0:100), INTENT(inout)                       :: prod100
285    ! annual release from the 10/100 year-turnover pool compartments
286    REAL(r_std),DIMENSION(npts,10), INTENT(inout)                       :: flux10
287    REAL(r_std),DIMENSION(npts,100), INTENT(inout)                      :: flux100
288    ! release during first year following land cover change
289    REAL(r_std),DIMENSION(npts), INTENT(inout)                          :: convflux
290    ! total annual release from the 10/100 year-turnover pool
291    REAL(r_std),DIMENSION(npts), INTENT(inout)                          :: cflux_prod10, cflux_prod100
292    ! harvest above ground biomass for agriculture
293    REAL(r_std), DIMENSION(npts), INTENT(inout)                       :: harvest_above
294
295    ! land cover change flag
296    LOGICAL, INTENT(in)                                                :: lcchange
297
298    ! Land cover change variables + EndOfYear
299    ! Do update of yearly variables? This variable must be .TRUE. once a year
300    LOGICAL, INTENT(in)                                                :: EndOfYear
301
302
303    ! 0.4 local
304
305    ! total conversion of biomass to litter (gC/(m**2)) / day
306    REAL(r_std), DIMENSION(npts,nvm)                                   :: tot_bm_to_litter
307    ! total living biomass (gC/(m**2))
308    REAL(r_std), DIMENSION(npts,nvm)                                   :: tot_live_biomass
309    ! biomass increase, i.e. NPP per plant part
310    REAL(r_std), DIMENSION(npts,nvm,nparts)                            :: bm_alloc
311    ! total turnover rate (gC/(m**2)) / day
312    REAL(r_std), DIMENSION(npts,nvm)                                   :: tot_turnover
313    ! total soil and litter carbon (gC/(m**2))
314    REAL(r_std), DIMENSION(npts,nvm)                                   :: tot_litter_soil_carb
315    ! total litter carbon (gC/(m**2))
316    REAL(r_std), DIMENSION(npts,nvm)                                   :: tot_litter_carb
317    ! total soil carbon (gC/(m**2))
318    REAL(r_std), DIMENSION(npts,nvm)                                   :: tot_soil_carb
319    ! crown area of individuals (m**2)
320    REAL(r_std), DIMENSION(npts,nvm)                               :: cn_ind
321    ! fraction that goes into plant part
322    REAL(r_std), DIMENSION(npts,nvm,nparts)                        :: f_alloc
323    ! space availability for trees
324    REAL(r_std), DIMENSION(npts)                                   :: avail_tree
325    ! space availability for grasses
326    REAL(r_std), DIMENSION(npts)                                   :: avail_grass
327
328    INTEGER                                                       :: j
329
330    ! total products remaining in the pool after the annual release
331    REAL(r_std),DIMENSION(npts)                                   :: prod10_total, prod100_total 
332    ! total flux from conflux and the 10/100 year-turnover pool
333    REAL(r_std),DIMENSION(npts)                                       :: cflux_prod_total
334
335    ! "maximal" coverage fraction of a PFT (LAI -> infinity) on ground
336    REAL(r_std),DIMENSION(npts,nvm)                                :: veget_max_old
337
338    REAL(r_std), DIMENSION(npts)                                   :: vartmp
339
340    REAL(r_std), DIMENSION(npts,nvm)                          :: histvar
341
342    ! =========================================================================
343
344    IF (bavard.GE.3) WRITE(numout,*) 'Entering stomate_lpj'
345
346    !
347    ! 1 Initializations
348    !
349
350    !
351    ! 1.1 set outputs to zero
352    !
353    co2_to_bm(:,:) = zero
354    co2_fire(:,:) = zero
355    npp_daily(:,:) = zero
356    turnover_daily(:,:,:) = zero
357    resp_maint(:,:) = zero
358    resp_growth(:,:) = zero
359    harvest_above(:) = zero
360
361    !
362    ! 1.2  initialize some variables
363    !
364
365    bm_to_litter(:,:,:) = zero
366    cn_ind(:,:) = zero
367    veget_max_old(:,:) = veget_max(:,:)
368
369    !
370    ! 1.3 Prescribe some vegetation characteristics if the vegetation is not dynamic
371    !     IF the DGVM is not activated, the density of individuals and their crown
372    !     areas don't matter, but they should be defined for the case we switch on
373    !     the DGVM afterwards.
374    !     At first call, if the DGVM is not activated, impose a minimum biomass for
375    !       prescribed PFTs and declare them present.
376    !
377
378    CALL prescribe (npts, &
379         veget_max, PFTpresent, everywhere, when_growthinit, &
380         biomass, leaf_frac, ind, cn_ind)
381
382    !
383    ! 2 climatic constraints for PFT presence and regenerativeness
384    !   call this even when DGVM is not activated so that "adapted" and "regenerate"
385    !   are kept up to date for the moment when the DGVM is activated.
386    !
387
388    CALL constraints (npts, dt_days, &
389         t2m_month, t2m_min_daily, when_growthinit, &
390         adapted, regenerate)
391
392    !
393    ! 3 PFTs in and out, based on climate criteria
394    !
395
396    IF ( control%ok_dgvm ) THEN
397
398       !
399       ! 3.1 do introduction/elimination
400       !
401
402       CALL pftinout (npts, dt_days, adapted, regenerate, &
403            neighbours, veget, veget_max, &
404            biomass, ind, age, leaf_frac, npp_longterm, lm_lastyearmax, senescence, &
405            PFTpresent, everywhere, when_growthinit, need_adjacent, RIP_time, &
406            co2_to_bm, &
407            avail_tree, avail_grass)
408
409       !
410       ! 3.2 reset attributes for eliminated PFTs.
411       !     This also kills PFTs that had 0 leafmass during the last year. The message
412       !     "... after pftinout" is misleading in this case.
413       !
414
415       CALL kill (npts, 'pftinout  ', lm_lastyearmax, &
416            ind, PFTpresent, cn_ind, biomass, senescence, RIP_time, &
417            lai, age, leaf_age, leaf_frac, &
418            when_growthinit, everywhere, veget, veget_max, bm_to_litter)
419
420       !
421       ! 3.3 calculate new crown area and maximum vegetation cover
422       !
423
424       CALL crown (npts, PFTpresent, &
425            ind, biomass, &
426            veget_max, cn_ind, height)
427
428    ENDIF
429
430    !
431    ! 4 phenology
432    !
433    CALL histwrite (hist_id_stomate, 'WHEN_GROWTHINIT', itime, when_growthinit, npts*nvm, horipft_index)
434    CALL histwrite (hist_id_stomate, 'TIME_LOWGPP', itime, time_lowgpp, npts*nvm, horipft_index)
435
436    WHERE(PFTpresent)
437       histvar=un
438    ELSEWHERE
439       histvar=zero
440    ENDWHERE
441    CALL histwrite (hist_id_stomate, 'PFTPRESENT', itime, histvar, npts*nvm, horipft_index)
442
443    WHERE(gdd_midwinter.EQ.undef)
444       histvar=val_exp
445    ELSEWHERE
446       histvar=gdd_midwinter
447    ENDWHERE
448    CALL histwrite (hist_id_stomate, 'GDD_MIDWINTER', itime, histvar, npts*nvm, horipft_index)
449
450    WHERE(ncd_dormance.EQ.undef)
451       histvar=val_exp
452    ELSEWHERE
453       histvar=ncd_dormance
454    ENDWHERE
455    CALL histwrite (hist_id_stomate, 'NCD_DORMANCE', itime, histvar, npts*nvm, horipft_index)
456
457    CALL phenology (npts, dt_days, PFTpresent, &
458         veget_max, &
459         tlong_ref, t2m_month, t2m_week, gpp_daily, &
460         maxmoiavail_lastyear, minmoiavail_lastyear, &
461         moiavail_month, moiavail_week, &
462         gdd_m5_dormance, gdd_midwinter, ncd_dormance, ngd_minus5, &
463         senescence, time_lowgpp, time_hum_min, &
464         biomass, leaf_frac, leaf_age, &
465         when_growthinit, co2_to_bm, lai)
466
467    !
468    ! 5 allocation
469    !
470
471    CALL alloc (npts, dt_days, &
472         lai, veget_max, senescence, when_growthinit, &
473         moiavail_week, tsoil_month, soilhum_month, &
474         biomass, age, leaf_age, leaf_frac, rprof, f_alloc)
475
476    !
477    ! 6 maintenance and growth respiration. NPP
478    !
479
480    CALL npp_calc (npts, dt_days, &
481         PFTpresent, &
482         tlong_ref, t2m_daily, tsoil_daily, lai, rprof, &
483         gpp_daily, f_alloc, bm_alloc, resp_maint_part,&
484         biomass, leaf_age, leaf_frac, age, &
485         resp_maint, resp_growth, npp_daily)
486
487    IF ( control%ok_dgvm ) THEN
488
489       ! new provisional crown area and maximum vegetation cover after growth
490
491       CALL crown (npts, PFTpresent, &
492            ind, biomass, &
493            veget_max, cn_ind, height)
494
495    ENDIF
496
497    !
498    ! 7 fire.
499    !
500
501    CALL fire (npts, dt_days, litterpart, &
502         litterhum_daily, t2m_daily, lignin_struc, &
503         fireindex, firelitter, biomass, ind, &
504         litter, dead_leaves, bm_to_litter, black_carbon, &
505         co2_fire)
506
507    IF ( control%ok_dgvm ) THEN
508
509       ! reset attributes for eliminated PFTs
510
511       CALL kill (npts, 'fire      ', lm_lastyearmax, &
512            ind, PFTpresent, cn_ind, biomass, senescence, RIP_time, &
513            lai, age, leaf_age, leaf_frac, &
514            when_growthinit, everywhere, veget, veget_max, bm_to_litter)
515
516    ENDIF
517
518    !
519    ! 8 tree mortality. Does not depend on age, therefore does not change crown area.
520    !
521
522    CALL gap (npts, dt_days, &
523         npp_longterm, turnover_longterm, lm_lastyearmax, &
524         PFTpresent, biomass, ind, bm_to_litter)
525
526    IF ( control%ok_dgvm ) THEN
527
528       ! reset attributes for eliminated PFTs
529
530       CALL kill (npts, 'gap       ', lm_lastyearmax, &
531            ind, PFTpresent, cn_ind, biomass, senescence, RIP_time, &
532            lai, age, leaf_age, leaf_frac, &
533            when_growthinit, everywhere, veget, veget_max, bm_to_litter)
534
535    ENDIF
536
537    !
538    ! 9 calculate vcmax, vjmax and photosynthesis temperatures
539    !
540
541    CALL vmax (npts, dt_days, &
542         leaf_age, leaf_frac, &
543         vcmax, vjmax)
544
545    CALL assim_temp (npts, tlong_ref, t2m_month, &
546         t_photo_min, t_photo_opt, t_photo_max)
547
548    !
549    ! 10 leaf senescence and other turnover processes. New lai
550    !
551
552    CALL turn (npts, dt_days, PFTpresent, &
553         herbivores, &
554         maxmoiavail_lastyear, minmoiavail_lastyear, &
555         moiavail_week,  moiavail_month,tlong_ref, t2m_month, t2m_week, veget_max, &
556         leaf_age, leaf_frac, age, lai, biomass, &
557         turnover_daily, senescence,turnover_time)
558
559    !
560    ! 11 light competition
561    !
562
563    IF ( control%ok_dgvm ) THEN
564
565       !
566       ! 11.1 do light competition
567       !
568
569       CALL light (npts, dt_days, &
570            PFTpresent, cn_ind, lai, maxfpc_lastyear, &
571            ind, biomass, veget_lastlight, bm_to_litter)
572
573       !
574       ! 11.2 reset attributes for eliminated PFTs
575       !
576
577       CALL kill (npts, 'light     ', lm_lastyearmax, &
578            ind, PFTpresent, cn_ind, biomass, senescence, RIP_time, &
579            lai, age, leaf_age, leaf_frac, &
580            when_growthinit, everywhere, veget, veget_max, bm_to_litter)
581
582    ENDIF
583
584    !
585    ! 12 establishment of saplings
586    !
587
588    IF ( control%ok_dgvm ) THEN
589
590       !
591       ! 12.1 do establishment
592       !
593
594       CALL establish (npts, dt_days, PFTpresent, regenerate, &
595            neighbours, resolution, need_adjacent, herbivores, &
596            precip_lastyear, gdd0_lastyear, lm_lastyearmax, &
597            cn_ind, lai, avail_tree, avail_grass, &
598            leaf_age, leaf_frac, &
599            ind, biomass, age, everywhere, co2_to_bm, veget_max)
600
601       !
602       ! 12.2 calculate new crown area (and maximum vegetation cover)
603       !
604
605       CALL crown (npts, PFTpresent, &
606            ind, biomass, &
607            veget_max, cn_ind, height)
608
609    ENDIF
610
611    !
612    ! 13 calculate final LAI and vegetation cover.
613    !
614
615    CALL cover (npts, cn_ind, ind, biomass, &
616         veget_max, veget_max_old, veget, &
617         lai, litter, carbon)
618
619    !
620    ! 14 the whole litter stuff:
621    !    litter update, lignin content, PFT parts, litter decay,
622    !    litter heterotrophic respiration, dead leaf soil cover.
623    !    No vertical discretisation in the soil for litter decay.
624    !
625    ! added by shilong for harvest
626    IF(harvest_agri) THEN
627       CALL harvest(npts, dt_days, veget_max, veget, &
628            bm_to_litter, turnover_daily, &
629            harvest_above)
630    ENDIF
631
632    ! 15.1  Land cover change
633
634    !shilong adde turnover_daily
635    IF(EndOfYear) THEN
636       IF (lcchange) THEN
637          CALL lcchange_main (npts, dt_days, veget_max, veget_max_new, &
638               biomass, ind, age, PFTpresent, senescence, when_growthinit, &
639               everywhere, veget, & 
640               co2_to_bm, bm_to_litter, turnover_daily, bm_sapl, tree, cn_ind,flux10,flux100, &
641!!$               prod10,prod100,prod10_total,prod100_total,&
642!!$               convflux,cflux_prod_total,cflux_prod10,cflux_prod100,leaf_frac,&
643               prod10,prod100,convflux,cflux_prod10,cflux_prod100,leaf_frac,&
644               npp_longterm, lm_lastyearmax, litter, carbon)
645
646       ENDIF
647    ENDIF
648!MM déplacement pour initialisation correcte des grandeurs cumulées :
649    cflux_prod_total(:) = convflux(:) + cflux_prod10(:) + cflux_prod100(:)
650    prod10_total(:)=SUM(prod10,dim=2)
651    prod100_total(:)=SUM(prod100,dim=2)
652    !
653    ! 16 total heterotrophic respiration
654
655    tot_soil_carb=zero
656    tot_litter_carb=zero
657    DO j=2,nvm
658
659       tot_litter_carb(:,j) = tot_litter_carb(:,j) + (litter(:,istructural,j,iabove) + &
660            &          litter(:,imetabolic,j,iabove) + &
661            &          litter(:,istructural,j,ibelow) + litter(:,imetabolic,j,ibelow))
662
663       tot_soil_carb(:,j) = tot_soil_carb(:,j) + (carbon(:,iactive,j) + &
664            &          carbon(:,islow,j)+  carbon(:,ipassive,j))
665
666    ENDDO
667    tot_litter_soil_carb = tot_litter_carb + tot_soil_carb
668
669    tot_live_biomass = biomass(:,:,ileaf) + biomass(:,:,isapabove) + biomass(:,:,isapbelow) +&
670         &             biomass(:,:,iheartabove) + biomass(:,:,iheartbelow) + &
671         &             biomass(:,:,iroot)+ biomass(:,:,ifruit)+ biomass(:,:,icarbres)
672
673    tot_turnover = turnover_daily(:,:,ileaf) + turnover_daily(:,:,isapabove) + &
674         &         turnover_daily(:,:,isapbelow) + turnover_daily(:,:,iheartabove) + &
675         &         turnover_daily(:,:,iheartbelow) + turnover_daily(:,:,iroot) + &
676         &         turnover_daily(:,:,ifruit) + turnover_daily(:,:,icarbres)
677
678    tot_bm_to_litter = bm_to_litter(:,:,ileaf) + bm_to_litter(:,:,isapabove) +&
679         &             bm_to_litter(:,:,isapbelow) + bm_to_litter(:,:,iheartbelow) +&
680         &             bm_to_litter(:,:,iheartabove) + bm_to_litter(:,:,iroot) + &
681         &             bm_to_litter(:,:,ifruit) + bm_to_litter(:,:,icarbres)
682
683    !
684    ! 17 history
685    !
686
687    ! 2d
688
689    CALL histwrite (hist_id_stomate, 'RESOLUTION_X', itime, &
690         resolution(:,1), npts, hori_index)
691    CALL histwrite (hist_id_stomate, 'RESOLUTION_Y', itime, &
692         resolution(:,2), npts, hori_index)
693    CALL histwrite (hist_id_stomate, 'CONTFRAC', itime, &
694         contfrac(:), npts, hori_index)
695
696    CALL histwrite (hist_id_stomate, 'LITTER_STR_AB', itime, &
697         litter(:,istructural,:,iabove), npts*nvm, horipft_index)
698    CALL histwrite (hist_id_stomate, 'LITTER_MET_AB', itime, &
699         litter(:,imetabolic,:,iabove), npts*nvm, horipft_index)
700    CALL histwrite (hist_id_stomate, 'LITTER_STR_BE', itime, &
701         litter(:,istructural,:,ibelow), npts*nvm, horipft_index)
702    CALL histwrite (hist_id_stomate, 'LITTER_MET_BE', itime, &
703         litter(:,imetabolic,:,ibelow), npts*nvm, horipft_index)
704
705    CALL histwrite (hist_id_stomate, 'DEADLEAF_COVER', itime, &
706         deadleaf_cover, npts, hori_index)
707
708    CALL histwrite (hist_id_stomate, 'TOTAL_SOIL_CARB', itime, &
709         tot_litter_soil_carb, npts*nvm, horipft_index)
710    CALL histwrite (hist_id_stomate, 'CARBON_ACTIVE', itime, &
711         carbon(:,iactive,:), npts*nvm, horipft_index)
712    CALL histwrite (hist_id_stomate, 'CARBON_SLOW', itime, &
713         carbon(:,islow,:), npts*nvm, horipft_index)
714    CALL histwrite (hist_id_stomate, 'CARBON_PASSIVE', itime, &
715         carbon(:,ipassive,:), npts*nvm, horipft_index)
716
717    CALL histwrite (hist_id_stomate, 'T2M_MONTH', itime, &
718         t2m_month, npts, hori_index)
719    CALL histwrite (hist_id_stomate, 'T2M_WEEK', itime, &
720         t2m_week, npts, hori_index)
721
722    CALL histwrite (hist_id_stomate, 'HET_RESP', itime, &
723         resp_hetero(:,:), npts*nvm, horipft_index)
724
725    CALL histwrite (hist_id_stomate, 'BLACK_CARBON', itime, &
726         black_carbon, npts, hori_index)
727
728    CALL histwrite (hist_id_stomate, 'FIREINDEX', itime, &
729         fireindex(:,:), npts*nvm, horipft_index)
730    CALL histwrite (hist_id_stomate, 'LITTERHUM', itime, &
731         litterhum_daily, npts, hori_index)
732    CALL histwrite (hist_id_stomate, 'CO2_FIRE', itime, &
733         co2_fire, npts*nvm, horipft_index)
734    CALL histwrite (hist_id_stomate, 'CO2_TAKEN', itime, &
735         co2_to_bm, npts*nvm, horipft_index)
736    ! land cover change
737    CALL histwrite (hist_id_stomate, 'CONVFLUX', itime, &
738         convflux, npts, hori_index)
739    CALL histwrite (hist_id_stomate, 'CFLUX_PROD10', itime, &
740         cflux_prod10, npts, hori_index)
741    CALL histwrite (hist_id_stomate, 'CFLUX_PROD100', itime, &
742         cflux_prod100, npts, hori_index)
743    CALL histwrite (hist_id_stomate, 'HARVEST_ABOVE', itime, &
744         harvest_above, npts, hori_index)
745
746    ! 3d
747
748    CALL histwrite (hist_id_stomate, 'LAI', itime, &
749         lai, npts*nvm, horipft_index)
750    CALL histwrite (hist_id_stomate, 'VEGET', itime, &
751         veget, npts*nvm, horipft_index)
752    CALL histwrite (hist_id_stomate, 'VEGET_MAX', itime, &
753         veget_max, npts*nvm, horipft_index)
754    CALL histwrite (hist_id_stomate, 'NPP', itime, &
755         npp_daily, npts*nvm, horipft_index)
756    CALL histwrite (hist_id_stomate, 'GPP', itime, &
757         gpp_daily, npts*nvm, horipft_index)
758    CALL histwrite (hist_id_stomate, 'IND', itime, &
759         ind, npts*nvm, horipft_index)
760    CALL histwrite (hist_id_stomate, 'TOTAL_M', itime, &
761         tot_live_biomass, npts*nvm, horipft_index)
762    CALL histwrite (hist_id_stomate, 'LEAF_M', itime, &
763         biomass(:,:,ileaf), npts*nvm, horipft_index)
764    CALL histwrite (hist_id_stomate, 'SAP_M_AB', itime, &
765         biomass(:,:,isapabove), npts*nvm, horipft_index)
766    CALL histwrite (hist_id_stomate, 'SAP_M_BE', itime, &
767         biomass(:,:,isapbelow), npts*nvm, horipft_index)
768    CALL histwrite (hist_id_stomate, 'HEART_M_AB', itime, &
769         biomass(:,:,iheartabove), npts*nvm, horipft_index)
770    CALL histwrite (hist_id_stomate, 'HEART_M_BE', itime, &
771         biomass(:,:,iheartbelow), npts*nvm, horipft_index)
772    CALL histwrite (hist_id_stomate, 'ROOT_M', itime, &
773         biomass(:,:,iroot), npts*nvm, horipft_index)
774    CALL histwrite (hist_id_stomate, 'FRUIT_M', itime, &
775         biomass(:,:,ifruit), npts*nvm, horipft_index)
776    CALL histwrite (hist_id_stomate, 'RESERVE_M', itime, &
777         biomass(:,:,icarbres), npts*nvm, horipft_index)
778    CALL histwrite (hist_id_stomate, 'TOTAL_TURN', itime, &
779         tot_turnover, npts*nvm, horipft_index)
780    CALL histwrite (hist_id_stomate, 'LEAF_TURN', itime, &
781         turnover_daily(:,:,ileaf), npts*nvm, horipft_index)
782    CALL histwrite (hist_id_stomate, 'SAP_AB_TURN', itime, &
783         turnover_daily(:,:,isapabove), npts*nvm, horipft_index)
784    CALL histwrite (hist_id_stomate, 'ROOT_TURN', itime, &
785         turnover_daily(:,:,iroot), npts*nvm, horipft_index)
786    CALL histwrite (hist_id_stomate, 'FRUIT_TURN', itime, &
787         turnover_daily(:,:,ifruit), npts*nvm, horipft_index)
788    CALL histwrite (hist_id_stomate, 'TOTAL_BM_LITTER', itime, &
789         tot_bm_to_litter, npts*nvm, horipft_index)
790    CALL histwrite (hist_id_stomate, 'LEAF_BM_LITTER', itime, &
791         bm_to_litter(:,:,ileaf), npts*nvm, horipft_index)
792    CALL histwrite (hist_id_stomate, 'SAP_AB_BM_LITTER', itime, &
793         bm_to_litter(:,:,isapabove), npts*nvm, horipft_index)
794    CALL histwrite (hist_id_stomate, 'SAP_BE_BM_LITTER', itime, &
795         bm_to_litter(:,:,isapbelow), npts*nvm, horipft_index)
796    CALL histwrite (hist_id_stomate, 'HEART_AB_BM_LITTER', itime, &
797         bm_to_litter(:,:,iheartabove), npts*nvm, horipft_index)
798    CALL histwrite (hist_id_stomate, 'HEART_BE_BM_LITTER', itime, &
799         bm_to_litter(:,:,iheartbelow), npts*nvm, horipft_index)
800    CALL histwrite (hist_id_stomate, 'ROOT_BM_LITTER', itime, &
801         bm_to_litter(:,:,iroot), npts*nvm, horipft_index)
802    CALL histwrite (hist_id_stomate, 'FRUIT_BM_LITTER', itime, &
803         bm_to_litter(:,:,ifruit), npts*nvm, horipft_index)
804    CALL histwrite (hist_id_stomate, 'RESERVE_BM_LITTER', itime, &
805         bm_to_litter(:,:,icarbres), npts*nvm, horipft_index)
806    CALL histwrite (hist_id_stomate, 'MAINT_RESP', itime, &
807         resp_maint, npts*nvm, horipft_index)
808    CALL histwrite (hist_id_stomate, 'GROWTH_RESP', itime, &
809         resp_growth, npts*nvm, horipft_index)
810    CALL histwrite (hist_id_stomate, 'AGE', itime, &
811         age, npts*nvm, horipft_index)
812    CALL histwrite (hist_id_stomate, 'HEIGHT', itime, &
813         height, npts*nvm, horipft_index)
814    CALL histwrite (hist_id_stomate, 'MOISTRESS', itime, &
815         moiavail_week, npts*nvm, horipft_index)
816    CALL histwrite (hist_id_stomate, 'VCMAX', itime, &
817         vcmax, npts*nvm, horipft_index)
818    CALL histwrite (hist_id_stomate, 'TURNOVER_TIME', itime, &
819         turnover_time, npts*nvm, horipft_index)
820    ! land cover change
821    CALL histwrite (hist_id_stomate, 'PROD10', itime, &
822         prod10, npts*11, horip11_index)
823    CALL histwrite (hist_id_stomate, 'PROD100', itime, &
824         prod100, npts*101, horip101_index)
825    CALL histwrite (hist_id_stomate, 'FLUX10', itime, &
826         flux10, npts*10, horip10_index)
827    CALL histwrite (hist_id_stomate, 'FLUX100', itime, &
828         flux100, npts*100, horip100_index)
829
830    IF ( hist_id_stomate_IPCC > 0 ) THEN
831       vartmp(:)=SUM(tot_live_biomass*veget_max,dim=2)/1e3*contfrac
832       CALL histwrite (hist_id_stomate_IPCC, "cVeg", itime, &
833         vartmp, npts, hori_index)
834       vartmp(:)=SUM(tot_litter_carb*veget_max,dim=2)/1e3*contfrac
835       CALL histwrite (hist_id_stomate_IPCC, "cLitter", itime, &
836         vartmp, npts, hori_index)
837       vartmp(:)=SUM(tot_soil_carb*veget_max,dim=2)/1e3*contfrac
838       CALL histwrite (hist_id_stomate_IPCC, "cSoil", itime, &
839         vartmp, npts, hori_index)
840       vartmp(:)=(prod10_total + prod100_total)/1e3
841       CALL histwrite (hist_id_stomate_IPCC, "cProduct", itime, &
842         vartmp, npts, hori_index)
843       vartmp(:)=SUM(lai*veget_max,dim=2)*contfrac
844       CALL histwrite (hist_id_stomate_IPCC, "lai", itime, &
845         vartmp, npts, hori_index)
846       vartmp(:)=SUM(gpp_daily*veget_max,dim=2)/1e3/one_day*contfrac
847       CALL histwrite (hist_id_stomate_IPCC, "gpp", itime, &
848         vartmp, npts, hori_index)
849       vartmp(:)=SUM((resp_maint+resp_growth)*veget_max,dim=2)/1e3/one_day*contfrac
850       CALL histwrite (hist_id_stomate_IPCC, "ra", itime, &
851         vartmp, npts, hori_index)
852       vartmp(:)=SUM(npp_daily*veget_max,dim=2)/1e3/one_day*contfrac
853       CALL histwrite (hist_id_stomate_IPCC, "npp", itime, &
854         vartmp, npts, hori_index)
855       vartmp(:)=SUM(resp_hetero*veget_max,dim=2)/1e3/one_day*contfrac
856       CALL histwrite (hist_id_stomate_IPCC, "rh", itime, &
857         vartmp, npts, hori_index)
858       vartmp(:)=SUM(co2_fire*veget_max,dim=2)/1e3/one_day*contfrac
859       CALL histwrite (hist_id_stomate_IPCC, "fFire", itime, &
860         vartmp, npts, hori_index)
861       vartmp(:)=harvest_above/1e3/one_day*contfrac
862       CALL histwrite (hist_id_stomate_IPCC, "fHarvest", itime, &
863         vartmp, npts, hori_index)
864       vartmp(:)=cflux_prod_total/1e3/one_day*contfrac
865       CALL histwrite (hist_id_stomate_IPCC, "fLuc", itime, &
866         vartmp, npts, hori_index)
867       vartmp(:)=(SUM((gpp_daily-(resp_maint+resp_growth+resp_hetero)-co2_fire) &
868            &        *veget_max,dim=2)-cflux_prod_total-harvest_above)/1e3/one_day*contfrac
869       CALL histwrite (hist_id_stomate_IPCC, "nbp", itime, &
870         vartmp, npts, hori_index)
871       vartmp(:)=SUM(tot_bm_to_litter*veget_max,dim=2)/1e3/one_day*contfrac
872       CALL histwrite (hist_id_stomate_IPCC, "fVegLitter", itime, &
873         vartmp, npts, hori_index)
874       vartmp(:)=SUM(SUM(soilcarbon_input,dim=2)*veget_max,dim=2)/1e3/one_day*contfrac
875       CALL histwrite (hist_id_stomate_IPCC, "fLitterSoil", itime, &
876         vartmp, npts, hori_index)
877       vartmp(:)=SUM(biomass(:,:,ileaf)*veget_max,dim=2)/1e3*contfrac
878       CALL histwrite (hist_id_stomate_IPCC, "cLeaf", itime, &
879         vartmp, npts, hori_index)
880       vartmp(:)=SUM((biomass(:,:,isapabove)+biomass(:,:,iheartabove))*veget_max,dim=2)/1e3*contfrac
881       CALL histwrite (hist_id_stomate_IPCC, "cWood", itime, &
882         vartmp, npts, hori_index)
883       vartmp(:)=SUM(( biomass(:,:,iroot) + biomass(:,:,isapbelow) + biomass(:,:,iheartbelow) ) &
884            &        *veget_max,dim=2)/1e3*contfrac
885       CALL histwrite (hist_id_stomate_IPCC, "cRoot", itime, &
886         vartmp, npts, hori_index)
887       vartmp(:)=SUM(( biomass(:,:,icarbres) + biomass(:,:,ifruit))*veget_max,dim=2)/1e3*contfrac
888       CALL histwrite (hist_id_stomate_IPCC, "cMisc", itime, &
889         vartmp, npts, hori_index)
890       vartmp(:)=SUM((litter(:,istructural,:,iabove)+litter(:,imetabolic,:,iabove))*veget_max,dim=2)/1e3*contfrac
891       CALL histwrite (hist_id_stomate_IPCC, "cLitterAbove", itime, &
892         vartmp, npts, hori_index)
893       vartmp(:)=SUM((litter(:,istructural,:,ibelow)+litter(:,imetabolic,:,ibelow))*veget_max,dim=2)/1e3*contfrac
894       CALL histwrite (hist_id_stomate_IPCC, "cLitterBelow", itime, &
895         vartmp, npts, hori_index)
896       vartmp(:)=SUM(carbon(:,iactive,:)*veget_max,dim=2)/1e3*contfrac
897       CALL histwrite (hist_id_stomate_IPCC, "cSoilFast", itime, &
898         vartmp, npts, hori_index)
899       vartmp(:)=SUM(carbon(:,islow,:)*veget_max,dim=2)/1e3*contfrac
900       CALL histwrite (hist_id_stomate_IPCC, "cSoilMedium", itime, &
901         vartmp, npts, hori_index)
902       vartmp(:)=SUM(carbon(:,ipassive,:)*veget_max,dim=2)/1e3*contfrac
903       CALL histwrite (hist_id_stomate_IPCC, "cSoilSlow", itime, &
904         vartmp, npts, hori_index)
905       DO j=1,nvm
906          histvar(:,j)=veget_max(:,j)*contfrac(:)*100
907       ENDDO
908       CALL histwrite (hist_id_stomate_IPCC, "landCoverFrac", itime, &
909         histvar, npts*nvm, horipft_index)
910       vartmp(:)=(veget_max(:,3)+veget_max(:,6)+veget_max(:,8)+veget_max(:,9))*contfrac*100
911       CALL histwrite (hist_id_stomate_IPCC, "treeFracPrimDec", itime, &
912          vartmp, npts, hori_index)
913       vartmp(:)=(veget_max(:,2)+veget_max(:,4)+veget_max(:,5)+veget_max(:,7))*contfrac*100
914       CALL histwrite (hist_id_stomate_IPCC, "treeFracPrimEver", itime, &
915         vartmp, npts, hori_index)
916       vartmp(:)=(veget_max(:,10)+veget_max(:,12))*contfrac*100
917       CALL histwrite (hist_id_stomate_IPCC, "c3PftFrac", itime, &
918         vartmp, npts, hori_index)
919       vartmp(:)=(veget_max(:,11)+veget_max(:,13))*contfrac*100
920       CALL histwrite (hist_id_stomate_IPCC, "c4PftFrac", itime, &
921         vartmp, npts, hori_index)
922       vartmp(:)=SUM(resp_growth*veget_max,dim=2)/1e3/one_day*contfrac
923       CALL histwrite (hist_id_stomate_IPCC, "rGrowth", itime, &
924         vartmp, npts, hori_index)
925       vartmp(:)=SUM(resp_maint*veget_max,dim=2)/1e3/one_day*contfrac
926       CALL histwrite (hist_id_stomate_IPCC, "rMaint", itime, &
927         vartmp, npts, hori_index)
928       vartmp(:)=SUM(bm_alloc(:,:,ileaf)*veget_max,dim=2)/1e3/one_day*contfrac
929       CALL histwrite (hist_id_stomate_IPCC, "nppLeaf", itime, &
930         vartmp, npts, hori_index)
931       vartmp(:)=SUM(bm_alloc(:,:,isapabove)*veget_max,dim=2)/1e3/one_day*contfrac
932       CALL histwrite (hist_id_stomate_IPCC, "nppWood", itime, &
933         vartmp, npts, hori_index)
934       vartmp(:)=SUM(( bm_alloc(:,:,isapbelow) + bm_alloc(:,:,iroot) )*veget_max,dim=2)/1e3/one_day*contfrac
935       CALL histwrite (hist_id_stomate_IPCC, "nppRoot", itime, &
936         vartmp, npts, hori_index)
937
938       CALL histwrite (hist_id_stomate_IPCC, 'RESOLUTION_X', itime, &
939            resolution(:,1), npts, hori_index)
940       CALL histwrite (hist_id_stomate_IPCC, 'RESOLUTION_Y', itime, &
941            resolution(:,2), npts, hori_index)
942       CALL histwrite (hist_id_stomate_IPCC, 'CONTFRAC', itime, &
943            contfrac(:), npts, hori_index)
944
945    ENDIF
946
947    IF (bavard.GE.4) WRITE(numout,*) 'Leaving stomate_lpj'
948
949  END SUBROUTINE StomateLpj
950
951  SUBROUTINE harvest(npts, dt_days, veget_max, veget, &
952       bm_to_litter, turnover_daily, &
953       harvest_above)
954    ! 0.1 input
955
956    ! Domain size
957    INTEGER, INTENT(in)                                            :: npts
958
959    ! Time step (days)
960    REAL(r_std), INTENT(in)                                         :: dt_days
961
962    ! new "maximal" coverage fraction of a PFT (LAI -> infinity) on ground
963    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)                 :: veget_max
964    ! 0.2 modified fields
965
966    ! fractional coverage on natural/agricultural ground, taking into
967    !   account LAI (=grid-scale fpc)
968    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)                    :: veget
969
970    ! conversion of biomass to litter (gC/(m**2 of nat/agri ground)) / day
971    REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(inout)           :: bm_to_litter
972
973    ! Turnover rates (gC/(m**2 of ground)/day)
974    REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(inout)          :: turnover_daily
975    ! harvest above ground biomass for agriculture
976    REAL(r_std), DIMENSION(npts), INTENT(inout)                       :: harvest_above
977    ! 0.4 local
978
979    ! indices
980    INTEGER(i_std)                                                     :: i, j, k, l, m
981
982    ! biomass increase (gC/(m**2 of ground))
983    REAL(r_std)                                                        :: above_old
984    ! yearly initialisation
985    above_old             = zero
986    harvest_above         = zero
987    DO i = 1, npts
988       DO j=1, nvm
989          IF (j > 11) THEN
990             above_old = turnover_daily(i,j,ileaf) + turnover_daily(i,j,isapabove) + &
991                  &       turnover_daily(i,j,iheartabove) + turnover_daily(i,j,ifruit) + &
992                  &       turnover_daily(i,j,icarbres) + turnover_daily(i,j,isapbelow) + &
993                  &       turnover_daily(i,j,iheartbelow) + turnover_daily(i,j,iroot)
994
995             turnover_daily(i,j,ileaf) = turnover_daily(i,j,ileaf)*0.55
996             turnover_daily(i,j,isapabove) = turnover_daily(i,j,isapabove)*0.55
997             turnover_daily(i,j,isapbelow) = turnover_daily(i,j,isapbelow)*0.55
998             turnover_daily(i,j,iheartabove) = turnover_daily(i,j,iheartabove)*0.55
999             turnover_daily(i,j,iheartbelow) = turnover_daily(i,j,iheartbelow)*0.55
1000             turnover_daily(i,j,iroot) = turnover_daily(i,j,iroot)*0.55
1001             turnover_daily(i,j,ifruit) = turnover_daily(i,j,ifruit)*0.55
1002             turnover_daily(i,j,icarbres) = turnover_daily(i,j,icarbres)*0.55
1003             harvest_above(i)  = harvest_above(i) + veget_max(i,j) * above_old *0.45
1004
1005          ENDIF
1006       ENDDO
1007    ENDDO
1008
1009!!$    harvest_above = harvest_above
1010  END SUBROUTINE harvest
1011END MODULE stomate_lpj
Note: See TracBrowser for help on using the repository browser.