source: branches/ORCHIDEE_EXT/ORCHIDEE/src_stomate/stomate_npp.f90 @ 257

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

Externalized version merged with the trunk
File size: 15.8 KB
Line 
1! Npp: Maintenance and growth respiration
2! We calculte first the maintenance rspiration. This is substracted from the
3! allocatable biomass (and from the present biomass if the GPP is too low).
4! Of the rest, a part is lost as growth respiration, while the other part is
5! effectively allocated.
6!
7! $Header: /home/ssipsl/CVSREP/ORCHIDEE/src_stomate/stomate_npp.f90,v 1.14 2010/04/20 14:12:04 ssipsl Exp $
8! IPSL (2006)
9!  This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
10!
11MODULE stomate_npp
12
13  ! modules used:
14
15  USE ioipsl
16  USE stomate_data
17  USE constantes
18  USE pft_parameters
19
20  IMPLICIT NONE
21
22  ! private & public routines
23
24  PRIVATE
25  PUBLIC npp_calc,npp_calc_clear
26
27  ! first call
28  LOGICAL, SAVE                                              :: firstcall = .TRUE.
29
30CONTAINS
31
32  SUBROUTINE npp_calc_clear
33    firstcall=.TRUE.
34  END SUBROUTINE npp_calc_clear
35
36  SUBROUTINE npp_calc (npts, dt, &
37       PFTpresent, &
38       tlong_ref, t2m, tsoil, lai, rprof, &
39       gpp, f_alloc, bm_alloc, resp_maint_part,&
40       biomass, leaf_age, leaf_frac, age, &
41       resp_maint, resp_growth, npp)
42
43    !
44    ! 0 declarations
45    !
46
47    ! 0.1 input
48
49    ! Domain size
50    INTEGER(i_std), INTENT(in)                                        :: npts
51    ! time step (days)
52    REAL(r_std), INTENT(in)                                     :: dt
53    ! PFT exists
54    LOGICAL, DIMENSION(npts,nvm), INTENT(in)                  :: PFTpresent
55    ! long term annual mean 2 meter reference temperature
56    REAL(r_std), DIMENSION(npts), INTENT(in)                    :: tlong_ref
57    ! 2 meter temperature
58    REAL(r_std), DIMENSION(npts), INTENT(in)                    :: t2m
59    ! soil temperature (K)
60    REAL(r_std), DIMENSION(npts,nbdl), INTENT(in)               :: tsoil
61    ! leaf area index
62    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: lai
63    ! root depth (m)
64    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: rprof
65    ! gross primary productivity (gC/days/(m**2 of ground))
66    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)               :: gpp
67    ! fraction that goes into plant part
68    REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(in)        :: f_alloc
69    ! maintenance respiration of different plant parts (gC/day/m**2 of ground)
70    REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(in)             :: resp_maint_part
71
72    ! 0.2 modified fields
73
74    ! biomass (gC/(m**2 of ground))
75    REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(inout)     :: biomass
76    ! leaf age (days)
77    REAL(r_std), DIMENSION(npts,nvm,nleafages), INTENT(inout)  :: leaf_age
78    ! fraction of leaves in leaf age class
79    REAL(r_std), DIMENSION(npts,nvm,nleafages), INTENT(inout)  :: leaf_frac
80    ! age (years)
81    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)            :: age
82
83    ! 0.3 output
84
85    ! maintenance respiration (gC/day/m**2 of total ground)
86    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)              :: resp_maint
87    ! autotrophic respiration (gC/day/m**2 of total ground)
88    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)              :: resp_growth
89    ! net primary productivity (gC/day/m**2 of ground)
90    REAL(r_std), DIMENSION(npts,nvm), INTENT(out)              :: npp
91    ! biomass increase, i.e. NPP per plant part
92    REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(out)       :: bm_alloc
93
94    ! 0.4 local
95
96    ! soil levels (m)
97    REAL(r_std), SAVE, DIMENSION(0:nbdl)                        :: z_soil
98    ! root temperature (convolution of root and soil temperature profiles)
99    REAL(r_std), DIMENSION(npts,nvm)                           :: t_root
100    ! maintenance respiration coefficients at 0 deg C (g/g d**-1)
101    REAL(r_std), DIMENSION(npts,nvm,nparts)                    :: coeff_maint
102    ! temperature which is pertinent for maintenance respiration (K)
103    REAL(r_std), DIMENSION(npts,nparts)                         :: t_maint
104    ! integration constant for root profile
105    REAL(r_std), DIMENSION(npts)                                :: rpc
106    ! long term annual mean temperature, C
107    REAL(r_std), DIMENSION(npts)                                :: tl
108    ! slope of maintenance respiration coefficient (1/K)
109    REAL(r_std), DIMENSION(npts)                                :: slope
110    ! growth respiration of different plant parts (gC/day/m**2 of ground)
111    REAL(r_std), DIMENSION(npts,nparts)                         :: resp_growth_part
112    ! allocatable biomass (gC/m**2 of ground) for the whole plant
113    REAL(r_std), DIMENSION(npts,nvm)                           :: bm_alloc_tot
114    ! biomass increase
115    REAL(r_std), DIMENSION(npts)                                :: bm_add
116    ! new biomass
117    REAL(r_std), DIMENSION(npts)                                :: bm_new
118    ! leaf mass in youngest age class (gC/m**2 of ground)
119    REAL(r_std), DIMENSION(npts,nvm)                           :: leaf_mass_young
120    ! leaf mass after maintenance respiration
121    REAL(r_std), DIMENSION(npts,nvm)                           :: lm_old
122    ! biomass created when biomass<0 because of dark respiration (gC/m**2 of ground)
123    REAL(r_std), DIMENSION(npts,nvm)                           :: bm_create
124    ! maximum part of allocatable biomass used for respiration
125    REAL(r_std), DIMENSION(npts)                                :: bm_tax_max
126    ! biomass that remains to be taken away
127    REAL(r_std), DIMENSION(npts)                                :: bm_pump
128    ! Index
129    INTEGER(i_std)                                             :: i,j,k,l,m
130
131    ! =========================================================================
132
133    IF (bavard.GE.3) WRITE(numout,*) 'Entering npp'
134    !
135    ! 1 Initializations
136    !
137
138    !
139    ! 1.1 first call
140    !
141
142    IF ( firstcall ) THEN
143
144       ! 1.1.1 soil levels
145
146       z_soil(0) = zero
147       z_soil(1:nbdl) = diaglev(1:nbdl)
148
149       ! 1.1.2 messages
150
151       WRITE(numout,*) 'npp:'
152
153       WRITE(numout,*) '   > max. fraction of allocatable biomass used for'// &
154            ' maint. resp.:', tax_max
155
156       firstcall = .FALSE.
157
158    ENDIF
159
160    !
161    ! 1.2 set output to zero
162    !
163
164    bm_alloc(:,:,:) = zero
165    resp_maint(:,:) = zero
166    resp_growth(:,:) = zero
167    npp(:,:) = zero
168
169    !
170    ! 1.3 root temperature: convolution of root and temperature profiles
171    !     suppose exponential root profile.
172    !
173
174    DO j = 2,nvm
175
176       ! 1.3.1 rpc is an integration constant such that the integral of the root profile is 1.
177
178       rpc(:) = un / ( un - EXP( -z_soil(nbdl) / rprof(:,j) ) )
179
180       ! 1.3.2 integrate over the nbdl levels
181
182       t_root(:,j) = zero
183
184       DO l = 1, nbdl
185
186          t_root(:,j) = &
187               t_root(:,j) + tsoil(:,l) * rpc(:) * &
188               ( EXP( -z_soil(l-1)/rprof(:,j) ) - EXP( -z_soil(l)/rprof(:,j) ) )
189
190       ENDDO
191
192    ENDDO
193
194    !
195    ! 1.4 total allocatable biomass
196    !
197
198    bm_alloc_tot(:,:) = gpp(:,:) * dt
199
200    !
201    ! 2 define maintenance respiration coefficients
202    !
203
204    DO j = 2,nvm
205
206       !
207       ! 2.1 temperature which is taken for the plant part we are talking about
208       !
209
210       ! 2.1.1 parts above the ground
211
212       t_maint(:,ileaf) = t2m(:)
213       t_maint(:,isapabove) = t2m(:)
214       t_maint(:,ifruit) = t2m(:)
215
216       ! 2.1.2 parts below the ground
217
218       t_maint(:,isapbelow) = t_root(:,j)
219       t_maint(:,iroot) = t_root(:,j)
220
221       ! 2.1.3 heartwood: does not respire. Any temperature
222
223       t_maint(:,iheartbelow) = t_root(:,j)
224       t_maint(:,iheartabove) = t2m(:)
225
226       ! 2.1.4 reserve: above the ground for trees, below for grasses
227
228       IF ( tree(j) ) THEN
229          t_maint(:,icarbres) = t2m(:)
230       ELSE
231          t_maint(:,icarbres) = t_root(:,j)
232       ENDIF
233
234       !
235       ! 2.2 calculate coefficient
236       !
237
238       tl(:) = tlong_ref(:) - ZeroCelsius
239       slope(:) = maint_resp_slope(j,1) + tl(:) * maint_resp_slope(j,2) + &
240            tl(:)*tl(:) * maint_resp_slope(j,3)
241
242       DO k = 1, nparts
243
244          coeff_maint(:,j,k) = &
245               MAX( coeff_maint_zero(j,k) * &
246               ( un + slope(:) * (t_maint(:,k)-ZeroCelsius) ), zero )
247
248       ENDDO
249
250    ENDDO
251
252    !
253    ! 3 calculate maintenance and growth respiration.
254    !   NPP = GPP - maintenance resp - growth resp.
255    !
256    !
257    DO j = 2,nvm
258
259       !
260       ! 3.1 maintenance respiration of the different plant parts
261       !
262       !
263       ! 3.2 Total maintenance respiration of the plant
264       !     VPP killer:
265       !     resp_maint(:,j) = SUM( resp_maint_part(:,:), DIM=2 )
266       !
267
268       resp_maint(:,j) = zero
269
270       ! with the new calculation of hourly respiration, we must verify that
271       ! PFT has not been killed after calcul of resp_maint_part in stomate
272       DO k= 1, nparts
273          WHERE (PFTpresent(:,j))
274             resp_maint(:,j) = resp_maint(:,j) + resp_maint_part(:,j,k)
275          ENDWHERE
276       ENDDO
277       !
278       ! 3.3 This maintenance respiration is taken away from the newly produced
279       !     allocatable biomass. However, we avoid that no allocatable biomass remains.
280       !     If the respiration is larger than a given fraction of the allocatable biomass,
281       !     the rest is taken from the tissues themselves.
282       !     We suppose that respiration is not dependent on leaf age ->
283       !     do not change age structure.
284       !
285
286       ! maximum part of allocatable biomass used for respiration
287       bm_tax_max(:) = tax_max * bm_alloc_tot(:,j)
288
289       DO i = 1, npts
290
291          IF ( ( bm_alloc_tot(i,j) .GT. zero ) .AND. &
292               ( ( resp_maint(i,j) * dt ) .LT. bm_tax_max(i) ) ) THEN
293
294             bm_alloc_tot(i,j) = bm_alloc_tot(i,j) - resp_maint(i,j) * dt
295
296             !Shilong        ELSEIF ( resp_maint(i,j) .GT. zero ) THEN
297          ELSEIF ( resp_maint(i,j) .GT. min_stomate ) THEN
298
299             ! remaining allocatable biomass
300             bm_alloc_tot(i,j) = bm_alloc_tot(i,j) - bm_tax_max(i)
301
302             ! biomass that remains to be taken away from tissues
303             bm_pump(i) = resp_maint(i,j) * dt - bm_tax_max(i)
304
305             ! take biomass from tissues
306
307             biomass(i,j,ileaf) = biomass(i,j,ileaf) - &
308                  bm_pump(i) * resp_maint_part(i,j,ileaf) / resp_maint(i,j)
309             biomass(i,j,isapabove) = biomass(i,j,isapabove) - &
310                  bm_pump(i) * resp_maint_part(i,j,isapabove) / resp_maint(i,j)
311             biomass(i,j,isapbelow) = biomass(i,j,isapbelow) - &
312                  bm_pump(i) * resp_maint_part(i,j,isapbelow) / resp_maint(i,j)
313             biomass(i,j,iroot) = biomass(i,j,iroot) - &
314                  bm_pump(i) * resp_maint_part(i,j,iroot) / resp_maint(i,j)
315             biomass(i,j,ifruit) = biomass(i,j,ifruit) - &
316                  bm_pump(i) * resp_maint_part(i,j,ifruit) / resp_maint(i,j)
317             biomass(i,j,icarbres) = biomass(i,j,icarbres) - &
318                  bm_pump(i) * resp_maint_part(i,j,icarbres) / resp_maint(i,j)
319
320          ENDIF
321
322       ENDDO   ! Fortran95: WHERE - ELSEWHERE construct
323
324       !
325       ! 3.4 dispatch allocatable biomass
326       !
327
328       DO k = 1, nparts
329          bm_alloc(:,j,k) = f_alloc(:,j,k) * bm_alloc_tot(:,j)
330       ENDDO
331
332       !
333       ! 3.5 growth respiration of a plant part is a given fraction of the
334       !     remaining allocatable biomass.
335       !
336
337       resp_growth_part(:,:) = frac_growthresp * bm_alloc(:,j,:) / dt
338
339       bm_alloc(:,j,:) = ( un - frac_growthresp ) * bm_alloc(:,j,:)
340
341       !
342       ! 3.6 Total growth respiration of the plant
343       !     VPP killer:
344       !     resp_growth(:,j) = SUM( resp_growth_part(:,:), DIM=2 )
345       !
346
347       resp_growth(:,j) = zero
348
349       DO k = 1, nparts
350          resp_growth(:,j) = resp_growth(:,j) + resp_growth_part(:,k)
351       ENDDO
352
353    ENDDO
354
355    !
356    ! 4 update the biomass, but save the old leaf mass for later
357    !   "old" leaf mass is leaf mass after maintenance respiration
358    !
359
360    lm_old(:,:) = biomass(:,:,ileaf)
361
362    biomass(:,:,:) = biomass(:,:,:) + bm_alloc(:,:,:)
363
364    !
365    ! 5 biomass can become negative in some rare cases, as the GPP can be negative
366    !   (dark respiration).
367    !   In this case, set biomass to some small value. This creation of matter is taken into
368    !   account by decreasing the autotrophic respiration. In this case, maintenance respiration
369    !   can become negative !!!
370    !
371
372    DO k = 1, nparts
373
374       DO j = 2,nvm
375
376          WHERE ( biomass(:,j,k) .LT. zero )
377
378             bm_create(:,j) = min_stomate - biomass(:,j,k)
379             
380             biomass(:,j,k) = biomass(:,j,k) + bm_create(:,j)
381             
382             resp_maint(:,j) = resp_maint(:,j) - bm_create(:,j) / dt
383
384          ENDWHERE
385
386       ENDDO
387
388    ENDDO
389
390    !
391    ! 6 Calculate the NPP (gC/(m**2 of ground/day)
392    !
393
394    DO j = 2,nvm
395       npp(:,j) = gpp(:,j) - resp_growth(:,j) - resp_maint(:,j)
396    ENDDO
397
398    !
399    ! 7 leaf age
400    !
401
402    !
403    ! 7.1 Decrease leaf age in youngest class if new leaf biomass is higher than old one.
404    !
405
406    DO j = 2,nvm
407       leaf_mass_young(:,j) = leaf_frac(:,j,1) * lm_old(:,j) + bm_alloc(:,j,ileaf)
408    ENDDO
409
410    DO j = 2,nvm
411       WHERE ( ( bm_alloc(:,j,ileaf) .GT. zero ) .AND. &
412         ( leaf_mass_young(:,j) .GT. zero ) )
413
414          leaf_age(:,j,1) = MAX ( zero, &
415               & leaf_age(:,j,1) * &
416               & ( leaf_mass_young(:,j) - bm_alloc(:,j,ileaf) ) / &
417               & leaf_mass_young(:,j) )
418         
419       ENDWHERE
420    ENDDO
421
422    !
423    ! 7.2 new age class fractions (fraction in youngest class increases)
424    !
425
426    ! 7.2.1 youngest class: new mass in youngest class divided by total new mass
427
428    DO j = 2,nvm
429       WHERE ( biomass(:,j,ileaf) .GT. min_stomate )
430
431          leaf_frac(:,j,1) = leaf_mass_young(:,j) / biomass(:,j,ileaf)
432
433       ENDWHERE
434    ENDDO
435
436    ! 7.2.2 other classes: old mass in leaf age class divided by new mass
437
438    DO m = 2, nleafages
439
440       DO j = 2,nvm
441          WHERE ( biomass(:,j,ileaf) .GT. min_stomate )
442
443             leaf_frac(:,j,m) = leaf_frac(:,j,m) * lm_old(:,j) / biomass(:,j,ileaf)
444
445          ENDWHERE
446       ENDDO
447
448    ENDDO
449
450    !
451    ! 8 Plant age (years)
452    !
453
454    !
455    ! 8.1 Increase age at every time step
456    !
457
458    WHERE ( PFTpresent(:,:) )
459
460       age(:,:) = age(:,:) + dt/one_year
461
462    ELSEWHERE
463
464       age(:,:) = zero
465
466    ENDWHERE
467
468    !
469    ! 8.2 For grasses, decrease age
470    !     if new biomass is higher than old one.
471    !     For trees, age is treated in "establish" if vegetation is dynamic,
472    !     and in turnover routines if it is static (in this case, only take
473    !     into account the age of the heartwood).
474    !
475
476    DO j = 2,nvm
477
478       IF ( .NOT. tree(j) ) THEN
479
480          ! Only four compartments for grasses
481          ! VPP killer:
482          ! bm_new(:) = SUM( biomass(:,j,:), DIM=2 )
483          ! bm_add(:) = SUM( bm_alloc(:,j,:), DIM=2 )
484
485          bm_new(:) = biomass(:,j,ileaf) + biomass(:,j,isapabove) + &
486               biomass(:,j,iroot) + biomass(:,j,ifruit)
487          bm_add(:) = bm_alloc(:,j,ileaf) + bm_alloc(:,j,isapabove) + &
488               bm_alloc(:,j,iroot) + bm_alloc(:,j,ifruit)
489
490          WHERE ( ( bm_new(:) .GT. zero ) .AND. ( bm_add(:) .GT. zero ) )
491             age(:,j) = age(:,j) * ( bm_new(:) - bm_add(:) ) / bm_new(:)
492          ENDWHERE
493
494       ENDIF
495
496    ENDDO
497
498    !
499    ! 9 history
500    !
501
502    CALL histwrite (hist_id_stomate, 'BM_ALLOC_LEAF', itime, &
503         bm_alloc(:,:,ileaf), npts*nvm, horipft_index)
504    CALL histwrite (hist_id_stomate, 'BM_ALLOC_SAP_AB', itime, &
505         bm_alloc(:,:,isapabove), npts*nvm, horipft_index)
506    CALL histwrite (hist_id_stomate, 'BM_ALLOC_SAP_BE', itime, &
507         bm_alloc(:,:,isapbelow), npts*nvm, horipft_index)
508    CALL histwrite (hist_id_stomate, 'BM_ALLOC_ROOT', itime, &
509         bm_alloc(:,:,iroot), npts*nvm, horipft_index)
510    CALL histwrite (hist_id_stomate, 'BM_ALLOC_FRUIT', itime, &
511         bm_alloc(:,:,ifruit), npts*nvm, horipft_index)
512    CALL histwrite (hist_id_stomate, 'BM_ALLOC_RES', itime, &
513         bm_alloc(:,:,icarbres), npts*nvm, horipft_index)
514
515    IF (bavard.GE.4) WRITE(numout,*) 'Leaving npp'
516
517  END SUBROUTINE npp_calc
518
519END MODULE stomate_npp
Note: See TracBrowser for help on using the repository browser.