source: branches/ORCHIDEE_EXT/ORCHIDEE/src_stomate/stomate_lpj.f90 @ 335

Last change on this file since 335 was 335, checked in by didier.solyga, 13 years ago

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