Context Navigation

Back to Documentation/ORCHIDEE_DOC

Documentation/ORCHIDEE_DOC: diffuco_trans_co2.f90

File diffuco_trans_co2.f90, 28.1 KB (added by maignan, 12 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	!! ==============================================================================
37	SUBROUTINE 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(.15REAL(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_1x_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_4x_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_stdswdown(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_3x_3 - (4._r_std.7_r_stdx_2vj2(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
678	END

Download in other formats:

Original Format