source: branches/ORCHIDEE_EXT/ORCHIDEE/src_stomate/stomate_alloc.f90 @ 386

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

Correct2 wrong loops in slowproc when dgvm is activated. Replace PRINT instructions by WRITE instructions. Add a call to ipslerr in stomate alloc in case of wrong values of L0 and R0
File size: 18.4 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)) THEN
155       CALL ipslerr (3,'in module stomate_alloc', &
156            &           'Something wrong happened', &
157            &           'L0 negative or division by zero if S0 = 1', &
158            &           '(Check your parameters.)')
159    ENDIF
160
161    !
162    ! 1.1 first call
163    !
164
165    IF ( firstcall ) THEN
166
167       ! 1.1.1 soil levels
168
169       z_soil(0) = zero
170       z_soil(1:nbdl) = diaglev(1:nbdl)
171
172       ! 1.1.2 info about flags and parameters.
173
174       WRITE(numout,*) 'alloc:'
175
176       WRITE(numout,'(a,$)') '    > We'
177       IF ( .NOT. ok_minres ) WRITE(numout,'(a,$)') ' do NOT'
178       WRITE(numout,*) 'try to reach a minumum reservoir when severely stressed.'
179
180       WRITE(numout,*) '   > Time to put initial leaf mass on (d): ',tau_leafinit
181
182       WRITE(numout,*) '   > scaling depth for nitrogen limitation (m): ', &
183            z_nitrogen
184
185       WRITE(numout,*) '   > sap allocation above the ground / total sap allocation: '
186       WRITE(numout,*) '       trees:', alloc_sap_above_tree
187       WRITE(numout,*) '       grasses:', alloc_sap_above_grass
188
189       WRITE(numout,*) '   > standard root alloc fraction: ', R0
190
191       WRITE(numout,*) '   > standard sapwood alloc fraction: ', S0
192
193       WRITE(numout,*) '   > standard fruit allocation: ', f_fruit
194
195       WRITE(numout,*) '   > minimum/maximum leaf alloc fraction: ', min_LtoLSR,max_LtoLSR
196
197       WRITE(numout,*) '   > maximum time (d) during which reserve is used:'
198       WRITE(numout,*) '       trees:',reserve_time_tree
199       WRITE(numout,*) '       grasses:',reserve_time_grass
200
201       firstcall = .FALSE.
202
203    ENDIF
204
205    !
206    ! 1.2 initialize output
207    !
208
209    f_alloc(:,:,:) = zero
210    f_alloc(:,:,icarbres) = un
211    !
212    ! 1.3 Convolution of the temperature and humidity profiles with some kind of profile
213    !     of microbial density gives us a representative temperature and humidity
214    !
215
216    ! 1.3.1 temperature
217
218    ! 1.3.1.1 rpc is an integration constant such that the integral of the root profile is 1.
219    rpc(:) = un / ( un - EXP( -z_soil(nbdl) / z_nitrogen ) )
220
221    ! 1.3.1.2 integrate over the nbdl levels
222
223    t_nitrogen(:) = 0.
224
225    DO l = 1, nbdl
226
227       t_nitrogen(:) = &
228            t_nitrogen(:) + tsoil_month(:,l) * rpc(:) * &
229            ( EXP( -z_soil(l-1)/z_nitrogen ) - EXP( -z_soil(l)/z_nitrogen ) )
230
231    ENDDO
232
233    ! 1.3.2 moisture
234
235    ! 1.3.2.1 rpc is an integration constant such that the integral of the root profile is 1.
236    rpc(:) = un / ( un - EXP( -z_soil(nbdl) / z_nitrogen ) )
237
238    ! 1.3.2.2 integrate over the nbdl levels
239
240    h_nitrogen(:) = zero
241
242    DO l = 1, nbdl
243
244       h_nitrogen(:) = &
245            h_nitrogen(:) + soilhum_month(:,l) * rpc(:) * &
246            ( EXP( -z_soil(l-1)/z_nitrogen ) - EXP( -z_soil(l)/z_nitrogen ) )
247
248    ENDDO
249
250    !
251    ! 1.4 for light limitation: lai on natural part of the grid cell or lai of this
252    !       agricultural PFT
253    !
254
255    ! mask agricultural vegetation
256    ! mean LAI on natural part
257
258    natveg_tot(:) = zero
259    lai_nat(:) = zero
260
261    DO j = 2, nvm
262
263       IF ( natural(j) ) THEN
264          veget_max_nat(:,j) = veget_max(:,j)
265       ELSE
266          veget_max_nat(:,j) = zero
267       ENDIF
268
269       ! sum up fraction of natural space covered by vegetation
270       natveg_tot(:) = natveg_tot(:) + veget_max_nat(:,j)
271
272       ! sum up lai
273       lai_nat(:) = lai_nat(:) + veget_max_nat(:,j) * lai(:,j)
274
275    ENDDO
276
277    DO j = 2, nvm
278
279       IF ( natural(j) ) THEN
280          lai_around(:,j) = lai_nat(:)
281       ELSE
282          lai_around(:,j) = lai(:,j)
283       ENDIF
284
285    ENDDO
286
287    !
288    ! 1.5 LAI below which carbohydrate reserve is used
289    !
290
291    lai_happy(:) = lai_max(:) * lai_max_to_happy
292
293    !
294    ! 2 Use carbohydrate reserve
295    !   This time constant implicitly takes into account the dispersion of the budburst
296    !   data. Therefore, it might be decreased at lower resolution.
297    !
298
299    ! save old leaf mass
300
301    lm_old(:,:) = biomass(:,:,ileaf)
302
303    DO j = 2, nvm
304
305       !
306       ! 2.1 determine mass to be translocated to leaves and roots
307       !
308
309       ! determine maximum time during which reserve is used
310
311       IF ( tree(j) ) THEN
312          reserve_time = reserve_time_tree
313       ELSE
314          reserve_time = reserve_time_grass
315       ENDIF
316
317       ! conditions: 1/ plant must not be senescent
318       !             2/ lai must be relatively low
319       !             3/ must be at the beginning of the growing season
320
321       WHERE ( ( biomass(:,j,ileaf) .GT. zero ) .AND. & 
322            ( .NOT. senescence(:,j) ) .AND. &
323            ( lai(:,j) .LT. lai_happy(j) ) .AND. &
324            ( when_growthinit(:,j) .LT. reserve_time ) ) 
325
326          ! determine mass to put on
327          use_reserve(:) = &
328               MIN( biomass(:,j,icarbres), &
329               deux * dt/tau_leafinit * lai_happy(j)/ sla(j) )
330
331          ! grow leaves and fine roots
332
333          transloc_leaf(:) = L0/(L0+R0) * use_reserve(:)
334
335          biomass(:,j,ileaf) = biomass(:,j,ileaf) + transloc_leaf(:)
336          biomass(:,j,iroot) = biomass(:,j,iroot) + ( use_reserve(:) - transloc_leaf(:) )
337
338          ! decrease reserve mass
339
340          biomass(:,j,icarbres) = biomass(:,j,icarbres) - use_reserve(:)
341
342       ELSEWHERE
343
344          transloc_leaf(:) = zero
345
346       ENDWHERE
347
348       !
349       ! 2.2 update leaf age
350       !
351
352       ! 2.2.1 Decrease leaf age in youngest class.
353
354       leaf_mass_young(:) = leaf_frac(:,j,1) * lm_old(:,j) + transloc_leaf(:)
355
356       WHERE ( ( transloc_leaf(:) .GT. min_stomate ) .AND. ( leaf_mass_young(:) .GT. min_stomate ) )
357
358          leaf_age(:,j,1) = MAX( zero, leaf_age(:,j,1) * ( leaf_mass_young(:) - transloc_leaf(:) ) / &
359               leaf_mass_young(:) )
360
361       ENDWHERE
362
363       ! 2.2.2 new age class fractions (fraction in youngest class increases)
364
365       ! 2.2.2.1 youngest class: new mass in youngest class divided by total new mass
366
367       WHERE ( biomass(:,j,ileaf) .GT. min_stomate )
368
369          leaf_frac(:,j,1) = leaf_mass_young(:) / biomass(:,j,ileaf)
370
371       ENDWHERE
372
373       ! 2.2.2.2 other classes: old mass in leaf age class divided by new mass
374
375       DO m = 2, nleafages
376
377          WHERE ( biomass(:,j,ileaf) .GT. min_stomate )
378
379             leaf_frac(:,j,m) = leaf_frac(:,j,m) * lm_old(:,j) / biomass(:,j,ileaf)
380
381          ENDWHERE
382
383       ENDDO
384
385    ENDDO         ! loop over PFTs
386
387    !
388    ! 3 Calculate fractional allocation.
389    !   The fractions of NPP allocated to the different compartments depend on the
390    !   availability of light, water, and nitrogen.
391    !
392
393    DO j = 2, nvm
394
395       RtoLSR(:)=0
396       LtoLSR(:)=0
397       StoLSR(:)=0
398
399       ! for the moment, fixed partitioning between above and below the ground
400       ! modified by JO/NV/PF for changing partitioning with stand age
401       ! we could have alloc_sap_above(npts,nvm) but we have only
402       ! alloc_sap_above(npts) as we make a loop over j=2,nvm
403       !
404       IF ( tree(j) ) THEN
405
406          alloc_sap_above (:) = alloc_min(j)+(alloc_max(j)-alloc_min(j))*(1.-EXP(-age(:,j)/demi_alloc(j)))
407
408          !IF (j .EQ. 3) WRITE(*,*) '%allocated above = 'alloc_sap_above(1),'age = ',age(1,j)
409       ELSE
410          alloc_sap_above(:) = alloc_sap_above_grass
411       ENDIF
412
413       ! only where leaves are on
414
415       WHERE ( biomass(:,j,ileaf) .GT. min_stomate )
416
417          !
418          ! 3.1 Limiting factors: weak value = strong limitation
419          !
420
421          ! 3.1.1 Light: depends on mean lai on the natural part of the
422          !       grid box (light competition).
423          !       For agricultural PFTs, take its own lai for both parts.
424          !MM, NV
425          WHERE( lai_around(:,j) < 10 )
426             limit_L(:) = MAX( 0.1_r_std, EXP( -undemi * lai_around(:,j) ) )
427          ELSEWHERE
428             limit_L(:) = 0.1_r_std
429          ENDWHERE
430          ! 3.1.2 Water
431
432          limit_W(:) = MAX( 0.1_r_std, MIN( un, moiavail_week(:,j) ) )
433
434          ! 3.1.3 Nitrogen supply: depends on water and temperature
435          !       Agricultural PFTs can be limited by Nitrogen for the moment ...
436          !       Replace this once there is a nitrogen cycle in STOMATE !
437
438          ! 3.1.3.1 water
439
440          limit_N_hum(:) = MAX( undemi, MIN( un, h_nitrogen(:) ) )
441
442          ! 3.1.3.2 temperature
443
444          limit_N_temp(:) = 2.**((t_nitrogen(:)-ZeroCelsius-25.)/10.)
445          limit_N_temp(:) = MAX( 0.1_r_std, MIN( un, limit_N_temp(:) ) )
446
447          ! 3.1.3.3 combine water and temperature factors to get nitrogen limitation
448
449          limit_N(:) = MAX( 0.1_r_std, MIN( un, limit_N_hum(:) * limit_N_temp(:) ) )
450
451          ! 3.1.4 Among water and nitrogen, take the one that is more limited
452
453          limit_WorN(:) = MIN( limit_W(:), limit_N(:) )
454
455          ! 3.1.5 strongest limitation
456
457          limit(:) = MIN( limit_WorN(:), limit_L(:) )
458
459          !
460          ! 3.2 Ratio between allocation to leaves, sapwood and roots
461          !
462
463          ! preliminary root allocation
464
465          RtoLSR(:) = &
466               MAX( .15_r_std, &
467               R0 * trois * limit_L(:) / ( limit_L(:) + deux * limit_WorN(:) ) )
468
469          ! sapwood allocation
470
471          StoLSR(:) = S0 * 3. * limit_WorN(:) / ( 2. * limit_L(:) + limit_WorN(:) )
472
473          ! leaf allocation
474
475          LtoLSR(:) = un - RtoLSR(:) - StoLSR(:)
476          LtoLSR(:) = MAX( min_LtoLSR, MIN( max_LtoLSR, LtoLSR(:) ) )
477
478          ! roots: the rest
479
480          RtoLSR(:) = un - LtoLSR(:) - StoLSR(:)
481
482       ENDWHERE
483
484       ! no leaf allocation if LAI beyond maximum LAI. Biomass then goes into sapwood
485
486       WHERE ( (biomass(:,j,ileaf) .GT. min_stomate) .AND. (lai(:,j) .GT. lai_max(j)) )
487
488          StoLSR(:) = StoLSR(:) + LtoLSR(:)
489
490          LtoLSR(:) = zero
491
492       ENDWHERE
493
494       !
495       ! 3.3 final allocation
496       !
497
498       DO i = 1, npts
499
500          IF ( biomass(i,j,ileaf) .GT. min_stomate ) THEN
501
502             IF ( senescence(i,j) ) THEN
503
504                ! 3.3.1 senescent: everything goes into carbohydrate reserve
505
506                f_alloc(i,j,icarbres) = 1.0
507
508             ELSE
509
510                ! 3.3.2 in growing season
511
512                ! to fruits
513                f_alloc(i,j,ifruit) = f_fruit
514
515                ! allocation to the reserve is proportional to the leaf and root allocation.
516                ! Leaf, root, and sap allocation are rescaled.
517                ! No allocation to reserve if there is much biomass in it
518                !   (more than the maximum LAI: in that case, rescale=1)
519
520                IF ( ( biomass(i,j,icarbres)*sla(j) ) .LT. 2*lai_max(j) ) THEN
521                   carb_rescale(i) = un / ( un + ecureuil(j) * ( LtoLSR(i) + RtoLSR(i) ) )
522                ELSE
523                   carb_rescale(i) = un
524                ENDIF
525
526                f_alloc(i,j,ileaf) = LtoLSR(i) * ( 1.-f_alloc(i,j,ifruit) ) * carb_rescale(i)
527
528                f_alloc(i,j,isapabove) = StoLSR(i) * alloc_sap_above(i) * &
529                     ( un - f_alloc(i,j,ifruit) ) * carb_rescale(i)
530                f_alloc(i,j,isapbelow) = StoLSR(i) * ( un - alloc_sap_above(i) ) * &
531                     ( un - f_alloc(i,j,ifruit) ) * carb_rescale(i)
532
533                f_alloc(i,j,iroot) = RtoLSR(i) * ( 1.-f_alloc(i,j,ifruit) ) * carb_rescale(i)
534
535                ! this is equivalent to:
536                ! reserve alloc = ecureuil*(LtoLSR+StoLSR)*(1-fruit_alloc)*carb_rescale
537                f_alloc(i,j,icarbres) = ( un - carb_rescale(i) ) * ( 1.-f_alloc(i,j,ifruit) )
538
539             ENDIF  ! senescent?
540
541          ENDIF    ! there are leaves
542
543       ENDDO      ! Fortran95: double WHERE construct
544
545    ENDDO          ! loop over PFTs
546
547    !
548    ! 4 root profile
549    !
550
551
552    IF (bavard.GE.4) WRITE(numout,*) 'Leaving alloc'
553
554  END SUBROUTINE alloc
555
556
557END MODULE stomate_alloc
Note: See TracBrowser for help on using the repository browser.