source: branches/ORCHIDEE_EXT/ORCHIDEE/src_stomate/stomate_alloc.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: 18.2 KB
Line 
1! allocation to the roots, stems, leaves, "fruits" and carbohydrate reserve.
2! Reproduction: for the moment, this is simply a 10% "tax".
3! This should depend on the limitations that the plant experiences. If the
4! plant fares well, it will have fruits. However, this means that we should
5! also "reward" the plants for having grown fruits by making the
6! reproduction rate depend on the fruit growth of the past years. Otherwise,
7! the fruit allocation would be a punishment for plants that are doing well.
8! "calculates" root profiles (in fact, prescribes it for the moment).
9!
10! $Header: /home/ssipsl/CVSREP/ORCHIDEE/src_stomate/stomate_alloc.f90,v 1.10 2009/03/31 12:11:22 ssipsl Exp $
11! IPSL (2006)
12!  This software is governed by the CeCILL licence see ORCHIDEE/ORCHIDEE_CeCILL.LIC
13!
14MODULE stomate_alloc
15
16  ! modules used:
17
18  USE ioipsl
19  USE pft_parameters
20  USE stomate_data
21  USE constantes
22
23  IMPLICIT NONE
24
25  ! private & public routines
26
27  PRIVATE
28  PUBLIC alloc,alloc_clear
29
30  ! first call
31  LOGICAL, SAVE                                             :: firstcall = .TRUE.
32CONTAINS
33  SUBROUTINE alloc_clear
34    firstcall = .TRUE.
35  END SUBROUTINE alloc_clear
36
37  SUBROUTINE alloc (npts, dt, &
38       lai, veget_max, senescence, when_growthinit, &
39       moiavail_week, tsoil_month, soilhum_month, &
40       biomass, age, leaf_age, leaf_frac, rprof, f_alloc)
41
42    !
43    ! 0 declarations
44    !
45
46    ! 0.1 input
47
48    ! Domain size
49    INTEGER(i_std), INTENT(in)                                       :: npts
50    ! time step (days)
51    REAL(r_std), INTENT(in)                                    :: dt
52    ! Leaf area index
53    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)              :: lai
54    ! "maximal" coverage fraction of a PFT ( = ind*cn_ind )
55    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)              :: veget_max
56    ! is the plant senescent? (only for deciduous trees - carbohydrate reserve)
57    LOGICAL, DIMENSION(npts,nvm), INTENT(in)                 :: senescence
58    ! how many days ago was the beginning of the growing season
59    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)              :: when_growthinit
60    ! "weekly" moisture availability
61    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)              :: moiavail_week
62    ! "monthly" soil temperature (K)
63    REAL(r_std), DIMENSION(npts,nbdl), INTENT(in)              :: tsoil_month
64    ! "monthly" soil humidity
65    REAL(r_std), DIMENSION(npts,nbdl), INTENT(in)              :: soilhum_month
66    !  age (days)
67    REAL(r_std), DIMENSION(npts,nvm), INTENT(in)              :: age
68
69    ! 0.2 modified fields
70
71    ! biomass (gC/m**2)
72    REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(inout)    :: biomass
73    ! leaf age (days)
74    REAL(r_std), DIMENSION(npts,nvm,nleafages), INTENT(inout) :: leaf_age
75    ! fraction of leaves in leaf age class
76    REAL(r_std), DIMENSION(npts,nvm,nleafages), INTENT(inout) :: leaf_frac
77
78    ! 0.3 output
79
80    ! root depth. This will, one day, be a prognostic variable. It will be calculated by
81    ! STOMATE (save in restart file & give to hydrology module!). For the moment, it
82    ! is prescribed.
83    REAL(r_std), DIMENSION(npts,nvm), INTENT(inout)             :: rprof
84    ! fraction that goes into plant part
85    REAL(r_std), DIMENSION(npts,nvm,nparts), INTENT(out)      :: f_alloc
86
87    ! 0.4 local
88
89    ! below this lai, the carbohydrate reserve is used
90    REAL(r_std), DIMENSION(nvm)                               :: lai_happy
91    ! limiting factor light
92    REAL(r_std), DIMENSION(npts)                               :: limit_L
93    ! limiting factor nitrogen
94    REAL(r_std), DIMENSION(npts)                               :: limit_N
95    ! factors determining limit_N: 1/ temperature
96    REAL(r_std), DIMENSION(npts)                               :: limit_N_temp
97    ! factors determining limit_N: 2/ humidity
98    REAL(r_std), DIMENSION(npts)                               :: limit_N_hum
99    ! limiting factor water
100    REAL(r_std), DIMENSION(npts)                               :: limit_W
101    ! limiting factor in soil (nitrogen or water)
102    REAL(r_std), DIMENSION(npts)                               :: limit_WorN
103    ! limit: strongest limitation amongst limit_N, limit_W and limit_L
104    REAL(r_std), DIMENSION(npts)                               :: limit
105    ! soil temperature used for N parameterization
106    REAL(r_std), DIMENSION(npts)                               :: t_nitrogen
107    ! soil humidity used for N parameterization
108    REAL(r_std), DIMENSION(npts)                               :: h_nitrogen
109    ! integration constant for vertical profiles
110    REAL(r_std), DIMENSION(npts)                               :: rpc
111    ! ratio between leaf-allocation and (leaf+sapwood+root)-allocation
112    REAL(r_std), DIMENSION(npts)                               :: LtoLSR
113    ! ratio between sapwood-allocation and (leaf+sapwood+root)-allocation
114    REAL(r_std), DIMENSION(npts)                               :: StoLSR
115    ! ratio between root-allocation and (leaf+sapwood+root)-allocation
116    REAL(r_std), DIMENSION(npts)                               :: RtoLSR
117    ! rescaling factor for carbohydrate reserve allocation
118    REAL(r_std), DIMENSION(npts)                               :: carb_rescale
119    ! mass taken from carbohydrate reserve (gC/m**2)
120    REAL(r_std), DIMENSION(npts)                               :: use_reserve
121    ! mass taken from carbohydrate reserve and put into leaves (gC/m**2)
122    REAL(r_std), DIMENSION(npts)                               :: transloc_leaf
123    ! mass in youngest leaf age class (gC/m**2)
124    REAL(r_std), DIMENSION(npts)                               :: leaf_mass_young
125    ! old leaf biomass (gC/m**2)
126    REAL(r_std), DIMENSION(npts,nvm)                          :: lm_old
127    ! maximum time (d) during which reserve is used
128    REAL(r_std)                                                :: reserve_time
129    ! lai on natural part of the grid cell, or of this agricultural PFT
130    REAL(r_std), DIMENSION(npts,nvm)                          :: lai_around
131    ! vegetation cover of natural PFTs on the grid cell (agriculture masked)
132    REAL(r_std), DIMENSION(npts,nvm)                          :: veget_max_nat 
133    ! total natural vegetation cover on natural part of the grid cell
134    REAL(r_std), DIMENSION(npts)                               :: natveg_tot
135    ! average LAI on natural part of the grid cell
136    REAL(r_std), DIMENSION(npts)                               :: lai_nat
137    ! intermediate array for looking for minimum
138    REAL(r_std), DIMENSION(npts)                               :: zdiff_min
139    ! fraction of sapwood allocation above ground (SHOULD BE CALCULATED !!!!)
140    REAL(r_std), DIMENSION(npts)                               :: alloc_sap_above
141    ! soil levels (m)
142    REAL(r_std), SAVE, DIMENSION(0:nbdl)                       :: z_soil
143    ! Index
144    INTEGER(i_std)                                            :: i,j,l,m
145
146    ! =========================================================================
147
148    IF (bavard.GE.3) WRITE(numout,*) 'Entering alloc'
149   
150    !
151    ! 1 Initialization
152    !
153    L0 = 1. - R0 - S0  ! defined in constantes.f90
154    IF ((L0 .LT. zero) .OR. (S0 .EQ. un)) STOP 'L0 negative or division by zero if S0 = 1' 
155
156    !
157    ! 1.1 first call
158    !
159
160    IF ( firstcall ) THEN
161
162       ! 1.1.1 soil levels
163
164       z_soil(0) = zero
165       z_soil(1:nbdl) = diaglev(1:nbdl)
166
167       ! 1.1.2 info about flags and parameters.
168
169       WRITE(numout,*) 'alloc:'
170
171       WRITE(numout,'(a,$)') '    > We'
172       IF ( .NOT. ok_minres ) WRITE(numout,'(a,$)') ' do NOT'
173       WRITE(numout,*) 'try to reach a minumum reservoir when severely stressed.'
174
175       WRITE(numout,*) '   > Time to put initial leaf mass on (d): ',tau_leafinit
176
177       WRITE(numout,*) '   > scaling depth for nitrogen limitation (m): ', &
178            z_nitrogen
179
180       WRITE(numout,*) '   > sap allocation above the ground / total sap allocation: '
181       WRITE(numout,*) '       trees:', alloc_sap_above_tree
182       WRITE(numout,*) '       grasses:', alloc_sap_above_grass
183
184       WRITE(numout,*) '   > standard root alloc fraction: ', R0
185
186       WRITE(numout,*) '   > standard sapwood alloc fraction: ', S0
187
188       WRITE(numout,*) '   > standard fruit allocation: ', f_fruit
189
190       WRITE(numout,*) '   > minimum/maximum leaf alloc fraction: ', min_LtoLSR,max_LtoLSR
191
192       WRITE(numout,*) '   > maximum time (d) during which reserve is used:'
193       WRITE(numout,*) '       trees:',reserve_time_tree
194       WRITE(numout,*) '       grasses:',reserve_time_grass
195
196       firstcall = .FALSE.
197
198    ENDIF
199
200    !
201    ! 1.2 initialize output
202    !
203
204    f_alloc(:,:,:) = zero
205    f_alloc(:,:,icarbres) = un
206    !
207    ! 1.3 Convolution of the temperature and humidity profiles with some kind of profile
208    !     of microbial density gives us a representative temperature and humidity
209    !
210
211    ! 1.3.1 temperature
212
213    ! 1.3.1.1 rpc is an integration constant such that the integral of the root profile is 1.
214    rpc(:) = un / ( un - EXP( -z_soil(nbdl) / z_nitrogen ) )
215
216    ! 1.3.1.2 integrate over the nbdl levels
217
218    t_nitrogen(:) = 0.
219
220    DO l = 1, nbdl
221
222       t_nitrogen(:) = &
223            t_nitrogen(:) + tsoil_month(:,l) * rpc(:) * &
224            ( EXP( -z_soil(l-1)/z_nitrogen ) - EXP( -z_soil(l)/z_nitrogen ) )
225
226    ENDDO
227
228    ! 1.3.2 moisture
229
230    ! 1.3.2.1 rpc is an integration constant such that the integral of the root profile is 1.
231    rpc(:) = un / ( un - EXP( -z_soil(nbdl) / z_nitrogen ) )
232
233    ! 1.3.2.2 integrate over the nbdl levels
234
235    h_nitrogen(:) = zero
236
237    DO l = 1, nbdl
238
239       h_nitrogen(:) = &
240            h_nitrogen(:) + soilhum_month(:,l) * rpc(:) * &
241            ( EXP( -z_soil(l-1)/z_nitrogen ) - EXP( -z_soil(l)/z_nitrogen ) )
242
243    ENDDO
244
245    !
246    ! 1.4 for light limitation: lai on natural part of the grid cell or lai of this
247    !       agricultural PFT
248    !
249
250    ! mask agricultural vegetation
251    ! mean LAI on natural part
252
253    natveg_tot(:) = zero
254    lai_nat(:) = zero
255
256    DO j = 2, nvm
257
258       IF ( natural(j) ) THEN
259          veget_max_nat(:,j) = veget_max(:,j)
260       ELSE
261          veget_max_nat(:,j) = zero
262       ENDIF
263
264       ! sum up fraction of natural space covered by vegetation
265       natveg_tot(:) = natveg_tot(:) + veget_max_nat(:,j)
266
267       ! sum up lai
268       lai_nat(:) = lai_nat(:) + veget_max_nat(:,j) * lai(:,j)
269
270    ENDDO
271
272    DO j = 2, nvm
273
274       IF ( natural(j) ) THEN
275          lai_around(:,j) = lai_nat(:)
276       ELSE
277          lai_around(:,j) = lai(:,j)
278       ENDIF
279
280    ENDDO
281
282    !
283    ! 1.5 LAI below which carbohydrate reserve is used
284    !
285
286    lai_happy(:) = lai_max(:) * lai_max_to_happy
287
288    !
289    ! 2 Use carbohydrate reserve
290    !   This time constant implicitly takes into account the dispersion of the budburst
291    !   data. Therefore, it might be decreased at lower resolution.
292    !
293
294    ! save old leaf mass
295
296    lm_old(:,:) = biomass(:,:,ileaf)
297
298    DO j = 2, nvm
299
300       !
301       ! 2.1 determine mass to be translocated to leaves and roots
302       !
303
304       ! determine maximum time during which reserve is used
305
306       IF ( tree(j) ) THEN
307          reserve_time = reserve_time_tree
308       ELSE
309          reserve_time = reserve_time_grass
310       ENDIF
311
312       ! conditions: 1/ plant must not be senescent
313       !             2/ lai must be relatively low
314       !             3/ must be at the beginning of the growing season
315
316       WHERE ( ( biomass(:,j,ileaf) .GT. zero ) .AND. & 
317            ( .NOT. senescence(:,j) ) .AND. &
318            ( lai(:,j) .LT. lai_happy(j) ) .AND. &
319            ( when_growthinit(:,j) .LT. reserve_time ) ) 
320
321          ! determine mass to put on
322          use_reserve(:) = &
323               MIN( biomass(:,j,icarbres), &
324               deux * dt/tau_leafinit * lai_happy(j)/ sla(j) )
325
326          ! grow leaves and fine roots
327
328          transloc_leaf(:) = L0/(L0+R0) * use_reserve(:)
329
330          biomass(:,j,ileaf) = biomass(:,j,ileaf) + transloc_leaf(:)
331          biomass(:,j,iroot) = biomass(:,j,iroot) + ( use_reserve(:) - transloc_leaf(:) )
332
333          ! decrease reserve mass
334
335          biomass(:,j,icarbres) = biomass(:,j,icarbres) - use_reserve(:)
336
337       ELSEWHERE
338
339          transloc_leaf(:) = zero
340
341       ENDWHERE
342
343       !
344       ! 2.2 update leaf age
345       !
346
347       ! 2.2.1 Decrease leaf age in youngest class.
348
349       leaf_mass_young(:) = leaf_frac(:,j,1) * lm_old(:,j) + transloc_leaf(:)
350
351       WHERE ( ( transloc_leaf(:) .GT. min_stomate ) .AND. ( leaf_mass_young(:) .GT. min_stomate ) )
352
353          leaf_age(:,j,1) = MAX( zero, leaf_age(:,j,1) * ( leaf_mass_young(:) - transloc_leaf(:) ) / &
354               leaf_mass_young(:) )
355
356       ENDWHERE
357
358       ! 2.2.2 new age class fractions (fraction in youngest class increases)
359
360       ! 2.2.2.1 youngest class: new mass in youngest class divided by total new mass
361
362       WHERE ( biomass(:,j,ileaf) .GT. min_stomate )
363
364          leaf_frac(:,j,1) = leaf_mass_young(:) / biomass(:,j,ileaf)
365
366       ENDWHERE
367
368       ! 2.2.2.2 other classes: old mass in leaf age class divided by new mass
369
370       DO m = 2, nleafages
371
372          WHERE ( biomass(:,j,ileaf) .GT. min_stomate )
373
374             leaf_frac(:,j,m) = leaf_frac(:,j,m) * lm_old(:,j) / biomass(:,j,ileaf)
375
376          ENDWHERE
377
378       ENDDO
379
380    ENDDO         ! loop over PFTs
381
382    !
383    ! 3 Calculate fractional allocation.
384    !   The fractions of NPP allocated to the different compartments depend on the
385    !   availability of light, water, and nitrogen.
386    !
387
388    DO j = 2, nvm
389
390       RtoLSR(:)=0
391       LtoLSR(:)=0
392       StoLSR(:)=0
393
394       ! for the moment, fixed partitioning between above and below the ground
395       ! modified by JO/NV/PF for changing partitioning with stand age
396       ! we could have alloc_sap_above(npts,nvm) but we have only
397       ! alloc_sap_above(npts) as we make a loop over j=2,nvm
398       !
399       IF ( tree(j) ) THEN
400
401          alloc_sap_above (:) = alloc_min(j)+(alloc_max(j)-alloc_min(j))*(1.-EXP(-age(:,j)/demi_alloc(j)))
402
403          !IF (j .EQ. 3) WRITE(*,*) '%allocated above = 'alloc_sap_above(1),'age = ',age(1,j)
404       ELSE
405          alloc_sap_above(:) = alloc_sap_above_grass
406       ENDIF
407
408       ! only where leaves are on
409
410       WHERE ( biomass(:,j,ileaf) .GT. min_stomate )
411
412          !
413          ! 3.1 Limiting factors: weak value = strong limitation
414          !
415
416          ! 3.1.1 Light: depends on mean lai on the natural part of the
417          !       grid box (light competition).
418          !       For agricultural PFTs, take its own lai for both parts.
419          !MM, NV
420          WHERE( lai_around(:,j) < 10 )
421             limit_L(:) = MAX( 0.1_r_std, EXP( -undemi * lai_around(:,j) ) )
422          ELSEWHERE
423             limit_L(:) = 0.1_r_std
424          ENDWHERE
425          ! 3.1.2 Water
426
427          limit_W(:) = MAX( 0.1_r_std, MIN( un, moiavail_week(:,j) ) )
428
429          ! 3.1.3 Nitrogen supply: depends on water and temperature
430          !       Agricultural PFTs can be limited by Nitrogen for the moment ...
431          !       Replace this once there is a nitrogen cycle in STOMATE !
432
433          ! 3.1.3.1 water
434
435          limit_N_hum(:) = MAX( undemi, MIN( un, h_nitrogen(:) ) )
436
437          ! 3.1.3.2 temperature
438
439          limit_N_temp(:) = 2.**((t_nitrogen(:)-ZeroCelsius-25.)/10.)
440          limit_N_temp(:) = MAX( 0.1_r_std, MIN( un, limit_N_temp(:) ) )
441
442          ! 3.1.3.3 combine water and temperature factors to get nitrogen limitation
443
444          limit_N(:) = MAX( 0.1_r_std, MIN( un, limit_N_hum(:) * limit_N_temp(:) ) )
445
446          ! 3.1.4 Among water and nitrogen, take the one that is more limited
447
448          limit_WorN(:) = MIN( limit_W(:), limit_N(:) )
449
450          ! 3.1.5 strongest limitation
451
452          limit(:) = MIN( limit_WorN(:), limit_L(:) )
453
454          !
455          ! 3.2 Ratio between allocation to leaves, sapwood and roots
456          !
457
458          ! preliminary root allocation
459
460          RtoLSR(:) = &
461               MAX( .15_r_std, &
462               R0 * trois * limit_L(:) / ( limit_L(:) + deux * limit_WorN(:) ) )
463
464          ! sapwood allocation
465
466          StoLSR(:) = S0 * 3. * limit_WorN(:) / ( 2. * limit_L(:) + limit_WorN(:) )
467
468          ! leaf allocation
469
470          LtoLSR(:) = un - RtoLSR(:) - StoLSR(:)
471          LtoLSR(:) = MAX( min_LtoLSR, MIN( max_LtoLSR, LtoLSR(:) ) )
472
473          ! roots: the rest
474
475          RtoLSR(:) = un - LtoLSR(:) - StoLSR(:)
476
477       ENDWHERE
478
479       ! no leaf allocation if LAI beyond maximum LAI. Biomass then goes into sapwood
480
481       WHERE ( (biomass(:,j,ileaf) .GT. min_stomate) .AND. (lai(:,j) .GT. lai_max(j)) )
482
483          StoLSR(:) = StoLSR(:) + LtoLSR(:)
484
485          LtoLSR(:) = zero
486
487       ENDWHERE
488
489       !
490       ! 3.3 final allocation
491       !
492
493       DO i = 1, npts
494
495          IF ( biomass(i,j,ileaf) .GT. min_stomate ) THEN
496
497             IF ( senescence(i,j) ) THEN
498
499                ! 3.3.1 senescent: everything goes into carbohydrate reserve
500
501                f_alloc(i,j,icarbres) = 1.0
502
503             ELSE
504
505                ! 3.3.2 in growing season
506
507                ! to fruits
508                f_alloc(i,j,ifruit) = f_fruit
509
510                ! allocation to the reserve is proportional to the leaf and root allocation.
511                ! Leaf, root, and sap allocation are rescaled.
512                ! No allocation to reserve if there is much biomass in it
513                !   (more than the maximum LAI: in that case, rescale=1)
514
515                IF ( ( biomass(i,j,icarbres)*sla(j) ) .LT. 2*lai_max(j) ) THEN
516                   carb_rescale(i) = un / ( un + ecureuil(j) * ( LtoLSR(i) + RtoLSR(i) ) )
517                ELSE
518                   carb_rescale(i) = un
519                ENDIF
520
521                f_alloc(i,j,ileaf) = LtoLSR(i) * ( 1.-f_alloc(i,j,ifruit) ) * carb_rescale(i)
522
523                f_alloc(i,j,isapabove) = StoLSR(i) * alloc_sap_above(i) * &
524                     ( un - f_alloc(i,j,ifruit) ) * carb_rescale(i)
525                f_alloc(i,j,isapbelow) = StoLSR(i) * ( un - alloc_sap_above(i) ) * &
526                     ( un - f_alloc(i,j,ifruit) ) * carb_rescale(i)
527
528                f_alloc(i,j,iroot) = RtoLSR(i) * ( 1.-f_alloc(i,j,ifruit) ) * carb_rescale(i)
529
530                ! this is equivalent to:
531                ! reserve alloc = ecureuil*(LtoLSR+StoLSR)*(1-fruit_alloc)*carb_rescale
532                f_alloc(i,j,icarbres) = ( un - carb_rescale(i) ) * ( 1.-f_alloc(i,j,ifruit) )
533
534             ENDIF  ! senescent?
535
536          ENDIF    ! there are leaves
537
538       ENDDO      ! Fortran95: double WHERE construct
539
540    ENDDO          ! loop over PFTs
541
542    !
543    ! 4 root profile
544    !
545
546
547    IF (bavard.GE.4) WRITE(numout,*) 'Leaving alloc'
548
549  END SUBROUTINE alloc
550
551
552END MODULE stomate_alloc
Note: See TracBrowser for help on using the repository browser.