Documentation/ORCHIDEE_DOC: diffuco_trans_co2.f90

File diffuco_trans_co2.f90, 28.1 KB (added by maignan, 9 years ago)
Line 
1!! ==============================================================================\n
2!! SUBROUTINE   : diffuco_trans_co2
3!! AUTHOR       :
4!! CREATION DATE:
5!!
6!> BRIEF        : This subroutine computes carbon assimilation and stomatal
7!! conductance, following respectively Farqhuar et al. (1980) and Ball et al. (1987).\n
8!!
9!! DESCRIPTION (functional, design, flags):\n
10!! The equations are different depending on the photosynthesis mode (C3 versus C4).\n
11!! Assimilation and conductance are computed over 20 levels of LAI and then
12!! integrated at the canopy level.\n
13!! This routine also computes partial beta coefficient: transpiration for each
14!! type of vegetation.\n
15!! There is a main loop on the PFTs, then inner loops on the points where
16!! assimilation has to be calculated.\n
17!! This subroutine is called by diffuco_main only if photosynthesis is activated
18!! for sechiba (flag STOMATE_OK_CO2=TRUE), otherwise diffuco_trans is called.\n
19!!
20!! REFERENCES   :
21!! - Ball, J., T. Woodrow, and J. Berry (1987), A model predicting stomatal
22!! conductance and its contribution to the control of photosynthesis under
23!! different environmental conditions, Prog. Photosynthesis, 4, 221– 224.
24!! - Collatz, G., M. Ribas-Carbo, and J. Berry (1992), Coupled photosynthesis
25!! stomatal conductance model for leaves of C4 plants, Aust. J. Plant Physiol.,
26!! 19, 519–538.
27!! - Farquhar, G., S. von Caemmener, and J. Berry (1980), A biochemical model of
28!! photosynthesis CO2 fixation in leaves of C3 species, Planta, 149, 78–90.
29!!
30!! FLOWCHART   :
31!!
32!! REVISIONS    : N. de Noblet     2006/06
33!!                - addition of q2m and t2m as input parameters for the
34!!                calculation of Rveget
35!!                - introduction of vbeta23\n
36!! ==============================================================================
37SUBROUTINE diffuco_trans_co2 (kjpindex, dtradia, swdown, temp_air, pb, qair, q2m, t2m, rau, u, v, q_cdrag, humrel, &
38                                assim_param, ccanopy, &
39                                veget, veget_max, lai, qsintveg, qsintmax, vbeta3, rveget, rstruct, cimean, vbetaco2, vbeta23)
40!! INTERFACE DESCRIPTION
41
42!! INPUT SCALAR   
43    INTEGER(i_std), INTENT(in)                               :: kjpindex         !! Domain size (-)
44    REAL(r_std), INTENT (in)                                 :: dtradia          !! Time step (s)
45!! INPUT FIELSS   
46    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: swdown           !! Downwelling short wave flux (W/m^2)
47    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: temp_air         !! Air temperature (K)
48    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: pb               !! Lowest level pressure (Pa)
49    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: qair             !! Lowest level specific humidity (kg/kg)
50! N. de Noblet - 2006/06 - addition of q2m and t2m
51    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: q2m              !! 2m specific humidity (kg/kg)
52    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: t2m              !! 2m air temperature (K)
53! N. de Noblet - 2006/06 - addition of q2m and t2m
54    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: rau              !! air density (kg/m3)
55    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: u                !! Lowest level wind speed (m/s)
56    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: v                !! Lowest level wind speed (m/s)
57    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: q_cdrag          !! Surface drag (m/s)
58   
59    REAL(r_std),DIMENSION (kjpindex,nvm,npco2), INTENT (in)  :: assim_param      !! min+max+opt temps, vcmax, vjmax for photosynthesis (K, umol/m^2/s) 
60    REAL(r_std),DIMENSION (kjpindex), INTENT (in)            :: ccanopy          !! CO2 concentration inside the canopy (ppm)
61    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)        :: humrel           !! Soil moisture stress (-)
62    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)        :: veget            !! Type of vegetation fraction (fraction)
63    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)        :: veget_max        !! Max. vegetation fraction (LAI -> infty) (fraction)
64    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)        :: lai              !! Leaf area index (m^2/m^2)
65    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)        :: qsintveg         !! Water on vegetation due to interception (kg/m^2)
66    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)        :: qsintmax         !! Maximum water on vegetation (kg/m^2)
67! N. de Noblet - 2006/06
68    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (in)        :: vbeta23          !! Beta for fraction of wetted foliage that will transpire (mm/d)   
69! N. de Noblet - 2006/06
70
71!! OUTPUT FIELDS   
72    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out)       :: vbeta3           !! Beta for Transpiration (mm/d)
73    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out)       :: rveget           !! Surface resistance of vegetation (s/m)
74    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out)       :: rstruct          !! structural resistance (s/m)
75    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out)       :: cimean           !! mean intercellular CO2 concentration (umole/m^2/s)
76    REAL(r_std),DIMENSION (kjpindex,nvm), INTENT (out)       :: vbetaco2         !! beta for CO2 (mm/d)
77
78!! LOCAL VARIABLES
79    REAL(r_std),DIMENSION (kjpindex,nvm)  :: vcmax                               !! maximum rate of carboxylation (umol/m^2/s)
80    REAL(r_std),DIMENSION (kjpindex,nvm)  :: vjmax                               !! maximum rate of Rubisco regeneration (umol/m^2/s)
81    REAL(r_std),DIMENSION (kjpindex,nvm)  :: tmin                                !! minimum temperature for photosynthesis (K)
82    REAL(r_std),DIMENSION (kjpindex,nvm)  :: topt                                !! optimal temperature for photosynthesis (K)
83    REAL(r_std),DIMENSION (kjpindex,nvm)  :: tmax                                !! maximum temperature for photosynthesis (K)
84    INTEGER(i_std)                        :: ji, jv, jl                          !! indices (-)
85    REAL(r_std), DIMENSION(kjpindex)      :: leaf_ci_lowest                      !! intercellular CO2 concentration at the lowest LAI level (ppm)
86    INTEGER(i_std), DIMENSION(kjpindex)   :: ilai                                !! counter for loops on LAI levels (-)
87    REAL(r_std), DIMENSION(kjpindex)      :: zqsvegrap
88    REAL(r_std)                           :: speed
89    ! Assimilation
90    LOGICAL, DIMENSION(kjpindex)         :: assimilate                           !! where assimilation is to be calculated (-)
91    LOGICAL, DIMENSION(kjpindex)         :: calculate
92    INTEGER(i_std)                       :: nic,inic,icinic
93    INTEGER(i_std), DIMENSION(kjpindex)  :: index_calc
94    INTEGER(i_std)                       :: nia,inia,nina,inina,iainia
95    INTEGER(i_std), DIMENSION(kjpindex)  :: index_assi,index_non_assi
96
97    REAL(r_std), DIMENSION(kjpindex)      :: vc2                                 !! rate of carboxylation (at a specific LAI level) (umol/m^2/s)
98    REAL(r_std), DIMENSION(kjpindex)      :: vj2                                 !! rate of Rubisco regeneration (at a specific LAI level) (umol/m^2/s)
99    REAL(r_std), DIMENSION(kjpindex)      :: assimi                              !! assimilation (at a specific LAI level) (umol/m^2/s)
100    REAL(r_std)                           :: x_1,x_2,x_3,x_4,x_5,x_6
101    REAL(r_std), DIMENSION(kjpindex)      :: gstop                               !! stomatal conductance at topmost level (m/s)
102    REAL(r_std), DIMENSION(kjpindex)      :: gs                                  !! stomatal conductance (umol/m^2/s)
103    REAL(r_std), DIMENSION(kjpindex)      :: Kc                                  !! Michaelis-Menten constant for the rubisco enzyme catalytic activity for CO2 (ppm)
104    REAL(r_std), DIMENSION(kjpindex)      :: Ko                                  !! Michaelis-Menten constant for the rubisco enzyme catalytic activity for O2 (ppm)
105    REAL(r_std), DIMENSION(kjpindex)      :: CP                                  !! CO2 compensation point (ppm)
106    REAL(r_std), DIMENSION(kjpindex)      :: vc                                  !! rate of carboxylation (umol/m^2/s)
107    REAL(r_std), DIMENSION(kjpindex)      :: vj                                  !! rate of Rubisco regeneration (umol/m^2/s)
108    REAL(r_std), DIMENSION(kjpindex)      :: kt                                  !! temperature-dependent pseudo-first-order rate constant of assimilation response to CO2 (C4) (?)
109    REAL(r_std), DIMENSION(kjpindex)      :: rt                                  !! temperature-dependent non-photosynthetic respiration (C4) (?)
110    REAL(r_std), DIMENSION(kjpindex)      :: air_relhum                          !! air relative humidity (?)
111    REAL(r_std), DIMENSION(kjpindex)      :: water_lim                           !! water limitation factor (fraction)
112    REAL(r_std), DIMENSION(kjpindex)      :: temp_lim                            !! temperature limitation factor (fraction)
113    REAL(r_std), DIMENSION(kjpindex)      :: gstot                               !! total stomatal conductance at the canopy level (m/s)
114    REAL(r_std), DIMENSION(kjpindex)      :: assimtot                            !! total assimilation at the canopy level (umol/m^2/s)
115    REAL(r_std), DIMENSION(kjpindex)      :: leaf_gs_top                         !! stomatal conductance at topmost level (umol/m^2/s)
116    REAL(r_std), DIMENSION(nlai+1)        :: laitab                              !! tabulated LAI steps (m^2/m^2)
117    REAL(r_std), DIMENSION(kjpindex)      :: qsatt                               !! surface saturated humidity (g/g)
118    REAL(r_std), DIMENSION(nvm,nlai)      :: light                               !! fraction of light that gets through (fraction)
119    REAL(r_std), DIMENSION(kjpindex)      :: ci_gs
120    REAL(r_std)                           :: cresist                             !! coefficient for resistances (?)
121
122!! PARAMETER VALUES
123    REAL(r_std), PARAMETER                :: laimax = 12.                        !! maximum LAI (m^2/m^2)
124    REAL(r_std), PARAMETER                :: xc4_1 = .83
125    REAL(r_std), PARAMETER                :: xc4_2 = .93
126
127!> @defgroup Photosynthesis The Mysteries of Photosynthesis
128!> @{   
129    !> 1. Preliminary calculations\n
130    !
131    !> 1.1 Calculate LAI steps\n
132    !> The integration at the canopy level is done over nlai fixed levels.
133    !! \latexonly
134    !! \input{diffuco_trans_co2_1.1.tex}
135    !! \endlatexonly
136!> @}
137!> @codeinc
138    DO jl = 1, nlai+1
139      laitab(jl) = laimax*(EXP(.15*REAL(jl-1,r_std))-1.)/(EXP(.15*REAL(nlai,r_std))-un)
140    ENDDO
141!> @endcodeinc
142
143!> @addtogroup Photosynthesis
144!> @{   
145    !>
146    !> 1.2 Calculate light fraction for each LAI step\n
147    !> The available light follows a simple Beer extinction law.
148    !> The extinction coefficients (ext_coef) are PFT-dependant constants and are defined in constant_co2.f90.
149    !! \latexonly
150    !! \input{diffuco_trans_co2_1.2.tex}
151    !! \endlatexonly
152!> @}
153!> @codeinc
154    DO jl = 1, nlai
155      DO jv = 1, nvm
156        light(jv,jl) = exp( -ext_coef(jv)*laitab(jl) )
157      ENDDO
158    ENDDO 
159!> @endcodeinc
160    !
161    ! Photosynthesis parameters
162    !
163    ! temperatures in K
164    !
165    tmin(:,:) = assim_param(:,:,itmin)
166    tmax(:,:) = assim_param(:,:,itmax)
167    topt(:,:) = assim_param(:,:,itopt)
168    !
169    vcmax(:,:) = assim_param(:,:,ivcmax)
170    vjmax(:,:) = assim_param(:,:,ivjmax)
171
172!> @addtogroup Photosynthesis
173!> @{   
174    !>
175    !> 1.3 Estimate relative humidity of air\n
176    !! \latexonly
177    !! \input{diffuco_trans_co2_1.3.tex}
178    !! \endlatexonly
179!> @}
180    !
181! N. de Noblet - 2006/06 - We use q2m/t2m instead of qair.
182!    CALL qsatcalc (kjpindex, temp_air, pb, qsatt)
183!    air_relhum(:) = &
184!      ( qair(:) * pb(:) / (0.622+qair(:)*0.378) ) / &
185!      ( qsatt(:)*pb(:) / (0.622+qsatt(:)*0.378 ) )
186!> @codeinc
187    CALL qsatcalc (kjpindex, t2m, pb, qsatt)
188    air_relhum(:) = &
189      ( q2m(:) * pb(:) / (0.622+q2m(:)*0.378) ) / &
190      ( qsatt(:)*pb(:) / (0.622+qsatt(:)*0.378 ) )
191!> @endcodeinc
192! N. de Noblet - 2006/06
193    !
194    !> @addtogroup Photosynthesis
195    !> @{   
196    !> 2. Loop over vegetation types\n
197    !> @}
198    !
199    DO jv = 1,nvm
200      !
201      !> @addtogroup Photosynthesis
202      !> @{   
203      !>
204      !> 2.1 Initializations\n
205      !! \latexonly
206      !! \input{diffuco_trans_co2_2.1.tex}
207      !! \endlatexonly
208      !> @}     
209      !
210      ! beta coefficient for vegetation transpiration
211      !
212      rstruct(:,jv) = rstruct_const(jv)
213      rveget(:,jv) = undef_sechiba
214      !
215      vbeta3(:,jv) = zero
216      vbetaco2(:,jv) = zero
217      !
218      cimean(:,jv) = ccanopy(:)
219      !
220      ! mask that contains points where there is photosynthesis
221      ! For the sake of optimization, computations are done only for convenient points.
222      ! nia is the number of points where the assimilation is calculated and nina the number of points where photosynthesis is not calculated
223      ! (based on criteria on minimum or maximum values on LAI, vegetation fraction, shortwave incoming radiation, temperature and relative humidity).
224      ! For the points where assimilation is not calculated, variables are initialized to specific values.
225      ! The assimilate(kjpindex) array contains the logical value (TRUE/FALSE) relative to this photosynthesis calculation.
226      ! The index_assi(kjpindex) array indexes the nia points with assimilation, whereas the index_no_assi(kjpindex) array indexes the nina points with no
227      ! assimilation.
228      nia=0
229      nina=0
230      !
231      DO ji=1,kjpindex
232        !
233        IF ( ( lai(ji,jv) .GT. 0.01 ) .AND. &
234              ( veget_max(ji,jv) .GT. 1.E-8 ) ) THEN
235          IF ( ( veget(ji,jv) .GT. 1.E-8 ) .AND. &
236               ( swdown(ji) .GT. min_sechiba ) .AND. &
237               ( temp_air(ji) .GT. tmin(ji,jv) ) .AND. &
238               ( temp_air(ji) .LT. tmax(ji,jv) ) .AND. &
239               ( humrel(ji,jv) .GT. min_sechiba )       ) THEN
240             !
241             assimilate(ji) = .TRUE.
242             nia=nia+1
243             index_assi(nia)=ji
244             !
245          ELSE
246             !
247             assimilate(ji) = .FALSE.
248             nina=nina+1
249             index_non_assi(nina)=ji
250             !
251          ENDIF
252        ELSE
253          !
254          assimilate(ji) = .FALSE.
255          nina=nina+1
256          index_non_assi(nina)=ji
257          !
258        ENDIF
259        !
260      ENDDO
261      !
262      gstot(:) = zero
263      assimtot(:) = zero
264      !
265      zqsvegrap(:) = zero
266      WHERE (qsintmax(:,jv) .GT. min_sechiba)
267          zqsvegrap(:) = MAX(zero, qsintveg(:,jv) / qsintmax(:,jv))
268      ENDWHERE
269      !
270      WHERE ( assimilate(:) )
271        water_lim(:) = MIN( 2.*humrel(:,jv), un )
272      ENDWHERE
273      ! give a default value of ci for all pixel that do not assimilate
274      DO jl=1,nlai
275         DO inina=1,nina
276            leaf_ci(index_non_assi(inina),jv,jl) = ccanopy(index_non_assi(inina)) * .667_r_std
277         ENDDO
278      ENDDO
279      !
280      ilai(:) = 1
281      !
282      ! Here is calculated photosynthesis (Farqhuar et al., 1980)
283      ! and stomatal conductance (Ball et al., 1987)
284      !
285      ! Calculate temperature dependent parameters
286      !
287      IF ( is_c4(jv) )  THEN
288        !
289        !> @addtogroup Photosynthesis
290        !> @{   
291        !>
292        !> 2.2 Calculate temperature dependent parameters for C4 plants\n
293        !! \latexonly
294        !! \input{diffuco_trans_co2_2.2.tex}
295        !! \endlatexonly
296        !> @}           
297        !
298        IF (nia .GT. 0) then
299          DO inia=1,nia
300            !
301            !> @codeinc
302            x_1 = 0.177 * EXP( 0.069*(temp_air(index_assi(inia))-tp_00) ) 
303            ! = 2.0**(((temp_air(index_assi(inia))-tp_00)-25.0)/10.0)
304            !
305            kt(index_assi(inia)) = 0.7 * x_1 * 1.e6
306            rt(index_assi(inia)) = 0.8 * x_1 / &
307              ( un + EXP(1.3*(temp_air(index_assi(inia))-tmax(index_assi(inia),jv))) )
308            !
309            vc(index_assi(inia)) = vcmax(index_assi(inia),jv) & 
310                * 0.39 * x_1 * water_lim(index_assi(inia)) / &
311!                 * 0.39 * x_1  / &
312                ( (1.0+EXP(0.3*(tmin(index_assi(inia),jv)-temp_air(index_assi(inia))))) &
313                * (1.0+EXP(0.3*(temp_air(index_assi(inia))-topt(index_assi(inia),jv)))) )
314            !> @endcodeinc
315            !
316          ENDDO
317        ENDIF
318        !
319        IF (nina .GT. 0) then
320          ! For the points where assimilation is not calculated, variables are initialized to specific values.
321          DO inina=1,nina
322            kt(index_non_assi(inina)) = zero
323            rt(index_non_assi(inina)) = zero
324            vc(index_non_assi(inina)) = zero
325           ENDDO
326        ENDIF
327        !
328      ELSE
329        !
330        !> @addtogroup Photosynthesis
331        !> @{   
332        !>
333        !> 2.3 Calculate temperature dependent parameters for C3 plants\n
334        !! \latexonly
335        !! \input{diffuco_trans_co2_2.3.tex}
336        !! \endlatexonly
337        !> @}           
338        !
339        IF (nia .GT. 0) then
340          !
341          DO inia=1,nia
342            !
343            !> @codeinc
344            temp_lim(index_assi(inia)) = &
345              (temp_air(index_assi(inia))-tmin(index_assi(inia),jv)) * &
346              (temp_air(index_assi(inia))-tmax(index_assi(inia),jv))
347            temp_lim(index_assi(inia)) = temp_lim(index_assi(inia)) / &
348              (temp_lim(index_assi(inia))-(temp_air(index_assi(inia))-&
349               topt(index_assi(inia),jv))**2)
350            !
351            Kc(index_assi(inia)) = 39.09 * EXP(.085*(temp_air(index_assi(inia))-tp_00))
352            !
353            Ko(index_assi(inia)) = 2.412 * 210000. &
354                * EXP(.085*(temp_air(index_assi(inia))-tmin(index_assi(inia),jv))) / &
355                  (temp_air(index_assi(inia))-tmin(index_assi(inia),jv))
356            !
357            CP(index_assi(inia)) = 42. * EXP( 9.46*(temp_air(index_assi(inia))-tp_00-25.)/&
358                                                           temp_air(index_assi(inia)) )
359            !
360            vc(index_assi(inia)) = vcmax(index_assi(inia),jv) * &
361               temp_lim(index_assi(inia)) * water_lim(index_assi(inia))
362!               temp_lim(index_assi(inia))
363            vj(index_assi(inia)) = vjmax(index_assi(inia),jv) * &
364               temp_lim(index_assi(inia)) * water_lim(index_assi(inia))
365!                temp_lim(index_assi(inia))
366            !> @endcodeinc
367            !
368          ENDDO
369          !
370        ENDIF
371        !
372        IF (nina .GT. 0) then
373          ! For the points where assimilation is not calculated, variables are initialized to specific values.
374          DO inina=1,nina
375            temp_lim(index_non_assi(inina)) = zero
376            Kc(index_non_assi(inina)) = zero
377            Ko(index_non_assi(inina)) = zero
378            CP(index_non_assi(inina)) = zero
379            !
380            vc(index_non_assi(inina)) = zero
381            vj(index_non_assi(inina)) = zero
382          ENDDO
383        ENDIF
384      ENDIF      ! C3/C4
385      !
386      !> @addtogroup Photosynthesis
387      !> @{   
388      !>
389      !> 2.4 Loop over LAI steps to estimate assimilation and conductance\n
390      !> @}           
391      !
392      DO jl = 1, nlai
393        !
394        nic=0
395        calculate(:) = .FALSE.
396        !
397        IF (nia .GT. 0) then
398          DO inia=1,nia
399            calculate(index_assi(inia)) = (laitab(jl) .LE. lai(index_assi(inia),jv) )
400            IF ( calculate(index_assi(inia)) ) THEN
401              nic=nic+1
402              index_calc(nic)=index_assi(inia)
403            ENDIF
404          ENDDO
405        ENDIF
406        !
407        !> @addtogroup Photosynthesis
408        !> @{   
409        !>
410        !> 2.4.1 Vmax is scaled into the canopy due to reduction of nitrogen\n
411        !! \latexonly
412        !! \input{diffuco_trans_co2_2.4.1.tex}
413        !! \endlatexonly
414        !> @}           
415        !
416        x_1 = ( un - .7_r_std * ( un - light(jv,jl) ) )
417        !
418        IF ( nic .GT. 0 ) THEN
419          DO inic=1,nic
420            !> @codeinc
421            vc2(index_calc(inic)) = vc(index_calc(inic)) * x_1
422            vj2(index_calc(inic)) = vj(index_calc(inic)) * x_1
423            !> @endcodeinc
424          ENDDO
425        ENDIF
426        !
427        IF ( is_c4(jv) )  THEN
428          !
429          !> @addtogroup Photosynthesis
430          !> @{   
431          !>
432          !> 2.4.2 Assimilation for C4 plants (Collatz et al., 1991)\n
433          !! \latexonly
434          !! \input{diffuco_trans_co2_2.4.2.tex}
435          !! \endlatexonly
436          !> @}           
437          !
438          DO ji = 1, kjpindex
439            assimi(ji) = zero
440          ENDDO
441          !
442          IF (nic .GT. 0) THEN
443             DO inic=1,nic
444              !> @codeinc
445              x_1 = - ( vc2(index_calc(inic)) + 0.092 * 2.3* swdown(index_calc(inic)) * &
446                    ext_coef(jv) * light(jv,jl) )
447              x_2 = vc2(index_calc(inic)) * 0.092 * 2.3 * swdown(index_calc(inic)) * &
448                    ext_coef(jv) * light(jv,jl) 
449              x_3 = ( -x_1 - sqrt( x_1*x_1 - 4.0 * xc4_1 * x_2 ) ) / (2.0*xc4_1)
450              x_4 = - ( x_3 + kt(index_calc(inic)) * leaf_ci(index_calc(inic),jv,jl) * &
451                    1.0e-6 )
452              x_5 = x_3 * kt(index_calc(inic)) * leaf_ci(index_calc(inic),jv,jl) * 1.0e-6
453              assimi(index_calc(inic)) = ( -x_4 - sqrt( x_4*x_4 - 4. * xc4_2 * x_5 ) ) / (2.*xc4_2)
454              assimi(index_calc(inic)) = assimi(index_calc(inic)) - &
455                                                      rt(index_calc(inic))
456              !> @endcodeinc
457            ENDDO
458          ENDIF
459        ELSE
460          !
461          !> @addtogroup Photosynthesis
462          !> @{   
463          !>
464          !> 2.4.3 Assimilation for C3 plants (Farqhuar et al., 1980)\n
465          !! \latexonly
466          !! \input{diffuco_trans_co2_2.4.3.tex}
467          !! \endlatexonly
468          !> @}           
469          !
470          DO ji = 1, kjpindex
471            assimi(ji) = zero
472          ENDDO
473          !
474          IF (nic .GT. 0) THEN
475            DO inic=1,nic
476              !> @codeinc
477              x_1 = vc2(index_calc(inic)) * leaf_ci(index_calc(inic),jv,jl) / &
478                    ( leaf_ci(index_calc(inic),jv,jl) + Kc(index_calc(inic)) * &
479                    ( un + 210000._r_std / Ko(index_calc(inic)) ) )
480              x_2 = .8855_r_std*swdown(index_calc(inic))*ext_coef(jv)*light(jv,jl)
481              x_3 = x_2+vj2(index_calc(inic))
482              x_4 = ( x_3 - sqrt( x_3*x_3 - (4._r_std*.7_r_std*x_2*vj2(index_calc(inic))) ) ) / &
483                    (2._r_std*.7_r_std)
484              x_5 = x_4 * leaf_ci(index_calc(inic),jv,jl) / &
485                    ( 4.5_r_std *  leaf_ci(index_calc(inic),jv,jl) + &
486                    10.5_r_std*CP(index_calc(inic)) ) 
487              x_6 = MIN( x_1, x_5 )
488              assimi(index_calc(inic)) = x_6 * ( un - CP(index_calc(inic))/&
489                 leaf_ci(index_calc(inic),jv,jl) ) - .011_r_std * vc2(index_calc(inic))
490              !> @endcodeinc
491            ENDDO
492          ENDIF
493        ENDIF
494        !
495        IF (nic .GT. 0) THEN
496          !
497          DO inic=1,nic
498            !
499            !! 2.4.4 Estimate conductance (Ball et al., 1987)
500            !
501            icinic=index_calc(inic)
502!            gs(icinic) = water_lim(icinic) * &
503            gs(icinic) = &
504                 ( gsslope(jv) * assimi(icinic) * &
505                 air_relhum(icinic) / ccanopy(icinic) ) &
506                  + gsoffset(jv)
507            gs(icinic) = MAX( gs(icinic), gsoffset(jv) )
508          ENDDO
509          !
510          DO inic=1,nic
511            icinic=index_calc(inic)
512            !
513            ! the new ci is calculated with
514            ! dci/dt=(ccanopy-ci)*gs/1.6-A
515            ! ci=ci+((ccanopy(icinic)-ci)*gs/1.6-&
516            !    assimi(icinic))*dtradia
517            ! we verify that ci is not out of possible values
518            !
519            ci_gs(icinic) = MIN( ccanopy(icinic), MAX( CP(icinic), &
520                   ( ccanopy(icinic) - 1.6_r_std * assimi(icinic) / &
521                     gs(icinic) ) ) ) - leaf_ci(icinic,jv,jl)
522          ENDDO
523          DO inic=1,nic
524            icinic=index_calc(inic)
525            !to avoid some problem of numerical stability, the leaf_ci is bufferized
526            leaf_ci(icinic,jv,jl) = leaf_ci(icinic,jv,jl) + ci_gs(icinic)/6.
527          ENDDO
528          !
529          DO inic=1,nic
530            icinic=index_calc(inic)
531            !
532            ! this might be the last level for which Ci is calculated. Store it for
533            ! initialization of the remaining levels of the Ci array.
534            !
535            leaf_ci_lowest(icinic) = leaf_ci(icinic,jv,jl)
536          ENDDO
537          !
538          DO inic=1,nic
539            icinic=index_calc(inic)
540            !
541            !! 2.4.5 Integration at the canopy level
542            ! total assimilation and conductance
543            assimtot(icinic) = assimtot(icinic) + &
544              assimi(icinic) * (laitab(jl+1)-laitab(jl))
545            gstot(icinic) = gstot(icinic) + &
546              gs(icinic) * (laitab(jl+1)-laitab(jl))
547            !
548            ilai(icinic) = jl
549            !
550          ENDDO
551          !
552        ENDIF
553        !
554        ! keep stomatal conductance of topmost level
555        !
556        IF ( jl .EQ. 1 ) THEN
557          !
558          leaf_gs_top(:) = zero
559          !
560          IF ( nic .GT. 0 ) then
561            DO inic=1,nic
562              leaf_gs_top(index_calc(inic)) = gs(index_calc(inic))
563             ENDDO
564          ENDIF
565          !
566        ENDIF
567        !
568        IF (nia .GT. 0) THEN
569          !
570          DO inia=1,nia
571            !
572            IF ( .NOT. calculate(index_assi(inia)) ) THEN
573              !
574              !   a) for plants that are doing photosynthesis, but whose LAI is lower
575              !      than the present LAI step, initialize it to the Ci of the lowest
576              !      canopy level
577              !
578              leaf_ci(index_assi(inia),jv,jl) = leaf_ci_lowest(index_assi(inia))
579              !
580            ENDIF
581            !
582          ENDDO
583          !
584        ENDIF
585        !
586      ENDDO  ! loop over LAI steps
587      !
588      !! 2.5 Calculate resistances
589      !
590      IF (nia .GT. 0) THEN
591        !
592        DO inia=1,nia
593          !
594          iainia=index_assi(inia)
595          !
596          ! conversion from mmol/m^2/s to m/s
597          !
598          gstot(iainia) = .0244*(temp_air(iainia)/tp_00)*&
599             (1013./pb(iainia))*gstot(iainia)
600          gstop(iainia) = .0244 * (temp_air(iainia)/tp_00)*&
601             (1013./pb(iainia))*leaf_gs_top(iainia)*&
602              laitab(ilai(iainia)+1)
603          !
604          rveget(iainia,jv) = un/gstop(iainia)
605          !
606        ENDDO
607        !
608        DO inia=1,nia
609          !
610          iainia=index_assi(inia)
611          !
612          ! rstruct is the difference between rtot (=1./gstot) and rveget
613          !
614          ! Correction Nathalie - le 27 Mars 2006 - Interdire a rstruct d'etre negatif
615          !rstruct(iainia,jv) = un/gstot(iainia) - &
616          !     rveget(iainia,jv)
617          rstruct(iainia,jv) = MAX( un/gstot(iainia) - &
618               rveget(iainia,jv), min_sechiba)
619          !
620        ENDDO
621        !
622        DO inia=1,nia
623          !
624          iainia=index_assi(inia)
625          !
626          speed = MAX(min_wind, wind(index_assi(inia)))
627          !
628          ! beta for transpiration
629          !
630          ! Corrections Nathalie - le 28 mars 2006 - sur conseils Fred Hourdin
631          ! Introduction d'un potentiometre (rveg_pft) pour regler la somme rveg+rstruct
632          !vbeta3(iainia,jv) = veget_max(iainia,jv) * &
633          !  (un - zqsvegrap(iainia)) * &
634          !  (un / (un + speed * q_cdrag(iainia) * (rveget(iainia,jv) + &
635          !   rstruct(iainia,jv))))
636          cresist=(un / (un + speed * q_cdrag(iainia) * &
637               (rveg_pft(jv)*(rveget(iainia,jv) + rstruct(iainia,jv)))))
638
639          vbeta3(iainia,jv) = veget_max(iainia,jv) * &
640            (un - zqsvegrap(iainia)) * cresist + &
641!!$          ! Ajout Nathalie - Juin 2006
642!!$          vbeta3(iainia,jv) = vbeta3(iainia,jv) + &
643          ! Corrections Nathalie - le 09 novembre 2009 : veget => veget_max
644!               MIN( vbeta23(iainia,jv), veget(iainia,jv) * &
645               MIN( vbeta23(iainia,jv), veget_max(iainia,jv) * &
646!               zqsvegrap(iainia) * humrel(iainia,jv) * &
647               zqsvegrap(iainia) * cresist )
648          ! Fin ajout Nathalie
649          !
650          ! beta for assimilation. The difference is that surface
651          ! covered by rain (un - zqsvegrap(iainia)) is not taken into account
652          ! 1.6 is conversion for H2O to CO2 conductance
653          ! vbetaco2(iainia,jv) = veget_max(iainia,jv) * &
654          !   (un / (un + q_cdrag(iainia) * &
655          !   (rveget(iainia,jv))))/1.6
656          !
657          vbetaco2(iainia,jv) = veget_max(iainia,jv) * &
658             (un / (un + speed * q_cdrag(iainia) * &
659             (rveget(iainia,jv) + rstruct(iainia,jv)))) / 1.6
660          !
661          ! cimean is the "mean ci" calculated in such a way that assimilation
662          ! calculated in enerbil is equivalent to assimtot
663          !
664          cimean(iainia,jv) = ccanopy(iainia) - &
665             assimtot(iainia) / &
666             ( vbetaco2(iainia,jv)/veget_max(iainia,jv) * &
667             rau(iainia) * speed * q_cdrag(iainia))
668          !
669        ENDDO
670        !
671      ENDIF
672      !
673    END DO         ! loop over vegetation types
674    !
675    IF (long_print) WRITE (numout,*) ' diffuco_trans_co2 done '
676
677
678END