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

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

Import first version of ORCHIDEE_EXT
File size: 15.7 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) = 0.
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       rpc(:) = 1. / ( 1. - EXP( -z_soil(nbdl) / rprof(:,j) ) )
178
179       ! 1.3.2 integrate over the nbdl levels
180
181       t_root(:,j) = zero
182
183       DO l = 1, nbdl
184
185          t_root(:,j) = &
186               t_root(:,j) + tsoil(:,l) * rpc(:) * &
187               ( EXP( -z_soil(l-1)/rprof(:,j) ) - EXP( -z_soil(l)/rprof(:,j) ) )
188
189       ENDDO
190
191    ENDDO
192
193    !
194    ! 1.4 total allocatable biomass
195    !
196
197    bm_alloc_tot(:,:) = gpp(:,:) * dt
198
199    !
200    ! 2 define maintenance respiration coefficients
201    !
202
203    DO j = 2,nvm
204
205       !
206       ! 2.1 temperature which is taken for the plant part we are talking about
207       !
208
209       ! 2.1.1 parts above the ground
210
211       t_maint(:,ileaf) = t2m(:)
212       t_maint(:,isapabove) = t2m(:)
213       t_maint(:,ifruit) = t2m(:)
214
215       ! 2.1.2 parts below the ground
216
217       t_maint(:,isapbelow) = t_root(:,j)
218       t_maint(:,iroot) = t_root(:,j)
219
220       ! 2.1.3 heartwood: does not respire. Any temperature
221
222       t_maint(:,iheartbelow) = t_root(:,j)
223       t_maint(:,iheartabove) = t2m(:)
224
225       ! 2.1.4 reserve: above the ground for trees, below for grasses
226
227       IF ( tree(j) ) THEN
228          t_maint(:,icarbres) = t2m(:)
229       ELSE
230          t_maint(:,icarbres) = t_root(:,j)
231       ENDIF
232
233       !
234       ! 2.2 calculate coefficient
235       !
236
237       tl(:) = tlong_ref(:) - ZeroCelsius
238       slope(:) = maint_resp_slope(j,1) + tl(:) * maint_resp_slope(j,2) + &
239            tl(:)*tl(:) * maint_resp_slope(j,3)
240
241       DO k = 1, nparts
242
243          coeff_maint(:,j,k) = &
244               MAX( coeff_maint_zero(j,k) * &
245               ( 1. + slope(:) * (t_maint(:,k)-ZeroCelsius) ), zero )
246
247       ENDDO
248
249    ENDDO
250
251    !
252    ! 3 calculate maintenance and growth respiration.
253    !   NPP = GPP - maintenance resp - growth resp.
254    !
255    !
256    DO j = 2,nvm
257
258       !
259       ! 3.1 maintenance respiration of the different plant parts
260       !
261       !
262       ! 3.2 Total maintenance respiration of the plant
263       !     VPP killer:
264       !     resp_maint(:,j) = SUM( resp_maint_part(:,:), DIM=2 )
265       !
266
267       resp_maint(:,j) = zero
268
269       ! with the new calculation of hourly respiration, we must verify that
270       ! PFT has not been killed after calcul of resp_maint_part in stomate
271       DO k= 1, nparts
272          WHERE (PFTpresent(:,j))
273             resp_maint(:,j) = resp_maint(:,j) + resp_maint_part(:,j,k)
274          ENDWHERE
275       ENDDO
276       !
277       ! 3.3 This maintenance respiration is taken away from the newly produced
278       !     allocatable biomass. However, we avoid that no allocatable biomass remains.
279       !     If the respiration is larger than a given fraction of the allocatable biomass,
280       !     the rest is taken from the tissues themselves.
281       !     We suppose that respiration is not dependent on leaf age ->
282       !     do not change age structure.
283       !
284
285       ! maximum part of allocatable biomass used for respiration
286       bm_tax_max(:) = tax_max * bm_alloc_tot(:,j)
287
288       DO i = 1, npts
289
290          IF ( ( bm_alloc_tot(i,j) .GT. zero ) .AND. &
291               ( ( resp_maint(i,j) * dt ) .LT. bm_tax_max(i) ) ) THEN
292
293             bm_alloc_tot(i,j) = bm_alloc_tot(i,j) - resp_maint(i,j) * dt
294
295             !Shilong        ELSEIF ( resp_maint(i,j) .GT. zero ) THEN
296          ELSEIF ( resp_maint(i,j) .GT. min_stomate ) THEN
297
298             ! remaining allocatable biomass
299             bm_alloc_tot(i,j) = bm_alloc_tot(i,j) - bm_tax_max(i)
300
301             ! biomass that remains to be taken away from tissues
302             bm_pump(i) = resp_maint(i,j) * dt - bm_tax_max(i)
303
304             ! take biomass from tissues
305
306             biomass(i,j,ileaf) = biomass(i,j,ileaf) - &
307                  bm_pump(i) * resp_maint_part(i,j,ileaf) / resp_maint(i,j)
308             biomass(i,j,isapabove) = biomass(i,j,isapabove) - &
309                  bm_pump(i) * resp_maint_part(i,j,isapabove) / resp_maint(i,j)
310             biomass(i,j,isapbelow) = biomass(i,j,isapbelow) - &
311                  bm_pump(i) * resp_maint_part(i,j,isapbelow) / resp_maint(i,j)
312             biomass(i,j,iroot) = biomass(i,j,iroot) - &
313                  bm_pump(i) * resp_maint_part(i,j,iroot) / resp_maint(i,j)
314             biomass(i,j,ifruit) = biomass(i,j,ifruit) - &
315                  bm_pump(i) * resp_maint_part(i,j,ifruit) / resp_maint(i,j)
316             biomass(i,j,icarbres) = biomass(i,j,icarbres) - &
317                  bm_pump(i) * resp_maint_part(i,j,icarbres) / resp_maint(i,j)
318
319          ENDIF
320
321       ENDDO   ! Fortran95: WHERE - ELSEWHERE construct
322
323       !
324       ! 3.4 dispatch allocatable biomass
325       !
326
327       DO k = 1, nparts
328          bm_alloc(:,j,k) = f_alloc(:,j,k) * bm_alloc_tot(:,j)
329       ENDDO
330
331       !
332       ! 3.5 growth respiration of a plant part is a given fraction of the
333       !     remaining allocatable biomass.
334       !
335
336       resp_growth_part(:,:) = frac_growthresp * bm_alloc(:,j,:) / dt
337
338       bm_alloc(:,j,:) = ( 1. - frac_growthresp ) * bm_alloc(:,j,:)
339
340       !
341       ! 3.6 Total growth respiration of the plant
342       !     VPP killer:
343       !     resp_growth(:,j) = SUM( resp_growth_part(:,:), DIM=2 )
344       !
345
346       resp_growth(:,j) = zero
347
348       DO k = 1, nparts
349          resp_growth(:,j) = resp_growth(:,j) + resp_growth_part(:,k)
350       ENDDO
351
352    ENDDO
353
354    !
355    ! 4 update the biomass, but save the old leaf mass for later
356    !   "old" leaf mass is leaf mass after maintenance respiration
357    !
358
359    lm_old(:,:) = biomass(:,:,ileaf)
360
361    biomass(:,:,:) = biomass(:,:,:) + bm_alloc(:,:,:)
362
363    !
364    ! 5 biomass can become negative in some rare cases, as the GPP can be negative
365    !   (dark respiration).
366    !   In this case, set biomass to some small value. This creation of matter is taken into
367    !   account by decreasing the autotrophic respiration. In this case, maintenance respiration
368    !   can become negative !!!
369    !
370
371    DO k = 1, nparts
372
373       DO j = 2,nvm
374
375          WHERE ( biomass(:,j,k) .LT. zero )
376
377             bm_create(:,j) = min_stomate - biomass(:,j,k)
378             
379             biomass(:,j,k) = biomass(:,j,k) + bm_create(:,j)
380             
381             resp_maint(:,j) = resp_maint(:,j) - bm_create(:,j) / dt
382
383          ENDWHERE
384
385       ENDDO
386
387    ENDDO
388
389    !
390    ! 6 Calculate the NPP (gC/(m**2 of ground/day)
391    !
392
393    DO j = 2,nvm
394       npp(:,j) = gpp(:,j) - resp_growth(:,j) - resp_maint(:,j)
395    ENDDO
396
397    !
398    ! 7 leaf age
399    !
400
401    !
402    ! 7.1 Decrease leaf age in youngest class if new leaf biomass is higher than old one.
403    !
404
405    DO j = 2,nvm
406       leaf_mass_young(:,j) = leaf_frac(:,j,1) * lm_old(:,j) + bm_alloc(:,j,ileaf)
407    ENDDO
408
409    DO j = 2,nvm
410       WHERE ( ( bm_alloc(:,j,ileaf) .GT. zero ) .AND. &
411         ( leaf_mass_young(:,j) .GT. zero ) )
412
413          leaf_age(:,j,1) = MAX ( zero, &
414               & leaf_age(:,j,1) * &
415               & ( leaf_mass_young(:,j) - bm_alloc(:,j,ileaf) ) / &
416               & leaf_mass_young(:,j) )
417         
418       ENDWHERE
419    ENDDO
420
421    !
422    ! 7.2 new age class fractions (fraction in youngest class increases)
423    !
424
425    ! 7.2.1 youngest class: new mass in youngest class divided by total new mass
426
427    DO j = 2,nvm
428       WHERE ( biomass(:,j,ileaf) .GT. min_stomate )
429
430          leaf_frac(:,j,1) = leaf_mass_young(:,j) / biomass(:,j,ileaf)
431
432       ENDWHERE
433    ENDDO
434
435    ! 7.2.2 other classes: old mass in leaf age class divided by new mass
436
437    DO m = 2, nleafages
438
439       DO j = 2,nvm
440          WHERE ( biomass(:,j,ileaf) .GT. min_stomate )
441
442             leaf_frac(:,j,m) = leaf_frac(:,j,m) * lm_old(:,j) / biomass(:,j,ileaf)
443
444          ENDWHERE
445       ENDDO
446
447    ENDDO
448
449    !
450    ! 8 Plant age (years)
451    !
452
453    !
454    ! 8.1 Increase age at every time step
455    !
456
457    WHERE ( PFTpresent(:,:) )
458
459       age(:,:) = age(:,:) + dt/one_year
460
461    ELSEWHERE
462
463       age(:,:) = zero
464
465    ENDWHERE
466
467    !
468    ! 8.2 For grasses, decrease age
469    !     if new biomass is higher than old one.
470    !     For trees, age is treated in "establish" if vegetation is dynamic,
471    !     and in turnover routines if it is static (in this case, only take
472    !     into account the age of the heartwood).
473    !
474
475    DO j = 2,nvm
476
477       IF ( .NOT. tree(j) ) THEN
478
479          ! Only four compartments for grasses
480          ! VPP killer:
481          ! bm_new(:) = SUM( biomass(:,j,:), DIM=2 )
482          ! bm_add(:) = SUM( bm_alloc(:,j,:), DIM=2 )
483
484          bm_new(:) = biomass(:,j,ileaf) + biomass(:,j,isapabove) + &
485               biomass(:,j,iroot) + biomass(:,j,ifruit)
486          bm_add(:) = bm_alloc(:,j,ileaf) + bm_alloc(:,j,isapabove) + &
487               bm_alloc(:,j,iroot) + bm_alloc(:,j,ifruit)
488
489          WHERE ( ( bm_new(:) .GT. zero ) .AND. ( bm_add(:) .GT. zero ) )
490             age(:,j) = age(:,j) * ( bm_new(:) - bm_add(:) ) / bm_new(:)
491          ENDWHERE
492
493       ENDIF
494
495    ENDDO
496
497    !
498    ! 9 history
499    !
500
501    CALL histwrite (hist_id_stomate, 'BM_ALLOC_LEAF', itime, &
502         bm_alloc(:,:,ileaf), npts*nvm, horipft_index)
503    CALL histwrite (hist_id_stomate, 'BM_ALLOC_SAP_AB', itime, &
504         bm_alloc(:,:,isapabove), npts*nvm, horipft_index)
505    CALL histwrite (hist_id_stomate, 'BM_ALLOC_SAP_BE', itime, &
506         bm_alloc(:,:,isapbelow), npts*nvm, horipft_index)
507    CALL histwrite (hist_id_stomate, 'BM_ALLOC_ROOT', itime, &
508         bm_alloc(:,:,iroot), npts*nvm, horipft_index)
509    CALL histwrite (hist_id_stomate, 'BM_ALLOC_FRUIT', itime, &
510         bm_alloc(:,:,ifruit), npts*nvm, horipft_index)
511    CALL histwrite (hist_id_stomate, 'BM_ALLOC_RES', itime, &
512         bm_alloc(:,:,icarbres), npts*nvm, horipft_index)
513
514    IF (bavard.GE.4) WRITE(numout,*) 'Leaving npp'
515
516  END SUBROUTINE npp_calc
517
518END MODULE stomate_npp
Note: See TracBrowser for help on using the repository browser.